### Analysis file for Mule Deer 420 Project ###

#Tell R where to find the data
setwd('C:/Users/ranglackdh/Documents/Students/Undergrad/Bryan/Manuscript/Revision2')

#Read in the data file, call it df
df <- read.csv("analysis.csv", header = TRUE, sep = ",")

#Check data structure, set year as factor
str(df)
df$year <- as.factor(df$year)

names(df)

summary(df)

#Pseudothreshold functional form 
df2 <- log((df[, c(8:20)])+0.001)
colnames(df2) <- paste(colnames(df2), "ps", sep = "_")

df <- cbind(df, df2)

#standardize covariates 
df$elev_std <- ((df$elev -(mean(df$elev)))/(2*sd(df$elev)))
df$slope_std <- ((df$slope -(mean(df$slope)))/(2*sd(df$slope)))
df$rough_std <- ((df$rough -(mean(df$rough)))/(2*sd(df$rough)))
df$ag_std <- ((df$ag -(mean(df$ag)))/(2*sd(df$ag)))
df$canopy_std <- ((df$canopy -(mean(df$canopy)))/(2*sd(df$canopy)))
df$forest_std <- ((df$forest -(mean(df$forest)))/(2*sd(df$forest)))
df$ndviamp_std <- ((df$ndviamp -(mean(df$ndviamp)))/(2*sd(df$ndviamp)))
df$ndvitin_std <- ((df$ndvitin -(mean(df$ndvitin)))/(2*sd(df$ndvitin)))
df$develop_std <- ((df$develop -(mean(df$develop)))/(2*sd(df$develop)))
df$range_std <- ((df$range -(mean(df$range)))/(2*sd(df$range)))
df$riparian_std <- ((df$riparian -(mean(df$riparian)))/(2*sd(df$riparian)))
df$road_std <- ((df$road -(mean(df$road)))/(2*sd(df$road)))
df$urban_std <- ((df$urban -(mean(df$urban)))/(2*sd(df$urban)))
df$elev_ps_std <- ((df$elev_ps -(mean(df$elev_ps)))/(2*sd(df$elev_ps)))
df$slope_ps_std <- ((df$slope_ps -(mean(df$slope_ps)))/(2*sd(df$slope_ps)))
df$rough_ps_std <- ((df$rough_ps -(mean(df$rough_ps)))/(2*sd(df$rough_ps)))
df$ag_ps_std <- ((df$ag_ps -(mean(df$ag_ps)))/(2*sd(df$ag_ps)))
df$canopy_ps_std <- ((df$canopy_ps -(mean(df$canopy_ps)))/(2*sd(df$canopy_ps)))
df$forest_ps_std <- ((df$forest_ps -(mean(df$forest_ps)))/(2*sd(df$forest_ps)))
df$ndviamp_ps_std <- ((df$ndviamp_ps -(mean(df$ndviamp_ps)))/(2*sd(df$ndviamp_ps)))
df$ndvitin_ps_std <- ((df$ndvitin_ps -(mean(df$ndvitin_ps)))/(2*sd(df$ndvitin_ps)))
df$develop_ps_std <- ((df$develop_ps -(mean(df$develop_ps)))/(2*sd(df$develop_ps)))
df$range_ps_std <- ((df$range_ps -(mean(df$range_ps)))/(2*sd(df$range_ps)))
df$riparian_ps_std <- ((df$riparian_ps -(mean(df$riparian_ps)))/(2*sd(df$riparian_ps)))
df$road_ps_std <- ((df$road_ps -(mean(df$road_ps)))/(2*sd(df$road_ps)))
df$urban_ps_std <- ((df$urban_ps -(mean(df$urban_ps)))/(2*sd(df$urban_ps)))

#Quadratic Functional form
df3 <- (df[, c(34:37,40,41,43,45)]^2)
colnames(df3) <- paste(colnames(df3), "sq", sep = "_")

df <- cbind(df, df3)

#Export standarized data to file so you can start here next time 
write.csv(df, "analysis_2.csv")
df <- read.csv("analysis_2.csv", header = TRUE, sep = ",")

df$effort_ps <- log(df$effort+0.001)
df$effort_std <- ((df$effort -(mean(df$effort)))/(2*sd(df$effort)))
df$effort_ps_std <- ((df$effort_ps -(mean(df$effort_ps)))/(2*sd(df$effort_ps)))
df$effort_std_sq <- df$effort_std^2

df$year <- as.factor(df$year)

#Add a very small number to the harvest density so you can fit a Gamma link GLM
df$deer_100km2 <- df$deer_100km +0.0000001

#Look at correlations between covariate
#install.packages("GGally")
#install.packages("ggplot2")
library(GGally)
library(ggplot2)

#create Correlation matrix
ggcorr(df[c(35:47,69)],
       label=TRUE, label_round = 2,
       max_size = 6,
       size = 4,
       hjust = 0.75,
       angle = -45,
       palette = "PuOr") + labs(title = "Covariate Correlation")

#Fit models to determine functional form

#create text for elevation model statements

elev_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

elev_models$model=paste("m",1:3,"_elev",sep="")

elev_models$cov.string[1]="elev_std"
elev_models$cov.string[2]="elev_ps_std"
elev_models$cov.string[3]="elev_std + elev_std_sq"

#create text for effort model statements

effort_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

effort_models$model=paste("m",1:3,"_effort",sep="")

effort_models$cov.string[1]="effort_std"
effort_models$cov.string[2]="effort_ps_std"
effort_models$cov.string[3]="effort_std + effort_std_sq"

#create table of slope models

slope_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

slope_models$model=paste("m",1:3,"_slope",sep="")

slope_models$cov.string[1]="slope_std"
slope_models$cov.string[2]="slope_ps_std"
slope_models$cov.string[3]="slope_std + slope_std_sq"

#create table of roughness model statements

rough_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

rough_models$model=paste("m",1:3,"_rough",sep="")

rough_models$cov.string[1]="rough_std"
rough_models$cov.string[2]="rough_ps_std"
rough_models$cov.string[3]="rough_std + rough_std_sq"

#create table of ag models

ag_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

ag_models$model=paste("m",1:3,"_ag",sep="")

ag_models$cov.string[1]="ag_std"
ag_models$cov.string[2]="ag_ps_std"
ag_models$cov.string[3]="ag_std + ag_std_sq"

#create table of canopy models

canopy_models=data.frame(matrix(nrow=2,ncol=2,dimnames=list(1:2,c("model","cov.string"))))

canopy_models$model=paste("m",1:2,"_canopy",sep="")

canopy_models$cov.string[1]="canopy_std"
canopy_models$cov.string[2]="canopy_ps_std"

#create table of forest models

forest_models=data.frame(matrix(nrow=2,ncol=2,dimnames=list(1:2,c("model","cov.string"))))

forest_models$model=paste("m",1:2,"_forest",sep="")

forest_models$cov.string[1]="forest_std"
forest_models$cov.string[2]="forest_ps_std"

#create table of NDVI amplitude models

ndviamp_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

ndviamp_models$model=paste("m",1:3,"_ndviamp",sep="")

ndviamp_models$cov.string[1]="ndviamp_std"
ndviamp_models$cov.string[2]="ndviamp_ps_std"
ndviamp_models$cov.string[3]="ndviamp_std + ndviamp_std_sq"

#create table of Time Intergrated NDVI models

ndvitin_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

ndvitin_models$model=paste("m",1:3,"_ndvitin",sep="")

ndvitin_models$cov.string[1]="ndvitin_std"
ndvitin_models$cov.string[2]="ndvitin_ps_std"
ndvitin_models$cov.string[3]="ndvitin_std + ndvitin_std_sq"

#create table of development models

develop_models=data.frame(matrix(nrow=2,ncol=2,dimnames=list(1:2,c("model","cov.string"))))

develop_models$model=paste("m",1:2,"_develop",sep="")

develop_models$cov.string[1]="develop_std"
develop_models$cov.string[2]="develop_ps_std"

#create table of range models

range_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

range_models$model=paste("m",1:3,"_range",sep="")

range_models$cov.string[1]="range_std"
range_models$cov.string[2]="range_ps_std"
range_models$cov.string[3]="range_std + range_std_sq"

#create table of riparian models

riparian_models=data.frame(matrix(nrow=2,ncol=2,dimnames=list(1:2,c("model","cov.string"))))

riparian_models$model=paste("m",1:2,"_riparian",sep="")

riparian_models$cov.string[1]="riparian_std"
riparian_models$cov.string[2]="riparian_ps_std"

#create table of road density models

road_models=data.frame(matrix(nrow=3,ncol=2,dimnames=list(1:3,c("model","cov.string"))))

road_models$model=paste("m",1:3,"_road",sep="")

road_models$cov.string[1]="road_std"
road_models$cov.string[2]="road_ps_std"
road_models$cov.string[3]="road_std + road_std_sq"

#create table of urbanization models

urban_models=data.frame(matrix(nrow=2,ncol=2,dimnames=list(1:2,c("model","cov.string"))))

urban_models$model=paste("m",1:2,"_urban",sep="")

urban_models$cov.string[1]="urban_std"
urban_models$cov.string[2]="urban_ps_std"

#create a list to hold the fitted elevation models 

mod.list.elev=vector("list",3)

#Run elev suite of models 

mod.list.elev[[1]]<-glm(paste("deer_100km2 ~",elev_models$cov.string[1]), family = Gamma, data=df)
mod.list.elev[[2]]<-glm(paste("deer_100km2 ~",elev_models$cov.string[2]), family = Gamma, data=df)
mod.list.elev[[3]]<-glm(paste("deer_100km2 ~",elev_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

#install.packages("AICcmodavg")
library(AICcmodavg)
#install.packages("XLConnect")
library(XLConnect)
#install.packages("MuMIn")
library(MuMIn)

elev_aic <- aictab(cand.set = mod.list.elev, modnames = elev_models$cov.string, sort = TRUE)

elev_tab <- model.sel(mod.list.elev)

writeWorksheetToFile(file = "func_forms.xlsx", data = elev_aic, sheet = "elev_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = elev_tab, sheet = "elev_tab")

#create a list to hold the fitted effort models 

mod.list.effort=vector("list",3)

#Run effort suite of models 

mod.list.effort[[1]]<-glm(paste("deer_100km2 ~",effort_models$cov.string[1]), family = Gamma, data=df)
mod.list.effort[[2]]<-glm(deer_100km2 ~ effort_ps_std, family = Gamma, data=df, start = c(0.5, 0.5))
mod.list.effort[[3]]<-glm(paste("deer_100km2 ~",effort_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table
effort_aic <- aictab(cand.set = mod.list.effort, modnames = effort_models$cov.string, sort = TRUE)

effort_tab <- model.sel(mod.list.effort)

writeWorksheetToFile(file = "func_forms.xlsx", data = effort_aic, sheet = "effort_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = effort_tab, sheet = "effort_tab")


#create a list to hold the fitted slope models 

mod.list.slope=vector("list",3)

#Run slope suite of models 

mod.list.slope[[1]]<-glm(paste("deer_100km2 ~",slope_models$cov.string[1]), family = Gamma, data=df)
mod.list.slope[[2]]<-glm(paste("deer_100km2 ~",slope_models$cov.string[2]), family = Gamma, data=df)
mod.list.slope[[3]]<-glm(paste("deer_100km2 ~",slope_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

slope_aic <- aictab(cand.set = mod.list.slope, modnames = slope_models$cov.string, sort = TRUE)

slope_tab <- model.sel(mod.list.slope)

writeWorksheetToFile(file = "func_forms.xlsx", data = slope_aic, sheet = "slope_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = slope_tab, sheet = "slope_tab")

#create a list to hold the fitted rough models 

mod.list.rough=vector("list",3)

#Run rough suite of models 

mod.list.rough[[1]]<-glm(paste("deer_100km2 ~",rough_models$cov.string[1]), family = Gamma, data=df)
mod.list.rough[[2]]<-glm(paste("deer_100km2 ~",rough_models$cov.string[2]), family = Gamma, data=df)
mod.list.rough[[3]]<-glm(paste("deer_100km2 ~",rough_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

rough_aic <- aictab(cand.set = mod.list.rough, modnames = rough_models$cov.string, sort = TRUE)

rough_tab <- model.sel(mod.list.rough)

writeWorksheetToFile(file = "func_forms.xlsx", data = rough_aic, sheet = "rough_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = rough_tab, sheet = "rough_tab")

#create a list to hold the fitted ag models 

mod.list.ag=vector("list",3)

#Run ag suite of models 

mod.list.ag[[1]]<-glm(paste("deer_100km2 ~",ag_models$cov.string[1]), family = Gamma, data=df)
mod.list.ag[[2]]<-glm(paste("deer_100km2 ~",ag_models$cov.string[2]), family = Gamma, data=df)
mod.list.ag[[3]]<-glm(paste("deer_100km2 ~",ag_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

ag_aic <- aictab(cand.set = mod.list.ag, modnames = ag_models$cov.string, sort = TRUE)

ag_tab <- model.sel(mod.list.ag)

writeWorksheetToFile(file = "func_forms.xlsx", data = ag_aic, sheet = "ag_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = ag_tab, sheet = "ag_tab")

#create a list to hold the fitted canopy models 

mod.list.canopy=vector("list",2)

#Run canopy suite of models 

mod.list.canopy[[1]]<-glm(paste("deer_100km2 ~",canopy_models$cov.string[1]), family = Gamma, data=df)
mod.list.canopy[[2]]<-glm(paste("deer_100km2 ~",canopy_models$cov.string[2]), family = Gamma, data=df)

#generate AIC table

canopy_aic <- aictab(cand.set = mod.list.canopy, modnames = canopy_models$cov.string, sort = TRUE)

canopy_tab <- model.sel(mod.list.canopy)

writeWorksheetToFile(file = "func_forms.xlsx", data = canopy_aic, sheet = "canopy_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = canopy_tab, sheet = "canopy_tab")

#create a list to hold the fitted forest models 

mod.list.forest=vector("list",2)

#Run forest suite of models 

mod.list.forest[[1]]<-glm(paste("deer_100km2 ~",forest_models$cov.string[1]), family = Gamma, data=df)
mod.list.forest[[2]]<-glm(paste("deer_100km2 ~",forest_models$cov.string[2]), family = Gamma, data=df)

#generate AIC table

forest_aic <- aictab(cand.set = mod.list.forest, modnames = forest_models$cov.string, sort = TRUE)

forest_tab <- model.sel(mod.list.forest)

writeWorksheetToFile(file = "func_forms.xlsx", data = forest_aic, sheet = "forest_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = forest_tab, sheet = "forest_tab")

#create a list to hold the fitted ndvi amplitude models 

mod.list.ndviamp=vector("list",3)

#Run ndvi amplitude of models 

mod.list.ndviamp[[1]]<-glm(paste("deer_100km2 ~",ndviamp_models$cov.string[1]), family = Gamma, data=df)
mod.list.ndviamp[[2]]<-glm(paste("deer_100km2 ~",ndviamp_models$cov.string[2]), family = Gamma, data=df)
mod.list.ndviamp[[3]]<-glm(paste("deer_100km2 ~",ndviamp_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

ndviamp_aic <- aictab(cand.set = mod.list.ndviamp, modnames = ndviamp_models$cov.string, sort = TRUE)

ndviamp_tab <- model.sel(mod.list.ndviamp)

writeWorksheetToFile(file = "func_forms.xlsx", data = ndviamp_aic, sheet = "ndviamp_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = ndviamp_tab, sheet = "ndviamp_tab")

#create a list to hold the fitted Time intergrated ndvi models 

mod.list.ndvitin=vector("list",3)

#Run Time intergrated ndvi suite of models 

mod.list.ndvitin[[1]]<-glm(paste("deer_100km2 ~",ndvitin_models$cov.string[1]), family = Gamma, data=df)
mod.list.ndvitin[[2]]<-glm(paste("deer_100km2 ~",ndvitin_models$cov.string[2]), family = Gamma, data=df)
mod.list.ndvitin[[3]]<-glm(paste("deer_100km2 ~",ndvitin_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

ndvitin_aic <- aictab(cand.set = mod.list.ndvitin, modnames = ndvitin_models$cov.string, sort = TRUE)

ndvitin_tab <- model.sel(mod.list.ndvitin)

writeWorksheetToFile(file = "func_forms.xlsx", data = ndvitin_aic, sheet = "ndvitin_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = ndvitin_tab, sheet = "ndvitin_tab")

#create a list to hold the fitted development models 

mod.list.develop=vector("list",2)

#Run slope development of models 

mod.list.develop[[1]]<-glm(paste("deer_100km2 ~",develop_models$cov.string[1]), family = Gamma, data=df)
mod.list.develop[[2]]<-glm(paste("deer_100km2 ~",develop_models$cov.string[2]), family = Gamma, data=df)

#generate AIC table

develop_aic <- aictab(cand.set = mod.list.develop, modnames = develop_models$cov.string, sort = TRUE)

develop_tab <- model.sel(mod.list.develop)

writeWorksheetToFile(file = "func_forms.xlsx", data = develop_aic, sheet = "develop_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = develop_tab, sheet = "develop_tab")

#create a list to hold the fitted range models 

mod.list.range=vector("list",3)

#Run range suite of models 

mod.list.range[[1]]<-glm(paste("deer_100km2 ~",range_models$cov.string[1]), family = Gamma, data=df)
mod.list.range[[2]]<-glm(paste("deer_100km2 ~",range_models$cov.string[2]), family = Gamma, data=df)
mod.list.range[[3]]<-glm(paste("deer_100km2 ~",range_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

range_aic <- aictab(cand.set = mod.list.range, modnames = range_models$cov.string, sort = TRUE)

range_tab <- model.sel(mod.list.range)

writeWorksheetToFile(file = "func_forms.xlsx", data = range_aic, sheet = "range_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = range_tab, sheet = "range_tab")

#create a list to hold the fitted riparian models 

mod.list.riparian=vector("list",2)

#Run riparian suite of models 

mod.list.riparian[[1]]<-glm(paste("deer_100km2 ~",riparian_models$cov.string[1]), family = Gamma, data=df)
mod.list.riparian[[2]]<-glm(paste("deer_100km2 ~",riparian_models$cov.string[2]), family = Gamma, data=df)

#generate AIC table

riparian_aic <- aictab(cand.set = mod.list.riparian, modnames = riparian_models$cov.string, sort = TRUE)

riparian_tab <- model.sel(mod.list.riparian)

writeWorksheetToFile(file = "func_forms.xlsx", data = riparian_aic, sheet = "riparian_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = riparian_tab, sheet = "riparian_tab")

#create a list to hold the fitted road models 

mod.list.road=vector("list",3)

#Run road suite of models 

mod.list.road[[1]]<-glm(paste("deer_100km2 ~",road_models$cov.string[1]), family = Gamma, data=df)
mod.list.road[[2]]<-glm(paste("deer_100km2 ~",road_models$cov.string[2]), family = Gamma, data=df)
mod.list.road[[3]]<-glm(paste("deer_100km2 ~",road_models$cov.string[3]), family = Gamma, data=df)

#generate AIC table

road_aic <- aictab(cand.set = mod.list.road, modnames = road_models$cov.string, sort = TRUE)

road_tab <- model.sel(mod.list.road)

writeWorksheetToFile(file = "func_forms.xlsx", data = road_aic, sheet = "road_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = road_tab, sheet = "road_tab")

#create a list to hold the fitted urban models 

mod.list.urban=vector("list",2)

#Run urban suite of models 

mod.list.urban[[1]]<-glm(paste("deer_100km2 ~",urban_models$cov.string[1]), family = Gamma, data=df)
mod.list.urban[[2]]<-glm(paste("deer_100km2 ~",urban_models$cov.string[2]), family = Gamma, data=df)

#generate AIC table

urban_aic <- aictab(cand.set = mod.list.urban, modnames = urban_models$cov.string, sort = TRUE)

urban_tab <- model.sel(mod.list.urban)

writeWorksheetToFile(file = "func_forms.xlsx", data = urban_aic, sheet = "urban_aic")
writeWorksheetToFile(file = "func_forms.xlsx", data = urban_tab, sheet = "urban_tab")

#Combine top functional forms into top model

full <- glm(deer_100km2 ~ urban_std + elev_std + elev_std_sq + slope_std + slope_std_sq + slope_ps_std + rough_ps_std +
              ag_std + ag_std_sq + canopy_std + forest_std + ndviamp_std + ndviamp_std_sq +
              ndvitin_std + ndvitin_std_sq + develop_std + range_std + range_std_sq + riparian_std + riparian_ps_std + road_std +
              road_std_sq + year + effort_std + effort_std_sq + effort_ps_std, family = Gamma, data=df, na.action=na.fail)

mods <- dredge(full, subset = !(elev_std_sq && ndviamp_std) && !(elev_std_sq && ndvitin_std) &&
               !(elev_std_sq && ndvitin_std_sq) && !(elev_std_sq && range_std_sq) && !(slope_std && rough_ps_std) && 
               !(slope_std_sq && rough_ps_std) && !(rough_ps_std && ag_std_sq) && !(ag_std_sq && ndviamp_std) &&
               !(ag_std_sq && ndvitin_std) && !(ag_std_sq && ndvitin_std_sq) && !(ag_std_sq && range_std_sq) &&
               !(canopy_std && forest_std) && !(ndviamp_std && ndvitin_std) &&
               !(ndviamp_std && ndvitin_std_sq) && !(ndviamp_std && range_std_sq) && !(ndvitin_std && range_std_sq) &&
               !(ndvitin_std_sq && range_std_sq) && !(develop_std && road_std_sq) && !(develop_std && urban_std) &&
               !(range_std_sq && road_std_sq) && !(road_std_sq && urban_std) && !(elev_std_sq && canopy_std) &&
               !(slope_ps_std && slope_std) && !(slope_ps_std && rough_ps_std) && !(riparian_ps_std && riparian_std) &&
               !(effort_std_sq && effort_ps_std) &&
               c(dc(elev_std, elev_std_sq), dc(elev_std_sq, elev_std), dc(slope_std, slope_std_sq), dc(slope_std_sq, slope_std),
               dc(ag_std, ag_std_sq), dc(ag_std_sq, ag_std), dc(ndvitin_std, ndvitin_std_sq), dc(ndvitin_std_sq, ndvitin_std),
               dc(range_std, range_std_sq), dc(range_std_sq, range_std), dc(road_std, road_std_sq), dc(road_std_sq, road_std),
               dc(ndviamp_std, ndviamp_std_sq), dc(ndviamp_std_sq, ndviamp_std), dc(effort_std, effort_std_sq), 
               dc(effort_std_sq, effort_std)))

library(XLConnect)
options(java.parameters = "-Xmx")
write.csv(mods, file = "dredge.csv")

top <- glm(deer_100km2 ~ canopy_std + effort_ps_std + ndviamp_std + ndviamp_std_sq  + road_std + road_std_sq + 
             rough_ps_std, family = Gamma, data=df, na.action=na.fail)

df$ndviamp_sq <- df$ndviamp^2
df$road_sq <- df$road^2

top.orig <- glm(deer_100km2 ~ canopy + log(effort+0.001) + ndviamp + ndviamp_sq + road + road_sq + log(rough+0.001),
                family = Gamma, data=df, na.action=na.fail)

summary(top)
summary(top.orig)

##Plot the results
#Canopy Cover
dat.predict1=expand.grid(canopy = seq(min(df$canopy),max(df$canopy),length=100),
                         effort = mean (df$effort),
                         ndviamp = mean(df$ndviamp),
                         rough = mean(df$rough),
                         road = mean(df$road))

dat.predict1$ndviamp_sq <- dat.predict1$ndviamp^2
dat.predict1$road_sq <- dat.predict1$road^2

# predict X*Beta, i.e., linear predictor 
pred1=predict(top.orig, newdata = dat.predict1 ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
dat.predict1$fit=pred1$fit
dat.predict1$se.fit=pred1$se.fit
dat.predict1$lci=dat.predict1$fit-1.96*dat.predict1$se.fit
dat.predict1$uci=dat.predict1$fit+1.96*dat.predict1$se.fit

##Plotting
library(ggplot2)
require(cowplot)

fit.canopy <- ggplot(dat.predict1, aes(x=canopy,y=fit)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=lci, ymax=uci, linetype=NA), alpha=0.2, show.legend=FALSE) + 
  theme(legend.title=element_blank(), legend.position=c(.65,.1), legend.key.size = unit(.25, "cm")) +
  xlab("Canopy Cover (%)") + 
  ylab(bquote('Predicted MD Harvest/100km'^~2))

#NDVI
dat.predict2=expand.grid(canopy = mean(df$canopy),
                         effort = mean (df$effort),
                         ndviamp = seq(min(df$ndviamp),max(df$ndviamp),length=100),
                         rough = mean(df$rough),
                         road = mean(df$road))

dat.predict2$ndviamp_sq <- dat.predict2$ndviamp^2
dat.predict2$road_sq <- dat.predict2$road^2

# predict X*Beta, i.e., linear predictor 
pred2=predict(top.orig, newdata = dat.predict2 ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
dat.predict2$fit=pred2$fit
dat.predict2$se.fit=pred2$se.fit
dat.predict2$lci=dat.predict2$fit-1.96*dat.predict2$se.fit
dat.predict2$uci=dat.predict2$fit+1.96*dat.predict2$se.fit

##Plotting
fit.ndvi <- ggplot(dat.predict2, aes(x=ndviamp,y=fit)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=lci, ymax=uci, linetype=NA), alpha=0.2, show.legend=FALSE) + 
  theme(legend.title=element_blank(), legend.position=c(.65,.1), legend.key.size = unit(.25, "cm")) +
  xlab("NDVI Amplitude") + 
  ylab(bquote('Predicted MD Harvest/100km'^~2))

#effort
dat.predict3=expand.grid(canopy = mean(df$canopy),
                         ndviamp = mean(df$ndviamp),
                         rough = mean(df$rough),
                         effort = seq(min(df$effort),max(df$effort),length=100),
                         road = mean(df$road))

dat.predict3$ndviamp_sq <- dat.predict3$ndviamp^2
dat.predict3$road_sq <- dat.predict3$road^2

# predict X*Beta, i.e., linear predictor 
pred4=predict(top.orig, newdata = dat.predict3 ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
dat.predict3$fit=pred4$fit
dat.predict3$se.fit=pred4$se.fit
dat.predict3$lci=dat.predict3$fit-1.96*dat.predict3$se.fit
dat.predict3$uci=dat.predict3$fit+1.96*dat.predict3$se.fit

fit.effort <- ggplot(dat.predict3, aes(x=effort,y=fit)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=lci, ymax=uci, linetype=NA), alpha=0.2, show.legend=FALSE) + 
  theme(legend.title=element_blank(), legend.position=c(.65,.1), legend.key.size = unit(.25, "cm")) +
  xlab("Mean Hunter Effort (days)") + 
  ylab(bquote('Predicted MD Harvest/100km'^~2))


#rough
dat.predict4=expand.grid(canopy = mean(df$canopy),
                         effort = mean (df$effort),
                         ndviamp = mean(df$ndviamp),
                         rough = seq(min(df$rough),max(df$rough),length=100),
                         road = mean(df$road))

dat.predict4$ndviamp_sq <- dat.predict4$ndviamp^2
dat.predict4$road_sq <- dat.predict4$road^2

# predict X*Beta, i.e., linear predictor 
pred4=predict(top.orig, newdata = dat.predict4 ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
dat.predict4$fit=pred4$fit
dat.predict4$se.fit=pred4$se.fit
dat.predict4$lci=dat.predict4$fit-1.96*dat.predict4$se.fit
dat.predict4$uci=dat.predict4$fit+1.96*dat.predict4$se.fit

fit.rough <- ggplot(dat.predict4, aes(x=rough,y=fit)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=lci, ymax=uci, linetype=NA), alpha=0.2, show.legend=FALSE) + 
  theme(legend.title=element_blank(), legend.position=c(.65,.1), legend.key.size = unit(.25, "cm")) +
  xlab("Roughness (%)") + 
  ylab(bquote('Predicted MD Harvest/100km'^~2))


#Road
dat.predict5=expand.grid(canopy = mean(df$canopy),
                         effort = mean (df$effort),
                         ndviamp = mean(df$ndviamp),
                         rough = mean(df$rough),
                         road = seq(min(df$road),max(df$road),length=100))

dat.predict5$ndviamp_sq <- dat.predict5$ndviamp^2
dat.predict5$road_sq <- dat.predict5$road^2

# predict X*Beta, i.e., linear predictor 
pred5=predict(top.orig, newdata = dat.predict5 ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
dat.predict5$fit=pred5$fit
dat.predict5$se.fit=pred5$se.fit
dat.predict5$lci=dat.predict5$fit-1.96*dat.predict5$se.fit
dat.predict5$uci=dat.predict5$fit+1.96*dat.predict5$se.fit


fit.road <- ggplot(dat.predict5, aes(x=road,y=fit)) + 
  geom_line() + 
  geom_ribbon(aes(ymin=lci, ymax=uci, linetype=NA), alpha=0.2, show.legend=FALSE) + 
  theme(legend.title=element_blank(), legend.position=c(.65,.1), legend.key.size = unit(.25, "cm")) +
  xlab(bquote('Roads (m/km'^2*')')) + 
  ylab(bquote('Predicted MD Harvest/100km'^~2))

## Plot Grid
#Plot the plots
grid <- plot_grid(fit.canopy, fit.effort, fit.ndvi, fit.road, fit.rough, ncol=2, labels = "AUTO")


### VALIDATION OF THE MODEL
###2014
install.packages("plyr")
library(plyr)
library(ggplot2)

train <- subset(df, year != "2014")
test <- subset(df, year == "2014")

Fold <- glm(deer_100km2 ~ canopy + effort + ndviamp + ndviamp_sq + log(rough+0.001) +
              road + road_sq, family = Gamma, data=train, na.action=na.fail)
  
####  start
pred <- predict(Fold, newdata = test, se.fit = T, type="response",conf.type="mean")

result <- cor.test(pred$fit, test$deer_100km2, method="spearman", exact = FALSE)

r2014 <- result$estimate
p2014 <- result$p.value

counts <- as.data.frame(cbind(pred$fit, test$deer_100km2))

ggplot(counts, aes(y = counts[,1], x = counts[,2])) + geom_point()

####2015
train <- subset(df, year != "2015")
test <- subset(df, year == "2015")


####  start
pred <- predict(Fold, newdata = test, se.fit = T, type="response",conf.type="mean")

result <- cor.test(pred$fit, test$deer_100km2, method="spearman", exact = FALSE)

r2015 <- result$estimate
p2015 <- result$p.value

counts <- as.data.frame(cbind(pred$fit, test$deer_100km2))

ggplot(counts, aes(y = counts[,1], x = counts[,2])) + geom_point()


####2016
train <- subset(df, year != "2016")
test <- subset(df, year == "2016")

####  start
pred <- predict(Fold, newdata = test, se.fit = T, type="response",conf.type="mean")

result <- cor.test(pred$fit, test$deer_100km2, method="spearman", exact = FALSE)

r2016 <- result$estimate
p2016 <- result$p.value

counts <- as.data.frame(cbind(pred$fit, test$deer_100km2))

ggplot(counts, aes(y = counts[,1], x = counts[,2])) + geom_point()


kfold <- (r2014 +r2015 +r2016)/3
pave <- (p2014 + p2015 + p2016)/3

kfold
pave

# null H tests that there is no rank cor between the vectors.

county <- read.csv("County_Predict.csv", header = T, sep= ",")

county$ndviamp_sq <- county$ndviamp^2
county$road_sq <- county$road^2

# predict X*Beta, i.e., linear predictor for the county level
pred=predict(top.orig, newdata = county ,se.fit = T,
              type="response",conf.type="mean")

# add predictions and se's to prediction data.frame
county$fit=pred$fit
county$se.fit=pred$se.fit
county$lci=pred$fit-1.96*pred$se.fit
county$uci=pred$fit+1.96*pred$se.fit

write.csv(county, file = "county.csv")