priors for multivariate mixed model in MCMCglmm with random intercepts and slopes - R-SIG-mixed-models

Fri, Jan 23, 2015 10:49 AM #

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