二分结果变量的回归

Regression on binary outcomes

发布于

2025年8月11日

1 概率与odds

一个事件 (例如，感冒) 有两种结果：发生或者不发生。我们用变量Y表示该事件，Y的取值有两个：1或0。Y是一个二分变量。发生的概率表示为P，那么不发生的概率则为1-P。下表呈现了一项对法国滑雪运动员的随机对照试验，该试验旨在考察维他命C与感冒的关系，安慰剂组不服用维他命C，维他命C组服用维他命C，经过一段时间后调查其是否患上感冒 (Muthen et al., 2016)。表中的第1-3列为频次。第4列为感冒的概率，其中0.221与0.122为条件概率，条件概率指的是在某种条件下的概率。0.221表示在安慰剂这种条件下的感冒的概率，0.122表示在维他命C这种条件下的感冒的概率。

表 1: 维他命C与感冒
n 未感冒 感冒 概率 Odds
安慰剂组 140 109 31 0.221 0.284
维他命C组 139 122 17 0.122 0.139
总计 279 231 48 0.208 0.263

感冒的概率大，还是不感冒的概率大？统计学家常通过比值来比较概率的大小，这个比值被称为odds：

odds = P/(1-P)

Odds

odds ratio是概率的比值。

显然，当P = 1-P时，odds = 1；当P > 1-P时，odds > 1；当P < 1-P时，odds < 1。

可见，就总体而言(N = 279)，感冒的odds为0.263，该值小于1，表明总体而言感冒的概率小于未感冒的概率。

2 Odds Ratio (OR)

P与odds有怎样的关系？我们可以计算并作图观察：

P <- seq(0, 1, 0.01)
odds <- P/(1 - P)
plot(odds, P)
图 1

可见，odds与P的关系具有单调性，odds越大，P越大，尽管它们之间的关系不是线性的。 因此，我们可以用odds比较不同条件下事件发生的概率。同样地，我们通过作比来比较odds的大小，这个比值被称为odds ratio(常被译为优势比)。

Odds Ratio (OR)

Oodds Ratio (OR) 是概率的比值的比值。

在前文的案例中，与未服用维他命C的人相比，服用维他命C的人是否不易感冒？安慰剂组(未服用维他命C)的odds为0.284，维他命C组的odds为0.139，那么，

odds ratio = odds维他命C/odds安慰剂 = 0.139/0.284 = 0.489

可得，与安慰剂组相比，维他命C组感冒的概率更小。

3 结果为二分变量的逻辑回归

当结果是二分变量时，我们如何用回归表示自变量与结果变量的的关系呢？让我们看一下下面这个例子。下面的数据描述了英国煤矿工人呼吸急促(breathlessness)这一症状与年龄的关系，第2列是各年龄的人数，第3列是有这一症状的人数，P是相应的比例。

表 2
age n breathlessness f P odds logit
22 1952 1 16 0.0081967 0.0082645 -4.7957905
27 1791 1 32 0.0178671 0.0181922 -4.0067648
32 2113 1 73 0.0345480 0.0357843 -3.3302456
37 2783 1 169 0.0607258 0.0646519 -2.7387382
42 2274 1 223 0.0980651 0.1087275 -2.2189110
47 2393 1 357 0.1491851 0.1753438 -1.7410066
52 2090 1 521 0.2492823 0.3320586 -1.1024437
57 1750 1 558 0.3188571 0.4681208 -0.7590289
62 1136 1 478 0.4207746 0.7264438 -0.3195942
22 1952 0 1936 0.9918033 121.0000000 4.7957905
27 1791 0 1759 0.9821329 54.9687500 4.0067648
32 2113 0 2040 0.9654520 27.9452055 3.3302456
37 2783 0 2614 0.9392742 15.4674556 2.7387382
42 2274 0 2051 0.9019349 9.1973094 2.2189110
47 2393 0 2036 0.8508149 5.7030812 1.7410066
52 2090 0 1569 0.7507177 3.0115163 1.1024437
57 1750 0 1192 0.6811429 2.1362007 0.7590289
62 1136 0 658 0.5792254 1.3765690 0.3195942

首先，我们尝试以breathlessness为结果变量作线性回归，breathlessness是一个二分变量。这里需要补充一点背景知识：二分变量服从伯努利分布，二分变量Y的 (平均值) 期望为事件发生的概率P，二分变量Y的方差为P(1 - P)。当我们以breathlessness为结果变量作线性回归时，实际上是用年龄age预测breathlessness，该方程估计得到的是在相应年龄上的breathlessness的平均值，也就是breathlessness = 1 (有呼吸急促症状) 的概率。接下来我们看结果：

# 以breathlessness为结果变量的线性回归
lm.fit <- lm(breathlessness ~ age, weights = f, dat)
summary(lm.fit)
## 
## Call:
## lm(formula = breathlessness ~ age, data = dat, weights = f)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.009  -6.059   3.269  12.543  17.326 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.268336   0.287116  -0.935    0.364
## age          0.009794   0.006735   1.454    0.165
## 
## Residual standard error: 10.78 on 16 degrees of freedom
## Multiple R-squared:  0.1167, Adjusted R-squared:  0.06153 
## F-statistic: 2.115 on 1 and 16 DF,  p-value: 0.1652

可见breathlessness = -0.268336 + 0.009794*age。我们可以根据该回归方程，使用age预测breathlessness，估计得到相应年龄breathlessness的患病率。在下图中，我们将实际数据 (散点图) 与回归方程 (直线) 进行对比，可以发现实际上breathlessness与age的关系不是一条直线，而是一条曲线，线性回归方程与实际数据的偏差较大。

ggplot(dat[dat$breathlessness == 1,], aes(age, P)) +
  geom_point() +
  geom_abline(slope = 0.009794, intercept = -0.268336)
图 2

4 Logit

当变量间关系不满足线性关系时，我们可以对原始数据进行一定的的转换，使其更加接近线性关系，继而再进行线性拟合。这里，我们可以对odds进行自然对数转换，将odds转化为以e为底的幂次：

odds = e^logit^

logit = log(odds)

对数转换后的变量为上表中的logit。log转换同样具有单调性，odds越大，logit越大。接下来，我们以logit为结果变量进行线性回归：

# 以breathlessness的logit为结果变量的线性回归
lm.fit2 <- lm(logit ~ age, subset = (breathlessness == 1), dat)
summary(lm.fit2)
## 
## Call:
## lm(formula = logit ~ age, data = dat, subset = (breathlessness == 
##     1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25431 -0.07937  0.04203  0.11581  0.14768 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.968911   0.178466  -39.05 1.88e-09 ***
## age          0.110338   0.004062   27.17 2.35e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1573 on 7 degrees of freedom
## Multiple R-squared:  0.9906, Adjusted R-squared:  0.9893 
## F-statistic:   738 on 1 and 7 DF,  p-value: 2.349e-08

可见logit = -6.968911 + 0.110338*age。年龄的系数是正的，表明随着年龄的增长，logit会增加，即breathlessness的odds会增加，即breathlessness的患病概率会增加。

我们可以根据该回归方程，使用age预测breathlessness的logit，估计得到相应年龄的breathlessness的logit。在下图中，我们将实际数据的logit (散点图) 与回归方程 (直线) 进行对比，可以发现logit与age的关系近乎直线。

ggplot(dat[dat$breathlessness == 1,], aes(age, logit)) +
  geom_point() +
  geom_abline(slope = 0.110338, intercept = -6.96891)
图 3

然而，线性回归假定结果变量的残差服从正态分布，这对于logit并不完美适用。既然我们将Y转换为了logit，更准确的统计方法是直接基于logit的分布 (即logistic累积概率分布) 来估计回归系数，这就是逻辑回归。另外，线性回归与逻辑回归的拟合算法不同，线性回归的拟合目标是使误差的平方和最小，而逻辑回归的拟合目标是使预测的正确率最大。这决定了两种算法所估计出来的回归方程大概率是不同的。

逻辑回归

逻辑回归就是以“逻辑值” (logit) 为“中间”结果变量的回归。

# 以breathlessness为结果变量的逻辑回归
glm.fit <- glm(
  factor(breathlessness) ~ age,
  family = binomial(link = "logit"),
  data = dat,
  weights = f
)
summary(glm.fit)
## 
## Call:
## glm(formula = factor(breathlessness) ~ age, family = binomial(link = "logit"), 
##     data = dat, weights = f)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.564333   0.124152  -52.87   <2e-16 ***
## age          0.102492   0.002454   41.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14318  on 17  degrees of freedom
## Residual deviance: 11993  on 16  degrees of freedom
## AIC: 11997
## 
## Number of Fisher Scoring iterations: 6

可见，logit = -6.564333 + 0.102492*age，该结果与前文得到的结果类似。

5 逻辑回归结果的解读

逻辑回归结果的解读分成三个层次：

  1. 回归系数的正负与显著性,
  2. odds, odds ratio，
  3. Y的概率。

我们从这三个层次上解读上述结果：

  1. 年龄的逻辑回归系数是正的，且达到了显著，表明年龄对于呼吸急促症状的发生概率具有显著的正向预测作用；

  2. 我们对回归方程等号两侧进行幂运算。一般的，对于逻辑回归方程Y = b0 + b1 *X，我们有：

    \[odds = e^Y = e^{b0 + b1 \cdot X} = e^{b0} \cdot e^{b1 \cdot X}\]

    在我们的案例中，我们可以得到：

    \[odds = e^{logit} = e^{-6.564333 + 0.102492 \cdot age}\]

    当X = 0时，\(odds_{X = 0} = e^{-6.564333}\)

    当X = 1时，\(odds_{X = 1} = e^{-6.564333 + 0.102492} = e^{-6.564333} \cdot e^{0.102492}\)

    \({Odds Ratio} = \frac {odds_{X = 1}} {odds_{X = 0}} = e^{0.102492} \approx 1.108\)

    即，age每增加一个单位，患上呼吸急促症状的odds会提高约1.108倍。

  3. 我们可以进一步计算概率的变化量。

    一般的，对于逻辑回归方程Y = b0 + b1*X，我们有：

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

    在上式中，函数F为标准logistic累积概率分布函数，该函数左边连着Y=1的概率，右边连着关于X的线性回归方程(b0 + b1*X)，F的作用是对关于X的线性回归方程进行非线性转换，将其转换为Y=1的概率。这里的函数\[\frac {1}{1 + e^{formula}}\]被称为连接函数(link function)。

连接函数

连接函数就是连接概率与线性回归方程的非线性函数。

那么，我们可以基于连接函数，使用年龄估计相应年龄患上呼吸急促症状的概率。注意，P(Y = 1| b0 + b1*X)与X的关系不是线性关系，对于不同水平的X，X增加一个单位所导致的P的变化量是不同的。下图呈现了年龄age_fake从0到100时患上呼吸急促症状的概率P_est。

age_fake <- seq(0, 100, 1)
P_est <- 1/(1 + exp(-1*(-6.564333 + 0.102492*age_fake)))
plot(age_fake, P_est)
图 4

我们以32岁与37岁为例，计算并比较其患上呼吸急促症状的概率：

# 比较32岁与37岁的概率
1/(1 + exp(-1*(-6.564333 + 0.102492*c(32, 37))))
## [1] 0.03610368 0.05884899

32岁时患上呼吸急促症状的概率估计为0.036，37岁时为0.059，患上呼吸急促症状的概率随年龄而增加。

6 Logistic累积概率分布

逻辑回归用logistic累积概率分布函数拟合X与Y=1的概率的关系。一般的，logistic累积概率分布函数为：

\[P = \frac {1}{1 + e^{-1 \cdot ( \mu + s \cdot X)}}\]

当 \(\mu = 0，s = 1\)时，该函数是标准logistic累积概率分布函数。

\(\mu\)与\(s\)的变化会给上述概率曲线带来哪些变化？我们作图演示。在下图中，P1是标准logistic累积概率分布。与P1相比，P2、P3的\(\mu\)的变化会导致图中曲线截距的变化，\(\mu\)越大，截距越大。P3、P4的\(s\)的变化会导致曲线的陡峭程度变化，\(s\)越大越陡峭。

data_fake <- 
  tibble(X = seq(-6, 6, 0.1)) |> 
  mutate(
    P1 = 1/(1 + exp(-1*X)),
    P2 = 1/(1 + exp(-1*(X - 2))),
    P3 = 1/(1 + exp(-1*(X - 4))),
    P4 = 1/(1 + exp(-1*(0.5*X))),
    P5 = 1/(1 + exp(-1*(2*X)))
  ) |> 
pivot_longer(
  cols = c(P1:P5),
  names_to = "logit_formula",
  values_to = "P"
) 
ggplot(data_fake, aes(X, P, color = logit_formula, shape = logit_formula)) +
  geom_point() +
  scale_shape_manual(values = c(16, 15, 17, 16, 16))

7 Probit回归

前文图4可以看出logistic累积概率分布图与正态累积概率分布图近似。我们可以将其作图比较。在下图中，横坐标为y*，纵坐标为累积概率P，红色表示标准logistic分布，绿色表示正态分布。

y_star <- seq(-6, 6, 0.1)
P_logistic <- 1/(1 + exp(-y_star))
P_norm <- pnorm(y_star)
dat_logistic_norm <- 
  tibble(y_star,P_logistic,P_norm) |> 
  pivot_longer(cols = c(P_logistic, P_norm),
               names_to = "dist",
               values_to = "P")
ggplot(dat_logistic_norm, aes(y_star, P, color = dist)) +
  geom_point()
图 5

所有分布的累积概率都是从0变成1，是否可以用正态累积概率函数作为连接函数 (link function ) 来拟合Y=0变成Y=1的概率变化呢？可以，这就是probit回归。一般的，对于probit回归方程，我们有：

\[P(Y=1|b0+b1 \cdot X) = \Phi(probit) = \Phi(b0+b1 \cdot X) = \int_{-\infty}^{b0+b1 \cdot X} \phi(z; 0, 1) dz\]

在上式中，\(\Phi(b0+b1 \cdot X)\)是连接函数，该函数是标准正态累积概率分布函数，该函数以\(probit = b0+b1 \cdot X\)为横坐标，求相应的正态分布累积概率。注意：逻辑回归方程的结果变量实际上是logit，不是Y，其中logit是logistic累积概率分布的横坐标；而probit回归方程的结果变量实际上是probit，也不是Y，其中probit是标准正态累积概率分布的横坐标，通过对\(P(Y=1|b0+b1 \cdot X)\)作\(\Phi\)的逆运算，即\(\Phi^{-1}(P)\)，我们可以计算得到\(probit\)：

\[Probit(P) = \Phi^{-1}(P)\]

有时，logit与probit也表示为\(Y^*\)，\(Y^*\)被称为Y的潜在变量 (underlying variable) 。尽管Y是二分变量，取值只能是0或者1，但每个case为0或者1的倾向是不同的，潜在变量\(Y^*\)就反映了这种倾向，\(Y^*\)越大，Y的取值越接近1、越可能是1。我们为\(Y^*\)设定一个阈限\(\tau\)，当\(Y^* > \tau\)时，我们就预测Y=1，当\(Y^* \leq \tau\)时，我们就预测Y=0。

Probit回归

Probit回归就是以probit为“中间”结果变量的回归。

我们对前文的案例作probit回归：

# 以breathlessness为结果变量的probit回归
glm.fit <- glm(
  factor(breathlessness) ~ age,
  weights = f,
  family = binomial(link = "probit"),
  dat
)
summary(glm.fit)
## 
## Call:
## glm(formula = factor(breathlessness) ~ age, family = binomial(link = "probit"), 
##     data = dat, weights = f)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.584158   0.061762  -58.03   <2e-16 ***
## age          0.054844   0.001272   43.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14318  on 17  degrees of freedom
## Residual deviance: 11981  on 16  degrees of freedom
## AIC: 11985
## 
## Number of Fisher Scoring iterations: 5

可得，Probit = -3.584158 + 0.054844*age。

8 Probit回归结果的解读

我们可以根据该方程计算probit值，将probit值作为标准正态分布的横坐标，然后计算相应的累积概率P_probit。

dat$probit <- -3.584158 + 0.054844*dat$age
dat$P_probit <- pnorm(dat$probit)
dat[1:9, c("probit", "P_probit")]
##     probit    P_probit
## 1 -2.37759 0.008713095
## 2 -2.10337 0.017716719
## 3 -1.82915 0.033688570
## 4 -1.55493 0.059981371
## 5 -1.28071 0.100147773
## 6 -1.00649 0.157089960
## 7 -0.73227 0.232001893
## 8 -0.45805 0.323458260
## 9 -0.18383 0.427073410

9 Logistic与probit回归系数的估计方法

前文提到线性回归的拟合目标是使预测Y的误差平方和最小。Probit回归与Logistic回归一样，其拟合目标都是使预测Y的准确率（似然函数）最大，即，使下面的函数的值最大（下面的函数是似然函数的对数函数）：

\[logL = \sum_i^n \mu_i logP_i + (1-\mu_i)log(1-P_i)\]

其中，\(P_i=Prob(\mu_i = 1|x_i)\)。

回归系数为多少时上述函数的值最大？幸运地是，寻找最优解的这一过程将由计算机去完成。

10 阈限 \(\tau\)

lavaan与Mplus一样，当结果变量是分类变量时默认使用probit回归。与传统的逻辑回归或者probit回归不同，lavaan与Mplus的计算结果中没有报告截距，而是报告了阈限Thresholds：

model <- "
breathlessness ~ age
"
fit <- sem(model, dat, ordered = "breathlessness", sampling.weights = "f")
summary(fit)
## lavaan 0.6-19 ended normally after 2 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                         2
## 
##   Number of observations                            18
##   Sampling weights variable                          f
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 0.000       0.000
##   Degrees of freedom                                 0           0
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   breathlessness ~                                    
##     age               0.055    0.025    2.192    0.028
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     brethlssnss|t1    3.584    1.096    3.270    0.001
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .breathlessness    1.000

如何解释这种差异？事实上，阈限Thresholds = 截距*(-1)。

当结果变量Y是二分变量时，逻辑回归方程的结果变量是logit，logit是Y的潜在变量，是一个连续变量，当logit高于某个阈限\(\tau\)时，Y的预测值将为1，当logit低于这个阈限\(\tau\)时，Y的预测值将为1，logit通常是基于标准logistic分布，\(\tau\)通常被设为0。Probit回归与此类似，Probit通常基于标准正态分布，\(\tau\)通常也被设为0。可以看出，实际上\(\tau\)这个参数具有意义、易于理解。

lavaan与Mplus采用的参数设定过程可以理解为“参数重设”的过程：

第一步，基于标准logistic累积概率分布或标准正态累积概率分布，设定\(\tau = 0\)，估计得到逻辑回归方程或probit回归的截距为b0，b0非0；

第二步，设定\(\tau = -b0\)，此时估计得到的截距就是0！

观察图5，我们可以将上述过程想象为坐标系不动、将累积概率分布图向左平移b0个单位，或者累积概率分布图不动、将纵坐标向右平移b0个单位，使得分布中心由\(\tau = 0\)移动到新的位置\(\tau = -b0\)。

11 拟合程度的比较

下面，我们基于前文的线性回归、logistic回归、probit回归分别估计相应的概率P_linear(红色)，P_logistic（黄色），P_probit（绿色），并将其与实际的P（蓝色点）进行比较。可见，线性回归的拟合较差，logistic回归与probit回归对实际数据的拟合较好。logistic回归与probit回归的拟合程度可以通过AIC、BIC指标来比较，AIC、BIC越小，其拟合越好。通过对前文报告的AIC进行比较，我们可以看出probit对数据的拟合更好。一般情况下，probit回归对数据的拟合程度优于logistic回归。

# label: fig-6
dat$P_linear <- -0.268336 + 0.009794*dat$age
dat$P_logistic <- 1/(1 + exp(-1*(-6.564333 + 0.102492*dat$age)))
ggplot(dat[dat$breathlessness == 1, ], aes(x = age)) + 
  geom_path(aes(y = P_linear), color = "red", linewidth = 1) +
  geom_path(aes(y = P_logistic), color = "yellow", linewidth = 1) +
  geom_path(aes(y = P_probit), color = "green", linewidth = 1) +
  geom_point(aes(y = P), color = "blue", size = 3)

12 参考文献

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