R에서
multcomp 패키지의
glht를 이용하면 다중비교(multiple comparison)가 가능하다. 더 자세한 설명은
multcomp 매뉴얼을 참조하자.
기본
R에 기본으로 들어있는
warpbreaks 데이터세터를 예제로 한다. 이 데이터는 다음과 같이 연속형 반응변수
breaks와 범주형 설명변수
wool,
tension으로 구성되어있다.
> summary(warpbreaks)
breaks wool tension
Min. :10.00 A:27 L:18
1st Qu.:18.25 B:27 M:18
Median :26.00 H:18
Mean :28.15
3rd Qu.:34.00
Max. :70.00
다음과 같은 선형모형을 적합하였다고 하자.
> mod <- lm(breaks ~ wool + tension, data=warpbreaks)
그리고 나서
tension의 세 수준 L, M, H 사이의 다중 비교를 할 필요가 있다고 하자. 물론 우선
multcomp패키지를 로딩해야 하고
glht() (General Linear Hypothesis Test)를 이용한다.
> require(multcomp)
> mod.mc <- glht(mod, linfct=mcp(tension="Tukey"))
결과는 다음과 같이
summary() 함수를 이용해 볼 수 있다.
plot(mod.mc)으로 그래프를 그려 볼 수도 있다.
> summary(mod.mc)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = breaks ~ wool + tension, data = warpbreaks)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
M - L == 0 -10.000 3.872 -2.582 0.0336 *
H - L == 0 -14.722 3.872 -3.802 0.0011 **
H - M == 0 -4.722 3.872 -1.219 0.4474
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
교호작용이 있는 경우 다중비교
교호작용(interaction)이 있는 모형이 필요하다고 하자. 예를 들면 다음과 같다.
> mod1 <- lm(breaks ~ wool * tension, data=warpbreaks)
이 경우에 모든 가능한 조합은 6가지가 있다. 다음과 같이 6개의 셀을 평균을 비교하는 것이다.
> xtabs(~wool+tension, warpbreaks)
tension
wool L M H
A 9 9 9
B 9 9 9
이때 우리가 다중비교를 원하는 대상이
wool의 동일한 수준 내에서
tension을 비교하는 것이라고 하자. 우선 모형을 만들때 다음과 같이
interaction() 함수를 이용해 새로운 변수를 만들고 절편은 없이 하는 것이 생각하기 편리하다.
tw <- with(warpbreaks, interaction(tension, wool))
mod1 <- lm(breaks ~ tw-1, warpbreaks)
이렇게 했을 때 회귀계수는 각 셀의 평균이 된다. 다음과 같다.
> coef(mod1)
twL.A twM.A twH.A twL.B twM.B twH.B
44.55556 24.00000 24.55556 28.22222 28.77778 18.77778
이것과 순서를 맞추어 다음처럼 contrast matrix를 만들어주자.
contr <- rbind("A:M-L" = c(-1,1,0,0,0,0),
"A:H-L" = c(-1,0,1,0,0,0),
"A:H-M" = c(0,-1,1,0,0,0),
"B:M-L" = c(0,0,0,-1,1,0),
"B:H-L" = c(0,0,0,-1,0,1),
"B:H-M" = c(0,0,0,0,-1,1))
첫행을 보면 "A:M-L"이라고 이름을 붙였는데
wool=A에서
tension=M과
tension=L을 비교하겠다는 것이다. 위에 적합한 모형의 회귀계수 중 두번째 것에서 첫번째 것을 빼준 것이 이에 해당한다. 따라서
c(-1,1,0,0,0,0)이 된다. 이제 다음과 같은 결과를 얻을 수 있다.
> summary(glht(mod1, linfct=contr))
Simultaneous Tests for General Linear Hypotheses
Fit: lm(formula = breaks ~ tw - 1, data = warpbreaks)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
A:M-L == 0 -20.5556 5.1573 -3.986 0.00133 **
A:H-L == 0 -20.0000 5.1573 -3.878 0.00185 **
A:H-M == 0 0.5556 5.1573 0.108 0.99996
B:M-L == 0 0.5556 5.1573 0.108 0.99996
B:H-L == 0 -9.4444 5.1573 -1.831 0.30801
B:H-M == 0 -10.0000 5.1573 -1.939 0.25536
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
