12-13장: 일원분산분석(One-way ANOVA) & Kurskal-Wallis Test

[안재형]

##### One-way ANOVA

> PlantGrowth = read.csv("PlantGrowth.csv")
> with(PlantGrowth, tapply(weight, group, mean))
 ctrl  trt1  trt2
5.032 4.661 5.526
> with(PlantGrowth, tapply(weight,group,sd))
     ctrl      trt1      trt2
0.5830914 0.7936757 0.4425733

# boxplot

> boxplot(weight ~ group, col="red", data=PlantGrowth)

 # fit one-way anova

> out = lm(weight ~ group, data=PlantGrowth)
> anova(out)
Analysis of Variance Table

Response: weight
          Df  Sum Sq Mean Sq F value  Pr(>F) 
group      2  3.7663  1.8832  4.8461 0.01591 *
Residuals 27 10.4921  0.3886                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# 회귀진단

> par(mfrow=c(2,2))

> plot(out)

> shapiro.test(resid(out))

# Install "multcomp" library and call the library

install.packages("multcomp")

library(multcomp)

### Dunnett

> out = lm(weight ~ group, data=PlantGrowth)
> dunnett = glht(out,linfct=mcp(group="Dunnett"))
> summary(dunnett)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = weight ~ group, data = PlantGrowth)

Linear Hypotheses:
                 Estimate Std. Error t value Pr(>|t|)
trt1 - ctrl == 0  -0.3710     0.2788  -1.331   0.323
trt2 - ctrl == 0   0.4940     0.2788   1.772    0.153
(Adjusted p values reported -- single-step method)

> plot(dunnett)

### Tukey

> tukey = glht(out,linfct=mcp(group="Tukey"))
> summary(tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = weight ~ group, data = PlantGrowth)

Linear Hypotheses:
                 Estimate Std. Error t value Pr(>|t|) 
trt1 - ctrl == 0  -0.3710     0.2788  -1.331  0.3908 
trt2 - ctrl == 0   0.4940     0.2788   1.772   0.1979 
trt2 - trt1 == 0   0.8650     0.2788   3.103   0.0119 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

> plot(tukey) 

> out = lm(weight ~ group, data=PlantGrowth)
> weight = rnorm(n=nrow(PlantGrowth),mean=fitted(out),sd=summary(out)$sigma)

##### Kurskal-Wallis Test

> kruskal.test(weight~group, data=PlantGrowth)

        Kruskal-Wallis rank sum test

data:  weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

[아사달]

음.. 이해가 잘 안돼요. 복사해서 붙여봐도 에러메세지가...
> dunnett=glht(out,linfct=mcp(group='Dunnett'))
Error in mcp2matrix(model, linfct = linfct) :
Variable(s) ‘group’ have been specified in ‘linfct’ but cannot be found i$

[안재형]

저는 잘 되는데요. 이걸 한꺼번에 돌려도 안되나요?
install.packages("multcomp")
library(multcomp)
out = lm(weight ~ group, data=PlantGrowth)
dunnett = glht(out,linfct=mcp(group="Dunnett"))

[열정]

install.packages("multcomp") 부분은 RStudio나 R 버전에 따라서 실행이 안되기도 하네요

[DrAhn]

저도 책에서 install.packages("multcomp") 보고 입력했더니
Warning message:
package ‘mulcomp’ is not available (for R version 3.1.0)
라고 뜨네요.

그래서 바로 위에 써주신

install.packages("multcomp")
library(multcomp)
out = lm(weight ~ group, data=PlantGrowth)
dunnett = glht(out,linfct=mcp(group="Dunnett"))

복사하고 엔터 치니깐 되네요. 왜그런거죠?

[안재형]

이상하네요. 저는 3.1.0에서 해봤는데 잘 되는데요.
RStudio 버전은 상관없을텐데, R은 최신버전에서 해야할겁니다.

[NiKe]

안선생님 질문있습니다!
13장에 Kruskal-Wallis test 가 one-way ANOVA 의 비모수 방법이라고 하셨는데, 그럼 one-way ANOVA의 회귀진단으로부터 잔차의 가정이 만족되지 않는경우에 Kruskal-Wallis test를 쓰는게 아닌가요? 12장에서 PlantGrowth 데이터를 이용해서 weight~group로 lm하고 plot(out)과 shapiro.test(resid(out))으로 잔차의 가정을 체크했을 때 모두 만족되었었는데, 13장에서 왜 weight~group 으로 Kruskal-Wallis test를 하는 건가 해서요...답변 주시면 감사하겠습니다!

[안재형]

Kruskal-Wallis test에 적합한 데이터가 없어서 그냥 one-way anova에 썼던 데이터를 사용했습니다^^

[NiKe]

아하. 그랬군요. ^^ 그럼 원래대로라면 one-way ANOVA의 회귀진단으로 부터 가정이 만족되지 않았을 때 Kruskal-wallis test를 쓰는게 맞는거군요. ㅎㅎ 답변 감사합니다! 이제 책의 끝을 향해 달려가고 있습니다. 이론적으로만 배웠던 것들을 직접 응용하면서 해보니 정리가 잘 되는거 같아요!! XD

[안재형]

책에 설명을 좀 했었어야했는데요^^ 2판이 나가게 되면 적절한 데이터를 찾던지, 그런 얘기를 추가하겠습니다. 감사합니다~

[곽대현]

안녕하세요,
12.5에서
dunnett = glht(out,linfct=mcp(group="Dunnett"))
명령어를 실행하면

Error in glht(out, linfct = mcp(group = "Dunnett")) :
could not find function "glht"

Error in mcp2matrix(model, linfct = linfct) :
Variable(s) ‘group’ of class ‘character’ is/are not contained as a factor in ‘model’.

라고 오류가 나오는데요, 왜그런지 알 수 있을까요?
해결법도 궁금합니다.

[안재형]

install.packages("multcomp")
설치하셨나요?