[R-sig-ME] Fwd: priors for multivariate mixed model in MCMCglmm with random intercepts and slopes

[Malcolm Fairbrother]

Hi John,

This is not my substantive area of expertise at all, and I'm not completely
confident I can help. But, nobody else has responded, and I have fitted
some multivariate models with MCMCglmm. So, for what it's worth...

For starters, I'm not sure I understand the data. If you're saying each of
7 properties are observed on 38 fish each observed over the course of 124
months, where does 1021 come from? (38x124 is 4712...) Are you saying you
have repeated observations on fish over the course of 124 months, but for
each actual fish you have fewer observations than 124?

Next, two questions about your models:
Why do you allow the slope for Age to vary randomly across months?
And why model growth (so, say, % or absolute change in length) rather than
the property that is growing (length)? The latter would seem more
straightforward to me.

As regards the priors, if I understand you correctly, you've kept the same
prior specification (prior1) when fitting models (m3 and m4) with an
additional parameter at each higher level. So maybe try changing 7 to 8 at
those levels?

prior2  <- list(R=list(V=diag(7), nu=7), G=list(G1=list(V=diag(8), nu=8),
G2=list(V=diag(8), nu=8)))

Additionally, also for what it's worth, I've found parameter-expanded
priors to be more uninformative, and Jarrod has made them pretty easy to
use. So (assuming you want them to be uninformative) you might try:

prior2a  <- list(R=list(V=diag(7), nu=7.02), G=list(G1=list(V=diag(8),
nu=8.02, alpha.mu=rep(0,8), alpha.V=1000*diag(8)), G2=list(V=diag(8),
nu=8.02, alpha.mu=rep(0,8), alpha.V=1000*diag(8))))

Hope that helps...?

Don't feel bad. You're far from the first person to write to this list with
questions about priors for MCMCglmm!

Cheers,
Malcolm




> Date: Thu, 22 Jan 2015 22:23:26 +1100
> From: John Morrongiello <jrmorrongiello at gmail.com>
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] Fwd: priors for multivariate mixed model in
>         MCMCglmm with random intercepts and slopes
>
> Just wondering if anyone has any thoughts on this question about MCMCglmm
> priors for a multivariate mixed model with random intercepts and slopes?
>
> ---------- Forwarded message ----------
> From: John Morrongiello <jrmorrongiello at gmail.com>
> Date: Sun, Jan 18, 2015 at 3:56 PM
> Subject: priors for multivariate mixed model in MCMCglmm with random
> intercepts and slopes
> To: r-sig-mixed-models at r-project.org
>
> Hi
> I've got a data set with 7 response variables measured monthly through time
> across 38 individuals (FishID). There are 124 months (MonthCode)
> representing each sampling time, with the data set having 1021 observations
> in total. I have previously analysed this data using a series of univariate
> mixed models fitted in lme4 of the form
>
> m1<-lmer(AvGrowth~ Age + (Age|FishID) + (Age|MonthCode), data=alldata1)
> .....mx<-lmer(.....)
>
> Here, Age is a continuous covariate that models an age-dependent decline in
> growth; the slope of this relationship allowed to vary amongst individuals
> (FishID) and through time (MonthCode). To this model structure I have also
> added environmental effects like temperature etc.
>
> As the 7 response variables are all measured at the same time from each
> fish (they are estimates of growth and otolith microchemistry), I thought
> to fit a multivariate mixed model to estimate covariances and also
> succinctly ascertain the importance of various environmental effects.
>
> I have had success fitting m2 below where the overall model intercept is
> suppressed and the trait-dependent iAge effect is allowed to vary by
> individual and MonthCode (similar to m1):
>
>
> prior1<-list(R=list(V=diag(7),nu=7),G=list(G1=list(V=diag(7),nu=7),G2=list(V=diag(7),nu=7)))
>
>
> m2<-MCMCglmm(cbind(AvGrowth.std,NaCa.std,SrCa.std,logMgCa.std,logBca.std,logBaCa.std,logLiCa.std)
> ~(trait + trait:Age -1),
> random=~us(trait:Age):FishID + us(trait:Age):MonthCode,
> rcov=~us(trait):units,
> family=rep("gaussian",7), prior=prior1,
> nitt=60000,thin=25,burnin=10000, data=alldata1,verbose=FALSE)
>
> When I use posterior.mode(m2$VCV) however, I don't get estimates of the
> random intercept terms FishID and MonthCode, nor their covariance with Age,
> rather just variances dependent on the iAge slope (e.g.
> AvGrowth.std:Age):AvGrowth.std:Age).FishID). I therefore tried to
> explicitly code a correlated random intercept and slope using us(1+
> trait:(Age):FishID in m3:
>
>
> m3<-MCMCglmm(cbind(AvGrowth.std,NaCa.std,SrCa.std,logMgCa.std,logBca.std,logBaCa.std,logLiCa.std)
> ~(trait + trait:Age -1),
> random=~us(1+ trait:Age):FishID + us(1+ trait:Age):MonthCode,
> rcov=~us(trait):units,
> family=rep("gaussian",7), prior=prior1,
> nitt=60000,thin=25,burnin=10000, data=alldata1,verbose=FALSE)
>
> This returns the error:
> Error in priorformat(if (NOpriorG) { :  V is the wrong dimension for some
> prior$G/prior$R elements
>
> So I'm tipping something is wrong with my priors. m4, where I estimate an
> overall model intercept and specific random intercepts also returns the
> same error.
>
>
> m4<-MCMCglmm(cbind(AvGrowth.std,NaCa.std,SrCa.std,logMgCa.std,logBca.std,logBaCa.std,logLiCa.std)
> ~(trait + trait:Age +1),
> random=~us(1+ trait:Age):FishID + us(1+ trait:Age):MonthCode,
> rcov=~us(trait):units,
> family=rep("gaussian",7), prior=prior1,
> nitt=60000,thin=25,burnin=10000, data=alldata1,verbose=FALSE)
>
> Would someone have some tips about how to get the priors for m3 and m4
> working? I've worked through the course notes and read some email posts,
> but admit that I a little at sea in terms of understanding the specifics of
> what is being coded on the prior side of things.
>
> Cheers
>
> John
>

	[[alternative HTML version deleted]]