지난번에 이어 interaction plot을 ggplot2와 ggiraph 패키지를 사용하여 구현하였다. 오늘 소개할 함수는 다음과 같다.
자동차의 연비에 관한 mtcars데이타에서 연비(mpg)를 반응(종속)변수로, 공차중량(wt)과 마력(hp)를 영향(독립)변수로 하여 다중회귀모형을 만들어 보자. 이때 반응변수의 상호작용을 보기 위해 hp/wt를 입력한다. 참고로 A/B는 완전교차(complete crossing)를 의미하며 AB= A+B+A:B(상호작용)이다. AB*C는 A+B+C+A:B+A:C+B:C+A:B:C 를 의미한다.
fit=lm(mpg~hp*wt,data=mtcars)
summary(fit)
Call:
lm(formula = mpg ~ hp * wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-3.0632 -1.6491 -0.7362 1.4211 4.5513
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
hp -0.12010 0.02470 -4.863 4.04e-05 ***
wt -8.21662 1.26971 -6.471 5.20e-07 ***
hp:wt 0.02785 0.00742 3.753 0.000811 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.153 on 28 degrees of freedom
Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13회귀모형을 살펴보면 hp도 의미가 있고 wt도 의미가 있으며 hp:wt의 상호작용도 의미가 있다. 이와 같은 상호 작용을 시각화 하기 위해 ggEffect()함수를 만들었다. 사용방법은 다음과 같다.
source("ggAncova.R")
ggEffect(fit,interactive=TRUE)ggEffect()함수의 원형은 다음과 같다.
ggEffect=function(fit,
x=1,
probs=c(0.10,0.5,0.90),
point=TRUE,
xvalue=NULL,
se=FALSE,
use.rownames=FALSE,
interactive=FALSE){
...
}인수중 probs의 기본 값은 probs=c(0.10,0.5,0.90)이며 이는 영향변수 x축 변수가 아닌 변수(여기서는 hp)의 10,50,90퍼센타일에 대한 회귀선을 그려준다. 만일 5,25,50,75,95퍼센타일의 5개의 회귀선을 그리고자 한다면 다음과 같이 하면 된다.
ggEffect(fit,probs=c(0.05,0.25,0.5,0.75,0.95),interactive=TRUE)두번째 영향변수를 x축변수로 그림을 그리려면 x=2로 주면 되고 점의 tooltip을 데이타의 rownames로 바꾸려면 use.rownames인수를 TRUE로 변경하면 된다.
ggEffect(fit,x=2,use.rownames=TRUE,interactive=TRUE)저자가 만든 moonBook 패키지에 있는 radial데이타에서 동맥경화 정도(NTAV) 에 미치는 나이와 흡연과의 상호작용을 시각화해보면 다음과 같다.
require(moonBook)
fit2=lm(NTAV~age*smoking,data=radial)
ggEffect(fit2,interactive=TRUE)ANCOVA는 하나의 반응변수와 하나의 covariate가 있는 모형이다. mtcars 데이타에서 cyl는 엔진의 개수로 4기통, 6기통, 8기통으로 입력되어 있다. 이 데이타에서 연비(mpg)를 반응변수로 하고 공차중량(wt)를 공변량으로 하고 엔진수를 범주형 변수로 하여 ANCOVA를 시행해보면 다음과 같다.
fit=lm(mpg~wt+factor(cyl),data=mtcars)
summary(fit)
Call:
lm(formula = mpg ~ wt + factor(cyl), data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-4.5890 -1.2357 -0.5159 1.3845 5.7915
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 33.9908 1.8878 18.006 < 2e-16 ***
wt -3.2056 0.7539 -4.252 0.000213 ***
factor(cyl)6 -4.2556 1.3861 -3.070 0.004718 **
factor(cyl)8 -6.0709 1.6523 -3.674 0.000999 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.557 on 28 degrees of freedom
Multiple R-squared: 0.8374, Adjusted R-squared: 0.82
F-statistic: 48.08 on 3 and 28 DF, p-value: 3.594e-11위의 분석 결과에서 엔진수가 4인 경우, 6인 경우, 8인 경우 mpg와 wt사이의 회귀직선의 intercept가 각기 다르다. 이러한 ANCOVA 모형을 다음과 같이 시각화할 수 있다.
ggAncova(mtcars,"mpg","wt","cyl",interactive=TRUE)저자가 만든 moonBook패키지의 radial데이타는 동맥경화정도(NTAV)를 나이와 흡연유무에 따라 분석해보면 다음과 같다.
ggAncova(radial,"NTAV","age","smoking",interactive=TRUE)다음은 자료를 두 개의 범주형 변수로 나누어 평균과 표준편차/표준에러를 구하고 그 자료를 이용하여 그래프를 그려주는 ggCatepillarPlot()을 소개한다.
예를 들어 moonBook패키지의 acs데이타 중 환자의 나이를 성별과 진단명으로 각각 평균과 표준편차를 구하면 다음과 같이 할 수 있다.
require(plyr)
ddply(acs,.(sex,Dx),summarize,mean=mean(age,na.rm=T),sd=sd(age,na.rm=T)) sex Dx mean sd
1 Female NSTEMI 70.88000 11.35268
2 Female STEMI 69.10714 10.36214
3 Female Unstable Angina 67.71895 10.67366
4 Male NSTEMI 61.14563 11.56941
5 Male STEMI 59.42727 11.72396
6 Male Unstable Angina 61.44130 10.56934이 결과를 이용하여 다음과 같은 그림을 그릴수 있다.
ggCatepillarPlot(acs,"age","sex","Dx",interactive=TRUE)또 다른 예로 mtcars 데이타로 자동/수동변속기 및 엔진수에 따라 그래프를 그려보면 다음과 같다.
mtcars$transmission=ifelse(mtcars$am,"manual","automatic")
ggCatepillarPlot(mtcars,"mpg","transmission","cyl",interactive=TRUE)이 글에 사용된 ggAncova.R 파일은 웹R을 통하여 공개할 예정이다.