Search code examples
bayesiananovawinbugsjags

BUGS model for (nested?) repeated measures ANOVA


I was wondering if anyone has code for a BUGS/JAGS model for a repeated measures ANOVA? Basically, I have a response (y) that I want to model against Time of day, Day, and Treatment. I would also like to include two interaction terms, Treatment x Time of Day and Treatment x Day. There are about 20 individuals in the study, who were measured 4 times per day over about 1 week. I'm not entirely sure where to start, and I'm concerned that the Time of day covariate should also be nested within the Day covariate? If anyone has code for the likelihood portion of the BUGS/JAGS model, it would be greatly appreciated. I can take care of priors. Just can't seem to get off the ground with this one.


Solution

  • There are a few ambiguities in your question.

    1. Do you want Time of Day and Day to enter as continuous covariates or as discrete factors?

    2. Do you want individual identity to enter the model as a fixed or random effect?

    3. If either Day or Time of Day is a factor, do you want to include it as a fixed or random effect?

    4. You ask about whether Time of Day should be nested within Day. This is impossible to answer without knowing more about your data and your aims.

    Here's an example of code that assumes that you want to treat individuals as a random effect.

    Also assumed: Treatment, Time.of.day, and Day have constant slopes across all individuals. It would be straightforward to extend this model to a fixed- or random-slopes model where different individuals get separate modeled slopes. For example, for a random-slopes model, you'd just modify the beta parameters below to treat them in a manner similar to the alpha parameter.

    Following the OP's request, this is the likelihood portion only, and does not include the priors.

    for(i in 1:n.observations){
       y[i] ~ dnorm(alpha[individual[[i]] + beta1*Day[i] + beta2*Time.of.day[i] + beta3*Treatment[i] + beta4*Treatment[i]*Day[i] + beta5*Treatment[i]*Time.of.day[i], tau.obs)
    }
    # individual[i] contains the numerical index representing the individual that corresponds to observation i.
    
    for(j in 1:n.individuals){
       alpha[j] ~ dnorm(mu, tau)
    }