| Title: | Debiased Score Tests for Goodness of Fit and Model Comparison |
|---|---|
| Description: | Debiased (Neyman-orthogonalized) score tests for assessing whether a semiparametric or parametric regression model is well-specified and for comparing nested models. The test employs a hunt-and-test strategy: on a held-out hunt sample, it fits the null model and uses machine learning to find a direction in which the null model's score seems positive; on an independent test sample, it assesses the significance of the score in the hunted direction. The test employs orthogonalization to eliminate the bias from estimating the null model, yielding a test statistic that is asymptotically standard normal under the null without requiring a parametric form for the alternative. Methods are provided for 'glm', 'lm' and 'mgcv::gam' fits as well as for detecting heterogeneous treatment effects. The methodology is described in Dhawan, Guo and Shah (2026) <doi:10.48550/arXiv.2607.28861>. |
| Authors: | F. Richard Guo [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-2081-7398>), Aditya Dhawan [aut] |
| Maintainer: | F. Richard Guo <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.0 |
| Built: | 2026-09-19 10:46:46 UTC |
| Source: | https://github.com/richardkwo/dscoretest |
Compare two models using the debiased score test
compare_models(fit.0, ...)compare_models(fit.0, ...)
fit.0 |
A fitted null-model object. Methods are provided for
|
... |
Additional arguments passed to the dispatched method,
notably the alternative-model fit |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
compare_models.glm, compare_models.lm,
compare_models.gam,
gof_test, dScoreTest
set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 5 * exp(dat$x1 + dat$x3) + rnorm(n) * 3 fit.0 <- glm(y ~ x1 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 3)) fit.1 <- glm(y ~ x1 + x2 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 4)) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) compare_models(fit.0, fit.1, hunt.style="wls") anova(fit.0, fit.1) # test a misspecified null model: should be rejected fit.00 <- glm(y ~ x2, family = gaussian(link = "log"), data = dat, start = rep(1, 2)) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") plot(compare_models(fit.00, fit.1, hunt.style="wls")) anova(fit.00, fit.1)set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 5 * exp(dat$x1 + dat$x3) + rnorm(n) * 3 fit.0 <- glm(y ~ x1 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 3)) fit.1 <- glm(y ~ x1 + x2 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 4)) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) compare_models(fit.0, fit.1, hunt.style="wls") anova(fit.0, fit.1) # test a misspecified null model: should be rejected fit.00 <- glm(y ~ x2, family = gaussian(link = "log"), data = dat, start = rep(1, 2)) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") plot(compare_models(fit.00, fit.1, hunt.style="wls")) anova(fit.00, fit.1)
Debiased score test of the null fit.0 against the alternative
fit.1, both fitted with mgcv::gam. GLM fit.1 is
used to hunt for signal that fit.0 potentially misses.
## S3 method for class 'gam' compare_models( fit.0, fit.1, hunt.style = "optimal", trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), verbose = FALSE, ... )## S3 method for class 'gam' compare_models( fit.0, fit.1, hunt.style = "optimal", trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), verbose = FALSE, ... )
fit.0 |
The null model as a fitted |
fit.1 |
The alternative model as a fitted |
hunt.style |
Hunting algorithm with the following options.
|
trim.outlier.hunt |
If |
splits |
Numeric vector of length 2 or 3 giving the relative sizes
of the sample splits; rescaled internally to sum to one.
Default is |
verbose |
Default |
... |
Unused; present for S3 generic/method consistency. |
Nesting is checked by predictor-variable name only, not by basis span:
two models with the same predictors but different smooth specifications
(e.g. s(x, k = 5) vs s(x, k = 20)) will pass the check.
The test remains meaningful as long as fit.1's class contains the
relevant alternative directions. Factor predictors, offset() terms,
weights arguments, and multi-column responses are not supported.
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) dat <- mgcv::gamSim(eg=1, n=500, dist="normal", scale=1) dat <- dat[, 1:5] # test fit.0 against fit.1: well-specified (f3 = 0) and should not be rejected fit.0 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x2), data = dat) fit.1 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat) compare_models(fit.0, fit.1) plot(compare_models(fit.0, fit.1)) anova(fit.0, fit.1) # mis-specified model: drops f2 and should be rejected fit.00 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x3), data = dat) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") anova(fit.00, fit.1)set.seed(42) dat <- mgcv::gamSim(eg=1, n=500, dist="normal", scale=1) dat <- dat[, 1:5] # test fit.0 against fit.1: well-specified (f3 = 0) and should not be rejected fit.0 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x2), data = dat) fit.1 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat) compare_models(fit.0, fit.1) plot(compare_models(fit.0, fit.1)) anova(fit.0, fit.1) # mis-specified model: drops f2 and should be rejected fit.00 <- mgcv::gam(y ~ s(x0) + s(x1) + s(x3), data = dat) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") anova(fit.00, fit.1)
Debiased score test of the null model fit.0 against the alternative
fit.1. GLM fit.1 is used to hunt for signal that fit.0
potentially misses.
## S3 method for class 'glm' compare_models( fit.0, fit.1, hunt.style = "optimal", trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), verbose = FALSE, ... )## S3 method for class 'glm' compare_models( fit.0, fit.1, hunt.style = "optimal", trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), verbose = FALSE, ... )
fit.0 |
The null model as a fitted |
fit.1 |
The alternative model as a fitted |
hunt.style |
Hunting algorithm with the following options.
|
trim.outlier.hunt |
If |
splits |
Numeric vector of length 2 or 3 giving the relative sizes
of the sample splits; rescaled internally to sum to one.
Default is |
verbose |
Default |
... |
Unused; present for S3 generic/method consistency. |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 5 * exp(dat$x1 + dat$x3) + rnorm(n) * 3 fit.0 <- glm(y ~ x1 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 3)) fit.1 <- glm(y ~ x1 + x2 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 4)) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) compare_models(fit.0, fit.1, hunt.style="wls") anova(fit.0, fit.1) # test a misspecified null model: should be rejected fit.00 <- glm(y ~ x2, family = gaussian(link = "log"), data = dat, start = rep(1, 2)) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") plot(compare_models(fit.00, fit.1, hunt.style="wls")) anova(fit.00, fit.1)set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 5 * exp(dat$x1 + dat$x3) + rnorm(n) * 3 fit.0 <- glm(y ~ x1 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 3)) fit.1 <- glm(y ~ x1 + x2 + x3, family = gaussian(link = "log"), data = dat, start = rep(1, 4)) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) compare_models(fit.0, fit.1, hunt.style="wls") anova(fit.0, fit.1) # test a misspecified null model: should be rejected fit.00 <- glm(y ~ x2, family = gaussian(link = "log"), data = dat, start = rep(1, 2)) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") plot(compare_models(fit.00, fit.1, hunt.style="wls")) anova(fit.00, fit.1)
Debiased score test of the null fit.0 against the alternative
fit.1. Both are internally refit as Gaussian-family GLMs and the
call is dispatched to compare_models.glm.
## S3 method for class 'lm' compare_models(fit.0, fit.1, ...)## S3 method for class 'lm' compare_models(fit.0, fit.1, ...)
fit.0 |
The null model as a fitted |
fit.1 |
The alternative model as a fitted |
... |
Additional arguments passed to |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 1 + dat$x1 + 2 * dat$x3 + rnorm(n) fit.0 <- lm(y ~ x1 + x3, data = dat) fit.1 <- lm(y ~ x1 + x2 + x3, data = dat) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) anova(fit.0, fit.1) # misspecified model: should be rejected fit.00 <- lm(y ~ x1 + x2, data = dat) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") anova(fit.00, fit.1)set.seed(42) n <- 500 dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n)) dat$x3 <- dat$x3 + (dat$x1 + dat$x2) / 3 dat$y <- 1 + dat$x1 + 2 * dat$x3 + rnorm(n) fit.0 <- lm(y ~ x1 + x3, data = dat) fit.1 <- lm(y ~ x1 + x2 + x3, data = dat) # test fit.0 against fit.1: should not be rejected compare_models(fit.0, fit.1) anova(fit.0, fit.1) # misspecified model: should be rejected fit.00 <- lm(y ~ x1 + x2, data = dat) compare_models(fit.00, fit.1) plot(compare_models(fit.00, fit.1)) compare_models(fit.00, fit.1, hunt.style="wls") anova(fit.00, fit.1)
Transform a hunted function into a debiased function
, where is the
projection of onto the null model.
debias_standard( h.hat, X.debias, fit.debias, predict_fun, weight_fun, wls_method, arg.wls_method )debias_standard( h.hat, X.debias, fit.debias, predict_fun, weight_fun, wls_method, arg.wls_method )
h.hat |
A list as returned by one of |
X.debias |
Part of X for debiasing. |
fit.debias |
Null model fitted on the debiasing sample of X and y. |
predict_fun, weight_fun, wls_method, arg.wls_method
|
They must be compatible with
|
The projection is obtained by fitting the null model
(via wls_method, weighted by weight_fun(fit.debias, X.debias))
with the hunted values h.hat$h(X.debias) as response. This projection
uses all columns of X, even when the hunt itself is driven by
only a subset of covariates (h.hat$X.cols).
A list with elements:
m.h.fitThe null model fitted (over all columns of X)
to project and debias .
hThe debiased hunt function
with signature h(X).
Test whether a semiparametric (e.g., GAM) or parametric (e.g., glm) regression model is well-specified. The test is a debiased (Neyman-orthogonalized) score test computed via sample splitting: on a held-out hunt sample, the null model is fit and a flexible ML algorithm is used to hunt for a direction in which the null model's score seems positive; on an independent test sample, that direction's score is evaluated to assess the significance. The test employs orthogonalization to eliminate plug-in bias from estimating the null model, so the resulting test statistic is asymptotically standard normal under the null without requiring a parametric form for the alternative.
dScoreTest( y, X, score_fun, weight_fun, fit_method, wls_method, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, debias.method = "standard", debias_fun = NULL, trim.outlier.hunt = TRUE, X.cols.hunt = 1:ncol(X), splits = c(0.5, 0.5), arg.fit_method = NULL, arg.wls_method = NULL, arg.hunt_fun = NULL, predict_fun = stats::predict, predict_fun_hunt = NULL, verbose = FALSE )dScoreTest( y, X, score_fun, weight_fun, fit_method, wls_method, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, debias.method = "standard", debias_fun = NULL, trim.outlier.hunt = TRUE, X.cols.hunt = 1:ncol(X), splits = c(0.5, 0.5), arg.fit_method = NULL, arg.wls_method = NULL, arg.hunt_fun = NULL, predict_fun = stats::predict, predict_fun_hunt = NULL, verbose = FALSE )
y |
Numeric response vector of length n. |
X |
Numeric covariate matrix of dimension n x p. |
score_fun |
Function with signature |
weight_fun |
Function with signature |
fit_method |
Function with signature |
wls_method |
Function with signature |
hunt.style |
Hunting algorithm with the following options.
|
hunt.method |
Built-in method for hunting. Currently available:
When this is set to any other value, arguments |
hunt_fun |
Default |
debias.method |
Debiasing method. Currently available:
When set to any other value, |
debias_fun |
Default |
trim.outlier.hunt |
If |
X.cols.hunt |
Integer vector selecting which columns of
|
splits |
Numeric vector of length 2 or 3 giving the relative sizes
of the sample splits; rescaled internally to sum to one.
Default is |
arg.fit_method |
Named list of additional arguments passed to
|
arg.wls_method |
Named list of additional arguments passed to
|
arg.hunt_fun |
Extra arguments (default |
predict_fun |
Function with signature |
predict_fun_hunt |
Default |
verbose |
Default |
For most scenarios, use one of these methods instead:
Use gof_test to test whether a fitted model is
well-specified against a nonparametric alternative. S3 methods are
provided for glm (gof_test.glm), lm
(gof_test.lm) and mgcv::gam
(gof_test.gam).
Use compare_models to test a null model fit.0
against an alternative supermodel fit.1 in the same model class.
Similar to anova, method can be used to conduct a
significance test of one or more predictors.
In contrast with gof_test, this method targets the alternative
fit.1.
S3 methods are provided for glm (compare_models.glm),
lm (compare_models.lm) and mgcv::gam
(compare_models.gam).
Use hte_test_conditional to test treatment effect
heterogeneity.
Use dScoreTest directly for full control over the score,
weight, refit and hunt routines: this is the underlying engine that the
S3 methods wrap.
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
Maintainer: F. Richard Guo [email protected] (ORCID) [copyright holder]
Authors:
Aditya Dhawan [email protected]
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
Useful links:
Report bugs at https://github.com/richardkwo/dScoreTest/issues
plot.dScoreTest, summary.dScoreTest,
hunt_optimal, hunt_wls,
hunt_vanilla, new_dScoreTest
## An example for customizing a dScoreTest: ## Conditional mean independence: is E[Y | X] a function of X[, 1:3] alone, ## i.e. do X4 and X5 carry no further information once X1, X2, X3 are given? set.seed(1) n <- 500 X <- matrix(rnorm(n * 5), n, 5) y <- X[, 1] + X[, 2]^2 + sin(X[, 3]) + rnorm(n) # null TRUE y.alt <- y + X[, 4] * X[, 5] # null FALSE ## Null model class: an arbitrary function of X[, 1:3], fitted by a ## regression forest. The fitter and predict_fun subset to those columns, ## while the hunt still searches all five columns for a direction of ## misspecification. Note honesty = FALSE with tuning: an underfitted null ## model is itself misspecified, and the test then (correctly) rejects on ## that lack of fit rather than on any dependence on X4, X5. fit_method <- function(y, X, ...) grf::regression_forest(X[, 1:3, drop = FALSE], y, honesty = FALSE, tune.parameters = "all") wls_method <- function(y, X, w, ...) grf::regression_forest(X[, 1:3, drop = FALSE], y, sample.weights = w, honesty = FALSE, tune.parameters = "all") predict_fun <- function(fit, X, ...) predict(fit, X[, 1:3, drop = FALSE])$predictions ## Square loss l(f, y) = (y - f)^2 / 2 gives the score l'(f, y) = f - y ## (a negative residual) and the weight l''(f, y) = 1. score_fun <- function(fit, y, X, ...) predict_fun(fit, X) - y weight_fun <- function(fit, X, ...) rep(1, nrow(X)) ## Null holds: no evidence against it. dScoreTest(y, X, score_fun, weight_fun, fit_method, wls_method, predict_fun = predict_fun) ## Null fails: the dependence on X4 and X5 is detected. dScoreTest(y.alt, X, score_fun, weight_fun, fit_method, wls_method, predict_fun = predict_fun)## An example for customizing a dScoreTest: ## Conditional mean independence: is E[Y | X] a function of X[, 1:3] alone, ## i.e. do X4 and X5 carry no further information once X1, X2, X3 are given? set.seed(1) n <- 500 X <- matrix(rnorm(n * 5), n, 5) y <- X[, 1] + X[, 2]^2 + sin(X[, 3]) + rnorm(n) # null TRUE y.alt <- y + X[, 4] * X[, 5] # null FALSE ## Null model class: an arbitrary function of X[, 1:3], fitted by a ## regression forest. The fitter and predict_fun subset to those columns, ## while the hunt still searches all five columns for a direction of ## misspecification. Note honesty = FALSE with tuning: an underfitted null ## model is itself misspecified, and the test then (correctly) rejects on ## that lack of fit rather than on any dependence on X4, X5. fit_method <- function(y, X, ...) grf::regression_forest(X[, 1:3, drop = FALSE], y, honesty = FALSE, tune.parameters = "all") wls_method <- function(y, X, w, ...) grf::regression_forest(X[, 1:3, drop = FALSE], y, sample.weights = w, honesty = FALSE, tune.parameters = "all") predict_fun <- function(fit, X, ...) predict(fit, X[, 1:3, drop = FALSE])$predictions ## Square loss l(f, y) = (y - f)^2 / 2 gives the score l'(f, y) = f - y ## (a negative residual) and the weight l''(f, y) = 1. score_fun <- function(fit, y, X, ...) predict_fun(fit, X) - y weight_fun <- function(fit, X, ...) rep(1, nrow(X)) ## Null holds: no evidence against it. dScoreTest(y, X, score_fun, weight_fun, fit_method, wls_method, predict_fun = predict_fun) ## Null fails: the dependence on X4 and X5 is detected. dScoreTest(y.alt, X, score_fun, weight_fun, fit_method, wls_method, predict_fun = predict_fun)
Returns a fitted CATE where covariates has no effect
modification given the remaining covariates. It employs R-loss and cross fitting
to estimate
fit_CATE( y, X, w = rep(1, nrow(X)), S = 1:(ncol(X) - 1), folds.crossfit = 5, randomized = FALSE )fit_CATE( y, X, w = rep(1, nrow(X)), S = 1:(ncol(X) - 1), folds.crossfit = 5, randomized = FALSE )
y |
Numeric response vector of length n. |
X |
Matrix |
w |
Non-negative numeric weight vector of length n. Defaults to
|
S |
A subset of |
folds.crossfit |
An integer for the number of folds in cross fitting using the R-loss to estimate CATE. When it is 1, no cross fitting is used. |
randomized |
If |
An object of class "CATE":
control_mean_funCATE_fun,
which only depends on through .
SS as specified.
pNumber of covariates, which equals ncol(Z).
## A randomized trial in which the treatment effect depends on Z1 only. set.seed(2) n <- 500 Z <- matrix(rnorm(n * 2), n, 2, dimnames = list(NULL, c("Z1", "Z2"))) Tr <- rbinom(n, 1, 0.5) # randomized treatment tau <- Z[, "Z1"] # true CATE y <- Z[, "Z1"] + Z[, "Z2"] + Tr * tau + rnorm(n) ## S = NULL leaves the CATE unrestricted, so it may depend on Z1 and Z2. fit <- fit_CATE(y, cbind(Tr, Z), S = NULL, randomized = TRUE, folds.crossfit = 2) cor(fit$CATE_fun(Z), tau) # close to 1: the true CATE is recovered sd(fit$CATE_fun(Z)) ## S = 1 bars Z1 from modifying the effect. Because the true CATE depends ## on Z1 alone, the fitted CATE then collapses to nearly a constant. fit0 <- fit_CATE(y, cbind(Tr, Z), S = 1, randomized = TRUE, folds.crossfit = 2) sd(fit0$CATE_fun(Z)) # much smaller than above ## predict() gives the fitted outcome mean mu0(Z) + T * tau(Z). head(predict(fit, cbind(Tr, Z)))## A randomized trial in which the treatment effect depends on Z1 only. set.seed(2) n <- 500 Z <- matrix(rnorm(n * 2), n, 2, dimnames = list(NULL, c("Z1", "Z2"))) Tr <- rbinom(n, 1, 0.5) # randomized treatment tau <- Z[, "Z1"] # true CATE y <- Z[, "Z1"] + Z[, "Z2"] + Tr * tau + rnorm(n) ## S = NULL leaves the CATE unrestricted, so it may depend on Z1 and Z2. fit <- fit_CATE(y, cbind(Tr, Z), S = NULL, randomized = TRUE, folds.crossfit = 2) cor(fit$CATE_fun(Z), tau) # close to 1: the true CATE is recovered sd(fit$CATE_fun(Z)) ## S = 1 bars Z1 from modifying the effect. Because the true CATE depends ## on Z1 alone, the fitted CATE then collapses to nearly a constant. fit0 <- fit_CATE(y, cbind(Tr, Z), S = 1, randomized = TRUE, folds.crossfit = 2) sd(fit0$CATE_fun(Z)) # much smaller than above ## predict() gives the fitted outcome mean mu0(Z) + T * tau(Z). head(predict(fit, cbind(Tr, Z)))
Debiased score test for goodness of fit
gof_test(object, ...)gof_test(object, ...)
object |
A fitted model object. Methods are provided for |
... |
Additional arguments passed to the dispatched method. |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
gof_test.glm, gof_test.lm,
gof_test.gam,
compare_models, dScoreTest
set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) # log(E[y]) ~ X well-specified y0 <- 5 * exp(X[,1] + X[,3]) + rnorm(n) * 3 fit.0 <- glm(y0 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.0) # log(E[y]) ~ X misspecified y1 <- y0 + exp(6 * cos(X[,1]/6)^2) / sqrt(n) fit.1 <- glm(y1 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.1)set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) # log(E[y]) ~ X well-specified y0 <- 5 * exp(X[,1] + X[,3]) + rnorm(n) * 3 fit.0 <- glm(y0 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.0) # log(E[y]) ~ X misspecified y1 <- y0 + exp(6 * cos(X[,1]/6)^2) / sqrt(n) fit.1 <- glm(y1 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.1)
Debiased score test for goodness of fit of an mgcv::gam fit.
## S3 method for class 'gam' gof_test( object, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, trim.outlier.hunt = TRUE, X.cols.exclude = NULL, splits = c(0.5, 0.5), arg.hunt_fun = NULL, predict_fun_hunt = NULL, verbose = FALSE, ... )## S3 method for class 'gam' gof_test( object, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, trim.outlier.hunt = TRUE, X.cols.exclude = NULL, splits = c(0.5, 0.5), arg.hunt_fun = NULL, predict_fun_hunt = NULL, verbose = FALSE, ... )
object |
Fitted |
hunt.style |
Hunting algorithm with the following options.
|
hunt.method |
Built-in method for hunting. Currently available:
When this is set to any other value, arguments |
hunt_fun |
Default |
trim.outlier.hunt |
If |
X.cols.exclude |
Columns in |
splits |
Numeric vector of length 2 or 3 giving the relative sizes
of the sample splits; rescaled internally to sum to one.
Default is |
arg.hunt_fun |
Extra arguments (default |
predict_fun_hunt |
When a customized |
verbose |
Default |
... |
Unused; present for S3 generic/method consistency. |
Only the numeric predictors appearing in stats::model.frame(object)
are exposed to the hunt; X.cols.exclude indexes into these
predictor variables (not basis columns). Factor-by smooths and other
non-numeric predictors are not currently supported. Formulas using
offset() terms, a weights argument, or a multi-column
response (e.g. cbind(succ, fail) ~ ...) are also not supported.
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) dat <- mgcv::gamSim(eg=1, n=500, dist="normal", scale=2, verbose = FALSE) dat.0 <- dat[,1:5] # well-specified fit.0 <- mgcv::gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat.0) test.0 <- gof_test(fit.0) # f3=0, also well-specified fit.1 <- mgcv::gam(y~s(x0)+s(x1)+s(x2),data=dat.0) test.1 <- gof_test(fit.1) plot(test.1) # misspecified dat.1 <- dat.0 dat.1$y <- dat.1$y * dat$f0 fit.2 <- mgcv::gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat.1) test.2 <- gof_test(fit.2) plot(test.2)set.seed(42) dat <- mgcv::gamSim(eg=1, n=500, dist="normal", scale=2, verbose = FALSE) dat.0 <- dat[,1:5] # well-specified fit.0 <- mgcv::gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat.0) test.0 <- gof_test(fit.0) # f3=0, also well-specified fit.1 <- mgcv::gam(y~s(x0)+s(x1)+s(x2),data=dat.0) test.1 <- gof_test(fit.1) plot(test.1) # misspecified dat.1 <- dat.0 dat.1$y <- dat.1$y * dat$f0 fit.2 <- mgcv::gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat.1) test.2 <- gof_test(fit.2) plot(test.2)
Debiased score test for goodness of fit of GLM.
## S3 method for class 'glm' gof_test( object, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, trim.outlier.hunt = TRUE, X.cols.exclude = NULL, splits = c(0.5, 0.5), arg.hunt_fun = NULL, predict_fun_hunt = NULL, verbose = FALSE, ... )## S3 method for class 'glm' gof_test( object, hunt.style = "optimal", hunt.method = "grf", hunt_fun = NULL, trim.outlier.hunt = TRUE, X.cols.exclude = NULL, splits = c(0.5, 0.5), arg.hunt_fun = NULL, predict_fun_hunt = NULL, verbose = FALSE, ... )
object |
Fitted glm object. |
hunt.style |
Hunting algorithm with the following options.
|
hunt.method |
Built-in method for hunting. Currently available:
When this is set to any other value, arguments |
hunt_fun |
Default |
trim.outlier.hunt |
If |
X.cols.exclude |
Columns in |
splits |
Numeric vector of length 2 or 3 giving the relative sizes
of the sample splits; rescaled internally to sum to one.
Default is |
arg.hunt_fun |
Extra arguments (default |
predict_fun_hunt |
When a customized |
verbose |
Default |
... |
Unused; present for S3 generic/method consistency. |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) # log(E[y]) ~ X well-specified y0 <- 5 * exp(X[,1] + X[,3]) + rnorm(n) * 3 fit.0 <- glm(y0 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.0) # log(E[y]) ~ X misspecified y1 <- y0 + exp(6 * cos(X[,1]/6)^2) / sqrt(n) fit.1 <- glm(y1 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.1)set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) # log(E[y]) ~ X well-specified y0 <- 5 * exp(X[,1] + X[,3]) + rnorm(n) * 3 fit.0 <- glm(y0 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.0) # log(E[y]) ~ X misspecified y1 <- y0 + exp(6 * cos(X[,1]/6)^2) / sqrt(n) fit.1 <- glm(y1 ~ X, family = gaussian(link = "log"), start=rep(1,4)) gof_test(fit.1)
Debiased score test for goodness of fit of an lm. Internally
refits the model as a Gaussian-family GLM and dispatches to
gof_test.glm.
## S3 method for class 'lm' gof_test(object, ...)## S3 method for class 'lm' gof_test(object, ...)
object |
Fitted |
... |
Additional arguments passed to |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) X[,3] <- X[,3] + X[,1] + X[,2] / 2 y0 <- 1 + X %*% c(1,1,2) + rnorm(n) # well-specified fit.0 <- lm(y0 ~ X) test.0 <- gof_test(fit.0) plot(test.0) y1 <- y0 + cos(X[,1]) # mis-specified fit.1 <- lm(y1 ~ X) test.1 <- gof_test(fit.1) plot(test.1)set.seed(42) n <- 500 X <- matrix(rnorm(n * 3), nrow = n) X[,3] <- X[,3] + X[,1] + X[,2] / 2 y0 <- 1 + X %*% c(1,1,2) + rnorm(n) # well-specified fit.0 <- lm(y0 ~ X) test.0 <- gof_test(fit.0) plot(test.0) y1 <- y0 + cos(X[,1]) # mis-specified fit.1 <- lm(y1 ~ X) test.1 <- gof_test(fit.1) plot(test.1)
For binary treatment , covariates and outcome , it tests
the hypothesis
where is a subset of covariates and
is the CATE. When is the full set, this amounts to testing
is a constant; when is a proper subset, this amounts to testing
has no further effect modification while holding fixed.
hte_test_conditional( y, Tr, Z, S = 1:ncol(Z), hunt.style = "optimal", folds.crossfit = 5, trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), arg.hunt_grf = list(honesty = FALSE, tune.parameters = "all"), verbose = FALSE, randomized = FALSE )hte_test_conditional( y, Tr, Z, S = 1:ncol(Z), hunt.style = "optimal", folds.crossfit = 5, trim.outlier.hunt = TRUE, splits = c(0.5, 0.5), arg.hunt_grf = list(honesty = FALSE, tune.parameters = "all"), verbose = FALSE, randomized = FALSE )
y |
Numeric response vector of length n. |
Tr |
Binary (0/1) treatment vector of length n. |
Z |
Numeric covariate matrix of dimension n x p. |
S |
A subset of |
hunt.style |
One of |
folds.crossfit |
Number of cross-fitting folds passed to
|
trim.outlier.hunt, splits, verbose
|
Passed through to
|
arg.hunt_grf |
Arguments passed to |
randomized |
If |
An object of class "dScoreTest": a list whose key elements
are the debiased test statistic t.stat and the one-sided p-value
p.val (right tail of the standard normal), along with the test-set
score residuals, the hunted direction, and the call. It has
print,
summary and
plot methods.
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
set.seed(1) n <- 600 Z <- matrix(rnorm(n * 3), n, 3) Tr <- rbinom(n, 1, plogis(Z[, 1])) y <- Z[, 2] + Z[, 3] + Tr * (1 + Z[, 1]) + rnorm(n) # CATE varies with Z1 # allow modification by Z1 (S = {2,3}): well-specified, should not reject hte_test_conditional(y, Tr, Z, S = c(2, 3)) # forbid all modification (constant CATE): misspecified, should reject hte_test_conditional(y, Tr, Z, S = 1:3)set.seed(1) n <- 600 Z <- matrix(rnorm(n * 3), n, 3) Tr <- rbinom(n, 1, plogis(Z[, 1])) y <- Z[, 2] + Z[, 3] + Tr * (1 + Z[, 1]) + rnorm(n) # CATE varies with Z1 # allow modification by Z1 (S = {2,3}): well-specified, should not reject hte_test_conditional(y, Tr, Z, S = c(2, 3)) # forbid all modification (constant CATE): misspecified, should reject hte_test_conditional(y, Tr, Z, S = 1:3)
Hunt by fitting residuals on X optimally, trained by solving a weighted least squares. The hunted function is also pre-debiased. This is the procedure of Section 2.2.1 of Dhawan, Guo and Shah (2026).
hunt_optimal( wls_hunt_method, wls_method, score_fun, weight_fun, fit, y, X, X.cols = 1:ncol(X), binary.y = FALSE, trim.outlier = TRUE, arg.wls_hunt_method = NULL, arg.wls_method = NULL, predict_fun = stats::predict, predict_fun_hunt = stats::predict )hunt_optimal( wls_hunt_method, wls_method, score_fun, weight_fun, fit, y, X, X.cols = 1:ncol(X), binary.y = FALSE, trim.outlier = TRUE, arg.wls_hunt_method = NULL, arg.wls_method = NULL, predict_fun = stats::predict, predict_fun_hunt = stats::predict )
wls_hunt_method |
Function with signature |
wls_method |
Function with signature |
score_fun |
Function with signature |
weight_fun |
Function with signature |
fit |
Fitted null model. Must support |
y |
Response vector of length n. |
X |
Covariates of dim n x p. |
X.cols |
Subset of covariates to hunt. Default |
binary.y |
Set to |
trim.outlier |
If |
arg.wls_hunt_method |
Named list of additional arguments passed to
|
arg.wls_method |
Named list of additional arguments passed to
|
predict_fun |
Function with signature |
predict_fun_hunt |
Function with signature |
If y is binary, then set binary.y=TRUE. Meanwhile,
predict_fun(fit, X, type="response") must output the predicted probabilities.
An object of class "hunt", a list with elements:
hunt.fitThe fitted hunt model produced by
wls_hunt_method.
trim.boundsThe Tukey IQR trimming bounds, or
c(-Inf, Inf) when trim.outlier = FALSE.
predict_fun_huntThe prediction function for
hunt.fit, as supplied.
X.colsThe columns of X used for the hunt, as
supplied.
hA function with signature h(X) giving the
orthogonalized hunted signal.
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
Hunt by fitting residuals on X.
hunt_vanilla( fit_hunt_method, resids, X, X.cols = 1:ncol(X), trim.outlier = TRUE, arg.fit_hunt_method = NULL, predict_fun_hunt = stats::predict )hunt_vanilla( fit_hunt_method, resids, X, X.cols = 1:ncol(X), trim.outlier = TRUE, arg.fit_hunt_method = NULL, predict_fun_hunt = stats::predict )
fit_hunt_method |
Function with signature |
resids |
Residuals (i.e., negative scores) of length n from the null model. |
X |
Covariates of dim n x p. |
X.cols |
Subset of covariates to hunt. (Default: |
trim.outlier |
If |
arg.fit_hunt_method |
Named list of additional arguments passed to
|
predict_fun_hunt |
Function with signature |
An object of class "hunt", a list with elements:
hunt.fitThe fitted hunt model produced by
fit_hunt_method.
trim.boundsThe Tukey IQR trimming bounds, or
c(-Inf, Inf) when trim.outlier = FALSE.
predict_fun_huntThe prediction function for
hunt.fit, as supplied.
X.colsThe columns of X used for the hunt, as
supplied.
hA function with signature h(X) giving the hunted
signal.
Hunt by fitting residuals on X, trained by solving a weighted least squares. See Proposition 2 of Dhawan, Guo and Shah (2026).
hunt_wls( wls_hunt_method, resids, X, X.cols = 1:ncol(X), trim.outlier = TRUE, arg.wls_hunt_method = NULL, predict_fun_hunt = stats::predict )hunt_wls( wls_hunt_method, resids, X, X.cols = 1:ncol(X), trim.outlier = TRUE, arg.wls_hunt_method = NULL, predict_fun_hunt = stats::predict )
wls_hunt_method |
Function with signature |
resids |
Residuals (i.e., negative scores) of length n from the null model. |
X |
Covariates of dim n x p. |
X.cols |
Subset of covariates to hunt. (Default: |
trim.outlier |
If |
arg.wls_hunt_method |
Named list of additional arguments passed to
|
predict_fun_hunt |
Function with signature |
An object of class "hunt", a list with elements:
hunt.fitThe fitted hunt model produced by
wls_hunt_method.
trim.boundsThe Tukey IQR trimming bounds, or
c(-Inf, Inf) when trim.outlier = FALSE.
predict_fun_huntThe prediction function for
hunt.fit, as supplied.
X.colsThe columns of X used for the hunt, as
supplied.
hA function with signature h(X) giving the hunted
signal.
Dhawan, A., Guo, F. R. and Shah, R. D. (2026). The debiased score test: hunt-and-test for semiparametric hypotheses. arXiv:2607.28861. https://arxiv.org/abs/2607.28861
Internal worker that builds a dScoreTest object from a fixed
three-way (or two-way) sample split. Called by dScoreTest.
Splits, hunting algorithm, and predict semantics are taken as fully
resolved arguments — no defaults are inferred from the data.
new_dScoreTest( y, X, idx.hunt, idx.debias, idx.test, score_fun, weight_fun, fit_method, wls_method, hunt.style = "optimal", hunt.method = "customized", debias.method = "standard", debias_fun = debias_standard, fit_hunt_method = NULL, wls_hunt_method = NULL, X.cols.hunt = 1:ncol(X), binary.y = FALSE, trim.outlier.hunt = TRUE, predict_fun = stats::predict, predict_fun_hunt = stats::predict, arg.fit_method = NULL, arg.wls_method = NULL, arg.fit_hunt_method = NULL, arg.wls_hunt_method = NULL )new_dScoreTest( y, X, idx.hunt, idx.debias, idx.test, score_fun, weight_fun, fit_method, wls_method, hunt.style = "optimal", hunt.method = "customized", debias.method = "standard", debias_fun = debias_standard, fit_hunt_method = NULL, wls_hunt_method = NULL, X.cols.hunt = 1:ncol(X), binary.y = FALSE, trim.outlier.hunt = TRUE, predict_fun = stats::predict, predict_fun_hunt = stats::predict, arg.fit_method = NULL, arg.wls_method = NULL, arg.fit_hunt_method = NULL, arg.wls_hunt_method = NULL )
y |
Numeric response vector of length n. |
X |
Numeric covariate matrix of dimension n x p. |
idx.hunt |
Integer indices into |
idx.debias |
Integer indices into |
idx.test |
Integer indices into |
score_fun |
Function with signature |
weight_fun |
Function with signature |
fit_method |
Function with signature |
wls_method |
Function with signature |
hunt.style |
One of |
hunt.method |
String for the hunting method. |
debias.method |
String for the debiasing method, recorded in the
returned object's |
debias_fun |
Function performing the debiasing, with the same signature
as |
fit_hunt_method |
Required when |
wls_hunt_method |
Required when |
X.cols.hunt |
Integer or name vector selecting which columns of
|
binary.y |
Logical. When |
trim.outlier.hunt |
Logical. Passed to the chosen |
predict_fun |
Function with signature |
predict_fun_hunt |
Function with signature |
arg.fit_method, arg.wls_method, arg.fit_hunt_method, arg.wls_hunt_method
|
Named lists of additional arguments forwarded to the corresponding
fitter via |
A list of class "dScoreTest" with elements:
t.statDebiased test statistic
.
p.valOne-sided p-value (right tail of the standard normal).
residsScore residuals on the test subsample.
hOrthogonalized hunted direction on the test subsample.
h.rawHunted direction before the outer debias projection.
hunted_funThe debiased hunted function
, a function that can be applied to X.
DataList with X, y, and the three index
vectors.
CallNamed list of methods,
hunt.style, hunt.method, debias.method, both
predict functions, and the four arg.* lists.
dScoreTest, hunt_optimal,
hunt_wls, hunt_vanilla
Diagonotic plots for the test:
Histogram of , where .
against the index , where refers to the
-th observation in the full dataset. Only those 's in the
test split are drawn. The mean is drawn as a horizontal line.
Extremes values under the null can result in bad normal approximation.
In this case, consider setting trim.outlier.hunt=TRUE.
Residuals (negative scores) versus the hunted signal. A horizontal segment is drawn between each pair of raw hunted signal and the debiased hunted signal. If debiased gets higher, colored in red; otherwise colored in green. A regression line (blue) with a large positive slope indicates the model is misspecified.
Normalized drawn in order.
## S3 method for class 'dScoreTest' plot(x, ...)## S3 method for class 'dScoreTest' plot(x, ...)
x |
A |
... |
Further graphical parameters passed to underlying plotting functions. |
No return value; called for its side effect of producing the diagnostic plots described above.
Computes the outcome mean
for a "CATE" object returned by fit_CATE.
## S3 method for class 'CATE' predict(object, X, ...)## S3 method for class 'CATE' predict(object, X, ...)
object |
A |
X |
Matrix |
... |
Unused, for S3 consistency. |
Numeric vector of length nrow(X) giving
.
Print the score test
## S3 method for class 'dScoreTest' print(x, ...)## S3 method for class 'dScoreTest' print(x, ...)
x |
A |
... |
Unused, for S3 consistency. |
The input x, invisibly. Called for the side effect of
printing a summary of the test to the console.
Score for a fitted CATE model
score_fun(fit, y, X, ...) ## S3 method for class 'CATE' score_fun(fit, y, X, ...)score_fun(fit, y, X, ...) ## S3 method for class 'CATE' score_fun(fit, y, X, ...)
fit |
A fitted model object. |
y |
Numeric response vector. |
X |
Covariate matrix. |
... |
Additional arguments passed to methods. |
Numeric vector of scores.
score_fun(CATE): Score for a "CATE" object: the residual
, i.e. predicted value minus y.
Reports the headline statistic and p-value, the sample-split sizes, a
raw-vs-debiased comparison of the test statistic (so the effect of the
outer-projection debiasing step is visible), and summary
digests of the diagnostic vectors L, L.raw, h,
h.raw and resids.
## S3 method for class 'dScoreTest' summary(object, ...) ## S3 method for class 'summary.dScoreTest' print(x, ...)## S3 method for class 'dScoreTest' summary(object, ...) ## S3 method for class 'summary.dScoreTest' print(x, ...)
object |
A |
... |
Unused, for S3 consistency. |
x |
A |
A list of class "summary.dScoreTest".
Weight for a fitted CATE model
weight_fun(fit, X, ...) ## S3 method for class 'CATE' weight_fun(fit, X, ...)weight_fun(fit, X, ...) ## S3 method for class 'CATE' weight_fun(fit, X, ...)
fit |
A fitted model object. |
X |
Covariate matrix. |
... |
Additional arguments passed to methods. |
Numeric vector of weights.
weight_fun(CATE): Weight for a "CATE" object: constant 1 for
each row of X, since the outcome model has a squared-error loss.