如何计算随机效应各分组边际效应？模拟数据斜率无变化排查

我正在为项目使用多水平模型，流程是模拟数据、拟合随机效应模型，再生成每个Subject分组的边际效应。虽然模型设置和sleepstudy示例完全一致，但拟合后的模型无法输出有差异的边际效应斜率。

sleepstudy示例（正常运行） #

以下代码可以正常生成不同斜率的分组曲线：

library(tidyverse)
library(lme4)
library(ggeffects)

data("sleepstudy")
m <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
me <- ggpredict(m, terms = c("Days", "Subject [sample=8]"), type = "random")
plot(me)

运行后，图中会显示8个Subject对应的8条斜率不同的曲线。

自定义模拟数据的问题 #

我模拟了和sleepstudy结构近似的数据，但拟合模型后无法生成变化的斜率：

N = 18 # 受试者数量
Nt = 9 # 测试次数

d <- tibble(
  Subject = factor(sprintf("%03d", 1:N)),       # 生成受试者编号
  subj_b0 = rnorm(n = N, mean = 250, sd = 20),  # 生成随机截距
  subj_b1 = rnorm(n = N, mean = 10, sd = 6)     # 生成随机斜率       
) %>%                                           
  mutate(Days = list(0:Nt)) %>%
  unnest(Days) %>%
  mutate(
    Y = subj_b0 + Days*subj_b1,
    Y = Y + rnorm(n = nrow(.), sd = 15)          # 添加误差项
  )    

fit <- lmer(Y ~ 1 + Days + (1 + Days|Subject), data = d)
p <- ggpredict(fit, terms = c("Days", "Subject [sample=8]"), type = "random")
plot(p)

这段代码生成的图中，所有Subject的斜率完全一致。

模型摘要对比 #

模型摘要确认Days的斜率存在变异，且其他估计参数和sleepstudy近似，但ggpredict()返回的分组斜率始终一致：

自定义数据拟合模型的摘要 #

summary(fit)

Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ 1 + Days + (1 + Days | Subject)
   Data: d

REML criterion at convergence: 1553.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.57322 -0.57070 -0.01073  0.59463  3.02611 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 745.22   27.30        
          Days         37.95    6.16    0.12
 Residual             176.54   13.29        
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  254.517      6.693  38.030
Days           9.033      1.492   6.053

sleepstudy模型的摘要 #

summary(m)

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy

REML criterion at convergence: 1743.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4634  0.0231  0.4634  5.1793 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.10   24.741       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.825  36.838
Days          10.467      1.546   6.771

内容的提问来源于stack exchange，提问作者sometimes_sci