#' GetSummary
#'
#' Generate summary information for report
#' @export
GetSummary <- function(alpha){
common <- BayesCorrOcc::code_lookup$common[BayesCorrOcc::code_lookup$alpha == toupper(alpha)]
## Data, JAGS object, and occupancy estimates
occ <- readRDS(paste0('inst/output/', alpha, '/occ.rds'))
dat <- readRDS(paste0("inst/output/", alpha, "/bbs_data.rds"))
ppc <- readRDS(paste0("inst/output/", alpha, "/ppc.rds"))
years <- seq(from = dat$start_year, to = dat$end_year)
nYears <- dat$nYears
## Extract beta coef estimates and CI's
## Beta coefficients for psi, gamma, & epsilon
psi_beta_tab <- MakeBetatab(alpha)
## Beta coefficients for th, th0, p, & omega
p_beta_tab <- MakeBetatab(alpha, nuisance = TRUE)
## Indicator variable posterior means
psi_g_tab <- MakeGtab(alpha)
summ <- list(spp_name = common,
spp_alpha = alpha,
n.routes = sum(apply(dat$h, 1, max, na.rm = TRUE)),
n.buffer = dim(dat$h)[1] - sum(apply(dat$h, 1, max, na.rm = TRUE)),
psi.betas = psi_beta_tab,
p.betas = p_beta_tab,
g.tab = psi_g_tab,
p = ppc$p)
summ
}
#' MakeBetatab
#'
#' Make Beta table for including in report
#' @param nuisance If true, returns AIC table of theta, theta', p, & omega; if false, returns AIC table for psi, gamma & epsilon
MakeBetatab <- function(alpha, nuisance = FALSE){
sim_list <- readRDS(paste0("inst/output/", alpha, "/jags_fit.rds"))
if(!nuisance){
## Covert covs included in psi, gam, & eps models to factor, set levels
covs_use <- factor(c("tmp", "Twet", "Prec", "Pwarm", "dtr",
"sq_tmp", "sq_Twet", "sq_Prec", "sq_Pwarm", "sq_dtr"),
levels = c("tmp", "Twet", "Prec", "Pwarm", "dtr",
"sq_tmp", "sq_Twet", "sq_Prec", "sq_Pwarm", "sq_dtr"))
## Data frame containing beta coeffecients and se
beta_est <- matrix(NA, nrow = 10, ncol = 1)
colnames(beta_est) <- c("$\\beta$")#, "$\\gamma$", "$\\epsilon$")
## Fill in intercept values
# beta_est[1, 1] <- paste(trunc(sim_list$mean$b[1]*100)/100,
# " (", trunc(sim_list$q2.5$b[1]*100)/100,
# " -- ", trunc(sim_list$q97.5$b[1]*100)/100, ")", sep = "")
# beta_est[1, 2] <- paste(trunc(sim_list$mean$beta.gam0*100)/100,
# " (", trunc(sim_list$q2.5$beta.gam0*100)/100,
# " -- ", trunc(sim_list$q97.5$beta.gam0*100)/100, ")", sep = "")
# beta_est[1, 3] <- paste(trunc(sim_list$mean$beta.eps0*100)/100,
# " (", trunc(sim_list$q2.5$beta.eps0*100)/100,
# " -- ", trunc(sim_list$q97.5$beta.eps0*100)/100, ")", sep = "")
## Fill in coefficients for climate covariates
for(i in 1:10){
# Psi
beta_est[i, 1] <- paste(trunc(sim_list$mean$betaT[i]*100)/100,
" (", trunc(sim_list$q2.5$betaT[i]*100)/100,
" -- ", trunc(sim_list$q97.5$betaT[i]*100)/100, ")", sep = "")
# Gamma
# beta_est[i + 1, 2] <- paste(trunc(sim_list$mean$betaT.gam[i]*100)/100,
# " (", trunc(sim_list$q2.5$betaT.gam[i]*100)/100,
# " -- ", trunc(sim_list$q97.5$betaT.gam[i]*100)/100, ")", sep = "")
#
# Epsilon
# beta_est[i + 1, 3] <- paste(trunc(sim_list$mean$betaT.eps[i]*100)/100,
# " (", trunc(sim_list$q2.5$betaT.eps[i]*100)/100,
# " -- ", trunc(sim_list$q97.5$betaT.eps[i]*100)/100, ")", sep = "")
}
## Rename climate covariates
covs_use <- dplyr::recode(covs_use, tmp = "Temp", sq_tmp = "Temp^2",
Twet = "Temp, Wettest Qrt", sq_Twet = "Temp, Wettest Qrt^2",
Prec = "Precip", sq_Prec = "Precip^2",
Pwarm = "Precip, Warmest Qrt", sq_Pwarm = "Precip, Warmest Qrt^2",
dtr = "Diurnal temp range", sq_dtr = "Diurnal temp range^2")
## Covert to data frame, add intercept, covert to character, replace NA with "-"
beta_df <- as.data.frame(beta_est)
covs <- data.frame(cov = as.character(levels(covs_use[order(covs_use)]))[1:10])
beta_df <- dplyr::bind_cols(covs, beta_df)
# beta_df[, 2:4] <- as.character(unlist(beta_df[, 2:4]))
names(beta_df)[1] <- ""
}else{
## Data frame containing beta coeffecients and se for theta, theta', p, omega ----
beta_est <- matrix(NA, nrow = 5, ncol = 3)
colnames(beta_est) <- c( "$p$", "$\\theta$", "$\\theta'$")
# Fill in intercept values (if annual == TRUE, different intercept for each year so don't include)
beta_est[1, 1] <- paste(trunc(sim_list$mean$alpha0*100)/100,
" (", trunc(sim_list$q2.5$alpha0*100)/100, " -- ",
trunc(sim_list$q97.5$alpha0*100)/100,")", sep = "")
beta_est[1, 2] <- paste(trunc(sim_list$mean$xpsi[1]*100)/100,
" (", trunc(sim_list$q2.5$xpsi[1]*100)/100, " -- ",
trunc(sim_list$q97.5$xpsi[1]*100)/100,")", sep = "")
beta_est[1, 3] <- paste(trunc(sim_list$mean$xpsi[2]*100)/100,
" (", trunc(sim_list$q2.5$xpsi[2]*100)/100, " -- ",
trunc(sim_list$q97.5$xpsi[2]*100)/100,")", sep = "")
# Fill in novice observer & wind effects, random observer variance
beta_est[2, 1] <- paste(trunc(sim_list$mean$alpha1*100)/100,
" (", trunc(sim_list$q2.5$alpha1*100)/100, " -- ",
trunc(sim_list$q97.5$alpha1*100)/100,")", sep = "")
beta_est[3, 1] <- paste(trunc(sim_list$mean$alpha2*100)/100,
" (", trunc(sim_list$q2.5$alpha2*100)/100, " -- ",
trunc(sim_list$q97.5$alpha2*100)/100,")", sep = "")
beta_est[4, 1] <- paste(trunc(sim_list$mean$alpha3*100)/100,
" (", trunc(sim_list$q2.5$alpha3*100)/100, " -- ",
trunc(sim_list$q97.5$alpha3*100)/100,")", sep = "")
beta_est[5, 1] <- paste(trunc(sim_list$mean$sigma.obs*100)/100,
" (", trunc(sim_list$q2.5$sigma.obs*100)/100, " -- ",
trunc(sim_list$q97.5$sigma.obs*100)/100,")", sep = "")
## Covert to data frame, add intercept, covert to character, replace NA with "-"
beta_df <- as.data.frame(beta_est)
covs <- data.frame(cov = c("Intercept", "Novice Observer", "Wind", "Twedt", "/sigma^2_{Obs}"))
beta_df <- dplyr::bind_cols(covs, beta_df)
beta_df[, 2:4] <- as.character(unlist(beta_df[, 2:4]))
beta_df[is.na(beta_df)] <- "-"
names(beta_df)[1] <- ""
}
beta_df
}
#' MakeGtab
#'
#' Make indicator variable table for including in report
#' @param alpha 4-letter species code
MakeGtab <- function(alpha){
sim_list <- readRDS(paste0("inst/output/", alpha, "/jags_fit.rds"))
## Covert covs included in psi, gam, & eps models to factor, set levels
covs_use <- factor(c("tmp", "Twet", "Prec", "Pwarm", "dtr",
"sq_tmp", "sq_Twet", "sq_Prec", "sq_Pwarm", "sq_dtr"),
levels = c("tmp", "Twet", "Prec", "Pwarm", "dtr",
"sq_tmp", "sq_Twet", "sq_Prec", "sq_Pwarm", "sq_dtr"))
## Data frame containing beta coeffecients and se
beta_est <- matrix(NA, nrow = 10, ncol = 1)
colnames(beta_est) <- c("g")#, "$\\gamma$", "$\\epsilon$")
## Fill in coefficients for climate covariates
for(i in 1:10){
# Psi
beta_est[i, 1] <- trunc(sim_list$mean$g[i]*100)/100
# Gamma
# beta_est[i, 2] <- trunc(sim_list$mean$g.gam[i]*100)/100
# Epsilon
# beta_est[i, 3] <- trunc(sim_list$mean$g.eps[i]*100)/100
}
## Rename climate covariates
covs_use <- dplyr::recode(covs_use, tmp = "Temp", sq_tmp = "Temp^2",
Twet = "Temp, Wettest Qrt", sq_Twet = "Temp, Wettest Qrt^2",
Prec = "Precip", sq_Prec = "Precip^2",
Pwarm = "Precip, Warmest Qrt", sq_Pwarm = "Precip, Warmest Qrt^2",
dtr = "Diurnal temp range", sq_dtr = "Diurnal temp range^2")
## Covert to data frame, add intercept, covert to character, replace NA with "-"
beta_df <- as.data.frame(beta_est)
covs <- data.frame(cov = as.character(levels(covs_use[order(covs_use)]))[1:10])
beta_df <- dplyr::bind_cols(covs, beta_df)
# beta_df[, 2:4] <- as.character(unlist(beta_df[, 2:4]))
names(beta_df)[1] <- ""
beta_df
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.