Bayesian t-test hypothesis testing for two independent groups
14 Jul 2012 - Michael Tsikerdekis
Introduction
This method can be used in the same circumstances that one would use the regular independent t-test; when you want to statistically compare the means of two groups. Both groups should have their data normally distributed.
References
This method is based on the book "A Practical Course in Bayesian Graphical Modeling" by Michael Lee and Eric-Jan Wagenmakers. Additionally, a published scientific article can be found here. Either or both are good to cite when using this method. Some of the code may has been changed in order to make the application of the analysis easier. If you want to learn more about the model and the code you can read the book or the article.
For the procedure you need R and Openbugs.
Procedure highlights
Input
setwd("/the/directory/for/both/of/your/files/") #this will help the program find the location of the model
group1 = c(49,13,50,64,21,23,15,7,32,8,17,15,19,16,33) #your set of values for group 1
group2 = c(41,20,53,41,20,44,31,24,24,44,15,35,32,25,35) #your set of values for group 2
openbugsmodel = "Ttest_2.txt" #specify the location for the file that contains the OpenBugs code
priorforh1 = dcauchy(0) #prior for the case of delta<>0
priorforh2 = 2*dcauchy(0) #prior for the case of delta<0
priorforh3 = 2*dcauchy(0) #prior for the case of delta>0
#Advanced input
itterations = 30000 #Don't change these unless you know what you're doing
burnin = 3000 #Don't change these unless you know what you're doing
group1 = c(49,13,50,64,21,23,15,7,32,8,17,15,19,16,33) #your set of values for group 1
group2 = c(41,20,53,41,20,44,31,24,24,44,15,35,32,25,35) #your set of values for group 2
openbugsmodel = "Ttest_2.txt" #specify the location for the file that contains the OpenBugs code
priorforh1 = dcauchy(0) #prior for the case of delta<>0
priorforh2 = 2*dcauchy(0) #prior for the case of delta<0
priorforh3 = 2*dcauchy(0) #prior for the case of delta>0
#Advanced input
itterations = 30000 #Don't change these unless you know what you're doing
burnin = 3000 #Don't change these unless you know what you're doing
Set the first three lines according to your setup and data or feed the variables your own data. You also need to set the file that contains the Openbugs model which you can find at the end of this page. If you wish to change the priors you can, just remember that you need to adjust the prior for hypothesis 2 and hypothesis 3 in order apply only for positive or negative numbers. Also you can change the iterations and the burnin if you want to improve your results. These need to be reported in your paper later on.
Output
The hypotheses for this test are:
- H0 - δ = 0: there is no difference between the two groups
- H1 - δ <>0: there is a difference between the two groups
- H2 - δ < 0: group 2 has larger values than group 1
- H3 - δ > 0: group 1 has larger values than group 2
The output produces a set of results in text along with the probability distribution plots for each one of them. Both are useful for making a decision about your hypothesis. As an example, you can use the graph and determine if at the point δ=0 your posterior(your results after the data) are higher or lower than the prior(your initial belief). If the posterior is higher than the prior at δ=0 then it reinforces the fact that the null hypothesis(H0) is probably true. If the posterior is lower than the prior then the data weakens your belief that the null hypothesis is true. Bayes factors are automatically reported on the graphs.
-----------Results--------------
Bayes Factor(BF10) for H1 delta<>0 over H0 delta=0: 0.5346907
Bayes Factor(BF01) for H0 delta=0 over H1 delta<>0: 1.87024
---
Bayes Factor(BF20) for H2 delta<0 over H0 delta=0: 0.9432402
Bayes Factor(BF02) for H0 delta=0 over H2 delta<0: 1.060175
---
Bayes Factor(BF30) for H3 delta> over H0 delta=0: 0.1256324
Bayes Factor(BF03) for H0 delta=0 over H3 delta>0: 7.959729
---
Bayes Factor(BF10) for H1 delta<>0 over H0 delta=0: 0.5346907
Bayes Factor(BF01) for H0 delta=0 over H1 delta<>0: 1.87024
---
Bayes Factor(BF20) for H2 delta<0 over H0 delta=0: 0.9432402
Bayes Factor(BF02) for H0 delta=0 over H2 delta<0: 1.060175
---
Bayes Factor(BF30) for H3 delta> over H0 delta=0: 0.1256324
Bayes Factor(BF03) for H0 delta=0 over H3 delta>0: 7.959729
---
Please read the Bayes Factor page for how to interpret it.
In the first two cases the evidence is "Barely worth mentioning" for H0. But, the third result (7.67) is considered "Substantial" evidence in favor of H0, indicating that when we consider if group 1 has a bigger effect than group 2, there is substantial evidence to say that is unlikely(providing proof for H0).
When publishing you are going to have to report also the process you obtained your results and the numbers of itterations for the MCMC test along with the burnin value and the chains(in this case 3).
Code
R code
# clears workspace:
rm(list=ls(all=TRUE))
#---------------INPUT DATA------------------
setwd("/the/directory/for/both/of/your/files/") #this will help the program find the location of the modelthe model
group1 = c(1,5,4,2,1,2,3,1,2,3,1,2,5,4,7,8,9,7,5,3,2,3,4,6,7,6,4,7) #your set of values for group 1
group2 = c(5,4,3,3,4,6,6,8,7,6,5,7,7,6,5,6,7,5,4,5,6,7,8,5,6,7,8,9) #your set of values for group 2
openbugsmodel = "Ttest_2.txt"
priorforh1 = dcauchy(0) #prior for the case of delta<>0
priorforh2 = 2*dcauchy(0) #prior for the case of delta<0
priorforh3 = 2*dcauchy(0) #prior for the case of delta>0
#Advanced input
itterations = 30000 #Don't change these unless you know what you're doing
burnin = 3000 #Don't change these unless you know what you're doing
#-------------------------------------------
library(R2OpenBUGS)
n1 = length(group1)
n2 = length(group2)
# Rescale:
group2 = group2-mean(group1)
group2 = group2/sd(group1)
group1 = as.vector(scale(group1))
data=list("group1", "group2", "n1", "n2") # to be passed on to WinBUGS
inits=function()
{
list(delta=rnorm(1,0,1),mu=rnorm(1,0,1),sigma=runif(1,0,5))
}
# Parameters to be monitored
parameters=c("delta")
# The following command calls WinBUGS with specific options.
# For a detailed description see Sturtz, Ligges, & Gelman (2005).
samples = bugs(data, inits, parameters,
model.file =openbugsmodel,
n.chains=3, n.iter=itterations, n.burnin=burnin, n.thin=1,
DIC=T,
codaPkg=F, debug=F)
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
samples$summary # overview
# Please work through the analyses below one at a time
######################################################
# Analysis 1. H1: delta ~ Cauchy (unrestricted case)
######################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
#============ BFs based on logspline fit ===========================
library(polspline) # this package can be installed from within R
fit.posterior = logspline(delta.posterior)
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh1 # height of order--restricted prior at delta = 0
BF10 = prior/posterior
BF01 = posterior/prior
cat("-----------Results--------------\r\n")
cat("Bayes Factor(BF10) for H1 delta<>0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF01) for H0 delta=0 over H1 delta<>0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[1][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][1])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = -3
xhigh = 3
yhigh = 4
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(-4,-3,-2,-1,0,1,2,3,4), lab=c("-4","-3","-2","-1","0", "1", "2", "3", "4"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(-4,-3,-2,-1,0,1,2,3,4), lab=c("-4","-3","-2","-1","0", "1", "2", "3", "4"))
axis(2)
points(0, dcauchy(0), pch=19, cex=2)
text(-1,3.5, expression(H[0]: delta == 0),cex=2)
text(-1,3, expression(H[1]: delta != 0),cex=2)
text(-1,2.5, bquote(.(BFtext) == .(BFplot)),cex=2)
###########################################################################
# Analysis 2. H1: delta ~ Cauchy^- (restricted case, negative values only)
###########################################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
# selects only negative delta's:
delta.posterior = delta.posterior[which(delta.posterior<0)]
#============ BFs based on logspline fit ===========================
fit.posterior = logspline(delta.posterior,ubound=0) # NB. note the bound
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh2 # height of order--restricted prior at delta = 0
BF10 = prior/posterior
BF01 = posterior/prior
cat("Bayes Factor(BF20) for H2 delta<0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF02) for H0 delta=0 over H2 delta<0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[2][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][2])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = -3
xhigh = 0
yhigh = 12
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(-3,-2,-1,0), lab=c("-3","-2","-1","0"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) 2*dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(-3,-2,-1,0), lab=c("-3","-2","-1","0"))
axis(2)
points(0, 2*dcauchy(0), pch=19, cex=2)
text(-2,10, expression(H[0]: delta == 0),cex=2)
text(-2,8, expression(H[2]: delta < 0),cex=2)
text(-2,6, bquote(.(BFtext) == .(BFplot)),cex=2)
###########################################################################
# Analysis 3. H1: delta ~ Cauchy^+ (restricted case, positive values only)
###########################################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
# selects only positive delta's:
delta.posterior = delta.posterior[which(delta.posterior>0)]
#============ BFs based on logspline fit ===========================
fit.posterior = logspline(delta.posterior,lbound=0) # NB. note the bound
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh3 # height of order--restricted prior at delta = 0
BF10 = prior/posterior #this produces the decimal version. if you inverse this you get the odds
BF01 = posterior/prior
cat("Bayes Factor(BF30) for H3 delta>0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF03) for H0 delta=0 over H3 delta>0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[3][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][3])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = 0
xhigh = 3
yhigh = 12
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(0,1,2,3,4), lab=c("0", "1", "2", "3", "4"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) 2*dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(0,1,2,3,4), lab=c("0", "1", "2", "3", "4"))
axis(2)
points(0, 2*dcauchy(0), pch=19, cex=2)
text(2,10, expression(H[0]: delta == 0),cex=2)
text(2,8, expression(H[3]: delta > 0),cex=2)
text(2,6, bquote(.(BFtext) == .(BFplot)),cex=2)
rm(list=ls(all=TRUE))
#---------------INPUT DATA------------------
setwd("/the/directory/for/both/of/your/files/") #this will help the program find the location of the modelthe model
group1 = c(1,5,4,2,1,2,3,1,2,3,1,2,5,4,7,8,9,7,5,3,2,3,4,6,7,6,4,7) #your set of values for group 1
group2 = c(5,4,3,3,4,6,6,8,7,6,5,7,7,6,5,6,7,5,4,5,6,7,8,5,6,7,8,9) #your set of values for group 2
openbugsmodel = "Ttest_2.txt"
priorforh1 = dcauchy(0) #prior for the case of delta<>0
priorforh2 = 2*dcauchy(0) #prior for the case of delta<0
priorforh3 = 2*dcauchy(0) #prior for the case of delta>0
#Advanced input
itterations = 30000 #Don't change these unless you know what you're doing
burnin = 3000 #Don't change these unless you know what you're doing
#-------------------------------------------
library(R2OpenBUGS)
n1 = length(group1)
n2 = length(group2)
# Rescale:
group2 = group2-mean(group1)
group2 = group2/sd(group1)
group1 = as.vector(scale(group1))
data=list("group1", "group2", "n1", "n2") # to be passed on to WinBUGS
inits=function()
{
list(delta=rnorm(1,0,1),mu=rnorm(1,0,1),sigma=runif(1,0,5))
}
# Parameters to be monitored
parameters=c("delta")
# The following command calls WinBUGS with specific options.
# For a detailed description see Sturtz, Ligges, & Gelman (2005).
samples = bugs(data, inits, parameters,
model.file =openbugsmodel,
n.chains=3, n.iter=itterations, n.burnin=burnin, n.thin=1,
DIC=T,
codaPkg=F, debug=F)
# Now the values for the monitored parameters are in the "samples" object,
# ready for inspection.
samples$summary # overview
# Please work through the analyses below one at a time
######################################################
# Analysis 1. H1: delta ~ Cauchy (unrestricted case)
######################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
#============ BFs based on logspline fit ===========================
library(polspline) # this package can be installed from within R
fit.posterior = logspline(delta.posterior)
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh1 # height of order--restricted prior at delta = 0
BF10 = prior/posterior
BF01 = posterior/prior
cat("-----------Results--------------\r\n")
cat("Bayes Factor(BF10) for H1 delta<>0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF01) for H0 delta=0 over H1 delta<>0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[1][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][1])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = -3
xhigh = 3
yhigh = 4
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(-4,-3,-2,-1,0,1,2,3,4), lab=c("-4","-3","-2","-1","0", "1", "2", "3", "4"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(-4,-3,-2,-1,0,1,2,3,4), lab=c("-4","-3","-2","-1","0", "1", "2", "3", "4"))
axis(2)
points(0, dcauchy(0), pch=19, cex=2)
text(-1,3.5, expression(H[0]: delta == 0),cex=2)
text(-1,3, expression(H[1]: delta != 0),cex=2)
text(-1,2.5, bquote(.(BFtext) == .(BFplot)),cex=2)
###########################################################################
# Analysis 2. H1: delta ~ Cauchy^- (restricted case, negative values only)
###########################################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
# selects only negative delta's:
delta.posterior = delta.posterior[which(delta.posterior<0)]
#============ BFs based on logspline fit ===========================
fit.posterior = logspline(delta.posterior,ubound=0) # NB. note the bound
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh2 # height of order--restricted prior at delta = 0
BF10 = prior/posterior
BF01 = posterior/prior
cat("Bayes Factor(BF20) for H2 delta<0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF02) for H0 delta=0 over H2 delta<0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[2][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][2])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = -3
xhigh = 0
yhigh = 12
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(-3,-2,-1,0), lab=c("-3","-2","-1","0"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) 2*dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(-3,-2,-1,0), lab=c("-3","-2","-1","0"))
axis(2)
points(0, 2*dcauchy(0), pch=19, cex=2)
text(-2,10, expression(H[0]: delta == 0),cex=2)
text(-2,8, expression(H[2]: delta < 0),cex=2)
text(-2,6, bquote(.(BFtext) == .(BFplot)),cex=2)
###########################################################################
# Analysis 3. H1: delta ~ Cauchy^+ (restricted case, positive values only)
###########################################################################
# Collect posterior samples across all chains:
delta.posterior = samples$sims.list$delta
# selects only positive delta's:
delta.posterior = delta.posterior[which(delta.posterior>0)]
#============ BFs based on logspline fit ===========================
fit.posterior = logspline(delta.posterior,lbound=0) # NB. note the bound
# 95% confidence interval:
x0=qlogspline(0.025,fit.posterior)
x1=qlogspline(0.975,fit.posterior)
posterior = dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
prior = priorforh3 # height of order--restricted prior at delta = 0
BF10 = prior/posterior #this produces the decimal version. if you inverse this you get the odds
BF01 = posterior/prior
cat("Bayes Factor(BF30) for H3 delta>0 over H0 delta=0: ",BF10,"\r\n")
cat("Bayes Factor(BF03) for H0 delta=0 over H3 delta>0: ",BF01,"\r\n")
cat("---\r\n")
if (BF10>=BF01){
BFplot=BF10
BFtext = bquote(BF[3][0])
}else{
BFplot=BF01
BFtext = bquote(BF[0][3])
}
BFplot=round(BFplot,2)
#============ Plot Prior and Posterior ===========================
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
xlow = 0
xhigh = 3
yhigh = 12
Nbreaks = 80
y = hist(delta.posterior, Nbreaks, prob=T, border="white", ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=2, lty=1, ylab="Density", xlab=" ", main=" ", axes=F)
#white makes the original histogram -- with unwanted vertical lines -- invisible
lines(c(y$breaks, max(y$breaks)), c(0,y$intensities,0), type="S", lwd=2, lty=1)
axis(1, at = c(0,1,2,3,4), lab=c("0", "1", "2", "3", "4"))
axis(2)
mtext(expression(delta), side=1, line = 2.8, cex=2)
#now bring in log spline density estimation:
par(new=T)
plot(fit.posterior, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lty=1, lwd=1, axes=F)
points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
# plot the prior:
par(new=T)
plot ( function( x ) 2*dcauchy( x, 0, 1 ), xlow, xhigh, ylim=c(0,yhigh), xlim=c(xlow,xhigh), lwd=1, lty=1, ylab=" ", xlab = " ", axes=F)
axis(1, at = c(0,1,2,3,4), lab=c("0", "1", "2", "3", "4"))
axis(2)
points(0, 2*dcauchy(0), pch=19, cex=2)
text(2,10, expression(H[0]: delta == 0),cex=2)
text(2,8, expression(H[3]: delta > 0),cex=2)
text(2,6, bquote(.(BFtext) == .(BFplot)),cex=2)
OpenBugs code
Ttest_2.txt
model
{
for (i in 1:n1)
{
group1[i] ~ dnorm(muX,lambdaXY)
}
for (i in 1:n2)
{
group2[i] ~ dnorm(muY,lambdaXY)
}
lambdaXY <- pow(sigmaXY,-2)
delta ~ dnorm(0,lambdaDelta)
lambdaDelta ~ dchisqr(1)
sigma ~ dnorm(0,sigmaChi)
sigmaChi ~ dchisqr(1)
sigmaXY <- abs(sigma)
mu ~ dnorm(0,muChi)
muChi ~ dchisqr(1)
alpha <- delta*sigmaXY
muX <- mu + alpha*0.5
muY <- mu - alpha*0.5
}
{
for (i in 1:n1)
{
group1[i] ~ dnorm(muX,lambdaXY)
}
for (i in 1:n2)
{
group2[i] ~ dnorm(muY,lambdaXY)
}
lambdaXY <- pow(sigmaXY,-2)
delta ~ dnorm(0,lambdaDelta)
lambdaDelta ~ dchisqr(1)
sigma ~ dnorm(0,sigmaChi)
sigmaChi ~ dchisqr(1)
sigmaXY <- abs(sigma)
mu ~ dnorm(0,muChi)
muChi ~ dchisqr(1)
alpha <- delta*sigmaXY
muX <- mu + alpha*0.5
muY <- mu - alpha*0.5
}
Additional information
The model specifications are:
Group 1 ~ Normal distribution ( μ + α/2, σ^2) #These are the only known values
Group 2 ~ Normal distribution ( μ - α/2, σ^2) #These are the only known values
σ ~ Cauchy distribution(0,1)+
μ ~ Cauchy distribution(0,1)
α = σ * δ
δ ~ Cauchy distribution(0,1) # This is assumed for H1 for the model while H0 is considered δ=0
Group 2 ~ Normal distribution ( μ - α/2, σ^2) #These are the only known values
σ ~ Cauchy distribution(0,1)+
μ ~ Cauchy distribution(0,1)
α = σ * δ
δ ~ Cauchy distribution(0,1) # This is assumed for H1 for the model while H0 is considered δ=0