Modelos jerárquicos para comparaciones múltiples: contexto de resultados múltiples

Acabo de (re) leer el Por qué de Gelman (generalmente) no tenemos que preocuparnos por las comparaciones múltiples . En particular, la sección "Resultados múltiples y otros desafíos" menciona el uso de un modelo jerárquico para situaciones en las que hay múltiples medidas relacionadas de la misma persona / unidad en diferentes momentos / condiciones. Parece tener una serie de propiedades deseables.

Entiendo que esto no es necesariamente algo bayesiano. ¿Podría alguien mostrarme cómo construir adecuadamente un modelo multivariado multinivel usando rjags y / o lmer (JAGS y BUGS normales también deberían estar bien, así como otras bibliotecas de modelos mixtos, por ejemplo, MCMCglmm) para que pueda jugar con él para comparar y resultados de contraste? El tipo de situación para la que me gustaría un modelo se refleja en los datos del juguete a continuación (medidas repetidas y multivariadas):

set.seed(69)
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?

Creo que tengo una solución parcial razonable para el modelo jerárquico bayesiano. rjagsCódigo a continuación ...

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
writeLines(modelstring,con="model.txt")

Nuevamente, el guión de medidas repetidas bayesianas de Kruschke

Finalmente encontré una solución de literatura para mi problema. Modelos bayesianos para múltiples resultados anidados en dominios por Thurston et al. 2009. Proponen un modelo jerárquico para dominios únicos o múltiples que refleja la naturaleza dependiente del dominio de las variables. Incorpora efectos aleatorios para individuos e individuos a través de dominios (si hay múltiples dominios). También se puede extender fácilmente para incluir medidas repetidas o diseños longitudinales.
Nota: publicaré un modelo JAGS aquí para completar la respuesta pronto