Comparing variance components

Date: Tue, 16 Feb 2016 20:59:33 +0000
From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
To: Wen Huang <whuang.ustc at gmail.com>,  "Thompson,Paul"
        <Paul.Thompson at SanfordHealth.org>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] [R] Comparing variance components
Message-ID: <56C38DB5.2010709 at ed.ac.uk>
Content-Type: text/plain; charset=utf-8; format=flowed

Hi Wen,

The question sounds sensible to me, but you can't do what you want to do
in lmer because it does not allow heterogenous variances for the
residuals. You can do it in nlme:

model.lme.a<- lme(y~Exp, random=~1|G,  data=my_data)
model.lme.b<- lme(y~Exp, random=~0+Exp|G,
weights=varIdent(form=~1|Exp),  data=my_data)

or MCMCglmm (or asreml if you have it):

model.mcmc.a<- MCMCglmm(y~Exp, random=~G,  data=my_data)
model.mcmc.b<- MCMCglmm(y~Exp, random=~idh(Exp):G, rcov=~idh(Exp):units,
data=my_data)

The first model assumes common variances for each experiment, the second
allows the variances to differ. You can comapre model.lme.a and
model.lme.b using a likelihood ratio test (2 parameters) or you can
compare the posterior distributions in the Bayesian model.

Note that this assumes that the levels of the random effect differ in
the two epxeriments (and they have been given separate lables). If there
is overlap then an additional assumption of model.a is that the random
effects have a correlation of 1 between the two experiments when they
are associated with the same factor level.

Cheers,

Jarrod



On 16/02/2016 20:28, Wen Huang wrote:

Hi Paul,

Thank you. That is a neat idea. How would you implement that? Could you

write an example code on how the model should be fitted? Sorry for my
ignorance.

Thanks,
Wen

On Feb 16, 2016, at 1:18 PM, Thompson,Paul

<Paul.Thompson at SanfordHealth.org> wrote:

Are you computing two estimates of reliability and wishing to compare

them? One possible method is to set both into the same design, treat the
design effect (Exp 1, Exp 2) as a fixed effect, and compare them with a
standard F test.

-----Original Message-----
From: R-sig-mixed-models [mailto:

r-sig-mixed-models-bounces at r-project.org] On Behalf Of Wen Huang

Sent: Tuesday, February 16, 2016 11:57 AM
To: Doran, Harold
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] [R] Comparing variance components

Hi Harold,

Thank you for your input. I was not very clear. I wanted to compare the

sigma2_A?s from the same model fitted to two different data sets. The same
for sigma2_e?s. The motivation is when I did the same experiment at two
different times, whether the variance due to A (sigma2_A) is bigger at one
time versus another. The same for sigma2_e, whether the residual variance
is bigger for one experiment versus another.

Thanks,
Wen

On Feb 16, 2016, at 12:40 PM, Doran, Harold <HDoran at air.org> wrote:

(adding R mixed group). You actually do not want to do this test, and

there is no "shrinkage" here on these variances. First, there are
conditional variances and marginal variances in the mixed model. What you
are have below as "A" is the marginal variances of the random effects and
there is no shrinkage on these, per se.

The conditional means of the random effects have shrinkage and each

conditional mean (or BLUP) has a conditional variance.

Now, it seems very odd to want to compare the variance between A and

then what you have as sigma2_e, which is presumably the residual variance.
These are variances of two completely different things, so a test comparing
them seems strange, though I suppose some theoretical reason could exists
justifying it, I cannot imagine one though.

-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Wen
Huang
Sent: Tuesday, February 16, 2016 10:57 AM
To: r-help at r-project.org
Subject: [R] Comparing variance components

Dear R-help members,

Say I have two data sets collected at different times with the same
design. I fit a mixed model using in R using lmer

lmer(y ~ (1|A))

to these data sets and get two estimates of sigma2_A and sigma2_e

What would be a good way to compare sigma2_A and sigma2_e for these

two data sets and obtain a P value for the hypothesis that sigma2_A1 =
sigma2_A2? There is obvious shrinkage on these estimates, should I be
worried about the differential levels of shrinkage on these estimates and
how to account for that?

Thank you for your thoughts and inputs!



    [[alternative HTML version deleted]]