#####
#####
#####read in the appropriate libraries
library(AICcmodavg)
library(lme4)
library(lattice)
library(mgcv)
library(gamm4)
library(car)
library(cvTools)
library(influence.ME)
library(sp)
library(rgdal)
library(raster)
library(maptools)
library(ape)
library(geosphere)
library(geoR)
library(faraway)
library(PresenceAbsence)

#####
#####
#####read in data
Ponded <- read.csv("BinomialNonRowcropsExtremeChange3.csv", header = TRUE)
head(Ponded)

#####
#####
#####log transform predictor variables with outliers and no zero values
Ponded$log.Pitdist <- log(Ponded$Pitdist)
Ponded$log.NewfootPA <- log(Ponded$NewfootPA)
Ponded$log.lsprsPostChange <- log(Ponded$lsprsPostChange)
Ponded$log.wprsPostChange <- log(Ponded$wprsPostChange)

#####
#####
#####scale all predictor variables, including those that have been log-transformed
Ponded$scaled.log.Pitdist <- scale(Ponded$log.Pitdist, center = TRUE, scale = TRUE)
Ponded$scaled.log.NewfootPA <- scale(Ponded$log.NewfootPA, center = TRUE, scale = TRUE)
Ponded$scaled.espro2 <- scale(Ponded$espro2PostChange, center = TRUE, scale = TRUE)
Ponded$scaled.log.lsprs <- scale(Ponded$log.lsprsPostChange, center = TRUE, scale = TRUE)
Ponded$scaled.fatmaxm <- scale(Ponded$fatmaxmPostChange, center = TRUE, scale = TRUE)
Ponded$scaled.log.wprs <- scale(Ponded$log.wprsPostChange, center = TRUE, scale = TRUE)
Ponded$scaled.wtmaxm <- scale(Ponded$wtmaxmPostChange, center = TRUE, scale = TRUE)
Ponded$scaled.wtmaxsf <- scale(Ponded$wtmaxsfPostChange, center = TRUE, scale = TRUE)

#####
#####
#####read in nominal variables as factors
Ponded$Wetland = factor(Ponded$Wetland)
Ponded$Wettype = factor(Ponded$Wettype)
Ponded$Year = factor(Ponded$Year)

#####
#####
#####transform wetland type to binomial
Ponded$SemiPerm <- Ponded$Wettype == 'Semi-permanent'
Ponded$Seasonal <- Ponded$Wettype == 'Seasonal'
Ponded$Temporary <- Ponded$Wettype == 'Temporary'

for(i in 1:length(Ponded$SemiPerm)){
	Ponded$SemiPerm[i] <- recode(Ponded$SemiPerm[i], "c('TRUE') = 1; ('FALSE') = 0")
}
for(i in 1:length(Ponded$Seasonal)){
	Ponded$Seasonal[i] <- recode(Ponded$Seasonal[i], "c('TRUE') = 1; ('FALSE') = 0")
}
for(i in 1:length(Ponded$Temporary)){
	Ponded$Temporary[i] <- recode(Ponded$Temporary[i], "c('TRUE') = 1; ('FALSE') = 0")
}

#####
#####
#####Make manual predictions with best-supported generalized linear mixed model
#with original data
Ponded$Pred.Presence.ManualFixedOnly <- NA
for(i in 1:nrow(Ponded)){
	Ponded$Pred.Presence.ManualFixedOnly[[i]] = 1.96446+ 
	+ 0.71939 * Ponded$SemiPerm[[i]]+ 
	- 0.87143 * Ponded$Temporary[[i]]+
	- ((1.76324 * Ponded$log.NewfootPA[[i]] - 1.76324 * -4.20119) / 0.62257)+
	+ ((0.42099 * Ponded$log.lsprsPostChange[[i]] - 0.42099 * 5.38196) / 0.31275)+
	- ((0.60281 * Ponded$fatmaxmPostChange[[i]] - 0.60281 * 16.08833) / 1.12164)+
	+ ((0.16996 * Ponded$log.wprsPostChange[[i]] - 0.16996 * 4.23183) / 0.79371)+
	- ((0.92181 * Ponded$wtmaxmPostChange[[i]] - 0.92181 * 6.53900) / 1.42316)
}

Ponded$Pred.Presence.ManualFixedOnly <- ilogit(Ponded$Pred.Presence.ManualFixedOnly)

head(Ponded)

#####
#####
#####assign inundation/non-inundation at probability threshold of 0.73
Ponded$Pred.Presence.ManualFixedOnlyPresence <- NA
for(x in 1:nrow(Ponded)){
	if(Ponded$Pred.Presence.ManualFixedOnly[x] < 0.73){
		Ponded$Pred.Presence.ManualFixedOnlyPresence[x] <- 0
	}else
		Ponded$Pred.Presence.ManualFixedOnlyPresence[x] <- 1
}

head(Ponded)

#determine number of ponding events
nrow(subset(Ponded, Pred.Presence.ManualFixedOnlyPresence == 1))

#obtain proportion of wetlands ponded
nrow(subset(Ponded, Pred.Presence.ManualFixedOnlyPresence == 1)) / nrow(Ponded)

#obtain mean probability of inundation
mean(Ponded$Pred.Presence.ManualFixedOnly)

#####
#####
#####subset wetlands that were ponded under the Modest Change Scenario
PresentExtremeChange <- subset(Ponded, Pred.Presence.ManualFixedOnlyPresence == 1)

#write to csv file
write.csv(PresentExtremeChange, "BinomialNonRowcropsExtremeChangePresent3.csv")