[R-sig-ME] heterogeneous residual variances using rcov in MCMCglmm?
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Oct 27 11:04:36 CEST 2011
Hi Malcolm,
I would be very careful about fitting heterogenous residuals to more
than a few
groups, especially 107. There's probably very little information to
estimate each variance (at least in BTdata). Quantitative geneticists
have become interested in this type of model, but they usually use a
"double-hierachcial" approach, where the 107 residual variances come
from some distribution (often log-normal). You can't fit these type of
model in MCMCglmm, but you could of course in JAGS/WinBUGS.
Lee/Nelder's (controversial?) HGLM approach can also be extended to
such an approach, and can even be implemented using existing REML
machinery (ASReml, and perhaps with some work, lme too). See
Ronnengard (2010) Genetics Selection Evolution 42 8.
Cheers,
Jarrod
Quoting Malcolm Fairbrother <m.fairbrother at bristol.ac.uk> on Thu, 27
Oct 2011 09:09:03 +0100:
> Hi Dave,
>
> Thanks, I was able to figure out a workable prior:
>
> prior3 <- list(R = list(V = diag(length(unique(BTdata$dam))), nu =
> 0.002), G = list(G1 = list(V = 1, nu = 0.002)))
> m3 <- MCMCglmm(tarsus ~ sex, random = ~dam, rcov=~idh(dam):units,
> prior = prior3, data = BTdata)
>
> Yes, this means m3$VCV has 107 columns (a variance for the dam
> level, and a residual variance for each of the 106 dams). I agree
> this seems unparsimonious, and for this kind of data is probably
> overkill. (Though the DIC does come out lower for this model, for
> whatever it's worth. And all the results seem sensible.)
>
> The paper I referenced discusses an application of mixed models to a
> very different kind of data: panel data on countries. In that realm,
> it seems more reasonable to think that level-2 units differ
> sufficiently such that level-1 variance may vary across them. The
> paper suggests that in the context of heterogeneous residual
> variances, the SEs for the fixed effects can be problematic, whereas
> allowing for heterogeneous residual variances deals with the problem.
>
> That's all IF I have understood correctly, however, which I may not
> have. Further input from anyone on the substance, not just the
> MCMCglmm arguments, would still be most welcome.
>
> Cheers,
> Malcolm
>
>
>
> On 26 Oct 2011, at 11:00, r-sig-mixed-models-request at r-project.org wrote:
>
>> Date: Tue, 25 Oct 2011 10:16:38 -0700
>> From: David Atkins <datkins at u.washington.edu>
>> To: r-sig-mixed-models at r-project.org
>> Subject: Re: [R-sig-ME] heterogeneous residual variances using rcov in
>> MCMCglmm?
>> Message-ID: <4EA6EEF6.3070608 at u.washington.edu>
>> Content-Type: text/plain; charset=windows-1252; format=flowed
>>
>>
>> Malcolm--
>>
>> I haven't used MCMCglmm for this previously, though I suspect your
>> intuition is correct re. rcov. However, you do need to provide a prior
>> for R that matches the dimensions of your (heterogeneous) residual variance.
>>
>> For example, it looks to me like the following fits heterogeneous
>> residual variances by *sex* (below).
>>
>> However, fitting a residual variance by *dam*, which has 106 unique
>> levels seems... unparsimonious? (Caveat: I didn't look at the paper you
>> referenced, so perhaps that is precisely what is suggested.)
>>
>> Hope that helps.
>>
>> cheers, Dave
>>
>>> prior2 <- list(R = list(V = diag(3), nu = 1.002), G = list(G1 =
>> list(V = 1, nu = 0.002)))
>>
>>> m2 <- MCMCglmm(tarsus ~ sex, random = ~dam, rcov=~idh(sex):units,
>> prior = prior2, data = BTdata)
>>
>> MCMC iteration = 0
>> [snip]
>> MCMC iteration = 13000
>>> summary(m2)
>>
>> Iterations = 3001:12991
>> Thinning interval = 10
>> Sample size = 1000
>>
>> DIC: 2015.288
>>
>> G-structure: ~dam
>>
>> post.mean l-95% CI u-95% CI eff.samp
>> dam 0.2554 0.166 0.3539 886.7
>>
>> R-structure: ~idh(sex):units
>>
>> post.mean l-95% CI u-95% CI eff.samp
>> Fem.units 0.5891 0.4970 0.6780 830.7
>> Male.units 0.6367 0.5430 0.7309 1000.0
>> UNK.units 0.5125 0.2885 0.7435 1000.0
>>
>> Location effects: tarsus ~ sex
>>
>> post.mean l-95% CI u-95% CI eff.samp pMCMC
>> (Intercept) -0.39849 -0.51948 -0.26687 1000 <0.001 ***
>> sexMale 0.76884 0.65555 0.87882 1000 <0.001 ***
>> sexUNK 0.15959 -0.06393 0.39012 1000 0.166
>> ---
>> Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>>
>> --
>> Dave Atkins, PhD
>> Research Associate Professor
>> Department of Psychiatry and Behavioral Science
>> University of Washington
>> datkins at u.washington.edu
>>
>> Center for the Study of Health and Risk Behaviors (CSHRB)
>> 1100 NE 45th Street, Suite 300
>> Seattle, WA 98105
>> 206-616-3879
>> http://depts.washington.edu/cshrb/
>> (Mon-Wed)
>>
>> Center for Healthcare Improvement, for Addictions, Mental Illness,
>> Medically Vulnerable Populations (CHAMMP)
>> 325 9th Avenue, 2HH-15
>> Box 359911
>> Seattle, WA 98104
>> http://www.chammp.org
>> (Thurs)
>>
>> Dear list,
>>
>> Am I right in thinking that the "rcov" term in a call to MCMCglmm can
>> allow for heterogeneous residual (level-1) variances, where the
>> heterogeneity is across level-2 units? If so, how do I make use of rcov
>> in this way?
>>
>> Consider the "BTdata" dataset included in the MCMCglmm package. The
>> following worked fine to estimate a two-level model, with observations
>> nested within dams, with a single residual variance across all dams:
>>
>> prior <- list(R = list(V = 1, nu = 0.002), G = list(G1 = list(V = 1, nu
>> = 0.002)))
>> m1 <- MCMCglmm(tarsus ~ sex, random = ~dam, prior = prior, data = BTdata)
>>
>> However, the variance of tarsus varies quite a bit by dam:
>>
>> summary(by(BTdata, BTdata$dam, function(x) var(x$tarsus)))
>>
>> (The variance is NA for one dam because there is only one observation
>> for that dam.) How therefore would I allow for heterogeneous residual
>> variances? I tried the following to estimate a unique residual variance
>> for each dam:
>>
>> m2 <- MCMCglmm(tarsus ~ sex, random = ~dam, rcov=~idh(dam):units, prior
>> = prior, data = BTdata)
>>
>> However, this fails and asks for a different set of priors, which I have
>> not been able to figure out. So two questions:
>> (1) Does this call to MCMCglmm make sense? (Or should I be trying some
>> other approach? I don't believe this can be done using lme4 and (RE)ML.)
>> (2) And how should I be specifying the priors for the MCMCglmm call?
>>
>> If you're wondering why I want to do this, I am following the advice
>> given in: http://pan.oxfordjournals.org/content/15/2/165.
>>
>> Any assistance would be much appreciated.
>>
>> Cheers,
>> Malcolm
>
> _______________________________________________
> 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.
More information about the R-sig-mixed-models
mailing list