使用交互项对多重插补建模
Model multiple imputation with interaction terms
根据 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))
对于被动插补,无论您是使用被动插补变量,还是在调用 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
用作预测变量。如果没有这一步,被动插补就没有意义,因为插补模型会忽略乘积项,从而与分析模型不一致。
根据 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))
对于被动插补,无论您是使用被动插补变量,还是在调用 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
用作预测变量。如果没有这一步,被动插补就没有意义,因为插补模型会忽略乘积项,从而与分析模型不一致。