Section 13 Generalized linear mixed modeling (species richness, vegetation data and planting year)

In this script, we run generalized linear mixed models to test the association between first order jacknife scores and restoration type. In addition, we run generalized linear mixed models to test associations between species richness and habitat (vegetation structure) using site-pair name (actively restored and naturally regenerating were specified to be paired) and repeat visits as random effects. Lastly, we assess associations between year since restoration began and first order jackknife scores of birds.

13.1 Install required libraries

library(tidyverse)
library(dplyr)
library(stringr)
library(vegan)
library(ggplot2)
library(scico)
library(psych)
library(rcompanion)
library(multcomp)
library(lme4)
library(sjPlot)

# Source any custom/other internal functions necessary for analysis
source("code/01_internal-functions.R")

13.2 Load the necessary data for statistical modeling

# We load the subset data
datSubset <- read.csv("results/datSubset.csv")

# Load species-trait data to essentially check for associations by habitat type
trait_dat <- read.csv("data/species-trait-dat.csv")

# richness by Visit
# this data basically gives you richness per visit and adds a visit number for each consecutive visit to that site
richnessPerVisit <- datSubset %>%
  group_by(Site, Date, Restoration.type) %>%
  transform() %>% replace(is.na(.), 0) %>%
  summarise_at(.vars = vars(c("IP":"HSWP")),.funs = sum) %>% 
  mutate_at(vars(c("IP":"HSWP")),~ replace(., . > 0, 1)) %>%
  rowwise() %>% 
  mutate(richness = sum(c_across(IP:HSWP))) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  dplyr::select(Site, Restoration.type, Date, richness, visit, siteCode)

# combine the Detections dataframe with the trait dataset
nDetectionsTrait <- datSubset %>%
  group_by(Site, Restoration.type, Date) %>%
  transform() %>% replace(is.na(.), 0) %>%
  summarise_at(.vars = vars(c("IP":"HSWP")),.funs = sum) %>%
  pivot_longer(cols=IP:HSWP, names_to="Species_Code", values_to="count") %>%
  left_join(.,trait_dat, by=c("Species_Code"="species_annotation_codes")) %>%
  mutate(forRichness = case_when(count>0 ~ 1,count==0 ~ 0)) %>%
  rename(., nDetections = count)

# Load data from previous scripts for use in a GLM
vegData <-  read.csv("results/summaryVeg.csv") %>%
  filter(!str_detect(Site_ID, 'OLCAP5B'))
vegPcaScores <- read.csv("results/pcaVeg.csv") %>%
  filter(!str_detect(Site_ID, 'OLCAP5B'))
jackAll <- read.csv("results/jackAll.csv")
jackRainforest <- read.csv("results/jackRainforest.csv")
jackOpencountry <- read.csv("results/jackOpencountry.csv")

13.3 Getting data ready in a format for generalized linear modeling

# All birds
modelDataAll <- vegPcaScores %>%  
  rename(Site = Site_ID) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  rename(Restoration.type = Site_type) %>%
  mutate(across(Restoration.type, factor))  %>%
  add_column(jacknife = jackAll$jack1,
             year = vegData$plantingYear) %>% 
  mutate ("roundjk" = round(jacknife))

# rainforest birds
modelData_rainForest <- vegPcaScores %>%  
  rename(Site = Site_ID) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  rename(Restoration.type = Site_type) %>%
  mutate(across(Restoration.type, factor))  %>%
  add_column(jacknife = jackRainforest$jack1,
             year = vegData$plantingYear) %>% 
  mutate ("roundjk" = round(jacknife))

# open country birds
modelData_openCountry <- vegPcaScores %>%  
  rename(Site = Site_ID) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  rename(Restoration.type = Site_type) %>%
  mutate(across(Restoration.type, factor))  %>%
  add_column(jacknife = jackOpencountry$jack1,
             year = vegData$plantingYear) %>% 
  mutate ("roundjk" = round(jacknife))

13.4 Getting data ready for generalized linear mixed modeling with richness data as well as visits.

# data for the GLMM (overall richness)
glmmAll <- richnessPerVisit[,-2] %>%
  full_join(modelDataAll[,-8], by = c("Site","siteCode"))

# rainforest birds richness for glmm
glmmRainforest <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habitat == "RF") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelData_rainForest, by = c("Site","siteCode","Restoration.type"))
  
# open-country birds
glmmOpencountry <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habitat == "OC") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelData_openCountry, by = c("Site","siteCode","Restoration.type"))

# Let's look at species by foraging habit
# canopy birds
glmmCanopy <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habit == "CAN") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelDataAll[,-8], by = c("Site","siteCode","Restoration.type"))
  
# ground-feeding birds
glmmGround <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habit == "GRD") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelDataAll[,-8], by = c("Site","siteCode","Restoration.type"))

# mid-storey birds
glmmMidStorey <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habit == "MID") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelDataAll[,-8], by = c("Site","siteCode","Restoration.type"))

# understorey birds
glmmUnderStory <- nDetectionsTrait %>%
  group_by(Site, Restoration.type, Date) %>% 
  filter (habit == "UND") %>%
  summarise(richness = sum(forRichness)) %>%
  group_by(Site) %>% 
  mutate(visit = row_number()) %>%
  mutate(siteCode = str_extract(Site, pattern = "\\w+\\d+")) %>%
  mutate(siteCode = factor(siteCode)) %>%
  full_join(modelDataAll[,-8], by = c("Site","siteCode","Restoration.type"))

13.5 Running the generalized linear mixed models

We now run generalized linear mixed models (GLMM) assuming Poisson errors and using log link functions to examine the effects of restoration type (benchmark, actively restored and passively restored) on the jackknife estimates of bird species richness (for all, rainforest, and open-country species), followed by TukeyHSD multiple comparisons tests of means.

# all birds
glmm_alljk <- glmer(roundjk ~ Restoration.type +(1|siteCode), data = modelDataAll, family = poisson(link = log))
summary(glmm_alljk)

tukey_glmmAllJack <- summary(glht(glmm_alljk, linfct=mcp(Restoration.type ="Tukey")))
cld(tukey_glmmAllJack)

# The above result suggests that there is a significant difference in first order jacknife estimates for benchmark sites and passively restored site (but no difference between active-passive and active-benchmark).

# rainforest birds
glmm_rainForestJack <- glmer(roundjk ~ Restoration.type + +(1|siteCode), data = modelData_rainForest, family = poisson(link = log))
summary(glmm_rainForestJack)

tukey_glmmRainForestJack <- summary(glht(glmm_rainForestJack, linfct=mcp(Restoration.type ="Tukey")))
cld(tukey_glmmRainForestJack)

# The above result suggests no significant difference in means between any treatment types

# open country birds
glmm_openCountryJack <- glmer(roundjk ~ Restoration.type + +(1|siteCode), data = modelData_openCountry, family = poisson(link = log))
summary(glmm_openCountryJack)

tukey_glmmOpenCountryJack <- summary(glht(glmm_openCountryJack, linfct=mcp(Restoration.type ="Tukey")))
cld(tukey_glmmOpenCountryJack)

# For open country birds, there is a significant difference in first-order jacknife estimates between benchmark and passive sites and benchmark and active sites

13.6 Testing the role of habitat

Testing the role of habitat (vegetation structure) and foraging habit on species richness within a generalized linear modeling framework

# Testing the role of habitat first

# all bird species
glmm_allBirds <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmAll, family = poisson(link = log))
summary(glmm_allBirds)

plot_model(glmm_allBirds, type="pred", terms=c("PC1","PC2"))
report::report(glmm_allBirds)

# significant negative association with PC2

# rainforest birds
glmm_Rainforest <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmRainforest, family = poisson(link = log))
summary(glmm_Rainforest)

plot_model(glmm_Rainforest, type="pred", terms=c("PC1","PC2"))
report::report(glmm_Rainforest)

# Results above suggest PC2 is significantly negatively associated with richness of rainforest birds. 

glmm_openCountry <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmOpencountry, family = poisson(link = log))
summary(glmm_openCountry)

plot_model(glmm_openCountry, type="pred", terms=c("PC1","PC2"))
report::report(glmm_openCountry)

# statistically significant positive association with PC1 and significant negative association with PC2

# Testing the role of foraging habit

# canopy birds (no significant association)
glmm_Canopy <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmCanopy, family = poisson(link = log))
summary(glmm_Canopy)

# ground-feeding birds (Marginal association between richness of ground-feeding birds and PC2)
glmm_Ground <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmGround, family = poisson(link = log))
summary(glmm_Ground)

# mid-storey birds (Marginal association between PC1 and richness of mid-storey birds)
glmm_MidStorey <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmMidStorey, family = poisson(link = log))
summary(glmm_MidStorey)

# understory birds (no significant association)
glmm_Understory <- glmer(richness ~ PC1 + PC2 + (1|siteCode) + (1|visit), data = glmmUnderStory, family = poisson(link = log))
summary(glmm_Understory)

13.7 Effect of planting year

Lastly, running a generalized linear model to test the effect of year since restoration on jacknife estimates of species richness

# Let's look at overall species richness first
# filter only data for restored sites
allRestored <- glmmAll %>%
  filter(Restoration.type=="Active") %>%
  mutate(yearSinceRestoration = (2022-year))

glmmAllRest <- glmer(roundjk ~ yearSinceRestoration + (1|visit), data = allRestored, family = poisson(link = log))
summary(glmmAllRest)

plot_model(glmmAllRest, type="pred")
report::report(glmmAllRest)

# no significant association for overall richness

# rainforest bird richness
rainRestored <- glmmRainforest %>%
  filter(Restoration.type=="Active") %>%
  mutate(yearSinceRestoration = (2022-year))

glmmRain <- glmer(roundjk ~ yearSinceRestoration + (1|visit), data = rainRestored, family = poisson(link = log))
summary(glmmRain)

plot_model(glmmRain, type="pred")
report::report(glmmRain)

# no significant association for rainforest bird species richness

# open-country richness
openRestored <- glmmOpencountry %>%
  filter(Restoration.type=="Active") %>%
  mutate(yearSinceRestoration = (2022-year))

glmmOpen <- glmer(roundjk ~ yearSinceRestoration + (1|visit), data = openRestored, family = poisson(link = log))
summary(glmmOpen)

plot_model(glmmOpen, type="pred")
report::report(glmmOpen)

# The effect of yearSinceRestoration is statistically significant and positive for open country birds (beta = 0.02, 95% CI [4.05e-03, 0.04], p = 0.016; Std. beta = 0.06, 95% CI [0.01, 0.11])