2023年数维杯数学建模C题宫内节育器的生产求解全过程文档及程序

2023年数维杯数学建模

C题 宫内节育器的生产

原题再现：

  宫内节育器（IUD）是一种相对安全、有效、经济、可逆、简便，广大妇女易接受的节育器具，目前已成为我国育龄妇女的主要避孕措施。据悉，我国约70%妇女选用 IUD 作为避孕方法,占世界 IUD 避孕总人数的 80%。某公司研发了两种型号的 VCu 记忆型宫内节育器，分别为 VCu260 记忆型宫内节育器和VCu380 记忆型宫内节育器。该公司研发的 VCu 记忆型宫内节育器采用镍钛记忆合金丝支架，除具有独特形状记忆功能外，还具有抗腐蚀、耐磨损、超弹性和对身体的副作用较小等优点。但是节育器在给女性带来方便的同时，也可能会引起疼痛、不适、脱落或者出血等症状。为了探究这两种型号的节育器是否适合投入生产，特与已被临床应用的 MCu 功能性宫内节育器一起临床试验。一共有三组，MCu 功能性宫内节育器作为对照组（标记为 1 组），VCu260 与 VCu380 记忆型宫内节育器均为试验组（分别标记为 2、3 组），每组入组均为 525 例，且分别于第一、三、六、十二个月就主诉情况进行随访记录。本次试验选择两家医院并同时进行临床试验。节育器的理化指标、受试者的身体指标、随访时的主诉情况均在一定程度上反映节育器的质量。
  附件 1 提供了两个医院临床受试者的身体指标与节育器的理化指标，附件 2提供了两个医院随访时的主诉情况记录。请建立数学模型研究下列问题：
  1. 根据附件 1 与附件 2，分析两个医院的临床数据有无显著性差异，若存在显著性差异，对导致这种差异的因素进行分析。
  2. 结合附件 1 与附件 2，分析受试者的身体指标与随访主诉情况的联系，并说明受试者的身体指标是否是受试者出现不适状况的主要因素。
  3. 根据受试者的身体指标、节育器的理化指标与随访时的主诉情况，建立节育器质量模型，并分析 VCu260 与 VCu380 记忆型宫内节育器的质量哪个更优，更适合生产。
  4. 结合问题 3，根据建立的节育器质量模型，探究影响宫内节育器质量的决定因素。
  附录说明
  随访时只对出现不适情况的受试者进行记录，使用感受良好的受访者不予记录；部分受试者在随访时难以联系，造成失访，会有部分数据缺失。
  使用节育器的目的是防止怀孕，生产节育器的前提是节育器在很大程度上能够避免怀孕。
  若受试者在使用宫内节育器前期出现不适症状，但后期症状消失，可认为受试者适应该节育器。例如受试者在第一、三个月份随访时均主诉有不适症状，但是在第六、十二个月份随访时均表示没有不适症状，可认为前三个月为受试者的适应期，此时不考虑节育器的质量问题。
  附件
  附件 1：两个院临床受试者和节育器的基本数据
  附件 2：两个医院随访的节育器使用后主诉情况

整体求解过程概述(摘要)

  在生物医疗技术迅猛发展之下，宫内节育器（IUD）因其具有相对安全、有效和简便等优点逐渐受到广大妇女青睐，目前已成为我国育龄妇女的主要避孕措施。本文利用多种统计学方法对受试者的身体指标、节育器的理化指标和随访主诉情况进行建模分析，通过综合评价的手段在 Vcu260 型和 VCu380 型节育器中选择一个更优的产品进行生产，并且给出了影响节育器质量的相关因素。
  对于数据处理方面，对数据的缺失值进行了替换填充，并且对于数据特征中重复、异常的数据点进行了处理，利用题目所给的条件对附件 2 的数据进行了预处理。
  针对问题 1，观察数据发现身体指标为定量数据、理化指标和主诉情况为定性数据；首先通过数据形态验证了身体指标服从正态分布，并对身体指标进行了独立样本 t 检验，对理化指标和主诉情况进行了卡方检验，发现只有身体指标中的初潮年龄、月经周期和主诉情况中的脱落、因症取出…等 8 个指标没有显著性差异，其余的 9 个指标存在显著性差异，最后还分析了出现差异性可能的原因。
  针对问题 2，首先对数据进行了皮尔逊相关性分析，初步分析了身体指标与主诉情况的联系，然后利用典型相关分析模型再次衡量了身体指标与主诉情况的联系:发现宫颈深度使主诉情况倾向好的方面，月经经期使主诉情况倾向不好的方面，并且通过分析主诉情况发现节育环对怀孕的抑制是十分有效的，负面的主要影响是分泌物增多和经期/周期异常。最后通过比较身体指标与主诉情况、理化指标与主诉情况的典型相关系数大小判断出身体指标不是其主要原因，理化指标影响力比身体指标多 12.64%。
  针对问题 3，首先建立了二级指标评价体系，并将所有指标都转变为了正向型指标，接着利用熵权法确定各个指标的权重，然后对样本按照 VCu260 型和 VCu380 型节育器的标准进行了分类，利用 Topsis 评价模型对节育器质量进行多级评价，最后按加权平均的方法计算了每个样本对节育器质量的总评分，以两个类别的样本总评分的均值来判断哪个类型的节育器更优。最后计算得到 VCu260 型和 VCu380 型节育器均分为0.6921 和 0.7001，VCu380 型节育器优于 VCu260 型节育器，更适合生产。
  针对问题 4，首先利用问题 3 的权重模型初步分析出理化指标是影响节育器质量的决定因素；接着利用身体指标和理化性质对于是否出现不适症状进行了 Logistic 回归分析，同样可以发现节育器的理化指标是影响节育器质量的决定因素，其中使用大型节育器是负向因素，放置节育器时宫颈扩张情况是正向因素。
  在模型的验证中，利用斯皮尔曼相关系数对皮尔逊相关系数进行了误差分析，发现这两种方法计算出来的相关系数相差不大，最大误差不超过 3%，进而验证了模型的稳定性。

问题分析：

  问题 1 的分析
  对于问题 1，分析两个医院的临床数据有无显著性差异，并针对有显著差异的指标原因探究。观察附件一的数据可以发现，临床数据指标主要分为两个方面，一是受试者的身体指标包括年龄、初潮年龄、月经周期、月经经期和宫颈深度；二是节育器的理化指标包括既往节育器使用情况（IUD、无和其他）、节育器型号（小、中和大）和放置节育器时扩张情况，其中身体指标为定量变量，理化指标为定性变量。通常差异性分析可以通过 t 检验、卡方检验和非参数检验进行，这三种方法的适用范围不同，t检验要求变量为定量变量且服从正态分布，卡方检验主要用于定性变量，非参数检验主要用于不满足正态分布的数据。首先对数据进行预处理后，将数据分为定量和定性两类，定量的数据再接着进行正态分布检验，若服从正态分布，则采用 t 检验，反之则采用非参数检验；对于定性的数据采用卡方检验。经检验后，对于存在显著性差异的指标进行探究性因素分析，给出导致存在差异性可能的原因。
  问题 2 的分析
  对于问题 2，分析受试者的身体指标与随访主诉情况的联系，并说明受试者的身体指标是否是受试者出现不适症状的主要原因。首先，由于某些受试者的数据缺失，为此要对这些缺失值进行处理，并且如果第一、三个月有不适症状，但第六、十二没有不适症状，则认为该症状对该样本无影响。考虑到受试者的身体指标与随访主诉情况这两类指标里均包含多个子指标，而典型相关分析正好对应这种“多对多”的指标相关性分析，因此采用典型相关分析来分析这两类指标之间的联系。同样地，也可以将受试者的身体指标与随访主诉情况进行典型相关性分析，比较受试者的身体指标与随访主诉情况和节育器的理化指标与随访主诉情况典型相关性的值大小即可以判断受试者的身体指标是否是受试者出现不适症状的主要原因。
  问题 3 的分析
  对于问题 3，需要根据受试者的身体指标、节育器的理化指标与随访主诉情况，建立质量评价模型，来判断 VCu260 型节育器和 VCu380 节育器哪个更适合生产。首先，建立多级评价指标体系，并需要考虑指标的性质，将全部指标都转变为正向型指标；由于每一个指标的影响程度都不相同，因此需要确定各个指标的权重，本文采用熵权法确实权重；接着建立 Topsis 节育器质量评价体系，分别计算 3 个一级指标的评分，再对这 3 个一级指标评分按加权求和计算每一样本的总评分。最后，按 VCu260 型节育器和 VCu380 节育器的样本进行分类，计算这两类节育器的评分均分，均分高的则更适合生产。
  问题 4 的分析
  对于问题 4，需要结合问题 3 的节育器质量模型探究影响节育器质量的决定因素。由于问题 3 已经通过熵权法计算了每一个指标的权重，故可以依据这个指标权重来判断哪个指标对节育器质量的影响更大。为进一步验证指标的重要性，建立受试者的身体指标和节育器的理化指标与是否出现不适症状的 Logistic 回归模型，通过其指标的影响系数以及显著性进一步说明哪个指标对节育器质量影响更大。

模型假设：

  1) 假设数据来源真实有效；
  2) 假设除题目所给指标外，其他方面的指标对节育器的质量无影响；
  3) 假设两个医院的数据为独立数据，不相互影响。

论文缩略图：

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：

clc;clear;close all
set(0,'defaultfigurecolor','w');
%% loaddata
Type11 = '附件 1';
Path11 = '.\Data\附件 1：两个院临床受试者和节育器的基本数据.xlsx';
Sheet11 = "一院临床受试者和节育器的基本数据";
Data11 = loadData(Path11, Sheet11, Type11);
Type12 = '附件 1';
Path12 = '.\Data\附件 1：两个院临床受试者和节育器的基本数据.xlsx';
Sheet12 = "二院临床受试者和节育器的基本数据";
Data12 = loadData(Path12,Sheet12,Type12);
Type21 = '附件 2';
Path21 = '.\Data\附件 2：两个医院随访的节育器使用后主诉情况.xlsx';
Sheet21 = "一院随访的节育器使用后主诉情况";
Data21 = loadData(Path21,Sheet21,Type21);
Type22 = '附件 2';
Path22 = '.\Data\附件 2：两个医院随访的节育器使用后主诉情况.xlsx';
Sheet22 = "二院随访的节育器使用后主诉情况";
Data22 = loadData(Path22,Sheet22,Type22);
%% 预处理
Data11f = fillmissing(Data11,"constant",0);
Data12f = fillmissing(Data12,"constant",0);
Data21f = fillmissing(Data21,"constant",0);
Data22f = fillmissing(Data22,"constant",0);
Data22f(233:426,:) = [];
loss1 = sum(Data21f.loss1) + sum(Data22f.loss1);
loss3 = sum(Data21f.loss3) + sum(Data22f.loss3);
loss6 = sum(Data21f.loss6) + sum(Data22f.loss6);
loss12 = sum(Data21f.loss1) + sum(Data22f.loss12);
n_data = size(Data21f,1) + size(Data22f,1);
alpha1 = loss1 / n_data;
alpha3 = loss3 / n_data;
alpha6 = loss6 / n_data;
alpha12 = loss12 / n_data;
fprintf('alpha1=%.4f%%,\nalpha3=%.4f%%,\nalpha6=%.4f%%,\nalpha12=%.4f%%n',alpha1*100,alpha3*100,alpha6*100,alpha12*100)
%% 合并分组
Data11f.category = ones(size(Data11f,1),1);
Data12f.category = ones(size(Data12f,1),1)*2;
Data1All = [Data11f;Data12f];
Data21f.category = ones(size(Data21f,1),1);
Data22f.category = ones(size(Data22f,1),1)*2;
Data2All = [Data21f;Data22f(2:end,:)];
% DataAll = [Data1All,Data2All];
Data2AllTemp = Data2All;

Data2AllTemp = Data2AllTemp(2:end,7:42);
%% 再处理
outputNames2 = {'num','group','fallOff','takeOut','pregnant','blood','pain',... 'mQuantity','hypersecretion' ,'anomaly', 'sick','category'};
outputData2All = table('Size',[size(Data2All,1)-1,length(outputNames2)],'VariableTypes',{'double','double','double', ...
 
'double','double','double','double','double','double','double','double','double'},'VariableNames', outputNames2);
outputData2All.num = Data2All.num(2:end);
outputData2All.group = Data2All.group(2:end);
for i = 1:size(Data2AllTemp,2)/4
 DataTemp = Data2AllTemp(:,i*4-3:i*4);
 sum12 = sum(table2array(DataTemp(:,1:2)),2);
 sum34 = sum(table2array(DataTemp(:,3:4)),2);
 idx12 = find(sum12 > 0);
 idx34 = find(sum34 > 0);
 outputData2All.(outputNames2{i+2})(idx34) = 1;
end
outputData2All.category = Data2All.category(2:end);
outputNames1 = 
{'num','group','age','menarcheAge','menstrualCycle','menstrualPeriod','iud',... 'uterineDepth','iudModel' ,'expansion', 'category'};
outputData1All = table('Size',[size(Data1All,1),11],'VariableTypes',{'double','double','double', ...
'double','double','double','double','double','double','double','double'},'VariableNames', outputNames1);
outputData1All.num = Data1All.num;
outputData1All.group = Data1All.group;
outputData1All.age = Data1All.age;
outputData1All.menarcheAge = Data1All.menarcheAge;
outputData1All.menstrualCycle = Data1All.menstrualCycle;
outputData1All.menstrualPeriod = Data1All.menstrualPeriod;
outputData1All.uterineDepth = Data1All.uterineDepth;
outputData1All.expansion = Data1All.expansion;
outputData1All.category = Data1All.category;
% Data1AllTemp = Data1All;
idx1 = find(Data1All.iud == 1);
idx2 = find(Data1All.iudNone == 1);
idx3 = find(Data1All.iudOther == 1);
outputData1All.iud(idx1) = 1;
outputData1All.iud(idx2) = 2;
outputData1All.iud(idx3) = 3;
idx1 = find(Data1All.iudModelS == 1);
idx2 = find(Data1All.iudModelM == 1);
idx3 = find(Data1All.iudModelB == 1);
outputData1All.iudModel(idx1) = 1;
outputData1All.iudModel(idx2) = 2;
outputData1All.iudModel(idx3) = 3;
% 补充主诉情况的数据
outputData2AllNew = 
table('Size',[size(Data1All,1),length(outputNames2)],'VariableTypes',{'double','double','double', ...'double','double','double','double','double','double','double','double','double'},'VariableNames', outputNames2);
outputData2AllNew.num = outputData1All.num;
outputData2AllNew.group = outputData1All.group;
outputData2AllNew.category = outputData1All.category;
% [UU,idxUU1] = intersect(outputData2AllNew.num,outputData2All.num);
% idxNew = find(outputData2AllNew.num == outputData2All.num);
k = 1;
% aa = [];
for i = 1:size(outputData2AllNew,1)
 for j = 1:size(outputData2All,1)
 if outputData2AllNew.num(i) == outputData2All.num(j) && 
outputData2AllNew.group(i) == outputData2All.group(j)...
 && outputData2AllNew.category(i) == outputData2All.category(j)
 outputData2AllNew(i,:) = outputData2All(j,:);
 k = k + 1;
% aa = [aa,j];
 break
 end
 end
end
outputDataAll = [outputData1All,outputData2AllNew(:,3:end-1)];
outputDataAll1 = [Data1All,outputData2AllNew(:,3:end-1)];
outputDataAll1.sick = ~outputDataAll1.sick;
% aaa = sum(outputDataAll1{:,16:23},2);
% length(find(aaa>0))
%% 正态分布检验
Name = 
{'age','menarcheAge','menstrualCycle','menstrualPeriod','uterineDepth'};
Label = {'年龄/岁','初潮年龄/岁','月经周期/天','月经经期/天','宫腔深度/cm'};
Beta = [11.5,3.3,3.8,2.3,2.0];
nName = length(Name);
for i = 1:nName
 U = unique(Data1All.(Name{i}));
 C = zeros(length(U),1);
 Xticks = cell(1,length(U));
 for j = 1:length(U)
 C(j) = length(find(Data1All.(Name{i}) == U(j)));
 Xticks{j} = num2str(U(j));
 end
 figure(i)
 b = bar(U,C);
 b.FaceColor = [0.8500 0.3250 0.0980];
 b.EdgeColor = [0.3,0.3,0.3];
 hold on
 [miu,sigma] = normfit(Data1All.(Name{i}));
 xx = min(U):(max(U)-min(U))*0.005:max(U);
 yy = normpdf(xx,miu,sigma)*max(C)*Beta(i);
 plot(xx,yy,'-','Color',[0.2 0.2 0.2],'LineWidth',1.3)
 
 xlabel(Label{i})
 ylabel('频数')
 set(gca,'FontName','TimesSimSun','FontSize',11) %设置坐标轴刻度字体名
称，大小
 exportgraphics(gca,['.\Fig\bar',Name{i},'.png'], 'Resolution',600)
end
%% 绘制相关系数矩阵
idxP = [3:10,12,13,15:20];
[R,P] = corrcoef(table2array(outputDataAll(:,idxP)));
Label = outputDataAll.Properties.VariableNames(idxP);
scrsz = get(0,'ScreenSize'); %%%% 获取屏幕的尺寸
figure('Name','hotR','Position',[0 30 scrsz(3) scrsz(4)-95]); 
hot_figure = heatmap(Label,Label,R);
hot_figure.GridVisible = 'off';
colormap(gca,'jet')
set(gca,'FontName','Times New Roman','FontSize',12.5) %设置坐标轴刻度字体名
称，大小
exportgraphics(gca,'./Fig/Pearson.png', 'Resolution',600)
R1 = corr(table2array(outputDataAll(:,idxP)),'type' , 'Spearman');
figure('Name','hotR_P','Position',[0 30 scrsz(3) scrsz(4)-95]); 
hot_figure = heatmap(Label,Label,R1);
hot_figure.GridVisible = 'off';
colormap(gca,'jet')
set(gca,'FontName','Times New Roman','FontSize',12.5) %设置坐标轴刻度字体名
称，大小
exportgraphics(gca,'./Fig/Pilsman.png', 'Resolution',600)
Rl = tril(R);
Ru = triu(R1);
RR = Rl + Ru;
for i = 1:length(RR)
 RR(i,i) = 1;
end
f = figure('Name','hotRR','Position',[0 30 scrsz(3) scrsz(4)-95]); 
hot_figure = heatmap(Label,Label,RR);
hot_figure.GridVisible = 'off';
colormap(gca,'jet')
set(gca,'FontName','Times New Roman','FontSize',12.5) %设置坐标轴刻度字体名称，大小
exportgraphics(gca,'./Fig/PP.png', 'Resolution',600)
idxR = find(R ~= 1);
Rexpension = R(idxR);
idxR1 = find(R1 ~= 1);
Rexpension1 = R1(idxR1);
E = Rexpension - Rexpension1;
figure('Name','Rline')
p1 = plot(Rexpension,'-
ko','MarkerIndices',1:round(length(Rexpension)*0.05):length(Rexpension),
'LineWidth',1.0);
p1.Color = [0.3,0.3,0.3];
hold on 
p2 = plot(Rexpension1,'--
r*','MarkerIndices',2:round(length(Rexpension)*0.05):length(Rexpension1),'LineWidth',1.0);
set(gca,'FontName','Times New Roman','FontSize',12.5) %设置坐标轴刻度字体名称，大小
legend('Pearson','Spearman','Location','northwest','fontSize',15)
xlabel('Index')
ylabel('相关系数','FontName','宋体')
xlim([0,length(Rexpension)])
exportgraphics(gca,'./Fig/相关系数对比.png', 'Resolution',600)
figure('Name','error') 
% plot(E)
s = stem(E);
s.Color = [0.5,0.5,0.5];
s.LineWidth = 1.0;
s.MarkerFaceColor = [0.8,0.8,0.8];
set(gca,'FontName','Times New Roman','FontSize',12.5) %设置坐标轴刻度字体名
称，大小
xlim([0,length(E)])
xlabel('Index')
ylabel('Error')
exportgraphics(gca,'./Fig/相关系数误差.png', 'Resolution',600)
%% 矩阵散点图
f = figure('Name','figMatrix','Position',[0 30 scrsz(3) scrsz(4)-95]);
[Sf,AX,BigAx,H,HAx] = plotmatrix(table2array(outputDataAll(:,idxP)));
for i = 1:length(idxP)
 H(i).FaceColor = [0.3,0.3,0.3];
 H(i).EdgeColor = [0.1,0.1,0.1];
end
for i = 1:length(idxP)
 for j = 1:length(idxP)
 Sf(i,j).MarkerSize = 2;
 Sf(i,j).Color = 'k';
 Sf(i,j).Marker = 'o';
 Sf(i,j).MarkerFaceColor = 'k';
 set(AX(i,j),'FontName','Times New Roman','FontSize',9) %设置坐标轴刻度字体名称，大小
 end
end
exportgraphics(f,'./Fig/散点矩阵.jpg', 'Resolution',600)
%% % 保存数据
writetable(Data11f,'.\Data\附件 1：两个院临床受试者和节育器的基本数据改.xlsx','Sheet','一院临床受试者和节育器的基本数据
','WriteMode','overwritesheet');
writetable(Data12f,'.\Data\附件 1：两个院临床受试者和节育器的基本数据改.xlsx','Sheet','二院临床受试者和节育器的基本数据
','WriteMode','overwritesheet');
writetable(Data21f,'.\Data\附件 2：两个医院随访的节育器使用后主诉情况改.xlsx','Sheet','一院随访的节育器使用后主诉情况
','WriteMode','overwritesheet');
writetable(Data22f,'.\Data\附件 2：两个医院随访的节育器使用后主诉情况改.xlsx','Sheet','二院随访的节育器使用后主诉情况
','WriteMode','overwritesheet');
writetable(Data1All,'.\Data\附件 1 汇总.xlsx','Sheet','Sheet1','WriteMode','overwritesheet');
writetable(outputData1All,'.\Data\附件1汇总.xlsx','Sheet','Sheet1','WriteMode','overwritesheet');
writetable(outputData2All,'.\Data\附件2汇总.xlsx','Sheet','Sheet1','WriteMode','overwritesheet');
writetable(outputDataAll,'.\Data\附件12汇总.xlsx','Sheet','Sheet1','WriteMode','overwritesheet');
writetable(outputData,'.\Data\附件汇总无合并.xlsx','Sheet','Sheet1','WriteMode','overwritesheet');
fprintf('数据已保存至.\\Data\\中!\n')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可