机器算法验证 - 在没有相应固定效应的 lmer 模型中包含随机斜率项是否合理？ - 吾爱随笔录

在没有相应固定效应的 lmer 模型中包含随机斜率项是否合理？

机器算法验证 r 混合模式随机效应模型 lme4-nlme

2022-03-31 09:38:42

我有一个实验，我向参与者展示了多种刺激，并希望控制刺激的显示顺序。我很好奇是否可以只考虑随机效应中的顺序而不包括固定效应中的项（即，因为在这种情况下顺序是一个令人讨厌的变量）。问题是下面的第二个公式是否有效：

mod1 <- lmer(Response ~ Order + (Order|Subject), data = dat)
mod2 <- lmer(Response ~ 1 + (Order|Subject), data = dat)

就可重现的示例而言，可以使用 sleepstudy 数据：

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fm2 <- lmer(Reaction ~ 1 + (Days | Subject), sleepstudy)

同样的问题也适用：如果我不希望主体内天数的随机斜率取决于天数的固定效应，那么 fm2 是否合理拟合。换句话说，对于 fm2，除了固定效应截距之外，随机效应将解释与主题相关的数据的所有方面。使用 sjPlot 绘制随机斜率时，考虑到预期结果，结果似乎是合理的：

library(sjPlot)
sjp.lmer(fm1, type = "rs.ri")

sjp.lmer(fm2, type = "rs.ri")

library(ggplot2)
ggplot(sleepstudy) + geom_line(aes(x = Days, y = Reaction, color = Subject))

1个回答

我相信这个问题与经常想知道的“必须始终在线性回归中包含截距项”非常相似，对此商定的答案是“是的，除非您有充分的理由不这样做”。

在进行任何实验之前，我试图考虑没有固定效应项会发生什么。让我们详细写出你的两个模型。第一个，具有固定效应斜率，是

y \sim N (μ_{α} + α_{[i]} + (μ_{β} + β_{[i]}) x, σ)

$y \sim N(\mu_{\alpha} +\alpha_{[i]} + (\mu_{\beta} + \beta_{[i]}) x, \sigma)$

α \sim N (0, σ_{α})

$\alpha \sim N(0, \sigma_{\alpha})$

β \sim N (0, σ_{β})

$\beta \sim N(0, \sigma_{\beta})$

其中是天数，我们有一个随机截距和每个主题的随机斜率。在另一种情况下，没有固定斜率，模型是 $x$ $\alpha_{[i]}$ $\beta_{[i]}$

y \sim N (μ_{α} + α_{[i]} + β_{[i]} x, σ)

$y \sim N(\mu_{\alpha} + \alpha_{[i]} + \beta_{[i]} x, \sigma)$

α \sim N (0, σ_{α})

$\alpha \sim N(0, \sigma_{\alpha})$

β \sim N (0, σ_{β})

$\beta \sim N(0, \sigma_{\beta})$

不同之处在于，在第二个模型中，我们先验地假设随机斜率的平均值为零。这意味着，我们期望与各个主题相关的斜率均匀分布在斜率附近（例如，一半应该是负数，一半应该是正数）。 $0$

现在，在您的数据模型中，这似乎不是真的。在您的第二个图中，每个主题内的估计斜率都是正的。看起来此模型对您的数据无效。包含固定斜率包括作为自由度的主体斜率的平均值，在此图中，您可以看到随机斜率均匀地聚集在零附近，如您所愿。

至于从模型中的参数推断，我相信模型的这种错误陈述会导致以下参数估计存在偏差

主体斜率将偏向零，因为可能性中均值为零的假设会将它们拉向零。
随机斜率的估计标准偏差会太大，因为膨胀这个参数会让斜率聚集在它们的真实、非零均值周围，而不会受到如此严重的惩罚。

在这里，我将创建一些模拟数据，其中真实的主题平均斜率非零

library("lme4")
library("arm")
set.seed(154)

N_classes = 50
N_obs <- 10000

random_intercepts <- structure(
  rnorm(N_classes), names = as.character(1:N_classes)
)

random_slopes <- structure(
  rnorm(N_classes, mean = 1), names = as.character(1:N_classes)
)

classes <- sample(as.character(1:N_classes), size = N_obs, replace = TRUE)
x <- runif(N_obs)
y <- random_intercepts[classes] + random_slopes[classes] * x + rnorm(N_obs)

df <- data.frame(class = factor(classes), x = x, y = y)

第一个模型很好地估计了所有真实参数

> M <- lmer(y ~ x + (x | class), data = df)
> display(M)
lmer(formula = y ~ x + (x | class), data = df)
        coef.est coef.se
(Intercept) 0.01     0.15   
x           1.02     0.15   

Error terms:
 Groups   Name        Std.Dev. Corr 
 class    (Intercept) 1.03          
          x           1.01     0.19 
 Residual             1.00

看起来这里所有的参数都被很好地估计了，包括随机斜率的标准偏差。

这是没有固定斜率的模型

> N <- lmer(y ~ (x | class), data = df)
> display(N)
lmer(formula = y ~ (x | class), data = df)
coef.est  coef.se 
   -0.14     0.15 

Error terms:
 Groups   Name        Std.Dev. Corr 
 class    (Intercept) 1.04          
          x           1.43     0.24 
 Residual             1.00

随机斜率标准差的估计值为1.43，证实了我的直觉，即它会偏大。

模型中主体斜率的平均值M很好

> mean(fixef(M)["x"] + ranef(M)$class$x)
[1] 1.015418

在另一个模型上，我的直觉似乎不太正确

> mean(ranef(N)$class$x)
[1] 0.9858566

看起来模型比确保完全满足随机斜率假设的正态性更认真地拟合数据。总而言之，随机斜率标准差的膨胀似乎是最严重的问题。

其它你可能感兴趣的问题

上一篇中间奖励可以用于强化学习吗？下一篇聚类和 A/B 测试