Linear mixed effect model - 이번엔 제발 이해해 보자...

통계를 돌리다가 LMM을 쓸 일이 있는데

 

아직도 정확히 이해를 못하고 있어서 웹을 서칭해보다가 좋은 자료를 발견했다. 

 

https://bodowinter.com/tutorial/bw_LME_tutorial1.pdf

 

Bodo Winter라는 분이 10년 전쯤 작성한, R에서의 LME tutorial이고 2개의 pdf 파일로 이뤄져 있는데, 하나 하나 따라가며 이해해보도록 하자. 

위와 같은 자료를 보면 female이 male보다 voice pitch가 높은 것은 자명하다고 볼 수도 있으나

실제로는 비슷한데 sampling 때문에 저렇게 보일 수도 있다. 

즉, 여성이 남성보다 유의미하게 pitch가 높을 가능성에 대해서 생각하는 것이 모든 통계의 기본이다. 

 

여기서 linear model이 등장한다.

 

Linear model

우리가 관심있는 것은 pitch가 성별에 따라 다른가 이다.

picth ~ sex

이렇게 쓰고 "pitch predicted by sex" 혹은 "pitch as a function of sex"로 표현한다. 

동시에 좌변은 dependent variable, 우변은 independent variable 혹은 explanatory variable, 혹은 predictor 라고 부른다. 

 

하지만 다른 이름이 있는데 그게 바로 "fixed effect" 이다. 

 

 

허나 세상은 그렇게 완벽하지 않고, sex 이외에도 pitch에 영향을 미치는 변수들이 굉장히 많다. 

심지어 실험에서 우리가 control할 수 없는 변수들이 많기 때문에 이러한 "random" factors를 포함하면 식을 다음과 같이 쓸 수 있다. 다음의 epsilon은 error term이라고도 부른다.

picth ~ sex + ε

 

실제로 데이터를 만들어 보면, 

pitch = c(233, 204,242,130,112,142)
sex = c(rep("female",3),rep("male",3))
my.df = data.frame(sex, pitch)

여기에 linear model을 만들어 보면, 

xmdl = lm(pitch ~ sex, my.df)
summary(xmdl)

 

위와 같은 결과를 얻을 수 있다. 

결과에서는 내가 만든 formula를 보여주고, 

residuals (잔차)와 coefficients, residual standard error, multiple R-squared 들을 보여주게 된다. 

하나씩 자세히 알아보자. 

 

Multiple R-squared

R제곱이라고도 하는 이 인자는 0에서 1 사이의 값을 갖는다.

전에 다뤘던 ANOVA에서 "ANOVA는 회귀분석의 특별한 형태"라고 했던 기억이 어렴풋이 있다. 

ANOVA에서는 변동을 이용하여 두 군의 차이를 분석하는데, 

회귀모형에서

 

총변동 SST = 설명 가능한 변동 SSR + 설명 불가한 변동 SSE

 

로 표현할 수 있고, 

 

R squared = SSR / SST = 1 - SSE / SST 가 된다. 

 

즉 R 제곱이 1에 가까울수록 우리의 모델의 설명력이 높아진다고 생각하면 된다.

 

예를 들어 위의 경우 R제곱이 0.921이면, 데이터셋에서 92.1%가 우리의 모델로 설명된다는 뜻이 되겠다. 

 

 

Adjusted R squared

이는 위의 R제곱과는 약간 다른데, R제곱이 얼마나 설명할 수 있는지에만 초점을 맞췄다고 하면,

Adjusted R squared는 얼마나 많은 fixed effect를 설명을 위해 사용하였는가에도 관심이 있다. 

 

 

Overall p-value

위에서 모델의 p-value는 0.002407이다.

즉 성별이 pitch에 영향이 없을 가능성이 0.002407이므로 모델이 통계적으로 유의하다고 할 수 있다. 

우리는 이 모델에 1개의 fixed effect만을 사용했기 때문에 coefficient table의 sexmale 오른쪽에 있는 probability 또한 0.00241로 이 값과 같은 것을 볼 수 있다. 

 

 

Coefficient table

 

근데 왜 sex가 아닌 sexmale이라고 table에 표시되어 있을까?

intercept 부분을 보면 이 값이 226.33Hz인 것을 볼 수 있는데 이것은 바로 여성들의 pitch 평균값이다. 

 

또한 sexmale의 값이 -98.33으로 음수임을 볼 수 있는데 226.33에서 -98.33을 해주면 남성들의 pitch 평균값인 128 정도가 나옴을 알 수 있다.

 

즉 우리가 만든 linear model은 데이터를 다음과 같이 보고 있는 것이다. 

x축의 female과 y축의 female pitch average가 그래프의 원점이 된다고 생각하면 편리하다. 

 

 

왜 female이 기준이 되었는가? -> lm() function이 알파벳 순서로 인식하기 때문

어떻게 categorical variable을 저런식으로 선형으로 연결할 수 있는가? -> 두 카테고리의 차이가 두 카테고리 사이의 기울기와 비례하기 때문

 

 

연속변수로 바꿔보면

age = c(14,23,35,48,52,67)
pitch = c(252,244,240,233,212,204)
my.df = data.frame(age,pitch)
xmdl = lm(pitch ~ age, my.df)
summary(xmdl)

위와 같은 결과를 얻을 수 있다. 

 

위에서의 estimate은 age가 0일 때의 pitch 값이 되겠다. 

이 estimate을 좀더 의미 있는 숫자로 바꾸고 싶다면 variable을 centering하면 된다. 

my.df$age.c = my.df$age - mean(my.df$age)
xmdl = lm(pitch ~ age.c, my.df)
summary(xmdl)

age와 mean의 차이로 새로운 변수를 만들어서, 이 변수로 모델을 만들어 주면

slope은 똑같고 intercept만 바뀐 것을 볼 수 있는데, age.c = 0일때가 결국 age = mean일때 이므로

평균 나이일 때의 pitch인 230.8333을 편하게 볼 수 있다.

 

 

Linear model이 "model"인 이유

4가지의 가정(assumption)이 포함되어 있기 때문이다. 

1) Linearity

이렇게 residual들에 box를 쳐서 박스를 돌려보면?

이러한 residuals vs fitted value 그래프를 얻을 수 있다. 

위의 경우 데이터가 많지 않아서 패턴이 보이지 않지만, 다음과 같은 curvy pattern을 보이는 경우도 있다. 

즉 이러한 경우에는 model의 linerity가 떨어진다고 볼 수 있다. 

이러할 때에는

• 중요한 fixed effect를 빼먹지 않았는지 확인하고
• response(종속변수에) non-linear transform을 해보거나
• fixed effect에 non-linear transform을 해보거나
• 카테고리컬 데이터의 경우 logistic model로 변경한다.

 

 

2) Absence of Collinearity

2개의 서로 다른 fixed effects가 서로 연관이 있는 경우에 collinear라고 한다.

이러한 변수들이 모델에 함께 포함되면 굉장히 unstable해지는데, 서로가 서로의 설명력을 "훔치기" 때문이다. 

 

3) Homoskedasticity or Absence of heteroskedasticity

등분산성이라고 보면 될 거 같다.

잔차들이 비슷한 정도로 떨어져 있다는 가정인데, 

이런식으로 blob-like appearance를 띄게 된다. 

 

이게 만족하지 않는 경우를 보면, 

fitted value가 증가할 수록 residual도 증가하는 것을 보이기 때문에 heteroskedasticity라고 할 수 있다. 

 

 

4) Normality of residuals

잔차들이 정규성을 띄어야 한다.

이는 잔차 히스토그램이나 qqplot을 그려보면 된다. 

 

 

5) Influential data points

영향점이라고도 하며, 특별히 모델에 영향을 크게 주는 데이터를 말한다.

R의 dfbeta 함수를 사용해 보면, 

여기 dfbeta value들이 나오게 되는데, 

특정 data point가 제외됐을 때 coefficient가 얼마나 수정되는지를 뜻한다. 

 

자세히 보면, 원래 age 변수의 coefficient는 -0.9099 였는데, 

data point 1이 제외되면 coefficient가 0.0643 만큼 바뀐다, 즉 -0.9099 minus 0.06437573 만큼 바뀐다는 의미이다.

(slope이 음수면 dfbeta가 빼지고, slope이 양수면 더해진다.)

 

따라서 slope의 절반정도 이상을 변경하는 data point는 influential data point로 보고 잘 확인해야 한다.

이러한 영향점을 그냥 삭제하는 것은 연구 윤리 위반이며

 

포함한 분석과 제거한 분석을 모두 보고하는 것이 중요하다. 

혹은 명확한 에러가 있는 경우 제거할 수도 있다. 

 

6) Independencies

선형 모델의 가장 중요한 가정은 바로 독립성이다.

즉 모든 데이터들이 독립이어야 한다. = 각 데이터들은 다른 individual로부터 측정된 것이어야 한다는 점이다.

 

여기서 반복측정된 데이터를 위해서 나온 개념이 바로 mixed model이다. 

 

 

 

Mixed model

 

Bodo Winter 씨의 연구에서 실제로 사용한 모델을 살펴 보자. 

목소리 pitch가 어떻게 달라지는 지를 연구한 것이다.

 

pitch ~ politeness + sex + ε

이미 살펴보았듯

politeness, sex : fixed effect

ε : error term, 컨트롤할 수 없는 부분이자 "probabilistic" or "stochastic" part

 

만약 한 대상에서 반복 측정한다면? 위에서 살펴본 독립성 가정이 침해되기 때문에

"subject"에 대한 random effect를 추가함으로써 이를 해결할 수 있다는 idea이다. 

 

예를 들어 subject 1은 baseline pitch가 233Hz, subject 2는 210Hz 인 것 처럼 말이다. 

이러한 개인 간의 차이를 random intercepts를 포함함으로써 해결한다.

 

이제 왜 mixed effect model이라는 말이 붙었는지 이해가 간다. 

기존의 선형 모델은 fixed effect를 제외하고는 error term으로 간주하였으나

(이해할 수 있는, 혹은 systematic한 fixed effect vs. 이해할 수 없고 컨트롤할 수 없는 error term)

 

mixed model에서는 fixed effect에 1개 이상의 random effect를 추가하게 된다. 

 

pitch ~ politeness + sex + (1|subject) + ε

(1|subject) 라는 표현은 "assume an intercept that’s different for each subject" 과 같다. 

 

 

한 가지 더 추가해 보자.  대화의 시나리오에 따라 pitch가 어떻게 달라지는지 알아보기 위해 

• 교수님에게 부탁하는 상황 : polite condition
• 동료에게 부탁하는 상황 : informal condition

등으로 condition을 나누었다. 

상황은 총 7가지이다. (부탁 / 지각한 것에 대한 변명 / 등등..)

 

 

위에서 살펴본 subject의 반복측정 상황에서 by-subject variation이 있던 것처럼

지금도 by-item(by-situation) variation이 생기게 된다.

즉 다음과 같이 모델을 수정할 수 있다. 

 

pitch ~ politeness + sex + (1|subject) + (1|item) + ε

 

즉 item에 따라 intercept를 달리 가져감으로써 독립이 아닌 문제를 해결하였다. 

 

이제 개인별 분석을 할 때에는 by-item variation을 무시하면 되고

item별 분석을 할 때에는 by-subject variation을 무시하면 된다!

 

R에서는 LME를 위해 "lme4"라는 패키지를 사용한다.

저자가 제공하는 dataset을 이용해 본다. 

library(lme4)
politeness=
  read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
  
which(is.na(politeness$frequency))
which(!complete.cases(politeness))

데이터 무결성을 함께 체크하였다. 

다음과 같은 구조이며, 각 subject의 gender, scenario, attitude, frequency(Hz)가 포함되어 있다. 

 

boxplot으로 attitude, gender에 따른 frequency의 분포를 살펴보는 코드는 다음과 같다.

boxplot(frequency ~ attitude*gender,
        col=c("white","lightgray"),politeness)

(이렇게 간단하다니................)

양쪽 성에서 모두 polite한 경우에 frequency가 낮음이 보인다.

이제 mixed model을 만들어 확인해 보자. 

 

politeness.model = lmer(frequency ~ attitude 
                        + (1|subject) + (1|scenario), data=politeness)
                        
summary(politeness.model)

 

 

위의 REML criterion 등은 일단 넘어가고, random effect를 살펴 보자.

 

boxplot에서 볼 수 있듯이 item 간의 변동 보다는 개인간의 변동이 더 크고

위에서 variance 항목으로 나타난다.

 

Residual은 scenario와 subject 모두와 관련 없는 항인데 이것이 우리의 새로운 ε가 되는 것이다.

 

그 아래에 있는 Fixed effects table은 linear model에서와 유사하다. (attitudepol & sexmale)

t value는 estimate을 Std. error로 나눈 값이다.

 

또한 여기서의 Intercept는 어떤 의미인가? -> informal condition에서의 data의 평균이다!

 

 

gender을 추가해보자!

politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politeness)

 

새로운 fixed effect인 gender을 넣어주고 나면,

subject와 관련된 variance가 크게 감소한 것을 알 수 있다!

 

 

 

 

LMM의 significance - Likelihood ratio test

불행히도 LM처럼 직관적인 p-value는 없기 때문에 저자는 likelihood ratio test라는 것을 소개하고 있다. 

쉽게 생각하면, fixed effect를 하나씩 제거(null model)해보고 full model과 비교하여 얼마나 달라지는지를 비교하는 방법이다. 

 

1개의 fixed effect (attitude)를 제거한 null model을 만들어 보자. 

politeness.null = lmer(frequency ~ gender +
                         (1|subject) + (1|scenario), data=politeness, REML=FALSE)
politeness.model = lmer(frequency ~ attitude +
                          gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)

** RMEL은 FALSE로 했는데 일단 모델끼리 비교할 때에는 이렇게 해야된다고 알고 넘어가자

** 그 이유는 (Pinheiro & Bates, 2000; Bolker et al., 2009) 에 나와있다고 한다.

 

두 모델을 anova로 비교하면 된다.

** 이 이유는 "Wilk's theorem" 때문이라고 하나 기술적으로 복잡해서 넘어간다.

anova(politeness.null,politeness.model)

 

즉 다음과 같이 기술할 수 있다. 

 

“… politeness affected pitch (χ2 (1)=11.62, p=0.00065), lowering it by about 19.7 Hz ± 5.6 (standard errors) …”

 

 

이제는 gender의 영향도 궁금하기 때문에 다음 처럼 할 수도 있다. 

위의 null model은 fixed effect를 모두 제거하므로 frequency의 평균값을 담고 있다. (intercept only model)

이 두 모델이 유의미하게 다르다고 해도, 이게 gender 때문인지 attitude 때문인지 구분이 안된다.

 

추가적으로, 두 변수간에 상호작용을 확인하고 싶은 경우

(예를 들어 남성과 여성 간에 formal / informal하게 말하는 것이 반대 작용을 하거나

혹은 여성에서는 차이가 있고 남성에서는 차이가 없는 경우?)

 

다음과 같이 테스트할 수 있다. 

 

 

Random slopes vs Random intercepts

다음 명령어로 model의 coefficient를 살펴 보자.

coef(politeness.model)

각 item들과 subject 들이 다른 intercept를 가진 것을 알 수 있다.

동시에 attitude와 gender는 모두 같은 값들이 배정된 것을 볼 수 있다!! -> item, subject와 관련 없이 attitude와 gender가 미치는 영향은 모두 같다!

이게 우리가 이 모델을 Random intercept model이라고 부르는 이유이다.

 

하지만 이게 아닐 수도 있다.

예를 들면 어떤 item(상황)에서는 좀더 polite하게 될 수도 있고, 반대가 될 수도 있기 때문이다.

즉, politeness의 영향은 item에 따라, subject에 따라 다를 수 있다!

이때 필요한 것이 Random slope model이다. 

politeness.model = lmer(frequency ~ attitude + gender 
                        + (1+attitude|subject) 
                        + (1+attitude|scenario),
                        data=politeness,
                        REML=FALSE)

 

위의 “(1+attitude|subject)”의 뜻은 

• 1 : intercept가 subject에 따라 달라진다.
• attitude : attitude의 계수도 subject에 따라 달라진다.

coefficient를 확인해 보면

이제 item과 subject에 따라 intercept 뿐만 아니라 attitude의 계수도 달라진 것을 볼 수 있다!!!

 

 

이제 새롭게 구한 Random slope model의 p-value 비스무리 한 것을 구해보자. 

 

politeness.model = lmer(frequency ~ attitude + gender 
                        + (1+attitude|subject) 
                        + (1+attitude|scenario),
                        data=politeness,
                        REML=FALSE)
politeness.null = lmer(frequency ~ gender 
                       + (1+attitude|subject) 
                       + (1+attitude|scenario),
                       data=politeness, 
                       REML=FALSE)
                       
anova(politeness.null,politeness.model)

attitude에 대한 null model과 비교했을 때 유의미하게 다른 것을 알 수 있다!

 

 

이 random slope를 어디까지 적용해야 하느냐, 꼭 해야 하느냐에 대한 질문에 대해서

저자는 최대한 해라! 라고 대답하고 있다. 

 

 

 

추가적으로..

데이터 영향점 분석을 위한 dfbeta는 lme에는 사용 불가능 하므로 influence.ME 패키지를 사용하자.

 

 

Write up

이런 글을 써주기 위해서는

citation()함수를 이용하여

citation()
citation(lme4)

R과 패키지 버젼을 꼭 써주자!