Friberg myelosuppression model

library(nlmixr2)

wbc <- function() {
  ini({
    ## Note that the UI can take expressions
    ## Also note that these initial estimates should be provided on the log-scale
    log_CIRC0 <- log(7.21)
    log_MTT <- log(124)
    log_SLOPU <- log(28.9)
    log_GAMMA <- log(0.239)
    ## Initial estimates should be high for SAEM ETAs
    eta.CIRC0  ~ .1
    eta.MTT  ~ .03
    eta.SLOPU ~ .2
    ##  Also true for additive error (also ignored in SAEM)
    prop.err <- 10
  })
  model({
    CIRC0 = exp(log_CIRC0 + eta.CIRC0)
    MTT = exp(log_MTT + eta.MTT)
    SLOPU = exp(log_SLOPU + eta.SLOPU)
    GAMMA = exp(log_GAMMA)

    # PK parameters from input dataset
    CL = CLI
    V1 = V1I
    V2 = V2I
    Q = 204

    CONC = centr/V1

    # PD parameters
    NN = 3
    KTR = (NN + 1)/MTT
    EDRUG = 1 - SLOPU * CONC
    FDBK = (CIRC0 / circ)^GAMMA

    CIRC = circ

    prol(0) = CIRC0
    tr1(0) = CIRC0
    tr2(0) = CIRC0
    tr3(0) = CIRC0
    circ(0) = CIRC0

    d/dt(centr) = periph * Q/V2 - centr * (CL/V1 + Q/V1)
    d/dt(periph) = centr * Q/V1 - periph * Q/V2
    d/dt(prol) = KTR * prol * EDRUG * FDBK - KTR * prol
    d/dt(tr1) = KTR * prol - KTR * tr1
    d/dt(tr2) = KTR * tr1 - KTR * tr2
    d/dt(tr3) = KTR * tr2 - KTR * tr3
    d/dt(circ) = KTR * tr3 - KTR * circ

    CIRC ~ prop(prop.err)
  })
}

Fit model using saem

d3 <- read.csv("Simulated_WBC_pacl_ddmore_samePK_nlmixr.csv", na.strings = ".")

fit.S <- nlmixr(wbc, d3, est="saem", list(print=0), table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00

library(xpose.nlmixr2)

xpdb <- xpose_data_nlmixr(fit.S)

plot(fit.S)

print(dv_vs_pred(xpdb) +
      ylab("Observed Neutrophil Count (10^9/L)") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(dv_vs_ipred(xpdb) +
      ylab("Observed Neutrophil Count (10^9/L)") +
      xlab("Individual Predicted Neutrophil Count (10^9/L)"))


print(res_vs_pred(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(res_vs_idv(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Time (h)"))


print(prm_vs_iteration(xpdb))


print(absval_res_vs_idv(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Time (h)"))


print(absval_res_vs_pred(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(ind_plots(xpdb, nrow=3, ncol=4) +
      ylab("Predicted and Observed Neutrophil Count (10^9/L)") +
      xlab("Time (h)"))


print(res_distrib(xpdb) +
      ylab("Density") +
      xlab("Conditional Weighted Residuals"))


vpcPlot(fit.S, n=500, n_bins = 10, show=list(obs_dv=TRUE),
        ylab = "Neutrophil Count (10^9/L)", xlab = "Time (h)")
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04


vpcPlot(fit.S, n=500, bins = c(0,170,300,350,500,600,900,3000,4580),
        show=list(obs_dv=TRUE), ylab = "Neutrophil Count (10^9/L)", xlab = "Time (h)")
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04

Fit model using FOCEi

fit.F <- nlmixr(wbc, d3, est="focei", list(print=0), table=list(cwres=TRUE, npde=TRUE))
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:27 
#> done
#> [====|====|====|====|====|====|====|====|====|====] 0:00:02

FOCEi goodness of fit plots

xpdb <- xpose_data_nlmixr(fit.F)
plot(fit.F)

print(dv_vs_pred(xpdb) +
      ylab("Observed Neutrophil Count (10^9/L)") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(dv_vs_ipred(xpdb) +
      ylab("Observed Neutrophil Count (10^9/L)") +
      xlab("Individual Predicted Neutrophil Count (10^9/L)"))


print(res_vs_pred(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(res_vs_idv(xpdb) +
      ylab("Conditional Weighted Residuals") +
      xlab("Time (h)"))


print(absval_res_vs_idv(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Time (h)"))


print(absval_res_vs_pred(xpdb, res = 'IWRES') +
      ylab("Individual Weighted Residuals") +
      xlab("Population Predicted Neutrophil Count (10^9/L)"))


print(ind_plots(xpdb, nrow=3, ncol=4) +
      ylab("Predicted and Observed Neutrophil Count (10^9/L)") +
      xlab("Time (h)"))


print(res_distrib(xpdb) +
      ylab("Density") +
      xlab("Conditional Weighted Residuals"))


# 10 bins is slightly better than auto bin
vpcPlot(fit.F, n=500, n_bins = 10, show=list(obs_dv=TRUE),
        ylab = "Neutrophil Count (10^9/L)", xlab = "Time (h)")
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04


# specify bins
vpcPlot(fit.F, n=500, bins = c(0, 170, 300, 350, 500, 600, 900, 3000, 4580),
        show=list(obs_dv=TRUE),
        ylab = "Neutrophil Count (10^9/L)", xlab = "Time (h)")
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04

Yuan Xiong

Fit model using saem

Fit model using FOCEi

FOCEi goodness of fit plots