使用 Survey 包对堆叠插补中的观察值进行加权
Use the Survey package to weight observations in stacked imputations
我正在探索估算数据中的模型变量选择。
一种技术是以长格式堆叠插补（其中 M 个插补数据集中的 n 个观察创建一个 n x M 长的数据集），并使用加权回归按比例减少每个观察对插补数量的贡献。如果我们将堆叠数据集视为一个大数据集，标准误差会太小。
我正在尝试使用 svyglm 中的 weights 参数来解释堆叠数据，从而产生您期望的 SE，即 n 个观测值，而不是 n x 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)
长数据给出的 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 包允许在乘法估算数据集和模型池中进行逐步选择。然而，它是计算密集型的，并且对于大型数据集来说速度很慢，特别是如果需要引导过程（例如在内部验证期间）——因此需要不太密集的堆叠方法。
抱歉，不，这行不通。
要处理具有权重的堆叠插补数据，您需要 频率权重，因此 1/10 的权重意味着您有 1/10 的观察值。使用 svydesign，您指定 抽样权重 ，因此 1/10 的权重意味着您的观察代表总体中的 10 个观察。这些将（并且应该）给出不同的标准错误。当你实际上有插补时假装你有频率权重是一个聪明的技巧，可以避免让软件理解它在做什么，这很好但与 survey 不兼容，后者理解它在做什么并且正在做一些不同的事情。
目前，如果您想将 svyglm 与多重插补一起使用，您需要分别计算标准误差——最方便的是使用 mitools::MIcombine 的 Rubin 规则，它被设置为与调查包（请参阅 with.svyimputationList 和 withPV 的帮助）。
可能值得向 mitools 或 survey 开发人员提出功能请求（引用示例）以允许对插补进行堆叠分析，但这不仅仅是一个问题调整权重。
我正在探索估算数据中的模型变量选择。
一种技术是以长格式堆叠插补（其中 M 个插补数据集中的 n 个观察创建一个 n x M 长的数据集），并使用加权回归按比例减少每个观察对插补数量的贡献。如果我们将堆叠数据集视为一个大数据集，标准误差会太小。
我正在尝试使用 svyglm 中的 weights 参数来解释堆叠数据，从而产生您期望的 SE，即 n 个观测值，而不是 n x 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)
长数据给出的 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 包允许在乘法估算数据集和模型池中进行逐步选择。然而，它是计算密集型的，并且对于大型数据集来说速度很慢，特别是如果需要引导过程（例如在内部验证期间）——因此需要不太密集的堆叠方法。
抱歉，不，这行不通。
要处理具有权重的堆叠插补数据，您需要 频率权重，因此 1/10 的权重意味着您有 1/10 的观察值。使用 svydesign，您指定 抽样权重 ，因此 1/10 的权重意味着您的观察代表总体中的 10 个观察。这些将（并且应该）给出不同的标准错误。当你实际上有插补时假装你有频率权重是一个聪明的技巧，可以避免让软件理解它在做什么，这很好但与 survey 不兼容，后者理解它在做什么并且正在做一些不同的事情。
目前，如果您想将 svyglm 与多重插补一起使用，您需要分别计算标准误差——最方便的是使用 mitools::MIcombine 的 Rubin 规则，它被设置为与调查包（请参阅 with.svyimputationList 和 withPV 的帮助）。
可能值得向 mitools 或 survey 开发人员提出功能请求（引用示例）以允许对插补进行堆叠分析，但这不仅仅是一个问题调整权重。