선형모형에서 설계행렬(design matrix) \(\pmb X\)가 완전계수(full rank)행렬일 때 회귀계수의 추정치는 최소제곱법에서 구해진 정규방정식의 유일한 해로 구해진다.
\[ (\pmb X^t \pmb X ) \pmb \beta = \pmb X^t \pmb y \quad \Rightarrow \quad \hat {\pmb \beta} = (\pmb X^t \pmb X )^{-1} \pmb X^t \pmb y \]
그러나 여러 가지 실험이나 자료의 형태에서 설계행렬 \(\pmb X\)의 계수가 완전하지 않을 때가 있으며(less than full rank)
\[ rank(\pmb X) =r < p= \text{ number of columns in } \pmb X \]
이러한 경우에는 정규방정식에서 유일한 해가 존재하지 않는다. 이 장에서는 이러한 경우의 해결 방법을 알아 보고 일원 배치법에 어떻게 적용되는 자를 알아본다.
이제 일원배치법에 대한 통계적 모형에서 모수에 대한 추정을 생각해 보자.
\[ y_{ij} = \mu + \alpha_i + e_{ij} ,\quad i=1,2,\dots,a,~~j=1,2, \dots,r \tag{7.1}\]
추정해야할 모수는 전체 평균 \(\mu\)와 각 그룹의 처리 효과 \(\alpha_1,\alpha_2, \dots, \alpha_a\) 그리고 분산 \(\sigma_E^2\)이다. 전체 평균과 그룹의 효과는 오차제곱합(Sum of Square Error; SSE)을 최소로 하는 모수를 추정하는 최소제곱법(Least Square method; LS)으로 구할 수 있다.
\[ \min_{\mu, \alpha_1, \dots \alpha_a} \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \mu -\alpha_i)^2 =\min_{\mu, \alpha_1, \dots \alpha_a} SSE \tag{7.2}\]
위의 오차제곱합이 모든 모수에 대하여 미분가능한 이차식으므로 최소제곱 추정량은 제곱합을 모수에 대하여 미분하고 0 으로 놓아 방정식을 풀어서 얻을 수 있다.
오차제곱합을 모수 \(\mu\)와 \(\alpha_1,\alpha_2,\dots,\alpha_a\) 로 미분하여 0 으로 놓은 방정식은 다음과 같다.
\[ \begin{aligned} & \pardifftwo{}{\mu} SSE = -2 \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \mu -\alpha_i) = 0 \\ & \pardifftwo{}{\alpha_i} SSE = -2 \sum_{j=1}^r (y_{ij} - \mu -\alpha_i) = 0 , \quad i=1,2,\dots, a \end{aligned} \]
위의 방정식을 정리하면 다음과 같은 \(a+1\)개의 방정식을 얻는다.
\[ \begin{aligned} \mu +\frac{ \sum_{i=1}^a \alpha_i}{a} & = \bar {\bar y}\\ \mu + \alpha_1 & = \bar {y}_{1.} \\ \mu + \alpha_2 & = \bar {y}_{2.} \\ \cdots & \cdots \\ \mu + \alpha_a & = \bar {y}_{a.} \\ \end{aligned} \tag{7.3}\]
위의 방정식에서 첫 번째 방정식은 다른 \(a\)개의 방정식을 모두 합한 방정식과 같다. 따라서 모수는 \(a+1\)개이지만 실제 방정식의 개수는 \(a\)개이므로 유일한 해가 얻어지지 않는다. 따라서 유일한 해를 구하려면 하나의 제약조건이 필요하며 일반적으로 다음과 같은 두 개의 조건 중 하나를 사용한다.
첫 번째 효과 \(\alpha_1\)를 0으로 놓는 조건을 주는 것이다 (\(\alpha_1=0\)). set-to-zero 조건 하에서는 다음과 같은 추정량이 얻어진다.
\[ \hat \mu = \bar {y}_{1.}, \quad \hat \alpha_1=0, ~~ \hat \alpha_i = \bar {y}_{i.} -\bar {y}_{1.},~~i=2,\dots,a \tag{7.4}\]
처리들의 효과의 합은 0이라는 조건을 주는 것이다 ( \(\sum_{i=1}^a \alpha_i=0\)). sum-to-zero 조건에서는 계수의 추정치가 다음과 같이 주어진다.
\[ \hat \mu = \bar {\bar {y}}, \quad \hat \alpha_i = \bar {y}_{i.} -\bar {\bar {y}},~~i=1,2,\dots,a \tag{7.5}\]
여기서 유의할 점은 개별 모수들의 추정량은 조건에 따라서 달라지지만 집단의 평균을 나타내는 모수 \(\mu+ \alpha_i\) 에 대한 추정량은 언제나 같다.
\[ \widehat{\mu+ \alpha_i} = \hat \mu + \hat {\alpha}_i = \bar {y}_{i.} \]
만약에 자료를 아래와 같은 평균 모형으로 나타낼 경우에는 각 평균 \(\mu_i\) 는 각 그룹의 표본 평균으로 추정된다.
\[ y_{ij} = \mu_i + e_{ij} \]
평균 모형에서 각 그룹의 모평균에 대한 최소제곱 추정량은 \(\hat \mu_i = \bar {y}_{i.}\) 이며 이는 주효과 모형에서의 추정량과 동일하다.
또한 모형에 관계없이 오차항의 분산 \(\sigma_E^2\) 에 대한 추정량은 다음과 같이 주어진다.
\[ \hat \sigma_E^2 = \frac{ \sum_i \sum_j (y_{ij} - \hat \mu -\hat \alpha_i )^2}{a(r-1)} \]
일원배치 모형 식 7.1 를 다음과 같은 벡터를 이용한 선형모형(linear model, regression model) 형태로 나타내고자 한다.
\[ \pmb y = \pmb X \pmb \beta +\pmb e \tag{7.6}\]
위의 선형모형식의 요소 \(\pmb y\), \(\pmb X\), \(\pmb \beta\), \(\pmb e\)는 다음과 같은 벡터와 행렬로 표현된다.
\[ \begin{bmatrix} y_{11} \\ y_{12} \\ \vdots \\ y_{1r} \\ y_{21} \\ y_{22} \\ \vdots \\ y_{2r} \\ \vdots \\ y_{a1} \\ y_{a2} \\ \vdots \\ y_{ar} \\ \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & . & . & 0 \\ 1 & 1 & 0 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & 0 & . & . & 0 \\ 1 & 0 & 1 & . & . & 0 \\ 1 & 0 & 1 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & . & . & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & . & . & 1 \\ 1 & 0 & 0 & . & . & 1 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & . & . & 1 \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{1} \\ \alpha_{2} \\ \vdots \\ \alpha_{a} \\ \end{bmatrix} + \begin{bmatrix} e_{11} \\ e_{12} \\ \vdots \\ e_{1r} \\ e_{21} \\ e_{22} \\ \vdots \\ e_{2r} \\ \vdots \\ e_{a1} \\ e_{a2} \\ \vdots \\ e_{ar} \\ \end{bmatrix} \tag{7.7}\]
이제 위에서 논의한 최소제곱법을 선형 모형 식 7.6 에 적용하면 다음과 같이 표현할 수 있다.
\[ \min_{\mu, \alpha_1, \dots \alpha_a} \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \mu -\alpha_i)^2 = \min_{\pmb \beta } ( \pmb y - \pmb X \pmb \beta )^t( \pmb y - \pmb X \pmb \beta ) \tag{7.8}\]
최소제곱법의 기준을 만족하는 계수 \(\pmb \beta\)는 다음과 같은 정규방정식(normal equation)의 해(solution)이다.
\[ \pmb X^t \pmb X \pmb \beta = \pmb X^t \pmb y \tag{7.9}\]
정규방정식 식 7.9 을 일워배치의 선형모형식 식 10.16 에 나타난 \(\pmb y\), \(\pmb X\)로 이용하여 나타내면 다음과 같다.
\[ \begin{bmatrix} ar & r & r & \cdot & \cdot & r \\ r & r & 0 & \cdot & \cdot & 0 \\ r & 0 & r & \cdot & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ r & 0 & 0 & \cdot & \cdot & r \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{1} \\ \alpha_{2} \\ \cdot \\ \cdot \\ \alpha_{a} \\ \end{bmatrix} = \begin{bmatrix} ar \bar {\bar y} \\ r {\bar y}_{1.}\\ r \bar y_{2.}\\ \cdot \\ \cdot \\ r \bar y_{a.} \end{bmatrix} \tag{7.10}\]
정규방정식 식 7.10 는 위에서 구한 최소제곱법에서 유도된 방정식 식 7.3 과 같다.
여기서 유의할 점은 선형모형식 식 10.16 의 계획행렬 \(\pmb X\) 가 완전 계수(full rank) 행렬이 아니다. 계획행렬 \(\pmb X\)의 첫 번째 열은 다른 열을 합한 것과 같다. 또한 정규 방정식 식 7.10 에서 \(\pmb X^t \pmb X\) 행렬도 완전계수 행렬이 아니다. 따라서 \(\pmb X^t \pmb X\) 행렬의 역행렬은 존재하지 않는다.
이러한 이유로 모수에 대한 유일한 추정량이 존재하지 않기 때문에 앞에서 언급한 제약 조건을 고려해야 정규방정식을 풀 수 있다.
만약 Set-to-zero 조건을 가정한다면 모수에서 \(\alpha_1\)을 제외하고 선형모형식 식 10.16 를 다음과 같이 다시 표현할 수 있다.
효과 \(\alpha_1\)을 0 으로 놓는다는 것은 \(\alpha_1\)을 추정할 필요가 없으므로 모수벡터 \(\pmb \beta\) 에서 \(\alpha_1\)를 빼고 게획행렬에서도 대응하는 열을 제거하는 것이다.
\[ \begin{bmatrix} y_{11} \\ y_{12} \\ \vdots \\ y_{1r} \\ y_{21} \\ y_{22} \\ \vdots \\ y_{2r} \\ \vdots \\ y_{a1} \\ y_{a2} \\ \vdots \\ y_{ar} \\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & . & . & 0 \\ 1 & 0 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & . & . & 0 \\ 1 & 1 & . & . & 0 \\ 1 & 1 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & . & . & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & . & . & 1 \\ 1 & 0 & . & . & 1 \\ 1 & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & . & . & 1 \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{2} \\ \alpha_{3} \\ \vdots \\ \alpha_{a} \\ \end{bmatrix} + \begin{bmatrix} e_{11} \\ e_{12} \\ \vdots \\ e_{1r} \\ e_{21} \\ e_{22} \\ \vdots \\ e_{2r} \\ \vdots \\ e_{a1} \\ e_{a2} \\ \vdots \\ e_{ar} \\ \end{bmatrix} \tag{7.11}\]
이제 수정된 모형식 식 7.11 에 최소제곱법을 적용하여 정규방정식을 구하면 다음과 같은 방정식을 얻는다.
\[ \begin{bmatrix} ar & r & r & \cdot & \cdot & r \\ r & r & 0 & \cdot & \cdot & 0 \\ r & 0 & r & \cdot & \cdot & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ r & 0 & 0 & \cdot & \cdot & r \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{2} \\ \alpha_{3} \\ \cdot \\ \cdot \\ \alpha_{a} \\ \end{bmatrix} = \begin{bmatrix} ar \bar {\bar y} \\ r {\bar y}_{2.}\\ r \bar y_{3.}\\ \cdot \\ \cdot \\ r \bar y_{a.} \end{bmatrix} \tag{7.12}\]
위의 정규방정 식 7.12 를 풀면 위에서 언급한 sum-to-zero 조건에서 구해지는 모수의 추정량 식 7.4 를 얻을 수 있다.
이제 Sum-to-zero 조건에서 모수의 추정에 대해 알아보자. 조건 \(\sum_{i=1}^a \alpha_i =0\) 조건을 마지막 모수 \(\alpha_a\)에 대하여 표현하면 다음과 같다.
\[ \alpha_a = -\alpha_1 - \alpha_2 - \dots - \alpha_{a-1} \]
따라서 마지막 처리 \(\alpha_a\) 에 대한 관측값에 대한 모형은 다음과 같아 쓸 수 있다.
\[ y_{aj} = \mu + \alpha_a + e_{aj} = \mu +( -\alpha_1 - \alpha_2 - \dots - \alpha_{a-1}) + e_{ij} \]
이러한 결과를 모형방정식에 반영한다. 즉, 모수벡터 \(\pmb \beta\) 에서 \(\alpha_a\)를 제거하고 게획행렬에 위의 마지막 처리에 대한 효과식을 반영하면 다음과 같은 선형모형식을 얻는다.
\[ \begin{bmatrix} y_{11} \\ y_{12} \\ \vdots \\ y_{1r} \\ y_{21} \\ y_{22} \\ \vdots \\ y_{2r} \\ \vdots \\ y_{a1} \\ y_{a2} \\ \vdots \\ y_{ar} \\ \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & . & . & 0 \\ 1 & 1 & 0 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & 0 & . & . & 0 \\ 1 & 0 & 1 & . & . & 0 \\ 1 & 0 & 1 & . & . & 0 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & . & . & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & . & . & 1 \\ 1 & 0 & 0 & . & . & 1 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & . & . & 1 \\ 1 & 0 & 0 & . & . & 1 \\ 1 & -1 & -1 & . & . & -1 \\ 1 & -1 & -1 & . & . & -1 \\ 1 & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & -1 & -1 & . & . & -1 \\ 1 & -1 & -1 & . & . & -1 \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{1} \\ \alpha_{2} \\ \vdots \\ \alpha_{a-1} \\ \end{bmatrix} + \begin{bmatrix} e_{11} \\ e_{12} \\ \vdots \\ e_{1r} \\ e_{21} \\ e_{22} \\ \vdots \\ e_{2r} \\ \vdots \\ e_{a1} \\ e_{a2} \\ \vdots \\ e_{ar} \\ \end{bmatrix} \tag{7.13}\]
이제 수정된 모형식 식 7.13 에 최소제곱법을 적용하여 정규방정식을 구하면 다음과 같은 방정식을 얻는다.
\[ \begin{bmatrix} ar & 0 & 0 & \cdot & \cdot & 0 \\ 0 & 2r & r & \cdot & \cdot & r \\ 0 & r & 2r & \cdot & \cdot & r \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & r & r & \cdot & \cdot & 2r \\ \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_{1} \\ \alpha_{2} \\ \cdot \\ \cdot \\ \alpha_{a-1} \\ \end{bmatrix} = \begin{bmatrix} ar \bar {\bar y} \\ r {\bar y}_{1.}-r {\bar y}_{a.} \\ r \bar y_{2.}-r {\bar y}_{a.}\\ \cdot \\ \cdot \\ r \bar y_{a-1,.} -r {\bar y}_{a.} \end{bmatrix} \tag{7.14}\]
위의 정규방정 식 7.14 를 풀면 위에서 언급한 sum-to-zero 조건에서 구해지는 모수의 추정량 식 7.5 를 얻을 수 있다.
교과서 예제 7.3 은 4년제 대학교의 학년별 영어시험 점수 자료이다. 이 자료는 4개의 학년에 대하여 각각 6명의 학생들에 대한 영어시험 점수이다. 이 자료를 이용하여 일원배치법을 적용하고 각 학년의 평균 점수에 대한 추정량을 구해보자.
head(english1) score grade
1 81 1
2 75 1
3 69 1
4 90 1
5 72 1
6 83 1english1$grade <- factor(english1$grade)
summary(english1) score grade
Min. :62.00 1:6
1st Qu.:72.00 2:6
Median :79.00 3:5
Mean :77.33 4:4
3rd Qu.:81.00
Max. :94.00 이제 다음과 같이 일원배치 모형을 적합하고 결과를 보자.
fit_anova_set0 <- lm(score ~ grade, data=english1)
summary(fit_anova_set0)
Call:
lm(formula = score ~ grade, data = english1)
Residuals:
Min 1Q Median 3Q Max
-9.500 -5.500 0.600 4.667 11.667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.333 2.868 27.312 1.75e-15 ***
grade2 -3.833 4.056 -0.945 0.3579
grade3 -6.933 4.254 -1.630 0.1215
grade4 9.167 4.535 2.021 0.0593 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.025 on 17 degrees of freedom
Multiple R-squared: 0.4341, Adjusted R-squared: 0.3342
F-statistic: 4.347 on 3 and 17 DF, p-value: 0.01905함수 lm() 을 이용하여 일원배치 모형을 적합한 결과를 보면 default 로 grade 에 대하여 1학년 효과 grade1 는 0 으로 고정되고 grade2, grade3, grade4 에 대한 추정치는 각각 -3.83, - 6.933, 9.167 로 추정된다. 즉, 범주형 변수의 첫 번째 수준을 0으로 고정하고 이를 기준으로 다른 수준에 대한 효과를 추정한 것이다.
\[ \hat \mu = 78.333, \quad \hat \alpha_1 =0, \quad \hat \alpha_2 = -3.83, \quad \hat \alpha_3 = -6.933, \quad \hat \alpha_4 = 9.167 \] 사용된 계획행렬을 보면 다음과 같다.
model.matrix(fit_anova_set0) (Intercept) grade2 grade3 grade4
1 1 0 0 0
2 1 0 0 0
3 1 0 0 0
4 1 0 0 0
5 1 0 0 0
6 1 0 0 0
7 1 1 0 0
8 1 1 0 0
9 1 1 0 0
10 1 1 0 0
11 1 1 0 0
12 1 1 0 0
13 1 0 1 0
14 1 0 1 0
15 1 0 1 0
16 1 0 1 0
17 1 0 1 0
18 1 0 0 1
19 1 0 0 1
20 1 0 0 1
21 1 0 0 1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$grade
[1] "contr.treatment"이제 각 학년에 영어 평균 \(\mu + \alpha_i\) 에 대한 추정값 \(\widehat{ \mu + \alpha_i}\) 를 구해보자.
\[ \widehat{\mu + \alpha_i} = \hat {\mu} + \hat {\alpha}_i \]
emmeans::emmeans(fit_anova_set0, "grade") grade emmean SE df lower.CL upper.CL
1 78.3 2.87 17 72.3 84.4
2 74.5 2.87 17 68.4 80.6
3 71.4 3.14 17 64.8 78.0
4 87.5 3.51 17 80.1 94.9
Confidence level used: 0.95 함수 lm 에서 default 로 설정되는 set-to-zero 조건은 다음과 명령어로 지정할 수 있다.
options(contrasts=c("contr.treatment", "contr.poly"))이제 일원배치 모형에서 sum-to-zero 조건을 적용하여 모수를 추정해 보자. sum-to-zero 조건을 적용하려면 다음과 같은 명령어를 실행해야 한다.
options(contrasts=c("contr.sum", "contr.poly"))이제 sum-to-zero 조건을 적용하여 일원배치 모형을 적합하고 결과를 보자.
fit_anova_sum0 <- lm(score ~ grade, data=english1)
summary(fit_anova_sum0)
Call:
lm(formula = score ~ grade, data = english1)
Residuals:
Min 1Q Median 3Q Max
-9.500 -5.500 0.600 4.667 11.667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.933 1.554 50.135 <2e-16 ***
grade1 0.400 2.555 0.157 0.8775
grade2 -3.433 2.555 -1.344 0.1967
grade3 -6.533 2.711 -2.410 0.0276 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.025 on 17 degrees of freedom
Multiple R-squared: 0.4341, Adjusted R-squared: 0.3342
F-statistic: 4.347 on 3 and 17 DF, p-value: 0.01905sum-to-zero 조건하에서 계획행렬은 다음과 같다.
model.matrix(fit_anova_sum0) (Intercept) grade1 grade2 grade3
1 1 1 0 0
2 1 1 0 0
3 1 1 0 0
4 1 1 0 0
5 1 1 0 0
6 1 1 0 0
7 1 0 1 0
8 1 0 1 0
9 1 0 1 0
10 1 0 1 0
11 1 0 1 0
12 1 0 1 0
13 1 0 0 1
14 1 0 0 1
15 1 0 0 1
16 1 0 0 1
17 1 0 0 1
18 1 -1 -1 -1
19 1 -1 -1 -1
20 1 -1 -1 -1
21 1 -1 -1 -1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$grade
[1] "contr.sum"sum-to-zero 조건하에서 일원배치 모형을 적합한 결과를 보면 다음과 같다.
\[ \hat \mu = 77.933, \quad \hat \alpha_1 = 0.400, \quad \hat \alpha_2 = -3.433, \quad \hat \alpha_3 = - 6.533, \quad \hat \alpha_4 = -(\hat \alpha_1 + \hat \alpha_2 + \hat \alpha_3) = 9.566 \]
이렇게 set_to_zero 조건과 sum-to-zero 조건을 적용한 결과는 다르지만 각 학년의 평균에 대한 추정치는 동일하다.
emmeans::emmeans(fit_anova_sum0, "grade") grade emmean SE df lower.CL upper.CL
1 78.3 2.87 17 72.3 84.4
2 74.5 2.87 17 68.4 80.6
3 71.4 3.14 17 64.8 78.0
4 87.5 3.51 17 80.1 94.9
Confidence level used: 0.95 참고로 교과서에서 4학년을 기준으로 다른 학년들과 비교하기 위하여 변수 grade 의 수준(level) 의 순서를 설정할 때 다음과 같이 4학년을 첫 번째 수준으로 설정하면 된다.
options(contrasts=c("contr.treatment", "contr.poly"))
english1$grade <- factor(english1$grade, levels = c(4,1,2,3))
str(english1)'data.frame': 21 obs. of 2 variables:
$ score: int 81 75 69 90 72 83 65 80 73 79 ...
$ grade: Factor w/ 4 levels "4","1","2","3": 2 2 2 2 2 2 3 3 3 3 ...fit_anova_set0_1 <- lm(score ~ grade, data=english1)
summary(fit_anova_set0_1)
Call:
lm(formula = score ~ grade, data = english1)
Residuals:
Min 1Q Median 3Q Max
-9.500 -5.500 0.600 4.667 11.667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 87.500 3.513 24.910 8.06e-15 ***
grade1 -9.167 4.535 -2.021 0.05927 .
grade2 -13.000 4.535 -2.867 0.01069 *
grade3 -16.100 4.713 -3.416 0.00329 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.025 on 17 degrees of freedom
Multiple R-squared: 0.4341, Adjusted R-squared: 0.3342
F-statistic: 4.347 on 3 and 17 DF, p-value: 0.01905설계행렬 \(\pmb X\)의 계수가 완전하지 않을 때 회귀 계수를 추정하기 위한 방법으로서 다음과 같은 세 가지 방법이 있다.
\(\pmb X\)의 계수가 완전하지 않을 때 설계행렬의 열을 다시 구성하여 계수를 완전하게 하는 방법이 있다. 즉 \(\pmb X = (\pmb X_1, \pmb X_2)\)으로 표시하고 \(\pmb X_1\)을 \(n \times r~ (r < p)\)라고 하며 어떤 행렬 \(\pmb F\)가 존재하여 \(\pmb X_2 = \pmb X_1 \pmb F\)의 관계를 가진다고 가정하자. 이러한 관계는 \(\pmb X_2\)의 열들이 \(\pmb X_1\)의 열들의 선형결합으로 표현될 수 있다는 것을 의미한다. 이러한 경우에 선형모형은 다음과 같이 표현될 수 있다.
\[ \pmb y= \pmb X \pmb \beta + \pmb e = \pmb X_1 (\pmb I, \pmb F)\pmb \beta + \pmb e = \pmb X_1 \pmb \alpha +\pmb e \] 여기서 새롭게 조정된 계수 \(\pmb \alpha\)와 처음의 계수 \(\pmb \beta\)는 다음과 같은 관계가 있다.
\[ \pmb \alpha = (\pmb I, \pmb F)\pmb \beta = (\pmb \beta_1, \pmb \beta_2) \]
따라서 새롭게 구성된 선형모형 \(\pmb y=\pmb X_1 \pmb \alpha +\pmb e\)에서 새로운 계수의 추정치는 \(\hat {\pmb \alpha} = (\pmb X_1^t \pmb X_1 )^{-1} \pmb X_1^t \pmb y\) 이다.
회귀게수에 부가 조건(side condition)을 주면 유일한 계수의 추정치를 구할 수 있다. 즉 \((p-r) \times p\) 행렬 \(\pmb H\)를 고려하고 \(\pmb H \pmb \beta =0\)이라는 부가조건을 가정하자. 즉 모든 \(\pmb \eta = R(\pmb X)\)에 대하여 \(\pmb \eta=\pmb X \pmb \beta\)와 \(\pmb H \pmb \beta =0\)를 만족하는 \(\pmb \beta\)는 유일하게 존재한다.
이러한 부가 조건 \(\pmb H \pmb \beta =0\)과 정규방정식 \((\pmb X^t \pmb X ) \pmb \beta = \pmb X^t \pmb y\)를 동시에 만족하는 유일한 해를 구하고 이를 최소제곱추정량으로 한다. 이러한 부가 조건을 주는 방법은 분산분석을 이용하는 여러 가지 선형 모형 (예: 일원 배치법)에 자주 사용된다.
\(\pmb X\)의 계수가 완전하지 않을 때 일반화 역행렬(generalized inverse matrix)를 이용하면 회귀계수의 추정치를 구할 수 있다.
여기서 \(m \times n\) 행렬 \(\pmb A\)의 일반화 역행렬 \(\pmb A^{-}\)는 다음을 만족하는 행렬이다.
\[ \pmb A = \pmb A \pmb A^{-} \pmb A \]
일반화 역행렬은 일반적으로 유일하지 않다. \(\pmb A\)가 정방행렬이고 정칙행렬일 때 유일하게 존재하며 \(\pmb A^- = \pmb A^{-1}\)이다. 정규방정식 \(\pmb X^t \pmb y=\pmb X^t \pmb X \hat {\pmb \beta}\)의 양변에 \(\pmb X^t \pmb X (\pmb X^t \pmb X )^-\)를 곱하면
\[\pmb X^t \pmb X (\pmb X^t \pmb X )^- \pmb X \pmb y = \pmb X^t \pmb X (\pmb X^t \pmb X )^- \pmb X^t \pmb X \hat {\pmb \beta} = \pmb X^t \pmb X \hat {\pmb \beta} = \pmb X^t \pmb y \]
이므로 \(\hat {\pmb \beta} = (\pmb X^t \pmb X )^- \pmb X^t \pmb y\)는 정규방정식의 해가 된다. 앞에서 언급하였듯이 일반화 역함수를 이용한 계수의 추정량은 유일하지 않다. 그러나 반응변수의 추정량 \(\hat {\pmb y} = \pmb X \hat {\pmb \beta}\)는 추정된 계수에 관계없이 유일하다.
앞 절에서 보았듯이 일원배치법을 선형 모형식으로 표현하는 경우 평균에 대한 모수는 모두 \(a+1\) 개가 있다.
\[ \mu, \alpha_1, \alpha_2, \cdots, \alpha_a \]
하지만 모형식에서 계획행렬 \(\pmb X\)가 완전 계수 행렬이 아니기 때문에 1개의 제약 조건을 가정하고 모수를 추정하였다. 하지만 제약 조건이 달라지면 각 모수의 추정량이 달라지기 때문에 각 모수는 유일한 값으로 추정이 불가능하다.
이렇게 각 모수들은 제약 조건에 따라서 유일하게 추정이 불가능하지만 앞 절에서 보았듯이 \(\mu + \alpha_i\) 에 대한 추정량은 제약조건에 관계없이 표본 평균 \(\bar y_{i.}\)으로 동일하게 추정되어 진다.
그러면 어떤 모수들은 유일하게 추정이 불가능하고 어떤 모수들이 유일하게 추정이 가능할까?
이제 제약조건이 달라도 유일하게 추정이 가능한 모수들의 형태를 살펴보자.
선형모형 \(\pmb y =\pmb X \pmb \beta + \pmb e\) 에서 계획행렬 \(\pmb X\)의 계수가 완전하지 않으면 모수 벡터 \(\pmb \beta\)는 유일한 값으로 추정할 수 없다.
이제 임의의 벡터 \(\pmb c\)가 있을 때 모수들의 선형결합 \(\psi = \pmb c^t \pmb \beta\)를 고려하자.
예를 들어 일원배치 모형에서는 다음과 같은 모수들의 선형결합을 고려하는 것이다.
\[ \psi = \pmb c^t \pmb \beta = [ c_0~ c_1~ c_2~ \cdots~~c_a] \begin{bmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_a \end{bmatrix} =c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 + \cdots + c_a \alpha_a \]
위에서 본 것처럼 하나의 모수 \(\alpha_1\)에 대한 유일한 추정은 불가능하다.
\[ \alpha_1 = (0) \mu + (1) \alpha_1 + (0) \alpha_2 + \cdots + (0) \alpha_a \]
하지만 모수의 조합 \(\mu+ \alpha_2\) 은 유일한 추정이 가능하다.
\[ \mu + \alpha_1 = (1) \mu + (1) \alpha_1 + (0) \alpha_2 + \cdots + (0) \alpha_a \]
이제 문제는 선형조합 \(\psi= \pmb c^t \pmb \beta\) 에서 계수들 \(c_0, c_1, \dots, c_a\)가 어떤 값을 가지는 경우 유일한 추정이 가능한 지 알아내는 것이다.
이제 \(\psi = \pmb c^t \pmb \beta\) 에 대한 유일한 추정량 \(\hat \psi\) 이 있다고 가정하자. 선형 모형에서 추정량 \(\hat \psi\)의 형태는 관측값에 대한 선형함수가 되어야 한다. 따라서 추정량을 \(\hat \psi = \pmb a^t \pmb y\) 로 나타낼 수 있다. 이제 추정량 \(\hat \psi\)의 기대값은 \(\psi=\pmb c^t \pmb \beta\)이어야 하므로 다음이 성립해야 한다.
\[ E(\hat \psi| \pmb X) = E(\pmb a^t \pmb y| \pmb X) = \pmb a^t E(\pmb y| \pmb X) = \pmb a^t \pmb X \pmb \beta = \pmb c^t \pmb \beta \]
위의 식에서 가장 마지막 두 항의 관계를 보면 다음이 성립해야 한다.
\[ \pmb a^t \pmb X = \pmb c^t \quad \text{ equivalently }\quad \pmb c = \pmb X^t \pmb a \tag{7.15}\]
즉 추정가능한 모수의 조합 \(\psi = \pmb c^t \pmb \beta\)에서 계수 벡터 \(\pmb c\) 는 계획행렬에 있는 행들의 선형 조합으로 표시되어야 한다는 것이다. 이렇게 유일하게 추정이 가능한 모수의 조합을 추정가능한 함수(estimable function)이라고 한다.
2개의 수준이 있고 반복이 2번 있는 일원배치 \((a=2,r=2)\) 에 대한 선형모형 식 10.16 을 생각해보자. 이 경우 계획행렬 \(\pmb X\) 과 모수벡터 \(\pmb \beta\) 는 다음과 같다.
\[ \pmb X = \begin{bmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{bmatrix} \quad \pmb \beta = \begin{bmatrix} \mu \\ \alpha_1 \\ \alpha_2 \end{bmatrix} \]
이제 유일하게 추정 가능한 모수 조합 \(\psi\) 은 어떤 형태일까?
\[ \psi = \pmb c^t \pmb \beta = c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 \]
위의 식 7.15 에서 추정가능한 모수의 조합에 대한 계수 벡터 \(\pmb c\) 는 다음과 같은 조건을 만족해야 한다.
\[ \pmb c = {\pmb X}^t \pmb a \]
이제 임의의 벡터 \(\pmb a\) 에 대하여 \(\pmb c= \pmb X^t \pmb a\)의 형태를 보자.
\[ \begin{aligned} \pmb c &= \pmb X^t \pmb a \\ & = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \end{bmatrix} \\ & = a_1 \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + a_2 \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + a_3 \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} + a_4 \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \\ & = (a_1 + a_2) \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + (a_3 + a_4) \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \\ &= b_1 \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + b_2 \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \end{aligned} \tag{7.16}\]
이제 \(\pmb X^t \pmb a\) 는 계획행렬 \(\pmb X\)에 있는 유일한 행들의 선형조합임을 알 수 있다.
따라서 유일하게 추정 가능한 모수의 선형조합 \(\psi = \pmb c^t \pmb \beta\) 에 대한 계수 벡터 \(\pmb c =[ c_0 ~ c_1 ~ c_2]^t\) 는 계획행렬 \(\pmb X\)의 유일한 행들의 선형 조합으로 구성되어야 한다.
\[ \pmb c = \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = b_1 \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + b_2 \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \tag{7.17}\]
첫 번째 처리에 대한 효과 모수 \(\alpha_1\) 를 선형조합으로 나타내면
\[ \alpha_1 = c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 = (0) \mu + (1) \alpha_1 + (0) \alpha_2 \]
따라서 조건 식 7.17 에서 \(\pmb c^t = [0~1~0]\)을 만들수 있는 계수 \(b_1\)과 \(b_2\)를 찾아야 하는데 이는 불가능하다. 따라서 모수 \(\alpha_1\) 은 추정 불가능하다.
\[ \pmb c = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} = b_1 \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + b_2 \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \]
모수 조합 \(\mu + \alpha_1\) 를 선형조합으로 나타내면
\[ \mu + \alpha_1 = c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 = (1) \mu + (1) \alpha_1 + (0) \alpha_2 \]
따라서 조건 식 7.17 에서 \(\pmb c^t = [1~1~0]\)을 만들수 있는 계수는 \(b_1=1\)과 \(b_2=0\) 이므로 추정이 가능하다.
\[ \pmb c = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} = (1) \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + (0) \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \]
\[ \alpha_1 -\alpha_2= c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 = (0) \mu + (1) \alpha_1 + (-1) \alpha_2 \]
따라서 조건 식 7.17 에서 \(\pmb c^t = [0~1~-1]\)을 만들수 있는 계수는 \(b_1=1\)과 \(b_2=-1\) 이므로 추정이 가능하다.
\[ \pmb c = \begin{bmatrix} 0 \\ 1 \\ -1 \end{bmatrix} = (1) \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix} + (-1) \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix} \]
이제 효과모형 식 7.1 을 고려하면 집단들 사이에 차이가 있는지에 대한 가설을 다음과 같이 고려할수 있다. 집단에 대한 효과가 모두 0이 되면 집단 간의 평균에 대한 차이는 없다.
\[ H_0: \alpha_1 = \alpha_2 =\cdots=\alpha_a =0 \quad \text{vs.} \quad H_1: \text{ not } H_0 \tag{7.18}\]
이제 집단 간의 변동(각 집단의 평균의 차이가 얼마나 나는지에 대한 통계량)과 집단 내의 변동(각 집단내에서 관측값들의 퍼진 정도)를 측정하는 통계량을 찾아서 검정 통계량을 구성해 보자.
일단 각 집단의 반복 측정값의 횟수는 모두 같다고 가정하자(\(r_i=r\)). 전체 평균과 집단의 평균을 정의하자.
\[ \bar {y}_{..} = \frac{\sum_{i=1}^a \sum_{j=1}^r y_{ij}}{ar}, \quad \bar {y}_{i.} = \frac{\sum_{j=1}^r y_{ij}}{r} \]
이제 하나의 관측값 \(y_{ij}\)과 전체 평균 \(\bar {y}_{..}\) 간의 편차(deviation)를 다음과 같이 분해해 보자.
\[ \underbrace{ y_{ij} - \bar {y}_{..} }_{\text{total deviation}} = \underbrace{ ( y_{ij} - \bar {y}_{i.} )}_{\text{within-group deviation}} + \underbrace{(\bar {y}_{i.} - \bar {y}_{..} )}_{\text{between-group deviation}} \tag{7.19}\]
식 7.19 에서 집단 평균과 총 평균의 편차 (\(\bar {y}_{i.} - \bar {y}_{..}\))는 처리의 효과를 측정할 수 있는 통계량이다. 집단 간의 차이를 반영하는 양으로 처리 효과 \(\alpha_i\)들에 의하여 발생한다.
집단 내의 관측값과 집단 평균의 차이 (\(y_{ij} - \bar {y}_{i.}\))는 집단 내의 변동을 나타내는 통계량으로 측정 오차 \(e_{ij}\)에 의하여 발생한다.
식 7.19 의 각 편차들은 양수와 음수로서 부호를 가지기 때문에 이를 변동으로 표현하기 위하여 차이를 제곱하여 합친 제곱합(sum of squares)을 고려해 보자.
\[ \begin{aligned} \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \bar {y}_{..} )^2 & = \sum_{i=1}^a \sum_{j=1}^r \left [ ( y_{ij} - \bar {y}_{i.} ) + (\bar {y}_{i.} - \bar {y}_{..} ) \right ]^2 \\ & = \sum_{i=1}^a \sum_{j=1}^r ( y_{ij} - \bar {y}_{i.} )^2 + \sum_{i=1}^a \sum_{j=1}^r (\bar {y}_{i.} - \bar {y}_{..} )^2 + 2 \sum_{i=1}^a \sum_{j=1}^r ( y_{ij} - \bar {y}_{i.} ) (\bar {y}_{i.} - \bar {y}_{..} ) \\ & = \sum_{i=1}^a \sum_{j=1}^r ( y_{ij} - \bar {y}_{i.} )^2 + \sum_{i=1}^a r (\bar {y}_{i.} - \bar {y}_{..} )^2 + 0 (why?) \\ \end{aligned} \]
결과적으로 다음과 같은 변동의 분해를 제곱합의 형식으로 얻을 수 있다.
\[ \underbrace{ \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \bar {y}_{..} )^2 }_{\text{total variation}} = \underbrace{ \sum_{i=1}^a \sum_{j=1}^r ( y_{ij} - \bar {y}_{i.} )^2 }_{\text{within-group variation}} + \underbrace{\sum_{i=1}^a \sum_{j=1}^r (\bar {y}_{i.} - \bar {y}_{..} )^2 }_{\text{between-group variation}} \tag{7.20}\]
분해식 식 7.20 에서 나타난 각 제곱합에 대한 이름과 의미를 살펴보자.
\[ SST = \sum_{i=1}^a \sum_{j=1}^r (y_{ij} - \bar {y}_{..} )^2 \]
\[ SSE = \sum_{i=1}^a \sum_{j=1}^r ( y_{ij} - \bar {y}_{i.} )^2 \]
\[ SSA = \sum_{i=1}^a \sum_{j=1}^r (\bar {y}_{i.} - \bar {y}_{..} )^2 =\sum_{i=1}^a r (\bar {y}_{i.} - \bar {y}_{..} )^2 \]
이제 분해식 식 7.20 을 다음과 같이 나타낼수 있다.
\[ SST = SSA + SSE \tag{7.21}\]
위의 분해식에서 볼 수 있듯이 집단 간의 변동의 크기를 나타내는 처리제곱합이 커질수록, 또는 집단내의 변동의 크기를 나타내는 오차제곱합이 작아질수록 귀무가설에 반대되는(즉, 집단 간의 평균이 유의한 차이가 난다는) 증거가 강해진다.
이제 가설 식 7.18 을 검정하기 위한 통계량을 구성해 보자. 먼저 다음과 같은 제곱합들을 각 자유도로 나눈 평균제곱합(Mean Sum of Squares)를 정의한다.
\[ MS_A = \frac{SSA}{\phi_A}, \quad MSE =\frac{SSE}{\phi_E} \tag{7.22}\]
집단 간의 변동과 집단 내의 변동의 상대적 비율로 그룹 간의 차이를 검정할 수 있다는 개념을 확장하여 다음과 같은 F-통계량 \(F_0\) 를 만들어 보자.
\[ F_0 = \frac{MSA}{MSE} = \frac{\text{between-group variation}} {\text{within-group variation}} \tag{7.23}\]
위 식 7.23 에서 정의된 F-통계량은 그룹 간에 평균의 차이가 클수록, 그룹 내의 차이가 작을 수록 그 값이 커진다. 따라서 F-통계량의 값이 크면 클수록 귀무가설에 반대되는 증거가 강해진다.
이렇게 전체의 변동을 집단 간의 변동과 집단 내의 변동으로 나누어 집단 간의 평균의 차이를 추론하는 방법을 분산분석(Analysis of Variance, ANOVA)이라고 한다.
이제 식 7.23 에서 정의된 F-통계량을 이용하여 가설 @ref(eq:hypo1)를 검정하는 통계적 방법을 만들어 보자. 일단 두 제곱합의 통계적 성질은 다음과 같다.
잔차 제곱합을 오차항의 분산으로 나눈 통계량은 자유도가 \(\phi_E\) 를 가지는 카이제곱 분포를 따른다.
\[ \frac{SSE}{\sigma_E^2} \sim \chi^2(\phi_E) \]
귀무가설이 참인 경우 처리 제곱합을 오차항의 분산으로 나눈 통계량은 자유도가 \(\phi_A\) 를 가지는 카이제곱 분포를 따른다.
\[ \frac{SSA}{\sigma_E^2} \sim \chi^2(\phi_A) \quad \text{ under } H_0 \]
잔차 제곱합과 처리 제곱합은 서로 독립이다.
따라서 귀무가설이 참인 경우 F-통계량은 자유도가 \(\phi_A, \phi_E\)를 가지는 F-분포를 따른다.
\[ F_0 = \frac{MSA}{MSE} = \frac{ \tfrac{SSA/\sigma_E^2}{\phi_A}} {\tfrac{SSE/\sigma_E^2}{\phi_E }} \sim F(\phi_A, \phi_E) \quad \text{ under } H_0 \] {#eq-anova-ftest)
유의수준 \(\alpha\)에서 F-통계량이 기각역을 벗어나면 귀무가설을 기각한다.
\[ \text{ Reject } H_0 \text{ if } F_0 > F(1-\alpha, \phi_A, \phi_E) \]
또는 다음과 같이 계산된 p-값이 유의수준 \(\alpha\) 보다 작으면 귀무가설을 기각한다.
\[ p-value = P[F(\phi_A, \phi_E) > F_0 ] \]
F-통계량을 정의할 때 편리하고 유용하게 사용되는 것이 다음과 같은 분산분석표(ANOVA table)이다.
| 요인 | 제곱합 | 자유도 | 평균제곱합 | \(F_0\) | p-값 |
|---|---|---|---|---|---|
| 처리 | \(SSA\) | \(\phi_A = a-1\) | \(MSA=SSA/\phi_A\) | \(F_0=MSA/MSE\) | \(P[F(\phi_A, \phi_E) > F_0 ]\) |
| 잔차 | \(SSE\) | \(\phi_E=a(r-1)\) | \(MSE=SSE/\phi_E\) | ||
| 총합 | \(SST\) | \(\phi_T = ar-1\) |
분산분석표를 이용한 F-검정으로 귀무가설을 기각하면 모든 처리 수준의 평균이 같지 않다는 결론을 내리고 어떤 집단 간에 평균의 차이가 유의한지 더 분석해야 한다. 평균 차이에 대한 신뢰구간과 가설 검정은 아래와 같이 주어진다.
두 수준 평균의 차이 \(\delta_{ij} = \mu_i - \mu_j\) 에 대한 \(100(1-\alpha)\) % 신뢰구간은 다음과 같이 주어진다.
\[ ( \bar {x}_{i.} - \bar {x}_{j.}) \pm t(1-\alpha/2, \phi_E) \sqrt{ \frac{2MSE}{r}} \tag{7.24}\]
두 평균의 차이 \(\delta_{ij}\) 에 대한 가설을 검정은 유의 수준 \(\alpha\)에서 다음과 같은 조건을 만족하면 위의 귀무가설을 기각한다.
\[ \left | \bar {x}_{i.} - \bar {x}_{j.} \right | > t(1-\alpha/2, \phi_E) \sqrt{ \frac{2MSE}{r}} \tag{7.25}\]
식 7.25 에서 검정을 위한 조건의 우변을 최소유의차(least significant difference; LSD) 라고 부른다. 두 수준의 차이가 유의하려면 두 평균 차이의 절대값이 최소한 최소유의차의 값보다 커야한다.
\[ \text{LSD} =t(1-\alpha/2, \phi_E) \sqrt{ \frac{2MSE}{r}} \]
일원배치 계획에서 수준의 개수가 a 개 인 경우 처리 수준들의 차이에 대하여 비교를 한다면 \(a \choose 2\) 개의 가설검정을 수행해야 한다. 예를 들어 처리 수준이 3개 있는 경우 다음과 같이 3개의 조합에 대하여 가설 검정을 수행할 수 있다.
\[ H_{01}: \mu_1 = \mu_2, \quad H_{02}: \mu_2 = \mu_3, \quad H_{03}: \mu_3 =\mu_1 \tag{7.26}\]
가설검정에서 사용되는 유의수준(significance level, \(\alpha\))에 대하여 생각해 보자. 지금까지 가설검정을 수행할 때 유의수준 5% 라는 말을 사용해 왔는데 이것이 무슨 의미를 가지는지 알아보자.
유의수준 5%라는 것은 수행하는 가설검정에서 귀무가설이 옳은 경우에 기각하는 확률을 말한다. 예를 들어 식 7.26 의 3개의 검정에 대하여 각각 t-검정을 수행하는 경우 귀무가설이 옳은데 우연하게 자료가 극단적으로 나와서 귀무가설을 기각하고 대립가설을 채택하는 확률이 유의수준이며 보통 5%를 사용한다. 이러한 오류를 제 1종의 오류(Type I error; false discovery error;false positve error)라고 한다.
위 식 7.26 에서 처럼 3개의 가설 검정을 동시에 실시한다면 각각의 가설검정에서 제 1 종의 오류를 범할 확률은 5%이다. 그런데 3개의 가설 검정을 동시에 실행하므로 다음과 같이 3개의 검정을 합쳐서 다음과 같은 확률에 관심이 있을 수 있다.
3개의 가설검정을 동시에 수행할 때 제 1종의 오류가 최소한 1번 발생할 확률은 얼마인가?
세 개의 가설검정을 동시에 수행하는 경우 세 검정 모두 제 1 종의 오류를 범하거나 두 개 또는 하나의 검정에서 제 1 종의 오류를 범할 사건의 확률은 얼마나 될까? 5%보다 작을까 아니면 클까? 또는 5%인가? 간단한 확률 공식을 이용하여 알아보자.
일단 두 개의 검정 \(H_{01}\) 과 \(H_{02}\)을 각각 유의수준 \(\alpha=0.05\)로서 동시에 수행 한다고 가정하고 다음과 같은 사건을 정의한다.
각 검정에서 제 1 종의 오류를 범할 확률을 \(\alpha\)라고 가정하자.
\[ P( A_1 ) = P(A_2) = \alpha =0.05 \]
이제 두개의 가설검정을 동시에 수행하는 경우 제 1 종의 오류를 최소한 1번 범하는 사건은 \(P(A_1 \cup A_2)\) 이며 여사건의 확률공식을 이용하면 다음과 같이 나타낼 수 있다.
\[ P( A_1 \cup A_2 ) = 1- P(A_1^c \cap A^c_2 ) \]
여기서 우리는 \(P(A_1^c)=P(A_2^c)=1-0.05=0.95\)를 알 수 있지만 두 사건의 교집합에 대한 확률은 계산하기 쉽지 않다. 왜냐하면 두 사건 \(A_1\)과 \(A_2\)가 일반적으로 독립이 아니어서 두 확률의 곱으로 쉽게 나타낼 수 없다.
만약에 두 사건이 독립이라면 다음과 같은 결과가 나온다. 즉 두 개의 독립인 가설검정을 동시에 수행하는 경우 최소한 1번의 제 1 종의 오류를 범하는 사건의 학률은 0.0975로 5%의 두 배 정도가 된다.
\[ P( A_1 \cup A_2 ) = 1- P(A_1^c \cap A^c_2 ) =1-P(A_1^c)P(A^c_2 ) = 1-(1-0.05)^2 = 0.0975 > 0.05 \]
만약 \(k\) 개의 독립인 가설검정을 동시에 수행하는 경우 제 1 종의 오류를 최소한 1번 이라도 범하는 사건의 학률은 \(1-(1-0.05)^k\)으로 급격하게 증가한다. 예를 들어 \(k=6\)인 경우 26.5% 로 5%의 5 배가 된다. 여기서 유의할 점은 이러한 결과는 모든 가설검정이 독립이고 여러 개의 가설검정들을 동시에 고려하는 경우이다.
즉, 두 개 이상의 가설검정을 동시에 고려해서 제 1 종의 오류를 최소한 1번 범할 경우를 오류라고 한다면 그 확률은 고려하는 검정의 개수가 증가함에 따라 빠르게 커진다.
이렇게 두 개 이상의 가설검정을 동시에 고려해서 계산하는 오류의 확률을 실험단위 오류(Experiment-wise error 또는 Family-wise error)라고 하며 반대로 가설검정을 동시에 고려하지 않고 개별적로 생각하는 오류를 개별단위 오류(Individual-wise error)라고 한다.
임상실험에서 신약(처리 1)의 효과가 위약(처리 2)보다는 우월하다는 사실을 입증하는 것이 일반적이다. 그런데 기존의 약(처리 3)보다 우월하다는 사실을 동시에 입증하려고 하는 경우도 있다. 이러한 경우 다음과 같은 두 개의 가설을 동시에 수행해야 한다.
\[ H_{01}: \mu_1 = \mu_2, \quad H_{02}: \mu_1 = \mu_3 \]
3개의 수준(신약, 위약, 기본의 약)을 가진 일원배치법으로 실험을 수행한 경우 첫 번쨰 가설 \(H_{01}\)은 \({\bar x}_{1.} - {\bar x}_{2.}\)를 이용하고 두 번쨰 가설 \(H_{02}\)은 \({\bar x}_{1.} - {\bar x}_{3.}\)을 이용하여 가설검정을 한다.
이러한 경우 각 검정에 대하여 유의 수준을 5% (개별단위 오류를 범할 확률이 5%) 라고 해도 실험단위 오류를 범할 확률은 5% 보다 크다.
여기서 유의할 점은 두 개의 가설에 대한 검정 통계량 \({\bar x}_{1.} - {\bar x}_{2.}\)과 \({\bar x}_{1.} - {\bar x}_{3.}\) 는 독립이 아니므로(why?) 실험단위 오류를 범할 확률은 5% 보다 크고 9.75% 보다는 작다.
다시 실험 단위 오류의 계산으로 돌아가서 만약에 두 사건이 독립이 아닌 경우에 실험적 오류를 통제할 수 있는, 즉 5%보다 작거나 같게 하는 방법에 대해서 알아보자 두 사건이 독립이 아닌 일반적인 경우에 확률 공식을 이용하여 다음과 같은 부등식을 얻을 수 있다.
\[ P( A_1 \cup A_2 ) \le P( A_1 ) + P( A_2 ) = (2)(0.05) = 0.1 \]
위의 결과를 보면 만약에 두 개의 가설검정을 동시에 수행하는 경우 각 가설검정에 대한 개별단위의 제 1 종 오류에 대한 확률을 반으로 줄이면(0.05/2=0.025) 실험적 오류가 5%보다 작거나 같게 된다.
\[ P( A_1 \cup A_2 ) \le P( A_1 ) + P( A_2 ) = (2)(0.05/2) = 0.05 \]
위에서 보인 같은 논리로서 \(k\) 개의 가설검정을 동시에 수행하는 경우 각 가설검정에 대한 개별적 1종 오류의 확률을 \(k\)배 줄이면(\(0.05/k\)) 실험단위 오류가 5%보다 작거나 같게 된다.
\[ P( A_1 \cup A_2 \cup ... \cup A_k ) \le (k)(0.05/k) = 0.05 \]
여기서 한 가지 유의할 점은 만약 두 개의 가설이 완전히 종속이거나(\(A_1 = A_2\)) 거의 종속이면 실험적 오류는 거의 변하지 않는다. 따라서 개별단위 1종 오류에 대한 수정은 거의 필요하지 않다.
\[ P( A_1 \cup A_2 ) = 1- P(A_1^c \cap A^c_2 ) \approx 1-P(A_1^c) = 0.05 \]
이렇게 실험단위 오류를 통제하기 위하여(5%보다 작거나 같게) 각 가설에 대한 개별단위 1 종 오류의 확률(유의수준)를 보정하는 방법을 다중비교(mutiple comparison) 라고 한다.
위에서 제시한 개별단위 1종 오류를 \(k\)배로 줄이는(0.05/k) 방법을 특별하게 본페로니 수정(Bonferroni correction)이라고 부른다. 본페로니 수정은 가장 보수적인 수정(most conservative correction)이라고 불리는데 그 이유는 실험적 오류가 가질 수 있는 가장 큰 값을 가정하고 보정하기 때문에 각각 수정한 개별단위 오류에 대한 유의수준이 너무 작게 되어(\(0.05/k\)) 귀무가설의 기각이 매우 힘들기 떄문이다.
만약 \(k\)개의 가설 검정에 본페로니 수정을 적용한다면 신뢰구간과 가설검정은 다음과 같이 수정된다.
두 수준 평균의 차이 \(\delta_{ij} = \mu_i - \mu_j\) 에 대한 본페로니 수정 신뢰구간은 다음과 같이 주어진다.
\[ ( \bar {x}_{i.} - \bar {x}_{j.}) \pm t(1-\alpha/(2k), \phi_E) \sqrt{ \frac{2MSE}{r}} \tag{7.27}\]
두 평균의 차이 \(\delta_{ij}\) 에 대한 가설을 본페로니 수정 검정은 다음과 같은 조건을 만족하면 귀무가설을 기각한다.
\[ \left | \bar {x}_{i.} - \bar {x}_{j.} \right | > t(1-\alpha/(2k), \phi_E) \sqrt{ \frac{2MSE}{r}} \tag{7.28}\]
기각역에 본페로니 수정을 하는 것은 윈래의 p-값에 가설의 개수 \(k\) 를 곱하여 수정 p-값을 사용하는 것과 같다.
\[ \text{Bonferoni adjusted p-value } = k \times \text{ unadjusted p-value } \tag{7.29}\]
일반적으로 각 가설검정들은 완전히 독립도 아니고 또한 완전한 종속도 아니다. 따라서 실험단위 오류는 각 가설 검정들이 어떻게 확률적으로 관련되어 있느냐에 따라 매우 달라진다. 이러한 이유로 인하여 다중비교의 방법은 매우 다양하며, 선택한 방법에 따라서 검정의 결과도 매우 달라질 수있는 사실에 유의해야 한다. 다중비교의 방법을 선택하는 것은 매우 어려운 일이다.
가설이 2개 이상 있는 경우 실험단위의 오류의 확률을 제어해야 하는지에 대한 판단은 상황에 따라서 달라진다.
앞에서 살펴본 임상실험의 예와 같이 중요한 의사 결정을 동시에 수행하는 2개 이상의 검정 결과에 따라서 해야할 경우 주로 다중 비교를 적용한다.
또한 다중 비교 방법은 실험의 설계와 목적에 따라서 많은 방법들이 존재한다. 주어진 실험 계획과 목적에 부합하는 다중 비교법을 선택해야 한다.
반면 탐색적인 목적으로 여러 개의 가설 검정을 동시에 수행하는 경우에는 다중비교를 적용하지 않거나 다중 비교보다 더 유연한 False Discovery Rate 방법(참조) 을 사용한다.