일원배치 분산분석의 원리

❗️블로그 옮김:  https://www.taemobang.com

방태모

분산분석(ANOVA : analysis of variance)의 목적은 세 집단 이상의 모평균 비교이다. 단, 집단 간 평균의 대소 관계는 알 수 없다. 이를 위해서는 추가적인 다중 비교(multiple comparison)가 필요하다. 요인이 1개인 CRD(Competely Randomized Design : 완전랜덤화설계)에서 쓰이는 분산분석 방법을 일원배치 분산분석(One-way ANOVA)이라고 한다. 일원배치 분산분석과 CRD 두 용어를 혼용하여 쓸 때가 있는데, 정확하게 정의할 필요가 있다.

 

CRD는 실험단위들의 각 처리 그룹 배치와 실험 순서를 랜덤으로 정하여 실험을 설계하는 경우를 말한다(다만, 실험에 따라서 순서가 무의미하면 동시에 실험을 진행). $t$개의 처리그룹마다, $r$번씩 반복하는(즉 실험단위들의 총 갯수 $N = rt$) CRD의 자료구조는 다음과 같다(각 처리마다 반복수가 동일한 균형자료, 그리고 고정효과일때만 고려한다).

 

1	$\cdots$	$t$
$Y_{11}$	$\cdots$	$Y_{t1}$
$\vdots$	$\vdots$	$\cdots$
$Y_{1r}$	$\cdots$	$Y_{tr}$

 

위 자료구조를 세 가지 형태의 모형식으로 표현할 수 있다(모형식 마지막에 $i = 1, \cdots, t$; $j = 1,  \cdots, r$은 꼭 써줘야 하지만, 지금은 편의상 생략).

 

• $Y_{ij} \sim N\left (\mu_i, \sigma^2   \right )$   $\cdots$ (식 1)
• $Y_{ij} = \mu_i + \epsilon_{ij}$,  $\epsilon_{ij} \sim N\left (0, \sigma^2   \right )$   $\cdots$ (식 2)
• $Y_{ij} = \mu + \tau_i +\epsilon_{ij}$,  $\epsilon_{ij} \sim N\left (0, \sigma^2   \right )$   $\cdots$ (식 3)

세 번째 모형식이 가장 보편적인 일원배치 분산분석의 모형식이다. 여기서 $\tau_i$는 $i$번째 처리 효과로 다음과 같이 정의된다:

 

$\tau_i = \mu_i - \mu$

 

$\tau_i$들에 대해 summation을 취하면 다음이 항상 성립한다:

 

$\sum \tau_i = \sum \mu_i - t\mu = 0$

 

각 처리 그룹의 평균 $\mu_i$에서 총평균(Grand mean) $\mu$를 뺀 값의 합이 0이 됨은 당연하다. 이를 고정효과(Fixed effect)를 가지는 일원배치 분산분석의 Side condition이라 한다.

 

 (1) 필요한 가정

 독립성(각 처리 그룹은 서로 독립), 정규성(정규 모집단으로 부터 관측되었다고 가정), 등분산성(각 처리 그룹의 분산은 동일)으로 총 3가지 가정이 필요하다. 먼저, 각 처리 그룹간의 독립성이 위배되면 후에 분산분석의 원리 부분에서 언급할 정규방정식(normal equations)이 유일한(unique) 해를 갖는 것이 불가능해진다. 두번째로, 정규성을 위배하면 세 집단의 모평균이 유의한 차이가 있는지에 대한 검정 결과를 신뢰할 수 없게된다($F$-분포는 정규모집단으로부터의 표본분포이므로). 마지막으로 등분산 가정이 필요한 이유는 집단간 분산이 유의하게 다르면, 마찬가지로 집단간 평균 차이의 검정결과가 잘못된 결론으로 내려질 수 있기 때문이다. 그리고 애초에 분산의 차이가 유의하게 존재하는 집단들을 평균만으로 비교하는 것은 합리적이지 않다고 생각이든다. 분석 수행후, 가정의 확인이 필요한 이유는 과연 이 분석결과가 신뢰할만한 결과인지를 확인하기 위해서이다. 즉, 표본들에서 관측된 차이가 단순히 "운"때문인지 아니면 실제로 모집단들도 표본에서의 추정과 같이 유의미한 차이가 있는것인지 확인하는 절차라고 볼 수 있다.

 

 (2) 분산분석의 원리

 일원배치 분산분석의 가설은 2가지 형태로 쓸 수 있다.

 

$H_0 : \mu_1 = \cdots = \mu_t \; vs \; H_1 : not\;H_0$

$H_0 : \tau_1 = \cdots = \tau_t \; vs \; H_1 : not\;H_0$

 

가장 보편적인 모형식이 (식 3)임을 고려하면 두번째 가설 형태가 선호된다. 그럼 일원배치 분산분석의 모수 $\mu, \mu_i, \tau_i, \epsilon_{ij}$들에 대한 추정은 어떻게 할까? 직관을 기반으로 써보자.

 

• $\hat{\mu} = \frac{1}{N} \sum_i^t \sum_j^r Y_{ij} \equiv \bar{Y}_{\cdot \cdot }$ (Grand mean : 총평균)
• $\hat{\mu}_i = \frac{1}{r} \sum_j^r Y_{ij} \equiv \bar{Y}_{i \cdot }$ (Treatment mean : 처리평균)
• $\hat{\tau}_i = \hat{\mu}_i - \hat{\mu} =  \bar{Y}_{i \cdot } - \bar{Y}_{\cdot \cdot}$
• $\hat{\epsilon}_{ij} = Y_{ij} - \bar{Y}_{i \cdot}$

이렇게 직관을 기반으로 쓴 모수들의 각 추정량들은 오차제곱합을 최소화하는 아이디어를 가진 LSE(Least Squares Estimation : 최소제곱법)를 이용해 모수를 추정한 결과와 동일하다. LSE를 통해 추정되는 두 모수 $\mu_i, \tau_i$는 다음의 목적함수를 최소화한다.

 

$\underset{\mu, \tau_1, \cdots, \tau_t}{arg\; min} Q = \sum_i^t \sum_j^r \epsilon_{ij}^2 =  \sum_i^t \sum_j^r \left ( Y_{ij} - \mu - \tau_i \right )^2$

 

이에 따라 해를 구해야하는 연립방정식은:

 

$\left\{\begin{matrix} \frac{\partial Q}{\partial \mu} = 0
\\ \frac{\partial Q}{\partial \tau_i} = 0\;,\; i = 1, \cdots, t

\end{matrix}\right.$

 

이를 정규방정식(normal equations)이라고 부르며, 대략적인 형태를 써보면 다음과 같다.

 

$N \hat{\mu} + r\hat{\tau}_1 + \cdots + r\hat{\tau}_t = N \bar{Y}_{\cdot \cdot} \\ 
r \hat{\mu} + r\hat{\tau}_1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;  = r \bar{Y}_{1 \cdot} \\
r \hat{\mu} \;\;\;\;\;\;\;\;\;\;\; + r\hat{\tau}_2 \;\;\;\;\;\;\;\;\; =  r \bar{Y}_{2 \cdot} \\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots  \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\
r \hat{\mu} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + r\hat{\tau}_t =  r \bar{Y}_{t \cdot} \\$

 

위 구조를 행렬 Form으로 바꿔보면, $\boldsymbol{A}_{(t+1 \times t+1)} \cdot \boldsymbol{x}_{(t+1 \times 1)} = \boldsymbol{b}_{(t+1 \times 1)}$ 와 같다. 여기서 맨 위의 방정식은 나머지 방정식들의 합으로 주어지므로, 서로 독립인 방정식이 $t$개 밖에 존재하지 않는다. 즉, 행렬 $\boldsymbol{A}$의 행과 열을 $m$과 $n$으로 나타낸다하면, 사실상 $m<n$인 상황이다. 즉 무수히 많은 해가 존재하게 된다. 그러나, 이때 일원배치 분산분석의 side conditions으로 $\sum \tau_i = 0$이 보장되므로, 기존의 맨 위의 식이 $N \hat{\mu} = N \bar{Y}_{\cdot \cdot}$로 정리된다. 그에따라 $\boldsymbol{A}$는 Full-rank 행렬이 되어 역행렬이 존재하고, 유일한 해를 보장받는 정규방정식이 된다. 어쨌든, 이렇게 얻어지는 해는 우리가 앞서 직관적으로 구해본 해와 동일하다.

 

(식 3)을 관측된 자료를 기반으로 추정하는 식으로 다시쓰면 다음과 같이 쓸 수 있다.

 

$ Y_{ij} = \mu + \tau_i +\epsilon_{ij}\;,\; \epsilon_{ij} \sim N\left (0, \sigma^2   \right ) \\
\;\;\;\; = \bar{Y}_{\cdot \cdot} + \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right ) + \left ( Y_{i j} - \bar{Y}_{i \cdot} \right )$

 

이제 분산분석에서 변동의 분해에서 봐왔던 익숙한 식으로 바꿀 수 있다.

 

$Y_{ij} - \bar{Y}_{\cdot \cdot} = \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right ) + \left ( Y_{i j} - \bar{Y}_{i \cdot} \right )$

 

차례대로 총편차, 처리효과편차, 오차편차 라고 한다. 단일 관측값 $Y_{ij}$ 하나에 대해서만 해석을 내리는 것은 무의미하므로, 관측값들이 가진 모든 정보를 합쳐서 해석하기위해 식을 한번 더 바꾸자. 양변을 제곱하고 모든 관측값들에 대해 합을 취할 것이다:

 

$\sum_i^t \sum_j^r \left ( Y_{ij} - \bar{Y}_{\cdot \cdot} \right )^2 = \sum_i^t \sum_j^r \left [ \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right ) + \left ( Y_{i j} - \bar{Y}_{i \cdot} \right )  \right ]^2 \\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum_i^t \sum_j^r \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right )^2 + \sum_i^t \sum_j^r \left ( Y_{i j} - \bar{Y}_{i \cdot} \right )^2 $

 

우변의 전체 제곱식에서 나타나는 cross-product는 0으로 사라진다. 위와 같은 관계식을 제곱합의 분할(partition of SS)라 하며, 각 제곱합 항을 순서대로 총제곱합(SST : total sum of squares), 처리제곱합(SStrt : treatments sum of squares), 오차제곱합(SSE : Error sum of squares)라고 부른다(제곱이라는 용어대신 변동(variation)이라는 단어를 사용하기도 한다). SST의 대부분을 SStr이 차지한다면 귀무가설을 기각할 수 있게되고, SSE가 차지한다면 귀무가설을 기각할 수 없게 될것이다. 실제 분산분석에서 검정통계량은 이들을 자유도로 나눈 형태인 다음의 평균제곱(MS : mean squares)을 이용한다.

 

$MStrt = \frac{1}{t-1}SStrt$,  $MSE = \frac{1}{N-t}SSe$

(자유도의 정의와 계산에 대한 이해가 부족하다면 자유도를 참고하면 된다)

 

여기서 중요한 점은, $MSE$는 오차항의 분산 $\sigma^2$의 unbiased estimator(비편향 추정량)라는 것이다. 이는 간단하게 증명된다.

 

   pf) 관측값들에 대해서는 $Y_{ij} \sim N\left ( \mu_i, \sigma^2 \right )$으로 정규모집단을 가정했으므로, 모수 $\mu_i$와 $\sigma^2$의 불편추정량은 각각 $\bar{Y}_{i \cdot}$과 $s_i^2$ 이다. 이때 $\hat{\sigma}^2 = s_i^2 = \frac{1}{r-1} \sum_j^r \left ( Y_{ij} - \bar{Y}_{i \cdot} \right )^2 $ 이라하면, 등분산 가정에 의해 모분산에 대한 $t$개의 추정량이 존재하는 것이다. 이 경우 간단한 해결책은 가중평균을 취하는 것이다:

 

$s^2 = \frac{(r-1)s_1^2 + \cdots + (r-1)s_t^2}{(r-1)t} = \frac{\sum_i^t s_i^2 (r-1)}{tr-t} = \frac{\sum_i^t \sum_j^r \left ( Y_{ij} - \bar{Y}_{i \cdot} \right )^2}{N-t} \;\; \left ( \because s_i^2 = \frac{1}{r-1}\sum_j^r \left ( Y_{ij} - \bar{Y}_{i \cdot} \right ) \right )$

$\therefore E[s^2] = E[MSE] = \sigma^2$

 

일원배치 분석분석의 검정 통계량과 검정 통계량의 귀무가설 하의 분포는 다음과 같이 주어진다.

 

$F_o = \frac{SStrt/(t-1)}{SSE/(N-t)} = \frac{MStrt}{MSE} = \frac{\sum_i \sum_j \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right )^2 / (t-1) }{\sum_i \sum_j \left ( Y_{ij} - \bar{Y}_{i \cdot}  \right )^2 / (N-t) } \underset{H_0}{\sim } F(t-1, N-t)$

 

검정통계량이 F-분포로 정의되는 이유는 뭘까? F-분포는 정규모집단에 대한 표본의 분포로, 카이제곱분포의 비로 정의된다. 카이제곱분포와 F-분포를 정의하면 다음과 같다.

 

$\frac{(n-1)s^2}{\sigma^2} = \frac{\sum (X_i - \bar{X})^2}{\sigma^2} \sim \chi^2\left ( n-1 \right ),$  where  $s^2 = \frac{1}{n-1} \sum\left ( X_i - \bar{X} \right )^2$

 

$\frac{U/r_1}{V/r_2} \sim F(r_1, r_2),$  where  $ U \sim \chi^2(r_1), \;\; V \sim \chi^2(r_2),\;\; U\perp V $

 

$\perp$는 독립임을 의미한다. 카이제곱분포와 F-분포가 정의되는 형태를 기억하고, 우리의 검정통계량을 다시 써보자.

 

$\frac{\sum_i \sum_j \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right )^2 / (t-1) }{\sum_i \sum_j \left ( Y_{ij} - \bar{Y}_{i \cdot}  \right )^2 / (N-t) } = \frac{\frac{\sum_i \sum_j \left ( \bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} \right )^2}{\sigma^2}/ (t-1) }{\frac{\sum_i \sum_j \left ( Y_{ij} - \bar{Y}_{i \cdot}  \right )^2}{\sigma^2}/(N-t) } $

 

이제 검정통계량이 왜 F-분포를 따르는지 보인다. 분자와 분모 각각 자유도가 $(t-1)$, $(N-t)$인 카이제곱분포를 따르게 되고, 결국 검정통계량은 $F(t-1, N-t)$ 분포를 따르게 된다. 또한 검정통계량은 결국 분산의 비를 비교하는 형태를 띠고 있다. 이러한 이유에서, 세 집단이상의 모평균 비교가 목적임에도 불구하고 "분산분석"이라는 이름을 갖는 것이다. 검정통계량의 분자와 분모 각각은 편차제곱을 자유도로 나눈 형태이며, 자유도로 나눈다는것은 평균을 취하는 의미와 같다. 편차제곱의 평균은 뭔가? "분산"이다. 즉, 검정통계량은 결국 분산의 비라고 할 수 있다. 이러한 이유에서, 세 집단의 모평균 비교가 목적임에도 불구하고 "분산분석"이라는 이름을 갖는 것이다.

 

(3) 정리

이때까지 배운 일원배치 분산분석 모형에 대한 것들을 간략하게 정리하자.

 

CRD Model  $Y_{ij} = \mu + \tau_i + \epsilon_{ij}, \; \epsilon_{ij} \sim N(0, \sigma^2)$
$\;\;\;\;\;(i = 1, \cdots, t;\;  j = 1, \cdots, r)$

 

$H_0 : \tau_1 = \cdots = \tau_r = 0\;\; vs \;\; H_1 : not \;\; H_0$

 

유의수준 또는 유의확률(p-value)을 통해 검정을 수행하면 되고, 유의수준을 통해 검정을 한다면 :

 

$ \textrm{if}\;\; F_0 > F(t-1, N-t; 1-\alpha),\; \textrm{we can reject}\;\; H_0$

 

단, 여기서 $1-\alpha$는 오른쪽이 아닌 그래프의 왼쪽영역을 말한다.

 

CRD에 대한 분산분석표에 대해서도 알아둘 필요가 있다.

 

Source	$SS$	$df$	$MS$	$F$
처리(Treatments)	$SStrt$	$t-1$	$MStrt$	$MStrt/MSE$
오차(Error)	$SSE$	$N-t$	$MSE$
총	$SST$	$N-1$

 

(4) 신뢰구간

 

 ○ 처리 평균 $\mu_i$에 대한 $100(1-\alpha)$% 신뢰구간

 신뢰 구간의 추정을 위해, 먼저 $\bar{Y}_{i \cdot}$이 무슨 분포를 따르는지 써보자. 각 처리 그룹의 평균은 $\bar{Y}_{i \cdot} \sim N(\mu_i, \frac{\sigma^2}{r})$과 같은 분포를 따른다. 또한, $MSE$는 $\sigma^2$의 불편추정량임을 앞서 보였으므로, MSE를 $\sigma^2$의 추정량으로 사용할 수 있다. 즉, 처리 평균 $\mu_i$에 대한 $100(1-\alpha)$% 신뢰구간은 다음과 같이 주어진다(모분산을 모르는 상황이므로, t-분포를 이용해 쓴다).

 

$\bar{Y}_{i \cdot} \pm t_(N-t; 1-\alpha/2)\frac{\sqrt{MSE}}{\sqrt{r}}$

 

 ○ 처리 평균의 차 $(\mu_i - \mu_j)$에 대한 $100(1-\alpha)$% 신뢰구간 
 위의 연장선상이다. 똑같이 먼저 $\bar{Y}_{i \cdot} - \bar{Y}_{j \cdot}$가 무슨 분포를 따르는지 써보자. 처리 평균의 차는 $\bar{Y}_{i \cdot} - \bar{Y}_{j \cdot} \sim N(\mu_i - \mu_j, \frac{2}{r}\sigma^2)$과 같은 분포를 따른다. 즉, 처리평균의 차 $(\mu_i - \mu_j)$에 대한 $100(1-\alpha)$% 신뢰구간은 다음과 같이 주어진다.

 

$(\bar{Y}_{i \cdot} - \bar{Y}_{j \cdot}) \pm t(N-t; 1-\alpha/2)\frac{\sqrt{2MSE}}{\sqrt{r}}\; (i\neq j)$

 

 

참고문헌

성내경 (2012). 실험설계와 분석 2판. 자유아카데미