[R-sig-ME] Mean of random effects same as fixed effect?

[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:
> On Sun, Jan 8, 2012 at 1:55 AM, Jens Åström <jens.astrom at slu.se> wrote:
>> 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)
> 
> 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).
> 
>> 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
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models