Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Supported by

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)

    ```

Sign In or Register to comment.