使用ri2包结合glm二项式模型的随机化推断问题

问题背景 #

我在针对二分类模型使用ri2包执行随机化推断流程时，通过glm指定了binomial(link='logit')族，但ri2似乎未识别该族，得到的系数与未指定族（默认高斯族，等价于使用lm）的结果相近。想了解以下三个问题：

1. 这种情况下随机化推断对我指定的模型是否有效；
2. 能否通过ri2对比处理组系数的p值来检验glm中处理组变量的显著性；
3. 是否存在让ri2正确识别glm族参数的方法。

复现代码 #

library(dplyr)
library(lmtest)
library(ri2)
library(fastDummies) # 补充加载依赖包

set.seed(123)
n <- 100

y <- rbinom(100, size = 1, prob = 0.5)

ngroup = 4
nrep = 15

treat <- rep(c("control", "treat1", "treat2", "treat3"), each = nrep)
strata <- sample(1:3, 100, replace = TRUE)
cluster <- sample(1:5, 100, replace = TRUE)

test_1 <- as.data.frame(cbind(treat, strata, cluster, y)) %>%
  mutate(y = as.numeric(y),
         cluster = as.numeric(cluster)) %>%
  fastDummies::dummy_cols(select_columns = "treat") %>%
  fastDummies::dummy_cols(select_columns = "strata")

### glm with family = binomial(link = "logit")
m1 <- glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, 
          data = test_1, family = binomial(link = "logit"))

### glm with a default family
m2 <- glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, 
          data = test_1)

### results with ri2
declaration <- declare_ra(N = 100, clusters = test_1$cluster)

ri2 <- 
  conduct_ri(
    formula = glm(y~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3, 
                  data = test_1, family = binomial(link = "logit")),
    assignment = "treat_treat1",
    declaration = declaration,
    data = test_1,
    sims = 1000
  )

ri2
coeftest(m1)
coeftest(m2)

问题解答 #

1. 随机化推断对指定模型是否有效 #

ri2的核心逻辑是基于随机化分配的潜在结果框架，有效性完全依赖于处理分配的随机性，而非模型的分布假设。即使你在conduct_ri()的formula里指定了logit族，ri2默认会使用lm_estimator()（线性回归）来估计处理效应——这就是你看到系数和高斯族结果相近的原因。

这种情况下，随机化推断仍然有效：OLS估计量在随机化实验中对平均处理效应（ATE）是无偏的，哪怕结果是二分类的。但要明确，此时你得到的是线性概率模型（LPM）的处理效应（即概率变化量），而非logit模型输出的对数优势比。

2. 能否用ri2的p值检验glm处理组变量的显著性 #

可以，但需要区分两种场景：

• 如果你用ri2默认的OLS估计器，得到的p值是基于LPM处理效应的随机化分布，用来检验的是“处理对结果概率的线性影响是否为0”；
• 如果你想检验logit模型中处理变量的对数优势比是否为0，需要自定义估计器，让ri2在每次模拟中拟合logit模型并提取对应系数，再基于该系数的随机化分布计算p值。

3. 让ri2正确识别glm族参数的方法 #

ri2支持通过estimator参数自定义估计逻辑，你可以编写一个函数来拟合指定族的glm模型，并提取目标系数。具体实现如下：

步骤1：定义自定义logit估计器

# 自定义估计器：拟合logit模型并提取treat_treat1的系数
logit_estimator <- function(data, formula, ...) {
  model <- glm(formula, data = data, family = binomial(link = "logit"))
  coef(model)["treat_treat1"] # 返回目标处理变量的系数
}

步骤2：在conduct_ri中使用自定义估计器

ri2_logit <- conduct_ri(
  formula = y ~ treat_treat1 + treat_treat2 + treat_treat3 + strata_2 + strata_3,
  assignment = "treat_treat1",
  declaration = declaration,
  data = test_1,
  sims = 1000,
  estimator = logit_estimator # 指定自定义估计器
)

# 查看结果
ri2_logit

这样，ri2会在每次随机化模拟中重新拟合logit模型，提取处理变量的对数优势比，最终生成基于该系数的随机化分布，并计算对应的p值。如果需要同时检验多个处理变量的系数，可以修改估计器函数返回系数向量，比如：

logit_multi_estimator <- function(data, formula, ...) {
  model <- glm(formula, data = data, family = binomial(link = "logit"))
  coef(model)[c("treat_treat1", "treat_treat2", "treat_treat3")]
}

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