Packages needed

You will need to install and load the following packages for these exercises:

install.packages("BayesFactor")
install.packages("MASS")
install.packages("devtools") # if not yet installed
library(devtools)
install_github("nicebread/BFDA", subdir="package")

library(BFDA)
library(BayesFactor)
library(MASS)

Exercise 2: Simulating data using BFDA and selecting trial parameters

The Bayes factor design analysis (BFDA) package provides a simple way to simulate and analyse sequential bayesian analyses for trial planning. For more complex analyses you can create your own custom simulations (e.g. see the scripts at https://osf.io/8mxda/ or https://osf.io/b397u/  for examples)

Please see https://rawgit.com/nicebread/BFDA/master/package/doc/BFDA_manual.html for the BFDA manual and full instructions

set.seed(19112022)

#for the following I use the naming convention sim.Hd followed by the simulated effect size, e.g. sim.Hd0.5 is for simulating the hypothesis of a between-group effect size of d = 0.5

sim.Hd0.5 <- BFDA.sim(expected.ES=0.5, type="t.between",prior=list("Cauchy",list(prior.location=0, prior.scale=sqrt(2)/2)), n.min=10, n.max=80, alternative="greater", boundary=Inf, B=1000, verbose=TRUE, stepsize = 1)

sim.Hd0.0 <- BFDA.sim(expected.ES=0, type="t.between", prior=list("Cauchy", list(prior.location=0, prior.scale=sqrt(2)/2)), n.min=10, n.max=80, alternative="greater", boundary=Inf, B=1000, verbose=TRUE, stepsize = 1)

#If you like (for later) you can now also run some in-between or larger values, e.g.

sim.Hd0.2 <- BFDA.sim(expected.ES=0.2, type="t.between", prior=list("Cauchy", list(prior.location=0, prior.scale=sqrt(2)/2)), n.min=10, n.max=80, alternative="greater", boundary=Inf, B=1000, verbose=TRUE, stepsize = 1)

sim.Hd0.4 <- BFDA.sim(expected.ES=0.4, type="t.between", prior=list("Cauchy", list(prior.location=0, prior.scale=sqrt(2)/2)), n.min=10, n.max=80, alternative="greater", boundary=Inf, B=1000, verbose=TRUE, stepsize = 1)

sim.Hd0.6 <- BFDA.sim(expected.ES=0.6, type="t.between", prior=list("Cauchy", list(prior.location=0, prior.scale=sqrt(2)/2)), n.min=10, n.max=80, alternative="greater", boundary=Inf, B=1000, verbose=TRUE, stepsize = 1)

#etc...

#Saving a BFDA simulation object:
saveRDS(sim.Hd0.5,"sim.Hd0.5.dd.mm.yy.RDS")
saveRDS(sim.Hd0.0,"sim.Hd0.0.dd.mm.yy.RDS")
#to load again
sim.Hd0.5<-readRDS("sim.Hd0.5.dd.mm.yy.RDS")
sim.Hd0.0<-readRDS("sim.Hd0.0.dd.mm.yy.RDS")

#First plot for H1
dev.new() #sometimes the plot doesn't work if you don't make a new window first
plot(sim.Hd0.5, n.min=10, n.max=80, boundary=c(1/5, 5), n.trajectories = 60)

#Then plot for H0
dev.new() #sometimes the plot doesn't work if you don't make a new window first
plot(sim.Hd0.0, n.min=10, n.max=80, boundary=c(1/5, 5), n.trajectories = 60)

#There are other functions in BFDA to explore the simulations/sequential analyses but I find the plots the most informative. However, feel free to explore.

Overview and setup

In this exercise you will analyse some (simulated) data using the trial parameters you selected, and find out when you hit the BF boundaries.

First you will need to simulate the trial data - please run the code below (which I’m happy to explain but there’s no need to understand it for the exercise - it’s just a means to provide you with the data to analyse):

#Simulates 100 participants in each arm of trial where:
#There is a correlation of r = 0.35 between pre and post-data
#One arm is not superior to control, the other is
set.seed(19112022)
samples=100
r=0.35
data<-mvrnorm(n=samples,mu=c(1,0.8),Sigma=matrix(c(1,r,r,1),nrow=2),empirical=FALSE)
X = data[,1]
Y = data[,2]
X<-as.integer(X*5+20)
Y<-as.integer(Y*5+20)
data<-mvrnorm(n=samples,mu=c(1,0.15),Sigma=matrix(c(1,r,r,1),nrow=2),empirical=FALSE)
X1 = data[,1]
Y1 = data[,2]
X1<-as.integer(X1*5+20)
Y1<-as.integer(Y1*5+20)
data<-mvrnorm(n=samples,mu=c(1,0.9),Sigma=matrix(c(1,r,r,1),nrow=2),empirical=FALSE)
X2 = data[,1]
Y2 = data[,2]
X2<-as.integer(X2*5+20)
Y2<-as.integer(Y2*5+20)
cdiff<-Y-X
T1diff<-Y1-X1
T2diff<-Y2-X2
predata<-c(X,X1,X2)
postdata<-c(Y,Y1,Y2)
diffdata<-postdata-predata
group<-c(rep("C",samples),rep("Tx1",samples),rep("Tx2",samples))
dseq<-seq(1,samples*3)
pid<-as.character(seq(1,samples*3))
mydata<-(cbind(predata,postdata,diffdata,group))
mydata<-mydata[sample(1:nrow(mydata)), ]
mydata<-data.frame(cbind(pid,dseq,mydata))
mydata$dseq<-as.integer(mydata$dseq)
mydata$predata<-as.integer(mydata$predata)
mydata$postdata<-as.integer(mydata$postdata)
mydata$diffdata<-as.integer(mydata$diffdata)

You now have a dataframe (mydata) with pre and post outcome data for a pretend 3-arm trial of treatments for depression (i.e. a decrease in score on the outcome measure is good). There are three treatment arms, “C” (control condition, e.g. TA), “Tx1” (New treatment 1) and “Tx2” (New treatment 2). One arm is superior to control, one isn’t. Each arm includes 100 participants.

Dataframe columns are:

• pid = participant id (a string of a number from 1 to 300)
• dseq = sequence in which they provided outcome data for the trial (numerical, here the same as pid for convenience)
• predata = score on the depression outcome measure at pre-treatment
• postdata = score on the depression outcome measure at post-treatment
• diffdata = post-treatment score minus pre-treatment score
• group = which group (C, Tx1, Tx2)

Exploring the data

Take a look at the data (e.g. with View(mydata) ) to get a sense of it. You can then use the functions in the code chunk below to explore the potential outcome of applying sequential Bayesian analyses.

First, run the code below to load the two functions:

BFsnapshot<-function(BFdata,n=300,rs=(sqrt(2)/2)){

  if (length(BFdata$pid)<n){
    return("Error: n larger than number of participants");
  }
  
  tryCatch(
        {
          tdata<-BFdata[1:n,]
          cat("\n\n Total N = ",n,"\n\n")
          cat("Control: n = ",length(tdata[tdata$group=="C",]$pid),", mean change (post minus pre-treatment) = ", mean(tdata[tdata$group=="C",]$diffdata,na.rm=TRUE)," (SD = ", sd(tdata[tdata$group=="C",]$diffdata,na.rm=TRUE),")\n\n",sep="")
         cat("Tx1: n = ",length(tdata[tdata$group=="Tx1",]$pid),", mean change (post minus pre-treatment) = ", mean(tdata[tdata$group=="Tx1",]$diffdata,na.rm=TRUE)," (SD = ", sd(tdata[tdata$group=="Tx1",]$diffdata,na.rm=TRUE),")\n\n",sep="")
         #effect size calculation
         d<-effectsize::cohens_d(tdata[tdata$group=="C",]$diffdata,tdata[tdata$group=="Tx1",]$diffdata)
         cat("Effect size vs. control: d=",d$Cohens_d," 95% CIs [",d$CI_low,",",d$CI_high,"]\n\n",sep="")
         #Calculates directional Bayesian t-test (nullinterval=c(0,Inf))
         BF<-BayesFactor::ttestBF(x=tdata[tdata$group=="C",]$diffdata,y=tdata[tdata$group=="Tx1",]$diffdata,nullinterval=c(0,Inf),rscale=rs)
        cat("BF vs. control: BF=",exp(BF@bayesFactor$bf),"\n\n",sep="")

        cat("Tx2: n = ",length(tdata[tdata$group=="Tx2",]$pid),", mean change (post minus pre-treatment) = ", mean(tdata[tdata$group=="Tx2",]$diffdata,na.rm=TRUE)," (SD = ", sd(tdata[tdata$group=="Tx2",]$diffdata,na.rm=TRUE),")\n\n",sep="")
        #effect size calculation
         d<-effectsize::cohens_d(tdata[tdata$group=="C",]$diffdata,tdata[tdata$group=="Tx2",]$diffdata)
         cat("Effect size vs. control: d=",d$Cohens_d," 95% CIs [",d$CI_low,",",d$CI_high,"]\n\n",sep="")
         #Calculates directional Bayesian t-test (nullinterval=c(0,Inf))
         BF<-BayesFactor::ttestBF(x=tdata[tdata$group=="C",]$diffdata,y=tdata[tdata$group=="Tx2",]$diffdata,nullinterval=c(0,Inf),rscale=rs)
        cat("BF vs. control: BF=",exp(BF@bayesFactor$bf),"\n\n",sep="") 
          
        },
        error=function(cond) {
          return (paste0("Something went wrong! ",cond))
        })
}

seqBFs<-function(BFdata,rs=(sqrt(2)/2)){
  tryCatch(
        {
          alldata<-BFdata[1:length(BFdata$pid),]
          dsTx1<-vector()
          dsTx2<-vector()
          BFsTx1<-vector()
          BFsTx2<-vector()
          nsTx1<-vector()
          nsTx2<-vector()
          is<-vector()
          #start at index number 10 to avoid problems
          for (i in 10:length(alldata$pid)){
          tdata<-alldata[1:i,]
          Ns<-as.data.frame(cbind(table(tdata[1:i,]$group)))
          dsTx1<-c(dsTx1,(effectsize::cohens_d(tdata[tdata$group=="C",]$diffdata,tdata[tdata$group=="Tx1",]$diffdata))$Cohens_d)
          BFsTx1<-c(BFsTx1,exp((BayesFactor::ttestBF(x=tdata[tdata$group=="C",]$diffdata,y=tdata[tdata$group=="Tx1",]$diffdata,nullinterval=c(0,Inf),rscale=rs)@bayesFactor$bf)))
            nsTx1<-c(nsTx1,Ns["Tx1",])
            
            dsTx2<-c(dsTx2,(effectsize::cohens_d(tdata[tdata$group=="C",]$diffdata,tdata[tdata$group=="Tx2",]$diffdata))$Cohens_d)
            BFsTx2<-c(BFsTx2,exp((BayesFactor::ttestBF(x=tdata[tdata$group=="C",]$diffdata,y=tdata[tdata$group=="Tx2",]$diffdata,nullinterval=c(0,Inf),rscale=rs)@bayesFactor$bf)))
            nsTx2<-c(nsTx2,Ns["Tx2",])
            is<-c(is,i)
          
        }
          allresults<-data.frame(cbind(is,nsTx1,dsTx1,BFsTx1,nsTx2,dsTx2,BFsTx2))
          return(allresults)
        },
        error=function(cond) {
          return (paste0("Something went wrong! ",cond))
        })
}

You can now use these two functions to explore the effect of doing sequential Bayesian analyses (here for simplicity, t-test on change scores).