[R] MCMCglmm heteroscedasticity dependent on predictor

Ben Bolker bbolker at gmail.com
Thu Sep 15 14:47:48 CEST 2011

Atle Torvik Kristiansen <atletorvik <at> gmail.com> writes:

> I have a dataset where the residual variance decreases with on one of
> the predictors (population size).
> Currently, the full model looks like this:
> prior<-list(R=list(V=1e-16, nu=-2),G1=list(V=diag(2), nu=2))
> m<-MCMCglmm(response~poly(population size,2)*poly(other
> predictor,2)+time, random=~us(1+time):population, data=data,
> prior=prior)
> Basically, it's a random regression with multiple populations measured
> multiple times.
> I have limited knowledge of MCMC, so:
> 1)  Does the specification of the prior seem sensible?

  Reasonably so for the group variance,
but the residual variance looks a little funny:
why is the expected value (V) so close to zero, and
why is 'nu' negative? Is that a typo?

> 2) How do i specify rcov? Is e.g. rcov=~us(population size):units a
> good approach?

  As far as I can see, you can't specify a response of
residual variance to a _continuous_ covariate in MCMCglmm --
only to categorical (grouping) variables.  I could be wrong,
but a skim through the "CourseNotes" vignette doesn't show
any other kinds of examples.
> 3) If I would like to include the other predictor in the rcov
> specification. Is this a good approach, rcov=~us(other
> predictor:population size):units?

  I would worry about this after you deal with #2.
> I know I could easily do this in nlme, but I'm hoping to avoid it. One
> reason is that I understand MCMC methods make it straightforward to
> assess the relative contribution of each predictor to the response.

  Hmmm.  Maybe you could expand on that.  How do you want to do that?
It might be easier to find a way to do that assessment with the
results of nlme than to 

 Alternatively, you may be stuck learning either AD Model Builder
or WinBUGS (in which you can structure the model however you
want ... you might try the glmmBUGS package to get started ...

 I would strongly recommend that you send follow-ups to 
the r-sig-mixed-models at r-project.org mailing list (I would
redirect this myself if I weren't posting via Gmane ...)

> Atle Torvik Kristiansen

More information about the R-help mailing list