BayesCorrOcc: Fit Bayesian Correlated Occupancy Model

#' PlotRoutes
#'
#' Plot map with BBS route locations
#' @param alpha 4-letter alpha code for species of interest
#' @param psi Should occupancy probability (in year 1) be plotted? (*doesn't work yet*)
#' @export

PlotRoutes <- function(alpha){
  dat <-  readRDS(here::here(paste0("inst/output/", alpha, "/bbs_data.rds")))

  rts <- data.frame(Longitude = dat$lon, Latitude = dat$lat, z = apply(dat$h, 1, function(x) max(x, na.rm = TRUE)))
  obs <- dplyr::filter(rts, z == 1)
  buff <- dplyr::filter(rts, z == 0)
  usa <- ggplot2::map_data("state")
  canada <- ggplot2::map_data("world", "Canada")
  mexico <- ggplot2::map_data("world", "Mexico")


  xmin <- min(dat$lon) - 2
  xmax <- max(dat$lon) + 2

  ymin <- min(dat$lat) - 2
  ymax <- max(dat$lat) + 2


  p <- ggplot() + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA), axis.title = element_blank())
  p <- p + coord_map(projection = "lambert", lat0 = 25, lat1 = 50, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
  p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
    geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  p <- p + geom_point(data = buff, aes(x = Longitude, y = Latitude), color = "grey70", size = 2)
  p <- p + geom_point(data = obs, aes(x = Longitude, y = Latitude), color = "#c45b4d", size = 3)
  p <- p + scale_x_continuous(breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                           to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
  p <- p + scale_y_continuous(breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                           to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
  p
}


#' PlotLat
#'
#' Plot graph of changes in latitude indices
#' @param alpha 4-letter alpha code for species of interest
#' @export

PlotLat <- function(alpha, ci = FALSE){
  indices <- read.csv(here::here(paste0('inst/output/', alpha, '/indices.csv')))
  lat.indices <- dplyr::filter(indices, ind == "n.lat" | ind == "s.lat" |
                               ind == "n.core" | ind == "s.core" | ind == "avg.lat")


  p <- ggplot(lat.indices, aes(x = Year, y = Est, group = ind, color = ind))
  p <- p + geom_line(aes(linetype = ind))
  p <- p + scale_y_continuous("Latitude")
  p <- p + scale_linetype_manual(values = c("solid", "dashed", "longdash", "dashed", "longdash"))
  p <- p + scale_color_manual(values = c("black", "grey35", "grey50", "grey35", "grey50"))
  p <- p + scale_x_continuous(breaks = seq(from = min(lat.indices$Year),
                                           to = max(lat.indices$Year),
                                           by = 2))
  p <- p + theme(legend.position = "none")

  if(ci){
    p <- p + geom_ribbon(aes(ymin = LCI, ymax = UCI), alpha = 0.2)
    p
  }else{
    p
  }
}

#' PlotLon
#'
#' Plot graph of changes in longitude indices
#' @param alpha 4-letter alpha code for species of interest
#' @export

PlotLon <- function(alpha, ci = FALSE){
  indices <- read.csv(here::here(paste0('inst/output/', alpha, '/indices.csv')))
  lon.indices <- dplyr::filter(indices, ind == "w.lon" | ind == "e.lon" |
                                 ind == "w.core" | ind == "e.core" | ind == "avg.lon")


  p <- ggplot(lon.indices, aes(x = Est, y = Year, group = ind, color = ind))
  p <- p + geom_path(aes(linetype = ind))
  p <- p + scale_x_continuous("(West)               Longitude                 (East)")
  p <- p + scale_linetype_manual(values = c("solid", "dashed", "longdash", "dashed", "longdash"))
  p <- p + scale_color_manual(values = c("black", "grey35", "grey50", "grey35", "grey50"))
  p <- p + scale_y_continuous(breaks = seq(from = min(lon.indices$Year),
                                           to = max(lon.indices$Year),
                                           by = 2))
  p <- p + theme(legend.position = "none")

  if(ci){
    p <- p + geom_ribbon(aes(xmin = LCI, xmax = UCI), alpha = 0.2)
    p
  }else{
    p
  }
}

#' PlotPsi
#'
#' Plot graph of changes in range-wide occupancy
#' @param alpha 4-letter alpha code for species of interest
#' @export

PlotPsi <- function(alpha){
  indices <- read.csv(here::here(paste0('inst/output/', alpha, '/indices.csv')))
  psi.indices <- dplyr::filter(indices, ind == "avg.psi")
  y.min <- 0
  y.max <- plyr::round_any(max(psi.indices$UCI), 0.1, f = ceiling)

  p <- ggplot(psi.indices, aes(x = Year, y = Est, group = ind))
  p <- p + geom_point()
  p <- p + geom_line()
  p <- p + geom_point(color = "white", size = 6)
  p <- p + geom_point(size = 4)
  p <- p + geom_ribbon(aes(ymin = LCI, ymax = UCI), alpha = 0.2)
  p <- p + scale_y_continuous("Occupancy", limits = c(y.min, y.max))
  p <- p + scale_x_continuous(breaks = seq(from = min(psi.indices$Year),
                                           to = max(psi.indices$Year),
                                           by = 2))
  p
}

#' MapPsi
#'
#' Annual occupancy probability maps
#' @param alpha 4-letter alpha code for species of interest
#' @param proj Should maps be projected? (Probably but much slower)
#' @export


MapPsi <- function(alpha, proj = TRUE){
  dat <-  readRDS(here::here(paste0("inst/output/", alpha, "/bbs_data.rds")))
  psi <- readRDS(here::here(paste0('inst/output/', alpha, '/psi.rds')))
  indices <- read.csv(here::here(paste0('inst/output/', alpha, '/indices.csv')))
  limits <- dplyr::filter(indices, ind %in% c("s.lat", "n.lat"))
  core <- dplyr::filter(indices, ind %in% c("s.core", "n.core"))
  center <- dplyr::filter(indices, ind %in% c("avg.lat"))

  usa <- ggplot2::map_data("state")
  canada <- ggplot2::map_data("world", "Canada")
  mexico <- ggplot2::map_data("world", "Mexico")

  xmin <- min(psi$Longitude) - 1
  xmax <- max(psi$Longitude) + 1

  ymin <- min(psi$Latitude) - 1
  ymax <- max(psi$Latitude) + 1


  if(proj){
    p <- ggplot() + facet_wrap(~Year)
    p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
      geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
      geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
    p <- p + geom_tile(data = psi, aes(x = Longitude, y = Latitude, fill = Psi))
    p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA,  color="#909090") +
      geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090")
    p <- p + geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = NA, color="#909090")
    p <- p + coord_map(projection = "lambert", lat0 = 25, lat1 = 50, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
    p <- p + geom_hline(data = limits, aes(yintercept = Est, group = ind), linetype = "dashed", color = "grey70")
    p <- p + geom_hline(data = core, aes(yintercept = Est, group = ind), linetype = "longdash", color = "grey70")
    p <- p + geom_hline(data = center, aes(yintercept = Est), color = "grey70")
    p <- p + scale_fill_continuous(low = "#F0F0F1", high = "#ef6e0b")
    p <- p + scale_x_continuous("Longitude",breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                                         to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
    p <- p + scale_y_continuous("Latitude", breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                                         to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
    p <- p + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA), axis.title = element_blank())
    p
  }else{
    p <- ggplot() + facet_wrap(~Year) + coord_fixed(ratio = 1.3, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
    p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
      geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
      geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
    p <- p + geom_raster(data = psi, aes(x = Longitude, y = Latitude, fill = Psi)) + facet_wrap(~Year)
    p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA,  color="#909090") +
      geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = NA, color="#909090")
    p <- p + geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
    p <- p + geom_hline(data = limits, aes(yintercept = Est, group = ind), linetype = "dashed", color = "grey70")
    p <- p + geom_hline(data = core, aes(yintercept = Est, group = ind), linetype = "longdash", color = "grey70")
    p <- p + geom_hline(data = center, aes(yintercept = Est), color = "grey70")
    p <- p + scale_fill_continuous(low = "#F0F0F1", high = "#ef6e0b")
    p <- p + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA), axis.title = element_blank())
    p <- p + scale_x_continuous("Longitude", breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                                          to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
    p <- p + scale_y_continuous("Latitude", breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                                         to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
    p
  }


}


#' MapDiff
#'
#' Initial/final & difference occupancy maps
#' @param alpha 4-letter alpha code for species of interest
#' @export


MapDiff <- function(alpha){
  psi <- readRDS(here::here(paste0('inst/output/', alpha, '/psi.rds')))

  psi1 <- dplyr::filter(psi, Year == min(Year))
  psil <- dplyr::filter(psi, Year == max(Year))
  psid <- dplyr::mutate(psil, Diff = Psi - psi1$Psi)
  psid <- dplyr::select(psid, -Psi)
  psid <- dplyr::rename(psid, Psi = Diff)

  psi2 <- dplyr::bind_rows(psi1, psid)
  psi2 <- dplyr::mutate(psi2, Period = rep(c("Initial", "Difference"), each = nrow(psi1)))

  usa <- ggplot2::map_data("state")
  canada <- ggplot2::map_data("world", "Canada")
  mexico <- ggplot2::map_data("world", "Mexico")

  xmin <- min(psi$Longitude) - 1
  xmax <- max(psi$Longitude) + 1

  ymin <- min(psi$Latitude) - 1
  ymax <- max(psi$Latitude) + 1

  p <- ggplot()
  p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
    geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  p <- p + geom_tile(data = psi1, aes(x = Longitude, y = Latitude, fill = Psi))
  p <- p + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA,  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = NA, color="#909090")
  p <- p + geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  p <- p + coord_map(projection = "lambert", lat0 = 25, lat1 = 50, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
  p <- p + scale_fill_continuous(low = "#F0F0F1", high = "#ef6e0b",  limits = c(0, 1), name =  "Occupancy \nProb.")
  p <- p + scale_x_continuous("Longitude",breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                                       to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
  p <- p + scale_y_continuous("Latitude", breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                                       to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
  p <- p + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA),
                 axis.title = element_blank(), legend.title = element_text(size = 10))
  p <- p + labs(title = paste0(unique(psi1$Year), " (Initial)"))

  q <- ggplot()
  q <- q + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
    geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  q <- q + geom_tile(data = psil, aes(x = Longitude, y = Latitude, fill = Psi))
  q <- q + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA,  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = NA, color="#909090")
  q <- q + geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  q <- q + coord_map(projection = "lambert", lat0 = 25, lat1 = 50, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
  q <- q + scale_fill_continuous(low = "#F0F0F1", high = "#ef6e0b", limits = c(0, 1), name =  "Occupancy \nProb.")
  q <- q + scale_x_continuous("Longitude",breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                                       to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
  q <- q + scale_y_continuous("Latitude", breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                                       to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
  q <- q + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA),
                 axis.title = element_blank(), legend.title = element_text(size = 10))
  q <- q + labs(title = unique(psil$Year))


  r <- ggplot()
  r <- r + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "#F0F0F1",  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = "#F0F0F1", color="#909090") +
    geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  r <- r + geom_tile(data = psid, aes(x = Longitude, y = Latitude, fill = Psi))
  r <- r + geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = NA,  color="#909090") +
    geom_polygon(data = canada, aes(x=long, y = lat, group = group), fill = NA, color="#909090")
  r <- r + geom_polygon(data = mexico, aes(x=long, y = lat, group = group),  fill = "#F0F0F1", color="#909090")
  r <- r + coord_map(projection = "lambert", lat0 = 25, lat1 = 50, xlim = c(xmin, xmax), ylim = c(ymin, ymax))
  r <- r + scale_fill_gradient2(low = "#c45b4d", mid ="#F0F0F1", high = "#268bd2", name = "Difference")
  r <- r + scale_x_continuous("Longitude",breaks = seq(from = 10*(xmin%/%10 + as.logical(xmin%%10)),
                                                       to = 10*(xmax%/%10 + as.logical(xmax%%10)), by = 10))
  r <- r + scale_y_continuous("Latitude", breaks = seq(from = 5*(ymin%/%5 + as.logical(ymin%%5)),
                                                       to = 5*(ymax%/%5 + as.logical(ymax%%5)), by = 10))
  r <- r + theme(panel.background = element_rect(fill = "#CDD2D4", color = NA),
                 axis.title = element_blank(), legend.title = element_text(size = 10))

  s <- gridExtra::grid.arrange(gridExtra::arrangeGrob(p, q, ncol = 2), r, ncol = 1)
  grid::grid.draw(s)
}

#' md2html
#'
#' Function to covert .md file to html (ht https://github.com/richfitz/modeladequacy)
#' @export

md2html <- function(filename, dest = NULL) {
  if(is.null(dest)){ dest <- paste0(tools::file_path_sans_ext(filename), ".html")}
  opts <- setdiff(markdownHTMLOptions(TRUE), 'base64_images')
  markdownToHTML(filename, dest, options = opts)
}

#' PostCheck
#'
#' Plot posterior predictive checks
#' @export

PostCheck <- function(alpha){
  ppc <- readRDS(here::here(paste0("inst/output/", alpha, "/ppc.rds")))

  ggplot(ppc$fit, aes(x = fit, y = fit.new)) + geom_point(size = 2, alpha = 0.5) +
    geom_abline() +
    scale_x_continuous("Discrepancy actual data") +
    scale_y_continuous("Discrepancy replicate data") +
    annotate("text", x = max(ppc$fit$fit) * 0.9, y = max(ppc$fit$fit.new) * 0.9,
             label = paste("p = ", round(ppc$p, digits = 2)), size = 10)

  }



#' PlotIndex
#'
#' Function to plot multi-species indices of range change
#' @export

PlotIndex <- function(group_name, lat = TRUE, start.year = 1972){
  indices <- read.csv(here::here(paste0("inst/output/indices/", group_name, ".csv")))

  if(lat){
    indices <- dplyr::filter(indices, grepl("lat", ind) & Year >= start.year)
    indices$ind <- droplevels(indices$ind)

    indices$ind <- dplyr::recode(indices$ind, s.lat = "Southern range limit",
                                 avg.lat = "Mean breeding latitude",
                                 n.lat = "Northern range limit")
    indices$ind <- factor(indices$ind, levels = c("Northern range limit", "Mean breeding latitude",  "Southern range limit"))
    tot_index <- dplyr::filter(indices, Species == "Composite")
    spp_index <- dplyr::filter(indices, Species != "Composite")

    p <- ggplot()
    p <- p + geom_line(data = spp_index, aes(x = Year, y = value,
                   group = interaction(ind, Species)), alpha = 0.15, size = 0.7)
    p <- p + scale_color_manual(values = rep("black", 5), guide = FALSE)
    p <- p + geom_line(data = tot_index, aes(x = Year, y = value, group = ind),
                       color = "#cb4b16", size = 1.25)
    p <- p + scale_y_continuous("Latitude")
    p <- p + facet_wrap(~ind, ncol = 1, scales = "free")

  }else{
    indices <- dplyr::filter(indices, grepl("lon", ind) & Year >= start.year)
    indices$ind <- droplevels(indices$ind)

    indices$ind <- dplyr::recode(indices$ind, w.lon = "Western range limit", avg.lon = "Mean breeding longitude", e.lon = "Eastern range limit")
    indices$ind <- factor(indices$ind, levels = c("Western range limit", "Mean breeding longitude",  "Eastern range limit"))

    tot_index <- dplyr::filter(indices, Species == "Composite")
    spp_index <- dplyr::filter(indices, Species != "Composite")

    p <- ggplot()
    p <- p + geom_path(data = spp_index, aes(x = value, y = Year,
                                             group = interaction(ind, Species)), alpha = 0.15, size = 0.7)
    p <- p + scale_color_manual(values = rep("black", 5), guide = FALSE)
    p <- p + geom_path(data = tot_index, aes(x = value, y = Year, group = ind),
                       color = "#cb4b16", size = 1.25)
    p <- p + scale_x_continuous("Longitude")
    p <- p + facet_wrap(~ind, scales = "free")
  }

  p

}