1 들어가며

Ch.5 Overview 에서 다항식 추세 MRM 의 전체 그림을 짧게 살폈다. 이 sub-post 는 그 출발점인 이차 추세 (quadratic trend) 모형 만 깊이 다룬다.

모형 형태와 모수의 의미
시간 효과 미분 — 추세가 평탄해지는 시점 (\(t^*\)) 의 계산과 임상적 함의
곡선 4 가지 패턴 (가속·감속 × 양·음)
Reisby 정신과 임상 데이터 (Hedeker Table 5.1) 의 적합 결과 해석
모형이 적합한 분산-공분산 행렬 \(\widehat V(y) = Z \widehat\Sigma_\upsilon Z' + \widehat\sigma^2 I\) 의 수치 재현
“모집단 평균은 직선이지만 개인은 곡선” 패턴의 통계적·임상적 함의 (Figure 5.3 직관)

직교 다항식 표현으로 같은 모형을 재표현하는 절차는 § 5.3 sub-post (작성 예정) 에서 별도로 다룬다.

2 모형 정의

정의: 곡선형 추세 MRM (Curvilinear Trend Mixed-Effects Model)

Level-1 (피험자 \(i\) 의 시점 \(t_{ij}\) 에서):

\[ y_{ij} = b_{0i} + b_{1i} t_{ij} + b_{2i} t_{ij}^2 + \varepsilon_{ij} \tag{5.1} \]

Level-2 (개인 모수의 모집단 분포):

\[ b_{0i} = \beta_0 + \upsilon_{0i}, \quad b_{1i} = \beta_1 + \upsilon_{1i}, \quad b_{2i} = \beta_2 + \upsilon_{2i} \tag{5.2} \]

오차 \(\varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2)\) 는 시점 간 독립. 랜덤 효과 \(\upsilon_i = (\upsilon_{0i}, \upsilon_{1i}, \upsilon_{2i})^\top \sim \mathcal{N}_3(0, \Sigma_\upsilon)\), \(\Sigma_\upsilon\) 는 \(3 \times 3\) 분산-공분산 행렬.

이 식을 다음과 같이 묶어 쓰면 의미가 더 분명해진다.

\[ y_{ij} = b_{0i} + (b_{1i} + b_{2i} t_{ij}) \, t_{ij} + \varepsilon_{ij} \]

괄호 안 \(b_{1i} + b_{2i} t_{ij}\) 가 시간이 지남에 따라 변하는 순간 시간 효과 (instantaneous time effect) 다. 선형 모형에서 시간 효과는 상수 \(b_{1i}\) 이지만, 이차 항이 추가되면 시간 효과 자체가 시간의 1 차 함수가 된다.

2.1 직관: \(t^2\) 한 항이 만들어내는 변화

모형	시간 효과 \(\partial y / \partial t\)	운동학 비유
선형 (\(\beta_1 t\))	\(\beta_1\) — 상수	등속 운동
이차 (\(+ \beta_2 t^2\))	\(\beta_1 + 2\beta_2 t\) — 1 차 함수	등가속 운동
3 차 (\(+ \beta_3 t^3\))	\(\beta_1 + 2\beta_2 t + 3\beta_3 t^2\) — 2 차 함수	가속도가 변하는 운동

\(\beta_2\) 의 부호가 가속도의 부호다. \(\beta_2 > 0\) 이면 시간 효과가 점점 커지고 (\(+\) 가속), \(\beta_2 < 0\) 이면 점점 작아진다 (\(-\) 가속, 감속).

이 비유 덕분에 \(\beta_2\) 의 임상적 의미를 곧바로 해석할 수 있다.

우울증 점수 (\(y\)): \(\beta_1 < 0\) 이면 호전 추세, \(\beta_2 > 0\) 이면 호전 속도가 둔화
학습 점수: \(\beta_1 > 0\) 이면 향상 추세, \(\beta_2 < 0\) 이면 향상 속도가 떨어짐 (천장효과)

3 Flattening Point — 추세가 평탄해지는 시점

이차 모형에서 시간 효과가 0 이 되는 시점은 미분으로 직접 구한다.

\[ \frac{\partial y}{\partial t} = \beta_1 + 2\beta_2 t = 0 \quad \Longrightarrow \quad t^* = -\frac{\beta_1}{2\beta_2} \tag{5.3} \]

왜 \(t^*\) 가 임상적으로 중요한가

\(t^*\) 는 추세가 부호를 뒤집는 변곡점 이다.

\(t < t^*\): \(y\) 가 한 방향 (예: 호전)
\(t = t^*\): 시간 효과 = 0, 곡선이 평탄
\(t > t^*\): 반대 방향 (예: 악화)

이 한 줄의 수식이 의사결정에 직접 들어간다.

\(t^*\) 가 추적 기간 내부 → 처치 효과가 한 번 뒤집힘 → 처치 기간 재설계 필요
\(t^*\) 가 추적 기간 외부 → 곡선성은 있으나 부호 뒤집힘 없음 → 단조 가속/감속 패턴

특히 평정 척도 데이터에서 \(t^*\) 를 계산하지 않고 \(\beta_2\) 만 보고하면, “곡선이 있다” 는 사실은 알아도 언제 평탄화되는지 를 놓친다.

3.1 4 가지 곡선 패턴

이차 모형 하나로 표현 가능한 추세는 부호 조합으로 4 가지가 된다.

패턴	\(\beta_1\)	\(\beta_2\)	형태	임상 예
Decelerating positive	\(+\)	\(-\)	상승 후 평탄화	학습 곡선 (천장)
Accelerating positive	\(+\)	\(+\)	상승이 가속	지수형 성장 초기
Decelerating negative	\(-\)	\(+\)	하락 후 평탄화	약물 호전 (Reisby 형태)
Accelerating negative	\(-\)	\(-\)	하락이 가속	급격한 악화

여기에 \(|\beta_2|\) 가 충분히 크면 \(t^*\) 가 추적 기간 내부에 들어가 J 자 또는 U 자가 된다. 즉 \(\beta_0, \beta_1, \beta_2\) 세 모수만으로 단조 증가/감소, J 자, U 자, 역 J 자, 역 U 자를 모두 표현할 수 있다.

4 Reisby 정신과 데이터 적용 — Hedeker Table 5.1

Hedeker (2006) 가 우울증 임상시험 (Reisby et al., HDRS 점수 6 주 추적) 데이터에 식 (5.1)·(5.2) 를 적합한 결과는 다음과 같다.

Table 5.1 — 이차 추세 MRM 추정 결과 (원본 시간 척도)

모수	추정값	SE	\(Z\)	\(p <\)
\(\beta_0\) (절편)	23.76	0.55	43.04	.0001
\(\beta_1\) (선형)	-2.63	0.48	-5.50	.0001
\(\beta_2\) (이차)	0.05	0.09	0.58	.56
\(\sigma_{\upsilon_0}^2\)	10.44	3.59
\(\sigma_{\upsilon_0 \upsilon_1}\)	-0.92	2.41
\(\sigma_{\upsilon_1}^2\)	6.64	2.76
\(\sigma_{\upsilon_0 \upsilon_2}\)	-0.11	0.42
\(\sigma_{\upsilon_1 \upsilon_2}\)	-0.94	0.49
\(\sigma_{\upsilon_2}^2\)	0.19	0.09
\(\sigma^2\) (오차)	10.52	1.11

Note: \(-2 \log L = 2207.64\). (Hedeker, 2006, Ch.5, p.84)

4.1 Flattening Point 수치 계산

식 (5.3) 에 추정값을 대입한다.

\[ \hat{t}^* = -\frac{\hat\beta_1}{2\hat\beta_2} = -\frac{-2.63}{2 \times 0.05} = 26.3 \text{ 주} \]

연구 기간이 0~5 주이므로 \(\hat{t}^* = 26.3\) 은 추적 기간 훨씬 밖 이다. 모집단 평균에서는 실질적으로 단조 감속 음의 추세 (decelerating negative) 이고, 곡선성이 매우 미미하다. \(\beta_2\) 의 Wald 검정이 비유의 (\(p \approx .56\)) 인 것도 같은 결론이다.

4.2 “평균은 직선·개인은 곡선” 패턴

\(\beta_2\) 만 보면 곡선성이 없다고 결론 짓기 쉽지만, 랜덤 효과 분산 을 함께 보면 다른 그림이 나온다.

두 가지 LR 검정 비교

선형 모형 (이차 항 전부 제거) vs 이차 모형:

\[ H_0: \beta_2 = \sigma_{\upsilon_2}^2 = \sigma_{\upsilon_0 \upsilon_2} = \sigma_{\upsilon_1 \upsilon_2} = 0 \]

\(\Delta \text{deviance} = 11.4\), df = 4 → 유의 (p < 0.05, 분산 검정 보정 없이도)

이차 항 고정 효과만 (랜덤 효과 분산은 0): \(\beta_2 \neq 0\), \(\sigma_{\upsilon_2}^2 = \sigma_{\upsilon_0 \upsilon_2} = \sigma_{\upsilon_1 \upsilon_2} = 0\) vs 이차 모형 (랜덤 효과 포함):

\(\Delta \text{deviance} = 11.0\), df = 3

→ 모형 개선의 거의 전부 (11.4 중 11.0) 가 이차 항을 랜덤 효과로 포함 한 데서 온다. \(\beta_2\) 자체 (모집단 평균) 의 기여는 0.4 에 불과하다.

이 결과를 직관적으로 설명하면 Figure 5.3 의 그림이 된다.

   y │   Subject 2 (가속 양의 추세)
     │      ╱─
     │    ╱        Average (직선)
     │  ╱─────────────
     │╱        ╲─
     │            ╲─  Subject 1 (감속 양의 추세)
     │
     └─────────────── t

Subject 1 은 상승 후 둔화 (decelerating positive), Subject 2 는 가속 상승 (accelerating positive). 두 곡선이 부호가 다른 곡선성을 가져서 평균에서 cancel out 된다 — 결과적으로 평균 곡선은 거의 직선.

실무 함의 — “평균만 보면 놓치는 것”

모집단 평균 추세가 직선이라고 해서 모든 개인이 직선으로 변한다는 의미가 아니다. \(\sigma_{\upsilon_2}^2\) 가 유의하면:

개인별 반응 다양성이 임상적으로 실재함
평균 처치 효과 외에 누가 가속 호전, 누가 감속 호전인가 가 핵심 정보
의료 의사결정에서는 random effect 의 경험 베이즈 (EB) 추정으로 개인별 곡선을 그려야 함

5 분산-공분산 구조 \(\widehat V(y) = Z \widehat\Sigma_\upsilon Z' + \widehat\sigma^2 I\)

이차 추세 MRM 의 마진 분포는 다음과 같다.

\[ y_i \mid \beta, \Sigma_\upsilon, \sigma^2 \sim \mathcal{N}(X_i \beta, \, V(y_i)) \]

\[ V(y_i) = Z_i \Sigma_\upsilon Z_i^\top + \sigma^2 I_{n_i} \tag{5.4} \]

Reisby 처럼 6 시점 (\(t = 0, 1, \ldots, 5\)) 에서 모든 시점이 관측된 경우, 디자인 행렬은 다음과 같다.

\[ Z' = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 2 & 3 & 4 & 5 \\ 0 & 1 & 4 & 9 & 16 & 25 \end{bmatrix}, \qquad \widehat\Sigma_\upsilon = \begin{bmatrix} 10.44 & -0.92 & -0.11 \\ -0.92 & 6.64 & -0.94 \\ -0.11 & -0.94 & 0.19 \end{bmatrix} \]

오차 분산 \(\widehat\sigma^2 = 10.52\).

직관: 이 행렬이 의미하는 것

\(Z_i \Sigma_\upsilon Z_i^\top\) 부분은 랜덤 효과로 만들어지는 시간 간 공분산 이다. - 같은 사람의 두 시점 \(t_j, t_{j'}\) 가 모두 절편·선형·이차 랜덤 효과를 공유하기 때문에 공분산이 발생한다. - \(t_j\) 와 \(t_{j'}\) 가 멀어질수록 \(\upsilon_{1i} t_j \cdot \upsilon_{1i} t_{j'}\) 같은 항이 커지므로 후반 시점일수록 분산이 커지는 경향이 자연스럽다.

\(\sigma^2 I\) 부분은 시점-내 독립 측정 오차 다. 같은 사람의 두 시점에서 임의로 발생하는 흔들림.

두 성분을 더하면 관측된 반복측정의 전체 분산-공분산이 된다.

5.1 수치 재현 — 첫 행과 대각 성분

식 (5.4) 를 직접 계산하여 모형이 적합한 공분산 구조를 확인한다. 이는 § 4.4 — 행렬 공식화 에서 다룬 마진 분포 공식과 동일하다.

import numpy as np

# Reisby 6 시점, 이차까지
t = np.arange(6)
Z = np.column_stack([np.ones_like(t), t, t**2])  # 6 x 3

# Hedeker Table 5.1 추정값
Sigma_v = np.array([
    [10.44, -0.92, -0.11],
    [-0.92,  6.64, -0.94],
    [-0.11, -0.94,  0.19],
])
sigma2 = 10.52

# 마진 분산-공분산
V_hat = Z @ Sigma_v @ Z.T + sigma2 * np.eye(6)

print("적합된 V(y):")
print(np.round(V_hat, 2))

print("\n대각 (시점별 분산):")
print(np.round(np.diag(V_hat), 2))

# 시점 간 표준편차
print("\n시점별 표준편차:")
print(np.round(np.sqrt(np.diag(V_hat)), 2))

기대 출력 (Hedeker p.85 표시 부분과 일치):

\(\widehat{\text{Var}}(Y_0) \approx 20.96\)
\(\widehat{\text{Var}}(Y_1) \approx 23.86\)
\(\widehat{\text{Var}}(Y_2) \approx 31.07\)
\(\widehat{\text{Var}}(Y_3) \approx 38.31\) …

대각 성분이 시간이 지남에 따라 단조 증가한다. 이것이 평정 척도 데이터에서 흔한 시점 후반 분산 증가 패턴이다 — 호전 속도가 사람마다 달라서, 시간이 갈수록 개인 차가 벌어진다.

5.2 적합 분산 vs 관측 분산 비교

Hedeker (2006, Table 5.4) 는 관측 표준편차와 모형 기반 추정 표준편차를 시점별로 비교한다.

시점	관측 SD	랜덤 절편만	랜덤 선형	랜덤 이차	랜덤 3 차
Week 0	4.53	5.93	4.98	4.58	4.50
Week 1	4.70	5.93	4.91	4.88	4.73
Week 2	5.49	5.93	5.24	5.57	5.31
Week 3	6.41	5.93	5.92	6.19	6.37
Week 4	6.97	5.93	6.84	6.78	7.06
Week 5	7.22	5.93	7.91	7.69	7.32

표가 말해주는 것

랜덤 절편만: 모든 시점에서 SD = 5.93 (상수). 시점에 따른 분산 변화를 전혀 잡지 못함
랜덤 선형 추가: 분산이 시점에 따라 변하기 시작 — 후반 시점에서 관측에 가까워짐
랜덤 이차 추가: 곡선형으로 분산 증가를 더 잘 따라감
랜덤 3 차 추가: Week 0~2 의 미세한 비단조 변화까지 반영

랜덤 효과 차수가 높아질수록 공분산 구조의 적합도 가 좋아진다. 이것이 평균 추세는 직선처럼 보여도 랜덤 이차 항을 포함하는 강한 정당화 — 모형이 적합한 분산 구조가 데이터 분산 구조와 더 일치한다.

6 합성 데이터로 적합 절차 시연

Reisby 데이터 자체는 공개되지 않았으므로, 그와 비슷한 구조의 합성 데이터로 같은 흐름을 재현한다.

6.1 Step 1: 순수 Numpy 로 우도 직접 계산

랜덤 효과 분산을 적분 (사실상 marginal 분포 사용) 으로 처리하면 가우스 다변량의 로그 우도가 된다.

import numpy as np

# 50명, 6 시점 합성 데이터 (Reisby 비슷한 모수)
np.random.seed(2026)
n_subj, n_time = 50, 6
t = np.arange(n_time)
Z = np.column_stack([np.ones_like(t), t, t**2])  # 6 x 3

beta = 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)
mu = Z @ beta  # 평균 곡선 (모든 사람 동일)
y = np.zeros((n_subj, n_time))
for i in range(n_subj):
    y[i] = mu + Z @ v[i] + np.random.normal(0, np.sqrt(sigma2_true), n_time)

# 마진 로그 우도 — 단일 사람당 N(X beta, Z Sigma_v Z' + sigma^2 I)
def marginal_loglik(params, y, Z):
    beta = params[:3]
    sigma2 = np.exp(params[3])  # 양수 보장
    # Sigma_v: 3x3 대칭, Cholesky 모수화
    L = np.zeros((3, 3))
    L[np.tril_indices(3)] = params[4:10]
    L[np.diag_indices(3)] = np.exp(np.diag(L))
    Sigma_v = L @ L.T

    n_subj, n_time = y.shape
    V = Z @ Sigma_v @ Z.T + sigma2 * np.eye(n_time)
    sign, logdetV = np.linalg.slogdet(V)
    Vinv = np.linalg.inv(V)
    mu = Z @ beta
    ll = 0
    for i in range(n_subj):
        r = y[i] - mu
        ll += -0.5 * (n_time * np.log(2 * np.pi) + logdetV + r @ Vinv @ r)
    return ll

# 한 점에서의 로그 우도 평가 (실제 모수에서)
init = np.concatenate([
    beta,
    [np.log(sigma2_true)],
    np.array([np.log(np.sqrt(10.4)), -0.9/np.sqrt(10.4), 0,
              np.log(np.sqrt(6.6)), 0, np.log(np.sqrt(0.2))])
])
ll_at_truth = marginal_loglik(init, y, Z)
print(f"진실 모수에서 log-likelihood = {ll_at_truth:.2f}")

이 코드의 핵심: \(V_i = Z \Sigma_\upsilon Z^\top + \sigma^2 I\) 를 직접 만들고 다변량 정규의 로그 우도를 계산한다. 실제 ML/REML 추정은 이를 모수에 대해 최적화하는 것 — scipy.optimize.minimize 로 확장 가능하다.

6.2 Step 2: statsmodels 로 quadratic MRM 적합

실무에서는 프레임워크를 사용한다.

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

# 합성 데이터를 long format 으로 변환
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

# 이차 추세 MRM (랜덤 절편 + 랜덤 선형 + 랜덤 이차)
m_quad = smf.mixedlm(
    "y ~ week + week2",
    df,
    groups=df["subject"],
    re_formula="~week + week2",
).fit(reml=False)

print(m_quad.summary())

# Flattening point
beta1_hat = m_quad.fe_params["week"]
beta2_hat = m_quad.fe_params["week2"]
t_star = -beta1_hat / (2 * beta2_hat) if abs(beta2_hat) > 1e-6 else float("inf")
print(f"\nFlattening point t* = {t_star:.2f} weeks")
print(f"(연구 기간 0~5주 {'내부' if 0 <= t_star <= 5 else '외부'})")

해석 포인트:

summary() 의 week2 행이 비유의여도 그것만으로 “곡선성 없음” 으로 결론 짓지 않는다
랜덤 효과 분산 (Var(week2)) 를 별도 LR 검정으로 평가
\(t^*\) 가 연구 기간 내부에 있으면 추세 부호 뒤집힘 발생, 외부면 단조 가속/감속

R lme4 에서

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

I(week^2) 로 원본 척도를 유지한다. poly(week, 2) 를 쓰면 자동으로 직교 다항식이 적용되어 § 5.3 표현이 된다.

6.3 Step 3: 개인별 곡선과 평균 곡선 시각화

Figure 5.3 의 직관 — “평균은 직선·개인은 곡선” — 을 재현한다.

import matplotlib.pyplot as plt

# 경험 베이즈 (BLUP) 랜덤 효과 추정
re_blup = m_quad.random_effects  # dict {subject_id: random effects}

# 각 사람의 적합 곡선
t_grid = np.linspace(0, 5, 100)
fig, ax = plt.subplots(figsize=(7, 5))

for sid, re in re_blup.items():
    b0 = m_quad.fe_params["Intercept"] + re["Group"]
    b1 = beta1_hat + re.get("week", 0)
    b2 = beta2_hat + re.get("week2", 0)
    ax.plot(t_grid, b0 + b1 * t_grid + b2 * t_grid**2,
            color="gray", alpha=0.3, linewidth=0.7)

# 평균 곡선 (고정 효과만)
b0_mean = m_quad.fe_params["Intercept"]
ax.plot(t_grid, b0_mean + beta1_hat * t_grid + beta2_hat * t_grid**2,
        "b-", linewidth=2.5, label="Average (fixed effects)")

ax.set_xlabel("Week")
ax.set_ylabel("HDRS-like score")
ax.set_title("Individual quadratic trends vs average")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()

생성된 그래프에서 회색 곡선들은 다양한 곡선성을 보이지만 (가속 호전, 감속 호전, 일부는 V 자), 평균 곡선 (파란선) 은 거의 직선에 가깝게 나타난다. 이것이 Hedeker Figure 5.3 이 강조하는 패턴 — 개인 곡선성이 부호로 cancel out 되어 평균에서 사라짐.

7 코드 검증 포인트

적합이 잘 되었는지 확인하는 체크리스트

\(\widehat{\text{Var}}(y_{ij})\) 의 시점별 값이 관측 분산과 가까운가? (Hedeker Table 5.4 형태로 비교)
\(\widehat\Sigma_\upsilon\) 가 양정치 (positive definite) 인가? np.linalg.eigvalsh(Sigma_v_hat).min() > 0
수렴 경고가 없는가? statsmodels 의 m.converged, lme4 의 singular fit 메시지 확인
시간을 중심화 (\(t - \bar{t}\)) 하면 다공선성이 줄어 수치 안정 — 안 되면 직교 다항식 사용 (§ 5.3)

8 응용 — 곡선형 추세 모형이 적합한 상황

분야	상황	예상 곡선 패턴	\(\beta_2\) 부호
정신과 임상시험	항우울제 치료 추적	Decelerating negative	\(+\)
약동학	단회 투여 후 혈중 농도	상승 후 하강 (역 U)	\(-\) (강)
발달 심리	영아 인지 점수 성장	Accelerating positive	\(+\)
학습 심리	새 기술 학습 (천장 효과)	Decelerating positive	\(-\)
통증 척도	만성 통증의 시간 추이	다양 — 개인별 분산 큼	변동

곡선형 추세 모형 (이차) 의 핵심 가치는 세 모수 (\(\beta_0, \beta_1, \beta_2\)) 만으로 6 가지 이상의 패턴 을 모두 표현할 수 있다는 점이다. 선형 모형으로 두 모수 (\(\beta_0, \beta_1\)) 만 가지면 단조 증가/감소만 가능하다.

9 다음 단계

주제	다룰 내용	위치
§ 5.3 직교 다항식	다공선성 해결, Cholesky 인수분해, 모수 변환, SE 변환 (\(G^+, vec/vech\))	§ 5.3 sub-post (작성 예정)
§ 5.3.5 cubic	3 차 추세 적합, 분산 기여 비율, S 자 곡선 해석	§ 5.3.5 sub-post (작성 예정)
Ch.6 CPM	랜덤 효과 없이 시간 공분산을 직접 모형화	미작성

10 관련 주제

선행 지식

Ch.5 Overview — 다항식 추세 MRM — 본 sub-post 의 통합 관점
§ 4.3 — 랜덤 절편·추세 MRM — 선형 추세 MRM 의 구조
§ 4.4 — 행렬 공식화 — \(V(y) = Z \Sigma_\upsilon Z' + \sigma^2 I\) 유도
§ 4.5 — MRM 추정론 — ML/REML 로 \(\Sigma_\upsilon\) 추정