##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##   Closed-population N-mixture model analysis for GRSP abundance
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## 7,230 6-min surveys conducted during breeding seasons of 2015 (n = 2,807) and 2016 (n = 4,423) in eastern Kansas. We counted 3,238 GRSPs
## (2015, n = 1,364; 2016, n = 1,874) during 1,887 surveys (2015, n = 772; 2016, n = 1,115). We 1,420 sites twice in 2015, and 
## 2,249 sites twice in 2016 (includes all sites from 2015); however, we occasionally could not visit sites due to poor road conditions
## or deteriorating weather. We visited sites in rounds. Round 1 occurred between 13 May and 20 June in 2015, and 20 May and 24 June
## in 2016; round 2 occurred between 17 June and 20 July in 2015, and 27 June and 29 July in 2016. Visits to the same site were always 
## seperated by >2 weeks. Some birds could have moved between primary sampling periods; therefore, we interpret abundance (lambda) 
## as the total number of individuals that were ever associated with a given site, rather than the number of individuals that permanently
## resided at a site during a breeding season (Chandler and Royle, 2015). We combined (i.e., stacked) encounter histories from each
## season and considered sites to be independent between years. 

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## PART I. Bring in data, view summary statistics and correlation matrices, create unmarked frame, and standardize predictor variables. 
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Bookkeeping 
rm(list = ls())

## Set working directory 
setwd("C:/Users/strixhunter/Documents/Research/HESP_GRSP_UPSA_GRPC/GRSP_NMix")

## Read in ch file
data <- scan("grsp_nmix_ch_final.csv", sep = ",", skip = 1L)

## Create matrix for encounter histories, which will accept missing data (NAs)
R <- 3669 #number of sites (1,420 in 2015 + 2,249 in 2016)
J <- 2 #number of visits
y <- matrix(data, nrow=R, ncol=J, byrow=TRUE)

## Bring in the site and obs covariates for full encounter histories
sitecovs = read.csv("grsp_nmix_sitecovs_final.csv")
obscovs = read.csv("grsp_nmix_obscovs_final.csv")

## Bring in files for assessing multicollinearity.

## This file only includes one record for each of 2,249 unique sites.
sitecovs2249 = read.csv("grsp_nmix_sitecovs_2249_uniquesites_for_corr_matrix.csv")

## View headers of input files
head(sitecovs)
head(obscovs)
head(sitecovs2249)

## Define categorical variables as factors for all files used during subsequent analyses
sitecovs$year <- as.factor(sitecovs$year)
obscovs$obs1 <- as.factor(obscovs$obs1)
obscovs$obs2 <- as.factor(obscovs$obs2)
obscovs$round1 <- as.factor(obscovs$round1)
obscovs$round2 <- as.factor(obscovs$round2)


## View summary statistics of unstandardized values associated with unique sites
summary(sitecovs2249)

## Create scatterplot and cov.coeff for predictor variables, based on 2,249 unique sites
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  # correlation coefficient
  r <- cor(x, y)
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste("r= ", txt, sep = "")
  text(0.5, 0.6, txt)
  
  # p-value calculation
  p <- cor.test(x, y)$p.value
  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
  txt2 <- paste("p= ", txt2, sep = "")
  if(p<0.01) txt2 <- paste("p= ", "<0.01", sep = "")
  text(0.5, 0.4, txt2)
}

## Site covs
corr.matrix.200 <- pairs(~ ta200 + ca200 + nwa200 + trees200, lower.panel = panel.cor, data=sitecovs2249)
corr.matrix.400 <- pairs(~ ta400 + ca400 + nwa400 + trees400, lower.panel = panel.cor, data=sitecovs2249)
corr.matrix.800 <- pairs(~ ta800 + ca800 + nwa800 + trees800, lower.panel = panel.cor, data=sitecovs2249)
corr.matrix.1600 <- pairs(~ ta1600 + ca1600 + nwa1600 + trees1600, lower.panel = panel.cor, data=sitecovs2249)

## Load packages for analyses
library(unmarked)
library(AICcmodavg)

#create unmarkedFramePCount dataframe for N-mixture models
grsp.umf <- unmarkedFramePCount(y=y, 
                                siteCovs=sitecovs[,c("year", "ta200", "ca200", "nwa200", "trees200", "ta400", "ca400", "nwa400", "trees400", 
                                                     "ta800", "ca800", "nwa800", "trees800", "ta1600", "ca1600", "nwa1600", "trees1600")],
                                obsCovs=list(season=obscovs[,c("round1", "round2")], observer=obscovs[,c("obs1","obs2")],
                                             wind=obscovs[,c("wind1","wind2")], tss=obscovs[,c("tss1","tss2")],temp=obscovs[,c("temp1", "temp2")]))

head(grsp.umf)
summary(grsp.umf)

## Standardize continuous predictor variables (after making the UMF)
siteCovs(grsp.umf)[,2:17] <- scale(siteCovs(grsp.umf)[,2:17])
obsCovs(grsp.umf)$wind <- scale(obsCovs(grsp.umf)$wind)
obsCovs(grsp.umf)$tss <- scale(obsCovs(grsp.umf)$tss)
obsCovs(grsp.umf)$temp <- scale(obsCovs(grsp.umf)$temp)

## View summary of UMF with standardized variables
summary(grsp.umf)

##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Part II. Fit models
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


## Identify the best mixture for modeling the raw encounter histories (constant model)
(fm1 <-  pcount(~1 ~1, grsp.umf, mixture = "P", se=TRUE, K=100)) 
(fm2 <-  pcount(~1 ~1, grsp.umf, mixture = "NB", se=TRUE, K=100))
(fm3 <-  pcount(~1 ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))

cand.mods <- list(fm1, fm2, fm3) 
mod.names <- c("fm1", "fm2", "fm3") 
aictable <- aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE, c.hat=1.15)
aictable
## Zero-inflated Poisson mixture is best, so proceed with it for subsequent models.

(global <-  pcount(~year+season+wind+temp+I(temp^2) ~(ca400+I(ca400^2))*trees400, mixture = "ZIP", grsp.umf, K=100)) 
## Goodness-of-Fit test for the top fitting model. Beware, with 1,000 simulations, this takes ~16 hours. 
## I have multiple equally-complex 'global' models because I am considering multiple spatial scales. I use the most
## complex model for lambda at the 400 m scale, which is the best scale. This is the top model. 
Nmix.gof.test(fm38, nsim = 500, plot.hist = TRUE) #c-hat = 1.15, p < 0.05


## Model detection probability
(fm4 <-  pcount(~year ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100)) 
(fm5 <-  pcount(~season ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))
(fm6 <-  pcount(~tss ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100)) ## tss is correlated with temp; temp performs better, so omit tss subsequently
(fm7 <-  pcount(~temp ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))
(fm8 <-  pcount(~temp+I(temp^2) ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))
(fm9 <-  pcount(~wind ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))
#(fm10 <- pcount(~observer ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100)) ## fails to converge, likely bc there are many observers, so I omit this effect from the global model

## Use confint(fm9, type="det") to check beta coeff confidence intervals for a given model
## Next, consider whether combinations of main effects of the informative variables above 
## provide a better fit than temp^2 alone. TSS is left out because it is correlated with temperature. 

cand.mods <- list(fm1, fm2, fm3, fm4, fm5, fm6, fm7, fm8, fm9) 
mod.names <- c("fm1", "fm2", "fm3", "fm4", "fm5", "fm6", "fm7", "fm8", "fm9") 
aictable <- aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE, c.hat=1.15)
aictable

(fm11 <- pcount(~(temp+I(temp^2))+wind ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100))  
    #discarded bc wind is an uninformative parameter 
#(fm12 <- pcount(~(temp+I(temp^2))+season ~1, grsp.umf, mixture = "ZIP", se=TRUE, K=100)) #discarded bc season is an uninformative parameter when temp^2 is included
    # I reconsider both wind and season as predictors of detectability after modeling lambda to see if they improve the final model

## Akaike weights for the four models that include temp+temp^2 sum to 1.0. 
## Moreover, beta slope cofficients overlap zero in fm11, fm12, and fm13, when they are combined with temp+temp^2.
## Thus, a non-linear effect of temperature is clearly best for predicting detection probability, so we 
## proceed with it for all subsequent models. 
  
## Now, run all possible models for abundance (lambda). But first, test whether abundance varied by year. 

(fm14 <-  pcount(~temp+I(temp^2) ~year, mixture = "ZIP", grsp.umf, K=100)) 

## Year is an uninformative paramater with beta slope coefficient CI's overlapping zero.

(fm15 <-  pcount(~temp+I(temp^2) ~ta200, mixture = "ZIP", grsp.umf, K=100)) 
(fm16 <-  pcount(~temp+I(temp^2) ~ta200+trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm17 <-  pcount(~temp+I(temp^2) ~ta200*trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm18 <-  pcount(~temp+I(temp^2) ~ta200+I(ta200^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm19 <-  pcount(~temp+I(temp^2) ~(ta200+I(ta200^2))+trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm20 <-  pcount(~temp+I(temp^2) ~(ta200+I(ta200^2))*trees200, mixture = "ZIP", grsp.umf, K=100)) 

(fm21 <-  pcount(~temp+I(temp^2) ~ca200, mixture = "ZIP", grsp.umf, K=100)) 
(fm22 <-  pcount(~temp+I(temp^2) ~ca200+trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm23 <-  pcount(~temp+I(temp^2) ~ca200*trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm24 <-  pcount(~temp+I(temp^2) ~ca200+I(ca200^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm25 <-  pcount(~temp+I(temp^2) ~(ca200+I(ca200^2))+trees200, mixture = "ZIP", grsp.umf, K=100)) 
(fm26 <-  pcount(~temp+I(temp^2) ~(ca200+I(ca200^2))*trees200, mixture = "ZIP", grsp.umf, K=100)) 

(fm27 <-  pcount(~temp+I(temp^2) ~ta400, mixture = "ZIP", grsp.umf, K=100)) 
(fm28 <-  pcount(~temp+I(temp^2) ~ta400+trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm29 <-  pcount(~temp+I(temp^2) ~ta400*trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm30 <-  pcount(~temp+I(temp^2) ~ta400+I(ta400^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm31 <-  pcount(~temp+I(temp^2) ~(ta400+I(ta400^2))+trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm32 <-  pcount(~temp+I(temp^2) ~(ta400+I(ta400^2))*trees400, mixture = "ZIP", grsp.umf, K=100)) 

(fm33 <-  pcount(~temp+I(temp^2) ~ca400, mixture = "ZIP", grsp.umf, K=100)) 
(fm34 <-  pcount(~temp+I(temp^2) ~ca400+trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm35 <-  pcount(~temp+I(temp^2) ~ca400*trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm36 <-  pcount(~temp+I(temp^2) ~ca400+I(ca400^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm37 <-  pcount(~temp+I(temp^2) ~(ca400+I(ca400^2))+trees400, mixture = "ZIP", grsp.umf, K=100)) 
(fm38 <-  pcount(~temp+I(temp^2) ~(ca400+I(ca400^2))*trees400, mixture = "ZIP", grsp.umf, K=100)) 

(fm39 <-  pcount(~temp+I(temp^2) ~ta800, mixture = "ZIP", grsp.umf, K=100)) 
(fm40 <-  pcount(~temp+I(temp^2) ~ta800+trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm41 <-  pcount(~temp+I(temp^2) ~ta800*trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm42 <-  pcount(~temp+I(temp^2) ~ta800+I(ta800^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm43 <-  pcount(~temp+I(temp^2) ~(ta800+I(ta800^2))+trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm44 <-  pcount(~temp+I(temp^2) ~(ta800+I(ta800^2))*trees800, mixture = "ZIP", grsp.umf, K=100)) 

(fm45 <-  pcount(~temp+I(temp^2) ~ca800, mixture = "ZIP", grsp.umf, K=100)) 
(fm46 <-  pcount(~temp+I(temp^2) ~ca800+trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm47 <-  pcount(~temp+I(temp^2) ~ca800*trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm48 <-  pcount(~temp+I(temp^2) ~ca800+I(ca800^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm49 <-  pcount(~temp+I(temp^2) ~(ca800+I(ca800^2))+trees800, mixture = "ZIP", grsp.umf, K=100)) 
(fm50 <-  pcount(~temp+I(temp^2) ~(ca800+I(ca800^2))*trees800, mixture = "ZIP", grsp.umf, K=100)) 

(fm51 <-  pcount(~temp+I(temp^2) ~ta1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm52 <-  pcount(~temp+I(temp^2) ~ta1600+trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm53 <-  pcount(~temp+I(temp^2) ~ta1600*trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm54 <-  pcount(~temp+I(temp^2) ~ta1600+I(ta1600^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm55 <-  pcount(~temp+I(temp^2) ~(ta1600+I(ta1600^2))+trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm56 <-  pcount(~temp+I(temp^2) ~(ta1600+I(ta1600^2))*trees1600, mixture = "ZIP", grsp.umf, K=100)) 

(fm57 <-  pcount(~temp+I(temp^2) ~ca1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm58 <-  pcount(~temp+I(temp^2) ~ca1600+trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm59 <-  pcount(~temp+I(temp^2) ~ca1600*trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm60 <-  pcount(~temp+I(temp^2) ~ca1600+I(ca1600^2), mixture = "ZIP", grsp.umf, K=100)) 
(fm61 <-  pcount(~temp+I(temp^2) ~(ca1600+I(ca1600^2))+trees1600, mixture = "ZIP", grsp.umf, K=100)) 
(fm62 <-  pcount(~temp+I(temp^2) ~(ca1600+I(ca1600^2))*trees1600, mixture = "ZIP", grsp.umf, K=100)) 

## Make candidate list for all models
cand.mods <- list(fm1, fm2, fm3, fm4, fm5, fm6, fm7, fm8, fm9, fm15, fm16, fm17, fm18, fm19, fm20,
                  fm21, fm22, fm23, fm24, fm25, fm26, fm27, fm28, fm29, fm30, fm31, fm32, fm33, fm34, fm35, fm36, fm37, fm38, fm39,
                  fm40, fm41, fm42, fm43, fm44, fm45, fm46, fm47, fm48, fm49, fm50, fm51, fm52, fm53, fm54, fm55, fm56, fm57, fm58, 
                  fm59, fm60, fm61, fm62)
mod.names <- c("fm1", "fm2", "fm3", "fm4", "fm5", "fm6", "fm7", "fm8", "fm9", "fm15", "fm16",
               "fm17", "fm18", "fm19", "fm20", "fm21", "fm22", "fm23", "fm24", "fm25", "fm26", "fm27", "fm28", "fm29", "fm30", "fm31",
               "fm32", "fm33", "fm34", "fm35", "fm36", "fm37", "fm38", "fm39", "fm40", "fm41", "fm42", "fm43", "fm44", "fm45", "fm46",
               "fm47", "fm48", "fm49", "fm50", "fm51", "fm52", "fm53", "fm54", "fm55", "fm56", "fm57", "fm58", "fm59", "fm60", "fm61",
               "fm62") 
aictable <- aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE, c.hat=1.15)
aictable

## The single-best best model includes an interaction between core grassland area (quadratic effect) and woody area at the 400 m radius scale.

## Last, we'll reconsider the effect of wind on detectability to see if it is more informative now that we 
## modeled variation in abundace. 

(fm63 <-  pcount(~temp+I(temp^2)+wind ~(ca400+I(ca400^2))*trees400, mixture = "ZIP", grsp.umf, K=100)) 

## Make candidate list for all models
cand.mods <- list(fm1, fm2, fm3, fm4, fm5, fm6, fm7, fm8, fm9, fm15, fm16, fm17, fm18, fm19, fm20,
                  fm21, fm22, fm23, fm24, fm25, fm26, fm27, fm28, fm29, fm30, fm31, fm32, fm33, fm34, fm35, fm36, fm37, fm38, fm39,
                  fm40, fm41, fm42, fm43, fm44, fm45, fm46, fm47, fm48, fm49, fm50, fm51, fm52, fm53, fm54, fm55, fm56, fm57, fm58, 
                  fm59, fm60, fm61, fm62, fm63)
mod.names <- c("fm1", "fm2", "fm3", "fm4", "fm5", "fm6", "fm7", "fm8", "fm9", "fm15", "fm16",
               "fm17", "fm18", "fm19", "fm20", "fm21", "fm22", "fm23", "fm24", "fm25", "fm26", "fm27", "fm28", "fm29", "fm30", "fm31",
               "fm32", "fm33", "fm34", "fm35", "fm36", "fm37", "fm38", "fm39", "fm40", "fm41", "fm42", "fm43", "fm44", "fm45", "fm46",
               "fm47", "fm48", "fm49", "fm50", "fm51", "fm52", "fm53", "fm54", "fm55", "fm56", "fm57", "fm58", "fm59", "fm60", "fm61",
               "fm62", "fm63") 
aictable <- aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE, c.hat=1.15)
aictable

## Prepare AIC table to be exported
aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE)
## Output model selection table to csv
write.csv(aictab(cand.set=cand.mods, modnames=mod.names, second.ord = FALSE), file = 'grsp_aictab_1Apr17.csv')


### mean abundance per site ###
nhat <- (sum(bup(ranef(fm63, K=100))))/3669*(1-plogis(-0.755))
   
### SE of mean abundance per site ###
set.seed(345)
(pb.N <- parboot(fm63, Nhat, nsim=200, report=10))
parboot
## SE = StdDev/sqrt(n)

## Estimates for Figure 4.
## Get point predictions for lambda at core area = 20 ha (z-scale = -0.249) and 40 ha (z-scale = 1.0) with either 
## woody area = 0 (z-scale = -0.72926), 1.76 (z-scale = -0.3936), or 5.35 (z-scale = 0.2945), which will then be used in figures.
## I chose woody area = 0, 1.76, or 5.35 based on the min, median, and 3rd quartile observered among all unique sites in our study
## at the 400 m radius scale. 

## 20 ha core area, 0 ha trees (min)
ppred.newdata.min <- data.frame(ca400 = -0.249, trees400 = -0.72926) 
ppred.min <- predict(fm63, type="state", newdata=ppred.newdata.min, backTran=FALSE, appendData=TRUE)
ppred.abun.min=exp(ppred.min$Predicted); ppred.abun.min
ppredSE.min <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.newdata.min)

## 20 ha core area, 1.76 trees (median)
ppred.newdata.md <- data.frame(ca400 = -0.249, trees400 = -0.3936) 
ppred.md <- predict(fm63, type="state", newdata=ppred.newdata.md, backTran=FALSE, appendData=TRUE)
ppred.abun.md=exp(ppred.md$Predicted)
ppredSE.md <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.newdata.md)

## 20 ha core area, 5.35 trees (3rd quartile)
ppred.newdata.3rd <- data.frame(ca400 = -0.249, trees400 = 0.2945) 
ppred.3rd <- predict(fm63, type="state", newdata=ppred.newdata.3rd, backTran=FALSE, appendData=TRUE)
ppred.abun.3rd=exp(ppred.3rd$Predicted)
ppredSE.3rd <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.newdata.3rd)

## 40 ha core area, 1.76 trees (median)
ppred.newdata.md <- data.frame(ca400 = 1.0, trees400 = -0.3936) 
ppred.md <- predict(fm63, type="state", newdata=ppred.newdata.md, backTran=FALSE, appendData=TRUE)
ppred.abun.md=exp(ppred.md$Predicted)
ppredSE.md <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.newdata.md)

## 40 ha core area, 5.35 trees (3rd quartile)
ppred.newdata.3rd <- data.frame(ca400 = 1.0, trees400 = 0.2945) 
ppred.3rd <- predict(fm63, type="state", newdata=ppred.newdata.3rd, backTran=FALSE, appendData=TRUE)
ppred.abun.3rd=exp(ppred.3rd$Predicted)
ppredSE.3rd <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.newdata.3rd)


## Estimates for Figure 5.
## Get point predictions for specific values associated with landscapes used in plate figure
## Landscapes used for plates, from left to right, top to bottom, are:
## Site ID 3893302; TGA = 25.4, CGA = 20.6 (z-scale = -0.21025), CGA/TGA = 0.81, WA = 0 (z-scale = -0.72926)
## Site ID 3893306; TGA = 35.5, CGA = 29.7 (z-scale = 0.353407), CGA/TGA = 0.83, WA = 0 (z-scale = -0.72926)
## Site ID 3831102; TGA = 25.5, CGA = 10.5 (z-scale = -0.83957), CGA/TGA = 0.41, WA = 0 (z-scale = -0.72926)
## Site ID 3891920; TGA = 35.9, CGA = 18.1 (z-scale = -0.36894), CGA/TGA = 0.50, WA = 0 (z-scale = -0.72926)
## Site ID 3831004; TGA = 25.0, CGA = 10.2 (z-scale = -0.86146), CGA/TGA = 0.41, WA = 5.5 (z-scale = 0.328042)
## Site ID 3891008; TGA = 35.5, CGA = 17.0 (z-scale = -0.43462), CGA/TGA = 0.48, WA = 5.4 (z-scale = 0.294485)

## Site ID 3893302
ppred.3893302 <- data.frame(ca400 = -0.21025, trees400 = -0.72926) 
ppred.3893302 <- predict(fm63, type="state", newdata=ppred.3893302, backTran=FALSE, appendData=TRUE)
ppred.abun.3893302=exp(ppred.3893302$Predicted)
ppredSE.3893302 <- predictSE(fm63, type = "response", backTran=TRUE, parm.type = "lambda", newdata = ppred.3893302)
ppred.abun.3893302

## Site ID 3893306
ppred.3893306 <- data.frame(ca400 = 0.353407, trees400 = -0.72926) 
ppred.3893306 <- predict(fm63, type="state", newdata=ppred.3893306, backTran=FALSE, appendData=TRUE)
ppred.abun.3893306=exp(ppred.3893306$Predicted)
ppredSE.3893306 <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.3893306)
ppred.abun.3893306

## Site ID 3831102
ppred.3831102 <- data.frame(ca400 = -0.83957, trees400 = -0.72926) 
ppred.3831102 <- predict(fm63, type="state", newdata=ppred.3831102, backTran=FALSE, appendData=TRUE)
ppred.abun.3831102=exp(ppred.3831102$Predicted)
ppredSE.3831102 <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.3831102)
ppred.abun.3831102

## Site ID 3891920
ppred.3891920 <- data.frame(ca400 = -0.36894, trees400 = -0.72926) 
ppred.3891920 <- predict(fm63, type="state", newdata=ppred.3891920, backTran=FALSE, appendData=TRUE)
ppred.abun.3891920=exp(ppred.3891920$Predicted)
ppredSE.3891920 <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.3891920)
ppred.abun.3891920

## Site ID 3831004
ppred.3831004 <- data.frame(ca400 = -0.86148, trees400 = 0.328042) 
ppred.3831004 <- predict(fm63, type="state", newdata=ppred.3831004, backTran=FALSE, appendData=TRUE)
ppred.abun.3831004=exp(ppred.3831004$Predicted)
ppredSE.3831004 <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.3831004)
ppred.abun.3831004

## Site ID 3891008
ppred.3891008 <- data.frame(ca400 = -0.43462, trees400 = 0.294485) 
ppred.3891008 <- predict(fm63, type="state", newdata=ppred.3891008, backTran=FALSE, appendData=TRUE)
ppred.abun.3891008=exp(ppred.3891008$Predicted)
ppredSE.3891008 <- predictSE(fm63, type = "response", backTran=FALSE, parm.type = "lambda", newdata = ppred.3891008)
ppred.abun.3891008