Bayesian stable isotope mixing models

stable isotopes

Bayesian

diet

MCMC

ecology tutorial

Build a Bayesian stable isotope mixing model from scratch in R with a Metropolis sampler, and read the wide, correlated posterior that two tracers produce.

Author

Tidy Ecology

Published

2026-07-19

The geometry post showed the problem: with two tracers and four sources, a consumer’s isotope signature is consistent with a whole band of diets, not one. A point estimate picks an arbitrary member of that band and hides the rest. A Bayesian mixing model does the honest thing instead: it reports the full band as a posterior distribution, so the uncertainty is visible rather than swept into a single number. This is what SIAR, MixSIR, and MixSIAR do. Here we build a small version from scratch, in base R, to see exactly where the width and the correlations come from.

The model

Each source has a mean isotope value and a spread; real sources are not points, they are clouds of measured individuals. The model treats the consumer’s tracer values as normal around the mixture mean, with a variance that combines the (fraction-weighted) source variances and a residual term:

\[ y_t \sim \mathcal{N}\!\left(\textstyle\sum_j p_j s_{jt},\ \ \textstyle\sum_j p_j^2 \sigma_{jt}^2 + \omega_t^2\right). \]

The diet fractions p live on the simplex (non-negative, sum to one), so the natural prior is a Dirichlet. A flat Dirichlet(1, …, 1) is the usual “uninformative” choice, uniform over the simplex. That is the whole model: a normal likelihood for the tracers and a Dirichlet prior for the fractions. Fitting it means sampling the posterior.

A sampler in one function

Because the fractions are constrained to the simplex, an ordinary random walk would keep proposing invalid points. A cleaner move is a Dirichlet proposal centred on the current point: draw p* from a Dirichlet whose mean is the current p, with a concentration kappa that tunes the step size. The proposal is not symmetric, so the Metropolis-Hastings ratio needs the proposal densities in both directions; the flat prior cancels in the ratio. The whole sampler is short.

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)                       # mixture mean per tracer
    v  <- as.numeric(t(Ssd^2) %*% (p^2)) + res^2         # source + residual variance
    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)         # propose from Dirichlet(kappa * p)
    aps <- pmax(kappa * ps, 1e-3)
    lls <- loglik(ps) + ldir(ps, alpha)
    # Hastings ratio: likelihood + prior + reverse/forward proposal densities
    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)
}

The 1e-3 floor on the Dirichlet parameters and the 1e-12 floor inside ldir matter: when a fraction is pushed to the edge of the simplex, log(0) would otherwise produce NaN and stall the chain. Guarding the logarithms keeps the sampler stable near the boundary.

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)                                       # residual SD per tracer

Does the sampler work? A determined test first

Before trusting any sampler on a hard problem, run it on an easy one where the answer is known. Three sources and two tracers is determined (the geometry post solved it exactly), so a consumer generated from a known diet should give a posterior that concentrates near that diet. Generate data from a diet of 40% C3, 30% C4, 30% Algae and fit.

S3 <- Smu[c("C3", "C4", "Algae"), ]; Sd3 <- Ssd[c("C3", "C4", "Algae"), ]
p_true3 <- c(0.40, 0.30, 0.30)
mu3 <- as.numeric(t(S3) %*% p_true3)
set.seed(101)
y3 <- cbind(rnorm(30, mu3[1], res[1]), rnorm(30, mu3[2], res[2]))
fit3 <- mix_mh(y3, S3, Sd3, res, alpha = rep(1, 3), iter = 40000, burn = 8000, seed = 202)
cal_acc <- round(fit3$acc, 3)
cal_ess <- round(min(apply(fit3$draws, 2, ess)))
cal_means <- round(colMeans(fit3$draws), 3)
c(acceptance = cal_acc, min_ESS = cal_ess)

acceptance    min_ESS 
     0.393   1801.000

rbind(posterior_mean = cal_means, truth = p_true3)

                [,1]  [,2]  [,3]
posterior_mean 0.404 0.302 0.294
truth          0.400 0.300 0.300

The acceptance rate is 0.393 and the smallest effective sample size across the three fractions is 1801, so the chain mixes well. The posterior means (0.404 C3, 0.302 C4, 0.294 Algae) sit right on the true diet. The sampler recovers a known answer when the geometry allows it. Now push it into the case that does not.

Three narrow, well-separated posterior density curves for C3, C4 and Algae, each peaking at its true value. — Figure 1: Posterior distributions for the determined three-source problem. Each fraction concentrates tightly around its true value (dashed line). When the geometry pins the answer, the Bayesian model finds it.

The underdetermined posterior

Bring back the fourth source. The consumer is generated from a diet of 40% C3, 10% C4, 20% Algae, 30% SPOM, but now four unknowns face three constraints. The sampler runs fine; the posterior tells the real story.

p_true4 <- c(0.40, 0.10, 0.20, 0.30)
mu4 <- as.numeric(t(Smu) %*% p_true4)
set.seed(303)
y4 <- cbind(rnorm(20, mu4[1], res[1]), rnorm(20, mu4[2], res[2]))
fit4 <- mix_mh(y4, Smu, Ssd, res, iter = 60000, burn = 12000, seed = 404)
d <- fit4$draws
colnames(d) <- rownames(Smu)
u_acc <- round(fit4$acc, 3)
u_ess <- round(min(apply(d, 2, ess)))
post <- t(apply(d, 2, function(z) c(mean = mean(z), sd = sd(z),
                                    lo = quantile(z, 0.025), hi = quantile(z, 0.975))))
rownames(post) <- rownames(Smu)
sd_spom <- round(post["SPOM", "sd"], 3)
sd_c3 <- round(post["C3", "sd"], 3)
round(cbind(post, truth = p_true4), 3)

       mean    sd lo.2.5% hi.97.5% truth
C3    0.366 0.058   0.250    0.472   0.4
C4    0.134 0.054   0.038    0.243   0.1
Algae 0.163 0.065   0.031    0.279   0.2
SPOM  0.337 0.070   0.208    0.476   0.3

The chain is healthy (acceptance 0.464, minimum effective sample size 136), but the posterior is broad. The SPOM fraction has a posterior standard deviation of 0.07 and a 95% credible interval roughly a quarter of the whole simplex wide; C3 is similar. Every one of these fractions has a true value comfortably inside its interval, so the model is not wrong, it is honestly uncertain. This width is the band from the geometry post, now expressed as probability.

Four wide, overlapping posterior density curves for the four sources, each spanning a large part of the zero-to-one range. — Figure 2: Posterior distributions for the four-source problem. Each fraction is broad and its interval covers the true value (dashed line), but no fraction is pinned down. The width is the geometric ambiguity of two tracers and four sources, made explicit.

A single point estimate would hide all of this. The minimum-norm linear solution, for instance, returns one plausible-looking diet:

M <- rbind(t(Smu), 1)
b <- c(colMeans(y4), 1)
p_point <- as.numeric(t(M) %*% solve(M %*% t(M), b))     # minimum-norm solution
names(p_point) <- rownames(Smu)
round(p_point, 3)

   C3    C4 Algae  SPOM 
0.362 0.136 0.159 0.343

Those four numbers look like an answer, but they are one arbitrary point inside the posterior cloud. Reporting them alone would imply a precision the data do not support.

Why the fractions trade off

The fractions are not just individually uncertain, they are locked together. Look at the posterior correlations.

C <- cor(d); dimnames(C) <- list(rownames(Smu), rownames(Smu))
cor_sc <- round(C["SPOM", "C3"], 2)
round(C, 2)

         C3    C4 Algae  SPOM
C3     1.00 -0.91  0.92 -0.98
C4    -0.91  1.00 -0.98  0.90
Algae  0.92 -0.98  1.00 -0.94
SPOM  -0.98  0.90 -0.94  1.00

Every pair is strongly correlated, and the strongest is SPOM against C3 at -0.98. The reason is in the mixing space: C3 and SPOM both sit at the depleted end of the carbon axis, so swapping one for the other barely moves the mixture. The data constrain the sum of the depleted sources tightly while leaving the split between them almost free. That near-one-to-one exchange is the correlation, and it is the fingerprint of two sources the tracers cannot separate.

A dense scatter of posterior draws forming a narrow downward-sloping band, showing a strong negative correlation between the SPOM and C3 fractions. — Figure 3: Posterior draws of the SPOM fraction against the C3 fraction. The tight negative ridge means the data fix their sum but not their individual values: whenever one goes up, the other comes down by nearly the same amount.

What the data can actually resolve

If the tracers fix the sum of the two depleted sources, then that sum is the quantity the data determine, even though neither source alone is. Add the posterior draws for C3 and SPOM (and, symmetrically, for the two enriched sources C4 and Algae) and the picture sharpens dramatically.

g_depleted <- d[, "C3"] + d[, "SPOM"]                    # true sum: 0.70
g_enriched <- d[, "C4"] + d[, "Algae"]                   # true sum: 0.30
summ <- function(z) c(mean(z), sd(z), as.numeric(quantile(z, c(0.025, 0.975))))
agg <- rbind(summ(g_depleted), summ(g_enriched))
dimnames(agg) <- list(c("C3 + SPOM", "C4 + Algae"), c("mean", "sd", "lo", "hi"))
sd_ind <- round(mean(apply(d, 2, sd)), 3)
sd_agg <- round(mean(c(sd(g_depleted), sd(g_enriched))), 3)
gain <- round(sd_ind / sd_agg, 1)
round(agg, 3)

            mean    sd    lo    hi
C3 + SPOM  0.703 0.017 0.671 0.737
C4 + Algae 0.297 0.017 0.263 0.329

c(individual_sd = sd_ind, aggregate_sd = sd_agg, narrower_by = gain)

individual_sd  aggregate_sd   narrower_by 
        0.062         0.017         3.600

The depleted group sums to a posterior mean of 0.703 with a standard deviation of just 0.017, and its credible interval brackets the true 0.70 tightly. The individual sources had an average posterior standard deviation of 0.062; the grouped sum is 3.6 times narrower. The strong negative correlation that made the individuals ambiguous is exactly what makes their sum precise: the errors cancel. A mixing model that cannot tell C3 from SPOM can still tell you, with confidence, how much of the diet came from depleted sources as a whole.

Two broad pale density curves for C3 and SPOM overlaid with one narrow dark density curve for their sum, which is much more concentrated. — Figure 4: The two individual depleted sources (wide, pale) against their sum (narrow, dark). Neither source alone is resolvable, but the quantity that matters ecologically, the total depleted contribution, is pinned down. Aggregating similar sources turns an unanswerable question into an answerable one.

Does more data help?

If the width came from sampling noise, more consumer measurements would shrink it. Some of it does: raise the consumer sample from 20 to 200 and refit.

set.seed(505)
y4b <- cbind(rnorm(200, mu4[1], res[1]), rnorm(200, mu4[2], res[2]))
fit4b <- mix_mh(y4b, Smu, Ssd, res, iter = 60000, burn = 12000, seed = 606)
sd_small <- round(mean(apply(d, 2, sd)), 3)
sd_large <- round(mean(apply(fit4b$draws, 2, sd)), 3)
ratio <- round(sd_large / sd_small, 2)
c(sd_n20 = sd_small, sd_n200 = sd_large, ratio = ratio)

 sd_n20 sd_n200   ratio 
  0.062   0.025   0.400

Ten times as much data shrinks the individual standard deviations to 0.4 of their former size, short of the reduction to about a third that pure sampling noise would give (the square root of ten). Part of the uncertainty shrinks with data; part of it is the geometry, and no consumer sample size removes that. The aggregate was already tight at twenty measurements. The lesson is to spend effort where it pays: more tracers to break the geometric ambiguity, or grouping to report what is identifiable, rather than more consumer replicates chasing a limit that data cannot cross.

Where to go next

This posterior is only as trustworthy as its inputs. The source values here were corrected to the consumer’s trophic level before any fitting, and that correction, the trophic discrimination factor, turns out to move the answer more than the data do. Before reporting a fit, run the standard checks: is the consumer inside the mixing polygon, is the posterior actually informed by the data rather than the prior, and which contributions are resolvable as groups.

Honest limits

The sampler here is deliberately minimal: one consumer, fixed source means and spreads, a flat prior, and a hand-rolled Metropolis step. Production tools add hierarchical structure, covariates, informative priors, and better-mixing samplers, and you should use them for real analyses (SIAR, MixSIAR, and their relatives). What the toy version buys is transparency: the wide posterior and the strong correlations are not artefacts of a black box, they are the honest consequence of asking two tracers to resolve four sources. A fancier sampler would report the same width. The fix is more information, not more computation.

References

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

Parnell AC, Inger R, Bearhop S, Jackson AL 2010 PLoS ONE 5:e9672 (10.1371/journal.pone.0009672).

Parnell AC, Phillips DL, Bearhop S, Semmens BX, Ward EJ, Moore JW, Jackson AL, Grey J, Kelley DJ, Inger R 2013 Environmetrics 24:387-399.

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

--- title: "Bayesian stable isotope mixing models" description: "Build a Bayesian stable isotope mixing model from scratch in R with a Metropolis sampler, and read the wide, correlated posterior that two tracers produce." date: "2026-07-19 13:00" categories: [R, stable isotopes, Bayesian, diet, MCMC, ecology tutorial] image: thumbnail.png image-alt: "Posterior histograms for four diet sources, wide and overlapping, with vertical lines marking the true values." --- The [geometry post](../isotope-mixing-models-in-r/) showed the problem: with two tracers and four sources, a consumer's isotope signature is consistent with a whole band of diets, not one. A point estimate picks an arbitrary member of that band and hides the rest. A Bayesian mixing model does the honest thing instead: it reports the full band as a posterior distribution, so the uncertainty is visible rather than swept into a single number. This is what SIAR, MixSIR, and MixSIAR do. Here we build a small version from scratch, in base R, to see exactly where the width and the correlations come from. ## The model Each source has a mean isotope value and a spread; real sources are not points, they are clouds of measured individuals. The model treats the consumer's tracer values as normal around the mixture mean, with a variance that combines the (fraction-weighted) source variances and a residual term: $$ y_t \sim \mathcal{N}\!\left(\textstyle\sum_j p_j s_{jt},\ \ \textstyle\sum_j p_j^2 \sigma_{jt}^2 + \omega_t^2\right). $$ The diet fractions `p` live on the simplex (non-negative, sum to one), so the natural prior is a Dirichlet. A flat Dirichlet(1, ..., 1) is the usual "uninformative" choice, uniform over the simplex. That is the whole model: a normal likelihood for the tracers and a Dirichlet prior for the fractions. Fitting it means sampling the posterior. ```{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") # Dirichlet density (log) and sampler; a small effective-sample-size helper. 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) } ess <- function(x) { x <- x - mean(x); n <- length(x) r <- acf(x, lag.max = min(2000, n - 1), plot = FALSE)$acf[-1] k <- which(r < 0.05)[1]; if (is.na(k)) k <- length(r) rho <- r[1:max(1, k - 1)]; n / (1 + 2 * sum(rho)) } ``` ## A sampler in one function Because the fractions are constrained to the simplex, an ordinary random walk would keep proposing invalid points. A cleaner move is a Dirichlet proposal centred on the current point: draw `p*` from a Dirichlet whose mean is the current `p`, with a concentration `kappa` that tunes the step size. The proposal is not symmetric, so the Metropolis-Hastings ratio needs the proposal densities in both directions; the flat prior cancels in the ratio. The whole sampler is short. ```{r} #| label: sampler 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) # mixture mean per tracer v <- as.numeric(t(Ssd^2) %*% (p^2)) + res^2 # source + residual variance 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) # propose from Dirichlet(kappa * p) aps <- pmax(kappa * ps, 1e-3) lls <- loglik(ps) + ldir(ps, alpha) # Hastings ratio: likelihood + prior + reverse/forward proposal densities 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) } ``` The `1e-3` floor on the Dirichlet parameters and the `1e-12` floor inside `ldir` matter: when a fraction is pushed to the edge of the simplex, `log(0)` would otherwise produce `NaN` and stall the chain. Guarding the logarithms keeps the sampler stable near the boundary. ```{r} #| label: sources 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) # residual SD per tracer ``` ## Does the sampler work? A determined test first Before trusting any sampler on a hard problem, run it on an easy one where the answer is known. Three sources and two tracers is determined (the geometry post solved it exactly), so a consumer generated from a known diet should give a posterior that concentrates near that diet. Generate data from a diet of 40% C3, 30% C4, 30% Algae and fit. ```{r} #| label: calibration S3 <- Smu[c("C3", "C4", "Algae"), ]; Sd3 <- Ssd[c("C3", "C4", "Algae"), ] p_true3 <- c(0.40, 0.30, 0.30) mu3 <- as.numeric(t(S3) %*% p_true3) set.seed(101) y3 <- cbind(rnorm(30, mu3[1], res[1]), rnorm(30, mu3[2], res[2])) fit3 <- mix_mh(y3, S3, Sd3, res, alpha = rep(1, 3), iter = 40000, burn = 8000, seed = 202) cal_acc <- round(fit3$acc, 3) cal_ess <- round(min(apply(fit3$draws, 2, ess))) cal_means <- round(colMeans(fit3$draws), 3) c(acceptance = cal_acc, min_ESS = cal_ess) rbind(posterior_mean = cal_means, truth = p_true3) ``` The acceptance rate is `r cal_acc` and the smallest effective sample size across the three fractions is `r cal_ess`, so the chain mixes well. The posterior means (`r cal_means[1]` C3, `r cal_means[2]` C4, `r cal_means[3]` Algae) sit right on the true diet. The sampler recovers a known answer when the geometry allows it. Now push it into the case that does not. ```{r} #| label: fig-calibration #| echo: false #| fig-cap: "Posterior distributions for the determined three-source problem. Each fraction concentrates tightly around its true value (dashed line). When the geometry pins the answer, the Bayesian model finds it." #| fig-alt: "Three narrow, well-separated posterior density curves for C3, C4 and Algae, each peaking at its true value." cal_df <- data.frame(frac = c(fit3$draws), src = factor(rep(rownames(S3), each = nrow(fit3$draws)), levels = rownames(S3))) tr_df <- data.frame(src = factor(rownames(S3), levels = rownames(S3)), t = p_true3) ggplot(cal_df, aes(frac, fill = src, colour = src)) + geom_density(alpha = 0.35, linewidth = 0.5) + geom_vline(data = tr_df, aes(xintercept = t, colour = src), linetype = "dashed", linewidth = 0.6, show.legend = FALSE) + scale_fill_manual(values = te_pal) + scale_colour_manual(values = te_pal) + labs(x = "diet fraction", y = "posterior density", fill = NULL, colour = NULL, title = "Determined problem: the posterior finds the truth") + theme_te() + theme(legend.position = "top") ``` ## The underdetermined posterior Bring back the fourth source. The consumer is generated from a diet of 40% C3, 10% C4, 20% Algae, 30% SPOM, but now four unknowns face three constraints. The sampler runs fine; the posterior tells the real story. ```{r} #| label: underdetermined p_true4 <- c(0.40, 0.10, 0.20, 0.30) mu4 <- as.numeric(t(Smu) %*% p_true4) set.seed(303) y4 <- cbind(rnorm(20, mu4[1], res[1]), rnorm(20, mu4[2], res[2])) fit4 <- mix_mh(y4, Smu, Ssd, res, iter = 60000, burn = 12000, seed = 404) d <- fit4$draws colnames(d) <- rownames(Smu) u_acc <- round(fit4$acc, 3) u_ess <- round(min(apply(d, 2, ess))) post <- t(apply(d, 2, function(z) c(mean = mean(z), sd = sd(z), lo = quantile(z, 0.025), hi = quantile(z, 0.975)))) rownames(post) <- rownames(Smu) sd_spom <- round(post["SPOM", "sd"], 3) sd_c3 <- round(post["C3", "sd"], 3) round(cbind(post, truth = p_true4), 3) ``` The chain is healthy (acceptance `r u_acc`, minimum effective sample size `r u_ess`), but the posterior is broad. The SPOM fraction has a posterior standard deviation of `r sd_spom` and a 95% credible interval roughly a quarter of the whole simplex wide; C3 is similar. Every one of these fractions has a true value comfortably inside its interval, so the model is not wrong, it is honestly uncertain. This width is the band from the geometry post, now expressed as probability. ```{r} #| label: fig-underdetermined #| echo: false #| fig-cap: "Posterior distributions for the four-source problem. Each fraction is broad and its interval covers the true value (dashed line), but no fraction is pinned down. The width is the geometric ambiguity of two tracers and four sources, made explicit." #| fig-alt: "Four wide, overlapping posterior density curves for the four sources, each spanning a large part of the zero-to-one range." u_df <- data.frame(frac = c(d), src = factor(rep(rownames(Smu), each = nrow(d)), levels = rownames(Smu))) t4_df <- data.frame(src = factor(rownames(Smu), levels = rownames(Smu)), t = p_true4) ggplot(u_df, aes(frac, fill = src, colour = src)) + geom_density(alpha = 0.30, linewidth = 0.5) + geom_vline(data = t4_df, aes(xintercept = t, colour = src), linetype = "dashed", linewidth = 0.6, show.legend = FALSE) + scale_fill_manual(values = te_pal) + scale_colour_manual(values = te_pal) + labs(x = "diet fraction", y = "posterior density", fill = NULL, colour = NULL, title = "Four sources, two tracers: a wide posterior") + theme_te() + theme(legend.position = "top") ``` A single point estimate would hide all of this. The minimum-norm linear solution, for instance, returns one plausible-looking diet: ```{r} #| label: point-solve M <- rbind(t(Smu), 1) b <- c(colMeans(y4), 1) p_point <- as.numeric(t(M) %*% solve(M %*% t(M), b)) # minimum-norm solution names(p_point) <- rownames(Smu) round(p_point, 3) ``` Those four numbers look like an answer, but they are one arbitrary point inside the posterior cloud. Reporting them alone would imply a precision the data do not support. ## Why the fractions trade off The fractions are not just individually uncertain, they are locked together. Look at the posterior correlations. ```{r} #| label: correlations C <- cor(d); dimnames(C) <- list(rownames(Smu), rownames(Smu)) cor_sc <- round(C["SPOM", "C3"], 2) round(C, 2) ``` Every pair is strongly correlated, and the strongest is SPOM against C3 at `r cor_sc`. The reason is in the mixing space: C3 and SPOM both sit at the depleted end of the carbon axis, so swapping one for the other barely moves the mixture. The data constrain the *sum* of the depleted sources tightly while leaving the *split* between them almost free. That near-one-to-one exchange is the correlation, and it is the fingerprint of two sources the tracers cannot separate. ```{r} #| label: fig-correlation #| echo: false #| fig-cap: "Posterior draws of the SPOM fraction against the C3 fraction. The tight negative ridge means the data fix their sum but not their individual values: whenever one goes up, the other comes down by nearly the same amount." #| fig-alt: "A dense scatter of posterior draws forming a narrow downward-sloping band, showing a strong negative correlation between the SPOM and C3 fractions." cc_df <- data.frame(C3 = d[, 1], SPOM = d[, 4]) ggplot(cc_df, aes(C3, SPOM)) + geom_point(alpha = 0.08, colour = te_green, size = 0.6) + labs(x = "C3 fraction", y = "SPOM fraction", title = "The depleted sources trade off one-for-one") + theme_te() ``` ## What the data can actually resolve If the tracers fix the sum of the two depleted sources, then *that sum* is the quantity the data determine, even though neither source alone is. Add the posterior draws for C3 and SPOM (and, symmetrically, for the two enriched sources C4 and Algae) and the picture sharpens dramatically. ```{r} #| label: aggregate g_depleted <- d[, "C3"] + d[, "SPOM"] # true sum: 0.70 g_enriched <- d[, "C4"] + d[, "Algae"] # true sum: 0.30 summ <- function(z) c(mean(z), sd(z), as.numeric(quantile(z, c(0.025, 0.975)))) agg <- rbind(summ(g_depleted), summ(g_enriched)) dimnames(agg) <- list(c("C3 + SPOM", "C4 + Algae"), c("mean", "sd", "lo", "hi")) sd_ind <- round(mean(apply(d, 2, sd)), 3) sd_agg <- round(mean(c(sd(g_depleted), sd(g_enriched))), 3) gain <- round(sd_ind / sd_agg, 1) round(agg, 3) c(individual_sd = sd_ind, aggregate_sd = sd_agg, narrower_by = gain) ``` The depleted group sums to a posterior mean of `r round(mean(g_depleted), 3)` with a standard deviation of just `r sd_agg`, and its credible interval brackets the true 0.70 tightly. The individual sources had an average posterior standard deviation of `r sd_ind`; the grouped sum is `r gain` times narrower. The strong negative correlation that made the individuals ambiguous is exactly what makes their sum precise: the errors cancel. A mixing model that cannot tell C3 from SPOM can still tell you, with confidence, how much of the diet came from depleted sources as a whole. ```{r} #| label: fig-aggregate #| echo: false #| fig-cap: "The two individual depleted sources (wide, pale) against their sum (narrow, dark). Neither source alone is resolvable, but the quantity that matters ecologically, the total depleted contribution, is pinned down. Aggregating similar sources turns an unanswerable question into an answerable one." #| fig-alt: "Two broad pale density curves for C3 and SPOM overlaid with one narrow dark density curve for their sum, which is much more concentrated." agg_df <- rbind( data.frame(frac = d[, "C3"], grp = "C3 (alone)"), data.frame(frac = d[, "SPOM"], grp = "SPOM (alone)"), data.frame(frac = g_depleted, grp = "C3 + SPOM (sum)")) agg_df$grp <- factor(agg_df$grp, levels = c("C3 (alone)", "SPOM (alone)", "C3 + SPOM (sum)")) ggplot(agg_df, aes(frac, colour = grp, fill = grp)) + geom_density(alpha = 0.25, linewidth = 0.6) + geom_vline(xintercept = 0.70, linetype = "dashed", colour = "grey30") + scale_colour_manual(values = c("#7fa793", "#b3a06a", te_green)) + scale_fill_manual(values = c("#a7c4b5", "#cdbf9a", te_green)) + labs(x = "diet fraction", y = "posterior density", colour = NULL, fill = NULL, title = "The sum is resolvable even when the parts are not") + theme_te() + theme(legend.position = "top") ``` ## Does more data help? If the width came from sampling noise, more consumer measurements would shrink it. Some of it does: raise the consumer sample from 20 to 200 and refit. ```{r} #| label: more-data set.seed(505) y4b <- cbind(rnorm(200, mu4[1], res[1]), rnorm(200, mu4[2], res[2])) fit4b <- mix_mh(y4b, Smu, Ssd, res, iter = 60000, burn = 12000, seed = 606) sd_small <- round(mean(apply(d, 2, sd)), 3) sd_large <- round(mean(apply(fit4b$draws, 2, sd)), 3) ratio <- round(sd_large / sd_small, 2) c(sd_n20 = sd_small, sd_n200 = sd_large, ratio = ratio) ``` Ten times as much data shrinks the individual standard deviations to `r ratio` of their former size, short of the reduction to about a third that pure sampling noise would give (the square root of ten). Part of the uncertainty shrinks with data; part of it is the geometry, and no consumer sample size removes that. The aggregate was already tight at twenty measurements. The lesson is to spend effort where it pays: more tracers to break the geometric ambiguity, or grouping to report what is identifiable, rather than more consumer replicates chasing a limit that data cannot cross. ## Where to go next This posterior is only as trustworthy as its inputs. The source values here were corrected to the consumer's trophic level before any fitting, and that correction, the [trophic discrimination factor](../discrimination-and-concentration-in-mixing/), turns out to move the answer more than the data do. Before reporting a fit, run the [standard checks](../checking-an-isotope-mixing-model/): is the consumer inside the mixing polygon, is the posterior actually informed by the data rather than the prior, and which contributions are resolvable as groups. ## Honest limits The sampler here is deliberately minimal: one consumer, fixed source means and spreads, a flat prior, and a hand-rolled Metropolis step. Production tools add hierarchical structure, covariates, informative priors, and better-mixing samplers, and you should use them for real analyses (SIAR, MixSIAR, and their relatives). What the toy version buys is transparency: the wide posterior and the strong correlations are not artefacts of a black box, they are the honest consequence of asking two tracers to resolve four sources. A fancier sampler would report the same width. The fix is more information, not more computation. ## Related tutorials - [Stable isotope mixing models in R](../isotope-mixing-models-in-r/) - [Discrimination and concentration in mixing models](../discrimination-and-concentration-in-mixing/) - [Checking a stable isotope mixing model](../checking-an-isotope-mixing-model/) - [A Metropolis sampler from scratch](../metropolis-hastings-from-scratch/) ## References Moore JW, Semmens BX 2008 Ecology Letters 11:470-480. Parnell AC, Inger R, Bearhop S, Jackson AL 2010 PLoS ONE 5:e9672 (10.1371/journal.pone.0009672). Parnell AC, Phillips DL, Bearhop S, Semmens BX, Ward EJ, Moore JW, Jackson AL, Grey J, Kelley DJ, Inger R 2013 Environmetrics 24:387-399. Stock BC, Jackson AL, Ward EJ, Parnell AC, Phillips DL, Semmens BX 2018 PeerJ 6:e5096 (10.7717/peerj.5096).