Fix 2022update

[kellijohnson-NOAA]

Implements

• change colors in warmer waters to be red thanks to @chantelwetzel-noaa
• saves output with error message if SDs are not estimable
• saves original Table_for_SS3.csv in addition to Index
• closes Error in amend_output(fit = NULL) fit is missing #41

Still need to add

• Add tweedie distribution #43

[kellijohnson-NOAA]

@chantelwetzel-noaa I compared runs for sablefish to those done for the 2019 stock assessment and the biomass estimates were within 12%. I think this is pretty good because it was with a different kmeans file. I need to check I get the same results when using the same kmeans file and then I think we can merge this into main/master. What do you think?

[kellijohnson-NOAA]

@chantelwetzel-noaa I have been thinking about this and do we need to exponentiate the log(sd) and convert the units then take the log again because I for sure did NOT do that.

[chantelwetzel-noaa]

Thanks for the reminder about this pull request. Initially, I did not think we would need to but I just compared results ran in 2019 for petrale sole and now I am not so sure. The biomass estimates from the new run are consistently lower ~6% using which are not that concerning since a different kmeans file was use, but the standard errors in the "converted" results are a fair bit lower (30-50%). This makes me wonder if there is a conversion issue. Have you done a comparable sablefish run with the same kmeans file?

[kellijohnson-NOAA]

I think the best thing to do would be to roll the code for {VAST} back to a version that had output in mt and see how Jim originally did the conversion compared to what is being done now. Then mimic that rather than just trying to match output. I can do that today.

[chantelwetzel-noaa]

I agree. Thank you for working on this.

[kellijohnson-NOAA]

Here are the sections of code that are doing the calculations. {VAST} and {FishStatsUtils} assume that the units are kg unless otherwise stated. Thus we can use

subdata[, "Catch_mt"] <- subdata[, "Catch_KG"] / 1000
units(subdata[, "Catch_mt]) <- "t"

and change line 107 of VAST_do to

b_i = subdata[, "Catch_mt"],

This just scales the estimates in Index.csv by 1000 for Estimate and Std. Error of Estimate. Std. Error for ln(Estimate) is not changed at all.

Current James-Thorson-NOAA/FishStatsUtils@92a2cd9

Table = cbind(
  expand.grid(dimnames(par_hat[[index_name]])),
  "Units" = make_unit_label( u=units(par_hat[[index_name]]), lab="", parse=FALSE ),
  "Estimate" = as.vector(par_hat[[index_name]]),
  "Std. Error for Estimate" = as.vector(par_SE[[index_name]]),
  "Std. Error for ln(Estimate)" = as.vector(par_SE[[log_index_name]])
)
write.csv( Table, file=paste0(DirName,"/Index.csv"), row.names=FALSE)

James-Thorson-NOAA/FishStatsUtils@80d5cbe

Table = NULL
  for( cI in 1:TmbData$n_c ){
    Tmp = data.frame(
      "Year"=year_labels,
      "Unit"=1,
      "Fleet"=rep(strata_names,each=TmbData$n_t),
      "Estimate_metric_tons"=as.vector(Index_ctl[cI,,,'Estimate']),
      "SD_log"=as.vector(log_Index_ctl[cI,,,'Std. Error']),
      "SD_mt"=as.vector(Index_ctl[cI,,,'Std. Error'])
    )
    if (TmbData$n_c>1) {
      Tmp = cbind( "Category"=category_names[cI], Tmp)
    }
    Table = rbind( Table, Tmp )
  }
  if (!is.null(total_area_km2)) {
    Table = cbind(Table, "Naive_design-based_index"=Design_t)
  }
  write.csv( Table, file=paste0(DirName,"/Table_for_SS3.csv"), row.names=FALSE)