Comparing variance components
Hi, Whether you should focus on differences in the intra-class correlation or the raw variances depends on the question. If it is suspected that the difference is due to differences in scale then I guess testing the IC would make more sense, because a change in scale would increase both variance components proportionally and the IC should remain the same. There is actually a way of just modelling scale differences in MCMCglmm using path analysis. It sounds odd but jut regress the response against itself in one of the experimental blocks. This works because y = y*beta + eta after a little rearangement becomes y = eta/(1-beta) and so 1/(1-beta) where beta is the path coefficient is an estimate of the scale (or 1/(1-beta)^2 the change in variance). For example, lets say Vg=Ve=1 in the first block but Vg=Ve=4 in the second (i.e. the scale is 2). G<-gl(50,6) g<-rnorm(50) Exp<-gl(2,150) y<-rnorm(300)+g[G] y[which(Exp==2)]<-2*y[which(Exp==2)] my_data<-data.frame(y=y, Exp=Exp, G=G) model.mcmc.c<-MCMCglmm(y~sir(~units:at.level(Exp,2),~units:at.level(Exp,2)), random=~G, rcov=~units, data=my_data) plot(1/(1-model.mcmc.c$Lambda)) Note that you can't just add y as a predictor in standard (g)lm(m) software because the likelihood is not in standard form. Regarding the variances, I usually use scaled F(1,1) priors in binomial models too. However, the long tail can be a bit of a problem in these models if the data set is small. In extreme cases it can lead to numerical problems: check the latent variables don't lie outside the range -25/25 for logit models. Reducing the scale in the prior can help. Cheers, Jarrod
On 17/02/2016 16:15, Dean Castillo wrote:
Hi Jarrod,
I have been trying to test a similar hypothesis for a while with only
limited success, so I wanted to thank you for your answer to Wen's
question. I did have a few additional questions of clarification, some
specific to my own data analysis.
For model.mcmc.a I think it is pretty straightforward to compare the
posterior distributions. For model.mcmc.b, now that you explicitly
modeled heterogenous variances for the residuals is it more
informative to examine the IC correlation for G for each Exp rather
than the variances themselves?
For my specific problem I have binomial data. I have been modeling the
data as proportions (bounded by 0-1 but can logit or arcsinsqrt
transform) as well as modeling it using "multinomial2" on the raw data.
The issue I am running into is that the variance of G for one of the
experimental blocks is very close to the boundary condition, while the
other is larger, when modeled as a proportion, and the HPD are very
wide. I have been using inv-gamma priors and will play around with the
parameter expanded priors as suggested in your course notes.
For the multinomial2 model should the priors for the variance be the
same? The posterior means are not as close to the boundary condition
(Exp1:G=0.8, Exp2:G=0.3) but the HPD are still very wide (1e-17,1.13)
for Exp2:G.
Any help is greatly appreciated
Dean
Date: Tue, 16 Feb 2016 20:59:33 +0000
From: Jarrod Hadfield <j.hadfield at ed.ac.uk
<mailto:j.hadfield at ed.ac.uk>>
To: Wen Huang <whuang.ustc at gmail.com
<mailto:whuang.ustc at gmail.com>>, "Thompson,Paul"
<Paul.Thompson at SanfordHealth.org>
Cc: r-sig-mixed-models at r-project.org
<mailto: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
<mailto: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
<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
<mailto: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
<mailto: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
<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 <mailto: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]]
>>>
>>> ______________________________________________
>>> R-help at r-project.org <mailto:R-help at r-project.org> mailing
list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide >>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>> _______________________________________________
>> 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 >> ----------------------------------------------------------------------- >> Confidentiality Notice: This e-mail message, including any attachments, >> is for the sole use of the intended recipient(s) and may contain >> privileged and confidential information. Any unauthorized review, use, >> disclosure or distribution is prohibited. If you are not the intended >> recipient, please contact the sender by reply e-mail and destroy >> all copies of the original message. > _______________________________________________ > 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 -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
-------------- next part -------------- An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20160219/c6290bdc/attachment-0001.pl>