无序分类结果变量的回归

Regression on unordered categorical outcomes

发布于

2025年8月19日

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

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

1 数据与案例

Keng程序包中包含数据depress,我们将depress中的变量依恋回避avo依恋焦虑anx按照得分高低进行两两组合,将被试划入四种依恋风格style:安全型,回避型,焦虑型,恐惧型。style是一个无序分类变量,它的四个类别之间没有顺序差异。我们的问题是:我们是否可以使用问题定向应对方式pm1预测被试的依恋风格style

data("depress")
depress <- depress |> 
  mutate(
    zavo = Scale(avo, sd = 1),
    zanx = Scale(anx, sd = 1),
    style = case_when(
      (zavo <= 0) & (zanx <= 0) ~ "secure",
      (zavo > 0) & (zanx <= 0) ~ "avoidant",
      (zavo <= 0) & (zanx > 0) ~ "anxious",
      (zavo > 0) & (zanx > 0) ~ "fearful"
    )
  )
head(depress$style)
## [1] "secure"  "anxious" "secure"  "anxious" "secure"  "secure"
table(depress$style)
## 
##  anxious avoidant  fearful   secure 
##       15       18       23       38
prop.table(table(depress$style))
## 
##   anxious  avoidant   fearful    secure 
## 0.1595745 0.1914894 0.2446809 0.4042553

2 无序分类结果变量的logistic回归

对于二分结果变量Y的logistic回归,logistic回归方程的结果变量是\(Y^*\),(Y = 1)的概率为 :

\[ \begin{align} P(Y = 1| X) &= F(Y^*)\\ &= F(logit)\\ &= \frac {1}{1 + e^{-1 \cdot (b0 + b1 \cdot X)}}\\ &= \frac {e^{b0 + b1 \cdot X}}{1 + e^{b0 + b1 \cdot X}} \end{align} \]

无序分类结果变量的逻辑回归同样也是二分结果变量逻辑回归的拓展。要点有二:

  1. 在二分结果变量的二项逻辑回归中,logit反映了两个类别之间的对比。在无序分类结果变量的逻辑回归中,多项logistic回归一般默认最后一个类别k为参照组,将其他类别与最后一个类别作对比,总计进行k-1次对比,需要用k-1个logistic回归方程。对于Y的类别c的回归方程为\(Y_c = b0 + b1 \cdot X\)。参照组的b0与b1被设定为0,换言之,参照组的回归方程是已知的,不需要估计。

  2. 分类变量Y有k个类别,对于某个类别c,其概率为:

\[P(Y = c|X) = F(Y^*) = F(logit) = \frac {e^{b0 + b1 \cdot X}}{\sum_i^k e^{b0 + b1 \cdot X}}\]

注意分母的变化。根据上式将Y的所有类别的P的公式写出并相加,P的和为1。

如何解释\(b0 + b1 \cdot X\)的含义?根据上面的公式,P(Y=c|X)为,

\[P(Y=c|X) = \frac {e^{b_0 + b_1 \cdot X}}{\sum_i^k e^{b_{0k} + b_{1k} \cdot X}}\] 参照组的概率P(Y=k|X)为,

\[P(Y=k|X) = \frac {e^{0 + 0 \cdot X}}{\sum_i^k e^{b_{0k} + b_{1k} \cdot X}} = \frac {1}{\sum_i^k e^{b_{0k} + b_{1k} \cdot X}}\] 那么,P(Y=c|X)与P(Y=k|X)的odds为,

\[\frac{P(Y=c|X)}{P(Y=k|X)} = e^{logit} = \frac{e^{b0 + b1 \cdot X}}{e^{0+0 \cdot X}} = e^{b0 + b1 \cdot X}\]

那么logit为,

\[log(\frac{P(Y=c|X)}{P(Y=k|X)}) = logit = b0 + b1 \cdot X\]

即,\(logit = b0 + b1 \cdot X\)反映了Y=c相对于Y=k的概率。对于b0与b1的解释可参考二项逻辑回归对b0与b1的解释。

3 案例演示

我们使用pm1预测style

fit.nom <- multinom(style ~ pm1, depress)
## # weights:  12 (6 variable)
## initial  value 130.311670 
## iter  10 value 121.539857
## final  value 121.539839 
## converged
summary(fit.nom)
## Call:
## multinom(formula = style ~ pm1, data = depress)
## 
## Coefficients:
##          (Intercept)        pm1
## avoidant    2.056519 -0.5843319
## fearful     3.405225 -0.9559211
## secure      1.821685 -0.2719101
## 
## Std. Errors:
##          (Intercept)       pm1
## avoidant    1.711787 0.5189068
## fearful     1.625192 0.4998518
## secure      1.550626 0.4594417
## 
## Residual Deviance: 243.0797 
## AIC: 255.0797
confint(fit.nom)
## , , avoidant
## 
##                 2.5 %    97.5 %
## (Intercept) -1.298523 5.4115605
## pm1         -1.601370 0.4327066
## 
## , , fearful
## 
##                 2.5 %     97.5 %
## (Intercept)  0.219908 6.59054260
## pm1         -1.935613 0.02377041
## 
## , , secure
## 
##                 2.5 %   97.5 %
## (Intercept) -1.217486 4.860857
## pm1         -1.172399 0.628579

可见,上述结果中没有anxious,即上述分析以anxious为参照组。以secure为例,secure组的回归方程为\(Y_{secure} = 1.821685 - 0.2719101*pm1\),但截距与斜率的95%置信区间包含0,即pm1的预测作用不显著。仅从方向上来看,pm1的预测作用为负,即个体的pm1的水平越高,个体越不可能为secure依恋风格。从odds上来看,-0.2719101意味着pm1每增加一个单位,odds缩小为越来的\(e^-0.2719101\) = 0.7619228。就概率而言,我们可以计算style的四个类别的概率:

首先写出三个回归方程(anxious作为参照组是已知的):

\[\begin{align} Y_{avoidant} &= 2.056519 -0.5843319 \cdot X\\ Y_{fearful} &= 3.405225 -0.9559211 \cdot X\\ Y_{secure} &= 1.821685 -0.2719101 \cdot X \end{align}\]

之后计算P公式的分母:

\[\begin{align} denominator &= e^0 + e^{2.056519 -0.5843319 \cdot X} + e^{3.405225 -0.9559211 \cdot X} + e^{1.821685 -0.2719101 \cdot X}\\ &= 1 + e^{2.056519 -0.5843319 \cdot X} + e^{3.405225 -0.9559211 \cdot X} + e^{1.821685 -0.2719101 \cdot X} \end{align}\]

之后计算每个类别的概率,下面以avoidant为例:

\[\begin{align} P_{avoidant} &= \frac{e^{2.056519 -0.5843319 \cdot X}}{e^0 + e^{2.056519 -0.5843319 \cdot X} + e^{3.405225 -0.9559211 \cdot X} + e^{1.821685 -0.2719101 \cdot X}} \\ &= \frac{e^{2.056519 -0.5843319 \cdot X}}{1 + e^{2.056519 -0.5843319 \cdot X} + e^{3.405225 -0.9559211 \cdot X} + e^{1.821685 -0.2719101 \cdot X}}\\ \end{align}\]

最后,将X的值代入即可。R语言计算概率更加方便,下面定义了一个P(X)函数,输入XP(X)将计算出X取特定值时四种依恋风格的概率。

# 定义计算概率的函数
P <- function(X) {
  # 概率公式的分子
  numerator_avoidant = exp(2.056519 - 0.5843319 * X)
  numerator_fearful = exp(3.405225 - 0.9559211 * X)
  numerator_secure = exp(1.821685 - 0.2719101 * X)
  # 概率公式的分母
  denominator = 1 + numerator_avoidant + numerator_fearful + numerator_secure
  # 四种依恋风格的概率
  P_anxious = 1 / denominator
  P_avoidant = exp(2.056519 - 0.5843319 * X) / denominator
  P_fearful = exp(3.405225 - 0.9559211 * X) / denominator
  P_secure = exp(1.821685 - 0.2719101 * X) / denominator
  # 以向量输出
  c(pm1 = X, P_anxious = P_anxious, P_avoidant = P_avoidant, P_fearful = P_fearful, P_secure = P_secure)
}
# 计算X = 1时四种依恋风格的概率
P(X = 1)
##        pm1  P_anxious P_avoidant  P_fearful   P_secure 
## 1.00000000 0.04619055 0.20133341 0.53489962 0.21757641
# 计算X = 1时四种依恋风格的累积概率
sum(P(X = 1)[2:5])
## [1] 1
# 计算(X = 1-5)时四种依恋风格的概率,并存入data.frame `dat`
dat <- data.frame(
  matrix(P(seq(1, 5, 0.1)), 
       ncol =  5,
       dimnames = list(
         NULL,
         c("pm1","anxious", "avoidant", "fearful", "secure"))
       )
)
head(dat)
##   pm1    anxious  avoidant   fearful    secure
## 1 1.0 0.04619055 0.2013334 0.5348996 0.2175764
## 2 1.1 0.04945604 0.2033316 0.5205032 0.2267092
## 3 1.2 0.05289835 0.2051401 0.5059773 0.2359843
## 4 1.3 0.05652145 0.2067495 0.4913456 0.2453835
## 5 1.4 0.06032886 0.2081512 0.4766326 0.2548873
## 6 1.5 0.06432361 0.2093375 0.4618639 0.2644750
# 计算(X = 1-5)时四种依恋风格的累积概率
rowSums(dat[2:5])
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1 1
# 作图,以pm1为横轴,以四种依恋风格的概率为纵轴
dat_l <- dat |> 
  pivot_longer(
    cols = c(2:5),
    values_to = "P",
    names_to = "style"
  )
ggplot(dat_l, aes(pm1, P, color = style)) +
  geom_path()

注意,尽管secure对应的logistic回归方程中pm1的回归系数为负,这表示随着pm1的增大,secure相对于anxious的odds在缩小,但是anxioussecure的概率都在增大,而且anxious的概率小于secure的概率。odds在缩小意味着什么?这表明secure的概率的增速小于anxious的概率的增速。

另外,rowSums()的计算表明,随着pm1的增长,四个类别的概率之和始终为一。