Skip to contents

rxode2 in a pipeline

2024-04-23

Source: vignettes/rxode2-pipeline.Rmd
rxode2-pipeline.Rmd

Setting up the rxode2 model for the pipeline

In this example we will show how to use rxode2 in a simple pipeline.

We can start with a model that can be used for the different simulation workflows that rxode2 can handle:

library(rxode2)
#> rxode2 2.1.2.9000 using 2 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

Ribba2012 <- function() {
  ini({
    k = 100
    tkde = 0.24
    eta.tkde = 0
    tkpq = 0.0295
    eta.kpq = 0
    tkqpp = 0.0031
    eta.kqpp = 0
    tlambdap = 0.121
    eta.lambdap = 0
    tgamma = 0.729
    eta.gamma = 0
    tdeltaqp = 0.00867
    eta.deltaqp = 0
    prop.sd <- 0
    tpt0 = 7.13
    eta.pt0 = 0
    tq0 = 41.2
    eta.q0 = 0
  })
  model({
    kde ~ tkde*exp(eta.tkde)
    kpq ~ tkpq * exp(eta.kpq)
    kqpp ~ tkqpp * exp(eta.kqpp)
    lambdap ~ tlambdap*exp(eta.lambdap)
    gamma ~ tgamma*exp(eta.gamma)
    deltaqp ~ tdeltaqp*exp(eta.deltaqp)
    d/dt(c) = -kde * c
    d/dt(pt) = lambdap * pt *(1-pstar/k) + kqpp*qp -
      kpq*pt - gamma*c*kde*pt
    d/dt(q) = kpq*pt -gamma*c*kde*q
    d/dt(qp) = gamma*c*kde*q - kqpp*qp - deltaqp*qp
    ## initial conditions
    pt0 ~ tpt0*exp(eta.pt0)
    q0 ~ tq0*exp(eta.q0)
    pt(0) = pt0
    q(0) = q0
    pstar <- (pt+q+qp)
    pstar ~  prop(prop.sd)
  })
}

This is a tumor growth model described in Ribba 2012. In this case, we compiled the model into an R object Ribba2012, though in an rxode2 simulation pipeline, you do not have to assign the compiled model to any object, though I think it makes sense.

Simulating one event table

Simulating a single event table is quite simple:

  • You pipe the rxode2 simulation object into an event table object by et().
  • When the events are completely specified, you simply solve the ODE system with rxSolve().
  • In this case you can pipe the output to plot() to conveniently view the results.
  • Note for the plot we are only selecting the selecting following:
    • pt (Proliferative Tissue),
    • q (quiescent tissue)
    • qp (DNA-Damaged quiescent tissue) and
    • pstar (total tumor tissue)
Ribba2012 %>% # Use rxode2
    et(time.units="months") %>% # Pipe to a new event table
    et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
    et(0, 250, by=0.5) %>% # Add some sampling times (not required)
    rxSolve() %>% # Solve the simulation
    plot(pt, q, qp, pstar) # Plot it, plotting the variables of interest
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

Simulating multiple subjects from a single event table

Simulating with between subject variability

The next sort of simulation that may be useful is simulating multiple patients with the same treatments. In this case, we will use the omega matrix specified by the paper:

## Add CVs from paper for individual simulation
## Uses exact formula:

lognCv = function(x){log((x/100)^2+1)}

library(lotri)
## Now create omega matrix
## I'm using lotri to quickly specify names/diagonals
omega <- lotri(eta.pt0 ~ lognCv(94),
               eta.q0 ~ lognCv(54),
               eta.lambdap ~ lognCv(72),
               eta.kqp ~ lognCv(76),
               eta.kqpp ~ lognCv(97),
               eta.deltaqp ~ lognCv(115),
               eta.tkde ~ lognCv(70))

omega
#>               eta.pt0    eta.q0 eta.lambdap   eta.kqp  eta.kqpp eta.deltaqp
#> eta.pt0     0.6331848 0.0000000   0.0000000 0.0000000 0.0000000   0.0000000
#> eta.q0      0.0000000 0.2558818   0.0000000 0.0000000 0.0000000   0.0000000
#> eta.lambdap 0.0000000 0.0000000   0.4176571 0.0000000 0.0000000   0.0000000
#> eta.kqp     0.0000000 0.0000000   0.0000000 0.4559047 0.0000000   0.0000000
#> eta.kqpp    0.0000000 0.0000000   0.0000000 0.0000000 0.6631518   0.0000000
#> eta.deltaqp 0.0000000 0.0000000   0.0000000 0.0000000 0.0000000   0.8426442
#> eta.tkde    0.0000000 0.0000000   0.0000000 0.0000000 0.0000000   0.0000000
#>              eta.tkde
#> eta.pt0     0.0000000
#> eta.q0      0.0000000
#> eta.lambdap 0.0000000
#> eta.kqp     0.0000000
#> eta.kqpp    0.0000000
#> eta.deltaqp 0.0000000
#> eta.tkde    0.3987761

With this information, it is easy to simulate 3 subjects from the model-based parameters:

set.seed(1089)
rxSetSeed(1089)
Ribba2012 %>% # Use rxode2
    et(time.units="months") %>% # Pipe to a new event table
    et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
    et(0, 250, by=0.5) %>% # Add some sampling times (not required)
    rxSolve(nSub=3, omega=omega) %>% # Solve the simulation
    plot(pt, q, qp, pstar) # Plot it, plotting the variables of interest
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> Warning: multi-subject simulation without without 'omega'

Note there are two different things that were added to this simulation: - nSub to specify how many subjects are in the model - omega to specify the between subject variability.

Simulation with unexplained variability

You can even add unexplained variability quite easily:

Ribba2012 %>% # Use rxode2
  ini(prop.sd=0.05) %>% # change variability
  et(time.units="months") %>% # Pipe to a new event table
  et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
  et(0, 250, by=0.5) %>% # Add some sampling times (not required)
  rxSolve(nSub=3, omega=omega) %>%
  plot(pt, q, qp, sim) # Plot it, plotting the variables of interest
#> ℹ change initial estimate of `prop.sd` to `0.05`
#> Warning: multi-subject simulation without without 'omega'


# note that sim is the simulated pstar since this is simulated from the
# model with a nlmixr2 endpoint

In this case we only added the sigma matrix to have unexplained variability on the pstar or total tumor tissue.

You can even simulate with uncertainty in the theta omega and sigma values if you wish.

Simulation with uncertainty in all the parameters (by matrices)

If we assume these parameters came from 95 subjects with 8 observations apiece, the degrees of freedom for the omega matrix would be 95, and the degrees of freedom of the sigma matrix would be 95*8=760 because 95 items informed the omega matrix, and 760 items informed the sigma matrix.

Ribba2012 %>% # Use rxode2
  ini(prop.sd = 0.05) %>%
  et(time.units="months") %>% # Pipe to a new event table
  et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
  et(0, 250, by=0.5) %>% # Add some sampling times (not required)
  rxSolve(nSub=3, nStud=3, omega=omega,
          dfSub=760, dfObs=95) %>% # Solve the simulation
  plot(pt, q, qp, sim) # Plot it, plotting the variables of interest
#> ℹ change initial estimate of `prop.sd` to `0.05`
#> Warning: multi-subject simulation without without 'omega'

Often in simulations we have a full covariance matrix for the fixed effect parameters. In this case, we do not have the matrix, but it could be specified by thetaMat.

While we do not have a full covariance matrix, we can have information about the diagonal elements of the covariance matrix from the model paper. These can be converted as follows:

rseVar <- function(est, rse){
  return(est*rse/100)^2
}

thetaMat <- lotri(tpt0 ~ rseVar(7.13,25),
                  tq0 ~ rseVar(41.2,7),
                  tlambdap ~ rseVar(0.121, 16),
                  tkqpp ~ rseVar(0.0031, 35),
                  tdeltaqp ~ rseVar(0.00867, 21),
                  tgamma ~ rseVar(0.729, 37),
                  tkde ~ rseVar(0.24, 33)
                  )

thetaMat
#>            tpt0   tq0 tlambdap    tkqpp  tdeltaqp  tgamma   tkde
#> tpt0     1.7825 0.000  0.00000 0.000000 0.0000000 0.00000 0.0000
#> tq0      0.0000 2.884  0.00000 0.000000 0.0000000 0.00000 0.0000
#> tlambdap 0.0000 0.000  0.01936 0.000000 0.0000000 0.00000 0.0000
#> tkqpp    0.0000 0.000  0.00000 0.001085 0.0000000 0.00000 0.0000
#> tdeltaqp 0.0000 0.000  0.00000 0.000000 0.0018207 0.00000 0.0000
#> tgamma   0.0000 0.000  0.00000 0.000000 0.0000000 0.26973 0.0000
#> tkde     0.0000 0.000  0.00000 0.000000 0.0000000 0.00000 0.0792

Now we have a thetaMat to represent the uncertainty in the theta matrix, as well as the other pieces in the simulation. Typically you can put this information into your simulation with the thetaMat matrix.

With such large variability in theta it is easy to sample a negative rate constant, which does not make sense. For example:

Ribba2012 %>% # Use rxode2
ini(prop.sd = 0.05) %>%
et(time.units="months") %>% # Pipe to a new event table
et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
et(0, 250, by=0.5) %>% # Add some sampling times (not required)
rxSolve(nSub=2, nStud=2, omega=omega,
thetaMat=thetaMat,
dfSub=760, dfObs=95) %>% # Solve the simulation
plot(pt, q, qp, pstar) # Plot it, plotting the variables of interest

#> ℹ change initial estimate of `prop.sd` to `0.05`
#> unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
#> @(lsoda.c:750
#> Warning message:
#> In rxSolve_(object, .ctl, .nms, .xtra, params, events, inits, setupOnly = .setupOnly) :
#>  Some ID(s) could not solve the ODEs correctly; These values are replaced with NA.

To correct these problems you simply need to use a truncated multivariate normal and specify the reasonable ranges for the parameters. For theta this is specified by thetaLower and thetaUpper. Similar parameters are there for the other matrices: omegaLower, omegaUpper, sigmaLower and sigmaUpper. These may be named vectors, one numeric value, or a numeric vector matching the number of parameters specified in the thetaMat matrix.

In this case the simulation simply has to be modified to have thetaLower=0 to make sure all rates are positive:

Ribba2012 %>% # Use rxode2
  ini(prop.sd = 0.05) %>%
  et(time.units="months") %>% # Pipe to a new event table
  et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
  et(0, 250, by=0.5) %>% # Add some sampling times (not required)
  rxSolve(nSub=2, nStud=2, omega=omega,
          thetaMat=thetaMat,
          thetaLower=0, # Make sure the rates are reasonable
          dfSub=760, dfObs=95) %>% # Solve the simulation
  plot(pt, q, qp, sim) # Plot it, plotting the variables of interest
#> ℹ change initial estimate of `prop.sd` to `0.05`
#> Warning: multi-subject simulation without without 'omega'

Summarizing the simulation output

While it is easy to use dplyr and data.table to perform your own summary of simulations, rxode2 also provides this ability by the confint function.

## This takes a little more time; Most of the time is the summary
## time.

sim0 <- Ribba2012 %>% # Use rxode2
  ini(prop.sd=0.05) %>%
  et(time.units="months") %>% # Pipe to a new event table
  et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
  et(0, 250, by=0.5) %>% # Add some sampling times (not required)
  rxSolve(nSub=10, nStud=10, omega=omega,
          thetaMat=thetaMat,
          thetaLower=0, # Make sure the rates are reasonable
          dfSub=760, dfObs=95) %>% # Solve the simulation
  confint(c("pt","q","qp","sim"),level=0.90); # Create Simulation intervals
#> ℹ change initial estimate of `prop.sd` to `0.05`
#> Warning: multi-subject simulation without without 'omega'
#> ! in order to put confidence bands around the intervals, you need at least 2500 simulations
#> summarizing data...done

sim0 %>% plot() # Plot the simulation intervals

Simulating from a data-frame of parameters

While the simulation from matrices can be very useful and a fast way to simulate information, sometimes you may want to simulate more complex scenarios. For instance, there may be some reason to believe that tkde needs to be above tlambdap, therefore these need to be simulated more carefully. You can generate the data frame in whatever way you want. The internal method of simulating the new parameters is exported too.

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

Ribba2012 <- Ribba2012()
# Convert to classic rxode2 model with ini attached

r <- Ribba2012$simulationIniModel
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

pars <- rxInits(r)
pars <- pars[regexpr("(prop|eta)",names(pars)) == -1]
print(pars)
#>           k        tkde        tkpq       tkqpp    tlambdap      tgamma 
#>    1.00e+02    2.40e-01    2.95e-02    3.10e-03    1.21e-01    7.29e-01 
#>    tdeltaqp        tpt0         tq0 rxerr.pstar 
#>    8.67e-03    7.13e+00    4.12e+01    1.00e+00
## This is the exported method for simulation of Theta/Omega internally in rxode2
df <- rxSimThetaOmega(params=pars, omega=omega,dfSub=760,
                      thetaMat=thetaMat, thetaLower=0, nSub=60,nStud=60) %>%
    filter(tkde > tlambdap) %>% as_tibble()
## You could also simulate more and bind them together to a data frame.
print(df)
#> # A tibble: 2,040 × 17
#>        k  tkde   tkpq tkqpp tlambdap tgamma tdeltaqp  tpt0   tq0 rxerr.pstar
#>    <dbl> <dbl>  <dbl> <dbl>    <dbl>  <dbl>    <dbl> <dbl> <dbl>       <dbl>
#>  1   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  2   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  3   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  4   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  5   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  6   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  7   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  8   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#>  9   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#> 10   100 0.716 0.0295 0.529    0.161  0.887    0.412  7.82  41.9           1
#> # ℹ 2,030 more rows
#> # ℹ 7 more variables: eta.pt0 <dbl>, eta.q0 <dbl>, eta.lambdap <dbl>,
#> #   eta.kqp <dbl>, eta.kqpp <dbl>, eta.deltaqp <dbl>, eta.tkde <dbl>
## Quick check to make sure that all the parameters are OK.
all(df$tkde>df$tlambdap)
#> [1] TRUE
sim1 <- r %>% # Use rxode2
    et(time.units="months") %>% # Pipe to a new event table
    et(amt=1, time=50, until=58, ii=1.5) %>% # Add dosing every 1.5 months
    et(0, 250, by=0.5) %>% # Add some sampling times (not required)
    rxSolve(df)
## Note this information looses information about which ID is in a
## "study", so it summarizes the confidence intervals by dividing the
## subjects into sqrt(#subjects) subjects and then summarizes the
## confidence intervals
sim2 <- sim1 %>% confint(c("pt","q","qp","sim"),level=0.90); # Create Simulation intervals
#> ! in order to put confidence bands around the intervals, you need at least 2500 simulations
#> summarizing data...done
save(sim2, file = file.path(system.file(package = "rxode2"), "pipeline-sim2.rds"), version = 2)
sim2 %>% plot()