MSEtool: Management Strategy Evaluation Toolkit

Documented in SS2OM

# === OM specification using SS3 (Methot 2012) stock assessments ====================

#' Reads MLE estimates from Stock Synthesis file structure into an operating model using package r4ss.
#'
#'
#' @description A function that uses the file location of a fitted SS3 model including input files to population the
#' various slots of an operating model with MLE parameter estimates. The function mainly populates the Stock and Fleet portions
#' of the operating model; the user still needs to parameterize most of the observation and implementation portions of the operating model.
#' @param SSdir A folder with Stock Synthesis input and output files in it.
#' @param nsim The number of simulations to take for parameters with uncertainty (for OM@@cpars custom parameters).
#' @param proyears The number of projection years for MSE
#' @param reps The number of stochastic replicates within each simulation in the operating model.
#' @param maxF The maximum allowable F in the operating model.
#' @param seed The random seed for the operating model.
#' @param interval The interval at which management procedures will update the management advice in \link[DLMtool]{runMSE}, e.g., 1 = annual updates.
#' @param Obs The observation model (class Obs). This function only updates the catch and index observation error.
#' @param Imp The implementation model (class Imp). This function does not update implementation parameters.
#' @param import_mov Logical, whether to import movement matrix from the assessment.
#' @param gender An integer that indexes the sex for importing life history parameters (1 = usually female, 2 = usually male, 1:2 = mean across both sexes).
#' @param age_rec Integer for the age of recruitment. The default is 1 for DLMtool operating models. Generally, should not be changed.
#' @param silent Whether to silence messages to the console.
#' @param Name The name of the operating model
#' @param Source Reference to assessment documentation e.g. a url
#' @param Author Who did the assessment
#' @param report Logical, if TRUE, the function will run \link[DLMtool]{runMSE} to generate historical data from the operating model
#' to compare against SS output. A markdown report will be generated.
#' @param filename If \code{report = TRUE}, character string for the name of the markdown and HTML files.
#' @param dir If \code{report = TRUE}, the directory in which the markdown and HTML files will be saved.
#' @param open_file If \code{report = TRUE}, whether the HTML document is opened after it is rendered.
#' @param ... Arguments to pass to \link[r4ss]{SS_output}.
#' @note Currently supports versions of r4ss on CRAN (v.1.24) and Github (v.1.34-38).
#' @details The function generally uses values from the terminal year of the assessment for most life history parameters (maturity, M, etc). This function
#' does detect time-varying growth in the assessment and uses annual length/weight-at-age for historical years.
#' Selectivity is derived from the F-at-age matrix.
#' @return An object of class OM.
#' @author T. Carruthers and Q. Huynh
#' @export
#' @seealso \link{SS2Data}
#' @importFrom stats acf
#' @importFrom reshape2 melt
#' @importFrom dplyr summarise group_by pull left_join
SS2OM <- function(SSdir, nsim = 48, proyears = 50, reps = 1, maxF = 3, seed = 1, interval = 1,
                  Obs = DLMtool::Generic_Obs, Imp = DLMtool::Perfect_Imp,
                  import_mov = TRUE, gender = 1:2, age_rec = 1, silent = FALSE,
                  Name = "OM generated by SS2OM function", Source = "No source provided", Author = "No author provided",
                  report = FALSE, filename = "SS2OM", dir = tempdir(), open_file = TRUE, ...) {

  replist <- SS_import(SSdir, ...)

  season_as_years <- FALSE
  if(replist$nseasons == 1 && all(replist$seasduration < 1)) {
    if(!silent) message(paste("Season-as-years detected in SS model. There is one season in the year with duration of", replist$seasduration, "year."))
    season_as_years <- TRUE
    nseas <- 1/replist$seasduration
    if(!silent) message("DLMtool operating model is an annual model. Since the SS model is seasonal, we need to aggregate over seasons.\n")
  } else {
    nseas <- replist$nseasons
    if(nseas > 1) {
      if(!silent) message("DLMtool operating model is an annual model. Since the SS model is seasonal, we need to aggregate over seasons.\n")
    }
  }

  # Create OM object
  Stock <- new("Stock")
  Fleet <- new("Fleet")
  OM <- new("OM", Stock = Stock, Fleet = Fleet, Obs = Obs, Imp = Imp)
  OM@nsim <- nsim
  OM@proyears <- proyears
  OM@Name <- Name
  OM@interval <- interval

  OM@Source <- paste0(Source, ". Author: ", Author, ".")

  mainyrs <- replist$startyr:replist$endyr
  OM@nyears <- nyears <- ceiling(length(mainyrs) / ifelse(season_as_years, nseas, 1))

  seas1_yind_full <- expand.grid(nseas = 1:nseas, true_year = 1:nyears) # Group assessment years to true years
  seas1_yind_full$assess_year <- mainyrs
  seas1_yind <- which(seas1_yind_full$nseas == 1) # Assessment years that correspond to first season of a true year

  OM@maxF <- maxF
  OM@reps <- reps
  OM@seed <- seed
  OM@CurrentYr <- ifelse(season_as_years, nyears, replist$endyr)

  # === Stock parameters =======================================================
  if(!silent && replist$nsexes == 2) {
    message("2-sex SS model found.")
    if(length(gender) == 1 && gender == 1) message("Life history parameters for gender code 1 (usually female) will be used for the OM.")
    if(length(gender) == 1 && gender == 2) message("Life history parameters for gender code 2 (usually male) will be used for the OM.")
    if(all(gender == c(1,2))) message("Life history parameters will be averaged from both genders for the OM.")
  }

  # ---- Stock-recruit relationship ----
  # R0 is unfished recruits at age = age_rec
  if(all(match(c("Gender", "Year"), names(replist$natage_annual_1_no_fishery), nomatch = 0))) { # SS 3.24

    N_at_age <- reshape2::melt(replist$natage_annual_1_no_fishery, id.vars = c("Bio_Pattern", "Gender", "Year"),
                               variable.name = "Age", value.name = "N")
    N_at_age <- N_at_age[N_at_age$Year == mainyrs[1] & N_at_age$Age == age_rec * ifelse(season_as_years, nseas, 1), ]
    N_at_age <- N_at_age[vapply(N_at_age$Gender, "%in%", logical(1), table = gender), ]

  } else { # SS 3.30

    N_at_age <- reshape2::melt(replist$natage_annual_1_no_fishery, id.vars = c("Bio_Pattern", "Sex", "Yr"),
                               variable.name = "Age", value.name = "N")
    N_at_age <- N_at_age[N_at_age$Yr == mainyrs[1] & N_at_age$Age == age_rec * ifelse(season_as_years, nseas, 1), ]
    N_at_age <- N_at_age[vapply(N_at_age$Sex, "%in%", logical(1), table = gender), ]

  }

  R0 <- sum(N_at_age$N)

  # In season-as-year model, R0 is the seasonal rate of recruitment, must adjust for annual model
  OM@R0 <- R0 * ifelse(season_as_years, nseas, 1)
  if(!silent) {
    message("R0 = ", round(R0), " (unfished abundance at age ", age_rec , ")")
    if(replist$nsexes == 2 && all(gender == c(1,2))) message("R0 is the sum of abundances of both genders.")
  }

  # Steepness is now deterministic
  if(packageVersion("r4ss") == 1.24) {
    SR_ind <- match(mainyrs, replist$recruit$year)
    SSB <- replist$recruit$spawn_bio[SR_ind]
    SSB0 <- replist$derived_quants[replist$derived_quants$LABEL == "SPB_Virgin", 2]
  } else {
    SR_ind <- match(mainyrs, replist$recruit$Yr)
    SSB <- replist$recruit$SpawnBio[SR_ind]
    SSB0 <- replist$derived_quants[replist$derived_quants$Label == "SSB_Virgin", 2]
  }

  if(replist$SRRtype == 3 || replist$SRRtype == 6) { # Beverton-Holt SR
    SR <- "BH"
    OM@SRrel <- 1L
    hs <- replist$parameters[grepl("steep", rownames(replist$parameters)), ]$Value
    if(!silent) message("Beverton-Holt stock-recruit relationship found with steepness = ", hs, ".\n")
  } else if(replist$SRRtype == 2) {
    SR <- "Ricker"
    OM@SRrel <- 2L
    hs <- replist$parameters[grepl("SR_Ricker", rownames(replist$parameters)), ]$Value
    if(!silent) message("Ricker stock-recruit relationship found with steepness = ", steep$Value, ".\n")
  } else if(replist$SRRtype == 7) {
    SR <- "BH"
    OM@SRrel <- 1L

    s_frac <- replist$parameters$Value[replist$parameters$Label == "SR_surv_Sfrac"]
    Beta <- replist$parameters$Value[replist$parameters$Label == "SR_surv_Beta"]

    s0 <- 1/SpR0
    z0 <- -log(s0)
    z_min <- z0 * (1 - s_frac)

    hs <- 0.2 * exp(z0 * s_frac * (1 - 0.2 ^ Beta))

    if(!silent) {
      message("Survival-based stock-recruit relationship was detected with steepness = ", round(hs, 2), ".")
      message("As an approximation, a Beverton-Holt relationship is used with the same steepness value.")
    }

  } else {
    SR <- OM@SRrel <- 1L

    rec <- replist$recruit$pred_recr[SR_ind] # recruits to age 0
    SpR0 <- SSB0/(R0 * ifelse(season_as_years, nseas, 1))

    hs <- SRopt(1, SSB, rec, SpR0, plot = FALSE, type = ifelse(SR == 1, "BH", "Ricker"))

    if(!silent) {
      message("From r4ss, SRRtype = ", replist$SRRtype)
      message("Steepness value not found. By default, estimating Beverton-Holt steepness from R and SSB estimates.\n")
      message("Steepness = ", hs, ".\n")
    }
  }
  OM@h <- rep(hs, 2)

  # ---- Max-age ----
  growdat <- getGpars(replist)      # Age-specific parameters in endyr
  growdat <- do.call(rbind, growdat)
  if("int_Age" %in% names(growdat)) {
    ages <- unique(growdat$int_Age)
  } else {
    ages <- unique(growdat$Age)
  }

  if(replist$nsexes == 1) {
    growdat$Len_Mat[growdat$Len_Mat < 0] <- 1
    growdat$Age_Mat[growdat$Age_Mat < 0] <- 1
  }

  maxage <- floor(max(ages)/ifelse(season_as_years, nseas, 1))

  seas1_aind_full <- expand.grid(nseas = 1:nseas, true_age = 0:maxage)[1:length(ages), ] # Group assessment ages to true ages
  seas1_aind_full$assess_age <- ages
  seas1_aind <- which(seas1_aind_full$nseas == 1 & seas1_aind_full$true_age >= age_rec) # Age indices that correspond to season 1

  OM@maxage <- maxage
  OM@cpars$plusgroup <- rep(1L, OM@nsim)
  if(!silent) message("Max age is ", maxage, ".")

  # ---- Growth ----
  if("int_Age" %in% names(growdat)) { # SS 3.30?

    if(season_as_years) {
      growdat <- growdat[vapply(growdat$int_Age, "%in%", logical(1), ages[seas1_aind]), ]
    } else growdat <- growdat[growdat$int_Age >= age_rec, ]
    growdat <- growdat[vapply(growdat$Sex, "%in%", logical(1), gender), ]

    Len_age_terminal <- summarise(group_by(growdat, int_Age), LAA = mean(Len_Beg)) %>% pull(2)
    Wt_age_terminal <- summarise(group_by(growdat, int_Age), WAA = mean(Wt_Beg)) %>% pull(2)
    Mat_age_terminal <- summarise(group_by(growdat, int_Age), MAA = mean(Len_Mat[Len_Mat >= 0] * Age_Mat[Age_Mat >= 0])) %>% pull(2)

  } else {

    if(season_as_years) {
      growdat <- growdat[vapply(growdat$Age, "%in%", logical(1), ages[seas1_aind]), ]
    } else growdat <- growdat[growdat$Age >= age_rec, ]
    growdat <- growdat[vapply(growdat$Sex, "%in%", logical(1), gender), ]

    Len_age_terminal <- unlist(summarise(group_by(growdat, Age), LAA = mean(Len_Beg))[, 2])
    Wt_age_terminal <- unlist(summarise(group_by(growdat, Age), WAA = mean(Wt_Beg))[, 2])
    Mat_age_terminal <- unlist(summarise(group_by(growdat, Age), MAA = mean(Len_Mat[Len_Mat >= 0] * Age_Mat[Age_Mat >= 0]))[, 2])
  }


  if(!is.null(replist$growthvaries) && replist$growthvaries) { # For all years except terminal year?
    if(!silent) {
      message("Time-varying growth found.")
      message("Projection period growth assumed to be the same as that in terminal year.")
    }
    if(season_as_years) stop("Can't deal with season_as_years yet when growth is time-varying.")

    Len_age <- reshape2::melt(replist$growthseries, id.vars = c("Morph", "Yr", "Seas", "SubSeas"), variable.name = "Age",
                              value.name = "LAA")
    Len_age <- Len_age[vapply(Len_age$Yr, "%in%", logical(1), table = mainyrs), ] # Subset by year
    Len_age <- Len_age[as.numeric(Len_age$Age)-1 >= age_rec, ] # Subset by age >= age_rec
    Len_age <- Len_age[Len_age$Seas == 1 & Len_age$SubSeas == 1, ] # Subset by season 1

    Len_age <- Len_age[vapply(Len_age$Morph, "%in%", logical(1), table = gender), ] # Subset by gender

    Len_age <- summarise(group_by(Len_age, Yr, Age), LAA = mean(LAA))
    Len_age_timevarying <- reshape2::acast(Len_age, list("Age", "Yr"), value.var = "LAA")

    Len_age <- as.matrix(cbind(Len_age_timevarying, Len_age_terminal))
    Len_age2 <- array(NA, dim = c(maxage, nyears+proyears, nsim))
    Len_age2[, 1:nyears, ] <- Len_age
    Len_age2[, nyears + 1:proyears, ] <- Len_age[, nyears]

    OM@cpars$Len_age <- aperm(Len_age2, c(3, 1, 2))

  } else {
    Len_age2 <- array(NA, dim = c(maxage, nsim, nyears+proyears))
    Len_age2[, , 1:nyears] <- Len_age_terminal
    Len_age2[, , nyears + 1:proyears] <- Len_age2[, , nyears]

    OM@cpars$Len_age <- aperm(Len_age2, c(2, 1, 3))

  }
  if(!silent) message("Length-at-age found.")

  #if(replist$wtatage_switch) {
  #  stop("Found empirical weights at age")
  #  Wt_age_emp <- melt(replist$wtatage, id.vars = c("Yr", "Seas", "Sex", "Bio_Pattern", "BirthSeas", "Fleet"),
  #                     variable.name = "Age", value.name = "WAA")
  #  Wt_age_emp$Yr <- abs(Wt_age_emp$Yr)
  #}
  GP <- replist$Growth_Parameters   # Some growth parameters (presumably in endyr)
  if (!is.null(GP$Platoon)) {
    GP <- GP[GP$Platoon == 1, ]
  }
  muLinf <- mean(GP$Linf[gender], na.rm = TRUE)
  cvLinf <- mean(GP$CVmax[gender], na.rm = TRUE)



  if(cvLinf > 1) cvLinf <- cvLinf/muLinf
  OM@LenCV <- rep(cvLinf, 2)

  OM@Ksd <- OM@Kgrad <- OM@Linfsd <- OM@Linfgrad <- c(0, 0)
  OM@K <- OM@Linf <- OM@t0 <- c(0, 0) # not used - vB pars estimated from Len_age internally

  # Weight at age
  OM@a <- mean(replist$Growth_Parameters$WtLen1[gender], na.rm = TRUE)
  OM@b <- mean(replist$Growth_Parameters$WtLen2[gender], na.rm = TRUE)


  if(exists("Len_age_timevarying")) { # This is a better solution for weight at age when growth is time varying
    OM@cpars$Wt_age <- OM@a * OM@cpars$Len_age ^ OM@b
  } else {
    Wt_age2 <- array(NA, dim = c(maxage, nsim, nyears + proyears))
    Wt_age2[, , 1:nyears] <- Wt_age_terminal
    Wt_age2[, , nyears + 1:proyears] <- Wt_age2[, , nyears]
    OM@cpars$Wt_age <- aperm(Wt_age2, c(2, 1, 3)) # dims = nsim, max_age, nyears+proyears
  }
  if(!silent) message("Weight-at-age found.")

  # ---- Maturity ----
  Mat_age <- array(NA, dim = c(maxage, nsim, nyears+proyears))
  Mat_age[, , 1:nyears] <- Mat_age_terminal/max(Mat_age_terminal)
  Mat_age[, , nyears + 1:proyears] <- Mat_age[, , nyears]
  OM@cpars$Mat_age <- aperm(Mat_age, c(2, 1, 3))  # dims = nsim, max_age, nyears+proyears
  if(!silent) message("Maturity-at-age found.")

  OM@L50 <- OM@L50_95 <- c(0, 0) # calculated internally

  # ---- M-at-Age ----
  if(all(match(c("Gender", "Year"), names(replist$M_at_age), nomatch = 0))) { # SS 3.24

    M_at_age <- reshape2::melt(replist$M_at_age, id.vars = c("Bio_Pattern", "Gender", "Year"), variable.name = "Age",
                               value.name = "M")
    M_at_age <- M_at_age[as.numeric(M_at_age$Age)-1 >= age_rec * ifelse(season_as_years, nseas, 1), ]  # Subset by age >= age_rec
    M_at_age <- M_at_age[vapply(M_at_age$Year, "%in%", logical(1), mainyrs), ]  # Subset by year
    M_at_age <- M_at_age[vapply(M_at_age$Gender, "%in%", logical(1), table = gender), ] # Subset by gender

    if(season_as_years) { # Mean across sub-ages, then sum across seasons
      M_at_age$true_Age <- seas1_aind_full$true_age[match(M_at_age$Age, seas1_aind_full$assess_age)]
      M_at_age <- summarise(group_by(M_at_age, Gender, Year, true_Age), M = mean(M))

      M_at_age$true_Year <- seas1_yind_full$true_year[match(M_at_age$Year, seas1_yind_full$assess_year)]
      M_at_age <- summarise(group_by(M_at_age, Gender, true_Year, true_Age), M = sum(M))

      M_ageArray <- reshape2::acast(M_at_age, list("true_Age", "true_Year"), value.var = "M")
    } else {
      M_at_age <- summarise(group_by(M_at_age, Year, Age), M = mean(as.numeric(M), na.rm = TRUE))
      M_ageArray <- reshape2::acast(M_at_age, list("Age", "Year"), value.var = "M")
    }


  } else { # SS 3.30
    M_at_age <- reshape2::melt(replist$M_at_age, id.vars = c("Bio_Pattern", "Sex", "Yr"), variable.name = "Age", value.name = "M")
    M_at_age <- M_at_age[vapply(M_at_age$Yr, "%in%", logical(1), table = mainyrs), ] # Subset by year
    M_at_age <- M_at_age[vapply(M_at_age$Sex, "%in%", logical(1), table = gender), ] # Subset by gender

    M_at_age <- M_at_age[as.numeric(M_at_age$Age)-1 >= age_rec * ifelse(season_as_years, nseas, 1), ]  # Subset by age >= age_rec

    M_at_age <- summarise(group_by(M_at_age, Yr, Age), M = mean(as.numeric(M), na.rm = TRUE))
    M_ageArray <- reshape2::acast(M_at_age, list("Age", "Yr"), value.var = "M")
  }

  if(ncol(M_ageArray) < nyears) M_ageArray <- cbind(M_ageArray, M_ageArray[, ncol(M_ageArray)])
  if(all(is.na(M_ageArray[maxage, ]))) M_ageArray[maxage, ] <- M_ageArray[maxage - 1, ]

  projM <- matrix(M_ageArray[, nyears], nrow = maxage, ncol = OM@proyears)
  M_ageArray <- cbind(M_ageArray, projM)
  M_ageArray <- replicate(OM@nsim, M_ageArray)
  OM@cpars$M_ageArray <- aperm(M_ageArray, c(3, 1, 2))

  OM@M <- M_ageArray[, nyears, 1]
  OM@M2 <- OM@M + 1e-6
  OM@Mexp <- OM@Msd <- c(0, 0)  # Time varying M would be in cpars
  if(!silent) message("Natural mortality found.")

  # ---- Depletion ----
  InitF <- replist$parameters$Value[grepl("InitF", replist$parameters$Label)]
  R_offset <- replist$parameters$Value[grepl("SR_R1_offset", replist$parameters$Label)]

  if(any(InitF > 0, na.rm = TRUE) || any(R_offset != 0, na.rm = TRUE)) {
    initD <- SSB[1]/SSB0
    OM@cpars$initD <- rep(initD, nsim)
    if(!silent) message("Initial depletion: OM@cpars$initD = ", round(initD, 2))
  }

  OM@D <- rep(replist$current_depletion, 2)
  if(!silent) message("Current depletion: OM@D = ", round(replist$current_depletion, 2), "\n")

  # ---- Recruitment deviations ----
  if(season_as_years) {
    replist$recruit$true_Yr <- seas1_yind_full$true_year[match(replist$recruit$Yr, seas1_yind_full$assess_year)]
    recruit <- summarise(group_by(replist$recruit, true_Yr), exp_recr = sum(exp_recr), pred_recr = sum(pred_recr)) # Need to sum over season_as_years
    hist_dev <- c(rep(1, maxage - 1), recruit$pred_recr[!is.na(recruit$true_Yr)]/recruit$exp_recr[!is.na(recruit$true_Yr)])
    dev_for_AC <- log(hist_dev)

  } else {

    year_first_rec_dev <- mainyrs[1] - maxage
    rec_ind <- match(year_first_rec_dev:(max(mainyrs) - age_rec * ifelse(season_as_years, nseas, 1)), replist$recruit$Yr)
    hist_dev <- replist$recruit$pred_recr[rec_ind]/replist$recruit$exp_recr[rec_ind] # In normal space
    hist_dev[is.na(hist_dev)] <- 1
    dev_for_AC <- replist$recruit$dev[rec_ind] # In log-space

  }

  dev_for_AC <- dev_for_AC[!is.na(dev_for_AC)]
  if(all(dev_for_AC == 0)) {
    if(!silent) message("Note: no recruitment deviates appear to be estimated.")
    AC <- 0
  } else {
    AC <- acf(dev_for_AC[dev_for_AC != 0], plot = FALSE)$acf[2, 1, 1]
  }
  if(is.na(AC)) AC <- 0
  if(!silent) message("Recruitment autocorrelation for projection period is estimated from historical recruitment deviations. OM@AC = ", round(AC, 3), ".")

  procsd <- replist$sigma_R_in
  procmu <- -0.5 * procsd^2 * (1 - AC/sqrt(1 - AC^2)) # adjusted log normal mean with AC
  Perr_hist <- matrix(hist_dev, nrow = nsim, ncol = maxage + nyears - 1, byrow = TRUE) # Historical recruitment is deterministic
  Perr_proj <- matrix(rnorm(proyears * nsim, rep(procmu, each = nsim),
                            rep(procsd, each = nsim)), nrow = nsim, ncol = proyears) # Sample recruitment for projection

  if(AC != 0) { # Add autocorrelation to projection recruitment
    Perr_proj[, 1] <- AC * Perr_hist[, ncol(Perr_hist)] + Perr_proj[, 1] * sqrt(1 - AC^2)
    for(y in 2:ncol(Perr_proj)) Perr_proj[, y] <- AC * Perr_proj[, y-1] + Perr_proj[, y] * sqrt(1 - AC^2)
  }

  OM@cpars$Perr_y <- cbind(Perr_hist, exp(Perr_proj))
  OM@Perr <- rep(procsd, 2) # Point estimate from assessment MLE
  OM@AC <- rep(AC, 2)
  if(!silent) message("Recruitment deviates found and future deviates sampled.")

  # ---- Movement modelling ----
  OM@Frac_area_1 <- OM@Size_area_1 <- OM@Prob_staying <- rep(0.5, 2)
  if(import_mov && nrow(replist$movement) > 0) {
    movement <- replist$movement[replist$movement$Seas == 1 & replist$movement$Gpattern == 1, ]
    if(nrow(movement) == 0) movement <- replist$movement[replist$movement$Seas == 1 & replist$movement$GP == 1, ]

    nareas <- length(unique(movement$Source_area))
    if(!silent) message(nareas, " area model found. Parameterizing movement matrix.")

    full_movement <- movement[, grepl("age", names(movement)) & names(movement) != "minage" & names(movement) != "maxage"]

    nages <- ncol(full_movement)
    mov <- array(NA, c(nsim, nages, nareas, nareas))

    for(i in 1:nrow(full_movement)) {
      from <- movement$Source_area[i]
      to <- movement$Dest_area[i]

      for(j in 1:ncol(full_movement)) mov[1:nsim, j, from, to] <- full_movement[i, j]
    }
    mov[is.na(mov)] <- 0

    if(season_as_years) mov <- mov[, seas1_aind, , , drop = FALSE]

    OM@cpars$mov <- mov
  }

  # Fleet parameters ===========================================================
  # ---- Vulnerability ----
  if(all(match(c("Gender", "Year"), names(replist$Z_at_age), nomatch = 0))) { # SS 3.24

    Z_at_age <- reshape2::melt(replist$Z_at_age, id.vars = c("Bio_Pattern", "Gender", "Year"), variable.name = "Age", value.name = "Z")
    Z_at_age <- Z_at_age[as.numeric(Z_at_age$Age)-1 >= age_rec * ifelse(season_as_years, nseas, 1), ]  # Subset by age >= age_rec
    Z_at_age <- Z_at_age[vapply(Z_at_age$Year, "%in%", logical(1), table = mainyrs), ]  # Subset by year
    Z_at_age <- Z_at_age[vapply(Z_at_age$Gender, "%in%", logical(1), table = gender), ] # Subset by gender

    if(season_as_years) { # Mean across sub-ages, then sum across seasons
      Z_at_age$true_Age <- seas1_aind_full$true_age[match(Z_at_age$Age, seas1_aind_full$assess_age)]
      Z_at_age <- summarise(group_by(Z_at_age, Gender, Year, true_Age), Z = mean(Z))

      Z_at_age$true_Year <- seas1_yind_full$true_year[match(Z_at_age$Year, seas1_yind_full$assess_year)]
      Z_at_age <- summarise(group_by(Z_at_age, Gender, true_Year, true_Age), Z = sum(Z))

      Z_ageArray <- reshape2::acast(Z_at_age, list("true_Age", "true_Year"), value.var = "Z")
    } else {
      Z_at_age <- summarise(group_by(Z_at_age, Year, Age), Z = mean(Z))
      Z_ageArray <- reshape2::acast(Z_at_age, list("Age", "Year"), value.var = "Z")
    }

  } else { # SS 3.30
    Z_at_age <- reshape2::melt(replist$Z_at_age, id.vars = c("Bio_Pattern", "Sex", "Yr"), variable.name = "Age", value.name = "Z")
    Z_at_age <- Z_at_age[vapply(Z_at_age$Yr, "%in%", logical(1), table = mainyrs), ] # Subset by year
    Z_at_age <- Z_at_age[as.numeric(Z_at_age$Age)-1 >= age_rec, ]                    # Subset by age >= age_rec
    Z_at_age <- Z_at_age[vapply(Z_at_age$Sex, "%in%", logical(1), table = gender), ] # Subset by gender

    Z_at_age <- summarise(group_by(Z_at_age, Yr, Age), Z = mean(Z))
    Z_ageArray <- reshape2::acast(Z_at_age, list("Age", "Yr"), value.var = "Z")

  }

  if(ncol(Z_ageArray) < nyears) Z_ageArray <- cbind(Z_ageArray, Z_ageArray[, ncol(Z_ageArray)])

  F_ageArray <- Z_ageArray - M_ageArray[, 1:nyears, 1]
  if(all(is.na(F_ageArray[maxage, ]))) F_ageArray[maxage, ] <- F_ageArray[maxage - 1, ]
  F_ageArray[F_ageArray < 1e-8] <- 1e-8

  F_norm <- apply(F_ageArray, 2, function(x) x/max(x))

  V <- array(NA, dim = c(maxage, nyears + proyears, nsim))
  V[, 1:nyears, ] <- array(F_norm, dim = c(maxage, nyears, nsim))
  V[, nyears + 1:proyears, ] <- V[, nyears, ]

  Find <- apply(F_ageArray, 2, max, na.rm = TRUE) # get apical F

  OM@cpars$V <- aperm(V, c(3, 1, 2))
  if(!silent) message("Selectivity found.")

  OM@L5 <- OM@LFS <- OM@Vmaxlen <- c(0,0) # calculated internally
  OM@isRel <- "FALSE" # these are real lengths not relative to length at 50% maturity

  # ---- Fishing mortality rate index ----
  OM@cpars$Find <- matrix(Find, nsim, nyears, byrow = TRUE) # is only historical years
  if(!silent) message("Historical F found.")

  OM@Spat_targ <- rep(1, 2)

  if(season_as_years) OM@EffYears <- 1:nyears else OM@EffYears <- mainyrs

  OM@EffLower <- OM@EffUpper <- OM@cpars$Find[1, 1:nyears]
  OM@Esd <- OM@qcv <- c(0, 0)
  OM@qinc <- c(0, 0)

  OM@Period <- OM@Amplitude <- rep(NaN, 2)


  # Observation model parameters ==============================================================================
  CSD <- replist$catch_error
  if(all(is.na(CSD)) && packageVersion("r4ss") == 1.35) CSD <- replist$catch_se
  if(!all(is.na(CSD))) {
    OM@Cobs <- range(sqrt(exp(CSD[CSD > 0]^2) - 1), na.rm = TRUE)
    if(!silent) {
      message("\nRange of error in catch (OM@Cobs) based on range of catch standard deviations: ", paste(CSD, collapse = " "))
    }
  }

  # Index observations -------------------------------------------------------
  Ind <- SS2Data_get_index(replist, mainyrs, season_as_years, nseas, index_season = "mean")

  if(is.null(Ind)) {
    message("No indices found.")
    if(packageVersion("DLMtool") >= 5.4) {
      Data@AddInd <- Data@CV_AddInd <- Data@AddIndV <- array(NA, c(1, 1, 1))
    }
  } else {
    Data <- new("Data")
    message(length(Ind$Iname), " indices of abundance found:")
    message(paste(Ind$Iname, collapse = "\n"))

    if(packageVersion("DLMtool") >= "5.4.4") {

      Data@AddInd <- Ind$AddInd
      Data@CV_AddInd <- sqrt(exp(Ind$SE_AddInd^2) - 1)

      Data@AddIunits <- Ind$AddIunits
      Data@AddIndType <- Ind$AddIndType

      if(season_as_years) {
        AddIndV <- apply(Ind$AddIndV, 1, function(x) {
          xx <- data.frame(assess_age = as.numeric(names(x)), sel = x) %>% left_join(seas1_aind_full[, -1], by = "assess_age")
          xx_agg <- aggregate(xx$sel, by = list(age = xx$true_age), mean, na.rm = TRUE)
          xx_agg$x[xx_agg$age >= age_rec]
        }) %>% t()
      } else {
        AddIndV <- Ind$AddIndV
      }

      Data@AddIndV <- array(AddIndV, c(1, dim(AddIndV)))

      OM@cpars$Data <- Data

      message("Updated Data@AddInd, Data@CV_AddInd, Data@AddIndV.")
    } else {
      message("\n\n *** Update DLMtool to latest version (5.4.4+) in order to add indices to OM (via OM@cpars$Data). *** \n\n")
    }
  }

  if(report) {
    if(!silent) message("\nRunning historical simulations to compare SS output and OM conditioning...\n")
    Hist <- runMSE(OM, Hist = TRUE)

    rmd_file <- file.path(system.file(package = "MSEtool"), "rmarkdown_templates", "SS2OM.Rmd")
    rmd <- readLines(rmd_file)

    write(rmd, file = file.path(dir, paste0(filename, ".rmd")))

    if(!silent) message("Rendering markdown file to HTML: ", file.path(dir, paste0(filename, ".html")))

    rmarkdown::render(file.path(dir, paste0(filename, ".rmd")), "html_document", paste0(filename, ".html"), dir,
                      output_options = list(df_print = "paged"), quiet = TRUE)
    message("Rendering complete.")

    if(open_file) browseURL(file.path(dir, paste0(filename, ".html")))
  }

  return(OM)
}