---
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?


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)
fit2 <- glm(native.richness~NLCD500m, data = beedata, family=poisson)
summary(fit1)
summary(fit2)

```

**Let's look at the output from the glm function more closely**
#Look at the p-values for fit 1 and fit 2. What do they tell us?


#Now let's plot our data
```{r}
#The following code 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))

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

e_richness3 <- within(e_richness2, {exotic.richness <- fit
LL <- fit - 1.96 * se.fit
UL <- fit + 1.96 * se.fit
})
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
})
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:
#Compare the above plot to the abundance plot in the focal paper. Is there a difference between abundance and species richness? 

#What is the species richness plot showing us?

Conservation Management Questions:
#Looping back to the big picture question,do you think do conservation management should implement stategies to support native species or do should there be a shift in supporting pollinator populations as a whole, regardless of whether or not it is detrimental to native species? Think about the data we looked at and the data in the focal paper