Skip to content

heterogeneous residual variances using rcov in MCMCglmm?

2 messages · Malcolm Fairbrother, David Atkins

Original
Malcolm Fairbrother
# 
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
David Atkins
# 
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