Why a Beta distribution?

In a stage-structured matrix with \(k\) possible fates (including death), survival probabilities form a multinomial, and the natural conjugate prior is a Dirichlet distribution (the approach used by raretrans). In this particular age-structured model, however, an individual of a given age can only survive to the next age class or die—exactly two possible outcomes. The Binomial likelihood paired with a \(\text{Beta}(\alpha, \beta)\) prior is therefore sufficient.

The Beta distribution is bounded on \([0, 1]\), making it appropriate for probabilities. Its two parameters \(\alpha\) and \(\beta\) can be interpreted as pseudo-counts: \(\alpha\) “prior survivals” and \(\beta\) “prior deaths”. Their sum \(\alpha + \beta\) acts as an effective sample size that controls how concentrated the distribution is around its mean.

Moment matching: from mean and SD to \(\alpha\) and \(\beta\)

The mean and variance of a \(\text{Beta}(\alpha, \beta)\) distribution are:

\[ \mu = \frac{\alpha}{\alpha + \beta}, \qquad \sigma^2 = \frac{\alpha \beta}{\left(\alpha + \beta\right)^2\left(\alpha + \beta + 1\right)} \]

Beissinger (1995) reports \(\mu\) and \(\sigma\) (Table 2). We can recover approximate \(\alpha\) and \(\beta\) by moment matching—equating the distribution’s theoretical moments to the observed summary statistics. Defining the auxiliary quantity \(\nu\):

\[\begin{align*} \nu &= \frac{\mu\left(1-\mu\right)}{\sigma^2} \\ \alpha &= \nu\,\mu \\ \beta &= \left(1-\mu\right)\,\nu \end{align*}\]

This conversion is valid as long as \(\sigma^2 < \mu(1 - \mu)\), which is satisfied for all entries in Beissinger’s Table 2. The code below performs this conversion and displays the resulting \(\alpha\) and \(\beta\) values.

kite_parms <- b1995_table2 %>% 
  slice(-(1:2)) %>% # remove top two rows
  separate(col = range_ns, into = c("min_ns","max_ns"), sep = "-") %>% 
  separate(col = range_s, into = c("min_s","max_s"), sep = "-") %>% 
  mutate(across(.cols = 3:12, .fns = as.numeric)) %>% # convert to numeric
  mutate(nu_ns = mean_ns * (1 - mean_ns) / sd_ns^2,
         alpha_ns = mean_ns * nu_ns,
         beta_ns = (1-mean_ns)*nu_ns,
         nu_s = mean_s * (1 - mean_s) / sd_s^2,
         alpha_s = mean_s * nu_s,
         beta_s = (1-mean_s)*nu_s) %>% 
  # set infinite vales to missing
  mutate(across(where(is.numeric), ~ifelse(is.nan(.x), NA_real_, .x))) 
kite_parms %>% 
  select(env_state, age_class, alpha_ns, beta_ns, alpha_s, beta_s) %>% 
  mutate(across(where(is.numeric), ~format(.x, digits = 2))) %>% 
  bind_rows(tibble(env_state = c("","Environmental State"),
                   age_class = c("", "Age Class"),
                   alpha_ns = c("Nest Success", "$\\alpha$"),
                   beta_ns = c("", "$\\beta$"),
                   alpha_s = c("Survival", "$\\alpha$"),
                   beta_s = c("", "$\\beta$")), .) %>% 
  mutate(across(.cols = 3:6, ~case_when(grepl("NA",.x) ~ " ", 
                                        TRUE ~ .x))) %>% 
#  filter(age_class != "Fledgling") %>% 
  hux() %>% 
  set_escape_contents(row = 2, col = 3:6, value = FALSE) %>% 
  theme_article(header_col = FALSE) %>% 
  set_na_string(value = " ") %>% 
  set_position("left") %>% 
  set_caption("Values of $\\alpha$ and $\\beta$ for each age class and 
              environmental state. Fledglings have zero probability of 
              nest success in all environmental states.")

Interpreting the effective sample size

One of the key benefits of re-expressing beliefs as Beta parameters is transparency: the sum \(\alpha + \beta\) is the effective sample size, so a reader can immediately see how much confidence the analyst has placed in each estimate.

Looking at the table above, nest-success estimates carry less confidence in drought and lag years than in high water years, consistent with what Beissinger (1995) states qualitatively. For survival, however, two patterns emerge that are inconsistent with Beissinger’s intended beliefs:

1. Within each environmental state, fledgling survival has a larger effective sample size than subadult or adult survival. This happens because fledgling survival is lower (mean \(\approx 0.5\)–\(0.9\)), so \(\mu(1-\mu)\) is larger relative to a fixed \(\sigma^2\).
2. High water year estimates appear less confident (smaller \(n\)) than lag year estimates, despite being based on the best data.

Both artifacts arise because Beissinger fixed the standard deviation at 0.03 or 0.10 regardless of the mean. For a normally distributed parameter, the mean and standard deviation are independent, so a fixed SD implies a fixed level of precision. For a Beta distribution, however, the variance depends on both the mean and the concentration. Fixing \(\sigma\) while changing \(\mu\) implicitly changes the effective sample size in counter-intuitive ways.

This observation has a practical lesson: when converting from a normal parameterisation to a Beta, it is worth checking whether the implied effective sample sizes are scientifically reasonable. If they are not, the analyst may want to specify \(\alpha\) and \(\beta\) directly rather than relying on moment matching.

pdf_seq <- seq(-0.2,1.2,0.001)
plot_colors <- RColorBrewer::brewer.pal(3, "RdYlBu")
kite_parms2 <- kite_parms %>% 
  mutate(env_state = rep(env_state[c(1,4,7)], each = 3),
         n_ns = alpha_ns + beta_ns,
         n_s = alpha_s + beta_s,
         lwr_ns = pbeta(min_ns, alpha_ns, beta_ns),
         upr_ns = pbeta(max_ns, alpha_ns, beta_ns),
         lwr_s = pbeta(min_s, alpha_s, beta_s),
         upr_s = pbeta(max_s, alpha_s, beta_s),
         plot_ns = map2(alpha_ns, beta_ns, ~tibble(x_ns = pdf_seq,
                                                   d_ns = dbeta(x_ns, .x, .y))),
         plot_ns_norm = map2(mean_ns, sd_ns, ~tibble(x_nsn = pdf_seq,
                                                     d_ns_norm = dnorm(x_nsn, .x, .y),
                                                     d_ns_tnorm = dnorm(x_nsn, .x, .y)/
                                                       pnorm(0, .x, .y, lower.tail = FALSE))),
         plot_s = map2(alpha_s, beta_s, ~tibble(x_s = pdf_seq,
                                                d_s = dbeta(x_s, .x, .y))),
         plot_s_norm = map2(mean_s, sd_s, ~tibble(x_sn = pdf_seq,
                                                  d_s_norm = dnorm(x_sn, .x, .y),
                                                  d_s_tnorm = dnorm(x_sn, .x, .y)/
                                                    pnorm(1, .x, .y))),
  ) %>% 
  unnest(cols = c("plot_ns", "plot_s", "plot_ns_norm","plot_s_norm"), 
         names_repair = "unique") %>% 
  mutate(env_state = factor(env_state, levels = c("Drought year","Lag year", "High year")),
         age_class = factor(age_class, levels = c("Fledgling", "Subadult", "Adult")))

ggplot(data = kite_parms2) + 
  geom_area(mapping = aes(x = x_s, y = d_s), alpha = 0.75, fill = plot_colors[1], color=plot_colors[1], size = 1) + 
  geom_area(mapping = aes(x = x_sn, y = d_s_tnorm), alpha = 0.75, fill = plot_colors[3], color = plot_colors[3], size = 1) +
  facet_grid(env_state~age_class, scales = "free_x") +
  scale_x_continuous(limits = c(0,1)) + 
  labs(x = "Annual survival", y = "Density")+
  theme_classic()

Probability distributions of annual survival by stage and environmental state. Beta distributions are red, truncated normal distributions are blue. Areas of overlap appear gray.

Comparing Beta and normal priors

The figure above overlays the Beta distribution (red) with the truncated normal (blue) used in the original PVA. Several points are worth noting:

• For most stage–environment combinations the two distributions are nearly indistinguishable, confirming that the normal approximation works well when the mean is far from 0 or 1 and the standard deviation is small.
• The approximation breaks down for subadult and adult survival in high water years, where mean survival is 0.95. Here the normal distribution assigns roughly 5% probability to values above 1.0—an impossible range for a probability. Beissinger (personal communication) dealt with this by rejecting inadmissible draws and resampling, effectively producing a truncated normal.
• The truncation has little practical effect in this case, but in general it shifts the effective mean upward and can subtly bias simulation results.

The Beta distribution avoids these boundary issues entirely because it is defined on \([0, 1]\) by construction. This is one of the simplest arguments for preferring Beta (or Dirichlet) priors over normal distributions when modelling probabilities.

Components of reproduction

In the snail kite model, the fertility entry in column \(j\) of the matrix is the product of four components:

\[ F_j = (\text{proportion nesting})_j \;\times\; (\text{nest attempts per year})_j \;\times\; (\text{nest success})_j \;\times\; (\text{young per successful nest}) \]

Beissinger (1995) treated all components except nest success as fixed values. For nest success—the probability that a nest produces at least one fledgling—means, standard deviations, and ranges were provided and sampled from normal distributions during the simulation.

Nest success as a Beta variable

Nest success is inherently a probability (bounded between 0 and 1), so a Beta distribution is the natural choice. Snyder et al. (1989) provide raw counts of successful and failed nests (their Table 3) that can be converted directly into Beta parameters: \(\alpha\) = number of successful nests and \(\beta\) = number of failed nests. This is more informative than the moment- matching approach because it uses the actual sample sizes rather than summary statistics.

Young per successful nest

Snyder et al. (1989) report an average of approximately 2 young per successful nest. Their Table 5 provides the actual distribution: 55, 110, and 47 nests produced 1, 2, or 3 young, respectively. This count data could be represented as a multinomial variable with a Dirichlet prior—exactly the approach that raretrans::fill_transitions() uses for transition probabilities. Although the mean number of young does not vary across environmental states, the sample sizes do, and this should be reflected in the confidence placed in each distribution. For simplicity, we use a fixed value of 2 young per successful nest in what follows.

Comparing prior sources

The code below builds Beta priors for nest success from the Snyder et al. (1989) counts and compares them with the moment-matched Beta and truncated normal distributions implied by Beissinger (1995).

repro_parms <- tribble(
  # prior_young is vector of nests with 1, 2, or 3 young. 
  # alpha_ns and beta_ns are the number of successful and unsuccessful nests
  ~env_state, ~age_class, ~prior_young, ~alpha_ns, ~beta_ns,
  "Drought year", "Subadult", c(1, 3, 1), 3, 91, 
  "Drought year", "Adult", c(1, 3, 1), 3, 91,
  "Lag year", "Subadult", c(6, 25, 8), 11, 58,
  "Lag year", "Adult", c(6, 25, 8), 11, 58,
  "High year", "Subadult", c(48, 82, 38), 100, 236,
  "High year", "Adult", c(48, 82, 38), 100, 236
)

repro_parms2 <- repro_parms %>% 
  mutate(n_ns = alpha_ns + beta_ns,
         plot_ns = map2(alpha_ns, beta_ns, ~tibble(x = pdf_seq,
                                                   d_ns = dbeta(x, .x, .y)))) %>% 
  unnest(plot_ns) %>% 
  mutate(env_state = factor(env_state, levels = c("Drought year","Lag year", "High year")),
         age_class = factor(age_class, levels = c("Fledgling", "Subadult", "Adult")))

ggplot(data = filter(repro_parms2, age_class != "Fledgling"),
       mapping = aes(x = x, y = d_ns)) + 
  geom_area(fill = plot_colors[2], alpha = 0.75, color=plot_colors[2], size = 1) + #Synder etal 1989
  geom_area(data = filter(kite_parms2, age_class != "Fledgling"), #Beissinger 1995
       mapping = aes(x = x_ns, y = d_ns), fill = plot_colors[1], alpha = 0.75, color=plot_colors[1], size = 1) + 
  geom_area(data = filter(kite_parms2, age_class != "Fledgling"), #Beissinger 1995
       mapping = aes(x = x_ns, y = d_ns_tnorm), fill = plot_colors[3], alpha = 0.75, color=plot_colors[3], size = 1) +   
  facet_grid(env_state~age_class, scales = "free_y") + 
  scale_x_continuous(limits = c(0,0.6)) +
#  scale_y_continuous(limits = c(0,30)) +
  labs(x = "Nest success", y = "Density") +
  rlt_theme

Prior beliefs about nest success using Beta distributions parameterized directly from Synder et al. (yellow, 1989), Beta distributions parameterized from Beissinger (red, 1995), and normal distributions from Beissinger (blue, 1995).

Interpretation

The figure reveals an important discrepancy. The yellow densities (Snyder et al. counts) are much narrower than the red densities (Beissinger moment-matched Beta), meaning the raw count data imply substantially more confidence in nest success than the summary statistics suggest. This may reflect a deliberate choice: because the other components of reproduction (proportion nesting, nest attempts, young per nest) were treated as known constants, Beissinger may have inflated the variance on nest success to partially compensate for the missing uncertainty in those components.

The blue densities (truncated normal from Beissinger) highlight a second issue. For drought years, where mean nest success is only 0.03, the normal distribution with \(\sigma = 0.10\) assigns substantial probability to negative values. After truncation to \([0, 1]\), the effective mean shifts upward to 0.092—roughly three times the intended value of 0.03. This systematic bias means the original PVA likely overestimated reproduction in drought years.

The Beta distribution avoids this problem entirely: it is defined on \([0, 1]\), so no truncation or rejection sampling is needed, and the mean remains exactly where the analyst intends. This is a compelling practical reason to prefer Beta priors for probabilities, especially when the mean is close to the boundaries.