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

#Tell R where to find the data
setwd('E:/documents/420 Project')

#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 = ",")

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

#create Correlation matrix
ggcorr(df[c(35:47)],
       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 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_100km ~",elev_models$cov.string[1]), data=df)
mod.list.elev[[2]]<-glm(paste("deer_100km ~",elev_models$cov.string[2]), data=df)
mod.list.elev[[3]]<-glm(paste("deer_100km ~",elev_models$cov.string[3]), data=df)

#generate AIC table

install.packages(AICcmodavg)
library(AICcmodavg)
install.packages(XLConnect)
library(XLConnect)
install.packages(xlsx)
library(xlsx)
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)

install.packages(XLConnect)
library(XLConnect)

write.xlsx2(elev_aic,"func_forms.xlsx", sheetName = "elev_aic",)
writeWorksheetToFile(file = "func_forms.xlsx", data = elev_tab, sheet = "elev_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_100km ~",slope_models$cov.string[1]), data=df)
mod.list.slope[[2]]<-glm(paste("deer_100km ~",slope_models$cov.string[2]), data=df)
mod.list.slope[[3]]<-glm(paste("deer_100km ~",slope_models$cov.string[3]), 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_100km ~",rough_models$cov.string[1]), data=df)
mod.list.rough[[2]]<-glm(paste("deer_100km ~",rough_models$cov.string[2]), data=df)
mod.list.rough[[3]]<-glm(paste("deer_100km ~",rough_models$cov.string[3]), 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_100km ~",ag_models$cov.string[1]), data=df)
mod.list.ag[[2]]<-glm(paste("deer_100km ~",ag_models$cov.string[2]), data=df)
mod.list.ag[[3]]<-glm(paste("deer_100km ~",ag_models$cov.string[3]), 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_100km ~",canopy_models$cov.string[1]), data=df)
mod.list.canopy[[2]]<-glm(paste("deer_100km ~",canopy_models$cov.string[2]), 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_100km ~",forest_models$cov.string[1]), data=df)
mod.list.forest[[2]]<-glm(paste("deer_100km ~",forest_models$cov.string[2]), 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_100km ~",ndviamp_models$cov.string[1]), data=df)
mod.list.ndviamp[[2]]<-glm(paste("deer_100km ~",ndviamp_models$cov.string[2]), data=df)
mod.list.ndviamp[[3]]<-glm(paste("deer_100km ~",ndviamp_models$cov.string[3]), 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_100km ~",ndvitin_models$cov.string[1]), data=df)
mod.list.ndvitin[[2]]<-glm(paste("deer_100km ~",ndvitin_models$cov.string[2]), data=df)
mod.list.ndvitin[[3]]<-glm(paste("deer_100km ~",ndvitin_models$cov.string[3]), 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_100km ~",develop_models$cov.string[1]), data=df)
mod.list.develop[[2]]<-glm(paste("deer_100km ~",develop_models$cov.string[2]), 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_100km ~",range_models$cov.string[1]), data=df)
mod.list.range[[2]]<-glm(paste("deer_100km ~",range_models$cov.string[2]), data=df)
mod.list.range[[3]]<-glm(paste("deer_100km ~",range_models$cov.string[3]), 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_100km ~",riparian_models$cov.string[1]), data=df)
mod.list.riparian[[2]]<-glm(paste("deer_100km ~",riparian_models$cov.string[2]), 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_100km ~",road_models$cov.string[1]), data=df)
mod.list.road[[2]]<-glm(paste("deer_100km ~",road_models$cov.string[2]), data=df)
mod.list.road[[3]]<-glm(paste("deer_100km ~",road_models$cov.string[3]), 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_100km ~",urban_models$cov.string[1]), data=df)
mod.list.urban[[2]]<-glm(paste("deer_100km ~",urban_models$cov.string[2]), 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_100km ~ urban_ps_std + elev_std + elev_std_sq + slope_std + slope_std_sq + rough_std +
              rough_std_sq + ag_std + ag_std_sq + canopy_ps_std + forest_std + forest_ps_std + ndviamp_std +
              ndvitin_std +ndvitin_std_sq + develop_ps_std + range_std + range_std_sq +riparian_ps_std + road_std +
              road_std_sq, data=df, na.action=na.fail)
library("MuMIn")
mods <- dredge(full, subset = !(forest_std && forest_ps_std) && !(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_std_sq) && 
               !(slope_std_sq && rough_std_sq) && !(rough_std_sq && 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_ps_std && forest_std) && !( canopy_ps_std && forest_ps_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_ps_std && road_std_sq) && !(develop_ps_std && urban_ps_std) &&
               !(range_std_sq && road_std_sq) && !(road_std_sq && urban_ps_std) && !(elev_std_sq && canopy_ps_std) &&
               c(dc(elev_std, elev_std_sq), dc(elev_std_sq, elev_std), dc(slope_std, slope_std_sq), dc(rough_std, rough_std_sq),
               dc(rough_std_sq, rough_std), dc(ag_std, ag_std_sq), dc(ag_std_sq, ag_std), dc(ndvitin_std, ndvitin_std_sq), 
               dc(range_std, range_std_sq), dc(range_std_sq, range_std), dc(road_std, road_std_sq), dc(road_std_sq, road_std)))

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

top <- glm(deer_100km ~ forest_std + ndvitin_std + ndvitin_std_sq + riparian_ps_std + rough_std + rough_std_sq +
             road_std + road_std_sq, data=df, na.action=na.fail)

df$ndvitin_sq <- df$ndvitin^2
df$rough_sq <- df$rough^2
df$road_sq <- df$road^2

top.orig <- glm(deer_100km ~ forest + ndvitin + ndvitin_sq + log(riparian+0.001) + rough + rough_sq +
                  road + road_sq, 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),
                         ndvitin = mean(df$ndvitin),
                         riparian = mean(df$riparian),
                         rough = mean(df$rough),
                         urban = mean(df$urban))

dat.predict1$ndvitin_sq <- dat.predict1$ndvitin^2
dat.predict1$rough_sq <- dat.predict1$rough^2

# predict X*Beta, i.e., linear predictor using across-strata means
# as reference 
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("Predicted Mule Deer Harvest/100km")

#NDVI
dat.predict2=expand.grid(canopy = mean(df$canopy),
                         ndvitin = seq(min(df$ndvitin),max(df$ndvitin),length=100),
                         riparian = mean(df$riparian),
                         rough = mean(df$rough),
                         urban = mean(df$urban))

dat.predict2$ndvitin_sq <- dat.predict2$ndvitin^2
dat.predict2$rough_sq <- dat.predict2$rough^2

# predict X*Beta, i.e., linear predictor using across-strata means
# as reference 
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=ndvitin,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("Time Integrated NDVI") + 
  ylab("Predicted Mule Deer Harvest/100km")


#riparian
dat.predict3=expand.grid(canopy = mean(df$canopy),
                         ndvitin = mean(df$ndvitin),
                         riparian = seq(min(df$riparian),max(df$riparian),length=100),
                         rough = mean(df$rough),
                         urban = mean(df$urban))

dat.predict3$ndvitin_sq <- dat.predict3$ndvitin^2
dat.predict3$rough_sq <- dat.predict3$rough^2

# predict X*Beta, i.e., linear predictor using across-strata means
# as reference 
pred3=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=pred3$fit
dat.predict3$se.fit=pred3$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

##Plotting
library(ggplot2)
require(cowplot)

fit.riparian <- ggplot(dat.predict3, aes(x=riparian,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("Riparian (%)") + 
  ylab("Predicted Mule Deer Harvest/100km")


#rough
dat.predict4=expand.grid(canopy = mean(df$canopy),
                         ndvitin = mean(df$ndvitin),
                         riparian = mean(df$riparian),
                         rough = seq(min(df$rough),max(df$rough),length=100),
                         urban = mean(df$urban))

dat.predict4$ndvitin_sq <- dat.predict4$ndvitin^2
dat.predict4$rough_sq <- dat.predict4$rough^2

# predict X*Beta, i.e., linear predictor using across-strata means
# as reference 
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

##Plotting
library(ggplot2)
require(cowplot)

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("Predicted Mule Deer Harvest/100km")


#Urban
dat.predict5=expand.grid(canopy = mean(df$canopy),
                         ndvitin = mean(df$ndvitin),
                         riparian = mean(df$riparian),
                         rough = mean(df$rough),
                         urban = seq(min(df$urban),max(df$urban),length=100))

dat.predict5$ndvitin_sq <- dat.predict5$ndvitin^2
dat.predict5$rough_sq <- dat.predict5$rough^2

# predict X*Beta, i.e., linear predictor using across-strata means
# as reference 
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

##Plotting
library(ggplot2)
require(cowplot)

fit.urban <- ggplot(dat.predict5, aes(x=urban,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("Urbanization (%)") + 
  ylab("Predicted Mule Deer Harvest/100km")

## Plot Grid
#Plot the plots
grid <- plot_grid(fit.canopy, fit.ndvi, fit.riparian, fit.rough, fit.urban, ncol=2)

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

dat.temp <- read.csv("E:/documents/420 Project/analysis_2.csv", header = TRUE)

dat.temp$ndvitin_sq <- dat.temp$ndvitin^2
dat.temp$rough_sq <- dat.temp$rough^2
dat.temp$road_sq <- dat.temp$road^2

train <- subset(dat.temp, year != "2014")
test <- subset(dat.temp, year == "2014")

Fold <- glm(deer_100km ~ forest + ndvitin + ndvitin_sq + log(riparian+0.001) + rough + rough_sq +
              road + road_sq, data=train, na.action=na.fail)
  
####  start
pred <- predict(Fold, newdata = test, se.fit = T, type="response",conf.type="mean")
f1pred <- pred$fit
f1pred[f1pred<0] <- 0

br = seq(min(f1pred)-0.01,max(c(max(f1pred),max(test$deer_100km))),by=(max(c(max(f1pred),max(test$deer_100km)))-(min(f1pred)-0.01))/10)
ranges = paste(head(br,-1), br[-1], sep="-")
freq = hist(f1pred, breaks=br, include.lowest=TRUE, plot=FALSE)
pred <- data.frame(bin = ranges, pred = freq$counts)

actual = test$deer_100km

freq2 = hist(actual, breaks=br, include.lowest=TRUE, plot=FALSE)
actual <- data.frame(bin = ranges, actual = freq2$counts)

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

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

counts <- cbind(pred, actual)

ggplot(counts, aes(y = pred, x = actual)) + geom_point()


####2015
train <- subset(dat.temp, year != "2015")
test <- subset(dat.temp, year == "2015")

Fold <- glm(deer_100km ~ forest + ndvitin + ndvitin_sq + log(riparian+0.001) + rough + rough_sq +
              road + road_sq, data=train, na.action=na.fail)

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

br = seq(min(f1pred)-0.01,max(c(max(f1pred),max(test$deer_100km))),by=(max(c(max(f1pred),max(test$deer_100km)))-(min(f1pred)-0.01))/10)
ranges = paste(head(br,-1), br[-1], sep="-")
freq = hist(f1pred, breaks=br, include.lowest=TRUE, plot=FALSE)
pred <- data.frame(bin = ranges, pred = freq$counts)

actual = test$deer_100km

freq2 = hist(actual, breaks=br, include.lowest=TRUE, plot=FALSE)
actual <- data.frame(bin = ranges, actual = freq2$counts)

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

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

counts <- cbind(pred, actual)

ggplot(counts, aes(y = pred, x = actual)) + geom_point()


####2016
train <- subset(dat.temp, year != "2016")
test <- subset(dat.temp, year == "2016")

Fold <- glm(deer_100km ~ forest + ndvitin + ndvitin_sq + log(riparian+0.001) + rough + rough_sq +
              road + road_sq, data=train, na.action=na.fail)

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

br = seq(min(f1pred)-0.01,max(c(max(f1pred),max(test$deer_100km))),by=(max(c(max(f1pred),max(test$deer_100km)))-(min(f1pred)-0.01))/10)
ranges = paste(head(br,-1), br[-1], sep="-")
freq = hist(f1pred, breaks=br, include.lowest=TRUE, plot=FALSE)
pred <- data.frame(bin = ranges, pred = freq$counts)

actual = test$deer_100km

freq2 = hist(actual, breaks=br, include.lowest=TRUE, plot=FALSE)
actual <- data.frame(bin = ranges, actual = freq2$counts)

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

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

counts <- cbind(pred, actual)

ggplot(counts, aes(y = pred, x = actual)) + geom_point()

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

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