Request for code used in ttestBF and BFttest.tstat in R
I'm new to Bayesian analysis and am finding it helpful to play with the BayesFactor package in R, in conjunction with reading the related articles. Thanks for these. I do want to be sure I understand what I'm doing, and it would be immensely helpful to be able to see the code for some of the functions.
Today I wrote a script to compare output from ttestBF and ttest.tstat, as I was confused by the differences, and it looks like the output from ttestBF is a raw odds ratio, whereas the output from ttest.tstat is a log odds ratio.
If I try to work out the odds ratios directly using the t distribution or Cauchy distribution (pt or pcauchy), I get values that are correlated with the Bayes factors from the package, but are very different in magnitude. I think this is almost certainly due to the fact that I am uncertain about how to incorporate the prior. If I could see how this is done in the package, this would be a great help to my understanding. I will append a short script (rmd) I have used to try and test this.
```{r generatedata}
#Generate 2 groups of 20, with mean difference m
makedat<-function(N,m){
mydat1<-rnorm((N/2),m,1)
mydat2<-rnorm((N/2),0,1)
myt<-t.test(mydat1,mydat2,paired=FALSE)$statistic
bf = ttestBF(mydat1,mydat2,mu=0,paired=FALSE,rscale="medium")
bf<-data.frame(bf)$bf # convert to data frame to access bf element
myresults<-c(myt,bf)
return(myresults)
}
```
## Create data frame to hold results for repeated runs
```{r dosums}
nrun<-10
sumdat<-data.frame(matrix(nrow=nrun,ncol=11))
colnames(sumdat)<-c('run','t','ttestBF','log.ttestBF','BFttest.tstat','p.t0','p.t1','logORt','p.cauchy0','p.cauchy1','logORcauchy')
N=40 #total N for both groups
m=.8 # true mean difference in population - can vary this
esteff=1 #effect size for H1 when computing BF
estcauch<-1
for (n in 1:nrun){
sumdat$run[n]<-n
myresults<-makedat(N,m)
sumdat$t[n]<-myresults[1]
sumdat$ttestBF[n]<-myresults[2]
sumdat$log.ttestBF[n]<-log(sumdat$ttestBF[n])
sumdat$BFttest.tstat[n]<-ttest.tstat(sumdat$t[n], N/2, N/2, rscale = 1,
complement = FALSE, simple = FALSE)$bf
sumdat$p.t0[n]<-1-pt(sumdat$t[n],df=(N-2),lower.tail=TRUE)
sumdat$p.t1[n]<-1-pt((esteff-sumdat$t[n]),df=(N-2),lower.tail=TRUE)
sumdat$logORt[n]<-log(sumdat$p.t1[n]/sumdat$p.t0[n])
sumdat$p.cauchy0[n]<-1-pcauchy(sumdat$t[n],lower.tail=TRUE)
sumdat$p.cauchy1[n]<-1-pcauchy((estcauch-sumdat$t[n]),lower.tail=TRUE)
sumdat$logORcauchy[n]<-log(sumdat$p.cauchy1[n]/sumdat$p.cauchy0[n])
}
plot(sumdat$log.ttestBF,sumdat$BFttest.tstat)
plot(sumdat$log.ttestBF,sumdat$logORt,pch=15)
lines(sumdat$log.ttestBF,sumdat$logORcauchy,type='p')
text(0,6,'Filled logOR.t \nand open logORcauchy')
```
Comments
Hi Dorothy,
All the code can be found on the github page.
`ttestBF()` can be found here: https://github.com/richarddmorey/BayesFactor/blob/master/pkg/BayesFactor/R/ttestBF.R
`ttest.tstat()` can be found here: https://github.com/richarddmorey/BayesFactor/blob/master/pkg/BayesFactor/R/ttest_tstat.R
The important "guts" of the `BayesFactor` t test functions can be found here (`meta.t.bf()` - a normal t test is treated as a meta-analysis with one t statistic): https://github.com/richarddmorey/BayesFactor/blob/master/pkg/BayesFactor/R/meta-ttest-utility.R
`ttestBF()` passes everything to `ttest.tstat()` which in turn passes everything to `meta.t.bf()`.
Basically, the "raw" functions like `ttest.stat` return a the logarithm of the Bayes factor because internally, that's how every Bayes factor is stored for numerical stability.
In order to see the relationship between the Cauchy prior and the Bayes factor, you can see the following code (this is similar to what you'd see in the `BayesFactor` code, but in the package a bunch of stability improvements have been made that make it more difficult to see the basics of what it is doing (e.g., corrections for huge t statistics, etc).
```
# Data setup
t_stat = 2
n1 = 10
n2 = 10
df = n1 + n2 - 2
effective_n = n1 * n2 / (n1 + n2)
# Prior
prior_scale = 1 # Cauchy prior scale
## Compute the marginal likelihood under the null, using Student's t
likelihood_null = dt(t_stat, df, ncp = 0)
## Compute the marginal likelihood under the alternative, using
## the noncentral t and Cauchy prior distribution over delta.
## Integrate over possible values of delta.
likelihood_alt =
integrate(
function(delta){
likelihood = dt(t_stat, df, ncp = delta * sqrt(effective_n))
prior = dcauchy(delta, scale = prior_scale)
return( likelihood * prior )
}, -Inf, Inf)[[1]]
# Bayes factor is the ratio of the two
bf_10 = like_alt / like_null
bf_10
## Do the same thing with BayesFactor package
## simple = TRUE yields a simpler version of the output
## with a Bayes factor (instead of log Bayes factor)
BayesFactor::ttest.tstat(t_stat, n1, n2, rscale = prior_scale, simple = TRUE)
```