ネストされた混合エフェクトモデルの出力を解釈する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 School:Classの係数は、全母集団平均から学校内のクラスの偏差のみを反映している - ではなく、の学校レベルの効果も反映されている。それは奇妙で間違っているかもしれませんが、（1）predict()（下記参照）と（2）lme4であなたが探しているものを得ることは本当に "ネスティング"の内部表現を持っていないので、与えられた係数セットにどのランダムな効果が含まれるべきかを一般的に決定する（これは説明であり、言い訳ではない）。何が価値があるために

、あなたはnlme::lmeを装備したモデルのために期待するようcoef()が仕事をして...

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) 

あなたは平均的なクラスの値を望んでいた場合は、英語と数学のスコアを平均化する必要があると思います。いくつかの退屈（およびtidyverseツールを使用

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 

- これはbの可能性..

nlme::lmeで同じモデルを同様に他の方法で行う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 

カラーポイントが予測されています。グレーポイントは係数（三角=クラスレベル、ダイヤモンド=学校レベル）です。