Answer Q1: To know if the model is correct we need more information. It depends on how you randomized the TreatNTreatN and TreatHTreatH combinations in each site (i.e. FactorSFactorS). If you assigned the TreatNTreatN and TreatHTreatH 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 locationdata = data)
because your replicate/block is nested within location.
Answer Q2: lmerlmer fits mixed-effect models and is a type of generalized linear mixed model with a guassianGaussian distribution. The function glmglm() can't fit random effects.
Answer Q3: You can't fit the model you specified using the glm()glm() function unless you treat FactorSFactorS as a fixed-effect,effect; you could use glmer()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 gaussianGaussian distribtion in glmerglmer, which is the same as doing the analysis in lmerlmer.
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?
