§ 6.4 — Bock WPSS 예시: 5 모형 적합 비교 + UN 임상 해석

교차 설계 우울증 데이터 · 직교 대비 · group × change of slope · conditional vs marginal 분산

Hedeker & Gibbons (2006) Ch.6 §6.4 의 자세한 풀이. Bock (1983b) 의 정신과 우울증 임상 데이터 (75 명, 6 주 추적, WPSS 척도) 를 사용해 5 가지 공분산 구조 (UN, Toeplitz, AR(1), CS) 를 적합하고 비교하는 케이스 스터디다. 두 치료 그룹의 cross-over 설계 (TCA-None vs None-TCA), 직교 대비 코딩 (linear trend, change of slope), group × change of slope 상호작용이 약물 효과 검정의 도구인 이유, UN 모형 적합 결과의 임상적 해석 (per-week change 계산, 약물·무약물 기간 호전 속도 비교) 까지 정리한다. 추정 conditional 분산 vs marginal 분산의 차이도 함께 다룬다.

Statistics
저자

Kwangmin Kim

공개

2026년 04월 30일

1 들어가며 — 왜 이 예시가 중요한가

Ch.6 의 5 가지 공분산 구조 — CS·AR(1), Toeplitz·UN, RE — 와 모형 선택 절차 를 모두 다뤘다. 이제 실제 임상 데이터에서 이 절차가 어떻게 작동하는지 본다.

§6.4 는 Bock (1983b) 의 우울증 임상 자료를 사용한 케이스 스터디다. 이 예시가 LDA 학습에서 차지하는 가치:

  1. cross-over 설계 — 두 그룹의 약물·무약물 순서가 정반대. group × time 상호작용이 약물 효과 검정의 정공법.
  2. 직교 대비 (orthogonal contrasts) — linear trend 와 change of slope 가 어떻게 6 시점을 분해하는지 명시.
  3. 5 구조 적합 비교 — 4 fully specified 구조 (UN, Toeplitz, AR(1), CS) 의 deviance·AIC 가 모두 같은 결론 (UN 채택) 을 가리킴.
  4. 임상 해석 — UN 모형의 회귀 계수가 per-week 변화율로 해석되어 “약물 효과” 라는 임상 결론으로 이어지는 흐름.
  5. conditional vs marginal 분산 — 추정 \(\widehat\Sigma\) 가 관측 \(\Sigma\) 보다 일관되게 작은 이유의 통계학적 직관.
한 줄 요약

“5 구조 모두 적합 → UN 만 데이터에 충분히 부합 → UN 모형의 group × change of slope 상호작용이 유의 → 두 그룹 모두 약물 기간에 호전 빠름 → 약물 효과 입증.”

이 한 줄을 수식·표·코드로 풀어내는 것이 이 sub-post 의 목표다.

2 데이터 — Bock 정신과 임상 자료

2.1 연구 설계

Bock (1983b) WPSS 데이터 요약
  • 표본: \(N = 75\) 우울증 환자
  • 추적 기간: 6 주, 매주 측정
  • 종속변수: WPSS (Weekly Psychiatric Status Scale) — 1~6 점 정신상태 평정
    • 1: 정상 / 2: 잔여 증상 / 3: 부분 회복 / 4: 뚜렷한 증상 / 5: 주요 우울 진단 기준 / 6: 극심한 손상
  • 두 그룹 (cross-over 설계):
    • TCA-None (\(n=46\)): 첫 3 주 약물 (TCA) → 다음 3 주 무약물
    • None-TCA (\(n=29\)): 첫 3 주 무약물 → 다음 3 주 약물 (TCA)
  • 결측: 모든 환자가 6 주 모두 관측 (완전 데이터)
cross-over 설계의 직관

이 설계의 핵심: 두 그룹이 약물·무약물 순서가 정반대.

Week:    1   2   3   4   5   6
TCA-None: [-- 약물 --] [-- 무약물 --]
None-TCA: [-- 무약물 -] [-- 약물 ---]

만약 약물이 정말 효과 있다면:

  • TCA-None: 첫 3 주 (약물) 에 호전 빠름, 다음 3 주 (무약물) 에 호전 느림.
  • None-TCA: 첫 3 주 (무약물) 에 호전 느림, 다음 3 주 (약물) 에 호전 빠름.

→ 두 그룹의 “전반-후반 기울기 변화” 패턴이 정반대 가 되어야 한다. 이것이 group × (전반-후반 기울기 차이) 상호작용이 약물 효과 검정의 도구인 이유.

2.2 Table 6.1 — 관측 통계

관측 평균 (mean), 표준편차 (SD), 상관 행렬

평균 (그룹별, 시점별):

그룹 \(N\) Wk1 Wk2 Wk3 Wk4 Wk5 Wk6
TCA-None 46 3.76 3.46 3.11 2.89 2.80 2.74
None-TCA 29 4.72 4.62 4.55 4.45 4.21 3.90

SD (시점별, 두 그룹 합산): 1.30, 1.40, 1.53, 1.61, 1.66, 1.65.

상관 행렬:

       Wk1   Wk2   Wk3   Wk4   Wk5   Wk6
Wk1   1.00
Wk2   0.91  1.00
Wk3   0.75  0.87  1.00
Wk4   0.68  0.82  0.91  1.00
Wk5   0.59  0.70  0.78  0.88  1.00
Wk6   0.60  0.68  0.72  0.84  0.96  1.00
데이터에서 읽히는 패턴

평균 패턴:

  • TCA-None: Wk1 → Wk6 에서 1.02 점 호전 (3.76 → 2.74). Wk1 ~ Wk3 (약물) 더 빠른 호전.
  • None-TCA: Wk1 → Wk6 에서 0.82 점 호전 (4.72 → 3.90). Wk4 ~ Wk6 (약물) 더 빠른 호전.

SD 패턴: 시점 후반으로 갈수록 SD 증가 (1.30 → 1.65). 약 27% 증가. → 시점별 분산 동일 이라는 CS·AR(1)·Toeplitz 가정 위반.

상관 패턴: lag 증가에 따라 감소하지만 매끄럽지 않음.

  • lag-1 상관: 0.91, 0.87, 0.91, 0.88, 0.96 (시점에 따라 변동, 후반 강해짐).
  • → 같은 lag 의 균질 상관 이라는 Toeplitz 가정 위반.

이 두 패턴이 미리 4 fully specified 구조 (CS·AR(1)·Toeplitz) 가 부적합할 것을 시사한다.

3 고정 효과 모형 — 직교 대비 코딩

3.1 선택된 contrasts

Bock 의 두 가지 시간 contrast

Linear trend — 6 주 전체의 평균 선형 추세:

Wk1 Wk2 Wk3 Wk4 Wk5 Wk6
Linear -5/2 -3/2 -1/2 1/2 3/2 5/2
  • Sum = 0 (week 3.5 에 중심).
  • 절편이 week 3.5 의 평균 의미.

Change of slope — 첫 3 주 vs 다음 3 주의 기울기 변화:

Wk1 Wk2 Wk3 Wk4 Wk5 Wk6
Change -1/2 0 1/2 1/2 0 -1/2
change of slope contrast 의 의미 — 직관

이 contrast 는 두 부분으로 나뉜다.

전반 3 주 (Wk1, Wk2, Wk3): -1/2, 0, 1/2 — 이 자체가 전반 기간의 선형 추세 contrast. 후반 3 주 (Wk4, Wk5, Wk6): 1/2, 0, -1/2 — 후반 기간의 선형 추세 contrast (부호 반전).

두 부분의 합이 “전반 - 후반 기울기 차이” 를 측정한다.

“TCA-None 그룹은 전반 3 주 (약물) 에 호전 가속 → linear trend 와 change of slope 가 모두 음수. None-TCA 그룹은 후반 3 주 (약물) 에 호전 가속 → change of slope 부호가 다른 방향.”

→ group × change of slope 상호작용이 약물 효과 검정의 정공법이 된다.

3.2 전체 고정 효과 명세

5 개 fixed effects + 절편:

  1. 절편 (\(\beta_0\)): week 3.5 평균 (TCA-None 그룹).
  2. Linear trend (\(\beta_1\)): TCA-None 의 6 주 전체 평균 선형 추세.
  3. Change of slope (\(\beta_2\)): TCA-None 의 전반-후반 기울기 차이.
  4. Group (\(\beta_3\)): None-TCA vs TCA-None 의 전체 평균 차이.
  5. Group × linear trend (\(\beta_4\)): 두 그룹의 전체 선형 추세 차이.
  6. Group × change of slope (\(\beta_5\)): 두 그룹의 기울기 변화 차이 → 약물 효과 검정.

4 5 구조 적합 — Table 6.2

ML 과 REML deviance
구조 \(q\) ML \(-2\log L\) REML \(-2\log L\)
UN 21 945.9 963.1
Toeplitz 6 988.9 1005.3
AR(1) 2 996.3 1013.0
CS 2 1185.8 1204.0

4.1 LR 검정 (UN 기준)

§ 6.3 모형 선택 sub-post 의 절차에 따라:

비교 \(\chi^2\) (ML) df 보정 p-value 결과
Toeplitz vs UN 988.9 - 945.9 = 43.0 21 - 6 = 15 0.000079 UN 채택
AR(1) vs UN 996.3 - 945.9 = 50.4 21 - 2 = 19 0.000057 UN 채택
CS vs UN 1185.8 - 945.9 = 239.9 21 - 2 = 19 < 0.0001 UN 채택

3 가지 fully specified 제약 모두 데이터와 통계적으로 부합 안 함. UN 사용이 강제된다.

왜 4 fully specified 모두 실패했나 — 데이터 패턴과의 충돌
가정 위반된 데이터 패턴
CS: 모든 lag 동일 상관 상관이 lag 에 따라 명확히 감소
AR(1): 지수 감쇠 + 정상성 SD 시점 증가 + lag 별 상관 비균질
Toeplitz: 시점 위치 무관 같은 lag 의 상관도 시점에 따라 변동
모두: 분산 동일 SD 시점 27% 증가

데이터의 두 가지 핵심 패턴 — (1) 시점별 분산 증가, (2) 같은 lag 의 상관 비균질 — 이 모두 fully specified 구조의 정상성 가정을 깬다.

→ UN 만이 이 두 패턴을 모두 표현 가능.

4.2 AIC 기준 일관 검증

ML deviance 기준 AIC = -2 log L + 2q:
- UN:        945.9 + 2(21) =  987.9   (가장 작음)
- Toeplitz:  988.9 + 2( 6) = 1000.9
- AR(1):     996.3 + 2( 2) = 1000.3
- CS:       1185.8 + 2( 2) = 1189.8

UN 의 AIC 가 압도적으로 작음 → AIC 도 UN 을 가리킴. LR + AIC 일치하므로 UN 채택의 강한 다중 증거.

5 UN 모형 — 고정 효과 추정 (Table 6.3)

UN 적합 회귀 계수
모수 ML 추정 SE \(Z\) \(p\)
Intercept 3.122 0.179 17.44 < .0001
Linear trend -0.198 0.036 -5.48 < .0001
Change of slope -0.255 0.102 -2.49 .015
Group 1.286 0.288 4.46 < .0001
Group × linear trend 0.017 0.058 0.28 .78
Group × change of slope 0.475 0.164 2.89 .005

Note: \(-2\log L = 945.9\).

5.1 결과의 4 단 해석

4 가지 핵심 결과

1. Linear trend 유의 (\(p < 0.0001\)):

  • \(\hat\beta_1 = -0.198\) — TCA-None 그룹의 6 주 전체 평균 호전 추세.
  • 매주 약 0.2 점씩 호전 (음수 = WPSS 점수 감소 = 증상 호전).

2. Change of slope 유의 (\(p = 0.015\)):

  • \(\hat\beta_2 = -0.255\) — TCA-None 그룹의 전반-후반 기울기 차이.
  • 음수 = 전반 기간 (약물) 의 기울기가 후반 (무약물) 보다 더 강한 음수 → 약물 기간에 호전 빠름.

3. Group 유의 (\(p < 0.0001\)):

  • \(\hat\beta_3 = 1.286\) — None-TCA 가 TCA-None 보다 평균 1.29 점 높음 (시간 평균).
  • 무작위 배정이 아닌 관찰 연구이므로 baseline 차이로 해석.

4. Group × change of slope 유의 (\(p = 0.005\)) ← 약물 효과 입증:

  • \(\hat\beta_5 = 0.475\) — 두 그룹의 기울기 변화 패턴이 정반대.
  • 양수 = TCA-None 의 음의 change of slope 와 None-TCA 의 양의 change of slope.
  • 즉 두 그룹 모두 약물 기간에 호전이 더 빠름 (cross-over 설계의 약물 효과).

5. Group × linear trend 비유의 (\(p = 0.78\)):

  • \(\hat\beta_4 = 0.017\) — 두 그룹의 전체 선형 추세는 비슷.
  • 두 그룹 모두 시간이 지나며 호전 → 자연 회복 + 약물 효과 혼재.

5.2 per-week 변화율 계산 — 임상 해석

TCA-None 그룹의 시점별 평균

추정 평균: \(\hat\mu_t = \hat\beta_0 + L_t \hat\beta_1 + C_t \hat\beta_2\)

여기서 \(L_t, C_t\) 는 각 시점의 linear, change of slope 대비 값.

Week 1 → Week 3 (전반, 약물 기간):

\[ \hat\mu_3 - \hat\mu_1 = (L_3 - L_1)\hat\beta_1 + (C_3 - C_1)\hat\beta_2 = 2(-0.198) + 1(-0.255) = -0.651 \]

per week change = \(-0.651 / 2 = -0.326\) (점/주).

Week 4 → Week 6 (후반, 무약물 기간):

\[ \hat\mu_6 - \hat\mu_4 = (L_6 - L_4)\hat\beta_1 + (C_6 - C_4)\hat\beta_2 = 2(-0.198) + (-1)(-0.255) = -0.141 \]

per week change = \(-0.141 / 2 = -0.071\) (점/주).

→ TCA-None 은 약물 기간에 약 4.6 배 빠른 호전 (\(-0.326 / -0.071 \approx 4.6\)).

None-TCA 그룹의 시점별 평균

None-TCA 의 절편·기울기는 TCA-None 절편·기울기 + 그룹 효과:

  • 절편: \(3.122 + 1.286 = 4.408\)
  • Linear trend: \(-0.198 + 0.017 = -0.181\)
  • Change of slope: \(-0.255 + 0.475 = 0.220\) (양수로 반전)

Week 1 → Week 3 (전반, 무약물 기간):

per week change = \(-0.181 + 0.5(0.220) = -0.071\) (점/주).

Week 4 → Week 6 (후반, 약물 기간):

per week change = \(-0.181 - 0.5(0.220) = -0.291\) (점/주).

→ None-TCA 도 약물 기간에 약 4.1 배 빠른 호전 (\(-0.291 / -0.071 \approx 4.1\)).

두 그룹의 약물 효과 일관성
그룹 무약물 기간 변화 약물 기간 변화 비율
TCA-None -0.071 (week 4-6) -0.326 (week 1-3) 4.6×
None-TCA -0.071 (week 1-3) -0.291 (week 4-6) 4.1×

놀랍게도 무약물 기간의 변화율 (-0.071) 이 두 그룹에서 동일. 약물 기간의 변화율도 비슷 (-0.29 ~ -0.33). cross-over 설계가 잘 작동했음을 보여주는 증거.

임상 결론: 두 그룹 모두 약물 기간 호전 속도가 무약물 기간보다 4 배 이상 빠름 → 약물 효과의 강한 증거.

6 추정 분산-공분산 행렬 (식 6.9, 6.10)

6.1 분산 — 시점별 변동

\(\widehat\Sigma\) 의 대각 
대각 (분산):  1.443  1.605  1.810  1.995  2.241  2.369
표준편차:     1.201  1.267  1.345  1.413  1.497  1.539

시점에 따라 분산이 매끄럽게 증가 — 표준편차가 6 주 만에 약 28% 증가.

“왜 시점 후반에 분산이 커지나” 의 임상 해석
  • 시점 초반 (Wk1, Wk2): 모두 우울증 환자 → 점수가 비슷한 범위에 몰림 → 분산 작음.
  • 시점 후반 (Wk5, Wk6): 호전된 환자 (낮은 점수) 와 호전 안 된 환자 (높은 점수) 가 분리 → 분산 큼.

이런 “fan-out” 패턴은 종단 임상 데이터에서 흔하다 — 처방 반응의 개인차가 시간이 지나며 누적.

6.2 추정 상관 — 식 (6.10)

       Wk1   Wk2   Wk3   Wk4   Wk5   Wk6
Wk1   1.000
Wk2   0.894 1.000
Wk3   0.699 0.836 1.000
Wk4   0.624 0.776 0.886 1.000
Wk5   0.522 0.641 0.725 0.847 1.000
Wk6   0.543 0.624 0.675 0.822 0.951 1.000

두 가지 패턴:

  1. lag 증가에 따라 상관 감소 — 거리가 멀수록 약함 (당연).
  2. 같은 lag 의 상관도 시점에 따라 변동 — 후반으로 갈수록 강해짐.
lag-1 상관의 시점별 변동 
Wk1-Wk2: 0.894
Wk2-Wk3: 0.836
Wk3-Wk4: 0.886
Wk4-Wk5: 0.847
Wk5-Wk6: 0.951

같은 lag-1 인데 0.836 ~ 0.951 범위. 이런 비균질성이 Toeplitz 의 가정 (같은 lag = 같은 상관) 을 깬다.

특히 Wk5-Wk6 이 0.951 로 매우 강한 — 마지막 두 주에 환자 상태가 안정되어 측정값이 거의 변하지 않음.

7 Conditional vs Marginal 분산 — 핵심 직관

7.1 관측값과 추정값의 차이

추정값 (식 6.9) vs 관측값 (Table 6.1) 
관측 SD (Table 6.1):  1.30  1.40  1.53  1.61  1.66  1.65
추정 SD (식 6.9):     1.20  1.27  1.35  1.41  1.50  1.54

추정 SD 가 일관되게 관측 SD 보다 작다.

7.2 왜 일관되게 작은가

conditional vs marginal 분산의 차이
  • 관측 SD: \(y\) 의 marginal 분산 — 그룹 차이·시간 효과·약물 효과를 모두 포함한 전체 변동.
  • 추정 \(\widehat\Sigma\): \(y \mid X\) 의 conditional 분산 — 공변량 (\(X\)) 효과를 빼고 남은 잔차 변동.

수식적으로:

\[ \text{Var}_{\text{marginal}}(y) = \text{Var}(\mu(X)) + E[\text{Var}_{\text{cond}}(y \mid X)] \]

(전체 분산 = 평균 변동의 분산 + 조건부 분산의 평균.)

→ marginal 분산 ≥ conditional 분산 이 항상 성립.

일상 비유 — 학년별 점수

학교 전체 학생의 시험 점수 분산 vs 같은 학년 안의 점수 분산을 비교하자.

  • 전체 분산 (marginal): 1 학년부터 6 학년까지 점수 차이를 모두 포함 → 큼.
  • 같은 학년 안 분산 (conditional): 학년 간 평균 차이는 빠지고 학년 내 변동만 → 작음.

학년 평균이 클수록 (학년 간 차이 큼) 두 분산의 차이도 커짐.

Bock 데이터에서 같은 원리:

  • 관측 SD: 그룹 차이 (1.29 점) + 시간 추세 (-0.198/주) + change of slope (-0.255) 효과를 모두 포함.
  • 추정 \(\widehat\Sigma\): 이 효과들을 빼고 남은 잔차 분산.

공변량 효과가 클수록 두 분산의 차이가 큼 — 본 사례에서는 약 8% 작아짐 (1.30 → 1.20 등).

7.3 일반 회귀에서의 동일 원리

이 차이는 다중 회귀에서도 정확히 같은 형태로 나타난다.

다중 회귀의 잔차 분산

다중 회귀 \(y = X\beta + \varepsilon\), \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\) 에서:

  • \(\text{Var}(y) = \text{Var}(X\beta) + \sigma^2\) (분산 분해).
  • 잔차 분산 \(\hat\sigma^2\) 는 항상 \(y\) 의 표본 분산보다 작거나 같다.
  • \(R^2\) 가 클수록 차이가 큼.

CPM 의 conditional vs marginal 분산은 이 직관의 시점별·다변량 일반화.

8 코드 재현 — Bock 분석

8.1 Step 1: 직교 대비 행렬 + 데이터 시뮬레이션

import numpy as np
import pandas as pd

# Bock contrasts
linear = np.array([-5/2, -3/2, -1/2, 1/2, 3/2, 5/2])
change = np.array([-1/2, 0, 1/2, 1/2, 0, -1/2])

# 두 contrast 의 직교성 검증 — sum of products = 0
print(f"sum(linear * change) = {np.sum(linear * change):.6f}")
# 0.0 — 두 contrast 가 직교

# 진짜 모수 (Hedeker Table 6.3)
beta = {
    "intercept": 3.122,
    "linear":   -0.198,
    "change":   -0.255,
    "group":     1.286,
    "g_lin":     0.017,
    "g_chg":     0.475,
}

# 시뮬레이션 분산-공분산 (식 6.9 의 핵심 구조)
sd = np.array([1.20, 1.27, 1.35, 1.41, 1.50, 1.54])
corr = np.array([
    [1.000, 0.894, 0.699, 0.624, 0.522, 0.543],
    [0.894, 1.000, 0.836, 0.776, 0.641, 0.624],
    [0.699, 0.836, 1.000, 0.886, 0.725, 0.675],
    [0.624, 0.776, 0.886, 1.000, 0.847, 0.822],
    [0.522, 0.641, 0.725, 0.847, 1.000, 0.951],
    [0.543, 0.624, 0.675, 0.822, 0.951, 1.000],
])
Sigma = np.outer(sd, sd) * corr

# 그룹 코딩
n_subj_TCA, n_subj_None = 46, 29
group = np.r_[np.zeros(n_subj_TCA), np.ones(n_subj_None)]

# 평균 곡선
mu = (beta["intercept"]
      + beta["linear"]   * linear[None, :]
      + beta["change"]   * change[None, :]
      + beta["group"]    * group[:, None]
      + beta["g_lin"]    * group[:, None] * linear[None, :]
      + beta["g_chg"]    * group[:, None] * change[None, :])

# 합성 데이터
np.random.seed(2026)
y = mu + np.random.multivariate_normal(np.zeros(6), Sigma, size=len(group))

# Long format
df = pd.DataFrame({
    "id":    np.repeat(np.arange(len(group)), 6),
    "week":  np.tile(np.arange(1, 7), len(group)),
    "y":     y.flatten(),
    "group": np.repeat(group, 6),
    "linear": np.tile(linear, len(group)),
    "change": np.tile(change, len(group)),
})
print(df.head(12))

8.2 Step 2: per-week 변화율 계산 함수

def per_week_change(beta_lin: float, beta_chg: float,
                    weeks: tuple[int, int]) -> float:
    """두 시점 사이의 평균 per-week 변화율

    공식: (mu_t2 - mu_t1) / (t2 - t1)
    여기서 mu_t = beta_0 + L_t * beta_lin + C_t * beta_chg.
    공변량 효과만 계산 (절편 상쇄).
    """
    L = np.array([-5/2, -3/2, -1/2, 1/2, 3/2, 5/2])
    C = np.array([-1/2, 0, 1/2, 1/2, 0, -1/2])
    t1, t2 = weeks[0] - 1, weeks[1] - 1  # 0-indexed
    delta_mu = (L[t2] - L[t1]) * beta_lin + (C[t2] - C[t1]) * beta_chg
    return delta_mu / (weeks[1] - weeks[0])


# TCA-None 그룹 (linear=-0.198, change=-0.255)
print(f"TCA-None Wk1→Wk3 (약물): {per_week_change(-0.198, -0.255, (1, 3)):.3f}")
# 약 -0.326

print(f"TCA-None Wk4→Wk6 (무약물): {per_week_change(-0.198, -0.255, (4, 6)):.3f}")
# 약 -0.071

# None-TCA 그룹 (linear=-0.198+0.017=-0.181, change=-0.255+0.475=0.220)
print(f"None-TCA Wk1→Wk3 (무약물): {per_week_change(-0.181, 0.220, (1, 3)):.3f}")
# 약 -0.071

print(f"None-TCA Wk4→Wk6 (약물): {per_week_change(-0.181, 0.220, (4, 6)):.3f}")
# 약 -0.291
검증 포인트

위 계산이 Hedeker Table 6.3 의 임상 해석 과 일치해야 한다.

  • TCA-None 약물: -0.326 (4.6× 빠름)
  • TCA-None 무약물: -0.071
  • None-TCA 무약물: -0.071 (TCA-None 무약물과 정확히 일치)
  • None-TCA 약물: -0.291 (4.1× 빠름)

두 그룹의 무약물 변화율이 정확히 같다는 사실이 cross-over 설계의 일관성을 확인한다.

8.3 Step 3: R 에서 5 모형 일괄 적합

library(nlme)

# 데이터 로드 (위 Python 시뮬레이션 결과)
df <- read.csv("bock_simulated.csv")

# 직교 대비 변수
df$linear <- with(df, -5/2 + (week - 1))  # -5/2, -3/2, -1/2, 1/2, 3/2, 5/2
df$change <- with(df, c(-1/2, 0, 1/2, 1/2, 0, -1/2)[week])

# UN — corSymm + varIdent
m_un <- gls(y ~ linear + change + group + group:linear + group:change,
            data = df,
            correlation = corSymm(form = ~ 1 | id),
            weights = varIdent(form = ~ 1 | week),
            method = "ML")

# Toeplitz — corARMA(p=0, q=5)
m_toep <- gls(y ~ linear + change + group + group:linear + group:change,
              data = df,
              correlation = corARMA(form = ~ week | id, p = 0, q = 5),
              method = "ML")

# AR(1) — corAR1
m_ar1 <- gls(y ~ linear + change + group + group:linear + group:change,
             data = df,
             correlation = corAR1(form = ~ week | id),
             method = "ML")

# CS — corCompSymm
m_cs <- gls(y ~ linear + change + group + group:linear + group:change,
            data = df,
            correlation = corCompSymm(form = ~ 1 | id),
            method = "ML")

# RE (랜덤 절편 + 기울기) — lme
library(nlme)
m_re <- lme(y ~ linear + change + group + group:linear + group:change,
            random = ~ week | id,
            data = df,
            method = "ML")

# 5 모형 비교
cat("Deviances (ML):\n")
cat(sprintf("UN:        %.1f\n", -2 * logLik(m_un)))
cat(sprintf("Toeplitz:  %.1f\n", -2 * logLik(m_toep)))
cat(sprintf("AR(1):     %.1f\n", -2 * logLik(m_ar1)))
cat(sprintf("CS:        %.1f\n", -2 * logLik(m_cs)))
cat(sprintf("RE(r=2):   %.1f\n", -2 * logLik(m_re)))

# AIC/BIC
AIC(m_un, m_toep, m_ar1, m_cs, m_re)
BIC(m_un, m_toep, m_ar1, m_cs, m_re)

# UN 의 고정 효과
summary(m_un)
RE 결과의 흥미

RE(\(r=2\), 랜덤 절편 + 기울기) 는 Hedeker Table 6.2 에 포함되어 있지 않다. 실제 적합하면 deviance 가 UN 보다 약간 크고 Toeplitz 보다 작은 위치에 놓일 가능성이 높다 — 시점별 분산 변동을 일부 표현 (랜덤 기울기) 하지만 같은 lag 상관 비균질성은 표현 못함.

실무에서는 RE 결과를 fully specified 4 구조와 함께 비교하는 것이 일반 절차.

9 5 모형의 데이터별 적합 패턴

9.1 그래픽 진단 — 어느 구조가 어디서 부족한가

시점별 SD 그림 — UN vs 다른 구조 
       Wk1   Wk2   Wk3   Wk4   Wk5   Wk6
관측:   1.30  1.40  1.53  1.61  1.66  1.65
UN:     1.20  1.27  1.35  1.41  1.50  1.54   (시점 변동 OK)
Toep:   1.43  1.43  1.43  1.43  1.43  1.43   (정상성 — 변동 못함)
AR(1):  1.43  1.43  1.43  1.43  1.43  1.43   (정상성 — 변동 못함)
CS:     1.43  1.43  1.43  1.43  1.43  1.43   (정상성 — 변동 못함)

UN 만이 시점별 SD 변동을 자연스럽게 표현. 나머지 3 구조는 일정한 SD 만 추정.

lag-1 상관 그림 
       Wk1-2 Wk2-3 Wk3-4 Wk4-5 Wk5-6
관측:   0.91  0.87  0.91  0.88  0.96
UN:     0.89  0.84  0.89  0.85  0.95   (시점 변동 OK)
Toep:   0.85  0.85  0.85  0.85  0.85   (lag-1 일정 — 변동 못함)
AR(1):  0.85  0.85  0.85  0.85  0.85   (lag-1 일정)
CS:     0.74  0.74  0.74  0.74  0.74   (모든 lag 일정 → 평균값)

여기서도 UN 만 시점별 lag-1 상관 변동 (특히 Wk5-Wk6 의 0.95) 을 표현.

9.2 CS 가 가장 나쁜 이유 정리

CS 의 deviance (1185.8) 가 다른 3 구조보다 압도적으로 큰 이유:

  1. 시점별 분산 변동 표현 못함 (다른 3 구조도 같음, 불리).
  2. lag 별 상관 변동 표현 못함 — AR(1)·Toeplitz 와 달리 모든 lag 동일 → 데이터의 “lag 따라 감소” 패턴 자체를 못 잡음.
  3. 추가 자유도가 적음 (\(q=2\), AR(1) 과 같지만 lag 정보 자체 사용 안 함).

→ CS = 가장 단순한 모형, 데이터가 가장 단순한 패턴 (모든 시점 동등) 일 때만 적합.

10 핵심 정리

한 페이지 요약
  1. 데이터: Bock (1983b) 75 명 우울증 환자, 6 주 추적, WPSS 척도, cross-over 약물 설계.
  2. 고정 효과: 직교 대비 (linear trend, change of slope) + group + 두 상호작용 = 6 모수.
  3. 5 구조 적합: UN deviance 945.9 압도적 작음 → 다른 4 구조 모두 LR 검정 유의 부족 (보정 p < 0.0001).
  4. 부적합 원인: 데이터의 (a) 시점별 분산 증가, (b) 같은 lag 상관 비균질성 — fully specified 의 정상성 가정 위반.
  5. UN 회귀 계수: linear trend 유의, group 유의, group × change of slope 유의 (\(p = 0.005\)) — 약물 효과 입증.
  6. per-week 변화율: 두 그룹 모두 약물 기간 약 -0.3 점/주, 무약물 기간 약 -0.07 점/주. 약 4 배 빠른 호전.
  7. conditional vs marginal: 추정 \(\widehat\Sigma\) < 관측 \(\Sigma\) — 공변량 효과 제거된 잔차 분산. 다중 회귀와 동일 원리.
  8. 추정 분산 패턴: SD 시점에 따라 1.20 → 1.54 (28% 증가) — 환자 반응 fan-out.
  9. 추정 상관 패턴: lag 증가하며 감소 + 후반 lag-1 강해짐 (Wk5-Wk6 0.951).
  10. 케이스 스터디 가치: 5 구조 비교 → 모형 선택 → 임상 해석 → 분산 구조 패턴 이해까지 한 사례에서 다 다룸.

Bock WPSS 예시는 Hedeker Ch.6 의 화룡점정 — 5 구조 이론을 실제 데이터에 적용하여 cross-over 임상 시험의 약물 효과를 통계적으로 입증하는 절차 전체를 보여준다.

11 Ch.6 마무리 — sub-post 시리즈 완성

이 sub-post 로 Ch.6 §6.1 ~ §6.4 의 모든 내용이 sub-post 시리즈로 정리됐다.

sub-post 다룬 섹션 핵심
Overview §6.1, §6.2 (5 구조 비교), §6.3 요약, §6.4 요약 “분산-공분산을 어떻게 모수화할까” — 5 구조 한 페이지
§ 6.2.1-6.2.2 CS, AR(1) “두 절약 구조 — lag 무시 vs 지수 감쇠”
§ 6.2.3-6.2.4 Toeplitz, UN “두 유연 구조 — 가정을 단계적으로 풀기”
§ 6.2.5 RE “메커니즘 기반 구조 — MRM 과의 다리”
§ 6.3 Model Selection “UN 기준 LR 검정과 2-step 절차”
이 포스트 § 6.4 Bock 예시 “5 구조의 데이터 적합 + 임상 해석”
다음 단계

§ 6.5 Summary 와 Ch.7 Overview 가 다음 자연스러운 흐름이다.

  • Ch.7 — RE + AC 결합 모형 (CPM + MRM 결합).
  • Ch.8 — GEE (CPM 의 비정규 일반화).

§6.5 는 SAS PROC MIXED 의 추가 구조 (NS-AR(1), 이질적 분산 CPM 등) 를 간략히 소개하는 짧은 마무리 섹션 — 별도 sub-post 작성 가치 낮음. Ch.7 §7.1 의 일부로 흡수하거나 본 시리즈에서 생략해도 무방.

12 관련 주제

선행 지식

관련

후속 주제

교재

  • Hedeker, D. & Gibbons, R. D. (2006). Longitudinal Data Analysis, Wiley, Ch.6 §6.4 (pp. 105-110)
  • Bock, R. D. (1983b). “The discrete Bayesian”, in H. Wainer & S. Messick (eds.), Modern Advances in Psychometric Research, Erlbaum — WPSS 데이터 출처
  • Wolfinger, R. (1993). “Covariance structure selection in general mixed models”, Communications in Statistics — Simulation and Computation 22, 1079-1106 — SAS PROC MIXED 의 추가 구조 (§6.5 보조)

Subscribe

Enjoy this blog? Get notified of new posts by email: