I want to test the effects of island area and land use, and the interaction between island area and land use on species richness. For land use, I have three groups, namely forest, farmland and mix. The data is based on transects on different islands, so the island ID is set as random effect.
My model looks like this:
#model = glmer(SR ~ Area + land_use + Area:land_use + (1|islandID))
#summary(model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: SR ~ Area + land_use + Area:land_use + (1 | islandID)
Data: transect_ZS
REML criterion at convergence: 184.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.66105 -0.56159 -0.00294 0.57259 1.72096
Random effects:
Groups Name Variance Std.Dev.
islandID (Intercept) 0.1524 0.3903
Residual 0.6805 0.8249
Number of obs: 70, groups: islandID, 34
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.9996 0.5187 57.0061 -1.927 0.05893 .
Area 0.9064 0.2834 40.9977 3.198 0.00267 **
land_useforest 0.6563 0.5569 62.0889 1.179 0.24309
land_usemix 0.9611 0.6373 55.3032 1.508 0.13723
Area:land_useforest -0.8318 0.3034 63.4045 -2.742 0.00793 **
Area:land_usemix -0.7756 0.4748 56.3692 -1.633 0.10795
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The results told me that island area and the interaction terms have significant effect on SR:
# > anova(model)
#Type III Analysis of Variance Table with Satterthwaite's method
# Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#Area 3.0359 3.03590 1 27.448 4.4615 0.04390 *
#land_use 1.5520 0.77601 2 57.617 1.1404 0.32679
#Area:land_use 5.1658 2.58288 2 60.935 3.7958 0.02795 *
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
And then I used lsmeans function to conduct Tukeys' pairwise comparison:
#lsmeans(model, pairwise ~ Area:land_use, adjust="tukey")
The results indicate that the species richness from farmland and forest is significantly different, right? I wonder if this difference should be seen as the significant difference of intercept of the species richness-area relationship between farmland and forest in this model? That is the species richness from farmland transects is higher than that from forest transects?
#$contrasts
contrast estimate SE df t.ratio p.value
1.19968425045037 farmland - 1.19968425045037 forest 3.4153 0.288 62.6 1.185 0.0466
1.19968425045037 farmland - 1.19968425045037 mix -0.0306 0.426 64.0 -0.072 0.9972
1.19968425045037 forest - 1.19968425045037 mix -0.3722 0.377 63.9 -0.987 0.5087
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates
But how to test if the slope of the species richness-area relationship between farmland and forest in this model is significant different? That is to testify if the species richness-area relationship from farmland transects is more steeper than that from forest transect?
I think you want
lstrends(model, pairwise ~ land_use, var = "Area", adjust="tukey")
The functions lsmeans and lstrends are in the emmeans package, in which they are equivalent to emmeans and emtrends respectively. So look at the documentation for those functions. The lsmeans package is just a front end.