##### ##### #####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("BinomialOnlyRowcropsModestChange3.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) Ponded$agtype <- factor(Ponded$agtype) ##### ##### #####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$Pivots <- Ponded$agtype == 1 Ponded$Dryland <- Ponded$agtype == 2 Ponded$Gravity <- Ponded$agtype == 3 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 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.82601+ + 0.78335 * Ponded$SemiPerm[[i]]+ - 1.73349 * Ponded$Temporary[[i]]+ - 0.33092 * Ponded$Dryland[[i]]+ + 0.45676 * Ponded$Gravity[[i]]+ - ((1.07370 * Ponded$log.NewfootPA[[i]] - 1.07370 * -4.00697) / 0.45949)+ + ((0.45429 * Ponded$log.lsprsPostChange[[i]] - 0.45429 * 5.40933) / 0.32265)+ - ((0.75777 * Ponded$fatmaxmPostChange[[i]] - 0.75777 * 16.01264) / 1.10196)+ + ((0.56492 * Ponded$log.wprsPostChange[[i]] - 0.56492 * 4.27221) / 0.76159)+ - ((0.98802 * Ponded$wtmaxmPostChange[[i]] - 0.98802 * 6.37463) / 1.41227) } Ponded$Pred.Presence.ManualFixedOnly <- ilogit(Ponded$Pred.Presence.ManualFixedOnly) head(Ponded) ##### ##### #####assign inundation/non-inundation at probability threshold of 0.55 Ponded$Pred.Presence.ManualFixedOnlyPresence <- NA for(x in 1:nrow(Ponded)){ if(Ponded$Pred.Presence.ManualFixedOnly[x] < 0.55){ 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 PresentModestChange <- subset(Ponded, Pred.Presence.ManualFixedOnlyPresence == 1) #write to csv file write.csv(PresentModestChange, "BinomialOnlyRowcropsModestChangePresent3.csv")