[R-sig-ME] multivariate mixed nested model

Thierry Onkelinx thierry.onkelinx at inbo.be
Thu Feb 2 10:12:57 CET 2017


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 op 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 op 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 op r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]



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