1 도입 — Ch.7 의 위치
Hedeker 책 Ch.4-5 에서 다룬 표준 혼합효과 회귀 모형 (MRM) 은 다음 두 가정을 깔고 있다:
- 랜덤 효과: \(v_i \sim N(0, \Sigma_v)\) — 피험자 간 이질성.
- 조건부 독립 오차: \(\varepsilon_i \sim N(0, \sigma^2 I_{n_i})\) — random effect 가 주어지면 같은 피험자의 오차들이 독립.
Ch.6 의 공분산 패턴 모형 (CPM) 은 다른 길로 — random effect 를 안 두고 \(V(y_i)\) 를 직접 (\(\Sigma_i\) = AR(1), Toeplitz 등) 모형화. 시간적 상관을 명시.
MRM (Ch.4-5): “피험자별 위치/기울기가 다르다.” Random effect 가 그 이질성을 잡는다.
CPM (Ch.6): “같은 피험자의 시간 점들이 상관 있다.” 공분산 구조가 그 상관을 잡는다.
Ch.7: 두 시점의 결합 — 피험자별 random effect (이질성) + 시간적 자기상관 (조건부 독립 위반). 가장 유연한 분산-공분산 모형.
식으로:
\[ V(y_i) = \underbrace{Z_i \Sigma_v Z_i'}_{\text{피험자 간 이질성 (MRM)}} \;+\; \underbrace{\sigma^2 \Omega_i}_{\text{조건부 자기상관 (CPM 와 결합)}} \]
\(\Omega_i\) 가 단위 행렬이면 표준 MRM (Ch.4-5), \(Z_i = 0\) 이면 CPM (Ch.6), 둘 다 일반이면 Ch.7.
본 포스트는 § 7.2 의 일반 framework + 다섯 자기상관 구조 (AR(1), MA(1), ARMA(1,1), Toeplitz, NS-AR(1)) 를 차례로 다룬다.
2 § 7.2 일반 Framework — 식 (7.1)-(7.6)
2.1 모형 정의 — 식 (7.1)
\(i\) 번째 피험자의 \(n_i \times 1\) 반응 벡터:
\[ y_i = X_i \beta + Z_i v_i + \varepsilon_i , \quad i = 1, \ldots, N . \tag{식 7.1} \]
각 항:
- \(y_i\): \(n_i \times 1\) 반응 벡터
- \(X_i\): \(n_i \times p\) 고정 효과 설계 행렬
- \(\beta\): \(p \times 1\) 고정 회귀 모수
- \(Z_i\): \(n_i \times r\) 랜덤 효과 설계 행렬
- \(v_i\): \(r \times 1\) 랜덤 효과 (\(v_i \sim N(0, \Sigma_v)\))
- \(\varepsilon_i\): \(n_i \times 1\) 오차 벡터
2.2 표준 MRM 의 V — 식 (7.2)
오차가 conditional independent (\(\varepsilon_i \sim N(0, \sigma^2 I_{n_i})\)) 이면:
\[ V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i} . \tag{식 7.2} \]
\(V(y_i)\) 의 두 부분:
\(Z_i \Sigma_v Z_i'\): 피험자 간 차이가 만드는 분산. 같은 피험자의 다른 시점은 같은 \(v_i\) 를 공유하므로 양의 공분산 발생. 이게 피험자 간 이질성.
\(\sigma^2 I_{n_i}\): 조건부 잔차 분산. 단위 행렬이라 다른 시점의 \(\varepsilon_{ij}, \varepsilon_{ik}\) 가 독립.
표준 MRM 은 같은 피험자의 시점들이 연관되는 유일한 통로가 random effect \(v_i\). 이게 너무 강한 가정.
2.3 AC 오차 추가 — 식 (7.3)
\(\varepsilon_i \sim N(0, \sigma^2 \Omega_i)\) 로 일반화 (\(\Omega_i\) 는 자기상관 행렬):
\[ V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 \Omega_i . \tag{식 7.3} \]
\(Z_i \Sigma_v Z_i'\): 피험자별 다른 위치/기울기. 같은 피험자 내에서 시점들이 연관되는 통로 1.
\(\sigma^2 \Omega_i\): 모든 피험자에서 공통인 시간적 자기상관 패턴. 같은 피험자 내 통로 2 — random effect 외에 추가로.
핵심 관찰: 첫 항은 개체 차이 (heterogeneity), 둘째 항은 시간 구조 (temporal structure). 두 종류의 의존성이 분리되어 있다.
예시: 우울증 점수의 시간 변화. - \(Z_i v_i\): 환자별로 baseline 우울 정도 + 회복 속도가 다름 (개인차). - \(\sigma^2 \Omega_i\): 한 환자 내에서 어제 우울했으면 오늘도 우울할 가능성 (단기 momentum) — 모든 환자에 공통.
표준 MRM 은 둘째 효과를 무시. Ch.7 은 둘 다 동시 모형화.
2.4 EB Estimator + Posterior Covariance — 식 (7.4)-(7.6)
\(y_i\) 와 \(v_i\) 의 결합 분포:
\[ \begin{pmatrix} y_i \\ v_i \end{pmatrix} \sim N\!\left( \begin{pmatrix} X_i \beta \\ 0 \end{pmatrix}, \begin{pmatrix} Z_i \Sigma_v Z_i' + \sigma^2 \Omega_i & Z_i \Sigma_v \\ \Sigma_v Z_i' & \Sigma_v \end{pmatrix} \right) . \tag{식 7.4} \]
랜덤 효과의 사후 (posterior) 평균 — Empirical Bayes (EB) 추정량:
\[ \hat{v}_i = \big[Z_i'(\sigma^2 \Omega_i)^{-1} Z_i + \Sigma_v^{-1}\big]^{-1} Z_i' (\sigma^2 \Omega_i)^{-1} (y_i - X_i \beta) . \tag{식 7.5} \]
사후 공분산:
\[ \Sigma_{v|y_i} = \big[Z_i'(\sigma^2 \Omega_i)^{-1} Z_i + \Sigma_v^{-1}\big]^{-1} . \tag{식 7.6} \]
Ch.4 의 EB estimator (표준 MRM):
\[ \hat{v}_i^{(\text{Ch.4})} = \big[Z_i'(\sigma^2 I_i)^{-1} Z_i + \Sigma_v^{-1}\big]^{-1} Z_i' (\sigma^2 I_i)^{-1} (y_i - X_i \beta) . \]
식 (7.5) 와 비교하면 단위 행렬 \(I_i\) 가 자기상관 행렬 \(\Omega_i\) 로 바뀐 것뿐. 식의 골격이 동일.
이는 자기상관 모형이 표준 MRM 의 자연스러운 일반화 임을 보여준다 — 모든 식이 \(I \to \Omega\) 로 일관되게 확장. 따라서 표준 MRM 의 추정 알고리즘 (EM, Fisher Scoring) 도 거의 그대로 적용 가능.
이제 다섯 자기상관 구조를 차례로 본다.
3 § 7.2.1 AR(1) — 1차 자기회귀 오차
3.1 정의 — 식 (7.7)
\[ \varepsilon_j = \rho \varepsilon_{j-1} + \xi_j , \tag{식 7.7} \]
\(\xi_j \sim N(0, \sigma^2)\) i.i.d., \(|\rho| < 1\) (정상성).
식 (7.7) 은 단순한 회귀 형식 — \(\varepsilon_j\) 를 \(\varepsilon_{j-1}\) 에 1차 회귀.
- \(\rho > 0\): 어제 양의 오차 → 오늘도 양의 오차 경향. 시간적 momentum.
- \(\rho < 0\): 어제 양 → 오늘 음. 진동 (oscillation).
- \(\rho = 0\): 독립 (표준 MRM).
\(\rho\) 가 1 에 가까울수록 인접 시점 간 상관 강함. \(|\rho| < 1\) 정상성 조건은 분산이 폭발하지 않도록 강제.
임상 예: 우울증 점수 — 어제 우울했으면 오늘도 우울할 가능성 (\(\rho > 0\)). 하루 사이의 short-term carryover.
3.2 분산과 공분산 — 식 (7.9)-(7.12)
식 (7.7) 에서 양변의 분산:
\[ V(\varepsilon_j) = \rho^2 V(\varepsilon_{j-1}) + \sigma^2 . \tag{식 7.9} \]
정상성 (\(V(\varepsilon_j) = V(\varepsilon_{j-1})\)) 가정 하에:
\[ V(\varepsilon_j) = \frac{\sigma^2}{1 - \rho^2} . \tag{식 7.10} \]
식 (7.9) 를 무한 과거로 풀면:
\[ V(\varepsilon_j) = \rho^2 V(\varepsilon_{j-1}) + \sigma^2 = \rho^2 [\rho^2 V(\varepsilon_{j-2}) + \sigma^2] + \sigma^2 = \cdots = \sigma^2 \sum_{k=0}^\infty \rho^{2k} . \]
이 무한 합이 수렴하려면 \(|\rho| < 1\) 필요 — 그렇지 않으면 분산이 발산. 수렴 합 = \(\sigma^2/(1-\rho^2)\), 정상성 분산.
\(\rho \to 1\) 의 극한에서 분산 폭발 → AR(1) 모형이 비현실적이 됨. 따라서 \(|\rho|\) 가 1 에 너무 가까우면 모형 진단 필요 (랜덤 워크 또는 비정상성 의심).
공분산은:
\[ \text{Cov}(\varepsilon_j, \varepsilon_{j-s}) = \frac{\rho^s \sigma^2}{1 - \rho^2} . \tag{식 7.12} \]
상관:
\[ \text{Corr}(\varepsilon_j, \varepsilon_{j-s}) = \rho^s . \]
lag-\(s\) 상관이 \(\rho^s\) — 시차에 따라 지수 감쇠.
3.3 분산-공분산 행렬 — 식 (7.13)
5 시점의 예 (\(n = 5\)):
\[ \Sigma = \frac{\sigma^2}{1 - \rho^2} \begin{pmatrix} 1 & \rho & \rho^2 & \rho^3 & \rho^4 \\ \rho & 1 & \rho & \rho^2 & \rho^3 \\ \rho^2 & \rho & 1 & \rho & \rho^2 \\ \rho^3 & \rho^2 & \rho & 1 & \rho \\ \rho^4 & \rho^3 & \rho^2 & \rho & 1 \end{pmatrix} . \]
식 (7.13) 의 행렬은 두 특성을 가진다:
- Toeplitz: 같은 시차의 원소가 같음. 대각선과 평행한 띠가 모두 동일 값.
- 지수 감쇠: 시차 \(s\) 의 상관이 \(\rho^s\) — 멀어질수록 빠르게 감소.
Toeplitz 는 시간적 일관성 (정상성), 지수 감쇠는 단기 momentum 의 두 가정을 한꺼번에 표현. AR(1) 의 우아함이 이 두 성질의 결합.
자유 모수: \(\sigma^2\) + \(\rho\) = 2 모수. \(n \times n\) 행렬 (\(n(n+1)/2\) 자유) 을 단 2 개로 압축. 매우 절약적.
3.4 Parameterization 차이 — 식 (7.14), (7.15)
CPM 의 식 (6.4) 에서 사용한 AR(1) parameterization:
\[ \sigma^{*2} \cdot \rho^{|j-k|} . \]
식 (7.13) 의 parameterization:
\[ \frac{\sigma^2}{1 - \rho^2} \cdot \rho^{|j-k|} . \]
관계:
\[ \sigma^{*2} = (1 - \rho^2) \sigma^2 \quad \Leftrightarrow \quad \sigma^2 = \frac{\sigma^{*2}}{1 - \rho^2} . \tag{식 7.14, 7.15} \]
식 (7.13) parameterization (시계열·계량경제학 표준): \(\sigma^2\) 는 innovation \(\xi_j\) 의 분산.
식 (6.4) parameterization (생물통계 표준): \(\sigma^{*2}\) 는 직접 관측 분산 (대각).
같은 모형의 다른 표기. 추정 결과 (분산 등) 는 다르지만 fitted values, predictions, 검정 통계량은 같다. 소프트웨어가 어느 쪽을 쓰는지 매뉴얼 확인 필수.
R nlme::corAR1 은 식 (6.4) 형태, SAS PROC MIXED 의 TYPE=AR(1) 도 비슷. Stata 의 xtreg, ar(1) 은 식 (7.13) 형태.
4 § 7.2.2 MA(1) — 1차 이동평균 오차
4.1 정의 — 식 (7.16)
\[ \varepsilon_j = \xi_j - \theta \xi_{j-1} , \tag{식 7.16} \]
\(\xi_j \sim N(0, \sigma^2)\) i.i.d.
식 (7.16) 은 AR(1) 과 다른 메커니즘:
- AR(1): 어제 오차 \(\varepsilon_{j-1}\) 가 오늘 오차에 영향.
- MA(1): 어제 disturbance \(\xi_{j-1}\) 의 일부가 오늘 오차에 들어옴.
차이는 미묘 — disturbance (\(\xi\)) 는 latent shock, error (\(\varepsilon\)) 는 observable quantity.
실용적 차이: AR(1) 은 무한 lag 까지 상관 (\(\rho^s\), 작아지지만 0 안 됨). MA(1) 은 lag-1 까지만 상관 — lag-2 이상은 정확히 0.
임상 예: 식이 일기 — 어제 무엇을 먹었나가 오늘에 약간 영향 (남은 음식 처리 등) 이지만 그저께는 영향 X. AR(1) 보다 짧은 memory.
4.2 분산-공분산 행렬
정상성 가정 하:
\[ \sigma^2 \Omega = \sigma^2 \begin{pmatrix} 1 + \theta^2 & -\theta & 0 & \cdots & 0 \\ -\theta & 1 + \theta^2 & -\theta & \cdots & 0 \\ 0 & -\theta & 1 + \theta^2 & \cdots & 0 \\ \vdots & & & \ddots & -\theta \\ 0 & 0 & \cdots & -\theta & 1 + \theta^2 \end{pmatrix} . \]
대각: \((1 + \theta^2) \sigma^2\), 첫째 off-diagonal: \(-\theta \sigma^2\), 그 외: 0.
\(\varepsilon_j = \xi_j - \theta \xi_{j-1}\), \(\varepsilon_{j+1} = \xi_{j+1} - \theta \xi_j\).
\(\text{Cov}(\varepsilon_j, \varepsilon_{j+1}) = \text{Cov}(\xi_j - \theta \xi_{j-1}, \xi_{j+1} - \theta \xi_j) = -\theta \sigma^2\) (\(\xi_j\) 가 양쪽에 등장).
\(\varepsilon_j\) 와 \(\varepsilon_{j+2} = \xi_{j+2} - \theta \xi_{j+1}\): 공통 disturbance 없음 → 공분산 0.
따라서 lag-2 이상 정확히 0. AR(1) 의 지수 감쇠와 다른 hard cutoff.
자유 모수: \(\sigma^2 + \theta\) = 2 모수 (AR(1) 과 같음).
MA(1) 형태는 CPM (Ch.6) 의 직접 분산-공분산 모형에서는 거의 안 쓰임 — 이유: 관측 데이터의 marginal 분산-공분산 매트릭스에서 lag-1 만 상관 인 패턴은 비현실적 (긴 시간에도 약하지만 상관 잔존).
그러나 MRM 의 conditional 오차 (random effect 제거 후 잔차) 에서는 합리적 — random effect 가 장기 패턴을 잡고, MA(1) 가 단기 momentum 만 추가.
이 점이 § 7.2 가 Ch.6 와 다른 이유 — MRM 과 결합한 AC 모형이 단독 CPM 보다 더 다양한 자기상관 구조를 자연스럽게 수용.
5 § 7.2.3 ARMA(1,1) — AR + MA 결합
5.1 정의 — 식 (7.17)
\[ \varepsilon_j = \rho \varepsilon_{j-1} + \xi_j - \theta \xi_{j-1} . \tag{식 7.17} \]
AR(1) 과 MA(1) 를 모두 포함. 자유 모수: \(\rho, \theta, \sigma^2\) = 3 모수.
ARMA(1,1) 의 분산-공분산 행렬에서:
- 대각: \(\sigma^2 (1 + \theta^2 - 2\rho\theta) / (1-\rho^2)\)
- 첫 lag 공분산: \(\sigma^2 (1-\rho\theta)(\rho-\theta) / (1-\rho^2)\)
- 그 후 lag \(s\): \(\rho^{s-1}\) 에 비례 (지수 감쇠)
lag-1 의 상관이 다른 lag 보다 더 강한 hump 패턴 을 모형화.
언제 적합한가: - AR(1) 만으로는 lag-1 상관이 충분히 안 잡힐 때. - 즉 인접 시점 간 매우 강한 상관 + 그 이후 일반 AR(1) 감쇠. - 임상에서 약물 약효의 carryover 와 같은 패턴.
6 § 7.2.4 Toeplitz — 일반 띠 구조
6.1 정의
각 시차 \(s = 1, \ldots, n-1\) 에 별도의 \(\rho_s\) 모수:
\[ \sigma^2 \Omega = \sigma^2 \begin{pmatrix} 1 & \rho_1 & \rho_2 & \cdots & \rho_{n-1} \\ \rho_1 & 1 & \rho_1 & \cdots & \rho_{n-2} \\ \rho_2 & \rho_1 & 1 & \cdots & \rho_{n-3} \\ \vdots & & & \ddots & \vdots \\ \rho_{n-1} & \rho_{n-2} & \rho_{n-3} & \cdots & 1 \end{pmatrix} . \]
자유 모수: \(\sigma^2 + (n-1)\) 개 \(\rho_s\) = \(n\) 모수.
AR(1) 은 lag-\(s\) 상관이 \(\rho^s\) 로 강제. Toeplitz 는 각 lag 가 독립적으로 추정 — 자유도 높지만 모수도 많음.
언제 적합한가: - 시차에 따라 상관이 비단조 (예: lag-1 상관 강 → lag-2 약 → lag-3 다시 강) 일 때. - 주기적 패턴 (week 단위, monthly 등).
자유 모수: AR(1) 의 2 개 → Toeplitz 의 \(n\) 개. \(n = 5\) 면 5 모수, \(n = 10\) 이면 10 모수. 시점 수가 많아질수록 모수 폭발.
\(s\)-차 Toeplitz: 처음 \(s\) 개 lag 만 비zero, 나머지는 0. 모수 절약.
3 시점의 경우, MRM 의 random intercept + Toeplitz(1) 오차 와 CPM 의 full Toeplitz(3) 가 동치이다.
CPM: [θ₁ θ₂ θ₃] Random Intercept + Toeplitz(1):
[θ₂ θ₁ θ₂] = σ²_v · J + σ²[1 ρ₁ 0]
[θ₃ θ₂ θ₁] [ρ₁ 1 ρ₁]
[0 ρ₁ 1]여기서 \(J\) 는 모두 1 인 행렬.
관계: - \(\theta_1 = \sigma_v^2 + \sigma^2\) (대각: 랜덤 절편 분산 + 오차 분산) - \(\theta_2 = \sigma_v^2 + \rho_1 \sigma^2\) (lag-1: 랜덤 절편 + AR 첫 항) - \(\theta_3 = \sigma_v^2\) (lag-2: 랜덤 절편만)
두 모형은 재모수화 (reparameterization) 일 뿐. 동일 fitted values + likelihood. 따라서 CPM 과 MRM-AC 가 같은 데이터 적합도를 줄 수 있고, 어느 쪽 표현이 자연스러운지가 선택 기준.
7 § 7.2.5 비정상 (Nonstationary) AR(1)
7.1 정의 — 식 (7.18)-(7.19)
식 (7.9) 와 같지만 정상성 가정 풀고 \(V(\varepsilon_0) = 0\):
\[ V(\varepsilon_j) = \rho^2 V(\varepsilon_{j-1}) + \sigma^2 , \tag{식 7.18} \]
이로부터:
\[ V(\varepsilon_1) = \sigma^2 , \quad V(\varepsilon_2) = (1 + \rho^2) \sigma^2 , \quad V(\varepsilon_3) = (1 + \rho^2 + \rho^4) \sigma^2 , \ldots \]
정상 AR(1) 은 모든 시점의 분산이 \(\sigma^2/(1-\rho^2)\) 로 같음. NS-AR(1) 은:
- 시점 1: \(\sigma^2\) (작음)
- 시점 \(j\): \(\sigma^2 \sum_{k=0}^{j-1} \rho^{2k}\) (증가)
- \(j \to \infty\): \(\sigma^2/(1-\rho^2)\) 로 수렴 (정상 AR(1) 의 분산)
즉 NS-AR(1) 은 시작 시점의 분산이 작고 시간에 따라 정상 AR(1) 의 분산으로 수렴하는 transition.
언제 적합한가: - 임상 시험의 baseline 측정 직후 (분산이 점차 증가). - 학습 효과 (시간이 지날수록 개인 차이가 누적). - 처치 시작 후 효과 누적의 분산이 점차 커지는 패턴.
7.2 Cholesky 인수분해 — 식 (7.19)
\(\Omega = \Upsilon \Upsilon'\) 로 분해:
\[ \Upsilon = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ \rho & 1 & 0 & \cdots & 0 \\ \rho^2 & \rho & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{n-1} & \rho^{n-2} & \rho^{n-3} & \cdots & 1 \end{pmatrix} . \tag{식 7.19} \]
\(\Omega\) 직접 계산 대신 lower triangular \(\Upsilon\) 만 다루면 됨:
- 역수 \(\Omega^{-1} = (\Upsilon')^{-1} \Upsilon^{-1}\), \(\Upsilon\) 의 역수가 또 lower triangular 라 빠름.
- 우도 계산에서 \(\log \det \Omega = 2 \log \det \Upsilon = 0\) (\(\Upsilon\) 의 대각이 모두 1 이므로!).
따라서 NS-AR(1) 은 계산상 매우 효율적. MIXREG 등의 소프트웨어 (Hedeker & Gibbons, 1996b) 가 Cholesky 형태로 직접 구현.
자유 모수: \(\sigma^2 + \rho\) = 2 모수 (정상 AR(1) 과 같음).
8 다섯 구조 비교 — 자유 모수와 적합 시점
| 구조 | 자유 모수 | 핵심 가정 | 적합한 패턴 |
|---|---|---|---|
| AR(1) | 2 (\(\sigma^2, \rho\)) | 정상성, 지수 감쇠 | 인접 시점 강한 상관, 그 이후 빠른 감쇠 |
| MA(1) | 2 (\(\sigma^2, \theta\)) | 정상성, lag-1 만 상관 | Hard cutoff, 단기 momentum 만 |
| ARMA(1,1) | 3 (\(\sigma^2, \rho, \theta\)) | 정상성, lag-1 강 + 후속 지수 | 복약 carryover 등 |
| Toeplitz | \(n\) (\(\sigma^2 + (n-1) \rho_s\)) | 정상성, lag 별 독립 추정 | 비단조·주기적 상관 |
| NS-AR(1) | 2 (\(\sigma^2, \rho\)) | 비정상, 시간 따라 분산 증가 | Baseline → 효과 누적 |
1. AR(1) 가 default: 시계열 분석에서 검증된 표준. 자유 모수 절약 (2개), 해석 직관적.
2. lag-1 매우 강 + 빠른 감쇠 → ARMA(1,1): lag-1 의 hump 가 보일 때.
3. lag-1 만 상관, 그 외 0 → MA(1): 거의 드물지만 short-memory 가 명확할 때.
4. 비단조 패턴 → Toeplitz: 데이터 탐색 후 자유 모수가 충분할 때 (큰 표본).
5. 분산 증가 → NS-AR(1): 학습/누적 효과가 명확.
진단 절차: 1. 표준 MRM 적합 → 잔차 시계열 검토. 2. 잔차의 ACF/PACF 도표로 패턴 확인. 3. 후보 구조 1-2개 적합 → AIC/BIC 또는 LRT 비교. 4. fixed effects 추정의 변화 검토 (구조에 강건한지).
9 자기상관 모수 검정 — LRT
표준 MRM (without AC) 와 MRM + AC 의 LRT:
\[ \chi^2_{LR} = -2[\log L_{\text{MRM}} - \log L_{\text{MRM+AC}}] . \]
자유도 = AC 모수 수 (AR(1) 은 1, ARMA(1,1) 은 2).
경계 문제 주의: \(\rho = 0\) 이 모수공간 경계 (\(|\rho| < 1\)) 와 멀어 일반적으로 standard \(\chi^2\) 적용 OK. 그러나 \(\sigma_v^2 = 0\) 같은 경계 문제는 별도 (mixture \(\chi^2\)).
10 코드 예시
10.1 Step 1 — R: nlme 의 corAR1 / corARMA / corCompSymm
library(nlme)
# 데이터: depression 점수 longitudinal
# subject, time, score, treatment
# 1. 표준 MRM (random intercept + slope, 조건부 독립 오차)
fit_mrm <- lme(score ~ time * treatment,
random = ~ time | subject,
data = dep)
# 2. MRM + AR(1) 오차
fit_ar1 <- lme(score ~ time * treatment,
random = ~ time | subject,
correlation = corAR1(form = ~ time | subject),
data = dep)
# 3. MRM + ARMA(1,1)
fit_arma <- lme(score ~ time * treatment,
random = ~ time | subject,
correlation = corARMA(form = ~ time | subject,
p = 1, q = 1),
data = dep)
# 4. MRM + Toeplitz (lag s=2)
# nlme 는 Toeplitz 직접 안 지원, gls 의 corSymm 또는 SAS 사용 권장
# 5. LRT 비교
anova(fit_mrm, fit_ar1)
anova(fit_ar1, fit_arma)
# AIC 비교
AIC(fit_mrm, fit_ar1, fit_arma)
# AR(1) 모수 추출
fit_ar1$modelStruct$corStruct
# rho 값 출력10.2 Step 2 — Python: statsmodels MixedLM (제한적)
import statsmodels.api as sm
import statsmodels.formula.api as smf
# statsmodels 의 MixedLM 은 random effect + 표준 오차만 지원.
# AR(1) 등 자기상관 오차는 직접 지원 안 함.
# 표준 MRM
model_mrm = smf.mixedlm("score ~ time * treatment", data=df,
groups=df["subject"], re_formula="~time")
fit_mrm = model_mrm.fit(method="lbfgs")
print(fit_mrm.summary())
# AC 오차를 원하면 R 의 nlme 또는 lme4 + lmerTest 추천.
# 또는 statsmodels 의 GLSAR (단순 AR(1) 회귀, random effects 없음).10.3 Step 3 — SAS PROC MIXED (참고)
PROC MIXED DATA=dep METHOD=REML;
CLASS subject treatment;
MODEL score = time treatment time*treatment / SOLUTION;
RANDOM intercept time / SUBJECT=subject TYPE=UN;
REPEATED / SUBJECT=subject TYPE=AR(1);
RUN;
* TYPE=AR(1) → AR(1) 오차;
* TYPE=ARMA(1,1) → ARMA(1,1);
* TYPE=TOEP(s) → s-차 Toeplitz;
* TYPE=AR(1)(NS) → 비정상 AR(1) (Hedeker MIXREG);11 핵심 요약
표준 MRM 의 조건부 독립 가정 (\(\varepsilon_i \sim N(0, \sigma^2 I)\)) 을 풀어 \(\varepsilon_i \sim N(0, \sigma^2 \Omega_i)\) 로 일반화하면 \(V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 \Omega_i\) (식 7.3) — 피험자 간 이질성 (random effect) 과 시간적 자기상관 (Ω) 의 두 성분이 결합된다. 다섯 자기상관 구조 — AR(1) (식 7.7, 지수 감쇠, 2 모수), MA(1) (식 7.16, lag-1만, 2 모수), ARMA(1,1) (식 7.17, lag-1 강+후속 감쇠, 3 모수), Toeplitz (lag 별 독립, \(n\) 모수), NS-AR(1) (식 7.18-7.19, 비정상 분산 증가, 2 모수) — 가 자유 모수와 패턴의 trade-off 를 제공한다. EB estimator (식 7.5) 와 posterior covariance (식 7.6) 는 표준 MRM 의 식에서 단위 행렬 \(I\) 를 \(\Omega\) 로 바꾼 형태 로 깔끔히 일반화된다.
| 구조 | 식 | 자유 모수 | 적합 시점 |
|---|---|---|---|
| AR(1) | (7.7) | 2 | 인접 강한 상관 + 지수 감쇠 |
| MA(1) | (7.16) | 2 | lag-1 만 상관 |
| ARMA(1,1) | (7.17) | 3 | lag-1 강 + 후속 감쇠 |
| Toeplitz | — | \(n\) | 비단조·주기적 |
| NS-AR(1) | (7.18-7.19) | 2 | 분산 시간 따라 증가 |
1. AR(1) 가 default: 자유 모수 2 개, 해석 명확. 시계열·계량경제 표준.
2. 잔차 ACF/PACF 로 후보 선정: 표준 MRM 적합 후 잔차의 자기상관 패턴 검토 → AR / MA / ARMA / Toeplitz 중 선택.
3. AIC vs BIC 차이 인식: BIC 가 Toeplitz 같은 모수 많은 모형을 더 강하게 페널라이즈. 보고 시 둘 다 제시.
4. Fixed effects 의 강건성 점검: 자기상관 구조 변경 시 \(\hat{\beta}\) 도 변하는지 확인. 큰 변화 → 자기상관 모형이 추정에 결정적 → 신중한 선택.
5. CPM 과의 동치 인식: Toeplitz, AR(1) 등은 CPM 으로도 표현 가능. MRM-AC 표현이 random effect 분리 해석 유리, CPM 표현이 분산 직접 모형화 명확. 같은 모형의 두 시각.
12 관련 주제
Hedeker 시리즈 (Ch.4-7)
- § 30 — 혼합효과 회귀모형 개요 — Ch.4 도입
- § 32 — 랜덤 절편 MRM — Ch.4 random intercept
- § 33 — 랜덤 절편·추세 MRM — Ch.5
- § 34 — 행렬 공식화 — Ch.6 matrix formulation
- § 35 — MRM 추정론 — ML/REML/EM/Fisher Scoring
관련 개념
- Ch.6 — 공분산 패턴 모형 (CPM, 미작성) — 본 포스트의 한 측면 (시간 구조) 만 사용
- Ch.8 — GEE — Marginal 접근의 또 다른 길
- Generalized Estimating Equations
시계열 / 계량경제 응용
13 참고 문헌
- Hedeker, D., & Gibbons, R. D. (2006). Longitudinal Data Analysis. Wiley. § 7.2.
- Chi, E. M., & Reinsel, G. C. (1989). Models for longitudinal data with random effects and AR(1) errors. JASA, 84(406), 452-459. (MRM + AR(1) 의 핵심 reference)
- Mansour, H., Nordheim, E. V., & Rutledge, J. J. (1985). Maximum likelihood estimation of variance components in repeated measures designs assuming autoregressive errors. Biometrics, 41(1), 287-294. (NS-AR(1) 의 원전)
- Hedeker, D., & Gibbons, R. D. (1996b). MIXREG: A computer program for mixed-effects regression analysis with autocorrelated errors. Computer Methods and Programs in Biomedicine, 49(3), 229-252. (NS-AR(1) Cholesky 구현)
- Jones, R. H., & Boadi-Boateng, F. (1991). Unequally spaced longitudinal data with AR(1) serial correlation. Biometrics, 47(1), 161-175. (불균등 시점 확장)
- Verbeke, G., & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer. (대안 reference)