6. 공분산분석
MK
2021 3 13
6.1 공분산분석(ANCOVA) 예제1
- moonBook 패키지에 포함되어 잇는 acs 데이터를 이용하여 공분산분석을 해보고자 합니다. 급성관상동맥증후군의 발병나이와 흡연의 가왕력을 비교하기 위해 먼저 boxplot과 densityplot을 그려봅시다.
require(moonBook)
data(acs)
acs$smoking=factor(acs$smoking, levels=c("Never","Ex-smoker","Smoker"))
plot(age~smoking,data=acs,color=c("red","green","blue"))
- 흡연의 기왕력에 따른 나이의 차이를 분산분석(ANOVA) 해봅시다.
## Analysis of Variance Table
##
## Response: age
## Df Sum Sq Mean Sq F value Pr(>F)
## smoking 2 12164 6082.2 49.497 < 2.2e-16 ***
## Residuals 854 104939 122.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## lm(formula = age ~ smoking, data = acs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.745 -7.458 0.491 8.491 32.542
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.5090 0.6084 109.322 <2e-16 ***
## smokingEx-smoker -0.7639 0.9861 -0.775 0.439
## smokingSmoker -8.0511 0.8677 -9.279 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.09 on 854 degrees of freedom
## Multiple R-squared: 0.1039, Adjusted R-squared: 0.1018
## F-statistic: 49.5 on 2 and 854 DF, p-value: < 2.2e-16- ANOVA에서 흡연자는 비흡연자에 비해 보다 어린 나이인 것을 알 수 있고, 흡연자에서 어린 나이에 acs가 발병하는 것을 알 수 있습니다.
- 여기에 비만 정도를 나타내는 BMI를 추가하여 공분산분석(ANCOVA)를 시행해봅시다.
##
## Call:
## lm(formula = age ~ smoking + BMI, data = acs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.817 -7.410 0.157 8.165 29.919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 85.0774 2.9141 29.195 < 2e-16 ***
## smokingEx-smoker -1.1299 1.0352 -1.091 0.275
## smokingSmoker -8.0736 0.8948 -9.023 < 2e-16 ***
## BMI -0.7669 0.1176 -6.521 1.28e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.86 on 760 degrees of freedom
## (93 observations deleted due to missingness)
## Multiple R-squared: 0.147, Adjusted R-squared: 0.1436
## F-statistic: 43.65 on 3 and 760 DF, p-value: < 2.2e-16colors=rainbow(acs$smoking)
colors=colors[c(3,2,1)]
plot(age~BMI, col=colors, data=acs)
legend("topright", legend=levels(acs$smoking), pch=1, col=colors,lty=1)
curve(85.077-0.767*x, col=colors[1],add=TRUE)
curve(85.077-0.767*x-1.130, col=colors[2],add=TRUE)
curve(85.077-0.767*x-8.074, col=colors[3],add=TRUE)
- smoking의 세 그룹을 각각 비교하기 위해 tukey 방법으로 다중비교를 해보면 다음과 같습니다.
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = age ~ smoking + BMI, data = acs)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Ex-smoker - Never == 0 -1.1299 1.0352 -1.091 0.518
## Smoker - Never == 0 -8.0736 0.8948 -9.023 <1e-04 ***
## Smoker - Ex-smoker == 0 -6.9437 1.0399 -6.677 <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)- 담배를 끊은 사람과 비흡연자 사이에 급성관상동맥증후군의 발병연령에는 차이가 없으나 흡연자는 비흡연자에 비해 8년, 담배를 끊은 사람에 비해 6.9년 일찍 발생합니다.
6.2 공분산분석(ANCOVA) 예제2
- moonBook 패키지의 radial 데이터가 포함되어 있습니다. 여기서 동맥경화 정도를 나타내는 것은 NTAV(normalized total atheroma volume)입니다. 이 연구에서는 NTAV의 변화를 나이와 성별을 예측변수로 ANCOVA 모형을 만드는 것입니다.
##
## Call:
## lm(formula = NTAV ~ age, data = radial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.231 -14.626 -4.803 9.685 100.961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.3398 14.6251 3.032 0.00302 **
## age 0.3848 0.2271 1.694 0.09302 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.94 on 113 degrees of freedom
## Multiple R-squared: 0.02477, Adjusted R-squared: 0.01614
## F-statistic: 2.87 on 1 and 113 DF, p-value: 0.09302plot(NTAV~age, data=radial, col=ifelse(radial$sex=="M","blue","red"))
abline(out,col="red",lwd=2)
title(expression(italic(NTAV==0.385%*%Age+44.34)),family="Times")
- 남자는 파란색, 여자는 빨간으로 구별하여 plot()을 해보면 남자가 여자보다 동맥경화 정도가 심한 것을 알 수 있습니다.
- 하지만 이와 같이 단순회귀모형을 적용할 때 나이가 0세일 때도 동맥경화가 있다는 것은 모순입니다. 따라서 이런 경우에는 y절편이 없는 non-intercept 모형을 선택해야 합니다. Non-intercept 모델을 적용하려면 lm 함수를 다음과 같이 사용합니다. lm(y~x-1)
##
## Call:
## lm(formula = NTAV ~ age - 1, data = radial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.235 -13.726 -4.613 10.168 98.343
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age 1.06532 0.03589 29.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.78 on 114 degrees of freedom
## Multiple R-squared: 0.8854, Adjusted R-squared: 0.8844
## F-statistic: 881 on 1 and 114 DF, p-value: < 2.2e-16plot(NTAV~age, data=radial, col=ifelse(radial$sex=="M","blue","red"))
abline(out1,col="red",lwd=2)
title(expression(italic(NTAV==1.065%*%Age)),family="Times")
- 하지만 그래프에서는 나이가 많을수록 NTAV의 분포가 넓어지는 것을 볼 수 있습니다. 이 회귀모형의 회귀분석의 가정을 만족시키는지 보기 위해 회귀진단을 해보겠습니다.
- 회귀진단 결과 역시 나이가 많을수록 잔차가 증가하고 Normal Q-Q plot에서도 벗어나는 것을 알 수 있습니다. 이와 같은 결과가 나타나는 이유는 회귀진단의 기본가정인 반응변수가 정규분포해야한다는 가정을 위반했기 때문입니다.
- NTAV의 분포는 왼쪽으로 심하게 쏠려 있습니다. 이러한 경우는 제곱근변환이나 로그변환이 필요한 경우 입니다. car 패키지의 powerTransformation() 함수를 써봅시다.
## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## radial$NTAV -0.0533 0 -0.4229 0.3164
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 0.07953273 1 0.77793
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 29.48666 1 5.6303e-08- 결과를 보면 NTAV를 정규화시키려면 NTAV^(-0.053)을 사용할 것을 제시하고 있습니다. -0.053은 0에 가깝기 때문에 모형을 정규성 적합시키기 위해 로그변환을 시도해보겠습니다. 또한 성별을 회귀진단에 포함시켜 ANCOVA를 해보면 다음과 같습니다.
##
## Call:
## lm(formula = log(NTAV) ~ (age - 1) + male, data = radial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.29896 -0.26180 0.05829 0.42898 1.93835
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age 0.059381 0.001116 53.19 < 2e-16 ***
## male 0.655927 0.101216 6.48 2.48e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5681 on 113 degrees of freedom
## Multiple R-squared: 0.9819, Adjusted R-squared: 0.9816
## F-statistic: 3072 on 2 and 113 DF, p-value: < 2.2e-16##
## Shapiro-Wilk normality test
##
## data: resid(out2)
## W = 0.99116, p-value = 0.671- 회귀진단 결과 성별과 나이가 모두 유의하며 shapiro.test() 결과 잔차의 정규분포 가정 또한 만족합니다. 로그변환된 모형의 adjusted R^2는 0.982로 반응변수 분산의 98.2%를 설명하고 있음을 알 수 있습니다.
plot(log(NTAV)~age, data=radial, col=ifelse(radial$sex=="M","blue","red"))
curve(0.0594*x, col='red', lty=2, add=TRUE)
curve(0.0594*x+0.656, col='blue', lty=1, add=TRUE)
legend("topleft", legend=c("Male","Female"), col=c("blue","red"),pch=1,lty=1:2)
- 반응 변수의 로그를 없앤 상태로 회귀선을 그려보면 다음과 같습니다.
plot(NTAV~age, data=radial, col=ifelse(radial$sex=="M","blue","red"))
curve(exp(0.0594*x), col='red', lty=2, add=TRUE)
curve(exp(0.0594*x+0.656), col='blue', lty=1, add=TRUE)
title(expression(NTAV==e^(0.0594%*%age+0.656%*%male)), family="Times")
legend("topleft", legend=c("Male","Female"), col=c("blue","red"),pch=1,lty=1:2)
- 남여의 표시를 각각 파란 M 빨간 F로 표시하려면 pch 부분만 추가하면 됩니다.
plot(NTAV~age, data=radial, col=ifelse(radial$sex=="M","blue","red"), pch=as.character(sex))
curve(exp(0.0594*x), col='red', lty=2, add=TRUE)
curve(exp(0.0594*x+0.656), col='blue', lty=1, add=TRUE)
title(expression(NTAV==e^(0.0594%*%age+0.656%*%male)), family="Times")
legend("topleft", legend=c("Male","Female"), col=c("blue","red"),pch=1,lty=1:2)