[R-sig-ME] heterogeneous residual variances using rcov in MCMCglmm?

[Malcolm Fairbrother]

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