Skip to content

Mean of random effects same as fixed effect?

3 messages · Douglas Bates, Jens Åström

Original
Jens Åström
# 
Hi all,

A couple of weeks ago I posted a question but got no answers. Here goes
a second attempt, now shorter and more general.


Are the following two model specifications interchangeable, or is there
a statistical reason for why it is not OK to express model 1 in the form
of model 2?

Model 1)
y=fixed.intercept+fixed.slope*x+random.intercept+random.slope*x

Model 2)
y=random.intercept*x+random.slope*x
fixed.intercept=mean(random.intercept)
fixed.slope=mean(random.slope)



The reasons for my asking is that I have trouble getting convergence
with model specification 1, when the random intercepts and random slopes
are correlated, but specifying it as model 2 seemed to work. This is me
trying to implement some standard mixed models in BUGS/JAGS. Original
post with complete working example is here:
http://markmail.org/message/vhqeq4j3kldttlt5



I'm happy for any comments, with or without BUGS/JAGS code.

/Jens Astrom
1 day later
Douglas Bates
#
On Sun, Jan 8, 2012 at 1:55 AM, Jens ?str?m <jens.astrom at slu.se> wrote:
You would need at least equal group sizes and identical values of the
covariate with respect to which you have a random slope to be able to
count on this.  Even then I'm not entirely sure it would work.

Generally the unconditional distribution of the random effects is
defined to have a mean of zero.  I don't know how you are defining
yours (and prefer not to wade through BUGS/JAGS model specifications
to find out).
5 days later
Jens Åström
# 
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
On 01/09/2012 09:02 PM, Douglas Bates wrote: