앞에서 계속 random effect model 이야기를 했으니 이제 예시를 하나 들 때도 된 것 같군요. 이 글에서 다루는 예시는 아주 유명한 데이터셋인데, 미국에 있는 8개의 고등학교에서 새로 도입한 교수학습법이 SAT 점수에 미친 효과를 측정한 자료입니다. [1] 통칭 “eight-schools example” 이라고 부르는데, 유명한 베이지안통계 교과서인 Bayesian Data Analysis에 실리면서 더욱 유명해진 감이 있습니다. 이 글의 내용 전개도 이 책의 서술 순서를 대체로 따라갑니다.
자료는 아주 간단합니다. 다음과 같습니다:
학교 효과 표준오차
A 28.39 14.9
B 7.94 10.2
C -2.75 16.3
D 6.82 11.0
E -0.64 9.4
F 0.63 11.4
G 18.01 10.4
H 12.16 17.6
각각의 학교에서 점수 상승(하강) 효과와 그 표준오차를 기록한 것을 알 수 있습니다. 이제 물음을 던집니다: 새로 도입한 교수학습법의 효과는 얼마인가? 어떻게 계산해야 할까요? 제 1감으로 떠오르는 것은 단순하게 모든 학교에서 관찰된 효과의 평균을 내는 것입니다. 그러면 8.82라는 숫자를 얻습니다.
그런데 이 숫자를 효과 값으로 사용하는 데에는 다소 문제가 있습니다. 일단 이 숫자는 학교별로 관찰된 편차를 깡그리 무시합니다. 다시 말해, 관측된 효과 값들 중에는 8.82라는 숫자와 상당히 동떨어져 보이는 값들이 있습니다. 이를테면 A 학교에서 관찰된 효과 값인 28.39는 8.82와는 상당히 동떨어진 값입니다. 그러면 A 학교에서 관측된 값은 어쩌면 그냥 우연히 좀 ‘튀어’ 보이는 값이 나온 것은 아니었을까요? 조금 의문스럽습니다.
그런데 사실 이 값 - 8.82 - 은 효과 크기들에 대한 fixed effect model의 (가중치를 사용하지 않은) 추정치입니다. 여기서 fixed effect model은 다음과 같이 정의됩니다 (\(y_i\)는 학교별 관측된 효과, \(\mu\)는 효과크기의 참값, \(e_i\)는 학교별 오차):
다시 말해 fixed effect model에서는 교수학습법의 효과의 참값은 모든 학교에서 똑같다고 간주합니다 (그래서 fixed effect model 입니다). 다르게는 “complete pooling” 방식이라고도 부릅니다. 그리고 이 방식을 적용하면 효과의 참값인 \(\mu\)에 대한 추정치는 그냥 앞에서 구한 (산술)평균이 됩니다. 그런데 그렇게 보기에는 학교 간 편차가 심해 보이고, 하나 더 문제가 되는 것은 사용할 수 있었던 다른 정보, 즉 표준오차가 무시됐다는 것입니다. 그런데 표준오차는 관측된 효과 크기들의 정밀도를 나타내는 척도라는 점에서 우려스럽습니다. 이를테면 H 학교에서 관측된 효과의 표준오차는 17.6인데, 이는 H 학교에서 관찰한 12.16이라는 값이 다른 학교에서 관측된 값들보다 - 더 큰 표준오차로 인해 - 더 불확실함을 의미합니다. 그런데 앞에서 평균을 구할 때 각 학교에 주어진 가중치는 1/8로 동일했습니다. 이게 정말 합당한 방식일까요? H 학교의 관측된 값에 대해서는 보다 적은 가중치가 주어져야 하지 않을까요? [2] 정리하자면, 이 fixed effect model은 좀 문제가 있어 보입니다.
그래서 다르게 생각해 보자면, 사실 학교마다 교수학습법의 효과가 갖는 참값은 다르며, 거기에 약간의 오차가 덧붙여져 실제로 관측된 값이 얻어졌다고 보는 게 더 타당할 수도 있겠습니다. (사실 학교마다 새로운 교수학습법을 받아들일 준비가 된 정도도 당연히 다르겠죠.) 이제 학교별 효과 크기가 다르다고 가정하고, 각각의 효과크기를 \(\theta_i\)라고 표기합시다. 그러면 새로운 통계모형은 다음과 같이 주어집니다:
(2)의 구조를 잘 보시면 \(y_i\)는 \(\theta_i\) (절편) 과 \(e_i\) (오차) 만으로 구성된 모델인데, 절편 자체도 분포를 갖는 모형이라는 것을 알 수 있습니다. 따라서 이 모형은 random effect 모델입니다. (좀 더 정확하게는 random intercept model 이라고 볼 수도 있겠죠.) 여기서 \(\sigma_i\) 들은 알고 있다고 가정됩니다. 즉, 위 자료에서 표준오차들이 \(\sigma_i\)의 노릇을 합니다. 이렇게 만든 모형을 핏팅하면 관측된 정보를 모두 사용하면서, 보다 현실적으로 말이 되는 값을 추정할 수 있습니다. [3] 실제로 모형을 추정할 때 관심의 대상이 되는 것은 \(\mu\) (참 효과크기들의 평균), \(\tau^2\) (학교별로 효과의 참값들이 퍼진 정도), 그리고 \(\theta_i\) (각 학교별 효과의 참값) 입니다.
이제 이 random effect model을 자료에 핏팅하겠습니다. 여기서는 R의 bayesmeta 패키지를 이용하겠습니다. 이 패키지에는 위에서 말한 자료가 Rubin1981 이라는 이름으로 들어있습니다. (참고로 코드는 [4] 에서 가져왔습니다):
> library(bayesmeta)
> schools <- bayesmeta(y=Rubin1981[,"effect"], sigma=Rubin1981[,"stderr"],
labels=Rubin1981[,"school"],
tau.prior=function(x){return(dhalfcauchy(x, scale=25))})
결과를 summary() 함수로 가져옵니다. 여기서는 아웃풋의 일부만 보여드립니다:
> summary(schools)
marginal posterior summary:
tau mu theta
mode 0.000000 7.997318 7.924963
median 4.768766 8.038269 7.995878
mean 5.899193 8.067574 8.067574
sd 4.886720 5.033346 9.164662
95% lower 0.000000 -1.810487 -10.597268
95% upper 15.207747 17.977655 26.996611
여기서는 \(\tau\), \(\mu\) 의 추정치를 보겠습니다. \(\tau\)의 mode는 0인데 이는 maximum a posteriori에 해당되며, \(\tau\)의 가장 그럴듯한 값은 0이라는 의미입니다. 하지만 이는 \(\tau\)의 값에 대한 불확실성을 전혀 고려하지 않은 수치죠. [5] 그래서 이를 고려한 값인 posterior mean을 봅니다. 그러면 \(\tau\)의 posterior mean은 5.90 가량인데, 이는 학교 간 효과 참값의 변동성이 대략 5.9 정도 된다는 의미입니다. 이 값이 위에 말한 random effect model에서의 \(\tau\)에 대한 추정치가 되겠습니다. 그리고 같은 요령으로 평균 효과크기인 \(\mu\)에 대해서도 보면, \(\mu\)의 가장 그럴듯한 값은 8.00이며 사후평균은 8.07로 서로 상당히 비슷함을 알 수 있습니다.
이제 각 학교에서의 효과의 참값에 대한 추정치들을 좀 들여다보겠습니다:
> schools$theta
A B C D E F G H
y 28.39 7.94 -2.75 6.82 -0.64 0.63 18.01 12.16
sigma 14.90 10.20 16.30 11.00 9.40 11.40 10.40 17.60
mode 8.91 7.94 7.48 7.83 6.82 7.28 8.99 8.08
median 10.16 7.97 6.99 7.78 6.13 6.78 9.85 8.33
mean 11.17 7.99 6.56 7.75 5.71 6.43 10.40 8.54
sd 7.92 6.17 7.45 6.37 6.28 6.61 6.64 7.48
95% lower -3.18 -4.40 -9.23 -5.15 -7.30 -7.35 -2.11 -6.50
95% upper 28.29 20.45 21.15 20.57 17.66 19.27 24.32 24.22
다 보기는 귀찮으니 (posterior) mean 부분만 보겠습니다. 예를 들어 A학교의 경우 관측된 효과 값은 28.39인데 참값에 대한 posterior mean은 11.16 입니다. 전체 평균인 8에 상당히 가깝게 당겨졌죠? 역시 28이라는 값은 우연히 관측된 값으로, 높은 표준오차 (14.9) 를 감안할 때 그렇게 믿을만하지 않다는 직관을 잘 반영해 주는 분석결과입니다. 반대로 C 학교의 경우를 보면, 관측된 값은 -2.75이나 추정된 효과 참값은 6.56으로 역시 전체 평균에 꽤 가깝게 당겨졌습니다. 이와 같이 hierarchical model 에서 추정치들이 전체 평균에 가깝게 당겨지는 현상은 매우 일반적인 것으로, 통계학에서는 shrinkage 라고 부릅니다. shrinkage는 특히 위에서 말한 complete pooling 방식과, 각 학교의 추정치들을 모두 각자의 효과의 참값으로 간주하는 “no pooling” 방식이라는 양 극단을 피해, 적절히 조절된 값으로 참값을 추정해 준다는 이점을 갖고 있습니다. (그런 의미에서 이를 “partial pooling” 이라고 부르기도 합니다.)
지금까지 가장 간단한 random effect 모형을 R 패키지를 사용하여 핏팅하고 그 결과를 살펴보았습니다. 한 가지 주의할 점은 \(\tau\)에 대한 사전분포에 따라 shrinkage의 강도가 상당히 달라질 수 있다는 것입니다. \(\tau\)의 값이 0에 가깝다는 사전 믿음이 강할수록 학교 간 차이가 적다는 편향이 생겨서, shrinkage가 꽤 강하게 들어갈 것이며, 반대로 0보다 크다는 사전믿음이 강할수록 학교 간 차이가 클 것이라는 편향이 생겨서 shrinkage가 약하게 작용할 것입니다. 따라서 prior의 강도에 따라 추론 결과가 어떻게 달라지는지 보는 sensitivity analysis가 사후에 이루어지는 것이 바람직합니다.
참고문헌
[2] 대개 inverse variance weighting이라는 방식으로 해결합니다. 이는 분산의 역수를 가중치로 한 평균을 내는 것을 의미합니다.
[3] 사실 이 모형은 random-effect meta-analysis 모형으로도 불립니다.
[4] 예제
[5] 빈도주의 메타분석 모형은 이런 이유로 \(\tau\)를 0으로 간주한 분석 결과를 제시합니다. 이것이 불만족스러워서 여기서 사용하지 않았습니다.
[6] 더 읽을거리