使用交互项对多重插补建模

Question

根据mice 包的文档，如果我们想在对交互项感兴趣时插补数据，我们需要使用被动插补。这是通过以下方式完成的。

library(mice)
nhanes2.ext <- cbind(nhanes2, bmi.chl = NA)
ini <- mice(nhanes2.ext, max = 0, print = FALSE)

meth <- ini$meth
meth["bmi.chl"] <- "~I((bmi-25)*(chl-200))"

pred <- ini$pred
pred[c("bmi", "chl"), "bmi.chl"] <- 0

imp <- mice(nhanes2.ext, meth = meth, pred = pred, seed = 51600, print = FALSE)

据说

Imputations created in this way preserve the interaction of bmi with chl

这里，在原始数据集中创建了一个名为bmi.chl 的新变量。 meth 步骤说明需要如何从现有变量中估算此变量。 pred 步骤表示我们不想从 bmi.chl 预测 bmi 和 chl。但是现在，如果我们想应用一个模型，我们该如何进行呢？ "~I((bmi-25)*(chl-200))"定义的乘积是否只是一种控制主效应估算值的方法，即bmi和chl？

如果我们要拟合的模型是glm(hyp~chl*bmi, family="binomial")，那么从估算数据中指定该模型的正确方法是什么？ fit1 或 fit2?

fit1 <- with(data=imp, glm(hyp~chl*bmi, family="binomial"))
summary(pool(fit1))

或者我们是否必须以某种方式使用创建的新变量的估算值，即bmi.chl？

fit2 <- with(data=imp, glm(hyp~chl+bmi+bmi.chl, family="binomial"))
summary(pool(fit2))

Answer 1

对于被动插补，您是使用被动插补变量还是在调用glm 时重新计算乘积项都没有关系。 fit1 和 fit2 在您的示例中产生不同结果的原因是不是只是对乘积项进行被动插补。

相反，您在相乘之前转换了两个变量（即，您计算 bmi-25 和 chl-100）。因此，被动插补变量bmi.chl 不代表乘积项bmi*chl，而是代表(bmi-25)*(chl-200)。

如果您只计算乘积项，那么 fit1 和 fit2 会产生应有的相同结果：

library(mice)
nhanes2.ext <- cbind(nhanes2, bmi.chl = NA)
ini <- mice(nhanes2.ext, max = 0, print = FALSE)

meth <- ini$meth
meth["bmi.chl"] <- "~I(bmi*chl)"

pred <- ini$pred
pred[c("bmi", "chl"), "bmi.chl"] <- 0
pred[c("hyp"), "bmi.chl"] <- 1

imp <- mice(nhanes2.ext, meth = meth, pred = pred, seed = 51600, print = FALSE)

fit1 <- with(data=imp, glm(hyp~chl*bmi, family="binomial"))
summary(pool(fit1))
# > round(summary(pool(fit1)),2)
#                est    se     t    df Pr(>|t|)   lo 95 hi 95 nmis  fmi lambda
# (Intercept) -23.94 38.03 -0.63 10.23     0.54 -108.43 60.54   NA 0.41   0.30
# chl           0.10  0.18  0.58  9.71     0.58   -0.30  0.51   10 0.43   0.32
# bmi           0.70  1.41  0.49 10.25     0.63   -2.44  3.83    9 0.41   0.30
# chl:bmi       0.00  0.01 -0.47  9.67     0.65   -0.02  0.01   NA 0.43   0.33

fit2 <- with(data=imp, glm(hyp~chl+bmi+bmi.chl, family="binomial"))
summary(pool(fit2))
# > round(summary(pool(fit2)),2)
#                est    se     t    df Pr(>|t|)   lo 95 hi 95 nmis  fmi lambda
# (Intercept) -23.94 38.03 -0.63 10.23     0.54 -108.43 60.54   NA 0.41   0.30
# chl           0.10  0.18  0.58  9.71     0.58   -0.30  0.51   10 0.43   0.32
# bmi           0.70  1.41  0.49 10.25     0.63   -2.44  3.83    9 0.41   0.30
# bmi.chl       0.00  0.01 -0.47  9.67     0.65   -0.02  0.01   25 0.43   0.33

这并不奇怪，因为mice 中的~I(bmi*chl) 和glm 中的bmi*chl 执行完全相同的操作。他们只是计算两个变量的乘积。

备注：

请注意，我添加了一行，说明在估算 hyp 时应使用 bmi.chl 作为预测变量。如果没有这一步，被动插补就没有任何意义，因为插补模型会忽略乘积项，从而与分析模型不一致。