---
title: "Proportions and success rates"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Proportions and success rates}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(drmTMB)
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

Proportion data are common in ecology and evolution, but the measurement
process matters. Some proportions are counts of successes out of known trials:
germinated seeds out of seeds planted, infected hosts out of hosts checked, or
occupied quadrats out of quadrats surveyed. Other proportions are continuous
rates already measured on the unit interval, such as leaf damage, vegetation
cover, or the fraction of time an animal spends foraging. This tutorial
assumes the location-scale reading pattern from [When variance carries
signal](location-scale.html) and extends it to a bounded response.

This article stays inside the currently implemented `drmTMB` bounded-response
surface:

- use `stats::binomial(link = "logit")` for ordinary 0/1 events or counted
  successes out of known trials with binomial sampling variation only;
- use `beta_binomial()` for counted successes out of known trials when
  extra-binomial variation is part of the scientific question;
- use `beta()` for continuous proportions strictly inside `(0, 1)`;
- use `zero_one_beta()` for continuous proportions on `[0, 1]` when exact 0 or
  1 values are structural outcomes rather than trial-count outcomes.

The source motivation comes from Nakagawa et al. (2026), who treat
overdispersed proportions as part of the same location-scale framework as
continuous and count responses. They distinguish discrete proportions from
continuous proportions, recommend beta-binomial models for successes out of
trials when extra-binomial variation is present, and warn that exact 0 or 1
values in continuous beta responses need explicit zero-one processes rather
than being hidden by the model.

## Choose The Proportion Route

Start from how the response was measured:

| Response | Example | First `drmTMB` family | Boundary status |
|---|---|---|---|
| 0/1 event indicator | early NDVI onset, host infected or not, quadrat occupied or not | `stats::binomial(link = "logit")` | rows are ordinary Bernoulli trials; no modelled `sigma` |
| Counted successes out of trials with ordinary binomial variation | germinated seeds out of planted seeds | `stats::binomial(link = "logit")` | successes can be 0 or all trials through the binomial count process |
| Counted successes out of trials with extra-binomial variation | germinated seeds out of planted seeds, with tray-to-tray probability scatter | `beta_binomial()` | adds modelled `sigma` for variation beyond binomial sampling |
| Continuous proportion strictly inside `(0, 1)` | leaf-area damage, vegetation cover, foraging time fraction | `beta()` | exact 0 or 1 is not supported by the strict beta likelihood |
| Continuous proportion with structural 0 or 1 values | 0% disease, 100% cover, complete absence or saturation | `zero_one_beta()` | `zoi` models exact boundary mass; `coi` models exact ones among boundary observations |

The distinction is not cosmetic. A denominator-aware model knows that 3
successes out of 6 trials and 30 successes out of 60 trials carry different
sampling information. A strict beta model assumes the observation is already a
continuous rate and that the boundaries are not part of the observed support.
For denominator-aware responses, write `cbind(successes, failures)`. Do not
write `weights = trials`, `successes / trials`, or
`cbind(successes, trials)`: top-level weights are likelihood weights, and the
second column of a two-column binomial response is failures, not trial totals.

## Plain Binomial Equation And Syntax

Use `stats::binomial(link = "logit")` when the response is an event
probability and ordinary binomial sampling variation is enough:

\[
\begin{aligned}
Y_i \mid n_i, \mu_i
  &\sim \operatorname{Binomial}(n_i, \mu_i),\\
\operatorname{logit}(\mu_i)
  &= \beta_0 + \beta_1 \text{canopy}_i + \beta_2 \text{NDVI}_i,\\
E[Y_i / n_i]
  &= \mu_i,\\
\operatorname{Var}(Y_i / n_i)
  &= \mu_i(1 - \mu_i) / n_i.
\end{aligned}
\]

For a 0/1 response, \(n_i = 1\):

```r
drmTMB(
  bf(early_onset ~ canopy + ndvi),
  family = stats::binomial(link = "logit"),
  data = onset_data
)
```

For success counts, keep both the successes and failures:

```r
drmTMB(
  bf(cbind(germinated, failed) ~ treatment + moisture),
  family = stats::binomial(link = "logit"),
  data = seed_trials
)
```

The fitted `mu` is the event probability. `fitted(fit)` and
`predict(fit, dpar = "mu")` return probabilities; multiply by
`successes_i + failures_i` only when the target summary is the expected number
of successes.

Here is a small 0/1 event example. Suppose each site is scored for whether a
vegetation index crosses an early-season threshold. The model asks how canopy
condition and spring NDVI change that event probability:

```{r}
set.seed(196)
n_onset <- 240
onset_data <- data.frame(
  canopy = factor(
    rep(c("open", "closed"), each = n_onset / 2),
    levels = c("open", "closed")
  ),
  ndvi = as.numeric(scale(runif(n_onset, 0.15, 0.85)))
)

closed <- as.numeric(onset_data$canopy == "closed")
mu_onset <- plogis(-0.85 - 0.35 * closed + 1.10 * onset_data$ndvi)
onset_data$early_onset <- rbinom(n_onset, size = 1, prob = mu_onset)

fit_onset <- drmTMB(
  bf(early_onset ~ canopy + ndvi),
  family = stats::binomial(link = "logit"),
  data = onset_data
)

coef(fit_onset, "mu")
```

Predicted probabilities are already on the response scale:

```{r}
new_onset <- data.frame(
  canopy = factor(c("open", "closed"), levels = levels(onset_data$canopy)),
  ndvi = c(0, 0)
)

data.frame(
  canopy = new_onset$canopy,
  ndvi = new_onset$ndvi,
  early_onset_probability = predict(fit_onset, newdata = new_onset, dpar = "mu")
)
```

The original plain-binomial slice was fixed-effect only. The current R engine
also fits ordinary unlabelled `mu` random intercepts and independent numeric
slopes at recovery grade. It still rejects non-logit links, factor responses,
proportions with trial weights, correlated or labelled slopes, structured
effects, bivariate or mixed responses, and `engine = "julia"`.

## Beta-Binomial Equation And Syntax

Use `beta_binomial()` when the response is still successes out of known trials
but the probability varies more than ordinary binomial sampling can explain.
The model uses

\[
\begin{aligned}
Y_i \mid n_i, p_i
  &\sim \operatorname{Binomial}(n_i, p_i),\\
p_i \mid \mu_i, \sigma_i
  &\sim \operatorname{Beta}(\alpha_i, \beta_i),\\
\operatorname{logit}(\mu_i)
  &= \beta_0
     + \beta_1 \text{sheltered}_i
     + \beta_2 \text{moisture}_i,\\
\log(\sigma_i)
  &= \gamma_0 + \gamma_1 \text{sheltered}_i,\\
\phi_i
  &= 1 / \sigma_i^2,\\
\alpha_i
  &= \mu_i\phi_i,\\
\beta_i
  &= (1 - \mu_i)\phi_i,\\
E[Y_i / n_i]
  &= \mu_i,\\
\operatorname{Var}(Y_i / n_i)
  &= \frac{\mu_i(1 - \mu_i)(1 + n_i\sigma_i^2)}
          {n_i(1 + \sigma_i^2)}.
\end{aligned}
\]

The matching `drmTMB` syntax is:

```r
drmTMB(
  bf(cbind(germinated, failed) ~ treatment + moisture, sigma ~ treatment),
  family = beta_binomial(),
  data = seed_trials
)
```

If trays, plots, or subjects are repeatedly observed, the first fitted
mixed-model slice is an ordinary `mu` random intercept, optionally paired with
an independent numeric `mu` slope:

```r
drmTMB(
  bf(cbind(germinated, failed) ~ treatment + moisture +
       (1 | tray) + (0 + moisture | tray),
     sigma ~ treatment),
  family = beta_binomial(),
  data = seed_trials
)
```

The random effects change the logit success-probability predictor. They do not
make `sigma` a random effect, and they do not open correlated slopes or
structured bounded-response effects.

Read each parameter before interpreting the fitted model:

| Symbol or syntax | Meaning | In the seed-germination example |
|---|---|---|
| \(Y_i\), `germinated` | number of successes | seeds that germinated in tray \(i\) |
| \(n_i\), `germinated + failed` | known number of trials | viable seeds planted in tray \(i\) |
| \(\mu_i\), `mu` | expected success probability | germination probability for a seed in that tray |
| \(\boldsymbol{\beta}_{\mu}\) | location coefficients on the logit probability scale | treatment and moisture effects on germination probability |
| \(\sigma_i\), `sigma` | extra-binomial scale | tray-to-tray probability variation beyond binomial sampling |
| \(\phi_i\) | beta precision | `phi_i = 1 / sigma_i^2`, so larger `sigma` means lower precision |
| \(\boldsymbol{\beta}_{\sigma}\) | scale coefficients on the log `sigma` scale | treatment effects on extra-binomial variation |

If `gamma_1` is the sheltered-treatment coefficient, then

\[
\frac{\sigma_\text{sheltered}}{\sigma_\text{open}}
  = \exp(\gamma_1),
\qquad
\frac{\phi_\text{sheltered}}{\phi_\text{open}}
  = \exp(-2\gamma_1).
\]

That second ratio reverses direction because \(\phi = 1 / \sigma^2\). In
`drmTMB`, report `sigma` when you want the direction to mean more or less
modelled variation.

## Seed Germination Example

Suppose seeds are planted in trays under open and sheltered microsites during a
restoration trial. Shelter may increase average germination by reducing heat
stress, and it may also make tray outcomes more predictable by buffering
microclimate. Each tray has a known number of seeds, so the denominator belongs
in the model.

This transparent simulation gives us known structure before fitting:

```{r}
set.seed(197)
n <- 360
seed_trials <- data.frame(
  treatment = factor(
    rep(c("open", "sheltered"), each = n / 2),
    levels = c("open", "sheltered")
  ),
  moisture = as.numeric(scale(runif(n, 0.1, 0.9))),
  trials = sample(18:32, n, replace = TRUE)
)

sheltered <- as.numeric(seed_trials$treatment == "sheltered")
mu_seed <- plogis(-0.55 + 0.70 * sheltered + 0.45 * seed_trials$moisture)
sigma_seed <- exp(-1.15 - 0.35 * sheltered)
phi_seed <- 1 / sigma_seed^2
tray_probability <- rbeta(
  n,
  shape1 = mu_seed * phi_seed,
  shape2 = (1 - mu_seed) * phi_seed
)

seed_trials$germinated <- rbinom(n, size = seed_trials$trials, prob = tray_probability)
seed_trials$failed <- seed_trials$trials - seed_trials$germinated

head(seed_trials)
```

Fit the beta-binomial location-scale model:

```{r}
fit_seed <- drmTMB(
  bf(cbind(germinated, failed) ~ treatment + moisture, sigma ~ treatment),
  family = beta_binomial(),
  data = seed_trials
)
```

Run diagnostics before interpreting coefficients:

```{r}
check_drm(fit_seed)
```

The `mu` coefficients describe the log-odds of germination. The `sigma`
coefficients describe extra-binomial variation in the latent tray-level
success probability:

```{r}
coef(fit_seed, "mu")
coef(fit_seed, "sigma")

sigma_ratio <- exp(coef(fit_seed, "sigma")["treatmentsheltered"])
c(
  sigma_ratio_sheltered_vs_open = sigma_ratio,
  precision_ratio_sheltered_vs_open = sigma_ratio^(-2)
)
```

For a response-scale summary, predict the expected success probability and the
proportion-level standard deviation for trays with the same number of seeds:

```{r}
new_seed_trays <- data.frame(
  treatment = factor(c("open", "sheltered"), levels = levels(seed_trials$treatment)),
  moisture = c(0, 0),
  trials = c(24, 24)
)

mu_hat <- predict(fit_seed, newdata = new_seed_trays, dpar = "mu")
sigma_hat <- predict(fit_seed, newdata = new_seed_trays, dpar = "sigma")
prop_var <- mu_hat * (1 - mu_hat) *
  (1 + new_seed_trays$trials * sigma_hat^2) /
  (new_seed_trays$trials * (1 + sigma_hat^2))

data.frame(
  treatment = new_seed_trays$treatment,
  trials = new_seed_trays$trials,
  expected_probability = mu_hat,
  expected_successes = new_seed_trays$trials * mu_hat,
  sigma = sigma_hat,
  phi = 1 / sigma_hat^2,
  proportion_sd = sqrt(prop_var)
)
```

This table keeps the two parts of the question separate: expected germination
probability and extra-binomial scatter around that probability.

```{r beta-binomial-tray-figure, eval=requireNamespace("ggplot2", quietly = TRUE), fig.width=7.2, fig.height=4.4, fig.cap="Beta-binomial tray summary for the seed-germination example. Faint points are observed tray proportions; overlaid points are fitted expected germination probabilities; vertical bars show plus or minus one fitted proportion-level standard deviation, not confidence intervals.", fig.alt="Jittered point plot of observed germination proportions for open and sheltered trays. Overlaid larger points show fitted expected probabilities, and vertical bars show one fitted proportion-level standard deviation for each treatment."}
library(ggplot2)

seed_trials$observed_proportion <- seed_trials$germinated / seed_trials$trials
seed_plot_summary <- data.frame(
  treatment = new_seed_trays$treatment,
  expected_probability = mu_hat,
  proportion_sd = sqrt(prop_var)
)
seed_plot_summary$lower <- pmax(
  0,
  seed_plot_summary$expected_probability - seed_plot_summary$proportion_sd
)
seed_plot_summary$upper <- pmin(
  1,
  seed_plot_summary$expected_probability + seed_plot_summary$proportion_sd
)

ggplot(seed_trials, aes(treatment, observed_proportion, colour = treatment)) +
  geom_jitter(width = 0.12, height = 0, alpha = 0.18, size = 0.9) +
  geom_errorbar(
    data = seed_plot_summary,
    aes(
      y = expected_probability,
      ymin = lower,
      ymax = upper
    ),
    width = 0.12,
    linewidth = 0.8
  ) +
  geom_point(data = seed_plot_summary, aes(y = expected_probability), size = 3) +
  scale_colour_manual(values = c("open" = "#0072B2", "sheltered" = "#009E73")) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    title = "Denominator-aware proportions can still show raw trays",
    subtitle = "Bars show fitted tray-level scatter, not confidence intervals",
    x = "Microsite treatment",
    y = "Germinated proportion",
    colour = "Treatment"
  ) +
  guides(colour = "none") +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(colour = "grey30")
  )
```

## Strict Continuous Proportions

Use `beta()` when the response is a continuous proportion strictly inside the
open interval. For example, vegetation cover measured from image analysis may
record 0.03, 0.41, or 0.82 rather than successes out of a known number of
quadrats. The beta model uses

\[
\begin{aligned}
Y_i \mid \mu_i, \sigma_i
  &\sim \operatorname{Beta}(\alpha_i, \beta_i),\\
\operatorname{logit}(\mu_i)
  &= \beta_0
     + \beta_1 \text{grazed}_i
     + \beta_2 \text{moisture}_i,\\
\log(\sigma_i)
  &= \gamma_0 + \gamma_1 \text{grazed}_i,\\
\phi_i
  &= 1 / \sigma_i^2,\\
\operatorname{Var}(Y_i)
  &= \frac{\mu_i(1 - \mu_i)\sigma_i^2}{1 + \sigma_i^2}.
\end{aligned}
\]

The fitted syntax is the same one-formula-per-parameter pattern:

```r
drmTMB(
  bf(cover ~ grazing + moisture, sigma ~ grazing),
  family = beta(),
  data = cover_data
)
```

A small transparent example:

```{r}
set.seed(198)
n_cover <- 300
cover_data <- data.frame(
  grazing = factor(
    rep(c("ungrazed", "grazed"), each = n_cover / 2),
    levels = c("ungrazed", "grazed")
  ),
  moisture = as.numeric(scale(runif(n_cover, 0.05, 0.95)))
)

grazed <- as.numeric(cover_data$grazing == "grazed")
mu_cover <- plogis(0.35 - 0.75 * grazed + 0.35 * cover_data$moisture)
sigma_cover <- exp(-1.10 + 0.45 * grazed)
phi_cover <- 1 / sigma_cover^2
cover_data$cover <- rbeta(
  n_cover,
  shape1 = mu_cover * phi_cover,
  shape2 = (1 - mu_cover) * phi_cover
)

fit_cover <- drmTMB(
  bf(cover ~ grazing + moisture, sigma ~ grazing),
  family = beta(),
  data = cover_data
)

check_drm(fit_cover)
```

Here, grazing changes both expected cover and the variation of cover among
sampling units:

```{r}
coef(fit_cover, "mu")
coef(fit_cover, "sigma")

new_cover <- data.frame(
  grazing = factor(c("ungrazed", "grazed"), levels = levels(cover_data$grazing)),
  moisture = c(0, 0)
)

mu_cover_hat <- predict(fit_cover, newdata = new_cover, dpar = "mu")
sigma_cover_hat <- predict(fit_cover, newdata = new_cover, dpar = "sigma")

data.frame(
  grazing = new_cover$grazing,
  expected_cover = mu_cover_hat,
  sigma = sigma_cover_hat,
  phi = 1 / sigma_cover_hat^2,
  cover_sd = sqrt(
    mu_cover_hat * (1 - mu_cover_hat) *
      sigma_cover_hat^2 / (1 + sigma_cover_hat^2)
  )
)
```

Do not use this strict beta route if the observed response includes exact 0 or
1 values. A sprayed plot with 0% damage or a quadrat with 100% cover is a
boundary-generating process, not an interior beta observation. Use
`zero_one_beta()` when those endpoints are structural outcomes.

## Continuous Proportions With Exact Boundaries

Use `zero_one_beta()` when the response is a continuous proportion on `[0, 1]`
and exact 0 or 1 values are generated by a separate boundary process. The
interior observations still use the beta mean-scale model. The two extra
probability formulas describe the boundary process:

\[
\begin{aligned}
\Pr(Y_i = 0) &= zoi_i(1 - coi_i),\\
\Pr(Y_i = 1) &= zoi_i coi_i,\\
\Pr(0 < Y_i < 1) &= 1 - zoi_i,\\
\operatorname{logit}(\mu_i) &= \beta_0 + \beta_1 \text{grazed}_i,\\
\log(\sigma_i) &= \gamma_0 + \gamma_1 \text{grazed}_i,\\
\operatorname{logit}(zoi_i) &= \delta_0 + \delta_1 \text{drought}_i,\\
\operatorname{logit}(coi_i) &= \kappa_0 + \kappa_1 \text{canopy}_i.
\end{aligned}
\]

```r
drmTMB(
  bf(cover ~ grazing, sigma ~ grazing, zoi ~ drought, coi ~ canopy),
  family = zero_one_beta(),
  data = cover_data
)
```

`predict(fit, dpar = "mu")` describes the mean among interior observations.
`predict(fit, dpar = "zoi")` is the probability of an exact boundary outcome,
and `predict(fit, dpar = "coi")` is the probability that a boundary outcome is
exactly 1. `fitted(fit)` returns the unconditional mean
`(1 - zoi) * mu + zoi * coi`, which is usually the response-scale summary an
applied reader expects.

## Current Boundary

The implemented bounded-response path is univariate. Use
`stats::binomial(link = "logit")` for ordinary fixed-effect event probability
models, `beta_binomial()` for overdispersed counted successes and failures,
`beta()` for a single continuous proportion strictly between 0 and 1, and
`zero_one_beta()` for a single continuous proportion with structural exact 0 or
1 values. `beta()` and `beta_binomial()` support ordinary unlabelled `mu`
random intercepts and independent numeric slopes as first mixed-model slices:

```r
drmTMB(
  bf(cbind(successes, failures) ~ treatment + dose +
       (1 | id) + (0 + dose | id),
     sigma ~ treatment),
  family = beta_binomial(),
  data = dat
)
```

```r
drmTMB(
  bf(proportion ~ treatment + dose + (1 | id) + (0 + dose | id),
     sigma ~ treatment),
  family = beta(),
  data = dat
)
```

Do not teach the following as fitted proportion examples yet:

- correlated random slopes, labelled covariance blocks, or random effects in
  `sigma`, `zoi`, or `coi`;
- `sd(group) ~ ...` random-effect scale models;
- `meta_V(V = V)`, or deprecated `meta_known_V(V = V)` for compatibility, with
  beta or beta-binomial responses;
- `phylo()`, `spatial()`, `animal()`, or `relmat()` bounded-response models,
  with one exception: `beta()` already accepts a single unlabelled `animal()`
  term -- an intercept or one-slope on `mu`, or an intercept on `sigma`, one
  endpoint at a time -- but only at recovery grade (trust the point estimate,
  not the interval), so it is held back here on tier grounds, not for want of
  likelihood code;
- bivariate or mixed-response bounded models such as
  `family = c(beta(), gaussian())`;
- ordered beta or beta-binomial zero-inflation;
- denominator shorthand such as `successes / trials`, `weights = trials`, or
  `cbind(successes, trials)` as a replacement for
  `cbind(successes, failures)`.

Most of those are useful future routes that still need likelihood code,
simulation recovery, diagnostics, and [source-map](source-map.html) evidence
before they become tutorial syntax.