The Robust Design

The problem we are stuck with

Right, where are we? Two models, and each one solves half the problem.

Lincoln-Petersen gives us population size \(N\), but only if we hold the population still with closure. Nothing changes between your two visits. That’s fine over a day or two. Stretch it over weeks or months and it’s pure fiction.

Cormack-Jolly-Seber gives us apparent survival \(\phi\), but it needs the population to be open. Animals are allowed to die between occasions, which is the whole point, because dying is the thing we’re trying to measure. The price is that CJS only looks at marked individuals and throws away \(N\) completely.

So we’re stuck with two half-answers. One knows about \(N\) but goes quiet on survival. The other knows about survival but can’t count heads. And of course what we actually want is both, out of the same study, at the same time.

This isn’t some minor annoyance we can shrug off. Almost every real conservation question needs both numbers at once. How many animals are there, and are enough of them surviving to keep that number up? Ask one without the other and you’ve barely said anything.

The robust design, which Kenneth Pollock came up with in 1982, fixes all of this with an idea so simple it’s almost annoying once it clicks.

The key insight: two time scales

LP and CJS fight because they want opposite things from the population. LP wants it frozen. CJS wants it changing. The robust design ends the fight by running at two different time scales at once, and handing each model exactly the thing it was asking for.

Here’s the shape of it:

library(ggplot2)
library(dplyr)

# Primary periods
primary <- data.frame(
  x_start = c(1, 4, 7, 10),
  x_end   = c(3, 6, 9, 12),
  y = 0.5,
  label = paste("Primary period", 1:4)
)

# Secondary occasions within each primary period
secondary <- data.frame(
  x = c(1.25, 2, 2.75,
         4.25, 5, 5.75,
         7.25, 8, 8.75,
         10.25, 11, 11.75),
  y = 0.5,
  primary = rep(1:4, each = 3)
)

# Open arrows between primary periods
arrows_df <- data.frame(
  x    = c(3.05, 6.05, 9.05),
  xend = c(3.95, 6.95, 9.95),
  y = 0.5
)

ggplot() +
  # Primary period boxes
  geom_rect(data = primary,
            aes(xmin = x_start, xmax = x_end,
                ymin = 0.15, ymax = 0.85),
            fill = "#d19527", alpha = 0.15, colour = "#d19527", linewidth = 1) +
  # Primary period labels
  geom_text(data = primary,
            aes(x = (x_start + x_end) / 2, y = 0.92, label = label),
            size = 3.2, colour = "#d19527", fontface = "bold") +
  # Secondary occasion points
  geom_point(data = secondary,
             aes(x = x, y = y),
             size = 5, colour = "#FF5733") +
  geom_text(data = secondary,
            aes(x = x, y = y),
            label = "s", size = 3, colour = "white", fontface = "bold") +
  # Open arrows between primary periods
  geom_segment(data = arrows_df,
               aes(x = x, xend = xend, y = y, yend = y),
               arrow = arrow(length = unit(0.25, "cm"), ends = "both"),
               linewidth = 0.8, colour = "grey40", linetype = "dashed") +
  # Labels for arrows
  annotate("text", x = c(3.5, 6.5, 9.5), y = 0.62,
           label = "Open\n(\u03d5)", size = 3, colour = "grey40") +
  # Closed label inside boxes
  annotate("text", x = c(2, 5, 8, 11), y = 0.22,
           label = "Closed\n(N, p)", size = 3, colour = "#d19527") +
  scale_x_continuous(limits = c(0.5, 13)) +
  scale_y_continuous(limits = c(0, 1.1)) +
  labs(x = NULL, y = NULL,
       title = "The robust design sampling structure",
       subtitle = "Orange dots = secondary occasions (s). Gold boxes = primary periods.") +
  theme_dark_site() +
  theme(axis.text = element_blank(),
        panel.grid = element_blank())

The primary periods are spaced far enough apart that the population has time to change between them. Animals die. New ones get recruited. This is the open part, run on CJS-style logic, and it’s what gives us \(\phi\).

The secondary occasions are the repeated visits inside a single primary period, packed close enough together in time that we can happily pretend the population is closed for that little window. No births, no deaths, nobody moving in or out for good. This is the closed part, run on LP-style logic, and it’s what gives us \(N\).

The robust design just nests one inside the other. Closed within a primary period. Open between them. LP and CJS stop treading on each other because they’re now working at different time scales.

That’s it. That’s the whole idea. Everything from here is details and bookkeeping.

The data structure

The data falls out straight from the way you sampled. Each animal now carries a capture history that works on two levels at once.

At the between-period level, each primary period gets boiled down to a single yes/no: did we see this animal at all during that period? That coarse summary is what feeds the open-population half of the model, the bit that estimates survival.

At the within-period level, you keep the full blow-by-blow of which secondary occasions the animal turned up on inside each primary period. That fine-grained detail feeds the closed-population half, the bit that estimates abundance and detection.

set.seed(7)

# 6 animals, 3 primary periods, 3 secondary occasions each
animals <- paste0("Animal ", 1:6)
periods <- paste0("P", rep(1:3, each = 3), "_s", rep(1:3, times = 3))

ch_matrix <- matrix(
  c(1,1,0, 1,0,1, 0,0,1,
    1,0,1, 0,0,0, 1,1,0,
    0,1,1, 1,1,0, 0,1,1,
    1,1,1, 0,1,0, 1,0,1,
    0,0,1, 1,0,1, 0,0,0,
    1,0,0, 1,1,1, 0,1,0),
  nrow = 6, byrow = TRUE
)

colnames(ch_matrix) <- periods
rownames(ch_matrix) <- animals

ch_df <- as.data.frame(ch_matrix) |>
  tibble::rownames_to_column("animal") |>
  tidyr::pivot_longer(-animal, names_to = "occasion", values_to = "detected") |>
  mutate(
    primary = as.integer(substr(occasion, 2, 2)),
    secondary = as.integer(substr(occasion, 5, 5)),
    detected_label = if_else(detected == 1, "Detected", "Not detected")
  )

ggplot(ch_df, aes(x = secondary, y = animal, fill = detected_label)) +
  geom_tile(colour = "white", linewidth = 1.2) +
  facet_grid(~ paste("Primary period", primary),
             switch = "x") +
  scale_fill_manual(values = c("Detected" = "#d19527",
                                "Not detected" = "#FF5733")) +
  scale_x_continuous(breaks = 1:3, labels = paste("s", 1:3)) +
  labs(x = "Secondary occasion", y = NULL, fill = NULL,
       title = "Robust design capture histories",
       subtitle = "Each primary period contains three secondary occasions") +
  theme_dark_site() +
  theme(legend.position = "bottom",
        panel.grid = element_blank(),
        strip.placement = "outside")

Every red tile still hides the same old ambiguity. The difference now is that the ambiguity shows up at two levels, and the model pulls them apart and deals with each on its own.

Inside a primary period, a red tile can only mean one thing: the animal was alive and present and you just missed it. We’ve assumed the population is closed for that window, so death and emigration are off the table. Imperfect detection is the only suspect left.

Between primary periods it’s murkier. An animal that never showed up across a whole primary period might have died before it started, or it might have been alive the entire time and dodged every single secondary occasion. That’s the ambiguity the open-population part of the model soaks up, exactly like CJS did a page ago.

Building the likelihood

The robust design likelihood splits into two neat halves. They get fitted together in one go, but you can understand them one at a time, so that’s exactly what we’ll do. One half, then the other.

Part 1: Within primary periods (closed population)

Inside a single primary period the population is closed. We’ve got a handful of repeated secondary occasions, and we want two things out of them: the detection probability \(p\), and the population size \(N_t\) for that period.

This is really just a closed-population mark-recapture model. The simplest flavour, sometimes called \(M_0\), assumes every animal has the same detection probability \(p\) on every secondary occasion. Write an animal’s within-period history as \(\mathbf{x}_{it} = (x_{it1}, x_{it2}, \ldots, x_{its})\) across \(s\) secondary occasions. Given the animal was actually there, the probability of seeing that exact history is:

\[\Pr(\mathbf{x}_{it} \mid \text{present}) = \prod_{j=1}^{s} p^{x_{itj}} (1-p)^{1-x_{itj}}\]

Don’t let the symbols spook you, it’s a Bernoulli product and nothing scarier. Each secondary occasion, you either caught the animal (\(x = 1\), which costs you a \(p\)) or you didn’t (\(x = 0\), which costs you a \(1-p\)). Multiply those across all the occasions and you’ve got the probability of the whole history.

Now flip it around. The probability of missing the animal completely across all \(s\) occasions in a period is every occasion going wrong in a row:

\[1 - \tilde{p}_t = (1-p)^s\]

So the probability of catching it at least once somewhere in primary period \(t\) is just one minus that:

\[\tilde{p}_t = 1 - (1-p)^s\]

That \(\tilde{p}_t\) is your detection probability for the whole period, and it’s the quiet little hinge that joins the within-period model to the between-period one. And here’s the payoff worth noticing: even if \(p\) on any single day is fairly rubbish, stack up enough secondary occasions and \(\tilde{p}_t\) climbs fast. Three cracks at 40% each and you’re already catching the large majority of animals that are actually there. That’s the whole reason we bother going back multiple times within a season.

To actually get \(N_t\) we lean on a Huggins-style conditional likelihood (you really don’t need to remember that name). The trick is that we never go at \(N_t\) head-on. Instead we condition on the animals we did catch at least once, estimate \(p\) from their histories, and then back out \(\hat{N}_t\) afterwards as:

\[\hat{N}_t = \frac{n_t}{\hat{\tilde{p}}_t}\]

where \(n_t\) is the number of distinct animals we caught in period \(t\). If that looks familiar, it should: it’s the same Lincoln-Petersen move in fancier clothes. We saw \(n_t\) animals, we reckon each one had a \(\hat{\tilde{p}}_t\) chance of being seen, so the population that must have been sitting out there is \(n_t\) divided by that probability.

Part 2: Between primary periods (open population)

Zoom out to the gaps between primary periods and the CJS machinery takes over. For each animal we now collapse its history down to one summary per primary period: \(\omega_{it} = 1\) if we caught it on any secondary occasion in period \(t\), and 0 if we never saw it that whole period.

The between-period likelihood is the exact CJS likelihood from the last page, with one swap. The detection probability feeding into it isn’t the raw single-day \(p\) any more. It’s the effective per-period detection \(\tilde{p}_t\) we just built back in Part 1.

\[\Pr(\omega_{it} = 1 \mid \text{alive at } t) = \tilde{p}_t\]

\[\Pr(\omega_{it} = 0 \mid \text{alive at } t) = 1 - \tilde{p}_t\]

Survival \(\phi_t\) works between primary periods, just like it did in CJS:

\[z_{i,t} \sim Bernoulli(\phi_{t-1} \times z_{i,t-1})\]

And what we actually observe, given whether the animal is alive, is:

\[\omega_{i,t} \sim Bernoulli(\tilde{p}_t \times z_{i,t})\]

The probability of never clapping eyes on the animal again from period \(t\) onward, \(\chi_t\), follows the same backwards recursion as before, except now \(\tilde{p}_t\) stands in for the old single-occasion \(p\):

\[\chi_T = 1\]

\[\chi_t = (1 - \phi_t) + \phi_t (1 - \tilde{p}_{t+1}) \chi_{t+1}\]

The full likelihood

Bolt the two halves together and the full robust design likelihood is just their product:

\[\mathcal{L} = \mathcal{L}_{\text{open}}(\phi, \tilde{p}) \times \prod_{t=1}^{T} \mathcal{L}_{\text{closed},t}(p_t)\]

The closed half gets worked out within each primary period and hands you \(\hat{p}_t\) and \(\hat{N}_t\). The open half takes the \(\tilde{p}_t\) built from those \(\hat{p}_t\) and hands you \(\hat{\phi}_t\). The whole thing hangs together because the two halves gossip to each other through \(\tilde{p}_t\). That single quantity is the bridge between the two time scales.

What we can now estimate

Here’s the haul. The robust design hands you a set of parameters that neither LP nor CJS could give you on their own:

Parameter	Symbol	Source
Population size at each primary period	\(N_t\)	Closed component
Detection probability (single occasion)	\(p\)	Closed component
Effective detection probability (per primary period)	\(\tilde{p}_t\)	Derived from \(p\) and \(s\)
Apparent survival between primary periods	\(\phi_t\)	Open component
Population growth rate	\(\lambda_t = N_{t+1}/N_t\)	Derived from \(N_t\)

\(\lambda_t\) is probably the single most useful thing on that list for day-to-day work. It tells you straight out whether the population is growing, holding steady, or sliding downhill between primary periods. And crucially it’s an honest number, because it has already accounted for the animals you missed, both at the single-occasion level and across a whole primary period. It won’t flatter you the way a raw count would.

A small simulation

Let’s build a robust design dataset from scratch and watch how the estimates behave when we already know the true answers. We’ll use six primary periods with three secondary occasions each, and keep it deliberately simple with a constant \(\phi\) and \(p\). Six primary periods gives the model five survival intervals to chew on, which is plenty to pin down \(\phi\).

Quick heads up: the modelling here uses a package called RMark. I’ll properly get into it later, but for now just know it’s a front end for a separate bit of software called program MARK, which has to be installed for any of this to run. If you want to grab it now, it’s here.

library(RMark)

set.seed(1988)

N_true      <- 300    # Initial population size
phi_true    <- 0.80   # Survival between primary periods
p_true      <- 0.40   # Detection on each secondary occasion
n_primary   <- 6
n_secondary <- 3

# Effective detection probability per primary period
p_tilde <- 1 - (1 - p_true)^n_secondary

cat("True N:", N_true,
    "\nTrue phi:", phi_true,
    "\nTrue p (single occasion):", p_true,
    "\nEffective p per primary period:", round(p_tilde, 3))

True N: 300 
True phi: 0.8 
True p (single occasion): 0.4 
Effective p per primary period: 0.784

# Simulate capture histories
sim_robust <- function(N, phi, p, n_prim, n_sec, seed = 1988) {
  set.seed(seed)
  
  # True alive state across primary periods
  alive <- matrix(0, nrow = N, ncol = n_prim)
  alive[, 1] <- 1
  for (t in 2:n_prim) {
    alive[, t] <- rbinom(N, 1, alive[, t - 1] * phi)
  }
  
  # Detection history: N x (n_prim * n_sec)
  det <- matrix(0, nrow = N, ncol = n_prim * n_sec)
  for (t in 1:n_prim) {
    for (s in 1:n_sec) {
      col <- (t - 1) * n_sec + s
      det[, col] <- rbinom(N, 1, alive[, t] * p)
    }
  }
  
  # Keep only animals detected at least once
  det[rowSums(det) > 0, ]
}

ch_data <- sim_robust(N_true, phi_true, p_true, n_primary, n_secondary)

ch_strings <- apply(ch_data, 1, paste, collapse = "")
rd_df <- data.frame(ch = ch_strings, stringsAsFactors = FALSE)

cat("Number of individuals detected at least once:", nrow(rd_df), "\n")

Number of individuals detected at least once: 281 

cat("First 10 capture histories:\n")

First 10 capture histories:

head(ch_strings, 10)

 [1] "101101100010111110" "000001101111100001" "111111000000011000"
 [4] "001100000000000000" "110101011110011101" "100010000000000000"
 [7] "010011110000000000" "000101010001000110" "100000000000000000"
[10] "110100100100000000"

# Process data for robust design in RMark
# time.intervals: 0 = within primary period (closed), 1 = between primary periods (open)
# 6 primary periods x 3 secondary occasions = 18 columns
# Intervals: 0 0 | 1 | 0 0 | 1 | 0 0 | 1 | 0 0 | 1 | 0 0 | 1 | 0 0
time_intervals <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0)

rd_processed <- process.data(rd_df,
                              model = "Robust",
                              time.intervals = time_intervals)

rd_ddl <- make.design.data(rd_processed)

# Fit model with constant phi and p, no temporary emigration
# GammaPrime and GammaDoublePrime fixed to -10 on logit scale (effectively zero)
# We will deal with temporary emigration properly on the next page
rd_fit <- mark(rd_processed, rd_ddl,
               model.parameters = list(
                 S              = list(formula = ~ 1),
                 p              = list(formula = ~ 1),
                 GammaPrime     = list(formula = ~ 1, fixed = -10),
                 GammaDoublePrime = list(formula = ~ 1, fixed = -10)
               ),
               output = FALSE,
               silent = TRUE)

beta <- rd_fit$results$beta

phi_est <- plogis(beta["S:(Intercept)", "estimate"])
phi_lcl <- plogis(beta["S:(Intercept)", "estimate"] -
                    1.96 * beta["S:(Intercept)", "se"])
phi_ucl <- plogis(beta["S:(Intercept)", "estimate"] +
                    1.96 * beta["S:(Intercept)", "se"])

p_est <- plogis(beta["p:(Intercept)", "estimate"])
p_lcl <- plogis(beta["p:(Intercept)", "estimate"] -
                  1.96 * beta["p:(Intercept)", "se"])
p_ucl <- plogis(beta["p:(Intercept)", "estimate"] +
                  1.96 * beta["p:(Intercept)", "se"])

cat("Survival estimate (true =", phi_true, "):",
    round(phi_est, 3),
    "  95% CI:", round(phi_lcl, 3), "to", round(phi_ucl, 3))

Survival estimate (true = 0.8 ): 0.781   95% CI: 0.75 to 0.809

cat("\nDetection estimate (true =", p_true, "):",
    round(p_est, 3),
    "  95% CI:", round(p_lcl, 3), "to", round(p_ucl, 3))

Detection estimate (true = 0.4 ): 0.025   95% CI: 0.024 to 0.026

Those estimates should land pretty close to the true values we picked. They won’t be bang on, because we’re working from a single simulated sample and samples wobble, but the true values should sit comfortably inside the confidence intervals.

What is still missing

The model we just fitted quietly pinned temporary emigration to zero. That was a deliberate cheat to keep things clean. In the real world, animals wander off out of your study area between primary periods and come swanning back later. They aren’t dead and they aren’t gone for good, but while they’re away you can’t detect them no matter how hard you look. And yes, that’s yet another flavour of ambiguous non-detection to throw on the pile.

That’s temporary emigration, and it’s handled by two parameters, \(\gamma'\) and \(\gamma''\) (gamma prime and gamma double prime), that we’ve swept under the rug for now. The whole next page is about dragging them back out.

For now, the thing to hang onto is the core structure. Closed within a primary period gives you \(N_t\) and \(p\). Open between primary periods gives you \(\phi\). The two time scales pull together through \(\tilde{p}_t\). That, right there, is the robust design.