• 出血性脑卒中临床智能诊疗建模

    先说下数据，随访流水号是患者的后续诊断号码，表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值越小显著性越大，确定显著影响关系后，给一些建议即可。

    代码（用的简单方法实现，大概就这么多内容，请自行更换算法，适当包装下）

    1. clear
    2. clc
    3.  
    4. %读取数据
    5. [~,~,data1]=xlsread('表1-患者列表及临床信息.xlsx','A1:W161');
    6. data1=string(data1);
    7. data1(find(ismissing(data1)==1))="";
    8. data1=[data1(:,1:15),[["血压最大值","血压最小值"];split(data1(2:end,16),"/")],data1(:,17:end)];
    9. data1(find(data1=="男"))="1";
    10. data1(find(data1=="女"))="2";
    11. [~,~,data2]=xlsread('表2-患者影像信息血肿及水肿的体积及位置.xlsx','A1:GZ161');
    12. data2=string(data2);
    13. data2(find(ismissing(data2)==1))="";
    14. [~,~,data31]=xlsread('表3-患者影像信息血肿及水肿的形状及灰度分布.xlsx','ED','B1:AG577');
    15. data31=string(data31);
    16. data31(find(ismissing(data31)==1))="";
    17. [~,~,data32]=xlsread('表3-患者影像信息血肿及水肿的形状及灰度分布.xlsx','Hemo','B1:AG546');
    18. data32=string(data32);
    19. data32(find(ismissing(data32)==1))="";
    20. [~,~,data4]=xlsread('附表1-检索表格-流水号vs时间.xlsx','C2:AB161');
    21. data4=string(data4);
    22. data4(find(ismissing(data4)==1))="";
    23.  
    24. %整理数据（这里将表1-3的信息全部整理，有缺失的取近似的样本取平均进行补充，每次检测的时间以每个患者发病时的时刻为0时刻进行计算）
    25. Table=["患者","流水号","时间/小时",data1(1,[5:14,16:24]),data2(1,3:24),"ED_"+data31(1,2:32),"Hemo_"+data32(1,2:32)];%记录指标数据
    26. for i=1:size(data1,1)-1
    27.     delta_t0=double(data1(i+1,15));%从发病到第一次诊断时的时间长度
    28.     a=find(data4==data1(i+1,4));
    29.     a1=mod(a,160);
    30.     if a1==0
    31.         a1=160;
    32.     end
    33.     a2=ceil(a/160);
    34.     t1=data4(a1,a2-1);%第一次诊断日期
    35.     table1=data1(i+1,[5:14,16:24]);
    36.     table2=data2(i+1,2:208);
    37.     %识别检查了多少次
    38.     a=find(table2=="");
    39.     if length(a)==0
    40.         a=length(table2);
    41.     else
    42.         a=a(1);
    43.     end
    44.     n=(a-1)/23;
    45.     %统计每次检查的数据，如果有缺失则不记录
    46.     for j=1:n
    47.         a=find(data4==table2((j-1)*23+1));
    48.         a1=mod(a,160);
    49.         if a1==0
    50.             a1=160;
    51.         end
    52.         a2=ceil(a/160);
    53.         t2=data4(a1,a2-1);%后续诊断日期
    54.         time=max((datenum(t2)-datenum(t1))*24+delta_t0,0);%以刚发病时刻为0，计算时间轴上第j次随诊距离发病的时间长度
    55.         if length(find(data31(:,1)==table2((j-1)*23+1))>0)
    56.             x1=[data1(i+1,1),table2((j-1)*23+1),time,table1,table2((j-1)*23+2:j*23),data31(find(data31(:,1)==table2((j-1)*23+1)),2:end)];
    57.         else%如果缺失则取样本差距前三个数据集的附件3指标取平均值补充
    58.             x1=[double(table2((j-1)*23+2:j*23));double(Table(2:end,23:44))];
    59.             x1=mapminmax(x1',0,1)';%数据标准化
    60.             d=pdist2(x1(1,:),x1(2:end,:));
    61.             [~,o]=sort(d);
    62.             o=o(1:min(3,length(o)));
    63.             x1=[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)];
    64.         end
    65.         if length(find(data32(:,1)==table2((j-1)*23+1))>0)
    66.             x2=[data32(find(data32(:,1)==table2((j-1)*23+1)),2:end)];
    67.         else%如果缺失则取样本差距前三个数据集的附件3指标取平均值补充
    68.             x2=[double(table2((j-1)*23+2:j*23));double(Table(2:end,23:44))];
    69.             x2=mapminmax(x2',0,1)';%数据标准化
    70.             d=pdist2(x2(1,:),x2(2:end,:));
    71.             [~,o]=sort(d);
    72.             o=o(1:min(3,length(o)));
    73.             x2=[mean(double(Table((o+1),45:75)),1)];
    74.         end
    75.         Table=[Table;x1,x2];
    76.     end
    77. end
    78.  
    79. save data Table
    1. %第一问
    2. clear
    3. clc
    4. warning off
    5. load data
    6.  
    7. %% a 是否要进行指标降维自行添加，降维后肯定能减少运算量
    8. Num=unique(Table(2:end,1),'stable');
    9. A=[];%患者第一次诊断时的指标，问题b也会用到
    10. B=[];%记录拟合参数
    11. Y=[];%是否发病
    12. T=zeros(length(Num),1);%静脉扩张时间
    13. for i=1:length(Num)
    14.     a=find(Table(:,1)==Num(i));
    15.     A=[A;double(Table(a(1),4:end))];
    16.     t=double(Table(a,3));
    17.     HM_volume=double(Table(a,23));
    18.     if length(HM_volume>0)
    19.         b=regress(HM_volume,[ones(length(t),1),t]);%线性回归，可以自己换算法
    20.     else%如果只有一个数据，则用样本数据匹配最接近的其他样本对应的参数
    21.         AA=mapminmax(A',0,1)';
    22.         d=pdist2(AA(end,:),AA(1:end-1,:));
    23.         [~,o]=min(d);
    24.         b=B(:,o);
    25.     end
    26.     B=[B,b];
    27.     t1=[t(1):0.1:48]';%将时间网格化，精度为0.1
    28.     HM_volume_48=[ones(length(t1),1),t1]*b;%时间设为第一次诊断到发病第48小时内
    29.     a1=find(HM_volume_48./HM_volume(1)>1.33);
    30.     a2=find(HM_volume_48-HM_volume(1)>6000);
    31.     aa=union(a1,a2);
    32.     if length(aa)>0
    33.         Y(i,1)=1;
    34.         T(i,1)=t1(aa(1));
    35.     else
    36.         Y(i,1)=0;
    37.     end
    38. end
    39. result1=[Num,Y,T];
    40. %% b 如果要用神经网络算法，可以加一个优化算法对网络结构的寻优
    41. %变量
    42. In=mapminmax(A',0,1);
    43. Out=Y';
    44. %bp神经网络
    45. Xtrain = In(:,1:100);
    46. Ytrain = Out(:,1:100);
    47. Xtest1 = In(:,101:130);
    48. Ytest1= Out(:,101:130);
    49. Xtest2 = In(:,131:160);
    50. Ytest2= Out(:,131:160);
    51. % 1. 创建网络
    52. net = newff(Xtrain,Ytrain,[fix(size(Xtrain,1)/2),fix(size(Xtrain,1)/4)],{'tansig','tansig'});
    53. % 2. 设置训练参数
    54. net.trainParam.epochs = 1000;%最大迭代次数，到达最大迭代次数则终止
    55. net.trainParam.goal = 1e-100;%训练误差，达到目标误差则终止
    56. net.trainParam.min_grad = 1e-100;%性能函数的最小梯度
    57. net.trainParam.lr = 1e-5;%学习速率
    58. net.trainParam.max_fail=100;%最大确认失败次数，终止条件之一
    59. % 3. 训练网络
    60. net = train(net,Xtrain,Ytrain);
    61. % 4. 仿真测试
    62. t_sim = max(sim(net,Xtrain),0);
    63. t_sim1 = max(sim(net,Xtest1),0);
    64. t_sim2 = max(sim(net,Xtest2),0);
    65. T_sim = round(t_sim);
    66. T_sim1 = round(t_sim1);
    67. T_sim2 = round(t_sim2);
    68. %混淆矩阵
    69. figure
    70. plotconfusion(categorical(Ytrain),categorical(T_sim));
    71. title('训练集')
    72. figure
    73. plotconfusion(categorical(Ytest1),categorical(T_sim1));
    74. title('测试集1')
    75. figure
    76. plotconfusion(categorical(Ytest2),categorical(T_sim2));
    77. title('测试集2')
    78. figure
    79. hold on
    80. auc=plot_roc(T_sim,Ytrain);
    81. auc1=plot_roc(T_sim1,Ytest1);
    82. auc2=plot_roc(T_sim2,Ytest2);
    83. legend(['训练集 auc=',num2str(round(auc,2))],['测试集1 auc=',num2str(round(auc1,2))],['测试集2 auc=',num2str(round(auc2,2))])
    84.  
    85. result1=[["患者","是否发生血肿扩张","血肿扩张时间","血肿扩张预测概率"];[result1,[t_sim';t_sim1';t_sim2']]];
    86. disp('结果见矩阵result1')
    1. %绘制roc曲线
    2. function auc=plot_roc(predict,ground_truth)
    3.  
    4. %根据该数目可以计算出沿x轴或者y轴的步长
    5. x_step = 1.0/length(predict);
    6. y_step = 1.0/length(predict);
    7. %首先对ground_truth中的分类器输出值按照从小到大排列
    8. [ground_truth,index] = sort(ground_truth);
    9. predict = predict(:,index);
    10.  
    11. for ii=1:size(predict,1)
    12.     %对predict中的每个样本分别判断他们是FP或者是TP
    13.     %遍历ground_truth的元素，
    14.     %若ground_truth[i]=1,则TP增加了1，往y轴方向上升y_step
    15.     %若ground_truth[i]=0,则FP增加了1，往x轴方向上升x_step
    16.     %初始点为(0.0,0.0)
    17.     x = 0.0;
    18.     y = 0.0;
    19.     X=[0];
    20.     Y=[0];
    21.     for i=1:length(ground_truth)
    22.         if ground_truth(i) == predict(ii,i)
    23.             y = y + y_step;
    24.         else
    25.             x = x + x_step;
    26.         end
    27.         X(i+1)=x;
    28.         Y(i+1)=y;
    29.     end
    30.     X(end+1)=x;
    31.     Y(end+1)=1;
    32.     X(end+1)=1;
    33.     Y(end+1)=1;
    34. end
    35. %画出图像     
    36. plot(X,Y,'-','LineWidth',1);
    37. %trapz，matlab自带函数，计算小矩形的面积,返回auc
    38. auc(ii) = trapz(X,Y);        
    1. %第二问
    2. clear
    3. clc
    4.  
    5. %% a
    6. load data
    7. Num=unique(Table(2:end,1),'stable');
    8. T=[];%记录复查时间
    9. ED=[];%记录水肿数据
    10. A=[];%记录除此诊断信息
    11. for i=1:length(Num)
    12.     a=find(Table(:,1)==Num(i));
    13.     A=[A;double(Table(a(1),4:end))];
    14.     T=[T;double(Table(a,3))];
    15.     ED=[ED;double(Table(a,34))];
    16. end
    17. figure
    18. plot(T,ED,'*')
    19. xlim([1,2000])
    20. %高斯模型
    21. hold on
    22. gaussModel = fit(T, ED, 'gauss1')
    23. plot(gaussModel)
    24. xlabel('时间')
    25. ylabel('水肿/10^-3ml')
    26. %计算残差
    27. ED_fit=gaussModel(T);
    28. Error=[];
    29. for i=1:length(Num)
    30.     a=find(Table(:,1)==Num(i))-1;
    31.     Error(i,1)=mean(abs(ED_fit(a)-ED(a)));
    32. end
    33. %% b
    34. %个体差异用Tabel第4-15列的指标来进行聚类
    35. cluster_n=5;%聚类中心
    36. [center, U, obj_fcn] = fcm(double(A(:,1:12)), cluster_n);
    37. figure%目标函数变化值
    38. plot(obj_fcn)
    39. xlabel('iteration')
    40. ylabel('obj.fcn_value')
    41. title('FCM聚类')
    42. [~,u]=max(U);%所属亚类
    43. ED_fit2=[];
    44. for i=1:cluster_n
    45.     a=find(u==i);
    46.     b=find(ismember(Table(2:end,1),Num(a))==1);
    47.     figure
    48.     plot(T(b),ED(b),'*')
    49.     xlim([1,2000])
    50.     %高斯模型
    51.     hold on
    52.     disp('亚类1：')
    53.     gaussModel = fit(T(b), ED(b), 'gauss1')
    54.     plot(gaussModel)
    55.     xlabel('时间')
    56.     ylabel('水肿/10^-3ml')
    57.     title(['亚类',num2str(i)])
    58.     ED_fit2(b,1)=gaussModel(T(b));
    59. end
    60. %计算残差
    61. Error2=[];
    62. for i=1:length(Num)
    63.     a=find(Table(:,1)==Num(i))-1;
    64.     Error2(i,1)=mean(abs(ED_fit2(a)-ED(a)));
    65. end
    66. result2=[["患者","残差(全体)","残差(亚类)","亚类"];[Num,Error,Error2,u']];
    67. %% c
    68. %计算水肿指标的变化率，在不同治疗方法下，改善效率
    69. K=[];
    70. for i=1:length(Num)
    71.     a=find(Table(:,1)==Num(i));
    72.     if length(a)>=5%至少五次检查数据才进行计算
    73.         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)));
    74.         kk=find(k<0);
    75.         if length(kk)==0
    76.             K(i,1)=0;
    77.         else
    78.             K(i,1)=mean(k(kk));
    79.         end
    80.     else
    81.         K(i,1)=NaN;
    82.     end
    83. end
    84. c=setdiff([1:length(K)],find(isnan(K)==1));
    85. G=double(Table(2:end,16:22));
    86. Z=Table(1,16:22);
    87. for i=1:length(Z)
    88.     [p(i),anovatab{i},stats{i}]=anova1(K(c),G(c,i)+1,'off');%单因素方差分析
    89.     fa=finv(0.95,anovatab{i}{2,3},anovatab{i}{3,3});%计算fa
    90.     F=anovatab{i}{2,5};%F值
    91.     if p(i)<=0.01 && F>fa
    92.         disp([Z(i)+"对水肿体积进展模式影响非常显著"])
    93.         fprintf('p值为%.4f<0.01，F值为%.2f>%.2f\n',p(i),F,fa)
    94.     elseif p(i)<=0.05 && F>fa
    95.         disp([Z(i)+"对水肿体积进展模式影响显著"])
    96.         fprintf('p值为%.4f<0.05，F值为%.2f>%.2f\n',p(i),F,fa)
    97.     else
    98.         disp([Z(i)+"对水肿体积进展模式影响不显著"])
    99.         fprintf('p值为%.4f，F值为%.2f\n',p(i),F)
    100.     end
    101. end
    102. disp('不同治疗对水肿进展模式的影响大小为：')
    103. [~,q]=sort(p);
    104. str=Z(q)+">";
    105. disp(str)
    106. %% d
    107. %同上述步骤求血肿体积的
    108. %计算水肿指标的变化率，在不同治疗方法下，改善效率
    109. K2=[];
    110. for i=1:length(Num)
    111.     a=find(Table(:,1)==Num(i));
    112.     if length(a)>=5%至少五次检查数据才进行计算
    113.         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)));
    114.         kk=find(k<0);
    115.         if length(kk)==0
    116.             K2(i,1)=0;
    117.         else
    118.             K2(i,1)=mean(k(kk));
    119.         end
    120.     else
    121.         K2(i,1)=NaN;
    122.     end
    123. end
    124. c=setdiff([1:length(K2)],find(isnan(K2)==1));
    125. for i=1:length(Z)
    126.     [p2(i),anovatab2{i},stats2{i}]=anova1(K2(c),G(c,i)+1,'off');%单因素方差分析
    127.     fa=finv(0.95,anovatab2{i}{2,3},anovatab2{i}{3,3});%计算fa
    128.     F=anovatab2{i}{2,5};%F值
    129.     if p2(i)<=0.01 && F>fa
    130.         disp([Z(i)+"对血肿肿体积进展模式影响非常显著"])
    131.         fprintf('p值为%.4f<0.01，F值为%.2f>%.2f\n',p2(i),F,fa)
    132.     elseif p2(i)<=0.05 && F>fa
    133.         disp([Z(i)+"对血肿体积进展模式影响显著"])
    134.         fprintf('p值为%.4f<0.05，F值为%.2f>%.2f\n',p2(i),F,fa)
    135.     else
    136.         disp([Z(i)+"对血肿体积进展模式影响不显著"])
    137.         fprintf('p值为%.4f，F值为%.2f\n',p2(i),F)
    138.     end
    139. end
    140. disp('不同治疗对血肿进展模式的影响大小为：')
    141. [~,q2]=sort(p2);
    142. str2=Z(q2)+">";
    143. disp(str2)
    144. %在算下血肿指标与水肿指标的相关性
    145. x0=double(Table(2:end,23));
    146. y0=double(Table(2:end,34));
    147. theta=x0'*y0/(norm(x0)*norm(y0));%余弦相似度
    148. fprintf('血肿指标与水肿指标的相关度为：%.4f\n',theta)
    149.  
    150. disp('结果见矩阵result2')
    1. %第三问
    2. clear
    3. clc
    4. %是否要进行指标降维自行添加
    5. load data
    6. Num=unique(Table(2:end,1),'stable');
    7. A=[];%首检
    8. Y=[];%首检Rms
    9. B=[];%随检
    10. N1=[];
    11. N2=[];
    12. Z=[Table(1,[4:end])];
    13. for i=1:length(Num)
    14.     a=find(Table(:,1)==Num(i));
    15.     A=[A;double(Table(a(1),[4:end]))];
    16.     N1=[N1;Table(a(1),1:2)];
    17.     B=[B;double(Table(a(2:end),[4:end]))];
    18.     N2=[N2;Table(a(2:end),1:2)];
    19. end
    20. %% a
    21. In=A;
    22. Y=xlsread('表1-患者列表及临床信息.xlsx','B2:B101');
    23. Out=Y;
    24. %bp神经网络
    25. Xtrain = In(1:100,:);
    26. Ytrain = Out;
    27. Xtest1 = In(101:130,:);
    28. Xtest2 = In(131:160,:);
    29. %%  训练模型
    30. trees = 100;                                      % 决策树数目
    31. leaf  = 5;                                        % 最小叶子数
    32. OOBPrediction = 'on';                             % 打开误差图
    33. OOBPredictorImportance = 'on';                    % 计算特征重要性
    34. Method = 'classification';                            % 分类还是回归
    35. net = TreeBagger(trees, Xtrain, Ytrain, 'OOBPredictorImportance', OOBPredictorImportance,...
    36.       'Method', Method, 'OOBPrediction', OOBPrediction, 'minleaf', leaf);
    37. %预测
    38. Y=string(predict(net,Xtrain));
    39. Y1=string(predict(net,Xtest1));
    40. Y2=string(predict(net,Xtest2));
    41. %混淆矩阵
    42. figure
    43. plotconfusion(categorical(Ytrain),categorical(Y));
    44. title('训练集')
    45. figure
    46. hold on
    47. auc=plot_roc(double(Y)',Ytrain');
    48. legend(['训练集 auc=',num2str(round(auc,2))])
    49. %% b
    50. Xyuce = B;
    51. Y3 = string(predict(net,Xyuce));
    52. result3="";
    53. for i=1:length(Num)
    54.     a=find(N2==Num(i));
    55.     for j=1:length(a)
    56.         result3(i,j)=Y3(a(j));
    57.     end
    58. end
    59. result3(ismissing(result3)==1)="";
    60. z=[];
    61. for i=1:size(result3,2)
    62.     z=[z,"第"+num2str(i)+"次随访"];
    63. end
    64. result3=[["患者","首次检测";Num,[Y;Y1;Y2]],[z;result3]];
    65. %% c
    66. %相关性分析
    67. X=double(Table(2:end,[4:end]));
    68. P=[];
    69. for i=1:size(X,2)
    70.     for j=1:size(X,2)
    71.         P(i,j)=X(:,i)'*X(:,j)/(norm(X(:,i))*norm(X(:,j)));
    72.     end
    73. end
    74. figure%热图
    75. set(gca,'position',[0.05 0.05 0.9 0.9])
    76. h=heatmap(P,'ColorbarVisible', 'on');
    77. h.FontSize = 8;
    78. h.CellLabelFormat = '%0.2g';
    79. resultp=["",Z;Z',P];%自行分析相关性结果
    80.  
    81. disp('预测结果见result3，相关性结果见resultp')

  • 相关阅读:
    PMP考试可以延缓考吗？解答来了!
    生成式人工智能能否使数字孪生在能源和公用事业行业成为现实？
    洛谷P3092 状压dp
    c# HashSet vs SortedSet
    【故障公告】今年第五次：数据库服务器 CPU 100%
    操作系统实验（4）—— 处理机调度算法模拟实现与比较
    华为云排名仅为第二？市场占有率营收被阿里双碾压
    软件测试开发从0到1（自用，更新中）
    北京有什么好的创业孵化器?中移产孵中心助创业
    SSH Remote Linux Servers By Public Key Authentication
  • 原文地址：https://blog.csdn.net/qq_39899679/article/details/133189476