1 들어가며

§ 5.2 곡선형 추세 MRM 에서 이차 추세 모형의 정의와 적합을 다뤘다. 이때 시간 변수 \(t\) 와 그 제곱 \(t^2\) 을 직접 회귀자로 사용했다.

이 단순한 표현에는 두 가지 실무적 문제가 있다.

다공선성 (collinearity): \(t\) 와 \(t^2\) 의 표본 상관이 1 에 가까워 추정량의 표준오차가 부풀려지고 수치적으로 불안정해진다.
척도 차이 (scale heterogeneity): 차수가 높아질수록 회귀자의 분산이 폭발적으로 커져서, 회귀 계수의 절대값이 차수 간 직접 비교 불가능하다.

직교 다항식 (orthogonal polynomial) 은 이 두 문제를 동시에 해결한다. 같은 모형의 재표현 (reparameterization) 일 뿐이지만, 수치적 안정성과 해석 편의가 크게 향상된다.

이 sub-post 는 다음을 다룬다.

다공선성 문제의 수치적 시연
시간 중심화 (\(t - \bar{t}\)) 의 부분 해법과 한계
Bock (1975) 의 Cholesky 인수분해 절차 (5 단계)
원본 척도와 직교 척도의 모형 표현 비교
모수 변환 공식 (\(\beta \leftrightarrow \gamma\), \(\Sigma_\upsilon \leftrightarrow \Sigma_\theta\))
표준오차 변환에 사용되는 \(G^+\) · vec/vech 연산
Reisby 데이터의 직교 적합 결과 (Table 5.2) 와 절편 의미 변화

2 다공선성 문제 — 수치 시연

원본 척도에서 \(t\) 와 \(t^2\) 의 상관이 얼마나 큰지 직접 확인한다.

import numpy as np

# 6 시점, 등간격
t = np.arange(6)  # 0, 1, 2, 3, 4, 5
t2 = t**2          # 0, 1, 4, 9, 16, 25

corr = np.corrcoef(t, t2)[0, 1]
print(f"corr(t, t^2) = {corr:.4f}")
# 출력: corr(t, t^2) = 0.9747

# 3 시점만 보면 더 극단적
t3 = np.arange(3)
print(f"corr(t, t^2) [3 timepoints] = {np.corrcoef(t3, t3**2)[0, 1]:.4f}")
# 출력: 약 0.98

표본 상관이 약 0.97 — 거의 완벽한 양의 상관이다.

다공선성이 만드는 문제

회귀에서 두 회귀자가 거의 완전 상관이면:

\(X^\top X\) 의 조건수가 폭발 → 행렬 역산이 수치적으로 불안정
회귀 계수의 표준오차 (\(\sqrt{(X^\top X)^{-1}}\) 의 대각) 가 부풀려짐
추정값 \(\hat\beta_1, \hat\beta_2\) 사이에 강한 음의 상관 → 한 모수가 0 에 가까울 때 다른 모수가 비정상적으로 크게 추정될 수 있음

이는 단순히 “고급 기법” 의 문제가 아니라, 이차 추세 적합이 실패하거나 수렴하지 않는 직접적 원인이다.

2.1 시간 중심화 (\(t - \bar{t}\)) — 부분 해법

시간을 중심화하면 다공선성이 줄어든다.

t_centered = t - t.mean()
print(f"corr(t_c, t_c^2) = {np.corrcoef(t_centered, t_centered**2)[0, 1]:.4f}")
# 출력: 0.0 (등간격·균형 데이터에서 정확히 0)

중심화의 효과와 한계

효과: 등간격이고 모든 시점에서 같은 수의 관측이 있으면, 중심화된 시간과 그 제곱의 상관이 정확히 0 이 된다. 불등간격이거나 결측이 있어도 상관이 크게 줄어든다.

한계: 중심화는 상관만 0 으로 만들 뿐이다. 두 회귀자의 분산은 여전히 다르다 — \(t_c\) 의 분산보다 \(t_c^2\) 의 분산이 훨씬 작다 (또는 큼). 회귀 계수의 절대값이 차수 간 비교 불가능한 상태는 그대로 남는다.

부수 효과: 절편의 의미가 바뀐다. 원본에서 \(\beta_0\) = “첫 시점 (\(t = 0\)) 의 평균” 이지만, 중심화 후에는 “시간 중간점 (\(t = \bar{t}\)) 의 평균” 이 된다.

직교 다항식은 중심화를 한 단계 더 밀고 나간 것 — 상관 0 + 단위 분산 까지 만든다.

3 Cholesky 인수분해 — 일반 절차

직교 다항식 표는 등간격일 때만 직접 사용 가능 (Pearson-Hartley 1976). 임의의 시간 구조 (불등간격, 결측 포함) 에서 직교 다항식을 얻으려면 Bock (1975) 의 Cholesky 인수분해 절차를 따른다.

직교 다항식 행렬 계산 — 5 단계

입력: 시간 행렬 \(T\), 즉 \(T = [\mathbf{1}, t, t^2, \ldots, t^{p-1}]\), \(n \times p\)

1 단계: \(T^\top T\) 계산. 이는 \(p \times p\) 대칭 행렬이며, 다항식 항들의 “유사 공분산 행렬” 역할을 한다.

2 단계: Cholesky 인수분해 \[ T^\top T = SS^\top \] \(S\) 는 하삼각, \(S^\top\) 는 그 전치 (상삼각).

3 단계: \(S^\top\) 의 역행렬 \((S^\top)^{-1}\) 계산.

4 단계: 직교 다항식 행렬 \[ T_{\text{ortho}} = T (S^\top)^{-1} \]

5 단계: 검증 — \(T_{\text{ortho}}^\top T_{\text{ortho}} = I_p\) 가 성립.

3.1 직관: 왜 Cholesky 가 직교화를 만드는가

이 절차는 사실 whitening 과 같은 원리다.

분산 행렬이 \(\Sigma\) 인 다변량 정규를 \(\Sigma^{-1/2}\) 로 변환하면 단위 분산이 된다. Cholesky 분해 \(\Sigma = LL^\top\) 에서 \(L^{-1}\) 이 그 역할을 한다.

여기서 다항식 항들의 “공분산” 역할을 하는 것이 \(T^\top T\) 다. \(T^\top T = SS^\top\) 의 \(S^{-\top}\) 가 그 역할 — \(T(S^\top)^{-1}\) 의 열들이 서로 직교하고 단위 분산을 갖는다.

한 줄 요약

직교 다항식 = “다항식 항들의 데이터 기반 whitening”. Cholesky 인수분해는 그 whitening 행렬을 효율적으로 계산하는 도구다.

3.2 6 시점 등간격 예시 — 수치 재현

Hedeker (2006, p.88) 의 6 시점 (\(t = 0, \ldots, 5\)) 이차 다항식 예시를 재현한다.

import numpy as np

# 시간 행렬: 절편, 선형, 이차
t = np.arange(6).astype(float)
T = np.column_stack([np.ones_like(t), t, t**2])

print("T =")
print(T)
# [[1, 0, 0],
#  [1, 1, 1],
#  [1, 2, 4],
#  [1, 3, 9],
#  [1, 4, 16],
#  [1, 5, 25]]

# 1 단계: T'T
TtT = T.T @ T
print("\nT'T =")
print(TtT)
# [[6, 15, 55],
#  [15, 55, 225],
#  [55, 225, 979]]

# 2 단계: Cholesky 분해 (S 하삼각)
S = np.linalg.cholesky(TtT)
Sp = S.T  # 상삼각
print("\nS' =")
print(np.round(Sp, 4))
# [[2.4495, 6.1237, 22.4537],
#  [0,      4.1833, 20.9165],
#  [0,      0,       6.1101]]

# 3 단계: (S')^-1
Sp_inv = np.linalg.inv(Sp)
print("\n(S')^-1 =")
print(np.round(Sp_inv, 4))
# [[0.4082, -0.5976, 0.5455],
#  [0,      0.2390, -0.8183],
#  [0,      0,      0.1637]]

# 4 단계: 직교 다항식 행렬
T_ortho = T @ Sp_inv
print("\nT_ortho =")
print(np.round(T_ortho, 4))

# 5 단계: 검증 — 단위 행렬에 가까운가
print("\nT_ortho' T_ortho =")
print(np.round(T_ortho.T @ T_ortho, 4))
# 단위 행렬 I_3 에 매우 가까움

기대 결과:

\(T_{\text{ortho}}\) 의 첫 열은 모두 \(1/\sqrt{6} \approx 0.4082\) (절편 표준화)
둘째 열은 \(\{-5, -3, -1, 1, 3, 5\}/\sqrt{70}\) (선형 직교)
셋째 열은 \(\{5, -1, -4, -4, -1, 5\}/\sqrt{84}\) (이차 직교)

이는 Pearson-Hartley 표 (Hedeker p.87) 와 정확히 일치한다.

4 모형 표현 — 두 좌표계 비교

원본 시간 척도와 직교 다항식 척도는 같은 모형의 두 좌표 표현 이다.

	원본 척도 (식 5.4)	직교 척도 (식 5.5)
디자인	\(X_i, Z_i\) (원본 시간)	\(X_i (S^\top)^{-1}, Z_i (S^\top)^{-1}\)
고정효과	\(\beta\)	\(\gamma\)
랜덤효과	\(\upsilon_i\)	\(\theta_i\)
모형	\(y_i = X_i \beta + Z_i \upsilon_i + \varepsilon_i\)	\(y_i = X_i (S^\top)^{-1} \gamma + Z_i (S^\top)^{-1} \theta_i + \varepsilon_i\)
마진 분산	\(V(y_i) = Z_i \Sigma_\upsilon Z_i^\top + \sigma^2 I\)	\(V(y_i) = Z_i (S^\top)^{-1} \Sigma_\theta S^{-1} Z_i^\top + \sigma^2 I\)

핵심: 두 모형은 수치적으로 동일

같은 데이터에 두 모형을 적합하면 로그 우도 \(-2 \log L\) 가 정확히 같다 (Hedeker Table 5.1·5.2 모두 2207.64). 표본의 마진 분포 \(f(y_i; \theta)\) 가 같기 때문이다 — 모수만 좌표 변환되었을 뿐 모형은 동일하다.

선택은 해석 편의의 문제이지, 모형 적합도의 문제가 아니다.

5 모수 변환 공식

두 좌표계의 모수는 일대일 대응한다.

5.1 회귀 계수 (고정 효과)

\[ \beta = (S^\top)^{-1} \gamma \qquad\Longleftrightarrow\qquad \gamma = S^\top \beta \tag{5.14} \]

랜덤 효과 (개인별 이탈) 도 같은 변환:

\[ \upsilon_i = (S^\top)^{-1} \theta_i \qquad\Longleftrightarrow\qquad \theta_i = S^\top \upsilon_i \tag{5.15} \]

5.2 분산-공분산 행렬

랜덤 효과 분산 행렬은 양쪽에서 곱하는 변환을 받는다.

\[ \Sigma_\upsilon = (S^\top)^{-1} \Sigma_\theta S^{-1} \qquad\Longleftrightarrow\qquad \Sigma_\theta = S^\top \Sigma_\upsilon S \tag{5.16} \]

직관: 왜 양쪽에서 곱하는가

분산 행렬은 이차 형식 (quadratic form) 의 계수다. \(\text{Var}(A x) = A \, \text{Var}(x) \, A^\top\) 같은 일반 공식에서, 변환 행렬 \(A\) 가 양쪽에 한 번씩 들어간다.

\(\upsilon = (S^\top)^{-1} \theta\) 이므로 \(A = (S^\top)^{-1}\) 을 대입하면: \[ \Sigma_\upsilon = (S^\top)^{-1} \Sigma_\theta \big[(S^\top)^{-1}\big]^\top = (S^\top)^{-1} \Sigma_\theta S^{-1} \]

5.3 회귀 계수의 표준오차 변환

추정량 \(\hat\beta\) 의 분산 행렬은 같은 변환을 따른다.

\[ V(\hat\beta) = (S^\top)^{-1} V(\hat\gamma) S^{-1} \qquad\Longleftrightarrow\qquad V(\hat\gamma) = S^\top V(\hat\beta) S \tag{5.18-19} \]

표준오차는 대각 성분의 제곱근.

# Hedeker Table 5.2 의 V(gamma_hat)
V_gamma = np.array([
    [1.8708, 0.5823, -0.1402],
    [0.5823, 0.7470,  0.0214],
    [-0.1402, 0.0214, 0.2914],
])

# 원본 척도의 V(beta_hat)
V_beta = Sp_inv @ V_gamma @ Sp_inv.T
print("V(beta_hat) =")
print(np.round(V_beta, 4))
# 대각 성분 제곱근 → SE
print("\nSE(beta) =", np.round(np.sqrt(np.diag(V_beta)), 4))
# 출력: [0.5521, 0.4790, 0.0883] — Table 5.1 의 SE 와 일치

5.4 분산 행렬 추정량의 표준오차 변환

\(\widehat\Sigma_\theta\) 같은 행렬형 추정량의 표준오차 변환은 한 단계 더 복잡하다. McCulloch (1982) 의 vec/vech 연산을 사용한다.

vec 와 vech — 행렬을 벡터로 펼치기

\(r \times r\) 행렬 \(\Sigma\) 에 대해:

\(\text{vec}\Sigma\): 모든 원소를 열 방향으로 쌓아 \(r^2 \times 1\) 벡터
\(\text{vech}\Sigma\): 대각 + 상삼각 (또는 하삼각) 의 고유 원소만 쌓아 \(r^* \times 1\) 벡터, \(r^* = r(r+1)/2\)

대칭 행렬은 \(r^* = r(r+1)/2\) 개의 고유 모수만 가지므로 vech 가 더 효율적이다.

예: \(r = 3\) 이면 \(r^2 = 9\), \(r^* = 6\).

vec 와 vech 사이에는 변환 행렬 \(G\) 와 그 Moore-Penrose 의사역 \(G^+\) 가 존재한다.

\[ \text{vec} \, \Sigma = G \, \text{vech} \, \Sigma \qquad \text{vech} \, \Sigma = G^+ \, \text{vec} \, \Sigma \]

이를 이용하면 분산 행렬 추정량의 표준오차 변환은 다음과 같다 (Hedeker 식 5.23-24).

\[ V(\text{vech}\,\widehat\Sigma_\theta) = G^+ (S^\top \otimes S^\top) G \cdot V(\text{vech}\,\widehat\Sigma_\upsilon) \cdot G^\top (S \otimes S) G^{+\top} \tag{5.23} \]

여기서 \(\otimes\) 는 Kronecker 곱.

직관: 벡터화·Kronecker 곱이 왜 필요한가

행렬을 벡터로 보는 이유: 분산 추정량의 분산을 다루려면 “분산의 분산” 같은 4 차원 텐서가 필요하다. Vec 으로 펼치면 이 4 차원 구조를 평범한 행렬로 다룰 수 있다.

Kronecker 곱: \(\Sigma_\theta = S^\top \Sigma_\upsilon S\) 같은 행렬 변환을 vec 형태로 쓰면 \(\text{vec}\Sigma_\theta = (S^\top \otimes S^\top) \text{vec} \Sigma_\upsilon\) 이 되기 때문.

vech ↔︎ vec 변환이 추가되는 이유: 대칭 행렬의 고유 모수만 다루기 위함. \(G^+\) 가 중복 원소를 제거.

이 모든 표기는 표준오차 한 줄 변환을 위한 도구일 뿐이다 — 통계 패키지가 자동 처리한다.

6 Reisby 직교 분석 — Hedeker Table 5.2

§ 5.2 에서 본 같은 데이터를 직교 다항식으로 재적합한 결과:

Table 5.2 — 이차 추세 MRM (직교 다항식 척도)

모수	추정값	SE	\(Z\)	\(p <\)
\(\gamma_0\) (절편)	43.24	1.37	31.61	.0001
\(\gamma_1\) (선형)	-9.94	0.86	-11.50	.0001
\(\gamma_2\) (이차)	0.31	0.54	0.58	.56
\(\sigma_{\theta_0}^2\)	111.91	21.60
\(\sigma_{\theta_0 \theta_1}\)	37.99	10.92
\(\sigma_{\theta_1}^2\)	37.04	8.90
\(\sigma_{\theta_0 \theta_2}\)	-10.14	6.19
\(\sigma_{\theta_1 \theta_2}\)	-0.82	3.80
\(\sigma_{\theta_2}^2\)	7.23	3.50
\(\sigma^2\)	10.52	1.11

Note: \(-2 \log L = 2207.64\) — Table 5.1 과 동일. (Hedeker, 2006, p.90)

6.1 두 표 비교 — 무엇이 같고 무엇이 다른가

동일한 것

로그 우도 \(-2 \log L = 2207.64\)
오차 분산 \(\widehat\sigma^2 = 10.52\)
곡선의 모양 — 두 모수 집합은 정확한 일대일 변환

6.2 다른 것

회귀 계수의 절대값과 SE
절편의 해석 의미
랜덤 효과 분산의 절대값
랜덤 효과 간 상관의 부호

6.3 절편 의미 변화 — 가장 중요한 해석 차이

직교 변환에서 가장 헷갈리는 지점

원본 척도의 \(\beta_0 = 23.76\) → 첫 시점 (\(t = 0\), baseline) 의 평균 반응

직교 척도의 \(\gamma_0 = 43.24\) → 모든 시점에 걸친 평균 반응 (전체 평균)

같은 데이터에서 절편이 23.76 vs 43.24 로 크게 다른 것은 정상이다.

검증: \(\gamma_0 = S^\top_{(1,\cdot)} \beta = 2.4495 \times 23.76 + 6.1237 \times (-2.63) + 22.4537 \times 0.05 \approx 43.24\)

6.4 랜덤 효과 상관 부호 변화

원본 척도와 직교 척도에서 랜덤 절편-선형의 상관 부호가 뒤집힌다.

	\(\widehat{\sigma}_{\cdot 0 \cdot 1}\)	표준화 상관
원본 (\(\upsilon\))	-0.92	\(-0.11\)
직교 (\(\theta\))	37.99	\(+0.59\)

왜 부호가 뒤집히는가

원본 척도에서 \(\upsilon_{0i}\) 는 “사람 \(i\) 의 baseline (첫 시점) 이탈” 이고, \(\upsilon_{1i}\) 는 “사람 \(i\) 의 선형 추세 이탈” 이다.

상관 = \(-0.11\) 의 의미: baseline 이 평균보다 높은 사람이 약간 더 음의 추세 (호전 빠름).

직교 척도에서 \(\theta_{0i}\) 는 “사람 \(i\) 의 평균 수준 (전체 시간 평균) 이탈” 이고, \(\theta_{1i}\) 는 “사람 \(i\) 의 선형 추세 이탈” 이다.

상관 = \(+0.59\) 의 의미: 평균 수준이 평균보다 높은 사람이 더 음의 추세 (호전 빠름).

수학적으로는 다음 관계 — baseline 이 높을수록 추세가 음 + 추세가 음일수록 평균 수준이 baseline 보다 낮음 → 평균 수준과 추세는 양의 상관 (둘 다 “초기 상태 + 회복 속도” 의 다른 측면).

같은 통계적 사실이 좌표계에 따라 다르게 보일 뿐이다. 어느 한쪽이 옳은 것이 아니라 두 표현이 동일한 모형의 다른 절단면이다.

6.5 분산 비율 — 직교 척도에서 더 직관적

직교 척도에서 랜덤 효과 분산은 같은 단위 (분산 1 의 회귀자) 라 직접 비교 가능하다.

\[ \frac{\widehat\sigma_{\theta_0}^2}{\widehat\sigma_{\theta_0}^2 + \widehat\sigma_{\theta_1}^2 + \widehat\sigma_{\theta_2}^2} = \frac{111.91}{111.91 + 37.04 + 7.23} = \frac{111.91}{156.18} \approx 71.7\% \]

세 분산의 상대 비율: 71.7% (절편) / 23.7% (선형) / 4.6% (이차)

분산 비율의 의미

이는 개인 이질성 (between-subject variation) 이 어떤 추세 차원에 분포하는지 보여준다.

71.7% — 절편: 사람마다 평균 우울 수준 (baseline 이라기보다 전체 평균) 이 다름
23.7% — 선형: 사람마다 호전 속도가 다름
4.6% — 이차: 사람마다 호전 가속/감속 패턴이 다름

원본 척도의 분산 (10.44, 6.64, 0.19) 으로는 이런 비교가 불가능하다 — 회귀자의 분산이 다르므로. 직교 변환의 가장 큰 해석 이득은 이 분산 비율의 직접 비교 다.

7 코드 예시 — 직교 변환과 검증

7.1 Step 1: Cholesky 절차 직접 구현 + 변환 검증

import numpy as np

# 6 시점, 이차까지
t = np.arange(6).astype(float)
T = np.column_stack([np.ones_like(t), t, t**2])

# Cholesky
S = np.linalg.cholesky(T.T @ T)
Sp = S.T
Sp_inv = np.linalg.inv(Sp)

# 직교 다항식 행렬
T_ortho = T @ Sp_inv

# 검증 1: 직교성
print("T_ortho' T_ortho =")
print(np.round(T_ortho.T @ T_ortho, 6))  # 단위 행렬

# 검증 2: 모수 변환 — Hedeker Table 5.2 의 gamma 를 Table 5.1 의 beta 로
gamma_hat = np.array([43.24, -9.94, 0.31])
beta_hat = Sp_inv @ gamma_hat
print(f"\nbeta_hat = {np.round(beta_hat, 2)}")
# 출력: [23.76, -2.63, 0.05] — Table 5.1 과 일치

# 검증 3: SE 변환
V_gamma = np.array([
    [1.8708, 0.5823, -0.1402],
    [0.5823, 0.7470,  0.0214],
    [-0.1402, 0.0214, 0.2914],
])
V_beta = Sp_inv @ V_gamma @ Sp_inv.T
SE_beta = np.sqrt(np.diag(V_beta))
print(f"\nSE(beta) = {np.round(SE_beta, 4)}")
# 출력: [0.5521, 0.4790, 0.0883] — Table 5.1 의 SE 와 일치

# 검증 4: 랜덤 효과 분산 변환
Sigma_theta = np.array([
    [111.91, 37.99, -10.14],
    [37.99, 37.04, -0.82],
    [-10.14, -0.82, 7.23],
])
Sigma_v = Sp_inv @ Sigma_theta @ Sp_inv.T
print("\nSigma_v =")
print(np.round(Sigma_v, 2))
# 출력: Table 5.1 의 [10.44, -0.92, -0.11; -0.92, 6.64, -0.94; -0.11, -0.94, 0.19] 와 거의 일치
# (Hedeker 책의 반올림으로 미세 차이 가능)

7.2 Step 2: statsmodels — 두 척도로 적합하여 동치 확인

import pandas as pd
import statsmodels.formula.api as smf

# § 5.2 의 합성 데이터 재사용 (위 코드에서 만든 df)
np.random.seed(2026)
n_subj, n_time = 50, 6
t = np.arange(n_time)

beta_true = np.array([24.0, -2.6, 0.05])
Sigma_v_true = np.array([
    [10.4, -0.9, -0.1],
    [-0.9, 6.6, -0.9],
    [-0.1, -0.9, 0.2],
])
sigma2_true = 10.5

v = np.random.multivariate_normal(np.zeros(3), Sigma_v_true, size=n_subj)
y = np.zeros((n_subj, n_time))
mu = beta_true[0] + beta_true[1] * t + beta_true[2] * t**2
for i in range(n_subj):
    indiv = (beta_true[0] + v[i, 0]) + (beta_true[1] + v[i, 1]) * t \
            + (beta_true[2] + v[i, 2]) * t**2
    y[i] = indiv + np.random.normal(0, np.sqrt(sigma2_true), n_time)

df = pd.DataFrame({
    "subject": np.repeat(np.arange(n_subj), n_time),
    "week": np.tile(t, n_subj),
    "y": y.flatten(),
})
df["week2"] = df["week"] ** 2

# 직교 다항식 변수 생성
T_full = np.column_stack([np.ones(n_time), t, t**2])
S_full = np.linalg.cholesky(T_full.T @ T_full)
Sp_inv_full = np.linalg.inv(S_full.T)
T_ortho_full = T_full @ Sp_inv_full  # n_time x 3

# long format 에 직교 다항식 변수 추가
df["cons"] = T_ortho_full[df["week"].values, 0]  # 모두 1/sqrt(6)
df["lin"] = T_ortho_full[df["week"].values, 1]
df["quad"] = T_ortho_full[df["week"].values, 2]

# 모형 1: 원본 척도
m_raw = smf.mixedlm(
    "y ~ week + week2",
    df,
    groups=df["subject"],
    re_formula="~week + week2",
).fit(reml=False)

# 모형 2: 직교 척도 (cons 가 절편 역할이라 별도 -1 처리)
m_ortho = smf.mixedlm(
    "y ~ -1 + cons + lin + quad",
    df,
    groups=df["subject"],
    re_formula="~-1 + cons + lin + quad",
).fit(reml=False)

print(f"raw  -2*logL = {-2 * m_raw.llf:.4f}")
print(f"ortho -2*logL = {-2 * m_ortho.llf:.4f}")
# 두 값이 동일 (수치 오차 한도 내)

# 변환 검증: gamma_hat 을 beta_hat 으로
gamma_hat = m_ortho.fe_params.values
beta_from_ortho = Sp_inv_full @ gamma_hat
print(f"\nbeta from ortho transform = {np.round(beta_from_ortho, 4)}")
print(f"beta from raw fit         = {np.round(m_raw.fe_params.values, 4)}")
# 두 값이 일치

R 사용자 — poly() 한 줄로

R 의 poly(week, degree = 2) 는 기본으로 직교 다항식을 반환한다 (Cholesky 인수분해 자동 적용).

library(lme4)
m_ortho <- lmer(y ~ poly(week, 2) + (poly(week, 2) | subject),
                data = df, REML = FALSE)
m_raw   <- lmer(y ~ week + I(week^2) + (week + I(week^2) | subject),
                data = df, REML = FALSE)

# 두 모형의 deviance 비교
deviance(m_ortho)  # 동일
deviance(m_raw)    # 동일

7.3 Step 3: 수렴 안정성 비교

직교 다항식이 수치적으로 더 안정적임을 시연한다.

import time

# 같은 데이터에서 두 모형 적합 시간 측정
times_raw, times_ortho = [], []
for seed in range(5):
    np.random.seed(seed)
    df_jitter = df.copy()
    df_jitter["y"] += np.random.normal(0, 0.5, len(df))

    t0 = time.time()
    smf.mixedlm("y ~ week + week2", df_jitter,
                groups=df_jitter["subject"],
                re_formula="~week + week2").fit(reml=False)
    times_raw.append(time.time() - t0)

    t0 = time.time()
    smf.mixedlm("y ~ -1 + cons + lin + quad", df_jitter,
                groups=df_jitter["subject"],
                re_formula="~-1 + cons + lin + quad").fit(reml=False)
    times_ortho.append(time.time() - t0)

print(f"raw  적합 시간 평균: {np.mean(times_raw):.3f}s")
print(f"ortho 적합 시간 평균: {np.mean(times_ortho):.3f}s")

직교 척도가 일반적으로 더 빠르고 수렴 경고가 적게 나온다. 시점 수가 많거나 차수가 높은 모형 (cubic, quartic) 일수록 차이가 두드러진다.

8 언제 직교 다항식을 써야 하는가

결정 가이드

원본 척도가 적합한 경우

시점 수가 적음 (예: 3~4 시점)
모형이 quadratic 까지만
이차 항의 부호와 절대 크기 자체에 임상적 의미가 있음 (예: 약물 효과의 가속도)
\(t^*\) (flattening point) 같은 직접 해석이 필요함

직교 다항식이 적합한 경우

시점 수가 많음 (예: 6 시점 이상)
cubic 이상의 고차 모형
수치적 수렴 문제 발생
분산 기여 비율을 직접 비교하고 싶음
baseline 보다 평균 수준 해석이 더 의미 있음 (예: 평균 우울 수준)

실무적으로는 둘 다 적합한 뒤 모수 변환으로 비교 하는 것이 가장 안전하다. 모형은 동일하므로 적합도는 같고, 해석에 더 적합한 좌표계를 골라 보고하면 된다.

9 핵심 정리

한 페이지 요약

다공선성 문제: \(t\) 와 \(t^2\) 의 상관이 1 에 가까워 추정 불안정 → 시간 중심화로 부분 해결, 직교 다항식으로 완전 해결
직교 다항식 = 다항식 항의 whitening: \(T(S^\top)^{-1}\) 로 직교 + 단위 분산 동시 달성
Cholesky 절차 5 단계: \(T \to T^\top T \to SS^\top \to (S^\top)^{-1} \to T(S^\top)^{-1}\)
두 모형은 수치 동일: \(-2 \log L\) 같음. 모수 변환은 \(S^\top\) 행렬로 결정
모수 변환: \(\beta = (S^\top)^{-1}\gamma\), \(\Sigma_\upsilon = (S^\top)^{-1}\Sigma_\theta S^{-1}\), \(V(\hat\beta) = (S^\top)^{-1}V(\hat\gamma)S^{-1}\)
분산 행렬 SE 변환: vec/vech 와 Kronecker 곱 (식 5.23) — 통계 패키지 자동 처리
절편 의미: 원본 = 첫 시점 평균, 직교 = 전체 시간 평균
랜덤 효과 상관 부호: 좌표계에 따라 뒤집힐 수 있음 — 같은 통계 사실의 두 측면
분산 기여 비율 비교: 직교 척도에서만 의미 있음 (회귀자가 같은 단위 분산)

10 다음 단계

주제	다룰 내용	위치
§ 5.3.4 고차 다항식	\(n - 1\) 차 한계, 분산-공분산 모수 한도	§ 5.3.5 sub-post (작성 예정)
§ 5.3.5 cubic 적합	Reisby Table 5.3 — 3 차 항 LR 검정, 분산 비율 (68.7/22.5/5.7/3.1), S 자 곡선	같은 sub-post
Ch.6 CPM	랜덤 효과 없이 시간 공분산 직접 모형화	미작성

11 관련 주제

선행 지식

Ch.5 Overview — 다항식 추세 MRM — 통합 관점
§ 5.2 — 곡선형 추세 MRM — 원본 척도 이차 모형
§ 4.4 — 행렬 공식화 — \(V(y) = Z\Sigma_\upsilon Z' + \sigma^2 I\) 유도
§ 4.5 — MRM 추정론 — ML/REML 의 일반 원리

관련 (다른 맥락의 직교 다항식)

§ 3.2 — MANOVA 일표본 (직교 다항식 P 행렬) — MANOVA 시간 효과 분해에서의 직교 다항식
§ 3.4 — Bock 수치 계산

후속 주제

§ 5.3.5 sub-post — Cubic Trend (작성 예정)
Ch.6 — Covariance Pattern Models (미작성)
Ch.7 Overview — MRM with Autocorrelated Errors

교재

Hedeker, D. & Gibbons, R. D. (2006). Longitudinal Data Analysis, Wiley, Ch.5 §5.3 (pp. 86-94)
Bock, R. D. (1975). Multivariate Statistical Methods in Behavioral Research, McGraw-Hill — Cholesky 인수분해 절차 원전
Pearson, E. S. & Hartley, H. O. (1976). Biometrika Tables for Statisticians — 등간격 직교 다항식 표
McCulloch, C. E. (1982). “Symmetric matrix derivatives with applications”, JASA — vec/vech 변환
Magnus, J. R. (1988). Linear Structures, Oxford University Press — \(G^+\) Moore-Penrose 의사역