---
title: "Demographic and environmental variance in an IPM"
description: "Separating demographic from environmental variance in population growth with an integrated model, why it depends on abundance, and when it is identifiable."
date: "2026-05-21 11:00"
categories: [R, population dynamics, capture-recapture, demography, ecology tutorial]
image: thumbnail.png
image-alt: "Two variance curves against population size: a declining demographic curve and a flat environmental line crossing at a threshold, with the observed abundance range shaded."
---
A population's numbers wobble from year to year for two quite different reasons. One is chance in a finite population: even with fixed vital rates, the number of survivors is a binomial draw and the number of recruits a Poisson draw, so the realised growth rate scatters around its expectation. This demographic noise shrinks as the population grows, because averages over more individuals are steadier. The other reason is the environment: good years and bad years shift survival and recruitment themselves, and this variation does not fade with population size. The two have opposite relationships with abundance, and telling them apart matters for any projection or extinction risk.
Counts alone cannot make the split, because the size of the demographic component depends on the vital rates, which the counts do not identify. An integrated population model (IPM) solves this. The capture-recapture and productivity data pin survival and recruitment, so the demographic variance is known; whatever fluctuation in the counts remains, beyond that known demographic scatter, is environmental. This post fits a state-space count model whose process variance carries both terms and recovers the environmental variance separately, then shows why abundance decides which term dominates and when the environmental part can be estimated at all.
## Two kinds of noise, two scalings
Write the realised per-capita growth rate as its expectation plus scatter. The demographic contribution to its variance is the binomial-plus-Poisson variance divided by population size: survival adds `phi(1 - phi)` and recruitment adds `f`, so the demographic variance of the growth rate is `phi(1 - phi) + f` over `N`. The environmental contribution is the variance of the rates themselves across years, call it `sigma_e` squared, and it does not depend on `N`. Setting the two equal gives a crossover abundance, `N*` equal to `phi(1 - phi) + f` divided by `sigma_e` squared: below it the population is ruled by demographic chance, above it by the environment.
```{r}
#| label: setup
#| code-fold: true
#| code-summary: "Setup: packages, colours and plot theme"
#| results: hide
#| message: false
#| warning: false
library(ggplot2)
paper <- "#f5f4ee"; ink <- "#16241d"; body <- "#2c3a31"; forest <- "#275139"
labelc <- "#46604a"; sage <- "#93a87f"; linec <- "#dad9ca"; faint <- "#5d6b61"
gold <- "#c9b458"; rust <- "#b5534e"
theme_te <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
plot.background = element_rect(fill = paper, colour = NA),
panel.background = element_rect(fill = paper, colour = NA),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = linec, linewidth = 0.3),
axis.title = element_text(colour = body),
axis.text = element_text(colour = faint),
plot.title = element_text(colour = ink, face = "bold"),
plot.subtitle = element_text(colour = labelc),
legend.title = element_text(colour = body),
legend.text = element_text(colour = body),
strip.text = element_text(colour = ink)
)
}
```
```{r}
#| label: simulate
set.seed(133)
Tyr <- 30; phi_bar <- 0.52; f_bar <- 0.48; sig_e <- 0.06; sy_true <- 8; N1_true <- 110
# environmental shocks perturb the rates each year; demographic draws on top
N <- numeric(Tyr); N[1] <- N1_true
for (t in 1:(Tyr - 1)) {
e <- rnorm(1, 0, sig_e)
phit <- plogis(qlogis(phi_bar) + e / (phi_bar + f_bar) * 2)
ft <- f_bar * exp(e / (phi_bar + f_bar))
N[t + 1] <- rbinom(1, N[t], phit) + rpois(1, ft * N[t])
}
C <- round(rnorm(Tyr, N, sy_true))
# capture-recapture and productivity: pin the mean rates (hence demographic size)
p_true <- 0.55; Kcr <- 10; Rrel <- rep(35, Kcr - 1)
recaps <- vector("list", Kcr - 1); never <- integer(Kcr - 1)
for (i in 1:(Kcr - 1)) {
alive <- rep(TRUE, Rrel[i]); fr <- integer(Rrel[i]); rc <- numeric(Kcr)
for (t in (i + 1):Kcr) {
alive <- alive & (runif(length(alive)) < phi_bar)
got <- alive & (runif(length(alive)) < p_true) & (fr == 0); fr[got] <- t
}
for (t in (i + 1):Kcr) rc[t] <- sum(fr == t)
recaps[[i]] <- rc; never[i] <- Rrel[i] - sum(fr > 0)
}
Jyr <- 12; prodB <- rep(45, Jyr); prodJ <- rpois(Jyr, prodB * f_bar)
```
The population runs from `r N[1]` to `r N[Tyr]` over the `r Tyr` years.
## A process variance with two terms
The count model is a Kalman filter, as in the earlier posts, but the one-step process variance now carries both parts: `N` times the demographic factor `phi(1 - phi) + f`, plus `N` squared times the environmental variance `sigma_e` squared. The demographic factor is fixed by the survival and recruitment that the capture-recapture and productivity data estimate, so the only new quantity is `sigma_e`.
```{r}
#| label: fit
marray_ll <- function(phi, p) {
ll <- 0
for (i in 1:(Kcr - 1)) {
ts <- (i + 1):Kcr; q <- phi^(ts - i) * (1 - p)^(ts - i - 1) * p
ll <- ll + sum(recaps[[i]][ts] * log(q)) + never[i] * log(1 - sum(q))
}
ll
}
kf <- function(phi, f, se2, sy, N1, demog = TRUE) {
lam <- phi + f; dv <- phi * (1 - phi) + f; a <- N1; P <- (2 * sy)^2; ll <- 0
for (t in seq_along(C)) {
if (t > 1) {
ap <- a; a <- lam * ap
Q <- (if (demog) dv * max(ap, 1) else 0) + se2 * max(ap, 1)^2
P <- lam^2 * P + Q
}
Fv <- P + sy^2; ll <- ll + dnorm(C[t], a, sqrt(Fv), log = TRUE)
K <- P / Fv; a <- a + K * (C[t] - a); P <- (1 - K) * P
}
ll
}
nll <- function(par, demog = TRUE) {
-(kf(plogis(par[1]), exp(par[2]), exp(par[3]), exp(par[4]), exp(par[6]), demog) +
marray_ll(plogis(par[1]), plogis(par[5])) +
sum(dpois(prodJ, prodB * exp(par[2]), log = TRUE)))
}
st <- c(qlogis(.5), log(.5), log(0.003), log(9), qlogis(.5), log(C[1]))
of <- optim(st, nll, demog = TRUE, method = "BFGS", hessian = TRUE,
control = list(maxit = 2000, reltol = 1e-11))
stopifnot(of$convergence == 0) # check before inverting the Hessian
Vf <- solve(of$hessian)
phi_h <- plogis(of$par[1]); f_h <- exp(of$par[2]); se2_h <- exp(of$par[3])
sig_e_h <- sqrt(se2_h)
ci_lo <- sqrt(exp(of$par[3] - 1.96 * sqrt(Vf[3, 3])))
ci_hi <- sqrt(exp(of$par[3] + 1.96 * sqrt(Vf[3, 3])))
dv_h <- phi_h * (1 - phi_h) + f_h
Nstar <- dv_h / se2_h
```
The model recovers the mean rates (survival `r round(phi_h, 3)`, recruitment `r round(f_h, 3)`) and estimates the environmental standard deviation of the growth rate at `r round(sig_e_h, 3)`, against a true `r sig_e`. The confidence interval runs from `r round(ci_lo, 3)` to `r round(ci_hi, 3)`: this thirty-year series does not pin the environmental term down. The crossover abundance implied by the fit is `r round(Nstar)`. Where that number falls relative to the abundances the population actually visited is what tells you which term ruled the series.
```{r}
#| label: fig-components
#| fig-cap: "The two contributions to growth-rate variance against population size. The demographic term (rust) falls as one over abundance; the environmental term (green) is flat. They cross at the threshold abundance (dashed), below which chance in a finite population dominates and above which year-to-year environmental variation does. The shaded band is the range the study population actually occupied."
#| fig-alt: "A plot with population size on the x-axis and variance of the per-capita growth rate on the y-axis: a rust curve for the demographic term descending as a hyperbola, a flat green line for the environmental term that the rust curve crosses, a dashed vertical line at that crossover abundance, and a pale shaded band spanning the abundances the population occupied."
#| fig-width: 7
#| fig-height: 5
Ng <- seq(min(N) * 0.8, max(N) * 1.15, length.out = 200)
comp <- rbind(
data.frame(N = Ng, v = dv_h / Ng, part = "demographic"),
data.frame(N = Ng, v = se2_h, part = "environmental")
)
ggplot(comp, aes(N, v, colour = part)) +
annotate("rect", xmin = min(N), xmax = max(N), ymin = -Inf, ymax = Inf,
fill = sage, alpha = 0.13) +
geom_vline(xintercept = Nstar, linetype = "dashed", colour = faint) +
geom_line(linewidth = 0.9) +
annotate("text", x = Nstar, y = max(comp$v) * 0.95,
label = paste0("N* = ", round(Nstar)), colour = faint, hjust = -0.1, size = 3.4) +
scale_colour_manual(values = c(demographic = rust, environmental = forest), name = NULL) +
labs(x = "population size (N)", y = "variance of per-capita growth rate",
title = "Which noise dominates depends on abundance",
subtitle = "demographic variance falls with N, environmental variance does not") +
theme_te() + theme(legend.position = "top")
```
## The cost of ignoring the demographic term
A common shortcut treats all the process variation as environmental, fitting a single per-capita process variance and dropping the demographic term. That single term must now absorb the demographic scatter as well, so it comes out larger than the environmental variance from the full decomposition.
```{r}
#| label: env-only
oe <- optim(st, nll, demog = FALSE, method = "BFGS",
control = list(maxit = 2000, reltol = 1e-11))
sig_e_env <- sqrt(exp(oe$par[3]))
```
Dropping the demographic term pushes the environmental estimate to `r round(sig_e_env, 3)`, about `r round(sig_e_env / sig_e_h, 2)` times the value from the correct decomposition. The point is structural: demographic noise mistaken for environmental noise makes a population look more at the mercy of the environment than it is. Tracing the two absolute variance terms along the trajectory shows how their shares change as the population grew.
```{r}
#| label: fig-decomp
#| fig-cap: "Absolute process variance split into its demographic and environmental parts along the fitted trajectory, ordered by abundance. The demographic part grows in proportion to abundance and the environmental part with its square, so the environmental share of the total rises from left to right. The stack is the total one-step process variance the model attributes to each year."
#| fig-alt: "A stacked area chart with population size on the x-axis and one-step process variance on the y-axis: a rust band for the demographic part and a green band for the environmental part, stacked together to give the total, the stack rising from left to right and the green band taking a larger share of it towards the right."
#| fig-width: 7
#| fig-height: 5
lam <- phi_h + f_h
ap <- Pp <- af <- Pf <- numeric(Tyr); a <- exp(of$par[6]); P <- (2 * exp(of$par[4]))^2
for (t in 1:Tyr) {
if (t > 1) { a <- lam * af[t - 1]; P <- lam^2 * Pf[t - 1] + dv_h * max(af[t - 1], 1) + se2_h * max(af[t - 1], 1)^2 }
ap[t] <- a; Pp[t] <- P
Fk <- P + exp(of$par[4])^2; K <- P / Fk
af[t] <- a + K * (C[t] - a); Pf[t] <- (1 - K) * P
}
as_ <- af
for (t in (Tyr - 1):1) { J <- Pf[t] * lam / Pp[t + 1]; as_[t] <- af[t] + J * (as_[t + 1] - ap[t + 1]) }
dd <- data.frame(N = as_,
demographic = dv_h * as_,
environmental = se2_h * as_^2)
dd <- dd[order(dd$N), ]
dl <- rbind(data.frame(N = dd$N, v = dd$demographic, part = "demographic"),
data.frame(N = dd$N, v = dd$environmental, part = "environmental"))
ggplot(dl, aes(N, v, fill = part)) +
geom_area(alpha = 0.85, colour = paper, linewidth = 0.2) +
scale_fill_manual(values = c(demographic = rust, environmental = forest), name = NULL) +
labs(x = "population size (N)", y = "one-step process variance",
title = "The noise budget shifts as the population grows",
subtitle = "environmental term grows with N squared, demographic share fades") +
theme_te() + theme(legend.position = "top")
```
## Takeaways
Demographic and environmental variability are not interchangeable: one is the sampling noise of a finite population and fades as one over abundance, the other is the variation of the vital rates and holds steady, and which matters depends on how large the population is. An integrated model separates them because the capture-recapture and productivity data fix the vital rates, and so fix how much demographic scatter to expect; the leftover count variation is the environmental part. Two cautions travel with this. The environmental variance is estimated from the spread of a few dozen years and is never sharp. And the separation only works if abundance varies enough to expose the different scalings; a population held near constant size over a short series leaves the two terms confounded. The final post in the series turns from estimation to checking: how to tell whether the data sets an integrated model combines actually agree.
```{r}
#| label: thumb
#| include: false
th <- ggplot(comp, aes(N, v, colour = part)) +
annotate("rect", xmin = min(N), xmax = max(N), ymin = -Inf, ymax = Inf,
fill = sage, alpha = 0.13) +
geom_vline(xintercept = Nstar, linetype = "dashed", colour = faint) +
geom_line(linewidth = 1.1) +
scale_colour_manual(values = c(demographic = rust, environmental = forest), guide = "none") +
labs(x = NULL, y = NULL) + theme_te()
ggsave("thumbnail.png", th, width = 7.5, height = 5.4, dpi = 100, bg = paper)
```
## References
Besbeas, P., Freeman, S. N., Morgan, B. J. T. and Catchpole, E. A. 2002. Biometrics 58(3):540-547 (10.1111/j.0006-341X.2002.00540.x).
Engen, S., Bakke, O. and Islam, A. 1998. Biometrics 54(3):840-846 (10.2307/2533838).
Lande, R., Engen, S. and Saether, B.-E. 2003. Stochastic Population Dynamics in Ecology and Conservation. Oxford University Press (ISBN 978-0-19-852525-7).
Saether, B.-E. and Engen, S. 2015. Trends in Ecology and Evolution 30(5):273-281 (10.1016/j.tree.2015.03.007).
Schaub, M. and Abadi, F. 2011. Journal of Ornithology 152(S1):227-237 (10.1007/s10336-010-0632-7).
Kery, M. and Schaub, M. 2012. Bayesian Population Analysis using WinBUGS. Academic Press (ISBN 978-0-12-387020-9).
## Related tutorials
- [Integrated population models from scratch](../integrated-population-model-basics/)
- [Implicit information in integrated models](../ipm-implicit-information/)
- [Stochastic population growth](../stochastic-population-growth/)
- [The Gompertz state-space model](../gompertz-state-space-model/)
