교호 작용이 있는 모형에서 다중 비교

R에서 multcomp 패키지의 glht를 이용하면 다중비교(multiple comparison)가 가능하다. 더 자세한 설명은 multcomp 매뉴얼을 참조하자.

• http://cran.r-project.org/web/packages/multcomp/index.html

기본

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)