Describe 레만투_혼합효과모형 here

간략한 이론
혼합 효과 모형(Mixed Effect Model)은 고정 효과(fixed effect)와 무선 효과(random effect)을 이용하여 종속변수를 설명한다. 아래에서 A는 고정 효과, B는 무선 효과이다.

Y = AX + BZ + e

HLM 또는 다층 모형(Multilevel Model)은 집단과 개인, 개인과 시행처럼 위계적인 형태를 띄고 있는 자료를 분석하기 위한 모형이다. 아래의 모형을 보면 절편(β0)과 회귀계수(β1)가 또다른 회귀식의 형태를 띄고 있다.

Y = β0 + β1X + e
β0 = γ00 + γ01Z + u0
β1 = γ10 + γ11Z + u1

위의 HLM은 혼합 효과 모형으로 다시 나타낼 수 있다. 맨 위의 회귀식에 아래의 두 회귀식을 대입하면 아래와 같다.

Y = (γ00 + γ01Z + u0) + (γ10 + γ11Z + u1)X + e

이것을 전개하면 아래와 같다.

Y = γ00 + γ10X + γ01Z + γ11ZX + u0 + u1X + e

이를 다시 혼합 모형의 표기법으로 바꾸면 아래와 같다.

Y = a0 + a1X + a2Z + a3ZX + b0 + b1X + e

lme4
아래에서 소개하는 nlme 패키지의 lme 함수보다 lme4 패키지의 lmer 함수를 사용할 것을 권한다. 자세한 설명은 추후에. 일단은 아래 링크를 참고하라.

참고: http://www.stat.columbia.edu/~cook/movabletype/archives/2006/01/fitting_multile.html

nlme
nlme 패키지는 비선형 혼합 효과 모형(nonlinear mixed effect model)을 지원한다. 여기서는 nlme를 이용하여 위계적 선형 모형((Hierarchical Linear Model:HLM)을 다루는 방법을 알아보겠다. 아래는 HLM을 선형 혼합 효과 모형으로 해석하는 방법을 간략히 소개한 것이다.

예제
먼저 nlme 패키지를 풀러들인다.

> library(nlme)

자료 만들기
nlme 패키지에 포함된 샘플 자료를 불러들인다. MathAchieve는 학생 자료이고, MathAchSchool은 학교 자료이다. Bryk와 Rodenbush(1992)의 HLM 책에 있는 예제 자료이다.

> data(MathAchieve)
> data(MathAchSchool)

Bryk와 Rodenbush의 HLM 프로그램과 달리 nlme에서는 각 수준의 자료를 따로 만들지 않는다. 먼저 분석할 자료를 만들어보자. 우선 학생 자료에서 학교 번호, 사회경제적 수준(socioeconomical status:SES), 수학성적을 뽑아서 HLM 분석을 위한 데이터프레임으로 만든다. 이름은 Bryk라고 하자.

> Bryk = as.data.frame(MathAchieve[c("School","SES","MathAch")])

tapply 함수를 이용하여 사회경제적 수준의 학교별 평균을 계산한다.

> attach(MathAchieve)
> mses = tapply(SES, School, mean)
> detach(MathAchieve)

Bryk 데이터 프레임에 학교별 평균을 덧붙여준다. 아래에서 첫째줄은 학생마다 그 학생이 속한 학교의 평균 SES를 덧붙여준다. 그리고 둘째줄은 각 학교의 평균 SES에 대한 학생의 SES의 편차를 구한다.

> Bryk$mses = mses[as.character(Bryk$School)]
> Bryk$cses = Bryk$SES - Bryk$mses

비슷한 방법으로 학교 자료에서 학교의 종류(공립학교, 카톨릭계 학교)를 뽑아서 Bryk 데이터프레임에 덧붙인다.

> sector = MathAchSchool$Sector
> names(sector) = row.names(MathAchSchool)
> Bryk$sector = sector[as.character(Bryk$School)]

학교의 분류는 이산변수이므로 factor로 바꿔준다.

> Bryk$sector = factor(Bryk$sector, levels = c("Public", "Catholic"))

contrast coding을 하기 위해 contrasts 함수를 사용한다. 이제 회귀식에서 Bryk$sector는 Public은 0, Catholic은 1로 코딩한다.

> contrasts(Bryk$sector)
         Catholic
Public          0
Catholic        1

분석
Bryk 데이터 프레임에 필요한 자료는 다 들어갔으니 분석을 실시한다. 분석할 모형은 아래와 같다. 간단히 말하면 학교의 평균적인 SES와 학교의 종류에 따라서 학생 SES와 수학 성적의 관계가 변화하는 것을 나타내는 모형이다.

(MathAch) = β0 + β1(cse) + e
β0 = γ00 + γ01(mses) + γ02(sector) + u0
β1 = γ10 + γ11(mses) + γ12(sector) + u1

위의 모형을 혼합 효과 모형으로 바꾸면 아래와 같다.

(MathAch) = γ00 + γ01(mses) + γ02(sector) + u0 + (γ10 + γ11(mses) + γ12(sector) + u1)(cses) + e
= γ00 + γ01(mses) + γ02(sector) + γ10(cses) + γ11(mses)(cses) + γ12(sector)(cses) + u1(cses) + u0 + e

위의 모형에서 고정 효과 부분만 R의 관계식으로 다시 쓰면 아래와 같다.

MathAch ~ 1 + mses + sector + cses + mses:cses + sector:cses

절편은 명시적으로 써주지 않아도 기본으로 포함되므로 빼고, 상호작용항은 줄여서 다시 쓰면 아래처럼 된다.

MathAch ~ mses*cses + sector*cses

무선 효과는 R의 관계식으로 아래처럼 된다. "| School"이 붙는 이유는 무선 효과가 학교 수준이기 때문이다.

~ 1 + cses | School

선형 혼합 효과 모형 분석은 lme 함수를 이용한다.

> bl = lme(MathAch ~ mses*cses + sector*cses, random=~cses|School, data=Bryk)

해석
다른 분석과 마찬가지로 summary 함수를 이용하면 분석 결과를 볼 수 있다.

> summary(bl)
Linear mixed-effects model fit by REML
 Data: Bryk 
       AIC      BIC    logLik
  46523.66 46592.45 -23251.83

Random effects:
 Formula: ~cses | School
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.5426150 (Intr)
cses        0.3182015 0.391 
Residual    6.0597955       

Fixed effects: MathAch ~ mses * cses + sector * cses 
                        Value Std.Error   DF  t-value p-value
(Intercept)         12.127931 0.1992919 7022 60.85510   0e+00
mses                 5.332875 0.3691684  157 14.44564   0e+00
cses                 2.945041 0.1556005 7022 18.92694   0e+00
sectorCatholic       1.226579 0.3062733  157  4.00485   1e-04
mses:cses            1.039230 0.2988971 7022  3.47688   5e-04
cses:sectorCatholic -1.642674 0.2397800 7022 -6.85076   0e+00
 Correlation: 
                    (Intr) mses   cses   sctrCt mss:cs
mses                 0.256                            
cses                 0.075  0.019                     
sectorCatholic      -0.699 -0.356 -0.053              
mses:cses            0.019  0.074  0.293 -0.026       
cses:sectorCatholic -0.052 -0.027 -0.696  0.077 -0.351

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.1592608 -0.7231893  0.0170471  0.7544510  2.9582205 

Number of Observations: 7185
Number of Groups: 160 

위의 결과로 회귀식을 다시 쓰면 아래처럼 된다.

(MathAch) = β0 + β1(cse) + e
β0 = 12.13 + 5.33(mses) + 1.23(sector) + u0
β1 = 2.95 + 1.04(mses) - 1.64(sector) + u1

u0 ~ N(0, 1.54)
u1 ~ N(0, 0.32)
e ~ N(0, 6.06)