inst/doc/tboot_bmr.R

In tboot: Tilted Bootstrap

## ----sim, cache = FALSE-------------------------------------------------------
library(tboot)
set.seed(2020)
quant1  <- rnorm(200) + 1
bin1    <- ifelse( (.5*quant1 + .5*rnorm(200)) > .5, 1, 0)
bin2    <- ifelse( (.5*quant1 + .5*rnorm(200)) > .5, 1, 0)
simData <- data.frame(quant1, bin1, bin2)
head(simData)

## ----target, cache = FALSE----------------------------------------------------
marginal_active <-  list(quant1=rnorm(5000, mean=.5, sd=.2),
                         bin1=rbeta(5000, shape1 = 50,shape2=50),
                         bin2=rbeta(5000, shape1 = 60,shape2=40))
marginal_pbo <-  list(quant1=rnorm(5000, mean=.2, sd=.2),
                      bin1=rbeta(5000, shape1 = 20,shape2=80),
                      bin2=rbeta(5000, shape1 = 30,shape2=70))

## ----weights, cache = FALSE---------------------------------------------------
bmr_active <- tweights_bmr(dataset = simData, marginal = marginal_active)
bmr_pbo <- tweights_bmr(dataset = simData, marginal = marginal_pbo)

## ----post_bmr, cache = FALSE, fig.height = 6, fig.width = 6, fig.align = "center"----
samples <- rbind(data.frame(trt="active", post_bmr(nsims=1e3, bmr_active)),
                 data.frame(trt="pbo", post_bmr(nsims=1e3, bmr_pbo)))
head(samples)

## ----pairsgraph, cache = FALSE, fig.height = 6, fig.width = 6, fig.align = "center"----
pairs(samples[,-1], col=ifelse(samples=="active","red","blue"), pch='.', cex=.5)

## ----marginal dist, cache = FALSE, fig.height = 6, fig.width = 6, fig.align = "center"----
library(ggplot2)
pltdta=do.call(rbind, lapply(c("quant1","bin1", "bin2"),
                             function(nm) {
                               rbind(data.frame(type="BMR", var=nm, trt=samples$trt,
                                                val=samples[[nm]]), 
                                     data.frame(type="marginal", var=nm,
                                                trt="active", 
                                                val=marginal_active[[nm]]),
                                     data.frame(type="marginal", var=nm, 
                                                trt="pbo", 
                                                val=marginal_pbo[[nm]]))
                             }))

  ggplot(pltdta, aes(fill=type, x=val)) + 
  geom_density(alpha=.3) + facet_grid(var~trt, scales = "free")



## ----tboot_bmr, cache = FALSE-------------------------------------------------
active_sample=tboot_bmr(nrow=100, weights_bmr=bmr_active)
head(active_sample)
                   

## ----attr, cache = FALSE------------------------------------------------------
attr(active_sample, "post_bmr")
                   

## ----sim_trial, cache = FALSE-------------------------------------------------
#Manage any errors by assuming the pvalue failed to reach statistical
#significance (i.e. pvalue is 1) but keep track of such errors. 
errorTrackGlobal=list()
manageError=function(expr)  {
  tryCatch(eval(quote(expr)), error=function(e){
    errorTrackGlobal[[length(errorTrackGlobal)+1]] <<- e$message
    return(1)
  }) 
}


#create function to simulate and analyze one virtual trial
sim_and_analyze=function() {
  active_sample=tboot_bmr(100, bmr_active)
  pbo_sample=tboot_bmr(100, bmr_pbo)
  data.frame(
    p_quant1=manageError(t.test(active_sample$quant1,pbo_sample$quant1)$p.value),
    p_bin1=manageError(fisher.test(active_sample$bin1,pbo_sample$bin1)$p.value),
    p_bin2=manageError(fisher.test(active_sample$bin2,pbo_sample$bin2)$p.value)
  )
}
#Simulate Pvalues
p_sim=do.call(rbind, replicate(100, sim_and_analyze(), simplify = FALSE))
head(errorTrackGlobal)
head(p_sim)