R软件-医学统计分析-有序logistic回归

统计学进展/累积 logisitic 回归模型/曾庆 1 累积 logisitic 回归模型 一、理论 （一）模型定义 假设结局变量 Y 有 J 个有序分类， 其自然结局顺序表示为 Y=1,2,., J， 每个分类 （结局） 的对应发生概率是 1, 2,., J，则其有序分类小于等于 j 的累积发生概率表示为 P(Yj)= 1+2+.+ j，因此可以通过指定累积概率 P(Yj)的阈值将整个结局变量 Y 的 J 个有序分 类从该指定阈值点截断， 使之成为二项结局分类。 另外有p个自变量记为 X=(x1, x2， ， xp)， 表示相应的影响因素（定性、定量均可） 。由此定义累积 logit(Yj) 函数： Jjj j xjyP xjyP xjyP xjyP jyPit + + = = = + . . ln )|( )|( ln )|(1 )|( ln)(log 21 21 该累积 logit(Yj)函数是两个累积概率比的对数值，它测定了结局变量 Y 小于等于分类 j 或者大于分类 j 的可能性大小。因为有 J 个有序分类结局，所以实际上可以写出至多 J-1 个的累积 logit 函数。 + = = = + = = = + J J Jj xJyP xJyP xJyP xJyP JyPit xyP xyP xyP xyP yPit J 121 2 1 . ln )| 1( )| 1( ln )| 1(1 )| 1( ln)1(log . ln )| 1( )| 1( ln )| 1(1 )| 1( ln)1(log 1L个 将每个累积 logit 函数用线性函数的形式表示为 += += ppjjj pp XXJyPit XXyPit J , 111 , 10 , 1 111110 .)1(log .)1(log 1 L个 该模型就是累积 logistic 回归模型(cumulative logit model)。对应的概率模型形式是 ).exp(1 ).exp( )( 110 110 pjpjj pjpjj XX XX jyP + + = 为简化上述模型，进一步假定对于所有 J-1 个累积 logit 函数，各个自变量 Xi所对应的 回归系数 i都是相等的。即对每一个累积 logit 函数各有一个不同的截距 j0，然而对所有的 累积 logit 函数，自变量 Xi却有一个相同的 i。在此假设条件下，不同累积 logit 的回归线相 互平行，只是截距 j0不同，这称为成比例发生比假设或平行线假设。满足平行线假设的模 型简化后的是 = +=+= p i iijppj XXXjyPit 1 0110 .)(log 写成系列模型则是 += += ppJ pp XXJyPit XXyPit .)1(log .)1(log 110 , 1 1110 L 该简化后的模型就称为成比例比数比累积 logit 模型(proportional-odds cumulative logit 统计学进展/累积 logisitic 回归模型/曾庆 2 model)。 该模型有 J 1 截距， p 个回归系数， 共 J + p 1 个回归系数。 一般认为累积 logistic 回归模型就是成比例比数比累积 logit 模型。对应的概率模型形式是 ).(exp1 1 ).exp(1 ).exp( )( 110110 110 ppjppj ppj XXXX XX jyP + = + + = 使用成比例比数比累积 logit 模型首先需要对平行线假设进行检验。如果平行假设被拒 绝，便说明自变量 Xi对不同的 logit 有不同的 i，因而说明成比例比数比累积 logit 模型不适 合，需要采用其他模型来进行资料的分析。 （二）似然函数和对数似然函数 对于有 n 个独立观察对象的样本， 第i个观察对象 Xi出现获得 Y=j 分类结局的概率记作 Pj=P(Y=j|Xi)，它是累积概率函数的差，即 Pj=P(Y=j|Xi)= P(Yj|Xi)P(Yj-1|Xi) 。由此构 造似然函数 L 为： = = n i J j ii n i J j y ij n i y iJ y i y i ijij iJii XjYPXjYPPPPPL 11 y 111 21 )| 1-(-)|()()( 21 L 式中 yij表示第i个观察对象自变量（探索变量）取值 Xi时已观察到结局变量 Y=j 等级 时所对应的编码，它满足 = = J j ij y 1 1，即只有一个 yij取值为 1（属于该类） ，其余均为 0。 相应的对数似然函数 LL 为 = = = n i J j jiijij n i J j ijij PitPityPyLL 11 1, 11 )(log)(lnlog)ln( 二、实例分析 【实例实例 1】 现有资料如表 1 所示， 试使用累积 Logistic 回归模型考察工作满意度与性别、 收入水平之间的关系。 表表 1 不同性别及收入水平人群不同性别及收入水平人群的的满意度情况满意度情况 满意度（人）（y） 性别 （gender） 收入水平 （income） 非常不满 稍微满意 比较满意 非常满意 女性 25000 0 2 4 2 男性 25000 0 1 9 6 在【实例 1】中，结果变量工作满意度 y 是多值有序变量，两个原因变量性别 gender 定性变量和收入水平 income 是有序变量。传统的分析方法包含分层的列联表 2 检验、高 统计学进展/累积 logisitic 回归模型/曾庆 3 维列联（244 表）的 CMH 2 检验，模型分析可以采用对数线性模型、累积 Logistic 回 归模型、多项 logistic 回归模型进行分析。累积 Logistic 回归模型进行分析同时考察性别与 收入水平对工作满意度的影响。下面采用累积 Logistic 回归模型进行分析。 在进行分析计算时， 数据的输入形式可以采用原始数据和频数数据形式。 原始数据形式 采用一行一例，各个指标（因素）对应一个变量的格式输入，而频数表形式在在上述格式基 础上再增加一个频数变量。 本实例数据表是频数表形式， 所以输入时课采用频数表格式输入。 其中的满意度 y 的结 果值是非常不满，稍微满意，比较满意，非常满意 4 个等级，而表内的数据值 1、3、11、2 等是各个等级对应的频数，因为表 1 是复合表的形式，它的 2 个因素变量(gender、income)、 1 个结果变量(y)的结果值具有规律性，适合程序生成，所以决定首先输入频数 f，然后再采 用程序分别完成其他各个因素的输入。 #输入所有的频数，输入完数据后，必须多一个空行 f-scan() 1 3 11 2 2 3 17 3 0 1 8 5 0 2 4 2 1 1 2 1 0 3 5 1 0 0 7 3 0 1 9 6 #数据因素变量 y-rep(1:4,8) gender-rep(c(F,M),c(4*4,4*4) income - rep( rep(c(25000), each=4), 2) #将数据合并为数据框(集)的形式 od - data.frame(y,gender,income,f, stringsAsFactors=F) 上述程序中采用 scan()函数进行频数输入，生成向量 f，需要注意的是，scan()函数的数 据行最后以空行表示数据输入完成。然后采用 rep()函数生成规律的其他因素的值。满意度 结果值用 14 表示，每一行顺序分别是 1、2、3、4 取值，总共 8 行；gender 语句中 4*4 表 示前面的数据结果值 F 、 M 每个都是 4 列 4 行共 16 个数据值； 而 income 语句中 each=4 表示每个数据值 4 个，外面的 rep 函数表示数据值重复了 2 次。 最后采用 data.frame()函数将所有向量组织成数据框（集）形式保存，形成数据整体。R 在转换为数据框的形式时，会自动将字符向量转换为因子变量，选项 stringsAsFactors=F 则 取消这个自动转换功能，因为自动转换功能对 gender 和 income 都会进行转换，而且 income 的转换顺序是25000、500015000、1500025000，这不符合数据的实际等级 情况。 #显示数据结构 str(od) 数据框 od 数据结构共 32 例， 4 个变量， 其中变量 y 和 income 是数值型 （整型和实型） ， gender 和 income 则是字符型。下面是数据框（集）的前 12 个数据，实际输入时也可以将数 据整理成下面的形式然后用 read.table()或者 scan()函数来输入。 #显示数据头 12 个 统计学进展/累积 logisitic 回归模型/曾庆 4 下面对输入数据进行整理，整理的结果保存在数据框 od2 中。函数 within()则指定整理 的数据框为 od。用 factor()函数将 gender 因子化；用 ordered()函数将 y、income 按正确的等 级形式生成等级变量， 其中 levels 后面的向量指定了等级的由大到小的顺序； 函数 eger() 则将等级变量的编码值保存在数值变量 income2 中。 以上整理语句由于不止 1 条所以需要写 在复合语句标志符号 中。 #整理数据 od2 - within(od, gender-factor(gender) income-ordered(income,levels=c(25000) ) eger(income) yy-ordered(y) ) #整理后的数据结构 整理后的数据框 od2 的结构显示共 32 例，6 个变量，除了原来的变量外新增 2 个变量 yy 和 income2，其中 yy 是有序等级变量，income2 是整型数值变量。显示的前 12 个数据值 可见，经过数据转化后，变量 gender、income 的屏幕可视形式没有变化，income2 中的数值 是 income 变量的对应等级编码。 #整理后的数据框（集）的头 12 个数据 统计学进展/累积 logisitic 回归模型/曾庆 5 上述数据框中，有频数为 0 的数据例，这对分析没有意义，必须删除，否则影响后续的 分析。可以采用子集函数 subset()，也可以采用数据框的下标运算来进行删除数据行。操作 完毕后的数据结果显示，原来的数据行中频数为 0 的被删除了。例如原来的 9 行是 0，现在 是删除后没有编号 9 的数据行了。 #删除频数为 0 的数据 od20) #显示头 12 个数据 head(od2,12) 数据整理好后，用 xtabs()函数进行列联表计算，可以将原来的数据表正确的显示出来。 xtabs(fincome+y+gender)。用 ftable()则可以用复合表的形式完全显示数据表 1。 统计学进展/累积 logisitic 回归模型/曾庆 6 #显示频数表 with(od2, ftable(xtabs(fgender+income+y) ) R 有多个软件包都能进行累计 logistics 回归模型的计算，这里介绍两个包，分别是 VGAM 和 MASS。其中 VGAM 包能进行所有类型的 logistic 回归模型的计算，而 MASS 仅 能计算累计 logisitic 回归模型。在使用相应的软件包在计算前必须首先用 library()或者 require()函数加载相应的软件包，其中 MASS 包是默认的安装包，而其他的包还需要用 install.package()先进行安装。 VGAM 除了能进行累计 logistic 回归模型的计算外还能进行累计 logistic 回归模型的平 行性假设的检验，而其他的软件包则不能；MASS 包虽然不能进行平行性检验，但是它可以 用 step()函数进行逐步的回归分析。 VGAM包使用vglm()函数来拟合广义线性模型， 通过family选项指定需要拟合的logistic 回归模型的类型。下面的模型语句中 y income2 指定 y 为结果变量，income2 为探索变量 （自变量） ， weights=f 指定频数变量 f （如果为原始数据则频数为 1， 则不需要使用） ， data=od2 指定所用数据框为 od2。模型中 family=cumulative(parallel=T) 指定拟合累计 logistic 回归模 型，而且 parallel=T 指定模型按平行性假定进行拟合。该选项可以简写为 family=propodds， 但是这样写的时候， 结果变量 y 的函数形式刚好与 family=cumulative(parallel=T)的形式相反。 所有计算结果保存在模型对象 om 中。 #加载 VGAM 包 library(VGAM) 统计学进展/累积 logisitic 回归模型/曾庆 7 #建立模型 om-vglm(y income2, weights=f, data=od2, family=cumulative(parallel=T) ) 为了检验模型 om 满足累积 logistic 回归模型的平行性假设，必须建立另外一个非平行 模型的假设，其模型的语句与上述形式基本相同，除了选项 family=cumulative(parallel=T)改 变为 family=cumulative(parallel=F)，这里 parallel=F 指定模型不满足平行性假设，这是实际 上是拟合多项 logistic 回归模型。 所有计算结果保存在模型对象 om2 中。 然后使用函数 lrtest() 进行似然比检验，检验两个模型（即 om 和 om2）是否有统计学差异。该检验的 H0：平行 性模型与非平行性模型的拟合效果一致（即满足平行性假设） ，H1：平行性模型与非平行性 模型拟合效果不一致。 #检验平行性假设的非平行模型 om2-vglm(y income2, weights=f, data=od2, family=cumulative(parallel=F) ) #似然比检验检验平行性假设 lrtest(om,om2) 上述似然比检验中平行模型 om 的对数似然函数值 LL 为-103.71， 而非平行模型 om2 的 对数似然函数值为-103.02， 两模型的对数似然值差为-0.6888884， 卡方值为 1.3778(即 2 倍对 数似然值的差)，p=0.5021，结果表明满足平行性假设。 整体模型的似然比检验的零假设是 H0：所有的回归系数值i 为 0，采用 lrtest()函数来 进行检验。下面的似然比检验结果表明，卡方值为 7.3627，p=0.0067，拒绝 H0，说明并非 所有的回归系数值i 都为 0 ，整体模型成立。 VGAM:lrtest(om) #显示结果 summary(om) 统计学进展/累积 logisitic 回归模型/曾庆 8 用 summary()函数显示结果。计算结果表明，经过 4 次迭代后模型就获得了结果。模型 的对数似然函数值 LL 为-103.7095，自由度 df 为 74，离差 deviance 为 207.4189，自由度 df 为 74。同时结果中指定建立的线性预测函数共 3 个，累积 logit 函数形式是从低等级 1 向高 等级 4 累积概率。 模型所得的回归系数值 b 列在 Coefficients：段。各个的回归系数的检验采用 wald Z 检 验，零假设是 H0：i=0。各个自变量的回归系数的 wald Z 检验值列在最后一列，如将 wald Z 值平方则是其他软件包（如 SAS、SPSS）中的 wald 卡方统计量。wald Z 值大于 1.96，则 p0.05，Z 值大于 2.58，则 p0.05，那么对应的自变量 没有统计学意义。上述结果中，income2 的回归系数 b 是-0.5134789，其标准误是 0.19186， wald Z检验的Z值是2.6763， 相应的wald 卡方检验值为7.162582， 对应p0.01， 说明income2 的回归系数不等于 0，有统计学意义，即收入对结果变量满意度有影响。由于只有一个自变 量 income2，所以该结果和模型总的似然比检验结果( 2 =7.3627，p=0.0067)一致。 所得累计 logistic 回顾模型的 logit 线性形式如下： income2 0.5134789- 2.6132914 )|3(log income2 0.5134789- 0.3762620- )|2(log income2 0.5134789- -2.0688262)| 1(log = = = Xypit Xypit Xypit 记满意度的 4 个等级非常不满、稍微满意、比较满意、非常满意的发生概率分别是 P1=p(y=1|X)、P2=p(y=2|X)、P3=p(y=3|X)、P4=p(y=4|X)，则 P1=p(y1|X)，P1+ P2=p(y 2|X)，P1+ P2+ P3=p(y3|X)，P4= p(y4|X)=1- P1- P2- P3。对应的概率模型形式如下： 统计学进展/累积 logisitic 回归模型/曾庆 9 ) income2) 0.5134789 914exp(2.61321/(1 )|4( income2) 0.5134789 2.6132914exp(1/(1 )|3( ) income2) 0.5134789 0.3762620exp(1/(1 )|3( income2) 0.5134789 620exp(0.37621/(1 )|2( ) income2) 0.5134789 2.0688262exp(/(11)|2( ) income2) 0.5134789 2.0688262exp(/(11 ) income2) 0.5134789- -2.0688262exp(1income2)/( 0.5134789- -2.0688262exp()| 1( += += += += += += += Xyp Xyp Xyp Xyp Xyp Xyp 或者 或者 或者 由概率模型形式可以知道，回归系数i 的估计值小于 0，说明自变量 X 的取值越大结 果变量低等级的发生概率越小，反之回归系数i 的估计值大于 0，自变量 X 的取值越大结 果变量高等级的发生概率也越大。本例中 income2 的回归系数值-0.5134789，小于 0 ，表示 收入越多，结果变量满意度的非常不满的发生概率 p(y1|X)越小，而非常满意发生的概率 p(y4|X)也越大。 可以进一步用 coef()函数提取回归系数计算各个自变量的 OR 值及其 95%可信区间。下 面的计算程序使用了 wald 可信区间近似方法来计算 OR 的可信区间。 #OR 计算, wald 近似 95%可信区间 beta-coef(summary(om) beta=|Z|=( 1-pnorm(abs(beta,3)*2, OR=exp(beta,1), OR95%CIL =exp(beta,1 - 1.96*beta,2), OR95%CIU=exp(beta,1 + 1.96*beta,2) beta 相似的，可以建立包含两个自变量 income2、gender 的累积 logistic 回归模型。其结果 如下： om6-vglm(y income2+gender, data=od2, weights=f, family=cumulative(parallel=T) ) summary(om6) 统计学进展/累积 logisitic 回归模型/曾庆 10 从结果中可见 gender 变量的 Z 值小于 1.96，表明该变量的回归系数没有统计学意义， 应该从模型中删除。该模型和包含单个变量的 income2 模型可以采用似然比检验进行比较， 用检验结果来选择更优的模型。 VGAM:lrtest(om,om6) 比较结果表明，两个模型没有显著性差别( 2 =0.017，p=0.9675)，应该采用模型参数更 少的模型，也就是只有单个自变量 income2 的模型 om。 此外选择更优模型还可以比较两个模型的信息统计量 AIC 或者 BIC，分别使用函数 AIC()、BIC()，信息统计量小的模型更优。本例中 om 模型的 AIC 是 215.4189，而 om6 模型 的 AIC 是 217.4173，也说明了应该选择单变量模型 om。 统计学进展/累积 logisitic 回归模型/曾庆 11 模型计算完毕后，使用 predict()函数可以计算连接函数的值，预测在给定 X 的情况下结 果变量的发生概率 p， 按照发生概率 p 的大小可以进一步确定在给定 X 的情况下， 结果变量 的分级结果。VGAM 包默认输出连接函数的值即 type=link，各个等级的发生概率用 type=response来计算，最可能的分级用 max.col()函数来确定。 #预测概率值 pp-predict(om,type=response) od3-cbind(od2,pp) head(od3,12) 上面的计算中，计算了给定 X 时，结果变量的可能发生概率，保存在 pp 中，最后 pp 与原数据框 od2 合并为新的数据框 od3。数据框 od3 的结果值中的1 、 2 、 3 、 4列 表示满意度各个等级的发生概率 p。 发生概率的最大等级就是估计 （预测） 的满意度的分级。 #预测分类 od3$class-max.col(pp) head(od3,12) 由估计的满意度分级和实际的满意度分级可以列出双向有序表， 求出估计的正确频数和 正确率。 由下面的计算结果可以看到， 本例的所有估计的满意度分级都是 3 级， 即比较满意， 这样得出该模型的估计正确的频数仅为 63 例，正确率仅有 60.58%，远没有达到好的水平。 这从一个方面说明了模型的拟合效果实际是不好的。 #预测分类正确率 with(od3,xtabs(fy+class) ) 统计学进展/累积 logisitic 回归模型/曾庆 12 prop.table( with(od3,xtabs(fy+class) ) ) 模型的拟合效果除了可以用分级正确率来衡量外， 还可以用模型的拟合卡方检验， 计算 广义决定系数来衡量。本例残差离差 deviance=207.4189 ，自由度 df=74，p0.001，说明模 型实际拟合不理想。 #拟合优度，广义 R2 om.null-vglm(y 1, data=od2, weights=f, family=cumulative(parallel=T) ) LLf - logLik(om) LL0 - logLik(om.null) N-sum(od$f) #McFadden pseudo-R2, Cox & Snell -R2, Nagelkerke -R2 cbind( McFadden R2= 1 - (LLf / LL0) , Cox & Snell R2=1 - exp(2/N) * (LL0 - LLf) , Nagelkerke R2=(1 - exp(2/N) * (LL0 - LLf) / (1 - exp(LL0)(2/N) ) 上面的程序计算广义 R2。广义 R2 的结果都远远的远离 1，接近 0，说明模型的实际拟 合效果并不好。 当变量数目比较多时，可以考虑采用自动逐步变量筛选。在 R 中自动逐步变量筛选用 step()函数。但是 VGAM 包的 vglm()函数却不能使用 step()函数来进行自变量筛选。要进行 自变量的筛选必须使用 MASS 包的 polr()函数才能进行。 下面是用 MASS 包 polr()来进行 om 模型的计算。 #MASS 包计算，首先加载 MASS 包 library(MASS) om3-polr(yyincome2,weights=f, method=logistic, data=od2) summary(om3) cbind(beta=coef(om3), t( confint(om3) ) ) exp( cbind(OR=coef(om3), t( confint(om3) ) ) ) 统计学进展/累积 logisitic 回归模型/曾庆 13 上述计算结果中，常数项（截距）的计算结果值与 vglm()函数的结果完全一致，但是自 变量 income2 的回归系数的值 0.5135 则刚好与 vglm()函数的 income2 的结果值-0.5135 符号 相反，这是因为两个函数定义的模型形式有差异所致。MASS 包的 polr()函数的累积 logistic 模型定义是etaXjypit j = 0 )|(log，采用的是减法形式，而一般的累积 logistic 回 归，如 vglm()，的函数定义形式是etaXjypit j += 0 )|(log，是加法形式，所以两 种定义计算出来的符号刚好相反。 MASS 包 polr()函数比较方便的是建立的模型可以直接使用 confint()函数来计算其参数 的可信区间，进而计算出 OR 的可信区间。 下面建立两变量的累积 logistic 回归模型。然后使用 drop1()函数进行自编量的似然比检 验。似然比检验显示自变量 income2 有统计学意义，而 gender 变量则没有统计学意义，可 以考虑删除。 om4-polr(yyincome2+gender,weights=f, method=logistic, data=od2) summary(om4) drop1(om4,test=Chisq) 统计学进展/累积 logisitic 回归模型/曾庆 14 在 R 中使用 step()函数进行逐步回归分析，采用的是信息量 AIC 选择自变量。R 的逐步 法函数 step()通过使用选项 direction=backward、 direction=forward、direction=both分别 指定了向后逐步法、向前逐步法和双向逐步法。默认采用向后法。如果是采用向后法必须首 先建立包含所有自变量的模型即全模型； 如果采用向前法、 逐步法除了必须要建立包含所有 变量的全模型作为逐步法的终点外还需建立最小变量模型作为逐步法起点， 一般最小模型是 不包含自变量的空模型。逐步法的结果输出中引入自变量用符号+ ，删除自变量用符号 - ，同时列出对应的 AIC 变化。 上面的分析已经建立了包含 2 个自变量的模型 om4，下面语句采用向后法筛选对模型 om4 进行自变量筛选。向后逐步计算结果保存在对象 ss1 中。 ss1-step(om4) #逐步法-向后法，ss1-step(om4,direction=backward) 统计学进展/累积 logisitic 回归模型/曾庆 15 逐步分析的对象 ss1 中包含了逐步法的简要过程，其分量名字是 anova，可以以摘要结 果的显示出来。 ss1$anova 最后，逐步法最后建立的模型可以用 summary()函数显示。 下面是指定双向逐步法的语句。 模型 om4.0 是不包含自变量的空模型， 作为逐步法的选 择模型的起点，而全模型 om4 则作为选择模型的终点。最后显示选择摘要和最终的模型结 果。 om4.0-polr(yy1,weights=f, method=logistic, data=od2) ss2-step(om4.0,formula(om4), direction=both) ss2$anova summary(ss2) 上面结果中，Start: AIC=220.78 开始时模型的 AIC，此时模型形式是 yy 1，模型外 的自变量是income2和gender， 由于引入income2后模型的AIC减少， 所以首先引入income2； 下一步 AIC=215.42，删除模型内的 income2 和增加模型外的 gender 都使 AIC 的值增大，所 以，逐步过程到此结束。 统计学进展/累积 logisitic 回归模型/曾庆 16 #逐步回归分析必须使用 MASS 包建立模型 library(MASS) om3-polr(yyincome2,weights=f, data=od2) summary(om3)