R에서 3개 이상 집단의 분산을 비교하는 방법(+사후검정)
세개 이상 집단의 분산을 비교하는 검정들은 levene's test 와 bartlett's test 등이 있다.
이들은 모든 집단의 분산이 같다는 귀무가설에 대한 p 값만 알려줄 뿐, 어느 집단 간에 차이가 있었는지는 알려주지 않습니다.
어느 집단 간에 차이가 있는지 사후검정을 해야하는데요. levene's test와 bartlett's test에 적용할 수 있는 사후검정 함수를 R에서는 찾지 못했습니다.
다른 방법을 통해 해결했는데요. 분산이 결국 편차제곱의 평균이므로, 각 원소의 편차제곱을 계산하고 이들의 평균비교를 해주면 됩니다.
아래 예시를 통해 알아봅시다.
a,b,c 세 그룹을 생성했습니다.
set.seed(1)
a=sample(0:100,30,replace=TRUE)
b=sample(10:90,30,replace=TRUE)
c=sample(20:80,30,replace=TRUE)
> a
[1] 26 37 57 91 20 90 95 66 63 6 20 17 69 38 77 50 72
[18] 100 38 78 94 21 65 12 26 38 1 38 87 34
> b
[1] 49 58 49 25 77 64 74 18 68 43 76 62 73 54 52 73 11 48 69 66 48 79
[23] 45 29 15 18 35 52 63 42
> c
[1] 75 37 48 40 59 35 49 66 25 73 40 71 41 40 49 74 72 43 67 78 46 63
[23] 44 39 66 32 63 27 34 28
각 그룹의 평균과 분산을 계산합니다.
m_a=mean(a)
m_b=mean(b)
m_c=mean(c)
각 그룹 원소의 편차제곱을 계산합니다.
ds_a=(a-m_a)^2
ds_b=(b-m_b)^2
ds_c=(c-m_c)^2
상자수염그림을 한번 그려봅시다.
> boxplot(ds_a,ds_b,ds_c,names=c("ds_a","ds_b","ds_c"))
이제 각 그룹에 분산분석과 사후분석을 적용하면 됩니다.
먼저 편차제곱들을 아래와 같이 1열은 그룹이름, 2열은 값 형태로 데이터프레임에 저장합니다.
x=c(rep("ds_a",length(ds_a)),
rep("ds_b",length(ds_b)),
rep("ds_c",length(ds_c)))
y=c(ds_a,ds_b,ds_c)
data=data.frame(x,y)
> head(data,5)
x y
1 ds_a 232.05444
2 ds_a 17.92111
3 ds_a 248.58778
4 ds_a 2476.72111
5 ds_a 450.85444
분산분석을 적용합시다.
> maov=aov(y~x,data)
> maov
Call:
aov(formula = y ~ x, data = data)
Terms:
x Residuals
Sum of Squares 5455403 41439423
Deg. of Freedom 2 87
Residual standard error: 690.1559
Estimated effects may be unbalanced
summary 함수를 통해 p값을 알아탭시다.
> summary(maov)
Df Sum Sq Mean Sq F value Pr(>F)
x 2 5455403 2727701 5.727 0.00461 **
Residuals 87 41439423 476315
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
p값이 0.004 이므로, 세 집단의 분산이 동일하다는 귀무가설이 기각됩니다.
이번에는 사후검정을 해봅시다. 본페로니 검정을 하겠습니다.
> my_bf=LSD.test(maov,"x",p.adj="bonferroni",group=F)
> my_bf
$statistics
MSerror Df Mean CV t.value MSD
476315.2 87 610.053 113.1305 2.44115 435.0069
$parameters
test p.ajusted name.t ntr alpha
Fisher-LSD bonferroni x 3 0.05
$means
y std r LCL UCL Min Max
ds_a 953.1167 996.9635 30 702.6688 1203.5646 10.4544444 3453.521
ds_b 387.0056 434.5934 30 136.5577 637.4535 0.6944444 1613.361
ds_c 490.0367 496.1229 30 239.5888 740.4846 0.0100000 1672.810
Q25 Q50 Q75
ds_a 232.05444 576.0544 1333.6544
ds_b 40.19444 240.6944 592.8611
ds_c 48.31000 358.2100 710.4100
$comparison
difference pvalue signif. LCL UCL
ds_a - ds_b 566.1111 0.0062 ** 131.10425 1001.1180
ds_a - ds_c 463.0800 0.0330 * 28.07314 898.0869
ds_b - ds_c -103.0311 1.0000 -538.03797 331.9757
$groups
NULL
attr(,"class")
[1] "group"
$comparison 을 보시면 a그룹과 b그룹, a그룹과 c그룹의 분산 유의차가 있고, b그룹과 c그룹의 분산 유의차는 없다는 것을 알 수 있습니다.
댓글