중첩 된 혼합 효과 모델 출력을 해석하는 중 R

Question

중첩 혼합 효과 모델에서 기준 계수를 해석하는 데 문제가 있습니다. Class.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) 

작업을 수행합니다. ..

nlme::lme 동일 모델 : 약간의 지루함 (및 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 

-이 ㄱ 수 뿐만 아니라 다른 방법으로 수행 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 

컬러 포인트가 예측이다; 회색 점은 계수입니다 (삼각형 = 클래스 수준, 다이아몬드 = 학교 수준).