---
title: "Does Urbanzation Favor Exotic Species?"
author: "Gabby Giacomangeli"
date: "2023-05-08"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown
Loading Packages
**If library() is not working, you need to install the package. Use install.package("nameofpackage") to install**

```{r}
#For this lesson, we will need to import the following packages:

#For Data Wrangling and Transformation:
library(dplyr)

#For Data Visualization and Exploration:
library(ggplot2)
library(AER)
library(MASS)
library(visreg)
library(AICcmodavg)

#For Data Import and Management
library(readr)
```

Importing Data Frame
```{r}
#Let's import our the data set:
beedata <- read.csv("Bee_Data.csv")

#Take a look at the data:
head(beedata)

```

```{r}
#Now we need the site codes, but we have to merge them to our bee data set. Look at the data in the environment and notice how "beedata" has 41 obs. of 32 variables

#Import site code data set:
sitecodes <- read.csv("SiteCodes.csv")

#Merge temperature data set to exotic bee data set:
beedata <- merge(beedata, sitecodes, by = "Locality", all.x=T)


#Take another look at the data in your environment. Notice how the "beedata" went from 41 obs of 32 variables to 41 obs. of 33 variables. This means we successfully added the new site codes column.

#Let's take a look at the data set to confirm that the site codes column is there:
head(beedata)

```

**At what buffer size did the authors use for their analysis? How did the authors of the focal paper determine that this km buffer is the best predictor of exotic bee richness? Let's find out**
```{r}

#We will use AIC. The AIC is designed to find the model that explains the most variation in the data. Once you've fit several regression models, you can compare the AIC value of each model. The lower the AIC, the better the model fit.

#Creating our models at each buffer:
r500m <- glm(exotic.richness~NLCD500m + Year, data = beedata, family = poisson)
r1km <- glm(exotic.richness~NLCD1km + Year, data = beedata, family = poisson)
r1.5km <- glm(exotic.richness~NLCD1.5km + Year, data = beedata, family = poisson)
r2km <- glm(exotic.richness~NLCD2km + Year, data = beedata, family = poisson)

#To compare these models and find which one is the best fit for the data, you can put them together into a list and use the aictab() command to compare all of them at once. Earlier, we loaded the library(AICcmodavg) which allows use to use the function aictab()

#Comparing our models:
cands <- list(r500m, r1km, r1.5km, r2km)
cand.names <- c("500m", "1km", "1.5km", "2km")
aictab(cand.set = cands, modnames = cand.names)

#Look at the ouput and determine which model is the better fit 
```

#Which buffer is the best predictor for exotic richness? How did you determine this?
**The 500m buffer is the best predictor for exotic richness. This was determined by comparing the AIC for each model. 500m had the lowest AIC, so this is the selected model. However, the Akaike’s rule of thumb is two models are essentially indistinguishable if the difference of their AICs is less than two. In the table above, the difference between 500m and 1km is less than 2, so while the 500m AIC is the lowest, it’s not distinguishable from the others.**


At this point, we are ready to perform our Poisson model analysis using the glm function. We fit the model and store it in the objects fit 1 & fit 2 and get a summary of both models

```{r}
#Does exotic richness increase with urbanization? How does this effect native richness? We will plot some data, but first we need to create models for our poisson regression

#Generalized linear models are fit using the glm( ) function. The glm() function makes it easy to perform other regression types since the only difference is the “family” argument. 


fit1<-glm(exotic.richness~NLCD500m,data = beedata,family = poisson(link="log"))

fit2 <- glm(native.richness~NLCD500m, data = beedata, family=poisson(link="log"))

summary(fit1)
summary(fit2)


```


#Now let's plot our data
```{r}
#The following code in this chunk was provided by the original R code (with minor changes) from Fitch, et al. 2019 that was not included in their final paper:

e_richness <- data.frame(NLCD500m=seq(0,100)/100,exotic.richness=rep(1,101))

e_richness2 <- cbind(e_richness, predict.glm(fit1, e_richness, type='response', se.fit=TRUE))

e_richness3 <- within(e_richness2, {exotic.richness <- fit
LL <- fit - 1.96 * se.fit
UL <- fit + 1.96 * se.fit
})

#What does predict.glm() do? It estimates standard errors of those predictions from a fitted generalized linear model object

n_richness <- data.frame(NLCD500m=seq(0,100)/100,native.richness=rep(1,101))
n_richness2 <- cbind(n_richness, predict.glm(fit2, n_richness, type='response', se.fit=TRUE))


n_richness2 <- within(n_richness2, {native.richness <- fit
LL <- fit - 1.96 * se.fit
UL <- fit + 1.96 * se.fit
})

#Now we use ggplot() to plot exotic and native species richness in response to proportional impervious surface:

ggplot(e_richness3) +
  theme_bw() +
  geom_point(data=beedata ,aes(x=NLCD500m, y=exotic.richness, col='Exotic'), size=2) +
  geom_ribbon(aes(x = NLCD500m, ymin = LL, ymax = UL), alpha = .25) +
  geom_line(aes(NLCD500m, exotic.richness, col='Exotic'), size = 2) +
  geom_point(data=beedata, aes(x=NLCD500m, y=native.richness, col='Native'), size=2) +
  geom_ribbon(data=n_richness2, aes(x = NLCD500m, ymin = LL, ymax = UL), alpha = .25) +
  geom_line(data=n_richness2, aes(NLCD500m, native.richness, col='Native'),size = 2) +
  labs(x = "Proportional impervious surface cover", y = "Species richness") + theme(axis.title = element_text(size=rel(1.5))) + theme(axis.title = element_text(size=rel(1.5)))
```

Questions:
Put back in R student version:
#Compare the above plot to the abundance plot in the focal paper. Is there a difference between abundance and species richness? 
**Similar to the abundance plot, richness of exotic bees increases with urbanization, but richness of native bees is largely unchanged.**

#What is the species richness plot showing us?
**This figure shows the relationship between urbanization and exotic bee proportional abundance. According to the figure, the absolute and proportional abundance and richness of exotic bees significantly increased with urbanization (Fitch, et al. 2019). Native bees showed no response to urbanization**