R cross-check — a flexible link for the Bliss beetle data¶

classic links (glm) + estimated flexible links (glmx Gosset, Aranda–Ordaz)¶

Companion to mixprobit_python.ipynb. There is no off-the-shelf Geweke–Keane mixture-of-normals probit in R (that is the from-scratch mixprobit_full_gibbs.py). But the question it answers — is the symmetric probit/logit link wrong, and does a flexible/asymmetric link fit better? — has a classic frequentist counterpart: estimate the link.

  • glm — the textbook probit / logit / cloglog benchmark.
  • glmx::gosset — a symmetric t-shaped link with an estimated tail parameter (does fatter-but-symmetric help?).
  • Aranda–Ordaz (1981) — an asymmetric link family (λ=1 → logit, λ→0 → cloglog), the method designed for exactly this dose–response problem.
In [1]:
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(glmx))

beetle <- read.csv("bliss_beetle.csv")            # Bliss (1935) grouped dose-mortality
dose <- beetle$dose
ntot <- beetle$ntot
kill <- beetle$kill
y <- cbind(kill, ntot - kill)
cat(sprintf("%d beetles, %d killed\n\n", sum(ntot), sum(kill)))

cat("classic link deviances (residual df = 6):\n")
for (lk in c("probit", "logit", "cloglog"))
  cat(sprintf("  %-8s %.2f\n", lk, deviance(glm(y ~ dose, binomial(link = lk)))))

481 beetles, 291 killed


classic link deviances (residual df = 6):

  probit   10.12
  logit    11.23
  cloglog  3.45

Symmetric vs. asymmetric flexibility¶

The probit and logit are symmetric; the beetle mortality curve is skewed. So the question is whether the link needs fatter tails (symmetric) or asymmetry.

  • Gosset (t-link): estimate the degrees of freedom — symmetric, so it can fatten tails but not skew.
  • Aranda–Ordaz: estimate λ — an asymmetric family. We fit (β₀, β₁, λ) jointly by maximum likelihood.
In [2]:
# Gosset: symmetric flexible link, profile the df parameter nu
dev_gosset <- function(nu) tryCatch(deviance(glm(y ~ dose, binomial(link = gosset(nu)))),
                                    error = function(e) 1e6)
og <- optimize(dev_gosset, c(0.2, 60))
cat(sprintf("Gosset (symmetric):   best nu=%.1f  deviance=%.2f  -> symmetric flexibility does NOT help\n",
            og$minimum, og$objective))

# Aranda-Ordaz asymmetric link: P = 1-(1+lam*exp(eta))^(-1/lam)  (lam=1 logit, lam->0 cloglog)
dc <- dose - mean(rep(dose, ntot))
negll <- function(par) {
  eta <- par[1] + par[2] * dc; lam <- exp(par[3]); base <- 1 + lam * exp(eta)
  if (any(base <= 0)) return(1e8)
  P <- pmin(pmax(1 - base^(-1 / lam), 1e-10), 1 - 1e-10)
  -sum(kill * log(P) + (ntot - kill) * log(1 - P))
}
fit <- optim(c(0, 5, 0), negll, method = "Nelder-Mead", control = list(maxit = 5000))
lam <- exp(fit$par[3])
obs <- kill / ntot
ll_sat <- sum(ifelse(kill > 0, kill * log(obs), 0) + ifelse(ntot - kill > 0, (ntot - kill) * log(1 - obs), 0))
dev_ao <- 2 * (ll_sat - (-fit$value))
cat(sprintf("Aranda-Ordaz (asym):  best lambda=%.3f  deviance=%.2f  -> data drives link to the cloglog (asymmetric) end\n",
            lam, dev_ao))

Gosset (symmetric):   best nu=60.0  deviance=10.24  -> symmetric flexibility does NOT help

Aranda-Ordaz (asym):  best lambda=0.000  deviance=3.44  -> data drives link to the cloglog (asymmetric) end
In [3]:
par(mfrow = c(1, 2))

## (a) dose-response curves
dg <- seq(min(dose) - 0.01, max(dose) + 0.01, length = 200)
pf <- glm(y ~ dose, binomial(link = "probit"))
eta_ao <- fit$par[1] + fit$par[2] * (dg - mean(rep(dose, ntot)))
P_ao <- 1 - (1 + lam * exp(eta_ao))^(-1 / lam)
plot(dose, kill / ntot, pch = 19, ylim = c(0, 1), xlab = "log dose (CS2)", ylab = "P(killed)",
     main = "Beetle dose-response")
lines(dg, predict(pf, data.frame(dose = dg), type = "response"), col = "blue", lwd = 2, lty = 2)
lines(dg, P_ao, col = "red", lwd = 2.4)                       # Aranda-Ordaz (= cloglog, since lambda->0)
legend("topleft", bty = "n", lwd = 2, lty = c(NA, 2, 1), pch = c(19, NA, NA),
       col = c("black", "blue", "red"),
       legend = c("observed", "probit (dev 10.1)", "Aranda-Ordaz = cloglog (dev 3.4)"))

## (b) the link's IMPLIED error density (in the y=1[eta+eps>0] convention), standardised
## probit error = N(0,1); the Aranda-Ordaz/cloglog error f_eps(t) = F_AO'(-t)
fAO  <- function(eta) exp(eta) * (1 + lam * exp(eta))^(-(1 / lam + 1))
feps <- function(t) fAO(-t)
gg <- seq(-8, 8, length = 4000); fg <- feps(gg); fg <- fg / sum(fg * (gg[2] - gg[1]))
mu <- sum(gg * fg) * (gg[2] - gg[1]); sdv <- sqrt(sum((gg - mu)^2 * fg) * (gg[2] - gg[1]))
z <- seq(-4, 4, length = 300)
plot(z, dnorm(z), type = "l", col = "blue", lwd = 2, lty = 2, ylim = c(0, 0.55),
     xlab = "disturbance (standardised)", ylab = "density",
     main = "Implied error: the link IS a distribution")
lines(z, sdv * feps(mu + sdv * z), col = "red", lwd = 2.4)    # standardised AO error
abline(v = 0, col = "gray", lwd = 0.6)
legend("topright", bty = "n", lwd = 2, lty = c(2, 1), col = c("blue", "red"),
       legend = c("normal (probit)", sprintf("Aranda-Ordaz (skew %.2f)",
                  sum(((gg - mu) / sdv)^3 * fg) * (gg[2] - gg[1]))))
par(mfrow = c(1, 1))
No description has been provided for this image

Conclusion — R confirms the Python result¶

model deviance
logit 11.23
conventional probit 10.12
Gosset (symmetric flexible) ≈10.2
cloglog 3.45
Aranda–Ordaz (asymmetric flexible, est. λ→0) 3.44
Python mixture-of-normals probit (m=2) 3.40

Two independent flexible approaches agree with the from-scratch Bayesian mixture:

  1. Symmetric flexibility is not enough — the Gosset t-link, free to fatten the tails, stays at ≈10.2 (probit-level). The problem is skew, not tail weight. (This is the frequentist echo of why a scale mixture fails and a full mixture is needed.)
  2. Asymmetric flexibility solves it — Aranda–Ordaz, given freedom over λ, drives the link to the cloglog (asymmetric) end, deviance 3.44 ≈ cloglog 3.45 ≈ the Python mixture 3.40.

So whether you estimate the link (Aranda–Ordaz, frequentist) or estimate the error distribution (Geweke–Keane mixture of normals, Bayesian), the data tells the same story: the beetle response is asymmetric, and accounting for that cuts the deviance from ~10 to ~3.4. In a binary model the error distribution is the link — two languages for one idea.