Bayesian Approach to Analyse Oneway RM ANOVA

Last post I analysed the fake data set using frequentists approach.  This time, I will use the Bayesian method to analyse the same data set by using Kruschkes code (which you may be able to download from Kruschkes blog).  The result will cross-examine against the frequentists approach to see if they show consistently result (supposedly we should get the same conclusion, because the difference is quite big).

First I loaded the data and did some pre-processes.  I modified Kruschkes code a bit, because my machine is operated by Linux.  Also remember to load rjags, so you can use BUGS and JAGS directly by calling it in R. 

tab14.3 <- read.table("tab14.3.txt", header=T)

#--------------------------------
# Bayesian RM-ANOVA
#--------------------------------
fileNameRoot="ANOVAonewayJagsWithinSubj" # for constructing output filenames
if ( .Platform$OS.type != "unix" ) { 
 X11 <- function( ... ) windows( ... ) 
}
require(rjags)

The if( .Platform@OS.type != "unix") statement checks the operating platform and makes sure that in the following codes, when the function, windows(), is issued, it will be converted to the function, X11(), because Linux uses X11(), instead of windows(), as the plotting tool.
   
Then, to write down the the model, I stored the content of the model into a variable, called modelstring, which is just a string variable in R. I listed all the important variables declaring in the model for the learning purpose.  It is very important to keep track of the vectors/matrices that are sent to BUGS from R.  This piece of information will be very helpful to debug especially when I was not very familiar with BUGS at the early stage of learning. 

#------------------------------------------------------------------------------
# THE MODEL.
# Key variables: Nx1Lvl, NSLvl, Ntotal, y, x1, a0, a1, aS, S, sigma, a1SDunabs, 

# aSSDunabs 
modelstring = "
model {
for ( i in 1:Ntotal ) {
 y[i] ~ dnorm( mu[i] , tau )
 mu[i] <- a0 + a1[x1[i]] + aS[S[i]]
}

# overall standard deviation
tau <- pow( sigma , -2 )
sigma ~ dunif(0, 10) # y values are assumed to be standardized

# grand mean
a0 ~ dnorm(0,0.001) # y values are assumed to be standardized

# independent variable
for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , a1tau ) }
a1tau <- 1 / pow( a1SD , 2 )
a1SD <- abs( a1SDunabs ) + .1
a1SDunabs ~ dt( 0 , 0.001 , 2 )

# subject variable
for ( jS in 1:NSLvl ) { aS[jS] ~ dnorm( 0.0 , aStau ) }
aStau <- 1 / pow( aSSD , 2 )
aSSD <- abs( aSSDunabs ) + .1
aSSDunabs ~ dt( 0 , 0.001 , 2 )
}
" # close quote for modelstring
writeLines(modelstring,con="modelonewayRM.txt")

The first for loop defines the lowest level of the model.  It declares quite a lot of variables, including i, Ntotal, y, mu, tau, a0, a1, x1, aS and S.  The "i" is the index for each row in the long format data frame.  In the case of this example data, the i goes from 1 to 45, which is 9 participants � 5 occasions of repeated measurements.  Note that i is a BUGS internal variable, so it is only valid inside BUGS code.  Also re-using this type of index variable inside BUGS code may not be a good idea.  

OK go back to the code.  Ntotal equals to 45.  The "y" variable with the i index represents each observational unit.  They are modelled by a normal density (dnorm) function with mean, mu and precision, tau.   BUGSs density function uses the parameter, precision, rather than standard deviation, to declare a normal distribution.  The standard deviation will be defined later.  Since the mu and tau have not been defined, the line comes after the ~ dnorm, defines what mu[i] is (each observational units mean is treated separately by the for loop).  This is the part you would usually see in a regression equation.  The symbol, <-, assigns the regression equation, beta0 (a0), beta1 (a1) � X1 (x1), and beta subject (aS)  � subject (S).  Beta 0 is the intercept, beta 1 and beta subject are slope for the treatment and subject factors, respectively.  In a one-factor RM-ANOVA, subject factor in fact is treated as a second factor with different subject treated as different level in that factor.  Thus, you could think of a one-way RM-ANOVA as a two-way factorial ANOVA without the interaction between the subject factor and the investigated factor.  The interaction cannot be modelled due to the fact that only one observation in each cell.   

After the modelled part was completed, the following lines outside the loop defines the unmodelled part, the precision, which was defined by using sigma, the standard deviation (i.e., the reversed 2 power function of the sigma is the precision).  Then the standard deviation was modelled by a uniform distribution with a range of 0 to 10.  Since the data to be feed into the Bayesian model will be standardised, the range of the sd is reasonable and relatively uninformative.

So now only three variables, a0, a1 and aS, need to be determined (y, x and S will be fed from the data loaded in R).  a0 is the grand mean and modelled by a normal density function with mean, 0 and standard deviation about 31.62 (equals to a precision of .001). 

a1 is the deflection (i.e., the differences between the grand mean and the means of each level of the independent variable).  j1 here represents each level of the tested independent variable, so it is from 1 to 5 (week 1 to week 5, Nx1Lv1=5) in this example.  The five deflections are all modelled by the same normal density, with mean, 0 and precision, a1tau, which later is defined by a truncated t distribution.  The "a1SDunabs" is firstly modelled by the t density function, with mean, 0, precision, .0001 and df, 2.  Then a1SDunabe is truncated by taking an absolute value, which then an .1 is added on it.  By doing this, the minimal value would not be very close to the zero.  The result is stored in a1SD, which then is assigned to a1tau by taking a reversed power of 2.

The subject factor basically was processed quite similar to a1.  Finally, the model description is sent to a txt file, called modelonewayRM.txt, which will be called by JAGS later. (to be continued)

* "Within-cell replications" refers to the repeated responses in a combination of the conditions. For example, in an experiment with 2 by 2 design, there are 4 cells. In a typical psychophysical experiment, one usually contributes a fixed number of responses for each cell, say 40. Thus, for this participant, a set of data will compose of 160 responses, 40 responses for each combination of condition (i.e., a cell).