1 왜 행렬 표현이 필요한가
앞선 포스트(31~33번)에서 MRM을 스칼라 수식으로 전개했다. 피험자 \(i\)의 \(j\)번째 관측값을 하나씩 쓰는 방식이다.
\[ y_{ij} = \beta_0 + \beta_1 t_{ij} + v_{0i} + v_{1i} t_{ij} + \varepsilon_{ij} \]
이 표현은 직관적이지만, 통계적 성질을 논할 때 불편하다.
- 우도함수 \(L(\beta, \Sigma_v, \sigma^2)\)는 개인 \(i\)의 모든 관측값 \(y_{i1}, \ldots, y_{in_i}\)의 결합 분포에서 나온다.
- 개인별 잔차 \(\varepsilon_{ij}\)들이 독립이어도, \(y_{ij}\)들은 \(v_{0i}\)와 \(v_{1i}\)를 통해 상관된다.
- 불완전 데이터(개인마다 다른 \(n_i\))를 일반적으로 처리하려면 개인 수준 행렬이 필요하다.
행렬 표현은 이 모든 것을 하나의 방정식으로 정리한다. 수식이 간결해질 뿐 아니라, 추정·검정·예측의 도출 과정이 명확해진다.
2 핵심 행렬 방정식
\[ \underbrace{y_i}_{n_i \times 1} = \underbrace{X_i}_{n_i \times p} \underbrace{\beta}_{p \times 1} + \underbrace{Z_i}_{n_i \times r} \underbrace{v_i}_{r \times 1} + \underbrace{\varepsilon_i}_{n_i \times 1} \tag{4.8} \]
\(i = 1, \ldots, N\) (피험자), \(j = 1, \ldots, n_i\) (관측값)
각 기호의 의미와 차원을 정리한다.
| 기호 | 차원 | 의미 |
|---|---|---|
| \(y_i\) | \(n_i \times 1\) | 피험자 \(i\)의 반응변수 벡터 |
| \(X_i\) | \(n_i \times p\) | 고정 효과 설계 행렬 |
| \(\beta\) | \(p \times 1\) | 고정 효과 모수 벡터 (모집단 공통) |
| \(Z_i\) | \(n_i \times r\) | 랜덤 효과 설계 행렬 |
| \(v_i\) | \(r \times 1\) | 개인별 랜덤 효과 벡터 |
| \(\varepsilon_i\) | \(n_i \times 1\) | 잔차 벡터 (within-subject 오차) |
\(p\): 고정 효과 수 (절편 + 기울기 + 공변수), \(r\): 랜덤 효과 수 (일반적으로 \(r \leq p\))
3 설계 행렬 해부
3.1 랜덤 절편·추세 모형
가장 기본적인 경우부터 시작한다. 절편과 기울기 모두 랜덤이면:
\[ y_i = \begin{pmatrix} y_{i1} \\ y_{i2} \\ \vdots \\ y_{in_i} \end{pmatrix}, \quad X_i = Z_i = \begin{pmatrix} 1 & t_{i1} \\ 1 & t_{i2} \\ \vdots & \vdots \\ 1 & t_{in_i} \end{pmatrix} \]
\[ \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix}, \quad v_i = \begin{pmatrix} v_{0i} \\ v_{1i} \end{pmatrix} \]
이 경우 \(X_i = Z_i\)다. 고정 효과에 포함된 모든 항이 랜덤 효과에도 대응된다.
직관: \(X_i \beta\)는 모집단 평균 궤적을 결정하고, \(Z_i v_i\)는 개인이 그 평균에서 얼마나 벗어나는지를 결정한다.
3.2 공변수를 포함한 모형
진단 집단 변수 \(DX_i\)와 \(DX_i \times t_{ij}\) 교호작용을 추가하면, \(X_i\)는 4열로 확장되지만 \(Z_i\)는 그대로다:
\[ X_i = \begin{pmatrix} 1 & t_{i1} & DX_i & DX_i \times t_{i1} \\ 1 & t_{i2} & DX_i & DX_i \times t_{i2} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & t_{in_i} & DX_i & DX_i \times t_{in_i} \end{pmatrix}, \quad Z_i = \begin{pmatrix} 1 & t_{i1} \\ 1 & t_{i2} \\ \vdots & \vdots \\ 1 & t_{in_i} \end{pmatrix} \]
\[ \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}, \quad v_i = \begin{pmatrix} v_{0i} \\ v_{1i} \end{pmatrix} \]
핵심 통찰: 행렬 표현에서 \(X_i\)는 수준-1 공변수(시변)와 수준-2 공변수(개인 수준)를 구분하지 않는다. \(DX_i\)는 모든 행에 걸쳐 같은 값이지만 \(X_i\)의 한 열로 들어간다. 다층(multilevel) 표현이 모형을 수준-1/수준-2로 분리해 설명하지만, 실제로는 단 하나의 모형이 추정된다.
4 분포 가정
\[ \varepsilon_i \sim N(0, \sigma^2 I_{n_i}) \]
\[ v_i \sim N(0, \Sigma_v) \]
잔차와 랜덤 효과는 서로 독립이다: \(v_i \perp \varepsilon_i\).
\(\Sigma_v\)는 랜덤 효과들의 분산-공분산 행렬이다. 랜덤 절편·추세 모형에서:
\[ \Sigma_v = \begin{pmatrix} \sigma^2_{v0} & \sigma_{v0v1} \\ \sigma_{v0v1} & \sigma^2_{v1} \end{pmatrix} \]
5 결합 다변량 정규분포
\(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 I_{n_i} & Z_i \Sigma_v \\ \Sigma_v Z_i' & \Sigma_v \end{pmatrix} \right) \tag{4.9} \]
이 결합 분포에서 각 블록의 의미를 읽어내면:
| 블록 | 위치 | 의미 |
|---|---|---|
| \(Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i}\) | 좌상 | \(y_i\)의 주변 공분산 행렬 |
| \(Z_i \Sigma_v\) | 우상 | \(y_i\)와 \(v_i\) 사이의 교차 공분산 |
| \(\Sigma_v Z_i'\) | 좌하 | 위의 전치 |
| \(\Sigma_v\) | 우하 | \(v_i\) 자체의 사전 공분산 |
도출 과정:
\(y_i = X_i \beta + Z_i v_i + \varepsilon_i\)에서:
\[ \text{Var}(y_i) = Z_i \text{Var}(v_i) Z_i' + \text{Var}(\varepsilon_i) = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i} \]
\[ \text{Cov}(y_i, v_i) = \text{Cov}(Z_i v_i + \varepsilon_i, v_i) = Z_i \text{Var}(v_i) = Z_i \Sigma_v \]
6 경험 베이즈 추정 도출
결합 정규분포 (4.9)에서 \(y_i\)를 관측한 후 \(v_i\)의 사후 분포를 구하면:
\[ v_i \mid y_i \sim N(\hat{v}_i,\; \Sigma_{v|y_i}) \]
사후 평균 \(\hat{v}_i\) (경험 베이즈 추정량):
\[ \hat{v}_i = \left[ Z_i' (\sigma^2 I_{n_i})^{-1} Z_i + \Sigma_v^{-1} \right]^{-1} Z_i' (\sigma^2 I_{n_i})^{-1} (y_i - X_i \hat{\beta}) \]
동치 형식:
\[ \hat{v}_i = \Sigma_v Z_i' V_i^{-1} (y_i - X_i \hat{\beta}) \]
사후 공분산:
\[ \Sigma_{v|y_i} = \left[ Z_i' (\sigma^2 I_{n_i})^{-1} Z_i + \Sigma_v^{-1} \right]^{-1} \tag{4.11} \]
6.1 공식 해부
\[ \hat{v}_i = \underbrace{\left[ \underbrace{Z_i' (\sigma^2 I)^{-1} Z_i}_{\text{데이터 정보}} + \underbrace{\Sigma_v^{-1}}_{\text{사전 정밀도}} \right]^{-1}}_{\text{정밀도 역수}} \underbrace{Z_i' (\sigma^2 I)^{-1}}_{\text{데이터 가중}} \underbrace{(y_i - X_i \hat{\beta})}_{\text{잔차}} \]
각 항의 역할:
- \(Z_i' (\sigma^2 I)^{-1} Z_i\): 개인 \(i\)의 데이터에서 랜덤 효과를 추정하는 정보량. 관측 수 \(n_i\)가 클수록, 잔차 분산 \(\sigma^2\)가 작을수록 커진다.
- \(\Sigma_v^{-1}\): 사전 정밀도. 모집단 분포 \(N(0, \Sigma_v)\)를 “당기는 힘”이다.
- \((y_i - X_i \hat{\beta})\): 고정 효과 평균을 뺀 잔차 벡터. 이 안에 개인 편차 정보가 담겨 있다.
6.2 Shrinkage: 데이터와 사전 정보의 타협
\(\hat{v}_i\)는 두 극단의 중간 어딘가다:
- 완전 데이터 극단 (\(n_i \to \infty\)): 데이터 정보 >> 사전 정보 → OLS로 추정한 개인 효과
- 데이터 없는 극단 (\(n_i = 0\)): 사전 정보만 → \(\hat{v}_i = 0\) (모집단 평균)
\[ \hat{v}_i = (\text{신뢰 가중치}) \times \text{OLS 개인 추정} + (1 - \text{신뢰 가중치}) \times 0 \]
관측이 적을수록 모집단 평균 쪽으로 “수축(shrink)”된다. 이 수축이 극단적인 개인 추정을 안정화하고, 불완전 데이터에서도 합리적인 예측을 가능하게 한다.
7 주변 공분산 구조 V(y_i)
\[ V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i} \]
7.1 수치 전개 예시: 3시점, 랜덤 절편·추세
\(n_i = 3\), 시점 \(t = 0, 1, 2\)이면:
\[ Z_i = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{pmatrix} \]
\[ Z_i \Sigma_v Z_i' = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} \sigma^2_{v0} & \sigma_{v0v1} \\ \sigma_{v0v1} & \sigma^2_{v1} \end{pmatrix} \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 2 \end{pmatrix} \]
전개하면:
\[ = \begin{pmatrix} \sigma^2_{v0} & \sigma^2_{v0} + \sigma_{v0v1} & \sigma^2_{v0} + 2\sigma_{v0v1} \\ \sigma^2_{v0} + \sigma_{v0v1} & \sigma^2_{v0} + 2\sigma_{v0v1} + \sigma^2_{v1} & \sigma^2_{v0} + 3\sigma_{v0v1} + 2\sigma^2_{v1} \\ \sigma^2_{v0} + 2\sigma_{v0v1} & \sigma^2_{v0} + 3\sigma_{v0v1} + 2\sigma^2_{v1} & \sigma^2_{v0} + 4\sigma_{v0v1} + 4\sigma^2_{v1} \end{pmatrix} \]
\(\sigma^2 I_3\)를 더하면 대각 원소에만 \(\sigma^2\)가 추가된다.
완성된 \(V(y_i)\):
\[ V(y_i) = \begin{pmatrix} \sigma^2_{v0} + \sigma^2 & \sigma^2_{v0} + \sigma_{v0v1} & \sigma^2_{v0} + 2\sigma_{v0v1} \\ \sigma^2_{v0} + \sigma_{v0v1} & \sigma^2_{v0} + 2\sigma_{v0v1} + \sigma^2_{v1} + \sigma^2 & \sigma^2_{v0} + 3\sigma_{v0v1} + 2\sigma^2_{v1} \\ \sigma^2_{v0} + 2\sigma_{v0v1} & \sigma^2_{v0} + 3\sigma_{v0v1} + 2\sigma^2_{v1} & \sigma^2_{v0} + 4\sigma_{v0v1} + 4\sigma^2_{v1} + \sigma^2 \end{pmatrix} \]
읽는 법:
- \((j,j)\) 원소: \(t = j-1\)에서의 분산 → 시점이 지날수록 증가 (if \(\sigma^2_{v1} > 0\))
- \((j,k)\) 원소 (\(j \ne k\)): 시점 \(j, k\)의 공분산 → 두 시점의 절대 거리와 위치에 의존
- 복합대칭과 달리 분산과 공분산이 모두 시점 의존적이다
8 Reisby 데이터: 공분산 구조 적합 (§4.5.1)
6시점 데이터(\(t = 0, 1, 2, 3, 4, 5\))에서 관측 공분산 행렬(하삼각)과 모형 추정 공분산 행렬을 비교한다.
8.1 관측 공분산 행렬
\[ V_{\text{obs}}(y) = \begin{pmatrix} 20.55 \\ 10.50 & 22.07 \\ 10.20 & 12.74 & 30.09 \\ 9.69 & 12.43 & 25.96 & 41.15 \\ 7.17 & 7.39 & 18.25 & 36.54 & 48.59 \\ 6.02 & 10.10 & 25.56 & 32.93 & 52.12 & 62.60 \end{pmatrix} \]
패턴을 읽으면: - 대각 원소(분산): Week 0에서 20.55 → Week 5에서 62.60으로 3배 이상 증가 - 비대각 원소(공분산): 대각선 근처 원소들이 크고 멀어질수록 작아지는 경향
이 패턴은 복합대칭(모든 분산 동일, 모든 공분산 동일)으로는 설명이 불가능하다.
8.2 모형 추정 공분산 행렬
Table 4.5 추정값: \(\hat{\sigma}^2_{v0} = 12.63\), \(\hat{\sigma}_{v0v1} = -1.42\), \(\hat{\sigma}^2_{v1} = 2.08\), \(\hat{\sigma}^2 = 12.22\)와
\[ Z' = \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 2 & 3 & 4 & 5 \end{pmatrix} \]
를 \(Z\hat{\Sigma}_v Z' + \hat{\sigma}^2 I_6\)에 대입하면:
\[ \hat{V}(y) = \begin{pmatrix} 24.85 \\ 11.21 & 24.08 \\ 9.79 & 12.52 & 27.48 \\ 8.37 & 13.18 & 18.00 & 35.03 \\ 6.95 & 13.84 & 20.73 & 27.63 & 46.74 \\ 5.53 & 14.50 & 23.47 & 32.44 & 41.41 & 62.60 \end{pmatrix} \]
8.3 적합 평가
21개 원소(하삼각)를 단 4개 모수 \((\hat{\sigma}^2_{v0}, \hat{\sigma}_{v0v1}, \hat{\sigma}^2_{v1}, \hat{\sigma}^2)\)로 추정했다.
| 특징 | 관측 | 모형 |
|---|---|---|
| 분산 증가 추세 | Week 0: 20.55 → Week 5: 62.60 | Week 0: 24.85 → Week 5: 62.60 |
| 대각 근처 공분산 크기 | 비교적 큼 | 비교적 큼 |
| 대각에서 먼 공분산 | 작아지는 경향 | 더 체계적으로 감소 |
모형은 분산의 시간 증가를 잘 포착하고, 대각 가까운 공분산도 합리적으로 추정한다. 단, 실제 데이터에 비해 패턴이 더 단조롭다 — 관측 행렬의 비규칙적 변동을 4개 모수로 완전히 재현할 수는 없다.
이것이 모형의 타협점이다: 모수가 적어 효율적이고 일반화에 강하지만, 세부 패턴을 모두 포착하지는 못한다.
9 공분산 구조의 유연성
\(V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i}\)가 포괄하는 구조들:
| 모형 | \(Z_i\) 구성 | \(\Sigma_v\) | 공분산 구조 |
|---|---|---|---|
| 랜덤 절편만 | \(\mathbf{1}_{n_i}\) (1벡터) | \([\sigma^2_{v0}]\) | 복합대칭 (CS) |
| 랜덤 절편·추세 | \([\mathbf{1}, \mathbf{t}]\) | \(2\times2\) | 시점 의존적 |
| 일반 \(r\)개 랜덤 효과 | \(r\)열 | \(r\times r\) | 점점 일반적 |
| + 자기상관 오차 | \(r\)열 | \(r\times r\) | \(\sigma^2 I \to \sigma^2 \Omega_i\) |
오른쪽으로 갈수록 유연하지만 모수가 늘어난다. ANOVA(복합대칭 강제)와 MANOVA(완전 비구조적)의 중간 어딘가에 MRM이 위치한다.
MRM은 이 두 극단 사이의 효율적인 절충점을 찾는다. 모수 수에 비해 설명하는 구조의 풍부함이 크다.
10 시변 공변수 확장 (§4.5.2)
지금까지의 공변수(\(t_{ij}\), \(DX_i\))는 시간에 따라 변하거나 개인에게 고정되어 있었다. 이제 시점마다 값이 달라지는 공변수 — 시변 공변수(time-varying covariates) — 를 다룬다.
10.1 Reisby 약물 농도 데이터
Reisby 연구에서 Week 2~5에 이미프라민(IMI)과 대사체(DMI) 혈중 농도를 측정했다. 이 약물 농도는 피험자마다 다르고 (\(2 \sim 312 \mu g/L\)), 시점마다도 다르다.
처리 방법: 1. 극단값 영향을 줄이기 위해 자연로그 변환: \(\ln \text{IMI}_{ij}\), \(\ln \text{DMI}_{ij}\) 2. 개인 간 비교를 위해 전체 평균 중심화: \(I_{ij} = \ln \text{IMI}_{ij} - \overline{\ln \text{IMI}}\), \(D_{ij} = \ln \text{DMI}_{ij} - \overline{\ln \text{DMI}}\)
10.2 수준-1 모형 (4.13)
\[ y_{ij} = b_{0i} + b_{1i} t_{ij} + b_{2i} I_{ij} + b_{3i} D_{ij} + \varepsilon_{ij} \]
10.3 수준-2 모형 (4.14)
\[ b_{0i} = \beta_0 + v_{0i}, \quad b_{1i} = \beta_1 + v_{1i}, \quad b_{2i} = \beta_2, \quad b_{3i} = \beta_3 \]
IMI, DMI 효과는 일단 고정으로만 추정한다 (랜덤 기울기 확장 가능).
행렬 표현에서 \(X_i\):
\[ X_i = \begin{pmatrix} 1 & t_{i,2} & I_{i,2} & D_{i,2} \\ 1 & t_{i,3} & I_{i,3} & D_{i,3} \\ 1 & t_{i,4} & I_{i,4} & D_{i,4} \\ 1 & t_{i,5} & I_{i,5} & D_{i,5} \end{pmatrix} \]
\(I_{ij}\)와 \(D_{ij}\)는 행마다 다른 값 — 시변 공변수는 \(X_i\)에 열로 포함될 뿐, 구조적으로 다르지 않다.
10.4 Table 4.9: 현재 약물 농도 → 현재 우울 점수
| 모수 | 추정값 | SE | \(z\) | \(p\) |
|---|---|---|---|---|
| \(\beta_0\) (절편) | 18.17 | 0.71 | 25.70 | < .001 |
| \(\beta_1\) (시간) | -2.03 | 0.28 | -7.15 | < .001 |
| \(\beta_2\) (\(\ln\) IMI) | 0.60 | 0.85 | 0.71 | .48 |
| \(\beta_3\) (\(\ln\) DMI) | -1.20 | 0.63 | -1.90 | .06 |
\(-2\ln L = 1502.5\)
약물 농도가 당시 HDRS 점수와 관련되는지 분석했다. 두 약물 모두 유의하지 않다.
10.5 Table 4.10: 약물 농도 → HDRS 변화량
\[ (y_{ij} - y_{i0}) = b_{0i} + b_{1i} t_{ij} + b_{2i} I_{ij} + b_{3i} D_{ij} + \varepsilon_{ij} \tag{4.15} \]
| 모수 | 추정값 | SE | \(z\) | \(p\) |
|---|---|---|---|---|
| \(\beta_0\) (절편) | -5.18 | 0.66 | -7.87 | < .001 |
| \(\beta_1\) (시간) | -1.97 | 0.29 | -6.90 | < .001 |
| \(\beta_2\) (\(\ln\) IMI) | 0.63 | 0.82 | 0.77 | ns |
| \(\beta_3\) (\(\ln\) DMI) | -1.97 | 0.60 | -3.26 | < .001 |
\(-2\ln L = 1498.8\)
DMI가 HDRS 변화량에 대해 강하게 유의하다 (\(p < .001\), \(\hat{\beta}_3 = -1.97\)). ln DMI가 1단위 증가할 때 HDRS 변화량이 약 1.97점 더 감소 — 더 많이 개선된다.
반면 IMI는 여전히 유의하지 않다. IMI와 DMI는 서로 상관(\(r \approx 0.18\)~\(0.23\))되므로, IMI 계수는 DMI를 통제한 후의 순수 효과를 나타낸다.
10.6 왜 현재 점수가 아닌 변화량이 중요한가
Table 4.11 상관표를 보면:
| Week 2 | Week 3 | Week 4 | Week 5 | |
|---|---|---|---|---|
| ln DMI vs HDRS 점수 | -0.177 | -0.075 | -0.246 | -0.293 |
| ln DMI vs HDRS 변화량 | -0.366 | -0.281 | -0.363 | -0.361 |
변화량과의 상관이 일관되게 크다. 약물 농도는 “지금 얼마나 우울한가”보다 “얼마나 나아졌는가”와 더 관련이 있다. 이는 약물 농도가 치료 반응의 생물학적 기전과 연결된다는 해석을 지지한다.
10.7 시변 공변수의 Within-Between 분해
시변 공변수 \(I_{ij}\)를 포함할 때 암묵적으로 가정하는 것이 있다: 개인 내(within-subject) 효과와 개인 간(between-subject) 효과가 같다는 것이다.
이를 확인하기 위해 분해한다:
\[ I_{ij} = \underbrace{\bar{I}_i}_{\text{개인 평균 (between)}} + \underbrace{(I_{ij} - \bar{I}_i)}_{\text{개인 내 편차 (within)}} \]
두 항을 \(X_i\)에 별도 열로 넣으면 within 효과와 between 효과를 독립적으로 추정할 수 있다. 만약 두 추정값이 유사하다면, 단일 계수 가정이 타당하다. 다르다면, 생태학적 오류(ecological fallacy)가 개입되었을 가능성이 있다.
11 구현 코드
11.1 Python (행렬 계산 직접 구현)
import numpy as np
# ============================================================
# 1. 행렬 설정 (n_i=6, 랜덤 절편·추세)
# ============================================================
t = np.array([0, 1, 2, 3, 4, 5])
n_i = len(t)
# 랜덤 효과 설계 행렬 Z_i
Z_i = np.column_stack([np.ones(n_i), t]) # (6,2)
# Reisby 추정값 (Table 4.5)
sigma2_v0 = 12.63
sigma_v0v1 = -1.42
sigma2_v1 = 2.08
sigma2 = 12.22
# Σ_v (랜덤 효과 공분산)
Sigma_v = np.array([[sigma2_v0, sigma_v0v1],
[sigma_v0v1, sigma2_v1]])
# V(y_i) = Z_i Σ_v Z_i' + σ² I
V_model = Z_i @ Sigma_v @ Z_i.T + sigma2 * np.eye(n_i)
print("모형 추정 공분산 행렬 V(y_i):")
print(np.round(V_model, 2))
# 대각: [24.85, 24.08, 27.48, 35.03, 46.74, 62.60]
# ============================================================
# 2. 경험 베이즈 추정 수식 확인
# ============================================================
# 특정 피험자의 랜덤 효과 추정
sigma2_I = sigma2 * np.eye(n_i)
# [Z'(σ²I)^{-1}Z + Σ_v^{-1}]^{-1}
precision_data = Z_i.T @ np.linalg.inv(sigma2_I) @ Z_i # 데이터 정밀도
precision_prior = np.linalg.inv(Sigma_v) # 사전 정밀도
posterior_cov = np.linalg.inv(precision_data + precision_prior) # Σ_{v|y}
print("\n사후 공분산 Σ_{v|y_i}:")
print(np.round(posterior_cov, 4))
# 가상 잔차 벡터로 EB 추정
beta_hat = np.array([23.58, -2.38]) # 고정 효과 추정값
X_i_simple = Z_i.copy()
y_i_example = np.array([28, 25, 20, 16, 12, 9]) # 예시 관측값
residual = y_i_example - X_i_simple @ beta_hat
v_hat = posterior_cov @ Z_i.T @ np.linalg.inv(sigma2_I) @ residual
print(f"\n경험 베이즈 랜덤 효과 추정:")
print(f"v̂₀ᵢ (절편 편차) = {v_hat[0]:.3f}")
print(f"v̂₁ᵢ (기울기 편차) = {v_hat[1]:.3f}")
print(f"개인 절편 b₀ᵢ = {beta_hat[0] + v_hat[0]:.3f}")
print(f"개인 기울기 b₁ᵢ = {beta_hat[1] + v_hat[1]:.3f}")
# ============================================================
# 3. 복합대칭 vs 랜덤 절편·추세 공분산 비교
# ============================================================
# 복합대칭 (랜덤 절편만)
sigma2_ri = 16.15 # Table 4.3 값
sigma2_e = 19.04
V_cs = sigma2_ri * np.ones((n_i, n_i)) + sigma2_e * np.eye(n_i)
print("\n복합대칭 모형 대각 (분산):")
print(np.round(np.diag(V_cs), 2)) # 모두 동일
print("\n랜덤 절편·추세 모형 대각 (분산):")
print(np.round(np.diag(V_model), 2)) # 시간에 따라 증가11.2 R (lme4 + 공분산 구조 시각화)
library(lme4)
library(Matrix)
library(ggplot2)
library(reshape2)
# ============================================================
# 1. 랜덤 절편·추세 모형 피팅
# ============================================================
reisby <- read.csv("reisby.csv")
fit <- lmer(hdrs ~ week + (1 + week | id), data = reisby, REML = FALSE)
summary(fit)
# 분산 모수 추출
vc <- VarCorr(fit)
Sigma_v <- as.matrix(vc$id) # 2×2 Σ_v
sigma2 <- sigma(fit)^2 # σ²
cat("Σ_v:\n"); print(round(Sigma_v, 2))
cat("σ²:", round(sigma2, 2), "\n")
# ============================================================
# 2. 모형 추정 공분산 행렬 V(y_i) 계산
# ============================================================
t_vals <- 0:5
Z <- cbind(1, t_vals) # 6×2
V_model <- Z %*% Sigma_v %*% t(Z) + sigma2 * diag(6)
rownames(V_model) <- colnames(V_model) <- paste0("Week", 0:5)
cat("\n모형 추정 공분산 행렬:\n")
print(round(V_model, 2))
# ============================================================
# 3. 관측 공분산 행렬
# ============================================================
reisby_wide <- reshape(reisby[c("id", "week", "hdrs")],
timevar = "week", idvar = "id", direction = "wide")
V_obs <- cov(reisby_wide[, -1], use = "pairwise.complete.obs")
cat("\n관측 공분산 행렬:\n")
print(round(V_obs, 2))
# ============================================================
# 4. 공분산 구조 시각화 (lag별 covariance plot)
# ============================================================
get_cov_by_lag <- function(V, name) {
n <- nrow(V)
result <- data.frame()
for (lag in 0:(n-1)) {
for (j in 1:(n-lag)) {
k <- j + lag
result <- rbind(result, data.frame(
lag = lag,
week = j - 1,
cov = V[j, k],
model = name
))
}
}
result
}
df_obs <- get_cov_by_lag(V_obs, "관측")
df_model <- get_cov_by_lag(V_model, "모형 추정")
df_all <- rbind(df_obs, df_model)
ggplot(df_all, aes(x = lag, y = cov, color = factor(week), linetype = model)) +
geom_line() +
geom_point() +
labs(
title = "Reisby 데이터: 관측 vs 모형 공분산 (Lag별)",
x = "Lag (주 차이)", y = "공분산",
color = "시작 Week", linetype = "출처"
) +
theme_minimal()
# ============================================================
# 5. 경험 베이즈 개인 추정
# ============================================================
re <- ranef(fit)$id
colnames(re) <- c("v0i", "v1i")
fe <- fixef(fit)
re$b0i <- fe["(Intercept)"] + re$v0i
re$b1i <- fe["week"] + re$v1i
re$id <- rownames(re)
# 사후 공분산 Σ_{v|y_i} 계산 (완전 데이터 피험자 기준)
ni <- 6 # 6시점
Z_mat <- cbind(1, 0:5)
sigma2_I_inv <- (1/sigma2) * diag(ni)
Sigma_v_inv <- solve(Sigma_v)
posterior_cov <- solve(t(Z_mat) %*% sigma2_I_inv %*% Z_mat + Sigma_v_inv)
cat("\n사후 공분산 Σ_{v|y_i} (완전 데이터 기준):\n")
print(round(posterior_cov, 4))
cat("사후 SD v0i:", round(sqrt(posterior_cov[1,1]), 3), "\n")
cat("사후 SD v1i:", round(sqrt(posterior_cov[2,2]), 3), "\n")
cat("사전 SD v0i:", round(sqrt(Sigma_v[1,1]), 3), "\n")
cat("사전 SD v1i:", round(sqrt(Sigma_v[2,2]), 3), "\n")
# 사후 SD < 사전 SD → 데이터가 불확실성을 줄임
# ============================================================
# 6. 시변 공변수 모형 (§4.5.2)
# ============================================================
# 약물 농도 데이터가 있을 경우
# reisby_drug: id, week, hdrs, log_imi, log_dmi (grand-mean centered)
# reisby_drug_sub <- subset(reisby_drug, week >= 2) # Week 2~5
# reisby_drug_sub$t_local <- reisby_drug_sub$week - 2 # 0~3 재코딩
# 모형 (4.13): 현재 HDRS ~ 약물 농도
# fit_drug <- lmer(
# hdrs ~ t_local + log_imi + log_dmi + (1 + t_local | id),
# data = reisby_drug_sub, REML = FALSE
# )
# 모형 (4.15): HDRS 변화량 ~ 약물 농도
# reisby_drug_sub$hdrs_change <- reisby_drug_sub$hdrs - reisby_drug_sub$baseline_hdrs
# fit_change <- lmer(
# hdrs_change ~ t_local + log_imi + log_dmi + (1 + t_local | id),
# data = reisby_drug_sub, REML = FALSE
# )12 CS vs MRM 공분산 구조 비교 요약
| 항목 | 복합대칭 (CS) | 랜덤 절편·추세 MRM |
|---|---|---|
| 분산 | 모든 시점 동일 | 시점 의존적 |
| 공분산 | 모든 시점 쌍 동일 | 시점 쌍 의존적 |
| 모수 수 | 2 (\(\sigma^2_{v0}\), \(\sigma^2\)) | 4 (\(\sigma^2_{v0}\), \(\sigma_{v0v1}\), \(\sigma^2_{v1}\), \(\sigma^2\)) |
| 적용 가능한 구조 | 단 하나 | 폭넓은 범위 |
| Reisby 적합도 | 나쁨 (LRT \(\chi^2=66.15\)) | 양호 |
13 정리
MRM의 행렬 표현 \(y_i = X_i \beta + Z_i v_i + \varepsilon_i\)는 단순한 표기 편의가 아니다.
- \(X_i\)와 \(Z_i\)의 분리는 고정/랜덤 효과의 역할을 명확히 하고, 공변수 확장 방법을 체계화한다.
- 결합 분포 (4.9)는 경험 베이즈 추정 (4.10)과 사후 공분산 (4.11)을 단계적으로 도출하는 기반이다.
- 공분산 행렬 (4.12) \(V(y_i) = Z_i \Sigma_v Z_i' + \sigma^2 I_{n_i}\)는 복합대칭부터 시점 의존적 구조까지 포괄하며, 모수 수 대비 높은 표현력을 제공한다.
- 시변 공변수는 \(X_i\)에 열로 포함되며, within-between 분해로 맥락 효과를 추가로 탐색할 수 있다.
이 행렬 공식화는 Ch.5(ML/REML 추정)와 Ch.6~7(공분산 구조 선택)의 논의를 위한 통합 언어다.
14 관련 주제
선행 지식
후속 주제
- ML/REML 추정 — 우도 최대화, 분산 추정 편향 보정
- 공분산 구조 선택 — AR(1)·비구조적·정보기준 (Ch.6 CPM, 미작성)