tdmore-dev / tdmore

In silico evaluation of precision dosing
https://tdmore-dev.github.io/tdmore/dev
GNU Affero General Public License v3.0
3 stars 1 forks source link

predict.tdmore guesses the output column to predict, predict.tdmorefit does not #24

Closed rfaelens closed 6 years ago

rfaelens commented 6 years ago

When using a tdmorefit object, you need to specify the output column:

data <- predict(pred, se.fit=TRUE, level=0.95, newdata = data.frame(TIME=seq(0, 24, by = 0.1), CONC=NA))

Can we not guess it as well?

data <- predict(pred, se.fit=TRUE, level=0.95, newdata = seq(0, 24, by = 0.1))

nluyckx commented 6 years ago

If I am not wrong, predict.tdmore does not guess anything but will output all the lhs and state variables. Should we do the same? Or should we guess the output columns from the 'observed' dataframe (like we do in the plot function)?

rfaelens commented 6 years ago

Test scenario in attachment

rfaelens commented 6 years ago
library(nlmixr)
#Model taken from literature: Soulele, K., et al. "Population pharmacokinetics of fluticasone propionate/salmeterol using two different dry powder inhalers." European Journal of Pharmaceutical Sciences 80 (2015): 33-42.
modelCode <- function(){
  ini({
    TVKa <- 3.87
    TVCL <- 659 #L/h
    TVV1 <- 56900 #L
    TVV2 <- 5550 #L
    TVQ <- 259 #L/h

    EKa ~ 0.04507129 #0.2123**2
    ECL ~ 0.1535856 #0.3919**2
    EV1 ~ 0.09223369 #0.3037**2
    EV2 ~ 0.208301 #0.4564**2
    EQ ~ 0.1015697# 0.3187**2

    EPS_ADD <- 1.91 #
    EPS_PROP <- 0.117
  })
  model({
    Ka <- TVKa * exp(EKa)
    CL <- TVCL * exp(ECL)
    V1 <- TVV1 * exp(EV1)
    V2 <- TVV2 * exp(EV2)
    Q <- TVQ * exp(EQ)
    K12 <- Q/V1
    K21 <- Q/V2

    d/dt(center) = - CL/V1 * center - K12*center + K21 * periph
    d/dt(periph) = K12*center - K21 * periph

    CONC = center / V1 * 1000
    CONC ~ prop(EPS_PROP) + add(EPS_ADD)
  })
}
nlmixrModel <- nlmixrUI(modelCode)

library(tdmore)
m1 <- tdmore(nlmixrModel)

regimen <- data.frame(
  TIME=seq(0, by=24, length.out=30),
  AMT=500 # 500ug standard dose
)

## Let tdmore figure it out
## TDMore cannot figure it out, because we have no observations)
predict(m1, regimen=regimen, newdata=NULL) 
#indeed, newdata is a required argument

## You get all the variables
predict(m1, regimen=regimen, newdata=seq(0, 30))

## This is an error!
predict(m1, regimen=regimen, newdata=data.frame(TIME=seq(0, 30)))

pred <- estimate(m1, regimen=regimen)

## TDMore guesses what you want
## This behaviour is very badly defined....
## Actually, we should simply show the proediction for the 'observed' data
## If observed data is empty, we should give back an empty data.frame
predict(pred, regimen=regimen, newdata=NULL)
predict(pred, regimen=regimen, newdata=NULL, se.fit=TRUE)

## Great, we get all the variables
predict(pred, regimen=regimen, newdata=seq(0, 30), se.fit=FALSE)

## Oh, we only get TIME...
## We should get all of the parameters
predict(pred, regimen=regimen, 
        newdata=seq(0, 30),
        se.fit=TRUE)

## Ok, also an error, same as for predict.tdmore
predict(pred, regimen=regimen, 
        newdata=data.frame(TIME=seq(0, 30)),
        se.fit=TRUE)

## Now we get TIME and CONC
predict(pred, regimen=regimen, 
        newdata=data.frame(TIME=seq(0, 30), CONC=NA),
        se.fit=TRUE)