[R-sig-ME] multivariate mixed nested model

Thu Feb 2 10:52:00 CET 2017

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

-- 
Paul V. Debes
Postdoctoral researcher

Email: paul.debes at utu.fi
Phone: +358 4659 11055