Hierarchical & latent-state Bayesian inference on real data

Is there extra variation in ant richness? An observation-level random effect

The data. ants.csv records ant species richness (the number of ant species found) at 22 sites across New England. At each site ants were sampled twice, once in a bog and once in the forest next to it, so the file has 44 rows. The covariates are forest (1 for the forest sample, 0 for the bog sample), latitude, and elevation. This is the dataset from section A of the lab.

Why this problem. A Poisson regression assumes the variance of the counts equals their mean. Raw richness does not look like that: its mean is 7.0 and its variance is 18.0. Some of that extra spread is expected, though, because habitat, latitude, and elevation shift the mean from row to row. The real question is whether any extra variation is left over once the covariates are in the model. If there is, a plain Poisson regression is overconfident: its standard errors and credible intervals come out too narrow.

One way to find out, and to fix it if needed, is to give every row its own small random nudge on the log mean, $\epsilon_i$. You met this device in Day 3, Problem 5, as an observation-level random effect: "observation-level" because there is one $\epsilon_i$ per row of data, not one per site or group. The nudges are drawn from a Normal distribution with standard deviation $\sigma_\epsilon$, and $\sigma_\epsilon$ is the answer to the question. A $\sigma_\epsilon$ near 0 means the covariates already account for the spread and plain Poisson was fine; a clearly positive $\sigma_\epsilon$ measures how much extra variation there is.

ants <- read.csv(here::here("Day5_BayesII", "data", "ants.csv"))
ants$forest    <- as.numeric(scale(ants$forest))
ants$latitude  <- as.numeric(scale(ants$latitude))
ants$elevation <- as.numeric(scale(ants$elevation))

1. Fit the model in NIMBLE with weakly informative priors ($\beta \sim \mathcal{N}(0, \sigma^2{=}4)$, $\sigma_\epsilon \sim \text{Uniform}(0,5)$). Report posterior medians and 95% CIs for $\beta_1,\beta_2,\beta_3$ and for $\sigma_\epsilon$.
2. Look at the posterior of $\sigma_\epsilon$. Does it point to meaningful leftover variation, or is it consistent with values close to 0? What does that tell you about overdispersion once the covariates are in the model?
3. Now fit the plain Poisson regression without the random effect, glm(richness ~ forest + latitude + elevation, family = poisson), using the scaled covariates.
   1. Compare your $\beta$ posterior medians and intervals to the glm() estimates. Where do they agree or differ, and why?
   2. Is there evidence of overdispersion from the glm() fit? Divide its residual deviance by its residual degrees of freedom; a ratio near 1 means no overdispersion, and a ratio well above 1 means the counts vary more than Poisson allows (Day 3, Problem 5).
   3. Does that agree with what you inferred from $\sigma_\epsilon$ in the Bayesian model in (b)? Was the random effect needed?

What does a "vague" log-scale prior actually assert?

Practice Exercise M1 had you draw from a vague log-scale prior and look at the multiplicative effect it implies across the latitude range; Practice Exercise M2 showed what that same prior does to the chains when the predictors are unscaled. This problem picks the M1 draw back up and takes it two steps further, onto a per-unit effect you can state in a sentence, and into an argument you would actually have with a reviewer. On a log link a prior on $\beta$ is a prior on a multiplicative effect $e^\beta$, and that translation is where "uninformative" priors stop being uninformative.

1. Draw $10^5$ values from $\mathcal{N}(0,\ \sigma^2 = 10000)$, the classic "vague" choice. Report the fraction for which $e^{\beta}$ exceeds a million-fold, and the fraction for which the effect exceeds a million-fold in either direction. State the second number in a sentence beginning "before seeing any data, this prior says…".
2. Repeat with $\mathcal{N}(0,\ \sigma^2 = 4)$. Report the same fraction, and also the central 95% range of multiplicative effects this prior allows (i.e. $e^{\pm 1.96 \times 2}$). Is that range wide enough for any ecological effect you have ever measured?
3. A reviewer objects that $\mathcal{N}(0, 4)$ is "informative" and demands $\mathcal{N}(0, 10000)$ for objectivity. Using your two numbers, write the two-sentence reply you would actually send.

The pooling spectrum for soil N₂O emissions

The data. Nitrous oxide (N₂O) is a potent greenhouse gas, and fertilized farm soils are a major source. N2OEmission.csv compiles 563 emission measurements from 107 sites. Each row is one measurement: a fertilizer applied at some rate at one site, and the N₂O emission that followed. The question running through Problems 3 and 4 is how emissions respond to the amount of nitrogen applied, and how much that differs among sites and fertilizers.

• emission: the N₂O emission measured. The response, used as $y_i = \log(\text{emission})$.
• n.input: the amount of nitrogen fertilizer applied. The predictor, used as $x_i = \log(\text{n.input})$ minus its mean. Because both are logged, the slope $\beta$ is unit-free: $\beta = 1$ means doubling the fertilizer doubles the emission.
• group.index: which of the 107 sites the row comes from, numbered 1–107 (group holds the original site IDs). This is the site index $j$, so row $i$ belongs to site $j[i]$. Sites have between 1 and 25 measurements.
• fert.index: which of 10 fertilizer types was applied, numbered 1–10 (fertilizer holds the original letter codes). This is the fertilizer index $k$. It is not used until Problem 4(b).
• carbon: soil organic carbon at the site, as a percentage. It is not used until Problem 4.
• reps: a replicate count carried over from the source compilation. Not used.

Model $\log(\text{emission})$ as a function of centred $\log(\text{n.input})$, with sites indexed by group.index. Inside the model you will write the site intercept as alpha[group[i]]: group[i] is row $i$'s site number, so this picks out that row's own intercept from the vector of 107. Section C of the practice lab explains that notation and has a fuller description of the data.

N2O   <- read.csv(here::here("Day5_BayesII", "data", "N2OEmission.csv"))
y     <- log(N2O$emission)
x     <- log(N2O$n.input) - mean(log(N2O$n.input))  # centred predictor
group <- N2O$group.index   # site number of each row; use it as alpha[group[i]]
nsites <- length(unique(group))

1. Fit the complete-pooling model $\mu_i = \alpha + \beta x_i$ and the partial-pooling random-intercept model $\alpha_j \sim \mathcal{N}(\mu_\alpha, \sigma_\alpha^2)$, $\mu_i = \alpha_{j[i]} + \beta x_i$. Report the posterior of $\beta$ from each.
2. Also fit the no-pooling model (independent $\alpha_j$). Plot the per-site $\alpha_j$ posterior medians from no-pooling (x-axis) vs. partial-pooling (y-axis) with the 1:1 line. Which sites shrink the most toward $\mu_\alpha$, and why?
3. Report the posterior of $\sigma_\alpha$. In plain language, what does its magnitude say about how different the sites are once N input is accounted for?

A site-level covariate on the random intercept, and a new random slope by fertilizer type

Every Bayesian hierarchical model so far today, the lab's N₂O models and Problem 3, has had random intercepts only. In Problem 3 each site got its own intercept $\alpha_j$, but all 563 measurements shared one N-input slope $\beta$: every site and every fertilizer was assumed to respond to added nitrogen at the same rate. Problem 3's model also says that sites differ, through $\sigma_\alpha$, but not why. This problem changes both, in one model that you fit once. It adds two pieces of structure at the same time, and they use two different columns, so keep them separate:

• A site-level covariate on the random intercept, to explain site differences with soil carbon. The site-level covariate is $w_j$, one value per site, built from the carbon column. Average carbon within each site to get $c_j$, the site's percent soil organic carbon (36 sites report more than one value), then take its logit, $w_j = \log\big(c_j/(100-c_j)\big)$. State how you scaled $w_j$ (for example, standardized to mean 0 and SD 1).
• A random slope, so the N-input slope differs by fertilizer type. This is new for today. Until now only intercepts varied by group; here the slope becomes a random effect too. You fit random slopes with lmer() on Day 3; this is the same idea written out in NIMBLE. The single $\beta$ is replaced by ten fertilizer-specific slopes $\beta_k$, each drawn from a shared distribution $\mathcal{N}(\mu_\beta, \sigma_\beta^2)$, so they are partially pooled toward the average slope $\mu_\beta$ exactly the way site intercepts are pooled toward their mean. The grouping for the slope is fert.index, not the site and not $w_j$. Fertilizer type is recorded per measurement, not per site (47 sites tested more than one type), so each row $i$ gets the slope of its own type, $\beta_{k[i]}$ with $k[i]$ = fert.index.

The full model, with both additions, written out with $j[i]$ = group.index and $k[i]$ = fert.index for row $i$:

• $\kappa$ (kappa) and $\eta$ (eta) are the intercept and slope of a regression of the site intercepts on soil carbon. $\kappa$ is the expected $\alpha_j$ for a site with $w_j = 0$ (average carbon, once $w$ is standardized), and it replaces Problem 3's $\mu_\alpha$. $\eta$ is how much a site's typical log-emission changes per unit of $w_j$.
• $\sigma_\alpha$ is the scatter of site intercepts around that line, the site variation carbon does not explain.
• $\mu_\beta$ is the average N-input slope across fertilizer types, and $\sigma_\beta$ is how much the ten type-specific slopes $\beta_k$ differ from it.

Fit this model once. Part (a) interprets the site-level effect of soil carbon from that fit, and part (b) interprets the random slope from the same fit.

1. Fit the full model above in NIMBLE, with both the soil-carbon covariate on the site intercepts and the random slope by fertilizer type. Site-level effect of soil carbon: report the posterior median, 95% credible interval, and $\Pr(\eta > 0)$ for $\eta$. Is there evidence that soil carbon affects baseline emissions, and in which direction? Say how big the effect is in plain terms.
2. Random slope by fertilizer type, from the same fit: report the posterior of $\sigma_\beta$ and rank the ten types by posterior median slope. Then test the ranking before you believe it: across your draws, count how often the top-ranked type actually holds the largest $\beta_k$ of the ten, and report that probability together with the type's 95% CI and the number of measurements it rests on. On this evidence, would you write “fertilizer type X produces the steepest emissions response” in a paper? One or two sentences.

Island lizard occupancy vs. perimeter:area ratio

Load the island-lizard data and fit logistic regression of presence on the scaled perimeter:area ratio:

liz      <- read.csv(here::here("Day5_BayesII", "data", "IslandsLizards.csv"))
x.scaled <- as.numeric(scale(liz$perimeterAreaRatio))
y        <- as.integer(liz$presence)

1. Fit in NIMBLE with the weakly informative prior $a, b \sim \mathcal{N}(0, \sigma^2{=}2.71)$ (i.e. dnorm(0, 1/2.71), chosen because it implies a near-uniform prior on $p$). Report R-hat for $a$ and $b$ (below 1.05 is the target; up to about 1.1 is acceptable). Use starting values near 0 (why is that a good idea here?).
2. What is the posterior probability that the slope $b$ is negative (i.e., more edge per unit area $\Rightarrow$ lower occupancy)?
3. A management-relevant difference. Add derived quantities inside the model, so MCMC delivers their full posterior:
   $p_{10} = \text{ilogit}(a + b\, x_{10}), \quad p_{20} = \text{ilogit}(a + b\, x_{20}), \quad \text{diff} = p_{20} - p_{10}$
   where $x_{10}$ and $x_{20}$ are the scaled versions of raw perimeter:area ratios 10 and 20 (scale them with the same mean and SD you used for the data). $p_{10}$ is the probability that an island with a ratio of 10 is occupied, so its posterior is the posterior distribution of that occupancy probability. It is not a posterior predictive distribution for $y$: that would be a new 0/1 outcome, drawn as y.new ~ dbern(p10). Report the posterior medians and 95% CIs of $p_{10}$, $p_{20}$ and $\text{diff}$. Does the interval for $\text{diff}$ exclude 0?
4. State in one sentence why computing $\text{diff}$ as a node inside nimbleCode gives correct uncertainty, whereas plugging point estimates of $a$ and $b$ into the formula would not.

Beverton-Holt fit with management quantities

Simulate a spawner-recruit time series from a known truth, then fit the Beverton-Holt model in NIMBLE and define the management quantities as deterministic nodes so their posteriors come for free.

set.seed(7)
Nyears <- 25; a.t <- 5; b.t <- 5000; sig.t <- 0.3
S <- rlnorm(Nyears, 7.5, 0.8)
R <- numeric(Nyears)
for (t in 1:Nyears)
  R[t] <- (a.t*S[t])/(1 + a.t/b.t*S[t]) * exp(rnorm(1, 0, sig.t))

1. Fit, with the reference points inside the model. Monitor $a,b,\sigma,N_\text{eq},S_\text{msy},H_\text{msy}$. Before you look at the output, work out the three true reference points by hand, substituting $a=5$ and $b=5000$ into the formulas above. Doing that arithmetic yourself is how you sanity-check any assessment. Report posterior medians and 95% CIs for all six quantities next to your hand-computed truths. Does every interval cover?
2. A harvest rate with a risk statement attached. The posterior median of $H_\text{msy}$ is, by construction, the harvest rate with a 50% chance of exceeding the true one. Report the 5th, 10th, 25th and 50th percentiles of the posterior. Which would you advise, and write the one-sentence risk statement that has to travel with it.
3. An escapement target. Report $S_\text{msy}$'s posterior median and 95% CI, and compute $\Pr(S_\text{msy} > 1000)$, $\Pr(S_\text{msy} > 1500)$ and $\Pr(S_\text{msy} > 2000)$. A manager asks you for “the” escapement target. Answer them in two sentences that use these numbers rather than hiding them.

Why you cannot build a derived quantity from summary numbers

You have posterior draws of $a$ and $b$ from Problem 6. Management quantities are functions of those parameters, and there are two ways to produce an interval for them, one right and one that looks reasonable and is not.

1. Compute $S_\text{msy} = b\sqrt{1/a} - b/a$ once per draw, then take the 5th, 50th and 95th percentiles of the result. Now do it the other way. Take the 5th, 50th and 95th percentiles of $a$ and of $b$ separately, and substitute those numbers into the formula. Report both sets of three numbers and the percentage difference at the 5th percentile.
2. Explain why the two disagree. What does cor(a, b) across your draws have to do with it, and why does summarising each parameter first throw away exactly the information the derived quantity needs?
3. Now repeat the comparison for $H_\text{msy} = 1 - \sqrt{1/a}$. The two approaches agree here, to the digit. What is different about this quantity, and what does that tell you about when the shortcut is safe? State the general rule in one sentence.

Repeat-count abundance and a goodness-of-fit check

An N-mixture model estimates the abundance of unmarked animals when you cannot count every individual. Repeated visits to the same closed site let it separate true abundance ($\lambda$) from detection probability ($p$), which a raw count confounds. (It is the abundance analog of the occupancy model; see the practice lab for a fuller explanation.) Here you will simulate repeated counts, fit the Royle (2004) N-mixture model in NIMBLE, and then check, step by step, whether data simulated from the fitted model look like the real counts.

library(MASS)
set.seed(9)
sites <- 150; reps <- 3; lambda <- 6; p.true <- 0.5
N <- rnegbin(sites, mu = lambda, theta = 2)   # NB, not Poisson
y <- matrix(NA, sites, reps)
for (i in 1:sites) y[i, ] <- rbinom(reps, N[i], p.true)

1. Fit the basic (constant $\lambda$, constant $p$) N-mixture. Remember that NIMBLE's binomial is dbinom(prob, size), the reverse of R's argument order. Give the latent abundances starting values with N = apply(y, 1, max). Here y is the sites × visits matrix of counts, the 1 means "work across each row" (each site), and max is R's ordinary max() function, so this returns the largest count seen at each site on any visit. The purpose: a binomial count can never exceed its size, so each site's true abundance $N_i$ must start at least as large as every count recorded there. Start it lower and that count is impossible, its log-probability is $-\infty$, and the sampler cannot begin from that point. (See the apply family tutorial for apply().) Report the posteriors of $\lambda$ and $p$.

2. Does the model reproduce the data? A posterior predictive check, step by step. The idea is simple. If the fitted model is a good description of how these counts arose, then fake datasets simulated from it should look like the real one. If the real data have some feature the fake datasets almost never show, the model is missing something.

   1. Simulate fake datasets inside the model. Add a second block to your nimbleCode that invents a brand-new set of sites, the same size as the real data, from the same model:

      for (i in 1:Nsites) {
        N.new[i] ~ dpois(lambda)
        for (j in 1:Nreps) {
          y.new[i, j] ~ dbinom(prob = p, size = N.new[i])
        }
      }

      N.new and y.new have no data attached, so at every MCMC iteration NIMBLE simply simulates them from the current draw of $\lambda$ and $p$. Each saved iteration therefore gives you one complete fake dataset of 150 sites × 3 visits, and the spread across iterations includes your uncertainty about $\lambda$ and $p$. This collection is the posterior predictive distribution of the data. Monitor y.new (thinning to every 10th iteration keeps the output manageable) and give it starting values, for example y.new = y with N.new started at N's starting values.
   2. Choose features you can read. Summarize any dataset of counts with three familiar numbers: the variance-to-mean ratio of all the counts (about 1 for Poisson-like counts, well above 1 when counts are more spread out than Poisson allows), the proportion of counts that are zero, and the largest count. Compute all three for the real data. Then compute them for every fake dataset: pull the y.new columns out of as.matrix(samples), so that each row is one fake dataset, and use apply(..., 1, ...) to summarize each row.
   3. Compare. For each feature, draw a histogram of its values across the fake datasets and mark the real data's value with a vertical line. Report the real value, the middle 95% of the fake values, and the posterior predictive p-value: the proportion of fake datasets whose value is at least as large as the real one. A p-value near 0.5 means the real data look typical; a value near 0 (or near 1) means the real data are more extreme than almost anything the model produces.
3. Which features does the model fail to reproduce, and in which direction? What does that tell you about the assumption that abundance is Poisson? Propose a fix, in words or as a modified nimbleCode block (for example, a random effect on $\log\lambda_i$ for each site, a Poisson-lognormal model), and say how you would expect it to change the three checks.

4. A standardized check: the chi-square-type sum. The three features in (b) were chosen with this problem in mind. A more general discrepancy, one that works the same way for any count model, is the chi-square-type sum over every count:

   $T(y, \theta) = \sum_{i}\sum_{j} \dfrac{\big(y_{ij} - \mathbb{E}[y_{ij} \mid \theta]\big)^2}{\mathbb{E}[y_{ij} \mid \theta]}$

   Each term is the squared distance between a count and the count the model expects, divided by that expected count. The division is what standardizes it: for counts the variance is roughly the mean, so a miss of 3 when 2 were expected counts for much more than a miss of 3 when 30 were expected. Because the sum needs nothing but each observation's expected value, the same recipe applies to any model, only $\mathbb{E}[y_{ij} \mid \theta]$ changes. That makes it a standard way to compute a Bayesian p-value and to compare how well different models fit the same data, such as this Poisson N-mixture and the fix you proposed in (c).

   1. The expected count. In this model every count has the same expected value, $\mathbb{E}[y_{ij} \mid \theta] = p\,\lambda$: on average a site holds $\lambda$ animals and each is detected with probability $p$.
   2. Compute it at every draw, twice. Unlike the features in (b), $T$ depends on the parameters, so compute it separately at each posterior draw $m$. Use that draw's $p^{(m)}\lambda^{(m)}$ to compute $T$ for the real counts, $T(y, \theta^{(m)})$, and for that draw's fake dataset, $T(y^\text{rep}_m, \theta^{(m)})$. Both use the same expected value, so the only difference between them is real data versus model-generated data.
   3. Compare. The Bayesian p-value is the proportion of draws at which the fake dataset's $T$ is at least as large as the real data's. Report it and the median of each set of $T$ values, and plot $T(y^\text{rep})$ against $T(y)$ with a 1:1 line: a model that fits scatters its points evenly on both sides of the line. Does this single number agree with what the three features in (b) told you?

Multi-season occupancy: recover persistence and colonization

Simulate a dynamic (multi-season) occupancy dataset and recover the year-specific persistence $\phi_t$ and colonization $\gamma_t$ from replicated detection/non-detection data with imperfect detection.

set.seed(10)
nS <- 50; nY <- 6; nR <- 3; p.t <- 0.5
phi.t <- c(0.70, 0.70, 0.80, 0.85, 0.85)  # TRUE persistence, one per transition
gam.t <- c(0.05, 0.15, 0.25, 0.30, 0.30)  # TRUE colonization, rising over time
z <- matrix(NA, nS, nY); z[, 1] <- rbinom(nS, 1, 0.2)
for (i in 1:nS) for (t in 2:nY)
  z[i,t] <- rbinom(1, 1, z[i,t-1]*phi.t[t-1] + (1-z[i,t-1])*gam.t[t-1])
y <- array(NA, c(nS, nR, nY))
for (t in 1:nY) for (i in 1:nS) for (j in 1:nR)
  y[i,j,t] <- rbinom(1, 1, p.t) * z[i,t]

1. Fit the model in NIMBLE, using the z-init trick (zst <- matrix(1L, nS, nY)) that the walkthrough works through. Report $\hat R$ for $p$ and for the $\gamma_t$.
2. Plot the posterior medians of $\gamma_t$ against year with their 95% intervals. Does the rising colonization signal come through? Report the mean width of those five intervals, and say what that width is telling you. Think about how many sites actually contribute information to a single year's transition.
3. Two ways to tighten them. Refit twice more, changing only how $\gamma_t$ is built:
   1. a hierarchical prior, $\text{logit}(\gamma_t) \sim \mathcal{N}(\mu_\gamma,\ \sigma_\gamma^2)$, so the years borrow strength from each other;
   2. a trend, $\text{logit}(\gamma_t) = g_0 + g_1 \cdot \text{year}_t$, which is a reasonable thing to try given how you simulated the data.
   For all three fits report two numbers for the five $\gamma_t$: the mean width of their 95% intervals (a measure of precision), and the RMSE of their posterior medians (a measure of accuracy).

   What RMSE is. Because you simulated the data, you know the true colonization probabilities: they are the vector gam.t in the simulation code above, one value per transition. The root mean squared error measures how far the estimates sit from that truth, on average:

   $\text{RMSE} = \sqrt{\dfrac{1}{5}\sum_{t=1}^{5}\big(\hat\gamma_t - \gamma_t^{\text{true}}\big)^2}$

   Take each year's error (posterior median minus truth), square it so misses in either direction count and big misses count extra, average the five squares, and take the square root to return to the original units. An RMSE of 0.05 means the estimates are typically about 0.05 away from the true probabilities. Smaller is better, and 0 would be perfect. In R, with gam.med holding the five posterior medians:

   sqrt(mean((gam.med - gam.t)^2))

   You can only compute RMSE here because the truth is known; with real data it is not. Which structure narrowed the intervals most? Which estimated $\gamma_t$ most accurately? Are they the same model?