有序结果变量的回归

Regression on ordinal outcomes

发布于

2025年8月17日

在读本文之前，请先读《二分结果变量的回归》一文。

library(tidyverse)
library(Keng)
library(MASS)

1 数据与案例

Keng程序包中的depress数据中包含了我们调查的学生第一个时间点的问题定向应对方式cope1i1p(1 = very seldom, 5 = very often)与第二个时间点的抑郁depr2i1(1 = very seldom, 4 = always)。我们预测问题定向应对方式的水平越高，学生的抑郁水平越低。我们首先对这两个变量作频次分析：

data("depress")
# 频次表
table(depress$cope1i10p, depress$depr2i1)
##    
##      1  2  3  4
##   1  0  4  2  3
##   2  0  6  3  1
##   3  8  6  2  2
##   4 13 17  2  3
##   5  8  8  4  2
# 条件概率表
dat <- data.frame(prop.table(table(depress$cope1i1p, depress$depr2i1), margin = 1))
names(dat) <- c("cope1i1p", "depr2i1", "P")
dat$cope1i1p <- as.numeric(dat$cope1i1p)
dat
##    cope1i1p depr2i1         P
## 1         1       1 0.0000000
## 2         2       1 0.2758621
## 3         3       1 0.3448276
## 4         4       1 0.3333333
## 5         5       1 0.5000000
## 6         1       2 0.3333333
## 7         2       2 0.4482759
## 8         3       2 0.3793103
## 9         4       2 0.5416667
## 10        5       2 0.3333333
## 11        1       3 0.3333333
## 12        2       3 0.1034483
## 13        3       3 0.1724138
## 14        4       3 0.1250000
## 15        5       3 0.0000000
## 16        1       4 0.3333333
## 17        2       4 0.1724138
## 18        3       4 0.1034483
## 19        4       4 0.0000000
## 20        5       4 0.1666667

严格地讲，我们需要将depr2i1作为有序类别变量进行统计分析。depr2i1有4个等级，等级越高表示抑郁症状越频繁。

2 底层逻辑

要考察问题应对对抑郁的预测作用，我们可以进行logistic回归或者probit回归。与二分结果变量不同，这里的depr2i1有多个类别，这里使用的logistic回归或者probit回归是对二分结果变量logistic回归与probit回归的拓展。如何拓展呢？要点有二：

  1. 由于结果变量有多个类别，我们可以通过“合并类别”的方式将结果变量转换为二分变量，然后进行logistic回归、probit回归。我们可以将depr2i1 = 1、depr2i1 = 2、depr2i1 = 3、depr2i1 = 4四个类别合并为depr2i1 = 1与depr2i1 = (2|3|4)两组，现在结果变量Y变成了二分变量，我们可以进行logistic回归与probit回归。

  2. 对于有k个分类水平的分类变量，我们需要进行总计k - 1次“合并”与回归分析。在本例中，我们需要进行3次回归分析。

    1. depr2i1 = 1 vs depr2i1 = (2|3|4)两组

    2. depr2i1 = (1|2) vs depr2i1 = (3|4)两组

    3. depr2i1 = (1|2|3) vs depr2i1 = 4两组

  3. 对于上面的三次回归分析，我们假定它们线性回归部分自变量的斜率相等，但截距不相等。

3 Logistic回归

接下来，我们使用MASS程序包中的polr()函数进行logistic回归：

fit.polr <- polr(
  factor(depr2i1) ~ cope1i1p,
  data = depress,
  Hess = TRUE,
  method = "logistic"
)
summary(fit.polr)
## Call:
## polr(formula = factor(depr2i1) ~ cope1i1p, data = depress, Hess = TRUE, 
##     method = "logistic")
## 
## Coefficients:
##           Value Std. Error t value
## cope1i1p -0.449     0.1891  -2.375
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.1804  0.6264    -3.4809
## 2|3 -0.2163  0.5817    -0.3719
## 3|4  0.7723  0.6034     1.2798
## 
## Residual Deviance: 229.1354 
## AIC: 237.1354
confint(fit.polr)
## Waiting for profiling to be done...
##       2.5 %      97.5 % 
## -0.82594404 -0.08155889

结果显示，Y* = -0.449*cope1i1p。该结果表明，(1)cope1i1p越大，抑郁水平升级（由1升为2，或由2升为3，或由3升为4）的概率越小，B = -0.449，95% CI = [-0.826, -0.082]。(2)cope1i1p每增加一个单位，抑郁水平升高一级的odds就变为原有水平的\(e^{-0.449}\) = 0.6382661倍。

4 基于概率的解释

对于logistic回归与probit回归结果的最高水平的解释是从概率上进行解释。接下来我们就根据回归分析的结果估算各个类别的概率。从上面的结果中，我们可以得到[1 vs (2|3|4)]、[(1|2) vs (3|4)]、[(1|2|3) vs 4]的回归方程：

Y1 = -2.1804 - 0.449*cope1i1p

Y2 = -0.2163 - 0.449*cope1i1p

Y3 = 0.7723 - 0.449*cope1i1p

通过上述三个方程，我们可以基于logistic累积概率分布依次计算cope1i1p取不同值时，depr2i1>=2、depr2i1>=3、depr2i1 =4的概率：

\[P_{geq2} = P(depr2i1>=2) = \frac {1} {1 + e^{-1 \cdot (-2.1804 - 0.449 \cdot cope1i1p)}} \]

\[P_{geq3} = P(depr2i1>=3) = \frac {1} {1 + e^{-1 \cdot (-0.2163 - 0.449 \cdot cope1i1p)}} \]

\[P_{eq4} = P(depr2i1=4) = \frac {1} {1 + e^{-1 \cdot (0.7723 - 0.449 \cdot cope1i1p)}} \]

下面，我们绘制出上述三个回归方程的logistic累积概率分布曲线。在下图中，横坐标为cope1i1p，纵坐标为概率，三条线分别表示上面的三个概率，绿色矩形圈出了cope1i1p的有效取值范围1-5。

plogis(q)

在R中，我们可以直接使用plogis()函数计算概率。逻辑回归基于标准logistic分布，因此plogis()的参数location = 0，scale = 1。我们很少与location、scale这两个参数打交道，我们通常只跟q打交道。实际上，q就是logistic回归方程的结果变量\(Y^*\)，也就是logit。

dat_dist <- data.frame(cope1i1p = seq(-15, 15, 1)) |> 
  mutate(
    P_geq2 = 1/(1 + exp(-1*(-2.1804 - 0.449*cope1i1p))),
    P_geq3 = 1/(1 + exp(-1*(-0.2163 - 0.449*cope1i1p))),
    # 下面的运算结果 = `plogis(0.7723 - 0.449*cope1i1p)`
    P_eq4 = 1/(1 + exp(-1*( 0.7723 - 0.449*cope1i1p))),
  ) |> 
  pivot_longer(
    cols = c(2:4),
    values_to = "P_est",
    names_to = "P"
  )
dat_area <- data.frame(
  cope1i1p = seq(1, 5, 1),
  P_est = 1
)
ggplot(dat_dist, aes(x = cope1i1p, y = P_est)) +
  geom_path(aes(colour = P)) +
  geom_area(fill = "seagreen", alpha = 0.2, data = dat_area) +
  scale_x_continuous(limits = c(-15, 15), 
                     breaks = c(-15, -10, -5, 0, 1, 5, 10, 15),
                     minor_breaks = NULL)

实际上，由于线性回归方程中的斜率相等，上述三条logistic概率曲线的形状是相同的。线性回归方程中的截距不同，这使得三条曲线沿着横坐标轴进行水平平移。

基于上面计算的概率，我们可以通过简单的加减运算计算得到depr2i1不同类别的概率：

\[P(depr2i1=1) = 1 - P(depr2i1>=2)\]

\[P(depr2i1=2) = P(depr2i1>=2) - P(depr2i1>=3)\]

\[P(depr2i1=3) = P(depr2i1>=3) - P(depr2i1>=2)\]

在MASS中，我们有一个更简便的计算结果变量不同类别的概率的方法：predict()。下面使用predict()基于前面得到的回归方程估计depr2i1不同取值的概率P_est，并以cope1i1p为横坐标，以depr2i1不同取值的概率为纵坐标绘制出曲线，depr2i1不同取值的实际概率P用散点表示。

dat_est_wide <- cbind(data.frame(cope1i1p = seq(-10, 10, 1)),
                      predict(fit.polr, 
                              newdata = data.frame(cope1i1p = seq(-10, 10, 1)), 
                              type = "p"))
dat_est_long <- pivot_longer(dat_est_wide,
                        cols = c(2:5),
                        values_to = "P_est",
                        names_to = "depr2i1")
dat_est_long$depr2i1 <- factor(dat_est_long$depr2i1)
dat_area <- data.frame(
  cope1i1p = seq(1, 5, 1),
  P_est = 1
)
ggplot(mapping = aes(x = cope1i1p)) +
  geom_path(aes(y = P_est, colour = depr2i1), data = dat_est_long) +
  geom_point(aes(y = P, colour = depr2i1, shape = depr2i1), size = 2, data = dat) +
  geom_area(aes(y = P_est), fill = "seagreen", alpha = 0.2, data = dat_area) +
  scale_x_continuous(limits = c(-10, 10), 
                     breaks = c(-10, -5, 0, 1, 5, 10),
                     minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

可见，P_est对P的拟合还是不错的。

5 Probit回归

fit2.polr <- polr(
  factor(depr2i1) ~ cope1i1p,
  data = depress,
  Hess = TRUE,
  method = "probit"
)
summary(fit2.polr)
## Call:
## polr(formula = factor(depr2i1) ~ cope1i1p, data = depress, Hess = TRUE, 
##     method = "probit")
## 
## Coefficients:
##            Value Std. Error t value
## cope1i1p -0.2694      0.111  -2.427
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -1.3154  0.3631    -3.6231
## 2|3 -0.1140  0.3468    -0.3286
## 3|4  0.4445  0.3509     1.2668
## 
## Residual Deviance: 228.9292 
## AIC: 236.9292
confint(fit2.polr)
## Waiting for profiling to be done...
##       2.5 %      97.5 % 
## -0.48801568 -0.05273406

Probit结果的解释方法与logistic回归类似。Logistic回归Residual Deviance = 229.1354，AIC = 237.1354，Probit回归Residual Deviance = 228.9292，AIC = 236.9292，皆小于Logistic回归，该结果表明Probit回归的拟合更好。

6 参考文献

Muthen, B.O, Muthen, L.K., & Asparouhov, T. (2016). Regression and Mediation Analysis using Mplus. Los Angeles, CA: Muthen & Muthen.