Skip to content

Commit

Permalink
Had issue w/ remotes install (probably cuz I broke GOA pollock)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jimianelli
jimianelli committed Feb 15, 2024
1 parent 254c5b9 commit b4901b1
Show file tree
Hide file tree
Showing 2 changed files with 8 additions and 202 deletions.
203 changes: 6 additions & 197 deletions content/AFSC-BSAI_AtkaMackerel.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,8 +333,9 @@ obj <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
# fitting the model
opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr,
control = list(eval.max = 8000, iter.max = 800))
#opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr,
# control = list(eval.max = 8000, iter.max = 800))
# method = "BFGS",
# control = list(maxit=1000000, reltol = 1e-15))
#print(opt)
Expand All @@ -346,200 +347,8 @@ opt <- nlminb(start=obj$par, objective=obj$fn, gradient=obj$gr,
#max(abs(obj$gr()))
sdr <- TMB::sdreport(obj)
sdr_fixed <- summary(sdr, "fixed")
report <- obj$report()
#sdr <- TMB::sdreport(obj)
#sdr_fixed <- summary(sdr, "fixed")
#report <- obj$report()
### Plotting
mycols <- c("FIMS" = "blue", "ASAP" = "red", "ASAP_orig" = "darkgreen")
for (i in 1:rdat$parms$nindices){
index_results <- data.frame(
survey = i,
year = years,
observed = rdat$index.obs[[i]],
FIMS = report$exp_index[[rdat$parms$nfleet+i]],
ASAP = rdat$index.pred[[i]]
)
if (i==1){
allinds_results <- index_results
}else{
allinds_results <- rbind(allinds_results, index_results)
}
}
#print(allinds_results)
comp_index <- ggplot(allinds_results, aes(x = year, y = observed)) +
geom_point() +
geom_line(aes(x = year, y = FIMS), color = "blue") +
geom_line(aes(x = year, y = ASAP), color = "red") +
facet_wrap(~survey, scales = "free_y", nrow = 2) +
xlab("Year") +
ylab("Index") +
ggtitle("Blue=FIMS, Red=ASAP") +
theme_bw()
#print(comp_index)
catch_results <- data.frame(
observed = fishing_fleet_index$index_data,
FIMS = report$exp_index[[1]],
ASAP = as.numeric(rdat$catch.pred[1,])
)
#print(catch_results)
comp_catch <- ggplot(catch_results, aes(x = years, y = observed)) +
geom_point() +
xlab("Year") +
ylab("Catch (mt)") +
geom_line(aes(x = years, y = FIMS), color = "blue") +
geom_line(aes(x = years, y = ASAP), color = "red") +
ggtitle("Blue=FIMS, Red=ASAP") +
theme_bw()
#print(comp_catch)
pop_results <- data.frame(
Year = c(years, max(years)+1, years, years, years, years, max(years)+1, years),
Metric = c(rep("SSB", 2*nyears+1), rep("F_mort", 2*nyears), rep("Recruitment", 2*nyears+1)),
Model = c(rep("FIMS", nyears+1), rep("ASAP", nyears), rep(c("FIMS", "ASAP"), each=nyears),
rep("FIMS", nyears+1), rep("ASAP", nyears)),
Value = c(report$ssb[[1]], rdat$SSB, report$F_mort[[1]], rdat$F.report, report$recruitment[[1]], as.numeric(rdat$N.age[,1]))
)
#print(pop_results)
# ggplot(filter(pop_results, Year <=2019), aes(x=Year, y=Value, color=Model)) +
# geom_line() +
# facet_wrap(~Metric, ncol=1, scales = "free_y") +
# theme_bw() +
# scale_color_manual(values = mycols)
orig_years <- seq(orig$parms$styr, orig$parms$endyr)
orig_pop_results <- data.frame(
Year = rep(orig_years, 3),
Metric = rep(c("SSB", "F_mort", "Recruitment"), each = length(orig_years)),
Model = "ASAP_orig",
Value = c(orig$SSB, orig$F.report, as.numeric(orig$N.age[,1]))
)
pop_results_3 <- rbind(pop_results, orig_pop_results)
#print(pop_results_3)
# ggplot(filter(pop_results_3, Year <=2019), aes(x=Year, y=Value, color=Model)) +
# geom_line() +
# facet_wrap(~Metric, ncol=1, scales = "free_y") +
# theme_bw() +
# scale_color_manual(values = mycols)
comp_FRSSB3 <- ggplot(pop_results_3, aes(x=Year, y=Value, color=Model)) +
geom_line() +
facet_wrap(~Metric, ncol=1, scales = "free_y") +
theme_bw() +
scale_color_manual(values = mycols)
#print(comp_FRSSB3)
FIMS_naa_results <- data.frame(
Year = rep(c(years, max(years)+1), each = nages),
Age = rep(ages, nyears+1),
Metric = "NAA",
Model = "FIMS",
Value = report$naa[[1]]
)
ASAP_naa_results <- data.frame(
Year = rep(years, each = nages),
Age = rep(ages, nyears),
Metric = "NAA",
Model = "ASAP",
Value = as.numeric(t(rdat$N.age))
)
orig_naa_results <- data.frame(
Year = rep(orig_years, each = nages),
Age = rep(ages, length(orig_years)),
Metric = "NAA",
Model = "ASAP_orig",
Value = as.numeric(t(orig$N.age))
)
naa_results <- rbind(FIMS_naa_results, ASAP_naa_results, orig_naa_results)
#print(naa_results)
# ggplot(filter(naa_results, Year <= 2019), aes(x=Year, y=Value, color=Model)) +
# geom_line() +
# facet_wrap(~Age, ncol=1, scales = "free_y") +
# ylab("NAA") +
# theme_bw() +
# scale_color_manual(values = mycols)
comp_naa2 <- ggplot(filter(naa_results, Year <= 2019, Model %in% c("ASAP", "FIMS")), aes(x=Year, y=Value, color=Model)) +
geom_line() +
facet_wrap(~Age, ncol=1, scales = "free_y") +
ylab("NAA") +
theme_bw() +
scale_color_manual(values = mycols)
#print(comp_naa2)
# ggplot(filter(naa_results, Year == 1973, Model %in% c("ASAP", "FIMS")), aes(x=Age, y=Value, color=Model)) +
# geom_line() +
# ylab("NAA in Year 1") +
# theme_bw() +
# scale_color_manual(values = mycols)
saveplots <- TRUE
if(saveplots){
ggsave(filename = "figures/NEFSC_YT_compare_index.png", plot = comp_index, width = 4, height = 4, units = "in")
ggsave(filename = "figures/NEFSC_YT_compare_catch.png", plot = comp_catch, width = 4, height = 4, units = "in")
ggsave(filename = "figures/NEFSC_YT_compare_FRSSB3.png", plot = comp_FRSSB3, width = 5, height = 6.5, units = "in")
ggsave(filename = "figures/NEFSC_YT_compare_NAA2.png", plot = comp_naa2, width = 5, height = 6.5, units = "in")
}
```
## Comparison figures
![Catch](figures/NEFSC_YT_compare_catch.png){width=4in}
![Indices](figures/NEFSC_YT_compare_index.png){width=4in}
![FRSSB](figures/NEFSC_YT_compare_FRSSB3.png){width=4in}
![NAA](figures/NEFSC_YT_compare_NAA2.png){width=4in}

## Comparison table

The likelihood components from FIMS and ASAP for the same data are shown in the table below. Note that the ASAP file had to turn on the use likelihood constants option to enable this comparison (this option should not be used when recruitment deviations are estimated).

```{r}
jnlltab <- data.frame(Component=c("Total","Index","Age Comp", "Rec"),
FIMS = c(report$jnll, report$index_nll, report$age_comp_nll, report$rec_nll),
ASAP = c(rdat$like$lk.total,
(rdat$like$lk.catch.total + rdat$like$lk.index.fit.total),
(rdat$like$lk.catch.age.comp + rdat$like$lk.index.age.comp),
rdat$like$lk.Recruit.devs))
print(jnlltab)
```

## What was your experience using FIMS? What could we do to improve usability?

Relatively easy to use by following the vignette. Creating wrappers for data input would help so that each element did not need to be assigned directly.

## List any issues that you ran into or found

Please [open an issue](https://github.com/NOAA-FIMS/FIMS/issues/new/choose) if you found something new.

* SSB calculations in FIMS assume 0.5 multiplier, which differs from ASAP [Issue #521](https://github.com/NOAA-FIMS/FIMS/issues/521).
* Output all derived values (this is mostly done)
* Fix recruitment estimation [Issue #364](https://github.com/NOAA-FIMS/FIMS/issues/364)
* Handle missing data, especially surveys [Issue #502](https://github.com/NOAA-FIMS/FIMS/issues/502)
* Weights at age that change over time
* Separate weights at age for catch, SSB, Jan-1 population, indices, etc.
* Fishery selectivity blocks or random effects
* Allow time-varying CVs and ESS (or alternative functions)
* Option for Index in numbers
* Timing of Index and SSB calculations within the year
* One-step-ahead residuals
* Reference points, projections, pushbutton retro

## What features are most important to add based on this case study?

* Missing values, would allow inclusion of the other 3 indices (too many missing years to fill for this example)

```{r}
# Clear C++ objects from memory
clear()
```
7 changes: 2 additions & 5 deletions content/AFSC-GOA-pollock.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,9 @@ installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages], repos = "http://cran.us.r-project.org")
}
<<<<<<< HEAD
remotes::install_github("NOAA-FIMS/FIMS")
=======
#remotes::install_github("NOAA-FIMS/FIMS")
#remotes::install_github("NOAA-FIMS/FIMS", ref='isNA')
>>>>>>> 84cf899 (Jim added an initial atka case)
remotes::install_github("kaskr/TMB_contrib_R/TMBhelper")
#remotes::install_github("kaskr/TMB_contrib_R/TMBhelper")
library(ggplot2)
library(tidyr)
Expand Down

0 comments on commit b4901b1

Please sign in to comment.