[R-sig-ME] [R] Comparing variance components
Jarrod Hadfield
j.hadfield at ed.ac.uk
Fri Feb 19 07:53:23 CET 2016
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>
More information about the R-sig-mixed-models
mailing list