#### 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: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.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.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_jan18.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_jan2018.csv")