[R-sig-ME] heterogeneous residual variances using rcov in MCMCglmm?
David Atkins
datkins at u.washington.edu
Tue Oct 25 19:16:38 CEST 2011
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
More information about the R-sig-mixed-models
mailing list