使用Survey包为堆叠多重插补数据加权校正标准误问题咨询

研究背景 #

我正在开展插补数据的模型变量选择相关研究，其中一种实现方案是将多重插补结果转换为长格式堆叠数据：M个插补数据集的n条原始观测会生成规模为n×M的长数据集，此时需要通过加权回归按插补次数比例降低单条观测的贡献，否则直接对堆叠数据集建模会得到偏小的标准误。
我尝试使用svyglm的weights参数处理堆叠数据，期望输出对应n条原始观测而非n×M条观测的标准误，具体实现过程如下：

复现代码 #

library(mice)

### create data
set.seed(42)
n <- 50
id <- 1:n
var1 <- rbinom(n,1,0.4)
var2 <- runif(n,30,80)
var3 <- rnorm(n, mean = 12, sd = 5)
var4 <- rnorm(n, mean = 100, sd = 20)
prob <- (((var1*var2)+var3)-min((var1*var2)+var3)) / (max((var1*var2)+var3)-min((var1*var2)+var3))
outcome <- rbinom(n, 1, prob = prob)
data <- data.frame(id, var1, var2, var3, var4, outcome)

### Add missingness
data_miss <- ampute(data)
patt <- data_miss$patterns
patt <- patt[2:5,]
data_miss <- ampute(data, patterns = patt)
data_miss <- data_miss$amp

## create 5 imputed datasets
nimp <- 5
imp <- mice(data_miss, m = nimp)

## Stack data
data_long <- complete(imp, action = "long")

## Generate model in stacked data (SEs will be too small)
modlong <- glm(outcome ~ var1 + var2 + var3 + var4, family = "binomial", data = data_long)
summary(modlong)

直接建模输出（标准误偏小） #

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.906417   0.965090  -3.012   0.0026 ** 
var1         2.221053   0.311167   7.138 9.48e-13 ***
var2        -0.002543   0.010468  -0.243   0.8081    
var3         0.076955   0.032265   2.385   0.0171 *  
var4         0.006595   0.008031   0.821   0.4115    

加权回归尝试代码 #

data_long$weight <- 1/nimp

library(survey)
des <- svydesign(ids = ~1, data = data_long, weights = ~weight)
mod_svy <- svyglm(formula = outcome ~ var1 + var2 + var3 + var4, family = quasibinomial(), design = des)

summary(mod_svy)

加权回归输出（标准误仍偏小） #

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.906417   1.036691  -2.804  0.00546 ** 
var1         2.221053   0.310906   7.144 1.03e-11 ***
var2        -0.002543   0.010547  -0.241  0.80967    
var3         0.076955   0.030955   2.486  0.01358 *  
var4         0.006595   0.008581   0.769  0.44288    

添加rescale = F参数禁止权重重缩放至样本量总和后，结果无变化：

mod_svy <- svyglm(formula = outcome ~ var1 + var2 + var3 + var4, family = quasibinomial(), design = des, rescale = F)
summary(mod_svy)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.906417   1.036688  -2.804  0.00546 ** 
var1         2.221053   0.310905   7.144 1.03e-11 ***
var2        -0.002543   0.010547  -0.241  0.80967    
var3         0.076955   0.030955   2.486  0.01358 *  
var4         0.006595   0.008581   0.769  0.44288 

预期参考结果 #

加权后得到的标准误应与单个插补数据集建模结果接近，单个插补数据集建模输出如下：

## Assess SEs in single imputation
mod_singleimp <-  glm(outcome ~ var1 + var2 + var3 + var4, family = "binomial", data = complete(imp,1))
summary(mod_singleimp)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.679589   2.116806  -1.266  0.20556   
var1         2.476193   0.761195   3.253  0.00114 **
var2         0.014823   0.025350   0.585  0.55874   
var3         0.048940   0.072752   0.673  0.50114   
var4        -0.004551   0.017986  -0.253  0.80026  

方案需求 #

• 希望获得该标准误偏小问题的解决思路
• 也欢迎提供其他可实现相同目标的轻量化方案

现有可选方案中，psfmi包支持多重插补数据集的逐步选择与模型合并，但在大数据集场景下计算负载高、运行速度慢，尤其是需要执行自助法内部验证时耗时更长，因此不满足使用需求。

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