[R-sig-ME] multivariate mixed nested model

[Claudio]

Thanks Paul,
I will try.
all the best
Claudio

Il giorno gio, 02/02/2017 alle 11.52 +0200, Paul Debes ha scritto:
> Hi Claudio
> 
> to make the model work in MCMCglmm, you need to explicitly specify
> the  
> interaction for each trait by adding "trait". This specification
> extends  
> to the dimensions of the prior(s) and distribution(s). If you wish
> to  
> specify covariances between traits (here for the residuals), your
> model  
> could be something like:
> 
> prior <- list(R = list(V = diag(3), nu = rep(0.002,3)))
> 
> m <- MCMCglmm(cbind(a,b,d) ~ -1 + trait*s,
> 	rcov = ~us(trait):units,
> 	family = rep("gaussian",3),
> 	data = df,
> 	prior = prior,
> 	verbose = FALSE,
> 	pl = TRUE)
> 
> summary(m)
> 
> Best,
> Paul
> 
> 
> 
>   On Thu, 02 Feb 2017 11:12:57 +0200, Thierry Onkelinx  
> <thierry.onkelinx at inbo.be> wrote:
> 
> > 
> > Dear Claudio,
> > 
> > It looks like you have only one observation (of each trait) per  
> > individual.
> > That is not sufficient to fit (0 + trait | individual). Hence the
> > error
> > message.
> > 
> > Using (1 | individuals) implies that each individual has a effect
> > which  
> > is
> > common for all traits. Whether that is a relevant assumption
> > depends on
> > your experiment. Tip: write the model as an equation for each
> > trait. Then
> > see if the random intercept for individual makes sense.
> > 
> > I'm not proficient with MCMCglmm. So I can help you with that.
> > 
> > Best regards,
> > 
> > 
> > 
> > ir. Thierry Onkelinx
> > Instituut voor natuur- en bosonderzoek / Research Institute for
> > Nature  
> > and
> > Forest
> > team Biometrie & Kwaliteitszorg / team Biometrics & Quality
> > Assurance
> > Kliniekstraat 25
> > 1070 Anderlecht
> > Belgium
> > 
> > To call in the statistician after the experiment is done may be no
> > more
> > than asking him to perform a post-mortem examination: he may be
> > able to  
> > say
> > what the experiment died of. ~ Sir Ronald Aylmer Fisher
> > The plural of anecdote is not data. ~ Roger Brinner
> > The combination of some data and an aching desire for an answer
> > does not
> > ensure that a reasonable answer can be extracted from a given body
> > of  
> > data.
> > ~ John Tukey
> > 
> > 2017-02-01 20:01 GMT+01:00 Claudio <oppela at gmail.com>:
> > 
> > > 
> > > Dear Thierry,
> > > thanks a lot, it is exactly the kind of suggestion I was looking
> > > for!
> > > However, when using (0 + traits|individuals) for the random
> > > effect of
> > > individuals I got the message:
> > > "Errore: number of observations (=1431) <= number of random
> > > effects
> > > (=1431) for term (0 + traits | individuals); the random-effects
> > > parameters and the residual variance (or scale parameter) are
> > > probably
> > > unidentifiable"
> > > The models work when I use for the individual part
> > > (1|individuals).
> > > 
> > > A second further question: in order to include in the model a
> > > continuous fixed covariate, am I doing the right thing when
> > > scripting:
> > > 
> > > mumo2 <- lmer(value~0 + traits + traits:species +
> > > traits:covariate (0 +
> > > traits|species:populations) + (1|individuals), mdata, REML=FALSE)
> > > 
> > > About cbind and MCMCglmm, below there is an example which causes
> > > the
> > > message:
> > > "Error in `[<-.data.frame`(`*tmp*`, , response.names, value =
> > > c(1, 2,
> > > 3,  : missing values are not allowed in subscripted assignments
> > > of data
> > > frames"
> > > 
> > > a = c(1, 2, 3, 5, 5, 7, 8, 9, 4, 5, 7, 8, 9, 4, 1, 2, 3, 5, 5, 4,
> > > 5, 7,
> > > 8, 9, 4, 1)
> > > b = c(8, 9, 4, 5, 7, 8, 9, 4, 1, 2, 3, 1, 2, 3, 5, 5, 7, 4, 1, 2,
> > > 3, 1,
> > > 2, 3, 6, 6)
> > > d = c(8, 9, 4, 5, 5, 7, 8, 9, 4, 5, 7, 8, 9, 4, 1, 2, 3, 2, 3, 1,
> > > 2, 3,
> > > 5, 5, 7, 6)
> > > s = c("m", "m", "f", "f", "m", "m", "m", "f", "m", "f", "f", "f",
> > > "m",
> > > "f", "f", "f", "m", "m", "m", "m", "m", "f", "m", "f", "f", "f")
> > > df = data.frame(a, b, d, s)
> > > 
> > > y<-cbind(a, b, d)
> > > prior <- list(R = list(V = 1, nu = 0.002))
> > > m <- MCMCglmm(y ~ s, family = "gaussian" , data = df, prior =
> > > prior,
> > > verbose = FALSE, pl = TRUE)
> > > summary(m)
> > > 
> > > all the best (and thanks again)
> > > Claudio
> > > 
> > > Il giorno lun, 30/01/2017 alle 16.02 +0100, Thierry Onkelinx ha
> > > scritto:
> > > > 
> > > > Dear Claudio,
> > > > 
> > > > I this you need to add the interaction with traits to all the
> > > > fixed
> > > > and random effects. Otherwise you assume that these have the
> > > > same
> > > > effect for each trait. Note that 0 + traits is identical to
> > > > traits -
> > > > 1.
> > > > 
> > > > mumo1 <- lmer(value~0 + traits + traits:species + (0 +
> > > > traits|species:populations) + (0 + traits|individuals), mdata,
> > > > REML=FALSE)
> > > > 
> > > > Your second question needs a reproducible example.
> > > > 
> > > > Best regards,
> > > > 
> > > > ir. Thierry Onkelinx
> > > > Instituut voor natuur- en bosonderzoek / Research Institute for
> > > > Nature and Forest
> > > > team Biometrie & Kwaliteitszorg / team Biometrics & Quality
> > > > Assurance
> > > > Kliniekstraat 25
> > > > 1070 Anderlecht
> > > > Belgium
> > > > 
> > > > To call in the statistician after the experiment is done may be
> > > > no
> > > > more than asking him to perform a post-mortem examination: he
> > > > may be
> > > > able to say what the experiment died of. ~ Sir Ronald Aylmer
> > > > Fisher
> > > > The plural of anecdote is not data. ~ Roger Brinner
> > > > The combination of some data and an aching desire for an answer
> > > > does
> > > > not ensure that a reasonable answer can be extracted from a
> > > > given
> > > > body of data. ~ John Tukey
> > > > 
> > > > 2017-01-28 16:43 GMT+01:00 Claudio <oppela at gmail.com>:
> > > > > 
> > > > > Hi all.
> > > > > I collected six body features (bf1-bf6)from three populations
> > > > > of a
> > > > > salamander and from two populations of another sister species
> > > > > of
> > > > > salamander.
> > > > > I would evaluate how the species (fixed) and population
> > > > > belonging
> > > > > (random) affect the body features, by comparing models built
> > > > > with
> > > > > lme4.
> > > > > For some models, I also want to include bf6 as covariate.
> > > > > Thus, in
> > > > > case
> > > > > of univariate analyses, some models, for example, could be:
> > > > > mo1<-lmer(bf1~species+(1|species:population), data,
> > > > > REML=FALSE)
> > > > > mo2<-lmer(bf1~species+bf6+(1|species:population), data,
> > > > > REML=FALSE)
> > > > > 
> > > > > However, I want to fit multivariate models, and my post is
> > > > > about
> > > > > this.
> > > > > First, I melted the data:
> > > > > mdata<-melt(data, id.vars = c("species", "population",
> > > > > "bf6"),
> > > > > measure.vars = c("bf1", "bf2","bf3","bf4","bf5"),
> > > > > variable.name =
> > > > > "traits)
> > > > > 
> > > > > Now the question.
> > > > > 1) Are the multivariate versions of the models mo1 and mo2
> > > > > above
> > > > > mumo1<-lmer(value~traits -1 + species +
> > > > > (1|species:populations) +
> > > > > (1|individuals), mdata, REML=FALSE)
> > > > > mumo1<-lmer(value~traits -1 + species + bf6 +
> > > > > (1|species:populations) +
> > > > > (1|individuals), mdata, REML=FALSE)
> > > > > 
> > > > > A secondary question, which in case I will move to a new
> > > > > post:
> > > > > it seemed to me that building multivariate models with
> > > > > MCMCglmm is
> > > > > easier. However, cbind did not work, even without missing
> > > > > values:
> > > > > to
> > > > > your knowledge, is there any issue?
> > > > > 
> > > > > thanks in advance
> > > > > Claudio
> > > > > 
> > > > > _______________________________________________
> > > > > R-sig-mixed-models at r-project.org mailing list
> > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > 	[[alternative HTML version deleted]]
> > 
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>