使用 MICE 的纵向多级插补模型中的随机效应

Question

我正在尝试使用纵向设计来估算数据集中的数据。有两个预测变量（实验组和时间）和一个结果变量（分数）。聚类变量为 id。

这是玩具数据

set.seed(345)
A0 <- rnorm(4,2,.5)
B0 <- rnorm(4,2+3,.5)
A1 <- rnorm(4,6,.5)
B1 <- rnorm(4,6+2,.5)
A2 <- rnorm(4,10,.5)
B2 <- rnorm(4,10+1,.5)
A3 <- rnorm(4,14,.5)
B3 <- rnorm(4,14+0,.5)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- rep(1:8,times = 4, length = 32)
time <- rep(0:3, each = 8, length = 32)
group <- rep(c("A","B"), times =2, each = 4, length = 32)
df <- data.frame(id = id, group = group, time = time,  score = score)

# plots
(ggplot(df, aes(x = time, y = score, group = group)) + 
    stat_summary(fun.y = "mean", geom = "line", aes(linetype = group)) +
    stat_summary(fun.y = "mean", geom = "point", aes(shape = group), size = 3) +
    coord_cartesian(ylim = c(0,18)))

# now place some NAs
df[sample(1:nrow(df), 10, replace = F),"score"] <- NA

df

如果我正确理解this post，在预测矩阵中我应该用-2 指定id 聚类变量，用1 指定两个固定预测变量time 和group。像这样

library(mice)

(ini <- mice(df, maxit=0))
(pred <- ini$predictorMatrix)
(pred["score",] <- c(-2, 1, 1, 0))
(imp <- mice(df, 
            method = c("", "", "", "2l.pan"),
            pred = pred, 
            maxit = 1, 
            seed = 71152))

我想知道的是：

1. 这是纵向随机截距插补模型吗？将 id 变量指定为 -2 将其指定为“类”变量，但在 this mice primer 中，它建议对于多级模型，您应该在数据框中创建一个包含所有 1 的变量作为常量，然后通过预测矩阵中的2 指定为随机截距。但是，这是基于2l.norm 函数而不是2l.pan 函数，所以我不确定我在哪里。 2l.pan 函数是否不需要此列，或者随机效应的规范？
2. 有什么方法可以指定纵向随机斜率模型，如果可以，如何指定？

Answer 1

这个答案对你来说可能有点晚了，但它可能能够帮助一些将来阅读这篇文章的人：

如何使用2l.pan

以下是有关使用mice 指定多级插补模型的一些详细信息。因为应用程序是纵向的，所以我使用术语“人员”来指代 2 级的单位。这些是 2l.pan 最相关的参数，如 mice 文档中所述：

type

长度为ncol(x) 的向量标识随机变量和类变量。 随机效应由2 标识。组变量（只有一个 是允许的）编码为-2。随机效应还包括固定 影响。如果对于协变量 X1 组均值应计算和 包括作为进一步的固定效果选择3。除了 3 中的效果，规范 4 还包括随机效果 X1.

对于用2l.pan 估算的变量，您可以在预测矩阵中使用 5 种不同的代码。人员标识符编码为-2（这与2l.norm 不同）。为了包含具有固定或随机效应的预测变量，这些变量分别用1 或2 编码。如果编码为2，则自动包含对应的固定效果。

此外，2l.pan 提供代码3 和4，它们与1 和1 和2 具有相似的含义，但将包含对该变量的人均值的附加固定效应。如果您尝试对随时间变化的预测变量的个人内部和个人之间的影响进行建模，这将非常有用。

intercept

判断是否自动添加拦截的逻辑。

默认情况下，2l.pan 包含截距作为固定和随机效果。因此，不需要在预测矩阵中包含常数项。如果设置了intercept=FALSE，则会更改此行为，并将截距从插补模型中删除。

groupcenter.slope

如果TRUE，在组均值的情况下（type 是3 或4）组均值 在进行插补之前对这些预测变量进行居中。 默认为FALSE。

使用此选项，可以将预测变量集中在人的平均值周围，而不是“按原样”包括预测变量（即，不居中）。这仅适用于编码为3 或4 的变量。对于编码为3 的预测变量，这不是很重要，因为有和没有居中的模型是相同的。

但是，当预测变量编码为4（即具有随机斜率）时，居中会改变随机效应的含义，因此随机斜率不再“按原样”应用于变量，而是应用于该变量的人内偏差。

---

在您的示例中，您可以为time 添加一个简单的随机斜率，如下所示：

library(mice)
ini <- mice(df, maxit=0)

# predictor matrix (following 'type')
pred <- ini$predictorMatrix
pred["score",] <- c(-2, 1, 2, 0)

# imputation method
meth <- c("", "", "", "2l.pan")

imp <- mice(df, method=meth, pred=pred, maxit=10, m=10)

在本例中，将time 编码为3 或4 没有多大意义，因为time 的人员均值对于所有人来说都是相同的。但是，如果您希望将随时间变化的协变量作为预测变量包含在插补模型中，3 和 4 会很有用。

intercept 和groupcenter.slope 等附加参数可以直接在对mice() 的调用中指定，例如：

imp <- mice(df, ..., groupcenter.slope=TRUE)

---

关于您的问题

所以，按照帖子中的说明回答您的问题：

1. 是的，2l.pan 提供了一个多级（或者更确切地说是两级）插补模型。默认情况下，截距既作为固定效应也作为随机效应包含（可以使用intercept=FALSE 更改），并且不需要在预测矩阵中指定（这与2l.norm 形成对比）。

2. 是的，您可以使用2l.pan 指定随机斜率。为此，具有随机斜率的预测变量在预测矩阵中编码为2 或4。如果编码 如2，包含随机斜率。如果编码为4，则包括随机斜率以及该变量的个人均值的附加固定效应。如果编码为4，则可以通过使用groupcenter.slope=TRUE来改变随机斜率的含义（见上文）。

本文还包括一些关于如何使用 2l.pan 和其他函数进行多重插补的工作示例：[Link]

Answer 2

pan 库不需要拦截词。

你可以使用

library(pan)
?pan

也就是说mice 使用了一个名为mice.impute.2l.pan 的pan 包装器并加载了mice 库，您可以查看该函数的帮助。它声明：它有一个名为intercept 的参数，即[a] Logical [and] determin[es] whether the intercept is automatically added. 默认情况下为TRUE。默认情况下，这被定义为随机截距。在浏览了鼠标包装器的 R 代码后发现了这一点：

if (intercept) { x <- cbind(1, as.matrix(x)) type <- c(2, type) }

其中pan 函数参数type 是Vector of length ncol(x) identifying random and class variables。默认情况下添加截距并定义为随机效果。

他们确实提供了一个例子，就像你所说的那样，在预测矩阵中用 1 表示“x”，以获得固定效果。

它还声明2l.norm、The random intercept is automatically added in mice.impute.2l.norm().

它有一些带有描述的示例。 The CRAN documentation for pan might help you.