---
title: "Warfarin PKPD"
author: "Nick Holford, Tomoo Funaki" 
output:
  html_document:
    toc: true
    toc_float: true
    number_sections: true
    self_contained: false
    theme: united
    highlight: tango
---

```{r include, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  out.width = "100%"
  )

```

# Libraries

The first step is to load the libraries we need - in this case, `nlmixr` and `xpose.nlmixr`. (We need `RxODE` and `xpose` as well, but those are loaded automatically by nlmixr.)

```{r libraries,echo=T,message=T,warnings=T,error=T,cache=F}

#Parameters
issim=3 # flag for data source: original data (0, simulated data lag time (1), simulated data transit (2) and (3)
ispds=F # flag to use SAEM (T) or FOCEI (F) for EBEs of PK parameters
isplot=0 # flag to show nlmixr diagnostic and residual plots (1) or not (0)

## Getting started
#Advice from nlmixr installation-notes.rtf (nlmixr_1.1.1-3)
#Make sure there is only one library path
.libPaths()
.libPaths(.libPaths()[regexpr("nlmixr",.libPaths())!=-1])
srclib=.libPaths() #srclib="M:\\Apps\\R\\R-4.0.5\\library"

## requires patchwork and xpose.nlmixr packages for markdown. 

library(ggplot2,lib=srclib)
library(nlmixr,lib=srclib)
library(RxODE,lib=srclib)
library(xpose,lib=srclib)
library(xpose.nlmixr,lib=srclib)
library(patchwork,lib=srclib)

#Show versions
packageVersion("vpc")
packageVersion("RxODE")
nlmixrver=packageVersion("nlmixr")
nlmixrver

```

# Local functions

```{r functions,echo=T,message=F,warnings=T,error=T,cache=T}

myvpc=function(nlmixr.fit,ylab) {
   nlmixr::vpc(nlmixr.fit,n=500, show=list(obs_dv=T,sim_median=T,pi=T), bins = bintimes,
       ylab = ylab, xlab = "Time (hours)", stratify="CMT", scales = "free_y")
}

# Goodness of fit plots using xpose
mygof=function(nlmixr.fit) {
   xp.fit <- xpose_data_nlmixr(nlmixr.fit)
   dvpred  <- dv_vs_pred(xp.fit) + labs(title = NULL, subtitle = NULL, caption = NULL)
   dvipred <- dv_vs_ipred(xp.fit) + labs(title = NULL, subtitle = NULL, caption = NULL)
   cwresidv <- res_vs_idv(xp.fit) + labs(title = NULL, subtitle = NULL, caption = NULL)   # default is CWRES
   cwrespred <- res_vs_pred(xp.fit) + labs(title = NULL, subtitle = NULL, caption = NULL) # default is CWRES  
   return (dvpred + dvipred + cwresidv + cwrespred + plot_layout(nrow=2) + plot_annotation(title=NULL))
}

```

# Data

Original warfarin PKPD data (O'Reilly et al. 1963, 1968)
#O'Reilly RA, Aggeler PM, Leong LS. Studies of the coumarin anticoagulant drugs: The pharmacodynamics of warfarin in man. J Clin Invest. 1963;42(10):1542-51.
#O'Reilly RA, Aggeler PM. Studies on coumarin anticoagulant drugs. Initiation of warfarin therapy without a loading dose. Circulation. 1968;38:169-77.

PKPD and turnover model for warfarin. (Holford 1986)
#Holford NH. Clinical pharmacokinetics and pharmacodynamics of warfarin. Understanding the dose-effect relationship. Clin Pharmacokinet. 1986;11(6):483-504.

NONMEM PKPD model used to simulate data (Funaki et al. 2018)
#Funaki T, Holford N, Fujita S. Population PKPD analysis using nlmixr and NONMEM.  Population Approach Group Japan2018.


```{r original-data,echo=T,message=F,warnings=T,error=T,cache=T}

# Read data
datadir="U:\\Research\\Pharmacometrics\\Users\\Tomoo Funaki\\warfarin\\data"
bintimes=c(seq(0,10,1),seq(12,144,12)) # bins for VPC
if (issim==0) {
   datafile="warfarin_dat.csv"
   rdata_ext=""
}
if (issim==1) {
   datafile="ka1_to_emax1_simln_items.csv"
   rdata_ext="_sim"
}
if (issim==2) {
   datafile="ka1_to_emax1_simtr_10.csv"
   rdata_ext="_simtr"
}

if (issim==3) {
   datafile="ka1_to_emax1_simtr2_10.csv"
   rdata_ext="_simtr2"
}

if (ispds) {
   rdata_ext=paste(rdata_ext,"_PDS",sep="")
} else {
   rdata_ext=paste(rdata_ext,"_PDF",sep="")
}

data <- read.csv(file = paste(datadir,datafile,sep="/"))
data$dv <- as.numeric(as.character(data$dv))
names(data) <- toupper(names(data))
data2 <- data
data2$DV <- ifelse(data2$TIME==0 & data2$DVID==1, NA, data2$DV)
```
## Plot data
Plot the PK and INR data
```{r spaghetti PKPD,echo=T,message=F,warnings=T,error=T,cache=T}

ggplot(data2, aes(x=TIME, y=DV)) +
  geom_line(aes(group=ID), col="red") +
  scale_x_continuous("Time (h)") +
  scale_y_continuous("Concentration | Effect") +
  facet_wrap(~DVID, scales = "free_y") +
  labs(title="Warfarin single-dose", subtitle="Concentration | Effect vs. time by individual")

```

# PK data

Create PK object from input data

```{r PK-data,echo=T,message=F,warnings=T,error=T,cache=T}
# Selecting only PK data
data.pk <- data[data$DVID==1, ]
data.pk <- data.pk[!is.na(data.pk$DV)|data.pk$AMT!=0, ]
data.pk$CMT <- ifelse(data.pk$AMT!=0, 1, 3)
save(data.pk,file=paste(datadir,"/","datapk",rdata_ext,".RData",sep=""))
#load(paste(datadir,"/",datapk",rdata_ext,".RData",sep=""))


```

# PK model
## pk.ka1t PK 

Let's start by fitting the PK model using a 1-compartment model with lag time, first-order absorption and elimination.

```{r pk.ka1t,echo=T,message=F,warnings=T,error=T,cache=T}

pk.ka1t <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    cp = center / v
    cp ~ prop(eps.pkprop) + add(eps.pkadd)
  })
}

```

### pk.fit.ka1t SAEM

```{r pk.fit.ka1tS,echo=T,message=F,warnings=T,error=T,cache=T}
pk.fit.ka1tS <- nlmixr(pk.ka1t, data.pk,  est="saem", control=list(print=0),table=list(cwres=T))
pk.fit.ka1tS
if (isplot) plot(pk.fit.ka1tS)
#bootstrapFit(pk.fit.ka1tS, nboot = 100, restart = TRUE) # for illustration but this takes a very long time

```

### pk.fit.ka1t VPC and GOF

The PK VPC using the original data does not describe the peak very well. Hard to see a problem with standard diagnostics.

```{r pk.fit.ka1t-vpc, message=F, warning=FALSE, error=F, echo=T, cache=T}

vpc=myvpc(pk.fit.ka1tS,ylab="PK Concentration (mg/L) saem")

gof=mygof(pk.fit.ka1tS)

vpc/gof

```

### pk.fit.ka1t FOCEI

```{r pk.fit.ka1tF,echo=T,message=F,warnings=T,error=T,cache=T}
pk.fit.ka1tF <- nlmixr(pk.ka1t, data.pk,  est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
pk.fit.ka1tF
if (isplot) plot(pk.fit.ka1tF)
```

### pk.fit.ka1t VPC & GOF

The PK VPC using the original data does not describe the peak very well. Hard to see a problem with standard diagnostics.

```{r pk.fit.ka1tF-vpc, message=F, warning=FALSE, error=F, echo=T, cache=T}

vpc=myvpc(pk.fit.ka1tF,ylab="PK Concentration (mg/L) focei")

gof=mygof(pk.fit.ka1tF)

vpc/gof

```

# PD data

## PD data object with PK EBEs

The PK model empirical Bayes estimates are combined with the effect data to allow a sequential PKPD model.

```{r pd-data,message=F,echo=T,warning=F,error=F,cache=T}

# Make PD data with SAEM (ispds=1) or FOCEI (ispds=0) EBEs for PK parameters
if (ispds) {
   fit.ipara <- pk.fit.ka1tS[pk.fit.ka1tS$TIME==24,]
} else {
   fit.ipara <- pk.fit.ka1tF[pk.fit.ka1tF$TIME==24,]
}

fit.ipara <- fit.ipara[c("ID", "ktr", "ka", "cl", "v")]
names(fit.ipara) <- c("ID", "IKTR", "IKA", "ICL", "IV")

data.pd <- data
data.pd <- merge(data.pd, fit.ipara, by="ID", all=T)

data.pd$CMT <- ifelse(data.pd$AMT!=0, 1,
                ifelse(data.pd$DVID==2 & data.pd$AMT==0, 4, 3))

#Create MDV=1 for dose record or when DV is recorded as zero
data.pd$MDV <- ifelse(data.pd$AMT>0 | data.pd$DV==0, 1, 0)
#select dosing records and PCA effect records
data.pd <- data.pd[data.pd$AMT!=0 | data.pd$DVID==2, ]

save(data.pd,file=paste(datadir,"/","datapd",rdata_ext,".RData",sep=""))
#load(file=paste(datadir,"/","datapd",rdata_ext,".RData",sep=""))

```

## PKPD data object
The PKPD data object includes both concentration and effect observations.

```{r pkpd-data,message=F,echo=T,warning=F,error=F,cache=T}

data.pkpd <- read.csv(file = paste(datadir,datafile,sep="/"))
names(data.pkpd) = toupper(names(data.pkpd))
data.pkpd$CMT <- data.pkpd$DVID + 2 # DVID 0->2 (center);DVID 2->3 (center); DVID 3->4 (effect)
data.pkpd$CMT[data.pkpd$AMT != 0] <- 1 # CMT 2->1
data.pkpd$EVID[data.pkpd$AMT != 0] <- 101

save(data.pkpd,file=paste(datadir,"/","datapkpd",rdata_ext,".RData",sep=""))
#load(file=paste(datadir,"/","datapkpd",rdata_ext,".RData",sep=""))


```

# Sequential PD Models

## pd.cp.logit Immediate PD logit

Sequential immediate effect model - logit for Emax
```{r pd.cp.logit, echo=T, message=F, error=F, warning=F,cache=T }
# define the models
pd.cp.logit <- function() {
  ini({
    poplogit <- 7.5
    tc50 <- log(0.5)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    effect=e0*(1-emax*cp/(c50+cp))

    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.CPS SAEM

```{r pd.fit.CPS,echo=F, message=F, error=F, warning=F,cache=T }
################ Immediate PD model SAEM ###########
pd.fit.CPS <- nlmixr(pd.cp.logit, data.pd,  est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pd.fit.CPS)
pd.fit.CPS
```

### pd.fit.CPS VPC & GOF

The VPC and GOF plots shows the immediate effect model describes the data poorly.

```{r pd.fit.CPS-vpc,echo=T, message=F, error=F, warning=F,cache=T }

vpc=myvpc(pd.fit.CPS,ylab="CPS Response (PCA) saem")

gof=mygof(pd.fit.CPS)

vpc/gof

```

## pd.ce.logit Effect CMT PD logit

Sequential effect compartment model - logit for Emax
```{r pd.ce.logit, echo=T, message=F, error=F, warning=F,cache=T }
pd.ce.logit <- function() {
  ini({
    poplogit <- 7.5
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    d/dt(ce) = kout*(cp - ce)

    effect=e0*(1-emax*ce/(c50+ce))

    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.CES SAEM

```{r pd.fit.CES, message=FALSE, warning=FALSE, error=F,cache=T}

pd.fit.CES <- nlmixr(pd.ce.logit, data.pd,  est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pd.fit.CES)
pd.fit.CES

```

### pd.fit.CES VPC & GOF

The VPC and GOF show the effect compartment model does not describe the data well.

```{r pd.fit.CES-vpc, echo=T,message=F, warning=FALSE,error=F,cache=T}

vpc=myvpc(pd.fit.CES,ylab="CES Response (PCA) saem")

gof=mygof(pd.fit.CES)

vpc/gof

```

## pd.turnover.logit Turnover PD logit

Sequential turnover model - logit for Emax

```{r pd.turnover.logit, echo=T, message=F, error=F, warning=F,cache=T }
pd.turnover.logit <- function() { # The logit transformation is used to constrain Emax in SAEM
  ini({
    poplogit <- 7.5
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    PD=1-emax*cp/(c50+cp)

    effect(0) = e0
    kin = e0*kout

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    d/dt(effect) = kin*PD -kout*effect


    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.TOS SAEM

Sequential turnover model logit

```{r pd.fit.TOS,echo=T, message=F, error=F, warning=F,cache=T }

################ Run nlmixr with PD and PKPD models ###########################
pd.fit.TOS <- nlmixr(pd.turnover.logit, data.pd,  est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pd.fit.TOS)
pd.fit.TOS

```

### pd.fit.TOS VPC & GOF

The VPC and GOF show the turnover model describes the PD data well.

```{r pd.fit.TOS-vpc, echo=T,message=F, error=F, warning=FALSE,cache=T}

vpc=myvpc(pd.fit.TOS,ylab="pd TOS Response (PCA) saem")

gof=mygof(pd.fit.TOS)

vpc/gof

```

## pd.cp.emax Immediate PD Emax

Sequential immediate effect model Emax<1

```{r pd.cp.emax, echo=T, message=F, error=F, warning=F,cache=T }
pd.cp.emax <- function() { 
  ini({
    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    poplogit = log(popemax/(1-popemax))
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    effect=e0*(1-emax*cp/(c50+cp))

    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.CPF FOCEI

```{r pd.fit.CPF,echo=F, message=F, error=F, warning=F,cache=T }
################ Immediate PD model FOCEI ###########
pd.fit.CPF <- nlmixr(pd.cp.emax, data.pd,  est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pd.fit.CPF)
pd.fit.CPF
```

### pd.fit.CPF VPC & GOF

The VPC and GOF shows the immediate effect model does not describe the PD data well.

```{r pd.fit.CPF-vpc, echo=T,message=F, error=F, warning=FALSE,cache=T}

vpc=myvpc(pd.fit.CPF,ylab="pd CPF Response (PCA) focei")

gof=mygof(pd.fit.CPF)

vpc/gof

```

## pd.ce.emax effect CMT Emax

Sequential effect compartment model Emax<1

```{r pd.ce.emax, echo=T, message=F, error=F, warning=F,cache=T }
pd.ce.emax <- function() { #clearer emax parameter constraint method which works with FOCEI
  ini({
    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    poplogit = log(popemax/(1-popemax))
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    d/dt(ce) = kout*(cp - ce)

    effect=e0*(1-emax*ce/(c50+ce))

    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.CEF FOCEI

```{r pd.fit.CEF, echo=T, message=FALSE, warning=FALSE, cache=T}
################ Effect compartment PD model FOCEI ###########
#   pd.fit.CEF <- nlmixr(pd.ce.emax, data.pd,  est="focei", control=foceiControl(print=100000)) # Default opt method fails
pd.fit.CEF <- nlmixr(pd.ce.emax, data.pd,  est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pd.fit.CEF)
pd.fit.CEF
```

### pd.fit.CEF VPC & GOF

```{r pd.fit.CEF-vpc, echo=T,message=F, error=F, warning=FALSE,cache=T}

vpc=myvpc(pd.fit.CEF,ylab="pd CEF Response (PCA) focei")

gof=mygof(pd.fit.CEF)

vpc/gof

```

## pd.turnover.emax Turnover PD Emax

Sequential turnover model Emax<1

```{r pd.turnover.emax, echo=T, message=F, error=F, warning=F,cache=T }
pd.turnover.emax <- function() {
  ini({
    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    poplogit = log(popemax/(1-popemax))
    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    # PK parameters from input dataset
    ktr = IKTR
    ka = IKA
    cl = ICL
    v = IV

    cp = center/v

    PD=1-emax*cp/(c50+cp)

    effect(0) = e0
    kin = e0*kout

    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center

    d/dt(effect) = kin*PD -kout*effect

    effect ~ add(eps.pdadd)
  })
}
```

### pd.fit.TOF FOCEI

```{r pd.fit.TOF, message=FALSE, warning=FALSE, echo=T, cache=T}

################ Turnover (true) PD model FOCEI ###########
pd.fit.TOF <- nlmixr(pd.turnover.emax, data.pd,  est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pd.fit.TOF)
pd.fit.TOF
```

### pd.fit.TOF VPC & GOF

Quitting from lines 704-711 (warf_simtr_PKPD_vpc_ui_saem_foce_B-18-v3_PDF.Rmd)
Error in foceiFit.data.frame0(data = .dat, inits = .inits, PKpars = .uif$theta.pars,  :
  Could not fit data.
Calls: <Anonymous> ... getVarCov.nlmixrFitCoreSilent -> .setCov -> foceiFit.data.frame0

```{r pd.fit.TOF-vpc, echo=T,message=F, error=F, warning=FALSE,cache=T}

#vpc=myvpc(pd.fit.TOF,ylab="pd TOF Response (PCA) focei")

#gof=mygof(pd.fit.TOF)

#vpc/gof

```

# Joint PKPD models

We can jointly fit the PK and PD data with a single model.

## pkpd.cp.logit Imm PKPD logit

Joint turnover model Emax<1

```{r pkpd.cp.logit, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.cp.logit <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    poplogit <- 7.5
    tc50 <- log(0.5)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    effect=e0*(1-emax*cp/(c50+cp))

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | ce
  })
}

```

### pkpd.fit.CPS SAEM

```{r pkpd.fit.CPS,echo=F, message=F, error=F, warning=F,cache=T }
################ Immediate PKPD logit model SAEM ###########
nlmixr(pkpd.cp.logit) # check multiple endpoint cmt and dvid assignments
pkpd.fit.CPS <- nlmixr(pkpd.cp.logit, data.pkpd, est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pkpd.fit.CPS)
pkpd.fit.CPS
```

### pkpd.fit.CPS VPC & GOF

```{r pkpd.fit.CPS-vpc, echo=T,message=F, error=F, warning=FALSE,cache=T}

vpc=myvpc(pkpd.fit.CPS,ylab="pkpd CPS Response (PCA) saem")

gof=mygof(pkpd.fit.CPS)

vpc/gof

```

## pkpd.ce.logit Effect CMT PKPD logit

Joint effect CMT logit

```{r pkpd.ce.logit, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.ce.logit <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    poplogit <- 7.5
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    d/dt(ce) = kout*(cp - ce)

    effect=e0*(1-emax*ce/(c50+ce))

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | ce
  })
}

```

### pkpd.fit.CES SAEM

Effect compartment model using the combined PKPD data.

```{r pkpd.fit.CES,echo=F,message=F,error=F,warning=F,cache=T}

################ Effect compartment PKPD logit model SAEM ###########
nlmixr(pkpd.ce.logit) # Show initial estimates and model
pkpd.fit.CES <- nlmixr(pkpd.ce.logit, data.pkpd, est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pkpd.fit.CES)
pkpd.fit.CES

```

### pkpd.fit.CES VPC & GOF

VPC medians do not agree. Note the time trends in the GOF plots (e.g. CPRED vs CWRES).

```{r pkpd.fit.CES-vpc,echo=T,message=F,error=F,warning=F,cache=T}

vpc=myvpc(pkpd.fit.CES,ylab="pkpd CES Response (PCA) saem")

gof=mygof(pkpd.fit.CES)

vpc/gof

```

## pkpd.turnover.logit TO PKPD logit

Joint turnover model logit

```{r pkpd.turnover.logit, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.turnover.logit <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    poplogit <- 7.5
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    logit=poplogit+eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    effect(0) = e0
    kin = e0*kout
    PD=1-emax*cp/(c50+cp)
    d/dt(effect) = kin*PD -kout*effect

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | effect
  })
}

```

### pkpd.fit.TOS SAEM

```{r pkpd.fit.TOS, message=FALSE, warning=FALSE, echo=T, cache=T}
################ Turnover (true) PKPD logit model SAEM ###########
nlmixr(pkpd.turnover.logit) # Show initial estimates and model
pkpd.fit.TOS <- nlmixr(pkpd.turnover.logit, data.pkpd, est="saem", control=list(print=0),table=list(cwres=T))
if (isplot) plot(pkpd.fit.TOS)
pkpd.fit.TOS

```

### pkpd.fit.TOS VPC

VPCs look acceptable and  no problem visible in the standard GOF plots.

```{r pkpd.fit.TOS-vpc, echo=T, message=FALSE, warning=FALSE,error=F, cache=T}

vpc=myvpc(pkpd.fit.TOS,ylab="pkpd TOS Response (PCA) saem")

gof=mygof(pkpd.fit.TOS)

vpc/gof

```

## pkpd.cp.emax Imm Effect PKPD Emax

Joint immediate effect model Emax<1

```{r pkpd.cp.emax, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.cp.emax <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    poplogit = log(popemax/(1-popemax))
    logit=poplogit + eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    effect=e0*(1-emax*cp/(c50+cp))

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | effect
  })
}

```

### pkpd.fit.CPF FOCEI

```{r pkpd.fit.CPF,echo=T, message=F, error=T, warning=F,cache=T }
################ Immediate PKPD emax model FOCEI ###########
pkpd.fit.CPF <- nlmixr(pkpd.cp.emax, data.pkpd, est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pkpd.fit.CPF)
pkpd.fit.CPF

```

### pkpd.fit.CPF VPC & GOF

Compiling VPC model...done
done (14.89 sec)
Quitting from lines 1036-1043 (warf_simtr2_PKPD_vpc_ui_saem_foce_B-18-v3_PDF.Rmd)
Error in as.character.factor(df[[col]]) : malformed factor
Calls: <Anonymous> ... unique -> rbind_dfs -> as.character -> as.character.factor
In addition: Warning message:
Undefined 'dvid' integer values in data: 1, 2

```{r pkpd.fit.CPF-vpc, echo=T, message=FALSE, warning=FALSE,error=F, cache=T}

#vpc=myvpc(pkpd.fit.CPF,ylab="pkpd CPF Response (PCA) saem")

#gof=mygof(pkpd.fit.CPF)

#vpc/gof

```

## pkpd.ce.emax Effect CMT PKPD Emax

Joint effect compartment model Emax<1

```{r pkpd.ce.emax, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.ce.emax <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    poplogit = log(popemax/(1-popemax))
    logit=poplogit + eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    d/dt(ce) = kout*(cp - ce)

    effect=e0*(1-emax*ce/(c50+ce))

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | ce
  })
}

```

### pkpd.fit.CEF FOCEI

```{r pkpd.fit.CEF, echo=F, message=F, error=F, warning=F,cache=T }

################ Effect compartment PKPD emax model FOCEI ###########
pkpd.fit.CEF <- nlmixr(pkpd.ce.emax, data.pkpd, est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pkpd.fit.CEF)
pkpd.fit.CEF

```

### pkpd.fit.CEF VPC & GOF

```{r pkpd.fit.CEF-vpc, echo=T, message=FALSE, warning=FALSE,error=F, cache=T}

vpc=myvpc(pkpd.fit.CEF,ylab="pkpd CEF Response (PCA) saem")

gof=mygof(pkpd.fit.CEF)

vpc/gof

```

## pkpd.turnover.emax TO PKPD Emax

Joint turnover model Emax<1

```{r pkpd.turnover.emax, echo=T, message=F, error=F, warning=F,cache=T }
pkpd.turnover.emax <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)

    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    eps.pkprop <- 0.1
    eps.pkadd <- 0.1

    popemax <- c(0,0.9,0.999)
    tc50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)

    eta.emax ~ .5
    eta.c50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5

    eps.pdadd <- 100
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    poplogit = log(popemax/(1-popemax))
    logit=poplogit + eta.emax
    emax = 1/(1+exp(-logit))
    c50 =  exp(tc50 + eta.c50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)

    cp = center / v
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl * cp

    effect(0) = e0
    kin = e0*kout
    PD=1-emax*cp/(c50+cp)
    d/dt(effect) = kin*PD -kout*effect

    cp ~ prop(eps.pkprop) + add(eps.pkadd) | center
    effect ~ add(eps.pdadd) | effect
  })
}

```

### pkpd.fit.TOF FOCEI

```{r pkpd.fit.TOF, message=FALSE, warning=FALSE, echo=T, cache=T}
################ Turnover (true) PKPD emax model FOCEI ###########
pkpd.fit.TOF <- nlmixr(pkpd.turnover.emax, data.pkpd, est="focei", control=list(print=0,outerOpt="bobyqa",table=list(cwres=T)))
if (isplot) plot(pkpd.fit.TOF)
pkpd.fit.TOF

```

### pkpd.fit.TOF VPC

```{r pkpd.fit.TOF-vpc, message=FALSE, warning=FALSE, echo=T, cache=T}

vpc=myvpc(pkpd.fit.TOF,ylab="pkpd TOF Response (PCA) saem")

gof=mygof(pkpd.fit.TOF)

vpc/gof

```

# Rmarkdown 

Output created: warf_simtr_PKPD_vpc_ui_saem_foce_B-18-v3_PDF.html

```{r rmd-warning, message=FALSE, warning=FALSE, echo=F, cache=F}

#%Rmethod%\R -e "rmarkdown::find_pandoc(dir='M:/Apps/R/Rstudio_1.4.1106/bin/pandoc'); rmarkdown::render('%script%.Rmd',output_file='%script%.html')"
warnings()

```