부록 B — 일원배치 모형과 최소제곱법

B.1 최소제곱법과 제약조건

이제 일원배치법에 대한 통계적 모형에서 모수에 대한 추정을 생각해 보자.

\[ y_{ij} = \mu + \alpha_i + e_{ij} \tag{B.1}\]

추정해야할 모수는 전체 평균 \(\mu\)와 각 그룹의 처리 효과 \(\alpha_i\) 그리고 분산 \(\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{B.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{B.3}\]

위의 방정식에서 첫 번째 방정식은 다른 \(a\)개의 방정식을 모두 합한 방정식과 같다. 따라서 모수는 \(a+1\)개이지만 실제 방정식의 개수는 \(a\)개이므로 유일한 해가 얻어지지 않는다. 따라서 유일한 해를 구하려면 하나의 제약조건이 필요하며 일반적으로 다음과 같은 두 개의 조건 중 하나를 사용한다.

B.1.1 set-to-zero condition

첫 번째 효과 \(\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{B.4}\]

B.1.2 sum-to-zero condition

처리들의 효과의 합은 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{B.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\) 에 대한 추정량은 다음과 같이 주어진다.

\[\begin{equation*} \hat \sigma_E^2 = \frac{ \sum_i \sum_j (y_{ij} - \hat \mu -\hat \alpha_i )^2}{a(r-1)} \end{equation*}\]

B.2 선형모형과 제약 조건

일원배치 모형 식 B.1 를 다음과 같은 벡터를 이용한 선형모형(linear regression model) 형태로 나타내고자 한다.

\[ \pmb y = \pmb X \pmb \beta +\pmb e \tag{B.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{B.7}\]

이제 위에서 논의한 최소제곱법을 선형 모형 식 B.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{B.8}\]

최소제곱법의 기준을 만족하는 계수 \(\pmb \beta\)는 다음과 같은 정규방정식(normal equation)의 해(solution)이다.

\[ \pmb X^t \pmb X \pmb \beta = \pmb X^t \pmb y \tag{B.9}\]

정규방정식 식 B.9 을 일워배치의 선형모형식 식 B.7 에 나타난 \(\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{B.10}\]

정규방정식 식 B.10 는 위에서 구한 최소제곱법에서 유도된 방정식 식 B.3 과 같다.

여기서 유의할 점은 선형모형식 식 B.7 의 계획행렬 \(\pmb X\) 가 완전 계수(full rank) 행렬이 아니다. 계획행렬 \(\pmb X\)의 첫 번째 열은 다른 열을 합한 것과 같다. 또한 정규 방정식 식 B.10 에서 \(\pmb X^t \pmb X\) 행렬도 완전계수 행렬이 아니다. 따라서 \(\pmb X^t \pmb X\) 행렬의 역행렬은 존재하지 않는다.

이러한 이유로 모수에 대한 유일한 추정량이 존재하지 않기 때문에 앞에서 언급한 제약 조건을 고려해야 정규방정식을 풀 수 있다.

B.2.1 Set-to-zero 조건에서의 모형과 최소제곱 추정량

만약 Set-to-zero 조건을 가정한다면 모수에서 \(\alpha_1\)을 제외하고 선형모형식 식 B.7 를 다음과 같이 다시 표현할 수 있다.
효과 \(\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{B.11}\]

이제 수정된 모형식 식 B.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{B.12}\]

위의 정규방정 식 B.12 를 풀면 위에서 언급한 sum-to-zero 조건에서 구해지는 모수의 추정량 식 B.4 를 얻을 수 있다.

B.2.2 Sum-to-zero 조건에서의 모형과 최소제곱 추정량

이제 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{B.13}\]

이제 수정된 모형식 식 B.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{B.14}\]

위의 정규방정 식 B.14 를 풀면 위에서 언급한 sum-to-zero 조건에서 구해지는 모수의 추정량 식 B.5 를 얻을 수 있다.

B.3 추정 가능한 함수

B.3.1 일원배치법에 추정가능한 모수

앞 절에서 보았듯이 일원배치법을 선형 모형식으로 표현하는 경우 평균에 대한 모수는 모두 \(a+1\) 개가 있다.

\[ \mu, \alpha_1, \alpha_2, \cdots, \alpha_a \]

하지만 모형식에서 계획행렬 \(\pmb X\)가 완전 계수 행렬이 아니기 때문에 1개의 제약 조건을 가정하고 모수를 추정하였다. 하지만 제약 조건이 달라지면 각 모수의 추정량이 달라지기 때문에 각 모수는 유일한 값으로 추정이 불가능하다.

이렇게 각 모수들은 제약 조건에 따라서 유일하게 추정이 불가능하지만 앞 절에서 보았듯이 \(\mu + \alpha_i\) 에 대한 추정량은 제약조건에 관계없이 표본 평균 \(\bar y_{i.}\)으로 동일하게 추정되어 진다.

그러면 어떤 모수들은 유일하게 추정이 불가능하고 어떤 모수들이 유일하게 추정이 가능할까?

이제 제약조건이 달라도 유일하게 추정이 가능한 모수들의 형태를 살펴보자.

B.3.2 추정가능한 모수의 함수

선형모형 \(\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{B.15}\]

즉 추정가능한 모수의 조합 \(\psi = \pmb c^t \pmb \beta\)에서 계수 벡터 \(\pmb c\) 는 계획행렬에 있는 행들의 선형 조합으로 표시되어야 한다는 것이다. 이렇게 유일하게 추정이 가능한 모수의 조합을 추정가능한 함수(estimable function)이라고 한다.

B.3.3 예제

2개의 수준이 있고 반복이 2번 있는 일원배치 \((a=2,r=2)\) 에 대한 선형모형 ?eq-lm2을 생각해보자. 이 경우 계획행렬 \(\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 \]

위의 식 식 B.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{B.16}\]

이제 \(\pmb X^t \pmb a\) 는 계획행렬 \(\pmb X\)에 있는 유일한 행들의 선형조합임을 알 수 있다.

위의 식 식 B.16 에서 유의할 점은 벡터 \(\pmb a=[a_1 ~a_2~a_3~a_4]^t\)는 임의로 주어진 벡터이다.

식 B.16 에서 \(a_1=1\), \(a_2=1\) 인 경우는 \(a_1=2\), \(a_2=0\) 인 경우와 동일하다.

따라서 유일하게 추정 가능한 모수의 선형조합 \(\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{B.17}\]

  • 처리의 효과를 나타내는 모수 \(\alpha_i\)는 추정이 불가능하다.

첫 번째 처리에 대한 효과 모수 \(\alpha_1\) 를 선형조합으로 나타내면

\[ \alpha_1 = c_0 \mu + c_1 \alpha_1 + c_2 \alpha_2 = (0) \mu + (1) \alpha_1 + (0) \alpha_2 \]

따라서 조건 식 B.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_i\)는 추정이 가능하다.

모수 조합 \(\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 \]

따라서 조건 식 B.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\)는 추정이 가능하다.

\[ \alpha_1 -\alpha_2= c_0 \alpha_0 + c_1 \alpha_1 + c_2 \alpha_2 = (0) \pmb u + (1) \alpha_1 + (-1) \alpha_2 \]

따라서 조건 식 B.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} \]

B.4 R 실습

B.4.1 예제 3.1

4개의 서로 다른 원단업체에서 직물을 공급받고 있다. 공급한 직물의 긁힘에 대한 저항력을 알아보기 위하여 각 업체마다 4개의 제품을 랜덤하게 선택하여 (\(a=4\), \(r=4\)) 일원배치법에 의하여 마모도 검사을 실시하였다.

B.4.2 자료의 생성

company<- as.factor(rep(c(1:4), each=4))
response<- c(1.93, 2.38, 2.20, 2.25,
             2.55, 2.72, 2.75, 2.70,
             2.40, 2.68, 2.32, 2.28,
             2.33, 2.38, 2.28, 2.25)
df31<- data.frame(company=company, response= response)
df31
   company response
1        1     1.93
2        1     2.38
3        1     2.20
4        1     2.25
5        2     2.55
6        2     2.72
7        2     2.75
8        2     2.70
9        3     2.40
10       3     2.68
11       3     2.32
12       3     2.28
13       4     2.33
14       4     2.38
15       4     2.28
16       4     2.25

각 수준에 대한 표보 평균을 구해보자.

df31s <- df31 %>% group_by(company)  %>%  summarise(mean=mean(response), median= median(response), sd=sd(response), min=min(response), max=max(response))
df31s
# A tibble: 4 × 6
  company  mean median     sd   min   max
  <fct>   <dbl>  <dbl>  <dbl> <dbl> <dbl>
1 1        2.19   2.22 0.189   1.93  2.38
2 2        2.68   2.71 0.0891  2.55  2.75
3 3        2.42   2.36 0.180   2.28  2.68
4 4        2.31   2.30 0.0572  2.25  2.38

B.4.3 선형모형의 적합(set-to-zero)

이제 자료를 다음과 같은 선형 모형으로 적합해 보자. 선형 모형의 적합은 lm() 함수를 사용한다.

\[ y_{ij} = \mu + \alpha_i + e_{ij} \]

여기서 선형식의 모수와 R의 변수는 다음과 같은 관계를 가진다,

선형식의 모수 R의 변수
\(\mu\) (Intercept)
\(\alpha_1\) company1
\(\alpha_2\) company2
\(\alpha_3\) company3
\(\alpha_4\) company4
fit1 <- lm(response~company,data=df31)
summary(fit1)

Call:
lm.default(formula = response ~ company, data = df31)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2600 -0.0700  0.0150  0.0625  0.2600 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.19000    0.07050  31.062 7.79e-13 ***
company2     0.49000    0.09971   4.914 0.000357 ***
company3     0.23000    0.09971   2.307 0.039710 *  
company4     0.12000    0.09971   1.204 0.251982    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.141 on 12 degrees of freedom
Multiple R-squared:  0.6871,    Adjusted R-squared:  0.6089 
F-statistic: 8.785 on 3 and 12 DF,  p-value: 0.002353

위에서 적합한 결과를 보면 평균 \(\mu\)와 4개의 처리 \(\alpha_1\), \(\alpha_2\), \(\alpha_3\), \(\alpha_4\) 가 모형에 있지만 모수의 추정량은 평균(intercept)과 3개의 모수(company2, company3, company4)만 추정량이 주어진다.

R 에서 옵션을 지정하지 않고 함수 lm()으로 선형모형을 적합하는 경우 set-to-zero 조건을 적용하며 자료에 나타난 처리의 수준들 중 순위가 가장 낮은 수준의 효과를 0으로 지정한다 (company1=0 ). set-to-zero 조건을 강제로 지정하려면 다음과 같은 명령문을 먼저 실행한다.

options(contrasts=c("contr.treatment", "contr.poly"))

위의 결과를 보면 (Intercept)에 대한 추정량이 첫 번째 처리 company1의 평균과 같은 것을 알 수 있다.

set-to-zero 조건에서의 계획행렬은 다음과 같이 볼 수 있다.

model.matrix(fit1)
   (Intercept) company2 company3 company4
1            1        0        0        0
2            1        0        0        0
3            1        0        0        0
4            1        0        0        0
5            1        1        0        0
6            1        1        0        0
7            1        1        0        0
8            1        1        0        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
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$company
[1] "contr.treatment"

이제 각 처리 평균에 대한 추정값 \(\widehat{\mu+ \alpha_i}\)을 구해보자.

emmeans(fit1, "company")
 company emmean     SE df lower.CL upper.CL
 1         2.19 0.0705 12     2.04     2.34
 2         2.68 0.0705 12     2.53     2.83
 3         2.42 0.0705 12     2.27     2.57
 4         2.31 0.0705 12     2.16     2.46

Confidence level used: 0.95 

이 경우 처리 평균에 대한 추정값은 산술 평균과 동일하게 나온다.

B.4.4 선형모형의 적합 (sum-to-zero)

이제 일원배치 모형에서 sum-to-zero 조건을 적용하여 모수를 추정해 보자. sum-to-zero 조건을 적용하려면 다음과 같은 명령어를 실행해야 한다.

options(contrasts=c("contr.sum", "contr.poly"))

이제 다시 선형모형을 적합하고 추정결과를 보자.

fit2 <- lm(response~company,data=df31)
summary(fit2)

Call:
lm.default(formula = response ~ company, data = df31)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2600 -0.0700  0.0150  0.0625  0.2600 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.40000    0.03525  68.081  < 2e-16 ***
company1    -0.21000    0.06106  -3.439 0.004901 ** 
company2     0.28000    0.06106   4.586 0.000626 ***
company3     0.02000    0.06106   0.328 0.748892    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.141 on 12 degrees of freedom
Multiple R-squared:  0.6871,    Adjusted R-squared:  0.6089 
F-statistic: 8.785 on 3 and 12 DF,  p-value: 0.002353

이제 sum-to-zero 조건에 따라서 위의 set-to-zero 결과와 모수의 추정값이 다르게 나타나는 것을 알 수 있다. 마지막 모수 company4(\(\alpha_4\))는 sum-to-zero 조건을 이용하여 다음과 같은 관계를 이용하여 구할 수 있다.

\[ \alpha_4 = -(\alpha_1 + \alpha_2 + \alpha_3) \]

sum-to-zero 조건에서의 계획행렬은 다음과 같이 볼 수 있다.

model.matrix(fit2)
   (Intercept) company1 company2 company3
1            1        1        0        0
2            1        1        0        0
3            1        1        0        0
4            1        1        0        0
5            1        0        1        0
6            1        0        1        0
7            1        0        1        0
8            1        0        1        0
9            1        0        0        1
10           1        0        0        1
11           1        0        0        1
12           1        0        0        1
13           1       -1       -1       -1
14           1       -1       -1       -1
15           1       -1       -1       -1
16           1       -1       -1       -1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$company
[1] "contr.sum"

이제 각 처리 평균에 대한 추정값 \(\widehat{\mu+ \alpha_i}\)을 구해보면 set-to-zero 조건에서의 추정값과 동일함을 알 수 있다.

emmeans(fit2, "company")
 company emmean     SE df lower.CL upper.CL
 1         2.19 0.0705 12     2.04     2.34
 2         2.68 0.0705 12     2.53     2.83
 3         2.42 0.0705 12     2.27     2.57
 4         2.31 0.0705 12     2.16     2.46

Confidence level used: 0.95 

B.4.5 분산분석

분산분석의 결과는 어떠한 제약 조건에서도 동일하다.

res1 <- anova(fit1)
res1
Analysis of Variance Table

Response: response
          Df Sum Sq  Mean Sq F value   Pr(>F)   
company    3 0.5240 0.174667  8.7846 0.002353 **
Residuals 12 0.2386 0.019883                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res2<- anova(fit2)
res2
Analysis of Variance Table

Response: response
          Df Sum Sq  Mean Sq F value   Pr(>F)   
company    3 0.5240 0.174667  8.7846 0.002353 **
Residuals 12 0.2386 0.019883                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1