[R-sig-ME] comparing posterior means

[Doran, Harold]

Doug:

I grab the variance/covariance matrix of the random effects in a way that I think will make you cringe, but will share it here and am interested in learning how it can be done more efficiently. Keep in mind basically two principles. First, my lme.eiv function (which stands for linear mixed model error-in-variables) uses henderson’s equations as typically described and then stands on your shoulders and uses almost all of the functions in the Matrix package for sparse matrices.

BTW, though my function is not available in a package, I am happy to share it with you. It has complete technical documentation and is written using S3 methods with various common extractor functions typically used (and more). The function is intended to be used when the variables on the RHS are measured with error. Otherwise, my function simply matches lmer’s output.

So, let me use the following to represent the linear model as

Xb = y

Assume X is the leftmost matrix in henderson’s equation (so it is a big 2x2 blocked matrix), b is a vector holding both the fixed and random effects, and y is the outcome. The matrix X is big (and very sparse in my applications) and so I find its Cholesky decomposition and then simply solve the triangular systems until a solution is reached.

After a solution is reached (here is where you will cringe), I find the inverse of the big matrix X as

mat1.inv <- solve(L, I)

where L is the Cholesky factor of X and I is a conformable identity matrix. This, in essence, finds the inverse of the big matrix X, and then I can grab the variances and covariances of everything I need from this after the solution is reached.

Harold



From: dmbates at gmail.com [mailto:dmbates at gmail.com] On Behalf Of Douglas Bates
Sent: Wednesday, October 15, 2014 12:37 PM
To: Doran, Harold
Cc: Ben Pelzer; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] comparing posterior means

On Wed, Oct 15, 2014 at 9:30 AM, Doran, Harold <HDoran at air.org<mailto: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 http://arxiv.org/abs/1406.5823

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 sparsity.)

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> [mailto: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<mailto: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<mailto: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<mailto:R-sig-mixed-models at r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]