In luke-a-rogers/mmmstan: Fit Markov Movement Models (MMM) via CmdStan and CmdStanR

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(mmmstan)
fit_models <- FALSE
# fit_models <- TRUE

Load data

if (fit_models) {
  # Define tag paths
  path_released <- "~/github/sablefish-data/data/tags_released.rda"
  path_recovered <- "~/github/sablefish-data/data/tags_recovered.rda"
  # Assign values
  tags_released <- read_from_path(path = path_released)
  tags_recovered <- read_from_path(path = path_recovered)
  # Tags
  tag_data <- dplyr::left_join(
    x = tags_released,
    y = tags_recovered,
    by = c(
      "tag_id", 
      "date_released", 
      "region_released_3", 
      "region_released_6", 
      "region_released_8", 
      "size_released", 
      "source"
    )
  )

  # Region arguments -----------------------------------------------------------

  # # Eight regions
  # list_regions_8 <- list(ai = 1, bs = 2, wg = 3, cg = 4, eg = 5, se = 6, bc = 7, cc = 8)
  # movement_pattern <- 2
  # movement_allow <- matrix(c(5, 7, 7, 5), nrow = 2, byrow = TRUE)
  # movement_disallow <- NULL

  # Three regions
  list_regions_3 <- list(ak = 1, bc = 2, cc = 3)
  movement_pattern <- 2
  movement_allow <- NULL
  movement_disallow <- NULL

  # Step arguments -------------------------------------------------------------

  step_interval <- "quarter"
  # step_duration_max <- NULL
  step_duration_max <- 12

  # Days argument --------------------------------------------------------------

  days_duration_min <- 90

  # Size arguments -------------------------------------------------------------

  # For fit mean
  list_sizes <- list(sml = 400:800)
  # For fit size
  list_sizes_small_large <- list(s = 400:549, l = 550:800)

  # Year arguments -------------------------------------------------------------

  year_start <- 1979
  year_end <- 2017

  # Column name arguments ------------------------------------------------------

  colname_date_released <- "date_released"
  colname_date_recovered <- "date_recovered"
  colname_region_released <- "region_released_3"
  colname_region_recovered <- "region_recovered_3"
  colname_size_released <- "size_released"

  # Arguments with defaults ----------------------------------------------------

  # Movement step priors
  mu_movement_step_diag <- NULL
  sd_movement_step_diag <- NULL
  # Fishing rate priors
  mu_fishing_rate <- NULL
  cv_fishing_rate <- NULL
  # Selectivity priors
  mu_selectivity <- NULL
  sd_selectivity <- NULL
  # Fishing weight priors
  mu_fishing_weight <- NULL
  sd_fishing_weight <- NULL
  # Natural mortality rate priors
  mu_natural_mortality_rate <- NULL
  sd_natural_mortality_rate <- NULL
  # Fractional (per tag) reporting rate priors
  mu_reporting_rate <- NULL
  sd_reporting_rate <- NULL
  # Fractional (per tag) initial loss rate priors
  mu_initial_loss_rate <- 0.1
  sd_initial_loss_rate <- 0.01
  # Instantaneous ongoing loss rate priors
  mu_ongoing_loss_rate <- 0.02
  sd_ongoing_loss_rate <- 0.001
  # Selectivity
  mu_selectivity <- NULL
  cv_selectivity <- NULL
  # Dispersion priors
  mu_dispersion <- 1.0
  sd_dispersion <- 0.5
  # Tolerance values
  tolerance_expected <- 1e-12
  tolerance_fishing <-  1e-12
  # CmdStanR
  data <- NULL
  chains <- 1
  step_size <- 0.01
  adapt_delta <- 0.95
  iter_warmup <- 250
  iter_sampling <- 750
  max_treedepth <- 10
  threads_per_chain <- parallel::detectCores()/(2*chains)

  # Use reduce sum -------------------------------------------------------------

  use_reduce_sum <- FALSE
  # use_reduce_sum <- TRUE
  refresh <- 10
}

Fit mean

if (fit_models) {
  # Fit mean
  fit_mean <- mmmstan(
    tag_data = tag_data,
    # Tag arguments
    list_regions = list_regions_3,
    list_sizes = list_sizes,
    year_start = year_start,
    year_end = year_end,
    step_interval = step_interval,
    step_duration_max = step_duration_max,
    days_duration_min = days_duration_min,
    colname_date_released = colname_date_released, # "date_released",
    colname_date_recovered = colname_date_recovered, # "date_recovered",
    colname_region_released = colname_region_released, # "region_released",
    colname_region_recovered = colname_region_recovered, # "region_recovered",
    colname_size_released = colname_size_released, # "size_released",
    # Movement index
    movement_pattern = movement_pattern,
    movement_allow = movement_allow,
    movement_disallow = movement_disallow,
    # Use reduce sum
    use_reduce_sum = use_reduce_sum,
    refresh = refresh
  )
  # View
  tibble::view(fit_mean$summary$movement_rate)
  tibble::view(fit_mean$summary$fishing_rate)
  # Shiny
  shinystan::launch_shinystan(fit_mean$rstanfit)
  # Plot
  ggplot2::ggplot(
    data = fit_mean$summary$movement_rate,
    mapping = ggplot2::aes(
      x = 1,
      y = .data$mean
    )
  ) +
    ggplot2::geom_col() +
    ggplot2::facet_grid(
      rows = ggplot2::vars(.data$x),
      cols = ggplot2::vars(.data$y)
    ) +
    ggsidekick::theme_sleek()
}

Fit size

if (fit_models) {
  # Fit size
  fit_size <- mmmstan(
    tag_data = tag_data,
    # Tag arguments
    list_regions = list_regions_3,
    list_sizes = list_sizes_small_large,
    year_start = year_start,
    year_end = year_end,
    step_interval = step_interval,
    step_duration_max = step_duration_max,
    colname_date_released = colname_date_released, # "date_released",
    colname_date_recovered = colname_date_recovered, # "date_recovered",
    colname_region_released = colname_region_released, # "region_released",
    colname_region_recovered = colname_region_recovered, # "region_recovered",
    colname_size_released = colname_size_released, # "size_released",
    # Movement index
    movement_pattern = movement_pattern,
    movement_allow = movement_allow,
    movement_disallow = movement_disallow,
    # Selectivity
    # mu_selectivity = array(0.9, c(2, 3)),
    # Use reduce sum
    use_reduce_sum = use_reduce_sum,
    refresh = refresh
  )
  # View
  tibble::view(fit_size$summary$movement_rate)
  tibble::view(fit_size$summary$fishing_rate)
  tibble::view(fit_size$summary$selectivity)
  # Plot size
  ggplot2::ggplot(
    data = fit_size$summary$movement_rate,
    mapping = ggplot2::aes(
      x = .data$l,
      y = .data$mean
    )
  ) +
    ggplot2::geom_col() +
    ggplot2::facet_grid(
      rows = ggplot2::vars(.data$x),
      cols = ggplot2::vars(.data$y)
    ) +
    ggsidekick::theme_sleek()
}