R/results.R

In sihoward/mossiemodel-dev: Temperature-dependent mosquito population model

#' Population size time series
#'
#' Plot time series of population size using model outputs.
#'
#' @param resdf Results dataframe generated by mosqmod::runModel().
#' @param selectPopn Character vector of population component to plot (options are "L", "L_1", "L_2", "L_3", "L_4", "L_5" or "M").
#' @param include_temp Include a plot facet of temperature time series (default = TRUE).
#' @param include_plasmod Include (default = TRUE) shaded regions on plot when
#'   temperatures are suitable for plasmodium development.
#' @param MTT Display minimum temperature threshold for mosquito development.
#' @param MTD Display minimum temperature threshold for plasmodium sporogeny.
#'
#' @importFrom rlang .data
#' @importFrom ggplot2 aes
#'
#' @export
plot_popn <- function(resdf,
                      selectPopn = c("L", "L_1", "L_2", "L_3", "L_4", "L_5", "M"),
                      include_temp = TRUE, include_plasmod = TRUE, MTT = NULL, MTD = NULL){

  # set global vars (to pass package checks)
  line <- xmax <- xmin <- yintercept <- NULL

  if(include_temp){
    # add temperature to columns in resdf to plot
    selectPopn <- c(selectPopn, "Tmean")

    if(!is.null(MTT) | !is.null(MTD)){

      # values available?
      ind <- c(!is.null(MTT), !is.null(MTD))

      # dataframe combining MTT and MTD
      mintemp_df <-
        data.frame(name = rep(selectPopn, 2),
                   yintercept = c(ifelse(selectPopn == "Tmean" & !is.null(MTT), MTT, NA),
                                  ifelse(selectPopn == "Tmean" & !is.null(MTD), MTD, NA)),
                   line = rep(c("mtt", "mtd"), each = length(selectPopn)))

      # create horizontal lines in temperature panel
      gg_mintempLine <-
        ggplot2::geom_hline(aes(yintercept = yintercept, linetype = line),
                            data = mintemp_df, show.legend = TRUE)
      gg_mintempScale <-
        ggplot2::scale_linetype_manual(name = "Min temperature",
                                       values = c('mtt' = "dashed", "mtd" = "dotted")[ind],
                                       labels = c('mtt' = "Plasmodium sporogeny",
                                                  'mtd' = "Mosquito development")[ind],
                                       guide = ggplot2::guide_legend(order = 3))
    } else {
      gg_mintempLine <- NULL
      gg_mintempScale <- NULL
    }

  }

  # stack selected columns
  d <- tidyr::pivot_longer(resdf, cols = dplyr::all_of(selectPopn))
  # respect plotting order of selectPopn
  d$name <- factor(d$name, levels = unique(d$name))
  # create factor for recorded vs projected temperature estimates
  d$Temperature <- factor(ifelse(d$source %in% "projected", "projected", "recorded"),
                          levels = c("recorded", "projected"))

  # rescale y-axis values for population sizes
  ## rows with population sizes
  popnRows <- !d$name %in% "Tmean"
  ## get max value
  maxval <- max(d$value[popnRows])
  ## potential y scale denominators
  yscale <- c(1, 100,1000,10000, 100000, 1000000, 10000000)
  ## pick denominator closest to population sizes
  ind <- (maxval - yscale > 0) & (maxval - yscale) == min(maxval - yscale[maxval - yscale > 0])
  ## divide values by denominator
  d$value[popnRows] <- d$value[popnRows] / yscale[ind]
  ## scale label
  if(yscale[ind] > 1){
    ylab_scaled <- paste0("Population size\n(",
                          format(yscale[ind], scientific = FALSE, big.mark = ","),
                          "s per square km)")
  } else {
    ylab_scaled <- paste0("Population size\n(number per square km)")
  }

  if(include_plasmod){
    # get temperature time series
    Tseq <- resdf[c("Date", "Tmean")]
    # check for non-contiguous dates and gaps
    if(any(diff(Tseq$Date) != 1)) stop("non-contiguous dates in temperature series passed to plasmod_devel()")
    # are thermal requirements met on each date?
    thermal_req <- plasmod_devel(Tseq$Tmean)[["thermal_req"]]

    # find start and stop dates when requirements are met
    startdate <- Tseq$Date[diff(c(FALSE, thermal_req, FALSE)) == 1]
    enddate <- Tseq$Date[diff(c(thermal_req, FALSE)) == -1]

    # combine development windows and repeat for each plot panel
    devel_windows <-
      data.frame(name = rep(selectPopn, each = length(startdate)),
                 xmin = rep(startdate, times = length(selectPopn)),
                 xmax = rep(enddate, times = length(selectPopn)),
                 ymin = -Inf, ymax = Inf)

    # store geoms to draw deveopment windows
    gg_DEVRect <-
      ggplot2::geom_rect(data = devel_windows,
                         ggplot2::aes(xmin = xmin, xmax = xmax, ymax = Inf, ymin = -Inf,
                                      fill = "reqmet"),
                         inherit.aes = FALSE, alpha = 0.15)

    gg_DEVscale <-
      ggplot2::scale_fill_manual(name = NULL,
                                 values = c('reqmet' = "tomato"),
                                 labels = c('reqmet' = "Avian malaria risk"),
                                 guide = ggplot2::guide_legend(order = 4, override.aes = list(linetype = 0)))
  } else {
    # store NULL when !include_plasmod
    gg_DEVRect <- NULL
    gg_DEVscale <- NULL
  }

  gg <-
    ggplot2::ggplot(d, ggplot2::aes(y = .data$value, x = .data$Date,
                                    color = .data$name)) +
    gg_DEVRect + gg_DEVscale +
    ggplot2::geom_line(linejoin = "round") +
    ggplot2::geom_line(ggplot2::aes(size = .data$Temperature), linejoin = "round") +
    ggplot2::scale_color_discrete(name = "Population",
                                  labels = function(x) c(L = "Larvae", M = "Adults", Tmean = "Temperature")[x],
                                  guide = ggplot2::guide_legend(order = 1)) +
    ggplot2::scale_size_discrete(range = c(1,0.5), guide = ggplot2::guide_legend(order = 2))


  if(!include_temp){
    gg <- gg + ggplot2::labs(y = ylab_scaled)
  } else {
    gg <- gg + ggplot2::facet_grid(name == "Tmean" ~ . , scale = "free_y",
                                   labeller = ggplot2::as_labeller(c('TRUE' = "Temperature\n(degrees C)",
                                                                     'FALSE' = ylab_scaled)),
                                   switch = "y") +
      ggplot2::labs(y = NULL) +
      gg_mintempLine + gg_mintempScale
  }

  # format x breaks and labels
  gg <-
    gg + ggplot2::scale_x_date(date_breaks = "6 months",
                               date_minor_breaks = "1 month",
                               label = function(x) format(x, "%b %Y"))

  return(gg)
}


#' Population size time series overlay
#'
#' Population size time series of model outputs with years overlaid.
#'
#' @param resdf Results dataframe generated by mosqmod::runModel().
#' @param selectPopn Character vector of population component to plot (options are "L", "L_1", "L_2", "L_3", "L_4", "L_5" or "M").
#'
#' @importFrom rlang .data
#'
#' @export
plot_popn_years <- function(resdf,
                       selectPopn = c("L", "L_1", "L_2", "L_3", "L_4", "L_5", "M")){

  # append last recorded temperature (removes gap when drawing geom_line)
  if(any(resdf$source != "projected")){
    last.recorded <- resdf[max(which(resdf$source != "projected")),]
    last.recorded$source <- "projected"
    resdf <- rbind(resdf, last.recorded)
  }

  # stack selected columns
  d <- tidyr::pivot_longer(resdf, cols = dplyr::all_of(selectPopn))
  # respect plotting order of selectPopn
  d$name <- factor(d$name, levels = unique(d$name))
  # calendar year
  d$Year <- as.numeric(format(d$Date, "%Y"))
  # Start years in July
  d$Year <- d$Year - (as.numeric(format(d$Date, "%m")) < 7)
  # drop any leap days
  d <- d[!grepl("02-29", format(d$Date, "%m-%d")),]
  # align start of each year with first date in sequence
  d <- dplyr::mutate(dplyr::group_by(d, .data$Year),
                     Date = min(d$Date) + (.data$Date - min(.data$Date)))
  # make year discrete
  d$Year <- factor(d$Year, levels = sort(unique(d$Year), decreasing = TRUE))
  # create factor for recorded vs projected temperature estimates
  d$Temperature <- factor(ifelse(d$source %in% "projected", "projected", "recorded"),
                          levels = c("recorded", "projected"))

  # rescale y-axis values for population sizes
  ## get max value
  maxval <- max(d$value)
  ## potential y scale denominators
  yscale <- c(1, 100,1000,10000, 100000, 1000000, 10000000)
  ## pick denominator closest to population sizes
  ind <- (maxval - yscale > 0) & (maxval - yscale) == min(maxval - yscale[maxval - yscale > 0])
  ## divide values by denominator
  d$value <- d$value / yscale[ind]
  ## scale label
  if(yscale[ind] > 1){
    ylab_scaled <- paste0("Population size\n(",
                          format(yscale[ind], scientific = FALSE, big.mark = ","),
                          "s per square km)")
  } else {
    ylab_scaled <- paste0("Population size\n(number per square km)")
  }

  # dataframe ordered by date
  d <- d[order(d$Date),]

  gg <- ggplot2::ggplot(d, ggplot2::aes(y = .data$value, x = .data$Date,
                                        color = .data$Year,
                                        size = .data$Temperature)) +
    ggplot2::geom_path(linejoin = "round") +
    ggplot2::labs(y = ylab_scaled) +
    ggplot2::scale_x_date(breaks = seq(min(d$Date), by = "3 month", length.out = 5),
                          date_minor_breaks = "1 month",
                          date_label = "%b") +
    ggplot2::scale_color_brewer(palette = "Dark2",
                                guide = ggplot2::guide_legend(override.aes = list(linewidth = 1))) +
    ggplot2::scale_size_discrete(breaks = c("recorded", "projected"), range = c(1,0.5))


  if(length(selectPopn) > 1){
    gg <- gg + ggplot2::facet_wrap(~name, scales = "free_y")
  }

  return(gg)
}