## ----include=FALSE------------------------------------------------------------ library(SAMprior) library(knitr) knitr::opts_chunk$set( fig.width = 1.62*4, fig.height = 4 ) ## setup up fast sampling when run on CRAN is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "true") ## NOTE: for running this vignette locally, please uncomment the ## following line: ## is_CRAN <- FALSE .user_mc_options <- list() if (is_CRAN) { .user_mc_options <- options(RBesT.MC.warmup=250, RBesT.MC.iter=500, RBesT.MC.chains=2, RBesT.MC.thin=1, RBesT.MC.control=list(adapt_delta=0.9)) } ## ----results="D_h",echo=FALSE------------------------------------------------- set.seed(123) std <- function(x) sd(x)/sqrt(length(x)) df_1 <- rnorm(40, 0, 3); df_2 <- rnorm(50, 0, 3); df_3 <- rnorm(60, 0, 3); dat <- data.frame(study = c(1,2,3), n = c(40, 50, 60), mean = round(c(mean(df_1), mean(df_2), mean(df_3)), 3), se = round(c(std(df_1), std(df_2), std(df_3)), 3)) kable(dat) ## ----message=FALSE------------------------------------------------------------ sigma = 3 # load R packages library(ggplot2) theme_set(theme_bw()) # sets up plotting theme set.seed(22) map_mcmc <- gMAP(cbind(mean, se) ~ 1 | study, weights=n,data=dat, family=gaussian, beta.prior=cbind(0, sigma), tau.dist="HalfNormal",tau.prior=cbind(0,sigma/2)) map_automix <- automixfit(map_mcmc) map_automix plot(map_automix)$mix ## ----message=FALSE------------------------------------------------------------ set.seed(234) sigma <- 3 ## Standard deviation in the current trial data.crt <- rnorm(35, mean = 0.4, sd = sigma) wSAM_LRT <- SAM_weight(if.prior = map_automix, delta = 0.5 * sigma, data = data.crt) cat('SAM weight: ', wSAM_LRT) ## ----message=FALSE------------------------------------------------------------ wSAM_PPR <- SAM_weight(if.prior = map_automix, delta = 0.5 * sigma, method.w = 'PPR', prior.odds = 3/7, data = data.crt) cat('SAM weight: ', wSAM_PPR) ## ----echo=FALSE, message=FALSE, warning=FALSE--------------------------------- weight_grid <- seq(-3, 3, by = 0.3) weight_res <- lapply(weight_grid, function(x){ res <- c() for(s in 1:300){ data.control <- rnorm(n = 35, mean = x, sd = sigma) res <- c(res, SAM_weight(if.prior = map_automix, delta = 0.5 * sigma, data = data.control)) } mean(res) }) df_weight <- data.frame(grid = weight_grid, weight = unlist(weight_res)) qplot(grid, weight, data = df_weight, geom = "line", main= "SAM Weight") + xlab('Sample mean from control trial')+ ylab('Weight') + geom_vline(xintercept = summary(map_automix)['mean'], linetype = 2, col = 'blue') ## ----message=FALSE------------------------------------------------------------ unit_prior <- mixnorm(nf.prior = c(1, summary(map_automix)['mean'], sigma)) SAM.prior <- SAM_prior(if.prior = map_automix, nf.prior = unit_prior, weight = wSAM_LRT, sigma = sigma) SAM.prior ## ----message=FALSE------------------------------------------------------------ # weak_prior <- mixnorm(c(1, summary(map_automix)[1], 1e4)) TypeI <- get_OC(if.prior = map_automix, ## MAP prior from historical data nf.prior = unit_prior, ## Unit-informative prior for mixture prior prior.t = mixnorm(c(1,0,1000)), ## Vague prior for treatment arm delta = 0.5*sigma, ## CSD for SAM prior ## Method to determine the mixture weight for the SAM prior method.w = 'LRT', n = 35, n.t = 70, ## Sample size for control and treatment arms if.rMAP = TRUE, ## Output robust MAP prior for comparison weight = 0.5, ## Weight for robust MAP prior ## Trial settings alternative = "greater", ## Direction of the posterior decision. Must be one of "greater" or "less" margin = 0, ## Clinical margin ## Treatment effect for control and treatment arms theta = c(summary(map_automix)['mean'], 0, -0.2, 2), theta.t = c(summary(map_automix)['mean'], -0.1, -0.2, 2)) kable(TypeI) ## ----message=FALSE------------------------------------------------------------ Power <- get_OC(if.prior = map_automix, ## MAP prior from historical data nf.prior = unit_prior, ## Unit-information prior for mixture prior prior.t = mixnorm(c(1,0,1000)), ## Vague prior for treatment arm delta = 0.5 * sigma, ## CSD for SAM prior ## Method to determine the mixture weight for the SAM prior method.w = 'LRT', n = 35, n.t = 70, ## Sample size for control and treatment arms if.rMAP = TRUE, ## Output robust MAP prior for comparison weight = 0.5, ## Weight for robust MAP prior ## Trial settings alternative = "greater", ## Direction of the posterior decision. Must be one of "greater" or "less" margin = 0, ## Clinical margin ## Treatment effect for control and treatment arms theta = c(summary(map_automix)['mean'], 0.1, 0.5, -2), theta.t = c(summary(map_automix)['mean'], 1.1, 2.0, -0.5)) kable(Power) ## ----------------------------------------------------------------------------- ## Calibrate the posterior probability cutoff PPC <- calibrate_cutoff_2arm(if.prior = map_automix, ## MAP prior from historical data nf.prior = unit_prior, ## Unit-information prior for mixture prior prior.t = mixnorm(c(1,0,1000)), ## Vague prior for treatment arm target = 0.05, ## Targeted Type I error rate n.t = 70, n = 35, ## Sample size for treatment and control arms, respectively theta.t = summary(map_automix)['mean'], ## The true effect for treatment arm theta = summary(map_automix)['mean'], ## The true effect for control arm sigma.t = sigma, sigma = sigma, ## Standard deviation in the treatment and control arms, respectively ## Method to determine the mixture weight for the SAM prior method = 'SAM', ## Method delta = 0.2, ## CSD for SAM prior ## Trial settings alternative = "greater", ## Direction of the posterior decision. Must be one of "greater" or "less". margin = 0. ## Clinical margin. ) ## ----include=F---------------------------------------------------------------- cat("Calibrated posterior probability cutoff:", round(PPC$cutoff, 4), "\n") ## ----------------------------------------------------------------------------- ## Simulate data for treatment arm data.trt <- rnorm(70, mean = 3, sd = sigma) ## first obtain posterior distributions... post_SAM <- postmix(priormix = SAM.prior, ## SAM Prior data = data.crt) post_trt <- postmix(priormix = unit_prior, ## Non-informative prior data = data.trt) ## Define the decision function decision = decision2S(PPC$cutoff, 0, lower.tail=FALSE) ## Decision-making decision(post_trt, post_SAM) ## ----------------------------------------------------------------------------- sessionInfo() ## ----include=FALSE------------------------------------------------------------ options(.user_mc_options)