Nothing
## ---- echo = FALSE------------------------------------------------------------
knitr::opts_chunk$set(fig.width = 6, fig.height = 4)
## ---- message = FALSE---------------------------------------------------------
library(dreamer)
library(dplyr)
library(ggplot2)
set.seed(888)
data <- dreamer_data_quad(
n_cohorts = c(20, 20, 20, 20, 20), # number of subjects in each cohort
doses = c(0, .5, 1, 1.5, 2), # dose administered to each cohort
b1 = 5, # intercept
b2 = 5,
b3 = - 2,
sigma = 2 # standard deviation
)
head(data)
ggplot(data, aes(dose, response)) + geom_point()
## -----------------------------------------------------------------------------
model_quad(
mu_b1 = 0,
sigma_b1 = 10,
mu_b2 = 0,
sigma_b2 = 10,
mu_b3 = 0,
sigma_b3 = 10,
shape = 1,
rate = .001,
w_prior = 1 / 2
)
## -----------------------------------------------------------------------------
# Bayesian model averaging
output <- dreamer_mcmc(
data = data,
n_adapt = 1e3,
n_burn = 1e3,
n_iter = 1e4,
n_chains = 2,
silent = TRUE,
# the argument name "mod_quad" is arbitrary and is chosen by the user
mod_quad = model_quad(
mu_b1 = 0,
sigma_b1 = 10,
mu_b2 = 0,
sigma_b2 = 10,
mu_b3 = 0,
sigma_b3 = 10,
shape = 1,
rate = .001,
w_prior = 1 / 2
),
# the argument name "mod_emax" is arbitrary and is chosen by the user
mod_emax = model_emax(
mu_b1 = 0,
sigma_b1 = 10,
mu_b2 = 0,
sigma_b2 = 10,
mu_b3 = 0,
sigma_b3 = 2,
mu_b4 = 1,
sigma_b4 = 5,
shape = 1,
rate = .001,
w_prior = 1 / 2
)
)
## -----------------------------------------------------------------------------
output
## -----------------------------------------------------------------------------
summary(output)
## -----------------------------------------------------------------------------
posterior(output)
# posterior Pr(effect of dose - effect of reference dose)
posterior(x = output, reference_dose = 0)
## -----------------------------------------------------------------------------
plot(output)
## -----------------------------------------------------------------------------
plot(output, data = data)
## -----------------------------------------------------------------------------
plot(output, reference_dose = 0) # adjust relative to dose = 0
## -----------------------------------------------------------------------------
plot(output, data = data, predictive = 10)
## -----------------------------------------------------------------------------
plot(output, reference_dose = 0, predictive = 10)
## -----------------------------------------------------------------------------
plot(output$mod_emax, data = data)
## -----------------------------------------------------------------------------
plot_comparison(output)
## -----------------------------------------------------------------------------
plot_comparison(
mod_emax = output$mod_emax,
mod_quad = output$mod_quad
)
## -----------------------------------------------------------------------------
# absolute: pr(mean at dose 0.50 > 1 | data)
pr_eoi(output, eoi = 1, dose = 0.50)
# relative: pr( (mean at dose 0.50) - (mean at dose 0.25) > 0.55 | data )
pr_eoi(output, eoi = 0.55, dose = 0.50, reference_dose = 0.25)
# vectorized
n_doses <- length(output$doses)
pr_eoi(
output,
eoi = rep(.55, n_doses),
dose = output$doses,
reference_dose = rep(0, n_doses)
)
## -----------------------------------------------------------------------------
# from the quadratic model
pr_eoi(
output$mod_quad,
eoi = 1.5,
dose = .75,
reference_dose = 0
)
## -----------------------------------------------------------------------------
pr_med(output, csd = 8)
# placebo adjusted
pr_med(
output,
csd = 3,
reference_dose = 0
)
# relative to placebo grp (dose = 0) for just the EMAX model
pr_med(
output$mod_emax,
csd = 3,
reference_dose = 0
)
## -----------------------------------------------------------------------------
# looking for smallest dose with 95% of the maximum efficacy
pr_medx(output, ed = 95)
## -----------------------------------------------------------------------------
post_medx(output, ed = 95)
## -----------------------------------------------------------------------------
post_perc_effect(
output,
dose = c(.05, .5)
)
## -----------------------------------------------------------------------------
diagnostics(output)
# single model
diagnostics(output$mod_emax)
## ---- fig.height = 8, fig.width = 6-------------------------------------------
# single model
plot_trace(output$mod_quad)
# traceplot for all parameters for each model using: plot_trace(output)
## -----------------------------------------------------------------------------
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 = 1
)
)
## -----------------------------------------------------------------------------
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
)
)
## -----------------------------------------------------------------------------
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
)
)
## -----------------------------------------------------------------------------
output_independent <- dreamer_mcmc(
data = data,
n_adapt = 1e3,
n_burn = 1e3,
n_iter = 1e4,
n_chains = 2,
silent = TRUE, # make rjags be quiet,
# this model has the same prior on the mean for each dose
mod_indep = model_independent(
mu_b1 = 0,
sigma_b1 = 1,
shape = 1,
rate = .001
)
)
# prior is too strong!
plot(output_independent, data = data)
output_independent2 <- dreamer_mcmc(
data = data,
n_adapt = 1e3,
n_burn = 1e3,
n_iter = 1e4,
n_chains = 2,
silent = TRUE, # make rjags be quiet,
# this model has the different priors on the mean for each dose
mod_indep = model_independent(
mu_b1 = c(0, 1, 2, 3, 4),
sigma_b1 = c(10, 10, 20, 20, 30),
shape = 1,
rate = .001,
doses = c(0, 0.5, 1, 1.5, 2)
)
)
plot(output_independent2, data = 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
)
ggplot(data_long, aes(time, response, group = dose, color = factor(dose))) +
geom_point()
## -----------------------------------------------------------------------------
# Bayesian model averaging
output_long <- dreamer_mcmc(
data = data_long,
n_adapt = 1e3,
n_burn = 1e3,
n_iter = 1e4,
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
)
)
)
## -----------------------------------------------------------------------------
plot(output_long, data = data_long)
# plot dose response at last time point
plot(output_long, times = 52, data = data_long)
## -----------------------------------------------------------------------------
posterior(output_long)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.