#####
#####
#####read in 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(gstat)
library(MuMIn)

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

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

#####
#####
#####eliminate duplicate values among years
Ponded2 <- subset(Ponded, !duplicated(Ponded$Wetland))
nrow(Ponded2)
colnames(Ponded2)

#obtain mean values for weather variables among years
Wetlands <- unique(Ponded$Wetland)
for(i in 1:nrow(Ponded2)){
	Wetland <- subset(Ponded, Wetland == Wetlands[i])
	Ponded2$espro2[i] <- mean(Wetland$espro2)
	Ponded2$lsrps[i] <- mean(Wetland$lsprs)
	Ponded2$fatmaxm[i] <- mean(Wetland$fatmaxm)
	Ponded2$wprs[i] <- mean(Wetland$wprs)
	Ponded2$wtmaxm[i] <- mean(Wetland$wtmaxm)
	Ponded2$wtmaxsf[i] <- mean(Wetland$wtmaxsf)
} 

#####
#####
#####write new data frames and csv files with assumed climate change scenarios
#Modest Change Scenario
PondedModest <- Ponded2
PondedModest$espro2PostChange <- PondedModest$espro2 + 0.2 * PondedModest$espro2
PondedModest$lsprsPostChange <- PondedModest$lsprs
PondedModest$fatmaxmPostChange <- PondedModest$fatmaxm + 2
PondedModest$wprsPostChange <- PondedModest$wprs + 0.1 * PondedModest$wprs
PondedModest$wtmaxmPostChange <- PondedModest$wtmaxm + 3
PondedModest$wtmaxsfPostChange <- PondedModest$wtmaxsf - 0.1 * PondedModest$wtmaxsf

#Extreme Change Scenario
PondedExtreme <- Ponded2
PondedExtreme$espro2PostChange <- PondedExtreme$espro2 + 0.4 * PondedExtreme$espro2
PondedExtreme$lsprsPostChange <- PondedExtreme$lsprs - 0.2 * PondedExtreme$lsprs
PondedExtreme$fatmaxmPostChange <- PondedExtreme$fatmaxm + 4
PondedExtreme$wprsPostChange <- PondedExtreme$wprs + 0.3 * PondedExtreme$wprs
PondedExtreme$wtmaxmPostChange <- PondedExtreme$wtmaxm + 5
PondedExtreme$wtmaxsfPostChange <- PondedExtreme$wtmaxsf - 0.2 * PondedExtreme$wtmaxsf

write.csv(PondedModest, "GaussianNoRowcropsModestChange3.csv")
write.csv(PondedExtreme, "GaussianNoRowcropsExtremeChange3.csv")

#####
#####
#####check for outliers in variable
boxplot(Ponded$Squaremeters)
dotchart(Ponded$Squaremeters)
dotchart(Ponded$Squaremeters, group = Ponded$Wettype)
dotchart(Ponded$Squaremeters, group = Ponded$agtype)
dotchart(Ponded$Squaremeters, group = Ponded$Year)
#outliers in ponded square meters

boxplot(Ponded$Pitdist)
dotchart(Ponded$Pitdist)
dotchart(Ponded$Pitdist, group = Ponded$Wettype)
dotchart(Ponded$Pitdist, group = Ponded$agtype)
dotchart(Ponded$Pitdist, group = Ponded$Year)
#outliers in pit distance

boxplot(Ponded$NewfootPA)
dotchart(Ponded$NewfootPA)
dotchart(Ponded$NewfootPA, group = Ponded$Wettype)
dotchart(Ponded$NewfootPA, group = Ponded$agtype)
dotchart(Ponded$NewfootPA, group  = Ponded$Year)
#outliers in perimeter to area ratio, with semi-permanent wetlands having
#much lower ratios

boxplot(Ponded$espro2)
dotchart(Ponded$espro2)
dotchart(Ponded$espro2, group = Ponded$Wettype)
dotchart(Ponded$espro2, group = Ponded$agtype)
dotchart(Ponded$espro2, group = Ponded$Year)

boxplot(Ponded$esprs)
dotchart(Ponded$esprs)
dotchart(Ponded$esprs, group = Ponded$Wettype)
dotchart(Ponded$esprs, group = Ponded$agtype)
dotchart(Ponded$esprs, group = Ponded$Year)

boxplot(Ponded$estmaxm)
dotchart(Ponded$estmaxm)
dotchart(Ponded$estmaxm, group = Ponded$Wettype)
dotchart(Ponded$estmaxm, group = Ponded$agtype)
dotchart(Ponded$estmaxm, group = Ponded$Year)

boxplot(Ponded$estminm)
dotchart(Ponded$estminm)
dotchart(Ponded$estminm, group = Ponded$Wettype)
dotchart(Ponded$estminm, group = Ponded$agtype)
dotchart(Ponded$estminm, group = Ponded$Year)

boxplot(Ponded$esvpdm)
dotchart(Ponded$esvpdm)
dotchart(Ponded$esvpdm, group = Ponded$Wettype)
dotchart(Ponded$esvpdm, group = Ponded$agtype)
dotchart(Ponded$esvpdm, group = Ponded$Year)
#there are outliers in spring mean vapor pressure deficit

boxplot(Ponded$lspro2)
dotchart(Ponded$lspro2)
dotchart(Ponded$lspro2, group = Ponded$Wettype)
dotchart(Ponded$lspro2, group = Ponded$agtype)
dotchart(Ponded$lspro2, group = Ponded$Year)

boxplot(Ponded$lsprs)
dotchart(Ponded$lsprs)
dotchart(Ponded$lsprs, group = Ponded$Wettpye)
dotchart(Ponded$lsprs, group = Ponded$agtype)
dotchart(Ponded$lsprs, group = Ponded$Year)
#outliers in summer precipitation sum

boxplot(Ponded$lstmaxm)
dotchart(Ponded$lstmaxm)
dotchart(Ponded$lstmaxm, group = Ponded$Wettype)
dotchart(Ponded$lstmaxm, group = Ponded$agtype)
dotchart(Ponded$lstmaxm, group = Ponded$Year)

boxplot(Ponded$lsvpdm)
dotchart(Ponded$lsvpdm)
dotchart(Ponded$lsvpdm, group = Ponded$Wettype)
dotchart(Ponded$lsvpdm, group = Ponded$agtype)
dotchart(Ponded$lsvpdm, group = Ponded$Year)

boxplot(Ponded$fapro1)
dotchart(Ponded$fapro1)
dotchart(Ponded$fapro1, group = Ponded$Wettype)
dotchart(Ponded$fapro1, group = Ponded$agtype)
dotchart(Ponded$fapro1, group = Ponded$Year)
#outliers are present in number of major autumn precipitation events

boxplot(Ponded$faprs)
dotchart(Ponded$faprs)
dotchart(Ponded$faprs, group = Ponded$Wettype)
dotchart(Ponded$faprs, group = Ponded$agtype)
dotchart(Ponded$faprs, group = Ponded$Year)
#outliers in fall precipitation sum

boxplot(Ponded$fatmaxm)
dotchart(Ponded$fatmaxm)
dotchart(Ponded$fatmaxm, group = Ponded$Wettype)
dotchart(Ponded$fatmaxm, group = Ponded$agtype)
dotchart(Ponded$fatmaxm, group = Ponded$Year)

boxplot(Ponded$fatminm)
dotchart(Ponded$fatminm)
dotchart(Ponded$fatminm, group = Ponded$Wettype)
dotchart(Ponded$fatminm, group = Ponded$agtype)
dotchart(Ponded$fatminm, group = Ponded$Year)

boxplot(Ponded$favpdm)
dotchart(Ponded$favpdm)
dotchart(Ponded$favpdm, group = Ponded$Wettype)
dotchart(Ponded$favpdm, group = Ponded$agtype)
dotchart(Ponded$favpdm, group = Ponded$Year)

boxplot(Ponded$wpro1)
dotchart(Ponded$wpro1)
dotchart(Ponded$wpro1, group = Ponded$Wettype)
dotchart(Ponded$wpro1, group = Ponded$agtype)
dotchart(Ponded$wpro1, group = Ponded$Year)
#outliers are present in number of major winter precipitation events

boxplot(Ponded$wprs)
dotchart(Ponded$wprs)
dotchart(Ponded$wprs, group = Ponded$Wettype)
dotchart(Ponded$wprs, group = Ponded$agtype)
dotchart(Ponded$wprs, group = Ponded$Year)
#outliers do exist in winter precipitation sum

boxplot(Ponded$wtmaxm)
dotchart(Ponded$wtmaxm)
dotchart(Ponded$wtmaxm, group = Ponded$Wettype)
dotchart(Ponded$wtmaxm, group = Ponded$agtype)
dotchart(Ponded$wtmaxm, group = Ponded$Year)

boxplot(Ponded$wtmaxsf)
dotchart(Ponded$wtmaxsf)
dotchart(Ponded$wtmaxsf, group = Ponded$Wettype)
dotchart(Ponded$wtmaxsf, group = Ponded$agtype)
dotchart(Ponded$wtmaxsf, group = Ponded$Year)

boxplot(Ponded$wtminm)
dotchart(Ponded$wtminm)
dotchart(Ponded$wtminm, group = Ponded$Wettype)
dotchart(Ponded$wtminm, group = Ponded$agtype)
dotchart(Ponded$wtminm, group = Ponded$Year)

boxplot(Ponded$wtminsf)
dotchart(Ponded$wtminsf)
dotchart(Ponded$wtminsf, group = Ponded$Wettype)
dotchart(Ponded$wtminsf, group = Ponded$agtype)
dotchart(Ponded$wtminsf, group = Ponded$Year)

boxplot(Ponded$ytminaf)
dotchart(Ponded$ytminaf)
dotchart(Ponded$ytminaf, group = Ponded$Wettype)
dotchart(Ponded$ytminaf, group = Ponded$agtype)
dotchart(Ponded$ytminaf, group = Ponded$Year)

boxplot(Ponded$ytminsf)
dotchart(Ponded$ytminsf)
dotchart(Ponded$ytminsf, group = Ponded$Wettype)
dotchart(Ponded$ytminsf, group = Ponded$agtype)
dotchart(Ponded$ytminsf, group = Ponded$Year)

#####
#####
#####log transform variables with outliers
Ponded$log.Meters <- log(Ponded$Squaremeters)
Ponded$log.Pitdist <- log(Ponded$Pitdist)
Ponded$log.NewfootPA <- log(Ponded$NewfootPA)
Ponded$log.esvpdm <- log(Ponded$esvpdm)
Ponded$log.lsprs <- log(Ponded$lsprs)
#Ponded$log.fapro1 <- log(Ponded$fapro1 + 0.5)
Ponded$log.faprs <- log(Ponded$faprs)
#Ponded$log.wpro1 <- log(Ponded$wpro1 + 0.5)
Ponded$log.wprs <- log(Ponded$wprs)

#####
#####
#####scale all predictor variables, including those that have already 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$espro2, center = TRUE, scale = TRUE)
Ponded$scaled.esprs <- scale(Ponded$esprs, center = TRUE, scale = TRUE)
Ponded$scaled.estmaxm <- scale(Ponded$estmaxm, center = TRUE, scale = TRUE)
Ponded$scaled.estminm <- scale(Ponded$estminm, center = TRUE, scale = TRUE)
Ponded$scaled.log.esvpdm <- scale(Ponded$esvpdm, center = TRUE, scale = TRUE)
Ponded$scaled.lspro2 <- scale(Ponded$lspro2, center = TRUE, scale = TRUE)
Ponded$scaled.log.lsprs <- scale(Ponded$log.lsprs, center = TRUE, scale = TRUE)
Ponded$scaled.lstmaxm <- scale(Ponded$lstmaxm, center = TRUE, scale = TRUE)
Ponded$scaled.lsvpdm <- scale(Ponded$lsvpdm, center = TRUE, scale = TRUE)
Ponded$scaled.fapro1 <- scale(Ponded$fapro1, center = TRUE, scale = TRUE)
Ponded$scaled.log.faprs <- scale(Ponded$log.faprs, center = TRUE, scale = TRUE)
Ponded$scaled.fatmaxm <- scale(Ponded$fatmaxm, center = TRUE, scale = TRUE)
Ponded$scaled.fatminm <- scale(Ponded$fatminm, center = TRUE, scale = TRUE)
Ponded$scaled.favpdm <- scale(Ponded$favpdm, center = TRUE, scale = TRUE)
Ponded$scaled.wpro1 <- scale(Ponded$wpro1, center = TRUE, scale = TRUE)
Ponded$scaled.log.wprs <- scale(Ponded$log.wprs, center = TRUE, scale = TRUE)
Ponded$scaled.wtmaxm <- scale(Ponded$wtmaxm, center = TRUE, scale = TRUE)
Ponded$scaled.wtmaxsf <- scale(Ponded$wtmaxsf, center = TRUE, scale = TRUE)
Ponded$scaled.wtminm <- scale(Ponded$wtminm, center = TRUE, scale = TRUE)
Ponded$scaled.wtminsf <- scale(Ponded$wtminsf, center = TRUE, scale = TRUE)
Ponded$scaled.wvpdm <- scale(Ponded$wvpdm, center = TRUE, scale = TRUE)
Ponded$scaled.ytminaf <- scale(Ponded$ytminaf, center = TRUE, scale = TRUE)
Ponded$scaled.ytminsf <- scale(Ponded$ytminsf, center = TRUE, scale = TRUE)

#####
#####
#####insert script for setting up fancy pairplots
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}

panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

#####
#####
#####construct pairplots to look at correlations between variables in same season
pairs(~ log.Meters + scaled.log.Pitdist + scaled.log.NewfootPA, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.espro2 + scaled.esprs + scaled.estmaxm + scaled.estminm + scaled.log.esvpdm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.espro2 + scaled.estmaxm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.lspro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.lsvpdm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.log.lsprs + scaled.lstmaxm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.fapro1 + scaled.log.faprs + scaled.fatmaxm + scaled.fatminm + scaled.favpdm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.fatmaxm + scaled.log.faprs, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.wpro1 + scaled.log.wprs + scaled.wtmaxm + scaled.wtmaxsf + scaled.wtminm + scaled.wtminsf + scaled.wvpdm, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.wtmaxsf + scaled.log.wprs, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.espro2 + scaled.estmaxm + scaled.log.lsprs + scaled.lstmaxm + scaled.log.faprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)
pairs(~ log.Meters + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded)

#####
#####
#####compare alternative random effects structures as a set to get delta AICc scores and weights
random.set <- list((log.Meters ~ Wettype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year)),
(log.Meters ~ Wettype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland)),
(log.Meters ~ Wettype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Year)))

test.random <- lapply(random.set, lmer, REML = TRUE, data = Ponded)
random.table <- aictab(test.random, modnames = as.character(random.set), sort = TRUE, second.ord = TRUE)
random.table

top.random <- lmer(log.Meters ~ Wettype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland), REML = TRUE, data = Ponded)
summary(top.random)

#####
#####
#####use backwards selection to remove unnecessary terms
Backwards1 <- drop1(top.random, test = "Chisq", trace = TRUE, REML = FALSE)
Backwards1

top.model <- lmer(log.Meters ~ Wettype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.fatmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland), REML = TRUE, data = Ponded)
summary(top.model)

#####
#####
#####obtain pseudo r-squared values for the best-supported model
r.squaredGLMM(top.model)

#####
#####
#####use plots to identify linear model assumption violations
hist(resid(top.model))
plot(fitted(top.model), resid(top.model))
xyplot(resid(top.model) ~ fitted(top.model), type = "smooth")
plot(resid(top.model) ~ Ponded$Wettype)
plot(resid(top.model) ~ Ponded$Year)
plot(resid(top.model) ~ Ponded$scaled.log.NewfootPA)
plot(resid(top.model) ~ Ponded$scaled.log.lsprs)
plot(resid(top.model) ~ Ponded$scaled.fatmaxm)
plot(resid(top.model) ~ Ponded$scaled.log.wprs)
plot(resid(top.model) ~ Ponded$scaled.wtmaxsf)

#####
#####
#####plot residuals against variables not included in the model, to identify any relationships
plot(resid(top.model) ~ Ponded$scaled.esprs)
plot(resid(top.model) ~ Ponded$scaled.estmaxm)
plot(resid(top.model) ~ Ponded$scaled.estminm)
plot(resid(top.model) ~ Ponded$scaled.log.esvpdm)
plot(resid(top.model) ~ Ponded$scaled.lstmaxm)
plot(resid(top.model) ~ Ponded$scaled.lsvpdm)
plot(resid(top.model) ~ Ponded$scaled.fapro1)
plot(resid(top.model) ~ Ponded$scaled.log.faprs)
plot(resid(top.model) ~ Ponded$scaled.fatminm)
plot(resid(top.model) ~ Ponded$scaled.favpdm)
plot(resid(top.model) ~ Ponded$scaled.wpro1)
plot(resid(top.model) ~ Ponded$scaled.wtmaxm)
plot(resid(top.model) ~ Ponded$scaled.wtminm)
plot(resid(top.model) ~ Ponded$scaled.wtminsf)
plot(resid(top.model) ~ Ponded$scaled.wvpdm)
plot(resid(top.model) ~ Ponded$scaled.ytminsf)
plot(resid(top.model) ~ Ponded$scaled.ytminaf)

#####
#####
#####experiment with computing confidence intervals
#manually
Coefficients <- coef(summary(top.model))
Lowers <- Coefficients[,1] - (1.96 * Coefficients[,2])
Uppers <- Coefficients[,1] + (1.96 * Coefficients[,2])
ConfInts <- rbind(Coefficients[,1], Lowers, Uppers)
ConfInts

#####
#####
#####make predictions with the best supported linear mixed model
#with only fixed effects considered
Ponded$Pred.Log.Meters.LMM.FixedOnly <- predict(top.model, re.form = ~0, type = 'response')
Ponded$Pred.Meters.LMM.FixedOnly <- exp(Ponded$Pred.Log.Meters.LMM.FixedOnly)
cor(Ponded$Squaremeters, Ponded$Pred.Meters.LMM.FixedOnly, method = "pearson")
xyplot(Ponded$Squaremeters ~ Ponded$Pred.Meters.LMM.FixedOnly, type = "smooth")

#####
#####
#####generate statistics
#Overall predicted
mean(Ponded$Pred.Meters.LMM.FixedOnly)
sum(Ponded$Pred.Meters.LMM.FixedOnly)

#Overall observed
mean(Ponded$Squaremeters)
sum(Ponded$Squaremeters)

#Divided by years
Year2004 <- subset(Ponded, Year == 2004)
mean(Year2004$Pred.Meters.LMM.FixedOnly)
sum(Year2004$Pred.Meters.LMM.FixedOnly)
mean(Year2004$Squaremeters)
sum(Year2004$Squaremeters)

Year2006 <- subset(Ponded, Year == 2006)
mean(Year2006$Pred.Meters.LMM.FixedOnly)
sum(Year2006$Pred.Meters.LMM.FixedOnly)
mean(Year2006$Squaremeters)
sum(Year2006$Squaremeters)

Year2007 <- subset(Ponded, Year == 2007)
mean(Year2007$Pred.Meters.LMM.FixedOnly)
sum(Year2007$Pred.Meters.LMM.FixedOnly)
mean(Year2007$Squaremeters)
sum(Year2007$Squaremeters)

Year2008 <- subset(Ponded, Year == 2008)
mean(Year2008$Pred.Meters.LMM.FixedOnly)
sum(Year2008$Pred.Meters.LMM.FixedOnly)
mean(Year2008$Squaremeters)
sum(Year2008$Squaremeters)

Year2009 <- subset(Ponded, Year == 2009)
mean(Year2009$Pred.Meters.LMM.FixedOnly)
sum(Year2009$Pred.Meters.LMM.FixedOnly)
mean(Year2009$Squaremeters)
sum(Year2009$Squaremeters)

#####
#####
#####make predictions with the best supported linear mixed model
#with fixed and random effects considered
Ponded$Pred.Log.Meters.LMM.FixedAndRandom <- predict(top.model, re.form = ~(1|Wetland) + (1|Year), type = 'response')
Ponded$Pred.Meters.LMM.FixedAndRandom <- exp(Ponded$Pred.Log.Meters.LMM.FixedAndRandom)
cor(Ponded$Squaremeters, Ponded$Pred.Meters.LMM.FixedAndRandom, method = "pearson")
xyplot(Ponded$Squaremeters ~ Ponded$Pred.Meters.LMM.FixedAndRandom, type = "smooth")

#####
#####
#####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 predictions manually, for comparative purposes
Ponded$Pred.Log.Meters.ManualFixedOnly <- NA
for(i in 1:nrow(Ponded)){
	Ponded$Pred.Log.Meters.ManualFixedOnly[[i]] = 8.25802+
	+ 1.22951 * Ponded$SemiPerm[[i]]+
	- 0.57479 * Ponded$Temporary[[i]]+
	- ((0.72171 * Ponded$log.NewfootPA[[i]] - 0.72171 * -4.38943) / 0.59353)+
	+ ((0.10567 * Ponded$espro2[[i]] - 0.10567 * 0.27358) / 0.46648)+
	+ ((0.26324 * Ponded$log.lsprs[[i]] - 0.26324 * 5.42392) / 0.29450)+
	- ((0.53600 * Ponded$fatmaxm[[i]] - 0.53600 * 15.97654) /  1.10901)+
	+ ((0.17197 * Ponded$log.wprs[[i]] - 0.17197 * 4.24774) / 0.79608)+
	+ ((0.63935 * Ponded$wtmaxsf[[i]] - 0.63935 *  25.65612) / 4.46749)
}

#####
#####
#####calculate the inverse log for predicted values
Ponded$Pred.Meters.ManualFixedOnly <- exp(Ponded$Pred.Log.Meters.ManualFixedOnly)
cor(Ponded$Squaremeters, Ponded$Pred.Meters.ManualFixedOnly, method = "pearson")
xyplot(Ponded$Squaremeters ~ Ponded$Pred.Meters.ManualFixedOnly, type = "smooth")

#####
#####
#####use k-fold cross validation function from Chris Jorgensen for validation
#with k = 10
source("cv.glmmFunction_v2.0.R")
Presence.LMM.CrossValidation <- cv.glmm(Ponded, top.model, k = 10, cvData = TRUE)
str(Presence.LMM.CrossValidation)
Presence.LMM.CrossValidation
head(Presence.LMM.CrossValidation$PredObsData)
summary(Presence.LMM.CrossValidation)
Presence.LMM.CrossValidation$delta
Presence.LMM.CrossValidation$CV

#read in validation data for all wetlands
ValidationAllWetlands <- read.csv("FloodedAreaValidationAllWetlands.csv", header = TRUE)
head(ValidationAllWetlands)

#get PCC for non-rowcrop-embedded wetlands
PredictedNonRowcrops <- Presence.LMM.CrossValidation$PredObsData[1]
ObservedNonRowcrops <- Presence.LMM.CrossValidation$PredObsData[2]
cor(PredictedNonRowcrops$pred, ObservedNonRowcrops$obs, method = "pearson")

#write validation data to a .csv file
ValidationDataNonRowcrops <- data.frame(Predicted = PredictedNonRowcrops$pred)
ValidationDataNonRowcrops$Observed <- ObservedNonRowcrops$obs
write.csv(ValidationDataNonRowcrops, "FloodedAreaValidationNonRowcrops.csv")

#plot validation results
#DummyLine <- data.frame(X = c(0,15))
#DummyLine$Y <- c(0,15)

par(mar = c(5, 5, 4, 2), oma = c(0, 1, 0, 0))
plot(PredictedNonRowcrops$pred, ObservedNonRowcrops$obs, pch = 6, cex = 0.005, col = "blue", frame.plot = TRUE, tck = 0.02, xlim = c(4, 15), ylim = c(-10, 15), xlab = expression(paste("Predicted log-transformed ponded area (", m^2, ")" )), ylab = expression(paste("Observed log-transformed ponded area (", m^2, ")" )), cex.lab = 1.25)
abline(lm(ObservedNonRowcrops$obs ~ PredictedNonRowcrops$pred), col = "blue", lwd = 2)
points(ValidationAllWetlands$Predicted, ValidationAllWetlands$Observed, pch = 6, cex = 0.005, col = "red")
abline(lm(ValidationAllWetlands$Observed ~ ValidationAllWetlands$Predicted), col = "red", lwd = 2, lty = 2)
#lines(DummyLine$X, DummyLine$Y, lty = 3, lwd = 2, col = "black")
legend(5, -3, c("Non-rowcrop-embedded (PCC = 0.58)", "All (PCC = 0.49)"), cex = 1.0, col = c("blue", "red"), lty = c(1, 2), lwd = c(2, 2), bty = 'n')

#####
#####
#####determine if the incorporation of smoothing terms via additive mixed modeling is appropriate
#not subset by any nominal variable
xyplot(resid(top.model) ~ Ponded$scaled.log.NewfootPA, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.log.lsprs, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.fatmaxm, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.wtmaxm, type = 'smooth')

#subset by wetland type
xyplot(resid(top.model) ~ Ponded$scaled.log.NewfootPA|Ponded$Wettype, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.log.lsprs|Ponded$Wettype, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.fatmaxm|Ponded$Wettype, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.wtmaxm|Ponded$Wettype, type = 'smooth')

#subset by year
xyplot(resid(top.model) ~ Ponded$scaled.log.NewfootPA|Ponded$Year, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.log.lsprs|Ponded$Year, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.fatmaxm|Ponded$Year, type = 'smooth')
xyplot(resid(top.model) ~ Ponded$scaled.wtmaxm|Ponded$Year, type = 'smooth')

#####
#####
#####use k-fold cross validation function from Chris Jorgensen for validation
#with k = 10
source("cv.glmmFunction_v2.0.R")
Presence.LMM.CrossValidation <- cv.glmm(Ponded, top.model, k = 10, cvData = TRUE)
str(Presence.LMM.CrossValidation)
head(Presence.LMM.CrossValidation$PredObsData)
summary(Presence.LMM.CrossValidation)
Presence.LMM.CrossValidation$delta
Presence.LMM.CrossValidation$CV

#get PCC between predicted and observed flooded acres
cor(Presence.LMM.CrossValidation$PredObsData[,2], Presence.LMM.CrossValidation$PredObsData[,1], method = "pearson")
xyplot(Presence.LMM.CrossValidation$PredObsData[,2] ~ Presence.LMM.CrossValidation$PredObsData[,1], type = "smooth")