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

[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

-- 
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