Calculate MCMC Diagnostics for Parameters

Description

Calculate MCMC diagnostics for individual parameters.

Usage

diagnostics(x)

Arguments

Value

A tibble listing the Gelman point estimates and upper bounds (obtained from coda::gelman.diag) and effective sample size (obtained from coda::effectiveSize) for each parameter within each model.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

# for all models
diagnostics(output)

# for a single model
diagnostics(output$mod_quad)

Generate Data from Dose Response Models

Description

See the model definitions below for specifics for each model.

Usage

dreamer_data_linear(
  n_cohorts,
  doses,
  b1,
  b2,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_linear_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_quad(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_quad_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_loglinear(
  n_cohorts,
  doses,
  b1,
  b2,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_loglinear_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_logquad(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_logquad_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_emax(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  b4,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_emax_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  b4,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_exp(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_exp_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_beta(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  b4,
  scale,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_beta_binary(
  n_cohorts,
  doses,
  b1,
  b2,
  b3,
  b4,
  scale,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_independent(
  n_cohorts,
  doses,
  b1,
  sigma,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

dreamer_data_independent_binary(
  n_cohorts,
  doses,
  b1,
  link,
  times,
  t_max,
  longitudinal = NULL,
  ...
)

Arguments

Value

A dataframe of random subjects from the specified model and parameters.

Functions

• dreamer_data_linear(): generate data from linear dose response.

• dreamer_data_linear_binary(): generate data from linear binary dose response.

• dreamer_data_quad(): generate data from quadratic dose response.

• dreamer_data_quad_binary(): generate data from quadratic binary dose response.

• dreamer_data_loglinear(): generate data from log-linear dose response.

• dreamer_data_loglinear_binary(): generate data from binary log-linear dose response.

• dreamer_data_logquad(): generate data from log-quadratic dose response.

• dreamer_data_logquad_binary(): generate data from log-quadratic binary dose response.

• dreamer_data_emax(): generate data from EMAX dose response.

• dreamer_data_emax_binary(): generate data from EMAX binary dose response.

• dreamer_data_exp(): generate data from exponential dose response.

• dreamer_data_exp_binary(): generate data from exponential binary dose response.

• dreamer_data_beta(): generate data from Beta dose response.

• dreamer_data_beta_binary(): generate data from binary Beta dose response.

• dreamer_data_independent(): generate data from an independent dose response.

• dreamer_data_independent_binary(): generate data from an independent dose response.

Linear

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * d

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Quadratic

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * d + b_3 * d^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Log-linear

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * log(d + 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Log-quadratic

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * log(d + 1) + b_3 * log(d + 1)^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

EMAX

y \sim N(f(d), \sigma^2)

f(d) = b_1 + (b_2 - b_1) * d ^ b_4 / (exp(b_3 * b_4) + d ^ b_4)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

b_4 \sim N(mu_b4, sigma_b4 ^ 2), (Truncated above 0)

1 / \sigma^2 \sim Gamma(shape, rate)

Here, b_1 is the placebo effect (dose = 0), b_2 is the maximum treatment effect, b_3 is the log(ED50), and b_4 is the hill or rate parameter.

Exponential

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * (1 - exp(- b_3 * d))

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), (truncated to be positive)

1 / \sigma^2 \sim Gamma(shape, rate)

Linear Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * d

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

Quadratic Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * d + b_3 * d^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

Log-linear Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * log(d + 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

Log-quadratic Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * log(d + 1) + b_3 * log(d + 1)^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

EMAX Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + (b_2 - b_1) * d ^ b_4 / (exp(b_3 * b_4) + d ^ b_4)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

b_4 \sim N(mu_b4, sigma_b4 ^ 2), (Truncated above 0)

Here, on the link(f(d)) scale, b_1 is the placebo effect (dose = 0), b_2 is the maximum treatment effect, b_3 is the log(ED50), and b_4 is the hill or rate parameter.

Exponential Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * (exp(b_3 * d) - 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), (Truncated below 0)

Independent

y \sim N(f(d), \sigma^2)

f(d) = b_{1d}

b_{1d} \sim N(mu_b1[d], sigma_b1[d] ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Independent Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_{1d}

b_{1d} \sim N(mu_b1[d], sigma_b1[d]) ^ 2

Longitudinal Linear

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + (t / t_max) * f(d)

a \sim N(mu_a, sigma_a)

Longitudinal ITP

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * ((1 - exp(- c1 * t))/(1 - exp(- c1 * t_max)))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

Longitudinal IDP

Increasing-Decreasing-Plateau (IDP).

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * (((1 - exp(- c1 * t))/(1 - exp(- c1 * d1))) * I(t < d1) + (1 - gam * ((1 - exp(- c2 * (t - d1))) / (1 - exp(- c2 * (d2 - d1))))) * I(d1 <= t <= d2) + (1 - gam) * I(t > d2))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

c2 \sim Uniform(a_c2, b_c2)

d1 \sim Uniform(0, t_max)

d2 \sim Uniform(d1, t_max)

gam \sim Uniform(0, 1)

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

head(data)

plot(data$dose, data$response)
abline(a = 1, b = 3)

# longitudinal data
set.seed(889)
data_long <- dreamer_data_linear(
  n_cohorts = c(10, 10, 10, 10), # number of subjects in each cohort
  doses = c(.25, .5, .75, 1.5), # dose administered to each cohort
  b1 = 0, # intercept
  b2 = 2, # slope
  sigma = .5, # standard deviation,
  longitudinal = "itp",
  times = c(0, 12, 24, 52),
  t_max = 52, # maximum time
  a = .5,
  c1 = .1
)

## Not run: 
  ggplot(data_long, aes(time, response, group = dose, color = factor(dose))) +
    geom_point()

## End(Not run)

Bayesian Model Averaging of Dose Response Models

Description

This function performs Bayesian model averaging with a selection of dose response models. See model for all possible models.

Usage

dreamer_mcmc(
  data,
  ...,
  n_adapt = 1000,
  n_burn = 1000,
  n_iter = 10000,
  n_chains = 4,
  silent = FALSE,
  convergence_warn = TRUE
)

Arguments

Details

The Bayesian model averaging approach uses data, multiple models, priors on each model's parameters, and a prior weight for each model. Using these inputs, each model is fit independently, and the output from the models is used to calculate posterior weights for each model. See Gould (2018) for details.

Value

A named list with S3 class "dreamer_bma" and "dreamer". The list contains the following fields:

• doses: a vector of the unique ordered doses in the data.

• times: a vector of the unique ordered times in the data.

• w_prior: a named vector with the prior probabilities of each model.

• w_post: a named vector with the posterior probabilities of each model.

• The individual MCMC fits for each model.

References

Gould, A. L. (2019). BMA-Mod: A Bayesian model averaging strategy for determining dose-response relationships in the presence of model uncertainty. Biometrical Journal, 61(5), 1141-1159.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
  data = data,
  n_adapt = 1e3,
  n_burn = 1e2,
  n_iter = 1e3,
  n_chains = 2,
  silent = TRUE,
  mod_linear = model_linear(
    mu_b1 = 0,
    sigma_b1 = 1,
    mu_b2 = 0,
    sigma_b2 = 1,
    shape = 1,
    rate = .001,
    w_prior = 1 / 2
  ),
  mod_quad = model_quad(
    mu_b1 = 0,
    sigma_b1 = 1,
    mu_b2 = 0,
    sigma_b2 = 1,
    mu_b3 = 0,
    sigma_b3 = 1,
    shape = 1,
    rate = .001,
    w_prior = 1 / 2
  )
)
# posterior weights
output$w_post
# plot posterior dose response
plot(output)

# LONGITUDINAL
library(ggplot2)
set.seed(889)
data_long <- dreamer_data_linear(
  n_cohorts = c(10, 10, 10, 10), # number of subjects in each cohort
  doses = c(.25, .5, .75, 1.5), # dose administered to each cohort
  b1 = 0, # intercept
  b2 = 2, # slope
  sigma = .5, # standard deviation,
  longitudinal = "itp",
  times = c(0, 12, 24, 52),
  t_max = 52, # maximum time
  a = .5,
  c1 = .1
)

## Not run: 
ggplot(data_long, aes(time, response, group = dose, color = factor(dose))) +
   geom_point()

## End(Not run)

output_long <- dreamer_mcmc(
   data = data_long,
   n_adapt = 1e3,
   n_burn = 1e2,
   n_iter = 1e3,
   n_chains = 2,
   silent = TRUE, # make rjags be quiet,
   mod_linear = model_linear(
      mu_b1 = 0,
      sigma_b1 = 1,
      mu_b2 = 0,
      sigma_b2 = 1,
      shape = 1,
      rate = .001,
      w_prior = 1 / 2, # prior probability of the model
      longitudinal = model_longitudinal_itp(
         mu_a = 0,
         sigma_a = 1,
         a_c1 = 0,
         b_c1 = 1,
         t_max = 52
      )
   ),
   mod_quad = model_quad(
      mu_b1 = 0,
      sigma_b1 = 1,
      mu_b2 = 0,
      sigma_b2 = 1,
      mu_b3 = 0,
      sigma_b3 = 1,
      shape = 1,
      rate = .001,
      w_prior = 1 / 2,
      longitudinal = model_longitudinal_linear(
         mu_a = 0,
         sigma_a = 1,
         t_max = 52
      )
   )
)

## Not run: 
# plot longitudinal dose-response profile
plot(output_long, data = data_long)
plot(output_long$mod_quad, data = data_long) # single model

# plot dose response at final timepoint
plot(output_long, data = data_long, times = 52)
plot(output_long$mod_quad, data = data_long, times = 52) # single model

## End(Not run)

Plot Prior

Description

Plot the prior over the dose range. This is intended to help the user choose appropriate priors.

Usage

dreamer_plot_prior(
  n_samples = 10000,
  probs = c(0.025, 0.975),
  doses,
  n_chains = 1,
  ...,
  times = NULL,
  plot_draws = FALSE,
  alpha = 0.2
)

Arguments

Value

The ggplot object.

Examples

# Plot prior for one model
set.seed(8111)
dreamer_plot_prior(
 doses = c(0, 2.5, 5),
 mod_quad_binary = model_quad_binary(
   mu_b1 = -.5,
   sigma_b1 = .2,
   mu_b2 = -.5,
   sigma_b2 = .2,
   mu_b3 = .5,
   sigma_b3 = .1,
   link = "logit",
   w_prior = 1
 )
)

# plot individual draws
dreamer_plot_prior(
 doses = seq(from = 0, to = 5, length.out = 50),
 n_samples = 100,
 plot_draws = TRUE,
 mod_quad_binary = model_quad_binary(
   mu_b1 = -.5,
   sigma_b1 = .2,
   mu_b2 = -.5,
   sigma_b2 = .2,
   mu_b3 = .5,
   sigma_b3 = .1,
   link = "logit",
   w_prior = 1
 )
)

# plot prior from mixture of models
dreamer_plot_prior(
 doses = c(0, 2.5, 5),
 mod_linear_binary = model_linear_binary(
   mu_b1 = -1,
   sigma_b1 = .1,
   mu_b2 = 1,
   sigma_b2 = .1,
   link = "logit",
   w_prior = .75
 ),
 mod_quad_binary = model_quad_binary(
   mu_b1 = -.5,
   sigma_b1 = .2,
   mu_b2 = -.5,
   sigma_b2 = .2,
   mu_b3 = .5,
   sigma_b3 = .1,
   link = "logit",
   w_prior = .25
 )
)

Posterior Plot of Bayesian Model Averaging

Description

Plots the posterior mean and quantiles over the dose range and plots error bars at the observed doses. If the data argument is specified, the observed means at each dose are also plotted.

Usage

## S3 method for class 'dreamer_mcmc'
plot(
  x,
  doses = attr(x, "doses"),
  times = attr(x, "times"),
  probs = c(0.025, 0.975),
  data = NULL,
  n_smooth = 50,
  predictive = 0,
  width = bar_width(doses),
  reference_dose = NULL,
  ...
)

Arguments

Value

Returns the ggplot object.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

plot(output)

## Not run: 
# with data
plot(output, data = data)

# predictive distribution
plot(output, data = data, predictive = 1)

# single model
plot(output$mod_linear)

## End(Not run)

Model Creation

Description

Functions which set the hyperparameters, seeds, and prior weight for each model to be used in Bayesian model averaging via dreamer_mcmc().

See each function's section below for the model's details. In the following, y denotes the response variable and d represents the dose.

For the longitudinal specifications, see documentation on model_longitudinal.

Usage

model_linear(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_quad(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_loglinear(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_logquad(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_emax(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  mu_b4,
  sigma_b4,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_exp(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  shape,
  rate,
  w_prior = 1,
  longitudinal = NULL
)

model_beta(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  mu_b4,
  sigma_b4,
  shape,
  rate,
  scale = NULL,
  w_prior = 1,
  longitudinal = NULL
)

model_independent(
  mu_b1,
  sigma_b1,
  shape,
  rate,
  doses = NULL,
  w_prior = 1,
  longitudinal = NULL
)

model_linear_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_quad_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_loglinear_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_logquad_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_emax_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  mu_b4,
  sigma_b4,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_exp_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_beta_binary(
  mu_b1,
  sigma_b1,
  mu_b2,
  sigma_b2,
  mu_b3,
  sigma_b3,
  mu_b4,
  sigma_b4,
  scale = NULL,
  link,
  w_prior = 1,
  longitudinal = NULL
)

model_independent_binary(
  mu_b1,
  sigma_b1,
  doses = NULL,
  link,
  w_prior = 1,
  longitudinal = NULL
)

Arguments

Value

A named list of the arguments in the function call. The list has S3 classes assigned which are used internally within dreamer_mcmc().

Linear

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * d

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Quadratic

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * d + b_3 * d^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Log-linear

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * log(d + 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Log-quadratic

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * log(d + 1) + b_3 * log(d + 1)^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

EMAX

y \sim N(f(d), \sigma^2)

f(d) = b_1 + (b_2 - b_1) * d ^ b_4 / (exp(b_3 * b_4) + d ^ b_4)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

b_4 \sim N(mu_b4, sigma_b4 ^ 2), (Truncated above 0)

1 / \sigma^2 \sim Gamma(shape, rate)

Here, b_1 is the placebo effect (dose = 0), b_2 is the maximum treatment effect, b_3 is the log(ED50), and b_4 is the hill or rate parameter.

Exponential

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * (1 - exp(- b_3 * d))

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), (truncated to be positive)

1 / \sigma^2 \sim Gamma(shape, rate)

Beta

y \sim N(f(d), \sigma^2)

f(d) = b_1 + b_2 * ((b3 + b4) ^ (b3 + b4)) / (b3 ^ b3 * b4 ^ b4) * (d / scale) ^ b3 * (1 - d / scale) ^ b4

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), Truncated above 0

b_4 \sim N(mu_b4, sigma_b4 ^ 2), Truncated above 0

1 / \sigma^2 \sim Gamma(shape, rate)

Note that scale is a hyperparameter specified by the user.

Independent

y \sim N(f(d), \sigma^2)

f(d) = b_{1d}

b_{1d} \sim N(mu_b1[d], sigma_b1[d] ^ 2)

1 / \sigma^2 \sim Gamma(shape, rate)

Independent Details

The independent model models the effect of each dose independently. Vectors can be supplied to mu_b1 and sigma_b1 to set a different prior for each dose in the order the doses are supplied to doses. If scalars are supplied to mu_b1 and sigma_b1, then the same prior is used for each dose, and the doses argument is not needed.

Linear Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * d

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

Quadratic Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * d + b_3 * d^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

Log-linear Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * log(d + 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

Log-quadratic Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * log(d + 1) + b_3 * log(d + 1)^2

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

EMAX Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + (b_2 - b_1) * d ^ b_4 / (exp(b_3 * b_4) + d ^ b_4)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2)

b_4 \sim N(mu_b4, sigma_b4 ^ 2), (Truncated above 0)

Here, on the link(f(d)) scale, b_1 is the placebo effect (dose = 0), b_2 is the maximum treatment effect, b_3 is the log(ED50), and b_4 is the hill or rate parameter.

Exponential Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_1 + b_2 * (exp(b_3 * d) - 1)

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), (Truncated below 0)

Beta Binary

y \sim Binomial(n, f(d)

link(f(d)) = b_1 + b_2 * ((b3 + b4) ^ (b3 + b4)) / (b3 ^ b3 * b4 ^ b4) * (d / scale) ^ b3 * (1 - d / scale) ^ b4

b_1 \sim N(mu_b1, sigma_b1 ^ 2)

b_2 \sim N(mu_b2, sigma_b2 ^ 2)

b_3 \sim N(mu_b3, sigma_b3 ^ 2), Truncated above 0

b_4 \sim N(mu_b4, sigma_b4 ^ 2), Truncated above 0

Note that scale is a hyperparameter specified by the user.

Independent Binary

y \sim Binomial(n, f(d))

link(f(d)) = b_{1d}

b_{1d} \sim N(mu_b1[d], sigma_b1[d]) ^ 2

Independent Binary Details

The independent model models the effect of each dose independently. Vectors can be supplied to mu_b1 and sigma_b1 to set a different prior for each dose in the order the doses are supplied to doses. If scalars are supplied to mu_b1 and sigma_b1, then the same prior is used for each dose, and the doses argument is not needed.

Longitudinal Linear

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + (t / t_max) * f(d)

a \sim N(mu_a, sigma_a)

Longitudinal ITP

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * ((1 - exp(- c1 * t))/(1 - exp(- c1 * t_max)))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

Longitudinal IDP

Increasing-Decreasing-Plateau (IDP).

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * (((1 - exp(- c1 * t))/(1 - exp(- c1 * d1))) * I(t < d1) + (1 - gam * ((1 - exp(- c2 * (t - d1))) / (1 - exp(- c2 * (d2 - d1))))) * I(d1 <= t <= d2) + (1 - gam) * I(t > d2))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

c2 \sim Uniform(a_c2, b_c2)

d1 \sim Uniform(0, t_max)

d2 \sim Uniform(d1, t_max)

gam \sim Uniform(0, 1)

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
  data = data,
  n_adapt = 1e3,
  n_burn = 1e2,
  n_iter = 1e3,
  n_chains = 2,
  silent = TRUE,
  mod_linear = model_linear(
    mu_b1 = 0,
    sigma_b1 = 1,
    mu_b2 = 0,
    sigma_b2 = 1,
    shape = 1,
    rate = .001,
    w_prior = 1 / 2
  ),
  mod_quad = model_quad(
    mu_b1 = 0,
    sigma_b1 = 1,
    mu_b2 = 0,
    sigma_b2 = 1,
    mu_b3 = 0,
    sigma_b3 = 1,
    shape = 1,
    rate = .001,
    w_prior = 1 / 2
  )
)
# posterior weights
output$w_post
# plot posterior dose response
plot(output)

# LONGITUDINAL
library(ggplot2)
set.seed(889)
data_long <- dreamer_data_linear(
  n_cohorts = c(10, 10, 10, 10), # number of subjects in each cohort
  doses = c(.25, .5, .75, 1.5), # dose administered to each cohort
  b1 = 0, # intercept
  b2 = 2, # slope
  sigma = .5, # standard deviation,
  longitudinal = "itp",
  times = c(0, 12, 24, 52),
  t_max = 52, # maximum time
  a = .5,
  c1 = .1
)

## Not run: 
ggplot(data_long, aes(time, response, group = dose, color = factor(dose))) +
   geom_point()

## End(Not run)

output_long <- dreamer_mcmc(
   data = data_long,
   n_adapt = 1e3,
   n_burn = 1e2,
   n_iter = 1e3,
   n_chains = 2,
   silent = TRUE, # make rjags be quiet,
   mod_linear = model_linear(
      mu_b1 = 0,
      sigma_b1 = 1,
      mu_b2 = 0,
      sigma_b2 = 1,
      shape = 1,
      rate = .001,
      w_prior = 1 / 2, # prior probability of the model
      longitudinal = model_longitudinal_itp(
         mu_a = 0,
         sigma_a = 1,
         a_c1 = 0,
         b_c1 = 1,
         t_max = 52
      )
   ),
   mod_quad = model_quad(
      mu_b1 = 0,
      sigma_b1 = 1,
      mu_b2 = 0,
      sigma_b2 = 1,
      mu_b3 = 0,
      sigma_b3 = 1,
      shape = 1,
      rate = .001,
      w_prior = 1 / 2,
      longitudinal = model_longitudinal_linear(
         mu_a = 0,
         sigma_a = 1,
         t_max = 52
      )
   )
)

## Not run: 
# plot longitudinal dose-response profile
plot(output_long, data = data_long)
plot(output_long$mod_quad, data = data_long) # single model

# plot dose response at final timepoint
plot(output_long, data = data_long, times = 52)
plot(output_long$mod_quad, data = data_long, times = 52) # single model

## End(Not run)

Model Creation: Longitudinal Models

Description

Assign hyperparameters and other values for longitudinal modeling. The output of this function is intended to be used as the input to the longitudinal argument of the dose response model functions, e.g., model_linear.

Usage

model_longitudinal_linear(mu_a, sigma_a, t_max)

model_longitudinal_itp(mu_a, sigma_a, a_c1 = 0, b_c1 = 1, t_max)

model_longitudinal_idp(
  mu_a,
  sigma_a,
  a_c1 = 0,
  b_c1 = 1,
  a_c2 = -1,
  b_c2 = 0,
  t_max
)

Arguments

Value

A named list of the arguments in the function call. The list has S3 classes assigned which are used internally within dreamer_mcmc().

Longitudinal Linear

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + (t / t_max) * f(d)

a \sim N(mu_a, sigma_a)

Longitudinal ITP

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * ((1 - exp(- c1 * t))/(1 - exp(- c1 * t_max)))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

Longitudinal IDP

Increasing-Decreasing-Plateau (IDP).

Let f(d) be a dose response model. The expected value of the response, y, is:

E(y) = g(d, t)

g(d, t) = a + f(d) * (((1 - exp(- c1 * t))/(1 - exp(- c1 * d1))) * I(t < d1) + (1 - gam * ((1 - exp(- c2 * (t - d1))) / (1 - exp(- c2 * (d2 - d1))))) * I(d1 <= t <= d2) + (1 - gam) * I(t > d2))

a \sim N(mu_a, sigma_a)

c1 \sim Uniform(a_c1, b_c1)

c2 \sim Uniform(a_c2, b_c2)

d1 \sim Uniform(0, t_max)

d2 \sim Uniform(d1, t_max)

gam \sim Uniform(0, 1)

Compare Posterior Fits

Description

Compare Posterior Fits

Usage

plot_comparison(..., doses, times, probs, data, n_smooth, width)

## Default S3 method:
plot_comparison(
  ...,
  doses = attr(list(...)[[1]], "doses"),
  times = NULL,
  probs = c(0.025, 0.975),
  data = NULL,
  n_smooth = 50,
  width = bar_width(doses)
)

## S3 method for class 'dreamer_bma'
plot_comparison(
  ...,
  doses = x$doses,
  times = NULL,
  probs = c(0.025, 0.975),
  data = NULL,
  n_smooth = 50,
  width = bar_width(doses)
)

Arguments

Details

If a Bayesian model averaging object is supplied first, all individual fits and the Bayesian model averaging fit will be plotted, with the model averaging fit in black (other model colors specified in the legend). Otherwise, named arguments must be supplied for each model, and only the models provided will be plotted.

Value

a ggplot object.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e3,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

plot_comparison(output)

# compare individual models
plot_comparison(linear = output$mod_linear, quad = output$mod_quad)

Traceplots

Description

Produces traceplots for each parameter for each model.

Usage

plot_trace(x)

Arguments

Value

No return value, called to create plots.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

# all parameters from all models
plot_trace(output)

# from a single model
plot_trace(output$mod_linear)

Posterior Distribution of Minimum X% Effective Dose

Description

Posterior Distribution of Minimum X% Effective Dose

Usage

post_medx(
  x,
  ed,
  probs,
  time,
  lower,
  upper,
  greater,
  small_bound,
  return_samples,
  ...
)

## S3 method for class 'dreamer_bma'
post_medx(
  x,
  ed,
  probs = c(0.025, 0.975),
  time = NULL,
  lower = min(x$doses),
  upper = max(x$doses),
  greater = TRUE,
  small_bound = 0,
  return_samples = FALSE,
  ...
)

## S3 method for class 'dreamer_mcmc'
post_medx(
  x,
  ed,
  probs = c(0.025, 0.975),
  time = NULL,
  lower = min(attr(x, "doses")),
  upper = max(attr(x, "doses")),
  greater = TRUE,
  small_bound = 0,
  return_samples = FALSE,
  index = 1:(nrow(x[[1]]) * length(x)),
  ...
)

Arguments

Details

The minimum X% effective dose is the dose that has X% of the largest effect for doses between lower and upper. When greater is TRUE, larger positive responses are considered more effective and vice versa. The X% response is calculated as small_bound + ed / 100 * (max_effect - small_bound) where "max_effect" is the maximum response for doses between lower and upper. The X% effective dose is the smallest dose which has X% response within the dose range. It is possible that for some MCMC samples, an X% effective dose may not exist, so probabilities are not guaranteed to sum to one.

Value

Posterior quantities and samples (if applicable), generally in the form of a list. The pr_edx_exists column gives the posterior probability that an EDX% effect exists.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e3,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

post_medx(output, ed = c(50, 90))

# from a single model
post_medx(output$mod_linear, ed = c(50, 90))

Calculate Posterior of a Dose's Percentage Effect

Description

Given a dose, the "percentage effect" is defined as (effect of the given dose - small_bound) / (maximum effect in dose range - small_bound). This function returns the posterior statistics and/or samples of this effect.

Usage

post_perc_effect(
  x,
  dose,
  probs,
  time,
  lower,
  upper,
  greater,
  small_bound,
  index,
  return_samples
)

## S3 method for class 'dreamer_bma'
post_perc_effect(
  x,
  dose,
  probs = c(0.025, 0.975),
  time = NULL,
  lower = min(x$doses),
  upper = max(x$doses),
  greater = TRUE,
  small_bound = 0,
  index = NA,
  return_samples = FALSE
)

## S3 method for class 'dreamer_mcmc'
post_perc_effect(
  x,
  dose,
  probs = c(0.025, 0.975),
  time = NULL,
  lower = min(attr(x, "doses")),
  upper = max(attr(x, "doses")),
  greater = TRUE,
  small_bound = 0,
  index = 1:(nrow(x[[1]]) * length(x)),
  return_samples = FALSE
)

Arguments

Value

A named list with the following components:

• stats: a tibble listing the dose, time (where relevant), probability a percentage effect exists, the average percentage effect, and the specified quantiles of the percentage effect.

• samps: a tibble with the posterior samples for each dose/time combination.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e3,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

post_perc_effect(output, dose = c(3, 5))

# from a single model
post_perc_effect(output$mod_linear, dose = c(3, 5))

Posterior Quantities from Bayesian Model Averaging

Description

Calculate posterior mean (and quantiles for specific doses for each MCMC iteration of the model.

Usage

posterior(
  x,
  doses,
  times,
  probs,
  reference_dose,
  predictive,
  return_samples,
  iter,
  return_stats
)

## S3 method for class 'dreamer_mcmc'
posterior(
  x,
  doses = attr(x, "doses"),
  times = attr(x, "times"),
  probs = c(0.025, 0.975),
  reference_dose = NULL,
  predictive = 0,
  return_samples = FALSE,
  iter = NULL,
  return_stats = TRUE
)

## S3 method for class 'dreamer_mcmc_independent'
posterior(
  x,
  doses = attr(x, "doses"),
  times = attr(x, "times"),
  probs = c(0.025, 0.975),
  reference_dose = NULL,
  predictive = 0,
  return_samples = FALSE,
  iter = NULL,
  return_stats = TRUE
)

## S3 method for class 'dreamer_bma'
posterior(
  x,
  doses = x$doses,
  times = x$times,
  probs = c(0.025, 0.975),
  reference_dose = NULL,
  predictive = 0,
  return_samples = FALSE,
  iter = NULL,
  return_stats = TRUE
)

Arguments

Value

A named list with the following elements:

• stats: a tibble the dose, posterior mean, and posterior quantiles.

• samps: the weighted posterior samples. Only returned if return_samples = TRUE.

Methods (by class)

• posterior(dreamer_mcmc): posterior summary for linear model.

• posterior(dreamer_mcmc_independent): posterior summary for independent model.

• posterior(dreamer_bma): posterior summary for Bayesian model averaging fit.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

posterior(output)

# return posterior samples of the mean
post <- posterior(output, return_samples = TRUE)
head(post$samps)

# from a single model
posterior(output$mod_quad)

# posterior of difference of doses
posterior(output, reference_dose = 0)

Calculate Probability of Meeting Effect of Interest (EOI)

Description

Calculate Pr(effect_dose - effect_reference_dose > EOI | data) or Pr(effect_dose > EOI | data).

Usage

pr_eoi(x, eoi, dose, reference_dose = NULL, time = NULL)

Arguments

Value

A tibble listing the doses, times, and Pr(effect_dose - effect_reference_dose > eoi) if reference_dose is specified; otherwise, Pr(effect_dose > eoi).

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

pr_eoi(output, dose = 3, eoi = 10)

# difference of two doses
pr_eoi(output, dose = 3, eoi = 10, reference_dose = 0)

# single model
pr_eoi(output$mod_linear, dose = 3, eoi = 10)

Pr(minimum efficacious dose)

Description

Calculates the posterior probability that each specified doses are the minimum effective dose in the set; i.e. the smallest dose that has a clinically significant difference (CSD).

Usage

pr_med(
  x,
  doses = attr(x, "doses"),
  csd = NULL,
  reference_dose = NULL,
  greater = TRUE,
  time = NULL
)

Arguments

Value

A tibble listing each dose and the posterior probability that each dose is the minimum efficacious dose.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e3,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

pr_med(output, csd = 10)

# difference of two doses
pr_med(output, csd = 3, reference_dose = 0)

# single model
pr_med(output$mod_quad, csd = 10)

Probability of minimum X% effective dose

Description

Calculate the probability a dose being the smallest dose that has at least X% of the maximum efficacy.

Usage

pr_medx(
  x,
  doses = attr(x, "doses"),
  ed,
  greater = TRUE,
  small_bound = 0,
  time = NULL
)

Arguments

Details

Obtaining the probability of a particular does being the minimum efficacious dose achieving ed% efficacy is dependent on the doses specified.

For a given MCMC sample of parameters, the 100% efficacy value is defined as the highest efficacy of the doses specified. For each posterior draw of MCMC parameters, the minimum ed% efficacious dose is defined as the lowest dose what has at least ed% efficacy relative to the 100% efficacy value.

The ed% effect is calculated as ed / 100 * (effect_100 - small_bound) + small_bound where effect_100 is the largest mean response among doses for a given MCMC iteration.

Value

A data frame with the following columns:

• dose: numeric dose levels.

• prob: Prob(EDX | data) for each dose. Note: these probabilities do not necessarily sum to 1 because the EDX may not exist. In fact, Pr(EDX does not exist | data) = 1 - sum(prob).

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = .1,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e3,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

pr_medx(output, ed = 80)

# single model
pr_medx(output$mod_linear, ed = 80)

Summarize Bayesian Model Averaging MCMC Output

Description

Summarize parameter inference and convergence diagnostics.

Usage

## S3 method for class 'dreamer_bma'
summary(object, ...)

Arguments

Value

Returns a named list with elements model_weights and summary containing the prior and posterior weights for each model and inference on parameters for each model as well as MCMC diagnostics.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

# all models (also show model weights)
summary(output)

# single model
summary(output$mod_linear)

Summarize Model Output

Description

Produces summaries for inference and diagnosing MCMC chains.

Usage

## S3 method for class 'dreamer_mcmc'
summary(object, ...)

Arguments

Value

A tibble with inference and diagnostics information for each parameter.

Examples

set.seed(888)
data <- dreamer_data_linear(
  n_cohorts = c(20, 20, 20),
  dose = c(0, 3, 10),
  b1 = 1,
  b2 = 3,
  sigma = 5
)

# Bayesian model averaging
output <- dreamer_mcmc(
 data = data,
 n_adapt = 1e3,
 n_burn = 1e3,
 n_iter = 1e4,
 n_chains = 2,
 silent = FALSE,
 mod_linear = model_linear(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 ),
 mod_quad = model_quad(
   mu_b1 = 0,
   sigma_b1 = 1,
   mu_b2 = 0,
   sigma_b2 = 1,
   mu_b3 = 0,
   sigma_b3 = 1,
   shape = 1,
   rate = .001,
   w_prior = 1 / 2
 )
)

# all models (also show model weights)
summary(output)

# single model
summary(output$mod_linear)