In PaulRegular/MSP: Multispecies production model incorporating covariance

Results {-}

Both cross-validation metrics (LOO-CV and Hindcast-CV scores) indicate that most multispecies production model formulations outperform a single-species production model when applied to seven species within three ecosystem production units (Northeast NL Shelf, Grand Bank, and Southern NL) off the east coast of Canada (Figure \@ref(fig:scores)). Focusing on overall scores, the performance of the single-species production model was similar to a multispecies formulation that assumes there is no correlation in the process errors across species or time. There tends to be more notable decreases in the scores as temporal and species correlations are introduced, indicating an improvement in the predictive ability of these models. The "shared correlation" and "just correlation" formulations, in particular, tended to receive the lowest scores, and dropping the species and temporal correlations in lieu of a shift covariate ("just shift" formulation) resulted in a deterioration of predictive ability. Scores were improved when temporal and species-to-species correlations were introduced along with the shift covariate ("full" formulation); however, the fit of the "full" model tended to be poorer than the "just correlation" formulation, which further indicates that the "shift" covariate degraded the predictive ability of the model. Subsequent plots focus on the best fitting formulation, "just correlation", to demonstrate model predictions.

The multispecies production model with unstructured species-to-species correlation and AR1 temporal correlation ("just correlation") offered an explanation of the trends in survey indices of focal species across three ecosystem production units with little signs of systematic bias (see residual plots included in model dashboards; Supplement 1). Predicted indices track observed values and, by estimating catchability parameters by species and survey, indices from temporally fragmented surveys are stitched together and their trends are used to inform a continuous underlying trend in biomass (Figure \@ref(fig:survey-trends)). The earlier Yankee and Engel eras of the Canadian surveys tended to receive lower catchability estimates than the Campelen era survey; indices since 1996 therefore tend to be closer, in relative terms, to the underlying estimates of biomass from the model.

Isolating residual changes in biomass not explained by reported fisheries landings or the production function (equation \@ref(eq:ms-logistic)) reveals substantive subtractions in the early 1990s across all three ecosystem production units (Figure \@ref(fig:pop-trends)). These residual subtractions exceed the absolute scale of landings taken in the late 1980s from the Northeast NL Shelf and Grand Bank ecosystem production units. The scale of changes derived from the production function and other processes are of a scale similar to landings reported in the Southern NL ecosystem production unit. Estimates of biomass exceeded the carrying capacity estimated for the Northeast NL Shelf through the 1980s, and all focal species displayed abrupt declines in the early 1990s. Comparing community composition in the early 1980s to the 2010s, there are no clear shifts in the relative biomass of the focal species on the Northeast NL Shelf. Elsewhere, total biomass exceeded or approached system carrying capacity in the late 1980s, after which all species declined. Estimates of total biomass have gradually increased on the Grand Bank and off Southern NL since the mid 1990s, largely due to increasing Redfish spp. biomass estimates.

regions <- c("2J3K", "3LNO", "3Ps")
rhos <- vector("list", 3)
names(rhos) <- regions
for (r in regions) {
    fits <- readRDS(here(paste0("analysis/NL_case_study/exports/fits_", r, ".rds")))
    rho <- c(multispic::inv_logit(fits$shared_cor$par$logit_rho),
             multispic::inv_logit(fits$shared_cor$par_lwr$logit_rho),
             multispic::inv_logit(fits$shared_cor$par_upr$logit_rho))
    rhos[[r]] <- paste0(formatC(rho[1], digits = 2, format = "f"), " (95% CI: ", 
                  formatC(rho[2], digits = 2, format = "f"), ", ", 
                  formatC(rho[3], digits = 2, format = "f"), ")")
}

Further inspection of the process errors reveals common patterns across all focal species across three ecosystem production units (Figure \@ref(fig:pe-cor)). Like the exponentiated and unstandardized process errors ("residual change") shown in Figure \@ref(fig:pop-trends), the standardized process errors highlight substantive subtractions in the early 1990s, representing time-series lows for 19 out of 21 populations. The standardized values also reveal parallel and periodic increases and decreases within each region. The species-to-species correlations in the process errors estimated by the "just correlation" model indicate that 84% (53 of 63) of pairs were positively correlated. Also note that the "shared correlation" model estimates of a common species-to-species correlation parameter were r rhos[["2J3K"]], r rhos[["3LNO"]], and r rhos[["3Ps"]] for the Northeast NL Shelf, Grand Bank, and Southern NL ecosystem production units, respectively. This model received similar cross-validation scores as the "just correlation" model which, taken together, further supports the inference that the process errors are primarily positively correlated across species within each production unit.

PaulRegular/MSP documentation built on Dec. 16, 2024, 1:59 p.m.