在R中解释嵌套混合效果模型输出

Question

我在解释嵌套混合效果模型中的基准系数时遇到问题。我已经安装了一个模型Test.Score〜Subject +（1 | School/Class），因为班级嵌套在学校内。当我看着系数但是使用COEF（模型），他们似乎直觉：

$`Class:School` 
     (Intercept) SubjectMaths 
1:A 82.73262 -4.108333 
1:B 83.98870 -4.108333 
1:C 82.26456 -4.108333 
2:A 82.25383 -4.108333 
2:B 78.22047 -4.108333 
2:C 80.18982 -4.108333 

$School 
(Intercept) SubjectMaths 
A 88.39636 -4.108333 
B 77.74404 -4.108333 
C 78.68460 -4.108333 

attr(,"class") 
[1] "coef.mer" 

学校内班级怎能舍弃值比那些刚内学校低得多？数据下面复制：

Test.Score = c(94,88,86,90,94,87,87,92,89,92,87,94,93,91,89,92,91, 
    91,95,91,82,84,90,81,92,89,85,94,88,94,94,94,86,94,93,84,82, 
    92,92,83,89,83,81,87,84,80,81,83,88,82,81,90,82,85,87,82,86, 
    84,87,88,82,91,95,77,88,87,79,75,91,77,82,91,95,92,89,83,79, 
    90,83,83,82,79,79,78,83,82,81,77,80,79,84,83,81,78,77,75,76, 
    76,84,75,78,78,71,79,70,75,75,78,76,71,76,76,73,71,80,70,71, 
    78,71,74,76,74,74,77,81,78,79,76,82,79,80,73,72,83,72,81,81, 
    72,79,74,67,75,71,66,65,71,73,69,65,67,71,72,68,73,65,65,74, 
    67,72,72,82,70,72,86,89,87,87,88,74,92,70,89,86,63,68,74,88, 
    71,88,91,76,86,75,79,76,69,86,71,78,67,67,73,69,81,79,78,80, 
    72,81,69,72,75,76,68,72,78,78,77,71,73,70,77,75,75,69,77,74, 
    76,68,78,76,75,68,74,69,78,76,70,79,78,67,65,86,88,65,88,73, 
    66,65,85) 
School = rep(c("A","B","C"), each = 80) 
Class = rep(c("1","2"), each = 20,6) 
Subject = rep(c("English","Maths"), each = 40, 3) 
data = data.frame(Test.Score, School, Class, Subject) 
data$Class = factor(data$Class) 
mod = lmer(Test.Score ~ Subject + (1|School/Class), REML = F, 
    data = data) 
coef(mod) 

Answer 1

的问题是，每个随机效应所列的系数包括只的特定随机效应的影响。特别是，2级的系数只反映学校内班级与总体平均值的偏差 - 而不是学校的影响。这可能看起来很奇怪或错误，但（1）你可以通过predict()（见下面）和（2）lme4得到你想要的内容，但它并不真正具有“嵌套”的内部表示，因此很难一般确定哪些随机效应应该包含在给定的一组系数中（这是一个解释，而不是借口）。

对于它的价值，coef()不工作，你期待装有nlme::lme机型...

library(lme4) 
## using sum-to-zero contrasts for convenience 
mod = lmer(Test.Score ~ Subject + (1|School/Class), REML = FALSE, 
      data = data, contrasts=list(Subject=contr.sum)) 
pframe <- with(data,expand.grid(School=levels(School), 
          Subject=levels(Subject), 
          Class=levels(Class))) 
pframe$Test.Score <- predict(mod,newdata=pframe) 

如果你想要的平均等级值，你需要平均的英语和数学成绩。 ..

在nlme::lme同样的模式：

mod2 = nlme::lme(Test.Score ~ Subject, random = ~ 1|School/Class, method="ML", 
    data = data, contrasts=list(Subject=contr.sum)) 
coef(mod2)   ## Class within School 
coef(mod2,level=1) ## School-level 

对于一些单调乏味（和tidyverse工具 - 这可以b E在其他方面也做了），重新排列系数策划：

rr2 <- tibble::rownames_to_column(coef(mod)[["Class:School"]]) 
rr2 <- dplyr::rename(rr2,Test.Score=`(Intercept)`) 
rr2 <- tidyr::separate(rowname,data=rr2,into=c("Class","School")) 
rr2$Subject <- NA 

rr3 <- tibble::rownames_to_column(coef(mod)[["School"]]) 
rr3 <- dplyr::rename(rr3,Test.Score=`(Intercept)`,School=rowname) 
rr3$Subject <- NA 
rr3$Class <- 1.5 

情节都在一起（数据，预测系数）：

library(ggplot2); theme_set(theme_bw()) 
ggplot(data,aes(Class,Test.Score,colour=Subject))+ 
    geom_boxplot()+ 
    geom_point(data=pframe,size=3,shape=16,position=position_dodge(width=0.75))+ 
    facet_wrap(~School,labeller=label_both)+ 
    geom_point(data=rr2,size=3,shape=17)+ 
    geom_hline(yintercept=fixef(mod)["(Intercept)"],lty=2)+ 
    geom_point(data=rr3,size=5,shape=18)+ 
    theme(panel.spacing=grid::unit(0,"lines")) ## cosmetic 

彩色点预测;灰色点是系数（三角形=班级;钻石=学校级别）。