###### ########################################################

rm(list=ls()) #will remove ALL objects

##############################################################

Bayes-Factor Calculations for T-tests

##############################################################

#Start of Settings

### Give a title for results output

Results.Title = ‘Normal(x,0,.5) N = 100 BS-Design, Obs.ES = 0′

### Criterion for Inference in Favor of H0, BF (H1/H0)

BF.crit.H0 = 1/3

### Criterion for Inference in Favor of H1

#set z.crit.H1 to Infinity to use Bayes-Factor, BF(H1/H0)

BF.crit.H1 = 3

z.crit.H1 = Inf

### Set Number of Groups

gr = 2

### Set Total Sample size

N = 100

### Set observed effect size

### for between-subject designs and one sample designs this is Cohen’s d

### for within-subject designs this is dz

obs.es = 0

### Set the mode of the alternative hypothesis

alt.mode = 0

### Set the variability of the alternative hypothesis

alt.var = .5

### Set the shape of the distribution of population effect sizes

alt.dist = 2 #1 = Cauchy; 2 = Normal

### Set the lower bound of population effect sizes

### Set to zero if there is zero probability to observe effects with the opposite sign

low = -3

### Set the upper bound of population effect sizes

### For example, set to 1, if you think effect sizes greater than 1 SD are unlikely

high = 3

### set the precision of density estimation (bigger takes longer)

precision = 100

### set the graphic resolution (higher resolution takes longer)

graphic.resolution = 20

### set limit for non-central t-values

nct.limit = 100

################################

# End of Settings

################################

# compute degrees of freedom

df = (N – gr)

# get range of population effect sizes

pop.es=seq(low,high,(1/precision))

# compute sampling error

se = gr/sqrt(N)

# limit population effect sizes based on non-central t-values

pop.es = pop.es[pop.es/se >= -nct.limit & pop.es/se <= nct.limit]

# function to get weights for Cauchy or Normal Distributions

get.weights=function(pop.es,alt.dist,p) {

if (alt.dist == 1) w = dcauchy(pop.es,alt.mode,alt.var)

if (alt.dist == 2) w = dnorm(pop.es,alt.mode,alt.var)

sum(w)

# get the scaling factor to scale weights to 1*precision

#scale = sum(w)/precision

# scale weights

#w = w / scale

return(w)

}

# get weights for population effect sizes

weights = get.weights(pop.es,alt.dist,precision)

#Plot Alternative Hypothesis

Title=”Alternative Hypothesis”

ymax=max(max(weights)*1.2,1)

plot(pop.es,weights,type=’l’,ylim=c(0,ymax),xlab=”Population Effect Size”,ylab=”Density”,main=Title,col=’blue’,lwd=3)

abline(v=0,col=’red’)

#create observations for plotting of prediction distributions

obs = seq(low,high,1/graphic.resolution)

# Get distribution for observed effect size assuming H1

H1.dist = as.numeric(lapply(obs, function(x) sum(dt(x/se,df,pop.es/se) * weights)/precision))

#Get Distribution for observed effect sizes assuming H0

H0.dist = dt(obs/se,df,0)

#Compute Bayes-Factors for Prediction Distribution of H0 and H1

BFs = H1.dist/H0.dist

#Compute z-scores (strength of evidence against H0)

z = qnorm(pt(obs/se,df,log.p=TRUE),log.p=TRUE)

# Compute H1 error rate rate

BFpos = BFs

BFpos[z < 0] = Inf

if (z.crit.H1 == Inf) z.crit.H1 = abs(z[which(abs(BFpos-BF.crit.H1) == min(abs(BFpos-BF.crit.H1)))])

ncz = qnorm(pt(pop.es/se,df,log.p=TRUE),log.p=TRUE)

weighted.power = sum(pnorm(abs(ncz),z.crit.H1)*weights)/sum(weights)

H1.error = 1-weighted.power

#Compute H0 Error Rate

z.crit.H0 = abs(z[which(abs(BFpos-BF.crit.H0) == min(abs(BFpos-BF.crit.H0)))])

H0.error = (1-pnorm(z.crit.H0))*2

# Get density for observed effect size assuming H0

Density.Obs.H0 = dt(obs.es,df,0)

# Get density for observed effect size assuming H1

Density.Obs.H1 = sum(dt(obs.es/se,df,pop.es/se) * weights)/precision

# Compute Bayes-Factor for observed effect size

BF.obs.es = Density.Obs.H1 / Density.Obs.H0

#Compute z-score for observed effect size

obs.z = qnorm(pt(obs.es/se,df,log.p=TRUE),log.p=TRUE)

#Show Results

ymax=max(H0.dist,H1.dist)*1.3

plot(type=’l’,z,H0.dist,ylim=c(0,ymax),xlab=”Strength of Evidence (z-value)”,ylab=”Density”,main=Results.Title,col=’black’,lwd=2)

par(new=TRUE)

plot(type=’l’,z,H1.dist,ylim=c(0,ymax),xlab=””,ylab=””,col=’blue’,lwd=2)

abline(v=obs.z,lty=2,lwd=2,col=’darkgreen’)

abline(v=-z.crit.H1,col=’blue’,lty=3)

abline(v=z.crit.H1,col=’blue’,lty=3)

abline(v=-z.crit.H0,col=’red’,lty=3)

abline(v=z.crit.H0,col=’red’,lty=3)

points(pch=19,c(obs.z,obs.z),c(Density.Obs.H0,Density.Obs.H1))

res = paste0(‘BF(H1/H0): ‘,format(round(BF.obs.es,3),nsmall=3))

text(min(z),ymax*.95,pos=4,res)

res = paste0(‘BF(H0/H1): ‘,format(round(1/BF.obs.es,3),nsmall=3))

text(min(z),ymax*.90,pos=4,res)

res = paste0(‘H1 Error Rate: ‘,format(round(H1.error,3),nsmall=3))

text(min(z),ymax*.80,pos=4,res)

res = paste0(‘H0 Error Rate: ‘,format(round(H0.error,3),nsmall=3))

text(min(z),ymax*.75,pos=4,res)

######################################################

### END OF Subjective Bayesian T-Test CODE

######################################################

### Thank you to Jeff Rouder for posting his code that got me started.

### http://jeffrouder.blogspot.ca/2016/01/what-priors-should-i-use-part-i.html