제 10강: 선형모형 1

데이터과학 입문

Author

Affiliation

원중호

서울대학교 통계학과

시작하기 전에

다음의 패키지가 설치되어 있지 않으면 설치한다.

# install.packages("MASS")
# install.packages("tidyverse")
library(MASS)
library(tidyverse)
sessionInfo()

R version 4.4.0 (2024-04-24)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Asia/Seoul
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4    
 [5] purrr_1.0.2     readr_2.1.5     tidyr_1.3.1     tibble_3.2.1   
 [9] ggplot2_3.5.1   tidyverse_2.0.0 MASS_7.3-60.2  

loaded via a namespace (and not attached):
 [1] gtable_0.3.5      jsonlite_1.8.8    compiler_4.4.0    tidyselect_1.2.1 
 [5] scales_1.3.0      yaml_2.3.8        fastmap_1.1.1     R6_2.5.1         
 [9] generics_0.1.3    knitr_1.46        htmlwidgets_1.6.4 munsell_0.5.1    
[13] pillar_1.9.0      tzdb_0.4.0        rlang_1.1.3       utf8_1.2.4       
[17] stringi_1.8.3     xfun_0.43         timechange_0.3.0  cli_3.6.2        
[21] withr_3.0.0       magrittr_2.0.3    digest_0.6.35     grid_4.4.0       
[25] rstudioapi_0.16.0 hms_1.1.3         lifecycle_1.0.4   vctrs_0.6.5      
[29] evaluate_0.23     glue_1.7.0        fansi_1.0.6       colorspace_2.1-0 
[33] rmarkdown_2.26    tools_4.4.0       pkgconfig_2.0.3   htmltools_0.5.8.1

선형 모형 (Venables and Ripley Ch. 6)

• 고전 통계학의 핵심
  
• 많은 통계 실무의 기초
  
• 다수의 현대적 모형 및 분석 기법은 선형 모형을 위해 개발된 방법론을 기반으로 함
  

An Analysis of Covariance Example

Whiteside 자료

1960년대에 영국 건축 연구소(UK Building Research Station)의 데릭 화이트사이드(Derek Whiteside)가 수집한 자료 (Hand et al., 1994).

• 잉글랜드 남동부에 있는 화이트사이드의 집에서 캐비티 벽 단열재를 설치하기 전과 후의 두 번의 “난방 시즌” 동안 주간 가스 소비량과 평균 외부 기온을 기록
  
• 실험 목적: 단열재가 가스 소비에 미치는 영향을 평가
  

---

as_tibble(whiteside)

# 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을 사용할 수 있다.

gasB <- lm(Gas ~ Temp, data = MASS::whiteside, subset = Insul == "Before")
# gasB <- lm(Gas ~ 1+Temp, data = MASS::whiteside, subset = Insul == "Before") # same result as above
gasA <- 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 객체

분석 결과 요약

summary(gasB) # p. 142

## 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

summary(gasA)

## 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

---

표본 분산

(varB <- deviance(gasB) / gasB$df.resid)  # direct calculation 

[1] 0.07914867

(varB <- summary(gasB)$sigma^2)           # alternative 

[1] 0.07914867

---

• 두 회귀 모형을 하나의 lm 객체에 적합:

gasBA <- 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

• a가 요인인 a/x 형식의 모형공식 항은 “a의 각 수준에서 1 + x 꼴의 개별 회귀 모형”으로 생각할 수 있다.

• 추정치는 동일하지만 표준 오차가 달라진 것을 볼 수 있다. 이는 gasBA가 분산을 풀링해서 추정하기 때문이다.

---

이차식 회귀

gasQ <- 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

• ‘identity’ 함수 I()는 산술적 의미 그대로 인수를 평가

• 개별 회귀 계수를 보면 After 그룹에만 2차항이 필요할 가능성이 있지만 그 증거가 압도적이지 않음.

• 단순성을 이유로 개별 선형회귀 모형을 유지

---

분산분석

• 더 간단한 모형으로 기울기가 같은 회귀모형을 생각할 수 있다.
• 이 모형을 적합하고 anova 함수를 사용하여 수준별 기울기를 갖는 회귀 모형 내에서 검정할 수 있다. (왜 그런가?)

# MASS p. 143
gasPR <- 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

---

• 분산분석 결과는 수준별 기울기가 실제로 필요하다는 것을 보여준다.

• 수준별 기울기를 갖는 선형모형을 다른 방식으로 기술하면 분석분석을 보다 간단하고 유익하게 수행할 수 있다.

  gasBA1 <- lm(Gas ~ Insul * Temp, data = MASS::whiteside)
  summary(gasBA1)$coef

                    Estimate Std. Error    t value     Pr(>|t|)
  (Intercept)      6.8538277 0.13596397  50.409146 7.997414e-46
  InsulAfter      -2.1299780 0.18009172 -11.827185 2.315921e-16
  Temp            -0.3932388 0.02248703 -17.487358 1.976009e-23
  InsulAfter:Temp  0.1153039 0.03211212   3.590665 7.306852e-04

  • 함수 options와 모형 공식인 formula가 회귀모형의 모수에 영향을 미침을 알 수 있다.

모형 공식과 모형 행렬

모형 공식

선형 모형은 응답 벡터 \(\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} \]

• 모형행렬 \(X\) (\(X_a\) 아님)의 랭크는 얼마인가?

• 이 모형은 식별가능한가?

---

식별가능 모형 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)
  dat <- data.frame(a = factor(rep(1:3, each=3)), 
                    y = rnorm(9, rep(2:4, each=3), 0.1))
  mod <- lm(y ~ a - 1, data = dat)
  coef(mod)

        a1       a2       a3 
  1.986017 2.893121 4.082296 

---

    model.matrix(mod)

  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"

---

식별가능 모형 2: 부호화 행렬

• 알려진 사실 \(X_a\mathbf{1}_k = \mathbf{1}_n\)로부터 (왜 그런가?) \[ \begin{split} [\mathbf{1}_n~X_a]\begin{bmatrix} \mu \\ \boldsymbol{\alpha} \end{bmatrix} &= \mathbf{1}_n\mu + X_a\boldsymbol{\alpha} \\ &= X_a\mathbf{1}_k\mu + X_a\boldsymbol{\alpha} = [\mathbf{1}_n~X_a]\begin{bmatrix} 0 \\ \mathbf{1}_k\mu + \boldsymbol{\alpha} \end{bmatrix} \end{split} \] 임을 알 수 있다.

• \(\mu\)의 값이 임의임도 알 수 있다.

---

• 차원 \(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이 다른 모든 수준을 비교하는 표준 수준으로 취급된다.

mod_treatment <- 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"

---

(alph.star <- coef(mod_treatment)) # mu^* = (Intercept)

(Intercept)          a2          a3 
  1.9860166   0.9071045   2.0962792 

(Ca <- contr.treatment(length(levels(dat$a))))

  2 3
1 0 0
2 1 0
3 0 1

(alph <- Ca %*% alph.star[-1])    # letting $mu = mu^*$

       [,1]
1 0.0000000
2 0.9071045
3 2.0962792

---

헬메르트 부호화 (Helmert coding)

S-plus의 기본 부호화. 각 열의 합이 0이 되도록 한다.

mod_helmert <- lm(y ~ a, data=dat, contrasts = list(a = "contr.helmert"))
(alph.star <- coef(mod_helmert))

(Intercept)          a1          a2 
  2.9871445   0.4535523   0.5475757 

(Ca <- contr.helmert(length(levels(dat$a))))

  [,1] [,2]
1   -1   -1
2    1   -1
3    0    2

(alph <- Ca %*% alph.star[-1])   # letting $mu = mu^*$

         [,1]
1 -1.00112792
2 -0.09402339
3  1.09515132

sum(alph)   # numerically zero. why?

[1] 0

---

직교다항식 부호화 (orthogonal polynomial coding)

R의 순서형 변수에 대한 기본 부호화.

dat$b <- as.ordered(dat$a)
mod_poly <- lm(y ~ b, data=dat)
(alph.star <- coef(mod_poly))

(Intercept)         b.L         b.Q 
  2.9871445   1.4822933   0.1151547 

(Ca <- 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

(alph <- Ca %*% alph.star[-1])     # letting $mu = mu^*$

            [,1]
[1,] -1.00112792
[2,] -0.09402339
[3,]  1.09515132

sum(alph)   # numerically zero

[1] 0

모수 해석

• 원모수 \(\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\)으로 설정)

Ca <- 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

D <- 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\).

Note

왜 앞 식별가능모형1과 식별가능모형2–처리부호화 결과가 같았는가?

---

예: 헬메르트 부호화

Ca <- 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

D <- solve(cbind(1, Ca)) 
D %>% 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} \]

---

Note

왜 앞 식별가능모형2–헬메르트부호화 결과에서 모수의 합이 (수치적으로) 0이었는가?

• 헬메르트 부호화행렬의 모든 열은 직교한다. 따라서 복호화행렬의 모든 행은 직교한다. (왜 그런가?)

  t(Ca) %*% Ca

       [,1] [,2] [,3] [,4]
  [1,]    2    0    0    0
  [2,]    0    6    0    0
  [3,]    0    0   12    0
  [4,]    0    0    0   20

---

직교다항식 부호화

• \(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\)의 열은 수직선에서 등간격으로 설정.

---

(Ca <- 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 

D <- solve(cbind(1, Ca)) 
D %>% 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\).

Note

왜 앞 식별가능모형2–직교부호화 결과에서 모수의 합이 (수치적으로) 0이었는가?

---

2원 배치

• a, b가 각각 \(r\), \(s\)개 수준을 갖는 범주형 변수(factor)라 하자.

  dat2 <- dat  # copy
  dat2$b <- factor(rep(1:3, each=3))

2원 배치 가법 모형

Note

    y ~ 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  0  0
  3           1  0  0  0  0
  4           1  1  0  1  0
  5           1  1  0  1  0
  6           1  1  0  1  0
  7           1  0  1  0  1
  8           1  0  1  0  1
  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  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")$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  0  0
  3  1  0  0  0  0
  4  0  1  0  1  0
  5  0  1  0  1  0
  6  0  1  0  1  0
  7  0  0  1  0  1
  8  0  0  1  0  1
  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  1  0  0  0  0
3  1  0  0  0  0
4  0  1  0  1  0
5  0  1  0  1  0
6  0  1  0  1  0
7  0  0  1  0  1
8  0  0  1  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원 배치 비가법모형

Note

    y ~ a + b + a:b
    또는
    y ~ a * b

• 이 모형은 \(rs\) 수준을 갖는 1원배치 모형과 사실상 같다.

• 대수적 기술 \[ y_{ijk} = \mu + \alpha_{i} + \beta_{j} + \gamma_{ij} + \epsilon_{ijk}, \quad i=1,2,\dotsc,r; ~ j=1,2,\dotsc, s; ~ k=1,2,\dotsc, n_{ij} \]

---

• 모형 행렬 : \(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  0  0     0     0     0     0
  3           1  0  0  0  0     0     0     0     0
  4           1  1  0  1  0     1     0     0     0
  5           1  1  0  1  0     1     0     0     0
  6           1  1  0  1  0     1     0     0     0
  7           1  0  1  0  1     0     0     0     1
  8           1  0  1  0  1     0     0     0     1
  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"

---

• 다음과 비교하자.

  X_star <- numeric()
  XaCa <- model.matrix(y ~ a, data=dat2)[, -1] 
  XbCb <- model.matrix(y ~ b, data=dat2)[, -1]
  for (i in seq_len(ncol(XbCb))) 
      for (j in seq_len(ncol(XaCa)))
          X_star <- cbind(X_star, XbCb[, i] * XaCa[, j]) 
  X_star

    [,1] [,2] [,3] [,4]
  1    0    0    0    0
  2    0    0    0    0
  3    0    0    0    0
  4    1    0    0    0
  5    1    0    0    0
  6    1    0    0    0
  7    0    0    0    1
  8    0    0    0    1
  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

Note

    y ~ a * b - 1
    또는
    y ~ a + b + a:b - 1

• 대수적 기술 \[ y_{ijk} = \alpha_{i} + \beta_{j} + \gamma_{ij} + \epsilon_{ijk}, \quad i=1,2,\dotsc,r; ~ j=1,2,\dotsc, s; ~ k=1,2,\dotsc, n_{ij} \]

• 모형 행렬 : \(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  0  0     0     0     0     0
3  1  0  0  0  0     0     0     0     0
4  0  1  0  1  0     1     0     0     0
5  0  1  0  1  0     1     0     0     0
6  0  1  0  1  0     1     0     0     0
7  0  0  1  0  1     0     0     0     1
8  0  0  1  0  1     0     0     0     1
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

Note

    y ~ - 1 + a + a:b

• 대수적 기술 \[ y_{ijk} = \alpha_{i} + \beta_{ij} + \epsilon_{ijk}, \quad i=1,2,\dotsc,r; ~ j=1,2,\dotsc, s; ~ k=1,2,\dotsc, n_{ij} \]

• 모형 행렬 : \(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     0     0     0     0     0     0
3  1  0  0     0     0     0     0     0     0
4  0  1  0     0     1     0     0     0     0
5  0  1  0     0     1     0     0     0     0
6  0  1  0     0     1     0     0     0     0
7  0  0  1     0     0     0     0     0     1
8  0  0  1     0     0     0     0     0     1
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"