An Analysis of Covariance Example
Whiteside 자료
1960년대에 영국 건축 연구소(UK Building Research Station)의 데릭 화이트사이드(Derek Whiteside)가 수집한 자료 (Hand et al., 1994 ).
잉글랜드 남동부에 있는 화이트사이드의 집에서 캐비티 벽 단열재를 설치하기 전과 후의 두 번의 “난방 시즌” 동안 주간 가스 소비량과 평균 외부 기온을 기록
실험 목적: 단열재가 가스 소비에 미치는 영향을 평가
# A tibble: 56 × 3
Insul Temp Gas
<fct> <dbl> <dbl>
1 Before -0.8 7.2
2 Before -0.7 6.9
3 Before 0.4 6.4
4 Before 2.5 6
5 Before 2.9 5.8
6 Before 3.2 5.8
7 Before 3.6 5.6
8 Before 3.9 4.7
9 Before 4.2 5.8
10 Before 4.3 5.2
# ℹ 46 more rows
Insul: Before와 After의 두 수준을 갖는 요인 변수Temp: 섭씨 단위의 주간 평균 외부 기온Gas: 1,000 입방피트 단위의 주간 가스 소비량
먼저 자료를 단열재 설치 전과 후로 나누어 두 개의 패널에 최소제곱 직선을 함께 그려 보자.
ggplot (data= MASS:: whiteside, aes (x= Temp, y= Gas), facets= ~ Insul) + stat_smooth (method = "lm" ) + geom_point () + facet_wrap (~ Insul) + labs (x = "Average external temperature (deg. C)" ,y = "Gas consumption (1000 cubic feet)" )
주어진 온도 범위 내에서는 직선 모형이 적절한 것으로 보인다.
그래프로 미루어 볼 때, 단열재가 동일한 외부 온도에서 가스 소비량을 줄이지만 직선의 기울기, 즉 외부 온도가 떨어질수록 가스 소비량이 증가하는 속도에도 영향을 미치는 것으로 보인다.
선형모형과 lm()
이러한 문제를 정량적으로 탐색하기 위해 lm을 기본 함수로 하는 선형 모형을 적용한다.
lm (formula, data, subset, weights, na.action)
formula: 윌킨슨-로저스 모형 공식 (유일한 필수 인수)data: 데이터 프레임 (선택 사항)weights: 양의 가중치 벡터 (균등하지 않은 가중치가 필요한 경우)subset: 사용할 데이터의 하위 집합을 지정하는 인덱스 벡터 (기본적으로 모든 항목이 사용됨)na.action: 결측치를 처리하는 방법을 지정하는 함수 (기본값은 na.omit, 즉 생략. 적합치에 NA가 포함되어야 하는 경우는 na.exclude 사용)
그래프의 직선을 각각 적합하기 위해 아래와 같이 lm을 사용할 수 있다.
<- lm (Gas ~ Temp, data = MASS:: whiteside, subset = Insul == "Before" )# gasB <- lm(Gas ~ 1+Temp, data = MASS::whiteside, subset = Insul == "Before") # same result as above <- update (gasB, subset = Insul == "After" ) # modify the fitted model
적합된 모형 객체는 lm 클래스에 속한다. 객체에 대한 추가 작업을 수행하는 generic 함수는 다음과 같다.
print: 간단한 출력summary: 통상적인 회귀분석 결과 출력,coef (or coefficients): 회귀계수 추출resid (or residuals): 잔차 추출
fitted (or fitted.values): 적합치 추출deviance: 잔차제곱합 추출anova: 분산표를 순차적 분석 또는 여러 계층적 모형을 비교predict: 새로운 데이터에 대한 평균 예측, 선택적으로 표준 오차 포함plot: 진단 플롯
lm 객체분석 결과 요약
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.85383 0.11842 57.88 <2e-16 ***
## Temp -0.39324 0.01959 -20.08 <2e-16 ***
## ---
##
## Residual standard error: 0.2813 on 24 degrees of freedom## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.72385 0.12974 36.41 < 2e-16 ***
## Temp -0.27793 0.02518 -11.04 1.05e-11 ***
## ---
##
## Residual standard error: 0.3548 on 28 degrees of freedom
표본 분산
<- deviance (gasB) / gasB$ df.resid) # direct calculation <- summary (gasB)$ sigma^ 2 ) # alternative
<- lm (Gas ~ Insul / Temp - 1 , data = MASS:: whiteside) summary (gasBA)
Call:
lm(formula = Gas ~ Insul/Temp - 1, data = MASS::whiteside)
Residuals:
Min 1Q Median 3Q Max
-0.97802 -0.18011 0.03757 0.20930 0.63803
Coefficients:
Estimate Std. Error t value Pr(>|t|)
InsulBefore 6.85383 0.13596 50.41 <2e-16 ***
InsulAfter 4.72385 0.11810 40.00 <2e-16 ***
InsulBefore:Temp -0.39324 0.02249 -17.49 <2e-16 ***
InsulAfter:Temp -0.27793 0.02292 -12.12 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.323 on 52 degrees of freedom
Multiple R-squared: 0.9946, Adjusted R-squared: 0.9942
F-statistic: 2391 on 4 and 52 DF, p-value: < 2.2e-16
이차식 회귀
<- lm (Gas ~ Insul/ (Temp + I (Temp^ 2 )) - 1 , data = MASS:: whiteside)summary (gasQ)$ coef
Estimate Std. Error t value Pr(>|t|)
InsulBefore 6.759215179 0.150786777 44.826312 4.854615e-42
InsulAfter 4.496373920 0.160667904 27.985514 3.302572e-32
InsulBefore:Temp -0.317658735 0.062965170 -5.044991 6.362323e-06
InsulAfter:Temp -0.137901603 0.073058019 -1.887563 6.489554e-02
InsulBefore:I(Temp^2) -0.008472572 0.006624737 -1.278930 2.068259e-01
InsulAfter:I(Temp^2) -0.014979455 0.007447107 -2.011446 4.968398e-02
분산분석
더 간단한 모형으로 기울기가 같은 회귀모형을 생각할 수 있다.
이 모형을 적합하고 anova 함수를 사용하여 수준별 기울기를 갖는 회귀 모형 내에서 검정할 수 있다. (왜 그런가?)
# MASS p. 143 <- lm (Gas ~ Insul + Temp, data = MASS:: whiteside)anova (gasPR, gasBA)
Analysis of Variance Table
Model 1: Gas ~ Insul + Temp
Model 2: Gas ~ Insul/Temp - 1
Res.Df RSS Df Sum of Sq F Pr(>F)
1 53 6.7704
2 52 5.4252 1 1.3451 12.893 0.0007307 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
모형 공식
선형 모형은 응답 벡터 \(\boldsymbol{y}\) 와 설명 변수 행렬 또는 모형 행렬 \(X\) 로 결정된다.
모형 공식은 좌변에 응답변수를, 우변에는 특정 규칙에 따라 모형 행렬을 생성하는 방법에 대한 지시사항을 적는다.
예제: 3변수 다중선형모형
모형 공식
y ~ x1 + x2 + x3대수적 기술
\[
y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \epsilon_i,
\quad
i= 1, 2, \dotsc, n
\]
모형 행렬
\[
X = [\mathbf{1} ~ \mathbf{x}_1 ~ \mathbf{x}_2 ~ \mathbf{x}_3]
\]
절편 \(\beta_0\) 는 모형공식에서 암묵적으로 가정된다 (y ~ 1 + x1 + x2 + x3).
y ~ x1 + x2 + x3 - 1로 절편을 제거할 수 있다.
예제: 1원배치 모형
설명변수 a가 \(k\) 개 수준을 갖는 범주형 변수 (factor).
대수적 기술
\[
y_{ij} = \mu + \alpha_j + \epsilon_{ij},
\quad
i=1,2,\dotsc,n_j;
~
j=1,2,\dotsc, k;
~
n = \sum_j n_j
\]
벡터형식 및 모형행렬
\[
\begin{split}
\mathbf{y} &= X\boldsymbol{\beta} + \boldsymbol{\epsilon},
\\
\boldsymbol{\beta} &= [\mu, \alpha_1, \alpha_2, \dotsc, \alpha_k]^T = [\mu, \boldsymbol{\alpha}^T]^T,
\\
\mathbf{y} &= [y_{11}, \dotsc, y_{n_1,1}, \dotsc, y_{1k}, \dotsc, y_{n_k,k}]^T,
\\
X &= [\mathbf{1} ~ X_a],
\end{split}
\]
\(X_a\) 는 \(n\times k\) 행렬로 각 행 이 a의 소속 수준에 해당하는 열에만 1이 있는 벡터(dummy variable, one-hot vector)로 구성된 이진행렬 (\(n = \underset{j=1} { \stackrel{k}{\sum} } n_j\) ).
예: \(n_1=2\) , \(n_2=1=n_3\) , \(k=3\) \[
X_a = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
\]
\[
X_a = \begin{bmatrix} \mathbf{1}_{n_1} & O & \cdots & O \\ O & \mathbf{1}_{n_2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & O \\ O & \cdots & O & \mathbf{1}_{n_k}\end{bmatrix}
\]
식별가능 모형 1: 제약식
\[
y_{ij} = \alpha_j + \epsilon_{ij},
\quad
i=1,2,\dotsc,n_j;
~
j=1,2,\dotsc, k
\] 즉, \(\beta_0 = \mu = 0\) 으로 제약. 벡터형식: \(\mathbf{e}_1^T \boldsymbol{\beta} = 0\) ; \(\mathbf{e}_1 = [1, 0, \dotsc, 0]^T\) .
모형행렬? 랭크?
모형 공식: y ~ a - 1
set.seed (2024-03-01 )<- data.frame (a = factor (rep (1 : 3 , each= 3 )), y = rnorm (9 , rep (2 : 4 , each= 3 ), 0.1 ))<- lm (y ~ a - 1 , data = dat)coef (mod)
a1 a2 a3
1.986017 2.893121 4.082296
a1 a2 a3
1 1 0 0
2 1 0 0
3 1 0 0
4 0 1 0
5 0 1 0
6 0 1 0
7 0 0 1
8 0 0 1
9 0 0 1
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
차원 \(k \times (k-1)\) 인 최대랭크(full rank) 행렬 \(C_a\) 에 대해 정방행렬 \([\mathbf{1}_k~C_a]\) 가 정칙(nonsingular)이라고 가정하면 \[
\mathbf{1}_k\mu + \boldsymbol{\alpha} = [\mathbf{1}_k~C_a]\begin{bmatrix} \mu^* \\ \boldsymbol{\alpha}^* \end{bmatrix}
\] 인 \(\mu^*\) , \(\boldsymbol{\alpha}^*\) 가 항상 존재한다.
따라서, \[
\begin{split}
[\mathbf{1}_n~X_a]\begin{bmatrix} 0 \\ \mathbf{1}_k\mu + \boldsymbol{\alpha} \end{bmatrix}
&= X_a(\mathbf{1}_k\mu + \boldsymbol{\alpha})
\\
&= X_a[\mathbf{1}_k~C_a]\begin{bmatrix} \mu^* \\ \boldsymbol{\alpha}^* \end{bmatrix}
= [\mathbf{1}_n ~ X_a C_a]\begin{bmatrix} \mu^* \\ \boldsymbol{\alpha}^* \end{bmatrix}
.
\end{split}
\]
이때 \([\mathbf{1}_n ~ X_a C_a] = X_a[\mathbf{1}_k~C_a]\) 은 최대랭크 행렬이므로 (왜 그런가?) 이를 새로운 모형행렬로 삼으면 \(\mu^*\) , \(\boldsymbol{\alpha}^*\) 은 최소제곱법에 의해 유일하게 결정된다.
위 조건을 만족하는 행렬 \(C_a\) 를 부호화행렬(coding matrix)이라 부른다. 모형에서 \(\mu\) 는 임의이므로 \(\mu = \mu^*\) 로 두면 \[
\boldsymbol{\alpha} = C_a\boldsymbol{\alpha}^*
\] 가 된다.
부호화 행렬의 종류
처리 부호화 (treatment coding)
R의 기본 부호화. 항상 \(\alpha_1 = 0\) 으로 한다. 이렇게 하면 수준 1이 다른 모든 수준을 비교하는 표준 수준으로 취급된다.
<- lm (y ~ a, data= dat)model.matrix (mod_treatment)
(Intercept) a2 a3
1 1 0 0
2 1 0 0
3 1 0 0
4 1 1 0
5 1 1 0
6 1 1 0
7 1 0 1
8 1 0 1
9 1 0 1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
<- coef (mod_treatment)) # mu^* = (Intercept)
(Intercept) a2 a3
1.9860166 0.9071045 2.0962792
<- contr.treatment (length (levels (dat$ a))))<- Ca %*% alph.star[- 1 ]) # letting $mu = mu^*$
[,1]
1 0.0000000
2 0.9071045
3 2.0962792
헬메르트 부호화 (Helmert coding)
S-plus의 기본 부호화. 각 열의 합이 0이 되도록 한다.
<- lm (y ~ a, data= dat, contrasts = list (a = "contr.helmert" ))<- coef (mod_helmert))
(Intercept) a1 a2
2.9871445 0.4535523 0.5475757
<- contr.helmert (length (levels (dat$ a))))
[,1] [,2]
1 -1 -1
2 1 -1
3 0 2
<- Ca %*% alph.star[- 1 ]) # letting $mu = mu^*$
[,1]
1 -1.00112792
2 -0.09402339
3 1.09515132
sum (alph) # numerically zero. why?
직교다항식 부호화 (orthogonal polynomial coding)
R의 순서형 변수 에 대한 기본 부호화.
$ b <- as.ordered (dat$ a)<- lm (y ~ b, data= dat)<- coef (mod_poly))
(Intercept) b.L b.Q
2.9871445 1.4822933 0.1151547
<- contr.poly (length (levels (dat$ b))))
.L .Q
[1,] -7.071068e-01 0.4082483
[2,] -9.073800e-17 -0.8164966
[3,] 7.071068e-01 0.4082483
<- Ca %*% alph.star[- 1 ]) # letting $mu = mu^*$
[,1]
[1,] -1.00112792
[2,] -0.09402339
[3,] 1.09515132
sum (alph) # numerically zero
모수 해석
원모수 \(\mu\) , \(\alpha\) (식별불가능)와 축소된 모수 \(\mu^*\) , \(\alpha^*\) (식별가능) 사이의 관계 \[
\begin{bmatrix} \mu^* \\ \boldsymbol{\alpha}^* \end{bmatrix}
=
[\mathbf{1}_k~C_a]^{-1}
(\mathbf{1}_k\mu + \boldsymbol{\alpha})
\] 를 상기하자.
\(D = [\mathbf{1}_k~C_a]^{-1}\) 로 놓고 \(D\) 의 행을 \(\mathbf{d}_0^T, \mathbf{d}_1^T, \dotsc, \mathbf{d}_{k-1}^T\) 라 하면 \[
I_k = D[\mathbf{1}_k~C_a]
=
\begin{bmatrix} \mathbf{d}_0^T \\ \mathbf{d}_1^T \\ \vdots \\ \mathbf{d}_{k-1}^T \end{bmatrix}
[\mathbf{1}_k~C_a]
=
\begin{bmatrix}
\mathbf{d}_0^T\mathbf{1}_k & \mathbf{d}_0^TC_a \\
\mathbf{d}_1^T\mathbf{1}_k & \mathbf{d}_1^TC_a \\
\vdots & \vdots \\
\mathbf{d}_{k-1}^T\mathbf{1}_k & \mathbf{d}_{k-1}^TC_a
\end{bmatrix}
\] 이므로 \[
\mathbf{d}_0^T C_a = \mathbf{0}^T,
\quad
\mathbf{d}_0^T\mathbf{1}_k = 1,
\quad
\mathbf{d}_i^T\mathbf{1}_k = 0,
~
i=1, \dotsc, k-1
\] 임을 알 수 있다.
이제 \(D\) 를 복호화행렬이라 부르자.
어떤 벡터 \(\mathbf{v}\) 가 \(\mathbf{v}^T\mathbf{1} = 1\) 이면 평균벡터, \(\mathbf{v}^T\mathbf{1} = 0\) 이면 대조벡터라고 부르기로 하면, \(\mathbf{d}_0\) 은 평균벡터, \(\mathbf{d}_1, \dotsc, \mathbf{d}_{k-1}\) 은 대조벡터임을 알 수 있다.
역으로 \(k\times k\) 정칙행렬 \(D\) 를 첫번째 행이 평균벡터, 나머지 행이 대조벡터가 되도록 고르면 그 역행렬 \(D^{-1}\) 은 적당한 부호화행렬 \(C_a\) 에 대해 \(D^{-1} = [\mathbf{1}_k~ C_a]\) 가 된다. (왜 그런가?)
따라서 부호화행렬 \(C_a\) 의 기술과 복호화행렬 \(D\) 의 기술은 동등하다.
이제 \[
\small
\begin{bmatrix} \mu^* \\ \alpha^*_1 \\ \vdots \\ \alpha^*_{k-1} \end{bmatrix}
=
\begin{bmatrix} \mu^* \\ \boldsymbol{\alpha}^* \end{bmatrix}
=
D(\mathbf{1}_k\mu + \boldsymbol{\alpha})
=
\begin{bmatrix} \mu + \mathbf{d}_0^T\boldsymbol{\alpha} \\ \mathbf{d}_1^T\boldsymbol{\alpha} \\ \vdots \\
\mathbf{d}_{k-1}^T\boldsymbol{\alpha} \end{bmatrix}
\] 이 되어 \(\mu^*\) 은 원모수 성분 \((\mu+\alpha_1, \dotsc, \mu+\alpha_k)\) 의 \(\mathbf{d}_0^T = (d_{01}, \dotsc, d_{0k})\) 에 의한 가중평균임을, \(\alpha^*_i\) 는 \(\boldsymbol{\alpha}\) 의 \(\mathbf{d}_i\) 에 의한 대조(contrast)임을 알 수 있다. 즉, \[
\alpha^*_i = \mathbf{d}_i^T\boldsymbol{\alpha}, \quad i=1, \dotsc, k-1
\] 으로 모수 사이의 관계를 해석할 수 있다.
이때 \(\mu\) 는 임의이므로 \(\mu = \mu^*\) 으로 놓으면 가중평균 \(\mathbf{d}_0^T\boldsymbol{\alpha} = 0\) 이 되어 이를 제약조건으로 두고 \((\mu, \boldsymbol{\alpha})\) 를 추정하면 추정값 \((\hat{\mu}, \hat{\boldsymbol{\alpha}})\) 이 유일하게 정해지고 \(\widehat{\mu^*} = \hat{\mu}\) , \(\widehat{\alpha^*_i} = \mathbf{d}_i^T\hat{\boldsymbol{\alpha}}\) 가 된다.
따라서 부호화(복호화)행렬을 사용한 식별 기법은 제약식을 사용한 식별 기법의 특수한 경우임을 알 수 있다.
예: 처리 부호화 (\(\alpha_{1}=0\) 으로 설정)
<- contr.treatment (letters[1 : 5 ]) %>% print
b c d e
a 0 0 0 0
b 1 0 0 0
c 0 1 0 0
d 0 0 1 0
e 0 0 0 1
<- solve (cbind (1 , Ca)) %>% print
a b c d e
1 0 0 0 0
b -1 1 0 0 0
c -1 0 1 0 0
d -1 0 0 1 0
e -1 0 0 0 1
따라서 \(\mu^* = \mu + \alpha_1\) , \(\alpha^*_1 = \alpha_2 - \alpha_1\) , \(\dotsc\) , \(\alpha_{k-1}^* = \alpha_k - \alpha_1\) .
왜 앞 식별가능모형1과 식별가능모형2–처리부호화 결과가 같았는가?
예: 헬메르트 부호화
<- contr.helmert (letters[1 : 5 ]) %>% print
[,1] [,2] [,3] [,4]
a -1 -1 -1 -1
b 1 -1 -1 -1
c 0 2 -1 -1
d 0 0 3 -1
e 0 0 0 4
<- solve (cbind (1 , Ca)) %>% fractions %>% print
a b c d e
[1,] 1/5 1/5 1/5 1/5 1/5
[2,] -1/2 1/2 0 0 0
[3,] -1/6 -1/6 1/3 0 0
[4,] -1/12 -1/12 -1/12 1/4 0
[5,] -1/20 -1/20 -1/20 -1/20 1/5
따라서 \[
\small
\begin{split}
\mu^* &= \mu + \frac{1}{k}\sum_{j=1}^k \alpha_j
\\
\alpha^*_1 &= \frac{1}{2}(\alpha_2 - \alpha_1)
\\
\alpha^*_2 &= \frac{1}{3}(\alpha_3 - \frac{1}{2}\sum_{j=1}^2\alpha_j)
\\
& \vdots
\\
\alpha^*_{k-1} &= \frac{1}{k}(\alpha_k - \frac{1}{k-1}\sum_{j=1}^{k-1}\alpha_j)
\end{split}
\]
왜 앞 식별가능모형2–헬메르트부호화 결과에서 모수의 합이 (수치적으로) 0이었는가?
직교다항식 부호화
\(j\) 차 다항식 \(p_j(x)\) 열이 다음을 만족: \[
\begin{split}
& \sum_{i=1}^n p_j^2(x_i) = 1, \quad j = 0, 1, \dotsc, q
\\
& \sum_{i=1}^n p_j(x_i)p_k(x_i) = 0, \quad j \neq k
.
\end{split}
\]
\(x_i\) 의 열은 수직선에서 등간격으로 설정.
<- contr.poly (letters[1 : 5 ])) %>% zapsmall %>% print # L=linear, Q=quadratic, etc.
.L .Q .C ^4
[1,] -0.6324555 0.5345225 -0.3162278 0.1195229
[2,] -0.3162278 -0.2672612 0.6324555 -0.4780914
[3,] 0.0000000 -0.5345225 0.0000000 0.7171372
[4,] 0.3162278 -0.2672612 -0.6324555 -0.4780914
[5,] 0.6324555 0.5345225 0.3162278 0.1195229
poly (seq (4 , 6 , 0.5 ), degree= 4 ) %>% zapsmall # coefficients of 4th-degree orthogonal polynomial (Kennedy & Gentle (1980, pp.343--4))
1 2 3 4
[1,] -0.6324555 0.5345225 -0.3162278 0.1195229
[2,] -0.3162278 -0.2672612 0.6324555 -0.4780914
[3,] 0.0000000 -0.5345225 0.0000000 0.7171372
[4,] 0.3162278 -0.2672612 -0.6324555 -0.4780914
[5,] 0.6324555 0.5345225 0.3162278 0.1195229
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5 5 5 5
attr(,"coefs")$norm2
[1] 1.00000000 5.00000000 2.50000000 0.87500000 0.22500000 0.03214286
attr(,"degree")
[1] 1 2 3 4
attr(,"class")
[1] "poly" "matrix"
t (Ca) %*% Ca %>% fractions # orthogonal?
.L .Q .C ^4
.L 1 0 0 0
.Q 0 1 0 0
.C 0 0 1 0
^4 0 0 0 1
<- solve (cbind (1 , Ca)) %>% zapsmall %>% print
[,1] [,2] [,3] [,4] [,5]
0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
.L -0.6324555 -0.3162278 0.0000000 0.3162278 0.6324555
.Q 0.5345225 -0.2672612 -0.5345225 -0.2672612 0.5345225
.C -0.3162278 0.6324555 0.0000000 -0.6324555 0.3162278
^4 0.1195229 -0.4780914 0.7171372 -0.4780914 0.1195229
따라서 \(\mu^* = \mu + \frac{1}{k} \underset{j=1} { \stackrel{k}{\sum} } \alpha_j\) .
왜 앞 식별가능모형2–직교부호화 결과에서 모수의 합이 (수치적으로) 0이었는가?
2원 배치
2원 배치 가법 모형
# model formula ~ a + b
대수적 기술 \[
y_{ijk} = \mu + \alpha_{i} + \beta_{j} + \epsilon_{ijk},
\]
모형행렬 (랭크 부족) \[
X = [\mathbf{1}, X_a, X_b],
\quad
\mathbf{y} = [\mathbf{1}, X_a, X_b] \begin{bmatrix} \mu \\ \boldsymbol{\alpha} \\ \boldsymbol{\beta} \end{bmatrix}
\] 여기서 \(\mathbf{y}=(y_{ijk})\) , \(\boldsymbol{\alpha}=(\alpha_i)\) , \(\boldsymbol{\beta}=(\beta_j)\) 인 벡터.
축소모형의 모형행렬 : \(X^*=[\mathbf{1}, X_a C_a, X_b C_b]\)
model.matrix (y ~ a + b, data = dat2)
(Intercept) a2 a3 b2 b3
1 1 0 0 0 0
2 1 0 0 1 0
3 1 0 0 0 1
4 1 1 0 0 0
5 1 1 0 1 0
6 1 1 0 0 1
7 1 0 1 0 0
8 1 0 1 1 0
9 1 0 1 0 1
attr(,"assign")
[1] 0 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
model.matrix (y ~ a, data= dat2)
(Intercept) a2 a3
1 1 0 0
2 1 0 0
3 1 0 0
4 1 1 0
5 1 1 0
6 1 1 0
7 1 0 1
8 1 0 1
9 1 0 1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
model.matrix (y ~ b, data= dat2)
(Intercept) b2 b3
1 1 0 0
2 1 1 0
3 1 0 1
4 1 0 0
5 1 1 0
6 1 0 1
7 1 0 0
8 1 1 0
9 1 0 1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$b
[1] "contr.treatment"
절편 항이 없는 축소모형(y ~ a + b - 1)의 모형 행렬: \(X^*=[X_a ~ X_b C_b]\) .
이 경우 모형행렬 \(X^*\) 는 a와 b에 대해 비대칭이고, 순서에 의존한다 (왜 그런가?)
model.matrix (y ~ a + b - 1 , data = dat2)
a1 a2 a3 b2 b3
1 1 0 0 0 0
2 1 0 0 1 0
3 1 0 0 0 1
4 0 1 0 0 0
5 0 1 0 1 0
6 0 1 0 0 1
7 0 0 1 0 0
8 0 0 1 1 0
9 0 0 1 0 1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
model.matrix (y ~ b + a - 1 , data = dat2)
b1 b2 b3 a2 a3
1 1 0 0 0 0
2 0 1 0 0 0
3 0 0 1 0 0
4 1 0 0 1 0
5 0 1 0 1 0
6 0 0 1 1 0
7 1 0 0 0 1
8 0 1 0 0 1
9 0 0 1 0 1
attr(,"assign")
[1] 1 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$b
[1] "contr.treatment"
attr(,"contrasts")$a
[1] "contr.treatment"
2원 배치 비가법모형
~ a + b + a: b# or ~ a * b
이 모형은 \(rs\) 수준을 갖는 1원배치 모형과 사실상 같다.
대수적 기술 \[
\begin{split}
y_{ijk} &= \mu + \alpha_{i} + \beta_{j} + \gamma_{ij} + \epsilon_{ijk},
\\
&
i=1,2,\dotsc,r;
~
j=1,2,\dotsc, s;
~
k=1,2,\dotsc, n_{ij}
\end{split}
\]
모형 행렬 : \(X=[\mathbf{1}_n, X_a, X_b, X_a \odot X_b]\)
여기서 \(A \odot B\) 는 \(A\) 와 \(B\) 의 한 열씩 뽑아 원소끼리 곱해서 새로운 열 하나를 만드는 연산. \[
\begin{split}
\mathbf{y} &= [\mathbf{1}_n, X_a, X_b, X_a\odot X_b]
\begin{bmatrix} \mu \\ \boldsymbol{\alpha} \\ \boldsymbol{\beta} \\ \boldsymbol{\gamma} \end{bmatrix}
\\
&= \mathbf{1}_n\mu + X_a\boldsymbol{\alpha} + X_b\boldsymbol{\beta} + (X_a \odot X_b)\boldsymbol{\gamma}
\end{split}
\] 여기서 \(\boldsymbol{\gamma} = (\gamma_{ijk})\) 인 벡터.
\(X_a \odot X_b\) 는 a와 b로 정의되는 하위범주에 대한 접속(incidence) 행렬.
축소모형의 모형 행렬 : \(X^*=[\mathbf{1}, X_a C_a, X_b C_b, (X_a C_a)\odot (X_b C_b)]\)
열 수: \(1 + (r-1) + (s-1) +(r-1)(s-1) = rs\)
model.matrix (y ~ a * b, data = dat2)
(Intercept) a2 a3 b2 b3 a2:b2 a3:b2 a2:b3 a3:b3
1 1 0 0 0 0 0 0 0 0
2 1 0 0 1 0 0 0 0 0
3 1 0 0 0 1 0 0 0 0
4 1 1 0 0 0 0 0 0 0
5 1 1 0 1 0 1 0 0 0
6 1 1 0 0 1 0 0 1 0
7 1 0 1 0 0 0 0 0 0
8 1 0 1 1 0 0 1 0 0
9 1 0 1 0 1 0 0 0 1
attr(,"assign")
[1] 0 1 1 2 2 3 3 3 3
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
다음과 비교하자.
<- numeric ()<- model.matrix (y ~ a, data= dat2)[, - 1 ] <- model.matrix (y ~ b, data= dat2)[, - 1 ]for (i in seq_len (ncol (XbCb))) for (j in seq_len (ncol (XaCa)))<- cbind (X_star, XbCb[, i] * XaCa[, j])
[,1] [,2] [,3] [,4]
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 1 0 0 0
6 0 0 1 0
7 0 0 0 0
8 0 1 0 0
9 0 0 0 1
다음을 상기하자. \[
\begin{split}
\mathbf{y} &= [\mathbf{1}, X_a C_a, X_b C_b, (X_a C_a)\odot (X_b C_b)]
\begin{bmatrix} \mu \\ \boldsymbol{\alpha}^* \\ \boldsymbol{\beta}^* \\ \boldsymbol{\gamma}^* \end{bmatrix}
\\
&= \mathbf{1}_n\mu + X_aC_a\boldsymbol{\alpha}^* + X_bC_b\boldsymbol{\beta}^* + (X_aC_a \odot X_bC_b)\boldsymbol{\gamma}^*
\end{split}
\]
크로네커 곱 \(\otimes\) 를 이용, \((X_a C_a)\odot (X_b C_b) = (X_a \odot X_b)(C_b \otimes C_a)\) 로 쓸 수 있고, 따라서 \(\boldsymbol{\gamma} = (C_b \otimes C_a) \boldsymbol{\gamma}^*\) .
\(\boldsymbol{\gamma}\) , \(\boldsymbol{\gamma}^*\) 를 각각 \(r \times s\) 행렬 \(\Gamma\) , \((r-1)\times (s-1)\) 행렬 \(\Gamma^*\) 의 벡터 표현으로 이해하면 (\(\boldsymbol{\gamma}=\operatorname{vec}(\Gamma)\) ), \(\Gamma=C_a \Gamma^* C_{b}^T\) 가 성립한다. (왜 그런가?)
식별가능 제약은 \(C_a^T\) 의 널벡터 \(\mathbf{d}_a\) (즉 \(\mathbf{d}_a^T C_a = \mathbf{0}^T\) )와 \(C_b^T\) 의 널벡터 \(\mathbf{d}_b\) 에 대해 \[
\mathbf{d}_{a}^T \Gamma = \mathbf{0}^T,
\quad
\Gamma \mathbf{d}_{b} = \mathbf{0}
\] 으로 쓸 수 있다. (왜 그런가?)
무절편 2원배치 비가법모형 1
~ a * b - 1 # or ~ a + b + a: b - 1
대수적 기술 \[
\begin{split}
y_{ijk} &= \alpha_{i} + \beta_{j} + \gamma_{ij} + \epsilon_{ijk},
\\
& i=1,2,\dotsc,r;
~
j=1,2,\dotsc, s;
~
k=1,2,\dotsc, n_{ij}
\end{split}
\]
모형 행렬 : \(X^*=[X_a, X_b C_b, X_{a}^{(-r)} \odot (X_b C_b)]\)
\(X_{a}^{(-r)}\) 은 \(X_{a}\) 에서 마지막 열을 뺀 행렬
model.matrix (y ~ a * b - 1 , data = dat2)
a1 a2 a3 b2 b3 a2:b2 a3:b2 a2:b3 a3:b3
1 1 0 0 0 0 0 0 0 0
2 1 0 0 1 0 0 0 0 0
3 1 0 0 0 1 0 0 0 0
4 0 1 0 0 0 0 0 0 0
5 0 1 0 1 0 1 0 0 0
6 0 1 0 0 1 0 0 1 0
7 0 0 1 0 0 0 0 0 0
8 0 0 1 1 0 0 1 0 0
9 0 0 1 0 1 0 0 0 1
attr(,"assign")
[1] 1 1 1 2 2 3 3 3 3
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
무절편 2원배치 비가법모형 2
대수적 기술 \[
\begin{split}
y_{ijk} &= \alpha_{i} + \beta_{ij} + \epsilon_{ijk},
\\
& i=1,2,\dotsc,r;
~
j=1,2,\dotsc, s;
~
k=1,2,\dotsc, n_{ij}
\end{split}
\]
모형 행렬 : \(X^*=[X_a, X_{a} \odot (X_b C_b)]\)
model.matrix (y ~ - 1 + a + a: b, data = dat2)
a1 a2 a3 a1:b2 a2:b2 a3:b2 a1:b3 a2:b3 a3:b3
1 1 0 0 0 0 0 0 0 0
2 1 0 0 1 0 0 0 0 0
3 1 0 0 0 0 0 1 0 0
4 0 1 0 0 0 0 0 0 0
5 0 1 0 0 1 0 0 0 0
6 0 1 0 0 0 0 0 1 0
7 0 0 1 0 0 0 0 0 0
8 0 0 1 0 0 1 0 0 0
9 0 0 1 0 0 0 0 0 1
attr(,"assign")
[1] 1 1 1 2 2 2 2 2 2
attr(,"contrasts")
attr(,"contrasts")$a
[1] "contr.treatment"
attr(,"contrasts")$b
[1] "contr.treatment"
