Learn / Inside the pipeline / 07 · models

Fitting the decline: probit regression that survives perfect data

Alongside the observed intervals, the workflow fits a curve to each accession×storage series with enough data (≥3 rows, ≥2 distinct years): a binomial GLM of germinated/failed against years, with a probit link. This is the same mathematical family as the Ellis–Roberts viability equation the rest of this site fits — the GLM intercept is Ki, and σ is −1/slope.

The trap is separation: a series that goes 100% then 0% (or clears/fails the same split perfectly) has no finite maximum-likelihood slope — ordinary fitting chases infinity and either errors or returns a meaninglessly steep curve. Small monitoring series hit this constantly.

The fix is bias-reduced estimation (Firth-type penalization): fit the ordinary GLM first, and if it warns of fitted probabilities at 0 or 1 or produces a non-finite likelihood, refit with the penalized method, which always yields finite, usable estimates. Record which method fit each series — that's provenance about estimation quality.

The lab's own fitter (viability.js) matches R's standard implementation to four decimal places on clean data; the bias-reduction fallback is the piece to add when you script this yourself against small, messy series.

Build it yourself
# ordinary probit GLM with a bias-reduced fallback (brglm2)
fit_probit_robust <- function(sub) {
  sub$nongerm <- sub$NumberSown - sub$NumberGermin
  sep <- FALSE
  m <- withCallingHandlers(
    glm(cbind(NumberGermin, nongerm) ~ YearsDifference,
        data = sub, family = binomial(link = "probit")),
    warning = function(w) {
      if (grepl("numerically 0 or 1", conditionMessage(w))) sep <<- TRUE
      invokeRestart("muffleWarning")
    })
  if (!sep && is.finite(logLik(m))) { attr(m, "method") <- "glm"; return(m) }
  m2 <- brglm2::brglm(cbind(NumberGermin, nongerm) ~ YearsDifference,
                      data = sub, family = binomial(link = "probit"),
                      type = "AS_mixed")
  attr(m2, "method") <- "brglm2"   # bias-reduced: finite even under separation
  m2
}
# Ki = coef(m)[1];  sigma_years = -1 / coef(m)[2]
The math behind this step: Fitting the decline: probits, σ, and Ki
The observed interval: no model, just the latest year that cleared the barObserved vs supported: two independent confidence checks