I'm a neuroscientist trying to move away from frequentist to bayesian statistics, please bear with me...

I'm after a hypothesis test on some of my data: e.g., let's say I have reaction times for two conditions from 10 subjects, 100 trials per subject. All participants do both conditions, i.e. it is repeated measures.

I want to get a test for the difference between the two conditions. Normally, I would average reaction times for each condition in each subject and then do a paired (1-sample) t-test across difference between the two condition means in each subject.

With a Bayesian approach I could also average across conditions and then use the BayesFactor R Package to get a BF or Kruschke's BEST R package to get a region of practical equivalence (ROPE)...

However, I find the preceding step of averaging across the trials in each participant unsatisfactory for a Bayesian analysis: For example, if I just wanted a posterior measure of the response times for each condition across the group I'd do a hierarchical analsyis (using JAGS) and include all the data (data -> individual priors -> group prior -> hyperprior on group -> fixed hyperprior on the hyperprior). This way I can get posterior estimates of the group response time for each condition. However, I don't know how to 'test' if these distributions are different. I believe that this is made more complicated by the fact that the response time between conditions is correlated at the level of subjects (i.e. some are just naturally faster regardless of condition than others).

My question: is there a way to calculate a Bayes Factor at the group level using a hierarhical model using all the (correlated) data? Or, am I just over complicating things and should just use the methods I already know which average across trials as mentioned above.

I hope that makes sense.

Many Thanks!

PS. I'm aware that some people are not a fan of Bayes Factors, and I'm certainly open to other tests, however, personally they are a useful statistic for quick digestion of a result in a paper. I also know that some people are sceptical of testing a null of no effect, but again, I find this the most intuitive way to think about it - open to more sophisticated priors.

  • 1
    $\begingroup$ To model these data properly is straight forward but not trivial. You'd want a hierarchical model. At the lowest level, each subject's data are described by some suitable pair of distributions for the two conditions; perhaps each is a Weibull. The distribution for the control condition would have a subject-specific central tendency parameter, and the central tendency of the other condition would be expressed by a subject-specific deflection from the baseline cen tend. Then the distribution of central tendencies would be under a higher-level distribution, and the distribution of deflections. $\endgroup$ Dec 14, 2016 at 20:45

1 Answer 1


Here's quick outline of a JAGS model. It's not implemented or tested, but I think the logical structure is right. Feel free to edit, of course.

Suppose that the data file is in "long" format. The row index is denoted rowIdx. Columns are the following:

y[rowIdx] = the response value for a trial

subj[rowIdx] = the subject index; subjects assumed to be consecutively numbered 1:nSubj

treatmentGroup[rowIdx] = 0 or 1 code for control group or treatment group

Then the core of the JAGS model could look something like this:

for ( rowIdx in 1:nRowsInDataFile ) {
  y[rowIdx] ~ dnorm( muControl[subj[rowIdx]] 
                     + treatmentGroup[rowIdx] *   treatmentEffect[subj[rowIdx]] , 
                     1/sigma^2 )
for ( subjIdx in 1:nSubj ) {
  muControl[subjIdx] ~ dnorm( muControlOverall , 1/sigmaMu^2 )
  treatmentEffect[subjIdx] ~ dnorm( treatmentEffectOverall , 1/sigmaEff^2 )
# prior on top level parameters:
muControlOverall ~ dnorm( 0 , 1/100^2 ) # or whatever
treatmentEffOverall ~ dnorm( 0 , 1/100^2 ) # or whatever
sigma ~ dunif( 0 , 1000 ) # or gamma
sigmaMu ~ dunif( 0 , 1000 ) # prefer gamma
sigmaEff ~ dunif( 0 , 1000 ) # prefer gamma
  • 1
    $\begingroup$ Above assumes normal distribution for y merely for convenience. If y is RT data, you'll want a skewed distribution... $\endgroup$ Dec 19, 2016 at 23:28
  • $\begingroup$ This is great - thanks! I've been playing around with the code today and it seems to be giving reasonable estimates on some simulated data so I think that it's working! For anyone else reading this and interested in this kind of problem - buy John's book, I did about a year ago it's fantastic, highly recommended! $\endgroup$
    – B. Critt
    Dec 20, 2016 at 18:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.