多预测变量逻辑回归中连续变量交互效应的可视化（R实现）

问题描述 #

我有一个包含4个预测变量（CAPE、SATmax、PC1、PC2）、交互项CAPE:SATmax与二分类响应变量extry的数据集，想要可视化CAPE的斜率随SATmax取值的变化。已知线性回归下有类似实现方法，现在需要在多预测变量的Logistic回归中仅关注该交互项的实现方式，使用R语言，希望能提供示例代码。

数据集子集如下：

test <- structure(list(extry = c(0, 0, 0, 0, 0, 0), CAPE = c(1.14374437306378, 
20.1164451497721, 4.50769841148758, 26.3733996612395, 2.28748874612802, 
0.269116323073831), SATmax = c(11.4, 8.9, 10.3, 12.3, 9.6, 11.7
), PC1 = c(-2.65976813816683, -2.5478670787521, -2.58556444360627, 
-2.39983714790594, -1.51378585899909, -0.695703681902304), PC2 = c(0.37508773845025, 
-0.347147505858787, -1.35864530847998, -1.66587526284572, -1.86195614101445, 
-1.75492325860976)), row.names = c(NA, 6L), class = "data.frame")

已拟合的模型：

fit_int <- glm(extry~CAPE+ SATmax + PC1 + PC2+ CAPE:SATmax, family=binomial("logit"), data = test)

实现方法 #

在Logistic回归中，交互项CAPE:SATmax意味着CAPE对logit(响应概率)的影响随SATmax的取值变化。我们可以通过计算不同SATmax取值下CAPE的斜率（即logit尺度上的边际效应），再将其可视化。

方法1：使用emmeans包快速计算并可视化 #

emmeans包可以便捷计算模型在不同协变量取值下的边际效应，步骤如下：

1. 安装并加载所需包：

install.packages("emmeans")
library(emmeans)
library(ggplot2)

2. 计算不同SATmax取值下CAPE的斜率（logit尺度）：
   选择覆盖SATmax数据范围的序列值，计算每个取值下CAPE的斜率及置信区间：

# 生成SATmax的序列值
sat_seq <- seq(min(test$SATmax), max(test$SATmax), length.out = 20)

# 计算每个SATmax下CAPE的斜率
slopes <- emtrends(fit_int, ~ SATmax, var = "CAPE", at = list(SATmax = sat_seq))

# 转换为数据框方便绘图
slopes_df <- as.data.frame(slopes)

3. 可视化斜率变化：

ggplot(slopes_df, aes(x = SATmax, y = CAPE.trend)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = 0.2, fill = "blue") +
  labs(title = "CAPE的斜率随SATmax的变化（Logit尺度）",
       x = "SATmax",
       y = "CAPE的斜率（Logit尺度）") +
  theme_minimal()

方法2：手动计算斜率（理解底层逻辑） #

Logistic回归的模型公式为：
$$\text{logit}(P(extry=1)) = \beta_0 + \beta_1 CAPE + \beta_2 SATmax + \beta_3 PC1 + \beta_4 PC2 + \beta_5 CAPE \times SATmax$$
CAPE的斜率为 $\beta_1 + \beta_5 \times SATmax$，可直接提取模型系数计算：

1. 提取模型系数与方差协方差矩阵：

coefs <- coef(fit_int)
beta1 <- coefs["CAPE"]
beta5 <- coefs["CAPE:SATmax"]
vcov_mat <- vcov(fit_int)

2. 计算斜率及置信区间：

# 计算每个SATmax对应的斜率方差
slope_var <- vcov_mat["CAPE", "CAPE"] + sat_seq^2 * vcov_mat["CAPE:SATmax", "CAPE:SATmax"] + 
  2 * sat_seq * vcov_mat["CAPE", "CAPE:SATmax"]

# 计算斜率与95%置信区间
slope <- beta1 + beta5 * sat_seq
lower_ci <- slope - 1.96 * sqrt(slope_var)
upper_ci <- slope + 1.96 * sqrt(slope_var)

# 构建数据框
slopes_manual_df <- data.frame(SATmax = sat_seq,
                               slope = slope,
                               lower_ci = lower_ci,
                               upper_ci = upper_ci)

3. 可视化：

ggplot(slopes_manual_df, aes(x = SATmax, y = slope)) +
  geom_line(color = "red") +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, fill = "red") +
  labs(title = "CAPE的斜率随SATmax的变化（手动计算，Logit尺度）",
       x = "SATmax",
       y = "CAPE的斜率（Logit尺度）") +
  theme_minimal()

补充：概率尺度的斜率可视化 #

如果需要观察CAPE对响应概率的边际效应随SATmax的变化，可使用marginaleffects包：

install.packages("marginaleffects")
library(marginaleffects)

# 计算不同SATmax下CAPE的边际效应（概率尺度，固定PC1、PC2为均值）
me <- marginaleffects(fit_int, 
                      variables = "CAPE",
                      newdata = datagrid(SATmax = sat_seq, 
                                         PC1 = mean(test$PC1), 
                                         PC2 = mean(test$PC2)))

# 可视化
ggplot(me, aes(x = SATmax, y = dydx)) +
  geom_line(color = "darkgreen") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, fill = "darkgreen") +
  labs(title = "CAPE对extry概率的边际效应随SATmax的变化",
       x = "SATmax",
       y = "CAPE的边际效应（概率尺度）") +
  theme_minimal()

说明 #

• Logit尺度的斜率表示CAPE每增加1单位，logit(响应概率)的变化量，该变化量随SATmax线性变化；
• 概率尺度的边际效应表示CAPE每增加1单位，响应概率的变化量，它是非线性的，且依赖于其他变量的取值；
• 示例使用的是数据集子集，实际使用时请替换为完整数据集（子集响应变量全为0，模型无法拟合出有效系数）。

内容的提问来源于stack exchange，提问作者lola