1、基于DEA分析法的艾滋病评价和预测的模型 摘要 本文通过对所给数据进行统计分析,利用多项式拟合及DEA分析法研究并完成了对艾滋病疗法的评价及疗效的预测,得到了一些合理的结果。对于问题一,我们首先按照一个能够刻画病人病情的指标HIV/CD4值将所给病人分成三类, 分别进行数据挖掘,然后采用多项式拟合进行预测,得出结果如下:轻度患者在43周后停止治疗,中度患者42周后停止治疗,重度患者大约40周后停止治疗。为了解决问题二、三,本文引入了运筹学的一个新领域数据包络分析(Data Envelopment Analysis, 简称DEA)模型,一种用于评价相同类型部门(或单位)中各成员间的相对有效性的模
2、型,来评价四种疗法的优劣并作出预测。在不考虑医疗费用的情况下,疗法四是较优的,可以用作较长时期的治疗。若考虑费用影响,则它们的优劣次序为疗法二、疗法一、疗法四、疗法三。关键词 数据包络 正态分布 多项式曲线拟合 模型 模型一、问题重述 根据世界卫生组织资料报道显示,自1981年发现和报告艾滋病以来,全球感染者累计6900万人,其中死亡人数占40,每年新增500多万HIV感染病例,目前HIV/AIDS患者总数达4200多万人,而目前还没有治愈的药物和方法,只能靠预防。因而战胜爱滋病是全人类的共同任务。艾滋病是由艾滋病毒(英文简称HIV),这种病毒破坏人体的免疫系统,使人丧失抵抗力,CD4细胞在抵
3、御HIV的入侵中起着重要的作用。当CD4被HIV病毒攻击的时候,其数量急剧减少,而HIV迅速增加,导致AIDS发作。艾滋病治疗的目的,是尽量减少人体内HIV的数量,同时产生更多的CD4,至少要有效地降低CD4减少的速度,以提高人体免疫能力。本文根据美国艾滋病医疗试验机构ACTG公布的两组数据(见表1)。表1:实验疗法药品和实验数据数据药品名药品名药品名疗法名称备注ACTG320(附件1)zidovudinelamivudineIndinavir193A(附件2)600mgzidovudine400mgdianosine疗法一600ngzidovudine2.25mgzalcitabine疗法二
4、两种药按月轮换服用600mgzidovudine400mgdianosine疗法三600mgzidovudine400mgdianosine400mgnevirapine疗法四注:193A 实验人数是1300人分成4组。 ACTG320 实验人数是300人。进行实验,对艾滋病疗法进行评价及疗效的预测。问题一:利用附件1的数据,预测继续治疗的效果,或者确定最佳治疗终止时间(继续治疗指在测试终止后继续服药,如果认为继续服药效果不好,则可选择提前终止治疗)。问题二:利用附件2的数据,评价4种疗法的优劣(仅以CD4为标准),并对较优的疗法预测继续治疗的效果,或者确定最佳治疗终止时间。问题三:艾滋病药品
5、的主要供给商对不发达国家提供的药品价格如表2:表2:价格表药品名药量(mg)费用(美元)zidovudine6001.60dianosine4000.85zalcitabine2.251.85nevirapine4001.20若考虑病人的医疗费用,问题二中的评价和预测会有什么改变。二、问题分析题目中给出了两组艾滋病患者的有关数据。经过对附件1中的数据进行初步的统计分析,根据医学临床实验结果,病情不同的病人治疗效果可能不同,所以我们根据病人初始病情严重程度把病人分成三类,在每一类病人中用HIV/CD4的数值来刻画病人病情。采取多项式拟合的方法,得出HIV/CD4关于时间变化的曲线,根据曲线的变化
6、趋势,对该疗法的疗效进行预测.通过对附件2中数据的初步处理,得出该数据中病人的初始病情以及年龄接近于正态分布,从而得出这两个因素对诸疗法的疗效评价结果影响不大,所以不再按此进行分类。由于测试时刻以小数形式表示,样本点太多,我们对其进行了处理,并进行了合理分段。我们以一段时间内CD4浓度对数的最大值及其对应的时刻、该时段CD4初始浓度对数作为输入指标,该时段末CD4浓度的对数为输出指标,以四种疗法作为决策单元应用数据包络分析,对四种疗法进行评价。然后用问题一的方法进行预测。问题三再加入价格作为输入指标,以二的方法对四种疗法进行类似的评价和预测。.三、问题假设影响病情的因素中只考虑HIV和CD4,
7、忽略年龄等因素的影响。用HIV/CD4作为病情的衡量标准。假设所给数据准确,能够真实反映病人的病情。问题二中仅以CD4作为评价疗效的标准。四、模型建立与求解问题一:题目中要求预测继续治疗的结果或者最佳治疗时间,因而要一个病情标准对预测治疗的结果或者最佳治疗时间进行确定,根据临床医学和附件1中的数据,影响艾滋病患者的主要因素有HIV和CD4的数目,因而考虑用HIV和CD4的一个表达式来确定健康标准。由于HIV和CD4是相互影响的,因而用HIV和CD4的比值来表示健康状态是比较合理的,而病人的最佳状态是HIV取最小值,CD4取最大值。为了减少误差取HIV/CD4作为健康标准,而不取CD4/HIV。
8、由于附件中病人初始病情的严重程度不一样,在相同的时间里治疗的结果也不相同,因而有必要对附件1中提供的病人进行分类。 将附件1中的病人HIV/CD4计算出来,并对初始病情画直方图,如下图:根据数据分布的聚集程度将病人分为三类:轻度患者中度患者重度患者0 HIV/CD4=0.280.28 HIV/CD40.58由于同类病人的HIV/CD4值服从正态分布,我们用其均值表示这一类某一时刻的HIV/CD4值,那末我们就求得了这一类病人在不同测试时刻的HIV/CD4值(附录:表1)。我们采用多项式模型。设逼近函数的形式为:其中为健康指标,即HIV/CD4,为 测试时刻。对所得的样本,满足根据最小二乘法原理
9、,记,就可以选择的估计使得由前面的分析,对所分的三类病人应该得到三个不同的拟合式。第一类:对于第一类,也就是病情较轻的病人,我们取(理由见后面的模型检验)拟合曲线如下图(拟合程序见附录:程序1) 分析图表可知,药物对该类病人的疗效呈现出一定的波动,开始时由于药物的作用HIV/CD4有所下降,但是经过一段时期后艾滋病病毒产生了耐药性,病情有所反弹。因为病人继续用药(甚至加大药量),这可以增加机体本身的免疫力,使得HIV/CD4仍然会下降。但全身免疫系统仍在受到艾滋病病毒的破坏,加之随着时间的增加药物不断在体内积累,可能对病人产生了副作用,直到免疫系统功能再也不能维持最低的防御能力时,病情会迅速恶
10、化,如图上所示。大约在43周后应该停止用药。第二类:我们取。曲线拟合图形如下:药物对于一般病人的疗效出现了与上类病人相类似的情况,但是其病情恶化的时间较前面提前了,大约在42周,此时应该停止治疗。第三类:仍取,曲线拟合图形如下:分析图表得到,药物在一开始疗效十分明显,以后渐趋平缓。由于病人的免疫系统十分脆弱,会对药物产生较大的依赖性,但是随着体内艾滋病毒耐药性的增加,病人的病情比其它两类病人更早的出现病情的反弹。大约在40周HIV/CD4达到最低并且开始恶化,所以在40周应该停止用药。问题二:为了对四种疗法进行评价,我们采用由美国著名运筹学家查恩斯,库伯和罗兹在最优生产前沿面的理论及“相对效率
11、评价”概念基础上进一步发展而来的,新的行之有效的系统分析方法12-数据包络(DEA)。数据包络(DEA)是使用数学规划模型比较决策单元之间的相对效率,对决策单元作出评价。我们采用多指标投入和多指标产出而对相同类型单元的相对效率进行综合评价的BCC模型。BCC模型的基本模型是符号说明: 第个决策单元对第种输入的数据, 第个决策单元对第种输入的数据,; 对第种输入的一种度量(或称权重) 对第种输入的一种度量(或称权重); 投入产品的松弛变量 输出产品的松弛变量用BCC模型进行求解,则投入变量是每种药品的药量和服药前患者的CD4数目,输出变量是服药后CD4的值,从而可以得出值,为对各种疗法进行有效的
12、分析,我们将各疗法数据按时间进行分类,相同时间点的作为一类;从而可以按照患者在不同的疗程内,对疗法进行比较。而用这种方法要处理的数据量太多,将数据用MATLAB进行处理(见附录1),得出在同一时刻不同疗法下的CD4浓度。将时间按4个值作为一个疗程,对疗法进行评价,在这里以CD4的浓度作为评价标准的话,要反映疗效好坏,我们观察在疗程内患者的CD4浓度的最大值(即患者的最好状态),及达到最佳状态的时刻和患者服药前的病情严重程度(即服药前的CD4浓度),将此三者做为投入量,将服药后CD4浓度作为输出量。在BCC模型中进行运算,得出,在相同疗程内,以CD4为评价指标的疗效;在lingo软件中输入程序(
13、见附录:程序3),min=p;3.4657*x1+4.1649*x2+3.239*x3+3.47*x4+s1=3.4657*p;2.97*x1+2.934*x2+2.90*x3+2.835*x4+s2=2.97*p;6*x1+4*x2+7*x3+5*x4+s3=6*p;2.36*x1+2.783*x2+2.63*x3+3.12*x4=2.36;得出p值,结果如下图所示:表2:输出结果疗程第一种疗法第二种疗法第三种疗法第四种疗法疗程一110.6451疗程二0.91840.96611疗程三10.90690.9511疗程四10.896510.97疗程五0.77120.994211疗程六0.73590
14、.796710.7968疗程七110.90490.96疗程八0.985710.90571疗程九0.87600.814710.9559结果分析:从结果来看,四种疗法没有明显的优劣之分,各有千秋。都有疗效比较好的时期。横向比较:在治疗初期,第三种药相对而言是疗效相对差一些,而其他三种的疗效旗鼓相当。随着疗程的进行,患者身体上产生了抗药性,而在这期间疗法三和疗法四的疗效明显好于其他两种疗法(需要画图这里!),特别是第六疗程,疗法三明显好于其他三种,而疗法三和疗法四,在治疗中期相差不大。在后期的治疗中,疗法一相对差一些,疗法二相对好一些,但是他们的相对系数差别不是很大,因而这里认为四种疗法无明显的差别
15、。表3:横向比较结果初期中期后期疗法四疗法一疗法二疗法三疗法二疗法四疗法三和疗法四疗法一疗法三疗法二疗法一注:中期疗法三和疗法四差别很小纵向比较:在整个过程中疗法四的波动较小,可以认为疗法四在本次实验中的大部分时间内相对其他三种疗效处于相对高效状态。综上所述,我们认为疗法四相对其他三种而言比较好。但是在不同的时期而言,有不同的最好的疗法。在Matlab中编程对疗法4进行曲线拟合预测(程序见附录4),运行结果,在图中画出为:对上图分析可得,疗法4的疗效比较平稳,可以用作长期的治疗,对维持CD4细胞的数目,起非常重要的作用.问题三:由于条件是在问题二基础上,所以大体分析和问题二相同,我们引入“成本
16、疗效比”这一概念,来反映在一定价格下疗法的疗效。现在通过数据包络(DEA)进行有效的分析,问题三中是根据主要供给商对不发达国家提供的药品价格,对疗法重新进行评价。那么,就要建立费用和疗效之间的关系,可以将费用作为投入量,以“成本疗效”比为评价标准,由于是将药品供给不发达国家,那么就要考虑价格和疗效之间的相对权值。因而采用数据包络中的模型进行处理,因为其在投入中考虑投入产品的相对权值。以价格,某疗程中CD4的最大值和取最大值此时的疗程数,为投入变量,以该疗程结束时CD4浓度为输出变量,进行求解。模型的基本模型是:其中: 第个决策单元对第种输入的数据, 第个决策单元对第种输入的数据,; 对第种输入
17、的一种度量(或称权重) 对第种输入的一种度量(或称权重); 将二中用MATLAB整理的数据代入模型中,在Lingo软件中编程:max=(3.3390*x1+38*x2+2.7755*x3+4.65*x4)/2.8644*p;2.677*x1+37*x2+2.374*x3+1.25*x4=2.345*p;3.0459*x1+35*x2+2.3869*x3+3.45*x4=2.2485*p;3.1542*x1+40*x2+2.514*x3+2.45*x4=3.154*p;3.3390*x1+38*x2+2.7755*x3+4.65*x4=2.8644*p;求出相应的效率评价指数:这里我们取的是的倒
18、数。在软件中运行的结果如下图所示表:在实验中的疗法的相对效率指数(单位)疗程名疗法一疗法二疗法三疗法四疗程一0.39520.49630.70710.4786疗程二0.6030.66280.60540.6404疗程三0.718990.84150.713360.6782疗程四0.7803990.9650.8615440.9595疗程五1.18131.147890.92741.0159疗程六1.581841.60391.15691.5784疗程七1.431471.34061.477491.4752疗程八1.5252971.618371.60941.521疗程九1.846521.95161.52561
19、.702433总和10.0640210.627969.58419410.04963 由于得到的值反映的是在某疗程下该疗法的相对效率指数,因而对其求和求出在整个实验中该疗法的相对疗效。结果分析:结果中疗法二的和最大,因而说明疗法二在问题所给的对不发达国家提供的价格下是相对最好的,虽然疗法四是疗效最好的,但是由于其价格昂贵,而影响了其相对疗效,不适合不发达国家患者服用。在开始时按照时间和CD4均值进行曲线拟合,得出疗法一是最好的,是因为疗法一中的药品价格比其他疗法有明显的价格优势,那么在“成本疗效比”中,即使疗效值较小,那么其“成本疗效比”也是一个较大的值,因而会在图象中显示出明显的优势。五、模型
20、检验与分析在问题一中采用多项式进行拟合时,首先要解决的问题有:模型的结构参数(多项式的阶次n)为多少?参数估计的计算估计参数的区间拟合的良好程度对于以上的问题可按以下步骤编程(附录:程序2)计算:计算原数据的标准差依次用 1 到n阶多项式去拟合计算拟合多项式系数用元胞数组记录不同阶次多项式的系数计算各系数的误差用元胞数组记录不同阶次多项式系数的误差记录自由度计算拟合多项式值的范围用元胞数组记录不同阶次拟合多项式的均值用元胞数组记录不同阶次拟合多项式的离差第一类拟合的分析:Q代表适当度,可以用于拟合的良好度的判断。它与自由度和量有关。Q值越大、越稳定说明拟合的程度越好。从上面的图可以看出,我们用
21、多项式进行的拟合还是令人满意的。上图反映了用七次的多项式进行拟合的情况,在误差允许的范围内也是可以接受的。下面关于第二类、第三类的拟合分析可以采用同样的方法。第二类的拟合分析: 第三类的拟合分析: 模型优缺点分析模型优点:数据包络(DEA)具有较高的灵敏度与可靠性,可以对无法价格化甚至难以轻易确定适当权重的指标进行分析,各测量指标能够以原来的面目出现,不必统一单位,无须进行量纲处理,大大简化了测量过程,保证了原始信息的完整;另一方面也避免了权重的人为因素影响。缺点:对数据误差与缺失十分敏感,在取中间数值时,我们要取小数点后四位,这无形中提高了调查数据收集的难度,它所评价的是相对疗效,而不是绝对
22、疗效。可能这些疗法都是很好的,只是程度不同而已。 参考文献1 Charmes A,Cooper W W,LiS.Using DEA to Evaluate Relative Efficiences in the Economic Performance of Chinese CitiesJ.Socio-Economic Planning Sciences,19892 Chilingerian JA.Evaluating physician deficiency in hospitals: A multivariate analysis of best practicesJ.European j
23、ouranal of operational research,1998,80:546-5483 韩 梅 DEA方法在国外医疗卫生系统效益评价中的应用J 系统工程理论与实践,1989,(2):55-684 陈祥华 邱枫林 王爱杰 数据包络分析在乡镇卫生院资源配置效率分析中的应用 中图分类号:R197。62 文献标识码:A 文章编号:1004-7778(2005)12-0023-02 5 陈仲生 基于Matlab7.0的统计信息处理 湖南出版社附录程序1:%三类:% 被拟合的原始数据a=0.9017956350.2218750.2355707070.2018689850.2335938140.0
24、614444440.179008430.0657404540.00250.05750.0919004210.3840464170.0252237580.0286165060.0333330.0271930580.0722191320.0387390.096491;y=a;x=0 3 4 5 6 7 8 9 12 22 23 24 25 27 40 41 42 44 46;dy=std(y); % 原数据 y 的标准差for n=1:10 % 依次用 1 到 6 阶多项式去拟合a,S=polyfit(x,y,n); % 计算拟合多项式系数An=a; % 用元胞数组记录不同阶次多项式的系数da=d
25、y*sqrt(diag(inv(S.R*S.R); % 计算各系数的误差DAn=da; % 用元胞数组记录不同阶次多项式系数的误差freedom(n)=S.df; % 记录自由度ye,delta=polyval(a,x,S); % 计算拟合多项式值的范围YEn=ye; % 用元胞数组记录不同阶次拟合多项式的均值Dn=delta; % 用元胞数组记录不同阶次拟合多项式的离差chi2(n)=sum(y-ye).2)/dy/dy; % 计算不同阶次的 量。endQ=1-chi2cdf(chi2,freedom); % 用于判断拟合良好度% 适当度的图示subplot(1,2,1),plot(1:10
26、,abs(chi2-freedom),b)xlabel( 阶次 ),title(chi 2 与自由度 )subplot(1,2,2),plot(1:10,Q,r,1:10,ones(1,10)*0.5)xlabel( 阶次 ),title(Q 与 0.5 线 ) %/clf,plot(x,y,b+);hold onerrorbar(x,YE7,D7,r);hold offtitle( 较适当的七阶拟合 )text(0.1,5.5,chi2= num2str(chi2(7) int2str(freedom(7)text(0.1,5,freedom= int2str(freedom(7)text(
27、0.6,1.7,Q= num2str(Q(7) 0.5)程序2min=p;3.4657*x1+4.1649*x2+3.239*x3+3.47*x4+s1=3.4657*p;2.97*x1+2.934*x2+2.90*x3+2.835*x4+s2=2.97*p;6*x1+4*x2+7*x3+5*x4+s3=6*p;2.36*x1+2.783*x2+2.63*x3+3.12*x4=2.36;程序3a4=2.835648788 0 2.786025 4 3.470433333 53.3322 6 3.037043478 7 3.186595455 83.326403333 9 3.125369231
28、 10 2.621671429 11 2.340525 12 4.01095 13 3.580975 14 3.2911125 15 3.329335659 16 3.208598387 172.842664286 18 2.873088889 19 2.527411111 20 2.7681 21 3.223083333 22 3.003425 23 3.055254118 242.964329268 25 3.230042857 26 2.34047 27 3.3670875 28 2.779344444 292.602575 303.282757143 312.905062069 323
29、.2299175 33 2.775527273 34 2.866623077 352.548366667 362.0138 37 3.33905 382.989107143 39 2.864474359 40;p4=polytool(a4(2,:),a4(1,:),9)程序4: a4=2.835648788 0 2.786025 4 3.470433333 53.3322 6 3.037043478 7 3.186595455 83.326403333 9 3.125369231 10 2.621671429 11 2.340525 12 4.01095 13 3.580975 14 3.29
30、11125 15 3.329335659 16 3.208598387 172.842664286 18 2.873088889 19 2.527411111 20 2.7681 21 3.223083333 22 3.003425 23 3.055254118 242.964329268 25 3.230042857 26 2.34047 27 3.3670875 28 2.779344444 292.602575 303.282757143 312.905062069 323.2299175 33 2.775527273 34 2.866623077 352.548366667 362.0
31、138 37 3.33905 382.989107143 39 2.864474359 40;a4(1,:)=a4(1,:)./4.65; p4=polyfit(a4(2,:),a4(1,:),9);x=0:.1:42;y=polyval(p4,x)plot(a4(2,:),a4(1,:),+,x,y) 2.956138889262.936635714272.819935714283.091422222292.53665 303.505565217312.752240964322.6478 332.514509524342.78212 352.903111111362.28418 372.91
32、665 382.182172727393.15425862140;a3(1,:)=a3(1,:)./2.45;%第四类a4=2.8356487880表1:各类病人不同时刻的统计表第一类:时间03456789HIV/CD40.90180.22190.23560.20190.23360.06140.17900.06571222232425274041420.00250.05750.09190.38400.02520.02860.03330.02720.072244460.03870.0965第二类:时间0234567HIV/CD40.07690.02110.02650.03410.02640.01
33、650.021089101116202122230.03500.02190.05470.01330.24440.00640.03230.06600.0208 2425262728293233370.01950.01800.04010.01550.03660.02930.04290.12260.011538394041424345460.04500.02200.01670.01700.04040.04310.02690.0295第三类:时间034578910HIV/CD40.40330.13650.15500.14490.05390.07650.19520.0218122123242527293
34、80.71430.03880.03980.05890.13670.06311.02000.129639404142440.06990.05080.01620.01220.0149表2:疗法一:CD40.85051.10341.20720.91500.83780.84970.8493时间02.00004.00005.00006.00007.00008.00000.88820.80690.81820.85520.80810.73100.81269.000010.000011.000012.000013.000014.000015.00000.86580.77210.75810.55990.8932
35、0.80510.737216.000017.000018.000019.000020.000021.000022.00000.80640.74760.76920.71910.60920.79090.475823.000024.000025.000026.000027.000028.000029.00000.78450.78490.76660.78490.74980.88290.643930.000031.000032.000033.000034.000035.000036.00000.67380.83300.81240.651737.000038.000039.000040.0000疗法二:C
36、D40.85051.10341.20720.91500.83780.84970.8493时间02.00004.00005.00006.00007.00008.00000.88820.80690.81820.85520.80810.73100.81269.000010.000011.000012.000013.000014.000015.00000.86580.77210.75810.55990.89320.80510.737216.000017.000018.000019.000020.000021.000022.00000.80640.74760.76920.71910.60920.7909
37、0.475823.000024.000025.000026.000027.000028.000029.00000.78450.78490.76660.78490.74980.88290.643930.000031.000032.000033.000034.000035.000036.00000.67380.83300.81240.651737.000038.000039.000040.0000疗法三:CD40.60980.59910.74630.71660.65310.68530.71540.6721时间04.00005.00006.00007.00008.00009.000010.00000
38、.56380.50330.86260.77010.70780.71600.69000.611311.000012.000013.000014.000015.000016.000017.000018.00000.61790.54350.59530.69310.64590.65700.63750.694619.000020.000021.000022.000023.000024.000025.000026.00000.50330.72410.59770.55970.70600.62470.69460.596927.000028.000029.000030.000031.000032.000033.
39、000034.00000.61650.54800.43310.71810.64280.616035.000036.000037.000038.000039.000040.0000疗法四:CD40.60980.59910.74630.71660.65310.68530.71540.6721时间04.00005.00006.00007.00008.00009.000010.00000.56380.50330.86260.77010.70780.71600.69000.611311.000012.000013.000014.000015.000016.000017.000018.00000.61790.54350.59530.69310.64590.65700.63750.694619.000020.000021.000022.000023.000024.000025.000026.00000.50330.72410.59770.55970.70600.62470.69460.596927.000028.000029.000030.000031.000032.000033.000034.00000.61650.54800.43310.71810.64280.616035.000036.000037.000038.000039.000040.0000