1 개요
일표본 반복측정 MANOVA는 단일 집단이 시간에 따라 어떻게 변화하는지를 검정한다. 그러나 실제 연구에서 핵심 질문은 흔히 “두 집단의 변화 패턴이 다른가?” 다.
- 약물 집단과 위약 집단의 우울 점수가 시간에 따라 다르게 변화하는가?
- 교육 방식 A와 B를 받은 학생들의 성취도 성장 궤적이 다른가?
- 저위험·고위험 환자의 바이오마커 변화 패턴이 다른가?
이것이 다표본(s-sample) 반복측정 MANOVA가 답하는 질문이다.
일표본 MANOVA가 “시간 효과”만 검정한다면, 다표본 MANOVA는 “집단 효과”, “시간 효과”, “집단 × 시간 교호작용” 세 가지를 동시에 검정한다. 그 중 가장 중요한 것이 집단 × 시간 교호작용이다.
2 다표본 MANOVA 모형
\(s\) 개 집단( \(h = 1, \ldots, s\) ), 집단 \(h\) 에 \(N_h\) 명, 총 피험자 \(N = \sum_{h=1}^s N_h\), 각 피험자는 \(n\) 시점( \(j = 1, \ldots, n\) )에 완전 측정된다고 가정한다.
\[ \mathbf{y}_{hi} = \boldsymbol{\mu} + \boldsymbol{\gamma}_h + \boldsymbol{\varepsilon}_{hi}, \quad h = 1,\ldots,s, \quad i = 1,\ldots,N_h \tag{3.8} \]
- \(\mathbf{y}_{hi} \in \mathbb{R}^n\): 집단 \(h\) 의 피험자 \(i\) 의 반복 측정 벡터
- \(\boldsymbol{\mu} \in \mathbb{R}^n\): 전체 시점별 평균 벡터 (집단 공통)
- \(\boldsymbol{\gamma}_h \in \mathbb{R}^n\): 집단 \(h\) 의 효과 벡터 (시점별 집단 효과)
- \(\boldsymbol{\varepsilon}_{hi} \sim N(\mathbf{0}, \boldsymbol{\Sigma})\): 집단 공통 오차 공분산 행렬
2.1 모형의 핵심: 집단 효과가 벡터다
일원분산분석에서 집단 효과는 스칼라 \(\gamma_h\) 다 — “집단이 얼마나 높고 낮은가”.
다표본 반복측정 MANOVA에서 집단 효과는 \(n \times 1\) 벡터 \(\boldsymbol{\gamma}_h\) 다.
\[ \boldsymbol{\gamma}_h = (\gamma_{h1}, \gamma_{h2}, \ldots, \gamma_{hn})' \]
\(\gamma_{hj}\) 는 시점 \(j\) 에서 집단 \(h\) 의 추가 효과다. \(\boldsymbol{\gamma}_h\) 가 집단별로 다르면 집단 × 시간 교호작용이 존재한다.
직관: 우울증 임상시험을 예로 들자. 약물 집단 \(\boldsymbol{\gamma}_{\text{drug}} = (-2, -4, -5, -6)\) 이고 위약 집단 \(\boldsymbol{\gamma}_{\text{placebo}} = (-1, -1, -1, -1)\) 이면, 두 집단은 1시점에서는 비슷하게 감소하지만, 4시점에서 큰 차이가 생긴다. 이것이 교호작용의 직관이다.
2.2 동분산 가정
모형은 집단 간 공분산 행렬 동질성(homogeneity of variance-covariance) 을 가정한다:
\[ \text{Cov}(\mathbf{y}_{hi}) = \boldsymbol{\Sigma} \quad \text{모든 집단 } h \text{에 대해 동일} \]
이 가정을 Box’s M 검정으로 확인할 수 있으나, 표본 크기가 크면 Box’s M 검정 자체가 민감해지므로 주의해야 한다.
3 직교 다항 변환
일표본과 마찬가지로 \(\mathbf{P}\) 로 모형을 변환한다:
\[ \mathbf{P}\mathbf{y}_{hi} = \mathbf{P}\boldsymbol{\mu} + \mathbf{P}\boldsymbol{\gamma}_h + \mathbf{P}\boldsymbol{\varepsilon}_{hi} = \boldsymbol{\theta} + \boldsymbol{\delta}_h + \boldsymbol{\varepsilon}_{hi}^* \tag{3.9} \]
여기서
- \(\boldsymbol{\theta} = \mathbf{P}\boldsymbol{\mu}\): 전체 평균의 직교 다항 계수 벡터
- \(\boldsymbol{\delta}_h = \mathbf{P}\boldsymbol{\gamma}_h\): 집단 \(h\) 효과의 직교 다항 계수 벡터
- \(\boldsymbol{\varepsilon}_{hi}^* \sim N(\mathbf{0}, \boldsymbol{\Sigma}^* = \mathbf{P}\boldsymbol{\Sigma}\mathbf{P}')\)
\(\boldsymbol{\delta}_h\) 의 각 원소가 갖는 의미:
| 원소 | 의미 |
|---|---|
| \(\delta_{h0}\) | 집단 \(h\) 의 전반적 평균 수준 차이 |
| \(\delta_{h1}\) | 집단 \(h\) 의 선형 시간 추세 차이 |
| \(\delta_{h2}\) | 집단 \(h\) 의 2차(가속/감속) 추세 차이 |
| \(\delta_{h, n-1}\) | 집단 \(h\) 의 \((n-1)\) 차 추세 차이 |
교호작용의 다항식 언어: 집단 × 시간 교호작용은 “\(\boldsymbol{\delta}_h\) 의 \(h\) 에 따른 차이”다. 특히 \(\delta_{h1}\) 들이 집단 간에 다르면 → 집단별 선형 성장률이 다르다(교호작용). \(\delta_{h2}\) 들이 다르면 → 집단별 가속/감속 패턴이 다르다.
4 세 SSCP 행렬 체계
다표본 MANOVA는 세 개의 SSCP 행렬로 전체 변동을 분해한다. 각 행렬이 담는 정보의 층이 다르다.
| 변동원 | df | SSCP | 설명 |
|---|---|---|---|
| 시간 (Time) | 1 | \(\mathbf{SST}^* = N\mathbf{P}\bar{\mathbf{y}}_{..}\bar{\mathbf{y}}_{..}'\mathbf{P}'\) | 전체 평균 궤적의 변동 |
| 집단 간 (Group) | \(s-1\) | \(\mathbf{SSG}^* = \mathbf{P}\!\left(\sum_h N_h \bar{\mathbf{y}}_{h.}\bar{\mathbf{y}}_{h.}' - N\bar{\mathbf{y}}_{..}\bar{\mathbf{y}}_{..}'\right)\!\mathbf{P}'\) | 집단 평균 궤적의 변동 |
| 잔차 (Subj. within Groups) | \(N-s\) | \(\mathbf{SSR}^* = \mathbf{P}\!\left(\sum_h\sum_i \mathbf{y}_{hi}\mathbf{y}_{hi}' - \sum_h N_h\bar{\mathbf{y}}_{h.}\bar{\mathbf{y}}_{h.}'\right)\!\mathbf{P}'\) | 집단 내 피험자 간 변동 |
\(\bar{\mathbf{y}}_{h.}\) 은 집단 \(h\) 의 시점별 평균 벡터, \(\bar{\mathbf{y}}_{..}\) 는 전체 시점별 평균 벡터다.
4.1 SSG* 행렬의 유도
SSG는 집단 간 변동(Between-Groups SS)의 행렬 버전이다.
\[ \mathbf{SSG} = \sum_{h=1}^s N_h(\bar{\mathbf{y}}_{h.} - \bar{\mathbf{y}}_{..})(\bar{\mathbf{y}}_{h.} - \bar{\mathbf{y}}_{..})' = \sum_h N_h \bar{\mathbf{y}}_{h.}\bar{\mathbf{y}}_{h.}' - N\bar{\mathbf{y}}_{..}\bar{\mathbf{y}}_{..}' \]
직교 변환을 적용하면
\[ \mathbf{SSG}^* = \mathbf{P}\mathbf{SSG}\mathbf{P}' = \mathbf{P}\!\left(\sum_h N_h \bar{\mathbf{y}}_{h.}\bar{\mathbf{y}}_{h.}'\right)\!\mathbf{P}' - N\mathbf{P}\bar{\mathbf{y}}_{..}\bar{\mathbf{y}}_{..}'\mathbf{P}' \]
4.2 대각원소의 의미 체계 (3.10)
세 SSCP 행렬의 대각원소가 각각 무엇을 의미하는지가 이 분석의 핵심이다.
\(\mathbf{SST}^*\) 대각원소:
\[ \text{diag}(\mathbf{SST}^*) = (sst_0, sst_1, sst_2, \ldots, sst_{n-1}) \]
- \(sst_0 = N \cdot n\bar{y}_{..}^2\): 전체 평균의 제곱합
- \(sst_k (k \geq 1)\): 전체 평균 궤적의 \(k\) 차 시간 추세 성분
\(\mathbf{SSG}^*\) 대각원소 — 이것이 핵심이다:
\[ \text{diag}(\mathbf{SSG}^*) = (ssg_0, ssg_1, ssg_2, \ldots, ssg_{n-1}) \]
- \(ssg_0\): 집단 주효과 — 집단 간 전반적 수준 차이
- \(ssg_1\): 집단 × 선형 교호작용 — 집단 간 선형 성장률 차이
- \(ssg_2\): 집단 × 2차 교호작용 — 집단 간 가속/감속 패턴 차이
- \(ssg_k (k \geq 1)\): 집단 × \(k\) 차 시간 교호작용
\(\mathbf{SSR}^*\) 대각원소:
\[ \text{diag}(\mathbf{SSR}^*) = (ssr_0, ssr_1, ssr_2, \ldots, ssr_{n-1}) \]
- \(ssr_0\): 집단 내 피험자 간(Subjects within Groups) 분산
- \(ssr_k (k \geq 1)\): \(k\) 차 추세의 피험자-내 잔차(피험자의 \(k\) 차 추세가 집단 평균에서 벗어나는 정도)
직관 요약:
| 행렬 | “누구의” 변동 | “무엇에 대한” 변동 |
|---|---|---|
| \(\mathbf{SST}^*\) | 전체 집단 | 시간 추세 |
| \(\mathbf{SSG}^*\) | 집단 간 | 집단 주효과 + 집단×시간 교호작용 |
| \(\mathbf{SSR}^*\) | 집단 내 피험자 간 | 피험자 평균 수준 + 피험자-내 추세 이질성 |
5 분할구획 ANOVA 결과 추출
구형성이 성립하면, 세 SSCP 행렬에서 분할구획 ANOVA의 모든 SS를 직접 추출한다.
5.1 SS 추출 공식
\[ \begin{aligned} SS_T &= \sum_{k=1}^{n-1} sst_k \quad &\leftarrow \mathbf{SST}^* \text{의 하 } n-1 \text{개 대각원소 합}\\ SS_G &= ssg_0 \quad &\leftarrow \mathbf{SSG}^* \text{의 첫 번째 대각원소} \\ SS_{GT} &= \sum_{k=1}^{n-1} ssg_k \quad &\leftarrow \mathbf{SSG}^* \text{의 하 } n-1 \text{개 대각원소 합} \\ SS_{S(G)} &= ssr_0 \quad &\leftarrow \mathbf{SSR}^* \text{의 첫 번째 대각원소} \\ SS_R &= \sum_{k=1}^{n-1} ssr_k \quad &\leftarrow \mathbf{SSR}^* \text{의 하 } n-1 \text{개 대각원소 합} \end{aligned} \]
왜 \(SS_G = ssg_0\) 인가?
\(\boldsymbol{\delta}_h\) 의 첫 번째 원소 \(\delta_{h0} = \mathbf{p}_0'\boldsymbol{\gamma}_h = \frac{1}{\sqrt{n}}\sum_j \gamma_{hj}\) 는 집단 \(h\) 의 평균 수준 차이를 나타낸다. \(ssg_0 = \sum_h N_h(\delta_{h0} - \bar{\delta}_{.0})^2\) 은 집단 간 이 평균 수준 차이의 제곱합이므로, 곧 집단 주효과 SS다.
왜 \(SS_{GT} = \sum_{k=1}^{n-1} ssg_k\) 인가?
\(ssg_k (k \geq 1)\) 은 “\(k\) 차 추세에서 집단 간 차이”이므로, 이것이 집단 × \(k\) 차 추세 교호작용 SS다. 이들의 합이 전체 집단 × 시간 교호작용 SS가 된다.
5.2 분할구획 ANOVA 분산분해표
| 변동원 | df | SS | MS | F | 분모 |
|---|---|---|---|---|---|
| 시간 | \(n-1\) | \(SS_T\) | \(MS_T = SS_T/(n-1)\) | \(MS_T/MS_R\) | \(MS_R\) |
| 집단 | \(s-1\) | \(SS_G\) | \(MS_G = SS_G/(s-1)\) | \(MS_G/MS_{S(G)}\) | \(MS_{S(G)}\) |
| 집단 × 시간 | \((s-1)(n-1)\) | \(SS_{GT}\) | \(MS_{GT} = SS_{GT}/[(s-1)(n-1)]\) | \(MS_{GT}/MS_R\) | \(MS_R\) |
| 집단 내 피험자 | \(N-s\) | \(SS_{S(G)}\) | \(MS_{S(G)} = SS_{S(G)}/(N-s)\) | \(MS_{S(G)}/MS_R\) | \(MS_R\) |
| 잔차 | \((N-s)(n-1)\) | \(SS_R\) | \(MS_R = SS_R/[(N-s)(n-1)]\) |
5.3 왜 집단 효과의 분모가 \(MS_{S(G)}\) 인가?
이 질문이 분할구획 설계의 핵심이다.
집단 효과 검정: 집단 간 차이를 측정하려면 피험자 간(between-subjects) 변동을 사용해야 한다. \(MS_{S(G)}\) 는 “집단 내 피험자들이 서로 얼마나 다른가”를 나타내므로, 올바른 분모다.
시간 효과, 집단 × 시간 교호작용 검정: 이는 피험자 내(within-subjects) 변동이다. 같은 피험자 내에서 시점들을 비교하므로, 피험자-내 잔차 \(MS_R\) 이 분모가 된다.
직관: 분할구획(split-plot) 설계에서는 두 층위의 오차가 존재한다. - 전포구(whole-plot) 오차: 피험자 간 차이 = \(MS_{S(G)}\) → 집단 효과의 분모 - 분할구(subplot) 오차: 피험자 내 시점 간 차이 = \(MS_R\) → 시간 및 교호작용의 분모
이 두 오차를 구분하지 않으면 집단 효과 검정이 부정확해진다.
6 다변량 검정
6.1 검정 대상: 세 가지 효과
다표본 MANOVA에서 세 가지 다변량 검정이 필요하다.
| 효과 | 귀무가설 | 검정 방법 |
|---|---|---|
| 집단 주효과 | \(H_0: ssg_0 = 0\) (집단 간 전반적 수준 같음) | 일변량 F-검정 (피험자-간 비교) |
| 시간 효과 | \(H_0: \boldsymbol{\tau} = \mathbf{0}\) (모든 시점 평균 같음) | 다변량 F-검정 (3.13) |
| 집단 × 시간 | \(H_0: \boldsymbol{\gamma}_1 = \boldsymbol{\gamma}_2 = \cdots = \boldsymbol{\gamma}_s\) | 다변량 F-검정 (3.11) |
집단 주효과는 피험자-간 효과이므로, 다변량 검정 대신 일변량 F-검정이 사용된다.
6.2 집단 × 시간 교호작용 검정 (가장 중요)
\(\mathbf{SSG}^*\) 와 \(\mathbf{SSR}^*\) 의 하 \((n-1) \times (n-1)\) 부분행렬을 추출한 후 행렬식 방정식을 푼다:
\[ |\mathbf{SSG}_{(n-1)}^* - \lambda\, \mathbf{SSR}_{(n-1)}^*| = 0 \tag{3.11} \]
이 방정식은 \(r = \min(s-1, n-1)\) 개의 비영(non-zero) 고유값 \(\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_r > 0\) 을 갖는다.
왜 고유값이 \(\min(s-1, n-1)\) 개인가?
- \(\mathbf{SSG}_{(n-1)}^*\) 의 랭크는 최대 \(s-1\) (집단 간 대비 수)
- 부분행렬의 차원은 \(n-1\)
- 따라서 비영 고유값의 수는 두 값 중 작은 것
예: \(s = 2\) 집단, \(n = 4\) 시점이면 \(\mathbf{SSG}^*\) 의 랭크 = 1. 고유값이 하나뿐이므로, 일표본과 동일하게 단일 고유값으로 검정.
\[ \text{고유값이 하나} \Rightarrow \Lambda = \frac{1}{1+\lambda_1} \quad (s=2 \text{ 특수 케이스}) \]
일반적으로 Wilks’ \(\Lambda\) 는
\[ \Lambda = \prod_{k=1}^{r} \frac{1}{1+\lambda_k} \tag{3.12} \]
직관: \(\Lambda\) 가 \(s-1\) 개 고유값의 곱으로 표현되는 이유는, 각 고유값이 독립적인 “분리 방향”에서의 집단 차이를 표현하기 때문이다. 모든 방향에서 집단이 분리되어 있어야 \(\Lambda\) 가 작아진다.
6.3 시간 효과 검정 (다표본)
\(H_0: \boldsymbol{\tau} = \mathbf{0}\) (시간 효과 없음). 일표본과 동일한 행렬식 방정식:
\[ |\mathbf{SST}_{(n-1)}^* - \lambda\, \mathbf{SSR}_{(n-1)}^*| = 0 \tag{3.13} \]
중요한 점: 다표본에서 잔차 행렬 \(\mathbf{SSR}^*\) 는 집단 내 변동( \(df = N-s\) )을 담는다. 일표본의 \(\mathbf{SSR}^*\) ( \(df = N-1\) )보다 자유도가 줄었지만, 집단 구분을 통해 불순물(집단 간 차이)을 제거했으므로 더 순수한 잔차다.
교호작용이 유의하지 않으면? 이상적으로는 교호작용을 잔차에 풀링(pooling)하여 시간 효과 검정의 검정력을 높이고 싶지만, 대부분의 MANOVA 소프트웨어는 이 풀링을 직접 지원하지 않는다. 따라서 실무에서는 교호작용이 비유의일 때도 (3.13)을 그대로 사용한다.
6.4 집단 × 시간 개별 추세 성분 검정
전체 교호작용 검정 후, 어느 추세 성분에서 집단 차이가 발생하는지 파악한다.
구형성 위반 시 (MANOVA 방식), 각 추세 성분에 고유한 분모를 사용:
\[ F_{GT_k} = \frac{ssg_k / (s-1)}{ssr_k / (N-s)}, \quad k = 1, 2, \ldots, n-1 \tag{3.14} \]
분모 자유도: \(N - s\).
구형성 충족 시 (RM ANOVA 방식), 공통 분모 \(MS_R\) 을 사용:
\[ F_{GT_k}^{(A)} = \frac{ssg_k / (s-1)}{MS_R}, \quad df_\text{분모} = (N-s)(n-1) \]
직관 — 왜 분모가 \(ssr_k\) 인가?
\(ssr_k\) 는 ” \(k\) 차 추세에서 피험자들의 개인 간 차이”를 나타낸다. 집단 × \(k\) 차 추세 교호작용( \(ssg_k\) )을 검정하려면, 같은 추세 차원에서 피험자들이 자연적으로 얼마나 다른지( \(ssr_k\) )를 분모로 써야 한다. 다른 차수의 오차(\(ssr_1\) 을 \(ssg_2\) 의 분모로 쓰는 것)는 차원이 다른 비교를 하는 것이다.
7 부분 추세 모형: \(q < n\)
\(n\) 시점이지만 \(q < n\) 차수의 다항식만 검정하고 싶을 때 (예: 4시점이지만 선형·2차만 의미있다고 판단하는 경우), \(n - q - 1\) 개의 하부 행렬만 추출하여 검정한다.
구체적으로, \(q = 2\) (선형·2차)만 고려하면 \(\mathbf{SSG}^*\) 와 \(\mathbf{SSR}^*\) 의 \((n-q-1) \times (n-q-1)\) 부분행렬을 비교한다.
이를 통해 이론적으로 의미 없는 고차 추세를 포함하지 않고 검정력을 높일 수 있다.
8 수치 예제: 임상시험 시뮬레이션
\(s = 2\) 집단(약물/위약), \(N = 60\) 명( \(N_{\text{drug}} = N_{\text{placebo}} = 30\) ), \(n = 4\) 시점(0주, 2주, 4주, 8주), 반응 변수: HAM-D 우울 척도.
8.1 설계 목표
- 약물 집단: 강한 선형 감소 + 약한 가속 감소
- 위약 집단: 완만한 선형 감소
- 기대 결과: 집단 × 선형 추세 교호작용 유의
8.2 모집단 설정
| 집단 | 0주 | 2주 | 4주 | 8주 |
|---|---|---|---|---|
| 약물 | 22.0 | 16.0 | 12.0 | 10.0 |
| 위약 | 22.0 | 20.0 | 19.0 | 18.0 |
| 차이 | 0 | -4 | -7 | -8 |
약물과 위약의 차이가 시간이 지날수록 커지는 패턴 → 집단 × 시간 교호작용 존재.
8.3 SSG* 대각원소 계산
직교 다항식 \(\mathbf{p}_1\) 로 집단 평균 벡터를 변환:
\[ \hat{\delta}_{\text{drug}, 1} = \mathbf{p}_1'\bar{\mathbf{y}}_{\text{drug},.} = \frac{-3(22) + (-1)(16) + 1(12) + 3(10)}{\sqrt{20}} = \frac{-28}{\sqrt{20}} = -6.26 \]
\[ \hat{\delta}_{\text{placebo}, 1} = \mathbf{p}_1'\bar{\mathbf{y}}_{\text{placebo},.} = \frac{-3(22) + (-1)(20) + 1(19) + 3(18)}{\sqrt{20}} = \frac{-9}{\sqrt{20}} = -2.01 \]
집단 × 선형 교호작용 SS:
\[ ssg_1 = N_{\text{drug}}(\hat{\delta}_{\text{drug},1} - \bar{\hat{\delta}}_{.1})^2 + N_{\text{placebo}}(\hat{\delta}_{\text{placebo},1} - \bar{\hat{\delta}}_{.1})^2 \]
두 집단이 같은 크기( \(N_h = 30\) )이면 \(\bar{\hat{\delta}}_{.1} = (-6.26 - 2.01)/2 = -4.14\) 이므로
\[ ssg_1 = 30(-6.26 + 4.14)^2 + 30(-2.01 + 4.14)^2 = 30(4.50) + 30(4.55) \approx 270.1 \]
선형 성장률에서 두 집단의 큰 차이가 포착된다.
9 코드 예시
9.1 R: 수동 구현 + 내장 함수
library(tidyverse)
library(MASS)
set.seed(42)
N_g <- 30 # 집단당 피험자 수
n <- 4 # 시점 수
s <- 2 # 집단 수
N <- N_g * s
# 공분산 행렬 (집단 공통)
Sigma <- matrix(c(
16, 8, 7, 6,
8, 14, 8, 7,
7, 8, 12, 8,
6, 7, 8, 11
), nrow = 4)
# 모집단 평균
mu_drug <- c(22, 16, 12, 10)
mu_placebo <- c(22, 20, 19, 18)
# 데이터 생성
Y_drug <- mvrnorm(N_g, mu_drug, Sigma)
Y_placebo <- mvrnorm(N_g, mu_placebo, Sigma)
Y <- rbind(Y_drug, Y_placebo)
colnames(Y) <- paste0("week", c(0, 2, 4, 8))
group <- factor(rep(c("drug", "placebo"), each = N_g))
# ──────────────────────────────────────────
# Step 1: 직교 다항식 행렬 P
# ──────────────────────────────────────────
times <- 1:n
T_mat <- outer(0:(n-1), times, function(d, t) t^d)
S <- t(chol(T_mat %*% t(T_mat)))
P <- solve(S) %*% T_mat
# ──────────────────────────────────────────
# Step 2: 집단별 평균, 전체 평균
# ──────────────────────────────────────────
y_bar_drug <- colMeans(Y_drug)
y_bar_placebo <- colMeans(Y_placebo)
y_bar_all <- colMeans(Y)
# 직교 다항 추정값
delta_drug <- P %*% y_bar_drug
delta_placebo <- P %*% y_bar_placebo
delta_all <- P %*% y_bar_all
cat("직교 다항 추정값 (약물):", round(delta_drug, 3), "\n")
cat("직교 다항 추정값 (위약):", round(delta_placebo, 3), "\n")
# ──────────────────────────────────────────
# Step 3: 세 SSCP 행렬 계산
# ──────────────────────────────────────────
SST_star <- N * P %*% outer(y_bar_all, y_bar_all) %*% t(P)
SSG_raw <- N_g * outer(y_bar_drug, y_bar_drug) +
N_g * outer(y_bar_placebo, y_bar_placebo) -
N * outer(y_bar_all, y_bar_all)
SSG_star <- P %*% SSG_raw %*% t(P)
SSY_raw <- t(Y) %*% Y
SSR_star <- P %*% (SSY_raw - N_g * outer(y_bar_drug, y_bar_drug) -
N_g * outer(y_bar_placebo, y_bar_placebo)) %*% t(P)
cat("\nSST* 대각원소 (시간 분해):\n")
print(round(diag(SST_star), 3))
cat("\nSSG* 대각원소 (집단 분해):\n")
cat("집단 주효과 ssg0:", round(diag(SSG_star)[1], 3), "\n")
cat("집단×선형 ssg1:", round(diag(SSG_star)[2], 3), "\n")
cat("집단×2차 ssg2:", round(diag(SSG_star)[3], 3), "\n")
cat("집단×3차 ssg3:", round(diag(SSG_star)[4], 3), "\n")
cat("\nSSR* 대각원소 (잔차 분해):\n")
cat("집단내 피험자 ssr0:", round(diag(SSR_star)[1], 3), "\n")
cat("피험자×선형 ssr1:", round(diag(SSR_star)[2], 3), "\n")
cat("피험자×2차 ssr2:", round(diag(SSR_star)[3], 3), "\n")
cat("피험자×3차 ssr3:", round(diag(SSR_star)[4], 3), "\n")
# ──────────────────────────────────────────
# Step 4: 분할구획 ANOVA SS 추출
# ──────────────────────────────────────────
SS_T <- sum(diag(SST_star)[-1]) # 시간 SS
SS_G <- diag(SSG_star)[1] # 집단 주효과 SS
SS_GT <- sum(diag(SSG_star)[-1]) # 집단×시간 교호작용 SS
SS_SG <- diag(SSR_star)[1] # 집단내 피험자 SS
SS_R <- sum(diag(SSR_star)[-1]) # 잔차 SS
df_T <- n - 1; df_G <- s - 1; df_GT <- (s-1)*(n-1)
df_SG <- N - s; df_R <- (N-s)*(n-1)
MS_T <- SS_T / df_T; MS_G <- SS_G / df_G
MS_GT <- SS_GT / df_GT; MS_SG <- SS_SG / df_SG; MS_R <- SS_R / df_R
F_T <- MS_T / MS_R # 시간: 분모 = MS_R
F_G <- MS_G / MS_SG # 집단: 분모 = MS_S(G) ← 핵심!
F_GT <- MS_GT / MS_R # 교호작용: 분모 = MS_R
p_T <- pf(F_T, df_T, df_R, lower.tail = FALSE)
p_G <- pf(F_G, df_G, df_SG, lower.tail = FALSE)
p_GT <- pf(F_GT, df_GT, df_R, lower.tail = FALSE)
cat("\n=== 분할구획 ANOVA (SST*, SSG*, SSR*에서 추출) ===\n")
cat(sprintf("시간 SS=%.2f, df=%d, MS=%.2f, F=%.2f, p=%.4f\n",
SS_T, df_T, MS_T, F_T, p_T))
cat(sprintf("집단 SS=%.2f, df=%d, MS=%.2f, F=%.2f, p=%.4f (분모=MS_S(G))\n",
SS_G, df_G, MS_G, F_G, p_G))
cat(sprintf("집단×시간 SS=%.2f, df=%d, MS=%.2f, F=%.2f, p=%.4f\n",
SS_GT, df_GT, MS_GT, F_GT, p_GT))
cat(sprintf("집단내 피험자 SS=%.2f, df=%d, MS=%.2f\n",
SS_SG, df_SG, MS_SG))
cat(sprintf("잔차 SS=%.2f, df=%d, MS=%.2f\n",
SS_R, df_R, MS_R))
# ──────────────────────────────────────────
# Step 5: 다변량 교호작용 검정 (3.11)
# ──────────────────────────────────────────
SSG_sub <- SSG_star[-1, -1] # 하 (n-1)×(n-1) 부분행렬
SSR_sub <- SSR_star[-1, -1]
E <- t(chol(SSR_sub))
A <- solve(E) %*% SSG_sub %*% t(solve(E))
eigs <- sort(Re(eigen(A)$values), decreasing = TRUE)
r <- min(s-1, n-1) # 비영 고유값 수
lambda_vals <- eigs[1:r]
Wilks_L <- prod(1 / (1 + lambda_vals))
# s=2이면 정확한 F 변환 가능
df1 <- n - 1; df2 <- N - s - n + 2
F_GT_multi <- (1 - Wilks_L) / Wilks_L * df2 / df1
p_GT_multi <- pf(F_GT_multi, df1, df2, lower.tail = FALSE)
cat("\n=== 다변량 집단×시간 교호작용 검정 (3.11) ===\n")
cat("고유값:", round(lambda_vals, 4), "\n")
cat(sprintf("Wilks' Λ = %.4f\n", Wilks_L))
cat(sprintf("F(%d, %d) = %.2f, p = %.4f\n", df1, df2, F_GT_multi, p_GT_multi))
# ──────────────────────────────────────────
# Step 6: 개별 추세 교호작용 F-검정 (3.14)
# ──────────────────────────────────────────
cat("\n=== 개별 추세 교호작용 F-검정 (MANOVA 방식 3.14) ===\n")
trends <- c("선형", "2차", "3차")
for (k in 1:(n-1)) {
ssg_k <- diag(SSG_star)[k+1]
ssr_k <- diag(SSR_star)[k+1]
F_k <- (ssg_k / (s-1)) / (ssr_k / (N-s))
p_k <- pf(F_k, s-1, N-s, lower.tail = FALSE)
cat(sprintf("집단×%s: ssg=%.2f, ssr=%.2f, F(1,%d)=%.2f, p=%.4f\n",
trends[k], ssg_k, ssr_k, N-s, F_k, p_k))
}
# ──────────────────────────────────────────
# Step 7: 내장 manova() 검증
# ──────────────────────────────────────────
cat("\n=== 내장 manova() 검증 ===\n")
fit <- manova(Y ~ group)
print(summary(fit, test = "Wilks"))
print(summary.aov(fit)) # 각 시점별 단변량 ANOVA9.2 Python: statsmodels MANOVA
import numpy as np
import pandas as pd
from scipy import linalg, stats
from statsmodels.multivariate.manova import MANOVA
np.random.seed(42)
N_g, n, s = 30, 4, 2
N = N_g * s
# 공분산 행렬
Sigma = np.array([
[16, 8, 7, 6],
[ 8, 14, 8, 7],
[ 7, 8, 12, 8],
[ 6, 7, 8, 11]
], dtype=float)
mu_drug = np.array([22, 16, 12, 10], dtype=float)
mu_placebo = np.array([22, 20, 19, 18], dtype=float)
L = linalg.cholesky(Sigma, lower=True)
rng = np.random.default_rng(42)
Y_drug = mu_drug + rng.standard_normal((N_g, n)) @ L.T
Y_placebo = mu_placebo + rng.standard_normal((N_g, n)) @ L.T
Y = np.vstack([Y_drug, Y_placebo])
# ──────────────────────────────────────────
# Step 1: 직교 다항식 행렬
# ──────────────────────────────────────────
times = np.arange(1, n + 1)
T_mat = np.array([times**d for d in range(n)])
S = linalg.cholesky(T_mat @ T_mat.T, lower=True)
P = linalg.solve_triangular(S, T_mat, lower=True)
# ──────────────────────────────────────────
# Step 2: 집단별·전체 평균 및 직교 다항 계수
# ──────────────────────────────────────────
y_bar_drug = Y_drug.mean(axis=0)
y_bar_placebo = Y_placebo.mean(axis=0)
y_bar_all = Y.mean(axis=0)
delta_drug = P @ y_bar_drug
delta_placebo = P @ y_bar_placebo
print("직교 다항 추정값 (약물):", np.round(delta_drug, 3))
print("직교 다항 추정값 (위약):", np.round(delta_placebo, 3))
# ──────────────────────────────────────────
# Step 3: 세 SSCP 행렬
# ──────────────────────────────────────────
SST_star = N * P @ np.outer(y_bar_all, y_bar_all) @ P.T
SSG_raw = (N_g * np.outer(y_bar_drug, y_bar_drug) +
N_g * np.outer(y_bar_placebo, y_bar_placebo) -
N * np.outer(y_bar_all, y_bar_all))
SSG_star = P @ SSG_raw @ P.T
SSY_raw = Y.T @ Y
SSR_star = P @ (SSY_raw - N_g * np.outer(y_bar_drug, y_bar_drug)
- N_g * np.outer(y_bar_placebo, y_bar_placebo)) @ P.T
print("\nSST* 대각:", np.round(np.diag(SST_star), 2))
print("SSG* 대각:", np.round(np.diag(SSG_star), 2))
print("SSR* 대각:", np.round(np.diag(SSR_star), 2))
# ──────────────────────────────────────────
# Step 4: 분할구획 ANOVA SS 추출
# ──────────────────────────────────────────
SS_T = np.sum(np.diag(SST_star)[1:])
SS_G = np.diag(SSG_star)[0]
SS_GT = np.sum(np.diag(SSG_star)[1:])
SS_SG = np.diag(SSR_star)[0]
SS_R = np.sum(np.diag(SSR_star)[1:])
df_T, df_G, df_GT = n-1, s-1, (s-1)*(n-1)
df_SG, df_R = N-s, (N-s)*(n-1)
MS_T = SS_T / df_T; MS_G = SS_G / df_G
MS_GT = SS_GT / df_GT; MS_SG = SS_SG / df_SG; MS_R = SS_R / df_R
F_T = MS_T / MS_R
F_G = MS_G / MS_SG # 집단 효과: 분모 = MS_S(G)
F_GT = MS_GT / MS_R
for label, SS, df, F, denom_df in [
("시간", SS_T, df_T, F_T, df_R),
("집단", SS_G, df_G, F_G, df_SG),
("집단×시간", SS_GT, df_GT, F_GT, df_R),
]:
p = 1 - stats.f.cdf(F, df, denom_df)
print(f"{label}: SS={SS:.2f}, df={df}, F={F:.2f}, p={p:.4f}")
# ──────────────────────────────────────────
# Step 5: 다변량 교호작용 검정
# ──────────────────────────────────────────
SSG_sub = SSG_star[1:, 1:]
SSR_sub = SSR_star[1:, 1:]
E = linalg.cholesky(SSR_sub, lower=True)
A_gT = linalg.solve_triangular(E, SSG_sub, lower=True)
A_gT = linalg.solve_triangular(E, A_gT.T, lower=True).T
eigs = np.sort(np.real(linalg.eigvals(A_gT)))[::-1]
r = min(s-1, n-1)
lam = eigs[:r]
Wilks_GT = np.prod(1 / (1 + lam))
df1 = n-1; df2 = N - s - n + 2
F_GT_m = (1 - Wilks_GT) / Wilks_GT * df2 / df1
p_GT_m = 1 - stats.f.cdf(F_GT_m, df1, df2)
print(f"\n다변량 교호작용: Λ={Wilks_GT:.4f}, F({df1},{df2})={F_GT_m:.2f}, p={p_GT_m:.4f}")
# ──────────────────────────────────────────
# Step 6: 개별 추세 교호작용 F-검정 (3.14)
# ──────────────────────────────────────────
print("\n개별 추세 교호작용 (3.14):")
for k, trend in enumerate(["선형", "2차", "3차"], 1):
ssg_k = np.diag(SSG_star)[k]
ssr_k = np.diag(SSR_star)[k]
F_k = (ssg_k / (s-1)) / (ssr_k / (N-s))
p_k = 1 - stats.f.cdf(F_k, s-1, N-s)
print(f" 집단×{trend}: F(1,{N-s})={F_k:.2f}, p={p_k:.4f}")
# ──────────────────────────────────────────
# Step 7: statsmodels로 검증
# ──────────────────────────────────────────
cols = ["w0", "w2", "w4", "w8"]
df_wide = pd.DataFrame(Y, columns=cols)
df_wide["group"] = ["drug"] * N_g + ["placebo"] * N_g
formula = "w0 + w2 + w4 + w8 ~ group"
maov = MANOVA.from_formula(formula, data=df_wide)
result = maov.mv_test()
print("\nstatsmodels MANOVA:")
print(result.summary())10 핵심 정리
| 항목 | 내용 |
|---|---|
| 모형 | \(\mathbf{y}_{hi} = \boldsymbol{\mu} + \boldsymbol{\gamma}_h + \boldsymbol{\varepsilon}_{hi}\), \(\boldsymbol{\gamma}_h\) 은 \(n \times 1\) 집단 효과 벡터 |
| 세 SSCP 행렬 | \(\mathbf{SST}^*\) (시간), \(\mathbf{SSG}^*\) (집단+교호작용), \(\mathbf{SSR}^*\) (피험자-내 잔차) |
| \(\mathbf{SSG}^*\) 분해 | \(ssg_0\) = 집단 주효과, \(ssg_k (k \geq 1)\) = 집단 × \(k\) 차 교호작용 |
| 집단 효과 분모 | \(MS_{S(G)}\) (피험자-간 오차) — 시간·교호작용의 분모 \(MS_R\) 과 다름 |
| 다변량 교호작용 | \(\|\mathbf{SSG}_{(n-1)}^* - \lambda\mathbf{SSR}_{(n-1)}^*\| = 0\) → \(r = \min(s-1, n-1)\) 고유값 |
| Wilks’ \(\Lambda\) | \(\prod_{k=1}^r 1/(1+\lambda_k)\) — 고유값이 여러 개일 수 있음 |
| 개별 추세 검정 | \(F_{GT_k} = [ssg_k/(s-1)] / [ssr_k/(N-s)]\) — 각자의 분모 사용 (3.14) |
| 교호작용 = 핵심 | 임상·사회과학에서 “약물이 효과 있는가”의 핵심 답은 교호작용에 있음 |
11 관련 주제
선행 지식
후속 주제
관련 개념