- Setup
- Exploratory data analysis
- The model
- Assemble a model estimation data set
- Estimate parameters: Approach I
- Estimate parameters: Approach II
library(mrgsolve)
library(tidyverse)
theme_set(theme_bw())We have plasma concentrations fairly early after the dose and running out to 48 hours:
conc <- readRDS("two_endpoints_plasma.RDS")
conc. # A tibble: 9 x 3
. ID time CP
. <dbl> <dbl> <dbl>
. 1 1 0.25 1.73
. 2 1 0.5 2.25
. 3 1 1 3.59
. 4 1 2 4.62
. 5 1 6 2.67
. 6 1 10 2.88
. 7 1 16 1.80
. 8 1 24 1.62
. 9 1 48 0.379
ggplot(conc, aes(time,CP)) +
geom_point() + geom_line() +
ylab("Plasma concentration") + ylim(0,5)
urine <- readRDS("two_endpoints_urine.RDS")
urine. # A tibble: 6 x 3
. ID time UR
. <dbl> <dbl> <dbl>
. 1 1 2 1.23
. 2 1 10 7.25
. 3 1 20 9.90
. 4 1 30 12.7
. 5 1 40 15.4
. 6 1 72 16.9
ggplot(urine, aes(time,UR)) +
geom_point() + geom_line() +
ylab("Cumulative amount in urine")
The data above were simulated with this model, including the parameter values you see below.
mod <- mread("two_endpoints.cpp")[ SET ] end = 72, delta = 0.25, add = 0.05
[ PARAM ] CLnr = 0.97, V = 22.3, KA = 1.9, CLr = 0.2, dvtype = 0,
sigma1 = 1, sigma2 = 1
[ CMT ] GUT CENT URINE
[ ODE ]
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - CLnr*(CENT/V) - CLr*(CENT/V);
dxdt_URINE = CLr*(CENT/V);
[ SIGMA ] 0 0
[ TABLE ]
capture CP = (CENT/V)*exp(EPS(1));
capture UR = URINE*exp(EPS(2));
capture DV = dvtype ==2 ? UR : CP;
capture sigma = dvtype==2? sigma2 : sigma1;Check out the [ TABLE ] block. We check the value of dvtype and
return concentration if dvtype==1 and urine amount if dvtype==2.
Note that there isn’t really a loss of generality here. You can still
simulate from the model as you wish. In order for DV to be meaningful,
you’ll have to use a data set with dvtype defined. But you could still
use this model to simulate any compartment / output at any time … just
ignore DV and access CP or UR.
First, add a DV column to both data sets and set a flag for the type
of each DV
conc <- mutate(conc, DV = CP, dvtype = 1)
urine <- mutate(urine, DV = UR, dvtype = 2)Then make a single data set sorted by time of observation
obs <- bind_rows(conc,urine) %>% arrange(time)No longer need CP or UR
obs <- select(obs,ID,time,dvtype,DV)And we’ll add some columns that to make mrgsolve happy; this is the observations part of the data set
obs <- mutate(obs, evid = 0, cmt = 0, amt = 0)
obs. # A tibble: 15 x 7
. ID time dvtype DV evid cmt amt
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0.25 1 1.73 0 0 0
. 2 1 0.5 1 2.25 0 0 0
. 3 1 1 1 3.59 0 0 0
. 4 1 2 1 4.62 0 0 0
. 5 1 2 2 1.23 0 0 0
. 6 1 6 1 2.67 0 0 0
. 7 1 10 1 2.88 0 0 0
. 8 1 10 2 7.25 0 0 0
. 9 1 16 1 1.80 0 0 0
. 10 1 20 2 9.90 0 0 0
. 11 1 24 1 1.62 0 0 0
. 12 1 30 2 12.7 0 0 0
. 13 1 40 2 15.4 0 0 0
. 14 1 48 1 0.379 0 0 0
. 15 1 72 2 16.9 0 0 0
Just like we did in the Indometh example, pull a single row and set up
the dosing
dose <- obs %>% slice(1)
dose. # A tibble: 1 x 7
. ID time dvtype DV evid cmt amt
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0.25 1 1.73 0 0 0
Now, set up this chunk to be doses (not observations)
dose <- mutate(dose, time = 0, amt = 100, evid = 1, cmt = 1)We make this dvtype 0 and set DV to NA (it is dosing record so we
don’t want this to contribute information to parameter estimates)
dose <- mutate(dose, dvtype = 0, DV = NA)Then bind it all together and arrange
data <- bind_rows(obs,dose) %>% arrange(time)
data. # A tibble: 16 x 7
. ID time dvtype DV evid cmt amt
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0 0 NA 1 1 100
. 2 1 0.25 1 1.73 0 0 0
. 3 1 0.5 1 2.25 0 0 0
. 4 1 1 1 3.59 0 0 0
. 5 1 2 1 4.62 0 0 0
. 6 1 2 2 1.23 0 0 0
. 7 1 6 1 2.67 0 0 0
. 8 1 10 1 2.88 0 0 0
. 9 1 10 2 7.25 0 0 0
. 10 1 16 1 1.80 0 0 0
. 11 1 20 2 9.90 0 0 0
. 12 1 24 1 1.62 0 0 0
. 13 1 30 2 12.7 0 0 0
. 14 1 40 2 15.4 0 0 0
. 15 1 48 1 0.379 0 0 0
. 16 1 72 2 16.9 0 0 0
ofv1 <- function(p, data, yobs, dvcol = "DV", pred=FALSE) {
p <- lapply(p, exp)
names(p) <- names(theta)
mod <- param(mod,p)
out <- mrgsim_d(mod, data, output="df")
if(pred) return(as_tibble(out))
y_hat <- out[[dvcol]]
#sum((y_hat-yobs)^2,na.rm=TRUE)
llike <- dnorm(yobs, y_hat, out[["sigma"]], log = TRUE)
return(-1*sum(llike, na.rm=TRUE))
}yobs <- data[["DV"]]
theta <- log(c(CLnr = 2, V=10, KA = 1, CLr = 1, sigma1=1, sigma2=1))
fit <- nloptr::newuoa(theta, ofv1, data=data, yobs=yobs)The way we set this up, the true parameters are in the model
as.numeric(param(mod)). CLnr V KA CLr dvtype sigma1 sigma2
. 0.97 22.30 1.90 0.20 0.00 1.00 1.00
exp(fit$par). [1] 0.9238533 22.9345229 1.8367364 0.1941290 0.3579444 0.5502006
And make the plots for concentration and urine cumulative amount; we use
the dose data set (above) to get the smooth line
library(nlme)
prd <- ofv1(fit$par,data = dose, pred = TRUE)
prd <- ofv1(fit$par, data = data, pred=TRUE)
data <- mutate(data, PRED = prd[["DV"]], RES=DV-PRED)
h <- fdHess(fit$par, ofv1, data = data, yobs=yobs)
h$Hessian %>% solve %>% diag %>% sqrt. [1] 0.02121075 0.05928664 0.14850905 0.03558959 0.23648334 0.29009604
ggplot() +
geom_point(data=conc, aes(time,CP), size = 4) +
geom_line(data=prd, aes(time,CP),col="firebrick", lwd=1) +
scale_y_log10() + ylab("Plasma concentration")
ggplot() +
geom_point(data= urine, aes(time,UR), size = 4) +
geom_line(data=prd, aes(time,UR),col="cornflowerblue", lwd=1) +
ylab("Cumulative amout in urine")
I still think approach I is simpler and easier to pull off without loss of flexibility for simulation.
But if you didn’t want to do the DV approach, here is another way.
First, find which rows are type 1 (concentration) and which are type 2 (urine):
iconc <- which(data$dvtype==1)
iur <- which(data$dvtype==2)So we can use this iur vector to grab just the urine amounts
iur. [1] 6 9 11 13 14 16
data[iur,]. # A tibble: 6 x 9
. ID time dvtype DV evid cmt amt PRED RES
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 2 2 1.23 0 0 0 1.20 0.0384
. 2 1 10 2 7.25 0 0 0 6.41 0.843
. 3 1 20 2 9.90 0 0 0 10.6 -0.734
. 4 1 30 2 12.7 0 0 0 13.2 -0.527
. 5 1 40 2 15.4 0 0 0 14.8 0.528
. 6 1 72 2 16.9 0 0 0 16.8 0.0897
Now, create a vector of all the concentrations followed by all the urine
amounts; the final product is yobs
yplas <- data$DV[iconc]
yur <- data$DV[iur]
yobs <- c(yplas,yur)The OFV function will now have to modified as well
ofv2 <- function(p,data,pred=FALSE) {
p <- lapply(p, exp)
names(p) <- names(theta)
mod <- param(mod,p)
out <- mrgsim_d(mod,data)
if(pred) return(as_tibble(out))
## The action is here:
yplas_hat <- out$CP[iconc]
yur_hat <- out$UR[iur]
y_hat <- c(yplas_hat,yur_hat)
##---------------------------
sum((y_hat-yobs)^2)
}The function is mostly the same, but now we’ll slice the simulation output to first grab the concentrations and then grab the urine concentrations. NOTE: this works because the simulation output has the exact same design / setup as the input data. So we slice both the input data and the output by the same indices so we can be sure the observed and simulated values match up.
Ok, now try it out:
theta <- log(c(CLnr = 2, V =12, KA = 1, CLr = 1))
fit <- minqa::newuoa(par=theta, fn=ofv2, data=data)as.numeric(param(mod)). CLnr V KA CLr dvtype sigma1 sigma2
. 0.97 22.30 1.90 0.20 0.00 1.00 1.00
exp(fit$par). [1] 0.9165441 23.1189119 1.8859076 0.1937382
And plot
prd <- ofv2(fit$par,data = dose, pred = TRUE)
ggplot() +
geom_point(data= data[iconc,], aes(time,DV), size = 4) +
geom_line(data=prd, aes(time,CP),col="firebrick", lwd=1) +
scale_y_log10() + ylab("Plasma concentration")
ggplot() +
geom_point(data= data[iur,], aes(time,DV), size = 4) +
geom_line(data=prd, aes(time,UR),col="cornflowerblue", lwd=1) +
ylab("Cumulative amout in urine")