def.gof: Directed Ebrahim-Farrington (DEF) Goodness-of-Fit Test

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

Directed Ebrahim-Farrington (DEF) Goodness-of-Fit Test

Description

Performs the Directed Ebrahim-Farrington (DEF) goodness-of-fit test for a fitted binary logistic regression model. DEF concentrates its power on a small set of calibration-curve "shape" directions by projecting the grouped standardized residuals onto a low-dimensional basis and testing the squared length of that projection.

Naming note: this test is published under the name EDGE (Efficient Directed Grouped Examination), and edge.gof is the primary interface going forward. def.gof() is retained, unchanged, as a fully supported legacy name.

Usage

def.gof(
  object,
  predicted_probs = NULL,
  X = NULL,
  G = 10,
  basis = c("poly3", "poly2", "stukel", "sym", "ensemble"),
  method = c("satterthwaite", "imhof"),
  weights = c("unit", "score"),
  external = FALSE
)

Arguments

object	A fitted binary logistic glm, or a binary (0/1) response vector y (then supply predicted_probs).
predicted_probs	Numeric predicted probabilities; required when object is a y vector, ignored when it is a glm.
X	Optional design matrix, used only with the y/predicted_probs form: it enables the exact estimation-adjusted (\Omega) calibration (logit working weights assumed). Without it the conservative \chi^2_k reference is used and a warning is issued. Ignored when object is a glm, and ignored (with a warning) when external = TRUE.
G	Integer number of equal-frequency groups (default 10; must be >= 3), or "auto" for max(10, ceiling(n / 25)), the partition rule of the EDGE paper.
basis	One of "poly3" (default), "poly2", "stukel", "sym", or "ensemble". "sym" is one column, \eta|\eta| at the logit \eta of each group's mean fitted risk: Stukel's (1988) symmetric direction, aimed at tails that are too heavy or too light on both sides, for example a probit or cauchit truth fitted by a logit.
method	One of "satterthwaite" (default) or "imhof". Ignored when weights = "score".
weights	"unit" (default) is the statistic as published, S = r'P_Z r referred to a weighted chi-squared law. "score" multiplies each column by the square root of its group's variance, which for a logit fit makes the statistic the score test for adding the grouped shape to the model (a score-type test for other links). It is referred to chi-squared on the rank of its information matrix, which is the number of columns unless one is redundant (see Details).
external	Logical, default FALSE. TRUE treats the predicted probabilities as frozen (external validation of a fixed model): \Omega = I, a constant column joins the basis, and the statistic is referred to \chi^2 on the number of basis columns (see Details of def.gof). Supply y and predicted_probs, or a glm whose fitted probabilities are then taken as frozen. Requires weights = "unit" and a basis other than "ensemble"; method is not used.

Details

The observations are sorted by predicted probability and split into G equal-frequency groups; the standardized grouped residual vector r is projected onto a basis matrix Z of smooth shapes, giving S = (Z'r)'(Z'Z)^{-1}(Z'r). Its null distribution is a weighted sum of \chi^2_1 variables with weights equal to the eigenvalues of (Z'Z)^{-1}Z'\Omega Z, where \Omega = I - U(X'WX)^{-1}U' is the estimation-adjusted covariance of the grouped residuals. The p-value uses a Satterthwaite scaled-\chi^2 approximation (default) or Imhof's method (if the CompQuadForm package is installed). Bases: "poly2", "poly3" (default), "stukel", "sym"; "ensemble" runs "poly2", "poly3" and "stukel" and combines them via def.ensemble.gof.

Equal-frequency groups split tied fitted risks by row order. With many ties, as with grouped data or a model on discrete covariates, the result can therefore depend on the order of the rows, and randomising the row order is advised.

With weights = "score" each column of Z is multiplied by \sqrt{V_g}, the square root of its group's variance, so that Z'r becomes \sum_g z_g (O_g - E_g): for a logit fit, the score for adding the grouped shape to the model as a step covariate. Its information after adjusting for the fitted coefficients is Z'\Omega Z, and the statistic u'I^{-1}u is referred to a \chi^2 law on the rank of that information (the number of columns unless one is redundant), read from it after scaling to a correlation matrix. For a logit fit this is the Rao score test for adding the grouped columns, and it agrees with anova(..., test = "Rao") up to glm's convergence tolerance; for other links it is a score-type test. A column whose information after the fit is below 10^{-10} times its information before the fit (Z_s'Z_s, with Z_s the weighted columns) is one the model already spans, as when the fitted logit is constant. It is left out, and when no column is left the p-value is NA, with a warning of class def_no_information.

Which weighting to use. The unit form is the statistic as published and is the recommended default. Which weighting is better depends on the basis, and the two bases go opposite ways:

• basis = "poly3" (and the other polynomial bases): use weights = "unit". Over the 32 power scenarios of the EDGE paper the score form was behind in 24 and ahead in 1, by a mean of 0.039.

• basis = "sym": the score form is usually the better choice, and at high discrimination it is decisively so. It was ahead in 21 of the same 32 scenarios; against a probit truth at high discrimination its power was 0.598 against the unit form's 0.096, and 0.834 against 0.506, at nominal size. Weighting each column by \sqrt{V_g} gives the extreme groups more weight, and the symmetric shape \eta|\eta| is carried by those groups; the cubic basis already spans them through its own columns.

Data quality overrides this. The score weighting is the more fragile of the two when a few records carry corrupted predictions, because it leans on exactly the extreme groups such records occupy. With 25 corrupted covariates in 1000 observations and a correct model otherwise, poly3 with weights = "score" raised its false-alarm rate to 0.894 where the unit form reached 0.344. When the data may contain corrupted predictors, use weights = "unit" whatever the basis.

External mode. With external = TRUE the predicted probabilities are taken as frozen, as when a published model is checked on new data with its coefficients fixed. Nothing is estimated from these data, so \Omega = I exactly, and no score equation absorbs the overall level, so a column of ones is added to the basis before redundant columns are dropped. The statistic S = r'P_Z r is then referred to \chi^2 on the number of columns kept: d + 1 for a d-column basis (4 for "poly3" and "stukel", 3 for "poly2", 2 for "sym"; one fewer for "stukel" when every group lies on one side of 0.5). The groups are the same equal-frequency groups as in the default mode. No "conservative" warning is given, since \Omega = I is exact here, and Method is "external". Only weights = "unit" is available: the score form is a different projection, and it has not been validated for frozen predictions. When object is a glm, its response and fitted probabilities are used as the frozen predictions; the fit is not otherwise used, so this is a test of those predictions and not the estimation-adjusted test of the model.

With fewer events (or fewer non-events) than groups, the grouped reference distribution is unreliable. The p-value is still returned, with a warning; a smaller G avoids it. With no event, or no non-event, the model has no maximum-likelihood fit, and the p-value is NA, with a warning of class def_degenerate.

Value

A one-row data.frame with columns Test, Basis, Test_Statistic (the statistic S), df, Method, and p_value. For weights = "score", Method is "score" and df is the integer rank the statistic is referred to. For external = TRUE, Method is "external" and df is the integer number of basis columns, constant included. When basis = "ensemble", the return is that of def.ensemble.gof.

Author(s)

Ebrahim Khaled Ebrahim [e-mail]

References

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, 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")}

See Also

ef.gof, def.ensemble.gof.

Examples

## gof_demo carries a documented smooth calibration misfit: the risk bends in age,
## and a model linear in age misses it. The point of a directed test is to see that.
data("gof_demo", package = "ebrahim.gof")
wrong <- glm(outcome ~ age + bmi + sex + treatment,
             data = gof_demo, family = binomial())
def.gof(wrong)                       # default poly3 basis
def.gof(wrong, basis = "stukel")     # tail-shape basis
def.gof(wrong, basis = "sym")        # symmetric tail direction, one column
def.gof(wrong, weights = "score")    # score form of the poly3 basis
def.gof(wrong, basis = "ensemble")   # combine poly2, poly3 and stukel (CCT)

## give the model the term it was missing, and the same test stands down
right <- glm(outcome ~ poly(age, 2) + bmi + sex + treatment,
             data = gof_demo, family = binomial())
def.gof(right)

## external validation: freeze the model fitted on one half, test it on the other
dev <- gof_demo[1:250, ]; val <- gof_demo[-(1:250), ]
frozen <- glm(outcome ~ age + bmi + sex + treatment, data = dev, family = binomial())
p_val <- predict(frozen, newdata = val, type = "response")
def.gof(val$outcome, predicted_probs = p_val, external = TRUE)

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