I'm new in the forum and also I'm a beginner on a mixed effect model using lme4 package in R. After reading lot of questions and answer in different forums and also several tutorials and book chapters, I still have some doubts about my formulation of my model. I've done an experiment with two fixed effects and two levels each (Pollution: Yes/No, and Temperature: 17°/22°) and one random effect with 5 levels (Population: 1, 2, 3, 4, 5) that is nested on Pollution effect. For each combination I have 6 replicates. I'm interested on the effect of Pollution and Temperature on the fertilisation success. So I was considering the following model:

y = Pollution(Population) + Temperature + Pollution(Population)*Temperature + err

Hoping that my model includes Population as a random effect nested in Pollution and accounting for interaction with temperature I've been able to piece together the following:

model1<-lmer(rto_Embryos~ Pollution + (1|Pollution:Population) + Pollution*Temp + (1|Replicate), data=mydata, REML=F)

Is this accurately expressed? Any suggestions for how to improve my code if it does not correctly depicted?

Thank you very much in advance...

  • $\begingroup$ Whats your goal with the model? Do you have a specific question you are attempting to answer? $\endgroup$ – Matthew Drury Aug 21 '16 at 15:33
  • $\begingroup$ Please give more information about your experimental design. $\endgroup$ – Robert Long Aug 22 '16 at 9:39
  • $\begingroup$ Hi @MatthewDrury ... I'm trying to see if populations of a same species developed under polluted conditions are better coping with thermal acute stress at the embrionary stage. To that, I took samples of different populations living in two different conditions (Polluted vs Non-polluted), I took different animals and did crossed fertilization trials and then I incubated them undet two different temperatures (17º -Control, and 22º- treatment. As a response I counted the amount of fertilized eggs, embryos and larvae. $\endgroup$ – Gaston Aug 22 '16 at 15:14
  • $\begingroup$ So my intention is to see the effect of Pollution history on the performance of reproduction and larvae production under a acute thermal event. I don`t think a two-way ANOVA is useful, since I have fixed and random effect and nesting among some of the factors $\endgroup$ – Gaston Aug 22 '16 at 15:16
  • $\begingroup$ Hi @RobertLong, regarding your request, I have give some details in the lines above. But, to be more specific, I took 10 individuals of each population (5 in total, being 3 under polluted condition and 2 under non-polluted condition). 5 out of 10 individuals acted as a females and the remaining as males. We crossed those gametes and then splitted them in 6 petri dishes and kept them in the two different temperatures and we followed the development until settlement, looking at the perfomance of each stage as a response. Is this enough information?? $\endgroup$ – Gaston Aug 22 '16 at 15:23

I think what you're after is:

model1 <- lmer(rto_Embryos ~ Pollution * Temperature + (1 | Population), data = mydata, REML = F)

This syntax reflects the hierarchical nature of your data while enabling you to assess if the effect of temperature varies in presence of pollution. Keeping pollution in both, the random and the fixed portions of the model is incorrect. It seems you were trying to interact the fixed term with the random term but this is not how it works in lme4, check this and this for useful examples on lme4 syntax.

One aspect of your model that may cause problems, as Robert Long points out in the comments, is the low number of levels you have on the random effects. According to the mixed-models faq a minimum of 5 or 6 levels is required to have a reliable estimation of the variance. Conceptually, in my opinion, you should use a hierarchical model, but in practice you might run into issues with such low number of levels (i.e. your variance may be estimated to be 0 and/or the distribution of your random effects might be different from 0), if that's the case, including Population as a fixed term could be a viable workaround.

Also, for comparison with other models through likelihood ratio tests, it's ok to keep REML = F but you should set it to REML = TRUE for your final model.

Hope this helps

  • $\begingroup$ I'm not so sure about that. Is Replicate actually a grouping variable? Either way, with only 5 levels of Population, this isn't really sufficient for fitting random intercepts. $\endgroup$ – Robert Long Aug 25 '16 at 20:18
  • $\begingroup$ @RobertLong, thanks for that comment, I think you're right about replicates, the first time I read the question I assumed there was multiple observations per Replicate, but it seems more likely that each Replicate, or petri dish yields just one observation. I also added a cautionary note on the number of levels but I still think Population should be modeled as a random effect. $\endgroup$ – donlelek Aug 26 '16 at 14:41
  • $\begingroup$ @donlelek thank you very much for your reply! As you say, Replicate is a single observation made on a petri dish (6 in total per each treatment combination). $\endgroup$ – Gaston Aug 30 '16 at 8:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.