marginaleffects风险差异与LPM概率变化是否等同及计算逻辑问询

问题背景与疑问 #

我希望预测泰坦尼克号沉没后不同社会阶层的存活概率。由于这并非随机对照试验（RCT），各组特征难以均衡，因此需要纳入性别、年龄因素以控制存活差异。我采用logit模型进行分析，但发现优势比难以向非专业人士解释，而风险差异和风险比更易理解。现提出以下问题：

1. 风险差异是否等同于概率差，且近似于线性概率模型（LPM）的系数（即概率变化）？
2. 若二者等同，marginaleffects包如何计算平均差异？是假设特定个体特征（如男性、特定年龄），还是取特征均值（如mean(Age)）并使用predict()函数？如何单独计算与PClass相关的概率？

解答 #

问题1：风险差异与LPM系数的关系 #

• 风险差异（Risk Difference, RD）就是两组事件发生概率的差值，比如上层阶级存活概率减去下层阶级存活概率，和你理解的“概率差”完全一致。
• 它和线性概率模型（LPM）的系数含义近似，但并非完全等同：
  • LPM直接在概率层面建模，系数代表自变量每变化一个单位，事件发生概率的固定变化值；
  • logit模型推导的风险差异是基于非线性转换后的概率计算的差值，这个差值会随个体特征（如年龄、性别）变化而不同。avg_comparisons()计算的是平均风险差异——对每个个体计算风险差异后取均值，这个结果和LPM的系数会比较接近，但当存活概率接近0或1时，logit模型的非线性会让二者出现明显差距，而且LPM可能出现预测概率超出0-1范围的问题，logit模型则不会。

问题2：marginaleffects包的计算逻辑与单独概率计算 #

• avg_comparisons()计算平均风险差异的步骤：
  1. 遍历数据集中的每一个个体，用该个体自身的性别、年龄等特征，通过logit模型的predict()函数，分别计算其在不同Pclass下的存活概率；
  2. 对每个个体，计算Pclass不同类别之间的概率差值；
  3. 将所有个体的差值取平均值，得到最终的平均风险差异。
     这个过程既不用假设特定特征，也不用取特征均值，而是基于全样本个体的真实特征逐一计算后平均，属于「平均边际效应（AME）」的范畴。
• 要单独计算各Pclass对应的存活概率，可以使用avg_predictions()函数，按Pclass分组计算平均预测概率：

  # 未调整模型的各阶层存活概率
  avg_predictions(lrm_unadj, variables = "Pclass")
  
  # 调整性别、年龄后的各阶层存活概率
  avg_predictions(lrm_adj, variables = "Pclass")
  

原始分析代码 #

library(tidyverse)
## Load Titanic library to get the dataset
library(titanic)
library(marginaleffects)

## Load the datasets
data("titanic_train")
data("titanic_test")

## Setting Survived column for test data to NA
titanic_test$Survived <- NA

## Combining Training and Testing dataset
complete_data <- rbind(titanic_train, titanic_test) %>% 
  filter(!Survived == "Unknown") %>% 
  mutate(Pclass = factor(Pclass),
         Sex = factor(Sex),
         Survived = factor(Survived)) 


lrm_unadj <- glm(Survived ~ Pclass,
                 data = complete_data, family = binomial(link = "logit"))

lrm_adj <- glm(Survived ~ Pclass + Sex + Age,
               data = complete_data, family = binomial(link = "logit"))


# Risk difference (difference in probabilities) - should this be similar to LPM?
avg_comparisons(lrm_unadj, variables = "Pclass")
avg_comparisons(lrm_adj, variables = "Pclass") 

lm(as.numeric(Survived) ~  Pclass, data = complete_data)%>% 
lmtest::coeftest(., vcov = sandwich::vcovHC, type = "HC3") %>% 
  broom::tidy(conf.int = T)

lm(as.numeric(Survived) ~  Pclass + Sex + Age, data = complete_data)%>% 
lmtest::coeftest(., vcov = sandwich::vcovHC, type = "HC3") %>% 
  broom::tidy(conf.int = T)

内容的提问来源于stack exchange，提问作者allen.joseph