← Past problems · 2012 set

2012 · Problem A — American Elk

Read the official problem page →

The prompt, restated

In 2001 and 2002 the National Park Service released a small founder population of elk (Cervus canadensis) into Great Smoky Mountains National Park as part of an experimental reintroduction. Several years of monitoring data — herd counts, sex ratio, calf-to-cow ratio, and observed mortality — are summarized in the prompt. Teams are asked: (1) Build a mathematical model that uses this data to project the herd forward and decide whether the population will persist or die out under current vital rates. (2) Test the sensitivity of the conclusion to the parameters that are least well measured (typically calf survival and adult female survival). (3) Propose a management plan — predator control, supplemental releases, habitat enlargement, or selective relocation — that demonstrably improves the long-run growth rate, and quantify how much each lever moves the outcome. (4) Deliver three artifacts: a one-page summary sheet, a complete technical report, and a letter to the Wildlife Commissioner translating the recommendation into plain language with a budget-aware ranking of interventions.

The data are small and noisy — fewer than 100 animals over a handful of years — so the modeling judgment is mostly about how to extract a defensible vital-rate matrix from sparse counts, how to propagate uncertainty into the persistence verdict, and how to map each candidate intervention back to a specific cell of the matrix.

Key modeling idea

The natural framework is a stage-structured population model — a Lefkovitch matrix $L$ with three or four life stages (calf, yearling, adult cow, optionally adult bull). The dynamics are $\mathbf{n}_{t+1} = L\, \mathbf{n}_t$, where each column of $L$ encodes survival to the next stage plus age-specific fecundity for adult females. The dominant eigenvalue $\lambda_1$ of $L$ is the asymptotic growth rate: if $\lambda_1 > 1$ the herd grows, if $\lambda_1 < 1$ it declines toward extinction. The right eigenvector gives the stable stage distribution; the left eigenvector gives reproductive value, which is what tells you which stage to protect first.

Because the herd is small (~25–60 animals over the observed window), the deterministic $\lambda_1$ is not the whole story — demographic stochasticity and occasional catastrophe years (disease, severe winter) can drive the herd extinct even when $\lambda_1 > 1$. The persistence verdict therefore comes from a Monte Carlo population viability analysis (PVA): simulate many replicates with vital rates drawn from their sampling distributions and report the probability of extinction within a chosen horizon (e.g., 50 years).

Suggested approach

• Step 1 — Estimate vital rates from the data. From the prompt's herd counts and calf-to-cow ratios, back out annual calf survival, yearling survival, adult cow survival, and adult fecundity. Use a simple maximum-likelihood / method-of- moments estimator and report a standard error for each rate (technique 4 — regression / parameter fitting).
• Step 2 — Build the Lefkovitch matrix and compute $\lambda_1$. Use a 3-stage model (calf, yearling, adult cow). Report the dominant eigenvalue, the stable stage distribution, and the sensitivity / elasticity of $\lambda_1$ to each matrix entry (Caswell 2001).
• Step 3 — Stochastic PVA. Draw vital rates each year from beta / lognormal distributions calibrated to the Step 1 estimates, run 5,000 replicate 50-year trajectories with demographic stochasticity (Poisson births, binomial survival), and report the extinction probability and quasi-extinction time (technique 10 — Monte Carlo, technique 5 — density-dependent growth for the carrying-capacity ceiling).
• Step 4 — Score interventions. Translate each candidate management lever into a matrix change: predator control raises calf survival by $\Delta s_c$; supplemental releases add $r$ animals to the adult class every $k$ years; habitat expansion raises the carrying capacity $K$. Re-run the PVA per lever and rank by the marginal reduction in extinction probability per dollar.
• Step 5 — Write the three deliverables. The summary sheet states the persistence verdict and one recommendation; the technical report carries the matrix, eigenanalysis, and PVA charts; the commissioner letter avoids notation and speaks in herd-size and budget terms.

Data sources to consider

Common pitfalls

• Using exponential growth without a carrying capacity. A pure $N_{t+1} = \lambda N_t$ model with $\lambda > 1$ predicts infinite elk. Add a density-dependent ceiling at the park's habitat-based carrying capacity (technique 5) — even if it doesn't bind for decades, the judges look for it.
• Treating $\lambda_1 > 1$ as "the herd persists". Small founder populations can wink out from stochasticity alone. Always report an explicit extinction probability with a horizon.
• Ignoring sex ratio. A two-sex model is overkill, but if the bull- to-cow ratio is heavily skewed by selective mortality you must show the cow class drives $\lambda_1$ and bulls can be aggregated.
• One-shot sensitivity. Vary each parameter alone, then jointly via Latin-hypercube sampling. Adult-cow survival typically dominates elasticity — if your ranking says calf survival dominates, recheck the matrix.
• Cost-blind interventions. Predator control and supplemental releases have very different price tags. Without a per-dollar effect estimate, the commissioner letter is not actionable.
• Confusing data-fitting with forecasting. The herd counts you have are autocorrelated; do not treat each annual count as an independent draw when computing uncertainty on $\lambda$.

Python sketch

Lefkovitch matrix + eigenanalysis + stochastic PVA with intervention comparison.

import numpy as np

rng = np.random.default_rng(0)

# --- Step 1: vital rates [illustrative, drawn from NPS GRSM reports] ---
s_c, s_y, s_a = 0.55, 0.80, 0.88      # calf, yearling, adult-cow annual survival
f_a           = 0.45                  # female calves per adult cow per year
K             = 250                   # park carrying capacity (cows-equivalent)

def lefkovitch(s_c, s_y, s_a, f_a):
    return np.array([
        [0.0, 0.0, f_a],              # calves from adult cows
        [s_c, 0.0, 0.0],              # calves -> yearlings
        [0.0, s_y, s_a],              # yearlings + adults -> adults
    ])

L = lefkovitch(s_c, s_y, s_a, f_a)
lam = np.linalg.eigvals(L).real.max()
print(f"deterministic lambda_1 = {lam:.3f}")  # [illustrative] ~1.05

# --- Step 3: stochastic PVA ---
def simulate(L_mean, n0, years=50, reps=5000, K=250,
             sigma=(0.10, 0.06, 0.04, 0.08)):
    extinct = 0
    final   = np.zeros(reps)
    for r in range(reps):
        n = np.array(n0, dtype=float)
        for _ in range(years):
            # environmental stochasticity on each rate (beta-clipped)
            sc = np.clip(rng.normal(s_c, sigma[0]), 0, 1)
            sy = np.clip(rng.normal(s_y, sigma[1]), 0, 1)
            sa = np.clip(rng.normal(s_a, sigma[2]), 0, 1)
            fa = max(0.0, rng.normal(f_a, sigma[3]))
            # density dependence on fecundity
            fa *= max(0.0, 1 - n.sum() / K)
            Lr = lefkovitch(sc, sy, sa, fa)
            # demographic stochasticity (Poisson births, binomial survival)
            births  = rng.poisson(Lr[0, 2] * n[2])
            surv_c  = rng.binomial(int(n[0]), sc)
            surv_y  = rng.binomial(int(n[1]), sy)
            surv_a  = rng.binomial(int(n[2]), sa)
            n = np.array([births, surv_c, surv_y + surv_a], dtype=float)
            if n.sum() < 2:
                extinct += 1
                break
        final[r] = n.sum()
    return extinct / reps, final.mean()

n0 = [5, 5, 15]   # rough founder stage distribution [illustrative]
p_ext, n_mean = simulate(L, n0)
print(f"baseline: P(extinct in 50 yr) = {p_ext:.2%}, E[N_50] = {n_mean:.0f}")

# --- Step 4: intervention sweep ---
scenarios = {
    "baseline":              (s_c,        s_a,  0),
    "+ predator control":    (s_c + 0.10, s_a,  0),
    "+ adult protection":    (s_c,        min(0.95, s_a + 0.04), 0),
    "+ release 5 cows/3yr":  (s_c,        s_a,  5),
}
for name, (sc, sa, release) in scenarios.items():
    Lx = lefkovitch(sc, s_y, sa, f_a)
    # (release handling omitted for brevity — add to adult class every 3 years)
    pe, _ = simulate(Lx, n0)
    print(f"{name:24s} P(extinct) = {pe:.2%}")

Sensitivity & validation checklist

• Compute eigen-elasticities of $\lambda_1$ — adult-cow survival typically has the largest elasticity (~0.5 [illustrative]); calf survival is usually second.
• Latin-hypercube sample all four vital rates within their 95% CIs and report the fraction of draws yielding $\lambda_1 > 1$. If less than 70%, the persistence verdict is fragile and the letter should say so.
• Vary the carrying capacity $K$ by ±50% — does the 50-year extinction probability move by more than 5 percentage points? If yes, name $K$ as a research priority.
• Cross-validate against the Pennsylvania or Kentucky elk reintroduction trajectories — your model should reproduce their observed 10-year growth within ~20%.
• Confirm the stable stage distribution from the right eigenvector matches the observed stage proportions in the latest census within sampling error.
• Run the PVA at 25-year and 100-year horizons — extinction probability should be monotonically non-decreasing; if it isn't, you have a bug.

Related pages

• Techniques — Logistic / density-dependent growth
• Techniques — Markov chains (matrix-population dynamics)
• Techniques — Monte Carlo simulation (PVA)
• Techniques — Parameter fitting
• Anatomy of a winning paper