if (require("RBesT")) {
###############################
## Normal outcome
###############################
data(biom)
## define shapes for which to calculate optimal contrasts
doses <- c(0, 0.05, 0.2, 0.6, 1)
modlist <- Mods(emax = 0.05, linear = NULL, logistic = c(0.5, 0.1),
linInt = c(0, 1, 1, 1), doses = doses)
## specify an informative prior for placebo, weakly informative for other arms
plc_prior <- mixnorm(inf = c(0.8, 0.4, 0.1), rob = c(0.2, 0.4, 10))
vague_prior <- mixnorm(c(1, 0, 10))
## i-th component of the prior list corresponds to the i-th largest dose
## (e.g. 1st component -> placebo prior; last component prior for top dose)
prior <- list(plc_prior, vague_prior, vague_prior, vague_prior, vague_prior)
m1 <- bMCTtest(dose, resp, biom, models=modlist, prior = prior)
## now supply a critical value (= threshold for maxmimum posterior probability)
m2 <- bMCTtest(dose, resp, biom, models=modlist, prior = prior, critV = 0.99)
####################################
## Binary outcome with covariates
####################################
if (FALSE) {
## generate data
logit <- function(p) log(p / (1 - p))
inv_logit <- function(y) 1 / (1 + exp(-y))
doses <- c(0, 0.5, 1.5, 2.5, 4)
## set seed and ensure reproducibility across R versions
set.seed(1, kind = "Mersenne-Twister", sample.kind = "Rejection", normal.kind = "Inversion")
group_size <- 100
dose_vector <- rep(doses, each = group_size)
N <- length(dose_vector)
## generate covariates
x1 <- rnorm(N, 0, 1)
x2 <- factor(sample(c("A", "B"), N, replace = TRUE, prob = c(0.6, 0.4)))
## assume approximately logit(10%) placebo and logit(35%) asymptotic response with ED50=0.5
prob <- inv_logit(emax(dose_vector, -2.2, 1.6, 0.5) + 0.3 * x1 + 0.3 * (x2 == "B"))
dat <- data.frame(y = rbinom(N, 1, prob),
dose = dose_vector, x1 = x1, x2 = x2)
## specify an informative prior for placebo (on logit scale), weakly informative for other arms
plc_prior <- mixnorm(inf = c(0.8, -2, 0.5), rob = c(0.2, -2, 10))
vague_prior <- mixnorm(c(1, 0, 10))
prior <- list(plc_prior, vague_prior, vague_prior, vague_prior, vague_prior)
## candidate models
mods <- Mods(emax = c(0.25, 1), sigEmax = rbind(c(1, 3), c(2.5, 4)), betaMod = c(1.1, 1.1),
placEff = logit(0.1), maxEff = logit(0.35)-logit(0.1),
doses = doses)
fit_cov <- glm(y~factor(dose) + 0 + x1 + x2, data = dat, family = binomial)
covariate_adjusted_estimates <- function(mu_hat, S_hat, formula_rhs,
doses, other_covariates, n_sim) {
## predict every patient under *every* dose
oc_rep <- as.data.frame(lapply(other_covariates, function(col) rep(col, times = length(doses))))
d_rep <- rep(doses, each = nrow(other_covariates))
pdat <- cbind(oc_rep, dose = d_rep)
X <- model.matrix(formula_rhs, pdat)
## average on probability scale then backtransform to logit scale
mu_star <- logit(tapply(inv_logit(X %*% mu_hat), pdat$dose, mean))
## estimate covariance matrix of mu_star
pred <- replicate(n_sim, logit(tapply(inv_logit(X %*% drop(mvtnorm::rmvnorm(1, mu_hat, S_hat))),
pdat$dose, mean)))
return(list(mu_star = as.numeric(mu_star), S_star = cov(t(pred))))
}
ca <- covariate_adjusted_estimates(coef(fit_cov), vcov(fit_cov), ~factor(dose)+0+x1+x2,
doses, dat[, c("x1", "x2")], 1000)
bMCTtest(doses, ca$mu_star, S = ca$S_star, type = "general", models = mods, prior = prior)
}
################################################
## example with contrasts handed over
################################################
data(biom)
## define shapes for which to calculate optimal contrasts
doses <- c(0, 0.05, 0.2, 0.6, 1)
modlist <- Mods(emax = 0.05, linear = NULL, sigEmax = c(0.5, 5),
linInt = c(0, 1, 1, 1), doses = doses)
## specify an informative prior for placebo, weakly informative for other arms
plc_prior <- mixnorm(inf = c(0.8, 0.4, 0.1), rob = c(0.2, 0.4, 10), sigma = 0.7)
vague_prior <- mixnorm(c(1, 0, 10), sigma = 0.7)
prior <- list(plc_prior, vague_prior, vague_prior, vague_prior, vague_prior)
## use prior effective sample sizes to calculate optimal contrasts
prior_ess <- unlist(lapply(prior, ess))
n_grp <- as.numeric(table(biom$dose))
weights <- n_grp + prior_ess
cmat <- optContr(modlist, w = weights)
bMCTtest(dose, resp, biom, models=modlist, prior = prior, contMat = cmat)
}
Run the code above in your browser using DataLab