Implementing the Robust Design in RMark

Author

Deon Roos

Published

July 24, 2026

Putting it all together

Up to now we’ve built the robust design up from first principles. You know what the model’s doing, why it needs two time scales, what each parameter means, and where all the ambiguous zeros are hiding. Time to actually point it at some data.

This page is the practical one. Less theory, more code. By the end you’ll be able to take a real robust design dataset, bully it into the right shape, fit a model with covariates, and draw figures that actually say something about what you found. This is the page you’ll be flicking back to when you’re elbow-deep in your own data and can’t remember which argument goes where.

We’re sticking with RMark the whole way. Quick reminder: RMark is just an R front end for program MARK, which you have to install separately. If you haven’t done that yet, grab it here.

What should influence your parameters?

Before you fit anything, stop and think hard about which covariates go where. This is a biological question, not a statistical one, and if you get it wrong here then it doesn’t matter how cleanly the model runs, it won’t be describing the real world.

The full robust design has four parameters, and each one means something different biologically. That biology is what should be steering your covariate choices, so let’s take them one at a time.

What influences survival ($\phi$)?

Survival is estimated between primary periods, so the covariates that belong here are things that vary between those periods and could plausibly tip an animal towards living or dying.

Body condition is the classic one. An animal going into winter in poor nick is less likely to see spring. Age matters for most species, with survival usually lowest in juveniles and the very old. Sex can matter if males and females are up against different predation or energetic demands. Conditions between periods get thrown in a lot too: how brutal the winter was, how bad the drought got, how much food was about. And if you ran any interventions, supplementary feeding, a predator cull, a disease treatment, those go here as well.

The question to keep asking yourself is: between two of my sampling occasions, what decided whether an animal was still alive at the next one?

What influences detection ($p$)?

Detection is estimated within primary periods, out of the pattern across your secondary occasions. So covariates here need to vary at the secondary occasion level, across the repeated visits inside a single primary period, and plausibly change how likely you are to catch or spot an animal that’s genuinely there.

Weather on the day is the obvious one: rain, temperature and wind all mess with trap success and how well an observer can see. Time of day matters for plenty of species. Effort belongs here too if it wobbled between occasions, like how many traps you set or how many hours you searched. And if different people did the sampling, who was out that day can go in as well.

The question here is: across my repeated visits within a season, what made me more or less likely to detect an animal that was actually sitting right there?

What influences availability ($\gamma'$ and $\gamma''$)?

Availability is about whether an animal is even in your study area during a primary period at all. So covariates here should be the things that push animals in and out of the study area between periods.

Season is the big one, especially for anything migratory. A covariate that just flags breeding versus non-breeding season will soak up a lot of predictable coming and going. Good habitat outside your study area might tempt animals away. And for some species reproductive stage counts: pregnant or lactating females often sit tighter to a site than animals that aren’t breeding.

The structure of covariates in RMark

RMark splits covariates into two kinds.

Individual covariates are fixed traits of an animal: sex, age at first capture, a body condition measured just the once. Each one goes straight into the data frame as its own column next to ch, a single value per row.

Time-varying covariates change from occasion to occasion: the temperature on each sampling day, the season at each primary period. For survival and availability these vary at the primary period level, so you add them as numbered columns (varname1, varname2, and so on, one per primary period). For detection they vary at the secondary occasion level, so you add one column per secondary occasion across the whole study. RMark reads the number on the end and lines each value up with the right row of the design matrix for you.

How should your data be organised?

The capture history string

The basic unit of data in RMark is the capture history string: one string of 1s and 0s per animal, a single digit for every sampling occasion.

A robust design with 3 primary periods of 3 secondary occasions each gives you a 9-character string. Six primary periods of 3 secondary occasions gives you 18. The order is all the secondary occasions of primary period 1, then all of primary period 2, and so on down the line.

Animal tagged in period 1, missed in period 2, detected in period 3:
Primary 1       Primary 2       Primary 3
s1 s2 s3        s1 s2 s3        s1 s2 s3
1  1  0         0  0  0         0  1  1

String: "110000011"

The data frame you hand to RMark needs, at the bare minimum, a column called ch holding these strings. Individual covariates ride alongside as extra columns. Time-varying covariates go in as those numbered column sets.

Let’s simulate a dataset to play with. We’ll use body condition as an individual covariate that drives survival, and temperature as a time-varying covariate that drives detection.

Code

library(RMark)
library(ggplot2)
library(dplyr)
library(tidyr)

set.seed(1988)

N_true      <- 400
phi_int     <- 0.75   # Baseline survival (logit scale intercept)
phi_cond    <- 0.50   # Effect of body condition on survival (logit scale)
p_int       <- -0.40  # Baseline detection (logit scale intercept)
p_temp      <- 0.60   # Effect of temperature on detection (logit scale)
gp_true     <- 0.40   # GammaPrime: prob of staying unavailable
gdp_true    <- 0.20   # GammaDoublePrime: prob of becoming unavailable
n_primary   <- 6
n_secondary <- 3
n_occasions <- n_primary * n_secondary

# Individual covariate: one value per animal
body_condition <- rnorm(N_true, mean = 0, sd = 1)

# Occasion-level covariate: one value per animal per secondary occasion
temp_matrix <- matrix(
  rnorm(N_true * n_occasions, mean = 0, sd = 1),
  nrow = N_true,
  ncol = n_occasions
)

# Simulate survival
alive <- matrix(0, nrow = N_true, ncol = n_primary)
alive[, 1] <- 1
for (t in 2:n_primary) {
  phi_i <- plogis(phi_int + phi_cond * body_condition)
  alive[, t] <- rbinom(N_true, 1, alive[, t - 1] * phi_i)
}

# Simulate availability
avail <- matrix(0, nrow = N_true, ncol = n_primary)
avail[, 1] <- 1
for (t in 2:n_primary) {
  for (i in 1:N_true) {
    if (alive[i, t] == 0) {
      avail[i, t] <- 0
    } else if (avail[i, t - 1] == 1) {
      avail[i, t] <- rbinom(1, 1, 1 - gdp_true)
    } else {
      avail[i, t] <- rbinom(1, 1, 1 - gp_true)
    }
  }
}

# Simulate detection incorporating temperature effect
det <- matrix(0, nrow = N_true, ncol = n_occasions)
for (t in 1:n_primary) {
  for (s in 1:n_secondary) {
    col <- (t - 1) * n_secondary + s
    p_col <- plogis(p_int + p_temp * temp_matrix[, col])
    det[, col] <- rbinom(N_true, 1, avail[, t] * p_col)
  }
}

# Keep only animals detected at least once
ever_seen  <- rowSums(det) > 0
det_seen   <- det[ever_seen, ]
cond_seen  <- body_condition[ever_seen]
temp_seen  <- temp_matrix[ever_seen, ]
n_seen     <- nrow(det_seen)

Code

# Build capture history strings
ch_strings <- apply(det_seen, 1, paste, collapse = "")

# Start data frame with ch and individual covariate
rd_df <- data.frame(
  ch             = ch_strings,
  body_cond   = cond_seen,
  stringsAsFactors = FALSE
)

# Add temperature as time-varying columns: temp1 through temp18
# RMark matches these to detection occasions by column number
temp_cols <- as.data.frame(temp_seen)
colnames(temp_cols) <- paste0("temp", 1:n_occasions)

rd_df <- bind_cols(rd_df, temp_cols)

The data frame now holds ch, the individual covariate body_cond (body condition), and 18 temperature columns, temp1 through temp18. When you drop temp into the detection formula, RMark uses that number on the end to match each temperature to the sampling occasion it belongs to.

Right, let’s actually look at the thing, because this is the exact shape your own data needs to be in before RMark will touch it. It’s too wide to sensibly print all 20 columns at once, so here are the first six animals with just the first three temperature columns (I’ve rounded the covariates to two decimals so the table stays readable):

Code

# First 6 animals: ch, the individual covariate, and the first 3 temp columns
preview <- head(rd_df[, c("ch", "body_cond", "temp1", "temp2", "temp3")])
preview$body_cond <- round(preview$body_cond, 2)
preview[, c("temp1", "temp2", "temp3")] <- round(preview[, c("temp1", "temp2", "temp3")], 2)
preview

                  ch body_cond temp1 temp2 temp3
1 000101101001000111      0.14  0.13 -0.39  2.82
2 100000000000000000     -0.34 -0.48 -0.35 -1.65
3 110000000000000000     -1.04  1.03  0.59  1.59
4 111000000000000000     -0.48  0.18  0.91  1.74
5 111100101100000000      1.26  0.10  1.76  0.71
6 100100110000000000      0.94  1.34  0.22 -0.87

And the full dimensions, so you can see how many animals made it through the “detected at least once” filter and the total number of columns:

Code

dim(rd_df)

[1] 365  20

Read one row at a time. The ch string is that animal’s whole 18-occasion detection history: six primary periods of three secondary occasions each, stitched end to end into one string of 1s and 0s. body_cond is its single body-condition value, and it’s the same number all the way across the row, because body condition is a property of the animal rather than of any one sampling day. Each temp column is the temperature on one specific secondary occasion, so those numbers do change across the row. Your own data wants to look exactly like this: one row per animal, one ch string, your individual covariates as single columns, and any time-varying covariate spread across one numbered column per occasion.

Setting up `time.intervals`

The time.intervals argument is how you tell RMark which gaps between consecutive occasions sit inside a primary period (closed, so 0) and which ones jump between primary periods (open, so 1). There’s one entry per gap, so 18 occasions means 17 gaps.

Occasion:   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
Gap:          1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
Interval:     0  0  1  0  0  1  0  0  1  0  0  1  0  0  1  0  0

Inside primary period 1 (occasions 1, 2, 3) gaps 1 and 2 are 0. The jump from primary period 1 to 2 is gap 3, so that’s a 1. Same pattern all the way along.

Code

build_time_intervals <- function(n_primary, n_secondary) {
  block     <- c(rep(0, n_secondary - 1), 1)
  intervals <- rep(block, n_primary)
  intervals[-length(intervals)]
}

time_intervals <- build_time_intervals(n_primary, n_secondary)
time_intervals

 [1] 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0

Always check the length: it has to come out to (n_primary * n_secondary) - 1.

Code

cat("Length of time_intervals:", length(time_intervals),
    "\nExpected:", n_primary * n_secondary - 1)

Length of time_intervals: 17 
Expected: 17

Writing out the model

Before you type a single line of R, write the whole model out by hand. Three reasons. It drags every assumption you’re making out into the open. It means your methods section is basically written already. And it catches the moments where your biological thinking and your actual code have quietly drifted apart from each other.

For our simulated data, here’s the model:

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

\[\text{logit}(\phi_{i,t}) = \beta_0 + \beta_1 \times \text{Condition}_i\]

\[a_{i,t} \mid z_{i,t} = 1 \sim \begin{cases} \text{Bernoulli}(1 - \gamma'') & \text{if } a_{i,t-1} = 1 \\ \text{Bernoulli}(1 - \gamma') & \text{if } a_{i,t-1} = 0 \end{cases}\]

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

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

\[\text{logit}(p_{i,t,s}) = \alpha_0 + \alpha_1 \times \text{Temperature}_{i,t,s}\]

where:

$z_{i,t}$ is the true alive/dead state of individual $i$ at primary period $t$
$\phi_{i,t}$ is apparent survival from primary period $t$ to $t+1$, varying by individual through body condition
$\beta_0$ is the survival intercept and $\beta_1$ is the effect of body condition
$a_{i,t}$ is the availability state of individual $i$ at primary period $t$, conditional on being alive
$\gamma''$ is the probability of becoming unavailable given currently available
$\gamma'$ is the probability of staying unavailable given currently unavailable
$\omega_{i,t}$ is 1 if individual $i$ was detected at least once during primary period $t$, 0 otherwise
$\tilde{p}_t$ is the effective detection probability across all $s$ secondary occasions
$p_{i,t,s}$ is detection probability on secondary occasion $s$ of primary period $t$, varying by temperature
$\alpha_0$ is the detection intercept and $\alpha_1$ is the effect of temperature

Once it’s laid out like that, turning it into RMark formula syntax is almost mechanical. $\beta_0 + \beta_1 \times \text{Condition}$ becomes ~ body_cond. $\alpha_0 + \alpha_1 \times \text{Temperature}$ becomes ~ temp. The gammas are held constant here, so they each just get ~ 1.

Fitting the model

Code

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

rd_ddl <- make.design.data(rd_processed)

rd_fit <- mark(
  rd_processed,
  rd_ddl,
  model.parameters = list(
    S                = list(formula = ~ body_cond),
    p                = list(formula = ~ temp),
    GammaPrime       = list(formula = ~ 1),
    GammaDoublePrime = list(formula = ~ 1)
  ),
  output = FALSE,
  silent = TRUE
)

Code

rd_fit$results$beta

                               estimate        se        lcl        ucl
S:(Intercept)                 0.6836572 0.0868179  0.5134941  0.8538204
S:body_cond                   0.5662978 0.0880560  0.3937080  0.7388876
GammaDoublePrime:(Intercept) -1.8358266 0.4107363 -2.6408698 -1.0307834
GammaPrime:(Intercept)       -1.7066375 1.6283436 -4.8981911  1.4849161
p:(Intercept)                -0.4947268 0.1207431 -0.7313832 -0.2580704
p:temp                        0.1953920 0.0539919  0.0895678  0.3012162
c:(Intercept)                -0.3779775 0.0634956 -0.5024289 -0.2535262
f0:(Intercept)                4.7518438 0.2117942  4.3367270  5.1669605
f0:session2                  -0.9039858 0.2115448 -1.3186137 -0.4893580
f0:session3                  -1.3291764 0.2481375 -1.8155260 -0.8428268
f0:session4                  -1.8007971 0.3008318 -2.3904275 -1.2111667
f0:session5                  -1.7028223 0.2828947 -2.2572959 -1.1483487
f0:session6                  -2.0558175 0.3328917 -2.7082852 -1.4033499

The beta table hands you every parameter estimate on the logit scale, along with standard errors and confidence intervals. The row names line up straight with the model equations: S:(Intercept) is $\beta_0$, S:body_cond is $\beta_1$, p:(Intercept) is $\alpha_0$, and p:temp is $\alpha_1$.

Visualising the results

Every predicted-relationship figure in this tutorial runs on the same four steps. Learn them once here and you can draw any parameter, from any model, for the rest of your career.

Step one: pull the coefficients out of the beta table. The table gives you an estimated intercept and slope for each parameter, on the logit scale. For a survival model with one covariate that’s the intercept ($\beta_0$) and the slope ($\beta_1$). For detection it’s the same pair, $\alpha_0$ and $\alpha_1$. Grab the standard errors for both while you’re there, you’ll want them in step three.

Step two: build a covariate sequence and predict on the logit scale. Make a vector of evenly spaced values covering the range of your covariate. For each value $x$, work out $\beta_0 + \beta_1 \times x$. That’s your prediction on the logit scale, right across the range.

Step three: get the uncertainty, using the delta method. You can’t just add and subtract 1.96 standard errors from the prediction, because the prediction is stitched together from two uncertain things at once, the intercept and the slope. The delta method deals with both. On the logit scale the variance of $\beta_0 + \beta_1 \times x$ is:

$\text{Var} = \text{SE}_{\beta_0}^2 + x^2 \times \text{SE}_{\beta_1}^2 + 2x \times \text{Cov}(\beta_0, \beta_1)$

That covariance term $\text{Cov}(\beta_0, \beta_1)$ comes out of the model’s variance-covariance matrix (rd_fit$results$beta.vcv), indexed by where the two parameters sit in the beta table.

Step four: squash everything back with plogis(). The prediction and both interval bounds are still on the logit scale, so run them through plogis() to turn them into probabilities before you plot. The lower bound is plogis(prediction - 1.96 * sqrt(variance)) and the upper is plogis(prediction + 1.96 * sqrt(variance)).

When you do this on your own data, the only things that change are the covariate name, where the parameters sit in the beta table, and the range of $x$ you predict over. The shape of the code stays exactly the same every single time.

Survival as a function of body condition

Code

beta <- rd_fit$results$beta
vcv  <- rd_fit$results$beta.vcv

phi_int_est  <- beta["S:(Intercept)",    "estimate"]
phi_cond_est <- beta["S:body_cond", "estimate"]
phi_int_se   <- beta["S:(Intercept)",    "se"]
phi_cond_se  <- beta["S:body_cond", "se"]

int_idx  <- 1
cond_idx <- 2

cond_seq <- seq(min(cond_seen), max(cond_seen), length.out = 100)

# Predicted survival on logit scale
logit_phi <- phi_int_est + phi_cond_est * cond_seq

# Delta method variance on logit scale
var_logit_phi <- phi_int_se^2 +
                 cond_seq^2 * phi_cond_se^2 +
                 2 * cond_seq * vcv[int_idx, cond_idx]

phi_df <- data.frame(
  condition = cond_seq,
  phi       = plogis(logit_phi),
  lwr       = plogis(logit_phi - 1.96 * sqrt(abs(var_logit_phi))),
  upr       = plogis(logit_phi + 1.96 * sqrt(abs(var_logit_phi)))
)

ggplot(phi_df, aes(x = condition)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#d19527", alpha = 0.2) +
  geom_line(aes(y = phi), colour = "#d19527", linewidth = 1) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(
    x        = "Body condition (standardised)",
    y        = expression(hat(phi))
  ) +
  theme_dark_site()

That confidence band comes from the delta method, carrying the uncertainty in both the intercept and the slope through the inverse-logit squash. What you’re looking for is whether the band would clear a flat horizontal line: if the lower edge sits clearly above, or the upper edge clearly below, some reference probability, that’s a sign the effect is real and not just noise.

Detection as a function of temperature

Code

p_int_est  <- beta["p:(Intercept)", "estimate"]
p_temp_est <- beta["p:temp",        "estimate"]
p_int_se   <- beta["p:(Intercept)", "se"]
p_temp_se  <- beta["p:temp",        "se"]

p_int_idx  <- 5
p_temp_idx <- 6

temp_seq <- seq(-3, 3, length.out = 100)

logit_p <- p_int_est + p_temp_est * temp_seq

var_logit_p <- p_int_se^2 +
               temp_seq^2 * p_temp_se^2 +
               2 * temp_seq * vcv[p_int_idx, p_temp_idx]

p_df <- data.frame(
  temp = temp_seq,
  p    = plogis(logit_p),
  lwr  = plogis(logit_p - 1.96 * sqrt(abs(var_logit_p))),
  upr  = plogis(logit_p + 1.96 * sqrt(abs(var_logit_p)))
)

ggplot(p_df, aes(x = temp)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#FF5733", alpha = 0.2) +
  geom_line(aes(y = p), colour = "#FF5733", linewidth = 1) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(
    x        = "Temperature (standardised)",
    y        = expression(hat(p))
  ) +
  theme_dark_site()

Population size estimates across primary periods

The per-period population size $\hat{N}_t$ is one of the headline numbers a robust design study exists to produce. RMark builds it from the f0 parameter, which is its estimate of how many animals were there but never once detected in a given primary period. Add that to the number you actually counted and you’ve got $\hat{N}_t$.

Code

f0_rows <- beta[grepl("^f0", rownames(beta)), ]

# Number of distinct individuals detected in each primary period
n_detected <- sapply(1:n_primary, function(t) {
  cols <- ((t - 1) * n_secondary + 1):(t * n_secondary)
  sum(rowSums(det_seen[, cols]) > 0)
})

# f0 is on log scale in RMark
f0_est <- exp(f0_rows[, "estimate"])
f0_lcl <- exp(f0_rows[, "lcl"])
f0_ucl <- exp(f0_rows[, "ucl"])

N_df <- data.frame(
  period = 1:n_primary,
  N_hat  = n_detected + f0_est,
  N_lcl  = n_detected + f0_lcl,
  N_ucl  = n_detected + f0_ucl
)

ggplot(N_df, aes(x = period, y = N_hat)) +
  geom_errorbar(aes(ymin = N_lcl, ymax = N_ucl),
                width = 0.2, colour = "#d19527", linewidth = 0.8) +
  geom_point(size = 3, colour = "#d19527") +
  geom_line(colour = "#d19527", linewidth = 0.6, linetype = "dashed") +
  scale_x_continuous(breaks = 1:n_primary) +
  labs(
    x        = "Primary period",
    y        = expression(hat(N))
  ) +
  theme_dark_site()

This plot is your population trend across the study, already corrected for imperfect detection. If it slopes downward, that’s not just you getting worse at spotting animals as the years go on. Because $p$ is estimated on its own, the $\hat{N}_t$ values track genuine changes in abundance. For most conservation work, this is the number everything else was building towards.

A note on your real data

Everything above ran on simulated data where we already knew the true answers. Once you’re working with your own data, a handful of practical things will bite you if you’re not careful.

Standardise your continuous covariates. Subtract the mean, divide by the standard deviation, before you fit. Two reasons. MARK gets numerically twitchy when covariates are on wildly different scales and can just fall over. And standardised coefficients are comparable in size, which really matters at write-up time: a coefficient of 0.3 on a standardised covariate is something you can talk about, whereas a coefficient of 0.003 on rainfall-in-millimetres is a number nobody, you included, will be able to make sense of.

Check for missing covariate values. RMark does not cope gracefully with NA in covariates. Either impute them or drop those individuals before you call process.data().

Get your time.intervals matching your data exactly. The single most common source of baffling RMark errors is the length of the capture history strings not matching the length of time.intervals. Use the build_time_intervals() helper from above instead of hand-typing the vector and hoping.

Keep your raw data and your RMark data frame apart. You’ll almost certainly want the original back later for plotting, checking, or sharing, so don’t overwrite it in place.

Implementing the Robust Design in RMark

Putting it all together

What should influence your parameters?

What influences survival (\(\phi\))?

What influences detection (\(p\))?

What influences availability (\(\gamma'\) and \(\gamma''\))?

The structure of covariates in RMark

How should your data be organised?

The capture history string

Setting up `time.intervals`

Writing out the model

Fitting the model

Visualising the results

Survival as a function of body condition

Detection as a function of temperature

Population size estimates across primary periods

A note on your real data

Putting it all together

What should influence your parameters?

What influences survival (\(\phi\))?

What influences detection (\(p\))?

What influences availability (\(\gamma'\) and \(\gamma''\))?

The structure of covariates in RMark

How should your data be organised?

The capture history string

Setting up time.intervals

Writing out the model

Fitting the model

Visualising the results

Survival as a function of body condition

Detection as a function of temperature

Population size estimates across primary periods

A note on your real data

Setting up `time.intervals`