如何修改R代码以复现Raudenbush&Bryk的Intercepts- and Slopes-as Outcomes模型正确随机效应？

我正在尝试用R复现Raudenbush和Bryk所著《Hierarchical Linear Models》中，基于High School and Beyond数据集的Intercepts- and Slopes-as Outcomes模型结果。使用的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)

代码运行后，结果与书中及HLM6复现的正确结果相近，但随机效应部分存在明显差异。

请问需要对代码做哪些修改，或是采用何种正确方式才能在R中获取与书中一致的准确结果？

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