[R-sig-ME] comparing posterior means

Douglas Bates bates at stat.wisc.edu
Wed Oct 15 18:36:57 CEST 2014

On Wed, Oct 15, 2014 at 9:30 AM, Doran, Harold <HDoran at air.org> wrote:

> Ben
> Yes, you can do this comparison of the conditional means using the
> variance of the linear combination AND there is in fact a covariance term
> between them. I do not believe that covariance term between BLUPs is
> available in lmer (I wrote my own mixed model function that does spit this
> out, however).
> Just to be didactic for a moment. Take a look at Henderson's equation(say
> at the link below)
> http://en.wikipedia.org/wiki/Mixed_model
> The covariance term between the blups that you would need comes from the
> lower right block of the leftmost matrix at the final solution. Lmer is not
> parameterized this way, so the comparison is not intended to show how that
> term would be extracted from lmer. Only to show that is does exist in the
> likelihood and can (conceivably) be extracted or computed from the terms
> given by lmer.

I would disagree, Harold, about the relationship between the formulation
used in lmer and that in Henderson's mixed model equations.  There is a
strong relationship, which is explicitly shown in

Also shown there is why the modifications from Henderson's formulation to
that in lmer lead to flexibility in model formulation and much greater
speed and stability in fitting such models.  Reversing the positions of the
random effects and fixed effects in the penalized least squares problem and
using a relative covariance factor instead of the covariance matrix allows
for the profiled log-likelihood or profiled REML criterion to be evaluated.
Furthermore, they allow for the sparse Cholesky decomposition to be used
effectively.  (Henderson's formulation does do as good a job of preserving

I believe you want the conditional variance-covariance matrix for the
random effects given the observed data and the parameter values.  The
sparse Cholesky factor L is the Cholesky factor of that
variance-covariance, up to the scale factor.  It is, in fact more stable to
work with the factor L than to try to evaluate the variance-covariance
matrix itself.

I'm happy to flesh this out in private correspondence if you wish.

> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:
> r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Pelzer
> Sent: Wednesday, October 15, 2014 8:56 AM
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] comparing posterior means
> Dear list,
> Suppose we have the following two-level null-model, for data from
> respondents (lowest level 1) living in countries (highest level 2):
> Y(ij) = b0j + eij = (b0 + u0j)  + eij
> b0j is the country-mean for country j
> b0 is the "grand mean"
> u0j is the deviation from the grand mean for country j, or the level-2
> residual eij is the level-1 residual
> The model is estimated by :  lmer(Y ~ 1+(1|country))
> My question is about comparing two particular posterior country-means.
> As for as I know, for a given country j, the posterior mean is equal to
> bb0 + uu0j, where bb0 is the estimate of b0 and uu0j is the posterior
> residual estimate of u0j.
> Two compare two particular posterior country means and test whether they
> differ significantly, would it be necessary to know the variance of
> bb0+uu0j for each of the two countries, or would it be sufficient to
> only know the variance of uu0j?
> The latter variance (of uu0j) can be extracted using
> rr <- ranef(modela, condVar=TRUE)
> attr(rr[[1]], "postVar")
> However, the variance of bb0+uu0j also depends on the variance of bb0
> and the covariance of bb0 and uu0j (if this covariance is not equal to
> zero, of course, which I don't know...).
> On the other hand, the difference between two posterior country means
> for country A and B say, would
> equal bb0 + u0A -(bb0 + u0B) = u0A - u0B meaning that I wouldn't need to
> worry about the variance of bb0.
> So my main question is about comparing and testing the difference
> between two posterior country means. Thanks for any help,
> Ben.
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

	[[alternative HTML version deleted]]

More information about the R-sig-mixed-models mailing list