A goodness-of-fit and calibration toolbox for logistic regression"

In ebrahim.gof: Goodness-of-Fit and Calibration Tests for Logistic Regression

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)
library(ebrahim.gof)

Why a toolbox

Before the inferences or predictions of a logistic regression model can be trusted, its goodness of fit has to be checked. Many tests exist for this, but in practice they are scattered across different packages with different interfaces, base R ships only the Hosmer--Lemeshow idea, and most of the classical tests lose power on sparse data --- the common situation where continuous covariates give almost every observation its own covariate pattern.

ebrahim.gof addresses this in one place. It:

• introduces the author's own sparse-data tests --- the omnibus Ebrahim--Farrington (EF) test, the directed EF / EDGE test that targets smooth calibration-shape departures, and a Cauchy-combination ensemble; and
• aggregates a wide range of classical and modern tests (Hosmer--Lemeshow, McCullagh, Osius--Rojek, le Cessie--van Houwelingen, Stute--Zhu, the binary-adaptive BAGofT test, and the givitiR calibration test) behind one call, run.all.gof(), so they can be compared head to head. Each aggregated test is obtained from its own package and is credited to its authors; they are provided for comparison, not claimed as original here.

The example data

The bundled gof_demo dataset was generated from a model whose true linear predictor is quadratic in age. If we fit a model that treats age linearly, it is therefore subtly misspecified through a smooth, low-dimensional calibration distortion --- exactly the misfit a directed test is built to detect and an omnibus test tends to miss.

data("gof_demo", package = "ebrahim.gof")
str(gof_demo)

fit <- glm(outcome ~ age + bmi + sex + treatment,
           data = gof_demo, family = binomial)

One call: the whole battery

run.all.gof() runs the battery and returns a tidy gof_battery table. Here we use the fast subset (include_slow = FALSE) and install = "no" so nothing is installed or resampled while the vignette builds; tests whose optional package is not present are simply skipped.

set.seed(1)
battery <- run.all.gof(fit, include_slow = FALSE, install = "no")
battery

Read the table by family. The global omnibus statistics (Pearson, deviance) spread their power over hundreds of covariate patterns and tend to miss this smooth misfit; on ungrouped data their sparse forms, and the F-test, are printed for comparison only and are not counted. The directed tests --- EDGE and EDGE.G10 (the cubic basis at two partitions), DEF.poly2, DEF.stukel, DEF.sym, and Stukel's score test --- concentrate their few degrees of freedom on the calibration-curve shape and detect it, as do the EF test and the ensemble. That contrast is the whole point of a directed test.

summary() condenses the table to one line per family, with the number of rejections and the smallest p-value in each:

summary(battery)

Sparse or grouped? Let the package show you both

Three of the battery's tests are pattern-sensitive: Pearson, Deviance, and McCullagh have a reference distribution that depends on whether they are computed per observation (the "sparse" form) or per distinct covariate pattern (the "grouped" form). With all-continuous covariates, as in gof_demo, every observation is its own pattern (800 patterns for 800 observations), so the two forms coincide; the package recognises this degenerate case and omits the redundant (grouped) rows, showing only the sparse form in the table above. (This is also the regime where the per-observation Pearson and deviance statistics lose their chi-squared reference entirely.) The paired (grouped) rows appear only when covariate patterns actually repeat --- as they do next.

The bundled gof_demo_grouped dataset makes the distinction real. It uses the same data-generating process, but age is recorded in 10-year bins and bmi as integers, so the 800 observations share only 328 covariate patterns. The battery detects this automatically and the two forms now disagree --- on this data the sparse deviance says "fine" (p about 0.17) while the grouped deviance rejects (p about 0.04), with the degrees of freedom dropping from 795 to 323:

data("gof_demo_grouped", package = "ebrahim.gof")
fit_g <- glm(outcome ~ age + bmi + sex + treatment,
             data = gof_demo_grouped, family = binomial)
set.seed(1)
run.all.gof(fit_g, include_slow = FALSE, install = "no")

Same model formula, same test family, opposite verdicts --- so handling the grouping wrong silently changes conclusions. The package's answer is to never make that choice for you invisibly: both forms are printed side by side whenever they differ, and in truly sparse data (where no grouped "correction" can rescue Pearson or the deviance) the partition and directed families below are the tools that remain valid.

The directed tests on their own

Each test is also callable directly. The EDGE test projects the grouped standardized residuals onto a small basis of smooth calibration shapes:

edge.gof(fit)                 # directed EDGE test (poly3 basis)
def.gof(fit, basis = "poly2") # a lower-order directed basis
ef.gof(fit)                   # the omnibus EF test

The Cauchy-combination ensemble pools the directed bases into a single p-value:

def.ensemble.gof(fit)

Where the misfit sits: the EDGE picture

plot() on an edge.gof() result draws the standardized residual of each risk group, the shape the test looked for (the curve), a pointwise 95% band (dashed) and a simultaneous band (shaded). Groups beyond the simultaneous band are drawn as triangles. In-sample the residuals are divided by their null standard deviations, which the fit makes smaller than one, and the legend calls the band approximate; for frozen predictions it is exact.

e <- edge.gof(fit)
plot(e)

The curve is below zero in the lowest- and highest-risk groups and above it in between: the model predicts too many events at both ends of the risk range and too few in the middle, the shape that the missing quadratic term in age leaves. scale = "rate" shows the same groups as observed against predicted event rates.

External validation of frozen predictions

A published model is checked on new patients without refitting it. Nothing is estimated from the validation outcomes, so the references are exact rather than approximate. run.all.external() is the battery for this setting. Here the first half of gof_demo plays the development data and the second half the validation data.

dev <- gof_demo[1:400, ]
val <- gof_demo[401:800, ]
covs <- c("age", "bmi", "sex", "treatment")
fit_dev <- glm(outcome ~ age + bmi + sex + treatment, data = dev, family = binomial)
p_val <- predict(fit_dev, newdata = val, type = "response")   # frozen
run.all.external(val$outcome, p_val)

On these 400 validation records none of the tests rejects at 0.05. The EDGE picture shows how far each group is from its prediction, on the event-rate scale:

e_ext <- edge.gof(y = val$outcome, predicted_probs = p_val, external = TRUE)
e_ext
plot(e_ext, scale = "rate")

Localizing the misfit

A rejection says that the model is wrong, not where. localize.gof() (a fitted model) and localize.external() (frozen predictions) split the misfit into parts --- the intercept, the slope, the link and the covariates --- test every combination of them by closed testing, and name the parts that carry misfit with the familywise error rate controlled at alpha. The result gives the update the named part calls for. In-sample the intercept and the slope are fitted, so they cannot be tested there.

The functions draw the verdict only in an interactive session; in a script or a report, call plot() explicitly.

loc <- localize.gof(fit, B = 199, seed = 1)
loc
plot(loc)

loc_ext <- localize.external(val$outcome, p_val, val[, covs], M = 199, seed = 1)
plot(loc_ext)

Monitoring a deployed model

edge.stream() keeps the external EDGE test of a deployed model up to date as outcomes arrive. It stores four sums per risk group, so each update costs the same however many records came before, and the statistic always equals the one computed on all records at once. The groups are fixed in advance from reference predictions.

set.seed(2)
s <- edge.stream(p_ref = p_val, G = 10)
for (day in 1:60) {                       # 100 patients a day
  lp <- rnorm(100, -1, 1)
  shift <- if (day > 30) 0.5 else 0      # from day 31 the event rate drifts upwards
  s <- update(s, rbinom(100, 1, plogis(lp + shift)), plogis(lp))
  if (day %in% c(30, 60)) print(s)
}

The test taken after 30 days sees nothing; the one after 60 days, half of them after the drift, rejects. One test at a fixed time is valid. Testing after every update and stopping at the first rejection needs the level spread over the looks; see ?edge.stream.

Beyond the battery

run.all.gof() selects among the classical tests and the directed ones. Several of the package's procedures sit outside it, because they need something the battery does not have: a penalty, a frozen set of weights, or a scorer of your own. They are called directly on the same fitted model.

Penalised (ridge) fits

Under a penalty the classical references are not merely weak but invalid: shrinkage biases the fitted probabilities, so the grouped residuals acquire a non-centrality and the usual chi-squared reference is wrong. Two functions apply the same correction and differ only in what they refer it to.

set.seed(4)
n <- 300; p <- 8
X <- matrix(rnorm(n * p), n, p)
y <- rbinom(n, 1, plogis(0.2 + X %*% c(0.9, -0.6, 0.4, rep(0, p - 3))))

# closed form: one fit, no resampling, and no Monte Carlo error in the p-value
calm.gof(X, y, lambda = 40)

shrink.gof() computes the same statistics and calibrates them by prepivoting instead, which costs B penalised refits. Prefer calm.gof() unless you need penalize to choose which columns are shrunk, or p >= n, which it refuses. Note the penalty scale: lambda here is n times a glmnet lambda, and calm.gof() takes a lambda_scale argument while shrink.gof() does not.

?calm.gof states where the reference is validated. It is not validated for strongly unbalanced outcomes once p/n is an appreciable fraction, and the function warns when you call it there.

Frozen-weight and pretrained tests

legoft() combines eleven statistics with weights fixed offline and shipped frozen, so two analysts obtain the same p-value; legoft.localize() then says which domain of evidence carries the misfit, with familywise error control. deepgof1() is a pretrained convolutional statistic whose level comes from your own parametric bootstrap rather than from the network, and plot() on its result draws the residual map the network read. None of these is retrained when you call it.

Building your own

gof.features() turns a fitted model into a feature vector of test statistics, and deploy.gof() refers a scorer of yours to a parametric bootstrap under the fitted model.

Going further

• The full battery. Set include_slow = TRUE to add the resampling/refit-based tests (le Cessie--van Houwelingen, Stute--Zhu, BAGofT, GAM-smooth tests) and calibration_plot = TRUE for the givitiR calibration test and belt. These need their respective optional packages; gof_install_suggests() installs them.
• Speed. The EDGE test is calibrated in closed form (no bootstrap, no refit), so it costs a few milliseconds even at large n.
• Choosing a test. For the structured misfit that dominates practice --- a wrong link or an omitted smooth term --- prefer a directed test (EDGE) over an omnibus one; keep an omnibus test in reserve for rough, high-frequency departures.

References

Farrington CP (1996). "On Assessing Goodness of Fit of Generalized Linear Models to Sparse Data." Journal of the Royal Statistical Society B, 58(2), 349--360.

Hosmer DW, Lemeshow S (1980). "Goodness of Fit Tests for the Multiple Logistic Regression Model." Communications in Statistics -- Theory and Methods, 9(10), 1043--1069.

ebrahim.gof documentation built on Oct. 11, 2026, 5:07 p.m.