3_habitat_characteristics_nb.Rmd

---
title: "Urban Habitat Characteristic Analysis"
author: "Brandi Pessman"
date: "`r Sys.Date()`"
output: html_document
---

## Setting the Working Directory and Global Code Chunk Options

I have set message and warning to FALSE for every chunk so the knitted document does not include messages and warnings. 

```{r setup}
knitr::opts_chunk$set(message = FALSE, warning = FALSE, echo = FALSE)
require("knitr")
opts_knit$set(root.dir = "/Users/bjpessman/Documents/phd_research_code/Agelenopsis_aggregation")
```

# Packages to Load

```{r packages}
library(tidyverse) # for making graphs
library(MASS) # negative binomial tests
library(emmeans) # for calculating predictions for graphing
source('functions/residual_plots_function.R') # produces residual plots
```

# Import Data

```{r import}
sites_wrangled <- readRDS('wrangled_data/sites_wrangled.RDS')
```

# Comparing Habitat Characteristics by Location

Here, we tested whether each predictor differs between the two Locations: urban forest and urban center. We build models, check assumptions of the model, and make graphs to represent differences. 

```{r tree cover}
tree100m <- glm(tree100m ~ Location, family = "poisson", data = sites_wrangled)
summary(tree100m)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
tree100m.nb <- glm.nb(tree100m ~ Location, data = sites_wrangled)
summary(tree100m.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(tree100m),
            se = plotrix::std.error(tree100m))

# get predictions
predictions <- summary(emmeans(tree100m.nb, ~Location),type = "response")

# graph raw data and predictions
ggplot(sites_wrangled, aes(x = Location, y = tree100m)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = response, color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = response - SE, 
                  ymax = response + SE,
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(-1, 100), breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  theme_classic() +
  ylab("Percent Tree Cover \n[100m radius]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=100, label="***", color="black", size = 10)

# Residual plots
residual_plots(tree100m.nb, sites_wrangled)
```

Tree cover is **significantly higher in the urban forest** than the urban center (z = 13.77, df = 1, 21, p < 0.001). The means and SE's are 47.15 ± 3.46 % for the urban forest and 1.58 ± 0.37 % for the urban center. 

```{r impervious cover}
imperv100m <- glm(imperv100m ~ Location, family = "poisson", data = sites_wrangled)
summary(imperv100m)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
imperv100m.nb <- glm.nb(imperv100m ~ Location, data = sites_wrangled)
summary(imperv100m.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(imperv100m),
            se = plotrix::std.error(imperv100m))

# get predictions
predictions <- summary(emmeans(imperv100m.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = imperv100m)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = response, color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = response - SE, 
                  ymax = response + SE,
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(-1, 100), breaks = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) +
  theme_classic() +
  ylab("Percent Impervious Cover \n[100m radius]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=100, label="***", color="black", size = 10)

# Residual plots
residual_plots(imperv100m.nb, sites_wrangled)
```

Impervious cover is **significantly higher in the urban center** than the urban forest (z = -14.49, df = 1, 21, p < 0.001). The means and SE's are 1.27 ± 0.36 % for the urban forest and 76.52 ± 2.63 % for the urban center. 

```{r road disturbance}
tdr <- glm(Traffic_Dist_Road ~ Location, family = "poisson", data = sites_wrangled)
summary(tdr)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
tdr.nb <- glm.nb(Traffic_Dist_Road ~ Location, data = sites_wrangled)
summary(tdr.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(Traffic_Dist_Road),
            se = plotrix::std.error(Traffic_Dist_Road))

# get predictions
predictions <- summary(emmeans(tdr.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = Traffic_Dist_Road)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = (response), color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = (response - SE), 
                  ymax = (response + SE),
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 1400), breaks = c(0, 200, 400, 600, 800, 1000, 1200, 1400)) +
  theme_classic() +
  ylab("Road Disturbance \n[vehicles/day/m]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=1400, label="***", color="black", size = 10)

# Residual plots
residual_plots(tdr.nb, sites_wrangled)
```

**Road disturbance is higher in the urban center** than the urban forest (z = -2.87, df = 1, 21, p = 0.004). The means and SE's are 104 ± 36 vehicles/day/m for the urban forest and 398 ± 125 vehicles/day/m for the urban center. 

```{r highway disturbance}
tdh <- glm(Traffic_Dist_Highway ~ Location, family = "poisson", data = sites_wrangled)
summary(tdh)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
tdh.nb <- glm.nb(Traffic_Dist_Highway ~ Location, data = sites_wrangled)
summary(tdh.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(Traffic_Dist_Highway),
            se = plotrix::std.error(Traffic_Dist_Highway))

# get predictions
predictions <- summary(emmeans(tdh.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = Traffic_Dist_Highway)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = (response), color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = (response - SE), 
                  ymax = (response + SE),
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 250), breaks = c(0, 50, 100, 150, 200, 250)) +
  theme_classic() +
  ylab("Highway Disturbance \n[vehicles/day/m]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=250, label="***", color="black", size = 10)

# Residual plots
residual_plots(tdh.nb, sites_wrangled)
```

**Highway disturbance is higher in the urban center** than the urban fores (z = -2.64, df = 1, 21, p = 0.008). The means and SE's are 35.6 ± 6.04 vehicles/day/m for the urban forest and 65.0 ± 9.82 vehicles/day/m for the urban center. 

```{r plant species richness}
plant <- glm(TotalSub ~ Location, family = "poisson", data = sites_wrangled)
summary(plant)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
plant.nb <- glm.nb(TotalSub ~ Location, data = sites_wrangled)
summary(plant.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(TotalSub),
            se = plotrix::std.error(TotalSub))
# get predictions
predictions <- summary(emmeans(plant.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = (TotalSub))) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = (response), color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = (response - SE), 
                  ymax = (response + SE),
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 25), breaks = c(0, 5, 10, 15, 20, 25)) +
  theme_classic() +
  ylab("Plant Species Richness \n[species/site]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=25, label="***", color="black", size = 10)

# Residual plots
residual_plots(plant.nb, sites_wrangled)
```

Plant species richness is **significantly higher in the urban forest** than the urban center (z = 6.61, df = 1, 21, p < 0.001). The means and SE's are 13.30 ± 1.56 species for the urban forest and 3.33 ± 0.58 species for the urban center. 

```{r spectral radiance}
spec <- glm(spec_rad ~ Location, family = "poisson", data = sites_wrangled)
summary(spec)

spec.nb <- glm.nb(spec_rad ~ Location, data = sites_wrangled)
summary(spec.nb)

# calcuate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(spec_rad),
            se = plotrix::std.error(spec_rad))

# get predictions
predictions <- summary(emmeans(spec.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = spec_rad)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = response, color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = response - SE, 
                  ymax = response + SE,
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(130, 165), breaks = c(130, 135, 140, 145, 150, 155, 160, 165)) +
  theme_classic() +
  ylab("Spectral Radiance \n[Watts/(m² * sr * µm)]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=165, label="**", color="black", size = 10)

# Residual plots
residual_plots(spec.nb, sites_wrangled)
```

Spectral radiance is **significantly higher in the urban center** than the urban forest (z = -3.079, df = 1, 21, p = 0.002). The means and SE's are 137 ± 3.71 Watts/(m² * sr * µm) for the urban forest and 153 ± 3.58 Watts/(m² * sr * µm) for the urban center. 

```{r light at night}
light <- glm(light_rad ~ Location, family = "poisson", data = sites_wrangled)
summary(light)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
light.nb <- glm.nb(light_rad ~ Location, data = sites_wrangled)
summary(light.nb)

# calculate mean and std. error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(light_rad),
            se = plotrix::std.error(light_rad))

# get predictions
predictions <- summary(emmeans(light.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = light_rad)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = response, color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = response - SE, 
                  ymax = response + SE,
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 175), breaks = c(0, 25, 50, 75, 100, 125, 150, 175)) +
  theme_classic() +
  ylab("Radiance [mcd/m²]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=175, label="***", color="black", size = 10)

# Residual plots
residual_plots(light.nb, sites_wrangled)
```

Light radiance is **significantly higher in the urban center** than the urban forest (z = -17.17, df = 1, 21, p < 0.001). The means and SE's are 5.89 ± 0.91 mcd/m² for the urban forest and 116.90 ± 9.39 mcd/m² for the urban center.

```{r patch area}
patch <- glm(patch_area_mm ~ Location, family = "poisson", data = sites_wrangled)
summary(patch)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
patch.nb <- glm.nb(patch_area_mm ~ Location, data = sites_wrangled)
summary(patch.nb)

# calculate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(patch_area_mm),
            se = plotrix::std.error(patch_area_mm))

# get predictions
predictions <- summary(emmeans(patch.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = (patch_area_mm))) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = (response), color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = (response - SE), 
                  ymax = (response + SE),
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 130310), breaks = c(0, 20000, 40000, 60000, 80000, 100000, 120000)) +
  theme_classic() +
  ylab("Patch Area [mm²]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=130310, label="***", color="black", size = 10)

# Residual plots
residual_plots(patch.nb, sites_wrangled)
```

Patch area is **significantly higher in the urban forest** than the urban center (z = 13.08, df = 1, 21, p < 0.001). The means and SE's are 46,295 ± 11,128 mm² for the urban forest and 654 ± 144 mm² for the urban center. 

```{r total road length}
road <- glm(road_length_m ~ Location, family = "poisson", data = sites_wrangled)
summary(road)

# data are likely overdispersed since the residual deviance / degrees of freedom are much larger than 1
road.nb <- glm.nb(road_length_m ~ Location, data = sites_wrangled)
summary(road.nb)

# calculate mean and standard error
sites_wrangled %>% 
  group_by(Location) %>% 
  summarize(mean = mean(road_length_m),
            se = plotrix::std.error(road_length_m))

# get predictions
predictions <- summary(emmeans(road.nb, ~Location),type = "response")

ggplot(sites_wrangled, aes(x = Location, y = road_length_m)) + 
  geom_jitter(color = "grey", width = 0.1, size = 1) +
  geom_point(aes(x = Location, y = response, color = Location), size = 2, data = predictions) + 
  geom_errorbar(aes(x = Location, 
                  ymin = response - SE, 
                  ymax = response + SE,
                  color = Location), data = predictions, inherit.aes = FALSE, width = 0.15, linewidth = 1.5) +
  scale_color_manual(values = c("Urban Center" = "#d95f02", "Urban Forest" = "#1b9e77")) +
  scale_x_discrete(labels=c("Urban Forest" = "Urban \nForest \nN = 10", "Urban Center" = "Urban \nCenter \nN = 12")) +
  scale_y_continuous(limits=c(0, 140), breaks = c(0, 20, 40, 60, 80, 100, 120, 140)) +
  theme_classic() +
  ylab("Total Length of Road \n[m, 100m radius]") + 
  theme(text = element_text(size = 18)) + 
  theme(axis.text.x=element_text(colour="black", size=18)) + 
  theme(axis.text.y=element_text(colour="black", size=18),
        axis.title.x = element_blank()) +
  theme(legend.position = "none") +
  annotate(geom="text", x=1.5, y=140, label="***", color="black", size = 10)

# Residual plots
residual_plots(road.nb, sites_wrangled)
```

Total road length is **significantly higher in the urban center** than the urban forest (z = -2.74, df = 1, 21, p = 0.006). The means and SE's are 15.6 ± 7.54 m for the urban forest and 93.0 ± 40.60 m for the urban center.