Dr. Kruschke,
We have a question about how to extend the LogisticOneWayAnovaJagsSTZ example from one to two factors. The difficulty is in specifying the correct nested array for a and b in dbeta(). Below is the pertinent code for the oneway case. We have tried many things, but none has worked. It would be greatly appreciated if you could explain how to expand it for two or more variables. Thank you. for ( i in 1:Ntotal ) { z[i] ~ dbin( theta[i] , n[i] ) theta[i] ~ dbeta( aBeta[x[i]] , bBeta[x[i]] ) T(0.001,0.999) } 
Administrator

Here's the idea. This is for a completely betweensubjects design: Each subject provides data for only one cell in a twofactor design. For subject i, there are n[i] trials and z[i] "successes". The data file consists of four columns, with one row per subject: z[i] n[i] x1[i] x2[i] (optional, unused: subject ID)      27 40 1 1 1 24 42 1 1 2 ... 35 45 1 2 19 36 40 1 2 20 ... 28 38 2 1 50 27 39 2 1 51 ... The hierarchical model estimates a probability of success for each individual, theta[i]. The individual probabilities in each cell come from a beta distribution with mean mu[x1[i],x2[i]] specific to that cell. The precision of the beta distribution is denoted kappa. For simplicity it is assumed to be the same for all cells. This is like the assumption of homogeneous variances in classical ANOVA. It can be relaxed, of course. The mean of each cell is a linear combination of baseline, additive factor deflections, and leftover interaction deflections, as in twofactor "ANOVA". Then the model statement could be something like the following. I have NOT tried this, so it may have minor glitches, but the gist is right. model { for ( i in 1:Ntotal ) { z[i] ~ dbin( theta[i] , n[i] ) theta[i] ~ dbeta( mu[x1[i],x2[i]] * kappa , (1mu[x1[i],x2[i]]) * kappa )T(0.001,0.999) } for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) { mu[j1,j2] < 1 / ( 1 + exp( ( a0 + a1[j1] + a2[j2] + a1a2[j1,j2] ) ) ) } } kappa ~ dgamma( 1.0 , 0.01 ) # or try dunif(2,1000) or whatever # The following is the usual prior on twofactor anova components: a0 ~ dnorm( 0 , 0.001 ) for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , 1/a1sd^2 ) } for ( j2 in 1:Nx2Lvl ) { a2[j1] ~ dnorm( 0.0 , 1/a2sd^2 ) } for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) { a1a2[j1,j2] ~ dnorm( 0.0 , 1/a1a2sd^2 ) } } } a1sd ~ dgamma(...,...) # or folded t or dunif a2sd ~ dgamma(...,...) # or folded t or dunif a1a2sd ~ dgamma(...,...) # or folded t or dunif # Convert a0,a1[],a2[],a1a2[] to sumtozero: ** SEE ANALOGOUS CODE IN TWO FACTOR ANOVA PROGRAMS (SEE BLOG) ** } Hope that helps! On Wed, Jul 3, 2013 at 9:51 AM, Shawn [via Doing Bayesian Data Analysis] <[hidden email]> wrote: Dr. Kruschke, 
Thank you very much for the lengthy response. We tried this, and its partial failure prompted the posting to the forum. The oneway versions with both factors give quite clean results. When the second factor is added in, something goes wrong and we end up with highly peaked marginals that are ridiculously wide, and the "main effect" estimates are off compared to the oneway cases. We will try a few other things and report back if a solution is found.

Administrator

Your description of the marginal distributions sounds like you are getting hypershrinkage. This is not "wrong", it's simply a consequence of the assumed hierarchical structure. When the withincell variance is very large compared to the betweencell variance, a good description is to estimate the cell means to be near each other. One way to check your program is to use fictitious data that have small withincell variance and relatively large betweencell variance. Be sure that there IS withincell variance so that the parameter estimates don't collapse. The program should recover the generating parameters well. On Fri, Jul 5, 2013 at 10:47 AM, Shawn [via Doing Bayesian Data Analysis] <[hidden email]> wrote: Thank you very much for the lengthy response. We tried this, and its partial failure prompted the posting to the forum. The oneway versions with both factors give quite clean results. When the second factor is added in, something goes wrong and we end up with highly peaked marginals that are ridiculously wide, and the "main effect" estimates are off compared to the oneway cases. We will try a few other things and report back if a solution is found. 
Thanks for this tidbit of useful information. After looking at your code more closely, I realized there was one small difference from our code pertaining to nested indexing that got us into trouble. Problem fixed.
Do you know of a good way to check the accuracy of the model created by JAGS? That is, is there a call that can return all node assignments? Such validation would be helpful for debugging. I have not found any such functionality yet in searching WinBUGS documentation. 
Administrator

Well, from the rjags manual, I find these functions that might be useful. First, compile the model, for example like this: > jagsModel = jags.model( "model.txt" , data=dataList , inits=initsList , + n.chains=nChains , n.adapt=adaptSteps ) Then you can examine the resulting jagsModel object using commands like these: variable.names(jagsModel) On Fri, Jul 5, 2013 at 1:39 PM, Shawn [via Doing Bayesian Data Analysis] <[hidden email]> wrote: Thanks for this tidbit of useful information. After looking at your code more closely, I realized there was one small difference from our code pertaining to nested indexing that got us into trouble. Problem fixed. 
Dr Kruschke,
I have one more question about twoway anovas. If I have a repeated measures design (still logistic), should I switch to the specification listed on p.434 in your book, or should be preceding model be expanded. If the latter, could you please describe how to include each participant in the mu[] model specification. Thanks for all of your input. The information provided by the posteriors is quite enlightening. 
Free forum by Nabble  Edit this page 