Primary Investigator: Althea ArchMiller
Collaborators: J. Fieberg, B. Dorazio, K. St. Clair
Description This model takes multiple years of operational survey data and “fixed” sightability model parameters (i.e., \(\Phi^{(s)}=\beta_{0g}^{(s)},\beta_{1g}^{(s)}\)) to estimate abundance for each year using temporal model with exchangeable random effects for lambda and a temporal spline for psi
# Documentation-related
library(devtools)
library(knitr)
library(ezknitr)
# Analysis-related
library(splines)
library(foreach)
library(reshape2)
# Remove environment and set seed
remove(list=ls())
set.seed(4587)
Augmented operational survey data
load("data/augdata.Rdata")
Years of analysis
year <- unique(operdata.augmented$year)
Posteriors from the model
load("data/output_data/posteriors_TEM_1stHalf.Rdata")
load("data/output_data/posteriors_TEM_2ndHalf.Rdata")
posteriors.TEM <- rbind(posteriors.TEM.1of2, posteriors.TEM.2of2)
Operational survey sampling characteristics
srates <- read.csv("data/srates0416.csv")
srates <- srates[srates$year %in% year,]
srates <- srates[srates$stratum!=4,]
Year-specific strata and plots from DNR
plotnumbs <- read.csv("data/plots_and_strata2016.csv")
colnames(plotnumbs)[16] <- "STRATA2016"
yearnames <- c("STRATA2004", "STRATA2009", "STRATA2012", "STRATA2013",
"STRATA2014", "STRATA2005", "STRATA2006", "STRATA2007",
"STRATA2008", "STRATA2010", "STRATA2011", "STRATA2015",
"STRATA2016")
plot.info <- reshape(data = plotnumbs,
varying = yearnames,
direction = "long", sep = "")
# Remove plots that were removed from the sampling strata
# after 2009 (5 6 26 27 28 29 30 31 288 289 290)
plot.info2 <- plot.info[plot.info$STRATA!=0 &plot.info$STRATA!=4,]
thin.index <- seq(from=20, to = nrow(posteriors.TEM), by = 20)
thin.posts.TEM <- as.data.frame(posteriors.TEM[thin.index,])
nl <- nrow(thin.posts.TEM)
Data frame to hold tauhats
tauhat.TEM <- NULL
Sampling characteristics
nyrL <- length(year)
nstrataL <- length(unique(operdata.augmented$h))
Dummy variables
dummy1h <- c(1,0,0)
dummy2h <- c(0,1,0)
dummy3h <- c(0,0,1)
Number of augmented for each stratum
B <- c(40, 60, 100)
for(yy in 1:nyrL){
# Pull operational survey data for year yy
sampled.plots <- unique(operdata.augmented$plotid[operdata.augmented$yr==yy])
plots.in.yy <- plot.info2$PLOTNUMBER[plot.info2$time==year[yy]]
unsampled.plots <- plots.in.yy[! plots.in.yy %in% sampled.plots]
# Define the number of unsampled plots
Nmiss <- length(unsampled.plots)
# Stratum to which each unsampled plot belongs
unsampled.h <- plot.info2$STRATA[plot.info2$time==year[yy]&
plot.info2$PLOTNUMBER %in% unsampled.plots]
# Number to augment for each unsampled plot based on stratum
naugMiss <- B[unsampled.h]
# Augment unsampled plotids and unsampled strata
unsamp.aug.plots <- rep(unsampled.plots[unsampled.h==1], each=B[1])
unsamp.aug.plots <- c(unsamp.aug.plots,
rep(unsampled.plots[unsampled.h==2], each=B[2]))
unsamp.aug.plots <- c(unsamp.aug.plots,
rep(unsampled.plots[unsampled.h==3], each=B[3]))
unsamp.aug.plots <- sort(unsamp.aug.plots)
unsamp.aug.strata <- unsampled.h[match(unsamp.aug.plots, unsampled.plots)]
# number of augmented groups
naug <- length(unsamp.aug.plots)
# Column names of lambda for year==yy
colname.lam <- c(
ifelse(yy<10, paste0("lam10",yy), paste0("lam1",yy)),
ifelse(yy<10, paste0("lam20",yy), paste0("lam2",yy)),
ifelse(yy<10, paste0("lam30",yy), paste0("lam3",yy))
)
posts.lam.TEM <- thin.posts.TEM[,colname.lam]
# Column names of a.psi for year==yy
colname.apsi <- c(
ifelse(yy<10, paste0("a.psi10",yy), paste0("a.psi1",yy)),
ifelse(yy<10, paste0("a.psi20",yy), paste0("a.psi2",yy)),
ifelse(yy<10, paste0("a.psi30",yy), paste0("a.psi3",yy))
)
posts.apsi.TEM <- thin.posts.TEM[,colname.apsi]
# Column names of b.psi for year==yy
colname.bpsi <- c(
ifelse(yy<10, paste0("b.psi10",yy), paste0("b.psi1",yy)),
ifelse(yy<10, paste0("b.psi20",yy), paste0("b.psi2",yy)),
ifelse(yy<10, paste0("b.psi30",yy), paste0("b.psi3",yy))
)
posts.bpsi.TEM <- thin.posts.TEM[,colname.bpsi]
# Column name for tau.sampled for year==yy
colname.tau <- ifelse(yy<6, paste0("tau.samp0",yy+4), paste0("tau.samp",yy+4))
######### for-loop
tauhat.TEM.temp <- foreach(ii=1:nl, .combine = "c") %do% {
# Psi for unsampled plots
psiM <- rbeta(
Nmiss,
as.numeric(posts.apsi.TEM[ii,unsampled.h]),
as.numeric(posts.bpsi.TEM[ii,unsampled.h])
)
# expand plot-level Psi values for each "group"
psiM.pergroup <- rep(psiM, times=naugMiss)
ynew <- rpois(naug, as.numeric(posts.lam.TEM[ii,unsamp.aug.strata]))
inpop <- rbinom(naug, 1, psiM.pergroup)
totM <- inpop*(ynew+1)
thin.posts.TEM[ii,colname.tau] + sum(totM)
}
# Merge with previous years' tauhats
tauhat.TEM <- c(tauhat.TEM, tauhat.TEM.temp)
}
Convert tauhats to data.frame and add years
tauhat.TEM <- as.data.frame(tauhat.TEM)
colnames(tauhat.TEM) <- "tauhat"
tauhat.TEM$Method <- "TEM"
tauhat.TEM$Year <- rep(year, each = nl)
save(tauhat.TEM, file = "data/output_data/tauhats_TEM.Rdata")
Session Information
devtools::session_info()
## Session info --------------------------------------------------------------
## setting value
## version R version 3.3.2 (2016-10-31)
## system x86_64, darwin13.4.0
## ui RStudio (1.0.136)
## language (EN)
## collate en_US.UTF-8
## tz America/Chicago
## date 2017-03-30
## Packages ------------------------------------------------------------------
## package * version date source
## abind 1.4-5 2016-07-21 CRAN (R 3.3.0)
## acepack 1.4.1 2016-10-29 CRAN (R 3.3.0)
## assertthat 0.1 2013-12-06 CRAN (R 3.3.0)
## backports 1.0.5 2017-01-18 CRAN (R 3.3.2)
## base64enc 0.1-3 2015-07-28 CRAN (R 3.3.0)
## boot 1.3-18 2016-02-23 CRAN (R 3.3.2)
## checkmate 1.8.2 2016-11-02 CRAN (R 3.3.0)
## cluster 2.0.5 2016-10-08 CRAN (R 3.3.2)
## coda * 0.19-1 2016-12-08 CRAN (R 3.3.2)
## codetools 0.2-15 2016-10-05 CRAN (R 3.3.2)
## colorspace 1.3-2 2016-12-14 CRAN (R 3.3.2)
## data.table 1.10.0 2016-12-03 CRAN (R 3.3.2)
## DBI 0.5-1 2016-09-10 CRAN (R 3.3.0)
## devtools * 1.12.0 2016-06-24 CRAN (R 3.3.0)
## digest 0.6.12 2017-01-27 CRAN (R 3.3.2)
## doBy * 4.5-15 2016-03-31 CRAN (R 3.3.0)
## dplyr * 0.5.0 2016-06-24 CRAN (R 3.3.0)
## evaluate 0.10 2016-10-11 CRAN (R 3.3.0)
## ezknitr * 0.6 2016-09-16 CRAN (R 3.3.0)
## foreach * 1.4.3 2015-10-13 CRAN (R 3.3.0)
## foreign 0.8-67 2016-09-13 CRAN (R 3.3.2)
## Formula * 1.2-1 2015-04-07 CRAN (R 3.3.0)
## ggplot2 * 2.2.1 2016-12-30 CRAN (R 3.3.2)
## ggthemes * 3.3.0 2016-11-24 CRAN (R 3.3.2)
## gridExtra * 2.2.1 2016-02-29 CRAN (R 3.3.0)
## gtable 0.2.0 2016-02-26 CRAN (R 3.3.0)
## highr 0.6 2016-05-09 CRAN (R 3.3.0)
## Hmisc * 4.0-2 2016-12-31 CRAN (R 3.3.2)
## htmlTable 1.9 2017-01-26 CRAN (R 3.3.2)
## htmltools 0.3.5 2016-03-21 CRAN (R 3.3.0)
## htmlwidgets 0.8 2016-11-09 CRAN (R 3.3.2)
## iterators 1.0.8 2015-10-13 CRAN (R 3.3.0)
## knitr * 1.15.1 2016-11-22 CRAN (R 3.3.2)
## labeling 0.3 2014-08-23 CRAN (R 3.3.0)
## lattice * 0.20-34 2016-09-06 CRAN (R 3.3.2)
## latticeExtra 0.6-28 2016-02-09 CRAN (R 3.3.0)
## lazyeval 0.2.0 2016-06-12 CRAN (R 3.3.0)
## magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
## markdown 0.7.7 2015-04-22 CRAN (R 3.3.0)
## MASS 7.3-45 2016-04-21 CRAN (R 3.3.2)
## Matrix 1.2-7.1 2016-09-01 CRAN (R 3.3.2)
## memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
## mime 0.5 2016-07-07 CRAN (R 3.3.0)
## munsell 0.4.3 2016-02-13 CRAN (R 3.3.0)
## nnet 7.3-12 2016-02-02 CRAN (R 3.3.2)
## plyr 1.8.4 2016-06-08 CRAN (R 3.3.0)
## R.methodsS3 1.7.1 2016-02-16 CRAN (R 3.3.0)
## R.oo 1.21.0 2016-11-01 CRAN (R 3.3.0)
## R.utils 2.5.0 2016-11-07 CRAN (R 3.3.0)
## R2jags * 0.5-7 2015-08-23 CRAN (R 3.3.0)
## R2WinBUGS 2.1-21 2015-07-30 CRAN (R 3.3.0)
## R6 2.2.0 2016-10-05 CRAN (R 3.3.0)
## RColorBrewer 1.1-2 2014-12-07 CRAN (R 3.3.0)
## Rcpp 0.12.9 2017-01-14 CRAN (R 3.3.2)
## reshape2 * 1.4.2 2016-10-22 CRAN (R 3.3.0)
## rjags * 4-6 2016-02-19 CRAN (R 3.3.0)
## rmarkdown 1.3 2016-12-21 CRAN (R 3.3.2)
## rpart 4.1-10 2015-06-29 CRAN (R 3.3.2)
## rprojroot 1.2 2017-01-16 CRAN (R 3.3.2)
## rstudioapi 0.6 2016-06-27 CRAN (R 3.3.0)
## scales 0.4.1 2016-11-09 cran (@0.4.1)
## SightabilityModel * 1.3 2014-10-04 CRAN (R 3.3.0)
## stringi 1.1.2 2016-10-01 CRAN (R 3.3.0)
## stringr 1.1.0 2016-08-19 CRAN (R 3.3.0)
## survival * 2.40-1 2016-10-30 CRAN (R 3.3.0)
## tibble 1.2 2016-08-26 CRAN (R 3.3.0)
## tidyr * 0.6.1 2017-01-10 CRAN (R 3.3.2)
## withr 1.0.2 2016-06-20 CRAN (R 3.3.0)
## yaml 2.1.14 2016-11-12 CRAN (R 3.3.2)
spun with: ezknitr::ezspin(“programs/ms_programs/e_calculate_TS_model_abundance.R”, out_dir = “output”, fig_dir = “figures”, keep_md = F) #knitr::opts_chunk$set(eval = F)