I am trying to work with a dataset (see below). The Locality includes Rapsfeld and Bluhflache fields, which are not independent. I am setting SiteID as the random factor. I ran the test with both lmer and nlme packages, as this site explains: https://rcompanion.org/handbook/G_03.html. However, I am not sure if I am doing it correctly. Also, I've never used these tests before (only have experience with independent variable tests), so which values from which test would be most relevant for data analysis?
I'm just a little lost and definitely inexperienced, here is my code so far:
Input = ("
year SiteID Locality SpeciesRichness SpeciesAbund
2016 3 Bluhflache 34 103
2016 9 Bluhflache 26 107
2016 12 Bluhflache 44 202
2016 13 Bluhflache 24 77
2016 14 Bluhflache 28 72
2016 27 Bluhflache 21 58
2016 29 Bluhflache 36 234
2016 37 Bluhflache 18 77
2016 39 Bluhflache 26 113
2016 56 Bluhflache 33 150
2016 61 Bluhflache 22 75
2016 62 Bluhflache 17 45
2017 3 Bluhflache 20 73
2017 3 Rapsfeld 21 131
2017 9 Bluhflache 12 40
2017 9 Rapsfeld 22 62
2017 12 Bluhflache 15 50
2017 12 Rapsfeld 18 126
2017 13 Bluhflache 16 181
2017 13 Rapsfeld 20 124
2017 27 Bluhflache 12 21
2017 27 Rapsfeld 21 220
2017 29 Bluhflache 27 164
2017 29 Rapsfeld 29 348
2017 32 Bluhflache 16 30
2017 37 Bluhflache 9 66
2017 37 Rapsfeld 5 40
2017 39 Bluhflache 18 97
2017 39 Rapsfeld 24 205
2017 67 Bluhflache 6 19
2017 67 Rapsfeld 21 108
")
SpeciesData = read.table(textConnection(Input),header=TRUE)
Each SiteID could be thought of has a block including a measurement for the Rapsfeld and a measurement for the Bluhflache. (is this correct??)
Define model and conduct analysis of variance (with lmer)
library(lme4)
library(lmerTest)
model = lmer(SpeciesRichness ~ Locality + (1|SiteID), data=SpeciesData, REML=TRUE)
1|SiteID indicates that SiteID is the random term.
anova(model)
output:
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Locality 15.216 15.216 1 19.992 0.2355 0.6327
Test the random effects in the model (with lmer)
The rand function from the lmerTest package will test the random effects in the model.
rand(model)
output for 1|SiteID:
ANOVA-like table for random-effects: Single term deletions
Model:
SpeciesRichness ~ Locality + (1 | SiteID)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 4 -106.19 220.39
(1 | SiteID) 3 -106.44 218.89 0.50097 1 0.4791
Mixed model with nlme:
Define model and conduct analysis of deviance:
library(nlme)
model2 = lme(SpeciesRichness ~ Locality, random=~1|SiteID, data=SpeciesData, method="REML")
anova(model2)
Output:
numDF denDF F-value p-value
(Intercept) 1 16 153.80039 <.0001
Locality 1 16 0.23553 0.634
Test the random effects in the model
model.fixed = gls(SpeciesRichness ~ Locality, data=SpeciesData, method="REML")
anova(model2,
model.fixed)
output:
Model df AIC BIC logLik Test L.Ratio
model2 1 4 220.3882 225.8574 -106.1941
model.fixed 2 3 218.8892 222.9911 -106.4446 1 vs 2 0.5009693
p-value
model2
model.fixed 0.4791
