无序分类结果变量的回归
Regression on unordered categorical outcomes
在读本文之前,请先读《二分结果变量的回归》一文。
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} \]
无序分类结果变量的逻辑回归同样也是二分结果变量逻辑回归的拓展。要点有二:
在二分结果变量的二项逻辑回归中,logit反映了两个类别之间的对比。在无序分类结果变量的逻辑回归中,多项logistic回归一般默认最后一个类别k为参照组,将其他类别与最后一个类别作对比,总计进行k-1次对比,需要用k-1个logistic回归方程。对于Y的类别c的回归方程为\(Y_c = b0 + b1 \cdot X\)。参照组的b0与b1被设定为0,换言之,参照组的回归方程是已知的,不需要估计。
分类变量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)
函数,输入X
,P(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在缩小,但是anxious
与secure
的概率都在增大,而且anxious
的概率小于secure
的概率。odds在缩小意味着什么?这表明secure
的概率的增速小于anxious
的概率的增速。
另外,rowSums()
的计算表明,随着pm1
的增长,四个类别的概率之和始终为一。