[R-sig-ME] Redirect from R-help Digest, Vol 103, Issue 15: [R] MCMCglmm heteroscedasticity dependent on predictor

[Ben Bolker]

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

> 
> ------------------------------
> 
> Message: 26
> Date: Thu, 15 Sep 2011 12:47:48 +0000
> From: Ben Bolker <bbolker at ...>
> To: <r-help at ...>
> Subject: Re: [R] MCMCglmm heteroscedasticity dependent on predictor
> Message-ID: <loom.20110915T143333-659 at ...>
> Content-Type: text/plain; charset="us-ascii"
> 
> 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 ... mailing list (I would
> redirect this myself if I weren't posting via Gmane ...)
> 
> > Atle Torvik Kristiansen
> > 
> >
> 
> Hello Ben,
> 
> I'm glad you replied.
> 
> 1) The residual prior is from Hadfield's course notes (Vignette) page
> 78, where he has specified a model with similar random effect. I haven't
> quite gotten to grips with the Inverse-Wishart, but I believe this is an
> improper flat prior. Not quite sure why I'm using it this way though?

  OK, from p. 26:

Although inverse-Wishart distributions with negative degree of belief 
parameters are not defined, the resulting posterior distribution can be 
defined if there is suficient replication. Specifying V=0 and n=-1 
is equivalent to a uniform prior for the standard deviation on the 
interval (0;1], and specifying V=0 and n=-2 is non-informative for a 
variance component.

  I can't say I've really come to grips with it either.
  You're setting a very uninformative prior for the residual
variance here, which seems fine (although I don't really see
why this specification is to be preferred to the default
(nu=0,V=1) ,,,
 
> 2) I've read through the course notes and overview, and I can't find
> anything on continous covariates either. I've also made several searches
> on Google with no result. However, I did some preliminary model
> selection on random effect structures, and end up with the same model,
> with or without the rcov specification. I am therefore guessing that
> MCMCglmm handles the residual heterogeneous variance in my dataset all
> by itself? I do end up with a very complex model with e.g. a cubic
> random effect, so I'm a bit concerned about overfitting.

   "handles" -- or "ignores".  I'm guessing the "rcov" is indeed
too complicated for you to fit, so things are just ending up
messing.  If you're lucky they're converging to zero.  Have you
looked at the trace plots?

> 3) Yes.
> 
> 4) Well, near the bottom of page 62 of the course notes, Hadfield
> describes a way of assessing the proportion of the total variance
> explained by one of the random effects (dams). It seems he divides the
> markov chain of the dams by all the other random effect markov chains.

  I think you mean the "variance of the dams" ?  I think that what
Hadfield is doing makes a bit more sense in the classical genetics
context -- partitioning genetic + environmental variance.  

> He also states somewhere else that in Bayesian models there is no
> technical difference between random and fixed effects. Thus, this
> approach should be applicable to a fixed effect also, i.e. doing this
> with model$Sol instead of model$VCV?  Seems reasonable to me, but I've
> been way wrong in the past.

  I think you might have to do it with mean-centered, squared values
of Sol.  I don't think this is necessarily going to be easy -- the
reason there's so much mess and controversy in the literature about
how to partition variance for multilevel models is that it's tricky
to define.

  On the other hand, I just poked around a bit and found

Gelman, Andrew, and Iain Pardoe. 2006. “Bayesian Measures of Explained
Variance and Pooling in Multilevel (Hierarchical) Models.”
Technometrics 48 (May):
241-251. doi:10.1198/004017005000000517. 
http://pubs.amstat.org/doi/abs/10.1198/004017005000000517.

which looks sensible and has R code at the end for application
to data from a WinBUGS run -- it might be 'straightforward'
(i.e. not necessarily easy but not conceptually difficult)
to adapt to MCMCglmm results.

> 
> Thank you for your help, it's been extremely helpful.

  I'm glad it's that kind of help and not the other kind :-)

> 
> Kind regards,
> 
> Atle Torvik Kristiansen
> 
>