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

统计学进展 累积 logisitic 回归模型 曾庆 1 累积 logisitic 回归模型 一 理论 一 模型定义 假设结局变量 Y 有 J 个有序分类 其自然结局顺序表示为 Y 1 2 J 每个分类 结局 的对应发生概率是 1 2 J 则其有序分类小于等于 j 的累积发生概率表示为 P Y j 1 2 j 因此可以通过指定累积概率 P Y j 的阈值将整个结局变量 Y 的 J 个有序分 类从该指定阈值点截断 使之成为二项结局分类 另外有p个自变量记为 X x1 x2 xp 表示相应的影响因素 定性 定量均可 由此定义累积 logit Y j 函数 Jjj j xjyP xjyP xjyP xjyP jyPit ln ln 1 ln log 21 21 该累积 logit Y j 函数是两个累积概率比的对数值 它测定了结局变量 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 模型 对应的概率模型形式是 exp 1 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 Y j Xi P Y j 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 ln log 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 维列联 2 4 4 表 的 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 函数生成规律的其他因素的值 满意度 结果值用 1 4 表示 每一行顺序分别是 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 5000 15000 15000 25000 这不符合数据的实际等级 情况 显示数据结构 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 后面的向量指定了等级的由大到小的顺序 函数 as integer 则将等级变量的编码值保存在数值变量 income2 中 以上整理语句由于不止 1 条所以需要写 在复合语句标志符号 中 整理数据 od2 within od gender factor gender income ordered income levels c 25000 income2 as integer 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 f income y gender 用 ftable 则可以用复合表的形式完全显示数据表 1 统计学进展 累积 logisitic 回归模型 曾庆 6 显示频数表 with od2 ftable xtabs f gender 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 则 p 0 05 Z 值大于 2 58 则 p0 05 那么对应的自变量 没有统计学意义 上述结果中 income2 的回归系数 b 是 0 5134789 其标准误是 0 19186 wald Z检验的Z值是2 6763 相应的wald 卡方检验值为7 162582 对应p 0 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 y 1 X P1 P2 p y 2 X P1 P2 P3 p y 3 X P4 p y 4 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 y 1 X 越小 而非常满意发生的概率 p y 4 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 f y class 统计学进展 累积 logisitic 回归模型 曾庆 12 prop table with od3 xtabs f y class 模型的拟合效果除了可以用分级正确率来衡量外 还可以用模型的拟合卡方检验 计算 广义决定系数来衡量 本例残差离差 deviance 207 4189 自由度 df 74 p 0 001 说明模 型实际拟合不理想 拟合优度 广义 R 2 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 yy income2 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 yy income2 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 yy 1 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 yy income2 weights f data od2 summary