使用 Survey 包对堆叠插补中的观察结果进行加权答案

【问题标题】：Use the Survey package to weight observations in stacked imputations使用 Survey 包对堆叠插补中的观察结果进行加权
【发布时间】：2021-11-09 17:02:08
【问题描述】：

我正在探索估算数据中的模型变量选择。

一种技术是以长格式堆叠插补（其中 M 个插补数据集中的 n 个观察创建一个长 n x M 的数据集），并使用加权回归按比例减少每个观察对插补数量的贡献。如果我们将堆叠的数据集视为一个大数据集，标准误差会太小。

我正在尝试使用 svyglm 中的 weights 参数来解释堆叠数据，从而产生您期望的具有 n 个观察而不是 n x M 观察的 SE。

举例说明：

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)

长数据导致 SE 过小，因为我们将数据集的大小增加了 5 倍

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)

加权回归给出了与未加权模型相似的 SE

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

我希望 SE 类似于在单个估算数据集中运行模型时获得的 SE

## 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 包允许在多重插补数据集和模型池中逐步选择。但是，对于大型数据集，它的计算量很大并且速度很慢，特别是如果需要引导流程（例如在内部验证期间） - 因此需要一种不太密集的堆叠方法。

【问题讨论】：

标签： glm survey imputation weighted

【解决方案1】：

抱歉，不，这行不通。

要处理具有权重的堆叠插补数据，您需要频率权重，因此 1/10 的权重意味着您有 1/10 的观察值。使用svydesign，您可以指定抽样权重，因此 1/10 的权重意味着您的观察代表总体中的 10 个观察。这些将（并且应该）给出不同的标准错误。假装你有频率权重，而实际上你有插补是一个聪明的技巧，以避免让软件理解它在做什么，这很好，但与 survey 不兼容，它知道它在做什么并且正在做不同的事情。

目前，如果您想将svyglm 与多个插补一起使用，您需要分别计算标准误——最方便的是使用鲁宾规则使用mitools::MIcombine，它设置为与调查包一起使用（请参阅帮助with.svyimputationList 和 withPV)。

可能值得向mitools 或survey 开发人员提出功能请求（并引用示例）以允许对插补进行堆叠分析，但这不仅仅是调整权重的问题。