05.txt

---
output:
  pdf_document: default
  html_document: default
always_allow_html: true
---

# 구획분석의 자료해석 {#ca-analysis}

\Large\hfill
한성필^[약동학 데이터를 분석하여 PK 파라미터를 찾는 업무를 하지 않는 독자들은 이 장을 읽지 않고 넘어가도 좋습니다.] 

\normalsize

```{r include=FALSE}
library(knitr)
library(tidyverse)
opts_chunk$set(echo=TRUE)
```

## 서론

약동학 구획분석을 한 문장으로 표현하자면 시간 농도(time-concentration data) 자료에 대한 약동학 파라미터 추정값(pharmacokinetic parameter estimates)을 구하는 과정이다. \index{약동학}
다시 말하자면, 인체의 여러 장기(organ) 또는 계(system)를 몇개의 구획으로 단순화하고 우리의 관심이 되는 특성을 파라미터화하여 각 구획들 간의 연관성을 파악하는 것이다. 
약동학 파라미터의 추정을 위해서 비선형 회귀분석(nonlinear regression analysis)이 사용 되어 왔고 현재 이 방식은 약동학 분석에서 가장 필수적인 방법론으로 자리 잡은 상태이다. \index{약동학}
추정에는 약동학 수식이 필요한데, 비선형 회귀분석 시 모델식(model formulation)으로 사용되며, 이 수식 안에 우리가 추정하고자 하는 약동학 파라미터가 모두 포함되어 있다. \index{약동학}

구획간의 약물 이동은 1차식 또는 선형 약동학(linear pharmacokinetics)을 따른다고 가정하게 된다. \index{약동학}
선형 약동학이란 간단히 말해서 용량이 두배가 되면 농도도 두 배가 되는 것을 말하며, D(용량) → C(농도)로 나타낼 경우 2D → 2C 의 관계가 성립하는 것을 말한다.\index{약동학}
1구획 모델에서는 구획들 간의 약물 분포가 균일하다는 가정 하에 인체를 단일구획으로 기술한다. 
1구획 모델에서는 약물이 작용부위로 신속히 이동한다고 가정하며 혈장을 자료 분석의 기준 구획으로 사용한다. 
1구획 모델의 로그 혈장농도값은 시간에 대해 하나의 기울기로 감소하는 양상을 나타낸다.
2구획 모델에서는 구획들 간의 약물 분포가 균일하지 않다고 보며 인체를 중심 구획 (central compartment) 과 말초구획 (peripheral compartment)로 나누어 기술한다. \index{말초구획}
중심구획은 관류(perfusion)가 빠른 기관(예: 간, 콩팥)을 나타내는데 이는 약물의 분포와 관계하고 α-phase로 나타낸다. \index{중심구획}
반면 말초구획은 관류가 느린 기관(예: 근육, 지방)을 나타내며 이는 약물의 소실과 관계하고 β-phase로 나타낸다. \index{말초구획}
2구획 모델의 혈장농도 로그값은 시간에 대해 기울기가 두 개인 이직선적으로 감소하는 양상을 나타내는데, 처음의 빠른 기울기는 분포기를 나타내고, 나중의 느린 기울기는 소실기를 나타내게 된다.

```{r fig5-1, fig.cap = "(ref:three-comp)", echo=FALSE, out.width = "75%"}
knitr::include_graphics("media-05/fig5-1.png")
```

(ref:three-comp) 3구획 약동학 모델에서 시간에 따른 농도 변화의 양상에 대한 모식도 [@han2014pharmacokinetic]

정맥주사 후 3구획 모델을 비유를 들어 설명하자면, 구획을 나타내는 원통이 그림 \@ref(fig:fig5-1)과 같이 세 개가 있다고 가정을 할 수 있다. 
약물을 정맥주사(i.v.bolus)하면 직후에 혈중 농도가 최대에 도달하였다가 시간이 지남에 따라 감소하게 되는데, 구획간의 분포 및 제거에 의해 약물의 농도가 감소를 하게 된다. 
여기서 V~1~은 중심구획, V~2~는 rapid peripheral compartment이고, V~3~은 slow peripheral compartment이다. \index{중심구획}
V~1~, V~2~, V~3~ 원통의 부피는 약물이 분포하는 용적이고, 높이는 농도에 해당한다. 
수평의 연결 파이프에 의해 약물이 분포 및 제거된다. 
정맥주사 후 초기에는 약물이 주로 V~1~에서 V~2~로 이동하면서 혈중농도가 급격하게 감소하고(rapid distribution phase), 그 이후에는 V~1~에서 V~3~로 약물이 이동하면서 동시에 V~2~에서는 거꾸로 V~1~으로 이동하며 전체적으로는 혈중농도가 감소(slow distribution phase)하게 된다.
그 후에는 V~2~, V~3~에서의 약물은 V~1~으로 이동해 나오면서 동시에 V~1~에서는 elimination이 일어나 농도가 아주 서서히 감소하게 된다. (terminal phase)

이렇게 약물의 이동에 대한 정량적인 표현을 위해서는 약동학 수식에 대해 익숙해 져야 할 필요가 있고, 수식을 잘 아는 것 만큼이나 자료를 해석하고 추정값을 실제로 구하는 작업을 충분히 경험할 필요가 있다. \index{약동학}

## 상용 소프트웨어를 이용한 구획분석 개론

```{r fig5-2, fig.cap = "다양한 상용 소프트웨어의 사용 예. A. GUI 기반 Phoenix NLME, B. GUI 기반 Lixoft Monolix, C. NONMEM의 제어구문", echo = FALSE}
knitr::include_graphics('media-05/fig5-2.png')
```

약동학 자료의 분석은 복잡하고 시간이 많이 걸리는 일이다. 모든 생물학 자료가 그러하듯 높은 변동성과 많은 노이즈가 PK 자료에 포함되어 있다. [@gabrielsson]\index{약동학}
따라서 비선형 회귀분석(Nonlinear regression analysis)이 필요한 구획분석을 위해서는 컴퓨터 소프트웨어의 사용이 필수적이다. 
이 때 Certara Phoenix NLME^[https://www.certara.com/software/phoenix-nlme/], Monolix^[https://lixoft.com/products/monolix/]와 같은 여러가지 소프트웨어들이 사용되고 있지만, 가장 널리쓰이는 도구는 ICON사의 NONMEM^[https://www.iconplc.com/innovation/nonmem/]이라고 하는 소프트웨어이다. (그림 \@ref(fig:fig5-2))
이것은 FORTRAN 코드로 프로그램 되어 있으며 2021년 현재 버전 7.5.0까지 지속적으로 업데이트 되고 있는 소프트웨어로, 구획분석 및 계량약리학 분야에서 가장 표준적인 도구로 쓰이고 있다. \index{계량약리학}
이에 대한 각종 자료를 이곳에서 다루기는 지면의 제약과 이론의 어려움 때문에 별도의 참고서적을 활용하는 것이 바람직할 것이다. [@nonmem;@gabrielsson;@basic]

상용 소프트웨어를 통한 구획분석 방법은 GUI 기반이고 많은 사용자가 확보되어 있다는 장점이 있으나,  상당한 비용을 매년 지불하여야 하며, 접근이 다소 쉽지 않다는 단점을 가지고 있다. 

## 구획분석에 활용할 수 있는 R package 소개

이를 극복하기 위해 통계 및 수학적 계산에 적합한 R 언어 기반의 오픈소스 소프트웨어가 활발히 개발되고 있다. 
즉, 비구획분석을 위한 R 패키지를 [3장](#nca-analysis)에서 소개한 것과 마찬가지로 구획분석용 R 패키지에 대한 다양한 시도가 일어나고 있는 것이다.
본 서에서는 무료 R 패키지 "wnl"를 사용한 자료해석과 예제에 대해 다루고자 한다. [@R-wnl]

NONMEM에서는 비선형 회귀분석을 위해서 최대 가능도 추정법(maximum likelihood estimation)을 사용한다. \index{가능도}
최대 가능도 추정법에서 얻어진 sample 데이터에서 Y가 현재 측정된 값, E(Y)는 모델에 의해 각각의 관측 지점에서의 예상되는 값, 그리고 V(Y)를 그 모델의 분산으로 표현한다. \index{가능도}
이 때 모델에 주어진 관측값에 대한 가능도(likelihood)는 다음과 같이 표현될 수 있다.\index{가능도}

\begin{equation}
L = \frac{1}{\sqrt{2\pi V(Y)}} \cdot e^{-\frac{1}{2V(Y)}\cdot(Y-E(Y))^2}
(\#eq:eq5-1)
\end{equation}

이때 이 일련의 관찰 값에 대한 가능도는 개개인 관찰값의 확률을 통해 나오게 된다. \index{가능도}
따라서 n개의 관찰값이 있으면 그 종합된 가능도는 n개의 개개인별 확률이 되고 subject에 대한 가능도는 다음과 같이 계산된다.\index{가능도}

\begin{equation}
L_i = \prod_{j=1}^{n_i}\frac{1}{\sqrt{2\pi V(Y_{ij})}} \cdot e^{-\frac{1}{2V(Y_{ij})}\cdot(Y_{ij}-E(Y_{ij}))^2}
(\#eq:eq5-2)
\end{equation}

이것은 -2 로그 가능도(-2 log likelihood)로 전환 될 수 있으며, 식 \@ref(eq:eq5-3)다음과 같이 표현된다.\index{가능도}

\begin{equation}
-2log(L_i) = n_i \cdot log(2\pi) + \sum_{j=1}^{n}(log(V(Y_{ij}))+ \frac{(Y_{ij} - E(Y_{ij}))^2}{V(Y_{ij})})
(\#eq:eq5-3)
\end{equation}

결국 NONMEM은 이 -2 로그 가능도의 값을 최소화(즉, 가능도의 값을 최대화) 해 주는 PK 파라미터 값들의 조합을 찾아준다. 
이 때, 오른쪽 수식에서 앞에있는 n~i~log(2π)는 상수이기 때문에 없애고, 실제 i번째 환자에서의 목적함수(objective function)는 다음과 같이 나타난다.\index{가능도}

\begin{equation}
OBJ_i = n_i \cdot logV(Y_i) + \sum_{j=1}^{n} \frac{(Y_{ij} - E(Y_{ij}))^2}{V(Y_{ij})}
(\#eq:eq5-4)
\end{equation}

실제 여기에서 derivation 을 정규 분포에만 적용한 것이 최대 가능도 측정법(maximum likelihood estimation method)^[최대 가능도 측정법(maximum likelihood estimation method)은 점 추정치를 찾는데까지는 오래 걸리나, 표준 오차(standard error of estimates)를 구하는게 유리하다. 최근에는 최소 제곱법(least square method) 보다 널리 쓰이고 있다.]이라면, 정규 분포라는 가정 없이도 적용한 것이 NONMEM이며 이 방식을 `확장된 최소 제곱법'(extended least square method)이라고 한다. \index{가능도}

이 책의 예제에서는 NONMEM의 확장된 최소 제곱법을 적용한 wnl 패키지의 `nlr()` 함수를 사용하였다.

각 모델에 맞는 적합한 모델에 따라 다음과 같은 script 형태로 입력해 주게 되며, 그에 대한 표준 오차(standard error of estimates), 공분산 행렬(covariance matrix of estimates), run test results on residuals, AIC, AICc, SBC를 내부 계산을 통하여 도출하게 된다.

이때 run test란 잔차(residual)들이 무작위적으로 분포하는지를 보기 위해서 각각의
잔차에서 최대 양수와 음수가 연속구간의 길이가 얼마나 되는지를 보는
방법이다. 이 때 계산된 p-value를 통해서 실제로 잔차가 무작위적으로 분포하는지, 편향되게 분포하는지를 고려해 볼 수 있다.

AIC, AICc, SBC는 Akaike information criterion, Akaike information
criterion corrected, Schwarz criterion을 의미하며 각각의 통계적 모델에서
구조가 적합한지를 확인할수 있도록 도와주는 지표이다.


## 1구획모델을 이용한 경구투여 후 약동학 자료 분석\index{약동학}

아래의 코드로 wnl 라이브러리를 불러오고 hands-on1.csv^[https://github.com/pipetcpt/pharmapk/blob/master/data.zip] 파일을 읽어올 수 있다. 

```{r}
library(wnl)
dPK02 = read.csv("data/hands-on1.csv", skip=1)
colnames(dPK02) = c("TIME", "DV")
plot(dPK02[,"TIME"], dPK02[,"DV"], xlim=c(0, 400), ylim=c(0, 2.5), 
     xlab="Time (min)", ylab="Concentration (ug/L)", pch=16)
```

먼저 NCA를 수행한다. T~max~는 40분, C~max~는 2.00 ug/L인 것을 알 수 있다.

```{r}
NonCompart::sNCA(dPK02[,"TIME"], dPK02[,"DV"], 
                 dose=100, doseUnit="ug", timeUnit="min")
```

아래에서 경구 복용 후 1차식^[1차식(first-order kinetics)에 의한 소실이란 체내에 남이 있는 특정시점의 약물량의 일정 분율만큼 약물이 단위시간당 소실되는 것을 말한다.]에 의해 소실되는 약물의 약동학 모델을 사용한다. \index{약동학}
이 경우에는 wnl 패키지 내부에 기본적 equation이 정의되어 있으므로 추가적인 equation을 입력하지 않고 각각의 파라미터의 값에 해당하는 THETA에 대해서만 입력하였다. 
필요한 THETA는 k(제거속도상수), k~a~(흡수속도상수), V(부피), t~lag~(흡수가 지연되는 시간)이다. 

### Compartmental analysis without Tlag

아래는 경구 투여 1구획, 지연 시간이 없는 모델이다. \index{지연 시간}
모델의 R 스크립트에서 THETA[1], THETA[2], THETA[3]는 각각 파라미터 중 k(제거속도상수), k~a~(흡수속도상수), V(부피)를 의미한다. 

```{r}
DOSE = 100

fPK02a = function(THETA) # Prediction function
{
  Ka = THETA[1]
  V  = THETA[2]
  K  = THETA[3]
  Cp = DOSE/V*Ka/(Ka - K)*(exp(-K*TIME) - exp(-Ka*TIME))
  return(Cp)
}
```

위 코드에서 `DOSE`는 Cp의 단위와 일치해야 하며, `fPK2a()`는 함수(function)로서, THETA를 입력 인자로 받고 Cp를 출력값으로 갖도록 작성된다. 
각각의 nonlinear regression을 위해서 THETA에 대한 초기값(initial estimate)을 아래와 같이 넣어준다. (`nlr()` 함수의 `IE()` 인자에 적용)

```{r}
TIME = dPK02[,"TIME"]
r1 = nlr(fPK02a, dPK02, pNames=c("ka", "V", "k"), IE=c(0.1, 30, 0.05))

r1$Est
```

### Compartmental analysis with T~lag~ 

경구 투여 1구획분석, 지연 시간이 있는 모델이다. k~a~, V, k에 T~lag~가 추가되어 농도를 나타낼 수 있다.\index{지연 시간}

\small 

```{r}
fPK02b = function(THETA) # Prediction function
{
  Ka   = THETA[1]
  V    = THETA[2]
  K    = THETA[3]
  tlag = THETA[4]

  Cp  = DOSE/V*Ka/(Ka - K)*(exp(-K*(TIME - tlag)) - exp(-Ka*(TIME - tlag))) 
  return(Cp)
}
```

```{r}
TIME = dPK02[,"TIME"]
r2 = nlr(fPK02b, dPK02, pNames=c("ka", "V", "k", "tlag"), 
         IE=c(0.1, 30, 0.05, 20))

r2$Est
```

\normalsize

### 모델링 결과

위 두개의 모델을 비교할 때 지연 시간이 있는 모델의 적합이 더 좋은 것을 알 수 있다. (그림 \@ref(fig:onecomp-lag)) C~max~ 부분을 주의깊게 살펴야 한다. \index{지연 시간}

```{r onecomp-lag, fig.cap = "(ref:onecomp-lag)", fig.width = 6, fig.height = 3.5}
plot(dPK02[,"TIME"], dPK02[,"DV"], xlim=c(0, 400), ylim=c(0, 2.5), 
     xlab="Time (min)", ylab="Concentration (ug/L)", pch=16)
TIME = 0:400
lines(TIME, fPK02a(r1$Est["PE", 1:3]), lty=2)
lines(TIME, fPK02b(r2$Est["PE", 1:4]))
```

(ref:onecomp-lag) wnl 패키지를 통해 얻어진 경구 투여 1구획분석 PK 파라미터로 400분까지 시뮬레이션한 결과. 모델에서 지연 시간(lag time)의 유무에 따라 C~max~ 부분의 차이가 두드러지는 것을 알 수 있음 (실선: 지연 시간이 있는 모델, 점선: 지연 시간이 없는 모델, 점: 관찰값)

## 1구획모델을 이용한 정맥주사 후 약동학 자료 분석\index{약동학}

wnl 라이브러리를 불러오고 자료를 읽어온다. 
이전 예제에서는 1명의 경구투약 후 농도 자료였고, 이번 예제는 4명의 i.v.투약 후 농도 자료이다. (그림 \@ref(fig:ggonecomp-iv))
`hands-on2.csv`^[https://github.com/pipetcpt/pharmapk/blob/master/data.zip]를 읽어 들인 후 그림을 통해 자료의 탐색을 수행한다.

```{r ggonecomp-iv, fig.cap = "1구획모델을 이용한 정맥주사 후 농도-시간 곡선 (N=4)"}
library(wnl)
dPK01 <-  read.csv("data/hands-on2.csv", skip=1)
colnames(dPK01) <- c("TIME", "DV", "ID")

library(ggplot2)
ggplot(dPK01, aes(TIME, DV, group = ID, color = as.factor(ID))) +
  geom_line() + geom_point() + scale_y_log10() +
  labs(color = "Subject No.")
```

- 4명의 대상자 모두 대략적으로 시간에 따른 농도 감소가 단항 지수함수적(monoexponential)인 것을 관찰할 수 있다.	
- 대상자 1과 2를 비교하면 2번 대상자가	AUC가 더 작으며, 따라서 청소율이 더 클 것이라\index{청소율}
예상할 수 있으며, Y	 절편이 거의 같은 것으로 보아 분포용적이 유사할 것으로 보인다.	
- 대상자 3과	4의 경우 곡선이 교차하는 형태로 눈으로는 어느 쪽이 AUC가 클지 알기 어렵고, Y절편에
해당하는 농도가 높은 쪽이 분포용적이 더 작을 것이라 예상할 수 있다.

이러한 탐색을 수치적인 해석으로 확장하기 위해 4명 자료의 NCA 분석을 수행한다. [3장](#nca-analysis)에서 다루었던 `tblNCA()` 함수를 사용해 계산할 수 있다.

```{r}
NonCompart::tblNCA(dPK01, key="ID", colTime="TIME", 
                   colConc="DV", dose=10, adm="Bolus")
```

### Compartmental analysis

V, K만 있으면 단항 지수함수적 농도 감소를 보이는 i.v.dosing의 농도를 나타낼 수 있으므로 아래와 같이 간단한 함수를 만들 수 있다.

```{r}
IDs = unique(dPK01[,"ID"])
nID = length(IDs)
DOSE = 10000 # ug

fPK01 = function(THETA) # Prediction function
{
  V  = THETA[1]
  K  = THETA[2]
  Cp = DOSE/V*exp(-K*TIME)  
  return(Cp)
}
```

Subject가 여러 명일 경우 각 대상자의 값을 구하기 위해서 `for`  함수를 사용하게 되는데, 다소 복잡해보이지만 `nlr` 함수를 여러번 사용하며 각 대상자의 약동학 파라미터 추정을 수행하는 것이 핵심이다.\index{약동학}

```{r}
Result = vector()
for (i in 1:nID) {
  cID = IDs[i]
  Data = dPK01[dPK01$ID == cID,]
  TIME = dPK01[dPK01$ID == cID,"TIME"]
  Res = nlr(fPK01, Data, pNames=c("V", "k"), IE=c(20, 0.2),
            SecNames=c("CL", "AUC", "AUMC" , "Thalf", "MRT"), 
            SecForms=c(~V*k, ~DOSE/V/k, ~DOSE/V/k/k, ~log(2)/k, ~1/k))
  Result = rbind(Result, cbind(ID=cID, Res$Est))
} ; Result
```

### 모델링 결과

정맥주사 후 약동학 자료에 대해 1구획모델로 구획분석을 수행하였고, 이때 구해진 PK 파라미터를 사용하여 150분까지 시뮬레이션 한 결과를 그림으로 나타내었다 (그림 \@ref(fig:onecomp-iv)). \index{약동학}\index{시뮬레이션}

\small

```{r onecomp-iv, fig.cap = "wnl 패키지를 통해 얻어진 정맥 투여 1구획모델 PK 파라미터로 150분까지 시뮬레이션한 결과. (실선: 예측값 곡선, 도형: 관찰값)", fig.width = 6, fig.height = 3.5}
plot(0, 1, type="n", xlim=c(0, 160), ylim=c(10, 1000), 
     log="y", xlab="Time (min)", ylab="Concentration (ug/L)")
for (i in 1:nID) {
  cID = IDs[i]
  TIME = dPK01[dPK01$ID == cID,"TIME"]
  points(TIME, dPK01[dPK01$ID == cID,"DV"], pch=14+i)
  cTHETA = Result[Result[,"ID"]==cID & rownames(Result)=="PE", c("V", "k")]
  lines(TIME, fPK01(cTHETA))   
}
```

\normalsize

## 2구획모델 경구 흡수 이용한 약동학 자료 분석\index{약동학}

wnl 라이브러리를 불러오고 자료를 읽어온다. 
`hands-on3.csv`^[https://github.com/pipetcpt/pharmapk/blob/master/data.zip]를 읽어 들인 후 모델링을 수행한다.\index{모델링}

```{r}
library(wnl)
dPK14 <- read.csv("data/hands-on3.csv", skip=1)
colnames(dPK14) = c("TIME", "DV") ; dPK14
```

### Compartmental analysis without Tlag

아래의 코드는 경구 투여 2구획, 지연 시간이 없는 모델이다. \index{지연 시간}
2차 지수로 감소하는 약물에서의 약동학 모델이며, 이때의 농도 값(Cp)을 구하기 위한 공식을 사용하였으며, K~a~(`Ka`), V~c~(`Vc`), k~21~(`k21`), alpha(`a`), beta(`b`)로 농도를 나타낼 수 있다.\index{약동학}

```{r}
Dpo = 23158

## without lag
fPK14a = function(THETA)
{
  Vc  = THETA[1]
  Ka  = THETA[2]
  k21 = THETA[3]
  a   = THETA[4] # alpha
  b   = THETA[5] # beta
  
  T1 = e$DATA[,"TIME"]
  Co = Ka*Dpo/Vc*((k21-a)/(Ka-a)/(b-a)*exp(-a*T1) + 
                    (k21-b)/(Ka-b)/(a-b)*exp(-b*T1) + 
                    (k21-Ka)/(a-Ka)/(b-Ka)*exp(-Ka*T1)) 
  
  return(Co)
}

r1 <- nlr(fPK14a, dPK14, 
          pNames=c("Vc/F", "Ka", "k21", "alpha", "beta"), 
          IE=c(350, 11, 1, 0.1, 0.01))

r1$Est
```

### Compartmental analysis with T~lag~ 

경구 투여 2구획, 지연 시간이 있는 모델이다. Ka, V, k21, alpha, beta 에 Tlag가 추가되어 농도를 나타낼 수 있다.\index{지연 시간}

```{r}
## with lag
fPK14b = function(THETA)
{
  Vc  = THETA[1]
  Ka  = THETA[2]
  k21 = THETA[3]
  a   = THETA[4] # alpha
  b   = THETA[5] # beta
  TL  = THETA[6] # Tlag
  
  T1 = e$DATA[,"TIME"]
  Co = Ka*Dpo/Vc*((k21-a)/(Ka-a)/(b-a)*exp(-a*(T1-TL)) + 
                    (k21-b)/(Ka-b)/(a-b)*exp(-b*(T1-TL)) + 
                    (k21-Ka)/(a-Ka)/(b-Ka)*exp(-Ka*(T1-TL))) 
  Co[Co < 0] = 0 # remove negative concentrations before tlag
  
  return(Co)
}

r2 <- nlr(fPK14b, dPK14, 
          pNames=c("Vc/F", "Ka", "k21", "alpha", "beta", "Tlag"), 
          IE=c(150, 11, 0.12, 0.1, 0.01, 0.05))
r2$Est
```

### 모델링 결과

1구획모델에서와 같이 지연 시간이 있는 모델의 적합이 더 좋은 것을 알 수 있다. 특히 C~max~ 부분에서의 적합 차이를 주의해서 살펴 볼 필요가 있다.\index{지연 시간}

```{r twocomp-lag, fig.cap = "(ref:twocomp-lag)", fig.width = 6, fig.height = 3.5}
plot(dPK14[,"TIME"], dPK14[,"DV"], xlim=c(0, 25), ylim=c(0, 250), 
     xlab="Time (min)", ylab="Concentration (ug/L)", pch=16)
TIME = dPK14[,"TIME"]
lines(TIME, fPK14a(r1$Est["PE", 1:5]), lty=2)
lines(TIME, fPK14b(r2$Est["PE", 1:6]))
```

(ref:twocomp-lag) wnl 패키지를 통해 얻어진 경구 투여 2구획모델 PK 파라미터로 25시간 시뮬레이션한 결과. 모델에서 지연 시간(lag time)의 유무에 따라 C~max~ 부분의 차이가 두드러지는 것을 알 수 있다. (실선: 지연 시간이 있는 모델, 점선: 지연 시간이 없는 모델, 점: 관찰값)

## Theoph 약동학 자료의 구획분석 {#theoph-wnl}

R에 내장된 Theoph 약동학 자료의 구획분석을 아래의 코드를 사용하여 수행한다. Minimization이 된 후 1구획 약동학 모델에서 언급된 파라미터 이외에 다른 파라미터에 대해서도 값을 구하기 위해서 다음과 같이 CL(제거 속도), t~1/2~(반감기), MRT(평균잔류시간)를 추가한 후, 입력한 파라미터값들을 구하기 위한 수식을 뒤쪽 부분 `SecForms' 인자에 입력하게 된다. 이런 방법을 통하여 우리가 원하는 파라미터들의 점추정치와(point estimate), 그리고 표준 오차(standard error of estimates)를 모두 구할 수 있다.\index{약동학}\index{반감기}

본문에 결과를 싣기에 지면이 부족하여 [별첨 B](#Theoph_nlr)에 수록하였다.

```{r eval=FALSE}
library(wnl)
tData = Theoph
colnames(tData) = c("ID", "BWT", "DOSE", "TIME", "DV")

fPK = function(THETA)     # Prediction function
{
  DOSE = 320000           # in microgram
  TIME = e$DATA[,"TIME"]  # use data in e$DATA

  K    = THETA[1]
  Ka   = THETA[2]
  V    = THETA[3]

  Cp   = DOSE/V*Ka/(Ka - K)*(exp(-K*TIME) - exp(-Ka*TIME))
  return(Cp)
}

IDs = unique(tData[,"ID"])
nID = length(IDs)
for (i in 1:nID) {
  Data = tData[tData$ID == IDs[i],]
  Theoph_nlr = nlr(fPK, Data, 
                   pNames=c("k", "ka", "V"), IE=c(0.1, 3, 500), 
                   SecNames=c("CL", "Thalf", "MRT"), 
                   SecForms=c(~V*k, ~log(2)/k, ~1/k))
  print(paste("## ID =", i, "##"))
  print(Theoph_nlr)
}
```

## 맺음말

이상으로 구획모델링(compartment modeling)에 대한 이론과 R코드를 사용한 자료해석을 알아보았다. \index{모델링}
이는 단순하지만 여러 현상을 묘사할 수 있는 유용한 도구 중 하나이다. 
즉, 간단한 수식으로 구성된 모델을 통해 현상을 분석하고 여러 변수를 연관시킨 후, 시뮬레이션을 통해 여러 상황들을 예측할 수 있다. \index{시뮬레이션}
여기에서는 가장 단순한 형태의 모델을 사용하여 공개 소프트웨어를 사용한 구획분석을 수행하였으나, 실제 자료 분석에는 보다 전문적인 소프트웨어를 사용하게 된다. 
그러나 본 교재에서 다룬 기초적인 내용이 전문적인 구획분석 도구를 이해하고 사용하는데 필요한 핵심 내용을 포함하고 있으므로 여전히 유용할 것이다.