Chapter 16 비선형 혼합 효과 모형 (Nonlinear Mixed Effects Model)

비선형 혼합 효과 모형을 처리할 수 있는 가장 효율적인 프로그램은 현재 NONMEM 이라는 software이다.

NONMEM을 설치하면 12명에게 theophylline 320mg을 경구 투약한 후 (한 사람당) 11회에 걸처 농도를 측정한 자료가 있다.

R에도 동일한 자료가 Theoph라는 이름으로 있으며 그림을 그리면 다음과 같다. 

DATA = Theoph
colnames(DATA) = c("ID", "BWT", "DOSE", "TIME", "DV")
DATA[,"ID"] = as.numeric(as.character(DATA[,"ID"]))

require(lattice)
xyplot(DV ~ TIME | as.factor(ID), data=DATA, type="b")
Lattice plot of theophylline data

Figure 16.1: Lattice plot of theophylline data

위의 자료는 약동학(pharmacokinetics)에서 경구 투약 one-compartment model로 적합시킬 수 있다. 경구 one-compartment model의 수식은 다음과 같다.

\[C(t) = \frac{\text{DOSE}}{V}\frac{K_a}{K_a - K}(e^{-K t} - e^{-K_a t})\]

위는 구조 모형(structural model)이고, 여기에 확률변수를 추가하여 표현하면 다음과 같다. 이제 \(C(t)\)를 \(F\)로 나타내기로 한다.

\[Y = F + F \cdot \epsilon_1 + \epsilon_2\]

\[\begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \end{bmatrix} = \overrightarrow{\epsilon} \sim MVN \left( \overrightarrow{0}, \begin{bmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \\ \end{bmatrix} \right)\]

위의 \(Y\)는 observation, \(F\)는 prediction, \(\epsilon_1\)는 proportional error, \(\epsilon_2\)는 additive error를 의미한다.

그런데, 위의 \(K_a\), \(K\), \(V\)는 약동학적 모수(parameter)이지만, 개인별로 비교적 고유한 값을 가지며, (더 정밀하게 보면 투약 때마다 매번 약간씩 달라질 수 있는데, 이것을 inter-occasional variation이라 부르며, 적어도 단기간에는 크게 변하지 않으므로, 여기서는 고려하지 않는다.) 그것이 전체 인구집단에서는 비교적 정규분포 또는 로그정규분포를 따른다고 가정한다. 이 예제에서는 로그정규분포를 가정한 것이다. (만약 정규분포 또는 로그정규분포를 따르지 않는다 하더라도 다른 분포를 적용하기가 쉽지 않다.)

이와 같이, 전체 인구에서 변하지 않는 부분을 \(\theta\)라는 상수로 표현하고, 변하는 부분을 \(\eta\)라는 확률변수로 표기하기로 한다.

\(\eta\) 들은 한 개인 안에서는 비교적 고유의 값을 가지지만, 전체 인구에서는 비교적 정규분포를 따른다고 가정하며, \(\eta\) 간에는 상관관계도 있다.

이를 수식으로 표현하면 다음과 같다.

\(K_a = \theta_1 e^{\eta_1}\)

\(V = \theta_2 e^{\eta_2}\)

\(K = \theta_3 e^{\eta_3}\)

\(\begin{bmatrix} \eta_1 \\ \eta_2 \\ \eta_3 \\ \end{bmatrix} = \overrightarrow{\eta} \sim MVN(\overrightarrow{0}, \Omega)\)

위의 자료를 분석하는 NONMEM control file은 다음과 같다. 

$PROB THEOPHYLLINE ORAL
$INPUT ID AMT=DROP TIME DV BWT
$DATA ../THEO
$PRED
  DOSE = 320
  KA   = THETA(1) * EXP(ETA(1))
  V    = THETA(2) * EXP(ETA(2))
  K    = THETA(3) * EXP(ETA(3))
  F    = DOSE/V*KA/(KA - K)*(EXP(-K*TIME) - EXP(-KA*TIME))
  IPRE = F
  Y    = F + F*EPS(1) + EPS(2)

$THETA (0, 2) (0, 50) (0, 0.1)
$OMEGA BLOCK(3)
  0.2
  0.1  0.2
  0.1  0.1  0.2
$SIGMA 0.1  0.1
$EST MAX=9999 METHOD=ZERO POSTHOC
$COV PRINT=ERS
$TAB ID TIME IPRE CWRES
     FILE=sdtab NOPRINT ONEHEADER
$TAB ID ETA(1) ETA(2) ETA(3)
     FILE=IndiEta.txt NOAPPEND ONEHEADER FIRSTONLY

NONMEM의 전통적인 추정 방법은 first order approximation (FO method), first order approximation with conditional estimation (FOCE method), Laplacian method 등 3가지가 있다.

FO method인 경우에는 $EST 부분에 METHOD=ZERO라고 쓰고, FOCE method인 경우에는 METHOD=COND라고 쓰며, Laplacian method인 경우에는 METHOD=COND LAPL이라고 쓰면 된다 (이하 LAPL로 표기).

Additive error가 아닌 경우에는 모두 INTERACTION이라는 option을 추가해야 한다.

NONMEM의 FO method에 의한 실행 결과를 요약하자면 다음과 같다.

점추정치(point estimate)는 다음과 같다. 아래에서 hat 기호는 일괄적으로 표시하지 않았다. 즉 아래에서 \(\theta_1\)은 \(\hat{\theta_1}\)를 의미한다.

\[\begin{bmatrix} \theta_1 \\ \theta_2 \\ \theta_3 \\ \end{bmatrix} = \begin{bmatrix} 3.17 \\ 38.3 \\ 0.105 \\ \end{bmatrix}\]

\[\begin{bmatrix} \omega_{11} & \omega_{12} & \omega_{13} \\ \omega_{21} & \omega_{22} & \omega_{23} \\ \omega_{31} & \omega_{32} & \omega_{33} \\ \end{bmatrix} = \begin{bmatrix} 1.20 & \omega_{12} & \omega_{13} \\ 0.137 & 0.0313 & \omega_{33} \\ 0.370 & 0.0434 & 0.251 \\ \end{bmatrix}\]

\[\begin{bmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \\ \end{bmatrix} = \begin{bmatrix} 0.0121 & 0 \\ 0 & 0.0543 \\ \end{bmatrix}\]

위의 점추정치에 대응하는 위치에 standard error(SE)를 표시하면 다음과 같다. 아래에서 hat 기호와 SE는 모두 생략했다. 즉, \(\theta_1\)은 \(\hat{SE(\hat{\theta_1})}\)를 의미한다.

\[\begin{bmatrix} \theta_1 \\ \theta_2 \\ \theta_3 \\ \end{bmatrix} = \begin{bmatrix} 0.641 \\ 1.68 \\ 0.0231\\ \end{bmatrix}\]

\[\begin{bmatrix} \omega_{11} & \omega_{12} & \omega_{13} \\ \omega_{21} & \omega_{22} & \omega_{23} \\ \omega_{31} & \omega_{32} & \omega_{33} \\ \end{bmatrix} = \begin{bmatrix} 0.421 & \omega_{12} & \omega_{13} \\ 0.0822 & 0.0198 & \omega_{33} \\ 0.340 & 0.0230 & 0.289 \\ \end{bmatrix}\]

\[\begin{bmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \\ \end{bmatrix} = \begin{bmatrix} 0.00358 & 0 \\ 0 & 0.0321 \\ \end{bmatrix}\]

\(\theta\)들의 SE는 상당히 작은 반면에(<50%), \(\omega\), \(\sigma\)들의 SE는 상대적으로 크다. 하지만, 분산 term들을 Wald type test에 의해 유의성을 판정하거나, 생략하지 않는다.

NONMEM 추정의 주요 step은 초기화(initialization), 추정(estimation), 공분산(covariance), 표(table) step으로 나눌 수 있다.