Hi Malcolm thanks for the suggestions regarding the priors. Changing the 7 to an 8 for the random effect priors worked a treat, as did your specification for uninformative priors. The model runs and results are sensible. In regards to your other questions about the data and model specification: I lied a bit in outlining the level of coverage. Some individuals have the 7 observations across each of 124 time periods, others less. This is partly because of different aged individuals and also because I ensured that I'm only analysing data from a given time point for an individual when all 7 response variables are present. I allow Age to randomly vary across months because the month term captures most of the extrinsic sources of variation in the model without specifying what this is (obviously this can be attributed when environmental terms are added to the model). It is highly likely thathow an individual responds to a change in the environment is age dependent, hence the random age slope on month (which just so happens to greatly improve model performance). Finally, I model growth as opposed to size at age as we only have information from the otolith on increment width which is a proxy for growth. I could do some form of back-calculation to convert increment measurements into fish size, but this process is not perfect and can introduce bias into the data. Also, I find model interpretation easier when I talk about growth (e.g. conditions in year x were good for growth, warm temperatures favour higher growth) rather than changes in fish size, which is very much conditional on the size of fish at the previous time step, and the range of fish that contribute to a given section of the chronology. Cheers John On Sat, Jan 24, 2015 at 5:49 AM, Malcolm Fairbrother <
M.Fairbrother at bristol.ac.uk> wrote:
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