2023年华中杯数学建模
A题 新型镇静药物临床实验疗效分析与预测
原题再现
临床研究是新药物研究中的关键环节。本题拟围绕一种新型镇静药物的临床实验数据分析展开。 尝试根据附件中提供的数据和相关材料，研究以下问题：
1. 关于术中、术后 24h 不良反应，新药组和原有药物组是否存在显著差异；能否建立一个有效的数学模型，根据患者基本信息和镇静药物种类，对患者术中、术后 24h 的不良反应进行预判。
2. 新药组和原有药物组在生命体征数据方面是否表现出显著差异；如果有显著差异，能否确定是由于新药造成，还是由其他因素造成。
3. 临床经验表明，用药后 3 分钟内的 IPI 数据对于病情后续恢复具有决定作用，尝试根据用药信息和患者信息对给药后 3 分钟以内的 IPI 数据进行预测。
4. 术后满意度与很多因素有关，包括护理、身体恢复程度等等，甚至有一些因素无法观测到。基于现有数据是否能够找出术后满意度与哪些因素有关？有怎样的关系。
附件说明：
1. 附件 1 提供了新药物“R 药”和原有药物“B 药”对比的临床实验研究数据，其中包括患者基本信息、镇静药及其诱导剂量，以及病患相关生命体征数据，包括 sbp, dbp, petco, RR, spO2, HR, IPI, moaas 等。
2. 附件 2 提供了病患信息采集表，以及相关指标说明。
3. 附件 3 提供了一份简短的相关背景信息阅读材料。
整体求解过程概述(摘要)
本题主要研究新型镇静药物在临床试验过程中的疗效并对药物使用进行可能后果预测，意在掌握新型镇静药物的疗效，对后续药物研发改进以及推广使用等具有重大的参考意义。
针对问题一，首先，我们使用 KS 检验，得到样本不服从正态分布。其次，基于样本为离散型变量，我们选择卡方检验对术中术后分别进行差异性显著分析。得到术中不良反应在两种药物使用中均有显著差异性；术后 24h 不良反应在两种药物使用中恶心、头昏、乏力、腹胀和腹痛在 99%的置信水平下存在显著性差异，呕吐、头晕、头痛、嗜睡和术后其他则不存在显著性差异。最后，我们采用 Logistic 回归模型，分别得到了术中、术后 24h 不良反应的相关系数，从而对这些不良反应进行较精准预判。
针对问题二，探究 R 药和 B 药在生命体征方面是否存在显著差异性，首先我们通过统计描述，得到样本分布类似于正态，但不是正态分布。其次，基于样本为连续型变量，我们选择 Mann-Whitney U 检验对整体进行差异性分析，得到除 moaas3 以外的生命体征数据在两种药物使用中在 99%的置信水平下有显著差异性。最后，由于探究相关因素的因果推断存在样本数据不均衡的情况，所以我们选择相关分析，利用Pearson 相关系数，得到生命体征数据与镇静药物相关性较弱，不能够确定差异是由新药造成。
针对问题三，根据用药信息和患者信息对给药后 3min 内的 IPI 进行预测。首先，我们考虑到全过程的 IPI 数据对三分钟以内的 IPI 数据的预测都会有影响，于是我们通过对十分钟以内的 IPI 数据进行加权平均，从而得到一个具有代表性的指标，对此进行偏最小二乘回归来对三分钟以内的 IPI 数据进行预测，得到酗酒对 IPI 的正向影响最大，IPI 越大，越处于高危状态。其次，吸烟、PONV、镇静药总剂量对 IPI 也有较大的正向影响。最后我们用 BP 神经网络模型和计算模型的均方误差两个方式对我们模型进行了评估，我们模型的拟合度更高，模型更为优良。
针对问题四，由于自变量个数过多，我们对自变量了进行降维处理。首先，我们利用最小二乘法（OLS）检验，得到自变量之间具有多重共线性。进而，我们选择用Lasso 回归将不重要变量的回归系数“惩罚”至零，得到术后满意度与年龄、镇静药总剂量和镇痛药总剂量负相关，与体重成正相关。
模型假设：
1. 假设各样本数据真实且具有代表性。
2. 连续的自变量与因变量的 logit 转换值之间存在线性关系。
3. 自变量患者的基本数据(身高、体重等）观测无误差；其他影响因子（如注射部位，对这两种特定药物的耐药性等）处于同一水平。
问题分析：
问题一的分析
第一小问首先需要我们判断对于术中、术后 24h 不良反应，新药组和原有药物组是否存在显著差异。我们知道进行差异性显著分析有多种办法，为了找到更适合我们样本的分析方法，我们初步使用 KS 检验来检验样本是否服从正态分布。在此基础上，又因为样本为离散型变量，我们最后选择卡方检验进行差异性显著分析。第二小问需要我们建立数学模型对患者术中、术后 24h 的不良反应进行预判。由于因变量的特点是二值变量，显然是一个二元决策问题。同时，分析各个检测指标对检测结果的影响，属于多元分析。为了较精准地对术中、术后的不良反应进行预判，我们对各种不良反应分别采用 Logistic 回归模型，从而提高预判的正确率。
问题二的分析
第一小问首先是需要判断新药组和原有药物组在生命体征数据方面是否表现出显著差异。在这里，我们将数据与第一问的数据进行了比较，发现这里的样本为连续型变量。进一步对数据进行分析，描述性统计的方法对样本是否服从正态分布进行了初步的判断。最后我们选择 Mann-Whitney U 检验对整体数据进行差异性显著分析。第二小问建立在第一问基础之上，确定差异是否由新药造成。判断新药是否是差异形成的原因，显然最理想的方法是因果推断。我们通过对样本进行分析，发现样本数量少且极度不均衡，通过查找资料，我们认为进行因果推断会造成较大的误差，从而难以进行准确地判定。所以我们采用相关分析，利用 Pearson 相关系数来说明各因素对差异造成的相关程度，从而以此作为结论。
问题三的分析
问题三需要我们根据用药信息和患者信息对给药后 3 分钟以内的 IPI 数据进行预测。首先，我们对手术全过程的 IPI 数据进行整理，发现 IPI15、IPI20 有大量的数据缺失，并且考虑到全过程的 IPI 数据对三分钟以内的 IPI 数据的预测都会有影响，如果我们只单独考虑给药后三分钟以内的 IPI 数据来进行预测，会因为丢失部分信息而造成预测结果的误差。于是我们在进行数据的清洗之后，通过对十分钟以内的 IPI数据进行加权平均，从而得到一个具有代表性的指标，在此基础之上，再通过偏最小二乘回归来对三分钟以内的 IPI 数据进行预测。最后我们还通过用 BP 神经网络进行机器学习的预测模型，和计算偏最小二乘回归模型的均方误差两个方式对我们模型进行了评估。
问题四的分析
问题四需要我们基于现有数据，找出与术后满意度相关的因素。通过分析，我们初步认为术后满意度与患者的基本信息、三分钟内的 IPI 数据、麻醉医生的满意程度、内镜医生的满意程度、术后的不良反应等都有相关性。由于自变量个数过多，我们需要对自变量进行降维处理。进一步，我们我们利用最小二乘法检验自变量之间是否具有多重共线性。最后，通过分析，我们选择用 Lasso 回归模型将不重要变量的回归系数“惩罚”至零。从而较精准地找出与术后满意度相关的因素，并且通过回归系数来说明各因素对术后满意度的相关程度。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
%导入处理后的数据 cleanedData
mu=mean(cleanedData);std=std(cleanedData);
rr=corrcoef(cleanedData); %求相关系数矩阵
A=(cleanedData-mu)./std; %数据标准化
X=A(:,[1:13]); %自变量
Y=A(:,14); %因变量
% 提出主成分对
[XL,YL,XS,YS,BETA,PCTVAR,MSE,stats] =plsregress(X,Y);
Xw=X\XS; %求自变量提出成分的系数，每列对应一个成分，这里 Xw 等于
stats.W
Yw=Y\YS; %求因变量提出成分的系数
a_0=PCTVAR(1,:);b_0=PCTVAR(2,:);% 自变量和因变量提出成分的贡献
率
a_1=cumsum(a_0);b_1=cumsum(b_0);% 计算累计贡献率
i=1;%初始值 i
while
((a_1(i)<0.9)&(a_0(i)>0.05)&(b_1(i)<0.9)&(b_0(i)>0.05)) % 提取
主成分的条件
i=i+1;
end
nc=i;% 选取的主成分对的个数
fprintf('%d 对成分分别为：\n',nc);% 打印主成分的信息
for i=1:nc
fprintf('第%d 对成分：\n',i);
fprintf('u%d=',i);
for k=1:size(X,2)%此处为变量 x 的个数
fprintf('+(%f*x_%d)',Xw(k,i),k);
end
fprintf('\n');
fprintf('v%d=',i);
fprintf('\n');
end
% 确定主成分后的回归分析
[XL2,YL2,XS2,YS2,BETA2,PCTVAR2,MSE2,stats2]
=plsregress(X,Y,nc);
n=size(X,2); m=size(Y,2);%n 是自变量的个数,m 是因变量的个数
BETA3(1,:)=mu(n+1:end)-mu(1:n)./std(1:n)*BETA2([2:end],:).
*std(n+1:end); %原始数据回归方程的常数项
BETA3([2:n+1],:)=(1./std(1:n))'*std(n+1:end).*BETA2([2:end
],:); %计算原始变量 x1,...,xn 的系数，每一列是一个回归方程
%贡献率画图
figure
percent_e1 = 100 * PCTVAR(1,:) / sum(PCTVAR(1,:));
pareto(percent_e1);
xlabel('主成分')
ylabel('贡献率(%)')
title('主成分对自变量的贡献率')
figure
percent_e2 = 100 * PCTVAR(2,:) / sum(PCTVAR(2,:));
pareto(percent_e2);
xlabel('主成分')
ylabel('贡献率(%)')
title('主成分对因变量的贡献率')
% 求预测值 yp 和真实值
yp=repmat(BETA3(1,:),[size(X,1),1])+A(:,[1:n])*BETA3([2:en
d],:)-3.6;
y0 = A(:,end-size(yp,2)+1:end);
bar(BETA2','k') % 画回归系数的直方图
% 三种方法检验网络性能
for i =1:size(Y,2)
yz = y0(:,i);% 真实值
yc = yp(:,i);% 预测值
perf = mse(y0,yp) % 均方误差
figure;
plotregression(yz,yc,['第',num2str(i),'个回归图'])% 第二
种方法，回归图
eps = yz-yc;
figure;
ploterrhist(eps,['第',num2str(i),'个误差直方图']) % 第三
种方法，误差直方图
end
