Introductory statistics with R: chapter 7

---

* ANOVA(Anaylsis of variance, 분산분석) and the Kruskal-Wallis test

 

* One-way ANOVA

두 집단 이상의 평균값이 차이를 보이는지 확인할 때 사용하는 분석방법으로, 개념적으로는 "집단 내 분산"과 "집단 간 분산"의 비율을 이용해 차이를 확인하는 것이며 F 분포를 이용한다. Null hypothesis는 "집단 간 차이가 존재하지 않는다" 이며, 따라서 p-value가 유의수준 이하인 경우 "어느 집단 간에는 차이가 존재한다"는 의미가 된다.

 

i개의 집단 간 비교를 시행하는 경우를 생각해보자. Xij는 i번째 집단의 j번째 관측치를 말한다. Xij 값은, 전체 평균 X bar에 집단 간 편차를 더하고, 거기에 i집단의 평균과 Xj와의 편차를 더한 값으로 표현할 수 있다. (즉, 전체 평균값에서 i 집단과의 편차를 더해 집단평균으로 범위를 좁히고, i집단 평균에서 개별 component와의 편차를 더해 개별 요소로 범위를 좁히는 개념이다)

그럼 위와 같이 정리가 된다. Mu는 전체 집단의 평균이며, alpha는 집단 간 편차, epsilon은 집단 내 편차이며 오차항이라고 한다 (error). 오차항은 독립적인 확률변수로, 평균이 0, 분산이 sigma^2인 정규분포를 따르며 alpha i값은, H0가 맞다면 0이 되어야 한다. 

그룹 내, 그룹 간, 그리고 전체 SSD(sum of squares)를 구해보자. 편차의 제곱이므로 분산과 유사한 개념이다.

위를 보면, 즉 전체 그룹의 sum of square는 그룹 간 SSD와 그룹 내 SSD의 합으로 표현된다. 또한 개념적으로, 집단 간 편차가 집단 내 편차(=랜덤)보다 유의미하게 커지게 된다면 이는 집단 간 유의미한 차이가 있음을 의미한다. 이 때 sum of square는 편차 제곱을 더한 값인데, 편차 제곱은 항상 양수이므로 집단 내 요소의 증가에 따라 무조건 증가하게 되어 설명력이 떨어진다 ("irrelevant grouping will always explain some part of the variation"). 그러므로 자유도(df, degree of freedom)으로 나눠줘서 이를 보정할 수 있음. 이 값이 바로 mean square(MS) 값이다.

** (n-1)S^2/sigma^2 ~ chisq(n-1) 인 것을 상기해보면, 표본 분산이 카이제곱 분포를 따르듯 MS 값 역시 카이제곱 분포를 따른다는 것을 직관적으로 이해할 수 있음.

 

전체 샘플의 크기가 N, 그룹의 개수가 k개 라고 한다면 그룹 '간' 편차에 대해서 df는 k-1, 그룹 '내' 편차에 대해서는 N-k, 전체에 대해서는 N-1이 된다. 그룹'간' 편차와 그룹'내' 편차의 df 값을 합하면 전체 df값이 된다.

MSw는 각 그룹 별 분산을 조합한 것이므로 전체 분산 sigma^2이 될 것이고, MSb의 분산은 H0가 참이라면 sigma^2이 될 것이 예상되지만, group effect(즉, 그룹에 따른 유의미한 편차)가 있다면 더 커질 것이다. 이 MSb와 MSw는 모두 카이제곱 분포를 따르므로, 카이제곱의 비율은 F분포를 따른다. 따라서 이를 통해 p-value를 구하는 것이 가능하다.

만약 그룹의 개수가 m이며, 각 그룹 별로 크기가 n으로 동일하다면, 위와 같이 적을 수 있다. 이 때 전체 샘플의 크기 N=mn 이며, 그룹 간 df는 m-1, 그룹 내 df는 m(n-1)이다. 그룹 내 df에 대해 조금 더 설명하면, 관측된 그룹 평균이 존재하므로 각 그룹의 df는 n-1이 된다. 이게 m개의 그룹에 대해 존재하므로 위와 같아짐.

F값은 그룹 간 편차가 없다면 ideally 1이 된다. 그리고 F test는 단측검정임을 기억하자.

 

* ANOVA의 세 가지 중요한 전제조건

1) 표본의 독립성: 각 표본들은 독립적으로 추출됨.

2) 정규성 가정: 각각의 그룹들이 추출된 모집단은 정규분포를 따름.

3) 등분산 가정: 각 그룹이 추출된 모집단의 분산은 동일하다.

  - var.test() 혹은 Bartlett's test로 확인이 가능하다.

 

* Bartlett's test: 등분산 가정을 검정

bartlett.test()를 이용해 모든 그룹 간 등분산가정을 만족하는지를 테스트할 수 있다. H0는 "모든 그룹 간 분산은 동일하다" 이므로 p-value > 0.05일 때, 이를 기각하지 못하므로 등분산을 가정한다.

 

anova()

R에서, one-way ANOVA는 regression analysis의 변형으로 존재한다. 따라서 lm()을 사용해 변수를 만들고, anova() 함수를 써서 분석을 한다.

R에서는 위의 분산분석표처럼 within, between group 으로 나타내는 것이 아니라 grouping factor(여기선 ventilation)와 잔차(residual)로 나타내게 된다. Grouping factor가 곧 between 이며, residual이 within을 의미한다.

 

주의사항

위의 분석표를 보고 p-value가 작으니 유의미한 결과가 나왔다고 해선 안된다!! 위의 tanner를 보면 df=1 이다. 실제 데이터에서는 tanner는 총 5개의 범주로 구성될 뿐 아니라, ANOVA는 보통 3그룹 이상의 비교에서 적용하므로 df=1일 수 없다. 이 경우는 tanner가 factor가 아닌 연속형 변수로 취급된 케이스이다. 따라서 이를 factor로 변환해서 test를 다시 돌려야 한다.

 

* Pairwise comparison and multiple testing

ANOVA를 통해 H0을 기각했다면, 이제 어느 그룹 간 차이가 존재하는지 확인해야 한다. 앞서 R에서는 ANOVA를 선형회귀의 변형으로 구현한다고 언급했다. 따라서 linear model의 summary에서 회귀계수를 확인해 단서를 찾아볼 수 있다.

여기서 절편(intercept)는 첫 번째 그룹을 의미하며, 다른 두 그룹이 2, 3번째 그룹과의 차이를 의미한다. 위 table에서 2, 3번째의 p-value는 바로 첫번째 그룹과 2, 3번째 그룹 간 차이를 검증한 p-value이다. 즉, 1-2번째 그룹 간 차이가 있다고 볼 수 있다(p<0.05). 하지만 이 경우 2-3번째 그룹 간의 차이는 알 수 없을 뿐더러 그룹이 더 많아진다면 더더욱 비교가 힘들다. 따라서 모든 그룹을 제대로 비교하기 위해선 multiple testing을 시행해야 한다.

 

Multiple testing

ANOVA로 그룹 간 차이가 있음이 확인된 이후, 차이가 나는 그룹을 확인하기 위해 여러 그룹을 서로 비교해보는 test로, 다중비교의 문제가 발생해 p-value exaggeration이 있다. 이를 간단히 설명하자면, p-value 0.05란 영가설이 참일 확률이 5% 이내라는 의미이다. 한 번의 test만을 놓고 보면 작은 확률처럼 보이지만, 100개의 test를 시행한다면 5번 정도는 잘못된 결론을 내리게 될 수 있다는 의미이기도 하다. 따라서 다중 검정을 시행하려면 p-value를 이에 맞게 보정해줘야 한다. 가장 잘 알려진 보정방법으로는 Bonferroni correction이 있으며, 이는 cutoff가 되는 p-value를 그룹 간 총 비교의 개수로 나눠주는 방법이다.

 

pairwise.t.test()

모든 종류의 two-group comparison을 해주는 함수로, adjustment를 해주면 다중 그룹 비교도 가능하다.

이 때 제시되는 p-value들은 이미 보정된 값들이다. 즉, group 간 비교의 개수가 3개라면 p-value에 3을 곱해 보정한 상태로 주워지는 것이다. 따라서 우리가 상정한 p-value(예를 들면 0.05)를 그냥 cut-off로 사용하면 된다.

 

p.adj를 따로 설정해주지 않는다면 default는 variant due to Holm 이라는 방법을 사용한다. 이 방법은 p-value의 값에 따라 서로 다른 보정치를 곱해준다. (가장 작은 값은 전체 비교의 개수 n, 두번째 작은 값은 n-1, 그 다음은 n-2, ....)

 

* Oneway test: 등분산 가정을 만족시키지 못한 경우

ANOVA는 비교하는 모든 그룹이 등분산일 것을 전제로 한다 (Bartlett test로 이를 확인한다). 이를 만족시키지 못했다면, 등분산 가정을 요구하지 않는 다른 test를 시행해야 할 것이다. 

 

oneway.test()

또한 사후검정으로 사용하는 pairwise.t.test 역시 argument를 통해 등분산가정이 F임을 설정해줘야 한다.

 

* Graphical presentation

본문 예제에서 stripchart()를 이용해 그룹 별 분포를 시각적으로 보여준다.

코드를 살펴보자.

tapply()란, 그룹별 처리를 위한 함수로 tapply(데이터, index(=factor), function) 의 형식으로 입력하면 factor에 따라 data에 function을 적용해주는 함수이다. 즉, xbar와 s, n은 모두 ventilation 그룹에 따라 mean, sd, length를 구한 값들이 된다. 그리고 sem은 s를 sqrt(n)으로 나눈, 표준오차값이 된다.

stripchart()를 이용해서 folate를 ventilation으로 describtion하되, jitter를 이용해 한줄이 아닌 흩어진 모양으로 출력했다. pch는 점의 종류를 말하는데, 검은색 dot을 말하며, vert=T를 통해 vertical한 stripchart를 만들었다.

arrow()로 화살표를 그렸는데, arrow() function의 형식은 arrow(x1,y1,x2,y2)이다. 위의 코드를 보면, 그룹별로 x좌표는 1, 2, 3 이렇게 할당되므로 1:3은 각 stripchart 그룹의 x좌표를 말한다고 볼 수 있다. 그리고 y좌표는 xbar, 즉 각 그룹의 평균값에서 각 그룹의 sem을 더하고 뺀 값을 각각 설정함으로써 수직으로 된 화살표를 그렸다. angle=90은 arrow의 화살 끝과 선이 이루는 각도로, 90도를 해놓으면 수직이 된다. code=3는 양측 화살임을 의미한다. 마지막으로 lines() 함수로 선을 그리는데, x좌표는 1:3, y좌표는 평균값이 되며 pch=4는 x 모양의 dot을 의미한다. type=b 는 점과 점을 연결하는 형태의 선을 의미하며 cex는 크기를 의미한다. 위 코드의 결과값을 보자.

위 그래프는 SEM값을 arrow()로 범위를 표시했지만, 이는 mean값이 모평균에 비해 어느 정도의 오차를 지니는지를 나타내주는 것이므로 때에 따라서는 confidence interval이나 SD가 더 적절할 수도 있다.

하지만 많은 경우 SEM값을 사용하는 이유는 여러 범주값 중 가장 작기 때문이다: 인용 - in many fields it appears to have become the tradition to use 1 SEM “because they are the smallest”; that is, it makes differences look more
dramatic.)

 

* Kruskal-Wallis test

One-way ANOVA의 nonparametic counterpart(즉, 정규성 가정 없이 사용할 수 있는 검정)이 바로 Kruskal-Wallis test 이다. 다른 비모수 검정들이 그렇듯이, 역시 group에 따라 rank를 매기고 그룹 간 rank의 sum of squares를 비교하는 과정을 거친다. 아래와 같이 kruskal.test()를 시행할 수 있다.

다른 비모수검정들과 마찬가지로, 비모수검정은 less-efficient하므로 ANOVA를 했을 때와는 달리 insignificant p-value가 계산되었다.

 

* Two-way ANOVA

하나의 factor에 대해서만 그룹을 나눠 분석한 것이 one-way ANOVA였다면, 여러 기준에 따라 cross-classified된 자료에 대해서도 ANOVA를 시행할 수 있다. 

이번에도 Xij에 대해 가정하자. One-way ANOVA때와는 달리, 관측값은 factor 당 한개씩이며 Xij는 m x n table의 i, j 번째 값을 의미한다. 그리고 row, column의 평균과 각 cell의 값에 대해 sum of square를 구할 수 있다. 이를 각각 SSDr, SSDc 라고 하고, 각 cell의 값에서 row, column에 따른 변동을 제거한 값의 sum of square를 SSDres라고 한다.

Xij에 대해 통계학적 모델을 세워보면, Xij는 전체평균값에 alpha와 beta 라는 요소의 영향을 더한 값 + 오차항으로 구성되며, 오차항은 독립변수이며 N(0,sigma^2)를 따른다. 

one-way ANOVA와 마찬가지로 df로 SSDr, SSDc, SSDres를 나누어 mean square를 구할 수 있으며 이 때 df는 n x m table에 대해 (n-1), (m-1), (m-1)(n-1)이 된다. 

F test는 SSDr, SSDc를 각각 SSDres로 나눈 값을 F값으로 하여 시행할 수 있다.

 

Two-way ANOVA에서 중요한 건 자료가 잘 'balanced'되어있어야 한다는 것이다: 마치 chi-square test를 시행하는 조건이 기대빈도 5 미만의 셀이 전체의 20%를 넘지 않아야 하는 전제조건이 있는 것 처럼?

 

본문에서 ISwR package의 heart.rate 이라는 dataset의 예제를 들어 설명하고 있음.

 

heart.rate 이라는 dataset의 외형은 아래와 같다

이를 만들기 위한 gl() 이라는 함수를 소개하고 있음. (데이터는 후략, 36번째 자료까지 있음)

 

gl()

generate level 함수로, 일정 규칙에 따라 factor를 부여해줌. 위의 데이터를 보면 subj는 1,2,3,4,5,6,7,8,9 가 반복되는 형태로 입력이 되어 있고, time은 0, 30, 60, 90이 각각 9개씩 있음. rep() 함수로 표시하자면 rep(c(1:9), 4)와 rep(c(0,30,60,120), c(9,9,9,9))와 같을 것이다. 이를 gl() 함수를 이용해 rep과 마찬가지로, 단 factor type으로 지정해서 생성할 수 있음.

 

gl(n,k,length,labels,ordered)

n: number of levels

k: number of replication(반복수)

length: length of the result(result의 총 길이)

labels: factor 이름을 지정 가능

ordered: a logical indicating whether the result should be ordered or not.

 

subj의 코드는 다음과 같다: subj <- gl(9,1,36)

time의 코드는 다음과 같다: time <- gl(4,9,36,labels=c(0,30,60,120))

 

factor들이 갖춰졌다면, two-way ANOVA 자체는 simple하게 실행 가능하다.

책에 따르면, 이 자료는 balanced design이므로 lm(hr~time+subj)와 같이 factor의 순서를 바꿔서 넣어도 정확히 같은 결과가 출력된다고 함. Unbalanced data에서는 factor 넣는 순서에 따라 값이 달라질 수 있다고 함.

 

* Graphics for repeated measures

위의 two-way ANOVA를 시행하는 경우, 각 factor에 대해 data가 어떻게 분포하는지 "spaghettigram"(면발을 늘어놓은 것 처럼 factor에 따른 값을 여러 직선으로 표시한 그림을 의미하는 듯)를 그려보는 것인 자료 이해에 도움이 된다. 이 그림은 interaction.plot()함수를 이용해 그릴 수 있다.

interaction.plot(x.factor, trace.factor, response, fun, ...)

x.factor: x축을 구성하는 factor 1

trace.factor: 스파게티 면발을 구성하는 factor 2

response: y축 데이터 (x, trace factor로 설명하려는 데이터)

fun: default는 mean으로, 위 데이터에서는 편의를 위해 각 factor 1, factor 2 에 해당하는 데이터가 한개 씩이었으나 데이터가 실제로 여러 개 존재하는 경우 어떤 식으로든 계산을 해서 하나의 값으로 만든 다음 그림을 그려야 한다. 이 때 대표값을 추출하는 함수를 fun argument에 입력함.

 

* The Friedman test

Two-way ANOVA에 대응하는 비모수적인 방법은 Friedman test이다. 이 때, 조건은 각 cell 별로 1개의 관측치만 존재해야 한다는 것. Friedman test는 각 row 내에서 관측치의 ranking을 매기고, column effect가 없다면 같은 row 내 분포는 유사할 것을 가정해서 검정을 진행한다.  Column sum of squares를 계산해 표준화하면 chi-square distributed test statistics를 얻을 수 있고, 이를 카이제곱 분포에 대입해 p-value를 얻어낸다.

* The ANOVA table in regression analysis

위에서 one-way, two-way ANOVA의 경우 variance table을 이용해 검정을 진행했다. 그리고 이 경우, linear model 함수를 이용해서 변수들을 연관지었다. 이 사실에서 엿볼 수 있듯이, variance table은 ANOVA 뿐만 아니라 모든 linear model에 대해 적용이 가능하다 (자세한건 12 챕터에서 다룬다고 함).

앞에 chapter6에서 thuesen에 대해 linear model을 만든 것을 생각해보자. 위는 thuesen의 linear model에 대해 one way AVOVA를 진행한 결과이다. 위의 ANOVA 결과와 lm(short.velocity~blood.glucose)의 summary한 결과를 비교해보면, 동일한 결과임을 확인할 수 있다. ANOVA에서 factor의 F 및 p-value는 lm summary의 F 및 p-value와 동일하며, residual standard error (0.2167)은 ANOVA의 residual mean Sq의 루트 값이다 (sqrt(0.04696)). R-squared는 SSDtotal에 대해 SSDb의 비율(0.1737 = 0.20727 / (0.20727+0.98610)), adjusted R^2는 다음과 같다: the adjusted R2 is the relative improvement of the residual variance, 0.1343 = (v − 0.04696)/v, where v = (0.2073 + 0.9861)/22 = 0.05425 is the variance of short.velocity if the glucose values are not taken into account.

 

(v값은 (SSD total)/(df total)이며, Adjusted R^2 = (v- mean sq of residuals)/v 가 된다. )