blrm_exnex {OncoBayes2}

Bayesian Logistic Regression Model for N-compounds with EXNEX

Description

Bayesian Logistic Regression Model (BLRM) for N compounds using EXchangability and NonEXchangability (EXNEX) modeling.

Usage

blrm_exnex(
  formula,
  data,
  prior_EX_mu_mean_comp,
  prior_EX_mu_sd_comp,
  prior_EX_tau_mean_comp,
  prior_EX_tau_sd_comp,
  prior_EX_corr_eta_comp,
  prior_EX_mu_mean_inter,
  prior_EX_mu_sd_inter,
  prior_EX_tau_mean_inter,
  prior_EX_tau_sd_inter,
  prior_EX_corr_eta_inter,
  prior_is_EXNEX_inter,
  prior_is_EXNEX_comp,
  prior_EX_prob_comp,
  prior_EX_prob_inter,
  prior_NEX_mu_mean_comp,
  prior_NEX_mu_sd_comp,
  prior_NEX_mu_mean_inter,
  prior_NEX_mu_sd_inter,
  prior_tau_dist,
  iter = getOption("OncoBayes2.MC.iter", 2000),
  warmup = getOption("OncoBayes2.MC.warmup", 1000),
  save_warmup = getOption("OncoBayes2.MC.save_warmup", TRUE),
  thin = getOption("OncoBayes2.MC.thin", 1),
  init = getOption("OncoBayes2.MC.init", 0.5),
  chains = getOption("OncoBayes2.MC.chains", 4),
  cores = getOption("mc.cores", 1L),
  control = getOption("OncoBayes2.MC.control", list()),
  backend = getOption("OncoBayes2.MC.backend", "rstan"),
  prior_PD = FALSE,
  verbose = FALSE
)

## S3 method for class 'blrmfit'
print(x, ..., prob = 0.95, digits = 2)

Arguments

formula

the model formula describing the linear predictors of the model. The lhs of the formula is a two-column matrix which are the number of occured events and the number of times no event occured. The rhs of the formula defines the linear predictors for the marginal models for each drug component, then the interaction model and at last the grouping and optional stratum factors of the models. These elements of the formula are separated by a vertical bar. The marginal models must follow a intercept and slope form while the interaction model must not include an interaction term. See the examples below for an example instantiation.

data

optional data frame containing the variables of the model. If not found in data, the variables are taken from environment(formula).

prior_EX_mu_mean_comp, prior_EX_mu_sd_comp

Mean and sd for the prior on the mean parameters \boldsymbol\mu_i = (\mu_{\alpha i}, \mu_{\beta i}) of each component. Two column matrix (intercept, log-slope) with one row per component.

prior_EX_tau_mean_comp, prior_EX_tau_sd_comp

Prior mean and sd for heterogeniety parameter \boldsymbol\tau_{si} = (\tau_{\alpha s i}, \tau_{\beta s i}) of each stratum. If no differential discounting is required (i.e. if there is only one stratum s = 1), then it is a two-column matrix (intercept, log-slope) with one row per component. Otherwise it is a three-dimensional array whose first dimension indexes the strata, second dimension indexes the components, and third dimension of length two for (intercept, log-slope).

prior_EX_corr_eta_comp

Prior LKJ correlation parameter for each component given as numeric vector. If missing, then a 1 is assumed corresponding to a marginal uniform prior of the correlation.

prior_EX_mu_mean_inter, prior_EX_mu_sd_inter

Prior mean and sd for population mean parameters \mu_{\eta k} of each interaction parameter. Vector of length equal to the number of interactions.

prior_EX_tau_mean_inter, prior_EX_tau_sd_inter

Prior mean and sd for heterogeniety parameter \tau_{\eta s k} of each stratum. Matrix with one column per interaction and one row per stratum.

prior_EX_corr_eta_inter

Prior LKJ correlation parameter for interaction given as numeric. If missing, then a 1 is assumed corresponding to a marginal uniform prior of the correlations.

prior_is_EXNEX_inter

Defines if non-exchangability is admitted for a given interaction parameter. Logical vector of length equal to the number of interactions. If missing FALSE is assumed for all interactions.

prior_is_EXNEX_comp

Defines if non-exchangability is admitted for a given component. Logical vector of length equal to the number of components. If missing TRUE is assumed for all components.

prior_EX_prob_comp

Prior probability p_{ij} for exchangability of each component per group. Matrix with one column per component and one row per group. Values must lie in [0-1] range.

prior_EX_prob_inter

Prior probability p_{\eta k j} for exchangability of each interaction per group. Matrix with one column per interaction and one row per group. Values must lie in [0-1] range.

prior_NEX_mu_mean_comp, prior_NEX_mu_sd_comp

Prior mean \boldsymbol m_{ij} and sd \boldsymbol s_{ij} = \mbox{diag}(\boldsymbol S_{ij}) of each component for non-exchangable case. Two column matrix (intercept, log-slope) with one row per component. If missing set to the same prior as given for the EX part. It is required that the specification be the same across groups j.

prior_NEX_mu_mean_inter, prior_NEX_mu_sd_inter

Prior mean m_{\eta k j} and sd s_{\eta k j} for each interaction parameter for non-exchangable case. Vector of length equal to the number of interactions. If missing set to the same prior as given for the EX part.

prior_tau_dist

Defines the distribution used for heterogeniety parameters. Choices are 0=fixed to it's mean, 1=log-normal, 2=truncated normal.

iter

number of iterations (including warmup).

warmup

number of warmup iterations.

save_warmup

save warmup samples (TRUE / FALSE). Only if set to TRUE, then all random variables are saved in the posterior. This substantially increases the storage needs of the posterior.

thin

period of saving samples.

init

positive number to specify uniform range on unconstrained space for random initialization. See stan.

chains

number of Markov chains.

cores

number of cores for parallel sampling of chains.

control

additional sampler parameters for NuTS algorithm.

backend

sets Stan backend to be used. Possible choices are "rstan" (default) or "cmdstanr".

prior_PD

Logical flag (defaults to FALSE) indicating if to sample the prior predictive distribution instead of conditioning on the data.

verbose

Logical flag (defaults to FALSE) controlling if additional output like stan progress is reported.

x

blrmfit object to print

...

not used in this function

prob

central probability mass to report, i.e. the quantiles 0.5-prob/2 and 0.5+prob/2 are displayed. Multiple central widths can be specified.

digits

number of digits to show

Details

blrm_exnex is a flexible function for Bayesian meta-analytic modeling of binomial count data. In particular, it is designed to model counts of the number of observed dose limiting toxicities (DLTs) by dose, for guiding dose-escalation studies in Oncology. To accommodate dose escalation over more than one agent, the dose may consist of combinations of study drugs, with any number of treatment components.

In the simplest case, the aim is to model the probability \pi that a patient experiences a DLT, by complementing the binomial likelihood with a monotone logistic regression

\mbox{logit}\,\pi(d) = \log\,\alpha + \beta \, t(d),

where \beta > 0. Most typically, d represents the dose, and t(d) is an appropriate transformation, such as t(d) = \log (d \big / d^*). A joint prior on \boldsymbol \theta = (\log\,\alpha, \log\,\beta) completes the model and ensures monotonicity \beta > 0.

Many extensions are possible. The function supports general combination regimens, and also provides framework for Bayesian meta-analysis of dose-toxicity data from multiple historical and concurrent sources.

For an example of a single-agent trial refer to example-single-agent.

Value

The function returns a S3 object of type blrmfit.

Methods (by generic)

• print: print function.

Combination of two treatments

For a combination of two treatment components, the basic modeling framework is that the DLT rate \pi(d_1,d_2) is comprised of (1) a "no-interaction" baseline model \tilde \pi(d_1,d_2) driven by the single-agent toxicity of each component, and (2) optional interaction terms \gamma(d_1,d_2) representing synergy or antagonism between the drugs. On the log-odds scale,

\mbox{logit} \,\pi(d_1,d_2) = \mbox{logit} \, \tilde \pi(d_1,d_2) + \eta \, \gamma(d_1,d_2).

The "no interaction" part \tilde \pi(d_1,d_2) represents the probability of a DLT triggered by either treatment component acting independently. That is,

\tilde \pi(d_1,d_2) = 1- (1 - \pi_1(d_1))(1 - \pi_2(d_2)).

In simple terms, P(no DLT for combination) = P(no DLT for drug 1) * P(no DLT from drug 2). To complete this part, the treatment components can then be modeled with monotone logistic regressions as before.

\mbox{logit} \, \pi_i(d_i) = \log\, \alpha_i + \beta_i \, t(d_i),

where t(d_i) is a monotone transformation of the doses, such as t(d_i) = \log (d_i \big / d_i^*).

The inclusion of an interaction term \gamma(d_1,d_2) allows DLT rates above or below the "no-interaction" rate. The magnitude of the interaction term may also be made dependent on the doses (or other covariates) through regression. As an example, one could let

\gamma(d_1, d_2) = \frac{d_1}{d_1^*} \frac{d_1}{d_2^*}.

The specific functional form is specified in the usual notation for a design matrix. The interaction model must respect the constraint that whenever any dose approaches zero, then the interaction term must vanish as well. Therefore, the interaction model must not include an intercept term which would violate this consistency requirement. A dual combination example can be found in example-combo2.

General combinations

The model is extended to general combination treatments consisting of N components by expressing the probability \pi on the logit scale as

\mbox{logit} \, \pi(d_1,\ldots,d_N) = \mbox{logit} \Bigl( 1 - \prod_{i = 1}^N ( 1 - \pi_i(d_i) ) \Bigr) + \sum_{k=1}^K \eta_k \, \gamma_k(d_1,\ldots,d_N),

Multiple drug-drug interactions among the N components are now possible, and are represented through the K interaction terms \gamma_k.

Regression models can be again be specified for each \pi_i and \gamma_k, such as

\mbox{logit}\, \pi_i(d_i) = \log\, \alpha_i + \beta_i \, t(d_i)

Interactions for some subset I(k) \subset \{1,\ldots,N \} of the treatment components can be modeled with regression as well, for example on products of doses,

\gamma_k(d_1,\ldots,d_N) = \prod_{i \in I(k)} \frac{d_i}{d_i^*}.

For example, I(k) = \{1,2,3\} results in the three-way interaction term

\frac{d_1}{d_1^*} \frac{d_2}{d_2^*} \frac{d_3}{d_3^*}

for drugs 1, 2, and 3.

For a triple combination example please refer to example-combo3.

Meta-analytic framework

Information on the toxicity of a drug may be available from multiple studies or sources. Furthermore, one may wish to stratify observations within a single study (for example into groups of patients corresponding to different geographic regions, or multiple dosing dose_info corresponding to different schedules).

blrm_exnex provides tools for robust Bayesian hierarchical modeling to jointly model data from multiple sources. An additional index j=1, \ldots, J on the parameters and observations denotes the J groups. The resulting model allows the DLT rate to differ across the groups. The general N-component model becomes

\mbox{logit} \, \pi_j(d_1,\ldots,d_N) = \mbox{logit} \Bigl( 1 - \prod_{i = 1}^N ( 1 - \pi_{ij}(d_i) ) \Bigr) + \sum_{k=1}^K \eta_{kj} \, \gamma_{k}(d_1,\ldots,d_N),

for groups j = 1,\ldots,J. The component toxicities \pi_{ij} and interaction terms \gamma_{k} are modelled, as before, through regression. For example, \pi_{ij} could be a logistic regression on t(d_i) = \log(d_i/d_i^*) with intercept and log-slope \boldsymbol \theta_{ij}, and \gamma_{k} regressed with coefficient \eta_{kj} on a product \prod_{i\in I(k)} (d_i/d_i^*) for some subset I(k) of components.

Thus, for j=1,\ldots,J, we now have group-specific parameters \boldsymbol\theta_{ij} = (\log\, \alpha_{ij}, \log\, \beta_{ij}) and \boldsymbol\nu_{j} = (\eta_{1j}, \ldots, \eta_{Kj}) for each component i=1,\ldots,N and interaction k=1,\ldots,K.

The structure of the prior on (\boldsymbol\theta_{i1},\ldots,\boldsymbol\theta_{iJ}) and (\boldsymbol\nu_{1}, \ldots, \boldsymbol\nu_{J}) determines how much information will be shared across groups j. Several modeling choices are available in the function.

• EX (Full exchangeability): One can assume the parameters are conditionally exchangeable given hyperparameters

  \boldsymbol \theta_{ij} \sim \mbox{N}\bigl( \boldsymbol \mu_{\boldsymbol \theta i}, \boldsymbol \Sigma_{\boldsymbol \theta i} \bigr),

  independently across groups j = 1,\ldots, J and treatment components i=1,\ldots,N. The covariance matrix \boldsymbol \Sigma_{\boldsymbol \theta i} captures the patterns of cross-group heterogeneity, and is parametrized with standard deviations \boldsymbol \tau_{\boldsymbol\theta i} = (\tau_{\alpha i}, \tau_{\beta i}) and the correlation \rho_i. Similarly for the interactions, the fully-exchangeable model is

  \boldsymbol \nu_{j} \sim \mbox{N}\bigl( \boldsymbol \mu_{\boldsymbol \nu}, \boldsymbol \Sigma_{\boldsymbol \nu} \bigr)

  for groups j = 1,\ldots, J and interactions k=1,\ldots,K, and the prior on the covariance matrix \boldsymbol \Sigma_{\boldsymbol \nu} captures the amount of heterogeneity expected in the interaction terms a-priori. The covariance is again parametrized with standard deviations (\tau_{\eta 1}, \ldots, \tau_{\eta K}) and its correlation matrix.

• Differential discounting: For one or more of the groups j=1,\ldots,J, larger deviations of \boldsymbol\theta_{ij} may be expected from the mean \boldsymbol\mu_i, or of the interactions \eta_{kj} from the mean \mu_{\eta,k}. Such differential heterogeneity can be modeled by mapping the groups j = 1,\ldots,J to strata through s_j \in \{1,\ldots,S\}, and modifying the model specification to

  \boldsymbol \theta_{ij} \sim \mbox{N}\bigl( \boldsymbol \mu_{\boldsymbol \theta i}, \boldsymbol \Sigma_{\boldsymbol \theta ij} \bigr),

  where

  \boldsymbol \Sigma_{\boldsymbol \theta ij} = \left( \begin{array}{cc} \tau^2_{\alpha s_j i} & \rho_i \tau_{\alpha s_j i} \tau_{\beta s_j i}\\ \rho_i \tau_{\alpha s_j i} \tau_{\beta s_j i} & \tau^2_{\beta s_j i} \end{array} \right).

  For the interactions, the model becomes

  \boldsymbol \nu_{j} \sim \mbox{N}\bigl( \boldsymbol \mu_{\boldsymbol \nu}, \boldsymbol \Sigma_{\boldsymbol \nu j} \bigr),

  where the covariance matrix \boldsymbol \Sigma_{\boldsymbol \nu j} is modelled as stratum specific standard deviations (\tau_{\eta 1 s_j}, \ldots, \tau_{\eta K s_j}) and a stratum independent correlation matrix. Each stratum s=1,\ldots,S then corresponds to its own set of standard deviations \tau leading to different discounting per stratum. Independent priors are specified for the component parameters \tau_{\alpha s i} and \tau_{\beta s i} and for the interaction parameters \tau_{\eta s k} for each stratum s=1,\ldots,S. Inference for strata s where the prior is centered on larger values of the \tau parameters will exhibit less shrinkage towards the the means, \boldsymbol\mu_{\boldsymbol \theta i} and \boldsymbol \mu_{\boldsymbol \nu} respectively.

• EXNEX (Partial exchangeability): Another mechansim for increasing robustness is to introduce mixture priors for the group-specific parameters, where one mixture component is shared across groups, and the other is group-specific. The result, known as an EXchangeable-NonEXchangeable (EXNEX) type prior, has a form

  \boldsymbol \theta_{ij} \sim p_{\boldsymbol \theta ij}\, \mbox{N}\bigl( \boldsymbol \mu_{\boldsymbol \theta i}, \boldsymbol \Sigma_{\boldsymbol \theta i} \bigr) +(1-p_{\boldsymbol \theta ij})\, \mbox{N}\bigl(\boldsymbol m_{\boldsymbol \theta ij}, \boldsymbol S_{\boldsymbol \theta ij}\bigr)

  when applied to the treatment-component parameters, and

  \boldsymbol \nu_{kj} \sim p_{\boldsymbol \nu_{kj}} \,\mbox{N}\bigl(\mu_{\boldsymbol \nu}, \boldsymbol \Sigma_{\boldsymbol \nu}\bigr)_k + (1-p_{\boldsymbol \nu_{kj}})\, \mbox{N}(m_{\boldsymbol \nu_{kj}}, s^2_{\boldsymbol \nu_{kj}})

  when applied to the interaction parameters. The exchangeability weights p_{\boldsymbol \theta ij} and p_{\boldsymbol \nu_{kj}} are fixed constants in the interval [0,1] that control the degree to which inference for group j is informed by the exchangeable mixture components. Larger values for the weights correspond to greater exchange of information, while smaller values increase robustness in case of outlying observations in individual groups j.

References

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

Neuenschwander, B., Wandel, S., Roychoudhury, S., & Bailey, S. (2016). Robust exchangeability designs for early phase clinical trials with multiple strata. Pharmaceutical statistics, 15(2), 123-134.

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.

Neuenschwander, B., Matano, A., Tang, Z., Roychoudhury, S., Wandel, S. Bailey, Stuart. (2014). A Bayesian Industry Approach to Phase I Combination Trials in Oncology. In Statistical methods in drug combination studies (Vol. 69). CRC Press.

Examples

## Setting up dummy sampling for fast execution of example
## Please use 4 chains and 100x more warmup & iter in practice
.user_mc_options <- options(OncoBayes2.MC.warmup=10, OncoBayes2.MC.iter=20, OncoBayes2.MC.chains=1,
                            OncoBayes2.MC.save_warmup=FALSE)

# fit an example model. See documentation for "combo3" example
example_model("combo3")

# print a summary of the prior
prior_summary(blrmfit, digits = 3)

# print a summary of the posterior (model parameters)
print(blrmfit)

# summary of posterior for DLT rate by dose for observed covariate levels
summ <- summary(blrmfit, interval_prob = c(0, 0.16, 0.33, 1))
print(cbind(hist_combo3, summ))

# summary of posterior for DLT rate by dose for new set of covariate levels
newdata <- expand.grid(
  stratum_id = "BID", group_id = "Combo",
  drug_A = 400, drug_B = 800, drug_C = c(320, 400, 600, 800),
  stringsAsFactors = FALSE
)
summ_pred <- summary(blrmfit, newdata = newdata, interval_prob = c(0, 0.16, 0.33, 1))
print(cbind(newdata, summ_pred))

# update the model after observing additional data
newdata$num_patients <- rep(3, nrow(newdata))
newdata$num_toxicities <- c(0, 1, 2, 2)
library(dplyr)
blrmfit_new <- update(blrmfit,
                      data = rbind(hist_combo3, newdata) %>%
                               arrange(stratum_id, group_id))

# updated posterior summary
summ_upd <- summary(blrmfit_new, newdata = newdata, interval_prob = c(0, 0.16, 0.33, 1))
print(cbind(newdata, summ_upd))
## Recover user set sampling defaults
options(.user_mc_options)