14.Rmd

# NONMEM의 \$COVARIANCE {#cov}
\index{\$COVARIANCE (\$COV, \$COVAR)}

\Large\hfill
배균섭
\normalsize

---

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 / 추정값의 표준오차}\index{standard errors of the estimates / 추정값의 표준오차}\index{\$COVARIANCE (\$COV, \$COVAR)}

## 실제 사례 {#cov-example}

다음은 공분산단계의 결과물을 이해하기 위한 예제 제어구문 파일의 내용이다. (코드 \@ref(exm:theoph-control-stream)) Theophylline 데이터셋은 12명의 사람에게 320mg의 theophylline을 1회 경구 투여한 자료이며 NONMEM 설치시에 `THEO` 라는 파일 이름으로 포함되어 있다. 그림은 \@ref(Theoph)절을 참고한다.\index{공분산 / covariance}\index{covariance / 공분산}
 
```{example, theoph-control-stream, echo=TRUE}
Theophylline의 경구 투여 약동학 자료 모델링을 위한 제어구문 파일
```
\vspace{-5ex} 
```perl
$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
```

위의 \$COV STEP에 의해 생성된 부분은 [별첨 A](#nonmem-output-file)과 같으며, 원래 output은 126열까지까지 있으나, 지면의 제약으로 인해 104열이후는 생략하였다.\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV에 의해 생성되는 첫 번째 부분(section)은 추정치의 표준오차 부분으로 출력의 직전 section인 FINAL PARAMETER ESTIMATE과 배치 모양이 똑같다. 이후의 출력들은 대개 이를 구하기 위해 필요한 중간계산 결과들이다.\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV 출력의 두 번째 부분은 추정값(estimate)의 분산-공분산행렬(variance-covariance matrix)이다. 이 행렬의 대각원소에다 제곱근을 취한 것이 추정값의 표준오차이다. 이는 뒤에 나오는 R matrix와 S matrix를 사용하여 \(R^{- 1}SR^{- 1}\)를 계산한 것이다.\index{공분산 / covariance}\index{covariance / 공분산}\index{variance-covariance matrix / 분산-공분산 행렬}\index{추정값의 표준오차 / standard errors of the estimates}\index{standard errors of the estimates / 추정값의 표준오차}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{standard errors of the estimates / 추정값의 표준오차}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV 출력의 세 번째 부분은 직전 부분인 분산-공분산행렬을 상관행렬(correlation matrix)로 변환한 것이다.  다만, 대각원소는 제곱근을 취하여 standard deviation(여기서는 standard error가 됨)을 표시하고 있다.\index{공분산 / covariance}\index{상관행렬 / correlation matrix}\index{correlation matrix / 상관행렬}\index{covariance / 공분산}\index{standard errors of the estimates / 추정값의 표준오차}\index{standard errors of the estimates / 추정값의 표준오차}\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV 출력의 네 번째 부분은 분산-공분산행렬의 역행렬이다. 내부적으로는 이 행렬을 먼저 구한 뒤, 이것의 역행렬을 분산-공분산행렬로 삼는다. 이것의 의미는 Fisher’s Information Matrix (NONMEM에서의 R matrix)와 유사하지만, 계산법은 더 보수적이어서 다르다.\index{공분산 / covariance}\index{covariance / 공분산}\index{R matrix / R 행렬}\index{R matrix / R 행렬}\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV output의 다섯 번째 부분은 상관행렬로부터 구한 고유치(eigenvalue)들을 작은 것부터 순서대로 나열한 것이다.  만약 상관행렬이 positive definite(양정치)라면 이것은 모두 양의 실수로 나와야 한다. 여기에 대한 추가 설명은 \@ref(eigenvalue)를 참고한다.\index{상관행렬 / correlation matrix}\index{correlation matrix / 상관행렬}\index{\$COVARIANCE (\$COV, \$COVAR)}

$COV 출력의 여섯 번째 부분은 R matrix이다. 이는 log likelihood를 parameter들로 이계미분한 것(hessian matrix)으로 Fisher’s Information Matrix(FIM)이라고도 한다. 일반적인 최대가능도추정법(maximum likelihood estimation, MLE)에서는 이것의 역행렬을 분산-공분산 행렬로 삼지만 NONMEM에서는 다음에 나오는 S matrix와 함께 sandwich estimator를 사용한다.\index{가능도 / likelihood}\index{공분산 / covariance}\index{분산-공분산 행렬 / variance-covariance matrix}\index{covariance / 공분산}\index{likelihood / 가능도}\index{variance-covariance matrix / 분산-공분산 행렬}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{\$COVARIANCE (\$COV, \$COVAR)}

$COV 출력의 마지막 부분은 S matrix이다. 이는 개별 대상자의 log likelihood를 파라미터들로 일계미분값(gradient vector)에 대한 cross product(n by 1 열벡터를 그것의 transpose(1 by n vector)와 곱한 것으로 n by n matrix가 된다)들을 모두 합한 것이다. MLE 이론에 따르면 이상적인 상황에서 S matrix와 R matrix의 기대값은 같다. 하지만, 실제로는 상당히 다를 수 있고, 정규분포 가정도 어긋날 수 있으므로 NONMEM은 \(R^{- 1}SR^{- 1}\) 를 분산-공분산행렬로 삼는다.\index{공분산 / covariance}\index{covariance / 공분산}\index{gradient / 기울기}\index{likelihood / 가능도}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{R matrix / R 행렬}\index{S matrix / S 행렬}\index{\$COVARIANCE (\$COV, \$COVAR)}

만약 어떤 원인(예를 들어, R이 singular)에 의해 분산-공분산행렬을 구할 수 없고, 따라서, 표준오차를 구할 수 없다면 $COV에서 `MATRIX=R` 또는 `MATRIX=S` 옵션으로 단순화하여 표준오차를 구할 수 있다. 이 경우에는 대개 표준오차가 더 작게 나온다. 어떤 행렬이 singular(비정칙)이어서 역행렬을 구할 수 없다면 NONMEM은 pseudo-inverse matrix를 구하게 되는데, 이는 유일한(unique) 해는 아니므로 해석에 유의해야 한다.\index{공분산 / covariance}\index{covariance / 공분산}\index{\$COVARIANCE (\$COV, \$COVAR)}

\$COV step이 실패하는 경우에는 먼저 파라미터 추정값에 이상이 없는지 살펴보고, 이상이 없어 보이면 붓스트랩이나 profile-likelihood등의 다른 방법을 이용하여 신뢰구간을 구하면 된다. 이 때 신뢰구간이 비대칭일 가능성이 크므로 표준오차는 제시하기 어려워진다. 혹자는 \$COV에 의한 표준오차 보다 붓스트랩에 의한 신뢰구간이 (시간은 훨씬 많이 걸리지만) 더 좋다고 한다. 붓스트랩을 할 때는 \$COV 절을 삭제하는 것이 시간을 절약하는 방법이다.\index{붓스트랩 / bootstrapping}\index{bootstrapping / 붓스트랩}\index{likelihood / 가능도}\index{\$COVARIANCE (\$COV, \$COVAR)}

## NONMEM document

NONMEM User’s guide와 help file에는 \$COV를 다음과 같이 설명한다. [@nonmem]\index{\$COVARIANCE (\$COV, \$COVAR)}

```
$COVARIANCE

 DISCUSSION:
 Optional.  Requests  that  the  NONMEM Covariance Step be implemented.
 This step outputs: standard errors, covariance matrix, inverse covari
 ance  matrix,  and the correlation form of the covariance matrix.

 OPTIONS:

 MATRIX=c
      Specifies that the covariance matrix will be different  from  the
      default  (R  sup  -1 S R sup -1).  MATRIX=R requests that 2 times
      the inverse R matrix be used.  MATRIX=S requests that 4 times the
      inverse  S  matrix be used.  (R and S are two matrices from  sta-
      tistical theory,  the Hessian and Cross-Product  Gradient  matri-
      ces,  respectively.)   With  MATRIX=R the standard errors will be
      more consistent with other nonlinear regression  software  imple-
      mentations.

 PRINT=[E][R][S]
      Additional outputs will be  printed  besides  the  defaults.   E:
      print  the  eigenvalues  of the correlation matrix.  R: print the
      matrix .5*R.  S: print the matrix .25*S.

 REFERENCES: Guide IV, section III.B.15
 REFERENCES: Guide V, section 9.4.2 , 10.6
```

\index{standard errors of the estimates / 추정값의 표준오차}\index{standard errors of the estimates / 추정값의 표준오차}
\index{standard errors of the estimates / 추정값의 표준오차}\index{standard errors of the estimates / 추정값의 표준오차}
\index{R matrix / R 행렬}\index{R matrix / R 행렬}
\index{correlation matrix / 상관행렬}\index{eigenvalues / 고유값}

## 이론적 배경 - 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 \left\{\frac{\partial^{2}}{\partial\theta^{2}}lnf(x;\theta) + \lbrack \frac{\partial}{\partial\theta}lnf(x;\theta) \rbrack^{2} \right\} 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
\end{split}
(\#eq:diff-diff)
\end{equation}

\begin{equation}
\begin{split}
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-2)
\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 행렬}\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 행렬}\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}

## R에서의 구현

\@ref(cov-example)의 예제를 R로 풀이하면 다음과 같다. 우선 CRAN에 올라와 있는 nmw package를 설치한다. [@R-nmw]

```r
if (!require(nmw)) install.packages("nmw")
```

`nmw::CovStep()` 함수를 보면 covariance step의 전반적인 개요를 알 수 있다.\index{covariance / 공분산}

```r
# nmw::CovStep
function () 
{
  e$STEP = "COV"
  StartTime = Sys.time()
  Rmat = hessian(Obj, e$FinalPara)/2
  Smat = CalcSmat()
  if (is.nan(Smat[1, 1])) 
    stop("Error - NaN produced")
  invR = solve(Rmat)
  Cov = invR %*% Smat %*% invR
  SE = sqrt(diag(Cov))
  Correl = cov2cor(Cov)
  InvCov = Rmat %*% solve(Smat) %*% Rmat
  EigenVal = sort(eigen(Correl)$values)
  RunTime = difftime(Sys.time(), StartTime)
  Result = list(RunTime, SE, Cov, Correl, InvCov, EigenVal, 
                Rmat, Smat)
  names(Result) = list("Time", "Standard Error", 
                       "Covariance Matrix of Estimates", 
                       "Correlation Matrix of Estimates", 
                       "Inverse Covariance Matrix of Estimates", 
                       "Eigen Values", 
                       "R Matrix", "S Matrix")
  return(Result)
}
```

S matrix를 계산하는 internal obj functions 함수 `nmw::CalcSmat()`는 다음과 같다.\index{S matrix / S 행렬}\index{S matrix / S 행렬}

```r
# nmw::CalcSmat
function () 
{
  Smat = matrix(rep(0, e$nPara * e$nPara), nrow = e$nPara, 
                ncol = e$nPara)
  for (i in 1:e$nID) {
    e$DATAi = e$DataRef[[i]]
    e$ETAi = e$EBE[i, e$EtaNames]
    e$nReci = e$Oi[i, "nRec"]
    if (e$METHOD == "ZERO") {
      gr = grad(OiS0, e$FinalPara)
    }
    else {
      gr = grad(OiS1, e$FinalPara)
    }
    Smat = Smat + gr %*% t(gr)
  }
  return(Smat/4)
}
```

## Theophylline 예제 데이터셋 {#Theoph}

R에는 기본적으로 nlme 패키지[@R-nlme]가 설치되어 있으며, NONMEM에도 있는 theophylline 데이터가 포함되어 있다. Theophylline 데이터를 그리면 다음과 같다.

(ref:theophconctime) Theophylline 경구 복용 후 농도-시간 곡선

```{r theophconctime, fig.cap="(ref:theophconctime)", out.width = '100%'}
require(nlme)
plot(Theoph)
```
위의 자료를 다음과 같이 경구 1 구획모델에 적합시켜 보기로 한다.
\begin{equation}
\begin{split}
  K_{a} & = \theta_{1}e^{\eta_{1}} \\
  V & = \theta_{2}e^{\eta_{2}} \\
  K & = \theta_{3}e^{\eta_{3}} \\
  C(t) & = \frac{\text{DOSE}}{V}\frac{K_{a}}{K_{a} - K}(e^{- Kt} - e^{- K_{a}t})
\end{split}
(\#eq:onecompfit)
\end{equation}

NONMEM 제어구문 파일의 내용은 \@ref(cov-example)에서 표시한 것과 같다.

R에서는 다음과 같은 준비 작업이 필요하다.
\small
```r
DataAll = Theoph
colnames(DataAll) = c("ID", "BWT", "DOSE", "TIME", "DV")
DataAll[,"ID"] = as.numeric(as.character(DataAll[,"ID"]))

nTheta = 3
nEta = 3
nEps = 2

THETAinit = c(2, 50, 0.1)
OMinit = matrix(c(0.2, 0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1, 0.2), 
                nrow=nEta, ncol=nEta)
SGinit = diag(c(0.1, 0.1))

LB = rep(0, nTheta) # Lower bound
UB = rep(1000000, nTheta) # Upper bound

FGD1 = deriv(~DOSE/(TH2*exp(ETA2))*TH1*exp(ETA1)/
             (TH1*exp(ETA1) - TH3*exp(ETA3))*
			 (exp(-TH3*exp(ETA3)*TIME)-exp(-TH1*exp(ETA1)*TIME)),
             c("ETA1","ETA2","ETA3"),
             function.arg=c("TH1", "TH2", "TH3", "ETA1", "ETA2", "ETA3", 
			                "DOSE", "TIME"),
             func=TRUE, hessian=FALSE)

H = deriv(~F + F*EPS1 + EPS2, c("EPS1", "EPS2"), 
          function.arg=c("F", "EPS1", "EPS2"), 
          func=TRUE)

PRED1 = function(THETA, ETA, DATAi)  # for FO and FOCE
{
  FGDres = FGD1(THETA[1], THETA[2], THETA[3], ETA[1], ETA[2], ETA[3], 
                DOSE=320, DATAi[,"TIME"])
  Gres = attr(FGDres, "gradient")
  Hres = attr(H(FGDres, 0, 0), "gradient")
  Res = cbind(FGDres, Gres, Hres)
  colnames(Res) = c("F", "G1", "G2", "G3", "H1", "H2")
  return(Res)
}
```
\normalsize

\index{1차 조건부 추정방법(FOCE) / first-order conditional estimation method(FOCE)}\index{1차추정법 / first-order method(FO)}\index{first-order conditional estimation method(FOCE) / 1차 조건부 추정방법(FOCE)}\index{first-order method(FO) / 1차추정법}\index{first-order conditional estimation method(FOCE) / 1차 조건부 추정방법(FOCE)}\index{first-order method(FO) / 1차추정법}\index{FOCE}\index{FO}

위 PRED 함수의 return 값 중에 F는 model prediction value이고, G들은 model predictiion function(structural model)을 각 η에 대하여 편미분한 값들이다. H들은 잔차에 대한 통계모형(model prediction에 noise가 섞여 관찰값을 이루는 모형)을 각 ε으로 편미분한 값들이다. F, G, H는 목적함수값을 구할 때 필요하다. NONMEM에서는 **nmtran**이 G와 H를 구하는 수식을 기호미분(symbolic differentiation)으로 구하나, R에서는 자동으로 해주는 것이 없으므로, 위와 같이 사용자가 직접 구성해 주어야 한다.\index{목적함수 / objective function}\index{objective function / 목적함수}\index{잔차 / residual error}\index{residual error / 잔차}\index{PRED}

이후에 다음과 같은 R script를 실행하면 NONMEM output과 형태는 다르지만 내용은 동일한 결과를 얻을 수 있다.\index{NONMEM 보고 파일 / NONMEM report file}\index{NONMEM report file / NONMEM 보고 파일}\index{NONMEM report file / NONMEM 보고 파일}
\small
```r
library(nmw)
# First Order Approximation Method
InitStep(DataAll, THETAinit=THETAinit, OMinit=OMinit, SGinit=SGinit, 
         LB=LB, UB=UB, Pred=PRED1, METHOD="ZERO")
(EstRes1 = EstStep())
## $`Initial OFV`
## [1] 141.3076
## 
## $Time
## Time difference of 2.999814 secs
## 
## $Optim
## $Optim$par
##  [1]  0.560417594 -0.167835388  0.148962362  0.995143048  0.056166719
##  [6]  0.151227211 -1.032468525  0.005776729  0.110936464 -0.956899772
## [11] -0.205559310
## 
## $Optim$value
## [1] 57.32106
## 
## $Optim$counts
## function gradient 
##       74       74 
## 
## $Optim$convergence
## [1] 0
## 
## $Optim$message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## 
## 
## $`Final Estimates`
##  [1]  3.16946754 38.25213460  0.10501808  1.19823325  0.13747849
##  [6]  0.03134899  0.37015671  0.04340042  0.25068582  0.01207782
## [11]  0.05427434
```
\normalsize
`InitStep()`과 `EstStep()`을 마친 후 `CovStep()`을 실행하면 아래와 같은 표준오차와 분산 행렬을 계산할 수 있다. 

\small
```r
(CovRes1 = CovStep())
## $Time
## Time difference of 1.0468 secs
## 
## $`Standard Error`
##  [1] 0.641076544 1.685217844 0.023072024 0.420617306 0.082197497
##  [6] 0.019812976 0.340273208 0.023052142 0.289524327 0.003576926
## [11] 0.032078283
## 
## $`Covariance Matrix of Estimates`
##                [,1]         [,2]          [,3]          [,4]          [,5]
##  [1,]  0.4109791347  0.339158144  5.746694e-03  0.2058089645  2.003772e-03
##  [2,]  0.3391581437  2.839959182  5.032613e-03  0.3376028346  3.490465e-02
##  [3,]  0.0057466939  0.005032613  5.323183e-04  0.0016294512 -1.041991e-03
##  [4,]  0.2058089645  0.337602835  1.629451e-03  0.1769189182  1.951490e-02
##  [5,]  0.0020037724  0.034904655 -1.041991e-03  0.0195149026  6.756428e-03
##  [6,] -0.0021925236  0.012804811 -2.503963e-04  0.0032072246  1.504690e-03
##  [7,]  0.1215890847  0.149089319  7.111900e-03  0.0575731487 -1.010198e-02
##  [8,]  0.0009971098  0.023865634  6.271266e-05  0.0042158445  8.584714e-04
##  [9,]  0.0669924083  0.057326151  6.226096e-03  0.0179862543 -1.309239e-02
## [10,]  0.0010500117  0.001807746  5.805488e-05  0.0005143569 -7.516774e-05
## [11,] -0.0049728997 -0.009950377 -4.790610e-04 -0.0010145003  9.532948e-04
##                [,6]          [,7]         [,8]          [,9]         [,10]
##  [1,] -2.192524e-03  0.1215890847 9.971098e-04  0.0669924083  1.050012e-03
##  [2,]  1.280481e-02  0.1490893190 2.386563e-02  0.0573261514  1.807746e-03
##  [3,] -2.503963e-04  0.0071119003 6.271266e-05  0.0062260963  5.805488e-05
##  [4,]  3.207225e-03  0.0575731487 4.215844e-03  0.0179862543  5.143569e-04
##  [5,]  1.504690e-03 -0.0101019780 8.584714e-04 -0.0130923877 -7.516774e-05
##  [6,]  3.925540e-04 -0.0028272756 2.326326e-04 -0.0032697941 -2.051327e-05
##  [7,] -2.827276e-03  0.1157858558 3.116262e-03  0.0940102394  9.767199e-04
##  [8,]  2.326326e-04  0.0031162617 5.314013e-04  0.0018656807  2.786064e-05
##  [9,] -3.269794e-03  0.0940102394 1.865681e-03  0.0838243357  8.055388e-04
## [10,] -2.051327e-05  0.0009767199 2.786064e-05  0.0008055388  1.279440e-05
## [11,]  1.806783e-04 -0.0038608273 2.199601e-04 -0.0033970159 -2.824858e-05
##               [,11]
##  [1,] -4.972900e-03
##  [2,] -9.950377e-03
##  [3,] -4.790610e-04
##  [4,] -1.014500e-03
##  [5,]  9.532948e-04
##  [6,]  1.806783e-04
##  [7,] -3.860827e-03
##  [8,]  2.199601e-04
##  [9,] -3.397016e-03
## [10,] -2.824858e-05
## [11,]  1.029016e-03
## 
## $`Correlation Matrix of Estimates`
##              [,1]       [,2]       [,3]        [,4]        [,5]       [,6]
##  [1,]  1.00000000  0.3139325  0.3885281  0.76325079  0.03802594 -0.1726174
##  [2,]  0.31393253  1.0000000  0.1294350  0.47628061  0.25198153  0.3835018
##  [3,]  0.38852814  0.1294350  1.0000000  0.16790689 -0.54943908 -0.5477629
##  [4,]  0.76325079  0.4762806  0.1679069  1.00000000  0.56444374  0.3848509
##  [5,]  0.03802594  0.2519815 -0.5494391  0.56444374  1.00000000  0.9239295
##  [6,] -0.17261745  0.3835018 -0.5477629  0.38485092  0.92392947  1.0000000
##  [7,]  0.55738714  0.2599936  0.9058832  0.40225837 -0.36117699 -0.4193635
##  [8,]  0.06747173  0.6143355  0.1179121  0.43479661  0.45306025  0.5093422
##  [9,]  0.36093637  0.1174929  0.9320626  0.14769593 -0.55014251 -0.5700142
## [10,]  0.45790382  0.2998965  0.7034659  0.34187510 -0.25566008 -0.2894510
## [11,] -0.24181804 -0.1840655 -0.6472826 -0.07518893  0.36154098  0.2842792
##             [,7]       [,8]       [,9]      [,10]       [,11]
##  [1,]  0.5573871 0.06747173  0.3609364  0.4579038 -0.24181804
##  [2,]  0.2599936 0.61433553  0.1174929  0.2998965 -0.18406548
##  [3,]  0.9058832 0.11791205  0.9320626  0.7034659 -0.64728263
##  [4,]  0.4022584 0.43479661  0.1476959  0.3418751 -0.07518893
##  [5,] -0.3611770 0.45306025 -0.5501425 -0.2556601  0.36154098
##  [6,] -0.4193635 0.50934216 -0.5700142 -0.2894510  0.28427925
##  [7,]  1.0000000 0.39727833  0.9542504  0.8024764 -0.35370524
##  [8,]  0.3972783 1.00000000  0.2795381  0.3378856  0.29745513
##  [9,]  0.9542504 0.27953807  1.0000000  0.7778421 -0.36576437
## [10,]  0.8024764 0.33788563  0.7778421  1.0000000 -0.24619292
## [11,] -0.3537052 0.29745513 -0.3657644 -0.2461929  1.00000000
## 
## $`Inverse Covariance Matrix of Estimates`
##              [,1]         [,2]        [,3]        [,4]        [,5]
##  [1,]   106.16085    -68.57396    6449.005    335.8698   -2554.409
##  [2,]   -68.57396     58.03937   -4878.746   -302.1420    2175.297
##  [3,]  6449.00514  -4878.74594  589180.809  26966.6055 -188642.065
##  [4,]   335.86981   -302.14199   26966.605   1681.5577  -11681.346
##  [5,] -2554.40932   2175.29716 -188642.065 -11681.3456   84767.297
##  [6,]  -386.87894    570.22260  -66147.099  -3404.8900   13635.511
##  [7,] -1202.16352    939.99684  -90186.464  -5086.8917   35747.140
##  [8,] 10794.57609  -8973.04621  795473.397  47387.2333 -336778.082
##  [9,]   -49.38187     87.68163  -10522.263   -442.6127    3308.451
## [10,] 11656.77324 -10122.84537  899033.055  53311.6422 -378718.161
## [11,] -1043.11500   1001.74635  -47225.438  -4879.5431   35063.038
##               [,6]         [,7]        [,8]         [,9]      [,10]
##  [1,]    -386.8789   -1202.1635   10794.576    -49.38187   11656.77
##  [2,]     570.2226     939.9968   -8973.046     87.68163  -10122.85
##  [3,]  -66147.0986  -90186.4639  795473.397 -10522.26321  899033.06
##  [4,]   -3404.8900   -5086.8917   47387.233   -442.61268   53311.64
##  [5,]   13635.5106   35747.1396 -336778.082   3308.45066 -378718.16
##  [6,]   72186.1449   10923.7488 -116902.668   2827.92008 -138707.39
##  [7,]   10923.7488   16640.0641 -149635.854    965.72182 -166637.08
##  [8,] -116902.6684 -149635.8536 1416416.077 -14025.69870 1587796.18
##  [9,]    2827.9201     965.7218  -14025.699    954.65511  -20047.21
## [10,] -138707.3931 -166637.0784 1587796.183 -20047.20949 2031529.82
## [11,]   15687.7641   14275.7793 -151936.736    935.29881 -170271.34
##              [,11]
##  [1,]   -1043.1150
##  [2,]    1001.7464
##  [3,]  -47225.4381
##  [4,]   -4879.5431
##  [5,]   35063.0376
##  [6,]   15687.7641
##  [7,]   14275.7793
##  [8,] -151936.7362
##  [9,]     935.2988
## [10,] -170271.3406
## [11,]   28036.5550
## 
## $`Eigen Values`
##  [1] 0.0002519304 0.0096729015 0.0108358602 0.0233184643 0.0520725533
##  [6] 0.2982375053 0.5047779131 0.9114702297 1.2088053283 3.2082379737
## [11] 4.7723193401
## 
## $`R Matrix`
##              [,1]         [,2]         [,3]         [,4]         [,5]
##  [1,]   17.924787   -1.3343223  -162.767654   -4.1309683    21.546405
##  [2,]   -1.334322    0.5507357    -7.672315    0.1118322    -1.462878
##  [3,] -162.767654   -7.6723148 34333.363150   86.0269293   433.962384
##  [4,]   -4.130968    0.1118322    86.026929   28.6263094  -177.270130
##  [5,]   21.546405   -1.4628778   433.962384 -177.2701302  1930.445843
##  [6,]   10.225928  -16.5210396    13.387686  272.9370786 -4270.878832
##  [7,]  -11.022690    2.9849069   -90.741373  -52.9261900   210.709300
##  [8,]   52.304346  -18.2457139   956.482064  164.3158075 -1421.957500
##  [9,]    7.044855   -2.2338946 -1350.939646   24.4536958   -43.763546
## [10,]  248.456482 -120.7991176 -7033.212482   50.2328789 -1013.856688
## [11,]   -1.752135   -5.2052276 -1992.414213    6.0120604   124.417556
##              [,6]        [,7]        [,8]         [,9]        [,10]
##  [1,]    10.22593  -11.022690    52.30435     7.044855    248.45648
##  [2,]   -16.52104    2.984907   -18.24571    -2.233895   -120.79912
##  [3,]    13.38769  -90.741373   956.48206 -1350.939646  -7033.21248
##  [4,]   272.93708  -52.926190   164.31581    24.453696     50.23288
##  [5,] -4270.87883  210.709300 -1421.95750   -43.763546  -1013.85669
##  [6,] 16610.43942 -139.814385  1113.59904    18.726078   4680.59998
##  [7,]  -139.81438  213.228947  -555.99366  -151.083275     96.25915
##  [8,]  1113.59904 -555.993663  4043.51428   130.794770   -555.76917
##  [9,]    18.72608 -151.083275   130.79477   236.875935    -20.42601
## [10,]  4680.59998   96.259149  -555.76917   -20.426010 192857.05263
## [11,]   -46.02961  -62.941133  -201.26760    92.656857   6568.90926
##              [,11]
##  [1,]    -1.752135
##  [2,]    -5.205228
##  [3,] -1992.414213
##  [4,]     6.012060
##  [5,]   124.417556
##  [6,]   -46.029614
##  [7,]   -62.941133
##  [8,]  -201.267605
##  [9,]    92.656857
## [10,]  6568.909257
## [11,]  3974.804398
## 
## $`S Matrix`
##               [,1]         [,2]         [,3]         [,4]        [,5]
##  [1,]    78.316509   -4.6468525  -1295.13192   -11.873085   142.72165
##  [2,]    -4.646852    0.7648878     64.36589     2.623533   -28.61925
##  [3,] -1295.131915   64.3658917 183632.39790  -230.636173   840.38211
##  [4,]   -11.873085    2.6235332   -230.63617    18.368716  -171.71679
##  [5,]   142.721653  -28.6192545    840.38211  -171.716794  2005.81552
##  [6,]  -145.835176   29.4905947   9000.10289   291.779615 -3809.95407
##  [7,]   -26.707401    0.2387057   3794.27704   -19.686952    51.76139
##  [8,]    44.375129   10.7614124 -10813.66435    84.841787  -765.19107
##  [9,]    13.946014   -4.4042212  -6396.75146     3.480210    87.90129
## [10,]  2039.647982 -397.4745827  -4148.02643 -1170.279733  8916.77585
## [11,]   279.500822  -47.3111189 -60483.51062   -22.729230   670.78875
##              [,6]          [,7]         [,8]         [,9]       [,10]
##  [1,]  -145.83518   -26.7074010     44.37513    13.946014   2039.6480
##  [2,]    29.49059     0.2387057     10.76141    -4.404221   -397.4746
##  [3,]  9000.10289  3794.2770370 -10813.66435 -6396.751456  -4148.0264
##  [4,]   291.77961   -19.6869516     84.84179     3.480210  -1170.2797
##  [5,] -3809.95407    51.7613883   -765.19107    87.901295   8916.7758
##  [6,] 12023.28652   188.5688359    667.62858  -711.894529  -3829.1367
##  [7,]   188.56884   129.3349739   -292.66398  -155.764410   1796.9713
##  [8,]   667.62858  -292.6639799   1121.03185   294.247258 -10631.8774
##  [9,]  -711.89453  -155.7644099    294.24726   327.282119   1812.2113
## [10,] -3829.13667  1796.9713202 -10631.87742  1812.211286 419517.6543
## [11,] -3489.01512 -1105.9231026   2773.71160  2358.454995  18067.4267
##              [,11]
##  [1,]    279.50082
##  [2,]    -47.31112
##  [3,] -60483.51062
##  [4,]    -22.72923
##  [5,]    670.78875
##  [6,]  -3489.01512
##  [7,]  -1105.92310
##  [8,]   2773.71160
##  [9,]   2358.45500
## [10,]  18067.42672
## [11,]  24042.66052
```
\normalsize

`PostHocEta()`, `TabStep()`의 실행화면은 아래와 같다.

\index{ID}
```r
(EBE1 = PostHocEta()) # Using e$FinalPara from EstStep()
##       ID       ETA1         ETA2        ETA3
##  [1,]  1 -0.6367109 -0.232258352 -0.73648224
##  [2,]  2 -0.5895843 -0.153341805 -0.06619115
##  [3,]  3 -0.3083755 -0.124816676 -0.21013190
##  [4,]  4 -1.0305984 -0.186821177 -0.21195510
##  [5,]  5 -0.8235560 -0.302352128 -0.24453948
##  [6,]  6 -1.0025271  0.068181532 -0.08745089
##  [7,]  7 -1.4316285 -0.097903076 -0.13802639
##  [8,]  8 -0.7541785 -0.039239022 -0.19621190
##  [9,]  9  0.7875803  0.010757282 -0.19937965
## [10,] 10 -1.4555649 -0.369057237 -0.40057582
## [11,] 11  0.1541451 -0.005061315 -0.08005791
## [12,] 12 -1.2863346 -0.388864841 -0.10134440
# (Tab1 = TabStep())
```

## 고유값(Eigenvalue) {#eigenvalue}
\index{고유값 / eigenvalues}\index{eigenvalues / 고유값}

공분산단계에서는 상관행렬(correlation matrix)의 고유값(eigenvalue)도 출력된다.\index{고유값 / eigenvalues}\index{공분산 / covariance}\index{상관행렬 / correlation matrix}\index{correlation matrix / 상관행렬}\index{covariance / 공분산}\index{eigenvalues / 고유값}

어떤 행렬 A를 스칼라 \(\lambda\)와 벡터u를 이용하여 다음과 같이 표시할 수 있을 때,

\begin{equation}
Au = \lambda u
(\#eq:au-lu)
\end{equation}

\(\lambda\)를 A의 고유값이라 하고 그 고유값에 해당하는 벡터 u를 고유벡터(eigenvector)라고 한다.\index{고유값 / eigenvalues}\index{eigenvalues / 고유값}

\(\lambda\)는

\begin{equation}
\left| A - \lambda I \right| = 0
(\#eq:a-li-zero)
\end{equation}

이라는 특성 방정식의 근들로 구할 수 있다.

Minimization에서 고유값을 출력하는 이유는 minimization 과정이 얼마나 어려웠는가를 표시하기 위해서이다. 조건수(condition number)란 고유값의 최대치/최소치에다 1/2승한 것인데, 이것이 크다면 그만큼 수치계산의 어려움을 겪었다는 것이다. 조건수가 크다고 모형이 잘못되었다는 것을 의미하는 것은 아니다. 혹시 minimization(fitting)에 실패하였다면 그 원인 중의 하나가 조건수가 크기 때문일 수 있다.\index{고유값 / eigenvalues}\index{eigenvalues / 고유값}\index{조건수 / condition number}\index{condition number / 조건수}

## 결론

NONMEM의 Covariance step을 이해하기 위해 R로 NONMEM의 output을 재현해 보았다. 그러나, R이 NONMEM을 대체할 수 없는 이유 두 가지는 다음과 같다.

1.  NMTRAN의 역할을 해줄 함수들이 없다.
2.  속도가 적어도 수 십에서 수 백 배 느리다.

나머지 FOCE나 Laplacian에 대해서는 nmw 패키지의 도움말^[https://cran.r-project.org/web/packages/nmw/nmw.pdf]을 보고 실행하면 알 수 있다. [@R-nmw]\index{1차 조건부 추정방법(FOCE) / first-order conditional estimation method(FOCE)}\index{first-order conditional estimation method(FOCE) / 1차 조건부 추정방법(FOCE)}\index{first-order conditional estimation method(FOCE) / 1차 조건부 추정방법(FOCE)}\index{FOCE}\index{FO}

<!--
참고문헌

1.  NONMEM User’s Guide I - VIII
2.  김우철. 수리통계학. 민영사. 2012
-->