艾滋病疗法的评价及疗效的预测
【摘要】本问题要求我们利用艾滋病服药治疗期间测试所得的数据预测继续治疗的效果并评价各种疗法的优劣。对三个问题我们都运用回归分析方法和方差分析方法来建立统计模型。
就问题一，首先对所有数据进行处理，建立了一元多项式回归模型，将多项式回归模型转化为多元线性回归模型运用Matlab软件求解，通过最小二乘法参数估计和统计检验，得到关于与测试时刻的比较准确的函数关系及其回归曲线，由此预测出药物治疗的效果，最终得出在药物治疗的第25.0000周停止用药最为合适的结论。然后提出分类的思想，即按不同病人在服药前期的CD4浓度大小不同将数据分成三组，对每组数据同样按前面的思路建立回归模型并求解，确定最佳治疗终止时间依次为第25.1295周，第24.0000周，第25.0000周。分类前后的预测结果大致相符。
问题二在增多了病人年龄因素的基础上对四种疗法的优劣进行了评价。首先，对每组治疗方案建立4种多元回归子模型进行定性分析，利用Matlab统计工具箱来处理，拟合出交互式图形，并分析剩余标准差s，筛选出每组治疗方案最佳模型。然后分析CD4浓度曲线，对4种疗法的优劣进行了排序，选择疗法四为最终方案，它的最佳终止服药时间为第17.8854周。然后用方差分析进行定量评价。定性和定量评价的结果相符。
问题三在问题二基础上增加了对治疗费用的考虑，建立多元回归模型，同样利用Matlab统计工具箱来处理。然后提出“疗效性价比”这一概念并建立模型，以单位费用的治疗的效果来评价4种治疗效果，在交互式图形中分别取治疗费用为0和200美元的两个时刻作为观测点，分析4种疗法的疗效性价比，最后得到第1，2，3，4种疗法依次较优，即疗法四为最优方案，在预测继续治疗的效果不佳的情况下，确定最佳终止时间为第17.8854周，终止治疗时共花费453.2464美元。
关键词：回归分析 治疗效果预测 方差分析 疗效性价比
1
一、问题的重述
艾滋病是当前人类社会最严重的瘟疫之一，是由艾滋病病毒(英文简称HIV)引起的，这种病毒破坏人的免疫系统，使人体丧失抵抗各种疾病的能力，从而严重危害人的生命。人类免疫系统的CD4细胞在抵御HIV的入侵中起着重要作用，当CD4被HIV感染而裂解时，其数量会急剧减少，HIV将迅速增加，导致AIDS发作。
艾滋病治疗的目的，是尽量减少人体内HIV的数量，同时产生更多的CD4，至少要有效地降低CD4减少的速度，以提高人体免疫能力。现在得到了美国艾滋病医疗试验机构ACTG公布的同时服用3种药物疗法和另外一组4种疗法的相关数据。
要求： 1）、预测继续治疗的效果，或者确定最佳治疗终止时间。 2）、评价4种疗法的优劣（仅以CD4为标准），并对较优的疗法预测继续治疗的效果，或者确定最佳治疗终止时间。
3）、如果病人需要考虑4种疗法的费用，找出对2）中的评价和预测（或者提前终止）有什么改变。
二、模型的假设
1、附件1、2数据所反映的规律基本符合AIDS病人药物治疗过程的病况规律，即能客观的反映现实情况。
2、在测试和预测期间AIDS病人不会死亡。 3、病人在服药测试初始CD4浓度相等。 4、问题三第一种疗法中第奇数个月（如第1个月）服用400mg didanosine （0.85美元/天），第偶数个月（如第2个月）服用600mg zidovudine（1.60美元/天）。 5、问题三治疗费用只由治疗药品费用组成。 6、每月按30天计算。
三、符号说明
符号 c h 表示意义 病人的CD4浓度（0.2个/ml） 病人的HIV浓度（单位不详） 第i类病人的CD4浓度（0.2个/ml），i?1,2,3 第i类病人的HIV浓度（单位不详），i?1,2,3 测试CD4和HIV浓度的时刻（单位：周） 问题二、三中第j种疗法的时刻（单位：周），j?1,2,3,4 对数处理后第j种疗法病人体内CD4浓度，j?1,2,3,4 ci hi d dj yj 2
符号 表示意义 第j种疗法的AIDS病人年龄，j?1,2,3,4 第j种疗法的拟合曲线的斜率的绝对值，j?1,2,3，4 第j种疗法的费用（单位：美元），j?1,2,3,4 停止治疗时为奇数月份时总的治疗费用 停止治疗时为偶数月份总的治疗费用 采用第j种疗法时费用增量，j?1,2,3,4 费用增量为?fj时CD4的浓度增量j?1,2,3,4 疗效性价比，j?1,2,3,4 aj kj fj s1 s2 ?fj ?yj pj
四、问题的分析、模型的建立及求解
4.1
问题一的分析及模型建立：
4.1.1 问题一的提出和数据分析
已知：同时服用3种药物（zidovudine，lamivudine，indinavir）的300多名病人每隔几周测试的CD4和HIV的浓度（每毫升血液里的数量）。
要求：预测继续治疗的效果，或者确定最佳治疗终止时间（继续治疗指在测试终止后继续服药，如果认为继续服药效果不好，则可选择提前终止治疗）。
分析：题目的附件1给出的CD4和HIV浓度是300多名病人每隔几周测试的，每个病人测试3到6次不等，加起来共有1665组数据。根据题目要求，我们应当寻找出药物治疗过程中CD4和HIV浓度随治疗时间变化而变化的规律图线，以此来判断和预测。
首先，以测试时刻（周）为横轴，分别以CD4、HIV的浓度为纵轴，作出浓度关于时刻的散点图进行观测，由于数据点太多，测试时刻相对较集中且每个病人测试次数只有几次，从散点图上不容易看出它们之间的规律。于是，再用回归分析法拟合出浓度和时刻的关系曲线并做统计检验。仔细观察分析数据我们发现，不同的病人在初始时刻CD4浓度（或HIV浓度）有所不同，即每个病人的病程是不一样的，于是我们提出按服药前（第0周或者第1周）CD4浓度的不同将病人分为三类，如表1所示：
3
表1 以CD4浓度为标准分类 第pi类病人 CD4浓度（0.2个/ml） p1 0－75 p2 75－150 p3 >150 然后对每一类数据各自分析讨论，同样用回归分析法拟合出浓度和时刻的关系曲线并做统计检验。以上两种方法同时用来解决问题一（即预测继续治疗的效果，确定最佳治疗终止时间），再比较两种方法的优劣。
4.1.2 模型的建立及求解
（一）模型一：分类前的回归模型
首先，记病人的CD4浓度（0.2个/ml）为c，HIV浓度（单位不详）为h（为了便于在同一个图像上作出CD4和HIV浓度与时刻的回归曲线，在数据时
d为测试CD4和HIV浓度的时刻均将HIV浓度数据扩大了50倍），（单位：周），
由附件1初步分析，随着服药时间的增加，CD4浓度先是增加，但是随着时间增加过多后c反而减少，而HIV浓度h先是减少然后又增加，所以不应拟合线性函数，可能是二次函数比较合适。
建立如下CD4浓度的回归模型：
c??0??1d??2d2?? （Ⅰ）
h??0??1d??2d2?? （Ⅱ）
其中?0～?2为回归系数, ?是相互独立的、期望为0、方差为?2、正态分布的随机变量，即??N(0,?2),?称（随机）误差。
1） 利用最小二乘法进行参数估计，确定系数?0～?2，得到预测方程； 2） 给定显著水平?＝0.05，对回归系数进行相关系数检验及F检验，确定回归
分析的效果，直到回归分析的效果好。
用Matlab软件命令regress求解,将回归系数的估计值如表2所示，并代入（Ⅰ）、（Ⅱ）式得到CD4浓度c和HIV浓度h的预测方程为
c?97.7277?6.6004d?0.1114d2 （Ⅲ）
h?439.7715?16.2617d?0.3289d2 （Ⅳ）
4
表2 模型一CD4和HIV浓度的计算结果 回归系数 回归系数估计值 CD4浓度 97.7277 6.6004 -0.1114 HIV浓度 439.7715 -16.2617 0.3289 回归系数置信区间 CD4浓度 [89.3095 106.1458] [5.2526 7.9481] [-0.1452 -0.0776] CD4浓度 [429.4523 450.0907] [-17.9164 -14.6071] [0.2874 0.3705] ?0 ?1 ?2 CD4浓度 HIV浓度 R2＝0. 1113 F＝98.1 p＝0 s2＝9990.4 R2＝0.2401 F＝247.3568 p＝0 s2＝1.5001×104 HIV and CD47006005004003002001000-100-10
由预测方程（Ⅲ）（Ⅳ）及图1上的回归曲线可以看出c是开口向下的抛物线而h是开口向上的抛物线，这符合AIDS患者用药后体内CD4浓度和HIV浓度随时间变化的规律。
回归系数?0～?2的置信区间比较短，而且都没有包含零点，表明回归模型有效。
HIV*50 or CD4? HIV? CD4010203040date(week)50607080图1 HIV浓度×50和CD4浓度的回归曲线
回归模型的统计检验 检验水平为?＝0.05 1） 相关系数检验法：相关系数R是反映随机变量X与Y呈线性程度的
一个度量指标，其取值范围R?1，当R接近于1时，变量X与Y呈密切线性相关，当R接近于0时，变量X与Y非线性相关。本模型CD4浓度和HIV浓度的计算结果中R2分别为0.1113和0.2401，都明显偏小，可见HIV和CD4浓度与检测时刻呈非线性相关，与我
们的模型假设是相符的。
5
2）
F检验法：若F?F1??(1,n?2)则回归效果好。用Matlab求得＝F0.05(1,1567)3.8474，由表二可看出，F的值都明显大于3.8474，表明回归方程的效果显著。
问题结果及分析
由图1得到曲线中CD4浓度最大值所对应的时刻为29.2317（周），HIV浓度最小值所对应的时刻为25.0000（周），也就是病人服药治疗期间，在第25.0000周体内HIV浓度降到最低，继续服药则HIV浓度又会反弹上升；同时CD4浓度还在上升，到了第29.2317周才上升到最大值，之后呈下降趋势。由于艾滋病治疗的目的，是尽量减少人体内HIV的数量，提高CD4的数量也是为了增加HIV，为了防止HIV浓度上升导致病情的恶化，在HIV浓度降到最低且尚未反弹时就应当停止服药。由此可以确定最佳治疗终止时间也就是HIV浓度降到最低的时间。
模型一得到的最佳治疗终止时间为第25.0000周。
（二）模型二：分类后的回归模型
将不同病人在初始时刻CD4浓度的不同进行分类后得出的数据分为三组，分别得到M1, M2, M3组矩阵。对每组矩阵数据再建立关于HIV、CD4浓度与测试时间的二次函数：
ci??0??1di??2di2??(i?1,2,3) （Ⅴ）
2hi??0??1di??2d?2, 3 ) （Ⅵ） i??(i1,用Matlab软件命令regress求解,将回归系数的估计值代入（Ⅴ）（Ⅵ）式得到CD4浓度ci和HIV浓度hi的预测方程
表3 回归模型表
浓度类型 数据类型 预测方程 R2 F Fa s2 M1 CD4浓度 c1?43.5456?5.9858d1?0.0906d12 2 c2?129.8522?7.0061d2?0.1241d20.1997 98.944 3.853 5351.5 0.1599 43.5 3.862 6734.1 M2 M3 2 c3?196.26?6.912d3?0.1132d30.1550 26.038 3.874 8083.6 HIV浓度 M1 M2 h1?458.1105?15.0890d1?0.3010d12 0.2029 100.91 3.853 16387 2 0.2963 96.229 3.862 1285.9 h2?425.5908?17.9384d2?0.3761d2 6
M3 2 h3?411.35?17.595d3?0.3466d30.3670 82.319 3.874 10527 p= 0 < 0.0001
HIV and CD47006005004003002001000-100-10HIV and CD4600500400HIV*50 or CD4HIV*50 or CD4300? HIV? HIV200? CD4100? CD40010203040date(week)50607080-100-10010203040date(week)50607080
图2 第一类回归曲线 图3 第二类回归曲线
HIV and CD4700600500400HIV*50 or CD4
? HIV3002001000-100-200-10? CD4010203040date(week)50607080
观察图2～图4，每一类CD4浓度和HIV浓度的回归曲线的变化趋势都大致相同，CD4浓度曲线都是先增大后减少，HIV浓度曲线都是先减少后上升。
图4 第三类回归曲线
问题结果及分析
从上面三个回归曲线图可以得到CD4浓度曲线最小值和HIV浓度曲线最大值所在的时刻如表四所示，如同模型一的分析，在HIV浓度降到最低且尚未反弹时就应当停止服药。由此可以判断不同类型病人都不宜在测试结束后还继续治疗，而是在服药到一定时间后终止，并且确定出最佳治疗终止时间分别为第25.1295周，第24.0000周，第25.0000周。
7
表4 曲线驻点的横纵坐标值 曲线类型 第一类 浓度值 时刻 33.0000 25.1295 第二类 浓度值 228.7633 106.0407 时刻 28.0000 24.0000 第三类 浓度值 301.0769 94.1240 时刻 30.0769 25.0000 CD4浓度曲线（min） 124.4511 HIV浓度曲线（max） 134.5057
（三）模型一与模型二对所得结果优劣评价
模型一中是对整体数据开始处理，初步得到治疗效果的评价，对所有病人都适用；模型二中是逐步细化将病人分类后得到三种评价结果，工作量稍微大点，但可作为模型一的参照检验，两种评价相差不大，这在一定程度上检验了结果的正确性。两种处理方法各有优缺点，相互补充。 4.2
问题二的分析及模型建立：
4.2.1 问题二的提出和数据分析
已知：将1300多名病人随机地分为4组，每组按下述4种疗法中的一种服药，大约每隔8周测试的CD4浓度。
四种疗法的日用药分别为：
疗法一：600mg zidovudine或400mg didanosine，这两种药按月轮换使用； 疗法二：600 mg zidovudine加2.25 mg zalcitabine； 疗法三：600 mg zidovudine加400 mg didanosine；
疗法四：600 mg zidovudine加400 mg didanosine，再加400 mg nevirapine。 要求：以CD4为标准评价4种疗法的优劣，并对较优的疗法预测继续治疗的效果，或者确定最佳治疗终止时间。
分析：由于采用了四种不同的用药方案，且不同年龄的病人身体免疫能力、产生CD4细胞的能力各有不同，这将直接影响CD4细胞浓度的变化趋势。由此得出影响CD4细胞浓度的因素有两个——病人的年龄差异和治疗方案的不同（具体表现为治疗的时间的长短），问题可归结为多元回归模型。采用多元回归来拟合出CD4的浓度与病人的年龄和治疗的时间的交互式图形，曲线的变化趋势表达了人体内CD4的浓度的变化趋势，曲线的斜率表达了CD4减少或增加的速率，进而根据据此判断治疗方案的优劣及确定最佳治疗时间。 4.2.2 模型的建立及求解 （一）建立多元回归模型：
4CD4的浓度做简化处理，处理为将治疗方案依次记为j?1,2,3,，
，并记做yj?log?CD4count?1?，将第j种治疗方案的病人的log??1?CD4count年龄记为aj(j?1,2,3,4)，治疗的时间记为dj?j?1,2,3,4?。
由于二元回归方程存在多种形式，为更加精确的拟合出变化曲线，我们采用
多元回归模型同时建立多个子模型拟合出CD4的浓度分别与病人的年龄和治疗的时间的交互式曲线，然后筛选出最佳模型作为它们各自的最终模型，作出最佳
8
模型对应的关系曲线[1]。
1.）purequadratic模型：包含线性项和纯二次项
2yj??0??1aj??2dj??3a2(j?1,2,3,4) j??4dj??2.）quadratic模型：包含线性项和完全二次项
2yj??0??1aj??2dj??3ajdj??4a2(j?1,2,3,4) j??5dj??3.）interaction模型：包含线性项和纯交互项
yj??0??1aj??2dj??3ajdj??4.）linear模型：只包含线性项
(j?1,2,3,4)
yj??0??1aj??2dj??(j?1,2,3,4)
（二）模型的求解：
用MATLAB统计工具箱中的rstool命令来处理，每个治疗方案都对应可输出四个交互式画面，其中每个交互式画面都在图左面提供了两个下拉式菜单，以便我们选择不同的子模型来分析曲线。将每种治疗方案的四个模型输出的回归系数?和剩余标准差s分别列入表5～表8中：
表5 治疗方案一与CD4浓度的模型输出
子模型名称 purequadratic quadratic interaction linear
表6 治疗方案二与CD4浓度的模型输出
子模型名称 purequadratic quadratic interaction linear
表7 治疗方案三与CD4浓度的模型输出
子模型名称 purequadratic quadratic interaction linear ?0 2.8293 2.8747 2.9047 2.8665 ?1 0.0041 0.0027 0.0018 0.0028 ?2 -0.0111 -0.0137 -0.0170 -0.0144 ?3 -0.0000 0.0001 0.0001 ?4 -0.0001 -0.0000 ?5 -0.0001 s 0.9310 0.9313 0.9307 0.9303 ?0 1.4525 1.5092 2.1844 2.1252 ?1 0.0565 0.0548 0.0212 0.0228 ?2 -0.0060 -0.0092 -0.0158 -0.0120 ?3 -0.0004 0.0001 0.0001 ?4 -0.0002 -0.0004 ?5 -0.0002 s 1.0321 1.0325 1.0326 1.0322 ?0 2.2972 2.3112 2.7134 2.6930 ?1 0.0271 0.0267 0.0073 0.0079 ?2 0.0041 0.0031 -0.0065 -0.0051 ?3 -0.0002 0.0000 0.0000 ?4 -0.0003 -0.0002 ?5 -0.0003 s 1.1215 1.1220 1.1220 1.1215 9
表8 治疗方案四与CD4浓度的模型输出 子模型名称 purequadratic quadratic interaction linear ?0 3.0105 3.0732 2.6306 2.5666 ?1 -0.0181 -0.0195 0.0099 0.0116 ?2 0.0331 0.0285 -0.0032 0.0008 ?3 0.0004 0.0001 0.0001 ?4 -0.0009 0.0004 ?5 -0.0009 s 1.1390 1.1394 1.1464 1.1460
筛选最佳模型：分析它们各自的剩余标准差s，选取s最小的一个，如有两项相同且最小，即取较为简单的一个方程，得到每种治疗方案的最终模型，
治疗方案一选取linear模型：
y1?2.8665?0.0028a1?0.0144d1 治疗方案二选取purequadratic模型：
22 y2?1.4525?0.0565a2?0.0060d2?0.0004a2?0.0002d2治疗方案三选取linear模型：
y3?2.6930?0.0079a3?0.0051d3 治疗方案四选取purequadratic模型：
22 y4?3.0105?0.0181a4?0.0331d4?0.0004a4?0.0009d4四种治疗方案按各自对应的最终模型输出的交互式画面如表5至表8下：
图5 疗法一的交互式画面 图6 疗法二的交互式画面
10
图7 疗法三的交互式画面 图8 疗法四的交互式画面
注：以上四图中虚曲线为因变量的置信区间，实曲线为自变量与因变量的关系曲线。
图形分析：
图形最左边给出了yj的预测值及其置信区间。
交互式画面左边一幅图形表示治疗时间dj固定时CD4浓度yj随病人年龄aj变化的变化曲线。
四幅图形中左半部分曲线均呈增长的总趋势，说明在经过同一治疗时间后病人体内的CD4浓度会随病人年龄的增长而增加。
交互式画面右半边一幅图形表示病人年龄aj固定时CD4浓度yj随治疗时间
dj变化而变化的曲线。
（三）四种疗法优劣的评价及预测结果
在图5～图8中，在年龄一定的情况下，只有第四种疗法的CD4浓度的变化曲线随治疗时间的延迟而产生先增高后降低的趋势，其他各疗法随着治疗的进行CD4的浓度在持续降低，由于艾滋病治疗的目的，是尽量使人体内产生更多的CD4，有效地降低CD4减少的速度，以提高人体免疫能力。所以第四种疗法最好，且最佳治疗终止时间应选在图形的最高点处，即确定d4?17.8854周为最佳治疗终止时间。
分析其他三种疗法曲线的走势，随着治疗的进行它们的CD4的浓度在持续降低，在此曲线斜率的绝对值表达了CD4降低的速率，斜率的绝对值越大，CD4减少得越快。设kj(j?1,2,3)为第j种疗法曲线斜率的绝对值，存在以下关系：
k3?k2?k1，说明随着治疗的进行，第三种疗法CD4减少的速度最慢，第二种疗法次之，第一种疗法减少的速度最快，从而确定优劣：疗法三>疗法二>疗法一。
综上可知，四种疗法的优劣依次如下（表9）（编号由优至劣如1号为最好，4号为最差）：
11
表9 四种疗法优劣比较结果 优劣序号 1. 2. 3. 4. 疗法名称 疗法四 疗法三 疗法二 疗法一 最佳治疗终止时间 第17.8854周 即选择疗法四：600 mg zidovudine加400 mg didanosine，再加400 mg nevirapine为最终方案，它的最佳终止时间为第17.8854周。
（四）方差分析
以上的回归分析是定性的判断四种疗法的优劣，下面我们来定量的讨论和检验四种不同的疗法对CD4浓度的影响，我们采用单因素方差分析方法处理。
由于附件2是将1300多名病人随机地分为4组，因此在此可以假设所有病人的CD4的初始浓度是相等的。由以上的回归模型得出四种中的前三种的CD4浓度几乎都是随治疗时间的增长而降低，只有第四种疗法的第17周附近出现了转折，因此对所有病人取具有代表性的第14至19周的CD4浓度作为观测值，再将所得观测值分别按对应的疗法分组，最后在Matlab里用命令anova1(X,group)对各疗法的优劣进行单因素方差分析并比较出优劣顺序。M文件见附录data203.m，结果如下[4]：
P?4.23923?10?5
并且得到两个图形界面（图9、图10）：
图一方差分析表中有6列: 第1列显示误差的来源；
第2列显示每一个误差来源的平方和(SS)；
第3列显示与每—个话中有误差来源相关的自由度(df)；
第4列显示均值平方和(MS)，它是误差来源平方和自由度的比值,即SS/df; 第5列显示F统计量，它是均值平方和的比值；
第6列显示p值,p值是F的函数(fcdf)；当F增加时P值减小。
图9 单因素方差分析表
图10 box图显示X的每组数据的箱形图。箱形图中心线上较大的差异对应于较大的F值和较小的P值。
12
因
P?4??5为
F0.05(2,930)= 3.0054,由方差分析表得F＝7.69 > F0.05(2,930),所以认为不同的疗法对病人CD4浓度的提高有显著影响。 各水平的效应值如下： ^??图10 ?1?x1?x＝2.8112－2.9883＝-0.1771 ?2?x2?x＝2.8597－2.9883＝-0.1286 ?3?x3?x＝3.0147－2.9883＝0.0264 ^??^???4?x4?x＝3.2566－2.9883＝0.2683 计算结果表明，效应值?4最大。说明第四种疗法对病人CD4浓度的提高值最大，因此第四种疗法?4为最优疗法。 4.3 问题三的分析及模型建立： ^^^??4.3.1 问题三的提出和数据分析 已知：艾滋病药品的主要供给商对不发达国家提供的各种药品价格不尽相同。具体如下：600mg zidovudine 1.60美元/天，400mg didanosine 0.85美元/天，2.25 mg zalcitabine 1.85美元/天，400 mg nevirapine 1.20美元/天。 要求：若病人需考虑4种疗法的费用，对（2）中的评价和预测（或者提前终止）将产生什么改变。 分析：在（2）中我们在讨论疗法的优劣及预测较优疗法继续治疗的效果时，只以CD4的浓度作为参考标准，在此应综合考虑治疗费用和CD4的浓度大小值，首先列出治疗费用与治疗时间的函数关系，然后采用多元回归模型来拟合出CD4的浓度与病人的年龄和治疗费用的交互式图形，再建立一个疗效性价比模型来进行分析。 在此提出一个疗效性价比的概念，定义： 疗效性价比＝CD4降低的浓度 / 所花治疗费用 13 为简化分析，我们采用比较在花费相同费用时四种疗法CD4浓度的降低程度来，判断四种疗法的优劣并做预测。 4.3.2 模型的建立及求解 （一）建立模型： 将第j种疗法的治疗费用记为fj(j?1,2,3,4)，第j种疗法的疗效性价比记为pj(j?1,2,3,4)，采用第j种疗法时费用增量为?fj(j?1,2,3,4)时CD4的浓度增量记为?yj(j?1,2,3,4)，治疗时间记为dj?j?1,2,3,4?。在采用第一种疗法时，记停止治疗时为奇数月份时总的治疗费用为s1，记停止治疗时为偶数月份总的治疗费用为s2。记： ?1该月为治疗的第奇数个月份?1该月为治疗的第偶数个月份，w2??。 w1???0该月为治疗的第偶数个月份?0该月为治疗的第奇数个月份第一步：根据病人日用药量、药品使用方案及药品价格确定治疗费用 fj(j?1,2,3,4)与治疗时间dj?j?1,2,3,4?的线性关系如下： 为使治疗费用尽可能少，假设采用第一种疗法时第一个月服用的是每400mg 0.85美元的didanosine，第二个月服用的是每600mg 1.60美元的zidovudine 。 ?7d??7d?)/2?)?0.85??(1)/2??30?1.60 1.）s1?(7d1?30??(13030向上取整向上取整??向下取整??向下取整?7d??7d?s2?(7d1?30??(1)/2?)?1.60??(1)/2??30?0.85 ?30向上取整?向下取整?30向上取整?向下取整w1?(7d1)/2 30向上取整求余w2?1?w1 f1?s1?w1?s2?w2 2．）f2?7d2?(1.60?1.85) 3．）f3?7d3?(1.60?0.85) 4．）f4?7d4?(1.60?0.85?1.20) 第二步：为简化计算，将CD4的浓度处理为log?CD4count?1?，并记做 yj?log?CD4count?1?。 14 我们以病人的年龄aj和病人的治疗费用fj作为自变量，以CD4的浓度yj作为因变量，利用多元回归分析拟合出CD4的浓度yj与病人的年龄aj和病人的治疗费用fj的交互式图形，输出每种治疗方案的四个模型输出的回归系数?和剩余标准差s，选取s最小的一个方程，如有两项最小，则取方程较为简单的一个，得到每种治疗方案的最终模型。（具体做法同问题二的处理）。 各种疗法的模型如下[1]： 治疗方案一选取purequadratic模型： y1?3.0105?0.0181a1?0.0013f1?0.0004a1f1 治疗方案二选取quadratic模型： y2?1.4525?0.565a2?0.0002f2?0.0004a2f2 治疗方案三选取linear模型： y3?2.8305?0.0042a3?0.0013f3 治疗方案四选取quadratic模型： y4?2.2972?0.0271a4?0.0002f4?0.0002a4f4 四种治疗方案对应的最终模型的交互式画面如下： 图11 治疗方案一的交互式画面（左：费用为0美元时；右：费用为200美元时） 图12 治疗方案二的交互式画面（左：费用为0美元时；右：费用为200美元时） 15 图13 治疗方案三的交互式画面（左：费用为0美元时；右：费用为200美元时） 图14 治疗方案四的交互式画面（左：费用为0美元时；右：费用为200美元时） 图15 治疗方案三、四CD4浓度达最大值的交互式画面（左：方案三；右：方案四） 第三步：根据疗效性价比的定义，建立一个关于四种不同疗法的疗效性价比模型：pj??yj?fj(j?1,2,3,4) 花费同样多的治疗费用，疗效性价比越高，单位费用的治疗的效果越明显，则治疗效果越理想，越明显。以此作为判断四种疗法优劣的依据。 第四步：判断治疗效果 16 首先分析四种疗法的交互式画面的左边，当治疗费用fj固定时，病人年龄dj与病人体内CD4浓度yj的变化曲线均呈增长的总趋势，说明在经过同一治疗时间后病人体内的CD4浓度会随病人年龄的增长而增加。 然后分析四种疗法的交互式画面的右边，前两种疗法的图形均呈下降趋势，第三、四种疗法的图形呈抛物线状，故当选择第三、四种疗法时，应在CD4浓度尚未达到最高或者达到最高的时刻停止治疗。第三、四种疗法CD4浓度达到最高值时费用分别为：126.6599和453.24美元，故至少应在费用达到453.24美元时停止治疗。 在治疗费用在0～453.24美元上分析： 1．读图分析：疗效性价比的正负； 2．分别取治疗费用为0和200美元的两个时刻（读图9～图12），记录CD4 的浓度，并计算这个区间内的疗效性价比。 表10 取0、200美元分析疗效性价比 方案名称 方案一 方案二 方案三 方案四 疗效性价比的正费用为0美元时负（读图判断） 刻CD4浓度 负 负 先正后负 先正后负 2.971 2.988 2.969 2.848 费用为200美元时刻CD4浓度 2.6306 2.926 2.979 3.015 疗效性价比（?10） -17.020 -3.100 0.500 0.835 ?4注：疗效性价比的正负要求图形上每一点都满足，而不是只有某两点符合。 注：图表中每种的第一、二个柱状图表示的分别为治疗费用为0、200美元时CD4的浓度（无单位），第三个柱状图表示的是此时的疗效性价比（?10）。 ?450-5-10-15-20CD4浓度（200美元）疗效性价比CD4浓度（0美元）方案一方案二方案三方案四图16 取0、200美元分析疗效性价比的直方图 分析：读直方图可知，治疗方案三、四优于方案一、二。 而对于前两种治疗方案，很明显方案一的疗效性价比要比方案二的小的多。 17 故第二种方案优于第一种方案。 在图上，疗效一和疗法二疗效性价比为负，说明采用这两种治疗方案治疗，不但没有使CD4的浓度升高，反倒使其降低，药物治疗起了负作用，或者带来负面影响。 以下在第三、四种疗法CD4浓度达到最高值时分别做疗效性价比分析。 表11 取最值点分析疗效性价比 费用为0美元时刻CD4浓度 2.969 2.848 费用为126.6599美元时刻CD4浓度 2.9847 疗效性价比（?10） 1.240 ?4方案名称 费用为453.24美元时刻CD4浓度 3.1428 疗效性价比（?10） 6.504 ?4方案三 方案四 由于在两种疗法的各自的最值点处疗效性价比6.504>1.240，故第四种方案优于第三种方案。 第五步：四种方案的优劣比较结果及相关信息，如表12所示。 表12 方案优劣表 方案优劣顺序 1. 2. 3. 4. 方案名称 方案四 方案三 方案二 方案一 最优方案终止治疗时间（周次） 17.8854 最优方案终止治疗时共花费治疗费用（美元） 453.2464 选择疗法四：600 mg zidovudine加400 mg didanosine，再加400 mg nevirapine为最终方案，最佳终止时间为第17.8854周，终止治疗时共花费453.2464美元。 五、模型的评价 5.1 模型的优点： 1．分析附件1 CD4和HIV浓度数据的时候，先从整体数据入手，对所有数据进行统一处理，然后再细化，将数据进行分类处理，由两种方法分别预测最佳治疗终止时间。这种逐步处理方法简洁易懂，给出不同的解法，使我们的评价预测具有很强的可比性和可读性。 2．评价问题二4种疗法的优劣时，建立了回归模型并用Matlab统计工具箱处理，定性的作出评价，同时我们还运用方差分析作了定量评价，两种评价得出的结果是一致的。 3．问题二三中利用Matlab统计工具箱命令来处理，建立了多个多元回归子模型，在交互式画面上选择变量，这样逐步筛选，使模型达到最优，提高了精确度。 18 4. 在问题三的解答中，我们提出了一个“性价比”的概念，建立“疗效性价比”模型，使抽象问题具体化、量化，既有效的完美的结合了CD4的浓度和治疗费用，又便于计算和比较。 5.2 模型的缺点： 由于所给两组都数据非常庞大，考虑到时间问题，为了避免复杂而繁琐的工作量，在作回归分析时，我们都没有对残差图中的一些异常点进行剔除处理。但是回归分析得到的结果都经过了各种统计检验，因而总的回归效果还是很好的。 六、模型的改进和推广 1．对问题作回归分析并建立回归模型时，若进行残差分析，作出残差及其置信区间分析图，对原始数据中的异常数据剔除处理后再计算，模型会有一定的改进。 2．我们在附件1中随机抽取一些病人数据分成几组，作出CD4浓度（或HIV浓度）与测试时刻的关系折线图，可以看到，折线图大致有两种不同的走势，一种是CD4浓度随着时间的增加而呈现上升的趋势，最后浓度稳定在最高点，另一种是CD4浓度随着时间的增加先是上升，过一段时间后转而下降。这说明不同的病人进行药物治疗后大致有两种不同的效果，一种是服药后病情都在好转，另一种是服药后前一段时间有效果，随后又变坏。若对问题一继续作深入讨论，可以对附件1 CD4和HIV浓度数据作更合理的调整分类，分类依据是药物治疗后CD4浓度变化趋势，将这两类病人用聚类分析的方法分离出来，得到两组相对比较类似的数据，再分别进行回归分析处理后，各自会预测出的治疗效果应该会不同。前一种预测应该是建议继续治疗，即在测试终止后继续服药，后一种预测应该是在测试期间当病人的CD4浓度值上升到最大时即终止服药。这种预测结论更加贴近现实。由于时间紧迫，我们没有具体进行操作，在此只作为模型的推广大致讨论了一下。 参考文献 [1] 姜启源，邢文训，谢金星,杨顶辉 大学数学实验，北京：清华大学出版社，2000 [2] 马知恩 等著，传染病动力学的数学建模与研究，北京:科学出版社，2004 [3] 傅鹂，龚劬，刘琼荪，何中市编著,数学实验，北京:科学出版社，2000 [4] （美）埃维森（Gudmund R.Iversen）等著；吴喜之 等译,统计学,北京: 高等教育出版社,2000 [5] 盛骤,谢式千,潘承毅,概率论与数理统计, 北京：高等教育出版社 ,2001 [6] 盛骤 等编，概率论与数理统计（第三版），北京：高等教育出版社，2001 [7] 苏金明等编，MATLAB工具箱应用 http://sshtm.ssreader.com/html/pages/showbook.asp?ssid=11193564， 2006年9月17日 19 八、 附录 8.1 程序代码： 8.1.1 问题一的程序代码 注：程序中data1为附件1数据的矩阵表示，由于篇幅过大，在此略去 n=size(data1); [ra,ir]=sort(data1(:,2)); W=[ra,ir]; k=1; for i=1:n(1) % if (data1(i,2)<=1)&data1(i,3)>75&data1(i,3)<=150 if (data1(i,2)<=1)&data1(i,3)>150 % if (data1(i,2)<=1)&data1(i,3)>=0&data1(i,3)<=75 s(k)=i; k=k+1; end end s'; M=data1(s',:); j=0; for i=2:n(1) if data1(i,1)-data1(i-1,1)==1 j=j+1; end end n2=size(M); k=1; for i=1:n2(1) for j=1:n(1) if data1(j,1)==M(i,1) ss(k)=j; k=k+1; end end end ss; data2=(data1(ss,:)); x1=data2(:,2); x2=x1.^2; y=data2(:,3); x=[x1,x2]; 20 X=[ones(size(x1),1),x1,x2]; %stepwise(x,y); [b1,bint1,r1,rint1,s1]=regress(y,X);%regress of cd4 after sort xx=-10:80; yy=b1(1)+b1(2)*xx+b1(3)*xx.^2;yyn=max(yy); xxn=(sqrt(b1(2)^2-4*b1(3)*(b1(1)-yyn))-b1(2))/(2*b1(3)); plot(xx,yy) max(yy) %plot(data2(:,2),data2(:,3),'.'); hold on x3=data2(:,4); x4=data2(:,4).^2; y2=data2(:,5)*50; %enlarge y2 to '50*y2' so it %will be clearer when ploted with yy; xs=[x3,x4]; Xs=[ones(size(x3),1),x3,x4]; %stepwise(x,y); [b2,bint2,r2,rint2,s2]=regress(y2,Xs);%regress of hiv after sort xx=-10:80; yy=b2(1)+b2(2)*xx+b2(3)*xx.^2;yyn=min(yy); xxn2=(sqrt(b2(2)^2-4*b2(3)*(b2(1)-yyn))-b2(2))/(2*b2(3)); plot(xx,yy) min(yy) grid text(60,320,'\\leftarrow HIV'); text(60,105,'\\leftarrow CD4'); title('HIV and CD4'); xlabel('date(week)'); ylabel('HIV*50 or CD4'); hold off xt1=data2(:,3);xt2=xt1.^2;xt3=1./(xt1+1);yt=data2(:,5); xt=[xt1,xt2]; Xt=[ones(size(xt1),1),xt1,xt2,xt3]; %stepwise(xt,yt); [b3,bint3,r3,rint3,s3]=regress(yt,Xt); xx=0.0000001:300; yyt=b3(1)+b3(2)*xx+b3(3)*xx.^2+b3(4)./(xx+1); figure(2); plot(xx,yyt) hold on 21 plot(data2(:,3),data2(:,5),'.') hold off xxn úte when cd4 reach max(cd4) xxn2 úte when hiv reach min(hiv) figure(3); bx1=data1(:,2);bx2=bx1.^2;by=data1(:,3); bX=[ones(size(bx1),1),bx1,bx2]; bx=[bx1,bx2]; %stepwise(bx,by) [b4,bint4,r4,rint4,s4]=regress(by,bX);%regress of cd4 with no sort xx=-10:80; yy1=b4(1)+b4(2)*xx+b4(3)*xx.^2; bx1=data1(:,4);bx2=bx1.^2;by=data1(:,5)*50; bX=[ones(size(bx1),1),bx1,bx2]; bx=[bx1,bx2]; %stepwise(bx,by) [b5,bint5,r5,rint5,s5]=regress(by,bX);%regress of hiv with no sort yy2=b5(1)+b5(2)*xx+b5(3)*xx.^2; %plot(xx,yy1,xx,yy2,bx1,by,'.') plot(xx,yy1,xx,yy2) grid text(60,320,'\\leftarrow HIV'); text(61,103,'\\leftarrow CD4'); title('HIV and CD4'); xlabel('date(week)'); ylabel('HIV*50 or CD4'); yyn=max(yy1) xxn3=(sqrt(b4(2)^2-4*b4(3)*(b4(1)-yyn))-b4(2))/(2*b4(3)); yyn=min(yy2) xxn4=(sqrt(b5(2)^2-4*b5(3)*(b5(1)-yyn))-b5(2))/(2*b5(3)); xxn3 xxn4 8.1.2 问题二回归分析的程序代码 size(data201); [ra,ir]=sort(data201(:,2)); D=[ra,ir]; for k=1:4 j(k)=0; for i=1:5036 if D(i,1)==k 22 j(k)=j(k)+1; end end end X1=data201(D(1:j(1),2),:); X2=data201(D(j(1)+1:j(1)+j(2),2),:); X3=data201(D(j(1)+j(2)+1:j(1)+j(2)+j(3),2),:); X4=data201(D(j(1)+j(2)+j(3)+1:j(1)+j(2)+j(3)+j(4),2),:); data1=X4; n=size(data1); [ra,ir]=sort(data1(:,2)); W=[ra,ir]; k=1; for i=1:n(1) %if (data1(i,4)<=1&data1(i,5)>3&data1(i,5)<4.5 if (data1(i,4)<=1)&data1(i,5)<3 %if (data1(i,4)<=1&data1(i,5)>4.5 s(k)=i; k=k+1; end end s'; M=data1(s',:); j=0; for i=2:n(1) if data1(i,1)-data1(i-1,1)==1 j=j+1; end end n2=size(M); k=1; for i=1:n2(1) for j=1:n(1) if data1(j,1)==M(i,1) ss(k)=j; k=k+1; end end end ss; 23 data2=(data1(ss,:)); x1=data1(:,3); x2=x1.^2; x3=data1(:,4); x4=x3.^2; x5=x1.*x3; y=data1(:,5); xx=[x1,x2,x3,x4,x5]; x=[x1,x3]; rstool(x,y,'purequadratic'); 8.1.3 问题二方差分析的程序代码 n=size(data201); k=1; for i=1:n(1) if data201(i,4)>=22&data201(i,4)<=26 data20(k,:)=data201(i,:); k=k+1; end end data20; n2=size(data20); %for j=1:4 k=1; for i=1:n2(1) if data20(i,2)==1 data21(k,:)=data20(i,:); k=k+1; end end X1=data21; clear data21; k=1; for i=1:n2(1) if data20(i,2)==2 data21(k,:)=data20(i,:); k=k+1; end end X2=data21; 24 clear data21; k=1; for i=1:n2(1) if data20(i,2)==3 data21(k,:)=data20(i,:); k=k+1; end end X3=data21; clear data21; k=1; for i=1:n2(1) if data20(i,2)==4 data21(k,:)=data20(i,:); k=k+1; end end X4=data21; n3=size(data21) %s(j)= sum(data21(:,5))/n3(1) X=[X1(:,5);X2(:,5);X3(:,5);X4(:,5)];group=[ones(140,1); ones(155,1)*2; ones(140,1)*3; ones(151,1)*4]; anova1(X,group) 8.1.3 问题三的程序代码 [ra,ir]=sort(data201(:,2)); D=[ra,ir]; for k=1:4 j(k)=0; for i=1:5036 if D(i,1)==k j(k)=j(k)+1; end end end 25 X1=data201(D(1:j(1),2),:); x14=X1(:,4); t=x14; a1=0.85; a2=1.60; yw1=(t*7-floor(ceil(t*7/30)/2)*30)*a1+floor(ceil(t*7/30)/2)*30*a2; yw2=(t*7-floor(ceil(t*7/30)/2)*30)*a2+floor(ceil(t*7/30)/2)*30*a1; w1=rem(ceil(t*7/30),2); w2=1-w1; yw=yw1.*w1+yw2.*w2; x11=X1(:,1);x12=X1(:,2);x13=X1(:,3);x14=X1(:,4);x15=X1(:,5);x16=yw; XX1=[x11,x12,x13,x14,x15,x16]; X2=data201(D(j(1)+1:j(1)+j(2),2),:); x21=X2(:,1);x22=X2(:,2);x23=X2(:,3);x24=X2(:,4);x26=X2(:,4)*(1.60+1.85)*7; x25=X2(:,5); XX2=[x21,x22,x23,x24,x25,x26]; X3=data201(D(j(1)+j(2)+1:j(1)+j(2)+j(3),2),:); x31=X3(:,1);x32=X3(:,2);x33=X3(:,3);x34=X3(:,4);x36=X3(:,4)*(1.60+0.85)*7; x35=X3(:,5); XX3=[x31,x32,x33,x34,x35,x36]; X4=data201(D(j(1)+j(2)+j(3)+1:j(1)+j(2)+j(3)+j(4),2),:); x41=X4(:,1);x42=X4(:,2);x43=X4(:,3);x44=X4(:,4); x46=X4(:,4)*(1.60+0.85+1.20)*7; x45=X4(:,5); XX4=[x41,x42,x43,x44,x45,x46]; data1=XX4; n=size(data1); j=1; ss(j)=data1(1,1); for i=2:n(1) if data1(i,1)-data1(i-1,1)>=1 j=j+1; ss(j)=data1(i,1); end end M=ss'; n2=size(M); x1=data1(:,3); x2=data1(:,6); y=data1(:,5); 26 x=[x1,x2]; rstool(x,y,'purequadratic') 27