Bayesian implementation of common dose response models using stan.
git clone https://github.com/maj-biostat/BayesDRM.git
cd BayesDRM
R CMD INSTALL .
You will need stan installed along with data.table
, ggplot2
, kableExtra
and loo
to run some of these.
Parameterised in terms of upper and lower asymptote. Emax parameter is derived.
Simulate data under this model.
emax_2p <- function(x, p0, b50){
pemax <- 1 - p0
p0 + pemax * x / (x + b50)
}
emax_2p_med <- function(pr_med, p0, b50){
b50 * (pr_med - p0) / (1-p0)
}
p0 <- 0.15
b50 <- 2
pr_med <- 0.8
x <- c(0, 5, 10, 15, 20)
p <- emax_2p(x, p0, b50)
med <- emax_2p_med(pr_med, p0, b50)
p
K <- length(p)
# 30 obs per dose
n <- 30
y <- rbinom(K, n, p)
plot(x, p, ylim = c(0, 1), type = "l")
points(x, y/n)
y <- rbinom(K, n, p)
ld <- list(N = K, y = y, n = rep(n, K), x = x,
pri_p0_a = 1, pri_p0_b = 1,
pri_b50_mu = 10, pri_b50_s = 4,
prior_only = F, pr_med = 0.8)
f1 <- BayesDRM::drm_emax2_bin(ld, refresh = 0)
f1
m <- as.matrix(f1, pars = c("p0", "b50", "p", "med", "yrep"))
# Leave one out cross validation:
log_lik_1 <- loo::extract_log_lik(f1, merge_chains = FALSE)
r_eff <- relative_eff(exp(log_lik_1), cores = 2)
loo_1 <- loo(log_lik_1, r_eff = r_eff, cores = 2)
print(loo_1)
Parameterised in terms of upper and lower asymptote. Emax parameter is derived.
Simulate data under this model.
emax_3p <- function(x, b0, bmax, b50){
b2 <- bmax - b0
eta <- b0 + b2 * x / (x + b50)
plogis(eta)
}
emax_3p_med <- function(pr_med, b0, bmax, b50){
eta_med <- qlogis(pr_med)
b2 <- bmax - b0
b50 * (eta_med - b0) / (b2 - eta_med + b0)
}
b0 <- qlogis(0.15)
bmax <- qlogis(0.95)
b50 <- 3
pr_med <- 0.8
x <- c(0, 5, 10, 15, 20)
p <- emax_3p(x, b0, bmax, b50)
med <- emax_3p_med(pr_med, b0, bmax, b50)
p
K <- length(p)
# 30 obs per dose
n <- 30
y <- rbinom(K, n, p)
plot(x, p, ylim = c(0, 1), type = "l")
points(x, y/n)
ld <- list(N = K, y = y, n = rep(n, K), x = x,
pri_b0_mu = -1.4, pri_b0_s = 2,
pri_b50_mu = 4, pri_b50_s = 3,
pri_bmax_mu = 3, pri_bmax_s = 0.5,
pr_med = 0.8, prior_only = F)
f1 <- BayesDRM::drm_emax3_bin(ld, refresh = 0)
f1
Smaller data set:
ld <- list(N = 5,
y = c(2, 4, 5, 9, 10),
n = c(10, 10, 10, 10, 10),
x = c(0, 3, 6, 12, 20),
pri_b0_mu = 0, pri_b0_s = 1.5,
pri_b50_mu = 4, pri_b50_s = 3,
# note - a stronger prior on the upper will also force
# down the lower bound.
pri_bmax_mu = 7, pri_bmax_s = 0.8,
pr_med = 0.8, prior_only = F)
f1 <- BayesDRM::drm_emax3_bin(ld, refresh = 0)
f1
log_lik_1 <- loo::extract_log_lik(f1, merge_chains = FALSE)
r_eff <- relative_eff(exp(log_lik_1), cores = 2)
loo_1 <- loo(log_lik_1, r_eff = r_eff, cores = 2)
print(loo_1)
dfig <- data.table(as.matrix(f1, pars = "p"))
dfig <- melt(dfig, measure.vars = names(dfig))
dfig <- dfig[, .(
mu = mean(value),
q_025 = quantile(value, prob = 0.025),
q_975 = quantile(value, prob = 0.975)
), keyby = variable]
dfig[, x := gsub("p[", "", variable, fixed = T)]
dfig[, x := as.integer(gsub("]", "", x, fixed = T))]
dfig[, x := ld$x[x]]
kableExtra::kable(dfig, format = "simple", digits = 2)
variable mu q_025 q_975 x
--------- ----- ------ ------ ---
p[1] 0.11 0.03 0.25 0
p[2] 0.46 0.27 0.66 3
p[3] 0.74 0.58 0.87 6
p[4] 0.92 0.85 0.97 12
p[5] 0.97 0.94 0.99 20
Parameterised in terms of upper and lower asymptote. Emax parameter is derived. Hill parameter is constrained by config, otherwise will usually be bimodal.
Simulate data under this model.
emax_4p <- function(x, b0, bmax, bh, b50){
b2 <- bmax - b0
eta <- b0 + b2 * x^bh / (x^bh + b50^bh)
plogis(eta)
}
b0 <- qlogis(0.15)
bmax <- qlogis(0.95)
b50 <- 7
bh <- 3
x <- c(0, 5, 10, 15, 20)
p <- emax_4p(x, b0, bmax, bh, b50)
p
K <- length(p)
# 30 obs per dose
n <- 30
y <- rbinom(K, n, p)
plot(x, p, ylim = c(0, 1), type = "l")
points(x, y/n)
ld <- list(N = K, y = y, n = rep(n, K), x = x,
pri_b0_mu = -1.4, pri_b0_s = 2,
pri_b50_mu = 4, pri_b50_s = 3,
pri_bmax_mu = 3, pri_bmax_s = 0.5,
pri_bh_mu = 1, pri_bh_s = 2,
# constrain hill to positive
pos_hill = 1, prior_only = F)
f1 <- BayesDRM::drm_emax4_bin(ld, refresh = 0)
f1