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

John Morrongiello jrmorrongiello at gmail.com
Sun Jan 18 05:56:16 CET 2015


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]]



More information about the R-sig-mixed-models mailing list