Checking a stable isotope mixing model

stable isotopes

Bayesian

diet

model checking

ecology tutorial

Three checks for a stable isotope mixing model in R: is the consumer inside the mixing polygon, is the posterior data-driven, and does source aggregation hold.

Author

Tidy Ecology

Published

2026-07-20

A stable isotope mixing model returns diet fractions for any input you give it. That is the problem: it returns them whether or not the input makes sense, whether or not the data inform the answer, and whether or not the quantities you are reading are the ones the data can actually resolve. A confident-looking posterior is not evidence that the fit means anything. This post collects three checks that separate a fit worth reporting from one that only looks like an answer, using the same four-source example as the rest of the mixing-model cluster.

Check one: is the consumer inside the mixing polygon?

A non-negative mixture of the sources can only land inside their convex hull. If the discrimination-corrected consumer falls outside that polygon, no valid diet exists, and yet the sampler will still return a posterior, often a narrow and confident one. The check is geometric and takes a moment; skipping it is how impossible fits get reported as findings.

inpoly <- function(pt, poly) {
  n <- nrow(poly); inside <- FALSE; j <- n
  for (i in 1:n) {
    xi <- poly[i, 1]; yi <- poly[i, 2]; xj <- poly[j, 1]; yj <- poly[j, 2]
    if (((yi > pt[2]) != (yj > pt[2])) &&
        (pt[1] < (xj - xi) * (pt[2] - yi) / (yj - yi) + xi)) inside <- !inside
    j <- i
  }
  inside
}
poly <- Smu[chull(Smu), ]
consumer_in  <- c(d13C = -24.3, d15N = 7.6)              # a plausible consumer
consumer_out <- c(d13C = -33,   d15N = 8)                # more depleted than any source
c(inside = inpoly(consumer_in, poly), outside = inpoly(consumer_out, poly))

 inside outside 
   TRUE   FALSE

The second consumer sits at d13C = -33, more depleted than every source, so it is outside the polygon: no non-negative diet can reach it. Fit it anyway and read what comes back.

set.seed(701)
y_out <- cbind(rnorm(20, consumer_out["d13C"], res[1]), rnorm(20, consumer_out["d15N"], res[2]))
fit_out <- mix_mh(y_out, Smu, Ssd, res, seed = 702)
do <- fit_out$draws; colnames(do) <- rownames(Smu)
spom_mean <- round(mean(do[, "SPOM"]), 3)
spom_sd <- round(sd(do[, "SPOM"]), 3)
round(rbind(mean = colMeans(do), sd = apply(do, 2, sd)), 3)

        C3    C4 Algae  SPOM
mean 0.032 0.001 0.001 0.966
sd   0.024 0.001 0.001 0.024

The fit is emphatic: it assigns 0.966 of the diet to SPOM with a posterior standard deviation of only 0.024. Taken at face value that reads as a precise, near-certain result. It is an artefact. The sampler cannot reach the consumer with any real mixture, so it collapses onto the single source nearest the consumer in the isotope plane and reports false confidence. Without the polygon check, you would report “SPOM dominates the diet” from a consumer that no combination of these sources can produce. A concentrated posterior on an infeasible consumer is a warning, not a result.

A d13C versus d15N biplot with four sources forming a polygon. One consumer sits inside; a red consumer sits to the left, clearly outside the polygon. — Figure 1: The mixing polygon with a feasible consumer (inside) and an infeasible one (outside, in red). The outside consumer is more depleted in carbon than any source, so no non-negative mixture reaches it. A mixing model fitted to it returns a confident but meaningless posterior.

Check two: is the posterior driven by the data or the prior?

A broad posterior might mean the data are weakly informative, or it might mean you are mostly seeing the prior. The flat Dirichlet is not flat on the individual fractions: for four sources, each fraction has a Beta(1, 3) prior with a standard deviation near 0.19. A useful diagnostic is how much the data shrink that: compare the prior standard deviation to the posterior standard deviation. A quantity the data barely touch stays close to its prior width; a quantity the data pin down shrinks a lot.

p_true <- c(0.40, 0.10, 0.20, 0.30); mu_in <- as.numeric(t(Smu) %*% p_true)
set.seed(303)
y_in <- cbind(rnorm(20, mu_in[1], res[1]), rnorm(20, mu_in[2], res[2]))
fit_in <- mix_mh(y_in, Smu, Ssd, res, seed = 404)
d <- fit_in$draws; colnames(d) <- rownames(Smu)
J <- 4
prior_sd_ind <- sqrt((J - 1) / (J^2 * (J + 1)))          # Beta(1, J-1): single source
prior_sd_agg <- sqrt(2 * (J - 2) / (J^2 * (J + 1)))      # Beta(2, J-2): sum of two sources
g <- d[, "C3"] + d[, "SPOM"]                             # the depleted group
shrink_ind <- round(prior_sd_ind / mean(apply(d, 2, sd)), 1)
shrink_agg <- round(prior_sd_agg / sd(g), 1)
c(prior_sd_individual = round(prior_sd_ind, 3), posterior_sd_individual = round(mean(apply(d, 2, sd)), 3),
  shrink_individual = shrink_ind)

    prior_sd_individual posterior_sd_individual       shrink_individual 
                  0.194                   0.062                   3.100

c(prior_sd_aggregate = round(prior_sd_agg, 3), posterior_sd_aggregate = round(sd(g), 3),
  shrink_aggregate = shrink_agg)

    prior_sd_aggregate posterior_sd_aggregate       shrink_aggregate 
                 0.224                  0.017                 13.100

The individual fractions shrink by about 3.1 times relative to their prior; the depleted group shrinks by 13.1 times. The data inform the group decisively and the individual sources only modestly. That is the quantitative version of the resolution problem: the numbers you can trust are the ones the data move far from the prior, and here that is the aggregate, not the individual sources. If a fraction’s posterior barely shrinks from its prior, you are reporting an assumption, not a measurement.

Two panels. Left: a broad dashed prior curve and a somewhat narrower solid posterior for one source. Right: a broad dashed prior and a very narrow solid posterior for the group, near the true value. — Figure 2: Prior (dashed) against posterior (solid) for a single source and for the depleted group. The individual source posterior is only modestly narrower than its broad prior. The group posterior is far tighter than its prior and sits near the true value, so the data, not the prior, are speaking.

Check three: is the resolvable quantity stable across aggregation choices?

The Bayesian post showed that summing the posterior draws for two similar sources gives a well-determined group contribution. A quick stability check is to reach the same quantity by a different route: instead of summing draws afterwards (a posteriori), merge the two sources into one before fitting (a priori) and refit with three sources. If the group contribution is real, the two routes should agree.

Smu3 <- rbind(Depleted = colMeans(Smu[c("C3", "SPOM"), ]), C4 = Smu["C4", ], Algae = Smu["Algae", ])
Ssd3 <- rbind(Depleted = sqrt(colMeans(Ssd[c("C3", "SPOM"), ]^2)), C4 = Ssd["C4", ], Algae = Ssd["Algae", ])
fit3 <- mix_mh(y_in, Smu3, Ssd3, res, alpha = rep(1, 3), seed = 808)
d3 <- fit3$draws; colnames(d3) <- rownames(Smu3)
post_hoc <- round(c(mean = mean(g), sd = sd(g)), 3)
a_priori <- round(c(mean = mean(d3[, "Depleted"]), sd = sd(d3[, "Depleted"])), 3)
rbind(a_posteriori = post_hoc, a_priori = a_priori, truth = c(mean = 0.70, sd = NA))

              mean    sd
a_posteriori 0.703 0.017
a_priori     0.703 0.014
truth        0.700    NA

Both routes put the depleted contribution at a posterior mean of about 0.70 with nearly the same spread, and both bracket the true value. The group contribution does not depend on whether you merge the sources before or after fitting, which is what makes it safe to report. Contrast that with the individual C3 and SPOM fractions, which change completely depending on arbitrary details of the fit: those are the quantities to leave out of the conclusions.

Putting the checks together

Run all three before believing a mixing model. Is the consumer inside the mixing polygon, so a valid diet exists at all? Does the posterior for each reported quantity shrink meaningfully from its prior, so you are reading data rather than assumption? And is the quantity stable to how sources are grouped, so it is not an artefact of one arbitrary aggregation? A fit that passes all three is worth reporting; a fit that fails any of them is telling you to fix the study design, add tracers, or report a coarser quantity, not to trust the numbers as they stand.

Honest limits

These are minimal, transparent versions of checks that mature tools formalise. The polygon test here uses source means; a fuller version simulates the mixing region accounting for source and discrimination uncertainty, which can shrink the feasible area. The prior-shrinkage diagnostic is a rule of thumb, not a formal model comparison. And agreement between the two aggregation routes confirms a quantity is stable, not that it is unbiased, since a wrong discrimination factor or a missing source would bias both routes together. The checks catch the common, fatal mistakes: an impossible consumer, a prior masquerading as a result, and an unresolvable fraction read as if it were pinned down. They do not substitute for a good sampling design and honest sources.

References

Phillips DL, Gregg JW 2003 Oecologia 136:261-269 (10.1007/s00442-003-1218-3).

Phillips DL, Inger R, Bearhop S, Jackson AL, Moore JW, Parnell AC, Semmens BX, Ward EJ 2014 Canadian Journal of Zoology 92:823-835 (10.1139/cjz-2014-0127).

Stock BC, Jackson AL, Ward EJ, Parnell AC, Phillips DL, Semmens BX 2018 PeerJ 6:e5096 (10.7717/peerj.5096).

Moore JW, Semmens BX 2008 Ecology Letters 11:470-480.

--- title: "Checking a stable isotope mixing model" description: "Three checks for a stable isotope mixing model in R: is the consumer inside the mixing polygon, is the posterior data-driven, and does source aggregation hold." date: "2026-07-20 11:00" categories: [R, stable isotopes, Bayesian, diet, model checking, ecology tutorial] image: thumbnail.png image-alt: "A mixing polygon with a consumer point sitting clearly outside it, flagged as infeasible." --- A stable isotope mixing model returns diet fractions for any input you give it. That is the problem: it returns them whether or not the input makes sense, whether or not the data inform the answer, and whether or not the quantities you are reading are the ones the data can actually resolve. A confident-looking posterior is not evidence that the fit means anything. This post collects three checks that separate a fit worth reporting from one that only looks like an answer, using the same four-source example as the rest of the [mixing-model cluster](../isotope-mixing-models-in-r/). ```{r} #| label: setup #| echo: false #| warning: false #| message: false library(ggplot2) theme_te <- function() { theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold", size = 13), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA), panel.grid.minor = element_blank()) } te_green <- "#275139" te_pal <- c("#275139", "#8c6d1f", "#7a1f1f", "#1f5c7a") ldir <- function(x, a) lgamma(sum(a)) - sum(lgamma(a)) + sum((a - 1) * log(pmax(x, 1e-12))) rdir <- function(a) { g <- pmax(rgamma(length(a), a, 1), 1e-300); g / sum(g) } mix_mh <- function(y, Smu, Ssd, res, alpha = rep(1, nrow(Smu)), kappa = 600, iter = 60000, burn = 12000, seed = 1) { n <- nrow(y); J <- nrow(Smu) loglik <- function(p) { mu <- as.numeric(t(Smu) %*% p); v <- as.numeric(t(Ssd^2) %*% (p^2)) + res^2 sum(dnorm(y, rep(mu, each = n), rep(sqrt(v), each = n), log = TRUE)) } set.seed(seed); p <- rep(1 / J, J); ll <- loglik(p) + ldir(p, alpha) keep <- matrix(NA_real_, iter, J); acc <- 0 for (t in 1:iter) { ap <- pmax(kappa * p, 1e-3); ps <- rdir(ap); aps <- pmax(kappa * ps, 1e-3) lls <- loglik(ps) + ldir(ps, alpha) if (log(runif(1)) < (lls - ll) + ldir(p, aps) - ldir(ps, ap)) { p <- ps; ll <- lls; acc <- acc + 1 } keep[t, ] <- p } list(draws = keep[(burn + 1):iter, , drop = FALSE], acc = acc / iter) } Smu <- rbind(C3 = c(-29, 4), C4 = c(-13, 6), Algae = c(-15, 12), SPOM = c(-28, 10)) Ssd <- rbind(C3 = c(1.0, 0.8), C4 = c(1.2, 1.0), Algae = c(0.9, 1.1), SPOM = c(1.1, 0.9)) colnames(Smu) <- colnames(Ssd) <- c("d13C", "d15N"); res <- c(0.5, 0.4) ``` ## Check one: is the consumer inside the mixing polygon? A non-negative mixture of the sources can only land inside their convex hull. If the discrimination-corrected consumer falls outside that polygon, no valid diet exists, and yet the sampler will still return a posterior, often a narrow and confident one. The check is geometric and takes a moment; skipping it is how impossible fits get reported as findings. ```{r} #| label: feasibility inpoly <- function(pt, poly) { n <- nrow(poly); inside <- FALSE; j <- n for (i in 1:n) { xi <- poly[i, 1]; yi <- poly[i, 2]; xj <- poly[j, 1]; yj <- poly[j, 2] if (((yi > pt[2]) != (yj > pt[2])) && (pt[1] < (xj - xi) * (pt[2] - yi) / (yj - yi) + xi)) inside <- !inside j <- i } inside } poly <- Smu[chull(Smu), ] consumer_in <- c(d13C = -24.3, d15N = 7.6) # a plausible consumer consumer_out <- c(d13C = -33, d15N = 8) # more depleted than any source c(inside = inpoly(consumer_in, poly), outside = inpoly(consumer_out, poly)) ``` The second consumer sits at `d13C = -33`, more depleted than every source, so it is outside the polygon: no non-negative diet can reach it. Fit it anyway and read what comes back. ```{r} #| label: outside-fit set.seed(701) y_out <- cbind(rnorm(20, consumer_out["d13C"], res[1]), rnorm(20, consumer_out["d15N"], res[2])) fit_out <- mix_mh(y_out, Smu, Ssd, res, seed = 702) do <- fit_out$draws; colnames(do) <- rownames(Smu) spom_mean <- round(mean(do[, "SPOM"]), 3) spom_sd <- round(sd(do[, "SPOM"]), 3) round(rbind(mean = colMeans(do), sd = apply(do, 2, sd)), 3) ``` The fit is emphatic: it assigns `r spom_mean` of the diet to SPOM with a posterior standard deviation of only `r spom_sd`. Taken at face value that reads as a precise, near-certain result. It is an artefact. The sampler cannot reach the consumer with any real mixture, so it collapses onto the single source nearest the consumer in the isotope plane and reports false confidence. Without the polygon check, you would report "SPOM dominates the diet" from a consumer that no combination of these sources can produce. A concentrated posterior on an infeasible consumer is a warning, not a result. ```{r} #| label: fig-feasibility #| echo: false #| fig-cap: "The mixing polygon with a feasible consumer (inside) and an infeasible one (outside, in red). The outside consumer is more depleted in carbon than any source, so no non-negative mixture reaches it. A mixing model fitted to it returns a confident but meaningless posterior." #| fig-alt: "A d13C versus d15N biplot with four sources forming a polygon. One consumer sits inside; a red consumer sits to the left, clearly outside the polygon." hull_df <- as.data.frame(rbind(poly, poly[1, ])) sd_df <- data.frame(Smu, lab = rownames(Smu)) ggplot() + geom_polygon(data = hull_df, aes(d13C, d15N), fill = te_green, alpha = 0.10, colour = te_green, linewidth = 0.6) + geom_point(data = sd_df, aes(d13C, d15N), colour = te_green, size = 3) + geom_text(data = sd_df, aes(d13C, d15N, label = lab), vjust = -1, size = 3.4) + geom_point(aes(consumer_in["d13C"], consumer_in["d15N"]), shape = 4, size = 4, stroke = 1.2) + annotate("text", x = consumer_in["d13C"], y = consumer_in["d15N"] - 0.8, label = "inside", size = 3.4) + geom_point(aes(consumer_out["d13C"], consumer_out["d15N"]), shape = 4, size = 4, stroke = 1.2, colour = "#7a1f1f") + annotate("text", x = consumer_out["d13C"], y = consumer_out["d15N"] - 0.8, label = "outside", size = 3.4, colour = "#7a1f1f") + labs(x = expression(delta^13 * C~"(permil)"), y = expression(delta^15 * N~"(permil)"), title = "Is the consumer inside the polygon?") + theme_te() ``` ## Check two: is the posterior driven by the data or the prior? A broad posterior might mean the data are weakly informative, or it might mean you are mostly seeing the prior. The flat Dirichlet is not flat on the individual fractions: for four sources, each fraction has a Beta(1, 3) prior with a standard deviation near 0.19. A useful diagnostic is how much the data shrink that: compare the prior standard deviation to the posterior standard deviation. A quantity the data barely touch stays close to its prior width; a quantity the data pin down shrinks a lot. ```{r} #| label: prior-sensitivity p_true <- c(0.40, 0.10, 0.20, 0.30); mu_in <- as.numeric(t(Smu) %*% p_true) set.seed(303) y_in <- cbind(rnorm(20, mu_in[1], res[1]), rnorm(20, mu_in[2], res[2])) fit_in <- mix_mh(y_in, Smu, Ssd, res, seed = 404) d <- fit_in$draws; colnames(d) <- rownames(Smu) J <- 4 prior_sd_ind <- sqrt((J - 1) / (J^2 * (J + 1))) # Beta(1, J-1): single source prior_sd_agg <- sqrt(2 * (J - 2) / (J^2 * (J + 1))) # Beta(2, J-2): sum of two sources g <- d[, "C3"] + d[, "SPOM"] # the depleted group shrink_ind <- round(prior_sd_ind / mean(apply(d, 2, sd)), 1) shrink_agg <- round(prior_sd_agg / sd(g), 1) c(prior_sd_individual = round(prior_sd_ind, 3), posterior_sd_individual = round(mean(apply(d, 2, sd)), 3), shrink_individual = shrink_ind) c(prior_sd_aggregate = round(prior_sd_agg, 3), posterior_sd_aggregate = round(sd(g), 3), shrink_aggregate = shrink_agg) ``` The individual fractions shrink by about `r shrink_ind` times relative to their prior; the depleted group shrinks by `r shrink_agg` times. The data inform the group decisively and the individual sources only modestly. That is the quantitative version of the resolution problem: the numbers you can trust are the ones the data move far from the prior, and here that is the aggregate, not the individual sources. If a fraction's posterior barely shrinks from its prior, you are reporting an assumption, not a measurement. ```{r} #| label: fig-prior #| echo: false #| fig-cap: "Prior (dashed) against posterior (solid) for a single source and for the depleted group. The individual source posterior is only modestly narrower than its broad prior. The group posterior is far tighter than its prior and sits near the true value, so the data, not the prior, are speaking." #| fig-alt: "Two panels. Left: a broad dashed prior curve and a somewhat narrower solid posterior for one source. Right: a broad dashed prior and a very narrow solid posterior for the group, near the true value." xx <- seq(0, 1, length.out = 400) prior_ind <- dbeta(xx, 1, J - 1); prior_agg <- dbeta(xx, 2, J - 2) pf <- rbind( data.frame(x = xx, dens = prior_ind, kind = "prior", panel = "single source (C3)"), data.frame(x = density(d[, "C3"], from = 0, to = 1)$x, dens = density(d[, "C3"], from = 0, to = 1)$y, kind = "posterior", panel = "single source (C3)"), data.frame(x = xx, dens = prior_agg, kind = "prior", panel = "depleted group (C3 + SPOM)"), data.frame(x = density(g, from = 0, to = 1)$x, dens = density(g, from = 0, to = 1)$y, kind = "posterior", panel = "depleted group (C3 + SPOM)")) pf$panel <- factor(pf$panel, levels = c("single source (C3)", "depleted group (C3 + SPOM)")) ggplot(pf, aes(x, dens, colour = kind, linetype = kind)) + geom_line(linewidth = 0.8) + facet_wrap(~panel) + scale_colour_manual(values = c(prior = "grey55", posterior = te_green)) + scale_linetype_manual(values = c(prior = "dashed", posterior = "solid")) + labs(x = "diet fraction", y = "density", colour = NULL, linetype = NULL, title = "How far does the data move each quantity from its prior?") + theme_te() + theme(legend.position = "top") ``` ## Check three: is the resolvable quantity stable across aggregation choices? The [Bayesian post](../bayesian-isotope-mixing-models/) showed that summing the posterior draws for two similar sources gives a well-determined group contribution. A quick stability check is to reach the same quantity by a different route: instead of summing draws afterwards (a posteriori), merge the two sources into one before fitting (a priori) and refit with three sources. If the group contribution is real, the two routes should agree. ```{r} #| label: aggregate-check Smu3 <- rbind(Depleted = colMeans(Smu[c("C3", "SPOM"), ]), C4 = Smu["C4", ], Algae = Smu["Algae", ]) Ssd3 <- rbind(Depleted = sqrt(colMeans(Ssd[c("C3", "SPOM"), ]^2)), C4 = Ssd["C4", ], Algae = Ssd["Algae", ]) fit3 <- mix_mh(y_in, Smu3, Ssd3, res, alpha = rep(1, 3), seed = 808) d3 <- fit3$draws; colnames(d3) <- rownames(Smu3) post_hoc <- round(c(mean = mean(g), sd = sd(g)), 3) a_priori <- round(c(mean = mean(d3[, "Depleted"]), sd = sd(d3[, "Depleted"])), 3) rbind(a_posteriori = post_hoc, a_priori = a_priori, truth = c(mean = 0.70, sd = NA)) ``` Both routes put the depleted contribution at a posterior mean of about 0.70 with nearly the same spread, and both bracket the true value. The group contribution does not depend on whether you merge the sources before or after fitting, which is what makes it safe to report. Contrast that with the individual C3 and SPOM fractions, which change completely depending on arbitrary details of the fit: those are the quantities to leave out of the conclusions. ## Putting the checks together Run all three before believing a mixing model. Is the consumer inside the mixing polygon, so a valid diet exists at all? Does the posterior for each reported quantity shrink meaningfully from its prior, so you are reading data rather than assumption? And is the quantity stable to how sources are grouped, so it is not an artefact of one arbitrary aggregation? A fit that passes all three is worth reporting; a fit that fails any of them is telling you to fix the study design, add tracers, or report a coarser quantity, not to trust the numbers as they stand. ## Honest limits These are minimal, transparent versions of checks that mature tools formalise. The polygon test here uses source means; a fuller version simulates the mixing region accounting for source and discrimination uncertainty, which can shrink the feasible area. The prior-shrinkage diagnostic is a rule of thumb, not a formal model comparison. And agreement between the two aggregation routes confirms a quantity is stable, not that it is unbiased, since a wrong discrimination factor or a missing source would bias both routes together. The checks catch the common, fatal mistakes: an impossible consumer, a prior masquerading as a result, and an unresolvable fraction read as if it were pinned down. They do not substitute for a good sampling design and honest sources. ## Related tutorials - [Stable isotope mixing models in R](../isotope-mixing-models-in-r/) - [Bayesian stable isotope mixing models](../bayesian-isotope-mixing-models/) - [Discrimination and concentration in mixing models](../discrimination-and-concentration-in-mixing/) - [Checking a community occupancy model](../checking-a-community-occupancy-model/) ## References Phillips DL, Gregg JW 2003 Oecologia 136:261-269 (10.1007/s00442-003-1218-3). Phillips DL, Inger R, Bearhop S, Jackson AL, Moore JW, Parnell AC, Semmens BX, Ward EJ 2014 Canadian Journal of Zoology 92:823-835 (10.1139/cjz-2014-0127). Stock BC, Jackson AL, Ward EJ, Parnell AC, Phillips DL, Semmens BX 2018 PeerJ 6:e5096 (10.7717/peerj.5096). Moore JW, Semmens BX 2008 Ecology Letters 11:470-480.