Solved – Nested mixed effects anova with lmer

generalized linear modellme4-nlme

enter image description here

I want to analyze some data from an experimental study with 3 factors:
treatN= nutrient manipulation 2levels N=nutrients added, A=ambient (fixed effect)

treatH= herbivory manipulation 3 levels C=ambient, H=moderate herbivory, HH=high herbivory (fixed effect)

Factor S= performed in 2 sites A and B(random effect)

I measured some response variables at the end of the experiment (area, necro…)
I want to perform a three way anova or any glm with mixed effects and with treat N and treatH nested in Factor S.

Q1:I used the package lmerTest, is this the right code? 

   m2 = lmer(area~ treatN*treatH + (1|FactorS) , data=data)

Q2:I assume that lmer is a glm and thus the errors don´t need to be normal? If not how can I test this?

Q3: I tried glm models 

    m4=glm(area ~ treatN*treatH + (1|FactorS), data=data, family=Gamma (link=identity))

but I get an error

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :contrasts can be applied only to factors with 2 or more levels
In addition: Warning message:
In Ops.factor(1, FactorS) : ‘|’ not meaningful for factors

how can I solve it?

Thanks!

Best Answer

Answer Q1: To know if the model is correct we need more information. It depends on how you randomized the TreatN and TreatH combinations in each site (i.e. FactorS). If you assigned the TreatN and TreatH combinations as a completely randomized design then I would say that yes your model is correct. If you randomized the the treatment combinations as a randomized complete block design I think the model should be: 

m2 = lmer(area ~ treatN*treatH + (1|FactorS/replicate), 
          data = data)

because your replicate/block is nested within location.

Answer Q2: lmer fits mixed-effect models and is a type of generalized linear mixed model with a Gaussian distribution. The function glm() can't fit random effects.

Answer Q3: You can't fit the model you specified using the glm() function unless you treat FactorS as a fixed-effect; you could use glmer() function doing the following:

m4 = glmer(area ~ treatN*treatH + (1|FactorS), 
           data = data, family = Gamma(link = "identity"))

If your data follows a normal distribution you can also use the Gaussian distribtion in glmer, which is the same as doing the analysis in lmer. 

m4 = glmer(area ~ treatN*treatH + (1|FactorS), 
           data = data, family = gaussian(link = "identity"))

(this will give you a warning saying you should just use lmer instead).

Here are a few links that can help you decide what distribution to use in your analysis:
When to use gamma GLMs?

Related Question