Mean of random effects same as fixed effect?

Thanks for the input! It is appreciated. In practice, my problem is
solved but there still remains some questions.

I was shown off list the conventional way of modelling this in a
Bayesian framework (below). It seems I was barking in the right general
direction, but not at the right tree. Prof. Bates objections may explain
the slight differences in results between my hack version and the more
conventional, but I am not sure.

I realise now that simply taking the mean of observations defined to be
normally distributed is not the same as estimating the mean of that
underlying distribution. For instance, there may be extreme values that
could be given undue influence by just taking the mean.

I believe that the conventional BUGS approach deals with this problem,
but I am curious if it is still subject to any of Prof. Bates
objections. The unconditional means of the random effects is here
defined as the fixed intercept and slope values. Remember, the initial
problem was that doing differently resulted in non-convergence.

When I manipulate the sleepstudy data set to produce a non-balanced data
set, the three methods (my hack version, the conventional BUGS, and
lmer) all gives slightly different answers, otherwise pretty much the
same. This leaves me curious if one can be considered more correct than
the other.

Any other suggestions of model specifications or other thoughts are
appreciated, otherwise I will settle with the conventional for now.

For future reference, this appears to be a conventional way to implement
models of the form
fm1<-lmer(Reaction~Days+(Days|Subject),data=sleepstudy)
in BUGS/JAGS (thanks again Kent!). (Inverse wishart is another way)

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)

   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] <- intercept
     B.hat[j,2] <- Days.par
   }
   intercept ~ dnorm(0, 1.0e-6)
   Days.par ~ dnorm(0, 1.0e-6)

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

   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 #set this at zero for non correlated random effects
}

/Jens

Mean of random effects same as fixed effect?

Thread (3 messages)