--- title: "A goodness-of-fit and calibration toolbox for logistic regression" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{A goodness-of-fit and calibration toolbox for logistic regression} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} 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. ```{r data} 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. ```{r battery} 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: ```{r battery-summary} 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: ```{r grouped} 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: ```{r directed} 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: ```{r ensemble} 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. ```{r edge-plot, fig.height = 3.6} 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. ```{r external} 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: ```{r external-edge, fig.height = 3.6} 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. ```{r localize, fig.height = 3.4} loc <- localize.gof(fit, B = 199, seed = 1) loc plot(loc) ``` ```{r localize-external, fig.height = 3.4} 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. ```{r stream} 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. ```{r} 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.