A re-edit to my previous post. I tried in cross validated but was closed due to coding inquiry.
I'm currently trying to find the best way to analyze whether habitat restoration treatments affect arthropod family richness on a sand dune ecosystem and am having trouble determining the random effects and fixed effects in my model and the proper model to run that will allow be to do a post-hoc test.
Two treatments (herbicide and mechanical removal) and a control were conducted on 9 randomized blocks within a preserve. Each block contains three plots (ex., 1H, 1S, 1C -> Block 1; 2H, 2S, 2C -> Block 2). And each block/plot was sampled twice in two different sessions (sess), Spring or Summer to account for the temporal difference in insect type activity, although I'm not really interested in between session information. There are a total of 54 data points where I received different numbers for my Families response variable.
Spring session x 9 blocks x 3 plots each = 27 plots, Summer session x 9 blocks x 3 plots each = 27 plots
Other factors that could be affecting arthropod composition could be abiotic factors like side (Side) of the sand dune preserve (East side - at the top of the ridge with larger sand particles or West side - lower elevation and similar to ground level with smaller sand particles). Also aspect (Aspect) in reference to sun direction. Some of the plots are facing N, others facing S, E, or W.
With my understanding, the response variable would be family richness (Families), the independent variable would be treatment (Treat), and random/confounding factors would be Block, Plot, and Aspect since those were not controlled? I'm not too sure about whether Side being fixed or random in my case.
I think Plot could be nested in Block.
Fixed effects: Treatment; Random effects: Block, Plot, Aspect; Side is unknown but I am leaning towards fixed effect.
I tried to do a linear model using lmer and following that with a post-hoc test, but I'm not too sure if this is the right result because of my statistical interpretation of the fixed and random effects.
lme <-lmer(Families~Treat + (1 |blck/plt ) + (1|Aspect) + (1|Side), data=df)
summary(lme)
library(multcomp)
summary(glht(lme, linfct = mcp(Treat = "Tukey")), test = adjusted("holm"))
Any help on whether the input is a correct model, or if I'm understanding my variables incorrectly would be greatly appreciated. Also would this be a proper post-hoc test to look at differences?
Thank you.
I tried an Anova, Lmer, lm and I'm still a bit confused on how to add multiple random effects in my code and get a proper pairwise comparison.
I think you want
Families ~ Treat + Aspect + Side + (1 |blck/plt )
Aspect and Side to be fixed rather than random are (1) conceptually, (for example) N/S/E/W are not exchangeable values, and you're not trying to generalize to other aspects; (2) practically, random effects with only 2 (Side) or 4 (Aspect) levels are hard to fit.1|blck/plt/Treat); however, since you (IUIC) a single observation of each treatment in each plot, this variation would be confounded with the residual variation anyway, so 1|blck/plt is appropriate.((H+S)/2 - C)H-S)
If your factor is ordered C, H, S I think this would be specified in emmeans as list(avgtrt = c(-1, 1/2, 1/2), trtdiff = c(0, 1, -1)).The advantages of choosing meaningful contrasts in advance are (1) by convention, people don't expect you to do multiplicity corrections; (2) it is slightly harder to fall into the trap of overinterpreting a large number of significant/NS outcomes, although in this case there are only three pairwise comparisons anyway, so it might not make as much difference.