Stochastic population forecasting

The vital package can be used for stochastic population forecasting with coherent components. This is based on the papers by Hyndman and Booth (2008) and Hyndman et al. (2013). Following Hyndman and Booth (2008), we use the following demographic growth-balance equations for each sex: \[ \begin{aligned} P_{t+1}(x+1) &= P_t(x) - D_t(x,x+1) + G_t(x,x+1), \qquad x = 0,1,\dots,p-2,\\ P_{t+1}(p+) &= P_t(p-1) + P_t(p+) - D_t(p-1+,p+) + G_t(p-1+,p+),\\ P_{t+1}(0) &= B_t - D_t(B,0) + G_t(B,0), \end{aligned} \] where

These are cohort rhomboids in the Lexis diagram. The cohort deaths are estimated from the survivorship ratios of the life table: \[ \begin{aligned} D_t(x,x+1) &= P_t(x)[1 - L(x+1)/L(x)], \qquad x = 0,1,\dots,p-2,\\ D_t(p-1+,p+) &= [P_t(p-1) + P_t(p+)][1 - T(p)/T(p-1)],\\ D_t(B,0) &= B_t[1 - L(0)/l(0)], \end{aligned} \] where \(L(x)\) denotes person-years lived at age \(x\), \(T(x)\) denotes person-years lived at age \(x\) and older, and \(l(0)\) is the radix of the life table. Net migration is then estimated as the residual from the growth-balance equations. The net_migration() function returns these estimates indexed by the age of each cohort at the end of the year, so that \(G_t(x,x+1)\) is at age \(x+1\), \(G_t(p-1+,p+)\) is at age \(p\), and \(G_t(B,0)\) is at age 0.

In the simulation of future populations, half of each cohort’s net migrants are added at the beginning of the year and half at the end. See Hyndman and Booth (2008) for details.

To simulate future births, deaths and net migrants, we develop three functional data models for fertility, mortality and migration. The models for mortality and migration use coherent components, so that the rates for males and females do not diverge over time.

The following example uses Norwegian mortality data up to 2023 and fertility data up to 2022, and produces simulated populations for the ten future years. Different models are used for each component to demonstrate the flexibility of the package, but other models can be used as well.

Mortality model

We use a coherent functional data model (Hyndman et al. 2013) for the log mortality rates. To ensure coherence, we compute the geometric mean of the sex-specific mortality rates and the corresponding ratios, using the make_pr() function. The data are smoothed first, and 6 components are used by default in each FDM, although only the first two are plotted here.

fit_mortality <- norway_mortality |>
  filter(Sex != "Total") |>
  smooth_mortality(Mortality) |>
  make_pr(.smooth) |>
  model(fdm = FDM(log(.smooth), coherent = TRUE))
autoplot(fit_mortality, 2)

Fertility model

For fertility, we use a functional mean model with a square root transformation, applied to the last 12 years of data. The plotted model shows the fitted values on the square root scale.

fit_fertility <- norway_fertility |>
  filter(Year > 2010) |>
  smooth_fertility(Fertility) |>
  model(fmean = FMEAN(sqrt(.smooth)))
autoplot(fit_fertility)

Net migration model

For net migration, we use a coherent functional data model. Because net migration values can be positive or negative, we can’t take products and ratios. Instead, we need to compute the means and corresponding differences using the make_sd() function.

netmig <- net_migration(
  norway_mortality |> filter(Sex != "Total"),
  norway_births
) |>
  make_sd(NetMigration)
fit_migration <- netmig |>
  model(fdm = FDM(NetMigration, coherent = TRUE))
autoplot(fit_migration)

Population simulation

The generate_population() function takes a starting population, and the three component models, and simulates future age-sex-specific population values. Here we produce 500 replicates of the future population.

pop <- norway_mortality |>
  filter(Sex != "Total", Year == max(Year))
future <- generate_population(
  starting_population = pop,
  mortality_model = fit_mortality,
  fertility_model = fit_fertility,
  migration_model = fit_migration,
  h = 10,
  n_reps = 500
)

One replicate is plotted below, along with the last few years of historical data.

future |>
  filter(.rep == "100") |>
  ggplot(aes(x = Age, y = Population, group = Year, color = Year)) +
  geom_line(
    data = norway_mortality |> filter(Year > 2010, Sex != "Total"),
    color = "grey",
    mapping = aes(group = Year)
  )  +
  geom_line() +
  scale_color_gradientn(colours = rainbow(10)[1:9]) +
  facet_grid(. ~ Sex)

The simulated populations can be used to compute any quantities that can be derived from populations numbers by sex and age. For example, the mean age of the population for the next 10 years

future |>
  group_by(Sex, .rep) |>
  summarise(mean_age = sum(Population * (Age + 0.5)) / sum(Population)) |>
  group_by(Sex) |>
  summarise(mean_age = mean(mean_age))

We can also plot population pyramids with prediction intervals. For example, here is the population pyramid for 2032 with a 95% prediction interval.

pyramid_2032 <- future |>
  filter(Year == 2032) |>
  mutate(Population = if_else(Sex == "Female", -Population, Population)) |>
  group_by(Age, Sex) |>
  summarise(
    lo = quantile(Population, 0.025),
    med = quantile(Population, 0.5),
    hi = quantile(Population, 0.975)
  )
pyramid_2032 |>
  ggplot(aes(x = Age)) +
  geom_ribbon(aes(ymin = lo, ymax = hi, colour = NULL),
    fill = "#c14b14", alpha = 0.2
  ) +
  geom_line(aes(y = med), color = "#c14b14") +
  facet_grid(. ~ Sex, scales = "free_x") +
  labs(y = "Population") +
  coord_flip() +
  guides(fill = "none", alpha = "none")

Stochastic population forecasting

Rob J Hyndman

Mortality model

Fertility model

Net migration model

Population simulation

References