I've just been (re-)reading Gelman's Why we (usually) don't have to worry about multiple comparisons. In particular the section "Multiple outcomes and other challenges" mentions using a hierarchical model for situations when there are multiple related measures from the same person/unit at different times/conditions. It appears to have a number of desirable properties.

I understand that this is not necessarily a Bayesian thing. Could somebody show me how to properly construct a multivariate multilevel model using rjags and/or lmer (regular JAGS and BUGS should be fine too, as well as other mixed model libraries e.g., MCMCglmm) so that I can play around with it to compare and contrast results? The type of situation I would like a model for is reflected in the toy data below (multivariate, repeated measures):

id     <- factor(rep(1:20, 2))                # subject identifier
dv1    <- c(rnorm(20), rnorm(20,  0.8, 0.3))  # dependent variable 1 data for 2 conditions
dv2    <- c(rnorm(20), rnorm(20,  0.3, 0.6))
dv3    <- c(rnorm(20), rnorm(20, -0.3, 0.8))
dv4    <- c(rnorm(20), rnorm(20,  0.2, 1  ))
dv5    <- c(rnorm(20), rnorm(20,  0.5, 4  ))
rmFac  <- factor(rep(c(1, 2), each=20))       # repeated measures factor
dvFac  <- factor(rep(1:5, each=40))           # dependent variable indicator

dfwide <- data.frame(id, dv1, dv2, dv3, dv4, dv5, rmFac)
dflong <- data.frame(id, dv = c(dv1, dv2, dv3, dv4, dv5), rmFac, dvFac) # just in case it's easier?
  • $\begingroup$ For me it is not clear what your question is... I'm missing that question mark :) $\endgroup$ – Rasmus Bååth Nov 2 '12 at 9:56
  • $\begingroup$ @RasmusBååth I agree, I've edited it to hopefully make it clearer what I would like. Thanks. $\endgroup$ – Matt Albrecht Nov 2 '12 at 11:27

I think I've got a reasonable partial solution for the hierarchical Bayesian model. rjags Code below....

dflong$dv <- scale(dflong$dv)[,1]
dataList = list(  
    y = dflong$dv, 
    rmFac  = dflong$rmFac ,
    dvFac  = dflong$dvFac ,
    id     = dflong$id ,
    Ntotal = length(dflong$dv) ,
    NrmLvl = length(unique(dflong$rmFac)),
    Ndep   = length(unique(dflong$dvFac)),
    NsLvl  = length(unique(dflong$id))

modelstring = "
model {
for( i in 1:Ntotal ) {
    y[i] ~ dnorm( mu[i] , tau[rmFac[i], dvFac[i]])
    mu[i] <- a0[ dvFac[i] ] + aS[id[i], dvFac[i]] + a1[rmFac[i] , dvFac[i]]
for (k in 1:Ndep){
    for ( j in 1:NrmLvl ) { 
        tau[j, k] <- 1 / pow( sigma[j, k] , 2 )
        sigma[j, k] ~  dgamma(1.01005,0.1005)
for (k in 1:Ndep) {
    a0[k] ~ dnorm(0, 0.001)
    for (s in 1:NsLvl){
        aS[s, k] ~ dnorm(0.0, sTau[k])
    for (j in 1:NrmLvl) {
        a1[j, k] ~ dnorm(0, a1Tau[k])
    a1Tau[k] <- 1/ pow( a1SD[k] , 2)
    a1SD[k]  ~ dgamma(1.01005,0.1005)

    sTau[k] <- 1/ pow( sSD[k] , 2)
    sSD[k]  ~ dgamma(1.01005,0.1005)
" # close quote for modelstring

Again, base Bayesian repeated measures script from Kruschke


I finally found a literature solution to my problem Bayesian models for multiple outcomes nested in domains by Thurston et al. 2009. They propose a hierarchical model for single or multiple domains that reflects the domain dependent nature of the variables. It incorporates random effects for individuals and individuals across domains (if there are multiple domains). It can also be easily extended to include repeated measures or longitudinal designs.
Note: I'll post a JAGS model on here to complete the answer soon


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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