1 종단 연구의 네 가지 장점
Hedeker & Gibbons(2006, Ch.1)는 종단 연구가 횡단 연구보다 우월한 네 가지 이유를 제시한다.
1.1 장점 1: 높은 검정력 (Statistical Power)
반복 측정한 관측치는 완전 독립이 아니므로, 같은 사람에게서 얻은 \(n_i\)회 측정은 \(n_i\)명의 독립 측정보다 더 많은 독립적 정보를 제공한다.
\[ \text{Var}(d_i) = 2\sigma^2(1 - \rho) \]
여기서 \(\rho\)는 같은 사람의 두 시점 간 상관이다. \(\rho > 0\)이면 차분의 분산이 줄어들어 동일 표본 크기에서 더 높은 검정력을 얻는다.
직관: 두 시점 혈압 변화를 측정할 때, 같은 사람의 두 측정치가 비슷할수록( \(\rho\) 높을수록) 실제 변화 신호가 노이즈 대비 더 선명하게 드러난다.
1.2 장점 2: 피험자가 자기 자신의 대조군
각 피험자가 처치 전·후를 모두 경험하므로, 유전적 소인·환경·측정되지 않은 개인 특성이 자동으로 통제된다.
\[ y_{i2} - y_{i1} = \text{처치 효과} + \text{개인 내 노이즈} \]
개인 간 이질성( \(\sigma_b^2\))이 오차에서 제거되므로 추정이 훨씬 효율적이다.
1.3 장점 3: 나이 효과와 코호트 효과의 분리
\[ \text{관측 변화} = \underbrace{f(\text{Age})}_{\text{개인 내 변화}} + \underbrace{h(\text{Cohort})}_{\text{출생 집단 효과}} \]
횡단 연구는 두 효과를 구분할 수 없다. 종단 연구는 동일인을 추적하여 시간에 따른 나이 효과(within-person)와 코호트 효과(between-person)를 분리한다.
어느 시점에 20~80세 1,000명을 한 번 측정하면, 고령일수록 혈압이 높다. 그런데 이 차이가 “나이가 들면 혈압이 오른다”(나이 효과)인가, 아니면 “고령 세대가 젊을 때부터 고염식 식문화를 가졌다”(코호트 효과)인가?
종단 연구는 동일인을 20세부터 추적하여 두 효과를 분리한다.
1.4 장점 4: 개인 변화에 대한 정보
횡단 데이터는 “집단 평균”만 추정한다. 종단 데이터는 개인별 변화 궤적을 추정하여 집단 내 이질성과 성장·변화의 결정 요인을 이해할 수 있다.
| 연구 유형 | 답할 수 있는 질문 |
|---|---|
| 횡단 | “지금 이 순간 집단 평균이 얼마인가?” |
| 종단 | “시간에 따라 각 개인이 어떻게 변하는가?” |
| 종단 | “변화 속도가 개인마다 다른가? 그 차이의 결정 요인은?” |
2 종단 데이터 분석의 도전
2.1 도전 1: 상관된 관측치
같은 피험자의 반복 측정치는 독립이 아니므로, 독립성을 가정하는 표준 회귀·ANOVA를 그대로 적용할 수 없다. 더 정교한 통계 방법(혼합효과 모형, GEE 등)이 필요하다.
2.2 도전 2: 결측치와 탈락(Dropout)
종단 연구에서 결측치는 세 가지 메커니즘을 갖는다:
| 메커니즘 | 정의 | 결과 |
|---|---|---|
| MCAR (완전 임의 결측) | 결측이 관측·미관측 변수 모두와 무관 | 단순 분석도 일치 추정 |
| MAR (임의 결측) | 관측된 변수로 결측 예측 가능 | 우도 기반 분석(LMM) 일치 |
| MNAR (비임의 결측) | 미관측 결과와 결측이 연관 | 모든 표준 방법이 편의 |
결측 처리의 흔한 오류 두 가지 (Hedeker & Gibbons, 2006, Ch.1):
- 완성자 분석(Completer analysis): 모든 시점을 완료한 피험자만 분석. 탈락 이유가 처치 효과와 혼재되면 편향 발생.
- LOCF (Last Observation Carried Forward): 마지막 관측값을 이후 결측에 대체. 탈락 전 노출 정도가 다른 피험자들을 동일하게 취급하는 문제.
혼합효과 회귀 모형은 MAR 가정 하에서 가용한 모든 데이터를 사용한다. 불완전 데이터를 가진 피험자도 분석에 포함되므로, 임의적 제외로 인한 편향이 줄어들고 검정력이 높아진다 (Hedeker & Gibbons, 2006, Ch.1).
2.3 도전 3: 시변 공변량 (Time-Varying Covariates)
결과변수뿐 아니라 예측변수도 시간에 따라 변할 수 있다. 예: 혈중 약물 농도(예측변수)와 건강 상태(결과변수)가 모두 시간에 따라 변한다.
이때 목표는 개인 내에서 두 변수의 동적 관계를 추정하는 것이다. 횡단 데이터는 이런 동적 관계를 추정할 수 없다.
2.4 도전 4: 이월 효과 (Carry-Over Effects)
교차 설계(crossover design)에서 피험자가 여러 처치를 순서대로 받을 때, 이전 처치의 효과가 다음 처치 기간까지 지속될 수 있다. 이월 효과를 처리하는 것은 통계적으로 복잡하며, 종단 연구의 가격(price)이다.
3 표기법과 데이터 배치
3.1 기본 표기법 (Hedeker & Gibbons, 2006, Ch.1.3)
\[ i = 1, \ldots, N \quad \text{(피험자 인덱스)} \]
균형 설계(balanced): \[ j = 1, \ldots, n \quad \text{(동일 측정 횟수)} \]
불균형 설계(unbalanced): \[ j = 1, \ldots, n_i \quad \text{(피험자마다 다른 측정 횟수)} \]
총 관측치 수: \[ \text{Total observations} = \sum_{i=1}^{N} n_i \]
피험자 \(i\)의 반응 벡터: \[ \mathbf{y}_i = (y_{i1}, y_{i2}, \ldots, y_{in_i})^\top \in \mathbb{R}^{n_i} \]
피험자 \(i\), 시점 \(j\)의 공변량 벡터 (\(p\)개, 절편 포함): \[ \mathbf{x}_{ij} = p \times 1 \]
시간 불변 공변량 (between-subject, 예: 성별): \(j = 1, \ldots, n_i\) 동안 일정 시변 공변량 (within-subject, 예: 나이): 피험자·시점별로 값이 다름
피험자 \(i\)의 공변량 행렬: \[ \mathbf{X}_i = n_i \times p \]
3.2 데이터 배치 (Long Format) (Hedeker & Gibbons, 2006, Ch.1.4)
종단 데이터는 long format (univariate layout)으로 표현한다: 측정값 하나가 한 행을 차지하고, 피험자당 \(n_i\)개의 행이 있다.
| Subject | Observation | Response | Covariate 1 | \(\cdots\) | Covariate \(p\) |
|---|---|---|---|---|---|
| 1 | 1 | \(y_{11}\) | \(x_{111}\) | \(\cdots\) | \(x_{11p}\) |
| 1 | 2 | \(y_{12}\) | \(x_{121}\) | \(\cdots\) | \(x_{12p}\) |
| \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\ddots\) | \(\vdots\) |
| 1 | \(n_1\) | \(y_{1n_1}\) | \(x_{1n_1 1}\) | \(\cdots\) | \(x_{1n_1 p}\) |
| \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\ddots\) | \(\vdots\) |
| \(N\) | 1 | \(y_{N1}\) | \(x_{N11}\) | \(\cdots\) | \(x_{N1p}\) |
| \(N\) | \(n_N\) | \(y_{Nn_N}\) | \(x_{Nn_N 1}\) | \(\cdots\) | \(x_{Nn_N p}\) |
이 구조는 다층 모형(multilevel model) 문헌에서 2-수준 설계라 부른다: - Level 1: 반복 관측치 (observations nested within subjects) - Level 2: 피험자 (subjects)
피험자가 다시 기관(병원, 학교 등) 안에 내포되면 3-수준 설계가 된다.
Wide format과의 비교:
| 형식 | 특징 | 용도 |
|---|---|---|
| Long (univariate) | 측정값 하나 = 한 행, 피험자당 \(n_i\)행 | LMM, GEE, 혼합효과 모형 |
| Wide (multivariate) | 피험자 하나 = 한 행, 시점별 변수가 열로 분리 | MANOVA, 전통적 반복측정 ANOVA |
4 핵심 문제: 독립성 가정의 위반과 OLS의 실패
4.1 왜 일반 OLS가 틀리는가
표준 선형 회귀모형은 다음을 가정한다:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}_N) \]
실제 종단 데이터의 오차 공분산 구조:
\[ \text{Cov}(\boldsymbol{\epsilon}) = \boldsymbol{\Sigma} = \text{diag}(\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2, \ldots, \boldsymbol{\Sigma}_N) \]
여기서 \(\boldsymbol{\Sigma}_i \neq \sigma^2 \mathbf{I}_{n_i}\) — 같은 피험자 내 관측치 간 상관이 존재한다.
4.2 OLS 표준오차의 편의
OLS 추정량 \(\hat{\boldsymbol{\beta}}_{OLS} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}\)는 불편이지만, 실제 분산:
\[ \text{Var}(\hat{\boldsymbol{\beta}}_{OLS}) = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \boldsymbol{\Sigma} \mathbf{X} (\mathbf{X}^\top \mathbf{X})^{-1} \]
OLS가 가정하는 분산 \(\hat{\sigma}^2 (\mathbf{X}^\top \mathbf{X})^{-1}\)와 다르므로, 표준오차가 과소추정된다.
| 문제 | 수식/설명 | 실무적 결과 |
|---|---|---|
| 편의된 SE | \(\widehat{\text{SE}}_{OLS} < \text{SE}_{true}\) | 좁은 신뢰구간, 거짓 유의성 (1종 오류 팽창) |
| 비효율적 추정 | \(\text{Var}(\hat{\boldsymbol{\beta}}_{OLS}) > \text{Var}(\hat{\boldsymbol{\beta}}_{GLS})\) | 불필요하게 큰 표본 요구 |
| 구조 누락 | \(\boldsymbol{\Sigma}_i\) 미모형화 | 개인 궤적, 시간 트렌드 파악 불가 |
직관: 30명을 4회 측정했는데 같은 사람의 측정치가 모두 비슷하다면, 실질적인 독립 정보는 120개가 아니라 30명 분량에 가깝다. OLS는 120개의 독립 관측치인 것처럼 취급하므로 표준오차를 과소추정한다.
4.3 올바른 추정량: GLS
일반화 최소제곱(GLS):
\[ \hat{\boldsymbol{\beta}}_{GLS} = (\mathbf{X}^\top \boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \mathbf{X}^\top \boldsymbol{\Sigma}^{-1} \mathbf{y}, \quad \text{Var}(\hat{\boldsymbol{\beta}}_{GLS}) = (\mathbf{X}^\top \boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \]
\(\boldsymbol{\Sigma}\)를 모르면 데이터로 추정해야 하며, 이것이 LMM·GEE의 핵심 작업이다.
5 변동의 분해와 ICC
5.1 ANOVA 기반 변동 분해
랜덤 절편 모형 \(y_{ij} = \mu + b_i + \epsilon_{ij}\) (\(b_i \sim N(0, \sigma_b^2)\), \(\epsilon_{ij} \sim N(0, \sigma_e^2)\), 독립)에서:
\[ \text{Var}(y_{ij}) = \sigma_b^2 + \sigma_e^2 \]
균형 설계 (\(n_i = T\))에서 일원배치 ANOVA를 통한 분해:
\[ SS_{Total} = SS_{Between} + SS_{Within} \]
\[ SS_{Between} = T\sum_{i=1}^N (\bar{y}_{i\cdot} - \bar{y}_{\cdot\cdot})^2, \quad SS_{Within} = \sum_{i=1}^N \sum_{j=1}^T (y_{ij} - \bar{y}_{i\cdot})^2 \]
평균 제곱의 기댓값:
\[ E[MS_{Between}] = \sigma_e^2 + T\sigma_b^2, \qquad E[MS_{Within}] = \sigma_e^2 \]
모적률 추정:
\[ \hat{\sigma}_e^2 = MS_{Within}, \qquad \hat{\sigma}_b^2 = \frac{MS_{Between} - MS_{Within}}{T} \]
직관: \(MS_{Between} \gg MS_{Within}\)이면 개인들의 기저 수준 차이가 크다. \(MS_{Between} \approx MS_{Within}\)이면 개인 간 차이보다 시점 내 변동이 지배한다.
5.2 급내상관계수 (ICC)
\[ \hat{\rho} = \frac{\hat{\sigma}_b^2}{\hat{\sigma}_b^2 + \hat{\sigma}_e^2} = \frac{MS_{Between} - MS_{Within}}{MS_{Between} + (T-1) \cdot MS_{Within}} \]
ICC는 또한 같은 피험자의 임의의 두 관측치 간 상관계수:
\[ \text{Corr}(y_{ij}, y_{ij'}) = \frac{\sigma_b^2}{\sigma_b^2 + \sigma_e^2} = \rho, \quad j \neq j' \]
| ICC | 해석 | 반복 측정 효율 |
|---|---|---|
| \(< 0.1\) | 개인 내 변동 지배 | 이득 거의 없음 |
| \(0.1 \sim 0.4\) | 중간 | 어느 정도 효율적 |
| \(0.4 \sim 0.7\) | 높음 | 효율적 |
| \(\geq 0.7\) | 매우 높음 | 큰 검정력 이득 |
\(\rho = 0.7\)이면 같은 사람의 관측치들은 서로 70%를 공유하고, 30%만 독립적인 정보다.
6 공분산 구조의 종류
\[ \boldsymbol{\Sigma}_i = \text{Cov}(\mathbf{y}_i) \in \mathbb{R}^{n_i \times n_i} \]
6.1 등분산 구조
1. 독립 (Independence) 모수 수: 1, OLS 가정, 종단 데이터에 부적합
\[ \boldsymbol{\Sigma} = \sigma^2\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \]
2. 복합 대칭 (Compound Symmetry, CS) 모수 수: 2, 반복측정 ANOVA의 구형성 가정과 동치
\[ \boldsymbol{\Sigma} = \sigma^2[(1-\rho)\mathbf{I} + \rho\mathbf{1}\mathbf{1}^\top] = \begin{pmatrix} \sigma^2 & \rho\sigma^2 & \rho\sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \sigma^2 & \rho\sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \rho\sigma^2 & \sigma^2 & \rho\sigma^2 \\ \rho\sigma^2 & \rho\sigma^2 & \rho\sigma^2 & \sigma^2 \end{pmatrix} \]
- 직관: “가까운 시점이든 먼 시점이든 같은 사람이면 똑같이 닮아 있다”
3. 1차 자기회귀 (AR(1)) 모수 수: 2, 시간이 지날수록 상관이 지수적으로 감소
\[ [\boldsymbol{\Sigma}]_{jk} = \sigma^2\rho^{|j-k|}, \quad \boldsymbol{\Sigma} = \sigma^2\begin{pmatrix} 1 & \rho & \rho^2 & \rho^3 \\ \rho & 1 & \rho & \rho^2 \\ \rho^2 & \rho & 1 & \rho \\ \rho^3 & \rho^2 & \rho & 1 \end{pmatrix} \]
- 직관: “지난주 측정이 오늘에 가장 영향 미치고, 먼 시점일수록 기하급수적으로 감소”
4. Toeplitz (TOEP) 모수 수: \(T\), AR(1)보다 유연 — 같은 간격이면 같은 상관, 단조성 비강제
\[ \boldsymbol{\Sigma} = \sigma^2\begin{pmatrix} 1 & \rho_1 & \rho_2 & \rho_3 \\ \rho_1 & 1 & \rho_1 & \rho_2 \\ \rho_2 & \rho_1 & 1 & \rho_1 \\ \rho_3 & \rho_2 & \rho_1 & 1 \end{pmatrix} \]
5. 비구조화 (Unstructured, UN) 모수 수: \(T(T+1)/2\), 가장 유연하나 소표본에서 불안정
\[ \boldsymbol{\Sigma} = \begin{pmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} & \sigma_{14} \\ \sigma_{12} & \sigma_{22} & \sigma_{23} & \sigma_{24} \\ \sigma_{13} & \sigma_{23} & \sigma_{33} & \sigma_{34} \\ \sigma_{14} & \sigma_{24} & \sigma_{34} & \sigma_{44} \end{pmatrix} \]
6.2 이분산 구조
6. 이분산 복합 대칭 (HCS) \([\boldsymbol{\Sigma}]_{jk} = \rho\sigma_j\sigma_k\;(j \neq k)\), \([\boldsymbol{\Sigma}]_{jj} = \sigma_j^2\) 모수 수: \(T+1\) — 각 시점의 분산이 다르지만 상관은 동일
7. 이분산 AR(1) (HAR(1)) \([\boldsymbol{\Sigma}]_{jk} = \rho^{|j-k|}\sigma_j\sigma_k\) 모수 수: \(T+1\) — 이분산 + AR(1) 상관
6.3 공분산 구조 선택 기준
| 구조 | 모수 수 | 적합 상황 |
|---|---|---|
| 독립 | 1 | 횡단 자료 (종단에 사용 금지) |
| CS | 2 | 구형성 가정 타당, 시점 간격 무관 |
| AR(1) | 2 | 등간격 시계열 측정 |
| Toeplitz | T | AR(1)보다 유연, 단조성 비강제 |
| HCS | T+1 | 시점별 분산 이질적 |
| 비구조화 | T(T+1)/2 | 표본 충분, 구조 가정 회피 |
REML 기반 AIC/BIC로 비교하고, 이론적 근거로 후보를 먼저 결정한다.
7 분석 방법 선택을 위한 고려사항
Hedeker & Gibbons(2006, Ch.1.5)는 분석 방법 선택 시 다섯 가지 요소를 고려해야 한다고 정리한다:
1. 결과변수의 형태
| 결과 형태 | 권장 분석 방법 |
|---|---|
| 연속형, 정규 분포 | 혼합효과 선형 회귀 (LMM) |
| 연속형, 비정규 (카운트 등) | 혼합효과 포아송/음이항 회귀 |
| 이진형 (yes/no) | 혼합효과 로지스틱/프로빗 회귀 |
| 순서형 | 혼합효과 비례 오즈 회귀 |
| 명목형 | 혼합효과 다항 회귀 |
2. 피험자 수 \(N\)
- 고급 혼합효과 모형은 대표본 이론에 기반 → \(N < 50\)이면 부적합 가능
- GEE의 샌드위치 추정량도 \(N\)이 작으면 하향 편의
3. 피험자당 관측 횟수 \(n_i\)
- \(n_i = 2\) (전·후 두 시점): 변화점수 분석 또는 ANCOVA로 충분
- \(n_i = n\) (균형): 반복측정 ANOVA/MANOVA 가능
- \(n_i\)가 피험자마다 다름: 혼합효과 회귀 모형이 필수
4. 공변량의 종류
- 시간 불변: 처치군·대조군, 성별, 인종 등
- 시변: 나이, 약물 농도, 컨디션 등 → 모형이 시변 공변량을 처리할 수 있어야 함
5. 분산-공분산 구조
- 동분산/이분산 여부
- 시간에 따른 공분산의 이질성(proximal occasions vs distal occasions)
- 잔차 자기상관
8 일반 분석 방법 분류
Hedeker & Gibbons(2006, Ch.1.6)는 종단 데이터 분석의 여섯 가지 일반 접근법을 정리한다:
8.1 접근 1: 유도 변수 방법 (Derived Variable Approach)
반복 측정을 하나의 요약 변수로 축소한 뒤 분석한다.
| 유도 변수 예 | 정의 |
|---|---|
| 시간 평균 | \(\bar{y}_i = n_i^{-1}\sum_j y_{ij}\) |
| 변화점수 | \(d_i = y_{i,n_i} - y_{i1}\) |
| 선형 추세 | 개인별 OLS 기울기 |
| AUC | 곡선 아래 면적 |
| LOCF | 마지막 관측값 이월 |
한계: 불균형 데이터에서 \(n_i\)가 다르면 유도 변수의 불확실성이 다르다 → 이분산성 위반. 시변 공변량 사용 불가. 검정력 손실.
8.2 접근 2: 반복측정 ANOVA
복합 대칭(CS) 가정 하에 시간 효과를 F-검정으로 평가한다. 한계: 결측 피험자를 완전 제외. 시변 공변량 불가. CS 가정이 실제로 성립하기 어렵다.
8.3 접근 3: MANOVA
Wide format의 반응 벡터를 다변량 반응으로 처리한다. 장점: CS 불필요, 비구조화 공분산 허용. 한계: 완전 데이터만 사용. 시변 공변량 불가. 시점 수가 많으면 공분산 모수 폭발.
8.4 접근 4: 혼합효과 회귀 모형 (Mixed-Effects Regression Models)
이 책의 주요 강조점이다. 정규 연속형부터 이진·순서·명목·카운트 결과까지 적용 가능하다. MAR 결측 하에서도 가용한 모든 데이터를 활용한다. 불균형 설계, 시변 공변량 처리 가능.
8.5 접근 5: 공분산 패턴 모형 (Covariance Pattern Models, CPM)
주변(marginal) 분포를 직접 모형화한다. 공분산 구조를 소수의 모수로 명세한다. 랜덤 효과를 따로 구분하지 않는 점에서 LMM과 다르다. 개인별 예측 불가.
8.6 접근 6: GEE (Generalized Estimating Equations)
준우도(quasi-likelihood) 기반의 주변 모형이다. 작동 상관이 잘못 지정되어도 \(\boldsymbol{\beta}\) 추정의 일치성을 보장한다. 한계: MAR이 아닌 MCAR 가정 필요. 소표본에서 샌드위치 SE 하향 편의.
9 가장 단순한 종단 분석
Hedeker & Gibbons(2006, Ch.1.7)는 2시점 설계에서의 단순 분석을 세 가지로 정리한다.
9.1 기본 설정
\(N\)명의 피험자, 사전(\(y_{i1}\))·사후(\(y_{i2}\)) 측정, 변화점수 \(d_i = y_{i2} - y_{i1}\).
전체 변화가 유의한지 검정: \(H_0: \mu_{y_1} = \mu_{y_2}\), 즉 \(H_0: \mu_{y_2} - \mu_{y_1} = 0\)
\[ t = \frac{\bar{d}}{s_d / \sqrt{N}} \overset{H_0}{\sim} t_{N-1} \]
동등하게 회귀 모형으로 표현:
\[ d_i = \beta_0 + e_i, \quad H_0: \beta_0 = 0 \text{ 검정} \]
9.2 1.7.1 변화점수 분석 (Change Score Analysis)
두 집단(처치군 \(x_i = 1\), 대조군 \(x_i = 0\))에서 변화량의 집단 간 차이를 검정:
\[ d_i = \beta_0 + \beta_1 x_i + e_i \]
- \(\beta_0\): 대조군의 평균 변화량
- \(\beta_1\): 두 집단의 평균 변화량 차이
가설 검정: - \(H_0: \beta_0 = 0\): 대조군의 평균 변화가 0인가? - \(H_0: \beta_1 = 0\): 두 집단의 평균 변화가 동일한가?
변화점수 분석은 사전 점수를 기울기 1로 고정한 ANCOVA와 동치임을 확인할 수 있다:
\[ d_i = \beta_0 + \beta_1 x_i + e_i \;\Leftrightarrow\; y_{i2} = y_{i1} + \beta_0 + \beta_1 x_i + e_i \]
그러나 사전-사후 관계의 기울기가 반드시 1일 필요는 없다 — 이것이 ANCOVA가 더 일반적인 이유다.
9.3 1.7.2 사후 점수의 ANCOVA (Analysis of Covariance of Post-test Scores)
사전 점수를 공변량으로 포함하여 사후 점수를 분석:
\[ y_{i2} = \beta_0 + \beta_1 x_i + \beta_2 y_{i1} + e_i \]
- \(\beta_2 \neq 1\): 사전-사후 관계가 일대일이 아닐 때 더 유연
- \(\beta_1\) 검정: 동일 사전 점수 조건에서 집단 간 사후 점수 차이 (조건부 비교)
- 무작위 배정 임상시험에서 변화점수 분석보다 더 효율적(검정력 높음)
비무작위 설계에서 집단 간 사전 점수가 다를 때:
- 변화점수 분석과 ANCOVA가 다른 결론을 낼 수 있다
- 두 방법이 서로 다른 질문에 답하기 때문:
- 변화점수: “두 집단의 평균 변화가 같은가?”
- ANCOVA: “같은 사전 점수를 가진 집단에서 사후 점수가 같은가?”
- 무작위 배정이면 사전 점수 분포가 동일 → 두 방법이 유사한 결론
9.4 1.7.3 변화점수의 ANCOVA (ANCOVA of Change Scores)
변화점수를 결과로, 사전 점수를 공변량으로:
\[ d_i = \beta_0 + \beta_1 x_i + \beta_2 y_{i1} + e_i \]
이를 전개하면:
\[ y_{i2} - y_{i1} = \beta_0 + \beta_1 x_i + \beta_2 y_{i1} + e_i \] \[ y_{i2} = \beta_0 + \beta_1 x_i + (1 + \beta_2) y_{i1} + e_i \]
이것은 사후 점수 ANCOVA (\(\beta_1\) 계수는 동일, \(y_{i1}\) 계수만 다름)와 처치 효과( \(\beta_1\) )에 대해 완전히 동일하다.
결론: 변화점수에 대한 ANCOVA와 사후 점수에 대한 ANCOVA는 처치 효과 검정에서 같은 결과를 낸다.
9.5 1.7.4 예시: THKS 흡연 예방 연구 (Hedeker & Gibbons, 2006, Ch.1.7.4)
Flay et al.(1988)의 Television School and Family Smoking Prevention and Cessation Project:
- 표본: 1,600명의 7학년 학생, 135개 학급, 28개 LA 학교
- 결과변수: 흡연 관련 지식 (THKS 척도)
- 시점: 중재 전(pre) · 중재 후(post)
- 설계: 학교를 네 집단에 무작위 배정
- CC (교실 커리큘럼) 단독
- TV (미디어 개입) 단독
- CC + TV 병합
- 대조군 (무처치)
전체 변화: 사전 평균 2.069 → 사후 평균 2.662 (증가량 0.59)
짝 \(t\)-검정: \(t_{1599} = 15.01\), \(p < .0001\) — 유의미한 지식 향상.
회귀 모형으로도 동일: \(\hat{\beta}_0 = 0.5925\), SE = 0.0395, \(t_{1599} = 15.01\).
집단별 기술 통계:
| 집단 | \(N\) | 사전 평균 | 사후 평균 | 차이 |
|---|---|---|---|---|
| CC=No, TV=No (대조) | 421 | 2.152 | 2.361 | 0.209 |
| CC=No, TV=Yes | 416 | 2.087 | 2.539 | 0.452 |
| CC=Yes, TV=No | 380 | 2.050 | 2.968 | 0.918 |
| CC=Yes, TV=Yes | 383 | 1.979 | 2.823 | 0.844 |
사후 점수 회귀 분석:
\[ \hat{y}_{i2} = 2.361 + 0.607\cdot CC + 0.177\cdot TV - 0.323\cdot CC \times TV \]
| 변수 | 추정값 | SE | \(t\) | \(p\) |
|---|---|---|---|---|
| 절편 | 2.3611 | 0.0665 | 35.52 | \(<.0001\) |
| CC | 0.6074 | 0.0965 | 6.29 | \(<.0001\) |
| TV | 0.1774 | 0.0943 | 1.88 | 0.0600 |
| CC × TV | −0.3234 | 0.1365 | −2.37 | 0.0180 |
사전 점수 조정 후 ANCOVA:
\[ \hat{y}_{i2} = 1.661 + 0.325\cdot y_{i1} + 0.641\cdot CC + 0.199\cdot TV - 0.322\cdot CC \times TV \]
ANCOVA 후 TV 효과가 유의해짐(\(p = .0273\)) — 사전 점수를 통제하면 검정력 향상.
해석: CC와 TV 중재가 각각 지식 향상에 기여하지만, 두 효과는 비가산적이다(음의 교호작용). CC가 없을 때 TV 효과가 더 크다.
10 LMM에서 \(\mathbf{y}_i\)의 주변 분포
10.1 주변 분포 유도
랜덤 절편 모형에서 \(b_i\)를 적분하면:
\[ \mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \mathbf{1}_{n_i} b_i + \boldsymbol{\epsilon}_i \]
\[ E[\mathbf{y}_i] = \mathbf{X}_i\boldsymbol{\beta}, \quad \text{Cov}(\mathbf{y}_i) = \sigma_b^2\mathbf{J}_{n_i} + \sigma_e^2\mathbf{I}_{n_i} \]
이 주변 공분산이 복합 대칭 구조와 정확히 일치한다.
\[ \mathbf{y}_i \sim N(\mathbf{X}_i\boldsymbol{\beta},\; \sigma_b^2\mathbf{J}_{n_i} + \sigma_e^2\mathbf{I}_{n_i}) \]
일반 LMM (\(\mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \mathbf{Z}_i\mathbf{b}_i + \boldsymbol{\epsilon}_i\), \(\mathbf{b}_i \sim N(\mathbf{0}, \mathbf{D})\)):
\[ \mathbf{y}_i \sim N(\mathbf{X}_i\boldsymbol{\beta},\; \mathbf{Z}_i\mathbf{D}\mathbf{Z}_i^\top + \sigma^2\mathbf{I}_{n_i}) \]
| \(\mathbf{Z}_i\) | 공분산 구조 |
|---|---|
| \(\mathbf{1}_{n_i}\) (랜덤 절편만) | CS |
| \((\mathbf{1}_{n_i}, \mathbf{t}_i)\) (랜덤 절편+기울기) | 시간 의존적 복잡 구조 |
직관: LMM의 랜덤 효과 구조가 자동으로 공분산 구조를 생성한다. \(\mathbf{Z}_i\)의 선택이 공분산 구조를 결정한다.
11 시간 효과의 분리
11.1 Within-Between 분해 (Mundlak 분해)
시변 공변량 \(x_{ij}\)를 두 성분으로 분리:
\[ x_{ij} = \underbrace{\bar{x}_{i\cdot}}_{\text{between}} + \underbrace{(x_{ij} - \bar{x}_{i\cdot})}_{\text{within}} \]
이를 모형에 명시하면:
\[ y_{ij} = \alpha + \beta_W(x_{ij} - \bar{x}_{i\cdot}) + \beta_B\bar{x}_{i\cdot} + b_i + \epsilon_{ij} \]
- \(\beta_W\): 동일 개인의 \(x\)가 평균보다 높을 때 \(y\)의 변화 (인과에 가까운 within 효과)
- \(\beta_B\): 평균 \(x\)가 높은 사람의 \(y\) 수준 차이 (between 효과)
\(\beta_W \neq \beta_B\)이면 내생성이 있음을 시사한다.
11.2 Hausman 검정
LMM(랜덤 효과)과 고정 효과 모형 중 선택:
\[ H = (\hat{\boldsymbol{\beta}}_{FE} - \hat{\boldsymbol{\beta}}_{RE})^\top [\widehat{\text{Var}}(\hat{\boldsymbol{\beta}}_{FE}) - \widehat{\text{Var}}(\hat{\boldsymbol{\beta}}_{RE})]^{-1} (\hat{\boldsymbol{\beta}}_{FE} - \hat{\boldsymbol{\beta}}_{RE}) \overset{H_0}{\sim} \chi^2_{df} \]
- 유의: 내생성 → 고정 효과 모형 선택
- 비유의: 랜덤 효과 모형이 더 효율적
12 검정력 이점과 표본 크기
12.1 표본 크기 비교
독립 두 표본:
\[ n_{indep} = \frac{2(z_{\alpha/2} + z_\beta)^2\sigma^2}{\delta^2} \]
반복 측정 (두 시점):
\[ n_{repeated} = \frac{(z_{\alpha/2} + z_\beta)^2 \cdot 2\sigma^2(1-\rho)}{\delta^2} \]
\[ \frac{n_{repeated}}{n_{indep}} = 1 - \rho \]
| \(\rho\) | 독립 설계 대비 표본 비율 | 절감 |
|---|---|---|
| 0.3 | 0.70 | 30% |
| 0.5 | 0.50 | 50% |
| 0.7 | 0.30 | 70% |
| 0.9 | 0.10 | 90% |
실제 설계에서는 예상 탈락률을 보정한다:
\[ n_{adjusted} = \frac{n_{repeated}}{1 - \text{dropout rate}} \]
13 주요 분석 방법 비교
13.1 방법별 핵심 수식
선형 혼합 모형 (LMM):
\[ \mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \mathbf{Z}_i\mathbf{b}_i + \boldsymbol{\epsilon}_i, \quad \mathbf{b}_i \sim N(\mathbf{0}, \mathbf{D}), \quad \boldsymbol{\epsilon}_i \sim N(\mathbf{0}, \sigma^2\mathbf{I}) \]
GEE:
\[ \sum_{i=1}^N \mathbf{D}_i^\top\mathbf{V}_i^{-1}(\mathbf{y}_i - \boldsymbol{\mu}_i) = \mathbf{0} \]
Sandwich 분산: \(\mathbf{M}_0^{-1}\left(\sum_i \mathbf{D}_i^\top\mathbf{V}_i^{-1}\hat{\boldsymbol{\epsilon}}_i\hat{\boldsymbol{\epsilon}}_i^\top\mathbf{V}_i^{-1}\mathbf{D}_i\right)\mathbf{M}_0^{-1}\)
고정 효과 모형 (within 변환):
\[ \tilde{y}_{ij} = \tilde{\mathbf{x}}_{ij}^\top\boldsymbol{\beta} + \tilde{\epsilon}_{ij}, \quad \tilde{y}_{ij} = y_{ij} - \bar{y}_{i\cdot} \]
13.2 방법 선택 가이드
| 기준 | LMM | GEE | 고정 효과 |
|---|---|---|---|
| 효과 해석 | Subject-specific | Population-average | Within-person |
| 정규성 가정 | 필요 | 불필요 | 불필요 |
| 결측 처리 | MAR | MAR + Robust SE | 시불변 변수 불가 |
| 개인 예측 (BLUP) | 가능 | 불가 | 불가 |
| 내생성 통제 | 불완전 | 불완전 | 강함 |
| 비선형 확장 | GLMM | 자연스러움 | 어려움 |
14 코드 예시
14.1 Step 1: 데이터 생성 및 ANOVA 기반 ICC 추정
import numpy as np
import pandas as pd
np.random.seed(42)
n_subjects = 30
n_times = 4
sigma_b, sigma_e = 2.0, 1.0
beta_0, beta_1 = 5.0, 0.8
b_i = np.random.normal(0, sigma_b, n_subjects)
records = [
{'subject': i, 'time': t,
'y': beta_0 + b_i[i] + beta_1 * t + np.random.normal(0, sigma_e)}
for i in range(n_subjects) for t in range(n_times)
]
df = pd.DataFrame(records)
# ANOVA 기반 ICC
grand_mean = df['y'].mean()
subj_means = df.groupby('subject')['y'].mean()
ss_between = n_times * ((subj_means - grand_mean) ** 2).sum()
ss_within = df.groupby('subject')['y'].apply(
lambda g: ((g - g.mean()) ** 2).sum()).sum()
ms_between = ss_between / (n_subjects - 1)
ms_within = ss_within / (n_subjects * (n_times - 1))
sigma2_b_hat = (ms_between - ms_within) / n_times
sigma2_e_hat = ms_within
icc = sigma2_b_hat / (sigma2_b_hat + sigma2_e_hat)
print(f"E[MS_Between] = σ_e² + T·σ_b² = {sigma_e**2 + n_times*sigma_b**2:.3f}, 추정 = {ms_between:.3f}")
print(f"E[MS_Within] = σ_e² = {sigma_e**2:.3f}, 추정 = {ms_within:.3f}")
print(f"ICC 추정 = {icc:.3f} (참: {sigma_b**2/(sigma_b**2+sigma_e**2):.3f})")14.2 Step 2: 변화점수 분석 및 ANCOVA 비교
import statsmodels.formula.api as smf
# 2시점 데이터: 처치군 vs 대조군
np.random.seed(0)
n = 50
group = np.repeat([0, 1], n) # 0=대조군, 1=처치군
y_pre = np.random.normal(10, 2, 2*n)
effect = np.where(group == 1, 2.0, 0.0)
y_post = y_pre + effect + np.random.normal(0, 1.5, 2*n)
df2 = pd.DataFrame({'y_pre': y_pre, 'y_post': y_post,
'diff': y_post - y_pre, 'group': group})
# 변화점수 분석
res_change = smf.ols("diff ~ group", df2).fit()
print("=== 변화점수 분석 ===")
print(res_change.summary().tables[1])
# 사후 점수 ANCOVA (사전 점수 조정)
res_ancova = smf.ols("y_post ~ group + y_pre", df2).fit()
print("\n=== ANCOVA (사후 점수 + 사전 점수 조정) ===")
print(res_ancova.summary().tables[1])
# group 계수: 두 분석에서 거의 동일14.3 Step 3: LMM 적합 (statsmodels)
import statsmodels.formula.api as smf
# 랜덤 절편 모형
model_ri = smf.mixedlm("y ~ time", df, groups=df["subject"])
result_ri = model_ri.fit(method="lbfgs")
print(result_ri.summary())
# Fixed: time ≈ beta_1, Random Group Var ≈ sigma_b^2, Residual ≈ sigma_e^2
# 랜덤 절편 + 기울기
model_rs = smf.mixedlm("y ~ time", df, groups=df["subject"], re_formula="~time")
result_rs = model_rs.fit(method="lbfgs")
# LRT: 랜덤 기울기 필요성 검정
from scipy.stats import chi2
lr_stat = -2 * (result_ri.llf - result_rs.llf)
p_val = chi2.sf(lr_stat, df=2)
print(f"LRT χ² = {lr_stat:.3f}, df=2, p = {p_val:.4f}")14.4 Step 4: R 코드 (lme4, geepack, plm)
library(lme4)
library(geepack)
library(plm)
set.seed(42)
n <- 30; T <- 4
df <- expand.grid(subject = factor(1:n), time = 0:(T-1))
b_i <- rnorm(n, 0, 2)
df$y <- 5 + b_i[as.integer(df$subject)] + 0.8*df$time + rnorm(n*T, 0, 1)
# ── LMM ──
fit_ri <- lmer(y ~ time + (1 | subject), data = df, REML = TRUE)
fit_rs <- lmer(y ~ time + (time | subject), data = df, REML = TRUE)
anova(fit_ri, fit_rs) # LRT: 랜덤 기울기 필요성
# ICC 추출
var_comp <- as.data.frame(VarCorr(fit_ri))
icc_lmer <- var_comp$vcov[1] / sum(var_comp$vcov)
cat("ICC from LMM:", round(icc_lmer, 3), "\n")
# ── GEE ──
fit_gee <- geeglm(y ~ time, data = df, id = subject, corstr = "exchangeable")
summary(fit_gee)
# ── Hausman Test: 고정 vs 랜덤 효과 모형 ──
pdata <- pdata.frame(df, index = c("subject", "time"))
fit_fe <- plm(y ~ time, data = pdata, model = "within")
fit_re <- plm(y ~ time, data = pdata, model = "random")
phtest(fit_fe, fit_re) # p<0.05: 고정 효과 선택, p≥0.05: 랜덤 효과 선택
# ── 변화점수 vs ANCOVA ──
set.seed(1)
n2 <- 50
df_2t <- data.frame(
group = rep(c(0, 1), each = n2),
y_pre = rnorm(2*n2, 10, 2)
)
df_2t$y_post <- df_2t$y_pre + ifelse(df_2t$group == 1, 2, 0) + rnorm(2*n2, 0, 1.5)
df_2t$diff <- df_2t$y_post - df_2t$y_pre
summary(lm(diff ~ group, df_2t)) # 변화점수 분석
summary(lm(y_post ~ group + y_pre, df_2t)) # ANCOVA
# group 계수가 거의 동일함을 확인15 분석 흐름 요약
종단 데이터 수집
↓
1. 데이터 구조 파악
N (피험자), n_i (시점), 균형 여부, 결측 패턴, 공변량 종류
↓
2. 탐색적 분석 (EDA)
- Spaghetti plot: 개인별 궤적
- Mean profile: 시점별 평균 ± SE
- 경험적 공분산/상관 행렬 시각화
- ANOVA 기반 ICC 추정
↓
3. 분석 방법 선택 (Hedeker 1.5 기준)
결과 형태 / N / n_i / 공변량 / 분산-공분산 구조 고려
↓
4. 공분산 구조 선택
CS / AR(1) / Toeplitz / UN → REML 기반 AIC/BIC 비교
↓
5. 분석 목적 판단
- Population-average 효과 → GEE
- Individual 예측 + 랜덤 효과 분리 → LMM/GLMM
- 내생성/시불변 혼란변수 통제 → 고정 효과 모형
- 내생성 검정 → Hausman test
↓
6. 모형 추정 및 진단
REML 추정, 잔차 정규성, Q-Q plot, QIC (GEE)
↓
7. 결과 해석
고정 효과 β, 랜덤 효과 σ_b², σ_e², ICC
↓
8. 민감도 분석
MAR vs MNAR 가정, 공분산 구조 변경 시 안정성 확인16 관련 주제
선행 지식
- GLM 개요 — 선형 회귀와 일반화 선형 모형의 기초
- 반복측정 ANOVA — 종단 분석의 전통적 접근
후속 주제
- 01 — 왜 Mixed Model인가 (GLM 한계, ICC, Fixed/Random)
- 02 — 모델 구조 (Random Intercept, Slope, 공분산)
- 03 — 추정과 모델 선택 (ML/REML, LRT, AIC/BIC)
- 08 — GEE (Marginal 효과, Working Correlation)
- 11 — LMM/GLMM/GEE/GAMM 비교 및 선택 가이드
관련 개념