Covariate sensitivity analysis via simulation

Evaluate the impact of continuous and categorical covariates on model parameters and exposure metrics. For each covariate the function perturbs its value to selected quantiles (continuous) or observed categories (categorical) while holding the remaining covariates at their reference values, simulates from the pharmacometric model under parameter uncertainty, and summarises the resulting percent change relative to the reference.

Usage

sg_covsens_sim(
  fpath_i = NULL,
  ds_parest = NULL,
  ds_covs = NULL,
  model,
  stimes,
  et,
  est_covmat = NULL,
  npop = 5,
  cont_cov_l,
  cat_cov_l,
  quantiles = c(0.1, 0.9),
  outputs = "Cc",
  aggr = c("min", "max", "mean", "auc"),
  seed = NULL,
  ci = 95
)

Arguments

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.

ds_parest

Data.frame. Parameter estimates table with columns parameter and value. Required when fpath_i is NULL; must be provided together with ds_covs. Default is NULL.

ds_covs

data.frame. Subject-level covariate dataset (one row per subject) containing both continuous and categorical covariate columns. Required when fpath_i is NULL; must be provided together with ds_parest. Default is NULL.

model

GMO. RxODE2 model to simulate from.

stimes

numeric vector. Sampling time points for steady-state simulation (e.g. generated by a cycling function over dosing intervals).

et

data.frame. Event (dosing) table. Must contain at least columns id, time, amt, evid and Regimen. Additional columns such as ii, addl, dur and Dose are used when present.

est_covmat

Data.frame. Parameter estimation covariance matrix. The first column (X1) must list parameter names; remaining columns (named identically) form the symmetric variance–covariance matrix.

npop

integer. Number of population-level simulation replicates drawn from the parameter uncertainty distribution. Default is 5.

cont_cov_l

A named list. Each element defines one continuous covariate and must itself be a list where par_vec is required; other components are optional/permissible. Name of each element should be equal the name of the continuous covariate:

par_vec: Character vector. Model parameter(s) affected by this covariate (e.g. c("CL")).
UTNAME: Character or NULL. Column name of the untransformed (back-transformed) covariate (e.g. "AGE"). If NULL or NA, defaults to name of the covariate. Default is NULL.
REF: Character or numeric. Reference value for the covariate. Use "median" to derive from data, or a numeric value. Default is "median".
NICENAME: Character or NULL. Display label for plots and tables (e.g. "Age, years"). Default is NULL.

cat_cov_l

A named list with. Each element defines one categorical covariate and must itself be a list where par_vec is required; other components are optional/permissible. Name of each element should be equal the name of the categorical covariate:

par_vec: Character vector. Model parameter(s) affected by this covariate (e.g. c("ka")).Required.
NICENAME: Character or NULL. Display label. Default is NULL.
REF: Character or NULL. Reference category value. If NULL, the first factor level (alphabetically) is used. Default is NULL.

quantiles

A numeric vector of length 2. Lower and upper quantiles of the continuous covariate distribution to test. Default is c(0.1, 0.9).

outputs

character vector. Name(s) of the model output variable(s) to evaluate (e.g. "Cc" or c("Cc", "Effect")). Default is "Cc".

aggr

character vector. Exposure aggregation metric(s) applied over the simulation time grid. Allowed values: "min" (Cmin), "max" (Cmax), "mean" (Cmean), "auc" (AUC).

seed

integer. Seed for the random number generator. Default is NULL.

ci

integer. Confidence interval (CI) level as a percentage. Allowed values: 95, 90, 80, 70, 50. Default is 95.

Value

A named list of four elements:

$PARSENS: data.frame. Full sensitivity results for model parameters: percent change statistics (mean, median, percentiles) for each covariate–parameter combination, with columns NICEN, VAR, KEY, LAB, Type and summary statistics (mean through P975).
$SUMPARSENS: data.frame. Compact summary table with columns Parameter, Covariate, Cov. percentile, Cov. value, Mean and 90%CI.
$EXPSENS: data.frame. Sensitivity of exposure metrics (Cmin, Cmax, Cmean, AUC) to covariates, structured identically to $PARSENS.
$COVREF: data.frame. Reference values for each covariate used in the analysis, with columns COV, NICEN, REF_VALUE and REF_SOURCE ("reference" for user-specified values, "median" when REF = "median" was used for a continuous covariate).

Details

Two data-source modes are supported (mutually exclusive):

File mode — supply fpath_i (path to a saved object with GFO tables such as SUMTAB, CATAB, and COTAB).
Table mode — supply both ds_parest (parameter estimates) and ds_covs (covariate dataset).

The analysis proceeds in two stages for each covariate type (continuous and categorical):

Parameter sensitivity — the covariate is set to its quantile or category value and the model is evaluated at time zero (no ODE integration) to quantify the direct effect on each parameter.
Exposure sensitivity — the full time-course is simulated over stimes and exposure metrics (aggr) are computed.

Uncertainty is propagated by sampling npop parameter vectors from a multivariate normal distribution parameterised by the population estimates and the estimation covariance matrix (est_covmat).

Examples

# \donttest{
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
library(rxode2)
#> rxode2 5.1.2 using 2 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

# --- Covariate definitions (only par_vec is required) ---
cont_cov_l <- list(
  LG_AGE = list(par_vec = c("CL")),
  LG_WEIGHT = list(
    UTNAME = "WEIGHT",
    NICENAME = "Weight, kg",
    par_vec = c("Vd")
  )
)

cat_cov_l <- list(
  SEX = list(par_vec = c("ka")),
  CYP2C9 = list(
    REF = "1",
    NICENAME = "CYP2C9 genotype",
    par_vec = c("CL")
  )
)

# --- Dosing ---
ev_t_input <- tribble(
  ~id, ~time, ~ii, ~amt, ~addl, ~dur, ~evid, ~Regimen,        ~Dose,
  1,   0,     336, 10,   21,    0.5,  1,     "0.3 mg/kg Q2W", 0.3
)
# --- Model ---
model <- RxODE({
  # Doses in mg
  # Time in hours

  ### Parameter values
  # Typical values
  ka_pop = 0.073;
  Vd_pop = 14.8;
  CL_pop = 0.347;

  # Random effects
  omega_ka = 0;
  omega_Vd = 0;
  omega_CL = 0;

  # Covariate effect
  # Continuous
  beta_CL_LG_AGE = 0.49990114;
  beta_Vd_LG_WEIGHT = 0.60529433;

  # Categorical
  beta_CL_CYP2C9_1_2 = -0.339;
  beta_CL_CYP2C9_1_3 = -0.574;
  beta_CL_CYP2C9_2_2 = -1.079;
  beta_CL_CYP2C9_2_3 = -0.745;
  beta_CL_CYP2C9_3_3 = -2.13;

  beta_ka_SEX_1 = -0.12198035;

  # Residual error
  Cc_b = 0;

  # Transformations
  ka_tv = exp(ka_pop);
  Vd_tv = exp(Vd_pop);
  CL_tv = exp(CL_pop);

  CL_multiplier = 1.0;  # Default/reference
  ka_multiplier = 1.0;

  if (SEX == "1") {ka_multiplier = exp(beta_ka_SEX_1)}

  if (CYP2C9 == "1") {
    CL_multiplier = exp(beta_CL_CYP2C9_1_2);
  } else if (CYP2C9 == "2") {
    CL_multiplier = exp(beta_CL_CYP2C9_1_3);
  } else if (CYP2C9 == "3") {
    CL_multiplier = exp(beta_CL_CYP2C9_2_2);
  } else if (CYP2C9 == "4") {
    CL_multiplier = exp(beta_CL_CYP2C9_2_3);
  } else if (CYP2C9 == "5") {
    CL_multiplier = exp(beta_CL_CYP2C9_3_3);
  }

  ka = ka_tv*ka_multiplier*exp(omega_ka);
  Vd = Vd_tv*exp(beta_Vd_LG_WEIGHT * LG_WEIGHT + omega_Vd); #Vd_tv*exp(omega_Vd);
  CL = CL_tv*CL_multiplier*exp(beta_CL_LG_AGE * LG_AGE + omega_CL);


  ### Explicit functions
  Cc = Ac/Vd;

  ### Initial conditions
  Ad(0) = 0;
  Ac(0) = 0;

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

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

# --- Estimation covariance (mock) ---
pnames     <- parest$parameter
npar       <- length(pnames)
set.seed(1)
m_cov      <- matrix(0.02, npar, npar)
diag(m_cov) <- 0.05 + runif(npar, 0, 0.05)
m_cov      <- (m_cov + t(m_cov)) / 2
est_covmat <- as_tibble(cbind(X1 = pnames, as.data.frame(m_cov)))
names(est_covmat)[-1] <- pnames

# --- Simulation times (steady-state cycle 10) ---
ss_cycle <- 10
stimes_ss <- c(
  ss_cycle * 4 * 7 * 24 + c(seq(0, 23.5, 0.5), seq(24, 335, 1)),
  ss_cycle * 4 * 7 * 24 + 2 * 7 * 24 + c(seq(0, 23.5, 0.5), seq(24, 335, 1))
)

# --- Run ---
result <- sg_covsens_sim(
  fpath_i=NULL, ds_parest = parest, ds_covs = ds_covval,
  model = model, stimes = stimes_ss, et = ev_t_input,
  est_covmat = est_covmat, npop = 10,
  cont_cov_l = cont_cov_l, cat_cov_l = cat_cov_l,
  quantiles = c(0.2, 0.8), aggr = c("max"),
  outputs = "Cc"
)
print(result[["PARSENS"]])
#> # A tibble: 10 × 25
#>    NICEN VAR   KEY    mean median    min   max    sd   P025    P05    P10    P15
#>    <chr> <chr> <chr> <dbl>  <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
#>  1 LG_A… CL    20th… 1.69   0.998 0.123  5.82   3.03 0.140  0.157  0.191  0.220 
#>  2 LG_A… CL    80th… 1.39   1.08  0.229  3.04   2.14 0.230  0.231  0.233  0.285 
#>  3 Weig… Vd    20th… 1.13   1.17  0.173  2.38   1.71 0.186  0.198  0.223  0.363 
#>  4 Weig… Vd    80th… 1.43   1.00  0.168  4.22   2.20 0.264  0.361  0.554  0.602 
#>  5 SEX   ka    NA    1.31   0.713 0.0627 4.09   2.52 0.0853 0.108  0.153  0.230 
#>  6 CYP2… CL    NA    2.63   1.59  0.541  8.72   3.63 0.563  0.585  0.630  0.685 
#>  7 CYP2… CL    NA    1.71   1.90  0.281  3.07   1.96 0.281  0.281  0.282  0.500 
#>  8 CYP2… CL    NA    0.838  0.828 0.105  1.96   1.57 0.141  0.177  0.250  0.301 
#>  9 CYP2… CL    NA    0.824  0.719 0.270  2.72   1.72 0.278  0.286  0.303  0.321 
#> 10 CYP2… CL    NA    0.220  0.195 0.0285 0.558  1.18 0.0315 0.0345 0.0406 0.0500
#> # ℹ 13 more variables: P25 <dbl>, P75 <dbl>, P85 <dbl>, P90 <dbl>, P95 <dbl>,
#> #   P975 <dbl>, Regimen <chr>, COVVAL <dbl>, BTR <chr>, BCOVVAL <dbl>,
#> #   LAB <fct>, Type <fct>, CATDES <chr>
print(result[["SUMPARSENS"]])
#> # A tibble: 10 × 6
#>    Parameter Covariate       `Cov. percentile` `Cov. value`  Mean `95%CI`      
#>    <chr>     <chr>           <chr>             <chr>        <dbl> <chr>        
#>  1 CL        LG_AGE          20th perc.        0.7971       1.69  0.14, 5.62   
#>  2 CL        LG_AGE          80th perc.        1.1304       1.39  0.23, 3.01   
#>  3 Vd        Weight, kg      20th perc.        63.3         1.13  0.186, 2.29  
#>  4 Vd        Weight, kg      80th perc.        84.1         1.43  0.264, 3.8   
#>  5 ka        SEX             Category          1            1.31  0.0853, 4.09 
#>  6 CL        CYP2C9 genotype Category          0            2.63  0.563, 8.03  
#>  7 CL        CYP2C9 genotype Category          2            1.71  0.281, 2.95  
#>  8 CL        CYP2C9 genotype Category          3            0.838 0.141, 1.83  
#>  9 CL        CYP2C9 genotype Category          4            0.824 0.278, 2.35  
#> 10 CL        CYP2C9 genotype Category          5            0.22  0.0315, 0.533
print(result[["EXPSENS"]])
#> # A tibble: 10 × 25
#>    NICEN     VAR   KEY    mean median    min   max    sd  P025   P05   P10   P15
#>    <chr>     <chr> <chr> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1 LG_AGE    Cc_C… 20th… 1.36   0.827 0.328   3.51  2.07 0.358 0.388 0.448 0.506
#>  2 LG_AGE    Cc_C… 80th… 1.58   0.694 0.158   8.50  3.53 0.164 0.171 0.183 0.188
#>  3 Weight, … Cc_C… 20th… 2.10   1.94  0.552   4.62  2.27 0.580 0.608 0.664 0.739
#>  4 Weight, … Cc_C… 80th… 2.02   1.37  0.783   5.28  2.44 0.795 0.806 0.830 0.857
#>  5 SEX       Cc_C… NA    1.45   1.52  0.548   3.22  1.80 0.586 0.624 0.700 0.743
#>  6 CYP2C9 g… Cc_C… NA    1.23   0.612 0.223   3.93  2.24 0.250 0.277 0.331 0.418
#>  7 CYP2C9 g… Cc_C… NA    0.913  0.497 0.196   3.32  1.99 0.205 0.214 0.231 0.240
#>  8 CYP2C9 g… Cc_C… NA    2.34   2.48  0.0990  6.86  2.86 0.275 0.451 0.802 0.917
#>  9 CYP2C9 g… Cc_C… NA    0.825  0.786 0.0831  2.26  1.62 0.112 0.141 0.199 0.286
#> 10 CYP2C9 g… Cc_C… NA    1.95   1.55  0.520   5.77  2.56 0.530 0.539 0.559 0.698
#> # ℹ 13 more variables: P25 <dbl>, P75 <dbl>, P85 <dbl>, P90 <dbl>, P95 <dbl>,
#> #   P975 <dbl>, Regimen <chr>, COVVAL <dbl>, BTR <chr>, BCOVVAL <dbl>,
#> #   LAB <fct>, Type <fct>, CATDES <chr>

# }