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 함수를 사용할 것을 권한다. 자세한 설명은 추후에. 일단은 아래 링크를 참고하라.
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)