01.Rmd

# 서론 {#intro}

비선형 혼합효과 모델링(Nonlinear mixed effect modeling)을 위해 가장 흔히 쓰이는 도구는 NONMEM®이라고 하는 소프트웨어이다. 현재 모델기반 약물개발 (Model-based drug development)를 위한 계량약리학 분야에서 가장 표준적인 도구로 생각되고 있으나 그 동작하는 원리를 이해하기는 쉽지 않다. 공개소프트웨어인 R을 통해서 NONMEM의 계산과 알고리듬을 구현하는 시도가 있었고, 성공적으로 결과값을 재현할 수 있음을 보였다. [@kim2015r; @bae2016r] 이 연구를 통해서 NONMEM의 동작원리를 단계별로 실행하며 학습할 수 있게 되었다. 더 나아가 저자들은 공개 소프트웨어인 nmw R 패키지를 개발 [@R-nmw], 배포하여 이러한 원리를 누구나 비용없이 접근할 수 있게 하였다. 

NONMEM의 실행과정은 대체로 initialization step, estimation step(추정단계), covariance step(공분산단계), table step(표단계)으로 나눌 수 있다. 공분산단계의 목적은 추정단계에서 추정된 파라미터(θ, Ω, Σ)들의 점 추정치(point estimate)에 대한 standard error(표준오차)를 구하는 것이다. 따라서, \$COV 절을 제어구문 파일에 포함시키지 않으면 표준오차 결과를 볼 수 없다. 공분산단계는 원칙적으로 추정단계가 성공했을 때만 의미가 있으므로 성공한 후에 실행되는데, 만약 추정 실패 시에도 종료시점에서의 값을 알고 싶으면 \$COV에 UNCOND 옵션을 추가해주면 된다. 공분산단계가 추정단계와 별도로 분리된 이유는 NONMEM이 사용자가 준 original scale에서 목적함수를 최소화하는 것이 아니라, unconstrained parameter (UCP)를 사용하기 때문이다. UCP는 과거에 scaled and transformed parameter (STP)라고도 불리었다.\index{공분산 / covariance}\index{목적함수 / objective function}\index{covariance / 공분산}\index{objective function / 목적함수}\index{standard errors of the estimates / 추정값의 표준오차}

## Maximum Likelihood Estimation

우도함수를 최대화하면서 모수를 추정하는 방법입니다. 
관측된 표본에 기초하여 관측 불가능한 파라미터(모수)를 추정하는 방법론 중 하나로, 표본들로부터 알려지지 않은 모집단 확률분포의 형태를 추정해가는 방법론입니다. (An estimation method that finds a combination of parameters which maximizes the likelihood of an observation dataset)

MLE는 주어진 Dataset(관측결과, D)의 발생 가능성(Likelihood, 우도)이 가장 높은 모수를 추정하고, 이를 바탕으로 최적의 f를 도출합니다. 이를 수식으로 표현하면 다음과 같습니다. 

정규분포의 경우 n~i~개의 관찰값을 갖는 대상 i에 대해서 -2 log likelihood (-2LL)은 다음과 같이 표현할 수 있습니다.

$$
-2LL_i = n_i \cdot log(2\pi) + \sum\limits_{j=1}^{n_i} log(\sigma_{ij}^2) + \sum\limits_{j=1}^{n_i} \frac{(y_{ij} - \mu_{ij})^2}{\sigma_{ij}^2}
$$

위 함수를 최소화 (우도를 최대화) 하는 모수의 조합은 무엇일까요? 그 결과가 maximum likelihood estimator/estimate (MLE)이고 그 목적함수는 다음과 같습니다.

$$
O_i = \sum\limits_{j=1}^{n_i} log(\sigma_{ij}^2) + \sum\limits_{j=1}^{n_i} \frac{(y_{ij} - \mu_{ij})^2}{\sigma_{ij}^2}
$$

MLE의 장점은 매우 직관적이라는 점입니다. 그리고 별다른 가정을 하지 않고 주어진 관측결과를 바탕으로 쉽계 계산할 수 있습니다. (물론, 관측결과가 이항분포를 따를 것이라는 가정이 들어갔지만, 다른 추정법에 비해서는 가정이 약한 편입니다.)

MLE는 빈도론을 기반으로 합니다. 관측결과에 의존해서 판단하기 때문에 관측결과의 크기(관측횟수)와 질(관측결과의 정확성)에 따라서 모수 추정의 성패가 크게 좌우됩니다.

## 이론적 배경

\index{MLE, maximum likelihood estimation}\index{likelihood / 가능도}

NONMEM의 수행은 대체로 초기화단계, 추정단계, 공분산단계, 표단계로 나눌 수 있다. 이중 공분산단계에 대하여 설명하고자 한다. 추정단계가 끝나면 최종 파라미터 추정값이 구해진다. NONMEM에서는 Y의 분포가 정규이든 아니든 상관없이 정규분포의 가능도 함수(엄밀히는 이것을 약간 변형한 것)를 목적함수로 한다. 즉, maximum likelihood estimation (MLE) 방법을 사용한다.\index{가능도 / likelihood}\index{공분산 / covariance}\index{목적함수 / objective function}\index{covariance / 공분산}\index{likelihood / 가능도}\index{objective function / 목적함수}

MLE에서는 loglikelihood가 목적함수인 경우에 Fisher’s Information Matrix의 역행렬을 이용하여 추정값들의 분산-공분산 행렬과 표준오차를 구할 수 있다. 표준오차는 분산-공분산 행렬의 대각원소를 1/2승 한 것이다.\index{공분산 / covariance}\index{목적함수 / objective function}\index{분산-공분산 행렬 / variance-covariance matrix}\index{covariance / 공분산}\index{likelihood / 가능도}\index{objective function / 목적함수}\index{variance-covariance matrix / 분산-공분산 행렬}

MLE 이론에 따르면 log likelihood (l)를 파라미터 추정값에서의 1계 편미분한 것의 제곱과 2계 편미분한 것이 부호만 반대인 같은 기대값을 가진다.\index{likelihood / 가능도}
\begin{equation}
\begin{split}
  \int & f(x;\theta)dx = 1 \\
  \int & exp(lnf(x;\theta))dx = 1
\end{split}
(\#eq:mle-theory)
\end{equation}

양변을 $\theta$에 대해 미분하면
\begin{equation}
\begin{split}
\int\left\{ \frac{\partial}{\partial\theta}lnf(x;\theta) \right\} exp(lnf(x;\theta))dx = 0 \\
E_{\theta}\left\lbrack \frac{\partial}{\partial\theta}lnf(X;\theta) \right\rbrack = 0
\end{split}
(\#eq:theta-differential)
\end{equation}

한 번 더 미분하면

\begin{equation}
\begin{split}
\int{\frac{\partial^{2}}{\partial\theta^{2}}lnf(x;\theta) + \lbrack \frac{\partial}{\partial\theta}lnf(x;\theta) \rbrack^{2} } exp(lnf(x;\theta))dx = 0 \\
E_{\theta}\left\{ \frac{\partial^{2}}{\partial\theta^{2}}lnf(X;\theta) + \left\lbrack \frac{\partial}{\partial\theta}lnf(X;\theta) \right\rbrack^{2} \right\} = 0 \\
I(\theta) = E_{\theta}\left\{ \frac{\partial^{2}}{\partial\theta^{2}}lnf(X;\theta) \right\} = E_{\theta}\left\{ - \left\lbrack \frac{\partial}{\partial\theta}lnf(X;\theta) \right\rbrack^{2} \right\} = Var_{\theta}\left\{ \frac{\partial}{\partial\theta}lnf(X;\theta) \right\} 
\end{split}
(\#eq:diff-diff)
\end{equation}

좀 더 복잡한 과정을 거치면 \(\sqrt{n}({\widehat{\theta}}^{\text{MLE}} - \theta)\)는 점근적으로 \(N(0,\lbrack I(\theta)\rbrack^{- 1})\) 를 따름을 증명할 수 있다. [@kim]

아주 이상적이지만(그러나 극히 드문) 상황에서는 목적함수의 1계 미분의 제곱과 2계 미분이 같겠지만, 실제 관찰값으로 계산하면\index{목적함수 / objective function}\index{objective function / 목적함수}
다르다. NONMEM에서는 1계 미분의 제곱행렬을 S matrix, 2계 미분 행렬(Hessian)을 R matrix라고\index{R matrix / R 행렬}\index{S matrix / S 행렬}
부른다. 대부분의 MLE software에서는 R matrix의 inverse만으로 estimate의 분산-공분산 행렬로 삼고\index{공분산 / covariance}\index{분산-공분산 행렬 / variance-covariance matrix}\index{covariance / 공분산}\index{variance-covariance matrix / 분산-공분산 행렬}\index{R matrix / R 행렬}
표준오차를 산출한다. NONMEM에서는 \(R^{- 1}SR^{- 1}\)을 estimate의 분산-공분산 행렬로 삼는다. 이것은\index{공분산 / covariance}\index{분산-공분산 행렬 / variance-covariance matrix}\index{covariance / 공분산}\index{variance-covariance matrix / 분산-공분산 행렬}
더 보수적인 방법이어서 신뢰구간이 더 넓게 나온다.

수식으로 표현하자면 다음과 같다. 여기에서 \(O_{i}\)는 개인(i)의 목적함수 값이며, \(O\)는 전체 목적함수 값이다.\index{목적함수 / objective function}\index{objective function / 목적함수}

\begin{equation}
\begin{split}
O = \sum_{i}^{}O_{i} \\
S = \sum_{i}^{}{\nabla_{O_{i}}{\nabla_{O_{i}}}^{T}} = \sum_{i}^{}\frac{\partial O_{i}}{\partial\theta}\frac{\partial O_{i}}{\partial\theta}^{T} \\
R = \frac{\partial^{2}O}{\partial\theta^{2}}
\end{split}
(\#eq:o-s-r)
\end{equation}

## First-order conditional estimation (FOCE) method 

FOCE는 FO 가정법보다 더 복잡합니다. EBE를 각 iteration 마다 추정하기 때문입니다.

> The FOCE is more complex than the first-order (FO) approximation method because it estimates the empirical Bayes estimate (EBE) for each iteration. By contrast, it is a further approximation of the Laplacian (LAPL) method, which uses second-order expansion terms. FOCE without INTERACTION can only be used for an additive error model, while FOCE with INTERACTION (FOCEI) can be used for any error model. The formula for FOCE without INTERACTION can be derived directly from the extension of the FO method, while the FOCE with INTERACTION method is a slight simplification of the LAPL method. Detailed formulas and R scripts are presented here for the reproduction of objective function values by NONMEM.