16.1 FO Method
FO method 결과를 R에서 재현하여 보면 다음과 같다.
16.1.1 초기화 단계 (Initialization Step)
먼저, \(\theta\), \(\eta\), \(\epsilon\)의 개수와 초기값, lower bound (LB), upper bound (UB), prediction function (PRED)등을 지정해야 한다.
PRED function은 구조 모형(약동학 모형)에 의한 예측값(\(F\))만 return하는 것이 아니라, 아래 G1, G2, G3로 표현된 \(\frac{\partial F}{\partial \eta_i}\)와 아래 H1, H2로 표현된 \(\frac{\partial Y}{\partial \epsilon_i}\)를 함께 return해야 한다.
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
FGD = 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=TRUE)
H = deriv(~F + F*EPS1 + EPS2, c("EPS1", "EPS2"), function.arg=c("F", "EPS1", "EPS2"),
func=TRUE)
PRED = function(THETA, ETA, DATAi)
{
FGDres = FGD(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")
if (e$METHOD == "LAPL") {
Dres = attr(FGDres, "hessian")
Res = cbind(FGDres, Gres, Hres, Dres[,1,1], Dres[,2,1], Dres[,2,2], Dres[,3,])
colnames(Res) = c("F", "G1", "G2", "G3", "H1", "H2", "D11", "D21", "D22",
"D31", "D32", "D33")
} else {
Res = cbind(FGDres, Gres, Hres)
colnames(Res) = c("F", "G1", "G2", "G3", "H1", "H2")
}
return(Res)
}위에서 PRED function을 만들 때 FO method 뿐 아니라, FOCE와 LAPL에 대비하여 만들었다. R에 내장된 deriv 함수는 첫 번째 인자로 주어진 수식을 두 번째 인자로 주어진 변수에 대해 기호 편미분(symbolic differentiation)하고, 1차 및 2차 편미분값을 돌려주는 함수를 만들어준다.
require(nmw)
InitStep(DATA, THETAinit=THETAinit, OMinit=OMinit, SGinit=SGinit, LB=LB, UB=UB,
Pred=PRED, METHOD="ZERO")위에서 METHOD에 “COND”나 “LAPL”로서 추정 방법을 다르게 지정할 수 있다. INTERACTION 옵션은 따로 지정할 필요 없이 기본적으로 포함된다.
16.1.2 추정 단계 (Estimation Step)
추정단계에서는 \(\theta\) , \(\Omega\), \(\Sigma\) 를 점추정(point estimation)하는 단계이다.
(EstRes = EstStep())$`Initial OFV`
[1] 141
$Time
Time difference of 1.27 secs
$Optim
$Optim$par
[1] 0.56042 -0.16784 0.14896 0.99514 0.05617 0.15123 -1.03247 0.00578 0.11094 -0.95690
[11] -0.20556
$Optim$value
[1] 57.3
$Optim$counts
function gradient
74 74
$Optim$convergence
[1] 0
$Optim$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
$`Final Estimates`
[1] 3.1695 38.2521 0.1050 1.1982 0.1375 0.0313 0.3702 0.0434 0.2507 0.0121 0.054316.1.3 공분산 단계 (Covariance Step)
Covariance step에서는 앞에서 구한 점추정치들의 분산-공분산 행렬을 추정하여, 이것으로 표준오차(standard error)를 구하게 된다. 선형 회귀나 단순한 이차 곡선 회귀이면 적합된 목적함수에서 바로 표준오차(SE)까지 계산이 가능하지만, 이번 비선형 회귀에서는 변수변환된 파라미터(scaled and transformed parameter, STP)를 사용하였기 때문에, 바로 표준오차(SE)를 구하는 것이 불가능하여 최종 추정값에서 다시 Hessian을 구하여 표준오차(SE)를 구한다.
(CovRes = CovStep())$Time
Time difference of 0.437 secs
$`Standard Error`
[1] 0.64108 1.68522 0.02307 0.42062 0.08220 0.01981 0.34027 0.02305 0.28952 0.00358 0.03208
$`Covariance Matrix of Estimates`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.410979 0.33916 5.75e-03 0.205809 2.00e-03 -2.19e-03 0.121589 9.97e-04 0.066992
[2,] 0.339158 2.83996 5.03e-03 0.337603 3.49e-02 1.28e-02 0.149089 2.39e-02 0.057326
[3,] 0.005747 0.00503 5.32e-04 0.001629 -1.04e-03 -2.50e-04 0.007112 6.27e-05 0.006226
[4,] 0.205809 0.33760 1.63e-03 0.176919 1.95e-02 3.21e-03 0.057573 4.22e-03 0.017986
[5,] 0.002004 0.03490 -1.04e-03 0.019515 6.76e-03 1.50e-03 -0.010102 8.58e-04 -0.013092
[6,] -0.002193 0.01280 -2.50e-04 0.003207 1.50e-03 3.93e-04 -0.002827 2.33e-04 -0.003270
[7,] 0.121589 0.14909 7.11e-03 0.057573 -1.01e-02 -2.83e-03 0.115786 3.12e-03 0.094010
[8,] 0.000997 0.02387 6.27e-05 0.004216 8.58e-04 2.33e-04 0.003116 5.31e-04 0.001866
[9,] 0.066992 0.05733 6.23e-03 0.017986 -1.31e-02 -3.27e-03 0.094010 1.87e-03 0.083824
[10,] 0.001050 0.00181 5.81e-05 0.000514 -7.52e-05 -2.05e-05 0.000977 2.79e-05 0.000806
[11,] -0.004973 -0.00995 -4.79e-04 -0.001015 9.53e-04 1.81e-04 -0.003861 2.20e-04 -0.003397
# Column 10, 11 not shown
$`Correlation Matrix of Estimates`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 1.0000 0.314 0.389 0.7633 0.038 -0.173 0.557 0.0675 0.361 0.458 -0.2418
[2,] 0.3139 1.000 0.129 0.4763 0.252 0.384 0.260 0.6143 0.117 0.300 -0.1841
[3,] 0.3885 0.129 1.000 0.1679 -0.549 -0.548 0.906 0.1179 0.932 0.703 -0.6473
[4,] 0.7633 0.476 0.168 1.0000 0.564 0.385 0.402 0.4348 0.148 0.342 -0.0752
[5,] 0.0380 0.252 -0.549 0.5644 1.000 0.924 -0.361 0.4531 -0.550 -0.256 0.3615
[6,] -0.1726 0.384 -0.548 0.3849 0.924 1.000 -0.419 0.5093 -0.570 -0.289 0.2843
[7,] 0.5574 0.260 0.906 0.4023 -0.361 -0.419 1.000 0.3973 0.954 0.802 -0.3537
[8,] 0.0675 0.614 0.118 0.4348 0.453 0.509 0.397 1.0000 0.280 0.338 0.2975
[9,] 0.3609 0.117 0.932 0.1477 -0.550 -0.570 0.954 0.2795 1.000 0.778 -0.3658
[10,] 0.4579 0.300 0.703 0.3419 -0.256 -0.289 0.802 0.3379 0.778 1.000 -0.2462
[11,] -0.2418 -0.184 -0.647 -0.0752 0.362 0.284 -0.354 0.2975 -0.366 -0.246 1.0000
$`Inverse Covariance Matrix of Estimates`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 106.2 -68.6 6449 336 -2554 -387 -1202 10795 -49.4 11657 -1043
[2,] -68.6 58.0 -4879 -302 2175 570 940 -8973 87.7 -10123 1002
[3,] 6449.0 -4878.7 589181 26967 -188642 -66147 -90186 795473 -10522.3 899033 -47225
[4,] 335.9 -302.1 26967 1682 -11681 -3405 -5087 47387 -442.6 53312 -4880
[5,] -2554.4 2175.3 -188642 -11681 84767 13636 35747 -336778 3308.5 -378718 35063
[6,] -386.9 570.2 -66147 -3405 13636 72186 10924 -116903 2827.9 -138707 15688
[7,] -1202.2 940.0 -90186 -5087 35747 10924 16640 -149636 965.7 -166637 14276
[8,] 10794.6 -8973.0 795473 47387 -336778 -116903 -149636 1416416 -14025.7 1587796 -151937
[9,] -49.4 87.7 -10522 -443 3308 2828 966 -14026 954.7 -20047 935
[10,] 11656.8 -10122.8 899033 53312 -378718 -138707 -166637 1587796 -20047.2 2031530 -170271
[11,] -1043.1 1001.7 -47225 -4880 35063 15688 14276 -151937 935.3 -170271 28037
$`Eigen Values`
[1] 0.000252 0.009673 0.010836 0.023318 0.052073 0.298238 0.504778 0.911470 1.208805 3.208238
[11] 4.772319
$`R Matrix`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 17.92 -1.334 -162.77 -4.131 21.55 10.2 -11.02 52.3 7.04 248.5
[2,] -1.33 0.551 -7.67 0.112 -1.46 -16.5 2.98 -18.2 -2.23 -120.8
[3,] -162.77 -7.672 34333.36 86.027 433.96 13.4 -90.74 956.5 -1350.94 -7033.2
[4,] -4.13 0.112 86.03 28.626 -177.27 272.9 -52.93 164.3 24.45 50.2
[5,] 21.55 -1.463 433.96 -177.270 1930.45 -4270.9 210.71 -1422.0 -43.76 -1013.9
[6,] 10.23 -16.521 13.39 272.937 -4270.88 16610.4 -139.81 1113.6 18.73 4680.6
[7,] -11.02 2.985 -90.74 -52.926 210.71 -139.8 213.23 -556.0 -151.08 96.3
[8,] 52.30 -18.246 956.48 164.316 -1421.96 1113.6 -555.99 4043.5 130.79 -555.8
[9,] 7.04 -2.234 -1350.94 24.454 -43.76 18.7 -151.08 130.8 236.88 -20.4
[10,] 248.46 -120.799 -7033.21 50.233 -1013.86 4680.6 96.26 -555.8 -20.43 192857.1
[11,] -1.75 -5.205 -1992.41 6.012 124.42 -46.0 -62.94 -201.3 92.66 6568.9
# Column 11 not shown
$`S Matrix`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 78.32 -4.647 -1295.1 -11.87 142.7 -145.8 -26.707 44.4 13.95 2040
[2,] -4.65 0.765 64.4 2.62 -28.6 29.5 0.239 10.8 -4.40 -397
[3,] -1295.13 64.366 183632.4 -230.64 840.4 9000.1 3794.277 -10813.7 -6396.75 -4148
[4,] -11.87 2.624 -230.6 18.37 -171.7 291.8 -19.687 84.8 3.48 -1170
[5,] 142.72 -28.619 840.4 -171.72 2005.8 -3810.0 51.761 -765.2 87.90 8917
[6,] -145.84 29.491 9000.1 291.78 -3810.0 12023.3 188.569 667.6 -711.89 -3829
[7,] -26.71 0.239 3794.3 -19.69 51.8 188.6 129.335 -292.7 -155.76 1797
[8,] 44.38 10.761 -10813.7 84.84 -765.2 667.6 -292.664 1121.0 294.25 -10632
[9,] 13.95 -4.404 -6396.8 3.48 87.9 -711.9 -155.764 294.2 327.28 1812
[10,] 2039.65 -397.475 -4148.0 -1170.28 8916.8 -3829.1 1796.971 -10631.9 1812.21 419518
[11,] 279.50 -47.311 -60483.5 -22.73 670.8 -3489.0 -1105.923 2773.7 2358.45 18067
# Column 11 not shown이 단계의 첫 번째 부분(section)은 추정치의 표준오차로 이전 단계의 FINAL PARAMETER ESTIMATE와 배치형태가 동일하다. 이후의 출력들은 대개 이것을 구하는데 필요한 중간 계산 결과들이다.
두 번째 부분은 추정값(estimate)의 분산-공분산 행렬(variance-covariance matrix)이다. 이 행렬의 대각원소에 제곱근을 취한 것이 추정값의 표준오차이다. 이는 뒤에 나오는 R matrix와 S matrix를 사용하여 \(R^{- 1}SR^{- 1}\)를 계산한 것이다.
세 번째 부분은 직전 부분인 분산-공분산 행렬을 상관 행렬(correlation matrix)로 변환한 것이다. NONMEM의 경우에는 대각원소에 제곱근을 취하여 표준 편차(standard deviation, 여기서는 standard error가 됨)를 표시하고 있다.
네 번째 부분은 분산-공분산 행렬의 역행렬이다. 내부적으로는 이 행렬을 먼저 구한 뒤, 이것의 역행렬을 분산-공분산 행렬로 삼는다. 이것의 의미는 Fisher’s Information Matrix (NONMEM에서의 R matrix)(의 추정량과) 유사하지만, 계산법은 더 보수적이어서 다르다.
다섯 번째 부분은 추정값들의 상관 행렬(correlation matrix of estimates)로부터 구한 고유치(eigenvalue)들을 작은 것부터 순서대로 나열한 것이다. 이런 상관 행렬은 양정치(positive definite)이기 때문에, 고유값들이 모두 양의 실수(positive real values)이다. 또한, 이 행렬은 실수 대칭 행렬(real symmetric matrix)이므로 실수인 고유값(real eigenvalues)을 갖는다. 추정된 상관계수들도 거의 항상 서로 다르다. 정리하자면, 추정값들의 상관 행렬은 full rank, 실수 대칭, 양정치 행렬(full-rank real symmetric positive definite matrix)이다. 만약 그렇지 않다면, NONMEM 등이 사용하는 BFGS 같은 계산 알고리듬(computation algorithm)상 최소화 실패(miniization failure)가 발생하게 된다. 따라서, 이 상관 행렬은 대각원소의 개수만큼 서로 다른 양의 실수 고유값(distinct positive real eigenvalues of rank)을 갖는다.
여섯 번째 부분은 R matrix이다. 이는 log likelihood를 parameter들로 이계 미분한 것(Hessian matrix)으로 Fisher’s Information Matrix(FIM)(의 추정량)이라고도 한다. 일반적인 최대 우도 추정법(maximum likelihood estimation, MLE)에서는 R matrix의 역행렬을 추정치의 분산공분산 행렬로 사용하지만 NONMEM에서는 다음에 나오는 S matrix를 함께 사용한다.
마지막 부분은 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}\)를 추정값의 분산-공분산 행렬로 사용한다.
만약 어떤 원인(예를 들어, R이 numerically singular)에 의해 분산공분산 행렬을 구할 수 없고, 따라서, 표준오차를 구할 수 없다면 $COV에서 MATRIX=R 또는 MATRIX=S 옵션으로 단순화하여 표준오차를 구할 수 있다.
이 경우에는 대개 표준오차가 더 작게 나온다. 어떤 행렬이 singular(비정칙)이어서 역행렬을 구할 수 없다면 NONMEM은 pseudo-inverse matrix를 구하게 되는데, 이는 유일한(unique) 해는 아니므로 해석에 유의해야 한다.
공분산 단계(Covariance step)가 실패하는 경우에는 먼저 파라미터 추정값에 이상이 없는지 살펴보고, 이상이 없어 보이면 붓스트랩이나 profile-likelihood 등의 다른 방법을 이용하여 신뢰구간을 구하면 된다. 이 때 신뢰구간이 비대칭일 가능성이 크므로, 표준 오차를 (산출하였다 하더라도) 제시하지 않는 경우가 많다. 혹자는 covariance step에 의한 표준오차보다 붓스트랩에 의한 신뢰구간이 (시간은 훨씬 많이 걸리지만) 더 좋다고 한다. 붓스트랩을 할 때는 공분산 단계(covariance step)와 표 단계(table step)를 삭제하면 실행 시간이 절약된다.
MLE 이론에 따르면 log likelihood를 파라미터 추정값에서 1계 편미분한 것의 제곱과 2계 편미분한 것의 기대값은 부호만 반대이면서 크기가 같다.
지금까지의 설명을 다시 정리하면 다음과 같다.
아주 이상적인 (그러나 극히 드문) 상황에서는 목적함수의 1계 미분의 제곱(S matrix)과 2계 미분(R matrix)이 같겠지만 (기대값은 같을 수도 있다), 실제 관찰값으로 계산하면 다르다. NONMEM에서는 1계 미분의 제곱행렬(cross product)을 S matrix, 2계 미분 행렬(Hessian)을 R matrix라고 부른다. 대부분의 MLE software에서는 R matrix의 inverse만으로 estimate의 분산-공분산 행렬로 삼고 표준오차를 산출한다. 반면, NONMEM에서는 \(R^{- 1}SR^{- 1}\)을 estimate의 분산-공분산 행렬로 삼는다. 이것은 더 보수적인 방법이어서 잔차(\(\epsilon\))와 interindividual random variability (\(\eta\))의 정규분포가정을 완화했다고 볼 수 있으며, 신뢰구간은 더 넓게 나온다.
수식으로 표현하자면 다음과 같다. 여기에서 \(O_{i}\)는 개인(i)의 목적함수 값이며, \(O\)는 전체 목적함수 값이다.
\[\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} \tag{16.1} \end{equation}\]
위의 설명은 가톨릭대학교 계량약리학연구소(PIPET) 대표저자 임동석의 ’계량약리학 워크샵 초급과정(2020. 주식회사 부크크)’의 14장에도 나와 있다.
16.1.4 Post Hoc Eta
FO method는 추정 과정(estimation step)에서는 모든 \(\eta\)를 0으로 고정하기 때문에 Empirical Bayes Estimate (EBE, 또는 realized eta, maximum a posteriori, post hoc eta)를 따로 구해주어야 한다. 반면 FOCE나 LAPL에서는 매 iteration마다 EBE를 추정하기 때문에, 나중에 따로 EBE를 추정할 필요는 없다. Estimation step에서 “HESSIAN OF POSTERIOR DENSITY IS NON-POSITIVE-DEFINITE DURING SEARCH”라는 error message를 보는 경우가 있는데, 이는 특정 subject에서 EBE추정에 실패했다는 뜻이다.
(EBE = PostHocEta()) # Using e$FinalPara from EstStep() ID ETA1 ETA2 ETA3
[1,] 1 -0.637 -0.23226 -0.7365
[2,] 2 -0.590 -0.15334 -0.0662
[3,] 3 -0.308 -0.12482 -0.2101
[4,] 4 -1.031 -0.18682 -0.2120
[5,] 5 -0.824 -0.30235 -0.2445
[6,] 6 -1.003 0.06818 -0.0875
[7,] 7 -1.432 -0.09790 -0.1380
[8,] 8 -0.754 -0.03924 -0.1962
[9,] 9 0.788 0.01076 -0.1994
[10,] 10 -1.456 -0.36906 -0.4006
[11,] 11 0.154 -0.00506 -0.0801
[12,] 12 -1.286 -0.38886 -0.1013EBE를 얻기 위한 objective function은 별도로 있는데, 다음에서 확인할 수 있다.
ObjEta = function(ETAi)
{
# External Variable : e$INTER, e$DATAi, e$THETA, e$invOM, e$SG, e$nEta, e$HNames
# External Function : e$PRED
FGHDi = e$PRED(e$THETA, ETAi, e$DATAi)
Ri = e$DATAi[,"DV"] - FGHDi[,"F"]
if (e$INTER == TRUE) {
Hi = FGHDi[, e$HNames, drop=FALSE]
} else {
FGHD0 = e$PRED(e$THETA, rep(0, e$nEta), e$DATAi)
Hi = FGHD0[, e$HNames, drop=FALSE]
}
Vi = diag(Hi %*% e$SG %*% t(Hi))
return(sum(log(Vi) + Ri*Ri/Vi) + t(ETAi) %*% e$invOM %*% ETAi)
}INTER option에 따라 V(Y)의 추정이 어떻게 달라지는지 주목해야 한다.
일반적으로 (적어도 NONMEM에서는) \(\theta\), \(\Omega\), \(\Sigma\)를 추정하기 위한 목적 함수(objective function)와 EBE를 추정하기 목적 함수(objective function)를 별도로 두고, 매 iteration마다 각각 최소화하는데, Youngjo Lee, Yudi Pawitan 등은 이 둘을 합해(likelihood 기준으로는 곱해) 하나의 목적 함수(objective function)를 구성하고 h-likelihood (extended likelihood의 특별한 case)라 하였다. (Yudi Pawitan. In All Likelihood. Oxford University Press. 2001)
하나의 objective function으로 구성하면 최적화 시간이 단축된다. 하지만, (연속형 Y 변수를 다루는 실제 자료에서는 일어날 것 같지 않지만) 일부 적용할 수 없는 경우들이 있어 완벅히 일반화할 수는 없다. (Lee Y, Nelder JA, Pawitan Y. Generalized Linear Models with Random Effects. 2e. Chapman & Hall/CRC. 2017)
16.1.5 Table Step
(sdtab = TabStep()) ID TIME DV PRED RES WRES G1 G2 G3 H1 H2
1 1 0.00 0.74 0.000 0.7400 3.1764 0.0000 0.000 0.0000 0.000 1
2 1 0.25 2.84 4.511 -1.6705 -0.1411 2.9495 -4.511 -0.0667 4.511 1
3 1 0.57 6.57 6.729 -0.1587 0.4537 2.3362 -6.729 -0.2572 6.729 1
4 1 1.12 10.50 7.444 3.0565 1.9984 0.6273 -7.444 -0.6497 7.444 1
5 1 2.02 9.66 6.984 2.6759 0.3614 -0.1475 -6.984 -1.2453 6.984 1
6 1 3.82 8.58 5.793 2.7871 -0.2128 -0.1979 -5.793 -2.1254 5.793 1
7 1 5.10 8.36 5.064 3.2957 0.2574 -0.1735 -5.064 -2.5389 5.064 1
8 1 7.03 7.47 4.135 3.3348 0.0668 -0.1417 -4.135 -2.9112 4.135 1
9 1 9.05 6.89 3.345 3.5452 0.5296 -0.1146 -3.345 -3.0643 3.345 1
10 1 12.12 5.94 2.423 3.5170 1.0898 -0.0830 -2.423 -3.0010 2.423 1
# Mid-lines are not shown here.
123 12 0.25 1.25 4.511 -3.2605 -0.2523 2.9495 -4.511 -0.0667 4.511 1
124 12 0.50 3.96 6.436 -2.4759 -0.6134 2.5903 -6.436 -0.2105 6.436 1
125 12 1.00 7.82 7.426 0.3939 -0.3457 0.8980 -7.426 -0.5636 7.426 1
126 12 2.00 9.72 6.998 2.7222 0.4872 -0.1429 -6.998 -1.2332 6.998 1
127 12 3.52 9.75 5.978 3.7717 1.7048 -0.2035 -5.978 -2.0051 5.978 1
128 12 5.07 8.57 5.080 3.4897 1.4530 -0.1741 -5.080 -2.5309 5.080 1
129 12 7.07 6.59 4.118 2.4721 -0.1349 -0.1411 -4.118 -2.9163 4.118 1
130 12 9.03 6.11 3.352 2.7582 0.8253 -0.1149 -3.352 -3.0637 3.352 1
131 12 12.05 4.57 2.441 2.1291 0.0660 -0.0836 -2.441 -3.0052 2.441 1
132 12 24.15 1.17 0.685 0.4850 -1.9998 -0.0235 -0.685 -1.7138 0.685 1위의 각 단계(step)를 이루는 함수의 내용은 CRAN의 source나, R console에서 함수 이름으로 살펴보면 알 수 있다.
더 자세하게 FO, FOCE, LAPL method 간에 목접 함수(objective function)의 차이는 다음의 두 논문을 참고한다.
Kim MG, Yim DS, Bae KS. R-based reproduction of the estimation process hidden behind NONMEM® Part 1: first-order approximation method. . 2015;23(1):01-07.
Bae KS, Yim DS. R-based reproduction of the estimation process hidden behind NONMEM® Part 2: First-order conditional estimation. . 2016;24(4):161-168.