run.all.gof: Run a Battery of Goodness-of-Fit Tests at Once

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

Run a Battery of Goodness-of-Fit Tests at Once

Description

Runs several goodness-of-fit tests for a binary logistic regression in one call and returns one tidy data.frame, one row per test. Pass a fitted glm to run the whole battery; pass (y, predicted_probs) to run the tests that need only predictions. Each test is wrapped so that a failure of one test never aborts the whole run.

Usage

run.all.gof(
  object,
  predicted_probs = NULL,
  X = NULL,
  tests = "all",
  G = 10,
  include_slow = TRUE,
  parallel = FALSE,
  ncores = NULL,
  calibration_plot = FALSE,
  install = c("ask", "no", "yes"),
  control = list()
)

Arguments

object	A fitted binary logistic glm, with one row per binary record and no prior weights, or a binary (0/1) response vector y (then supply predicted_probs). A grouped (cbind) or weighted fit is refused: expand it to one row per record. Another link gives a warning, since most tests assume the logit. A warning is also given when some fitted risks are 0 or 1 (separation) or when there are fewer than 10 events or non-events.
predicted_probs	Numeric predicted probabilities; required when object is a y vector.
X	Optional design matrix; lets the directed (DEF) and covariate-space tests run from the (y, predicted_probs) form. It should be the model matrix of the fit that gave predicted_probs: a data frame is coded with model.matrix(~ ., X), and a matrix without a constant column gets an intercept column in front.
tests	Either "all" (default) or a character vector of test names to run (e.g. c("EF","DEF.poly3","HL")).
G	Integer number of groups passed to the grouping tests (default 10), or "auto" for max(10, ceiling(n / 25)) (see def.gof). "auto" is resolved once, so every row, the ensemble rows included, uses the same number of groups; the directed rows record it in Note.
include_slow	Logical; when TRUE (the default) the full battery runs, including the slow tests: le Cessie-van Houwelingen smoothing (O(n^2)-O(n^3)), the GAM tests, Stute-Zhu, eHL, BAGofT, Projection, and GiViTI. Set FALSE for a quick run with the fast tests only. A one-time message notes this whenever slow tests are included.
parallel	Logical; when TRUE, the resampling loops of the slow bootstrap tests (Stute-Zhu and Lai-Liu-HL) are run on a local PSOCK cluster via parLapply (works on all platforms, including Windows). All other tests are unaffected. The default FALSE keeps every loop sequential, exactly as in previous versions.
ncores	Integer; the number of worker processes used when parallel = TRUE. The default NULL uses max(1, parallel::detectCores() - 1). Values below 2 fall back to the sequential path.
calibration_plot	Logical; when TRUE and GiViTI is among the tests, also compute and draw the GiViTI calibration belt and store it on the result (retrievable with plot()). Default FALSE.
install	One of "ask" (default), "no", or "yes", controlling what happens when a test in the run needs an optional package that is not installed. In an interactive session, "ask" lists the missing packages and asks before installing, and "yes" installs them without asking; "no" never installs (the test is just skipped with a note). In a non-interactive session (scripts, R CMD check) nothing is ever installed, regardless of this setting. See gof_install_suggests.
control	Optional named list of per-test options. Recognized entries: "Stute-Zhu" = list(B = ...) (bootstrap replicates); GiViTI = list(devel = "internal"/"external"); "Lai-Liu-HL" = list(n0 = ..., k = ..., alpha = ...); Stukel = list(form = "joint"/"lr"/"marginal") (the joint score test by default, the likelihood-ratio refit, or the pre-2.8.0 marginal sum); DEF.poly2, DEF.poly3, DEF.stukel and DEF.sym = list(weights = "unit"/"score", G = ...), where G may be "auto" (see def.gof); and BAGofT = list(...) which forwards to the binary adaptive test – nsim (resampling iterations; default 100), nsplits, ne (the estimation-split size), and the random-forest partitioner's tuning Kmax (maximum number of adaptive partition cells), ntree, nmin, mtry, maxnodes, and engine ("auto", "fast" or "package"; see bagoft.fast). Example: list(BAGofT = list(nsim = 200, Kmax = 8, ntree = 500)); and Projection = list(B = ..., scale = FALSE, max_n = 3000) (bootstrap refits, standardized covariates, and the sample size above which the row is skipped; see projection.gof).

Details

How to read the battery. Every test here answers the same question – does the fitted model describe the data it was fitted to – but they differ in the departure each is built to notice. That is why the panel is more informative than any single p-value: agreement across families is evidence of fit, and disagreement tells you what kind of misfit is present. A test that rejects points at the departure its own construction is sensitive to.

The one thing the panel cannot do is rescue an invalid reference distribution. Under sparse data – almost every covariate pattern unique – the classical chi-square references fail, which is the situation this package was written for.

Every row the battery returns carries a Family label, and the sections below are those labels: find the label in the output, then the section of the same name here.

Global and Standardized statistics (Family "Global", "Standardized"). These compare observed and fitted responses over the whole sample without grouping.

• Pearson – the sum of squared Pearson residuals. Its chi-square reference assumes many observations per covariate pattern. On sparse data that assumption fails and the test is unreliable in both directions; it is reported for completeness and comparison, not for use. It is printed without a significance flag and is not counted among the rejections of the header or of summary(); the Pearson (grouped) row, shown when covariate patterns repeat, is the one to read.

• Deviance – the likelihood-ratio statistic against the saturated model. On ungrouped binary data its expectation depends on the fitted risks alone and not on whether they agree with the responses, so it sits above n - p when the fitted risks are near one half and far below it when they are extreme; the direction of the error is set by the risk profile rather than by the fit. On a correctly specified model with mid-range risks it therefore rejects almost always. It is reported for completeness and comparison, not for use, and like Pearson it is not counted among the rejections.

• Osius-Rojek – rescues the Pearson statistic by standardizing it with its asymptotic mean and variance computed under increasing sample size rather than increasing cell counts, giving a normal reference that stays valid when patterns are unique; the deviate is referred two-sided. Fast, parameter-free, and a sensible default. It can be anti-conservative at small samples and conservative under a strongly skewed covariate.

• McCullagh – standardizes the Pearson statistic by its exact conditional moments rather than asymptotic ones (computed by the Kuss 2002 algorithm). Typically the most powerful of the global statistics. The deviate is referred to the upper tail of the normal, as in the GOFLOGIT macro, so a Pearson statistic below its conditional mean is not counted as misfit. The moments are computed from the p x p information matrix, so the cost is linear in the sample size.

• Copas-RSS – the unweighted sum of squares of the raw residuals, standardized by its own moments and referred two-sided to the normal (df is NA, as for the other normal deviates). Weighting each residual equally rather than by its variance makes it comparatively sensitive to misfit in the middle of the risk range.

• Information-Matrix – the White/Orme test. It compares two estimators of the information matrix that agree only if the model is correctly specified, so it is an omnibus check on the whole specification rather than on calibration alone.

Partition tests (Family "Partition"). These sort observations by fitted risk, group them, and compare observed with expected counts group by group.

• HL – the Hosmer-Lemeshow test, the field's default: G groups cut at percentiles of fitted risk (the "deciles of risk" when G = 10), referred to chi-square on G - 2 degrees of freedom. Note that this reference was established by simulation, not derivation. Familiar and cheap, but modest in power, and its result depends on the grouping.

• HL-equalwidth – the same statistic with groups cut at fixed probability intervals instead of percentiles. When fitted risks are concentrated in a narrow range, intervals can come out empty and the test may not be computable at all, which is why percentile grouping is usually preferred.

• Pigeon-Heyse – a variance-corrected Hosmer-Lemeshow statistic that accounts for the variability of fitted probabilities within each group. The correction is conservative in sparse designs, so a non-rejection carries less weight than the nominal level suggests.

• F-test – deviance residuals compared across the risk groups by a one-way analysis of variance. It is markedly liberal under sparsity and grows more so with sample size; it is included for comparison, printed without a significance flag and not counted among the rejections.

• EF, EF-normal – the Ebrahim-Farrington test, Farrington's correction on the Hosmer-Lemeshow risk groups. EF uses the chi-square reference on G - 2 degrees of freedom at the battery's G; EF-normal uses the plug-in normal reference, which is the reference for many groups, at G = "auto" (groups of about 25 records) unless control = list("EF-normal" = list(G = ...)) sets another number; it needs the model or X. See ef.gof. Tied fitted risks are put in random order.

• Lai-Liu-HL – Lai and Liu's procedure for using the Hosmer-Lemeshow test in large samples, where any test eventually rejects. It standardizes power to a reference sample size n0 and so returns no p-value: the statistic is the standardized power and the accept/reject decision appears in the Note column. Tune with control = list("Lai-Liu-HL" = list(n0 = ..., k = ...)).

• HL-largeN – Nattino, Pennell and Lemeshow's (2020) large-sample Hosmer-Lemeshow test. Instead of perfect fit it tests whether the misfit, measured by \epsilon = \sqrt{\lambda/n}, exceeds a tolerance \epsilon_0: the ordinary statistic is referred to a noncentral \chi^2_{G-2} with noncentrality \epsilon_0^2 n. By their convention \epsilon_0 is the misfit that would be just significant at n0 = 1e6; change it with control = list("HL-largeN" = list(n0 = ...)). The Note gives \epsilon_0 and the estimate \hat\epsilon. Meant for samples in the tens of thousands and above; in small samples it is conservative.

Directed tests (Family "Directed"). Rather than asking whether anything is wrong, these ask whether a particular shape of departure is present, which buys power when the guess is right.

• EDGE, EDGE.G10 – the EDGE test (the cubic basis of edge.gof) at its default partition, G = "auto", and at ten groups, whatever the battery's G. When the default partition has ten groups (n \le 250) the two are one test and only EDGE is shown. The statistic S is a weighted sum of chi-squares: the note says when its p-value comes from S/c \sim \chi^2_{df} (Satterthwaite), so Statistic and df alone do not give it.

• DEF.poly2, DEF.poly3, DEF.stukel, DEF.sym – the directed forms, each aiming the test at a smooth departure in the shape of the calibration curve: a quadratic or cubic drift in the linear predictor, Stukel's asymmetry-and-tail family, or Stukel's symmetric direction, tails too heavy or too light on both sides. Powerful when the misfit resembles the chosen basis, weaker when it does not. Each row takes weights and G through control, for example control = list(DEF.sym = list(weights = "score", G = "auto")); see def.gof. In the full battery DEF.poly3 is left out when it would repeat EDGE.G10 (the battery's G is 10 and control does not move either row); name it in tests to get it. On a sample with no event, or no non-event, there is no fitted model; these rows and the Stukel row are then NA, and Note says why.

• Stukel – a score test against Stukel's generalized logistic link, which nests the logit and lets the two tails bend independently. It is aimed squarely at link misspecification. The two tail directions are tested jointly on 2 degrees of freedom (1 when every fitted risk lies on one side of one half, which Note then says). Testing them jointly gives up a little power against one-sided (cloglog-type) departures. Up to 2.7.0 this row summed two marginal statistics and was liberal; see NEWS. control = list(Stukel = list(form = "lr")) gives the likelihood-ratio test for the same two directions, and form = "marginal" the old sum, for reproducing earlier results only. The likelihood-ratio refit can fail to converge under separation; the row is then NA, with a note. The joint form leaves out a direction whose information after the fit is below 10^{-10} times its information before the fit, as when the fitted logit is constant; when no direction is left the row is NA, with a note.

Covariate-space tests (Family "Covariate-space"). These partition the covariates themselves rather than the fitted risk, so they can see structure that risk-ordering averages away – an omitted interaction, for instance, need not disturb the marginal calibration at all.

• Tsiatis – a score test that adds indicator variables for regions of the covariate space and asks whether they improve the fit.

• Xie – clusters the covariate space and compares observed with expected counts per cluster, using the corrected degrees of freedom G - k/2 - 1 with k the number of predictors.

• Pulkstenis-Robinson – crosses a categorical covariate with risk groups, so it needs at least one categorical predictor. The package auto-detects one (any factor, character or logical, or a numeric with few distinct values, controlled by getOption("ebrahim.gof.pr.maxlev", 6)) and returns NA with a note when there is none.

Smoothing and GAM tests (Family "Smoothing", "GAM"). These replace grouping with a smoother, so nothing is lost to an arbitrary choice of bin edges.

• le-Cessie – the le Cessie-van Houwelingen score test, which smooths the residuals over the covariate space with a kernel and asks whether the smoothed surface is further from zero than chance allows. Sensitive to local structure that omnibus statistics average away, and one of the most expensive tests here, with cost growing sharply in the sample size.

• HL-GAM, PR-GAM, Xie-GAM – variants that fit a deliberately overfitted generalized additive model and use it to define the grouping, letting the data rather than the analyst choose where the boundaries fall. These need mgcv.

Resampling tests (Family "Bootstrap"). When a statistic has no usable closed-form reference, these build one by simulation.

• Stute-Zhu – a cumulative-residual test: residuals are accumulated along the fitted linear predictor and the largest excursion of that path is compared with a parametric bootstrap. It needs no binning and no bandwidth, and in practice it is the best-behaved test in the battery on size. The p-value is (1 + k)/(1 + B), with k the number of bootstrap statistics at least as large as the observed one, so it is never 0. Set the number of resamples with control = list("Stute-Zhu" = list(B = ...)).

• BAGofT – the binary adaptive test, which splits the data, uses one part to learn a partition that separates fitted from observed, and tests on the other. Its behaviour depends strongly on how many splits and resamples it is given; set them with control = list(BAGofT = list(nsim = ...)) and be aware that the published default is far more expensive than a single split. By default the row is computed by bagoft.fast, which returns the same p-value as the BAGofT package for the same seed and needs only randomForest (and dcov above five covariates); control = list(BAGofT = list(engine = "package")) calls BAGofT itself. When no simulated statistic is as extreme as the observed one, BAGofT returns 0; the row then reports 1/(nsim + 1), the smallest Monte Carlo p-value, and says so.

• Projection – the projection test of Escanciano (2006) as defined for logistic regression by Liu et al. (2024): the cumulative residual process is taken along every direction of the covariate space, not only along the fitted linear predictor as in Stute-Zhu, so it also sees departures such as an omitted interaction. Model-based bootstrap with B = 1000 refits by default; see projection.gof. Its weight matrix costs O(n^3) time and O(n^2) memory, so the row is skipped above n = 3000. Set control = list(Projection = list(B = ..., max_n = ...)).

Calibration tests (Family "Calibration"). These come from clinical prediction, and ask directly whether predicted risks match observed frequencies.

• GiViTI, GiViTI-external – the GiViTI polynomial calibration test, which fits a polynomial of the fitted risk and tests whether it departs from the identity, under the internal and external development assumptions respectively. The two are not interchangeable: on data used to fit the model the internal form is the appropriate one, so the full battery runs GiViTI only, and GiViTI-external runs when it is named in tests. For frozen predictions on new data use run.all.external. It also produces the calibration belt, which shows where on the risk scale a model drifts; see calibration_plot. Wraps givitiR in an isolated callr subprocess, so a failure inside its compiled dependencies returns NA instead of ending your session. Select with control = list(GiViTI = list(devel = "internal")).

• eHL – an e-value form of the Hosmer-Lemeshow test, reported as p = \min(1, 1/e). E-values are safe under optional stopping, which conventional p-values are not, but the conversion shown here is conservative.

• Cubic-LR – the cubic calibration likelihood-ratio test: the square and cube of the logit of the fitted risk are added to a logistic regression of the outcome on it, and the drop in deviance is referred to chi-squared on 2 degrees of freedom. It is the goodness-of-link idea of Pregibon (1980) with the cubic logistic model of Morgan (1985), and a fixed-degree relative of the GiViTI polynomial; fast, and needs only the predictions. See cubic.calib.gof, which documents its origin.

Combinations (Family "Ensemble"). Rather than choosing one test, these pool several.

• Ensemble.Vote(3DEF) and Ensemble.Univ(3DEF+EF) – Cauchy combinations of the directed tests, and of those plus the omnibus EF. The Cauchy combination is valid without knowing how the members correlate, which is what makes pooling dependent tests possible at all. The point is to avoid having to guess the departure in advance, at the cost of being slightly less powerful than the single best member would have been. Both rows always combine the unit form of DEF.poly2, DEF.poly3 and DEF.stukel at the battery's G. When control gives those rows other weights or another G, the ensemble rows do not follow, and their Note says "unit form".

• See also legoft, a pretrained combination whose weights are fixed offline and ship frozen, so two analysts running it on the same data obtain the same p-value.

Implementation notes. Tsiatis and Xie cluster the covariate space with k-means using a fixed internal seed, so results are reproducible and your own random stream is left untouched. The equal-frequency groups of HL, F-test, EF and the DEF rows split tied fitted risks by row order, so with many ties (grouped data, or a model on discrete covariates) their results can depend on the order of the rows; randomising the row order is advised. Every bundled test reproduces the implementation used in the original simulation study: Osius-Rojek follows LogisticDx's gof.glm, Copas-RSS follows the rms gof residual, and HL follows ResourceSelection::hoslem.test. The exception is Stukel: its default joint score statistic agrees with anova(..., test = "Rao") for the augmented model, up to glm's convergence tolerance, and only form = "marginal" reproduces LogisticDx (through statmod::glm.scoretest when statmod is installed).

Procedures the battery does not select. Some of the package's own methods are not part of the panel and are called directly on the fitted model, their p-values read beside it: deepgof1 and legoft are of that kind. Naming them in tests does not reach them.

For goodness of fit after penalized fitting, where none of the above references are valid because the coefficients are shrunk, use calm.gof or shrink.gof. These take the design, the response and the penalty rather than a fitted glm, so they are called separately from the battery.

Value

A data.frame (of class gof_battery) with columns Test, Family, Statistic, df, p_value, and Note, one row per test. A dedicated print method shows the rows grouped by family with formatted p-values and significance flags; the underlying columns remain available for programmatic use.

Note

Grouped vs sparse forms. Pearson, Deviance and McCullagh are reported in two forms: the default (sparse / one-trial) form and a "(grouped)" form computed on the distinct covariate patterns (each a Binomial(m_g, P_g)). The two are identical when every covariate pattern is unique (fully sparse data, as in the simulation) and differ only when patterns repeat (m_g > 1). To avoid clutter, the "(grouped)" row is shown only when it actually differs from the sparse form (i.e., when some pattern repeats); on fully sparse data it is a duplicate and is omitted. Osius-Rojek is always computed on covariate patterns, matching its classical (LogisticDx) definition.

Farrington vs EF. The original Farrington (1996) test is a grouped (covariate-pattern) test. The Ebrahim-Farrington (EF) test is its sparse-data counterpart: it does not group by covariate pattern but forms G data-dependent bins of the predicted risk, so it applies directly to fully sparse data. Use EF for sparse binary data; the grouped Farrington form is appropriate only when covariate patterns repeat.

Reproducibility of the parallel path. With parallel = TRUE the cluster's random-number streams are initialized with clusterSetRNGStream, seeded deterministically from the session's current RNG state. Two runs from the same set.seed state (and the same ncores) therefore give identical bootstrap p-values. Note that the parallel L'Ecuyer-CMRG streams necessarily differ from the serial RNG stream, so parallel = TRUE results differ (within Monte-Carlo error) from parallel = FALSE results at the same seed; this is standard and both are valid. Results also depend on ncores, because the replicates are split across workers.

Author(s)

Ebrahim Khaled Ebrahim [e-mail]

References

The aggregated tests are due to their original authors; they are provided here for comparison and credited as follows.

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1111/j.2517-6161.1996.tb02086.x")}

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/03610928008827941")}

McCullagh P (1985). "On the Asymptotic Distribution of Pearson's Statistic in Linear Exponential Family Models." International Statistical Review, 53(1), 61–67. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2307/1402880")}

Osius G, Rojek D (1992). "Normal Goodness-of-Fit Tests for Multinomial Models with Large Degrees of Freedom." Journal of the American Statistical Association, 87(420), 1145–1152. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.1992.10476271")}

le Cessie S, van Houwelingen JC (1991). "A Goodness-of-Fit Test for Binary Regression Models, Based on Smoothing Methods." Biometrics, 47(4), 1267–1282. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2307/2532385")}

Stukel TA (1988). "Generalized Logistic Models." Journal of the American Statistical Association, 83(402), 426–431. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.1988.10478613")}

Stute W, Zhu LX (2002). "Model Checks for Generalized Linear Models." Scandinavian Journal of Statistics, 29(3), 535–545. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1111/1467-9469.00304")}

Escanciano JC (2006). "A Consistent Diagnostic Test for Regression Models Using Projections." Econometric Theory, 22(6), 1030–1051. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1017/S0266466606060506")}

Liu H, Li X, Chen F, Haerdle W, Liang H (2024). "A Comprehensive Comparison of Goodness-of-Fit Tests for Logistic Regression Models." Statistics and Computing, 34, 175. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s11222-024-10487-5")}

Tsiatis AA (1980). "A Note on a Goodness-of-Fit Test for the Logistic Regression Model." Biometrika, 67(1), 250–251. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/67.1.250")}

Xie XJ, Pendergast J, Clarke W (2008). "Increasing the Power: A Practical Approach to Goodness-of-Fit Test for Logistic Regression Models with Continuous Predictors." Computational Statistics & Data Analysis, 52(5), 2703–2713. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.csda.2007.09.027")}

Pulkstenis E, Robinson TJ (2002). "Two Goodness-of-Fit Tests for Logistic Regression Models with Continuous Covariates." Statistics in Medicine, 21(1), 79–93. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/sim.943")}

Nattino G, Finazzi S, Bertolini G (2014). "A New Calibration Test and a Reappraisal of the Calibration Belt for the Assessment of Prediction Models Based on Dichotomous Outcomes." Statistics in Medicine, 33(14), 2390–2407. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/sim.6100")}

Pigeon JG, Heyse JF (1999). "An Improved Goodness of Fit Statistic for Probability Prediction Models." Biometrical Journal, 41(1), 71–82. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/(SICI)1521-4036(199903)41:1<71::AID-BIMJ71>3.0.CO;2-O")}

Copas JB (1989). "Unweighted Sum of Squares Test for Proportions." Journal of the Royal Statistical Society C, 38(1), 71–80. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2307/2347682")}

White H (1982). "Maximum Likelihood Estimation of Misspecified Models." Econometrica, 50(1), 1–25. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2307/1912526")}

Orme C (1988). "The Calculation of the Information Matrix Test for Binary Data Models." The Manchester School, 56(4), 370–376. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1111/j.1467-9957.1988.tb01339.x")}

Kuss O (2002). "Global Goodness-of-Fit Tests in Logistic Regression with Sparse Data." Statistics in Medicine, 21(24), 3789–3801. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/sim.1421")}

Lai X, Liu L (2018). "A Simple Test Procedure in Standardizing the Power of Hosmer-Lemeshow Test in Large Data Sets." Journal of Statistical Computation and Simulation, 88(13), 2463–2472. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/00949655.2018.1467912")}

Nattino G, Finazzi S, Bertolini G (2014). "A New Calibration Test and a Reappraisal of the Calibration Belt for the Assessment of Prediction Models Based on Dichotomous Outcomes." Statistics in Medicine, 33(14), 2390–2407. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/sim.6100")}

Zhang J, Ding J, Yang Y (2023). "Is a Classification Procedure Good Enough? A Goodness-of-Fit Assessment Tool for Classification Learning." Journal of the American Statistical Association, 118(542), 1115–1125. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.2021.1979010")}

Liu Y, Xie J (2020). "Cauchy Combination Test: A Powerful Test with Analytic p-Value Calculation under Arbitrary Dependency Structures." Journal of the American Statistical Association, 115(529), 393–402. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/01621459.2018.1554485")}

Hosmer DW, Hosmer T, le Cessie S, Lemeshow S (1997). "A Comparison of Goodness-of-Fit Tests for the Logistic Regression Model." Statistics in Medicine, 16(9), 965–980. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/(sici)1097-0258(19970515)16:9<965::aid-sim509>3.0.co;2-o")}

The methods introduced by this package, and the studies that evaluate them, are reported in the following. Reproduction materials for each are archived and citable.

Ebrahim EK, Khattab IG, El-Kotory A (2026). "A Modified Hosmer-Lemeshow Goodness-of-Fit Test for Asymmetric Links: Second-Order Power and Robustness." arXiv:2607.15454 [stat.ME]. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2607.15454")}

Ebrahim EK, El-Kotory A (2026). "Benchmarking Goodness-of-Fit and Calibration Algorithms for Logistic Regression Classifiers: A Large-Scale Simulation Study under Sparse Data." Journal of Intelligent Computing and Data Science, 2(2), 97-128. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.21608/jicds.2026.505341.1048")} Reproduction materials: \Sexpr[results=rd]{tools:::Rd_expr_doi("10.5281/zenodo.21286171")}

Ebrahim EK (2026). "Goodness-of-Fit Tests and Calibration Machine-Learning Algorithms for Logistic Regression with Sparse Data." M.Sc. thesis, Alexandria University. arXiv:2608.11140 [stat.ME]. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2608.11140")}

Ebrahim EK, Hussein OAE-A, El-Kotory A (2026). "A Grouped Calibration Test for Logistic Regression That Tolerates a Few Corrupted Records." arXiv:2608.20511 [stat.ME]. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2608.20511")} Reproduction materials (first version): \Sexpr[results=rd]{tools:::Rd_expr_doi("10.5281/zenodo.21247541")}

Ebrahim EK (2026). "EDGES: A Selection-Free Ensemble Goodness-of-Fit Test." Reproduction materials. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.5281/zenodo.21320865")}

Ebrahim EK (2026). "Detection Subspaces: A Theory of Goodness-of-Fit Tests." Reproduction materials. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.5281/zenodo.21687498")}

Ebrahim EK (2026). "Shrinkage Invalidates the Hosmer-Lemeshow Test: Goodness of Fit for Penalized Logistic Regression, with an Application to Glaucoma Diagnosis." arXiv:2609.06413 [stat.ME]. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2609.06413")} Reproduction materials for shrink.gof: \Sexpr[results=rd]{tools:::Rd_expr_doi("10.5281/zenodo.21900114")}

See Also

ef.gof, def.gof, def.ensemble.gof, gof_install_suggests.

Examples

set.seed(1)
n <- 500
x <- runif(n, -3, 3)
y <- rbinom(n, 1, 1 / (1 + exp(-(0.6 * x))))
fit <- glm(y ~ x, family = binomial())

## The fast tests. Every covariate pattern here is unique, so this is the
## sparse case the package is written for.
res <- run.all.gof(fit, include_slow = FALSE)
res

## The return value is a plain data.frame, so the panel can be read
## programmatically as well as printed.
res[res$p_value < 0.05, c("Test", "Family", "p_value")]   # what rejected
table(res$Family)                                        # coverage by family

## A correctly specified model: the panel should mostly agree, and any
## isolated rejection is the false positive you expect at the 5 percent level.
mean(res$p_value < 0.05, na.rm = TRUE)

## Now a model that is genuinely wrong -- the quadratic term is omitted.
y2  <- rbinom(n, 1, 1 / (1 + exp(-(0.6 * x + 0.5 * x^2))))
bad <- glm(y2 ~ x, family = binomial())
run.all.gof(bad, include_slow = FALSE)

## Pick specific tests, for instance one per family, which is the pairing the
## package recommends over relying on any single statistic.
run.all.gof(fit, tests = c("McCullagh", "HL", "Stukel", "Tsiatis"))

## Grouping tests take the number of groups; the choice is a convention
## rather than a derived optimum, so it is worth varying.
for (g in c(5, 10, 20))
  print(run.all.gof(fit, tests = "HL", G = g))


## The full battery (include_slow = TRUE by default). The slow tests need the
## suggested packages mgcv, randomForest, givitiR and callr; in an interactive
## session run.all.gof() offers to install any that are missing
## (install = "ask"). See also gof_install_suggests().
## The control= list forwards options to the individual tests; the reductions
## here keep the example quick without changing what it demonstrates.
run.all.gof(fit, install = "no",
            control = list("Stute-Zhu" = list(B = 50),
                           BAGofT = list(nsim = 20),
                           Projection = list(B = 99)))

## The GiViTI calibration belt shows WHERE on the risk scale a model drifts,
## which a single p-value cannot.
res2 <- run.all.gof(fit, tests = c("McCullagh", "GiViTI"),
                    calibration_plot = TRUE)
plot(res2)   # redraw the stored belt

## Run the bootstrap loops on a PSOCK cluster. Seeds are handled internally,
## so a parallel run reproduces a serial one.
set.seed(1)
run.all.gof(fit, tests = "Stute-Zhu", parallel = TRUE, ncores = 2,
            control = list("Stute-Zhu" = list(B = 50)))

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