Question: How do I interpret simulation results in multi-endpoint model?

[billdenney]

I'm running a multi-endpoint model. When simulating from that model, I do not know how to extract the simulation results with residual variable uncertainty. Specifically, there is a sim column in the output, but it is not clear which endpoint that represents.

Unless I'm missing something, I think that the ask here is create the output columns differently. Using the multiple endpoint model example from the vignette, I think that it would be preferable to remove the sim and ipred columns and add new columns named something like sim_cp, sim_effect, ipred_cp and ipred_effect. In that example, the right-hand-size of the column name would come from the left-hand-side of the ~ defining the error model.

library(nlmixr)

pk.turnover.emax3 <- 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
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    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
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err) | pca
  })
}

fit.TOS <- nlmixr(pk.turnover.emax3, warfarin, "saem", control=list(print=0),
                  table=list(cwres=TRUE, npde=TRUE))
#> i parameter labels from comments will be replaced by 'label()'
#> RxODE 1.0.7.1 using 4 threads (see ?getRxThreads)
#> Calculating covariance matrix
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> Compiling NPDE model...done
#> done

d_sim <- nlmixrSim(fit.TOS)
#> Compiling model...done

head(as.data.frame(d_sim))
#>   id time     ipred       sim         depot       gut    center    effect
#> 1  1  0.5  1.271356  1.397402  1.960791e+01 70.887933  9.459017 102.91348
#> 2  1  1.0  3.236853  3.358431  3.844697e+00 71.817791 24.082512 100.98672
#> 3  1  2.0  6.447958  8.203813  1.478168e-01 50.703242 47.973460  96.77507
#> 4  1  3.0  8.529106  9.845301  5.683104e-03 33.953178 63.457404  92.59923
#> 5  1  6.0 10.978957 11.101657  3.229628e-07 10.102235 81.684547  81.04046
#> 6  1  9.0 11.097592 10.063121 -8.917957e-10  3.005189 82.567202  71.04630

^{Created on 2021-04-17 by the reprex package (v2.0.0)}

[mattfidler]

In theory you can add addCov=TRUE to include the compartment information.

There is a minor bug in the simulation script that doesn't convert CMT to the factor it should be; Once that is fixed you have:

library(nlmixr)

pk.turnover.emax3 <- 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
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    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
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err) | pca
  })
}

fit.TOS <- nlmixr(pk.turnover.emax3, warfarin, "saem", control=list(print=0),
                  table=list(cwres=TRUE, npde=TRUE))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> RxODE 1.0.8 using 4 threads (see ?getRxThreads)
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> Compiling NPDE model...done
#> done

d_sim <- nlmixrSim(fit.TOS, addCov=TRUE)
#> Compiling model...done
head(as.data.frame(d_sim))
#>   id time     ipred       sim        depot       gut   center    effect CMT
#> 1  1  0.5  1.977727  2.049699 6.118439e+00 78.868113 14.94109 102.62801  cp
#> 2  1  1.0  4.362342  4.429278 3.743526e-01 66.308903 32.95609 100.62802  cp
#> 3  1  2.0  7.675170  8.543195 1.401397e-03 40.546346 57.98343  96.39803  cp
#> 4  1  3.0  9.573255 10.188081 5.246148e-06 24.642287 72.32286  92.23970  cp
#> 5  1  6.0 11.334601 11.387380 7.090154e-10  5.531205 85.62926  80.76960  cp
#> 6  1  9.0 11.097316 10.666202 1.298825e-11  1.241534 83.83665  70.86947  cp

^{Created on 2021-04-17 by the reprex package (v2.0.0)}