使用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

如果我正确理解了这篇文章，在预测矩阵中，我应该使用-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将其标识为“类”变量，但在这个mice入门指南中，建议对于多层级模型，您应该在数据框中创建一个全为1的变量作为常数，然后通过预测矩阵中的2指定为随机截距。然而，这是基于2l.norm函数而不是2l.pan函数，所以我真的不确定我在这里的位置。 2l.pan函数是否需要此列或随机效应的指定?
2. 是否有任何方法可以指定纵向随机斜率模型，如果有，如何指定？

- llewmills

Answer 1

这个答案可能对你来说有点晚了，但它可能会帮助一些将来阅读此内容的人：

如何使用2l.pan 下面是关于使用mice指定多级插补模型的一些细节。由于应用程序是纵向的，因此我使用术语“人员”来指代第二级别的单位。这些是mice文档中提到的与2l.pan相关的最相关参数：

type

向量，长度为ncol(x)，用于标识随机和类变量。 随机效应由2标识。组变量（仅允许一个）编码为-2。 随机效应还包括固定效应。如果要计算协变量X1的组均值并将其作为进一步的固定效应包括在内，则选择3。除了3中的效应之外，规范4还包括X1的随机效应。

在使用2l.pan进行变量插补时，您可以在预测矩阵中使用5种不同的代码。人员标识符编码为-2（这与2l.norm不同）。要包括具有固定或随机效应的预测变量，这些变量分别使用1或2进行编码。如果编码为2，则相应的固定效应会自动包括在内。

此外，2l.pan还提供了代码3和4，其含义与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 进行更改），并且不需要在预测变量矩阵中指定它。

2. 是的，你可以使用 2l.pan 指定随机斜率。为此，随机斜率的预测变量在预测变量矩阵中编码为 2 或 4。如果编码为 2，则随机斜率被包括在内。如果编码为 4，则随机斜率以及该变量的个人平均值的另一个固定效应也被包括在内。如果编码为 4，则可以利用 groupcenter.slope=TRUE 来改变随机斜率的含义（请参见上文）。

这篇文章还包括一些使用 2l.pan 和其他多层次插补函数的实例：[链接]

- SimonG

---

@SimonG，你的回答太棒了，让我开心了一整天。虽然前面的回答也不错，但是考虑到你的回答更加全面和准确，我必须把采纳答案给你。 - llewmills

Answer 2

pan库不需要截距项。

您可以使用以下函数进行详细了解。

library(pan)
?pan

话虽如此，mice使用围绕pan（即mice.impute.2l.pan）的包装器，如果已加载mice库，则可以查看该函数的帮助文档。它声明：它有一个名为intercept的参数，这是[a] Logical [and] determin[es] whether the intercept is automatically added.默认情况下为TRUE。这被默认定义为随机截距。在浏览mice包装器的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的示例，用于固定效应。

对于2l.norm，它还指出：The random intercept is automatically added in mice.impute.2l.norm().

它有一些带有描述的示例。 CRAN文档中的pan可能会对您有所帮助。

- Matt L.

---

谢谢@Matt L. 那么据你所知，没有随机斜率插补方法吗？ - llewmills

在 mice.impute.2l.norm() 的帮助文档中，有一段不寻常的内容：当前实现的算法无法处理指定为固定效应（type=1）的预测变量。使用 mice.impute.2l.norm() 时，当前的建议是将所有预测变量都指定为随机效应（type=2）。 因此，我的模型中所有的固定效应都应该被指定为随机效应，而传统的随机效应 id 则被指定为类变量。我必须承认我不太理解，但感谢您指引我方向。 - llewmills