有序逻辑回归的结果解读和假设检验

医学和生信笔记

发布于 2026-06-01 13:26:53

690

有序逻辑回归（ordinal logistic regression）适用于因变量为等级资料的数据。

使用孙振球《医学统计学》第4版例16-4的数据。

随机选取84例患者做临床试验，探讨性别和治疗方法对该病的影响。变量赋值为：性别（X1，男=0，女=1），治疗方法（X2，传统疗法=0，新型疗法=1），疗效（Y，无效=1，有效=2，痊愈=3）。

data16_4 <- read.csv("datasets/例16-04.csv",header = T)
head(data16_4)
##   X1 X2 Y
## 1  0  0 1
## 2  0  0 1
## 3  0  0 1
## 4  0  0 1
## 5  0  0 1
## 6  0  0 1
str(data16_4)
## 'data.frame':    84 obs. of  3 variables:
##  $ X1: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ X2: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Y : int  1 1 1 1 1 1 1 1 1 1 ...

为与流行病学上对优势比的解释保持一致（当某自变量X的OR值大于1时将其视作危险因素，小于1时将其视作保护因素），对有序逻辑回归因变量Y赋值时应将专业上最不利的等级赋予最小值，最有利的等级赋予最大值。如Y为疾病严重程度，则应按从“严重”到“轻”的顺序从低到高赋1，2，3。

# 因变量变为有序因子，注意赋值
data16_4$Y <- factor(data16_4$Y, levels = c(1,2,3),
               labels = c("无效","有效","痊愈"),
               ordered = T)

# 自变量变为无序因子
data16_4$X1 <- factor(data16_4$X1,levels = c(0,1), # 以男为参考
                      labels = c("男","女"))
data16_4$X2 <- factor(data16_4$X2,levels = c(0,1), # 以传统疗法为参考
                      labels = c("传统疗法","新型疗法"))

str(data16_4)
## 'data.frame':    84 obs. of  3 variables:
##  $ X1: Factor w/ 2 levels "男","女": 1 1 1 1 1 1 1 1 1 1 ...
##  $ X2: Factor w/ 2 levels "传统疗法","新型疗法": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Y : Ord.factor w/ 3 levels "无效"<"有效"<..: 1 1 1 1 1 1 1 1 1 1 ...

使用MASS::polr拟合有序逻辑回归：

library(MASS)

fit <- polr(Y ~ X1 + X2, data = data16_4,Hess = TRUE,method = "logistic")
summary(fit)
## Call:
## polr(formula = Y ~ X1 + X2, data = data16_4, Hess = TRUE, method = "logistic")
## 
## Coefficients:
##            Value Std. Error t value
## X1女       1.319     0.5381   2.451
## X2新型疗法 1.797     0.4718   3.809
## 
## Intercepts:
##           Value  Std. Error t value
## 无效|有效 1.8128 0.5654     3.2061 
## 有效|痊愈 2.6672 0.6065     4.3979 
## 
## Residual Deviance: 150.0294 
## AIC: 158.0294

结果解读

结果x1和x2的系数以及常数项的值与SPSS结果一致。不过这里给出的是t值，不是wald值，在样本量大时两个值基本一致。

R语言的polr函数的计算逻辑如下公式所示，它计算的是：落入低等级（坏结果）的累积概率。也就是下面这个公式中等号左侧的部分，如果回归系数β是负值，那么这个概率就会增加，如果β是正的，那么这个概率就减少。

logit P(Y ≤ k ∣ x) = ζk − β * X

本例X1女的系数是1.319，系数是正的，说明相比于参考组（男性），女性落在“疗效较差”等级（无效或有效）的对数几率降低了，优势比OR=exp(1.319)=3.73968，说明女性患者获得比男性更高疗效等级（例如从“无效”跳到“有效”，或从“有效”跳到“痊愈”）的累积优势比是男性的3.74倍。同理，新型疗法的疗效优于传统疗法。

Intercepts中，无效|有效的值是1.8128，这是参考组（男性+传统疗法）处于“无效”等级的对数几率。概率计算方法如下：

无效

说明对于使用传统疗法的男性患者，预测其疗效为“无效”的概率高达86.0%。

有效|痊愈的值是2.6672，这是参考组（男性+传统疗法）处于“无效或有效”等级的累积对数几率。计算可得：对于使用传统疗法的男性患者，预测其未能痊愈（即无效或有效）的概率约为93.5%。

有效

但是结果中没有给出P值，手动计算P值：

p <- pnorm(abs(coef(summary(fit))[, "t value"]),lower.tail = F)*2
p
##         X1女   X2新型疗法    无效|有效    有效|痊愈 
## 1.425572e-02 1.392807e-04 1.345300e-03 1.092866e-05

各个变量的系数以及OR值、可信区间等，与二分类逻辑回归计算方式一致：

# β值，也就是回归系数，两种方法都可以
coefficients(fit)
##       X1女 X2新型疗法 
##   1.318755   1.797300
#coef(fit)

# β值的95%可信区间
confint(fit)
## Waiting for profiling to be done...
##                2.5 %   97.5 %
## X1女       0.3022840 2.428170
## X2新型疗法 0.8995976 2.758261

# OR值，比值比（odds ratio）
exp(coef(fit))
##       X1女 X2新型疗法 
##   3.738765   6.033338

# OR值的95%的可信区间
exp(confint(fit))
## Waiting for profiling to be done...
##               2.5 %   97.5 %
## X1女       1.352945 11.33812
## X2新型疗法 2.458614 15.77240

假设检验

在拟合有序逻辑回归（因变量Y包括g个类别）时，需要对所拟合的g-1个方程对应的累计概率曲线的平行性进行检验，即检验各自变量在不同累计概率模型中的回归系数是否相同。当平行性假设未满足时，说明资料不适合有序logistic回归模型，应该采用多分类logistic回归模型。

Brant-Wald检验分别对每一个自变量进行Wald检验，判断该变量在不同切点（Cut-points）上的回归系数是否显著不同。

# 平行线检验(Brant-Wald test)
brant::brant(fit)
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      1.83    2   0.4
## X1女      1.59    1   0.21
## X2新型疗法       0.01    1   0.94
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

Omnibus(整体检验)是检验所有自变量是否同时满足平行线假设。然后下面是检验每一个自变量是否满足平行线假设。

P值>0.05，平行线检验通过，可以使用有序逻辑回归。

模型整体的显著性检验，可以借助似然比检验：

# 先构建一个只有截距的模型
fit0 <- polr(Y ~ 1, data = data16_4,Hess = TRUE,method = "logistic")

# 两个模型比较，似然比检验
anova(fit0, fit)
## Likelihood ratio tests of ordinal regression models
## 
## Response: Y
##     Model Resid. df Resid. Dev   Test    Df LR stat.     Pr(Chi)
## 1       1        82   169.9159                                  
## 2 X1 + X2        80   150.0294 1 vs 2     2  19.8865 4.80508e-05

P值<0.01，模型是有意义的。

计算伪R2，也就是书中的广义决定系数，可用于评价模型的拟合优度：

# 除了伪R2，还给出了更多统计量
DescTools::PseudoR2(fit, which = "all")
##        McFadden        CoxSnell      Nagelkerke   AldrichNelson VeallZimmermann 
##       0.1170373       0.2108068       0.2429443              NA              NA 
##           Efron McKelveyZavoina            Tjur             AIC             BIC 
##              NA              NA              NA     158.0294131     167.7526803 
##          logLik         logLik0              G2 
##     -75.0147065     -84.9579583      19.8865036

每一项的意义可以参考下面这张图：

对于有序逻辑回归，还可以使用LipsitzTest评价模型的拟合优度：

library(generalhoslem)
## Loading required package: reshape
lipsitz.test(fit)
## 
##  Lipsitz goodness of fit test for ordinal response models
## 
## data:  formula:  Y ~ X1 + X2
## LR statistic = 92.526, df = 9, p-value = 5.551e-16

p<0.05，说明模型拟合欠佳，结合伪R2也很小，说明这个模型拟合确实欠佳。

结果预测

获取模型预测的类别：

pred <- predict(fit, data16_4, type = "class")
head(pred)
## [1] 无效 无效 无效 无效 无效 无效
## Levels: 无效 有效 痊愈

获取模型预测的概率：

prob <- predict(fit, data16_4, type = "probs") # 或者使用 fitted(fit)
head(prob)
##        无效       有效       痊愈
## 1 0.8597003 0.07536263 0.06493706
## 2 0.8597003 0.07536263 0.06493706
## 3 0.8597003 0.07536263 0.06493706
## 4 0.8597003 0.07536263 0.06493706
## 5 0.8597003 0.07536263 0.06493706
## 6 0.8597003 0.07536263 0.06493706

oridnal

除了MASS之外，还可以用ordinal::clm做有序逻辑回归，它的输出结果中给的是z值而不是t值，但是这个模型不能直接提供给brant函数做平行检验。

# 简单演示下用法
library(ordinal)
model <- clm(Y ~ X1 + X2, data = data16_4, link = "logit")
summary(model)
## formula: Y ~ X1 + X2
## data:    data16_4
## 
##  link  threshold nobs logLik AIC    niter max.grad cond.H 
##  logit flexible  84   -75.01 158.03 5(0)  2.68e-07 4.5e+01
## 
## Coefficients:
##            Estimate Std. Error z value Pr(>|z|)    
## X1女         1.3188     0.5381   2.451 0.014256 *  
## X2新型疗法   1.7973     0.4718   3.809 0.000139 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##           Estimate Std. Error z value
## 无效|有效   1.8128     0.5654   3.206
## 有效|痊愈   2.6672     0.6065   4.398

nominal_test是一个针对比例优势假设的拟合优度检验，它通过检验“非平行模型是否显著优于平行模型”来间接判断平行线假设是否成立。

nominal_test(model)
## Tests of nominal effects
## 
## formula: Y ~ X1 + X2
##        Df  logLik    AIC     LRT Pr(>Chi)
## <none>    -75.015 158.03                 
## X1      1 -74.297 158.59 1.43620   0.2308
## X2      1 -74.936 159.87 0.15773   0.6913

两个变量的P值都是大于0.05的，说明两个变量都没有违背比例优势假设，即满足平行线假设。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-05-29，如有侵权请联系 cloudcommunity@tencent.com 删除

变量