R 中截距和斜率作为结果模型的随机效应

Question

我正在尝试使用统计软件 R 重现 Raudenbush 和 Bryk 所著的《分层线性模型》书中的截距和斜率结果模型的结果，该书涉及高中及以后的数据集。

hsb

数据集来自mlmhelpr包，我正在尝试的R代码是：

library(lme4)

library(mlmhelpr)

hsb=hsb

groupMeanses = aggregate(hsb$ses, list(hsb$id), FUN = mean, data=hsb)

names(groupMeanses) = c('id','groupMeanses')

hsb = merge(hsb, groupMeanses, by = c('id'))

hsb$ses_centered_a = (hsb$ses - hsb$groupMeanses)

hsb$ses_centered_b = (hsb$ses - hsb$meanses)

M3d = lmer(mathach ~ 1 + ses_centered_b * meanses + ses_centered_b * catholic 
      + (1+meanses|id) , REML = F, data=hsb, 
      control = lmerControl(optimizer = "bobyqa", boundary.tol = 1e-5))

summary(M3d)

及其结果：

如您所见，这两个结果（Raudenbush 和 Bryk 书中的结果以及我的 R 代码看起来非常接近，除了随机效应之外。我在互联网上找到了一个教程，该教程使用 HLM 6 适合相同的模型，可以正确再现结果，所以我一定是在 R 中做错了什么。

我的问题是

：我应该在代码中进行哪些修改，或者在 R 中获得正确结果的正确方法是什么？ 如有任何帮助，我们将不胜感激。

Answer 1

更详细一点：我用 lme4::lmer、

nlme::lme

和

glmmTMB::glmmTMB

拟合了相同的模型，它们都以（略有）不同的方式实现了混合模型拟合。他们都发现了“奇异拟合”，斜率标准差很小或最小（从 6e-6 到 0.05），斜率和截距之间的相关性高度可变（这是有道理的——当斜率 SD 趋于零时，斜率-截距相关性变为不明确的）。截距 SD 和残差 SD 的偏差和估计高度一致。

  package sd.(Intercept)           cor   sd.meanses sd.Observation deviance
1    lmer       1.520733 -1.0000000000 5.341325e-02       6.062259 46497.41
2     lme       1.520815 -0.0001895723 2.457569e-03       6.062254 46497.44
3 glmmTMB       1.520813  0.9723940794 5.643373e-06       6.062255 46497.44

我还使用 allFit() 函数在

lmer

 内尝试各种优化器 - 这些都给出了一致的答案。

下一个测试是从 HLM6 中进行估计，评估各种 R 包的对数似然（这对于

lmer

和

glmmTMB

来说相当简单），并看看我们是否恢复了 HLM6 偏差估计。不幸的是，HLM6 似乎没有报告截距斜率协方差（或相关性）。

library(lme4)
library(nlme)
library(glmmTMB)

M3d = lmer(mathach ~ 1 + ses_centered_b * meanses + ses_centered_b * catholic 
      + (1+meanses|id) , REML = FALSE, data=hsb, 
      control = lmerControl(optimizer = "bobyqa", boundary.tol = 1e-5))

M3e = glmmTMB(mathach ~ 1 + ses_centered_b * meanses + ses_centered_b * catholic 
      + (1+meanses|id) , REML = FALSE, data=hsb)

M3f = lme(mathach ~ 1 + ses_centered_b * meanses + ses_centered_b * catholic,
          random = ~ 1+meanses|id , method = "ML", data=hsb)


library(broom.mixed)
library(tidyverse)

fitlist <- list(lmer = M3d, lme = M3f, glmmTMB = M3e)

(fitlist
    |> purrr::map_dfr(tidy, effect = "ran_pars", .id = "package")
    |> select(package, term, estimate)
    |> mutate(across(term, ~gsub("_+",".", .)))
    |> mutate(across(term, ~gsub("cor.*", "cor", .)))
    |> pivot_wider(names_from = "term", values_from = "estimate")
    |> mutate(deviance = c(deviance(fitlist$lmer), deviance(fitlist$lme), 2*M3e$obj$fn()))
    |> as.data.frame()
)

aa <- allFit(M3d)

(tidy(aa, effect = "ran_pars")
    |> select(optimizer, term, estimate)
    |> pivot_wider(names_from = "term", values_from = "estimate")
    |> as.data.frame()
)