使用具有聚类数据的小鼠进行插补
Imputation using mice with clustered data
所以我使用 mice
包来估算缺失数据。我是插补的新手,所以我已经到了一定程度,但是 运行 进入了陡峭的学习曲线。举个玩具例子:
library(mice)
# Using nhanes dataset as example
df1 <- mice(nhanes, m=10)
如您所见,我使用大部分默认设置对 df1 进行了 10 次估算 - 我很乐意在回归模型、合并结果等中使用此结果。但是在我的现实生活数据中,我有来自不同国家/地区的调查数据。因此,缺失水平因国家/地区而异,特定变量的值(即年龄、教育水平等)也不同。因此,我想估算缺失值,允许按国家/地区进行聚类。所以我将创建一个没有缺失的分组变量(当然在这个玩具示例中缺少与其他变量的相关性,但在我的真实数据中它们存在)
# Create a grouping variable
nhanes$country <- sample(c("A", "B"), size=nrow(nhanes), replace=TRUE)
那么我如何告诉 mice()
这个变量与其他变量不同 - 即它是多级数据集中的一个级别?
您必须设置一个 predictorMatrix 来告诉老鼠使用哪个变量来估算另一个变量。一个快速的方法是使用 predictorM<-quickpred(nhanes)
然后将矩阵中的 1s 更改为 2,如果它是一个普通变量,如果它是不同国家的二级变量,则将其更改为 -2,并将其作为 predictorMatrix =predictorM
提交给 mice 命令。在方法命令中,您现在必须将方法设置为 2l.norm
(如果它是度量变量)或 2l.binom
(如果它是二进制变量)。对于后者,您需要 Sabine Zinn (https://www.neps-data.de/Portals/0/Working%20Papers/WP_XXXI.pdf) 编写的函数。不幸的是,我不知道世界上是否存在用于估算两级计数数据的方法。
请注意,对多级数据集进行估算会大大减慢该过程。根据我的经验,像 PMM 或 Baboon 包中的重采样方法在保持数据的层次结构方面效果很好,并且使用起来要快得多。
如果您将集群视为 "mixed-effects" 模型,那么您应该使用 mice
提供的用于集群数据的方法。这些方法可以在 manual 中找到,通常以 2l.something
.
为前缀
聚类数据方法的多样性在 mice
中有所限制,但我可以推荐使用 2l.pan
来处理较低级别单元中的缺失数据,并在聚类级别使用 2l.only.norm
.
作为混合效应模型的替代方案,您可以考虑使用虚拟指标来表示集群结构(即,每个集群一个虚拟变量)。当您从混合效应模型的角度考虑集群时,此方法并不理想。因此,如果您想进行混合效应分析,请尽可能坚持使用混合效应模型。
下面,我展示了这两种策略的示例。
准备:
library(mice)
data(nhanes)
set.seed(123)
nhanes <- within(nhanes,{
country <- factor(sample(LETTERS[1:10], size=nrow(nhanes), replace=TRUE))
countryID <- as.numeric(country)
})
案例 1:使用混合效应模型进行插补
本节使用 2l.pan
来估算具有缺失数据的三个变量。请注意,我通过在预测矩阵中指定 -2
来使用 clusterID
作为聚类变量。对于所有其他变量,我仅分配固定效应 (1
)。
# "empty" imputation as a template
imp0 <- mice(nhanes, maxit=0)
pred1 <- imp0$predictorMatrix
meth1 <- imp0$method
# set imputation procedures
meth1[c("bmi","hyp","chl")] <- "2l.pan"
# set predictor Matrix (mixed-effects models with random intercept
# for countryID and fixed effects otherwise)
pred1[,"country"] <- 0 # don't use country factor
pred1[,"countryID"] <- -2 # use countryID as cluster variable
pred1["bmi", c("age","hyp","chl")] <- c(1,1,1) # fixed effects (bmi)
pred1["hyp", c("age","bmi","chl")] <- c(1,1,1) # fixed effects (hyp)
pred1["chl", c("age","bmi","hyp")] <- c(1,1,1) # fixed effects (chl)
# impute
imp1 <- mice(nhanes, maxit=20, m=10, predictorMatrix=pred1, method=meth1)
案例 2:使用集群的虚拟指标 (DI) 进行插补
本节使用pmm
插补,聚类结构用"ad hoc"表示。也就是说,聚类不是由随机效应表示,而是由固定效应表示。这可能会夸大具有缺失数据的变量的集群级可变性,因此请确保您在使用它时知道自己在做什么。
# create dummy indicator variables
DIs <- with(nhanes, contrasts(country)[country,])
colnames(DIs) <- paste0("country",colnames(DIs))
nhanes <- cbind(nhanes,DIs)
# "empty" imputation as a template
imp0 <- mice(nhanes, maxit=0)
pred2 <- imp0$predictorMatrix
meth2 <- imp0$method
# set imputation procedures
meth2[c("bmi","hyp","chl")] <- "pmm"
# for countryID and fixed effects otherwise)
pred2[,"country"] <- 0 # don't use country factor
pred2[,"countryID"] <- 0 # don't use countryID
pred2[,colnames(DIs)] <- 1 # use dummy indicators
pred2["bmi", c("age","hyp","chl")] <- c(1,1,1) # fixed effects (bmi)
pred2["hyp", c("age","bmi","chl")] <- c(1,1,1) # fixed effects (hyp)
pred2["chl", c("age","bmi","hyp")] <- c(1,1,1) # fixed effects (chl)
# impute
imp2 <- mice(nhanes, maxit=20, m=10, predictorMatrix=pred2, method=meth2)
所以我使用 mice
包来估算缺失数据。我是插补的新手,所以我已经到了一定程度,但是 运行 进入了陡峭的学习曲线。举个玩具例子:
library(mice)
# Using nhanes dataset as example
df1 <- mice(nhanes, m=10)
如您所见,我使用大部分默认设置对 df1 进行了 10 次估算 - 我很乐意在回归模型、合并结果等中使用此结果。但是在我的现实生活数据中,我有来自不同国家/地区的调查数据。因此,缺失水平因国家/地区而异,特定变量的值(即年龄、教育水平等)也不同。因此,我想估算缺失值,允许按国家/地区进行聚类。所以我将创建一个没有缺失的分组变量(当然在这个玩具示例中缺少与其他变量的相关性,但在我的真实数据中它们存在)
# Create a grouping variable
nhanes$country <- sample(c("A", "B"), size=nrow(nhanes), replace=TRUE)
那么我如何告诉 mice()
这个变量与其他变量不同 - 即它是多级数据集中的一个级别?
您必须设置一个 predictorMatrix 来告诉老鼠使用哪个变量来估算另一个变量。一个快速的方法是使用 predictorM<-quickpred(nhanes)
然后将矩阵中的 1s 更改为 2,如果它是一个普通变量,如果它是不同国家的二级变量,则将其更改为 -2,并将其作为 predictorMatrix =predictorM
提交给 mice 命令。在方法命令中,您现在必须将方法设置为 2l.norm
(如果它是度量变量)或 2l.binom
(如果它是二进制变量)。对于后者,您需要 Sabine Zinn (https://www.neps-data.de/Portals/0/Working%20Papers/WP_XXXI.pdf) 编写的函数。不幸的是,我不知道世界上是否存在用于估算两级计数数据的方法。
请注意,对多级数据集进行估算会大大减慢该过程。根据我的经验,像 PMM 或 Baboon 包中的重采样方法在保持数据的层次结构方面效果很好,并且使用起来要快得多。
如果您将集群视为 "mixed-effects" 模型,那么您应该使用 mice
提供的用于集群数据的方法。这些方法可以在 manual 中找到,通常以 2l.something
.
聚类数据方法的多样性在 mice
中有所限制,但我可以推荐使用 2l.pan
来处理较低级别单元中的缺失数据,并在聚类级别使用 2l.only.norm
.
作为混合效应模型的替代方案,您可以考虑使用虚拟指标来表示集群结构(即,每个集群一个虚拟变量)。当您从混合效应模型的角度考虑集群时,此方法并不理想。因此,如果您想进行混合效应分析,请尽可能坚持使用混合效应模型。
下面,我展示了这两种策略的示例。
准备:
library(mice)
data(nhanes)
set.seed(123)
nhanes <- within(nhanes,{
country <- factor(sample(LETTERS[1:10], size=nrow(nhanes), replace=TRUE))
countryID <- as.numeric(country)
})
案例 1:使用混合效应模型进行插补
本节使用 2l.pan
来估算具有缺失数据的三个变量。请注意,我通过在预测矩阵中指定 -2
来使用 clusterID
作为聚类变量。对于所有其他变量,我仅分配固定效应 (1
)。
# "empty" imputation as a template
imp0 <- mice(nhanes, maxit=0)
pred1 <- imp0$predictorMatrix
meth1 <- imp0$method
# set imputation procedures
meth1[c("bmi","hyp","chl")] <- "2l.pan"
# set predictor Matrix (mixed-effects models with random intercept
# for countryID and fixed effects otherwise)
pred1[,"country"] <- 0 # don't use country factor
pred1[,"countryID"] <- -2 # use countryID as cluster variable
pred1["bmi", c("age","hyp","chl")] <- c(1,1,1) # fixed effects (bmi)
pred1["hyp", c("age","bmi","chl")] <- c(1,1,1) # fixed effects (hyp)
pred1["chl", c("age","bmi","hyp")] <- c(1,1,1) # fixed effects (chl)
# impute
imp1 <- mice(nhanes, maxit=20, m=10, predictorMatrix=pred1, method=meth1)
案例 2:使用集群的虚拟指标 (DI) 进行插补
本节使用pmm
插补,聚类结构用"ad hoc"表示。也就是说,聚类不是由随机效应表示,而是由固定效应表示。这可能会夸大具有缺失数据的变量的集群级可变性,因此请确保您在使用它时知道自己在做什么。
# create dummy indicator variables
DIs <- with(nhanes, contrasts(country)[country,])
colnames(DIs) <- paste0("country",colnames(DIs))
nhanes <- cbind(nhanes,DIs)
# "empty" imputation as a template
imp0 <- mice(nhanes, maxit=0)
pred2 <- imp0$predictorMatrix
meth2 <- imp0$method
# set imputation procedures
meth2[c("bmi","hyp","chl")] <- "pmm"
# for countryID and fixed effects otherwise)
pred2[,"country"] <- 0 # don't use country factor
pred2[,"countryID"] <- 0 # don't use countryID
pred2[,colnames(DIs)] <- 1 # use dummy indicators
pred2["bmi", c("age","hyp","chl")] <- c(1,1,1) # fixed effects (bmi)
pred2["hyp", c("age","bmi","chl")] <- c(1,1,1) # fixed effects (hyp)
pred2["chl", c("age","bmi","hyp")] <- c(1,1,1) # fixed effects (chl)
# impute
imp2 <- mice(nhanes, maxit=20, m=10, predictorMatrix=pred2, method=meth2)