tests/testthat/testPredict.R

In tdmore-dev/tdmore: In silico evaluation of precision dosing

library(nlmixr)
library(tdmore)
library(testthat)
library(RxODE)
library(dplyr)

expect_not_equal <- function(a, b) {
  expect_false(isTRUE(all.equal(a,b)))
}

m1 <- getModel("example2")
regimen <- data.frame(
  TIME=seq(0, by=24, length.out=10),
  AMT=500 # 500ug standard dose
)

#The below was generated by:
#predict(m1, newdata=c(0, 2, 4), regimen=regimen) %>% write.table(file="", row.names=F)
expected_output <- '
TIME Ka CL V1 V2 Q K12 K21 CONC center periph
0 3 700 60000 6000 250 0.004166666666666667 0.041666666666666664 8.333333333333334 500 0
2 3 700 60000 6000 250 0.004166666666666667 0.041666666666666664 8.076334766847998 484.5800860108799 3.9347334384664725
4 3 700 60000 6000 250 0.004166666666666667 0.041666666666666664 7.832422651165478 469.94535906992866 7.43479770669032
' %>% read.table(text=., header=T, row.names=c())

test_that("predict.tdmore() single profiles", {
  #empty newdata should still give all output columns
  expect_equal( tdmore:::predict.tdmore(m1, newdata=c()), subset(expected_output, FALSE))
  #no treatment regimen
  expect_equal( tdmore:::predict.tdmore(m1, newdata=c(0, 2)),
                expected_output %>% filter(TIME %in% c(0,2)) %>% mutate(CONC=0, center=0, periph=0) )
  #with treatment regimen
  expect_equal( tdmore:::predict.tdmore(m1, newdata=c(0, 2), regimen),
                expected_output %>% filter(TIME %in% c(0,2)) )
  #everything
  expect_equal( tdmore:::predict.tdmore(m1, newdata=c(0, 2, 4), regimen),
                expected_output )

  # only the TIME column
  expect_equal( tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=c(0, 2, 4)), regimen),
                expected_output %>% select(TIME) )
  # only the TIME and CONC column
  expect_equal( tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=c(0, 2, 4), CONC=NA), regimen),
                expected_output %>% select(TIME, CONC) )

  # If values in CONC, they are overwritten
  expect_equal( tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=c(0, 2, 4), CONC=runif(3)), regimen),
                expected_output %>% select(TIME, CONC) )

  # Also possible to ask for model outputs without any residual error model attached
  expect_equal( tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=c(0, 2, 4), periph=runif(3)), regimen),
                expected_output %>% select(TIME, periph) )

  tdmore:::predict.tdmore(m1, parameters=c(ECL=1),
          newdata=data.frame(TIME=c(0, 2, 4), periph=runif(3)),
          regimen)
})

#Generated by:
#predict(m2, newdata=c(0, 2, 4), regimen=regimen, covariates=c(WT=70)) %>% write.table(file="", row.names=F)
expected_output <- '
TIME KA CL V1 V2 Q K12 K21 CONC depot center periph
0 3.7 3.7 61 61 10 0.16393442622950818 0.16393442622950818 0 500 0 0
2 3.7 3.7 61 61 10 0.16393442622950818 0.16393442622950818 5.799123733991486 0.3056246732199422 353.74654777348064 102.0773668872091
4 3.7 3.7 61 61 10 0.16393442622950818 0.16393442622950818 4.292293518711172 1.8681218837224262e-4 261.82990464138146 157.54439316432845
' %>% read.table(text=., header=T)

m2 <- getModel("default")
test_that("predict.tdmore() single profiles with covariates model", {
  #missing covariates should give error
  expect_error( tdmore:::predict.tdmore(m2, newdata=c()), ".*covariate.*WT.*missing")

  #wrong format
  expect_error( tdmore:::predict.tdmore(m2, newdata=c(), covariates=list(WT=70)), ".*ovariates.*in wrong format")

  # covariate provided
  expect_equal( tdmore:::predict.tdmore(m2, newdata=c(), covariates=c(WT=70)), subset(expected_output, FALSE))

  #everything
  expect_equal( tdmore:::predict.tdmore(m2, newdata=c(0, 2, 4), regimen=regimen, covariates=c(WT=70)),
                expected_output )

  # time-varying covariate
  expect_equal( tdmore:::predict.tdmore(m2, newdata=c(0, 2), regimen=regimen,
                        covariates=data.frame(
                          TIME=c(0, 2.5),
                          WT=c(70, 700) #ridiculously high, but we should see difference
                          ),
                        covsInterpolation="locf"),
                expected_output %>% filter(TIME %in% c(0, 2)) )
  expect_not_equal( tdmore:::predict.tdmore(m2, newdata=c(0, 2), regimen=regimen,
                        covariates=data.frame(
                          TIME=c(0, 2.5),
                          WT=c(70, 700) #ridiculously high, but we should see difference
                        ),
                        covsInterpolation="linear"),
                expected_output %>% filter(TIME %in% c(0, 2)) )

  expect_not_equal( tdmore:::predict.tdmore(m2, newdata=c(0, 2), regimen=regimen,
                            covariates=data.frame(
                              TIME=c(0, 2.5),
                              WT=c(70, 700) #ridiculously high, but we should see difference
                            )),
                    expected_output %>% filter(TIME %in% c(0, 2, 4)) )
})


## TODO: The SE should be removed from predict.tdmore, see ticket #64
## Instead, these tests should be moved to testModelFrame
test_that("predict.tdmore() se", {
  #ask for SE
  prediction <- tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=0:2, CONC=NA), regimen, se=T, level=0.95)
  expect_setequal(names(prediction) , c("TIME", "CONC", "CONC.upper", "CONC.lower"))

  #SE ignores output without res_var
  prediction <- tdmore:::predict.tdmore(m1, newdata=data.frame(TIME=0:2, CONC=NA, periph=NA), regimen, se=T, level=0.95)
  expect_setequal(names(prediction) , c("TIME", "CONC", "periph", "CONC.upper", "CONC.lower"))
})

tdmore:::estimate(m2, regimen=regimen) #This does not give an error, because the model_predict routine is not even called
expect_error( tdmore:::estimate(m2, regimen=regimen, observed=data.frame(TIME=15, CONC=4)) ) #covariates not specified


pred <- tdmore:::estimate(m2, regimen=regimen, covariates=c(WT=70))
test_that("predict.tdmorefit() inherits parameters from the fit, even with no observed values", {
  expect_equal(
    ## TODO: should we create all output columns, or only TIME/CONC?
    predict(pred),
    expected_output %>% select(TIME,CONC) %>% filter(FALSE)
  )

  expect_equal(
    predict(pred, newdata=c(0, 2, 4)),
    expected_output
  )

  predict(pred, parameters=c(ECL=1), newdata=0:10)
})

ipred <- tdmore:::estimate(m2, regimen=regimen, covariates=c(WT=70), observed=data.frame(TIME=15, CONC=4))
# Generated by:
# predict(ipred, newdata=c(0, 2, 4, 15)) %>% write.table(file="", row.names=F)
expected_output <- '
"TIME" "KA" "CL" "V1" "V2" "Q" "K12" "K21" "CONC" "depot" "center" "periph"
0 3.7 2.71347592286047 51.0433635789244 51.0433635789244 10 0.195911854134334 0.195911854134334 0 500 0 0
2 3.7 2.71347592286047 51.0433635789244 51.0433635789244 10 0.195911854134334 0.195911854134334 6.76000923601633 0.305624668318032 345.053609230868 116.756194606358
4 3.7 2.71347592286047 51.0433635789244 51.0433635789244 10 0.195911854134334 0.195911854134334 5.03237148628253 0.000186812192765632 256.869167438531 173.817014761646
15 3.7 2.71347592286047 51.0433635789244 51.0433635789244 10 0.195911854134334 0.195911854134334 2.95423376266689 1.00853629280556e-12 150.79402804494 171.452154838672
' %>% read.table(text=., header=T)
test_that("predict.tdmorefit() with observed values", {
  expect_equal(
    predict(ipred),
    expected_output %>% filter(TIME==15) %>% select(TIME, CONC)
  )
  expect_equal(
    predict(ipred, newdata=15),
    expected_output %>% filter(TIME==15)
  )

  prediction <- predict(ipred, newdata=15, se.fit=T)
  vars <- expected_output %>% select(-TIME) %>% colnames(.)
  expect_setequal(
    colnames(prediction),
    c("TIME", paste0(vars, rep(c("", ".median", ".upper", ".lower"), each=length(vars))))
  )

  prediction <- predict(ipred, newdata=data.frame(TIME=15, CONC=NA), se.fit=T)
  expect_setequal(
    colnames(prediction),
    c("TIME", "CONC", "CONC.median", "CONC.lower", "CONC.upper")
  )

  # With no output values, should behave as expected
  expect_equal(
    predict(ipred, newdata=data.frame(TIME=15), se.fit=T),
    data.frame(TIME=15)
  )
})

## TODO: Predict.tdmorefit with se.fit
## TODO: predict.tdmorefit with se.fit and level=NA