code.qmd

---
title: "Stock Assessment Chapter Example"
format: pdf
---

## Outline
- Introduction 
    - Case Study Escanaba Lake Walleye
    - R packages (RTMB, TMB, stockassessment?)
    - Look at data
- Surplus production model
- Virtual Population Analysis
- SCAA (FSA)
- reference points?
- Model diagnostics and uncertainty

## Introduction

The analyses use the R programming environment (citation). 
All packages used in this section for Chapter X are listed below. 

```{R eval = FALSE}
library("RTMB") # (Kristensen 2024)
```

The data used in this section are from Walleye *Sander vitreus* in Escanaba Lake, Wisconsin. Escanaba Lake is a 119-ha, drainage lake in Wisconsin. The fishes of Escanaba Lake have been the basis of long-term research to determine effects of angling on the community (Petterson 1953; Churchill 1957; Kempinger and Churchill 1972; Kempinger et al. 1975; Kempinger and Carline 1977; Serns and Sempinger 1981; Serns 1982a, 1982b, 1986a, 1986b, 1987). Walleyes were introduced to the lake via stocking from 1933-1945, and were firmly established by 1948 (). 

For this example, we use Walleye data from 1981-2002, prior to the recreational regulatory changes in 2003. 
Fyke nets (1.2 x 1.8-m frame; 19.1-mm mesh) are used for the sampling data each year immediately after ice-out to mark adult Walleye for mark-recapture population estimates. 
- Creel harvest was used as recapture gear until 2003
- change in use of recapture gear (from creel harvest to electrofishing) and catchability of those gear has a big effect on the population estimates over time, and thus we only considered the data from 1981-2002 for this chapter.
- liberalized angling regulations - no bag limits, closed season, or minimum length limits
- gear, years, life history values recorded? aging info
- 

```{R}
load("data/dat.Rdata")
str(dat)
# years
years <- unique(dat$year)
print(years)
# ages
ages <- unique(dat$age)
print(ages)
# data types
print(dat$fleet_names)
```

```{R, echo = FALSE, fig.width=7, fig.height=8, fig.cap="Figure X. "}
# Recreational catch at age (numbers)
par(mfrow=c(2,1), mar = c(3,4,2,1))
matplot(years, dat$caa, type = "b", 
        xlab = "Years", ylab = "Catch at age (numbers)",
        main = "Recreational")
text(2001.5, max(dat$caa)*0.95, "A")
par(mar = c(4,4,1,1))
plot(years, dat$ct, type = "l", ylim = c(0, max(dat$ct)), 
    xlab = "Years", ylab = "Total Catch (numbers)")
text(2001.5, max(dat$ct)*0.95, "B")

# Fishery-independent survey at age (numbers/ha)
par(mfrow=c(2,1), mar = c(3,4,2,1))
matplot(years, dat$iaa, type = "b", 
        xlab = "Years", ylab = "Survey at age (numbers/ha)",
        main = "Fishery-independent survey")
text(2001.5, max(dat$iaa)*0.95, "A")
par(mar = c(4,4,1,1))
plot(years, dat$it, type = "l", ylim = c(0, max(dat$it)),
    xlab = "Years", ylab = "Total Catch")
text(2001.5, max(dat$it)*0.95, "B")
```

## RTMB basics

Data

Parameter objects are gathered in a list (`pars`) that also serves as initial guesses when fitting the model:

```{R, eval = FALSE}
pars <- list()
pars$a <- 1
pars$b <- 2
```

Note that only parameters of interest (ones that will be estimated) should be in the `pars` list.

The objective function `f` takes the parameter list as input, and serves as the...:

```{R, eval = FALSE}
f <- function(pars) {
  getAll(data, pars)

  obs <- OBS(obs)
  
  # initialize joint negative log likelihood
  jnll <- 0

  ...

  # Export
  REPORT(y)
  # Get predicted uncertainties 
  ADREPORT()

  # Return
  return(jnll)
}
```

This function calculates the negative log likelhood (or joint negative log likelihood).

The `getAll` function makes all the list elements of data and parameters visible inside the function so that one can write ...

The `OBS()` statement tells RTMB that ... is the response, which is needed to enable automatic simulation and residual calculations.

The `REPORT()` statement tells

The `ADREPORT()` statement tells RTMB that we want uncertainties for this calculation.

The objective function `f` is processed by RTMB using this call:

```{R, eval = FALSE}
obj <- MakeADFun(f, pars)
```

We optimize the model using `nlminb` (other gradient based optimizers in R can also be used):

```{R, eval = FALSE}
opt <- nlminb(obj$par, obj$fn, obj$gr)
```

Uncertainties (point estimate + standard error) can be calculated using:
```{R, eval = FALSE}
sd_rep <- sdreport(obj)
```

## Surplus Production Model



## Virtual Population Analysis
## Statistical Catch at Age Model
## Referemce Points
## Model Diagnostics and Uncertainty
## References