# code for bird response to tree removal
# sarah thompson, proofed for JAE sept 28, 2015
# analysis conducted using R version 3.1.2

# change to file location
setwd("C:/Users/sjthompson/Desktop/JAEonlinedataArchive") # The included file path was used by author, other users should replace it with their own working directory where the input csv file is located. 

# Import data
longnorm<-read.csv("datforarchive.csv", header=T)
head(longnorm)

# Load necessary packages
library(gdata) # for drop.levels()
library(plotrix) # for plotCI()
library(nlme)
library(lme4)
library(plyr)

##### definition of data fields #####
# stn.unit.year --- unique id (station, wpa/unit, year of survey)
# trt --- site was originally cut or control 
# wpa --- site (waterfowl production area)
# station  -- point count station number
# year ---- year & yearfac -- year as a factor
# spec --- species AOU code
# total100 -- total count of that species within 100 m
# ave_shrub --- average distance to shrub (average of 4 directions) (m)
# ave_slash --- average distance to slash pile (m)
# ave_tree --- average distance to tree (m)
# ave_rob_mean --- average of Robel readings (grass height and density, dm)
# ave_lit_mean --- average litter depth reading (cm)

# data set has a row for each SPECIES observed at the point, thus information like 
# robel and distance to shrub is repeated 

################################################################
#### data preparation for analysis #############################

# Subset into 2 parts, bird data and habitat data
# isolate columns with information on birds
speccols <- longnorm[, c("stn.unit.year", "spec", "total100")]

# separate from information on habitat
dat1 <- longnorm[, c("stn.unit.year","wpa","trt", "year", "observer",
                     "ave_rob_mean", "ave_lit_mean", "ave_shrub",
                     "ave_tree", "ave_slash" )]

# Reduce habitat information to 1 row for each stn.unit.yr/unique point count
covsonly1 <- dat1[!duplicated(dat1[]), ]

# reduce to only identifying columns
sites1 <- covsonly1[,c("stn.unit.year", "wpa",  "year", "trt", "observer" )]

# this shows the balanced design and 5 places w 19 counts 
table(sites1$wpa, sites1$year)

##### Function to create datasets for individual species #######
CreateDataset <- function(spec.i) {
  # creates data set for individual species, sums over points, fixes NAs
  #
  # Args:
  # only requires the name of the species
  # parses out single species, adds up all observations, adds formatting for 
  # models (treatment indicator etc)
  g.i <- drop.levels(subset(speccols, spec==spec.i, drop=TRUE))
  g.sum <- ddply(g.i, "stn.unit.year", summarise, total=sum(total100))
  dat <- merge(g.sum, sites1, by="stn.unit.year", all.x=TRUE, all.y=TRUE)
  dat$total[is.na(dat$total)] <- 0
  dat$trtnum <- as.numeric(dat$trt) - 1
  dat$year <- dat$year-2005
  dat$I.trt <- "0"
  dat$I.trt[dat$trt == "cut"] <- "1"
  dat$I.trt[dat$year == 0 & dat$trt == "cut"] <- "0"
  dat$I.trt <- as.factor(dat$I.trt)
  x <- order(dat$wpa)
  dat <- dat[x, ]
  nameout <- as.data.frame(dat)
  return(nameout)
}

# Data sub-dataset for each of 16 key species 
g.amgo <- CreateDataset("AMGO")
g.bhco <- CreateDataset("BHCO")
g.bobo <- CreateDataset("BOBO")
g.ccsp <- CreateDataset("CCSP")
g.cogr <- CreateDataset("COGR")
g.coye <- CreateDataset("COYE")
g.eaki <- CreateDataset("EAKI")
g.mawr <- CreateDataset("MAWR")
g.sosp <- CreateDataset("SOSP")
g.tres <- CreateDataset("TRES")
g.swsp <- CreateDataset("SWSP")
g.sewr <- CreateDataset("SEWR")
g.rwbl <- CreateDataset("RWBL")
g.yhbl <- CreateDataset("YHBL")
g.ywar <- CreateDataset("YWAR")
g.wifl <- CreateDataset("WIFL")


###### Create data sets for grouped species ##########
# create subsets
# woodland birds that nest in cavities
g.cavity<-drop.levels(subset(speccols, spec=="HOWR"|spec=="BCCH"|spec=="NOFL"|
                               spec=="GCFL"| spec=="WBNU"|spec=="DOWO"|spec=="HAWO"|
                               spec=="RBWO"| spec=="PIWO"|spec=="YBSA", drop=TRUE))    

# all woodland/parkland birds
g.woodall<-drop.levels(subset(speccols, spec=="OROR"|spec=="WAVI"|spec=="LEFL"|
                                spec=="CEDW"|spec=="BAOR"|spec=="RBGR"|
                                spec=="REVI"|spec=="BBCU"|spec=="BLJA"|
                                spec=="EAWP"|spec=="AMCR"|spec=="RTHU"|
                                spec=="INBU"|spec=="YTVI"|spec=="AMRE"|
                                spec=="OSFL"|spec=="CHSP"|spec=="EAPH"|
                                spec=="YBCU"|spec=="SCTA"|
                                spec=="NOCA"|spec=="GRCA"|spec=="MODO"|
                                spec=="AMRO"|spec=="HOFI"|spec=="WEKI"|
                                spec=="CHSP"|spec=="BRTH"|spec=="BRBL", drop=TRUE))


# combines all wetland (ground and floating nest) and waterfowl
g.wetall<-drop.levels(subset(speccols,  
                             spec=="CAGO"|spec=="RUDU"|spec=="RNDU"|
                               spec=="VIRA"|spec=="AMBI"|spec=="WISN"|
                               spec=="LEBI"|spec=="MALL"|spec=="BWTE"|
                               spec=="GADW"|spec=="NSHO"|spec=="BLTE"|
                               spec=="PBGR"|spec=="SORA"|spec=="AMCO"|
                               spec=="REDH"|spec=="CANV"|spec=="RNGR"|
                               spec=="EAGR"|spec=="FOTE"|spec=="CAGO"|
                               spec=="MALL"|spec=="BWTE"|spec=="GADW"|
                               spec=="NSHO"| spec=="RUDU"|spec=="CANV"|
                               spec=="RNDU"|spec=="HOME"|spec=="WODU"|
                               spec=="REDH", drop=TRUE))

# grassland birds, omitting common species
g.grass<-drop.levels(subset(speccols, spec=="SAVS"|spec=="UPSA"|
                              spec=="LCSP"|spec=="GRSP"|spec=="HOLA"| 
                              spec=="VESP"|spec=="WEME", drop=TRUE))

# all grass combined
g.grassall <-drop.levels(subset(speccols, spec=="SAVS"|spec=="UPSA"|
                                  spec=="LCSP"|spec=="GRSP"|spec=="HOLA"| 
                                  spec=="VESP"|spec=="WEME"
                                |spec=="BOBO"|spec=="CCSP"|spec=="SEWR"|
                                  spec=="NOHA", drop=TRUE))



##### Same as CreateDataset1, but for grouped species #####
CreateDataset2 <- function(group.i){
  # creates data set for grouped species, sums over points, fixes NAs
  #
  # Args:
  # only requires the name of the group, created above
  # 
  # Returns: 
  # the dataset
  g.sum <- ddply(group.i, "stn.unit.year", summarise, total = sum(total100))
  dat   <- merge(g.sum, sites1, by="stn.unit.year", all.x=TRUE, all.y=TRUE)
  dat$total[is.na(dat$total)] <- 0
  dat$trtnum <- as.numeric(dat$trt) - 1
  dat$year <- dat$year-2005
  dat$I.trt <- "0"
  dat$I.trt[dat$trt == "cut"] <- "1"
  dat$I.trt[dat$year == 0 & dat$trt == "cut"] <- "0"
  x <- order(dat$wpa)
  dat <- dat[x, ]
  nameout <- as.data.frame(dat)
  return(nameout)
}


# use CreateDataset2 to generate datasets
d.cavity   <- CreateDataset2(g.cavity)
d.woodall <- CreateDataset2(g.woodall)
d.grass    <- CreateDataset2(g.grass)
d.wetall  <- CreateDataset2(g.wetall)
d.grassall <- CreateDataset2(g.grassall)

##################################################################
###### Bird analysis with GEEs  ##################################
# required packages
library(geepack) # for geeglm()
library(MuMIn) # for QIC

# prepare lists that can be looped through analysis
# datasets
dat.list<-list(g.bhco, g.bobo, g.ccsp, g.coye, g.eaki, g.mawr, 
               g.sewr, g.swsp, g.ywar, d.cavity, d.woodall, 
               d.wetall, d.grassall, d.grass)
# names
spec.list<-c("g.bhco","g.bobo","g.ccsp", "g.coye","g.eaki","g.mawr",
             "g.sewr","g.swsp","g.ywar", "d.cavity", "d.woodall", 
             "d.wetall", "d.grassall", "d.grass")

# empty dataframe for outputting AIC values
dfAIC = data.frame(matrix(vector(),  length(dat.list), 7, 
                          dimnames=list(c(),  c("spec" ,"sum",  "null", "year", "trt", "yrptrt", "yrXtrt"))), 
                   stringsAsFactors=F)

# loop all through the analysis and output table
for (i in 1:13){
  dat<-dat.list[[i]]
  null <-geeglm(total ~ 1, data=dat, id=wpa, 
                family=poisson("log"), corstr="exchangeable", 
                na.action="na.omit")
  year <-geeglm(total ~ as.factor(year), data=dat, id=wpa, 
                family=poisson("log"), corstr="exchangeable", 
                na.action="na.omit")
  trt <-geeglm(total ~ as.factor(I.trt), data=dat, id=wpa, 
               family=poisson("log"), corstr="exchangeable", 
               na.action="na.omit")
  yrptrt <-geeglm(total ~ as.factor(year)+as.factor(I.trt), data=dat, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", 
                  na.action="na.omit")
  yrXtrt <-geeglm(total ~ as.factor(year)*trtnum, data=dat, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", 
                  na.action="na.omit")
  out<-QIC(null, year, trt, yrptrt, yrXtrt)
  minQ<-min(out[, 1])
  dQIC<-round(out[, 1]-minQ, 2)
  dfAIC[i, ]<-c( spec.list[i], sum(dat$total), dQIC)
}

### NOTE !!! grass rare wont work with yrXtrt, so must manually run this for #12
# and remove yrXtrt and put in an NA
i<-14
dat<-dat.list[[i]]
null <-geeglm(total ~ 1, data=dat, id=wpa, 
              family=poisson("log"), corstr="exchangeable", 
              na.action="na.omit")
year <-geeglm(total ~ as.factor(year), data=dat, id=wpa, 
              family=poisson("log"), corstr="exchangeable", 
              na.action="na.omit")
trt <-geeglm(total ~ as.factor(I.trt), data=dat, id=wpa, 
             family=poisson("log"), corstr="exchangeable", 
             na.action="na.omit")
yrptrt <-geeglm(total ~ as.factor(year)+as.factor(I.trt), data=dat, 
                id=wpa, family=poisson("log"), corstr="exchangeable", 
                na.action="na.omit")
#yrXtrt <-geeglm(total ~ as.factor(year)*trtnum, data=dat, 
#               id=wpa, family=poisson("log"), corstr="exchangeable", 
#              na.action="na.omit")
out<-QIC(null, year, trt, yrptrt)
minQ<-min(out[, 1])
dQIC<-round(out[, 1]-minQ, 2)
dfAIC[i, ]<-c( spec.list[i], sum(dat$total), dQIC, "NA")

# AIC table for birds only
dfAIC


############################################################################
############## plotting avian response #######################################

# yearXtrt mods for each species/group
gee.bhco <-geeglm(total ~ as.factor(year)*trtnum, data=g.bhco, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.bobo <-geeglm(total ~ as.factor(year)*trtnum, data=g.bobo, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.coye <-geeglm(total ~ as.factor(year)*trtnum, data=g.coye, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.eaki <-geeglm(total ~ as.factor(year)*trtnum, data=g.eaki, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.mawr<-geeglm(total ~ as.factor(year)*trtnum, data=g.mawr, 
                 id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.sewr<-geeglm(total ~ as.factor(year)*trtnum, data=g.sewr, 
                 id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.swsp<-geeglm(total ~ as.factor(year)*trtnum, data=g.swsp, 
                 id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.cavity <-geeglm(total ~ as.factor(year)*trtnum, data=d.cavity, 
                    id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.wetall <-geeglm(total ~ as.factor(year)*trtnum, data=d.wetall, 
                    id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.grassall <-geeglm(total ~ as.factor(year)*trtnum, data=d.grassall, 
                      id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.ywar <-geeglm(total ~ as.factor(year)*trtnum, data=g.ywar, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.ccsp <-geeglm(total ~ as.factor(year)*trtnum, data=g.ccsp, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")

gee.woodall <-geeglm(total ~ as.factor(year)*trtnum, data=d.woodall, 
                     id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.grass <-geeglm(total ~ as.factor(year)*trtnum, data=d.grass, 
                   id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
gee.grass2 <-geeglm(total ~ as.factor(year)+I.trt, data=d.grass, 
                    id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")

# predict, transform, plot output of GEE 
# function that plots year*trt output for GEE version of analysis
GeePlot<-function(mods.out, name.i, xax=FALSE){
  fit1<-mods.out
  ########### predict and transform #####################
  newdat <- expand.grid(year=(0:6),trtnum=(0:1),total = 0)
  mm.geeEX<-model.matrix(terms(fit1), newdat)
  newdat$pred.T = mm.geeEX %*% coef(fit1)
  tvar1.gee <- diag(mm.geeEX %*% tcrossprod(fit1$geese$vbeta, mm.geeEX))
  
  newdat$pred.lcl<- exp(newdat$pred.T-1.96*sqrt(tvar1.gee))
  newdat$pred.ucl <- exp(newdat$pred.T+1.96*sqrt(tvar1.gee))
  newdat$pred.T<-exp(newdat$pred.T)
  
  ######## Set up plot #######################
  y<-c(0,max(newdat$pred.ucl)+max(newdat$pred.ucl)*0.2)
  # control sites, grey open circles, solid line
  with (data = newdat, plot(year[trtnum==0]-0.01, pred.T[trtnum==0], 
                            col=gr, ylab="", xlab="", 
                            ylim=y, xaxt="n", xlim=c(0-.3, 6+.3),
                            cex=1, lty=1,lwd = 1.8,
                            cex.main=1, cex.axis=cx, type="l"))
  
  plotCI(newdat$year[newdat$trtnum==0]-.01, newdat$pred.T[newdat$trtnum==0],
         ui=newdat$pred.ucl[newdat$trtnum==0],
         li=newdat$pred.lcl[newdat$trtnum==0], 
         col=gr, add=TRUE, xaxt="n", pch=NA)
  with (data = newdat, points(year[trtnum==0]-.01, pred.T[trtnum==0], 
                              col=gr, bg="white",
                              pch=c(21), cex=1.4, lwd=2,
                              cex.main=.8, cex.axis=cx))
  
  ## option to add axis or just axis tick marks ################
  if (xax==TRUE) {
    axis(1, 0:6, c("0","1","2","3","4","5","6"), cex.axis=cx)
  }
  else {
    axis(1, 0:6, tck=-0.03, labels=FALSE)
  }
  title(main=name.i, font=4, cex.main=1.1)
  ##bird.tot<-sum(newdat$total)
  ##legend("topleft", paste("  n =", bird.tot), bty="n", cex=0.9)
  
  # cut sites, grey, squares
  with (data = newdat, points(year[trtnum==1]+0.1, pred.T[trtnum==1], 
                              col="black", ylab="", xlab="", cex=1.6, 
                              pch=18))
  plotCI(newdat$year[newdat$trtnum==1]+0.1, 
         newdat$pred.T[newdat$trtnum==1],
         ui=newdat$pred.ucl[newdat$trtnum==1],
         li=newdat$pred.lcl[newdat$trtnum==1],
         add=TRUE, scol="black", pch=NA)
  with (data = newdat, lines(year[trtnum==1]+0.1, 
                             pred.T[trtnum==1], 
                             col="black",lwd=1.2, lty=2))
}


###### GEE MULTI PANEL PLOT, IND SPECIES ##############
# prep for multi-panel plot of individual species
# Data in list
d.group<-list(g.bhco, g.bobo, g.ccsp, g.eaki, g.mawr,
              g.sewr, g.swsp, g.ywar)
# plot names
name.i<-c("A. Brown-headed cowbird",  "B. Bobolink", 
          "C. Clay-colored sparrow",    "D. Common yellowthroat",
          "E. Eastern kingbird",
          "F. Marsh wren" ,       "G. Sedge wren",
          "H. Swamp sparrow",   "I. Yellow warbler")

# model outputs in a list, all are year X trt
mods.out<-list(gee.bhco, gee.bobo, gee.ccsp, gee.coye, gee.eaki, 
               gee.mawr, gee.sewr, gee.swsp, gee.ywar)

# remove # to output high-res plot
# plot.new()
# tiff("Plotspecgee2.tiff", width = 6, height = 7, units = 'in', res = 300) ## turn on for tiff output

# begin formal plot, 3 rows by 2 cols
par(mar=c(1.5, 2.0, 1.5, 0.75))
par(mgp=c(2, 0.5, 0))
par(oma=c(3, 2, 0, 0))
par(mfrow=c(3,3))
cx<-0.9
gr<-"grey44"
cy<-1.4

#' BHCO
# summary(gee.bhco)
GeePlot(mods.out[[1]], name.i[[1]])
text(4,0.7,expression("--"), cex=cy)
text(5,0.7,expression("--"), cex=cy)
text(6,0.7,expression("--"), cex=cy)

#'  BOBO
# summary(gee.bobo)
GeePlot(mods.out[[2]], name.i[[2]])
text(2,1.1,expression("--"), cex=cy)

#' CCSP
# summary(gee.ccsp)
GeePlot(mods.out[[3]], name.i[[3]])
text(6.1,1.6,expression("+"), cex=cy)# no longer sig with dropped point

#' COYE
# summary(gee.coye)
GeePlot(mods.out[[4]], name.i[[4]])
text(1.1,1.35,expression("+"), cex=cy)

#' EAKI
# summary(gee.eaki)
GeePlot(mods.out[[5]], name.i[[5]])
text(3,0.44,expression("--"), cex=cy)
text(6.1,0.40,expression("+"), cex=cy) # no longer Sig.

#' MAWR
GeePlot(mods.out[[6]], name.i[[6]], xax=FALSE)
# summary(gee.mawr)# no sigs

#' SEWR
GeePlot(mods.out[[7]], name.i[[7]], xax=TRUE)
# summary(gee.sewr)
text(2,.9,expression("--"), cex=cy)
text(5.1,1.6,expression("+"), cex=cy)

#' SWSP
GeePlot(mods.out[[8]], name.i[[8]], xax=TRUE)
# summary(gee.swsp)# no sigs

#' YWAR
# summary(gee.ywar)
GeePlot(mods.out[[9]], name.i[[9]], xax=TRUE)
text(2.1, 1.2,expression("-"), cex=cy)
text(3, 1.2,expression("--"), cex=cy)
text(4, 1.35,expression("-"), cex=cy)
text(5, 1.35,expression("--"), cex=cy)
text(6, 1.35,expression("--"), cex=cy)

mtext("Years since start of study", side = 1, line = 0, outer = TRUE, padj=.5)
mtext("Birds per point (estimated)", side = 2, line = 0, outer = TRUE)
# dev.off()


###### MULTI PANEL PLOT GROUPS  GEE #########
name.i<-c("A. Cavity Nesting",
          "B. Tree/Shrub Nesting", 
          "C. All Wetland", 
          "D. All grassland combined", 
          "E. Rare grassland")

mods.out<-list(gee.cavity, gee.woodall,  gee.wetall, 
               gee.grassall, gee.grass)

# Multi- Panel Plot, GROUPS, GEE ######

# turn on if you want a high-res export
#plot.new()
#tiff("Plot_groups_gee2.tiff", width = 6, height = 5, units = 'in', res = 300)                
par(mfrow=c(3,2))
par(mar=c(1.5, 2.0, 1.5, 0.75))
par(mgp=c(2, 0.5, 0))
par(oma=c(3, 2, 0, 0))

layout(matrix(c(1,1,2,2,3,3,4,4,0,5,5,0), byrow=T, nrow=3))

cplus<-1.2 # cex of plus signs, I think

GeePlot(mods.out[[1]], name.i[[1]])
# summary(gee.cavity)
text(1, 0.4,expression("--"), cex=cplus)
text(2, 0.5,expression("--"), cex=cplus)
text(3, 0.5,expression("--"), cex=cplus)
text(4, 0.50,expression("--"), cex=cplus)
text(5, 0.80,expression("--"), cex=cplus)
text(6, 0.70,expression("--"), cex=cplus)
#legend("topright", "trt X year", bty="n")


GeePlot(mods.out[[2]], name.i[[2]])
# summary(gee.woodall)
text(3,1.8,expression("--"), cex=cplus)
text(4,1.8,expression("-"), cex=cplus)
text(5,2.2,expression("--"), cex=cplus)
text(6,1.9,expression("--"), cex=cplus)

GeePlot(mods.out[[3]], name.i[[3]], xax=TRUE)
# summary(gee.wetall)
text(2.1,1.1,expression("+"), cex=cplus)
text(4,  0.9,expression("+"), cex=cplus)
text(5.1,1.2,expression("++"), cex=cplus)
text(6.1,1.72,expression("++"), cex=cplus)

summary(gee.grassall)
GeePlot(mods.out[[4]], name.i[[4]], xax=TRUE)
text(6.1,2.5,expression("+"), cex=cplus)

# rare ground-grass  
fit1<-gee.grass

# specifically plot ground-grass rare, fixing errant estimate
newdat <- expand.grid(year=(0:6),trtnum=(0:1),total = 0)
mm.geeEX<-model.matrix(terms(fit1), newdat)
newdat$pred.T = mm.geeEX %*% coef(fit1)
tvar1.gee <- diag(mm.geeEX %*% tcrossprod(fit1$geese$vbeta, mm.geeEX))
newdat$pred.lcl<- exp(newdat$pred.T-1.96*sqrt(tvar1.gee))
newdat$pred.ucl <- exp(newdat$pred.T+1.96*sqrt(tvar1.gee))
newdat$pred.T<-exp(newdat$pred.T)
options(scipen=100)
newdat$pred.ucl # you can see one problematic value
newdat$pred.ucl[newdat$pred.ucl>1]<-0  # set to something reasonable for plot


y<-c(0, 0.4)
# control sites, grey open circles, solid line
with (data = newdat, plot(year[trtnum==0]-0.01, pred.T[trtnum==0], 
                          col=gr, ylab="", xlab="", 
                          ylim=y, xaxt="n", xlim=c(0-.3, 6+.3),
                          cex=1, lty=1,lwd = 1.8,cex.main=1, type="l"))

plotCI(newdat$year[newdat$trtnum==0]-.01, newdat$pred.T[newdat$trtnum==0],
       ui=newdat$pred.ucl[newdat$trtnum==0],
       li=newdat$pred.lcl[newdat$trtnum==0], 
       col=gr, add=TRUE, xaxt="n", pch=NA)
with (data = newdat, points(year[trtnum==0]-.01, pred.T[trtnum==0], 
                            col=gr, bg="white",pch=c(21), cex=1.4, lwd=2,
                            cex.main=.8, cex.axis=cx))

axis(1, 0:6, c("0","1","2","3","4","5","6"), cex.axis=cx)
title(main=name.i[[5]], font=4, cex.main=1.1)

# cut sites, grey, squares
with (data = newdat, points(year[trtnum==1]+0.1, pred.T[trtnum==1], 
                            col="black", ylab="", xlab="", cex=1.6, pch=18))
plotCI(newdat$year[newdat$trtnum==1]+0.1, 
       newdat$pred.T[newdat$trtnum==1],
       ui=newdat$pred.ucl[newdat$trtnum==1],
       li=newdat$pred.lcl[newdat$trtnum==1],
       add=TRUE, scol="black", pch=NA)
with (data = newdat, lines(year[trtnum==1]+0.1, pred.T[trtnum==1], col="black",
                           lwd=1.2, lty=2))
# end special plot for rare grassland birds
text(2,.25,expression("+"), cex=cplus)

mtext("Years since start of study", side = 1, line = 0, outer = TRUE, padj=.5)
mtext("Birds per point (estimated)", side = 2, line = 0, outer = TRUE)

#dev.off()

##############################################################
###### VEG MODELS with GLMMs #################################
covsonly1$year2 <- covsonly1$year
covsonly1$year <- covsonly1$year-2005 # don't run 2X!
covsonly1$I.trt <- "0"
covsonly1$I.trt[covsonly1$trt=="cut"] <- "1"
covsonly1$I.trt[covsonly1$year == 0 & covsonly1$trt=="cut"] <- "0"
covsonly1$I.trt <- as.factor(covsonly1$I.trt)
covsonly1$trtnum <- as.numeric(covsonly1$trt)-1
summary(covsonly1$year)


# runs through 5 veg covariates,  adds columns to a df 
# robel
null       <-lmer(covsonly1$ave_rob_mean ~ 1 + (1|wpa), data=covsonly1)
yearonly   <- lmer(covsonly1$ave_rob_mean ~ as.factor(year) +  (1|wpa), data=covsonly1)
treatment  <- lmer(covsonly1$ave_rob_mean ~ as.factor(I.trt) + (1|wpa), data=covsonly1)
yr_plus_trt<- lmer(covsonly1$ave_rob_mean ~ as.factor(year) + as.factor(I.trt) + (1|wpa),  data = covsonly1)
yr_x_trt  <- lmer(covsonly1$ave_rob_mean ~ as.factor(year)* trtnum + (1|wpa), data=covsonly1)

aic <- AIC(null, yearonly, treatment, yr_plus_trt, yr_x_trt)
rob <- aic[,2]-min(aic[,2])
aictab <- cbind(aic,rob)

# litter
null       <-lmer(covsonly1$ave_lit_mean ~ 1 +               (1|wpa) ,  data=covsonly1)
yearonly   <-lmer(covsonly1$ave_lit_mean ~ as.factor(year)+  (1|wpa) ,  data=covsonly1)
treatment  <-lmer(covsonly1$ave_lit_mean ~ as.factor(I.trt) + (1|wpa), data=covsonly1)
yr_plus_trt<-lmer(covsonly1$ave_lit_mean ~ as.factor(year)+ as.factor(I.trt) + (1|wpa),  data = covsonly1)
yr_x_trt   <-lmer(covsonly1$ave_lit_mean ~ as.factor(year)* as.factor(trt) + (1|wpa)  ,  data = covsonly1)

aic<-AIC(null, yearonly, treatment, yr_plus_trt, yr_x_trt)
lit<-aic[,2]-min(aic[,2])
aictab<-cbind(aictab,lit)

# shrub dist
null       <-lmer(covsonly1$ave_shrub ~ 1 +               (1|wpa) ,  data=covsonly1)
yearonly   <-lmer(covsonly1$ave_shrub ~ as.factor(year)+  (1|wpa) ,  data=covsonly1)
treatment  <-lmer(covsonly1$ave_shrub ~ as.factor(I.trt) + (1|wpa), data=covsonly1)
yr_plus_trt<-lmer(covsonly1$ave_shrub ~ as.factor(year)+ as.factor(I.trt) + (1|wpa),  data = covsonly1)
yr_x_trt   <-lmer(covsonly1$ave_shrub ~ as.factor(year)* as.factor(trt) + (1|wpa)  ,  data = covsonly1)

aic<-AIC(null, yearonly, treatment, yr_plus_trt, yr_x_trt)
shrubD<-aic[,2]-min(aic[,2])
aictab<-cbind(aictab,shrubD)

# tree distance
null      <-lmer(covsonly1$ave_tree ~ 1 +               (1|wpa) ,  data=covsonly1)
yearonly  <-lmer(covsonly1$ave_tree ~ as.factor(year)+  (1|wpa) ,  data=covsonly1)
treatment <-lmer(covsonly1$ave_tree ~ as.factor(I.trt) + (1|wpa), data=covsonly1)
yr_plus_trt<-lmer(covsonly1$ave_tree ~ as.factor(year)+ as.factor(I.trt) + (1|wpa),  data = covsonly1)
yr_x_trt  <-lmer(covsonly1$ave_tree ~ as.factor(year)* as.factor(trt) + (1|wpa)  ,  data = covsonly1)

aic<-AIC(null, yearonly, treatment, yr_plus_trt, yr_x_trt)
treeD<-aic[,2]-min(aic[,2])
aictab<-cbind(aictab,treeD)

# slash Distance
null       <-lmer(covsonly1$ave_slash ~ 1 +               (1|wpa) ,  data=covsonly1)
yearonly   <-lmer(covsonly1$ave_slash ~ as.factor(year)+  (1|wpa) ,  data=covsonly1)
treatment  <-lmer(covsonly1$ave_slash ~ as.factor(I.trt) + (1|wpa), data=covsonly1)
yr_plus_trt<-lmer(covsonly1$ave_slash ~ as.factor(year)+ as.factor(I.trt) + (1|wpa),  data = covsonly1)
yr_x_trt  <-lmer(covsonly1$ave_slash ~ as.factor(year)* as.factor(trt) + (1|wpa)  ,  data = covsonly1)

aic<-AIC(null, yearonly, treatment, yr_plus_trt, yr_x_trt)
slashD<-aic[,2]-min(aic[,2])

# AIC table, vegetation vars
aictabveg<-cbind(aictab,slashD)
vegtab<-t(aictabveg)
colnames(vegtab)<-c("null", "yearonly", "trtonly", "yr+trt", "yrxtrt")
vegtab

# top veg models##
# lmer and lme ~exact same outputs but lmer provides p-values for assessing
# significance of fixed effects
m.tree <- lme(ave_tree ~ as.factor(year)*trtnum, random=~1|wpa, data=covsonly1)
m.tree4 <- lmer(ave_tree ~ as.factor(year)*trtnum +(1|wpa), data=covsonly1)

m.shrub4 <- lmer(ave_shrub ~ as.factor(year)*trtnum + (1|wpa), data=covsonly1)
m.shrub <- lme(ave_shrub ~ as.factor(year)*trtnum, random=~1|wpa, data=covsonly1)

m.slash4 <- lmer(ave_slash ~ as.factor(year)*trtnum + (1|wpa), data=covsonly1)
m.slash <- lme(ave_slash ~ as.factor(year)*trtnum , random=~1|wpa, data=covsonly1)

m.lit4   <- lmer(ave_lit_mean ~ as.factor(year)*trtnum + (1|wpa), data=covsonly1)
m.lit <- lme(ave_lit_mean ~ as.factor(year)*trtnum , random=~1|wpa, data=covsonly1)

m.rob4   <- lmer(ave_rob_mean ~ as.factor(year)*trtnum + (1|wpa), data=covsonly1)
m.rob <- lme(ave_rob_mean ~ as.factor(year)*trtnum, random=~1|wpa, data=covsonly1)


# the plotting function is different bc of lmer vs glmer## prediction intervals 
VegPlot<-function(name.i, mods.out, xax=FALSE){
  fit1 <- mods.out
  g.all <- covsonly1
  # Model-based estimates from year model
  yrmod.x <- model.matrix(fit1)
  g.all$pred.yr <- yrmod.x%*%fixef(fit1) 
  g.all$pred.se.yr <- sqrt(diag(yrmod.x %*% tcrossprod(vcov(fit1),yrmod.x)))
  # back transform
  g.all$pred.T  <- g.all$pred.yr
  g.all$pred.lcl<- g.all$pred.yr-1.96*g.all$pred.se.yr
  g.all$pred.ucl<- g.all$pred.yr+1.96*g.all$pred.se.yr
  #Set up plot
  # order data, otherwise get spaghetti plot
  x<-order(g.all$trtnum, g.all$year)
  g.all <- g.all[x,]
  # this needs to be here, bc it needs to change with each iteration
  y<-c(0,max(g.all$pred.ucl)+max(g.all$pred.ucl)*0.2)
  # control sites, grey open circles, solid line
  with (data=g.all, plot(year[trtnum==0]-0.01, pred.T[trtnum==0], 
                         col=gr, ylab="", xlab="", 
                         ylim=y, xaxt="n", xlim=c(0-.3, 6+.3),
                         cex=1, lty=1,lwd = 1.8,
                         cex.main=1, cex.axis=cx, type="l"))
  plotCI(g.all$year[g.all$trtnum==0]-.01, g.all$pred.T[g.all$trtnum==0],
         ui=g.all$pred.ucl[g.all$trtnum==0],
         li=g.all$pred.lcl[g.all$trtnum==0], 
         col=gr, add=TRUE, xaxt="n", pch=NA)
  
  with (data=g.all, points(year[trtnum==0]-.01, pred.T[trtnum==0], 
                           col=gr, bg="white",
                           pch=c(21), cex=1.4, lwd=2,
                           cex.main=.8, cex.axis=cx))
  
  if (xax==TRUE) {
    axis(1, 0:6, c("0","1","2","3","4","5","6"), cex.axis=cx)
  }
  else {
    axis(1, 0:6, tck=-0.03, labels=FALSE)
  }
  
  title(main=name.i, font=4, cex.main=1.1)
  
  
  # cut sites, squares, black
  
  with (data=g.all, points(year[trtnum==1]+0.1, pred.T[trtnum==1], 
                           col="black", ylab="", xlab="", cex=1.6, 
                           pch=18))
  
  plotCI(g.all$year[g.all$trtnum==1]+0.1, 
         g.all$pred.T[g.all$trtnum==1],
         ui=g.all$pred.ucl[g.all$trtnum==1],
         li=g.all$pred.lcl[g.all$trtnum==1],
         add=TRUE, scol="black", pch=NA)
  with (data=g.all, lines(year[trtnum==1]+0.1, 
                          pred.T[trtnum==1], 
                          col="black",lwd=1.2, lty=2))
  
}

# make list of names and model outputs for figure panels
name.i<-list( "A. Dist. to tree (m)",  "B. Dist. to shrub (m)" ,
              "C. Dist. to slash (m)", "D. Litter depth (cm)" ,
              "E. Grass height (dm)")
mods.out<-list(m.tree4, m.shrub4, m.slash4, m.lit4, m.rob4)


#########################################################
##### MULTI-PANEL VEG PLOT ##############################
# dev.off() ## turn these on to output a high-quality tiff 
# plot.new()
# tiff("Plotvegmods.tiff", width = 5, height = 6, units = 'in', res = 300)                

par(mfrow=c(3, 2))
par(mar=c(1.5, 2.0, 1.5, 0.75))
par(mgp=c(2, 0.5, 0))
par(oma=c(3, 2, 0, 0))

layout(matrix(c(1,1,2,2,3,3,4,4,0,5,5,0), byrow=T, nrow=3))
#summary(m.tree)
VegPlot(name.i[[1]], mods.out[[1]])
text(1.1, 675, expression("++"), cex=cy)
text(2.1, 625, expression("++"), cex=cy)
text(3.1, 555, expression("++"), cex=cy)
text(4.1, 655, expression("++"), cex=cy)
text(5.1, 515, expression("++"),  cex=cy)
text(6.1, 455, expression("++"),  cex=cy)

VegPlot(name.i[[2]], mods.out[[2]])
#summary(m.shrub)
text(1.0, 78, expression("++"),  cex=cy)
text(2.1, 85, expression("++"),  cex=cy)
text(3.1, 88, expression("++"),  cex=cy)
text(4.0, 65, expression("++"),  cex=cy)
text(5.0, 60, expression("++"),  cex=cy)
text(6.0, 62, expression("++"),  cex=cy)

VegPlot(name.i[[3]], mods.out[[3]])
#summary(m.slash)
text(1, 797,expression("--"),  cex=cy)
text(2, 856,expression("--"),  cex=cy)
text(3, 825,expression("--"),  cex=cy)
text(4, 842,expression("--"),  cex=cy)
text(5, 857,expression("--"),  cex=cy)
text(6, 902,expression("--"),  cex=cy)

VegPlot(name.i[[4]], mods.out[[4]], xax=FALSE)
#summary(m.lit)
text(1.1, 6.9, expression("--"),  cex=cy)
text(2, 5.5, expression("--"), cex=cy)
text(3, 4.9, expression("--"),  cex=cy)
text(4, 5.4, expression("--"),  cex=cy)
text(5, 4.7, expression("--"), cex=cy)
text(6, 4.9, expression("--"),  cex=cy)

VegPlot(name.i[[5]], mods.out[[5]], xax=TRUE)
text(2, 6.1, expression("-"),  cex=cy+0.2)
text(5, 7.1, expression("--"),  cex=cy)

#summary(m.rob)


#' Legend
cth<-1.5 # y-height of cut legend
text(4.0, cth,  expression("cut"), cex=1.5, pos=4)
points(3.4, cth,  pch=18, cex=1.8)
segments(3.0, cth, 3.75, cth, lty=2)

ch<-0.7 # height of control legend
text(4,ch,  expression("control"), cex=1.5, pos=4)
segments(3, ch,3.75, ch, col=gr, lty=1, lwd=2)
points(3.4,ch,  pch=21, col=gr, lwd=2,bg="white", cex=1.8)

mtext("Years since start of study", side = 1, line = 0, outer = TRUE, padj=.5)
# dev.off()
### end veg plot ####


#################################################################
##### Create table of betas and SEs #################################
dat.list <- list(g.bhco, g.bobo, g.ccsp, g.coye,g.eaki,g.sewr, g.ywar, 
                 d.cavity, d.woodall, d.wetall, d.grassall, d.grass)
spec.list <- c("g.bhco","g.bobo","g.ccsp", "g.coye","g.eaki","g.sewr","g.ywar", 
               "d.cavity", "d.woodall", "d.wetall", "d.grassall", "d.grass")
df <- data.frame(matrix(vector(),  length(dat.list), 7, 
                        dimnames=list(c(),  c("spec", "yr1","yr2", "yr3", "yr4", "yr5", "yr6"))), 
                 stringsAsFactors=F)

# loop all through the analysis and output table
for (i in 1:11){
  dat<-dat.list[[i]]
  yrXtrt <-geeglm(total ~ as.factor(year)*trtnum, data=dat, 
                  id=wpa, family=poisson("log"), corstr="exchangeable", na.action="na.omit")
  est.se<-round(summary(yrXtrt)$coefficients$Std.err, 3)
  est<-round(summary(yrXtrt)$coefficients$Estimate, 3)
  
  tmp<-c(spec.list[i],
         paste(est[9]," (",est.se[9], ")", sep=""), 
         paste(est[10]," (",est.se[10],")", sep=""), 
         paste(est[11]," (",est.se[11],")", sep=""),
         paste(est[12]," (",est.se[12],")", sep=""),
         paste(est[13]," (", est.se[13],")", sep=""), 
         paste(est[14]," (",est.se[14],")", sep=""))
  df[i, ]<-tmp
}
df