冠心病与年龄

1、10.6冠心病与年龄问题 冠心病（Coionary Heart Disease,简称CHD）是一种常见的心脏疾病，严重地危 害着人类的健康。到目前为止，其病因尚未完全研究清楚，医学界普遍认同的、重要的易 患因素是、高血脂、高血压、糖尿病、动脉粥样硬化及家族史等。多项研究表明，冠心病 发病率随着年龄的增加而上升，在冠心病的流行病学研究中，年龄也是最常见的混杂因素 之一。为了更好的说明冠心病发病率与年龄的关系，医学家们对100名不同年龄的人进行观察, 表1中给出了这100名被观察者的年龄及他们是否患冠心病的数据（其中冠心病一栏中，1 代表被观察者患冠心病，0代表不患冠心病）。本节的目的是根据这些数

2、据建立数学模型，来分析冠心病发病率与年龄的关系，并进行统计预测闾表1100名被观察者的年龄与是否患冠心病的观察数据序号年龄冠心病序号年龄冠心病序号年龄冠心病序号年龄冠心病12002635051441765512230273505244177561324028360534507856142502936154451795615251303605546080570626031370564618157072603237157470825718280333705847083571928034380594718457110290353806048085571113003639061481865801230

3、0373916248187581133003840063490885811430039401644908959115300404106549190591163014141066500916001732042420675019260118320434206851093611193304442069520946212033045421705219562121340464307153196631223404743072531976402334148431735419864124340494407455099651253405044075551100691分析与假设 假设这100名被观察者是独立选取的

4、，记x为被观察者的年龄，Y为被 观察者患冠心病的情况（丫=1表示患冠心病，丫=0表示未患冠心病），显然丫是一个01变 量。利用表1的数据作出Y对x的散点图（见图1）。8.5从图1容易看出，直接对上述数据建立像前面几节那样的回归模型是行不通的，需要对 数据进行预处理。数据预处理的一种方法是将被观察者按年龄进行分组，并统计各年龄段中患冠心病的人 数，及患病人数占该年龄段总人数的比例（以下简称患病比例）。为方便起见，我们将年龄分 成8个年龄段，分段后的数据见表2。表2各年龄段的冠心病人数及比例年龄段年龄段中点人数患冠心病人数患病比例20 2924.51010.130 34321520.1335-39

5、371230.2540 444215503345-49471360.4650 5452850.6355-595717130.7660 6964.51080.80合计100430.43为考察患病比例与年龄的关系，首先根据表2数据做出患病比例对各年龄段中点的散点 图（见图2,为方便起见，散点的横坐标均简单地取各年龄段的中点，y20. 10： o 10 20 30,40,50,60从图2可以看出，冠心病患病比例随年龄的增大而递增，大致是一条介于0与1之间的S型曲线，这条曲线应该怎样用回归方程来确定呢？表1和表2中的患病比例实际上就是年 龄为x时（以下均取年龄段的中点）Y的平均值，用（条件）期望的符号

6、记作y = E（Y|x）患病比例y是年龄X的函数，其取值在区间0,1上。如果用普通的方法建立回归方程， 那么很容易求得其线性回归曲线或更接近与S型曲线的3次多项式回归曲线（分别见图3和 图4）,其回归模型的形式为V = Pq + /3X + 夕3X3 + （2）其中随机误差，服从均值为0的正态分布，特别地，当A =凡=0时为线性回归模型。然而在这个问题中，（2）式回归方程中y的取值不一定在0,1中，进一步说， 即使y的值在0,1中，由于在给定x时，误差项&也只能取0, 1两个值，显然5 不具有正态性，而且的方差依赖于x,具有异方差性，这些都违反了普通回归分 析的前提条件。因此，当Y为一个二分类

7、（或多分类）变量而不是连续变量时，用前几节介 绍的基于最小二乘法的回归分析是不合适的，需要用到新的回归模型。logit模型下面用万（X）表示年龄为X的被观察者患冠心病的概率，即乃（x） = P（Y = l|.r）（3）显然Y的（条件）期望为E （Y|x） =71 （x）,（条件）方差为D （Y|x）=4。）（1 一万）由 （1）式可知，兀（X）即为该年龄段的患病比例y.为了寻求患病比例加（x）与年龄x之间、形如图2的S型曲线的函数关系，并注意到万（x）在0,1区间取值，可以建立如在第5章、第7章多次用到的logistic模型 产心%（x） =hk（ 4 ）（4）的反函数写作In产乙=凡+2/（

8、5）1-4（X）（5）式左端可看作n（x）的交换，记作log（万（x） = ln,称为log”模型或logisfic回 1 - 7T（X）归模型。当n（X）在0,1取值时，1。8（小）取值为（8,+8）。在数据预处理上时，将被观察者的年龄分成k=8组，记第i组（i = l,2,女）年龄为七， 被观察人数为凡，患病人数为每位被观察者患病概率为多=?,/% ,这时log模型具 有如下形式：7t.log （乃i） = 111 - = PQ+ P,X；（6）1-%其中夕0,4是回归系数。合理地设？服从二项分布3（%,叫）,瓦），4可用最大似然法估计得 到网。模型求解log模型是一种广义线性模型（Gen

9、eralized Linear Model）,可利用MATLAB统计工具箱中的命令0疗/求解，通常的使用格式为：b= glmfit（x, y* distr / link）或b, dev, stats = glmfit（x, y, distf, link1）其中输入x为自变量数据矩阵，缺省时会自动添加一列1向量作为x的第一列；y为因变 量数据向量，而为估计系数时所用的分布，可以使历。疝出,poiss。等，缺省时为7。/；特别当取历nomial时，y可取一个2列矩阵，第一列为观察“成功”的次 数；link取历f模型见下面）等，缺省时为log输出b为回归系数的估 计值；dev为拟合偏差，是一般的残差

10、平方和的推广；stats输出一些统计指标，详见 MATLAB的帮助文件.用表2的数据输入并执行以下的程序：Age=24.5 32 37 42 47 52 57 64.5;C/?J = 1 2 3 5 6 5 13 8Total=10 15 12 15 13 8 17 10;proport = chd./Total,/?, dev, stats = ginifit (Age, ChdTotal, binomiallog it);log it Fit = gbnval(b. Age, log it),piot(Age, Pi obit Jo , Age Jogxlabel(Age ); ylabe

11、l(Pi: oportionoJCHD )b.bi = stats.se, dev得到log模型中的参数，夕。,4的最大似然估计值与它的标准差(见表3)拟合偏差为 0.524 2.图5给出了 logistic回归曲线和散点图。利用命令yhat.dylo, dyhi= glmvab, x/log if)还可以给出自变量为x时因变量y的预测值及置信度为95%的置信区间，结果如表4.表4患病比例的预测值与预测区间年龄段年龄X患病比例A 预测值乃（X）置信区间202924.50.10.07830.0282,0.19923034320.130.15740.0788,0.2898 35 39370.250

12、.24000.1461,0.368340 44420.330.34810.2458,0.466645 49470.460.47440.3625,0.589150 54520.630.60410.4720,0.7227 55 59570.760.72070.5668,0.8357606964.50.800.68550,6855,0.9366 模型评价与结果分析首先我们看到log模型(5)的右端是年龄x的线性函数，如果加入x的二次项后，是否能显著地提高模型的拟合程度呢，即考虑模型log if(乃(x) = 1+ Px +1-乃(x)_为了检验是否需要引入/项，在log isc回归中要用似然比检验统

13、计量网.执行程序历2,dev? = glnifit(geAge. lChd Total, binomial/)；b2, pval = 1 - chilcdf dev - dev2,l)计算得到夕0,4,凤的估计值分别为-5. 350 6, 0. 119 4和-0. 000 2, p值为0. 937 1, 表示模型中引入X2项并不能显著提高拟合程度。其次，处理这类问题的另一种广义线性模型是模型，其形式为万(X)=(几+*丫)(8)probit (兀(x)= t(4(x)=凡 + 0X(9)其中。是正态概率分布函数，它也是s型曲线，利用MATLAB统计工具箱中的命令g/矿 求解时，只需将log i

14、t改为probit即可。用表2数据计算得到模型系数的估计值分别为-2. 993 3和0. 624 4,标准 差分别为0.601 1和0. 12 8,拟合偏差为0.652 9.与log模型的结果比较见表5,拟合 曲线比较见图6,可以发现这两个模型的拟合程度不相上下。表5 欣/模型与log模型的比较年龄段年龄Z患病 比例预测值(logit)预测值(probit)置信区间(logit)202924. 50. 10. 07830. 07150, 028 2 ,0. 199 23034320. 130. 15740. 15950,078 8 ,0. 289 83539370. 250. 24000. 2

15、4690. 1461,0. 368 34044420. 330. 34810. 35480. 2458, 0. 466 64549470. 460. 47440. 17590. 3625, 0. 589 15054520. 630. 60410. 59940. 4720, 0. 751 055 59570. 760. 72070. 71350. 5668, 0. 831460-6961. 50. 800. 85010. 84890. 6854, 0. 943 0模型预测与进一步分析 通过上诉分析可知，对于我们的问题和观察数据，模型。Alog (X)= 一5.03 82+0.1050k(10)1

16、 7F(x)是一个合适的模型，从(10)式能够给出任何年龄的人患冠心病的概率及相应的置信区 间，例如，图7给出了年龄分别为20, 30, 40, 50, 60, 70, 80的人患冠心病的概率， 以及置信度为95%的置信区间.10g模型的另一个好处是其中的回归系数有很好的直观解释。10g模型与统计 中0dds (发生比或优势)的概念有密切的联系，所谓的odds就是事件的发生概率与不发生 概率之比。本节中，若记。dds(X)为年龄x的人患与不患冠心病概率之比，则Odds (x)二双(11)1-乃(工)于是10g模型可以表示为odds(x) = *+觥(12)当年龄增加一岁时。dds比(发生比率或

17、优势比)为” + 1)一： (13) odds(x) +侬于是C)(14)odds(x)即4为自变量增加1个单位时。dds比的对数。丹0时，e百1, x每增加1个单 位，。龙/s比会相应增加，且对任意整数k,有oddsx + k) = ekp- odds(x)(15)在模型(10)中注=0.1050，可以算出一个20岁的青年人患冠心病的概率仅为 A乃(20) = 0.0503 ,且发生比(患与不患冠心病的概率之比)为。c/c/s(20) = 0.0593,说明这个年龄的人患冠心病几乎是不大可能的。10年后，30岁人的发生比就变成 odds(30) = e10 xA x 0.0593 = 0.1

18、695 ,发生比(可解释作危险率)增大到20岁时的2. 8577 倍，而到 60 岁时，odds(60) = e40 xA x 0.0593 = 3.9545 ,危险率是 20 岁的 e40A = 66.6863倍。可见回归系数4在log”模型中有着重要的意义。这一点在peobit 模型终是无法体现的。最后，在log”模型中，人们常常感兴趣的是，x取何值时)(# = 0.5。由模型(10) 求解(1)=0-5,可得=47.98。这就是说，当你到48岁时，患冠心病的概率就会 大于不患冠心病的概率，要格外小心了！评注因变量是定性变量的回归分析作为一种有效的数据处理方法已被广泛应用， 尤其在医学、社会调查、生物信息处理等领域，这类回归模型属于广义线性模型的研究 范畴，本节只涉及因变量是0-1变量且只有一个自变量的情形，对多个变量反 匕”的 情形，log模型和模型的一般形式分别为(16)(16)p