Skip to contents

Evaluate data from oncology dose finding using a joint BLRM

scenario_jointBLRM.Rd

Evaluates data scenarios consisting of observations of one or more monotherapy or two-drug combination therapy dose-finding trials and computes posterior toxicities for a trial of interest and a set of doses of interest.

For an introduction to the use of this function, see the introduction vignette: vignette("intro_jointBLRM", package = "decider").

If multiple scenarios need to be evaluated, consider using the function scenario_list_jointBLRM() instead, which is a parallelized wrapper that processes a list of data scenarios within the same setting via scenario_jointBLRM().

A description of the underlying model and methods are given in the section Details.

Usage

scenario_jointBLRM(
   data=NULL,
   historical.data=NULL,
   doses.of.interest,
   dose.ref1,
   dose.ref2,
   trials.of.interest,
   types.of.interest=NULL,
   esc.rule=c("ewoc", "loss", "dynamic.loss"),
   dosing.intervals = c(0.16, 0.33, 0.6),
   ewoc.threshold = 0.25,
   loss.weights = c(1, 0, 1, 2),
   dynamic.weights = rbind(c(0.32, 0, 0.32, 0.36),
                           c(0.29, 0, 0.31, 0.4),
                           c(0.27, 0, 0.33, 0.4),
                           c(0.2,  0, 0.3,  0.5)
   ),
   prior.mu = list(mu_a1 =  c(logit(0.33), 2),
                   mu_b1 =  c(0,           1),
                   mu_a2 =  c(logit(0.33), 2),
                   mu_b2 =  c(0,           1),
                   mu_eta = c(0,           1.121)
   ),
   prior.tau = list(tau_a1 =  c(log(0.25),  log(2)/1.96),
                    tau_b1 =  c(log(0.125), log(2)/1.96),
                    tau_a2 =  c(log(0.25),  log(2)/1.96),
                    tau_b2 =  c(log(0.125), log(2)/1.96),
                    tau_eta = c(log(0.125), log(2)/1.96)
   ),
   saturating = FALSE,
   probs = c(0.025, 0.5, 0.975),
   iter = 26000,
   warmup = 1000,
   refresh = 0,
   adapt_delta = 0.8,
   max_treedepth = 15,
   chains = 4,
   digits = 5,
   seed=sample.int(.Machine$integer.max, 1),
   path = NULL,
   file.name = NULL,
   plot.decisions = FALSE,
   plot.combi.heatmap = TRUE,
   plot.int.probs.loss = FALSE,
   plot.return = FALSE,
   plot.file.format = "pdf",
   plot.width,
   plot.height,
   plot.unit,
   output.scen.config = FALSE
)

Arguments

data

List that contains the data scenario to be evaluated. Can be NULL if the prior shall be computed. The data list should have the following named entries, all of which need to be vectors of the same length (length should be the number of cohorts in the data):

  • data$dose1
    Numeric vector, each entry must be non-negative. Entry \(i\) should provide the dose level of compound 1 administered to cohort \(i\) in the data. Use 0 or NA to state that compound 1 was not used during treatment of a cohort. Note: For each cohort, either data$dose1 or data$dose2 must be positive.

  • data$dose2
    Numeric vector, each entry must be non-negative. Entry \(i\) should provide the dose level of compound 2 administered to cohort \(i\) in the data. Use 0 or NA to state that compound 2 was not used during treatment of a cohort. Note: For each cohort, either data$dose1 or data$dose2 must be positive.

  • data$n.pat
    Numeric vector, each entry must be a non-negative integer. Entry \(i\) should provide the number of patients in cohort \(i\) in the data. 0 is interpreted as the cohort not having been treated yet. If data contains solely cohorts with 0 patients, the function will sample from the prior.

  • data$n.dlt
    Numeric vector, each entry must be a non-negative integer. Entry \(i\) should provide the number of DLTs in cohort \(i\) in the data. In particular, the values need to be smaller or equal to the patient number of the corresponding cohort.

  • data$trial
    Numeric or character vector. The entries should be trial names, i.e. indicators for the trial to which the cohort belongs. These can either be numbers or strings (both will be converted to numbers internally). Note: As mixed vectors of numbers and strings will be converted to strings, entries such as 1 and "1" will be interpreted as the same trial.

historical.data

Optional named list, must have the same structure as data. It is equivalent to include observations as data or historical.data, the latter is included solely for convenience. Trial names across data and historical data must be consistent, in the sense that observations with the same entries in data$trial, respectively historical.data$trial are interpreted to belong to the same trial.

Note that the argument is only included for convenience, as it allows to perform multiple inferences using the same set of historical data with changing new data, without having to include the non-changing (historical) observations in the data argument.

doses.of.interest

Numeric matrix with two rows and non-negative entries. Each column gives a dose combination of interest, where the first row provides the dose level of compound 1, the second row of compound 2. NA entries will be interpreted as "dose was not administered", i.e. are internally converted to 0. Setting both doses to zero in one column is not a valid dose combination and will throw an error.

dose.ref1

Numeric, must be positive. Reference dose for compound 1.

dose.ref2

Numeric, must be positive. Reference dose for compound 2.

trials.of.interest

Optional vector of numerical or character trial numbers/names, for which the posterior is to be computed. If missing, all trials included in the data argument are used instead, or, if there are no trials in the data and historical.data, the priors for the doses of interest are computed. Trial names should be consistent with the numbers given in the data argument. If trial names are included that do not have observations in the data argument, the predictive DLT probabilities (based on a MAP approach) will be computed instead for these trials (i.e., the posterior for these trials will only depend on the data observed in other trials).

types.of.interest

Optional character vector with one entry for each entry of trials.of.interest which specifies the trial type of the corresponding trial of interest. Supported trial types are "mono1", "mono2", "combi", and "all". If types.of.interest is not provided, the trial types are inferred from the data given for the trial. If the type is neither specified by the user nor uniquely inferable from the data (e.g. trials which are not included in the data yet, or trials including observations from multiple types), the trial type is implicitly set to "all".

The available trial types specify which of the dose levels from doses.of.interest are included in the summary for the corresponding trial of interest. More specifically, if the type of a trial of interest is...

  • ..."mono1": Only doses of interest for which the second compound is 0 are included in the summary and (if activated) plots.

  • ..."mono2": Only doses of interest for which the fist compound is 0 are included in the summary and (if activated) plots.

  • ..."combi": Only doses of interest for which none of the compound is 0 are included in the summary and (if activated) plots.

  • ..."all": All doses of interest are included in the summary. If plots are activated, a plot is generated for each trial type detected within the doses of interest.

If the type is "all" and the escalation rule is "dynamic.loss", the dynamic Bayes risk (expected loss) is computed based on the reference dose of the trial type corresponding to the type of trial of interest. That is, if e.g. mono 1 and combination doses are included in the doses of interest and the type is "all", the dynamic Bayes risk will select for each dose whether it needs to use mono reference doses or the reference dose combination.

esc.rule

Optional character. Can be either "ewoc", "loss", "dynamic" or "dynamic.loss", where the latter two are treated synonymously. Note: If "loss" or "dynamic.loss" are chosen, the parameter dosing.intervals needs to contain three entries, otherwise it can either have two or three entries. If esc.rule is "ewoc", at most the first two entries of dosing.intervals are used. See the section Details for a description of the available methods.

dosing.intervals

Optional numeric with 1, 2, or 3 ascending positive entries. Must have three entries when the esc.rule is set to "loss", "dynamic.loss", or "dynamic", otherwise (i.e. when esc.rule is "ewoc") lengths 1, 2, or 3 are permitted.
If dosing interval has only one entry (requires that esc.rule is "ewoc"), the entry is interpreted as the boundary between DLT rates considered as target dosing and overdosing, i.e., no underdosing interval is considered in this case. Otherwise (including the default), entries 1 and 2 are interpreted as lower and upper boundary of the target dosing interval. A third entry is ignored when esc.rule is "ewoc", and for other values of esc.rule required as boundary between the dosing intervals containing the excessive and unacceptable DLT rates. The argument dosing.intervals defaults to c(0.16, 0.33, 0.6), which corresponds to using \([0.16, 0.33)\) as target interval, and \([0.33,0.6)\) as excessive toxicity interval for loss escalation. If dosing.intervals has only one entry, say x, the computation is equivalent to setting the argument to c(0, x), but the output tables and plots will be formatted slightly differently to reflect that the underdosing interval is not considered in this case.

ewoc.threshold

Optional numeric between 0 and 1. Overdosing thresholds for EWOC plots. Defaults to 0.25.

loss.weights

Optional numerical vector with four entries (which can be arbitrary numbers), the default is c(1,0,1,2). Specifies the interval weights/penalties that are used for static loss escalation. This is only needed when esc.rule is "loss", and the parameter is ignored otherwise. More precisely, loss.weights[1] is the weight of the underdosing interval, loss.weights[2] the weight of the target dosing interval, loss.weights[3] the weight of the excessively toxic dosing interval, and loss.weights[4] is the weight of the unacceptably toxic dosing interval. See the Details section below for explanations regarding the methodology.

dynamic.weights

Optional numerical matrix with four rows and four columns, and arbitrary numbers as entries. Specifies the interval weights/penalties that are used for dynamic loss escalation. This is only needed when esc.rule is "dynamic" or "dynamic.loss", and the parameter is ignored otherwise. Each row of dynamic.weights specifies one of the static loss weight vectors that are interpolated during dynamic loss escalation. See the Details section below for explanations regarding the methodology.

More precisely, dynamic.weights[1, ] is the static weight vector that is weighted by the reference probability of underdosing, dynamic.weights[2, ] is the static weight vector that is weighted by the reference probability of target dosing, dynamic.weights[3, ] is the static weight vector that is weighted by the reference probability of excessively toxic dosing, and dynamic.weights[4, ] is the static weight vector that is weighted by the reference probability of unacceptably toxic dosing. The default value for dynamic.weights is a matrix with the following four rows (in order):

  • c(0.32, 0, 0.32, 0.36),

  • c(0.29, 0, 0.31, 0.40),

  • c(0.27, 0, 0.33, 0.40),

  • c(0.20, 0, 0.30, 0.50).

prior.mu

Optional list that gives the prior distribution for the hyper means \(\mu\). The list must have named entries, which all need to be numeric vectors of length 2:

  • prior.mu$mu_a1
    Numeric with length two, defaults to c(logit(0.33), 2). Specifies mean and SD of the hypermean \(\mu_{\alpha_1}\) of the parameter \(log(\alpha_1)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.mu$mu_b1
    Numeric with length two, defaults to c(0, 1). Specifies mean and SD of the hypermean \(\mu_{\beta_1}\) of the parameter \(log(\beta_1)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.mu$mu_a2
    Numeric with length two, defaults to c(logit(0.33), 2). Specifies mean and SD of the hypermean \(\mu_{\alpha_2}\) of the parameter \(log(\alpha_2)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.mu$mu_b2
    Numeric with length two, defaults to c(0, 1). Specifies mean and SD of the hypermean \(\mu_{\beta_2}\) of the parameter \(log(\beta_2)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.mu$mu_eta
    Numeric with length two, defaults to c(0, 1.121). Specifies mean and SD of the hypermean \(\mu_{\eta}\) of the parameter \(\eta\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

prior.tau

Optional list that gives the prior distribution for the between-trial heterogeneities (hyper SD) \(\tau\). The list must have named entries, which all need to be numeric vectors of length 2:

  • prior.tau$tau_a1
    Numeric with length two, defaults to c(log(0.25), log(2)/1.96). Specifies mean and SD on log-scale of the heterogeneity \(\tau_{\alpha_1}\) of the parameter \(log(\alpha_1)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.tau$tau_b1
    Numeric with length two, defaults to c(log(0.125), log(2)/1.96). Specifies mean and SD on log-scale of the heterogeneity \(\tau_{\beta_1}\) of the parameter \(log(\beta_1)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.tau$tau_a2
    Numeric with length two, defaults to c(log(0.25), log(2)/1.96). Specifies mean and SD on log-scale of the heterogeneity \(\tau_{\alpha_2}\) of the parameter \(log(\alpha_2)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.tau$tau_b2
    Numeric with length two, defaults to c(log(0.125), log(2)/1.96). Specifies mean and SD on log-scale of the heterogeneity \(\tau_{\beta_2}\) of the parameter \(log(\beta_2)\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

  • prior.tau$tau_eta
    Numeric with length two, defaults to c(log(0.125), log(2)/1.96). Specifies mean and SD on log-scale of the heterogeneity \(\tau_{\eta}\) of the parameter \(\eta\) in the BLRM. The second entry must therefore be positive. See the Details section below for more detail.

saturating

Optional logical, defaults to FALSE. If TRUE, the BLRM will be using a saturating interaction term as described in OncoBayes2::blrm_formula_saturating(). Also refer to the Details section below.

probs

Optional numeric with arbitrary entries between 0 and 1. Provides the levels for the quantiles displayed in the output. Defaults to c(0.025, 0.5, 0.975).

iter

Optional integer, number of total MCMC iterations per chain. Defaults to 26000. Note: Number of warmup iterations is counted towards iter, i.e. of the iter many iterations, the first warmup many samples are not saved.

warmup

Optional integer, number of warmup iterations discarded from total MCMC iterations per chain. Defaults to 1000. Note: Number of warmup iterations is counted towards iter, i.e. of the iter many iterations, the first warmup many samples are not saved. Refer to rstan::sampling() from the rstan-package package for details.

refresh

Optional integer. Given to Stan's refresh argument for rstan::sampling(), defaults to 0. A positive value indicates that updates on sampling progress are printed every refresh many iterations. Set this to 0 to avoid output from Stan. Refer to rstan::sampling() from the rstan-package package for details.

adapt_delta

Optional numeric between 0.6 and 1, default is 0.8. Given to Stan's rstan::sampling() method in the control argument of Stan. Translates to target Metropolis acceptance probability during Hamiltonian Monte Carlo, respectively NUTS. Can be used to influence the stepsize Stan uses for leapfrog steps during the NUTS algorithm. See rstan::sampling(), rstan-package for details. Note: One can increase the value of this parameter in order to reduce the number of Stan warnings on "divergent transition". Larger values than 0.95 will slow down sampling. The function will not permit to set adapt_delta to values below 0.6 (the rstan-package default).

max_treedepth

Optional integer, defaults to 15. Should not be lower than 10. This is a parameter of the NUTS algorithm. Roughly speaking, NUTS constructs a search tree when generating a new proposal, until a stopping criterion (the NUTS criterion) is satisfied or, alternatively, until the maximum depth of the search tree is reached (to avoid endless looping). The argument max_treedepth allows to control the latter. The maximum treedepth will also automatically constrain the maximum number of leapfrog steps, and should therefore not be below 10. See also rstan::sampling() for details on max_treedepth. Note: The maximum treedepth should usually not be reached during sampling, otherwise Stan will print a warning. In this case, it is recommended to increase its value.

chains

Optional integer. Number of Markov chains constructed by Stan. Defaults to 4.

digits

Optional positive integer that specified the number of digits used in the outputs. Defaults to 5.

seed

Optional positive integer that specifies the seed to be used for the simulation. The default is sample.int(.Machine$integer.max, 1). Note that reproducibility can only be obtained when the function is executed on exactly the same computing architecture, using identical software versions (e.g. of the compiler, Stan, and R), and the same input specifications. This is due to the internal use of rstan-package for MCMC sampling, which is only reproducible under these restrictions (refer to the Stan reference manual from rstan-package's homepage for more detail).

path

Optional character that specified the path to save the resulting output files. If NULL (the default), no output is written (but still returned to R). Otherwise, it is checked whether path specifies a directory, and, if yes, all output is saved there.

file.name

Optional name for the output file. If NULL or missing, no output is saved (and also if no valid path is given). Results are only returned to R in this case. Note that plots will not be returned to R unless the additional argument plot.return is specified.

plot.decisions

Optional logical, defaults to FALSE. If TRUE, plots of escalation decisions according to the specified escalation rule are created. If plot.decisions is TRUE, it is recommended to supply a writable path and valid file.name, as plots will by default only be saved and not returned to R. Plots are only returned to R if the additional argument plot.return is specified. In particular, if plot.decisions is TRUE and plot.return is FALSE, the plots will not be accessible unless a valid file.name is given to enable saving the plots.

Depending on the value of esc.rule, the following plots are created:

  • If esc.rule=="ewoc", a diagram of the posterior probabilities that the DLT rate lies in the underdosing, target dosing and overdosing intervals is created, which indicates the boundary given in ewoc.threshold with red color coding (and, therefore, implicitly the decision of EWOC-based dose escalation rules). For trials of types "mono1" and "mono2" (cf. argument types.of.interest), a bar plot is created, and for trials of type "combi", either a heatmap of the overdosing probabilities can be created (if the value of plot.combi.heatmap is TRUE) or alternatively a bar plot in the same fashion as monotherapy trials.

  • If esc.rule=="loss", a diagram of the expected loss (Bayes risk) based on the weight choice in loss.weights is created. Further, the argument plot.int.probs.loss can be used to additionally create plots for the underdosing, target dosing and overdosing intervals. For trials of types "mono1" and "mono2" (cf. argument types.of.interest), a bar plot is created, and for trials of type "combi", either a heatmap of the expected loss can be created (if plot.combi.heatmap==TRUE) or alternatively a bar plot in the same fashion as monotherapy trials.

  • If esc.rule=="dynamic.loss", a diagram of the expected dynamic loss (dynamic Bayes risk) based on the weight choice in dynamic.weights is created. Further, the argument plot.int.probs.loss can be used to additionally create plots for the underdosing, target dosing and overdosing intervals. For trials of types "mono1" and "mono2" (cf. argument types.of.interest), a bar plot is created, and for trials of type "combi", either a heatmap of the expected dynamic loss can be created (if plot.combi.heatmap==TRUE) or alternatively a bar plot in the same fashion as monotherapy trials.

plot.combi.heatmap

Optional logical, defaults to TRUE. If the value is TRUE, combination therapy plots are created as heatmaps instead of bar plots. This affects all escalation rules.

plot.int.probs.loss

Optional logical, defaults to FALSE. Only has an effect if esc.rule is either "loss" or "dynamic.loss". In this case, if the value is TRUE, additional plots will be created that display the interval probabilities to complement the (always created) plots that display the resulting expected loss.

plot.return

Optional logical, defaults to FALSE. If set to TRUE, the functions return the created plots to R in the result list. In this case, the result list will have an additional entry, ...$plots, which contains the posterior plots as a named list. The returned plots will be given as a list of ggplot2::ggplot() objects, which can be displayed in R via the ggplot2::plot.ggplot method. Refer to the documentation of the ggplot2-package package for further detail.

plot.file.format

Optional character, defaults to "pdf". Can either be "pdf", "jpeg" or "png". Specifies the file format of the created output plots (if plots are saved). Note that "pdf" will use grDevices::cairo_pdf().

plot.width

Optional numerical value or vector, can have length 1 or 3 and must have positive entries. Provides the width of the output plots measured in the unit given in plot.unit. If plot.width is a single value, all plots will use the same width. Otherwise, if plot.width has length 3, the entry plot.width[1] is used as width for plots for trials of type "mono1", while plot.width[2] is used for trials of type "mono2" and plot.width[3] for trials of type "combi.".
plot.width has only an effect if all of the parameters plot.unit, plot.width, and plot.height are defined. In this case, The output plots will be created with the specified width and length (interpreted in the given unit of measurement). If one of these parameters is missing, the default device size is used.

plot.height

Optional numerical value or vector, can have length 1 or 3 and must have positive entries. Provides the height of the output plots measured in the unit given in plot.unit. If plot.height is a single value, all plots will use the same height. Otherwise, if plot.height has length 3, the entry plot.height[1] is used as height for plots for trials of type "mono1", while plot.height[2] is used for trials of type "mono2" and plot.height[3] for trials of type "combi.".
plot.height has only an effect if all of the parameters plot.unit, plot.width, and plot.height are defined. In this case, The output plots will be created with the specified width and length (interpreted in the given unit of measurement). If one of these parameters is missing, the default device size is used.

plot.unit

Optional character string, can be "in", "cm", or "mm". Provides the unit of measurements for plot.width and plot.height. Has only an effect if all of the parameters plot.unit, plot.width, and plot.height are defined. In this case, The output plots will be created with the specified width and length (interpreted in the given unit of measurement). If one of these parameters is missing, the default device size is used.

output.scen.config

Optional logical, defaults to FALSE. If TRUE, the results will contain the necessary specifications of the evaluated scenario, i.e., the data scenario, priors, escalation rule specification, and seed. Otherwise, only the posterior summaries are returned.

Value

List. The output list will have an entry for each trial of interest that provides a summary of the posterior toxicities. If output.scen.config

is TRUE, additional entries that give the input data are included.

If plot.return and plot.decisions are both TRUE, an additional entry is created that holds a list of all output plots as ggplot2::ggplot objects.

More precisely, the following list entries are always generated:

  • $trial-[...]
    Here, [...] denotes the given trial name in trials.of.interest. The entry gives a matrix that lists summary statistics and interval probabilities based on the posterior of the DLT rate of each of the doses of interest for a trial. Additionally, if loss escalation is active, the expected loss of each dose is listed.

If additionally output.scen.config is active, there will be the following additional entries:

  • $data
    Input data used to fit the joint BLRM. Includes all cohorts from both data and historical.data in a merged data matrix.

  • $prior
    Contains the specified (hyper-)prior distribution used by the joint BLRM.

  • $configuration
    Contains the remaining configurations, e.g. seed and escalation rule.

  • $'Stan options'
    Contains the arguments given to Stan, e.g. number of MCMC iterations and chains.

If the additional plot list is generated, there will be an entry:

  • $plots
    List that contains the output plots as ggplot2::ggplot objects. There will be one entry for each trial, namely

    • $plots$trial-[...]
      Here, [...] denotes the trial name of one of the trials of interest. The entry for a trial is either a single ggplot2::ggplot object or a list of such objects. The number of entries is determined by the remaining specifications for plots. E.g., if loss escalation is performed, there can be either just a plot of the expected loss for each dose, or additionally a second plot with the usual interval probabilities. Similarly, if the trial of interest is of type "all", there will be up to 3 plots, one for each trial type for which dose levels are included (mono 1, mono 2, combination).

Details

The joint BLRM is defined according to (Neuenschwander et al., 2014 and 2016). It allows to perform Bayesian logistic regression to estimate the dose-toxicity relationship of two different monotherapies and combination therapy with these compounds in a joint model, which includes a hierarchical prior for robust borrowing across trials. The following gives a brief introduction to the model, for more detail refer to the references provided at the end of this documentation page.

Model specification

The model assumes, first of all, that treatment with a dose or dose combination \(d\) is during trial \(j\) governed by its toxicity parameter \(\pi_j(d)\), so that $$r|\pi_j(d) \sim Binomial(n, \pi_j(d))$$ for the number of DLTs, \(r\), observed after treating a cohort of \(n\) patients at dose \(d\) during trial \(j\). The toxicity is then modeled based on the type of treatment. If \(d_i\) is a dose of compound \(i\), it is assumed that $$logit(\pi_j(d)) = log(\alpha_{ij}) + \beta_{ij}\cdot log\left(\frac{d_i}{d^*_i}\right)$$ for unknown parameters \(\alpha_{ij}>0\) and \(\beta_{ij}>0\) dependent on the trial \(j\) and compound \(i\). Moreover, \(d^*_i\) denotes a fixed reference dose for compound \(i\).

Further, if \(d=(d_1, d_2)\) is a dose combination of compounds 1 and 2, write \(\pi_j(d_i)\) for their toxicity according to the monotherapy model, and define the interaction-free toxicity in trial \(j\) as $$\pi_{0j}(d_1, d_2) = \pi_j(d_1) + \pi_j(d_2) - \pi_j(d_1)\cdot \pi_j(d_2).$$ Based on this, an additional interaction parameter, \(\eta_j\), is introduced, and the combination therapy toxicity with interaction is modeled as $$logit(\pi_j(d_1, d_2)) = logit(\pi_{0j}(d_1, d_2)) + \eta_j \cdot \frac{d_1}{d^*_1} \cdot \frac{d_2}{d^*_2}.$$ This concludes the definition of the likelihood of the model. In the following, we denote the interaction term as the function $$ g(\eta, d_1, d_2) = \eta \cdot \frac{d_1}{d^*_1} \cdot \frac{d_2}{d^*_2}.$$ This is the so-called linear (on logit scale) interaction term.

Saturating interaction term

Another possibility for modeling the interaction is the so-called saturating interaction term, which was proposed by the authors of the OncoBayes2 package, cf. OncoBayes2::blrm_formula_saturating(). This interaction model involves an additional scaling factor in dependence of the dose, with the aim of addressing potentially problematic behavior of the linear interaction model in case of large dose levels. More specifically, if one has \(\eta<0\) and considers the limit \(d_i\rightarrow \infty\), it holds $$ g(\eta, d_1, d_2) \rightarrow -\infty $$ and, consequently, \(\pi_j(d_1, d_2) \rightarrow 0\). Although this seems to be mostly affecting doses larger than the reference dose, this property may not be considered realistic, so that one may want to use modified models that do not show this behavior. For this aim, the functions in the decider package allow to use a saturating interaction term instead, which follows the ideas from OncoBayes2::blrm_formula_saturating(), to avoid the aforementioned phenomenon.

The saturating interaction term is based on the usual interaction term, but includes an additional scaling factor in dependence of the dose level, which ensures a finite limit when letting \(d\rightarrow \infty\). More precisely, the linear interaction term,\(g(\eta, d_1, d_2)\), is exchanged with the saturating term $$ \tilde g (\eta, d_1, d_2) = g(\eta, d_1, d_2) \cdot \frac{2}{1+\frac{d_1}{d^*_1} \cdot \frac{d_2}{d^*_2}} $$ The scaling factor is defined so that whenever \(d_i=d^*_i\) for both \(i\), the saturating interaction term is equal to the usual interaction term, while it holds \(\tilde g (\eta, d_1, d_2)\rightarrow 2\eta\) for \(d_i\rightarrow \infty\).

Prior structure and prior choice

Regarding the prior configuration of the joint BLRM, denote by $$\theta_j=(log(\alpha_{1j}), log(\beta_{1j}), log(\alpha_{2j}), log(\beta_{2j}), \eta_j)$$ the trial-specific parameter vector for trial \(j\). For the general hierarchical prior assume that $$\theta_j|\mu,\Sigma \sim Normal_5(\mu, \Sigma)$$ for a shared hyper mean vector \(\mu\) and a shared hyper covariance matrix \(\Sigma\), with entries $$\mu=(\mu_1, \mu_2, \mu_3, \mu_4, \mu_5)$$ and $$\Sigma = (\Sigma_{kl})$$ for \(k=1,...,5\), \(l=1,...,5\). The entries \(\Sigma_{kl}\) are defined as $$\Sigma_{kk}=\tau_k\cdot\tau_k$$ on the diagonal, and as $$\Sigma_{kl}=\rho_{kl}\cdot\tau_k\cdot\tau_l$$ for \(k\) not equal to \(l\). Here, \(\rho_{kl}\) are the correlation coefficients of parameters \(k\) and \(l\), and \(\tau_k\) the standard deviations of their respective parameters. As the \(\tau_k\) control the standard deviation across the parameters of different trials, these parameters are referred to as between-trial heterogeneities.

Regarding the prior for the model, it suffices to specify the distributions of \(\mu_k\), \(\tau_k\) and \(\rho_{kl}\). As a simplification, it is further assumed that all entries of \(\theta_j\) are uncorrelated, except for \(log(\alpha_{ij})\) and \(log(\beta_{ij})\) (i.e., the parameters belonging to the same monotherapy compound). That is, \(\rho_{kl}=0\) for all \(kl\), except for \(\rho_{12}\) and \(\rho_{34}\). As hyperprior distribution of the correlations, the two remaining correlation coefficients \(\rho_{12}\) and \(\rho_{34}\) are assumed to follow \(Uniform(-1, 1)\) distributions. Note that this cannot be changed in the functions scenario_jointBLRM(), sim_jointBLRM(), and fit_jointBLRM().

For the hyper means \(\mu_k\), it is assumed that $$\mu_k \sim Normal(m_{\mu_k}, s_{\mu_k}^2)$$ for fixed mean \(m_{\mu_k}\) and standard deviation \(s_{\mu_k}\). These values are specified in the prior.mu argument of scenario_jointBLRM() and sim_jointBLRM(). For simplicity, the entries of the list prior.mu are vectors of length two, which contain the fixed mean \(m_{\mu_k}\) and fixed SD \(s_{\mu_k}\). These entries are named for their corresponding parameters, as given in the following overview. For reasons of readability, we will use the corresponding parameter names as subscript by writing e.g. \(\mu_{\alpha_1}\) for the hypermean of \(log(\alpha_1)\). This is summarized in the following table.

Entrymu_a1mu_b1mu_a2mu_b2mu_eta
Parameter\(log(\alpha_{1j})\)\(log(\beta_{1j})\)\(log(\alpha_{2j})\)\(log(\beta_{2j})\)\(\eta_j\)
Hyper Mean\(\mu_1=\mu_{\alpha_1}\)\(\mu_2=\mu_{\beta_1}\)\(\mu_3=\mu_{\alpha_2}\)\(\mu_4==\mu_{\beta_2}\)\(\mu_{\eta}\)

For the hyper SDs \(\tau_k\), it is assumed that $$\tau_k \sim logNormal(m_{\tau_k}, s_{\tau_k}^2)$$ for fixed mean \(m_{\tau_k}\) and standard deviation \(s_{\tau_k}\), which are both to be given on log-scale (that is, \(\tau_k=exp(y)\) for a random variable \(y\) which has a \(Normal(m_{\tau_k}, s_{\tau_k}^2)\) distribution). These values are specified in the prior.tau argument of scenario_jointBLRM() and sim_jointBLRM(). For simplicity, the entries of the list prior.tau are vectors of length two, which contain the fixed mean \(m_{\tau_k}\) and fixed SD \(s_{\tau_k}\) on log-scale. These entries are named for their corresponding parameters, as given in the following overview.

Entrytau_a1tau_b1tau_a2tau_b2tau_eta
Parameter\(log(\alpha_{1j})\)\(log(\beta_{1j})\)\(log(\alpha_{2j})\)\(log(\beta_{2j})\)\(\eta_j\)
Hyper SD\(\tau_1=\tau_{\alpha_1}\)\(\tau_2=\tau_{\beta_1}\)\(\tau_3=\tau_{\alpha_2}\)\(\tau_4=\tau_{\beta_2}\)\(\tau_5=\tau_{\eta}\)

See (Neuenschwander et al., 2014 and 2016) for details regarding the prior choice.

Escalation rules

The function supports in general three different escalation rules and corresponding displays of the results. In this context, the term "escalation rule" refers to a method that is used to derive concrete dosing recommendations from the posterior computed by the BLRM. The supported escalation rules are the escalation with overdose control (EWOC) principle, static loss escalation, and dynamic loss escalation.

Escalation with overdose control (EWOC)

The EWOC principle goes back to Babb et al. (1998) and derives dosing recommendations based on the posterior probability for overdosing, where overdosing refers to DLT rates above the specified target interval. By default, the overdosing interval starts therefore at 0.33. Formally, EWOC considers the posterior distribution \(\pi_j(d)\) of the DLT rate for a dose \(d\) and demands that the posterior probability of \(\pi_j(d)\) lying in the overdosing interval is lower than a prespecified feasibility bound (which defaults to 0.25). As a formula, the EWOC condition is therefore for a dose \(d\) during trial \(j\) (using the default values 0.33 and 0.25 for lower boundary of the overdosing interval and feasibility bound): $$ Prob(\pi_j(d) \ge 0.33) < 0.25. $$ Only doses that satisfy this criterion are then considered as recommendations for the next dose. To obtain a concrete dose among the remaining levels, one typically demands an additional rule. Most commonly, the "optimal probability rule" (selecting the dose that has the largest probability of having a DLT rate within the target interval among the levels that satisfy EWOC) or the "maximal dose rule" (selecting the largest dose that satisfies EWOC) are chosen. For dose combinations, one of the compounds should additionally be selected to be maximized first in case of draws.

Loss function-based escalation

Static loss escalation goes back to Neuenschwander et al. (2008) and applies Bayesian decision theory to obtain dosing recommendations. That is, the selection of a dosing recommendation is viewed as a decision problem among different available dose levels and different decisions (dose levels) are evaluated using a formal loss function. To define the loss function, the range of DLT rates is divided in four disjoint dosing intervals, namely, underdosing (interval \(I_1\)), target dosing (\(I_2\)), excessively toxic dosing (\(I_3\)), and unacceptably toxic dosing (\(I_4\)). In contrast to EWOC, this corresponds to dividing the overdosing interval into two distinct intervals to achieve a finer differentiation.

To define a formal loss function based on these intervals, a pre-specified penalty value (or interval weight), \(w_i\), is assigned to each of the dosing intervals. The loss function is then defined to be \(w_i\) for the decision on a dose with DLT rate in interval \(I_i\). As it is desired to select doses with DLT rate in the target interval, one should typically assume that \(w_2=0\) and that the remaining penalties are positive numbers. The dosing recommendation is now achieved by calculating for each fixed dose level the expected value of the loss function, the so-called Bayes risk (or expected loss). By definition of the loss function, the Bayes risk of a dose \(d\) during trial \(j\) is $$ \sum_{i=1}^4 w_i \cdot Prob(\pi_j(d) \in I_i). $$ The recommended decision is then the dose with minimal Bayes risk among the possible doses for the next cohort.

In terms of the pre-specified weights \(w_i\), one should always set \(w_2 = 0\) for the target interval, but for the remaining weights manual tuning may be required to cover different situations. Values that often work are e.g. \(w_1 = 1\), \(w_3 = 1\), \(w_4 = 2\) for moderate to aggressive escalation decisions, or \(w_1 = 1\), \(w_3 = 2\), \(w_4 = 3\) for rather conservative escalation decisions. Note that \(w_1 = 1\), \(w_3 = 1\), \(w_4 = 1\) leads to a decision rule that simply maximizes the probability to lie in the target interval, which was found to lead to overly aggressive decisions in simulations according to Zhou et al. (2018).

Dynamic adjustment of loss function-based escalation

Dynamic loss escalation refers to a heuristic that aims to adapt the pre-specified loss weights of the method described in the previous subsection. The motivation for this are potential difficulties with respect to pre-specifying weights that ensure appropriate behavior in different situations (larger or lower true DLT rates), as well as the attempt to obtain an escalation rule that behaves more robustly across scenarios, in the sense that it obtains good performance in settings of very low DLT rates as well as in settings of very large ones. Static loss escalation often needs to be tuned to one of these settings to obtain good on-trial behavior, but may perform poorly in other settings due to this. Dynamic loss escalation tries to avoid this by adapting the loss weights to more aggressive or more conservative configurations using a heuristic based on the observed data and posterior.

To apply dynamic loss escalation, four fixed loss weight combinations $$W_k=(w_{k1}, w_{k2}, w_{k3}, w_{k4})$$ for \(k=1,...,4\) with different degrees of aggressiveness need to be selected. The weights should be normalized in the sense that the absolute values of each \(W_k\) must sum up to the same number (i.e., have identical \(l_1\) norms). As decisions from static loss functions are invariant under linear rescaling of the weights, one can assume without loss of generality that the weights sum up to 1.

For the weight adjustment, new loss weights are computed prior to every decision by interpolating between the fixed weights \(W_k\) based on the current posterior interval probabilities of the reference dose of the trial of interest. Denote the reference dose as \(d^*\), which can either be a monotherapy dose level (for monotherapy trials) or a dose combination. Now, for each dosing interval \(I_i\), the corresponding reference interval probabilities \(p^*_i\) are computed from the posterior of the BLRM, i.e., assuming that trial \(j\) is the trial of interest, compute for \(i=1,...,4\) $$ p^*_i = Prob(\pi_j(d^*) \in I_i). $$ Using these, define the current dynamic weight \(w^*\) as the interpolation $$ w^* = p^*_1 * W_1 + p^*_2 * W_2 + p^*_3 * W_3 + p^*_4 * W_4. $$ These dynamic weights can now used to define a formal loss function and arrive at a decision by computing the dose with the minimal expected loss based on the dynamic penalty weights for the intervals.

Intuitively, the heuristic assumes that the reference dose is a priori assumed to be a likely candidate for the MTD, and typically among the larger planned dose levels. That is, when the reference dose has a large probability of being an underdose, the dose-toxicity scenario was less toxic than expected a priori, and a more aggressive weight combination can be applied. Conversely, if the reference dose has a large probability of being overly toxic, the planned dose levels (and therefore, the investigational treatment) may likely be more toxic than expected a priori, so that a more conservative configuration of interval penalties (larger penalty for overdosing intervals) is appropriate. Due to this, the weights \(W_1\) should represent the most aggressive configuration, \(W_2\) and \(W_3\) should be moderate to conservative configurations, while \(W_4\) should be the most conservative among the prespecified weights. Refer to the parameter documentation of dynamic.weights for a proposed default configuration that was found to work well in many situations.

References

Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. https://mc-stan.org

Neuenschwander, B., Branson, M., & Gsponer, T. (2008). Critical aspects of the Bayesian approach to phase I cancer trials. Statistics in medicine, 27(13), 2420-2439, doi:10.1002/sim.3230.

Neuenschwander, B., Matano, A., Tang, Z., Roychoudhury, S., Wandel, S., & Bailey, S. (2014). A Bayesian Industry Approach to Phase I Combination Trials in Oncology. In: Zhao. W & Yang, H. (editors). Statistical methods in drug combination studies. Chapman and Hall/CRC, 95-135, doi:10.1201/b17965.

Neuenschwander, B., Roychoudhury, S., & Schmidli, H. (2016). On the use of co-data in clinical trials. Statistics in Biopharmaceutical Research, 8(3), 345-354, doi:10.1080/19466315.2016.1174149.

Babb, J., Rogatko, A., & Zacks, S. (1998). Cancer phase I clinical trials: Efficient dose escalation with overdose control. Statistics in medicine 17(10), 1103-1120.

Zhou, H., Yuan, Y., & Nie, L. (2018). Accuracy, safety, and reliability of novel phase I designs. Clinical Cancer Research, 24(21), 5483-5484 <doi: 10.1158/1078-0432.ccr-18-0168>.

See also

scenario_list_jointBLRM(), sim_jointBLRM(), rstan::sampling(), rstan-package, ggplot2::ggplot, ggplot2-package.

Examples

if (FALSE) {
result <- scenario_jointBLRM(
                   data=list(dose1 = c(1, 2, 4, 6, 8, 0,  0,  0,  1,  2),
                             dose2 = c(0, 0, 0, 0, 0, 10, 20, 30, 10, 10),
                             n.pat = c(3, 3, 3, 3, 3, 3,  6,  9,  3,  3),
                             n.dlt = c(0, 0, 0, 0, 1, 0,  0,  1,  0,  0),
                             trial = c(1, 1, 1, 1, 1, 2,  2,  3,  3,  3)
                             ),
                   trials.of.interest = c(1,       3),
                   types.of.interest =   c("mono1", "combi"),
                   doses.of.interest = rbind(
                     c(1, 2, 4, 6, 8, 12,  rep(c(1, 2, 4, 6, 8, 12),
                                               times=3)),
                     c(0, 0, 0, 0, 0, 0,   rep(c(10, 20, 30),
                                               each = 6 ))),
                   dose.ref1 = 12,
                   dose.ref2 = 30,
                   esc.rule = "dynamic.loss",
                   prior.mu = list(mu_a1 =  c(logit(0.33), 2),
                                   mu_b1 =  c(0,          1),
                                   mu_a2 =  c(logit(0.33), 2),
                                   mu_b2 =  c(0,          1),
                                   mu_eta = c(0,          1.121)),
                   prior.tau = list(tau_a1 =  c(log(0.25),  log(2)/1.96),
                                    tau_b1 =  c(log(0.125), log(2)/1.96),
                                    tau_a2 =  c(log(0.25),  log(2)/1.96),
                                    tau_b2 =  c(log(0.125), log(2)/1.96),
                                    tau_eta = c(log(0.125), log(2)/1.96)),
                   path = getwd(),
                   file.name = NULL,
                   iter=10000,
                   chains=4
                   )
}