Goodness-of-fit testing in high-dimensional generalized linear models

TL;DR

This paper introduces a flexible family of goodness-of-fit tests for high-dimensional generalized linear models, utilizing residual analysis and machine learning methods to detect model misspecification with controlled error rates.

Contribution

It presents a novel, adaptable testing framework that leverages residual signals and modern machine learning techniques for high-dimensional model validation.

Findings

Effective in detecting model misspecification in simulations

Controls type I error asymptotically under the null hypothesis

Demonstrates utility on real data examples

Abstract

We propose a family of tests to assess the goodness-of-fit of a high-dimensional generalized linear model. Our framework is flexible and may be used to construct an omnibus test or directed against testing specific non-linearities and interaction effects, or for testing the significance of groups of variables. The methodology is based on extracting left-over signal in the residuals from an initial fit of a generalized linear model. This can be achieved by predicting this signal from the residuals using modern flexible regression or machine learning methods such as random forests or boosted trees. Under the null hypothesis that the generalized linear model is correct, no signal is left in the residuals and our test statistic has a Gaussian limiting distribution, translating to asymptotic control of type I error. Under a local alternative, we establish a guarantee on the power of the test. We illustrate the effectiveness of the methodology on simulated and real data examples by testing goodness-of-fit in logistic regression models. Software implementing the methodology is available in the R package `GRPtests'.

Equations

$$\mathbb{E}(Y_i|x_i=x)=:m_0(x)=\mu (x^T{\beta }_0);$$

$$R_i:=\frac{Y_i-\mu (x_i^T\hat{\beta })}{\sqrt{{\mu }^{\prime }(x_i^T\hat{\beta })}},i=1,\dots ,n.$$

$$\hat{\beta }:=\underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{\frac{1}{n}\sum _{i=1}^n\rho (Y_i,x_i^T\beta )+\lambda \Vert \beta {\Vert }_1\right\}.$$

$$R_i=\frac{Y_i-\mu (x_i^T\hat{\beta })}{\sqrt{{\mu }^{\prime }(x_i^T\tilde{\beta })}},i=1,\dots ,n.$$

$$\underset{w\in W}{\sup }{w}^{T}R.$$

$$R_{\text{ora},i}:=\frac{Y_i-\mu (x_i^T\hat{\beta })}{D_{{\beta }_0,ii}},i=1,\dots ,n,$$

$$R_{\text{ora},i}$$

$$r_i\approx D_{{\beta }_0,ii}{x}_i^T({\beta }_0-\hat{\beta }).$$

$$R_{\text{ora}}\approx \epsilon +D_{{\beta }_0}X({\beta }_0-\hat{\beta }).$$

$$w^T{R}_{\text{ora}}\approx w^T\epsilon +w^T{D}_{{\beta }_0}X({\beta }_0-\hat{\beta })=:w^T\epsilon +\delta .$$

$$\Vert X^T{D}_{{\beta }_0}w{\Vert }_{\infty }\le C\sqrt{\log p},$$

$${\hat{\beta }}_{\text{ora-sq}}:=\underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{\frac{1}{\sqrt{n}}\Vert D_{{\beta }_0}(\hat{f}(X)-X\beta ){\Vert }_2+{\lambda }_{\mathrm{sq}}\Vert \beta {\Vert }_1\right\}.$$

$$w_{\text{ora}}:=\frac{D_{{\beta }_0}(\hat{f}(X)-X{\hat{\beta }}_{\text{ora-sq}})}{\Vert D_{{\beta }_0}(\hat{f}(X)-X{\hat{\beta }}_{\mathrm{ora}-\mathrm{sq}}){\Vert }_2},$$

$${\hat{\beta }}_{\mathrm{sq}}:=\underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{\frac{1}{\sqrt{n}}\Vert {\hat{D}}_A(\hat{f}(X)-X\beta ){\Vert }_2+{\lambda }_{\mathrm{sq}}\Vert \beta {\Vert }_1\right\}.$$

$${\hat{w}}_A:=\frac{{\hat{D}}_A(\hat{f}(X)-X{\hat{\beta }}_{\mathrm{sq}})}{\Vert {\hat{D}}_A(\hat{f}(X)-X{\hat{\beta }}_{\mathrm{sq}}){\Vert }_2}.$$

$${\hat{\gamma }}_j:=\underset{\gamma \in {\mathbb{R}}^{p-|G|}}{\mathrm{argmin}}\frac{1}{\sqrt{n}}\Vert \hat{D}(X_j-X_{-G}\gamma ){\Vert }_2+{\lambda }_{\text{nw}}\Vert \gamma {\Vert }_1,$$

$${\hat{w}}_j:=\frac{\hat{D}(X_j-X_{-G}{\hat{\gamma }}_j)}{\Vert \hat{D}(X_j-X_{-G}{\hat{\gamma }}_j){\Vert }_2}.$$

$$T^b:=\underset{j\in G}{\max }\left|\sum _{i=1}^n{\hat{w}}_{j,i}{\hat{R}}_{G,i}{e}_i^b\right|,$$

$$p_{\text{value}}:=\frac{1}{B+1}\biggl{(}1+\sum_{b=1}^{B}\mathds{1}_{\{T^{b}\geq T\}}\biggr{)}.$$

$$\Theta (\lambda ,{\beta }_0,X):=\left\{\vartheta \in {\mathbb{R}}^p:\Vert \vartheta -{\beta }_0{\Vert }_1\le s\lambda ,\Vert X(\vartheta -{\beta }_0){\Vert }_2^2/n\le s{\lambda }^2\right\},$$

$$\delta:=\mathbb{P}\bigl{(}\hat{\beta}\notin\Theta(\lambda,\beta_{0},X)\bigm{|}X\bigr{)}.$$

$$\left|\mathbb{P}\left(T\leq z|{Z_{A}}\right)-\Phi(z)\right|\leq\delta+C\bigl{\{}\lambda_{\mathrm{sq}}\sqrt{n}s\lambda+\|\hat{w}_{A}\|_{\infty}{s}(\lambda^{2}+\lambda_{A}^{2}){n}+K_{X}s\lambda_{A}+\|\hat{w}_{A}\|_{\infty}\bigr{\}},$$

$$\left|\mathbb{P}\left(T\le z|Z_A\right)-\Phi (z)\right|={\mathcal{O}}_P\left(\delta +\frac{s(K_X\log p+K_X^2\sqrt{\log p})}{\sqrt{n}}+\Vert {\hat{w}}_A{\Vert }_{\infty }s{K}_X^2\log p\right).$$

$$({\tilde{\beta }}_{\mathrm{sq}},{\tilde{\eta }}_{\mathrm{sq}}):=\underset{(\beta ,\eta )\in {\mathbb{R}}^p\times {\mathbb{R}}^n}{\mathrm{argmin}}\left\{\frac{1}{\sqrt{n}}\Vert {\hat{D}}_A(\hat{f}(X)-X\beta )-\eta {\Vert }_2+{\lambda }_{\mathrm{sq}}\Vert \beta {\Vert }_1+{\lambda }_{\eta }\Vert \eta {\Vert }_1\right\}$$

$${\tilde{w}}_A:=\frac{{\hat{D}}_A(\hat{f}(X)-X{\tilde{\beta }}_{\mathrm{sq}})-{\tilde{\eta }}_{\mathrm{sq}}}{\Vert {\hat{D}}_A(\hat{f}(X)-X{\tilde{\beta }}_{\mathrm{sq}})-{\tilde{\eta }}_{\mathrm{sq}}{\Vert }_2}.$$

$${\tau }_j^2:=\frac{1}{n}\mathbb{E}\Vert D_{{\beta }_0}(X_j-X_{-G}{\gamma }_{0,j}){\Vert }_2^2,$$

$${\gamma }_{0,j}:=\underset{\gamma \in {\mathbb{R}}^{p-|G|}}{\mathrm{argmin}}\frac{1}{n}\mathbb{E}\Vert D_{{\beta }_0}(X_j-X_{-G}\gamma ){\Vert }_2^2.$$

$$\varepsilon=D_{\beta_{0}}^{-1}\bigl{(}Y-\mu(X\beta_{0})\bigr{)}.$$

$${\Theta }_{-G}(\lambda ,{\beta }_0):=\left\{\vartheta \in {\mathbb{R}}^{p-|G|}:\Vert \vartheta -{\beta }_{0,-G}{\Vert }_1\le s\lambda ,\Vert X_{-G}(\vartheta -{\beta }_{0,-G}){\Vert }_2^2/n\le s{\lambda }^2\right\}.$$

$$\mathbb{E}(e^{{\eta }_i^2/c_1^2}|X)\le c_2\text{ and }\mathbb{E}({\eta }_i^2|X)\ge c_3,$$

Full text

1 Introduction

In recent years, there has been substantial progress in developing methodology for estimation in generalized linear models in high-dimensional settings, where the number of covariates in the model may be much larger than the number of observations. A standard technique for estimation is the Lasso for generalized linear models (Park and Hastie,, 2007), which has a fast implementation in the R package glmnet (Friedman et al.,, 2010) and is widely used. The Lasso enjoys good empirical and theoretical properties for estimation and variable selection, provided that we are searching for a sparse approximation to the regression coefficients in the generalized linear model.

Once a generalized linear model has been fitted to the high-dimensional data, it is important to assess the quality of the fit. Literature on testing goodness-of-fit in low-dimensional settings is extensive: we refer to Section 1.2 below for an overview. However, the methods typically rely on properties that only hold in low-dimensional settings such as asymptotic linearity and normality of the maximum likelihood estimator, for example. These may fail to hold with an increasing number of covariates in the model; as a consequence it is typically not possible to extend these approaches in an obvious way to the high-dimensional setting. This motivates us to develop a new method that may be used for detecting misspecification in the fit of a (potentially high-dimensional) generalized linear model.

To fix ideas, suppose we have data $(x_{i},Y_{i})_{i=1}^{n}$ formed of feature vectors $x_{i}\in\mathbb{R}^{p}$ and univariate responses $Y_{i}\in\mathcal{Y}\subseteq\mathbb{R}$. Let us write $X=[x_{1},\dots,x_{n}]^{T}=[X_{1},\dots,X_{p}]=(X_{ij})_{1\leq i\leq n,1\leq j\leq p}$ for the $n\times p$ design matrix and $Y=(Y_{1},\dots,Y_{n})^{T}$ for the vector of responses. Consider a generalized linear model (McCullagh and Nelder,, 1989) for the data. Specifically, consider the setting where the $Y_{i}$ are independent conditional on $X$, and the conditional distribution of $Y_{i}$ only depends on $X$ through the linear combination $x_{i}^{T}\beta_{0}$ for some coefficient vector $\beta_{0}\in\mathbb{R}^{p}$. In particular, for the conditional expectation, this implies a structure of the form

[TABLE]

the function $\mu(\cdot)$ is a known inverse link function and $\beta_{0}$ is unknown. Moreover, we assume that $\text{var}(Y_{i}|x_{i}=x)=\mu^{\prime}(x^{T}\beta_{0}).$ This structure of the variance arises in generalized linear models derived from exponential families with canonical link functions, such as logistic regression or Poisson log-linear models.

We will focus on the detection of misspecification in the conditional mean function. In a low-dimensional setting, we understand that the model is misspecified in the conditional mean when there does not exist a $\beta_{0}\in\mathbb{R}^{p}$ such that $m_{0}(x)=\mu(x^{T}\beta_{0})$. In a high-dimensional setting where $p\geq n$, this concept becomes more complicated at first sight; for example, with fixed design points $x_{1},\dots,x_{n}$, there always exists $\beta_{0}\in\mathbb{R}^{p}$ such that $m_{0}(x_{i})=\mu(x_{i}^{T}\beta_{0})$ for all $i=1,\dots,n$, meaning that the model can never be misspecified. However, in a high-dimensional setting, it is impossible to estimate $\beta_{0}$ consistently without additional structural assumptions. An assumption that is often used, and which we adopt in this paper, is sparsity of the model. Therefore, we address the question of whether a sparse model fits well to the observations, or whether a (sparse) non-linear model is more appropriate. If we restrict ourselves to sparse models, then misspecification can happen in the same way as in low-dimensional settings, even for fixed design. Some of the most important types of misspecification that are of interest in applications are missing nonlinear terms such as quadratic effects or interaction terms. Examples of generalized linear models that may be covered by our framework include logistic regression, Poisson regression, robust regression (Huber loss, Cauchy loss) and linear regression.

1.1 Overview of our contributions

We now briefly outline our strategy for goodness-of-fit testing; a more detailed description is given in Section 2. Let $\hat{\beta}$ be an estimate of $\beta_{0}$ derived from a Lasso-penalised generalized linear model (GLM Lasso) fit. Our starting point is the vector ${R}$ of Pearson residuals, with $i$-th coordinate

[TABLE]

Now consider taking as a test statistic the scalar product $w^{T}{R}$, for some (fixed) unit vector $w\in\mathbb{R}^{n}$. If the generalized linear model were correct, then $w^{T}{R}$ would be approximately an average of zero-mean random variables, and under reasonable conditions, should converge to a centred Gaussian random variable. On the other hand, if the model were misspecified, the residuals would contain some signal, and were $w$ to be positively correlated with this signal, the lack of fit should be exposed by the test statistic taking a large value.

In the alternative setting, the signal in the residuals may be picked up by more flexible regression methods, such as random forests (Breiman,, 2001) or boosted trees (Chen and Guestrin,, 2016). However using such flexible regressions to inform the choice of $w$ directly would make $w$ strongly dependent on ${R}$ even under the null; as such calibration of the resulting test statistic would be problematic. Our approach therefore is to construct $w$ based on an independent auxiliary dataset $(X_{A},Y_{A})$ (e.g. derived through sample splitting) in the following way. We first perform a GLM Lasso fit on the auxiliary dataset to obtain an additional set of residuals. Regressing these residuals back on to the explanatory variables $X_{A}$ using a flexible regression method, we obtain an estimated regression function $\tilde{f}:\mathbb{R}^{p}\to\mathbb{R}$ that aims to predict the signal in the residuals; we refer to the $n$-fold concatenation of such an $\tilde{f}$ as a residual prediction function $\hat{f}:\mathbb{R}^{n\times p}\rightarrow\mathbb{R}^{n}$. We may then choose $w$ proportional to $\hat{f}(X)$ to give a direction $w$ independent of $Y$.

One important issue that arises in the high-dimensional setting is that although components of ${R}$ are close to zero-mean under the null, their bias can drive a substantial shift in the mean of $w^{T}{R}$. To prevent this, we replace $w$ with the residuals from a particular weighted (square-root) Lasso regression of $w$ on to $X$. This final step ensures that $w$ is almost orthogonal to the bias in the residuals and as a consequence, the limiting distribution under the null is a centred Gaussian. A notable feature of our construction is that the asymptotic null distribution is essentially invariant to the residual prediction method used. This can therefore be as flexible as needed to detect the type of mean misspecification we would like to uncover.

We provide a software implementation of our methodology in the R package GRPtests (Janková et al.,, 2019).

1.2 Related literature

High dimensions. Our work is related to that of Shah and Bühlmann, (2018) who study goodness-of-fit tests for the linear model. They consider test statistics based on a proxy for the prediction error of a flexible regression method applied to the scaled residuals following a square-root Lasso fit to the data. It is shown that when a Gaussian linear model holds, these residuals depend only weakly on the unknown regression coefficients, motivating calibration of the tests via a parametric bootstrap. As there is no analogue of this result for other generalized linear models, it does not seem possible to extend this approach to our more general setting. Our methodology shares the idea of ‘predicting’ the residuals but, even when we specialize our approach to the Gaussian linear model, differs substantially in the construction of test statistics and the form of calibration.

In recent years, there has been much work on inference and testing in high-dimensional generalized linear models, particularly for the linear model. The work on significance testing includes flexible approaches based on (multiple) sample-splitting (Wasserman and Roeder,, 2009; Meinshausen et al.,, 2009; Meinshausen and Bühlmann,, 2010; Shah and Samworth,, 2013) which may be combined with other methods. Another line of work, initiated by Zhang and Zhang, (2014), proposes a method of de-biasing the Lasso that can be used for testing significance of variables in the linear regression. van de Geer et al., (2014) extend the methodology to generalized linear models; further developments include Javanmard and Montanari, (2014), Dezeure et al., (2017) and Yu et al., (2018); see also Belloni et al., (2014). General frameworks for testing low-dimensional hypotheses about the parameter can be based for example on Neyman orthogonality conditions (Chernozhukov et al.,, 2015, 2018) or on a profile likelihood testing framework (Ning et al.,, 2017). In recent work, Zhu and Bradic, (2017) propose a method for testing more general hypotheses about the parameter vector, such as the sparsity level of the model parameter and minimum signal strength. Javanmard and Lee, (2017) suggest a procedure to test similar hypotheses about the parameter in linear or logistic regression.

Low dimensions. There are numerous methods for testing goodness-of-fit of a model in low-dimensional settings, especially for the case of logistic regression, which is one of the focuses of this work. The most standard tests are residual deviance and Pearson’s chi-squared tests; however, they behave unsatisfactorily if the data contain only a small number of observations for each pattern of covariate values. There have been a number of strategies to circumvent this difficulty, mainly based on grouping strategies, residual smoothing or modifications of Pearson’s chi-squared test.

Hosmer and Lemeshow, (1980) proposed two methods of grouping based on ranked estimated logistic probabilities that form groups of equal numbers of subjects. The disadvantage of these tests (as noted in Le Cessie and Van Houwelingen, (1991)) is that as they are based on a grouping strategy in the space of responses, they lack power to detect departures from the model in regions of the covariate space that yield the same estimated probabilities. For example, a model with a quadratic term may have very different covariate values with the same estimated probability. Tsiatis, (1980) circumvents the difficulties faced by Hosmer–Lemeshow tests using a grouping strategy in the covariate space. However, different partitions of the space of covariates may still lead to substantially different conclusions.

Le Cessie and Van Houwelingen, (1991) introduced a test based on residual smoothing using nonparametric kernel methods. Smoothed residuals replace each residual with a weighted average of itself and other residuals that are close in the covariate space. If residuals close to each other are strongly correlated, smoothing does not affect the magnitude of the residuals strongly, while if they are not correlated smoothing will shrink the residuals towards zero. Su and Wei, (1991) proposed a goodness-of-fit test for the generalized linear model based on a cumulative sum of residuals, which was later adapted by Lin et al., (2002) and a weighted version was proposed in Hosmer and Hjort, (2002). Another approach based on modifications of Pearson’s chi-squared test was studied in Osius and Rojek, (1992) and Farrington, (1996) who derived a large-sample normal approximation for Pearson’s chi-squared test statistic.

1.3 Organization and notation

The rest of the paper is organised as follows. In Section 2 we motivate and present our goodness-of-fit testing methodology. In Section 3, we study its theoretical properties, providing guarantees on the type I error and power. In Section 4, we illustrate the empirical performance of the method on simulated and semi-real genomics data. A brief discussion is given in Section 5. Proofs are deferred to Section 6 and Appendix A.

For a vector $x\in\mathbb{R}^{d}$, we let $x_{j}$ denote its $j$-th entry and write $\|x\|_{p}:=(\sum_{j=1}^{d}|x_{j}|^{p})^{1/p}$ for $p\in\mathbb{N}$, $\|x\|_{\infty}:=\max_{j=1,\dots,d}|x_{j}|$ and $\|x\|_{0}$ for the number of non-zero entries of $x.$ For a matrix $A\in\mathbb{R}^{n\times p}$, we use the notation $A_{ij}$ or $(A)_{ij}$ for its $(i,j)$-th entry, $A_{j}$ to denote its $j$-th column and we let $\|A\|_{\infty}:=\max_{i,j}|A_{ij}|$. Letting $G\subseteq\{1,\dots,p\}$, we denote by $A_{G}$ the matrix containing only columns from $A$ whose indices are in $G$, and by $A_{-G}$ the columns of $A$ whose indices are in the complement of $G$. We use $\Lambda_{\min}(A)$ and $\Lambda_{\max}(A)$ to denote the minimum and maximum eigenvalue of a square matrix $A$.

For sequences of random variables $X_{n},Y_{n}$, we write $X_{n}=\mathcal{O}_{P}(Y_{n})$ if $X_{n}/Y_{n}$ is bounded in probability and $X_{n}=o_{P}(1)$ if $X_{n}$ converges to zero in probability. We write $a\lesssim b$ to mean that there exists $C>0$, which may depend on other quantities designated as constants in our assumptions, such that $a\leq Cb$. If $a\lesssim b$ and $b\lesssim a$, we write $a\asymp b$. Finally, for a function $f:\mathbb{R}\rightarrow\mathbb{R}$ and a vector $z=(z_{1},\dots,z_{n})\in\mathbb{R}^{n}$ we will use $f(z)$ to denote the coordinate-wise application of $f$ to $z$, that is $f(z)=(f(z_{1}),\dots,f(z_{n}))$.

2 Methodology: Generalized Residual Prediction tests

As discussed in Section 1.1, our Generalized Residual Prediction (GRP) testing methodology relies on an initial fit of the Lasso for generalized linear models, which is defined by

[TABLE]

Here $\rho:\mathcal{Y}\times\mathbb{R}\rightarrow\mathbb{R}$ is a loss function, usually derived from the negative log-likelihood associated with the model. Our general framework for goodness-of-fit testing will also assume we have available an auxiliary dataset $(X_{A},Y_{A})\in\mathbb{R}^{n_{A}\times p}\times\mathcal{Y}^{n_{A}}$ independent of $(X,Y)$, sharing the same conditional distribution structure as that of $(X,Y)$. In the rest of the paper, we take $n_{A}=n$ for simplicity, although this is not needed for our procedures. Consider the Pearson-type residuals

[TABLE]

Here $\tilde{\beta}\in\mathbb{R}^{p}$ is an additional estimate of $\beta_{0}$ that may be computed using the auxiliary dataset, or in certain circumstances may be taken as $\hat{\beta}$ itself: we discuss these two cases in the following sections. Given the vector ${R}$ of residuals, the basic form of our test statistic is $w^{T}{R}$; here $w\in\mathbb{R}^{n}$ is a direction typically derived using the auxiliary dataset. We describe in detail the construction of such a $w$ in Section 2.1, where the goal is general goodness-of-fit testing.

A further modification of the method can allow us to use multiple directions $w$ to test simultaneously for different departures from the null or to aggregate over different directions derived using flexible regression methods with different tuning parameters. Given a set $W\subseteq\mathbb{R}^{n}$ of direction vectors $w$, our proposed test statistic then takes the form

[TABLE]

We illustrate the use of this more general form of our test statistic for testing the significance of a group of variables. Such a problem may not immediately seem like goodness-of-fit testing, but is equivalent to testing the adequacy of a model not involving the group of variables in question. We explain how this may be addressed by our framework in Section 2.2 and describe a wild bootstrap procedure (Chernozhukov et al.,, 2013) to approximate the distribution of the test statistic under the null.

2.1 Goodness-of-fit testing

To motivate our general procedure for goodness-of-fit testing, consider the vector ${R}_{\text{ora}}$ of Pearson residuals with an oracle variance scaling, whose $i$-th component is given by

[TABLE]

where $D_{\beta_{0},ii}^{2}:={\mu^{\prime}(x_{i}^{T}\beta_{0})}$. We may decompose the residuals into noise and estimation error terms by writing

[TABLE]

where \varepsilon_{i}:=\bigl{\{}Y_{i}-\mu(x_{i}^{T}\beta_{0})\bigr{\}}/D_{\beta_{0},ii} and r_{i}:=\bigl{\{}\mu(x_{i}^{T}\beta_{0})-\mu(x_{i}^{T}\hat{\beta})\bigr{\}}/D_{\beta_{0},ii}. If the generalized linear model is correct, then $\mathbb{E}(\varepsilon_{i}|x_{i})=0$ and $\mathrm{Var}(\varepsilon_{i}|x_{i})=1$. Turning to the remainder term $r_{i}$, a first-order Taylor expansion of $\mu$ yields the approximation

[TABLE]

Writing $D_{\beta_{0}}$ for the diagonal matrix with entries $D_{\beta_{0},ii}$ for $i=1,\dots,n$, and $\varepsilon:=(\varepsilon_{1},\dots,\varepsilon_{n})$, we obtain the decomposition

[TABLE]

Consider a unit vector $w\in\mathbb{R}^{n}$; as discussed in Section 1.1, this will typically be constructed from an application of a residual prediction method on the auxiliary data. Our oracle $w^{T}{R}_{\text{ora}}$ then satisfies

[TABLE]

Under suitable conditions on $w$ and on the moments of the errors, the Berry–Esseen theorem should ensure that the pivot term $w^{T}\varepsilon$ is well approximated by a standardised Gaussian random variable. To keep the remainder term $\delta$ in (3) under control we can leverage the fact that under the null, we can expect $\|\hat{\beta}-\beta_{0}\|_{1}$ to be small. If $w$ satisfies a near-orthogonality condition

[TABLE]

for some $C>0$, then Hölder’s inequality will yield $|\delta|\leq C\sqrt{\log p}\|\hat{\beta}-\beta_{0}\|_{1}$, which asymptotically vanishes under suitable conditions on the sparsity of $\beta_{0}.$

To guarantee the near-orthogonality condition (4), we may use the square-root Lasso (Belloni et al.,, 2011; Sun and Zhang,, 2012): for $\lambda_{\mathrm{sq}}>0$, let

[TABLE]

The Karush–Kuhn–Tucker (KKT) conditions for the convex programme imply that the resulting vector of scaled residuals,

[TABLE]

satisfies the near-orthogonality property $\|X^{T}D_{\beta_{0}}w_{\text{ora}}\|_{\infty}\leq C\sqrt{\log p}$ when $\lambda_{\mathrm{sq}}=C\sqrt{\log p/n}$. Note that in performing this square-root Lasso regression, we are not assuming that $\hat{f}$ is well-approximated by a sparse linear combination of variables: we are simply exploiting the stationarity properties of the solution to the square-root Lasso optimisation problem111In principle, there is a possibility that we obtain a degenerate solution with ${\hat{f}}(X)=X\hat{\beta}_{\text{ora-sq}}$. However, we can obseve directly whether or not this occurs, and have never seen this happen in any of our numerical experiments..

From the reasoning above, we conclude that, under appropriate conditions, a simple test based on the asymptotic normality of $w_{\text{ora}}^{T}{R}_{\text{ora}}$ will keep the type I error under control. In order to create a version of the test statistic that does not require oracular knowledge of $D_{\beta_{0}}$, we may replace this quantity with variance estimates based either on $\hat{\beta}$ or on an estimate derived from the auxiliary data; we use the latter approach as this simplifies the analysis. The overall procedure is summarised in Algorithm 1 below.

Algorithm 1. Goodness-of-fit testing.

 

Input: sample $(X,Y)\in\mathbb{R}^{n\times p}\times\mathcal{Y}^{n}$; auxiliary sample $(X_{A},Y_{A})\in\mathbb{R}^{n\times p}\times\mathcal{Y}^{n};\lambda,\lambda_{A},\lambda_{\mathrm{sq}}>0.$

1. 1:

Estimation: Fit a GLM Lasso to $(X,Y)$ and $(X_{A},Y_{A})$ (with tuning parameters $\lambda,\lambda_{A}$ respectively) yielding estimators $\hat{\beta}$ and $\hat{\beta}_{A}$, respectively. 2. 2:

Residual prediction: Compute the residuals ${Y_{A}-{\mu}(X_{A}\hat{\beta}_{A})}$ and fit a flexible regression method of these residuals versus $X_{A}$ to obtain a prediction function ${\hat{f}}:\mathbb{R}^{n\times p}\rightarrow\mathbb{R}^{n}.$ 3. 3:

Near orthogonalization: Construct the diagonal weight matrix $\hat{D}_{A}^{2}:=\text{diag(}{\mu}^{\prime}(X\hat{\beta}_{A}))$ and compute an approximate projection of the prediction $\hat{D}_{A}{\hat{f}}(X)$ onto the column space of $\hat{D}_{A}X$:

[TABLE]

Define a direction

[TABLE] 4. 4:

Test statistic: Compute the residual vector $\hat{R}:=\hat{D}_{A}^{-1}(Y-{\mu}(X\hat{\beta}))$ and let $T:=\hat{w}_{A}^{T}\hat{R}$.

**Output: $p_{\text{value}}=1-\Phi(T)$

 **

In practice, the auxiliary dataset $(X_{A},Y_{A})$ would be obtained through sample splitting. The effect of the randomness induced by such a split can be mitigated using methods designed to aggregate over multiple sample splits, as studied for instance in Meinshausen et al., (2009).

2.2 Group testing

Our framework of residual prediction tests also encompasses significance testing of groups of regression coefficients in a generalized linear model. Suppose that we wish to test $H_{0}:\beta_{G}=0$ for a given group $G\subseteq\{1,\dots,p\}$. We first form the vector $\hat{R}_{G}$ of residuals based on a GLM Lasso fit of $Y$ on $X_{-G}$. Then, rather than constructing a single direction $w$ using an auxiliary dataset, we can use multiple directions given by the columns of $X_{G}$. Specifically, we use the test statistic $\max_{j\in G}|\hat{w}_{j}^{T}\hat{R}_{G}|$ where $\hat{w}_{j}$ is given by the scaled residuals of the weighted square-root Lasso regressions of $X_{j}$ on to $X_{-G}$.

Note that under the null, $X_{G}$ will be independent of the noise $\varepsilon$ (1), and so sample splitting is not necessary in this case to mitigate the potentially complicated dependence of the directions and residuals $\hat{R}_{G}$. The limiting distribution of the test however will not be Gaussian due to the maximisation over multiple directions. Instead, we argue that $\max_{j\in G}|\hat{w}_{j}^{T}\hat{R}_{G}|\approx\max_{j\in G}|\hat{w}_{j}^{T}\varepsilon|$ and then use a wild bootstrap procedure to approximate the distribution of this latter quantity. The overall procedure is summarised in Algorithm 2 below.

Algorithm 2. Group Test.

 

Input: Group $G\subseteq\{1,\dots,p\};$ sample $(X,Y)\in\mathbb{R}^{n\times p}\times\mathbb{R}^{n};$ $B\in\mathbb{N}$; $\lambda,\lambda_{\text{nw}}>0.$

1. 1:

Fit a GLM Lasso to $(X_{-G},Y)$ with a tuning parameter $\lambda$ to obtain an estimator $\hat{\beta}_{-G}\in\mathbb{R}^{p-|G|}$. Let $\hat{D}^{2}:=\text{diag}(\mu^{\prime}(X_{-G}\hat{\beta}_{-G})).$ Compute the vector of residuals $\hat{R}_{G}:=\hat{D}^{-1}(Y-{\mu}(X_{-G}\hat{\beta}_{-G})).$ 2. 2:

For each $j\in G,$ compute the nodewise regression estimator

[TABLE]

and let

[TABLE] 3. 3:

Evaluate the test statistic $T:=\max_{j\in G}|{\hat{w}_{j}^{T}\hat{R}_{G}}|.$ 4. 4:

For $b=1,\dots,B$ generate independent random variables $e_{1}^{b},\dots,e_{n}^{b}\sim\mathcal{N}(0,1)$ and let

[TABLE]

where $\hat{w}_{j,i}$ and $\hat{R}_{G,i}$ are the $i$-th entries of $\hat{w}_{j}$ and $\hat{R}_{G}$, respectively. 5. 5:

Calculate the p-value

[TABLE]

**Output: $p_{\text{value}}$

  **

Our Algorithm 2.2 is similar to the de-biased Lasso for generalized linear models (van de Geer et al.,, 2014). The main difference however is that the de-biased Lasso aims to ensure the directions $\hat{w}_{j}$ are almost orthogonal to $X_{-j}$, whereas we only impose near-orthogonality with respect to $X_{-G}$. Thus for large groups $G$, more of the direction of $X_{G}$ is preserved in the $\hat{w}_{j}$, which typically leads to better power of the test.

3 Theoretical guarantees

In this section we provide theoretical guarantees for the tests proposed in Algorithms 2.1 and 2.2. We consider an asymptotic regime with the sample size $n$ tending to infinity and the number of parameters $p=p_{n}$ growing as a function of $n$.

3.1 Size of the test

In the following sections, we show that under the null hypothesis, the size of the test is asymptotically correct. We explore goodness-of-fit testing in Section 3.1.1 and group testing in Section 3.1.2.

3.1.1 Goodness-of-fit testing

Here we show that our test statistic has a Gaussian limiting distribution and we establish a bound on the type I error of the test. In this section, we condition on the design and the auxiliary dataset. Our result makes use of the fact that when the model is well specified, the GLM Lasso performs well in terms of estimation. Specifically, under certain conditions, it holds with high probability that $\hat{\beta}\in\Theta(\lambda,\beta_{0},X)$, where $\Theta(\lambda,\beta_{0},X)$ is a local neighbourhood of $\beta_{0}$ defined by

[TABLE]

with $s:=\|\beta_{0}\|_{0}$ as the number of non-zero entries of $\beta_{0}$; see for example Bühlmann and van de Geer, (2011, Corollary 6.3). Sufficient conditions for this to occur include $\lambda\asymp\sqrt{\log p/n}$, $s=o(n/\log p)$ and further conditions on the tail behaviour of the errors $Y_{i}-\mu(x_{i}^{T}\beta_{0})$, the design matrix $X$ and the link function, as detailed below.

Condition 1**.**

Assume that $\mathbb{E}(Y_{i}|x_{i}=x)=\mu(x^{T}\beta_{0})$ and that $\emph{var}(Y_{i}|x_{i}=x)=\mu^{\prime}(x^{T}\beta_{0})$, that the inverse link function $u\mapsto\mu(u)$ is differentiable, $u\mapsto\mu^{\prime}(u)$ is Lipschitz with constant $L$, and that $\mu^{\prime}(u)>0$ for all $u\in\mathbb{R}$. Suppose moreover that the weights satisfy $\min_{i}D_{\beta_{0},ii}\geq d_{\min}$ for some constant $d_{\min}>0.$ Assume that \mathbb{E}\bigl{\{}|Y_{i}-\mu(x_{i}^{T}\beta_{0})|^{3}/D_{\beta_{0},ii}^{3}\bigm{|}X\bigr{\}}\leq C_{\varepsilon} for some constant $C_{\varepsilon}>0$, that $\max_{i=1,\dots,n}\|x_{i}\|_{\infty}\leq K_{X}$ for some $K_{X}\geq 1$ and that $12d_{\min}^{-2}LK_{X}s\lambda_{A}\leq 1.$

Condition 1 is satisfied for generalized linear models with canonical links under mild additional conditions. For example, in the case of logistic regression, the condition on the weights is satisfied if the class probability $\pi_{0}(x)=\mathbb{P}(Y_{i}=1|x_{i}=x)$ is bounded away from zero and one. Boundedness of the design (along with other conditions, including $12d_{\min}^{-2}LK_{X}s\lambda_{A}\leq 1$) guarantees that the weights can be consistently estimated. For our result below, it is convenient to introduce the shorthand notation ${Z_{A}}:=(X,X_{A},Y_{A})$.

Theorem 1**.**

Consider Algorithm 2.1 with tuning parameters $\lambda,\lambda_{A},\lambda_{\mathrm{sq}}>0$. Assume that Condition 1 is satisfied, that $\hat{\beta}_{A}\in\Theta(\lambda,\beta_{0},X_{A})$ and let

[TABLE]

Then there exists a constant222Here and below, the constants in the conclusions of our results may depend upon quantities introduced as constants in the relevant conditions for these results. $C>0$ such that whenever $Z_{A}$ satisfies $\hat{f}(X)\neq X\hat{\beta}_{\mathrm{sq}}$, we have for any $z\in\mathbb{R}$ that

[TABLE]

where $\Phi$ denotes the standard normal distribution function.

For the asymptotically optimal choice of tuning parameters $\lambda\asymp\lambda_{A}\asymp K_{X}\sqrt{\log p/n}$ and $\lambda_{\mathrm{sq}}\asymp\sqrt{\log p/n},$ the bound in Theorem 1 reduces to

[TABLE]

We now discuss the terms on the right-hand side of (7). The terms $\lambda_{\mathrm{sq}}\sqrt{n}s\lambda$ and $K_{X}s\lambda_{A}$ arise from bounding the bias term (near-orthogonalization step) and from bounding the weights, respectively. The presence of the $\|\hat{w}_{A}\|_{\infty}$ term in the bound stems from the contribution of each individual component $\hat{w}_{A,i}$ to the variance of the pivot term and that of the higher-order terms omitted in (2), which create a bias in the distribution of the test statistic.

To provide some intuition on the size of $\|\hat{w}_{A}\|_{\infty}$, recall that $\hat{w}_{A}$ is a vector in $\mathbb{R}^{n}$ with $\|\hat{w}_{A}\|_{2}=1$, so we may hope for $\|\hat{w}_{A}\|_{\infty}$ to be small; in fact, it can be shown in certain settings, and under additional technical conditions, that $\|\hat{w}_{A}\|_{\infty}=\mathcal{O}_{P}(\log n/\sqrt{n})$. We also remark that we observe $\hat{w}_{A}$, and if the size of its $\ell_{\infty}$-norm is a concern, then we can modify the square-root Lasso objective to control it explicitly. Indeed, consider setting

[TABLE]

and let

[TABLE]

Then the KKT conditions of the optimisation problem imply in particular both that a near-orthogonality condition similar to (4) is satisfied for suitable $\lambda_{\mathrm{sq}}$, and also that $\|\tilde{w}_{A}\|_{\infty}\leq\sqrt{n}\lambda_{\eta}$. Our empirical results in Section 4 however suggest that in practice $\|\hat{w}_{A}\|_{\infty}$ typically satisfies the necessary constraint and therefore we propose to use the simpler standard square-root Lasso without the above modifications.

3.1.2 Group testing

In this section, we derive theoretical properties for the group testing procedure proposed in Algorithm 2.2. Since we do not use sample splitting, we cannot directly apply the arguments of Theorem 1, as the direction $\hat{w}_{j}$ depends on $\hat{\beta}_{-G}$ via the weights $\hat{D}=D_{\hat{\beta}_{-G}}$. In order to understand this dependence, here we consider the setting of random bounded design and assume the response–covariate pairs $(Y_{i},x_{i})$ are all independent and identically distributed.

We aim to use the multiplier bootstrap procedure (Chernozhukov et al.,, 2013) to estimate the distribution of the test statistic $\max_{j\in G}|T_{j}|$ as described in Algorithm 2.2, but we first summarize a preliminary result which shows that, under appropriate conditions, $T_{j}$ can be asymptotically approximated by the zero-mean average $w_{j}^{T}\varepsilon$. Here we define $w_{j}:=D_{\beta_{0}}(X_{j}-X_{-G}\gamma_{0,j})/(\sqrt{n}\tau_{j})$, where

[TABLE]

and $\gamma_{0,j}$ is the population version of $\hat{\gamma}_{j}$ from Algorithm 2.2; i.e.

[TABLE]

Recall that

[TABLE]

In order to guarantee consistency of $\hat{\gamma}_{j}$ in Algorithm 2.2, we will introduce the additional requirement that $\gamma_{0,j}$ is sparse. Denote by $\beta_{0,-G}\in\mathbb{R}^{p-|G|}$ the subset of components of $\beta_{0}$ corresponding to indices in $G^{c}.$ We also define

[TABLE]

Condition 2**.**

• (i)

Let $\eta_{i}:=Y_{i}-\mu(x_{i}^{T}\beta_{0})$ and assume that there exist constants $c_{1},c_{2},c_{3}>0$ such that

[TABLE]

for all $i=1,\dots,n$.

• (ii)

There exists $K\geq 1$ such that $\|X\|_{\infty}\leq K$ and $\max_{j\in G}\|X_{-G}\gamma_{0,j}\|_{\infty}\leq K$.

• (iii)

*For some $\delta>0$ and all $\beta\in\mathbb{R}^{p}$ satisfying $\|\beta-\beta_{0}\|_{1}\leq\delta$ it holds that $c_{0}\leq\mu^{\prime}(x^{T}\beta)\leq C_{0}$ for some constants $c_{0},C_{0}>0$ and all $x\in\mathbb{R}^{p}$ with $\|x\|_{\infty}\leq K$. *

• (iv)

Denoting $\Sigma_{0}:=\mathbb{E}X^{T}D_{\beta_{0}}^{2}X/n$, we have $1/\Lambda_{\min}(\Sigma_{0})\leq C_{e}$ and $\|\Sigma_{0}\|_{\infty}\leq C_{e}$ for some constant $C_{e}>0.$

• (v)

We have $\max_{j\in G}\|\gamma_{0,j}\|_{0}\leq s$, $\|\beta_{0}\|_{0}\leq s$ and there exists a sequence $(a_{n})$ with $a_{n}\rightarrow 0$ and $K^{3}{s\log p}/{\sqrt{n}}\leq a_{n}$.

Proposition 1**.**

Assume that Conditions 1 and 2 are satisfied and assume that $\hat{\beta}_{-G}$ satisfies

[TABLE]

Consider $\hat{w}_{j},j\in G$ (assumed to be non-degenerate) as defined in Algorithm 2.2 with tuning parameters $\lambda\asymp\sqrt{\log p/n}$ and $\lambda_{\emph{nw}}\asymp K\sqrt{\log p/n}$. Assume that the null hypothesis $H_{0}:\beta_{0,G}=0$ holds. Then there exists a constant $C>0$ such that with probability at least $1-3/p$, we have

[TABLE]

Using Proposition 1 and the results of Chernozhukov et al., (2013), we can show that the quantiles of $\max_{j\in G}|T_{j}|$ can be approximated by the quantiles of $T^{b}:=\max_{j\in G}|\sum_{i=1}^{n}\hat{w}_{j,i}\hat{R}_{G,i}e_{i}^{b}|$ where $e_{1}^{b},\dots,e_{n}^{b}$ are independent $\mathcal{N}(0,1)$ random variables. We only need to guarantee that $T_{j}$ is well-approximated by $w_{j}^{T}\varepsilon$ and we pay a price of $\log|G|$ for testing $|G|$ hypothesis simultaneously, where $|G|$ denotes the cardinality of $G$.

Define the $\alpha$-quantile of $T^{b}$ conditional on $(x_{i},Y_{i})_{i=1}^{n}$ by

[TABLE]

where $\mathbb{P}_{e}$ is the probability measure induced by the multiplier variables $(e_{i}^{b})_{i=1}^{n}$ holding $(x_{i},Y_{i})_{i=1}^{n}$ fixed.

Theorem 2**.**

Assume the conditions of Proposition 1 and that there exist constants $C_{2},c_{2}>0$ such that

[TABLE]

Then there exist constants $c,C>0$ such that

[TABLE]

Theorem 2 shows that if the generalized linear model is correct and the null hypothesis $\beta_{0,G}=0$ holds, then Algorithm 2.2 produces an asymptotically valid p-value for testing this hypothesis.

3.2 Power analysis for goodness-of-fit testing

The choice of $w$ as postulated in Theorem 1 guarantees that the type I error for goodness-of-fit testing stays under control. We now provide guarantees on the local power of the test. To this end, let us suppose that the true model has conditional expectation function m_{0}(x)=\mu\bigl{(}x^{T}\beta_{0}+g_{0}(x)\bigr{)}. Here $g_{0}$ represents a small nonlinear perturbation of the linear predictor $x^{T}\beta_{0}$. Our aim here is to understand how this propogates through to the distribution of our test statistic. We will suppose that the perturbation is small enough that GLM Lasso estimates lie with high probability within local neighbourhoods of $\beta_{0}$. Let us first provide some intuition on the expected value of the test statistic under model misspecification. Writing $f_{0}(x)=x^{T}\beta_{0}+g_{0}(x)$, the expectation of the theoretical residuals \varepsilon=D^{-1}_{\beta_{0}}\bigl{(}Y-\mu(X\beta_{0})\bigr{)} is given by

[TABLE]

As argued in Section 2, the oracular test statistic $w^{T}R_{\text{ora}}$ can be approximated by the scalar product $w^{T}\varepsilon$. Therefore, in order to obtain good power properties, we should seek to construct a direction $w$ so as to maximize $\mathbb{E}w^{T}\varepsilon$. The oracular choice $w:=w_{0}/\|w_{0}\|_{2}$ yields by a Taylor expansion the approximation

[TABLE]

We therefore see that the test statistic behaves in expectation as a weighted $\ell_{2}$-norm of the nonlinear term $g_{0}(X)$. We now provide a theoretical justification which can be used for local asymptotic guarantees on the power of our method. We introduce the following conditions which are modifications of Condition 1 to account for the case when the model is misspecified.

Condition 3**.**

*Assume that $\mathbb{E}(Y_{i}|x_{i}=x)=\mu(f_{0}(x))$, the inverse link function $u\mapsto\mu(u)$ is differentiable, $u\mapsto\mu^{\prime}(u)$ is Lipschitz with constant $L$, and $\mu^{\prime}(u)>0$ for all $u\in\mathbb{R}$. Suppose moreover that the weights satisfy $\min_{i}D_{\beta_{0},ii}\geq d_{\min}$ for some constant $d_{\min}>0$. Assume that \mathbb{E}\bigl{(}|Y_{i}-\mu(f_{0}(x_{i}))|^{3}/D_{Y,ii}^{3}\bigm{|}X_{i}\bigr{)}\leq C_{\varepsilon} for a constant $C_{\varepsilon}>0$, where we denote $D_{Y}^{2}:=\emph{cov}(Y|X).$ Let $\max_{i=1,\dots,n}\|x_{i}\|_{\infty}\leq K_{X}$ for some $K_{X}\geq 1$ and assume that $12{d_{\min}^{-2}}{LK_{X}s\lambda}\leq 1$ and $|D_{Y,ii}^{2}D_{\beta_{0},ii}^{-2}-1|\leq 2d_{\min}^{-2}LK_{X}s\lambda$. *

Theorem 3**.**

Consider Algorithm 2.1 with tuning parameters $\lambda,\lambda_{A},\lambda_{\mathrm{sq}}$. Assume Condition 3, that $\hat{\beta}_{A}\in\Theta(\lambda,\beta_{0},X_{A})$ and let

[TABLE]

Then there exists a constant $C>0$ such that, whenever $Z_{A}$ is such that $\hat{f}(X)\neq X\hat{\beta}_{\mathrm{sq}}$, we have for any $z\in\mathbb{R}$ that

[TABLE]

where

[TABLE]

Under the null hypothesis, we have ${\bf\Delta}=0$ and $D_{Y}=D_{\beta_{0}}$; thus we recover the result of Theorem 1. The departure of the model from the null hypothesis is captured by ${\bf\Delta}$. Hence the theorem shows that to detect departures from the null, the direction $w$ must be “correlated” with the signal that remains in the residuals under misspecification, namely $\mu(f_{0}(X))-\mu(X\beta_{0})$. Under a local alternative, e.g. $f_{0}(X)=X\beta_{0}+g_{0}(X)/\sqrt{n}$, where $\|g_{0}(X)\|_{2}=1$, we have ${\bf\Delta}\asymp 1$.

Theorem 3 relies on rates of convergence of the “projected” estimator $\hat{\beta}$ in (10) when the model is misspecified. Oracle inequalities for Lasso-regularized estimators in high-dimensional settings have been well explored; we refer to Bühlmann and van de Geer, (2011) and the references therein. If there is misspecification, we hope that the projected estimator behaves as if it knows which variables are relevant for a linear approximation of the possibly nonlinear target $f_{0}.$ In Appendix A.4, we summarize how misspecification affects estimation of the best linear approximation, based on the approach of Bühlmann and van de Geer, (2011). These results guarantee that under a local alternative, the Lasso for generalized linear models still satisfies the condition

[TABLE]

3.3 Consequences for logistic regression

In this section we show how our general theory applies to the problem of goodness-of-fit testing for logistic regression models. We take $\mathcal{Y}=\{0,1\}$ and assume $(Y_{i}|x_{i}=x)\sim\text{Bernoulli}(\pi_{0}(x))$. Define

[TABLE]

that is, $\mathbb{E}(Y_{i}|x_{i}=x)=\pi_{0}(x)=\mu(f_{0}(x)),$ for the inverse link function $\mu(u)=1/(1+e^{-u}).$ The function $f_{0}$ may be potentially nonlinear in $x.$ The $\ell_{1}$-regularized logistic regression estimator is

[TABLE]

where $d(\xi):=\log(1+e^{\xi})$. Let $\beta_{0}\in\mathbb{R}^{p}$ be the best approximation obtained by a GLM (Bühlmann and van de Geer,, 2011, Section 6.3, p. 115). We define $S:=\{j:\beta_{0,j}\neq 0\}$ and $s:=|S|$.

Corollary 1 below follows by combining Theorem 3 with existing results on $\ell_{1}$-penalized logistic regression. We assume conditions, stated formally in Lemma 6 in Appendix A, which guarantee that with high probability this penalized estimator is sufficiently close to $\beta_{0}$.

Corollary 1**.**

Assume that the conditions of Lemma 6 in Appendix A hold and in addition assume that $\hat{\beta}_{A}\in\Theta(\lambda,\beta_{0},X_{A})$ and that $12c_{0}^{-2}Ks\lambda_{A}\leq 1$ where $c_{0}^{2}=(e^{\eta}/\epsilon_{0}+1)^{-2}$ and $K,\eta$ and $\epsilon_{0}$ are defined in Lemma 6. Suppose that $\lambda\asymp\lambda_{\rm{sq}}\asymp\lambda_{A}\asymp\sqrt{\log(2p)/n}.$ Then there exists a constant $C>0$ such that for any $z\in\mathbb{R}$, and whenever $Z_{A}$ is such that $\hat{f}(X)\neq X\hat{\beta}_{\mathrm{sq}}$,

[TABLE]

where

[TABLE]

4 Empirical results: Logistic regression

In this section we explore the empirical performance of the methods for goodness-of-fit testing and group testing in the setting of logistic regression. We begin by considering goodness-of-fit testing in low-dimensional settings in Section 4.1 and in high-dimensional settings in Section 4.2. Goodness-of-fit testing on semi-real data is investigated in Section 4.3, while in Section 4.4, we explore group testing in high-dimensional settings.

4.1 Low-dimensional settings

While for low-dimensional settings, there are numerous methods available for testing goodness-of-fit as discussed in Section 1.2, we show that our test from Algorithm 2.1 may be advantageous even here. We compare the performance of the our test (with residual prediction method being a random forest with default tuning parameter choices) against the Hosmer–Lemeshow $\hat{C}$ test, the Hosmer–Lemeshow $\hat{H}$ test (see Lemeshow and Hosmer Jr, (1982)) and the le Cessie–van Houwelingen–Copas–Hosmer unweighted sum of squares test (see Hosmer et al., (1997)). These tests are implemented in the function HLgof.test() in the R package MKmisc (Kohl,, 2018).

We simulated data from a logistic regression model with sample size $N=300$ and $p=10$ covariates according to

[TABLE]

where

[TABLE]

We considered different forms for the misspecification $g(\cdot)$:

• •

quadratic effect: (a) $g(u)=2u_{1}^{2}$, (b) $g(u)=2u_{5}^{2}$,

• •

interactions:

(c) $g(u)=u_{1}u_{2}$, (d) $g(u)=u_{1}u_{3}$, (f) $g(u)=u_{1}u_{4}$, (g) $g(u)=u_{4}u_{7}$.

Here $\sigma\geq 0$ measures the size of departure from the null hypothesis $H_{0}:\sigma=0$. Note that our GRP testing methodology requires an auxiliary sample of size $n_{A}$. We therefore randomly split the sample taking $n_{A}=n=N/2$, with $n$ being the number of observations in the main sample. The observation vectors $x_{i}$ follow a $\mathcal{N}_{10}(0,\Sigma_{0})$ distribution where

[TABLE]

is the Toeplitz matrix with correlation $\rho=0.6$. The results for the six settings above are shown in Figure 2. All methods maintain good control over type I error, but in most scenarios our GRP-test has significantly greater power compared with the other methods.

4.2 High-dimensional settings

Here we consider logistic regression models with two different types of misspecification from the presence of a pure quadratic effect and an interaction term. Specifically, we take the log-odds to be

[TABLE]

with

[TABLE]

and either (a) $g(u)=u_{1}u_{2}$ or (b) $g(u)=u_{1}^{2}/2$. The parameter $\sigma$ controls the degree of the misspecification and we look at $\sigma\in[0,3]$.

The observation vectors $x_{i}$ are independent with a Gaussian distribution $\mathcal{N}_{p}(0,\Sigma_{0}),$ where $\Sigma_{0}$ is the Toeplitz matrix (13) for a range of correlations $\rho\in\{0.4,0.6,0.8\}$. We consider a setting with $p=500$ and $N=800.$ The GRP-test requires sample splitting and we use split size $50\%$; thus the size of auxiliary sample is $n_{A}=400$. Again, we use a random forest as the residual prediction method.

In order to achieve better control of the type I error, we use a slight modification of the procedure proposed in Algorithm 1. Define $\hat{S}:=\{j:\hat{\beta}_{j}\neq 0\}$. Rather than calculating the direction $\hat{w}_{A}$ through $\hat{\beta}_{\mathrm{sq}}$, we instead use

[TABLE]

in its place within (6). In this way, we enforce that $\hat{w}_{A}$ is exactly orthogonal to $\hat{D}_{A}X_{\hat{S}}$. This helps to keep the remainder term arising from the asymptotic expansion of the test statistic under control, as can be seen from the following argument. Assume that a “beta-min condition” is satisfied, that is for all $j\in S:=\{k\in\{1,\dots,p\}:\beta_{0,k}\not=0\}$ it holds that $\beta_{0,j}\gtrsim\sqrt{s\log p/n}$, where $s:=|S|$. Then asymptotically, it holds that $\hat{S}\supseteq S$ with high probability (e.g. Bühlmann and van de Geer,, 2011, Corollary 7.6). On the event that this occurs, we have that the remainder term in (1) satisfies

[TABLE]

Even without such a beta-min condition, it is plausible that we will obtain a reduction in this bias term through this strategy of exact orthogonalization.

In the high-dimensional setting, there is no obvious method that we can use for comparison with our proposed GRP-test. Therefore as a theoretical benchmark, we consider an oracle GRP-test applied to a reduced design matrix containing only variables in the active set $S=\{1,\dots,5\}$, thereby reducing problem to a low-dimensional one. The results are reported in Tables 1 and 2, from which we see that the GRP-test does indeed control the type I error, and suffers only a relatively small loss in power compared with the oracle GRP-test.

4.3 Semi-real data example

We use a gene-expression dataset on lung cancer available from the NCBI database (Spira et al., (2007), https://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS2771) to illustrate the size and power performance of the goodness-of-fit test. We aim to detect if the model is a logistic regression, or if there are extra nonlinear effects. The full dataset contains airway epithelial gene expressions for 22215 genes from each of 192 smokers with (suspected) lung cancer, but this was reduced by taking the 500 genes with the largest variances. Having scaled the resulting variables, we fit a $\ell_{1}$-penalized logistic regression using cv.glmnet() from the package glmnet (Friedman et al.,, 2010) and obtained a parameter estimate $\hat{\beta}$ with its corresponding support set $\hat{S}$. We then fit a Gaussian copula model to the rows of the design matrix and generated a new, augmented design matrix $X$ by simulating a further $N=800$ observation vectors from this fitted model. Finally, we generated 800 new responses: Y_{i}|x_{i}=u\sim\text{Bern}\bigl{(}\hat{\pi}(u)\bigr{)} where

[TABLE]

$\bar{g}:=N^{-1}\sum_{i=1}^{N}g(x_{i})$ and $\hat{\sigma}_{g}^{2}:=(N-1)^{-1}\sum_{i=1}^{N}\{g(x_{i})-\bar{g}\}^{2}$, for the following three scenarios:

[TABLE]

where $j_{\ell}\in\hat{S}$, $\ell=1,\dots,4$ are uniformly sampled entries from $\hat{S}.$ We report rejection probabilities for all three scenarios from 100 repetitions in Table 3. In each case, the GRP-test is able to detect the misspecification relatively reliably, while keeping the type I error under control.

4.4 Group testing

Finally, we consider the problem of testing for the significance of groups of predictors using the methodology set out in Section 3.1.2. We compare the GRP-test (Algorithm 2.2) with the globaltest (Goeman et al.,, 2004) and the de-biased Lasso (van de Geer et al.,, 2014; Dezeure et al.,, 2015) for logistic regression. We consider logistic regression models with coefficient vector of the form

[TABLE]

for a range of values of $\theta\in[0,1.5]$, and look at testing the null hypothesis $\beta_{0,G}=0$ where $G=\{5,6,\ldots,p\}$. Thus larger values of $\theta$ correspond to more extreme violations of the null. Similarly to earlier examples, we use design matrices constructed via realisations of a random Gaussian design with Toeplitz covariance (13) and $\rho=0.6$. The results are reported in Figures 3 and 4. Both the GRP-test and the globaltest control the type I error at very close to the nominal level, while from Figure 3 the de-biased Lasso test is conservative; on the other hand, the GRP-test does very well in these examples in terms of power.

5 Discussion

In this work, we have introduced a new method for detecting conditional mean misspecification in generalized linear models based on predicting remaining signal in the residuals. For this task of prediction, we have a number of powerful machine learning methods at our disposal. Whilst these estimation performance of these methods is largely theoretically intractable, by employing sample-splitting and a careful debiasing strategy involving the square-root Lasso, our generalized residual prediction framework provides formal statistical tests with type I error control when used in conjunction with (essentially) arbitrary machine learning methods.

One requirement for these theoretical guarantees is that the sparsity of the true regression coefficient $\beta_{0}$ satisfies $s=o(\sqrt{n}/\log p)$, a condition that was also needed in related work on the de-biased Lasso (Zhang and Zhang,, 2014; van de Geer et al.,, 2014; Javanmard and Montanari,, 2014). It would be very interesting if this could be relaxed to $s=o(n/\log p)$ for instance, which would encompass settings where the GLM Lasso estimate satisfies $\|\hat{\beta}-\beta_{0}\|_{2}\to 0$ though $\|\hat{\beta}-\beta_{0}\|_{1}$ may be diverging.

Another interesting question is whether sample splitting can be completely avoided if we were able to obtain guarantees for an estimator $\hat{w}$ of a population direction $w$. Such alternatives to sample splitting could be particularly helpful for settings where there is dependence across the observations, such as in the case of generalized linear mixed effect models.

Acknowledgements: The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the ‘Statistical Scalability’ programme when work on this paper was undertaken, supported by EPSRC grant number EP/R014604/1. JJ is supported by a Swiss National Science Foundation fellowship. RDS is supported by an EPSRC First Grant and an EPSRC programme grant. PB is supported by the European Research Council under the grant agreement No. 786461 (CausalStats – ERC-2017-ADG). RJS is supported by an EPSRC fellowship and an EPSRC programme grant.

6 Proofs

6.1 Proofs for Section 3

6.1.1 Proofs for Section 3.1.1

Proof of Theorem 1.

This follows from the more general Theorem 3, noting that under the null hypothesis, ${\bf\Delta}=0$. The proof of Theorem 3 can be found in Section 6.1.3. ∎

6.1.2 Proofs for Section 3.1.2

Proof of Proposition 1.

For $j\in G$, let $u_{j}:=X_{j}-X_{-G}\gamma_{0,j}$ and define the sets

[TABLE]

[TABLE]

[TABLE]

where $C$ will be specified below. We first derive a high-probability bound for the set

[TABLE]

By Lemma 3 in Section A, there exists a constant $A>0$ (in the definition of $\mathcal{T}_{3,j}$) such that

[TABLE]

In the rest of the proof, we work on the event $\mathcal{T}.$ Define $D_{\beta_{0,-G}}:=\mu^{\prime}(X_{-G}\beta_{0,-G})$ (note that under $H_{0}$, it holds that $D_{\beta_{0,-G}}=D_{\beta_{0}}$). Now consider the decomposition

[TABLE]

where

[TABLE]

and where we write $D_{\hat{\beta}_{-G}}:=\hat{D}$. We first derive the rates of convergence for the estimator $\hat{\gamma}_{j}$ from Algorithm 2.2 and then proceed to bound the remainders. By Lemma 4 in Section A, there exist positive constants $C_{3}$ and $C_{4}$ such that on $\mathcal{T},$ we have

[TABLE]

where \hat{\tau}_{j}^{2}:=\bigl{\|}D_{\hat{\beta}_{-G}}(X_{j}-X_{-G}\hat{\gamma}_{j})\bigr{\|}_{2}^{2}/n.

Remainder $\text{rem}_{1,j}$: Let $k_{j}$ denote the index of column $X_{j}$ in the matrix $X_{j\cup G^{c}}$. For $j\in G$, we define

[TABLE]

and its population-level counterpart $\Gamma_{0,j}$ based on $\gamma_{0,j}$. First note that

[TABLE]

and similarly,

[TABLE]

Therefore, we obtain using Hölder’s inequality that

[TABLE]

On the set $\cap_{j\in G}\mathcal{T}_{3,j}\subseteq\mathcal{T},$ we have

[TABLE]

Next we bound $\|\hat{\Gamma}_{j}/\hat{\tau}_{j}-\Gamma_{0,j}/\tau_{j}\|_{1}.$ Firstly, we can decompose and bound

[TABLE]

Now note that

[TABLE]

Combining (21) with the fact that

[TABLE]

we obtain that

[TABLE]

where $e_{\ell}$ denotes the $\ell$-th standard basis vector in $\mathbb{R}^{p-|G|+1}$. Define $\Sigma_{\beta_{0},j\cup G^{c}}:=\mathbb{E}X_{j\cup G^{c}}^{T}D_{\beta_{0}}^{2}X_{j\cup G^{c}}/n$ and note from Condition 2(iv) that $\Sigma_{\beta_{0},j\cup G^{c}}$ is invertible. Consequently, $\Gamma_{0,j}/\tau_{j}^{2}=\Sigma_{\beta_{0},j\cup G^{c}}^{-1}e_{k_{j}}$, so $1/\tau_{j}^{2}=(\Sigma_{\beta_{0},j\cup G^{c}}^{-1})_{k_{j}k_{j}}.$ Further note that

[TABLE]

Therefore, by Condition 2(iv), we have

[TABLE]

Consequently, and by sparsity of $\gamma_{0,j}$ assumed in Condition 2(v), it follows that

[TABLE]

Moreover,

[TABLE]

By Condition 2(v), we can find $n_{0}\in\mathbb{N}$ such that $a_{n}\leq 1/(2C_{4}C_{e})$ for $n\geq n_{0}$. Then from (17) and (6.1.2), we obtain on $\mathcal{T}$ that for $n\geq n_{0}$,

[TABLE]

Using (20), we conclude that on $\mathcal{T}$,

[TABLE]

Consequently, combining (19) and (23), there exists a constant such that it holds on $\mathcal{T}$

[TABLE]

Remainder $\text{rem}_{2,j}$: By the mean value theorem, for each $i=1,\dots,n$, there exists $\alpha_{i}\in[0,1]$ such that

[TABLE]

where $\tilde{\beta}_{(i)}:=\alpha_{i}\hat{\beta}_{-G}+(1-\alpha_{i})\beta_{0,-G}.$ Consequently, using Hölder’s inequality and the KKT conditions from the optimization problem in Algorithm 2.2, we obtain

[TABLE]

where

[TABLE]

By Condition 2(ii), we have $\|X_{-G}\gamma_{0,j}\|_{\infty}\leq K$, so we obtain

[TABLE]

Since $\mu^{\prime}$ is Lipschitz and $\|x_{i}\|_{\infty}\leq K$, we obtain that on $\mathcal{T}$,

[TABLE]

Thus

[TABLE]

by Condition 2(iii) and Condition 2(v). Therefore, on $\mathcal{T}$, it follows that

[TABLE]

Consequently, on $\mathcal{T}$,

[TABLE]

Then using the result above, on $\mathcal{T}$, we obtain

[TABLE]

The result follows from (6.1.2), together with the bounds (14), (24), (25) and (6.1.2). ∎

Proof of Theorem 2.

We want to show that the quantiles of our test statistic for group testing,

[TABLE]

can be approximated by quantiles of its bootstrapped version

[TABLE]

where $(e_{i}^{b})_{i=1}^{n}$ is a sequence of independent and identically distributed $\mathcal{N}(0,1)$ random variables. We can apply Theorem 4 from Section A.3 together with Proposition 1. Adopting the notation of Theorem 4 we let

[TABLE]

and

[TABLE]

Note that the maxima above can be rewritten without the absolute values using the fact that for any $a\in\mathbb{R}$ it holds that $|a|=\max\{a,-a\}$. Thus for $i=1,\ldots,n$ and $j\in G$ we let $x_{ij}:=\sqrt{n}w_{j,i}\varepsilon_{i}$ and $\hat{x}_{ij}:=\sqrt{n}\hat{w}_{j,i}\hat{R}_{G,i}$. Moreover, for $i=1,\dots,n$ and $j\in G$ we also define $x_{i(j+|G|)}:=-x_{ij}$ and $\hat{x}_{i(j+|G|)}:=-\hat{x}_{ij}$. We will apply Theorem 4 with $x_{ij}$ and $\hat{x}_{ij}$ where $i=1,\dots,n$ and $j\in H:=G\cup\{j+|G|:j\in G\}.$

We now check that conditions (51), (52), (53), (54) and (55) needed for Theorem 4 are satisfied.

Checking condition (51):

First, by the tower property, we have

[TABLE]

By Condition 2(iii), it follows that $c_{0}\leq D_{\beta_{0},ii}^{2}\leq C_{0}$. Consequently, and using Condition 2(i), there exist constants $c,C$ such that $c\leq\mathbb{E}(\varepsilon_{i}^{2}|X)=\mathbb{E}(\eta_{i}^{2}/D_{\beta_{0},ii}^{2}|X)\leq C.$ It follows that

[TABLE]

Now recalling that

[TABLE]

we see that

[TABLE]

Therefore, combining (31) and (32), we obtain $c\leq\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}x_{ij}^{2}\leq C$ for $j\in H$, as required.

**Checking condition (52):

**Recall from (26) that we have the deterministic bound

[TABLE]

Using this bound, we will now check that for suitable $B_{n}\asymp K$,

[TABLE]

First observe that

[TABLE]

and

[TABLE]

Taking sufficiently large $B_{n}\asymp K$, we can therefore guarantee that (33) holds, as required.

**Checking condition (53):

**By Proposition 1, there exists a constant $C^{\prime}>0$ such that

[TABLE]

for $\zeta_{1}^{\prime}:=C^{\prime}K^{3}s{\log p}/\sqrt{n}$ and $\zeta_{2}^{\prime}:=3/p$. Next note that

[TABLE]

Now conditional on $(x_{ij})_{i=1}^{n}$ and $(\hat{x}_{ij})_{i=1}^{n}$ we have that $Z_{j}:=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(x_{ij}-\hat{x}_{ij})e_{i}^{b}\sim\mathcal{N}(0,\sigma_{j}^{2})$ where $\sigma_{j}^{2}:=\frac{1}{n}\sum_{i=1}^{n}(x_{ij}-\hat{x}_{ij})^{2}.$ Therefore,

[TABLE]

Then it follows by Borell’s inequality for any $t>0$,

[TABLE]

Taking $t:=\sqrt{2\log(2p+1)}\max_{j\in H}\sigma_{j}$ and noting that $|H|=2|G|\leq 2p,$ we obtain

[TABLE]

Denote $\Delta_{2}:=\max_{j\in H}\sigma_{j}^{2}.$ Then by Lemma 5 in Appendix A there exists a constant $C^{\prime\prime}>0$ such that

[TABLE]

Therefore, combining (34) and (35)

[TABLE]

So we can take

[TABLE]

and $\zeta_{2}:=4/p$ (in applying Theorem 4).

Checking conditions (54) and (55):

Finally, by assumption (9), there exist constants $C_{2}^{\prime},c_{2}>0$ such that

[TABLE]

and

[TABLE]

where $B_{n}\asymp K$. ∎

6.1.3 Proofs for Section 3.2

Proof of Theorem 3.

In this proof, it is convenient to write $\mathbb{P}_{A}(\cdot)$ and $\mathbb{E}_{A}(\cdot)$ as shorthand for $\mathbb{P}(\cdot|Z_{A})$ and $\mathbb{E}(\cdot|Z_{A})$ respectively. Consider the decomposition

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

There are three terms:

1. I.

The term $\phi$ is the pivot. By the Berry–Esseen theorem, we will show below that (after scaling) it is well approximated by a normal random variable. 2. II.

The term ${\bf\Delta}$ captures the deviation from the null hypothesis. If the null hypothesis is true, then ${\bf\Delta}=0$. 3. III.

The term $\text{rem}_{1}$ is a stochastic remainder term, for which we will develop a probabilistic bound below.

Let $\sigma:=\|D_{Y}\hat{D}_{A}^{-1}\hat{w}_{A}\|_{2}$. Then

[TABLE]

Now, for any $\epsilon>0$ and $z\in\mathbb{R}$,

[TABLE]

Similarly,

[TABLE]

Therefore, combining (36), (6.1.3) and (6.1.3), we find that

[TABLE]

where \mathrm{rem}_{0}:=\frac{1}{\sqrt{2\pi}}\bigl{(}\frac{1}{\sigma}-1\bigr{)}.

**Bound for the pivot.

**We apply the Berry–Esseen theorem for non-identically distributed summands to $Z_{0}:=\phi/\sigma$. Note that

[TABLE]

For $i=1,\dots,n$, denote U_{i}:=\hat{w}_{A,i}\hat{D}_{A,ii}^{-1}\bigl{\{}Y_{i}-\mu\bigl{(}f_{0}(x_{i})\bigr{)}\bigr{\}} and $\sigma_{i}^{2}:=\text{Var}(U_{i}|{Z_{A}})=\hat{w}_{A,i}^{2}\hat{D}_{A,ii}^{-2}D_{Y,ii}^{2}.$ Since $\mathbb{E}_{A}(U_{i})=0$, the Berry–Esseen theorem (Esseen,, 1942) yields that

[TABLE]

where $C_{1}>0$ is a universal constant. Hence using Condition 1,

[TABLE]

Bound for $\text{rem}_{0}$. To bound $|\sigma-1|$, we first bound $\|D_{Y}^{2}\hat{D}_{A}^{-2}-I\|_{\infty}$. First, since $\hat{\beta}_{A}\in\Theta(\lambda,\beta_{0},X_{A})$, $\mu^{\prime}$ is Lipschitz and $\|x_{i}\|_{\infty}\leq K_{X}$, we obtain

[TABLE]

Therefore, $\hat{D}_{A,ii}^{2}\geq D_{\beta_{0},ii}^{2}/2$ under the condition $6d_{\min}^{-2}LK_{X}s\lambda_{A}\leq 1/2.$ This then implies that

[TABLE]

Next, by assumption, we have

[TABLE]

(note that under $H_{0},$ it holds that $D_{Y}=D_{\beta_{0}}$, so $|\eta_{i,1}|=|D_{Y,ii}^{2}D_{\beta_{0},ii}^{-2}-1|=0$ and the required bound trivially holds). Then

[TABLE]

Finally, by assumption $6d_{\min}^{-2}LK_{X}s\lambda_{A}\leq 1/2$ and from the last display and Lemma 2, it follows that

[TABLE]

We also see from this calculation that under our conditions, $\sigma\geq 1/2$ and $\|D_{Y}\hat{D}_{A}^{-1}\|_{\infty}\leq 2$.

Bound for $\text{rem}_{1}$. A Taylor expansion of $\mu$ yields

[TABLE]

where $\tilde{\beta}_{(i)}=\alpha_{i}\beta_{0}+(1-\alpha_{i})\hat{\beta}$ for some $\alpha_{i}\in[0,1]$. Let $D_{\tilde{\beta}}$ denote a diagonal matrix with diagonal entries $D_{\tilde{\beta},{ii}}^{2}:=\mu^{\prime}(x_{i}^{T}\tilde{\beta}_{(i)}).$ Then

[TABLE]

where

[TABLE]

and

[TABLE]

Using Hölder’s inequality together with $\hat{\beta}\in\Theta(\lambda,\beta_{0},X)$ and the KKT conditions of the square-root Lasso (5), we have

[TABLE]

To bound the second term, $\text{rem}_{12}$, first by the Lipschitz property of $\mu^{\prime}$ we have

[TABLE]

Then, on the event that $\hat{\beta}\in\Theta(\lambda,\beta_{0},X)$,

[TABLE]

Therefore,

$\mathbb{P}_{A}(|\text{rem}_{1}|\geq r_{\text{rem}_{1}})\leq\delta$, where

[TABLE]

We conclude, using (39) and (40), and taking $\epsilon:=r_{\text{rem}_{1}}$ that

[TABLE]

The result follows. ∎

6.1.4 Proofs for Section 3.3

The logistic loss function is

[TABLE]

where d(\xi):=\log\bigl{(}1+e^{\xi}\bigr{)}, and we let $f_{\beta}(x):=x^{T}\beta$. We define the risk function

[TABLE]

and set $\beta_{0}:=\operatorname*{argmin}_{\beta\in\mathbb{R}^{p}}R(f_{\beta}|X)$.

Proof of Corollary 1.

We apply Theorem 3 to the case of logistic regression to obtain local guarantees on the power of the test. To this end, we need to bound $\delta$ in (10) and Condition 3 of Theorem 3.

To bound $\delta$, we note that by Lemma 6 in Section A with $t:=\log(2p)$, we have with probability at least $1-1/(2p)$ that

[TABLE]

In what follows we work on the event where this occurs. We next want to obtain a bound on $\|X(\hat{\beta}-\beta_{0})\|_{2}^{2}/n$. Note that the second derivative of the loss function is

[TABLE]

For $\|x\|_{\infty}\leq K$ and any $f$ with $\sup_{x:\|x\|_{\infty}\leq K}|f(x)-f_{0}(x)|\leq\eta$, we therefore have

[TABLE]

Note that for any $\tilde{\beta}$ on the line segment between $\beta_{0}$ and $\hat{\beta}$, we have

[TABLE]

Thus we can conclude using a Taylor expansion of the loss function that there exist $\{\tilde{\beta}_{(i)}:i=1,\ldots,n\}$, each on the line segment from $\beta_{0}$ to $\hat{\beta}$, such that

[TABLE]

We deduce that there exists a constant $\tilde{C}>0$ such that with $\lambda=\tilde{C}\sqrt{\log(2p)/n}$, we have \delta=\mathbb{P}\bigl{(}\hat{\beta}\notin\Theta(\lambda,\beta_{0},X)|X\bigr{)}\leq 1/(2p).

It remains to check that Condition 3 of Theorem 3 is satisfied. Firstly, the inverse link function $\mu(u)=1/(1+e^{-u})$ is differentiable and Lipschitz with constant $1$. Moreover, by (42), $D_{Y,ii}^{2}\geq d_{\min}^{2}:=c_{0}^{2}$ and also $D_{\beta_{0},ii}^{2}\geq c_{0}^{2}$. Finally, observe that \mathbb{E}\bigl{\{}|Y_{i}-\mu(f_{0}(x_{i}))|^{3}|X\bigr{\}}\leq 1. Moreover, $12d_{\min}^{-2}LK_{X}s\lambda=12c_{0}^{-2}LK_{X}s\lambda\leq 1$ by hypothesis. Therefore, Condition 3 is satisfied. ∎

Appendix A Appendix

A.1 Auxiliary lemmas

Lemma 1** (Hoeffding’s inequality for a maximum of $p$ averages).**

Suppose that for each $j=1,\ldots,p$, the random variables $Z_{1j},\ldots,Z_{nj}$ are independent with

[TABLE]

Then for all $t>0$

[TABLE]

Proof of Lemma 1.

Apply Corollary 17.1 in van de Geer, (2016) together with a union bound. ∎

Lemma 2**.**

Let $\tilde{A},\tilde{B}\in\mathbb{R}^{n\times n}$ be diagonal matrices and suppose that $\tilde{B}$ is invertible. Let $w\in\mathbb{R}^{n}$ satisfy $\|w\|_{2}=1$. Then

[TABLE]

Proof of Lemma 2.

[TABLE]

Hence

[TABLE]

as required. ∎

A.2 Auxiliary lemmas for Group Testing

Lemma 3**.**

Under the conditions of Proposition 1, we have

[TABLE]

Proof.

To obtain a probability bound for $\mathcal{T}_{1,j}$, we can apply Lemma 1, noting that

[TABLE]

with $Z_{ijk}:=u_{j,i}D_{\beta_{0},ii}X_{-G,i,k}$, where $X_{-G,i,k}$ is the $(i,k)$-th entry of the matrix $X_{-G}$ and $u_{j,i}$ is the $i$-th entry of $u_{j}.$ Note that by Condition 2 (ii), it follows that $|u_{j,i}|\leq 2K$ and by Condition 2 (iii), we have $|D_{\beta_{0},ii}|\leq C_{0}$. Therefore, $|Z_{ijk}|\leq c_{i}$ for $c_{i}:=2C_{0}K^{2}$ and for all $i,j,k$. Thus Lemma 1 implies that for all $t>0,$

[TABLE]

Therefore,

[TABLE]

For the set $\mathcal{T}_{3,j}$, by the sub-Gaussianity of $\eta_{i}=Y_{i}-\mu(x_{i}^{T}\beta_{0})$ from Condition 2 (i), there exists a constant $C>0$ such that

[TABLE]

Therefore,

[TABLE]

Using bounds (43) and (44), the fact that $|G|\leq p$ and the assumption $\mathbb{P}(\mathcal{T}_{2}^{c})\leq 1/p$, we obtain by a union bound that

[TABLE]

∎

Lemma 4**.**

Assume Conditions 1 and 2. For $j\in G$, we let $u_{j}:=X_{j}-X_{-j}\gamma_{0,j}$ and define the sets

[TABLE]

[TABLE]

Then there exist $\lambda_{\mathrm{nw}}\asymp\sqrt{\log p/n}$, and positive constants $C_{3}$ and $C_{4}$ such that on the set $\cap_{j\in G}(\mathcal{T}_{1,j}\cap\mathcal{T}_{2}),$ it holds that

[TABLE]

whenever $\hat{\tau}_{j}^{2}:=\|\hat{D}(X_{j}-X_{-G}\hat{\gamma}_{j})\|_{2}^{2}/n>0$.

Proof of Lemma 4.

To obtain rates of convergence for $\hat{\gamma}_{j}$ from Algorithm 2.2, we follow the arguments in the proof of Theorem 3.2 in van de Geer et al., (2014), which considers nodewise regression with random bounded design. The difference is that they define a nodewise regression program with design matrix $X_{-j}$ and use the Lasso, whereas we want to apply the nodewise regression with a smaller design matrix, $X_{-G}$, and we use the square-root Lasso. We also seek finite-sample, as opposed to asymptotic, bounds, but this requires only minor modifications. But since we assume that $\hat{\tau}_{j}>0$, the square-root Lasso program with penalty $\lambda_{\text{nw}}$ corresponds to the Lasso program with penalty $\lambda_{\text{Lasso}}:=\hat{\tau}_{j}\lambda_{\text{nw}}.$ We now check that the appropriate finite-sample analogues of conditions (D1)–(D5) of Theorem 3.2 from van de Geer et al., (2014) are satisfied for $\tilde{X}:=X_{j\cup{G^{c}}}$. Firstly, the analogues of (D1), (D2), (D4) are satisfied directly by the assumptions in Conditions 1 and 2. For (D3), first note that the smallest eigenvalue of $\Sigma_{\beta_{0},j\cup G^{c}}:=\mathbb{E}\tilde{X}^{T}D_{\beta_{0}}^{2}\tilde{X}/n$ is lower bounded by the smallest eigenvalue of $\Sigma_{0}=\mathbb{E}X^{T}D_{\beta_{0}}^{2}X/n$, which is in turn lower bounded by $1/C_{e}$. Similarly, $\|\Sigma_{\beta_{0},j\cup G^{c}}\|_{\infty}\leq\|\Sigma_{\beta_{0}}\|_{\infty}\leq C_{e},$. Finally, Condition (D5) is satisfied on $\mathcal{T}_{2}$.

As in the proof of Proposition 1, let $k_{j}$ denote the index of column $X_{j}$ in the matrix $X_{j\cup G^{c}}$. We write

[TABLE]

and $\Gamma_{0,j}$ for its analogy defined in terms of $\gamma_{0,j}.$ Then note that we can write $X_{j}-X_{-G}\hat{\gamma}_{j}=\tilde{X}\hat{\Gamma}_{j}$ and recall that $\hat{\tau}_{j}^{2}=\|\hat{D}\tilde{X}\hat{\Gamma}_{j}\|_{2}^{2}/{n}$ and $\tau_{j}^{2}=\mathbb{E}\|D_{\beta_{0}}\tilde{X}\Gamma_{0,j}\|_{2}^{2}/{n}$. By inspecting the proof of Theorem 3.2 of van de Geer et al., (2014), we conclude that there exist positive constants $C_{3}$ and $C_{4}$ such that on $\cap_{j\in G}(\mathcal{T}_{1,j}\cap\mathcal{T}_{2}),$ it holds that

[TABLE]

as required. ∎

The following lemma bounds a term $\Delta_{2}$ defined in the proof of Theorem 2.

Lemma 5**.**

Under the conditions of Theorem 2, there exists a constant $C^{\prime\prime}>0$ such that

[TABLE]

Proof of Lemma 5.

On the set $\mathcal{T}$ defined in the proof of Proposition 1, we have

[TABLE]

Now we bound $r_{1}.$ By similar arguments as in the proof of Proposition 1, we will now show that on $\mathcal{T},$

[TABLE]

First,

[TABLE]

Now we can use Condition 2(ii), (18), (27), (28) and (23) to bound the terms in the last display and obtain

[TABLE]

Moreover, by the sub-Gaussianity of $\eta_{i}$ and since $D^{2}_{\beta_{0},ii}\geq c_{0}$ (by Condition (2)(iii)), there exist constants $C,C^{\prime\prime\prime}>0$ such that with probability at least $1-1/p$, we have

[TABLE]

Therefore, with probability at least $1-1/p-\mathbb{P}(\mathcal{T}^{c})$,

[TABLE]

We now bound $r_{2}$. Using Condition 2(iii), together with the fact that $\hat{\beta}_{-G}\in\Theta_{-G}(\lambda,\beta_{0})$ on $\mathcal{T}$ we have that on this event,

[TABLE]

Then using the fact that $\hat{\beta}_{-G}\in\Theta_{-G}(\lambda,\beta_{0})$ on $\mathcal{T}$ and using that

[TABLE]

(which follows similarly as in the proof of Theorem 3) and using (47), we obtain that on $\mathcal{T}$,

[TABLE]

Then

[TABLE]

where we used $\|\hat{w}_{j}\|_{\infty}^{2}=\mathcal{O}(K^{2}/n)$ (which follows from (29)).

Overall, collecting (49) and (48) there exists a constant $C^{\prime\prime}$ such that

[TABLE]

∎

A.3 Multiplier bootstrap

We summarize Corollary 3.1 from Chernozhukov et al., (2013). To this end, we need the following condition.

Condition 4**.**

Let $(x_{i})_{i=1}^{n}$ be $n$ independent random vectors with values in $\mathbb{R}^{g}$ satisfying

[TABLE]

and

[TABLE]

*Define *

[TABLE]

Let $(e_{i})_{i=1}^{n}$ be a sequence of i.i.d. $\mathcal{N}(0,1)$ random variables independent of $(x_{i})_{i=1}^{n}$ and define

[TABLE]

Assume that there exist $\zeta_{1},\zeta_{2}\geq 0$ such that

[TABLE]

where $\mathbb{P}_{e}$ is the probability measure induced by the multiplier variables $(e_{i})_{i=1}^{n}$ holding $(x_{i})_{i=1}^{n}$ fixed.

Theorem 4** (Corollary 3.1 in Chernozhukov et al., (2013)).**

Suppose that Condition 4 is satisfied with

[TABLE]

and

[TABLE]

Then there exist constants $c,C>0$ depending only on $c_{1},C_{1},c_{2}$ and $C_{2}$ such that

[TABLE]

where $c_{W}(\alpha)$ is the $\alpha$-quantile of $W$ conditional on $(x_{i})_{i=1}^{n}$ given by

[TABLE]

A.4 Oracle inequalities for logistic regression under misspecification

We require a condition on the design matrix known as the compatibility condition (Bühlmann and van de Geer,, 2011).

Definition 1** (Compatibility constant).**

We say that the compatibility condition is met with constant $\phi>0$ if for all $\beta\in\mathbb{R}^{p}$ that satisfy $\|\beta_{S^{c}}\|_{1}\leq 3\|\beta_{S}\|_{1},$ it holds that

[TABLE]

For our final lemma, we use the notation of Sections 3.3 and 6.1.4.

Lemma 6**.**

Suppose that there exists a constant $K>0,$ such that

[TABLE]

and that $X$ satisfies the compatibility condition with constant $\phi>0.$ Take $t>0$ and let

[TABLE]

[TABLE]

Assume that there exist constants $\epsilon_{0},\eta$ such that

[TABLE]

[TABLE]

For some constant $M\geq 8$ take $\lambda$ satisfying $8\lambda_{0}\leq\lambda\leq M\lambda_{0}$ and $17\lambda s(e^{\eta}/\epsilon_{0}+1)^{2}/\phi^{2}\leq\eta/(2K)$, and further assume that

[TABLE]

Then with probability at least $1-e^{-t},$ it holds that

[TABLE]

Proof of Lemma 6.

The proof follows from Lemma 6.8 and Section 6.7 in Bühlmann and van de Geer, (2011). ∎