[R-sig-ME] Two correlated random effects in BUGS/JAGS

[Jens Åström]

Hi there modellers!

This is a question about how to specify a mixed model in BUGS/JAGS but
has application to mixed models in general. I've tried the BUGS list
without getting through, and I am confident several of the readers here
understands this better than me and should be able to answer.

Consider the example model in the lmer help with uncorrelated random
intercepts and slopes:
fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy)

I specified this as: (full working example further below)
##JAGS code##
mu[i]<-intercept+Days.par*Days[i]+subj.inter[subj[i]]+subj.Days[subj[i]]*Days[i]
intercept~dnorm(0.0,1.0E-06)
Days.par~dnorm(0.0,1.0E-06)
##

This seems to work fine as I get nice convergence and parameter
estimates that agree well with lmer. Note that I estimate both the fixed
effect of the intercept and slope explicitly, as well as the random
variations around those means (i.e. the random effects).

In contrast, the model with correlated random intercepts and slopes:
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

does not seem to work with such coding. Specifically, the fixed effects
estimates of the intercept and slope screw things up. No convergence
here (I hope I have tried everything). It does not seem to matter if I
use the multinomial normal distribution or the inverse Wishart for the
random effects either.

On the other hand, this seems to work:
##JAGS code##
mu[i] <-subj.inter[subj[i]] + subj.Days[subj[i]]*Days[i]
intercept<-mean(subj.inter[])
Days.par<-mean(subj.Days[])
##

This produces nice convergence and estimates that are close to those
from lmer (exept that the random effects estimates are not "scaled" as
deviations from the mean, i.e. the fixed effects). Note that I do not
estimate the fixed effects explicitly here, but calculate them as the
mean of the random effects. I am not sure this is OK. FYI, this does not
work in the uncorrelated example above.

That is the question. Is this correct procedure? It seems both plausible
and odd at the same time. Why the need for the different fixed effects
coding? Why should a more complex model estimate less parameters
(although the DIC estimates does not indicate that it does, oddly)?


IN SHORT: How would you code the model "fm1 <- lmer(Reaction ~ Days +
(Days|Subject), sleepstudy)" in BUGS/JAGS? This question extends to all
models with two correlated random effects I guess. It is also a question
whether the fixed effect can rightly be defined as the mean of the
random effects.


I would be grateful for any insights, suggestions or corrections.

/Jens Astrom



Full working example: (requires JAGS, produces two small text files)

##############################
require(lme4)
require(rjags)
load.module("glm") ##This module improves the performance of JAGS when
using glms and the like.

##Set up the data for JAGS
nr.subj<-length(levels(sleepstudy$Subject)) ## Number of levels of the
grouping factor Subject
subj<-sleepstudy$Subject ## Rename Subjects consecutively starting from 1
levels(subj)<-1:nr.subj
my.sleepstudy<-list("Reaction"=sleepstudy$Reaction,"Days"=sleepstudy$Days,"nr.subj"=nr.subj,"subj"=subj)
##list the variables that go in the JAGS model



fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),
sleepstudy) ##uncorrelated random intercepts and slopes
summary(fm2)

##Here is the JAGS version of fm2
cat(" ##Sleepstudy model with two, uncorrelated random effects. Simple
translation.
model{
for(i in 1:length(Reaction)){
Reaction[i]~dnorm(mu[i],tau.Reaction)
mu[i]<-intercept+Days.par*Days[i]+subj.inter[subj[i]]+subj.Days[subj[i]]*Days[i]
}

for(j in 1:nr.subj){
subj.inter[j]~dnorm(0.0,tau.subj.inter)
subj.Days[j]~dnorm(0.0,tau.subj.Days)
}

tau.subj.inter~dgamma(0.001,0.001)
sigma.subj.inter<-1/sqrt(tau.subj.inter)
tau.subj.Days~dgamma(0.001,0.001)
sigma.subj.Days<-1/sqrt(tau.subj.Days)
tau.Reaction~dgamma(0.001,0.001)
sigma.Reaction<-sqrt(1/tau.Reaction)

intercept~dnorm(0.0,1.0E-06)
Days.par~dnorm(0.0,1.0E-06)
}
",fill=T,file="fm2_simple.txt")



##Run this model in JAGS
sleep_uncorr<-jags.model("fm2_simple.txt",data=my.sleepstudy,n.chains=3)
update(sleep_uncorr,5000)
sleep_uncorr.res<-coda.samples(sleep_uncorr,c("intercept","Days.par","sigma.Reaction","sigma.subj.inter","sigma.subj.Days","subj.inter","subj.Days"),n.iter=20000,thin=10)
sleep_uncorr.par<-summary(sleep_uncorr.res)

sleep_uncorr.par ##These estimates agrees pretty well with lmer



##Now for the correlated random effects
fm1<-lmer(Reaction~Days+(Days|Subject),data=sleepstudy)

##JAGS version of fm1 (IS THIS TRICK OK?)
cat(" ##Sleepstudy model with two correlated random effects. Second attempt
model{

  for (i in 1:length(Reaction)){
    Reaction[i] ~ dnorm (mu[i], tau.Reaction)
    mu[i] <-subj.inter[subj[i]] + subj.Days[subj[i]]*Days[i]
  }
  tau.Reaction~dgamma(0.001,0.001)
  sigma.Reaction<-sqrt(1/tau.Reaction)

intercept<-mean(subj.inter[])
Days.par<-mean(subj.Days[])

  for (j in 1:nr.subj){
    subj.inter[j] <- B[j,1]
    subj.Days[j] <- B[j,2]
    B[j,1:2] ~ dmnorm (B.hat[j,], Tau.B[,])
    B.hat[j,1] <- mu.subj.inter
    B.hat[j,2] <- mu.subj.Days
  }
  mu.subj.inter ~ dnorm (0,1.0E-06)
  mu.subj.Days ~ dnorm (0,1.0E-06)

  Tau.B[1:2,1:2] <- inverse(Sigma.B[,])
  Sigma.B[1,1] <- pow(sigma.subj.inter, 2)
  sigma.subj.inter ~ dunif (0, 100)
  Sigma.B[2,2] <- pow(sigma.subj.Days, 2)
  sigma.subj.Days ~ dunif (0, 100)
  Sigma.B[1,2] <- rho*sigma.subj.inter*sigma.subj.Days
  Sigma.B[2,1] <- Sigma.B[1,2]
  rho ~ dunif (-1, 1)
  #rho<-0
}
",fill=T,file="fm1_simple2.txt")


##Now run this in JAGS

sleep_corr2<-jags.model("fm1_simple2.txt",data=my.sleepstudy,n.chains=3)
update(sleep_corr2,5000)
sleep_corr2.res<-coda.samples(sleep_corr2,c("intercept","Days.par","sigma.subj.inter","sigma.subj.Days","rho","subj.inter","subj.Days","sigma.Reaction"),n.iter=20000,thin=10)
sleep_corr2.par<-summary(sleep_corr2.res)

sleep_corr2.par ###These estimates are close to those from lmer. Is this
the way to specify the fixed effects when you have correlated random
effects?

#########################################