#### impact of drought analysis
### Reiter et al.
## date: January 2017; mod June 2017; mod December 2017
## This code is to run generalized additive mixed models to water data in the Central Valley
## this is for single models fit to all year types comnbined and statistical asseessment of
## overall differences between year types; the effect of the precipitation
## the effect of land ownership (in wetlands).
## required data: "wetlands_merge_scene_dates_ownership.csv"; "err2.csv";
## "Ag_merge_precip3.csv"
##########################################3
## 2000-2015 combined models
## addded in data to distinguish between protected and non protected wetlands
setwd("C:/water")
library(gamm4)
library(ggplot2)
## wetland data
wetwater<-read.csv("wetlands_merge_scene_dates_ownership.csv", header=T)
## read in bias correction data for reference
bc<-read.csv("C:/water/err2.csv", header=T)
## add in drought indicator
wetwater$drought<-ifelse(wetwater$water.year%in%c(2001:2002,2007:2009,2013:2015),1,0)
wetwater$ytp<-ifelse(wetwater$year<2013&wetwater$drought==1,1,
ifelse(wetwater$year>2012,2,0))
wetwater$ytpf<-ifelse(wetwater$year<2013&wetwater$drought==1,"1d",
ifelse(wetwater$year>2012,"2d","0d"))
## indicator for recent drought
wetwater$rd<-ifelse(wetwater$year>2011,1,0)
## standardize some covariates
wetwater$zday = (wetwater$yday-150)/100
wetwater$id = c(1:nrow(wetwater))
wetwater$sprecip4 = wetwater$precip.4wk/100
wetwater$sprecip2 = wetwater$precip.2wk/100
## correct for classification bias between sensors
wetwater$nfloodbc<-ifelse(wetwater$year<2012,
round(wetwater$n.flooded-wetwater$n.flooded*-.0127),
round(wetwater$n.flooded-wetwater$n.flooded*-.1065))
## some issueswith bias correction where the corrected number of flooded pixels is greater than the number sampled
## line below makes the correctd value equal to the number sampled.
wetwater$nfloodbc<-ifelse(wetwater$nfloodbc>wetwater$n.sampled,wetwater$n.sampled,wetwater$nfloodbc)
##
####### grab only seasonal wetlands
wetwater1seas<-subset(wetwater,wetwater$landcover==1 & wetwater$prop.sampled>=0.5|wetwater$landcover==10 & wetwater$prop.sampled>=0.5)
### make public private = 1,0
wetwater1seas$lc<-ifelse(wetwater1seas$landcover==10,1,0)
## models 2000-2015
wetwater1seasf1A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ ytpf+s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf1A$mer)
wetwater1seasf2A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ ytpf*lc + sprecip2 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf2A$mer)
wetwater1seasf3A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ytpf*lc+ sprecip4 +s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf3A$mer)
wetwater1seasf4A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd+ s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf4A$mer)
wetwater1seasf5A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd*lc+ sprecip2 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf5A$mer)
plot(residuals(wetwater1seasf1A$gam)~fitted(wetwater1seasf1A$gam))
wetwater1seasf6A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd*lc+ sprecip4 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1seas, weights=wetwater1seas$prop.sampled)
summary(wetwater1seasf6A$mer)
## need to get AIC
wetwater1seasAIC<-as.vector(c(summary(wetwater1seasf1A$mer)$AIC[1],
summary(wetwater1seasf2A$mer)$AIC[1],
summary(wetwater1seasf3A$mer)$AIC[1],
summary(wetwater1seasf4A$mer)$AIC[1],
summary(wetwater1seasf5A$mer)$AIC[1],
summary(wetwater1seasf6A$mer)$AIC[1]))
wetwater1seasR2<-as.vector(c(summary(wetwater1seasf1A$gam)$r.sq[1],
summary(wetwater1seasf2A$gam)$r.sq[1],
summary(wetwater1seasf3A$gam)$r.sq[1],
summary(wetwater1seasf4A$gam)$r.sq[1],
summary(wetwater1seasf5A$gam)$r.sq[1],
summary(wetwater1seasf6A$gam)$r.sq[1]))
wetwater1seasMod<-c("1A","2A","3A",
"4A","5A","6A")
## Model fit table
wetwater1seasTAB<-data.frame(cbind(wetwater1seasMod, wetwater1seasAIC,wetwater1seasR2))
write.csv(wetwater1seasTAB,"wetwater1seasTAB_apr2018.csv")
## semi-permanent wetland
wetwater1sp<-subset(wetwater,wetwater$landcover==2 & wetwater$prop.sampled>=0.5|wetwater$landcover==20 & wetwater$prop.sampled>=0.5)
### make public private = 1,0
wetwater1sp$lc<-ifelse(wetwater1sp$landcover==20,1,0)
## models 2000-2015
wetwater1spf1A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ ytpf+s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf1A$mer)
wetwater1spf2A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ ytpf*lc + sprecip2 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf2A$mer)
wetwater1spf3A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~ytpf*lc+ sprecip4 +s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf3A$mer)
wetwater1spf4A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd+ s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf4A$mer)
wetwater1spf5A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd*lc+ sprecip2 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf5A$mer)
plot(residuals(wetwater1spf1A$gam)~fitted(wetwater1spf1A$gam))
wetwater1spf6A<-gamm4(cbind(nfloodbc, n.sampled-nfloodbc)~rd*lc+ sprecip4 + s(zday, k=10),
random=~(1|water.year) + (1|id), family=binomial, data=wetwater1sp, weights=wetwater1sp$prop.sampled)
summary(wetwater1spf6A$mer)
## need to get AIC
wetwater1spAIC<-as.vector(c(summary(wetwater1spf1A$mer)$AIC[1],
summary(wetwater1spf2A$mer)$AIC[1],
summary(wetwater1spf3A$mer)$AIC[1],
summary(wetwater1spf4A$mer)$AIC[1],
summary(wetwater1spf5A$mer)$AIC[1],
summary(wetwater1spf6A$mer)$AIC[1]))
wetwater1spR2<-as.vector(c(summary(wetwater1spf1A$gam)$r.sq[1],
summary(wetwater1spf2A$gam)$r.sq[1],
summary(wetwater1spf3A$gam)$r.sq[1],
summary(wetwater1spf4A$gam)$r.sq[1],
summary(wetwater1spf5A$gam)$r.sq[1],
summary(wetwater1spf6A$gam)$r.sq[1]))
wetwater1spMod<-c("1A","2A","3A",
"4A","5A","6A")
## Model fit table
wetwater1spTAB<-data.frame(cbind(wetwater1spMod, wetwater1spAIC,wetwater1spR2))
write.csv(wetwater1spTAB,"wetwater1spTAB_apr2018.csv")
################################
## extract coefficients from models
seascoefs<-data.frame(rbind(summary(wetwater1seasf1A$mer)$coefficients,
summary(wetwater1seasf2A$mer)$coefficients,
summary(wetwater1seasf3A$mer)$coefficients,
summary(wetwater1seasf4A$mer)$coefficients,
summary(wetwater1seasf5A$mer)$coefficients,
summary(wetwater1seasf6A$mer)$coefficients))
seasmod<-c(rep("1a",4),rep("2a",8), rep("3a",8),
rep("4a",3),rep("5a",6), rep("6a",6))
cov<-c("int","d","rd","day",
"int","d","rd","lc","p2","lc*d","lc*rd","day",
"int","d","rd","lc","p4","lc*d","lc*rd","day",
"int","sd","day",
"int","sd","lc*sd","p2","lc","day",
"int","sd","lc*sd","p4","lc","day")
seascoefs2<-data.frame(cbind(seascoefs,seasmod,cov))
names(seascoefs2)<-c("est","se","zval","pval","model","covar")
## coefficient summary table
write.csv(seascoefs2,"seascoefs2_apr18.csv")
## semi-permanent wetlands
spcoefs<-data.frame(rbind(summary(wetwater1spf1A$mer)$coefficients,
summary(wetwater1spf2A$mer)$coefficients,
summary(wetwater1spf3A$mer)$coefficients,
summary(wetwater1spf4A$mer)$coefficients,
summary(wetwater1spf5A$mer)$coefficients,
summary(wetwater1spf6A$mer)$coefficients))
spmod<-c(rep("1a",4),rep("2a",8), rep("3a",8),
rep("4a",3),rep("5a",6), rep("6a",6))
cov<-c("int","d","rd","day",
"int","d","rd","lc","p2","lc*d","lc*rd","day",
"int","d","rd","lc","p4","lc*d","lc*rd","day",
"int","sd","day",
"int","sd","lc*sd","p2","lc","day",
"int","sd","lc*sd","p4","lc","day")
spcoefs2<-data.frame(cbind(spcoefs,spmod,cov))
names(spcoefs2)<-c("est","se","zval","pval","model","covar")
## coefficient summary table
write.csv(spcoefs2,"spcoefs2_apr2018.csv")