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

[안재형]

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

[NiKe]

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

[안재형]

NiKe 책에 설명을 좀 했었어야했는데요^^ 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")
설치하셨나요?