##### ##### #####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("PondedGaussianCombined.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) Ponded$agtype = factor(Ponded$agtype) ##### ##### #####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, "GaussianTotalModestChange3.csv") write.csv(PondedExtreme, "GaussianTotalExtremeChange3.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) 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.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.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.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.log.faprs + scaled.fatminm, 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.log.wprs + scaled.wtmaxsf, diag.panel = panel.hist, upper.panel = panel.smooth, lower.panel = panel.cor, data = Ponded) pairs(~ log.Meters + scaled.ytminaf + scaled.ytminsf, 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.fatminm + 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.lstmaxm + scaled.fatminm + 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 + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.fatminm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year)), (log.Meters ~ Wettype + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.fatminm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland)), (log.Meters ~ Wettype + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.fatminm + 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 + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.fatminm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year), 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.random.drop1 <- lmer(log.Meters ~ Wettype + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.lstmaxm + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year), REML = TRUE, data = Ponded) summary(top.random.drop1) Backwards2 <- drop1(top.random.drop1, test = "Chisq", trace = TRUE, REML = FALSE) Backwards2 top.random.drop2 <- lmer(log.Meters ~ Wettype + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year), REML = TRUE, data = Ponded) summary(top.random.drop2) Backwards3 <- drop1(top.random.drop2, test = "Chisq", trace = TRUE, REML = FALSE) Backwards3 top.model <- lmer(log.Meters ~ Wettype + agtype + scaled.log.NewfootPA + scaled.espro2 + scaled.log.lsprs + scaled.log.wprs + scaled.wtmaxsf + (1|Wetland) + (1|Year), REML = TRUE, data = Ponded) summary(top.model) ##### ##### #####obtain pseduo 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$agtype) 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.fatminm) 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.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.fatmaxm) 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.log.wprs) 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") } ##### ##### #####transform surrounding ag to binomial Ponded$NoAg <- Ponded$agtype == 0 Ponded$Pivots <- Ponded$agtype == 1 Ponded$Dryland <- Ponded$agtype == 2 Ponded$Gravity <- Ponded$agtype == 3 for(i in 1:length(Ponded$NoAg)){ Ponded$NoAg[i] <- recode(Ponded$NoAg[i], "c('TRUE') = 1; ('FALSE') = 0") } for(i in 1:length(Ponded$Pivots)){ Ponded$Pivots[i] <- recode(Ponded$Pivots[i], "c('TRUE') = 1; ('FALSE') = 0") } for(i in 1:length(Ponded$Gravity)){ Ponded$Gravity[i] <- recode(Ponded$Gravity[i], "c('TRUE') = 1; ('FALSE') = 0") } for(i in 1:length(Ponded$Dryland)){ Ponded$Dryland[i] = recode(Ponded$Dryland[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.25906+ + 1.32752 * Ponded$SemiPerm[[i]]+ - 0.76114 * Ponded$Temporary[[i]]+ - 0.28251 * Ponded$Pivots[[i]]+ + 0.36406 * Ponded$Dryland[[i]]+ - 0.56592 * Ponded$Gravity[[i]]+ - ((0.40144 * Ponded$log.NewfootPA[[i]] - 0.40144 * -4.23480) / 0.53077)+ + ((0.07812 * Ponded$espro2[[i]] - 0.07812 * 0.28077) / 0.47154)+ + ((0.30366 * Ponded$log.lsprs[[i]] - 0.30366 * 5.45610) / 0.29539)+ + ((0.26013 * Ponded$log.wprs[[i]] - 0.26013 * 4.31325) / 0.76437)+ + ((0.77571 * Ponded$wtmaxsf[[i]] - 0.77571 * 26.01538) / 4.40183) } ##### ##### #####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 #get PCC Predicted <- Presence.LMM.CrossValidation$PredObsData[1] Observed <- Presence.LMM.CrossValidation$PredObsData[2] cor(Predicted$pred, Observed$obs, method = "pearson") plot(Predicted$pred, Observed$obs, pch = 19, cex = 0.35, col = "black", frame.plot = FALSE, xlim = c(4, 15), ylim = c(-10, 15), xlab = "Predicted log-transformed ponded square meters", ylab = "Observed log-transformed ponded square meters") abline(lm(Observed$obs ~ Predicted$pred), col = "red", lwd = 2) #write validation data to a .csv file ValidationData <- data.frame(Predicted = Predicted$pred) ValidationData$Observed <- Observed$obs write.csv(ValidationData, "FloodedAreaValidationAllWetlands.csv") ##### ##### #####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.fatminm, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.wtmaxsf, 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.fatminm|Ponded$Wettype, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.wtmaxsf|Ponded$Wettype, type = 'smooth') #subset by agricultural field type xyplot(resid(top.model) ~ Ponded$scaled.log.NewfootPA|Ponded$agtype, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.log.lsprs|Ponded$agtype, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.fatminm|Ponded$agtype, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.wtmaxsf|Ponded$agtype, 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.fatminm|Ponded$Year, type = 'smooth') xyplot(resid(top.model) ~ Ponded$scaled.wtmaxsf|Ponded$Year, type = 'smooth')