Population Estimates

Andrew Edwards

Last rendered on 12 May, 2026

Estimates of animal populations

library(pacea)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following object is masked from 'package:testthat':
#> 
#>     matches
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

In pacea we include estimates of populations, including uncertainty, as calculated from recent stock assessments.

Currently, we have Pacific Hake and Pacific Herring, major components of the ecosystem, and Pacific Harbour Seals, the most abundant pinniped species in the Northeast Pacific. Click on links to jump to the relevant sections below.

NEW (since original release of pacea): we have updated the Pacific Hake estimates with those from the 2024 assessment, but retained the 2023 assessment’s estimates within the package (as hake_biomass_2023 etc.). See below and the help files for more details. Also added total biomass of age-1 hake and recruitment deviations. And Pacific Herring was added in May 2024.

Pacific Hake

Pacific Hake biomass

The hake biomass time series is saved as a tibble in pacea, as simply:

hake_biomass
#> # A tibble: 60 × 4
#>     year   low median  high
#>    <dbl> <dbl>  <dbl> <dbl>
#>  1  1966 0.560  0.947  1.85
#>  2  1967 0.600  0.961  1.88
#>  3  1968 0.604  0.972  1.94
#>  4  1969 0.715  1.11   2.23
#>  5  1970 0.818  1.27   2.60
#>  6  1971 0.826  1.31   2.68
#>  7  1972 0.929  1.49   3.09
#>  8  1973 1.10   1.76   3.60
#>  9  1974 1.07   1.70   3.43
#> 10  1975 0.908  1.47   2.95
#> # ℹ 50 more rows
tail(hake_biomass)
#> # A tibble: 6 × 4
#>    year   low median  high
#>   <dbl> <dbl>  <dbl> <dbl>
#> 1  2020 0.845  1.17   1.99
#> 2  2021 0.629  0.933  1.71
#> 3  2022 0.548  0.918  1.85
#> 4  2023 0.560  1.11   2.55
#> 5  2024 0.530  1.19   2.90
#> 6  2025 0.521  1.22   3.03

with the annual estimates given as median, low and high ends of the 95% credible intervals. Units are millions of tonnes of female spawning biomass at the start of the year, as documented in ?hake_biomass which also refers to ?hake_recruitment for full details and further background. For example, it is important to realise that the assessment considers a coastwide stock off the west coast of Canada and the United States.

To plot the time series of estimated coastwide biomass, simply do

plot(hake_biomass)

where the solid circles are the annual medians, with credible intervals given by the shaded regions (mostly reproducing Figure d in the latest assessment).

Note that the hake_biomass object also has class pacea_biomass:

class(hake_biomass)
#> [1] "pacea_biomass" "tbl_df"        "tbl"           "data.frame"

such that the above command plot(hake_biomass) first looks to use our tailored function plot.pacea_biomass(), which gives the style of plot shown above. We use this approach for hake_recruitment and various other objects in pacea, to enable easy default plotting. All functions have help files, for example ?plot.pacea_biomass.

We now also include the total biomass of age-1 hake, due to its potential influence on hake recruitment Vestfals et al. (2023). Note this is total female and male biomass (not spawning only like hake_biomass), and credible intervals are not current easily available.

hake_total_biomass_age_1
#> # A tibble: 60 × 2
#>     year median
#>    <dbl>  <dbl>
#>  1  1966  140. 
#>  2  1967  122. 
#>  3  1968  362. 
#>  4  1969  232. 
#>  5  1970   53.8
#>  6  1971  693. 
#>  7  1972   62.5
#>  8  1973   40.6
#>  9  1974  454. 
#> 10  1975   16.6
#> # ℹ 50 more rows
tail(hake_total_biomass_age_1)
#> # A tibble: 6 × 2
#>    year median
#>   <dbl>  <dbl>
#> 1  2020  15.7 
#> 2  2021 320.  
#> 3  2022 561.  
#> 4  2023   9.13
#> 5  2024  70.2 
#> 6  2025  70.5

for which units are thousands of tonnes, and our default plot is

plot(hake_total_biomass_age_1)

Pacific Hake recruitment

The hake_recruitment object contains annual estimates of recruitment of age-0 fish

hake_recruitment
#> # A tibble: 60 × 4
#>     year    low median  high
#>    <dbl>  <dbl>  <dbl> <dbl>
#>  1  1966 0.0521  1.65   9.95
#>  2  1967 0.223   4.89  15.2 
#>  3  1968 0.255   3.15   9.89
#>  4  1969 0.0455  0.731  4.16
#>  5  1970 4.77    9.43  23.0 
#>  6  1971 0.0726  0.845  3.08
#>  7  1972 0.0629  0.552  1.94
#>  8  1973 3.23    6.19  14.5 
#>  9  1974 0.0367  0.345  1.37
#> 10  1975 0.920   1.89   4.41
#> # ℹ 50 more rows
tail(hake_recruitment)
#> # A tibble: 6 × 4
#>    year    low median  high
#>   <dbl>  <dbl>  <dbl> <dbl>
#> 1  2020 1.59    3.40   9.14
#> 2  2021 2.89    7.06  19.3 
#> 3  2022 0.0141  0.131  1.01
#> 4  2023 0.474   0.879  1.85
#> 5  2024 0.476   0.882  1.87
#> 6  2025 0.476   0.884  1.88

for which values are billions of age-0 fish, and our default plot is

plot(hake_recruitment)

The dots represent the medians, with blue bars showing the 95% credible intervals. This reproduces Figure f in the latest stock assessment (without some extra details).

Note that, as mentioned in ?hake_recruitment the most recent few years are very uncertain as they are not informed by data. In particular, the final year represents recruitment at the start of the year but has no data to inform it, and so should definitely not be used for any analyses (we include it as it is shown in the stock assessment; it also essentially represents the prior distribution used for recruitment).

The estimated uncertainty of recruitments is large, and so we also include scaled values that make it easier to compare years. The 2010 recruitment is the second largest in the time series (below only 1980’s) and as such people have an intuition that it is large, and we can compare other years by looking at recruitment scaled by the 2010 recruitment:

hake_recruitment_over_2010
#> # A tibble: 60 × 4
#>     year     low median   high
#>    <dbl>   <dbl>  <dbl>  <dbl>
#>  1  1966 0.00337 0.110  0.543 
#>  2  1967 0.0154  0.336  0.738 
#>  3  1968 0.0166  0.214  0.508 
#>  4  1969 0.00296 0.0485 0.232 
#>  5  1970 0.407   0.626  1.02  
#>  6  1971 0.00510 0.0570 0.159 
#>  7  1972 0.00414 0.0373 0.102 
#>  8  1973 0.279   0.411  0.646 
#>  9  1974 0.00262 0.0231 0.0730
#> 10  1975 0.0757  0.126  0.207 
#> # ℹ 50 more rows
tail(hake_recruitment_over_2010)
#> # A tibble: 6 × 4
#>    year     low  median   high
#>   <dbl>   <dbl>   <dbl>  <dbl>
#> 1  2020 0.117   0.226   0.474 
#> 2  2021 0.228   0.466   1.01  
#> 3  2022 0.00100 0.00854 0.0616
#> 4  2023 0.0389  0.0576  0.0897
#> 5  2024 0.0391  0.0579  0.0905
#> 6  2025 0.0390  0.0581  0.0909
plot(hake_recruitment_over_2010)

This shows the medians and 95% credible intervals (red lines) of the recruitment divided by that in 2010, which clearly shows that, for example, the 2014 is clearly smaller than the 2010 recruitment (highlighted by the green dashed line), which is not intuitive in the above plot(hake_recruitment) plot (see ?hake_recruitment_over_2010 for more details and motivation).

This approach somewhat scales out some of the uncertainty in the absolute scale of \(R_0\), the mean unfished equilibrium recruitment. This is similarly seen by instead scaling recruitments by \(R_0\), as

hake_recruitment_over_R0
#> # A tibble: 60 × 4
#>     year    low median  high
#>    <dbl>  <dbl>  <dbl> <dbl>
#>  1  1966 0.0212  0.660 3.61 
#>  2  1967 0.0964  2.02  4.93 
#>  3  1968 0.103   1.30  3.40 
#>  4  1969 0.0184  0.299 1.46 
#>  5  1970 2.15    3.90  6.93 
#>  6  1971 0.0310  0.345 1.03 
#>  7  1972 0.0270  0.227 0.652
#>  8  1973 1.49    2.56  4.31 
#>  9  1974 0.0164  0.143 0.461
#> 10  1975 0.415   0.773 1.39 
#> # ℹ 50 more rows
tail(hake_recruitment_over_R0)
#> # A tibble: 6 × 4
#>    year     low median  high
#>   <dbl>   <dbl>  <dbl> <dbl>
#> 1  2020 0.621   1.39   3.26 
#> 2  2021 1.22    2.88   7.07 
#> 3  2022 0.00608 0.0530 0.364
#> 4  2023 0.302   0.364  0.382
#> 5  2024 0.303   0.366  0.385
#> 6  2025 0.303   0.367  0.386
plot(hake_recruitment_over_R0)

Also now available are the log-scale recruitment deviations from the stock-recruitment curve (see ?hake_recruitment_deviations), which are sometimes used to look for environmental drivers, e.g. Vestfals et al. (2023):

hake_recruitment_deviations
#> # A tibble: 57 × 4
#>     year    low  median  high
#>    <dbl>  <dbl>   <dbl> <dbl>
#>  1  1966 -2.81   0.637  2.30 
#>  2  1967 -1.30   1.75   2.61 
#>  3  1968 -1.23   1.31   2.23 
#>  4  1969 -2.98  -0.182  1.39 
#>  5  1970  1.82   2.37   2.91 
#>  6  1971 -2.45  -0.0493 1.02 
#>  7  1972 -2.61  -0.484  0.541
#>  8  1973  1.43   1.92   2.41 
#>  9  1974 -3.11  -0.960  0.195
#> 10  1975  0.161  0.740  1.29 
#> # ℹ 47 more rows
tail(hake_recruitment_deviations)
#> # A tibble: 6 × 4
#>    year    low median    high
#>   <dbl>  <dbl>  <dbl>   <dbl>
#> 1  2017 -0.217  0.388  0.984 
#> 2  2018 -3.25  -1.75  -0.670 
#> 3  2019 -3.11  -1.56  -0.429 
#> 4  2020  0.576  1.35   2.17  
#> 5  2021  1.29   2.10   2.97  
#> 6  2022 -4.04  -1.89   0.0314
plot(hake_recruitment_deviations)

Updates due to new stock assessments

Note that the Pacific Hake stock assessment is conducted annually at the start of each year, and these estimates will be updated in pacea by April each year. So estimates will change (particularly for the most recent years) as more data inform the stock assessment model. So hake_biomass, hake_recruitment etc. will be based on the most recent stock assessment (see ?hake_recruitment).

To ensure that estimates from previous assessments are still available in pacea (so you can ensure fully repeatable analyses that do not change when pacea is updated with the new estimates), we also save each assessment’s estimates with the year appended; for example, hake_biomass_2023 and hake_biomass_2024 for the spawning biomass estimates from the 2023 and 2024 assessments. Similarly, we also save hake_recruitment_2023, hake_recruitment_2024, hake_recruitment_over_2010_2023, hake_recruitment_over_2010_2024, hake_recruitment_over_R0_2023, and hake_recruitment_over_R0_2024, plus similar objects and later years once they become available. See ?hake_recruitment for details and references. The above examples will work on all these objects, for example, the values and plot shown in the previous code chunk, but based on the 2023 assessment results, are:

hake_recruitment_over_R0_2023
#> # A tibble: 58 × 4
#>     year    low median  high
#>    <dbl>  <dbl>  <dbl> <dbl>
#>  1  1966 0.0259  0.604 3.33 
#>  2  1967 0.103   1.81  4.50 
#>  3  1968 0.0923  1.17  3.05 
#>  4  1969 0.0171  0.265 1.34 
#>  5  1970 2.04    3.55  6.40 
#>  6  1971 0.0290  0.312 0.964
#>  7  1972 0.0247  0.208 0.613
#>  8  1973 1.40    2.37  4.08 
#>  9  1974 0.0177  0.137 0.438
#> 10  1975 0.389   0.728 1.32 
#> # ℹ 48 more rows
tail(hake_recruitment_over_R0_2023)
#> # A tibble: 6 × 4
#>    year    low median   high
#>   <dbl>  <dbl>  <dbl>  <dbl>
#> 1  2018 0.0735  0.249  0.719
#> 2  2019 0.0449  0.238  0.779
#> 3  2020 1.19    4.45  17.2  
#> 4  2021 0.0124  0.176  2.48 
#> 5  2022 0.0171  0.367  7.87 
#> 6  2023 0.0185  0.368  7.15
plot(hake_recruitment_over_R0_2023)

Note how hake_recruitment_over_R0_2023 contains no estimate for the 2024 recruitment, but hake_recruitment does (although it is not actually informed by data, so should not be used in any analysis). And see how the (highly uncertain) 2020 and 2021 recruitment estimates changed somewhat from 2023 assessment to the 2024 assessment, due to updated data:

filter(hake_recruitment_2023, year %in% c(2020, 2021))
#> # A tibble: 2 × 4
#>    year    low median  high
#>   <dbl>  <dbl>  <dbl> <dbl>
#> 1  2020 2.91   11.4   47.6 
#> 2  2021 0.0283  0.450  6.91
filter(hake_recruitment_2024, year %in% c(2020, 2021))
#> # A tibble: 2 × 4
#>    year   low median  high
#>   <dbl> <dbl>  <dbl> <dbl>
#> 1  2020  2.06   4.75  12.7
#> 2  2021  4.09  10.2   29.5

See page 13 of the 2024 assessment document for explanation.

Pacific Herring

Pacific Herring results from the latest assessment (DFO 2024; SR 2024/001) are saved in a similar manner to those for Pacific Hake. The spawning biomass time series is saved as a tibble in pacea, as simply:

herring_spawning_biomass
#> # A tibble: 365 × 5
#>     year region   low median  high
#>    <dbl> <fct>  <dbl>  <dbl> <dbl>
#>  1  1951 HG      7.20   9.86 14.4 
#>  2  1952 HG      3.71   6.55 11.5 
#>  3  1953 HG     15.8   20.7  28.6 
#>  4  1954 HG     42.7   51.9  65.5 
#>  5  1955 HG     55.9   63.8  72.9 
#>  6  1956 HG     11.7   14.7  18.0 
#>  7  1957 HG      3.57   4.93  6.60
#>  8  1958 HG      2.55   3.80  5.93
#>  9  1959 HG      4.37   7.75 13.3 
#> 10  1960 HG      8.45  12.8  19.4 
#> # ℹ 355 more rows
tail(herring_spawning_biomass)
#> # A tibble: 6 × 5
#>    year region   low median  high
#>   <dbl> <fct>  <dbl>  <dbl> <dbl>
#> 1  2018 WCVI    11.8   15.3  19.6
#> 2  2019 WCVI    11.7   15.4  20.2
#> 3  2020 WCVI    13.7   18.1  24.2
#> 4  2021 WCVI    17.8   24.0  32.8
#> 5  2022 WCVI    24.2   34.7  49.4
#> 6  2023 WCVI    24.4   41.2  67.5

with the annual estimates given as median, low and high ends of the 90% credible intervals. Units are thousands of tonnes of female spawning biomass each year, as documented in ?herring_spawning_biomass which also refers to ?herring_recruitment for full details and further background. For example, it is important to know that these are 90% credible intervals (not 95% like for hake), units are thousands of tonnes (not millions like for hake), and are for males and females combined (not females only like for hake). As usual, do read the help files. (We included spawning in the name herring_spawning_biomass as we may also include total biomass at some point if available).

Results are included for the five major herring stock assessment regions: Haida Gwaii (HG), Prince Rupert District (PRD), Central Coast (CC), Strait of Georgia (SOG), and West Coast of Vancouver Island (WCVI). These are given in the region column, and to see values for just one region it is easy to filter:

filter(herring_spawning_biomass,
       region == "HG")
#> # A tibble: 73 × 5
#>     year region   low median  high
#>    <dbl> <fct>  <dbl>  <dbl> <dbl>
#>  1  1951 HG      7.20   9.86 14.4 
#>  2  1952 HG      3.71   6.55 11.5 
#>  3  1953 HG     15.8   20.7  28.6 
#>  4  1954 HG     42.7   51.9  65.5 
#>  5  1955 HG     55.9   63.8  72.9 
#>  6  1956 HG     11.7   14.7  18.0 
#>  7  1957 HG      3.57   4.93  6.60
#>  8  1958 HG      2.55   3.80  5.93
#>  9  1959 HG      4.37   7.75 13.3 
#> 10  1960 HG      8.45  12.8  19.4 
#> # ℹ 63 more rows

Built-in plotting functions have been developed, with the default being to plot all regions:

plot(herring_spawning_biomass)

where the solid circles are the annual medians, with credible intervals given by the shaded regions, matching the style used above for hake.

We have made it easy to plot results for just one region:

plot(herring_spawning_biomass,
     region = "PRD")

See ?plot.pacea_biomass_herring for options, such as title = "short" to use the region’s acronym rather than full name.

The herring_recruitment object contains annual estimates of recruitment of age-2 fish

herring_recruitment
#> # A tibble: 355 × 5
#>     year region    low median  high
#>    <dbl> <fct>   <dbl>  <dbl> <dbl>
#>  1  1953 HG     0.683  1.04   1.70 
#>  2  1954 HG     0.0957 0.165  0.286
#>  3  1955 HG     0.0823 0.148  0.267
#>  4  1956 HG     0.108  0.176  0.267
#>  5  1957 HG     0.0885 0.138  0.211
#>  6  1958 HG     0.171  0.265  0.425
#>  7  1959 HG     0.0454 0.0888 0.176
#>  8  1960 HG     0.142  0.266  0.516
#>  9  1961 HG     0.240  0.418  0.741
#> 10  1962 HG     0.213  0.371  0.623
#> # ℹ 345 more rows
tail(herring_recruitment)
#> # A tibble: 6 × 5
#>    year region   low median  high
#>   <dbl> <fct>  <dbl>  <dbl> <dbl>
#> 1  2018 WCVI   0.313  0.429 0.595
#> 2  2019 WCVI   0.212  0.297 0.409
#> 3  2020 WCVI   0.578  0.803 1.11 
#> 4  2021 WCVI   0.480  0.677 0.977
#> 5  2022 WCVI   0.585  0.867 1.27 
#> 6  2023 WCVI   0.487  0.815 1.36

for which values are billions of age-2 fish, and our default plot is

plot(herring_recruitment)

This uses the same style of plot, with the dots represent the medians, with blue bars showing the 90% credible intervals. Again, it is easy (see ?plot.pacea_recruitment_herring) to plot just one of the regions:

plot(herring_recruitment,
     region = "SOG")

We aim to update herring_spawning_biomass and herring_recruitment annually as new assessments are conducted, and so have also saved the 2023 assessment’s results as herring_spawning_biomass_2023 and herring_recruitment_2023, which will not change going forward as pacea is updated.

Pacific Harbour Seals

Estimated abundances for seven regions (and coastwide) were calculated by DFO (2022; SAR 2022/034), and are included in pacea as:

harbour_seals
#> # A tibble: 640 × 5
#>    date       region   low  mean  high
#>    <date>     <fct>  <dbl> <dbl> <dbl>
#>  1 1965-01-01 SOG     593. 1384. 3231.
#>  2 1965-09-07 SOG     669. 1488. 3307.
#>  3 1966-05-14 SOG     755. 1599. 3387.
#>  4 1967-01-19 SOG     851. 1719. 3471.
#>  5 1967-09-25 SOG     959. 1848. 3562.
#>  6 1968-06-01 SOG    1080. 1988. 3661.
#>  7 1969-02-06 SOG    1214. 2140. 3772.
#>  8 1969-10-14 SOG    1364. 2305. 3895.
#>  9 1970-06-20 SOG    1531. 2485. 4036.
#> 10 1971-02-25 SOG    1715. 2683. 4195.
#> # ℹ 630 more rows

tail(harbour_seals)
#> # A tibble: 6 × 5
#>   date       region       low   mean    high
#>   <date>     <fct>      <dbl>  <dbl>   <dbl>
#> 1 2015-08-01 Coastwide 75072. 92401. 118166.
#> 2 2016-04-07 Coastwide 73146. 91001. 117663.
#> 3 2016-12-13 Coastwide 71189. 89638. 117317.
#> 4 2017-08-20 Coastwide 69222. 88313. 117113.
#> 5 2018-04-26 Coastwide 67259. 87026. 117031.
#> 6 2019-01-01 Coastwide 65316. 85774. 117053.

Estimates are for each region and for the coastwide population, reproducing Figure 3 of DFO (2022):

plot(harbour_seals)

Or just plot estimates for a single region:

plot(harbour_seals,
     region = "SOG")

For full details and reference see

?harbour_seals

We anticipate adding assessment output for other species soon, so stay tuned (check the NEWS). For zooplankton anomalies see the new zooplankton vignette: zooplankton.html.