보정 평균 (adjusted mean)은 성별, 연령 등 여러 가지 요소를 보정한 평균을 의미하며 각 그룹간 평균값의 차이를 더 분명하게 알려줄 수 있는 통계 분석 방법이라고 할 수 있습니다. R에서는 SAS의 LSMEANS과 같은 기능을 하는 패키지인 lsmeans 패키지가 널리 사용됩니다. 비록 이 기능은 emmeans이라는 새로운 패키지에 기능이 통합될 예정이지만 lsmeans 패키지를 오래 써왔고 기본 문법은 비슷한 것 같아 여기서는 lsmeans를 설명합니다. 예제는 앞서 소개한 moonBook 패키지의 acs 데이터를 사용합니다. 이 데이터는 심근 경색이나 협심증으로 병원을 방문한 환자의 데이터를 공개한 것입니다. 참고로 앞서 설명한 사후 검정 - Tukey's HSD test에서 이어지는 이야기라고 할 수 있습니다.
require(moonBook)
require(lsmeans)
str(acs)
앞서 포스트에서 설명했던 예제와 비슷하게 BMI에 따라 환자를 정상, 과체중, 비만 세 그룹으로 나눠 세 그룹간 나이의 차이가 있는지를 검증해 보겠습니다. 기본 아이디어는 과체중이나 비만인 경우 더 젊은 나이에 심근 경색이나 협심증이 생긴다는 것입니다. 이를 위해 각 그룹간 나이의 평균에 차이가 있는지를 분산분석 (ANOVA)로 검증합니다.
acs$obesity2[acs$BMI<23 span="">
acs$obesity2[acs$BMI>=23]=1
acs$obesity2[acs$BMI>=25]=2
table(acs$obesity2)
out=lm(age~factor(obesity2), data=acs)
anova(out)
out=aov(age~factor(obesity2), data=acs)
TukeyHSD(out)
> out=lm(age~factor(obesity2), data=acs)
> anova(out)
Analysis of Variance Table
Response: age
Df Sum Sq Mean Sq F value Pr(>F)
factor(obesity2) 2 3504 1752.0 13.133 2.469e-06 ***
Residuals 761 101517 133.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> out=aov(age~factor(obesity2), data=acs)
> TukeyHSD(out)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = age ~ factor(obesity2), data = acs)
$`factor(obesity2)`
diff lwr upr p adj
1-0 -2.012012 -4.531198 0.5071740 0.1464941
2-0 -4.964351 -7.253721 -2.6749811 0.0000013
2-1 -2.952339 -5.437984 -0.4666936 0.0149519
이 결과를 해석하면 정상 체중과 과체중 사이에는 유의한 차이가 없지만, 정상체중/과체중-비만 사이에는 유의한 차이가 있는 것으로 해석할 수 있습니다. 여기서 한 걸음 더 나아가 성별에 따른 차이를 보정하면 ANCOVA가 됩니다.
out=lm(age~factor(obesity2)+factor(sex), data=acs)
anova(out)
out=aov(age~factor(obesity2)+factor(sex), data=acs)
TukeyHSD(out)
> out=lm(age~factor(obesity2)+factor(sex), data=acs)
> anova(out)
Analysis of Variance Table
Response: age
Df Sum Sq Mean Sq F value Pr(>F)
factor(obesity2) 2 3504 1752.0 14.735 5.269e-07 ***
factor(sex) 1 11153 11152.6 93.798 < 2.2e-16 ***
Residuals 760 90364 118.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> out=aov(age~factor(obesity2)+factor(sex), data=acs)
> TukeyHSD(out)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = age ~ factor(obesity2) + factor(sex), data = acs)
$`factor(obesity2)`
diff lwr upr p adj
1-0 -2.012012 -4.390364 0.3663396 0.1161858
2-0 -4.964351 -7.125734 -2.8029677 0.0000003
2-1 -2.952339 -5.299024 -0.6056529 0.0090510
$`factor(sex)`
diff lwr upr p adj
Male-Female -8.097007 -9.739297 -6.454718 0
공분산분석 역시 비슷한 결과가 나왔습니다. 그런데 이렇게 보정했을 때 각 그룹간 나이의 차이는 어느 정도일지 눈으로 알기는 어렵습니다. 사후 검정을 포함해 보정 평균을 구하기 위해 lsmeans를 사용합니다.
out=lm(age~factor(obesity2)+factor(sex), data=acs)
marginal = lsmeans(out, ~ obesity2)
cld(marginal, alpha=0.05, sort = FALSE, Letters=letters, adjust="tukey") #compact letter display
> out=lm(age~factor(obesity2)+factor(sex), data=acs)
> marginal = lsmeans(out, ~ obesity2)
> cld(marginal, alpha=0.05, sort = FALSE, Letters=letters, adjust="tukey") #compact letter display
obesity2 lsmean SE df lower.CL upper.CL .group
0 66.69924 0.6719680 760 65.09121 68.30727 a
1 64.96966 0.7815843 760 63.09932 66.84000 a
2 62.02618 0.6576726 760 60.45236 63.60000 b
Results are averaged over the levels of: sex
Confidence level used: 0.95
Conf-level adjustment: sidak method for 3 estimates
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
성별로 보정했을 때 평균값은 66.69924, 64.96966, 62.02618로 나왔으며 lower.CL upper.CL 에 95% CI의 범위도 함께 나왔습니다. 이렇게 보니 정상/과체중 간의 차이보다 비만군과의 차이가 더 크다는 것을 알 수 있습니다. 만약 더 많은 변수를 보정하면 어떨까요?
out=lm(age~factor(obesity2)+factor(sex)+factor(smoking)+HDLC+factor(HBP)+factor(DM), data=acs)
marginal = lsmeans(out, ~ obesity2)
cld(marginal, alpha=0.05, sort = FALSE, Letters=letters, adjust="tukey") #compact letter display
> out=lm(age~factor(obesity2)+factor(sex)+factor(smoking)+HDLC+factor(HBP)+factor(DM), data=acs)
> marginal = lsmeans(out, ~ obesity2)
> cld(marginal, alpha=0.05, sort = FALSE, Letters=letters, adjust="tukey") #compact letter display
obesity2 lsmean SE df lower.CL upper.CL .group
0 66.55850 0.6859582 743 64.91691 68.20009 a
1 64.15678 0.7803798 743 62.28923 66.02433 b
2 61.10558 0.6589566 743 59.52861 62.68255 c
Results are averaged over the levels of: sex, smoking, HBP, DM
Confidence level used: 0.95
Conf-level adjustment: sidak method for 3 estimates
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
값에 차이가 생기면서 세 그룹 모두 유의한 차이가 있는 것으로 나타났습니다. (참고로 lm 대신 glm을 사용해도 같은 결과가 나옵니다.) 잘 보면 서로 95% CI 값이 겹치는 경우에도 차이가 있는 것으로 나오는 데 당연히 사후 검정 결과와 다를 수 있습니다.
lsmean 기능은 공분산분석이나 선형회귀 분석의 단점을 보완할 수 있습니다. 실제로는 매우 작은 차이인데 통계적으로는 샘플수가 많거나 측정이 정밀하면 마치 유의한 차이가 있는 것처럼 나타날 수 있습니다. 이럴 때 보정 평균을 구하면 실제값 사이에는 별 차이가 없는지 아니면 그래도 의미 있는 차이가 있는지를 한눈에 볼 수 있을 것입니다.
참고로 lsmean 기능은 plot 기능과 연결해서 사용할 수 있습니다.
out=lm(age~factor(obesity2)+factor(sex)+factor(smoking)+HDLC+factor(HBP)+factor(DM), data=acs)
marginal = lsmeans(out, ~ obesity2)
out2=cld(marginal, alpha=0.05, sort = FALSE, Letters=letters, adjust="tukey")
plot(out2)
이 역시 plot의 다양한 옵션과 함께 사용할 수 있으며 다른 분석 패키지와 연동할 수 있습니다. 좀 더 다양한 사용에 대해서는 아래 글을 참조해 주십시요.
댓글
댓글 쓰기