Perform simulations from an rxode2 model

Runs simulations using a GMO with explicit R objects for parameters (theta), variance–covariance (thetamat), omega, and sigma. Argument layout is described in GSI.

Usage

sg_sim(
  model,
  et,
  stimes = NULL,
  outputs = NULL,
  theta = NULL,
  omega = NULL,
  sigma = NULL,
  thetamat = NULL,
  covs = NULL,
  npop = 1,
  nsub = 1,
  aggr = NULL,
  addcov = TRUE,
  keep = NULL,
  scale = NULL,
  covint = "locf",
  inits = NULL,
  byID = NULL,
  byPOP = NULL,
  shared = NULL,
  ncores = 1,
  atol = 1e-08,
  rtol = 1e-06,
  fpath_i = NULL,
  ...
)

Arguments

model

GMO. RxODE2 model to simulate from.

et

Data.frame. Event table; a component of GSI.

stimes

Numeric vector. Sampling time points; a component of GSI. Default is NULL.

outputs

Vector of strings. Names of the model variables to output. If NULL, all variables returned. Default is NULL.

theta

Named vector or data frame. Population or baseline parameter values; a component of GSI. For sensitivity analysis, parameters listed in params are replaced by sampled values. Default is NULL.

omega

Named matrix or vector. Inter-individual variability matrix; a component of GSI.

sigma

Named matrix. Residual error variance matrix or its Cholesky decomposition; a component of GSI. Default is NULL.

thetamat

Named variance–covariance matrix for parameters on the normal scale; a component of GSI. Default is NULL.

covs

A list of covariate structures. Each element should be a list containing covariate relationship definitions with the following structure:

PAR - string, name of the parameter to which covariate effect is applied
COVNAME - string, name of the covariate column in the dataset
FUNC - string, functional form of the covariate effect
TRANS - string, transformation applied to covariate
REF - numeric, reference value for categorical covariates
INIT - numeric, initial value for the covariate effect parameter
EST - logical, whether to estimate the covariate effect Can be NULL, if no covariates are used. Default is NULL

npop

Integer. Number of population replicates. Default is 1.

nsub

Integer. Number of subjects sampled per population (omega/sigma matrices per ID). Default is 1.

aggr

A string that specifies the aggregation type. Set to NULL for no aggregation, ID for aggregation by time, population and ID, NPOP for aggregation by population and time, TIME for aggregation by time. Default is NULL.

addcov

Logical. If TRUE, columns with covariate values will be added to the resulting dataset. Default is TRUE.

keep

Vector of strings. Columns of event table to keep in the output data frame. Default is NULL.

scale

A numeric named vector. Scaling for ODE parameters of the system. The names must correspond to the parameter identifiers in the ODE specification. Each of the ODE variables will be divided by the scaling factor. Default is NULL.

covint

String. Specifies the interpolation method for time-varying covariates. When solving ODEs it often samples times outside the sampling time specified in event table. When this happens, the time varying covariates are interpolated. Currently this can be:

"linear" interpolation, which interpolates the covariate by solving the line between the observed covariates and extrapolating the new covariate value
"locf" last observation carried forward
"NOCB" next observation carried backward. This is the same method that NONMEM uses
"midpoint" last observation carried forward to midpoint; next observation carried backward to midpoint

Default is "locf"

inits

A named vector specifying initial conditions for model variables; Default is NULL.

byID

logical. If TRUE, replicate populations/subjects per event-table ID. Default NULL is treated as TRUE.

byPOP

logical. When both thetamat and omega are used: if TRUE, replicate subjects by population. Default NULL is treated as TRUE.

shared

logical. If TRUE, one shared population draw for all IDs. Default NULL is treated as TRUE.

ncores

Integer. Number of cores used for calculations. Default is 1.

atol

A numeric absolute tolerance used by the ODE solver to determine if a good solution has been achieved. This is also used in the solved linear model to check if prior doses do no add anything to the solution. Default is 1e-8.

rtol

Numeric. A numeric relative tolerance used by the ODE solver to determine if a good solution has been achieved. This is also used in the solved linear model to check if prior doses do not add anything to the solution. Default is 1e-6.

fpath_i

String, GFO, or a named list with a GFO component (and optionally GCO). If a string is given, the path to a .RData or .json file with a fit object is expected.

...

optional arguments passed to rxode2::rxSolve (e.g. method, maxsteps).

Value

GSO: a data frame with simulation results (long format).

Details

Index columns ID, POPN, and sim.id are included depending on what is used:

No uncertainty, no variability: Only ID (and TIME, VAR, VALUE)
Only uncertainty (thetamat): POPN and ID
Only variability (omega): sim.id and ID
Uncertainty and variability: POPN, sim.id, and ID

Other columns:

TIME: Timepoint of the simulations
VAR: Names of the output variables
VALUE: Optional. Simulated value. Returned when no aggregation applied
mean, median, ...: Optional. Aggregated statistic of simulated values. Returned when aggregation is applied
COV1, ...COVn: Optional. Covariates value. Returned when addcov is TRUE
KEEP1, ...KEEPn: Optional. Columns that were specified in the keep argument

Event table can use either id or ID. Parameters can be fixed per ID by including parameter names (e.g. *_pop, omega_*) as columns in the event table.

Examples

# \donttest{
library(rxode2)
library(dplyr)
library(tibble)

# Bundled [GMO] and [GSI] (see ?gmo_pk1c, ?gsi_pk1c)
sim_gmo <- do.call(sg_sim, c(list(model = gmo_pk1c), gsi_pk1c))
#> Warning: multi-subject simulation without without 'omega'
head(sim_gmo)
#> # A tibble: 6 × 5
#>    POPN    ID  TIME VAR   VALUE
#>   <dbl> <int> <dbl> <chr> <dbl>
#> 1     1     1 0     Cc    0    
#> 2     1     1 0.603 Cc    0.189
#> 3     1     1 1.21  Cc    0.279
#> 4     1     1 1.81  Cc    0.318
#> 5     1     1 2.41  Cc    0.331
#> 6     1     1 3.02  Cc    0.331

# Base 1-compartment PK model (no covariate)
mod <- rxode2::rxode2({
  ka_pop = 0.1;
  Vd_pop = 10;
  CL_pop = 0.5;

  omega_ka = 0;
  omega_Vd = 0;
  omega_CL = 0;

  Cc_b = 0;
  ka_tv = exp(ka_pop);
  Vd_tv = exp(Vd_pop);
  CL_tv = exp(CL_pop);

  ka = ka_tv * exp(omega_ka);
  Vd = Vd_tv * exp(omega_Vd);
  CL = CL_tv * exp(omega_CL);

  Cc = Ac / Vd;

  Ad(0) = 0;
  Ac(0) = 0;

  d/dt(Ad) = -ka * Ad;
  d/dt(Ac) = ka * Ad - CL * Cc;

  Cc_ResErr = Cc * (1 + Cc_b);
})
#>  
#>  

# Model with covariate: WTBL on Vd (Vd_tv = exp(Vd_pop) * (WTBL/70)^beta_WTBL_Vd_pop)

mod_cov <- rxode2::rxode2({
  ka_pop = 0.1;
  Vd_pop = 10;
  CL_pop = 0.5;

  beta_WTBL_Vd_pop = 0;

  omega_ka = 0;
  omega_Vd = 0;
  omega_CL = 0;

  Cc_b = 0;

  ka_tv = exp(ka_pop);
  Vd_tv = exp(Vd_pop) * (WTBL/70)^beta_WTBL_Vd_pop;
  CL_tv = exp(CL_pop);

  ka = ka_tv * exp(omega_ka);
  Vd = Vd_tv * exp(omega_Vd);
  CL = CL_tv * exp(omega_CL);

  Cc = Ac / Vd;

  Ad(0) = 0;
  Ac(0) = 0;

  d/dt(Ad) = -ka * Ad;
  d/dt(Ac) = ka * Ad - CL * Cc;

  Cc_ResErr = Cc * (1 + Cc_b);
})
#>  
#>  

et_test <- tibble::tribble(
  ~ID, ~TIME, ~EVID, ~CMT, ~AMT,
  1, 0, 1, 1, 10,
  2, 0, 1, 1, 50
)

stimes_test <- seq(0, 24, 0.1)
output_test <- c("Ac", "Ad", "Cc", "Cc_ResErr")
theta_test <- c(ka_pop = log(0.1), Vd_pop = log(10), CL_pop = log(0.5))

thetamat_test <- matrix(c(0.05, 0.01, 0, 0.01, 0.05, 0, 0, 0, 0.05), nrow = 3, byrow = TRUE)
rownames(thetamat_test) <- colnames(thetamat_test) <- c("ka_pop", "Vd_pop", "CL_pop")

omega_test <- matrix(c(0.2, 0.05, 0, 0.05, 0.2, 0, 0, 0, 0.2), nrow = 3, byrow = TRUE)
rownames(omega_test) <- colnames(omega_test) <- c("omega_ka", "omega_Vd", "omega_CL")

sigma_test <- matrix(0.1); rownames(sigma_test) <- colnames(sigma_test) <- "Cc_b"

# No variability
sim1 <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
               outputs = output_test)

# Population uncertainty, single scenario (ID)
sim2a <- sg_sim(model = mod, et = dplyr::filter(et_test, ID == 1),
                stimes = stimes_test, theta = theta_test, thetamat = thetamat_test,
                npop = 10)

# Population uncertainty, multiple IDs, byID = TRUE (populations are replicated
# for each scenario (ID)),
# shared = TRUE (the same populations are used for each scenario (ID))
sim2b <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
               thetamat = thetamat_test, npop = 10, byID = TRUE, shared = TRUE)
#> Warning: multi-subject simulation without without 'omega'

# Population uncertainty, multiple IDs, byID = FALSE (one solve over full event table;
# npop must equal n(ID); ID - virtual subject)
sim2c <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
               thetamat = thetamat_test, npop = nrow(et_test), byID = FALSE)
#> Warning: multi-subject simulation without without 'omega'


# Between-subject variability (BSV), multiple IDs, byID = TRUE (subjects are
# replicated for each scenario (ID)), shared = FALSE (separate subjects per ID)
sim3a <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                omega = omega_test, nsub = 10, byID = TRUE, shared = FALSE)


# Between-subject variability (BSV), byID = FALSE (one solve over full event table;
# nsub must equal n(ID); ID - virtual subject)
sim3b <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                omega = omega_test, nsub = nrow(et_test), byID = FALSE)


# BSV + uncertainty: byID = TRUE (populations/subjects are replicated for each
# scenario (ID)), byPOP = TRUE (subjects are replicated for each population),
# shared = FALSE (separate populations/subjects per ID)
sim4a <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                thetamat = thetamat_test, npop = 10, omega = omega_test, nsub = 5,
                byID = TRUE, byPOP = TRUE, shared = FALSE)

# BSV + uncertainty: byID = TRUE (populations/subjects are replicated for each
# scenario (ID)), byPOP = FALSE (subjects are not replicated between population;
# npop = nsub)
sim4b <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                thetamat = thetamat_test, npop = nrow(et_test), omega = omega_test,
                nsub = nrow(et_test), byID = TRUE, byPOP = FALSE)

# BSV + uncertainty: byID = FALSE (one solve over full event table; nsub and/or
# npop must equal n(ID); ID - virtual subject), byPOP = TRUE (subjects are
# replicated for each population), shared = FALSE (separate populations/subjects
# per ID)
sim4c <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                thetamat = thetamat_test, npop = 10, omega = omega_test, nsub = 5,
                byID = FALSE, byPOP = TRUE, shared = FALSE)
#> !Critical condition: nsub = n(ID) OR npop = n(ID); otherwise byID = TRUE

# BSV + uncertainty: byID = FALSE (one solve over full event table; nsub and/or
# npop must equal n(ID); ID - virtual subject), byPOP = FALSE (subjects are not
# replicated between population; npop = nsub)
sim4d <- sg_sim(model = mod, et = et_test, stimes = stimes_test, theta = theta_test,
                thetamat = thetamat_test, npop = nrow(et_test), omega = omega_test,
                nsub = nrow(et_test), byID = FALSE, byPOP = FALSE)

# Parameter override in event table
sim4 <- sg_sim(model = mod, et = dplyr::filter(et_test, ID == 1) %>%
                                 dplyr::mutate(Vd_pop = 2),
               stimes = stimes_test, theta = theta_test, thetamat = thetamat_test,
               npop = 1)


# Covariate (WTBL): pass covs and keep so WTBL is used and returned
theta_cov <- c(theta_test, beta_WTBL_Vd_pop = 0.5)
sim5 <- sg_sim(model = mod_cov, et = dplyr::filter(et_test, ID == 1) %>%
                                     dplyr::mutate(WTBL = 70),
               stimes = stimes_test, theta = theta_cov, covs = c("WTBL"),
               keep = c("WTBL"),
               thetamat = thetamat_test, npop = 5)
# }