1 개요
반복측정 MANOVA는 같은 피험자를 여러 시점에 반복 측정한 데이터를 다변량 반응 벡터로 취급하는 분석법이다. 각 피험자의 \(n\) 시점 관측값을 \(n \times 1\) 벡터로 보고, 이 벡터들의 공분산 구조를 자유롭게 추정한다.
이 포스트는 일표본(one-sample) 설계, 즉 집단 구분 없이 피험자들의 시간에 따른 변화를 분석하는 기본 구조를 다룬다. 다표본 설계는 별도 포스트에서 다룬다.
1.1 학습 지도
아래 질문에 답할 수 있으면 이 포스트를 이해한 것이다.
- 왜 MANOVA는 반복측정 ANOVA보다 더 유연한 공분산 가정을 허용하는가?
- 직교 다항식 행렬 \(\mathbf{P}\) 를 어떻게 구성하고, 왜 사용하는가?
- \(\mathbf{SST}^*\) 와 \(\mathbf{SSR}^*\) 의 대각원소가 각각 무엇을 의미하는가?
- \(\mathbf{SST}^*\) 와 \(\mathbf{SSR}^*\) 에서 반복측정 ANOVA 결과를 어떻게 추출하는가?
- 다변량 시간 효과 검정에서 행렬식 방정식과 Cholesky 분해의 역할은?
- 구형성 충족 시와 위반 시 개별 추세 F-검정이 어떻게 달라지고, 검정력은 왜 다른가?
2 일표본 MANOVA 모형
\(N\) 명의 피험자( \(i = 1, \ldots, N\) )가 모두 \(n\) 시점( \(j = 1, \ldots, n\) )에 완전 측정된다.
\[ \mathbf{y}_i = \boldsymbol{\mu} + \boldsymbol{\varepsilon}_i, \quad i = 1, \ldots, N \tag{3.1} \]
- \(\mathbf{y}_i = (y_{i1}, y_{i2}, \ldots, y_{in})' \in \mathbb{R}^n\): 피험자 \(i\) 의 반복 측정 벡터
- \(\boldsymbol{\mu} = (\mu_1, \mu_2, \ldots, \mu_n)' \in \mathbb{R}^n\): 시점별 모집단 평균 벡터
- \(\boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma})\): 오차 벡터, \(\boldsymbol{\Sigma}\) 는 \(n \times n\) 양정치 대칭 행렬
이 모형의 핵심은 \(\boldsymbol{\Sigma}\) 에 아무 제약이 없다는 점이다.
반복측정 ANOVA는 \(\boldsymbol{\Sigma} = \sigma_\pi^2 \mathbf{1}\mathbf{1}' + \sigma_e^2 \mathbf{I}\) 로 제약했지만, MANOVA에서 \(\boldsymbol{\Sigma}\) 의 \((j, k)\) 원소 \(\sigma_{jk}\) 는 시점 \(j\) 와 \(k\) 사이의 공분산을 아무 구조 없이 추정한다.
직관: 우울증 임상시험에서 “1주차-2주차 상관”과 “1주차-8주차 상관”이 같다고 가정하는 것이 반복측정 ANOVA다. MANOVA는 이 가정을 버리고, 데이터가 알아서 결정하게 한다. 일반적으로 인접한 시점일수록 더 강하게 상관되는데, MANOVA는 이 현실적인 패턴을 수용한다.
2.1 반복측정 ANOVA와의 연결
반복측정 ANOVA의 공분산 가정( \(\boldsymbol{\Sigma} = \sigma_\pi^2 \mathbf{1}\mathbf{1}' + \sigma_e^2 \mathbf{I}\) )을 MANOVA에 대입하면, 두 모형이 같아진다. 따라서 반복측정 ANOVA는 MANOVA의 특수 사례다.
이 관계 덕분에, 아래에서 볼 것처럼, MANOVA의 SSCP 행렬에서 반복측정 ANOVA의 제곱합(SS)을 직접 추출할 수 있다.
3 직교 다항식 성장곡선 모형
\(n\) 시점 평균 벡터 \(\boldsymbol{\mu} = (\mu_1, \ldots, \mu_n)'\) 를 시간의 함수로 모수화하는 방법이 필요하다. 가장 유연하고 해석하기 쉬운 방법이 직교 다항식 성장곡선 분석이다.
3.1 1단계: 시간 행렬 \(\mathbf{T}\) 구성
시점 값 \(t_1, t_2, \ldots, t_n\) 과 다항식 차수 \(q \leq n\) 에 대해 \(q \times n\) 시간 행렬 \(\mathbf{T}\) 를 다음과 같이 정의한다.
\[ \mathbf{T} = \begin{bmatrix} 1 & 1 & \cdots & 1 \\ t_1 & t_2 & \cdots & t_n \\ t_1^2 & t_2^2 & \cdots & t_n^2 \\ \vdots & & & \vdots \\ t_1^{q-1} & t_2^{q-1} & \cdots & t_n^{q-1} \end{bmatrix} \]
\(\mathbf{T}'\) 는 \(n \times q\) 행렬이고, 평균 벡터의 다항 모형은
\[ \boldsymbol{\mu} = \mathbf{T}' \boldsymbol{\beta} \]
\(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_{q-1})'\) 는 다항 계수 벡터다.
4시점 예시 ( \(t = 1, 2, 3, 4\), \(q = 4\) ):
\[ \mathbf{T} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 3 & 4 \\ 1 & 4 & 9 & 16 \\ 1 & 8 & 27 & 64 \end{bmatrix} \]
3.2 2단계: 직교화 — Cholesky 분해로 \(\mathbf{P}\) 유도
\(\mathbf{T}\) 를 그대로 사용하면 다공선성(collinearity) 문제가 생긴다. \(t_j\) 와 \(t_j^2\) 가 강하게 상관되므로, 회귀 계수 추정이 불안정해진다.
이를 해결하기 위해 \(\mathbf{T}\) 를 직교 다항식 행렬 \(\mathbf{P}\) 로 변환한다. 변환 방법은 \(\mathbf{T}\mathbf{T}'\) 의 Cholesky 분해다:
\[ \mathbf{T}\mathbf{T}' = \mathbf{S}\mathbf{S}', \quad \mathbf{P} = \mathbf{S}^{-1}\mathbf{T} \tag{Cholesky} \]
\(\mathbf{S}\) 는 \(q \times q\) 하삼각(lower triangular) 행렬이다. 이때 \(\mathbf{P}\) 는 다음 성질을 만족한다:
\[ \mathbf{P}\mathbf{P}' = \mathbf{S}^{-1}\mathbf{T}\mathbf{T}'(\mathbf{S}^{-1})' = \mathbf{S}^{-1}\mathbf{S}\mathbf{S}'(\mathbf{S}')^{-1} = \mathbf{I} \]
즉, \(\mathbf{P}\) 의 각 행은 단위 벡터이고 서로 직교한다.
직관: Cholesky 분해는 “\(\mathbf{T}\) 의 행들이 서로 겹치는 방향을 제거하는 회전”이다. Gram-Schmidt 직교화와 수치적으로 동치다. 직교화 후 각 행이 “순수하게 상수, 선형, 2차, 3차 성분만”을 담게 된다.
3.3 3단계: 4시점의 직교 다항식 행렬 \(\mathbf{P}_4\)
균등 간격 4시점( \(n = q = 4\) )의 경우, 직교 다항식은 통계학 교재(Pearson & Hartley, 1976)에 표로 수록되어 있다:
\[ \mathbf{P}_4 = \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\[6pt] \frac{-3}{\sqrt{20}} & \frac{-1}{\sqrt{20}} & \frac{1}{\sqrt{20}} & \frac{3}{\sqrt{20}} \\[6pt] \frac{1}{2} & \frac{-1}{2} & \frac{-1}{2} & \frac{1}{2} \\[6pt] \frac{-1}{\sqrt{20}} & \frac{3}{\sqrt{20}} & \frac{-3}{\sqrt{20}} & \frac{1}{\sqrt{20}} \end{bmatrix} = \begin{bmatrix} \mathbf{p}_0' \\ \mathbf{p}_1' \\ \mathbf{p}_2' \\ \mathbf{p}_3' \end{bmatrix} \tag{3.2} \]
각 행의 의미:
| 행 | 이름 | 패턴 | 검정 의미 |
|---|---|---|---|
| \(\mathbf{p}_0\) | 상수(constant) | 모두 \(+\) | 전체 평균 수준 |
| \(\mathbf{p}_1\) | 선형(linear) | \(-\)에서 \(+\) 로 증가 | 시간에 따른 선형 추세 |
| \(\mathbf{p}_2\) | 2차(quadratic) | \(+, -, -, +\) 패턴 | 가속·감속, U자·역U자형 |
| \(\mathbf{p}_3\) | 3차(cubic) | \(-, +, -, +\) 패턴 | S자 곡선형 변곡 |
직관: \(\mathbf{p}_1\) 계수 \((-3, -1, 1, 3)\) 은 “1시점이 제일 낮고 4시점이 제일 높은가?”를 묻는다. \(\mathbf{p}_2\) 계수 \((1, -1, -1, 1)\) 은 “중간 시점들이 끝 시점들보다 낮은가?”를 묻는다. 각 계수 벡터가 서로 직교하므로, 선형 추세 검정이 2차 추세 추정에 영향을 주지 않는다.
3.4 4단계: 직교 다항 변환된 모형
\(\mathbf{P}\) 로 모형 양변을 곱하면
\[ \mathbf{P}\mathbf{y}_i = \mathbf{P}\boldsymbol{\mu} + \mathbf{P}\boldsymbol{\varepsilon}_i \]
\[ \underbrace{\mathbf{P}\mathbf{y}_i}_{\text{변환된 반응}} = \underbrace{\boldsymbol{\theta}}_{\text{다항 계수 벡터}} + \underbrace{\boldsymbol{\varepsilon}_i^*}_{\text{변환된 오차}} \tag{3.3} \]
여기서
\[ \boldsymbol{\theta} = \mathbf{P}\boldsymbol{\mu}, \quad \boldsymbol{\varepsilon}_i^* = \mathbf{P}\boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma}^* = \mathbf{P}\boldsymbol{\Sigma}\mathbf{P}') \]
\(\hat{\boldsymbol{\theta}} = \mathbf{P}\bar{\mathbf{y}}\) 는 변환된 평균 벡터의 최소제곱 추정량이다. (\(\bar{\mathbf{y}}\) 는 \(n \times 1\) 시점별 표본 평균 벡터)
직관: 변환 전 모형은 “시점 \(j\) 의 평균이 얼마인가?”를 묻고, 변환 후 모형은 “선형 추세 계수가 얼마인가?”, “2차 추세 계수가 얼마인가?”를 묻는다. 직교 변환은 시점별 평균의 언어를 추세 성분의 언어로 바꾸는 것이다.
3.5 구형성 조건의 재해석
반복측정 ANOVA의 구형성(sphericity) 가정은 MANOVA의 언어로 표현하면:
\[ \boldsymbol{\Sigma}^*_{(n-1)} = \mathbf{P}_{(n-1)}\boldsymbol{\Sigma}\mathbf{P}_{(n-1)}' = \sigma_e^2 \mathbf{I}_{n-1} \]
즉, \(\boldsymbol{\Sigma}^*\) 의 하 \((n-1) \times (n-1)\) 부분행렬이 대각원소가 모두 \(\sigma_e^2\) 로 같고 비대각원소가 모두 0인 행렬일 조건이다.
이 조건이 성립하면 MANOVA에서 반복측정 ANOVA 결과를 추출하는 것이 정당화된다. MANOVA는 이 조건 없이도 타당한 검정을 제공한다.
4 SSCP 행렬 분해
직교 다항 변환된 모형에서 데이터의 변동을 SSCP(Sum of Squares and Cross-Products) 행렬로 분해한다.
| 변동원 | df | SSCP ( \(n \times n\) ) | \(E(\text{SSCP})\) |
|---|---|---|---|
| 시간(Time) | 1 | \(\mathbf{SST}^* = N\mathbf{P}\bar{\mathbf{y}}\bar{\mathbf{y}}'\mathbf{P}'\) | \(\mathbf{P}[\boldsymbol{\Sigma} + N\boldsymbol{\mu}\boldsymbol{\mu}']\mathbf{P}'\) |
| 잔차(Residual) | \(N-1\) | \(\mathbf{SSR}^* = \mathbf{P}(\mathbf{Y}'\mathbf{Y} - N\bar{\mathbf{y}}\bar{\mathbf{y}}')\mathbf{P}'\) | \((N-1)\mathbf{P}\boldsymbol{\Sigma}\mathbf{P}'\) |
| 전체(Total) | \(N\) | \(\mathbf{SSY}^* = \mathbf{P}\mathbf{Y}'\mathbf{Y}\mathbf{P}'\) |
여기서 \(\mathbf{Y}\) 는 \(N \times n\) 전체 데이터 행렬(피험자 × 시점), \(\bar{\mathbf{y}}\) 는 \(n \times 1\) 시점별 표본 평균 벡터다.
각 df는 피험자-간(between-subjects) 정보의 양을 반영한다. 시간 SSCP의 df가 1인 이유는 전체 평균 벡터가 1개의 피험자-간 정보 단위를 사용하기 때문이다.
4.1 SST* 행렬의 구조
\[ \mathbf{SST}^* = N\mathbf{P}\bar{\mathbf{y}}\bar{\mathbf{y}}'\mathbf{P}' = N(\mathbf{P}\bar{\mathbf{y}})(\mathbf{P}\bar{\mathbf{y}})' = N\hat{\boldsymbol{\theta}}\hat{\boldsymbol{\theta}}' \]
\(\mathbf{SST}^*\) 는 랭크-1 행렬이다. 즉, \(\hat{\boldsymbol{\theta}}\) 방향의 벡터 외적(outer product)이다.
대각원소의 의미 (3.3):
\[ \mathbf{SST}^* = \text{diag}\begin{pmatrix} \underbrace{N(\mathbf{p}_0'\bar{\mathbf{y}})^2}_{sst_0 \text{: 상수}} ,\; \underbrace{N(\mathbf{p}_1'\bar{\mathbf{y}})^2}_{sst_1 \text{: 선형}},\; \underbrace{N(\mathbf{p}_2'\bar{\mathbf{y}})^2}_{sst_2 \text{: 2차}},\; \ldots,\; \underbrace{N(\mathbf{p}_{n-1}'\bar{\mathbf{y}})^2}_{sst_{n-1} \text{: (n-1)차}} \end{pmatrix} \]
- \(sst_0 = N \cdot n\bar{y}^2_{..}\): 전체 평균의 제곱합 (상수항, 시간 효과와 무관)
- \(sst_k = N(\mathbf{p}_k'\bar{\mathbf{y}})^2 (k = 1, \ldots, n-1)\): \(k\) 차 시간 추세의 제곱합
직관: \(sst_1\) 은 “\(N\) 명의 평균 궤적이 선형으로 얼마나 기울어져 있는가”를 제곱으로 표현한 것이다. \(\bar{\mathbf{y}}\) 를 선형 방향 \(\mathbf{p}_1\) 에 투영한 거리의 제곱에 \(N\) 을 곱한 것이다.
4.2 SSR* 행렬의 구조
\[ \mathbf{SSR}^* = \mathbf{P}(\mathbf{Y}'\mathbf{Y} - N\bar{\mathbf{y}}\bar{\mathbf{y}}')\mathbf{P}' = \sum_{i=1}^N \mathbf{P}(\mathbf{y}_i - \bar{\mathbf{y}})(\mathbf{y}_i - \bar{\mathbf{y}})'\mathbf{P}' \]
마지막 등호: \(\mathbf{Y}'\mathbf{Y} - N\bar{\mathbf{y}}\bar{\mathbf{y}}' = \sum_i(\mathbf{y}_i - \bar{\mathbf{y}})(\mathbf{y}_i - \bar{\mathbf{y}})'\) (분산-공분산의 비편향 분자)
대각원소의 의미:
\[ \mathbf{SSR}^* = \text{diag}\begin{pmatrix} \underbrace{ssr_0}_{\text{피험자 SS}},\; \underbrace{ssr_1}_{\text{피험자}\times\text{선형}},\; \underbrace{ssr_2}_{\text{피험자}\times\text{2차}},\; \ldots,\; \underbrace{ssr_{n-1}}_{\text{피험자}\times(n-1)\text{차}} \end{pmatrix} \]
- \(ssr_0 = n\sum_i(\bar{y}_{i.} - \bar{y}_{..})^2\): 피험자 간 분산 (Subjects SS)
- \(ssr_k = \sum_i (\mathbf{p}_k'\mathbf{y}_i - \mathbf{p}_k'\bar{\mathbf{y}})^2\) (\(k = 1, \ldots, n-1\)): 피험자 내 \(k\) 차 추세의 개인 간 변동
직관: \(ssr_1\) 은 “각 피험자의 선형 기울기가 집단 평균 기울기에서 얼마나 벗어나는가”의 총합이다. \(ssr_0\) 은 “피험자들의 전반적 수준(절편)이 얼마나 다른가”를 나타낸다. 구형성이 성립하면 \(ssr_1, ssr_2, \ldots, ssr_{n-1}\) 이 모두 같은 기댓값을 가진다.
5 반복측정 ANOVA 결과 추출
\(\mathbf{SST}^*\) 와 \(\mathbf{SSR}^*\) 는 반복측정 ANOVA에 필요한 모든 SS를 담고 있다.
구형성이 성립할 때 추출 공식:
\[ SS_S = ssr_0 = \mathbf{SSR}^*[1,1] \quad (\text{피험자 SS}) \]
\[ SS_T = \sum_{k=1}^{n-1} sst_k = \text{하 } (n-1) \text{개 대각원소의 합} \quad (\text{시간 SS}) \]
\[ SS_R = \sum_{k=1}^{n-1} ssr_k = \text{하 } (n-1) \text{개 대각원소의 합} \quad (\text{잔차 SS}) \]
이로부터 반복측정 ANOVA 분산분해표:
| 변동원 | df | SS | MS | \(E(\text{MS})\) | F |
|---|---|---|---|---|---|
| 피험자 | \(N-1\) | \(SS_S = n\sum_i(\bar{y}_{i.}-\bar{y}_{..})^2\) | \(MS_S\) | \(\sigma_e^2 + n\sigma_\pi^2\) | \(MS_S/MS_R\) |
| 시간 | \(n-1\) | \(SS_T = N\sum_j(\bar{y}_{.j}-\bar{y}_{..})^2\) | \(MS_T\) | \(\sigma_e^2 + \frac{N}{n-1}\sum_j\tau_j^2\) | \(MS_T/MS_R\) |
| 잔차 | \((N-1)(n-1)\) | \(SS_R\) | \(MS_R\) | \(\sigma_e^2\) |
왜 “잔차”의 df가 \((N-1)(n-1)\) 인가?
\[ MS_R = \frac{SS_R}{(N-1)(n-1)} = \frac{\sum_{k=1}^{n-1} ssr_k}{(N-1)(n-1)} = \frac{(n-1) \cdot \overline{ssr}}{(N-1)(n-1)} = \frac{\overline{ssr}}{N-1} \]
\(\overline{ssr}\) 은 \(ssr_1, \ldots, ssr_{n-1}\) 의 평균이다. 따라서 \(MS_R\) 은 “\(n-1\) 개의 피험자-내 오차 분산 추정값의 평균” 이다. 구형성이 성립하면 이들의 기댓값이 같으므로, 평균이 정당화된다.
6 다변량 시간 효과 검정
6.1 검정 논리: 행렬의 비교
귀무가설 \(H_0: \boldsymbol{\mu}\) 의 모든 원소가 같다 (시간 효과 없음, \(\boldsymbol{\tau} = \mathbf{0}\)).
일변량 F-검정의 논리: \[ F = \frac{MS_T}{MS_R} \approx 1 \; (H_0 \text{하}) \quad \text{vs} \quad \frac{MS_T}{MS_R} \gg 1 \; (H_1 \text{하}) \]
다변량에서는 분수를 행렬로 일반화해야 한다. 두 \(n \times n\) 행렬의 “비”를 스칼라로 요약하는 방법이 행렬식 방정식(determinantal equation) 이다.
6.2 행렬식 방정식 (Determinantal Equation)
\(\mathbf{SST}^*\) 와 \(\mathbf{SSR}^*\) 의 하 \((n-1) \times (n-1)\) 부분행렬을 추출한다. (상수항 성분 \(sst_0, ssr_0\) 은 시간 효과 검정과 무관하므로 제외)
\[ |\mathbf{SST}_{(n-1)}^* - \lambda\, \mathbf{SSR}_{(n-1)}^*| = 0 \tag{3.4} \]
이 방정식을 \(\lambda\) 에 대해 풀면 고유값(eigenvalue) \(\lambda_1\) 을 얻는다.
일표본에서는 비영(non-zero) 고유값이 하나뿐이다. 이유: \(\mathbf{SST}_{(n-1)}^*\) 가 랭크-1 행렬이기 때문이다. (\(\mathbf{SST}^* = N\hat{\boldsymbol{\theta}}\hat{\boldsymbol{\theta}}'\) — 벡터의 외적은 랭크 1)
고유값의 직관:
- \(H_0\) 가 참이면: \(\mathbf{SST}_{(n-1)}^* \approx \mathbf{SSR}_{(n-1)}^*\) → \(|\mathbf{SST}^* - 1 \cdot \mathbf{SSR}^*| \approx 0\) → \(\lambda_1 \approx 1\)
- 시간 효과가 클수록: \(\mathbf{SST}_{(n-1)}^* \gg \mathbf{SSR}_{(n-1)}^*\) → \(\lambda_1 \gg 1\)
이것이 일변량 F-값과 다변량 고유값의 유사점이다: 둘 다 “신호/잡음”의 비를 나타낸다.
6.3 Cholesky 분해로 단순화
식 (3.4)는 두 행렬 고유값 문제(two-matrix eigenproblem) 다. 이를 더 계산하기 쉬운 한 행렬 고유값 문제로 변환한다.
\(\mathbf{SSR}_{(n-1)}^* = \mathbf{E}\mathbf{E}'\) (Cholesky 분해, \(\mathbf{E}\) 는 하삼각 행렬)로 놓으면
\[ |\mathbf{E}^{-1}\mathbf{SST}_{(n-1)}^*(\mathbf{E}^{-1})' - \lambda\mathbf{I}_{n-1}| = 0 \tag{3.5} \]
유도:
\[ |\mathbf{SST}^* - \lambda\mathbf{SSR}^*| = |\mathbf{SST}^* - \lambda\mathbf{E}\mathbf{E}'| = 0 \]
\(\mathbf{E}^{-1}\) 을 왼쪽에서, \((\mathbf{E}')^{-1} = (\mathbf{E}^{-1})'\) 을 오른쪽에서 곱하면
\[ |\mathbf{E}^{-1}\mathbf{SST}^*(\mathbf{E}^{-1})' - \lambda\mathbf{I}| = 0 \]
이제 표준 고유값 문제로, 대부분의 선형대수 라이브러리로 직접 계산할 수 있다.
직관: Cholesky 분해는 “\(\mathbf{SSR}^*\) 의”비율 척도”를 제거”하는 작업이다. \(\mathbf{E}^{-1}\) 는 잔차의 자연적 단위를 정규화하는 역할을 한다. 정규화 후, \(\mathbf{E}^{-1}\mathbf{SST}^*(\mathbf{E}^{-1})'\) 의 고유값이 잔차 대비 신호의 크기를 나타낸다.
6.4 다변량 검정 통계량
고유값 \(\lambda_1\) 으로부터 네 가지 검정 통계량을 계산한다. 일표본에서는 비영 고유값이 \(\lambda_1\) 하나뿐이므로, 네 통계량이 모두 같은 결론을 준다.
| 통계량 | 수식 | 해석 |
|---|---|---|
| Roy 최대근 | \(\theta_1 = \lambda_1\) | 주방향의 효과 크기 |
| Wilks’ \(\Lambda\) | \(\Lambda = \dfrac{1}{1 + \lambda_1}\) | 전체 분산 중 잔차 비율 |
| Hotelling-Lawley | \(T^2 = \lambda_1\) | 마할라노비스 거리의 일반화 |
| Pillai-Bartlett | \(V = \dfrac{\lambda_1}{1+\lambda_1} = 1 - \Lambda\) | 설명된 분산 비율 |
이 통계량들은 귀무가설 하에서 F-분포로 근사된다. 일표본에서 Wilks’ \(\Lambda\) 의 F-변환:
\[ F = \frac{1-\Lambda}{\Lambda} \cdot \frac{N-n+1}{n-1} \approx F(n-1, N-n+1) \quad \text{정확한 변환} \]
( \(n-1\) : 시간 자유도, \(N-n+1\) : 잔차 자유도)
각 통계량의 특성:
- Wilks’ \(\Lambda\): \(\Lambda = |\mathbf{SSR}^*| / |\mathbf{SST}^* + \mathbf{SSR}^*|\) 로도 쓸 수 있다. “오차 분산이 전체 분산에서 차지하는 비율”이므로, \(\Lambda = 1\) 은 효과 없음, \(\Lambda = 0\) 은 완벽한 분리.
- Roy 최대근: 단일 효과 방향이 하나일 때 가장 강력하나, 분포 근사가 어려움.
- Pillai-Bartlett: 귀무가설 위반 형태에 가장 견고(robust). 소표본·분포 이탈 시 선호.
- Hotelling-Lawley: 이변량 Hotelling \(T^2\) 검정의 일반화.
7 개별 시간 추세 성분 검정
전체 시간 효과 검정 후, 어떤 차수의 추세가 유의한지 확인한다. 구형성 충족 여부에 따라 두 가지 F-검정 방식이 존재한다.
7.1 방법 A: 구형성 충족 시 — 공통(averaged) 분모
구형성이 성립하면 \(ssr_1 = ssr_2 = \cdots = ssr_{n-1}\) (기댓값 동일). 따라서 이를 평균한 \(MS_R\) 을 공통 분모로 사용할 수 있다.
\[ MS_R = \frac{\sum_{k=1}^{n-1} ssr_k}{(n-1)(N-1)} \]
각 추세 성분의 F-검정 (3.6):
\[ F_k = \frac{sst_k}{\displaystyle\frac{\sum_{j=1}^{n-1} ssr_j}{(n-1)(N-1)}} = \frac{sst_k}{MS_R}, \quad k = 1, 2, \ldots, n-1 \tag{3.6} \]
자유도: \(df_\text{분자} = 1\), \(df_\text{분모} = (n-1)(N-1)\)
예: \(N = 64\), \(n = 4\) 이면 \(df_\text{분모} = 3 \times 63 = 189\).
7.2 방법 B: 구형성 위반 시 — 개별(separate) 분모
구형성이 성립하지 않으면 \(ssr_k\) 들의 기댓값이 서로 다를 수 있다. 이 경우 각 추세 성분에 고유한 오차 항 \(ssr_k\) 을 사용한다:
\[ F_k = \frac{sst_k}{ssr_k/(N-1)}, \quad k = 1, 2, \ldots, n-1 \tag{3.7} \]
자유도: \(df_\text{분자} = 1\), \(df_\text{분모} = N-1\)
예: \(N = 64\), \(n = 4\) 이면 \(df_\text{분모} = 63\).
7.3 검정력(Power) 비교
방법 A와 B의 핵심 차이는 분모 자유도다.
| 방법 | 분모 df | 기각역 ( \(\alpha = 0.05\) ) |
|---|---|---|
| A (구형성 충족, RM ANOVA) | \((n-1)(N-1)\) | \(F_{0.05}(1, 189) \approx 3.89\) |
| B (구형성 위반, MANOVA) | \(N-1\) | \(F_{0.05}(1, 63) \approx 3.99\) |
분모 자유도가 클수록 기각역이 작아져 검정력이 높아진다. 방법 A는 분모 df가 \((n-1)\) 배 더 크므로, 구형성이 성립하면 방법 B보다 더 강력하다.
직관: 방법 A는 \(n-1\) 개의 오차 추정값을 풀링(pooling) 하여 더 정밀한 오차 추정치를 만든다. “평균을 많이 낼수록 추정이 안정적”이라는 중심극한정리의 논리다. 구형성이 위반되면 이 풀링 자체가 부적절하므로, 방법 B처럼 각자 분모를 써야 한다.
- Mauchly 검정으로 구형성 평가
- 구형성 충족 ( \(p > 0.05\) ): 방법 A 사용 (더 강력)
- 구형성 위반 ( \(p \leq 0.05\) ): 방법 B(MANOVA) 또는 RM ANOVA + GG/HF 보정
- 소표본( \(N < 40\) )에서는 Mauchly 검정의 검정력이 낮으므로 보수적으로 방법 B 선택 권장
8 수치 예제: 어휘 성장 데이터
Bock(1975)의 어휘 성장 데이터를 이용하여 위 절차를 수치로 추적한다. \(N = 64\) 명, \(n = 4\) 시점(4~7학년), 어휘 점수.
8.1 기술통계
| 학년 | 평균 \(\bar{y}_{.j}\) |
|---|---|
| 4학년 (\(j=1\)) | 1.14 |
| 5학년 (\(j=2\)) | 2.54 |
| 6학년 (\(j=3\)) | 2.99 |
| 7학년 (\(j=4\)) | 3.47 |
| 전체 \(\bar{y}_{..}\) | 2.535 |
8.2 직교 다항 추정값 계산
\(\hat{\boldsymbol{\theta}} = \mathbf{P}_4 \bar{\mathbf{y}}\) 에서 \(\bar{\mathbf{y}} = (1.14, 2.54, 2.99, 3.47)'\):
\[ \begin{aligned} \hat{\theta}_0 &= \mathbf{p}_0'\bar{\mathbf{y}} = 0.5(1.14 + 2.54 + 2.99 + 3.47) = 5.070 \\ \hat{\theta}_1 &= \mathbf{p}_1'\bar{\mathbf{y}} = \frac{-3(1.14) - 1(2.54) + 1(2.99) + 3(3.47)}{\sqrt{20}} = \frac{9.35}{\sqrt{20}} = 1.666 \\ \hat{\theta}_2 &= \mathbf{p}_2'\bar{\mathbf{y}} = 0.5(1.14 - 2.54 - 2.99 + 3.47) = -0.461 \\ \hat{\theta}_3 &= \mathbf{p}_3'\bar{\mathbf{y}} = \frac{-1(1.14) + 3(2.54) - 3(2.99) + 1(3.47)}{\sqrt{20}} = 0.222 \end{aligned} \]
| 성분 | 추정값 | 해석 |
|---|---|---|
| 상수 \(\hat{\theta}_0\) | 5.070 | 전체 평균의 스케일 |
| 선형 \(\hat{\theta}_1\) | +1.666 | 강한 양의 선형 증가: 학년이 높아질수록 어휘 증가 |
| 2차 \(\hat{\theta}_2\) | -0.461 | 음의 2차: 증가 속도가 학년이 높아질수록 감소 |
| 3차 \(\hat{\theta}_3\) | 0.222 | 약한 3차, 유의하지 않을 가능성 |
8.3 SST* 대각원소 계산
\[ sst_k = N \hat{\theta}_k^2 = 64 \hat{\theta}_k^2 \]
\[ sst_1 = 64 \times (1.666)^2 = 64 \times 2.776 = 177.6 \]
\[ sst_2 = 64 \times (-0.461)^2 = 64 \times 0.213 = 13.6 \]
\[ sst_3 = 64 \times (0.222)^2 = 64 \times 0.049 = 3.2 \]
\[ SS_T = sst_1 + sst_2 + sst_3 = 177.6 + 13.6 + 3.2 = 194.4 \]
8.4 SSR* 대각원소 (교재 Table 3.1에서)
\[ ssr_0 = 116.2 \quad \text{(피험자 SS)} \]
\[ ssr_1 = 50.4, \quad ssr_2 = 44.0, \quad ssr_3 = 60.6 \]
\[ SS_R = 50.4 + 44.0 + 60.6 = 155.0 \]
8.5 다변량 시간 효과 검정
행렬식 방정식 (3.5)에서 Roy 최대근 \(\lambda_1 = 4.743\) (교재 계산값).
\[ \Lambda = \frac{1}{1 + 4.743} = \frac{1}{5.743} = 0.1741 \]
F-변환 (Finn, 1974 근사):
\[ F(3, 61) = \frac{(1-\Lambda)/3}{\Lambda/(N-n+1)} \approx 96.45, \quad p < 0.0001 \]
\(H_0\) 기각 — 학년에 따른 어휘 점수 변화가 매우 유의하다.
8.6 개별 추세 검정 비교
방법 A (구형성 충족 가정, RM ANOVA):
\[ MS_R^{(A)} = \frac{ssr_1 + ssr_2 + ssr_3}{(n-1)(N-1)} = \frac{155.0}{3 \times 63} = \frac{155.0}{189} = 0.820 \]
\[ F_1^{(A)} = \frac{177.6}{0.820} = 216.6, \quad df(1, 189) \]
방법 B (MANOVA, 각자 분모):
\[ F_1^{(B)} = \frac{177.6}{50.4/63} = \frac{177.6}{0.800} = 222.0, \quad df(1, 63) \]
| 추세 | \(sst_k\) | \(ssr_k\) | \(F^{(A)}\) (df=189) | \(F^{(B)}\) (df=63) | 해석 |
|---|---|---|---|---|---|
| 선형 | 177.6 | 50.4 | 216.6 *** | 222.0 *** | 강한 선형 증가 |
| 2차 | 13.6 | 44.0 | 16.6 *** | 19.5 *** | 유의한 감속 |
| 3차 | 3.2 | 60.6 | 3.9 * | 3.3 ns | 약한 3차, 경계선 |
방법 A의 \(F\) 값이 방법 B보다 낮지만, 기각역( \(F_{0.05}(1, 189) \approx 3.89\) )도 더 작다. 3차 추세는 방법 A에서 유의( \(p \approx 0.05\) ), 방법 B에서 비유의( \(p > 0.05\) )로 결론이 갈린다. 이 데이터에서 구형성이 기각되지 않으므로 방법 A가 더 적절하다.
결론: 어휘 점수는 선형으로 증가하되( \(p < 0.001\) ), 증가 속도가 둔화된다( \(p < 0.001\) ). 3차 변곡 패턴은 통계적으로 불분명하다.
9 코드 예시
9.1 R: 수동 구현 + 내장 함수
library(tidyverse)
library(MASS)
# 어휘 성장 데이터 재현 (N=64, n=4)
set.seed(42)
N <- 64
# 비구조적 공분산 행렬 (교재 Table 2.1 기반 근사)
Sigma <- matrix(c(
1.500, 0.800, 0.750, 0.700,
0.800, 1.800, 0.900, 0.850,
0.750, 0.900, 1.600, 0.950,
0.700, 0.850, 0.950, 1.700
), nrow = 4)
mu_true <- c(1.14, 2.54, 2.99, 3.47)
Y <- mvrnorm(N, mu = mu_true, Sigma = Sigma)
colnames(Y) <- paste0("grade", 4:7)
# ──────────────────────────────────────────
# Step 1: 직교 다항식 행렬 P4 구성 (원리 이해)
# ──────────────────────────────────────────
# 방법 1: 수동 Cholesky
n <- 4
times <- 1:n
# T 행렬 (4x4): 상수, 선형, 2차, 3차
T_mat <- outer(0:(n-1), times, function(d, t) t^d)
S <- t(chol(T_mat %*% t(T_mat))) # 하삼각 Cholesky 인수
P <- solve(S) %*% T_mat
cat("Cholesky로 유도한 P4:\n")
print(round(P, 4))
# 방법 2: 내장 poly() 활용
P_builtin <- t(poly(1:n, degree = n-1, raw = FALSE))
# (정규화 방식이 다를 수 있으므로 수동 구현을 권장)
# ──────────────────────────────────────────
# Step 2: 직교 다항 추정값 theta_hat
# ──────────────────────────────────────────
y_bar <- colMeans(Y) # 시점별 평균
theta_hat <- P %*% y_bar # 직교 다항 추정값
cat("\n직교 다항 추정값:\n")
cat("상수:", round(theta_hat[1], 4), "\n")
cat("선형:", round(theta_hat[2], 4), "\n")
cat("2차:", round(theta_hat[3], 4), "\n")
cat("3차:", round(theta_hat[4], 4), "\n")
# ──────────────────────────────────────────
# Step 3: SST*, SSR* 행렬 계산
# ──────────────────────────────────────────
SST_star <- N * P %*% outer(y_bar, y_bar) %*% t(P)
SSR_star <- P %*% (t(Y) %*% Y - N * outer(y_bar, y_bar)) %*% t(P)
cat("\nSST* 대각원소:\n")
cat("상수:", round(diag(SST_star)[1], 3), "\n")
cat("선형:", round(diag(SST_star)[2], 3), "\n")
cat("2차:", round(diag(SST_star)[3], 3), "\n")
cat("3차:", round(diag(SST_star)[4], 3), "\n")
cat("\nSSR* 대각원소:\n")
cat("피험자:", round(diag(SSR_star)[1], 3), "\n")
cat("피험자×선형:", round(diag(SSR_star)[2], 3), "\n")
cat("피험자×2차:", round(diag(SSR_star)[3], 3), "\n")
cat("피험자×3차:", round(diag(SSR_star)[4], 3), "\n")
# ──────────────────────────────────────────
# Step 4: 반복측정 ANOVA SS 추출
# ──────────────────────────────────────────
SS_S <- diag(SSR_star)[1] # 피험자 SS
SS_T <- sum(diag(SST_star)[-1]) # 시간 SS (선형+2차+3차)
SS_R <- sum(diag(SSR_star)[-1]) # 잔차 SS
df_S <- N - 1; df_T <- n - 1; df_R <- (N-1)*(n-1)
MS_S <- SS_S / df_S; MS_T <- SS_T / df_T; MS_R <- SS_R / df_R
F_T <- MS_T / MS_R
cat("\n=== 반복측정 ANOVA 결과 (SST*, SSR*에서 추출) ===\n")
cat(sprintf("피험자 SS=%.2f, df=%d, MS=%.4f\n", SS_S, df_S, MS_S))
cat(sprintf("시간 SS=%.2f, df=%d, MS=%.4f, F=%.2f\n", SS_T, df_T, MS_T, F_T))
cat(sprintf("잔차 SS=%.2f, df=%d, MS=%.4f\n", SS_R, df_R, MS_R))
# ──────────────────────────────────────────
# Step 5: 다변량 시간 효과 검정
# ──────────────────────────────────────────
# 하 (n-1)×(n-1) 부분행렬 추출
SST_sub <- SST_star[-1, -1]
SSR_sub <- SSR_star[-1, -1]
# Cholesky 분해 → 고유값
E <- t(chol(SSR_sub)) # 하삼각 인수
A <- solve(E) %*% SST_sub %*% t(solve(E)) # E^{-1} SST (E^{-1})'
eigenvalues <- sort(Re(eigen(A)$values), decreasing = TRUE)
lambda1 <- eigenvalues[1]
Wilks_Lambda <- 1 / (1 + lambda1)
Roy_root <- lambda1
HotellingLT <- lambda1
Pillai <- lambda1 / (1 + lambda1)
cat("\n=== 다변량 시간 효과 검정 ===\n")
cat(sprintf("Roy 최대근 (λ1): %.4f\n", Roy_root))
cat(sprintf("Wilks' Λ: %.4f\n", Wilks_Lambda))
cat(sprintf("Hotelling-Lawley: %.4f\n", HotellingLT))
cat(sprintf("Pillai-Bartlett: %.4f\n", Pillai))
# Wilks' Λ → F 변환 (일표본)
F_multi <- (1 - Wilks_Lambda) / Wilks_Lambda * (N - n + 1) / (n - 1)
cat(sprintf("F(%d, %d) = %.2f\n", n-1, N-n+1, F_multi))
# ──────────────────────────────────────────
# Step 6: 개별 추세 검정 — 방법 A vs 방법 B
# ──────────────────────────────────────────
cat("\n=== 개별 추세 성분 검정 ===\n")
trends <- c("선형", "2차", "3차")
for (k in 1:(n-1)) {
sst_k <- diag(SST_star)[k+1]
ssr_k <- diag(SSR_star)[k+1]
# 방법 A: 공통 분모 MS_R
F_A <- sst_k / MS_R
pval_A <- pf(F_A, 1, df_R, lower.tail = FALSE)
# 방법 B: 개별 분모
F_B <- sst_k / (ssr_k / (N-1))
pval_B <- pf(F_B, 1, N-1, lower.tail = FALSE)
cat(sprintf("%s: F_A=%.2f(p=%.4f) | F_B=%.2f(p=%.4f)\n",
trends[k], F_A, pval_A, F_B, pval_B))
}
# ──────────────────────────────────────────
# Step 7: 내장 함수로 검증
# ──────────────────────────────────────────
# MANOVA (Wide format)
fit_manova <- manova(Y ~ 1)
cat("\n내장 manova() 결과 (Pillai):\n")
print(summary(fit_manova, test = "Pillai"))
# 반복측정 ANOVA (Long format)
long_df <- as.data.frame(Y) |>
mutate(subject = factor(1:N)) |>
pivot_longer(-subject, names_to = "grade", values_to = "score") |>
mutate(time = factor(grade))
fit_aov <- aov(score ~ time + Error(subject/time), data = long_df)
cat("\n내장 aov() 결과 (RM ANOVA):\n")
print(summary(fit_aov))9.2 Python: 수동 구현 + statsmodels
import numpy as np
import pandas as pd
from scipy import linalg, stats
# 데이터 생성
np.random.seed(42)
N, n = 64, 4
Sigma = np.array([
[1.5, 0.8, 0.75, 0.70],
[0.8, 1.8, 0.90, 0.85],
[0.75, 0.90, 1.6, 0.95],
[0.70, 0.85, 0.95, 1.7]
])
mu_true = np.array([1.14, 2.54, 2.99, 3.47])
rng = np.random.default_rng(42)
L = linalg.cholesky(Sigma, lower=True)
Y = mu_true + rng.standard_normal((N, n)) @ L.T
# ──────────────────────────────────────────
# Step 1: 직교 다항식 행렬 P (Cholesky 방법)
# ──────────────────────────────────────────
times = np.arange(1, n + 1)
# T: (n x n) 멱행렬
T_mat = np.array([times**d for d in range(n)]) # (n, n)
S = linalg.cholesky(T_mat @ T_mat.T, lower=True)
P = linalg.solve_triangular(S, T_mat, lower=True) # S^{-1} T
print("직교 다항식 행렬 P:")
print(np.round(P, 4))
# ──────────────────────────────────────────
# Step 2: 직교 다항 추정값
# ──────────────────────────────────────────
y_bar = Y.mean(axis=0)
theta_hat = P @ y_bar
print("\n직교 다항 추정값:")
for label, val in zip(["상수", "선형", "2차", "3차"], theta_hat):
print(f" {label}: {val:.4f}")
# ──────────────────────────────────────────
# Step 3: SSCP 행렬
# ──────────────────────────────────────────
SST_star = N * P @ np.outer(y_bar, y_bar) @ P.T
SSR_star = P @ (Y.T @ Y - N * np.outer(y_bar, y_bar)) @ P.T
print("\nSST* 대각원소:", np.round(np.diag(SST_star), 3))
print("SSR* 대각원소:", np.round(np.diag(SSR_star), 3))
# ──────────────────────────────────────────
# Step 4: RM ANOVA SS 추출
# ──────────────────────────────────────────
SS_S = np.diag(SSR_star)[0]
SS_T = np.sum(np.diag(SST_star)[1:])
SS_R = np.sum(np.diag(SSR_star)[1:])
df_S, df_T, df_R = N-1, n-1, (N-1)*(n-1)
MS_T, MS_R = SS_T/df_T, SS_R/df_R
F_T = MS_T / MS_R
print(f"\n=== RM ANOVA 추출 ===")
print(f"SS_T={SS_T:.2f}, MS_T={MS_T:.4f}, F={F_T:.2f}")
print(f"SS_R={SS_R:.2f}, MS_R={MS_R:.4f}, df_R={df_R}")
# ──────────────────────────────────────────
# Step 5: 다변량 시간 효과 검정
# ──────────────────────────────────────────
SST_sub = SST_star[1:, 1:]
SSR_sub = SSR_star[1:, 1:]
# Cholesky → 고유값 문제
E = linalg.cholesky(SSR_sub, lower=True)
A = linalg.solve_triangular(E, SST_sub, lower=True)
A = linalg.solve_triangular(E, A.T, lower=True).T # E^{-1} SST (E^{-1})'
eigenvalues = np.real(linalg.eigvals(A))
eigenvalues = np.sort(eigenvalues[eigenvalues > 1e-10])[::-1]
lambda1 = eigenvalues[0]
Wilks_L = 1 / (1 + lambda1)
F_multi = (1 - Wilks_L) / Wilks_L * (N - n + 1) / (n - 1)
p_multi = 1 - stats.f.cdf(F_multi, n-1, N-n+1)
print(f"\n=== 다변량 검정 ===")
print(f"Roy 최대근: {lambda1:.4f}")
print(f"Wilks' Λ: {Wilks_L:.4f}")
print(f"F({n-1}, {N-n+1}) = {F_multi:.2f}, p = {p_multi:.6f}")
# ──────────────────────────────────────────
# Step 6: 개별 추세 F-검정
# ──────────────────────────────────────────
print("\n=== 개별 추세 검정 ===")
for k, trend in enumerate(["선형", "2차", "3차"], 1):
sst_k = np.diag(SST_star)[k]
ssr_k = np.diag(SSR_star)[k]
F_A = sst_k / MS_R
p_A = 1 - stats.f.cdf(F_A, 1, df_R)
F_B = sst_k / (ssr_k / (N-1))
p_B = 1 - stats.f.cdf(F_B, 1, N-1)
print(f"{trend}: F_A={F_A:.2f} (df={df_R}, p={p_A:.4f}) | "
f"F_B={F_B:.2f} (df={N-1}, p={p_B:.4f})")
# ──────────────────────────────────────────
# Step 7: statsmodels MANOVA로 검증
# ──────────────────────────────────────────
from statsmodels.multivariate.manova import MANOVA
df = pd.DataFrame(Y, columns=["g4", "g5", "g6", "g7"])
# 일표본 검정: intercept만 포함
maov = MANOVA.from_formula("g4 + g5 + g6 + g7 ~ 1", data=df)
result = maov.mv_test()
print("\nstatsmodels MANOVA 결과:")
print(result.summary())10 핵심 수식 요약
| 단계 | 수식 | 목적 |
|---|---|---|
| 모형 | \(\mathbf{y}_i = \boldsymbol{\mu} + \boldsymbol{\varepsilon}_i, \; \boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma})\) | 비구조적 공분산 허용 |
| P 유도 | \(\mathbf{P} = \mathbf{S}^{-1}\mathbf{T}, \; \mathbf{T}\mathbf{T}' = \mathbf{S}\mathbf{S}'\) | 직교 다항식 변환 |
| 변환 모형 | \(\mathbf{P}\mathbf{y}_i = \boldsymbol{\theta} + \boldsymbol{\varepsilon}_i^*\) | 추세 성분 분리 |
| SST* | \(N\mathbf{P}\bar{\mathbf{y}}\bar{\mathbf{y}}'\mathbf{P}'\) | 시간(추세) 제곱합 행렬 |
| SSR* | \(\mathbf{P}(\mathbf{Y}'\mathbf{Y}-N\bar{\mathbf{y}}\bar{\mathbf{y}}')\mathbf{P}'\) | 잔차 제곱합 행렬 |
| RM ANOVA 추출 | \(SS_T = \sum_{k=1}^{n-1}sst_k\), \(SS_R = \sum_{k=1}^{n-1}ssr_k\) | 구형성 충족 시 |
| 다변량 검정 | \(\|\mathbf{E}^{-1}\mathbf{SST}_{(n-1)}^*(\mathbf{E}^{-1})' - \lambda\mathbf{I}\| = 0\) | 고유값으로 검정 |
| Wilks’ Λ | \(\Lambda = 1/(1+\lambda_1)\) | F로 근사 변환 |
| 추세 검정 A | \(F_k = sst_k / MS_R\), df = \((n-1)(N-1)\) | 구형성 충족 (더 강력) |
| 추세 검정 B | \(F_k = sst_k / [ssr_k/(N-1)]\), df = \(N-1\) | 구형성 위반 |
11 관련 주제
선행 지식
후속 주제
관련 개념