先说下数据,随访流水号是患者的后续诊断号码,表3有对应的数据,首先需要做下数据整理,需要整理出每次诊断的指标(包括表1中人物信息、表2中的检查指标以及表3中的检查指标,表4中有对应的时间,以刚发病时为0时刻,基于流水号进行匹配整合)
这个题目主要涉及到机器学习算法和常规的评价算法,首先可以检查下数据是否有异常或缺失情况,采取必要处理,缺数据的是表3,可以通过求当前流水号的其他指标数据与集合中相同指标数据的差距,从而找到近似样本,对缺失数据进行补充,这一步骤和协同过滤算法类似。
第一问:
第一小问:是标记患者是否发生血肿扩张,看附件2中HM_volume的几次指标数值的变化是否出现体积增加≥6 mL或相对体积增加≥33%的情况,是则标记为1,流水号用附表1查下时间,以首次流水号上对应的时间减去发病到首次检测的间隔时间作为0时刻,依次计算出后续流水号对应的时刻,对时间和HM_volume指标值做插值拟合或者拟合关系式,并算出48小时的HM_volume值,然后判断是否发生血肿扩张,如果48小时内判断发生血肿扩张,则需要计算出刚好满足判断为血肿扩张的时间,注意HM_volume值单位是10^-3ml。
第二小问:基于患者首次检测的指标作为输入,是否48小时内是否发生水肿扩张为输出,用机器学习算法训练,训练集测试集划分见题目说明,模型训练好之后分别带入这几个数据集,算出对应的概率,这个概率别四舍五入了,一些算法最后输出的是整数,但是输出结构中能调出对应的概率。
第二问:
第一小问:建立时间与ED_volume值的关系模型,建议用高斯模型,并进行拟合检验,需要拿前100个患者的检测时间和ED_volume值数据来做分析,然后分别带入前100个患者的发病时间去计算下拟合值,并计算出残差平均值。
第二小问:先聚类,类别3-5个,然后分别按照第一小问的方式做就行。
第三小问:首先提取有5次检测以上的样本计算后一时刻与前一时刻的变化量,这里需要治疗对患者水肿的改善程度,只保留为负值的改变量,如果都为正值,那么就视为没有改善,是否选择不同的治疗作为两个组别,遍历每种治疗进行单因素方差分析,然后讨论不同治疗对水肿的进展模式的影响显著性。
第四小问:同样的方法对血肿进行分析,然后再分析下血肿指标与水肿之间的相关性。
第三问:
第一小问:建立指标与mRs的评分模型,还是采用机器学习算法做训练,然后预测,做训练前,可以做下主要指标提取的步骤,只预测101-160号患者首次检测的mRs值,训练集和测试集划分还是根据题目的来。
第二小问:预测随访的mRs值,患者的基础信息匹配附件1,连同检测指标,带入模型中进行预测。
第三小问:做下简单的相关性分析,这种医学数据直接用相关性算法肯定结果值很低,建议使用方差分析、t检验、F检验这类算法作显著性分析,p值越小显著性越大,确定显著影响关系后,给一些建议即可。
代码(用的简单方法实现,大概就这么多内容,请自行更换算法,适当包装下)
%预处理
clear
clc
%读取数据
[~,~,data1]=xlsread('表1-患者列表及临床信息.xlsx','A1:W161');
data1=string(data1);
data1(find(ismissing(data1)==1))="";
data1=[data1(:,1:15),[["血压最大值","血压最小值"];split(data1(2:end,16),"/")],data1(:,17:end)];
data1(find(data1=="男"))="1";
data1(find(data1=="女"))="2";
[~,~,data2]=xlsread('表2-患者影像信息血肿及水肿的体积及位置.xlsx','A1:GZ161');
data2=string(data2);
data2(find(ismissing(data2)==1))="";
[~,~,data3]=xlsread('表3-患者影像信息血肿及水肿的形状及灰度分布.xlsx','B1:AG577');
data3=string(data3);
data3(find(ismissing(data3)==1))="";
[~,~,data4]=xlsread('附表1-检索表格-流水号vs时间.xlsx','C2:AB161');
data4=string(data4);
data4(find(ismissing(data4)==1))="";
%整理数据
Table=["患者","流水号","时间/小时",data1(1,[5:14,16:24]),data2(1,3:24),data3(1,2:32)];%记录指标数据
for i=1:size(data1,1)-1
delta_t0=double(data1(i+1,15));%从发病到第一次诊断时的时间长度
a=find(data4==data1(i+1,4));
a1=mod(a,160);
if a1==0
a1=160;
end
a2=ceil(a/160);
t1=data4(a1,a2-1);%第一次诊断日期
table1=data1(i+1,[5:14,16:24]);
table2=data2(i+1,2:208);
%识别检查了多少次
a=find(table2=="");
if length(a)==0
a=length(table2);
else
a=a(1);
end
n=(a-1)/23;
%统计每次检查的数据,如果有缺失则不记录
for j=1:n
a=find(data4==table2((j-1)*23+1));
a1=mod(a,160);
if a1==0
a1=160;
end
a2=ceil(a/160);
t2=data4(a1,a2-1);%后续诊断日期
time=max((datenum(t2)-datenum(t1))*24+delta_t0,0);%以刚发病时刻为0,计算时间轴上第j次随诊距离发病的时间长度
if length(find(data3(:,1)==table2((j-1)*23+1))>0)
x=[data1(i+1,1),table2((j-1)*23+1),time,table1,table2((j-1)*23+2:j*23),data3(find(data3(:,1)==table2((j-1)*23+1)),2:end)];
else%如果缺失则取样本差距前三个数据集的附件3指标取平均值补充
x1=[double(table2((j-1)*23+2:j*23));double(Table(2:end,23:44))];
x1=mapminmax(x1',0,1)';%数据标准化
d=pdist2(x1(1,:),x1(2:end,:));
[~,o]=sort(d);
o=o(1:min(3,length(o)));
x=[data1(i+1,1),table2((j-1)*23+1),time,table1,table2((j-1)*23+2:j*23),mean(double(Table((o+1),45:75)),1)];
end
Table=[Table;x];
end
end
save data Table
%第一问
clear
clc
warning off
load data
%% a
Num=unique(Table(2:end,1),'stable');
A=[];%患者第一次诊断时的指标,问题b也会用到
B=[];%记录拟合参数
Y=[];%是否发病
T=zeros(length(Num),1);%静脉扩张时间
for i=1:length(Num)
a=find(Table(:,1)==Num(i));
A=[A;double(Table(a(1),4:end))];
t=double(Table(a,3));
HM_volume=double(Table(a,23));
if length(HM_volume>0)
b=regress(HM_volume,[ones(length(t),1),t]);%线性回归,可以自己换算法
else%如果只有一个数据,则用样本数据匹配最接近的其他样本对应的参数
AA=mapminmax(A',0,1)';
d=pdist2(AA(end,:),AA(1:end-1,:));
[~,o]=min(d);
b=B(:,o);
end
B=[B,b];
t1=[t(1):0.1:48]';%将时间网格化,精度为0.1
HM_volume_48=[ones(length(t1),1),t1]*b;%时间设为第一次诊断到发病第48小时内
a1=find(HM_volume_48./HM_volume(1)>1.33);
a2=find(HM_volume_48-HM_volume(1)>6000);
aa=union(a1,a2);
if length(aa)>0
Y(i,1)=1;
T(i,1)=t1(aa(1));
else
Y(i,1)=0;
end
end
result1=[Num,Y,T];
%% b
%变量
In=mapminmax(A',0,1);
Out=Y';
%bp神经网络
Xtrain = In(:,1:100);
Ytrain = Out(:,1:100);
Xtest1 = In(:,101:130);
Ytest1= Out(:,101:130);
Xtest2 = In(:,131:160);
Ytest2= Out(:,131:160);
% 1. 创建网络
net = newff(Xtrain,Ytrain,[fix(size(Xtrain,1)/2),fix(size(Xtrain,1)/4)],{'tansig','tansig'});
% 2. 设置训练参数
net.trainParam.epochs = 1000;%最大迭代次数,到达最大迭代次数则终止
net.trainParam.goal = 1e-100;%训练误差,达到目标误差则终止
net.trainParam.min_grad = 1e-100;%性能函数的最小梯度
net.trainParam.lr = 1e-5;%学习速率
net.trainParam.max_fail=100;%最大确认失败次数,终止条件之一
% 3. 训练网络
net = train(net,Xtrain,Ytrain);
% 4. 仿真测试
t_sim = max(sim(net,Xtrain),0);
t_sim1 = max(sim(net,Xtest1),0);
t_sim2 = max(sim(net,Xtest2),0);
T_sim = round(t_sim);
T_sim1 = round(t_sim1);
T_sim2 = round(t_sim2);
%混淆矩阵
figure
plotconfusion(categorical(Ytrain),categorical(T_sim));
title('训练集')
figure
plotconfusion(categorical(Ytest1),categorical(T_sim1));
title('测试集1')
figure
plotconfusion(categorical(Ytest2),categorical(T_sim2));
title('测试集2')
result1=[["患者","是否发生血肿扩张","血肿扩张时间","血肿扩张预测概率"];[result1,[t_sim';t_sim1';t_sim2']]];
%第二问
clear
clc
%% a
load data
Num=unique(Table(2:end,1),'stable');
T=[];%记录复查时间
ED=[];%记录水肿数据
A=[];%记录除此诊断信息
for i=1:length(Num)
a=find(Table(:,1)==Num(i));
A=[A;double(Table(a(1),4:end))];
T=[T;double(Table(a,3))];
ED=[ED;double(Table(a,34))];
end
figure
plot(T,ED,'*')
xlim([1,2000])
%高斯模型
hold on
gaussModel = fit(T, ED, 'gauss1')
plot(gaussModel)
xlabel('时间')
ylabel('水肿/10^-3ml')
%计算残差
ED_fit=gaussModel(T);
Error=[];
for i=1:length(Num)
a=find(Table(:,1)==Num(i))-1;
Error(i,1)=mean(abs(ED_fit(a)-ED(a)));
end
%% b
%个体差异用Tabel第4-15列的指标来进行聚类
cluster_n=5;%聚类中心
[center, U, obj_fcn] = fcm(double(A(:,1:12)), cluster_n);
figure%目标函数变化值
plot(obj_fcn)
xlabel('iteration')
ylabel('obj.fcn_value')
title('FCM聚类')
[~,u]=max(U);
ED_fit2=[];
for i=1:cluster_n
a=find(u==i);
b=find(ismember(Table(2:end,1),Num(a))==1);
figure
plot(T(b),ED(b),'*')
xlim([1,2000])
%高斯模型
hold on
disp('亚类1:')
gaussModel = fit(T(b), ED(b), 'gauss1')
plot(gaussModel)
xlabel('时间')
ylabel('水肿/10^-3ml')
title(['亚类',num2str(i)])
ED_fit2(b,1)=gaussModel(T(b));
end
%计算残差
Error2=[];
for i=1:length(Num)
a=find(Table(:,1)==Num(i))-1;
Error2(i,1)=mean(abs(ED_fit2(a)-ED(a)));
end
result2=[["患者","残差(全体)","残差(亚类)","亚类"];[Num,Error,Error2,u']];
%% c
%计算水肿指标的变化率,在不同治疗方法下,改善效率
K=[];
for i=1:length(Num)
a=find(Table(:,1)==Num(i));
if length(a)>=5%至少五次检查数据才进行计算
k=(double(Table(a(2:end),34))-double(Table(a(1:end-1),34)))./(double(Table(a(2:end),3))-double(Table(a(1:end-1),3)));
kk=find(k<0);
if length(kk)==0
K(i,1)=0;
else
K(i,1)=mean(k(kk));
end
else
K(i,1)=NaN;
end
end
c=setdiff([1:length(K)],find(isnan(K)==1));
G=double(Table(2:end,16:22));
Z=Table(1,16:22);
for i=1:length(Z)
[p(i),anovatab{i},stats{i}]=anova1(K(c),G(c,i)+1,'off');%单因素方差分析
fa=finv(0.95,anovatab{i}{2,3},anovatab{i}{3,3});%计算fa
F=anovatab{i}{2,5};%F值
if p(i)<=0.01 && F>fa
disp([Z(i)+"对水肿体积进展模式影响非常显著"])
fprintf('p值为%.4f<0.01,F值为%.2f>%.2f\n',p(i),F,fa)
elseif p(i)<=0.05 && F>fa
disp([Z(i)+"对水肿体积进展模式影响显著"])
fprintf('p值为%.4f<0.05,F值为%.2f>%.2f\n',p(i),F,fa)
else
disp([Z(i)+"对水肿体积进展模式影响不显著"])
fprintf('p值为%.4f,F值为%.2f\n',p(i),F)
end
end
disp('不同治疗对水肿进展模式的影响大小为:')
[~,q]=sort(p);
str=Z(q)+">";
disp(str)
%% d
%同上述步骤求血肿体积的
%计算水肿指标的变化率,在不同治疗方法下,改善效率
K2=[];
for i=1:length(Num)
a=find(Table(:,1)==Num(i));
if length(a)>=5%至少五次检查数据才进行计算
k=(double(Table(a(2:end),23))-double(Table(a(1:end-1),23)))./(double(Table(a(2:end),3))-double(Table(a(1:end-1),3)));
kk=find(k<0);
if length(kk)==0
K2(i,1)=0;
else
K2(i,1)=mean(k(kk));
end
else
K2(i,1)=NaN;
end
end
c=setdiff([1:length(K2)],find(isnan(K2)==1));
for i=1:length(Z)
[p2(i),anovatab2{i},stats2{i}]=anova1(K2(c),G(c,i)+1,'off');%单因素方差分析
fa=finv(0.95,anovatab2{i}{2,3},anovatab2{i}{3,3});%计算fa
F=anovatab2{i}{2,5};%F值
if p2(i)<=0.01 && F>fa
disp([Z(i)+"对血肿肿体积进展模式影响非常显著"])
fprintf('p值为%.4f<0.01,F值为%.2f>%.2f\n',p2(i),F,fa)
elseif p2(i)<=0.05 && F>fa
disp([Z(i)+"对血肿体积进展模式影响显著"])
fprintf('p值为%.4f<0.05,F值为%.2f>%.2f\n',p2(i),F,fa)
else
disp([Z(i)+"对血肿体积进展模式影响不显著"])
fprintf('p值为%.4f,F值为%.2f\n',p2(i),F)
end
end
disp('不同治疗对血肿进展模式的影响大小为:')
[~,q2]=sort(p2);
str2=Z(q2)+">";
disp(str2)
%在算下血肿指标与水肿指标的相关性
x0=double(Table(2:end,23));
y0=double(Table(2:end,34));
theta=x0'*y0/(norm(x0)*norm(y0));%余弦相似度
fprintf('血肿指标与水肿指标的相关度为:%.4f\n',theta)
clear
clc
load data
Num=unique(Table(2:end,1),'stable');
A=[];%首检
Y=[];%首检Rms
B=[];%随检
N1=[];
N2=[];
Z=[Table(1,[4:5,7:end])];
for i=1:length(Num)
a=find(Table(:,1)==Num(i));
A=[A;double(Table(a(1),[4:5,7:end]))];
Y=[Y;double(Table(a(1),6))];
N1=[N1;Table(a(1),1:2)];
B=[B;double(Table(a(2:end),[4:5,7:end]))];
N2=[N2;Table(a(2:end),1:2)];
end
%% a
[In,ps_in]=mapminmax(A',0,1);
Out=Y';
%bp神经网络
Xtrain = In(:,1:100);
Ytrain = Out(:,1:100);
Xtest1 = In(:,101:130);
Ytest1= Out(:,101:130);
Xtest2 = In(:,131:160);
Ytest2= Out(:,131:160);
% 1. 创建网络
net = newff(Xtrain,Ytrain,[fix(size(Xtrain,1)/2),fix(size(Xtrain,1)/4)],{'tansig','tansig'});
% 2. 设置训练参数
net.trainParam.epochs = 1000;%最大迭代次数,到达最大迭代次数则终止
net.trainParam.goal = 1e-100;%训练误差,达到目标误差则终止
net.trainParam.min_grad = 1e-100;%性能函数的最小梯度
net.trainParam.lr = 1e-5;%学习速率
net.trainParam.max_fail=100;%最大确认失败次数,终止条件之一
% 3. 训练网络
net = train(net,Xtrain,Ytrain);
% 4. 仿真测试
t_sim = max(sim(net,Xtrain),0);
t_sim1 = max(sim(net,Xtest1),0);
t_sim2 = max(sim(net,Xtest2),0);
T_sim = round(t_sim);
T_sim1 = round(t_sim1);
T_sim2 = round(t_sim2);
%混淆矩阵
figure
plotconfusion(categorical(Ytrain),categorical(T_sim));
title('训练集')
figure
plotconfusion(categorical(Ytest1),categorical(T_sim1));
title('测试集1')
figure
plotconfusion(categorical(Ytest2),categorical(T_sim2));
title('测试集2')
%% b
Xyuce = mapminmax(B','apply',ps_in);
t_sim3 = max(sim(net,Xyuce),0);
T_sim3 = round(t_sim3);
result3="";
for i=1:length(Num)
a=find(N2==Num(i));
for j=1:length(a)
result3(i,j)=string(T_sim3(a(j)));
end
end
result3=string(result3);
result3(ismissing(result3)==1)="";
z=[];
for i=1:size(result3,2)
z=[z,"第"+num2str(i)+"次随访"];
end
result3=[["患者","首次检测";Num,[T_sim';T_sim1';T_sim2']],[z;result3]];
%% c
%相关性分析
X=double(Table(2:end,[4:5,7:end]));
P=[];
for i=1:size(X,2)
for j=1:size(X,2)
P(i,j)=X(:,i)'*X(:,j)/(norm(X(:,i))*norm(X(:,j)));
end
end
figure%热图
set(gca,'position',[0.05 0.05 0.9 0.9])
h=heatmap(P,'ColorbarVisible', 'on');
h.FontSize = 8;
h.CellLabelFormat = '%0.2g';
resultp=["",Z;Z',P];%自行分析相关性结果
disp('预测结果见result3,相关性结果见resultp')