Simulation-calibration testing for inference in Lasso regressions
Matthieu Pluntz Cyril Dalmasso Pascale Tubert-Bitter Ismail Ahmed



The Wanderer: Charting WASP-77A b's Formation and Migration Using a System-Wide Inventory of Carbon and Oxygen Abundances
    Maleah Rhem
    September 9, 2024
=========================================================================================================================





We propose a test of the significance of a variable appearing on the Lasso path and use it in a procedure for selecting one of the models of the Lasso path, controlling the Family-Wise Error Rate.


Our null hypothesis depends on a set A of already selected variables and states that it contains all the active variables. We focus on the regularization parameter value from which a first variable outside A is selected. As the test statistic, we use this quantity’s conditional p-value, which we define conditional on the non-penalized estimated coefficients of the model restricted to A. We estimate this by simulating outcome vectors and then calibrating them on the observed outcome’s estimated coefficients. We adapt the calibration heuristically to the case of generalized linear models in which it turns into an iterative stochastic procedure. We prove that the test controls the risk of selecting a false positive in linear models, both under the null hypothesis and, under a correlation condition, when A does not contain all active variables.


We assess the performance of our procedure through extensive simulation studies. We also illustrate it in the detection of exposures associated with drug-induced liver injuries in the French pharmacovigilance database.


Keywords: Variable selection, high-dimensional regression, Lasso, empirical p-value, conditional p-value, FWER control.




§ INTRODUCTION

Variable selection in high-dimensional regressions is a classic problem in health data analysis. It aims to identify a limited number of factors associated with a given health event among a large number of candidate variables such as genetic factors or environmental or drug exposures. The Lasso regression <cit.> provides a series of sparse models which variables enter one after another as the regularization parameter λ decreases. To be used in variable selection, the Lasso requires a further procedure for choosing the value of λ and thus the associated model. By itself, the Lasso does not indicate whether the variables it selects are statistically significant.


<cit.> address this by adapting the paradigm of hypothesis testing to the Lasso. They develop a significance test for j_k, the k-th variable selected on the Lasso path, by measuring the impact of its inclusion in the model compared to the selection by the Lasso of the set A of the first k-1 variables only. Their test statistic is a difference of covariances between the Lasso results on two datasets — complete or restricted to the variables in A — at the regularization parameter λ_k+1 where the next variable enters the Lasso path. They determine the distribution of this statistic in the linear model, under the null hypothesis asserting that the first k-1 variables selected by the Lasso contain all the active variables. It is therefore possible to calculate the exact p-value associated with this test. However, in the binary model, the test statistic's exact distribution is not known.


We propose a testing procedure that addresses the same issue as the test by Lockhart et al.: produce, for each variable selected by the Lasso, a p-value that measures its significance while correcting from the drawbacks of this test: invalidity in nonlinear cases, and sometimes lack of power. We retain the idea of exploiting the value of the regularization parameter at which a given variable enters the Lasso path. All else being equal, the earlier on the Lasso path a variable is selected — that is, at a high value of λ — the stronger and more significant its association with the response. Like Lockhart et al., we propose to test in the framework of hypothesis testing, for each of the variables entering the Lasso path, whether its selection λ is higher than it should be if there were no link between this variable and the response variable.


To this end, we are inspired by the permutation selection of <cit.>. Their method aims to determine from the data an optimal value of λ, and then select the model given by the Lasso at this value. The idea is that the optimal λ is the one below which it is likely that the Lasso inadvertently selects variables that are in reality independent of Y. To find it, they generate a population of random permutations y^(l), l = 1, .., N of the observed response y. This imposes "independence"  between the permuted responses and each of the X_j (in the sense that all the y^(l)_i have the same distribution despite the variation in X_i,j). Then they perform the Lasso regression of each permuted response vector on X, and measure λ_0 (y^(l)), the smallest value of λ at which this Lasso regression selects no variable. The chosen regularization parameter is then the median of {λ_0 (y^(1)), .., λ_0 (y^(N)) }.


We adopt the idea of generating responses which follow the same distribution as the observed response, but which has been made artificially independent from the covariates. As in permutation selection, we perform the Lasso regression of the simulated response on the observed covariates and we focus on λ_0, the highest value of λ where a covariate, which by construction is in fact independent of the response, is selected. These λ_0 obtained from several simulated datasets make up a reference population which is representative of the case where the response and covariates are independent. Unlike permutation selection, where one applies the median λ_0 from this population to the data of interest, we compare the λ_0 obtained on the data of interest to this reference population in order to test the significance of this λ_0, estimating the p-value of the test by Monte Carlo from the reference population. This generalizes to quantities other than λ_0: the λ at which any variable enters the Lasso path, not just the first one. While Sabourin et al. retain only the distribution of the simulated response vectors while breaking their association with all the X_j through permutation, we instead retain the simulated responses' association with some of the covariates by applying a “calibration” to them. Therefore our algorithm is called simulation-calibration.


From a variable selection perspective, if the test concludes that the λ at which a variable enters the Lasso path is significant, we select the variable. By iterating the test, we obtain a complete model selection procedure.


Section <ref> presents our null hypothesis and basic notations. Section <ref> shows the need for, and defines, the conditional p-value which we which we will use as a test statistic in the linear case. Section <ref> presents the simulation-calibration algorithm estimating this statistic, and section <ref> shows that it is valid and permits control of the type 1 error. Section <ref> adapts this test and this algorithm to generalized linear models. Section <ref> presents iterative model selection procedures based on simulation-calibration. Section <ref> states a theorem which permits the model selection procedures to control the FWER (or FDR) under certain conditions in the linear case. Section <ref> shows extended simulation studies of the validity of the simulation-calibration test and of the iterative procedures' variable selection performances. Section <ref> shows an application to real-word pharmacovigilance data. 



§ PROBLEM STATEMENT AND NOTATIONS


Let X ∈ℝ^n × p be a matrix of p covariates and n observations and y ∈ℝ^n a response vector of size n. We consider the linear regression model:


    y = β_0 + Xβ + ϵ,     ϵ∼𝒩(0,σ^2 I_n )



where β_0 ∈ℝ, β∈ℝ^p, and σ > 0. Let β̂^Lasso (λ,y) be the Lasso estimator of β, which depends on the regularization parameter λ and the response y. It also depends on the covariate matrix X, but this is considered constant in the following.


Let A be a subset of 1, .., p. We assume that the covariates X_j, j ∈ A, are active, and we seek to determine if other covariates are active. We thus define the null and alternative hypotheses, which depend on A:

    H_0 (A) : ∀ j ∉ A,β_j = 0 
    
    H_1 (A) : ∃ j ∉ A,β_j ≠ 0



Although A is a priori an arbitrary subset of {1, .., p}, the idea, as in Lockhart et al., is to take A as the set of indices of the first k-1 covariates entering the Lasso path, for relatively small values of k (even if p can be large). The simplest particular case is A = ∅.


We consider the following statistic:

    λ_A(y) = sup{λ≥ 0 :  ∃ j ∉ A,  β̂_̂ĵ^Lasso (λ,y) ≠ 0 }



That is, the λ at which the first variable outside of A is selected, or the k-th variable selected if A includes the first k-1 variables. We will reject the null hypothesis if we observe a value of λ_A that is too high compared to what is expected under the null hypothesis. An abnormally high λ at which a variable enters the Lasso typically results from an association between the response and this variable. Therefore, we interpret a test of H_0 (A) based on λ_A as assessing the significance of the first variable selected by the Lasso outside of A, whose index we call j_A. 


However, H_0 (A) might be false while the variable j_A is inactive, if there exists an active variable not belonging to A that is selected "later" on the Lasso path, at a λ lower than λ_A. In this case, it is correct to reject H_0 (A) from a hypothesis testing perspective, but incorrect to conclude that the variable j_A is significant. In section <ref>, we show under which conditions the probability of this event is controlled.


Furthermore, it is possible that multiple variables enter the Lasso path simultaneously at the parameter λ_A, whether mathematically (for example, in the case of both X and y being binary, there is a non-zero probability of exact symmetry between the associations of y with two distinct covariates) or approximately if the difference between the λ at which two distinct variables are selected is too small to measure. In both cases, rejecting H_0 (A) leads to the selection of all these variables.



§ CONDITIONAL P-VALUE
 

If λ_A were used as the test statistic, its associated p-value would be, by definition:

    p_A^0 (y) = P_H_0 (A)( λ_A (Y) ≥λ_A (y) )


where Y is a random variable that follows the model <ref>.

In general, there is no simple explicit analytical expression for the result of a Lasso regression. This prevents determining the exact distribution of λ_A (y) and thus calculating p_A^0 (y) or other similar quantities. Therefore, we aim to estimate this probability using the Monte Carlo method. To do this, we need to simulate i.i.d. outcome vectors which follow model <ref> and satisfy H_0(A).

However, this not possible without additional information because even assuming H_0(A) is true, the parameters of model <ref> are not known: H_0(A) guarantees that β_j = 0 for j ∉ A but says nothing about β_j for j ∈ A or σ. A naive attempt to perform this simulation, where β_0, the β_j for j ∈ A, and σ are simply replaced by their estimates on the data, leads to a biased estimation of the p-value (see section <ref> and the supplementary material).

Therefore, instead of p_A^0 (y), we focus on a variant whose definition includes all the necessary information for its estimation by Monte Carlo. It requires conditioning on the parameters that remain unknown under H_0(A). Consider the following linear model:


    y = X_A β_A + ϵ_A,     ϵ_A ∼𝒩(0,σ_A^2 I_n)


where X_A is the matrix composed of a column vector of 1s and the columns of X whose indices belong to A, and β_A is the vector (β_j)_j ∈{0}∪ A. 

This is the reduced form that the model <ref> takes when H_0 (A) is true. Then we have σ = σ_A. Let θ_A = (β_A,σ_A) be the parameter vector of the model <ref>, Θ_A = ℝ^1+|A|×ℝ_+ the space of its possible values, and θ_A (y) the maximum likelihood estimator (unpenalized) of θ_A. We define the following test statistic:
 

    p_A (y) = P_H_0 (A)(λ_A (Y) ≥λ_A (y) | θ_A (Y) = θ_A (y) ).



This formula is similar to that of a p-value but with an added conditioning, making p_A (y) a conditional p-value.


p_A (y) is defined so that we can know its distribution as precisely as possible. We have the following result:

 
Under H_0 (A),

    ∀ t ∈ [0,1], P_H_0 (A)(p_A(Y) ≤ t ) ≤ t.



In other words, the null distribution of p_A(Y) stochastically dominates the uniform distribution on [0,1]. This phrase, used in decision theory, means that the cumulative distribution function of the dominant probability distribution is lower at every point than that of the dominated distribution, and thus the dominant distribution is systematically shifted towards higher values compared to the dominated distribution.


In practice, the distribution of λ_A (Y) conditional on θ_A (Y) = θ_A' is typically continuous for all parameter vectors with a non-zero standard deviation, i.e., any θ_A' ∈Θ_A^* where Θ_A^* = ℝ^|A|×ℝ_+^*. We call the assumption that this is true the continuity assumption. It implies the following simpler result:

 
Under the continuity assumption, p_A(Y) follows a uniform distribution on [0,1].

The proofs of Lemma <ref> and corollary <ref> are in the supplementary material.

Lemma <ref> allows us to use p_A (y) as a test statistic. This quantity is close to 0 when λ_A (y) is high, which is indicative of H_1 (A). Therefore, we wish to reject H_0 (A) if p_A (y) is sufficiently small. The lemma ensures that for all α∈ [0,1], rejecting H_0 (A) if and only if p_A (y) ≤α guarantees that the type 1 error is at most α. The p-value associated with p_A (y) is p_A(y) itself under the continuity assumption, and not larger than p_A(y) in the general case.




§ ALGORITHM ESTIMATING THE CONDITIONAL P-VALUE
 
Unlike p_A^0 (y), we are able to estimate p_A (y) with a Monte Carlo method. The problem with estimating p_A^0 (y) is that it is not possible to simulate Y following its distribution under H_0 (A), which is not known. To estimate p_A (y), due to the conditioning in its definition, we must instead simulate Y following its distribution under H_0 (A) conditional on θ_A (Y). The following Monte Carlo algorithm does this. It is parameterized by an integer N, the number of simulations, which controls estimation accuracy.



In practice, at step 2.1, we use θ_A^sim = θ_A (y), but other choices are possible because all values of θ_A^sim produce calibrated responses y^(l) that follow the same distribution.



§ PROPERTIES OF THE ALGORITHM
 


 §.§ Proof of consistency and properties for a given response vector



  §.§.§ Conditional distribution of the response vector

The algorithm's consistency is based on knowing the conditional distribution of Y under the null hypothesis. It is given by the following lemma:

 
Under H_0 (A), the distribution of Y conditional on θ_A (Y) = θ_A (y) is the uniform distribution on the set S ( θ_A (y) ) = { y' ∈ℝ^n |θ_A (y') =θ_A (y) }.

The proof is in the supplementary material.


The definition of p_A (y) (<ref>) is based on an arbitrary random vector Y which follows this conditional distribution. Therefore, Y can be replaced by any random variable following the same distribution, given by lemma <ref>:
 
Any random variable Y^u that follows a uniform distribution over S ( θ_A (y) ) satisfies P(λ_A (Y^u) ≥λ_A (y)) = p_A (y).




  §.§.§ Calibration and simulation under the conditional distribution

Corollary <ref> implies that to estimate p_A (y) by Monte Carlo, one just has to simulate random response vectors uniformly over S ( θ_A (y) ). Now, we show that step 2. of algorithm <ref> does this.

We use calibration functions, which we define as follows. For any parameter vectors θ_A^(1) = (β_A^(1), σ_A^(1)) ∈Θ_A^* and θ_A^(2) = (β_A^(2), σ_A^(2)) ∈Θ_A:

    cal_θ_A^(1)→θ_A^(2) : ℝ^n    →   ℝ^n 
    
         y'    →    X_A β_A^(2) + σ_A^(2)/σ_A^(1)(y' - X_A β_A^(1)).


They map one given value of θ_A to another, that is:
 
For all θ_A^(1)∈Θ_A^*, θ_A^(2)∈Θ_A, for all y^(1)∈ℝ^n satisfying θ_A (y^(1)) = θ_A^(1), we have:

    θ_A( cal_θ_A^(1)→θ_A^(2) (y^(1)) ) = θ_A^(2).


The proof is in the supplementary material.

Calibration functions are a way to "impose" on a response vector y' the condition that θ_A (y') = θ_A^(2) for any target parameter vector θ_A^(2). As a composition of a homothety and a translation, they can be seen informally as the simplest change one can make to y' while imposing this condition. In particular, the empirical correlation structure between y' and variables not belonging to A varies little between y' and its calibrated version: it is only modified to the extent that it is carried by the variables in A.


Calibration enables simulating the uniform distribution over S ( θ_A (y) ):

 
Step 2. of algorithm <ref> simulates each y^(l) following the uniform distribution over S ( θ_A (y) ).



By replacing Y and θ_A (y) in lemma <ref> with y^sim and θ_A (y^sim), we find that the distribution of y^sim conditional on θ_A (y^sim) is uniform over S (θ_A (y^sim)). Furthermore, conditional on θ_A (y^sim),

    cal_θ_A (y^sim) →θ_A (y): y' → X_A β_A (y) + σ_A (y)/σ_A (y^sim)(y' - X_A β_A (y^sim) )

is the composition of a homothety and a translation. These transformations preserve the uniformity of a probability distribution. Also, according to lemma <ref>, it maps S ( θ_A (y^sim) ) onto S ( θ_A (y) ). Therefore, conditional on θ_A (y^sim), y^(l) follows the uniform distribution over S ( θ_A (y) ). By integration, without conditioning y^(l) follows the uniform distribution over S ( θ_A (y) ).




  §.§.§ Distribution of the empirical conditional p-value

The preceding lemmas give us the distribution of p_A (y) produced by the algorithm. It is the binomial distribution:
 

    N p_A (y) ∼Bin (N, p_A (y)).



For each l, y^(l) follows the uniform distribution over S ( θ_A (y) ). Therefore, according to corollary <ref>, for each simulation:


    P ( λ_A (y^(l) ) ≥λ_A (y) ) = p_A (y).


N p_A (y) is thus the sum of N independent binary random variables with mean p_A (y).


p_A (y) is therefore an unbiased and, due to the law of large numbers, consistent estimator of p_A (y). Its variance is Var( p_A (y) ) = p_A (y)(1-p_A (y))/N, which converges to 0 as N goes to infinity.



 §.§ Properties for a random response vector
 

The above applies to a given response vector y. In the case of a random response vector Y that follows model <ref> (or model <ref> under the null hypothesis), it is interesting to know the distribution of p_A (Y). That is because we do not observe the theoretical test statistic p_A (Y): we only estimate it. The decision to select or not a variable is based on p_A (Y), and any control of the type I error must be inferred from the distribution of p_A (Y) under the null hypothesis. We determine this distribution exactly under the continuity assumption defined in <ref>, and then we infer a stochastic dominance of this distribution in the general case which guarantees the type I error control.


The relationship between p_A (y) and p_A (y) (lemma <ref>) allows us to translate results about the distribution of p_A (Y) to that of p_A (Y), moving from a continuous to a discrete distribution. We have two results with and without the continuity assumption:

 
Under H_0 (A) and under the continuity assumption, the distribution of p_A (Y) is the discrete uniform distribution over {0, 1/N, .., 1-1/N, 1}.


Under H_0 (A), the distribution of p_A (y) stochastically dominates the discrete uniform distribution over {0, 1/N, .., 1-1/N, 1}.

This lemmas are analogous to (and derived from) corollary <ref> and lemmas <ref> respectively. Their proofs are in the supplementary material; that of lemma <ref> uses lemma <ref> and a criterion for stochastic dominance.


Lemma <ref> allows us to control type I error: for all α∈ [0,1], producing a p_A (y) and then rejecting H_0 (A) if and only if p_A (y) ≤α ensures that the type I error is less than or equal to ⌊α N⌋ + 1/N + 1≤α + 1 - α/N + 1. Conversely, rejecting H_0 (A) if and only if p_A (y) ≤α - 1 - α/N ensures a type I error less than or equal to α.

The residual terms 1 - α/N + 1 and - 1 - α/N are very small (smaller than the granularity of estimation of p_A (y), which is 1/N) and can be controlled by the user since N can be chosen as large as desired (at the cost of computational time). It is also possible to adapt the definition of p_A to formally eliminate these residual terms. With:

    p_A^+ (y) = 1/N+1( 1 + ∑_l=1^N 1{λ_A (y^(l) ) ≥λ_A (y) }),

the rejection criterion p_A^+ (y) ≤α ensures error control at level α. However, p_A^+ (y) is a biased estimator of p_A (y).




§ CASE OF GENERALIZED LINEAR MODELS
 



 §.§ Problem statement
 
Our algorithm is valid in the the linear model. We propose an adaptation of it to two discrete generalized linear models: the binary model y|X ∼Bernoulli( f(X β) ), and the Poisson model y|X ∼Poisson( f(X β) ). The link function f is the logistic function in the binary model and the exponential in the Poisson model. Unlike the linear one, these two models do not have a standard deviation σ, therefore we write β_A or β_A instead of θ_A or θ_A.

As in the linear case, we assume that the variables X_j, j ∈ A are active, and we want to know whether any variable outside of A is active. The hypotheses H_0 (A), H_1 (A), and the statistic λ_A (y), are defined the same way. Again, ideally we would like to estimate p_A^0 (y) = P_H_0 (A)( λ_A (Y) ≥λ_A (y) ), which is not possible because the distribution of Y under H_0 (A) is not known. As in the linear case, we propose an algorithm that computes p_A (y), an estimation of a statistic that approximates p_A^0 (y). This algorithm is similar to the algorithm <ref>, differing only in its calibration method. However, the statistic estimated by p_A (y) in the generalized linear models is not the same as in the linear case and does not have its theoretical properties.

The solution we chose in the linear case was to estimate the conditional p-value, p_A (y), defined in <ref>. Unfortunately, in discrete models this quantity is no more interesting because the set E of values that y can take is much more restricted: it is finite in the binary model (E = {0,1}^n) and countable in the Poisson model (E = ℕ^n) while the set Θ_A = ℝ^|A| of possible parameters is uncountable. One can generally expect the estimator function β_A, which sends a finite or countable set to an uncountable set, to be injective: only one vector, y, should produce exactly the estimate β_A (y). Therefore, the condition β_A (Y) = β_A (y) implies Y = y hence λ_A (Y) = λ_A (y) so necessarily p_A (y) = 1.

However, if n is large enough, the number of values that y (and thus β_A (y)) can take is very high (2^n in the binary model). So even if only y' = y satisfies β_A (y') = β_A (y), a large number of distinct vectors y' can satisfy β_A (y') ≈β_A (y). Therefore, the exact conditioning can — informally — be replaced by a "conditioning by β_A (Y) ≈β_A (y)". More formally, this amounts to using a probability distribution P_β_A (y) on the set E, producing random vectors Y compatible with H_0(A) and which satisfy β_A (Y) ≈β_A (y). This last condition is quantifiable in terms of mean square error:

    MSE( P_β_A (y)) = E_Y ∼P_β_A (y)[ ( X_A β_A (Y) - X_A β_A (y) )^2 ]

which must be very close to 0. We will then estimate:

    p_A (y) = P_Y ∼P_β_A (y)(λ_A (Y) ≥λ_A (y) )

which is the approximate equivalent in discrete models of the test statistic p_A (y) of the linear model.

In the linear model, we simulated vectors y^(l) following the distribution produced by conditioning on λ_A (Y) = λ_A (y) (a distribution we explicitly determined) and then used the population of y^(l) to estimate p_A(y) via Monte Carlo. We simulated vectors following this distribution by using a calibration function on vectors that were first simulated following the linear model under the null hypothesis. Similarly, in discrete models, we propose to simulate y^(l) following a distribution P_β_A (y) that meets the conditions above but is not explicitly determined, and to use these y^(l) to estimate p_A (y) via Monte Carlo. The simulation of vectors following P_β_A (y) is also carried out by applying a calibration procedure to vectors simulated according to the generalized linear model under the null hypothesis. This calibration is, however, more complex than that of the linear model.




 §.§ Calibration in nonlinear models
 

Given two parameter vectors β_A^(1), β_A^(2)∈ℝ^n, we want a calibration procedure from β_A^(1) to β_A^(2) such that:
[Desirable properties of calibration] 
For every y^(1)∈ E, denoting y^(2) as its calibrated version:


  <ref>.1.y^(2) should be "close" to y^(1) to preserve as much as possible its correlation structure with variables outside of A.

  <ref>.2.If β_A (y^(1)) = β_A^(1) then
β_A ( y^(2) ) ≈β_A^(2).


Condition <ref>.2 is the informal adaptation of the formal lemma <ref> satisfied by the linear calibration. The equality between parameters, which is a property too strong in discrete models, is replaced by an approximation.

In linear models, calibration is a function cal_β_A^(1)→β_A^(2) : E → E involving multiplication by and addition of non-integer constants. In discrete models, this is not possible because the calibration must produce a vector of integer numbers. Our solution is to replace deterministic non-integers with random integers having the same expected value. Therefore, calibration in nonlinear models is a random procedure described by the conditional distribution of the calibrated vector given the initial vector: Y^(2) | Y^(1)∼P_β_A^(1)→β_A^(2).



  §.§.§ One-step calibration algorithm in nonlinear models

Below we present a one-step calibration procedure, which follows a conditional distribution denoted P^(1)_β_A^(1)→β_A^(2). It is a step of the complete procedure that simulates vectors following P_β_A^(1)→β_A^(2). In both generalized linear models, it is based on the prediction vectors produced by the initial and target parameters: e^(1) = f ( X β_A^(1)) and e^(2) = f ( X β_A^(2)).





  §.§.§ Properties of the one-step calibration algorithm

The two non-exclusive scenarios of the binary case are consistent because if e^(1)_i = e^(2)_i then Y^(2)_i = Y^(1)_i regardless of the formula used.  In both models, E[Y^(2) | Y^(1)] = Z and if β_A^(1) = β_A^(2) (the case where no calibration is necessary) then e^(1) = e^(2) and Y^(1) = Z = Y^(2).

If an individual i satisfies e^(1)_i ≈ e^(2)_i, then the probability that Y^(2)_i = Y^(1)_i is high. If β_A^(1)≈β_A^(2), this is usually true, thus, as desired (condition <ref>.1.), calibration does not modify much of Y^(1).

Furthermore, the calibration "transforms a vector following the parameters β_A^(1) into a vector following the parameters β_A^(2)" in the sense that, with Y^(1) being random:

    E[Y^(1)] = e^(1)E[Y^(2)] = E[Z ] = e^(2).


However, this property does not imply condition <ref>.2., which is about the estimates β_A (Y^(1)) and β_A (Y^(2)), or, equivalently, the prediction vectors they produce: ê^(m) = f ( β_A (Y^(m)) ), m = 1 or 2. The condition can be verified in the simple case where the model is reduced to an intercept, that is, A = ∅. The prediction vectors e^(1) and e^(2) are then constants.

 
In single-parameter models, if β_A (Y^(1)) = β_A^(1) that is, ê^(1) = e^(1), then:


  * E[ ê^(2) | Y^(1)] = e^(2)

  * Var( ê^(2) | Y^(1)) ≤1/n| e^(2) - e^(1)| with equality in the binary model.



We prove it in the supplementary material. In single-parameter models, P^(1)_β_A (Y^(1)) →β_A^(2) thus produces a distribution of Y^(2) where ê^(2) is distributed around e^(2) with low variance. Therefore, by applying f^-1, β_A (Y^(2)) is distributed around the target value β^(2) with low variance. This result, combined with the parameter transfer property from β^(1) to β^(2) (<ref>) which is valid in multi-parameter models, suggests that the distribution of β_A (Y^(2)) under P^(1)_β_A (Y^(1)) →β_A^(2) has a similar property in multi-parameter models.



  §.§.§ Iterative calibration algorithm

The  variance inequality of lemma <ref> indicates that, at least in single-parameter models, ê^(2) (or β_A (Y^(2))) approaches its target with greater precision when the initial vector (ê^(1) or β_A (Y^(1))) is already closer to it. It is therefore useful to reapply the calibration procedure following a distribution P^(1)_· →β_A^(2) to an already calibrated vector. We thus define of P_β_A^(1)→β_A^(cal), the conditional distribution describing the complete calibration procedure, via the following simulation algorithm.





 §.§ Simulation-calibration testing algorithm in generalized linear models


With the calibration algorithms <ref> and <ref>, we can propose the algorithm for generating the empirical conditional p-value p_A (y) in binary and Poisson generalized linear models. It is similar in form to algorithm <ref> for linear models, with two main differences:


  * theoretically, the quantity estimated by p_A (y) is not the same conditional p-value, p_A(Y) instead of p_A(Y) (see <ref>);

  * algorithmically, the calibration step is more complex, using the sub-algorithms <ref> and <ref>.




p_A (y) is intended to be interpreted the same way in linear and generalized linear models. However, the properties relating to its distribution under the null hypothesis (lemmas <ref> and <ref>) are only demonstrated in the linear case. The construction of p_A (y) in generalized linear models is designed so that these properties are approximately valid there.



§ VARIABLE SELECTION PROCEDURE
 

The simulation-calibration test measures the significance of a covariate selected by the Lasso. When no active variables are known a priori and a complete model needs to be selected, it is necessary to use the test iteratively.



 §.§ Notations and algorithm

For every k ∈{1, .., p}, let j_k be the index of the k-th variable selected by the Lasso and let A_k = { j_1, .., j_k } be the set of the first k variables selected, with A_0 = ∅. Each j_k is the first variable selected outside the set A_k-1, that is, following to the notation introduced in <ref>, j_k = j_A_k-1. This is a slight variant of the Lasso path where variables do not "reenter": if there is a λ such that β_j^Lasso (λ) ≠ 0 (which translates into a k such that j ∈ A_k), then by definition j belongs to all A_k', k' > k even when there are λ' < λ such that β_j^Lasso (λ') = 0.

The simulation-calibration test of H_0 (A_k-1) thus measures the significance of the variable j_k. Furthermore, since variables do not "reenter" the sequence of sets, the null hypotheses H_0 (A_k) become progressively weaker: if k < k', then A_k ⊂ A_k' and since the hypotheses concern the complements of these sets (on which they assert the nullity of β), H_0 (A_k) implies H_0 (A_k').





 §.§ Choice of the halting criterion


In this procedure, we carry out several hypothesis tests. However, since the procedure halts at the first test where the null hypothesis is not rejected, it cannot be considered a multiple testing procedure per se, in which a potentially high number of tests satisfy their null hypothesis. Therefore, it is not advisable to make the rejection threshold more stringent as done by procedures like Bonferroni and Benjamini and Hochberg <cit.>, which are used to control the risks of false positives deriving from the large number of p-values generated under the null hypothesis.

However, it is possible to adapt the halting criterion of the procedure to the fact that it evaluates a sequence of ordered tests. In the general problem of a sequence of ordered tests whose p-values p_1, .., p_m are measured, <cit.> proposed the ForwardStop criterion. This rejects the first k̂_F tests where:

    k̂_F    =   max{ k ∈{1, .., m } : p_k^FS≤α}
    
    p_k^FS   =    -1/k∑_i = 1^klog(1 - p_i).


This criterion has the advantage that, in the vast majority of cases, one can be compute it knowing only the first values of the sequence (p_k), because the sequence (p_k^FS) is calculated from p_i, i ≤ k and it generally increases with k. In practice, when we apply ForwardStop, we select the first k̂_F' variables, where:

    k̂_F' = min{ k ∈{1, .., m } :  p_k^FS > α} - 1

which amounts to replacing the condition to continue the algorithm at step 2 with p_k^FS≤α.

In contrast to ForwardStop, the simple criterion p_k ≤α is called thresholding. Thresholding controls both the FWER and the FDR at level α <cit.>, while ForwardStop, generally less conservative, controls the FDR at level α <cit.>.

These control results apply to the type I error, that is, incorrect rejection of the null hypothesis when it is true. However, in the sequential procedure, we often test H_0 (A) hypotheses that are not true: notably at the first step of the procedure, where A = ∅ and it is sufficient that there is an active variable for H_0 (A) to be incorrect. It is possible that there are active variables outside of A but that they are selected by the Lasso only at low values of λ and that the first variable selected outside of A is an inactive variable. In this case, retaining this variable constitutes a false positive from the perspective of variable selection, but not a type I error since H_0 (A) is not true. Section <ref> presents a result controlling this occurrence.



§ EXTENDED THEOREM: CONTROL OF THE SELECTION ERROR
 
We divide the set of covariates into three disjoint subsets: {1, .. ,p } = A ∪ B ∪ C with ∀ j ∈ C, β_j = 0.

We perform the simulation-calibration test of the hypothesis H_0 (A) with the goal of variable selection, that is, if the test rejects H_0 (A), we select the (almost certainly unique) variable j_A (Y) ∈ B∪ C such that β_j_A (Y)^Lasso (λ) ≠ 0 for λ near λ_A. This selection is a true positive if j_A (Y) ∈ B and a false positive if j_A (Y) ∈ C.
Since H_0 (A) is not true (because there may exist j ∈ B such that β_j ≠ 0), the results from the previous sections do not apply. However, we have the following result:

Assume that the active variables are orthogonal to the inactive variables, that is:

    X_C^T X_A∪ B=0.

Then, the simulation-calibration test of H_0 (A) at level α has a probability less than α + 1 - α/N+1 of selecting a false positive.

We prove it in the supplementary material.


Like under the null hypothesis (section <ref>), this result can also be understood as controlling at the level α the risk of selecting a false positive if the selection criterion is p_A (y) ≤α - 1 - α/N.


The theorem renders applicable the properties of control of the FWER and FDR by the variable selection procedure with thresholding and ForwardStop respectively, provided that the active variables are orthogonal to the inactive variables.




§ SIMULATION STUDIES
 
To assess the performance of the test by simulation-calibration, we conducted two simulation studies.

First, we measured the distribution of the p-value generated under H_0 (A) in a large number of different scenarios. By construction, this distribution is supposed to be the uniform on [0,1]. More precisely, we proved in Lemma <ref> that in the linear model, under the continuity assumption, p_A (Y) follows the uniform distribution on a discrete analogue of [0,1], and in Lemma <ref> that by relaxing the continuity assumption, in the linear model p_A (Y) stochastically dominates the uniform distribution. However, we do not have equivalent theorems in generalized linear models, although the properties of nonlinear calibration demonstrated in section <ref> suggest that a distribution close to the uniform is also to be expected in these models.

We also measured the performance of the variable selection procedures —– with both the thresholding and ForwardStop halting criteria —– using three usual metrics of variable selection: family-wise error rate (FWER), false discovery rate (FDR), and sensitivity. Due to the properties of thresholding and ForwardStop, we expect control of the FWER with the former and of the FDR with the latter, at least under the conditions where the extended theorem applies: in the linear model and with no correlation between covariates. Additionally, we compared these performances with those of an equivalent procedure based on the CovTest by <cit.>.



 §.§ Simulation plan

The following simulation plan is shared by the study of the p-value's distribution under the null and by that of the variable selection procedure. We simulate n_sim = 500 data sets for each of the 252 parameter sets (or scenarios). In all cases, the number of observations is n = 1000 and the number of covariates is p = 500. The parameters which vary across scenarios are:


  * The type of model: linear, binary with dense data, binary with sparse data, or Poisson. In binary models, dense data means that E[Y|X = 0] = 0.5 and sparse data that E[Y|X = 0] = 0.1, with the value of β_0 distinguishing these two cases.

  * The correlation matrix used to simulate the regressors: a Toeplitz matrix of coefficients ρ_(i,j)=ρ^| i - j |, with ρ=0, ρ=0.9, or ρ=0.99. These high values are used by <cit.> in their simulation study.

  * The number of active variables: 0, 1, 2, 5, or 10. They are drawn uniformly among the 500 covariates.

  * For scenarios with at least one active variable, the signal-to-noise ratio (SNR). This quantity is inspired by the SNR used by Sabourin et al. Unlike them, we use an empirical version which is defined for each generalised linear model. This is the ratio between the empirical variance of the signals (E_β[Y_i | X], i = 1, ... , n) and the empirical mean of the variances of the noises (Var_β(Y_i | X), i = 1, ... , n): E_β[Y_i | X], i = 1, ... , n):

    SNR(X, β) = 1/n-1Σ_i=1^n (E_β[Y_i | X] - 1/nΣ_i=1^n E_β[Y_i | X] )^2 /1/nΣ_i=1^n Var_β(Y_i | X).
 
A higher SNR means that the impact of each active variable on Y is more easily observable, which should make variable selection more effective. Scenarios with 0 active variables necessarily have a zero signal-to-noise ratio. We set SNR(X, β) = 1, 0.3, 0.1, 0.03,  or 0.01.




 §.§ p-value under the null hypothesis
 

To verify that the test follows its expected behavior under the null hypothesis, we assume that the set A of 0 to 10 active regressors is known and we test H_0 (A). For each scenario, in each of the n_sim = 500 data sets s characterized by A_s, X_s, and Y_s, we produce a p̂_s = p_A_s (Y_s, X_s) using the simulation-calibration algorithm based on N = 100 simulations of calibrated response vectors.


In one of the scenarios (linear model, ρ = 0.99, 1 active regressor, SNR = 1), as an example, we also produced in each of the 500 data sets a naive estimate of the unconditional p-value p_A^0 (y) (see section <ref>). It is obtained by Monte Carlo without the calibration step, i.e. by applying the algorithm <ref> without the calibration step (2.3.), taking y^(l) = y^sim with θ_A^sim = θ_A (y). This intends to illustrate the impact of calibration on the distribution of the p-values.


In each scenario, we observe whether the population of ( p̂_s )_1 ≤ s ≤ n_sim is distinguishable from a sample drawn the uniform distribution on [0,1]. Graphically, agreement or disagreement with the uniform distribution is observed on quantile-quantile (q-q) plots where for each s = 1, .., n_sim, the s-th smallest value p̂_(s) is represented at the coordinates ( s/n_sim, p̂_(s)).


The diagrams in the supplementary material illustrate the necessity of the calibration step for the empirical p-values produced to be valid. In this example (linear model, ρ = 0.99, 1 active variable, SNR = 0.1), p-values produced by simulation-calibration are compatible with the uniform distribution. In contrast, p-values produced without calibration deviate significantly, the smallest p-value among 500 being 0.13. At usual test levels, false positives are thus practically impossible instead of being possible with a controlled probability, suggesting a very low power of the test. At higher test levels, type I error is no longer controlled.


To systematically evaluate the adequacy of the p̂_s to the uniform distribution across all scenarios, we used the Kolmogorov-Smirnov (K-S) test on the population of ( p̂_s )_1 ≤ s ≤ n_sim in each scenario. We conducted the test in its two-sided version, where the null hypothesis is that the population considered is a sample from the target distribution (here the uniform distribution on [0,1]), and in one of its one-sided versions having a less strict null hypothesis: that the population is sampled from a distribution which has stochastically dominates the target distribution. We are interested in this one-sided test because, as seen in the conclusion of section <ref>, basing a test on an empirical p-value whose distribution dominates the uniform distribution on [0,1] allows control of its type I error.


Results vary across the model type: linear, dense binary, sparse binary, or Poisson. For each of the four types and each variant of the Kolmogorov-Smirnov test, we used the Bonferroni correction on all the p-values of the 63 K-S tests applied to the scenarios of that model type.





In the linear and dense binary models, the bilateral K-S tests do not reject the hypothesis that each population of p-values follows a uniform distribution over [0,1]. The smallest p-value is 0.0112 among the linear scenarios and 0.0546 among the balanced binary scenarios, which, combined with the Bonferroni correction for 63 tests, does not lead to rejection at usual test levels (63 × 0.0112 = 0.706). Figure <ref> illustrates this general adequacy with the uniform distribution across eight example scenarios.


In the sparse binary model, there are scenarios where the distribution of generated p-values significantly deviates from the uniform distribution and does not dominate it. At the 0.1 level and with the Bonferroni correction, the null hypothesis of the K-S test (unilateral as well as bilateral) is rejected in two out of 63 scenarios: those with maximum correlation among covariates (ρ = 0.99) and the number of active regressors (|A| = 10), and where the signal-to-noise ratio is high (SNR = 0.3 or SNR = 1). The p-values of the bilateral K-S test in these scenarios are respectively p = 2.86 × 10^-6 = 1.81 × 10^-4 / 63 and p = 2.55 × 10^-7 = 1.61 × 10^-5 / 63. Despite this non-adequacy, the q-q plots (Figure <ref>, which includes the SNR = 1 scenario where the deviation from the uniform distribution is most pronounced) show that P_H_0 (A)(p̂_A ≤α) ≤α at low values of α, the opposite occurring only when α exceeds about 0.4. Then, in practice, we have a more conservative control of the type I error than the nominal level at usual test levels. In the most divergent scenario, the diagram even shows a type I error significantly lower than its expected level, signaling a loss of test power. However, this effect is considerably less pronounced than in the uncalibrated simulation example.


In the Poisson model, the unilateral K-S tests do not reject the hypothesis that each population of p-values stochastically dominates the uniform distribution over [0,1] (the minimum p-value across unilateral tests is 0.0214 = 1.35/63). Thus, FWER control is at least as conservative as the nominal level at all test levels. However, at the 0.1 level and with the Bonferroni correction, the bilateral K-S tests reject the hypothesis of a uniform distribution of p-values in 5 out of 63 scenarios. These scenarios present a different profile from those where deviations from the uniform distribution were observed in the sparse binary case: 4 out of 5 are scenarios with zero correlation, all have an SNR of 0.1 or less, and the minimum p-value (p = 1.79 × 10^-6 = 1.13 × 10^-4 / 63) is reached at zero correlation, minimum signal-to-noise ratio (SNR = 0.01), and |A| = 10 known active regressors. Moreover, as seen in Figure <ref>, where the scenario with the minimum K-S test p-value appears, even then the deviation from the uniform is small, especially at the small quantiles. The actual level of first type error is therefore very close to its nominal level at its usual values.



 §.§ Variable selection procedure


In this simulation study, unlike in section <ref>, we did not assume the set of active regressors to be known. This set is estimated by the sequential procedure described in algorithm <ref>. The number of simulations in the calculation of empirical p-values is higher, N = 500, which allows for a more precise estimation of p-values.


At first, we applied the algorithm with the thresholding halting criterion at α = 0.95, a high level that produces on each dataset a relatively complete sequence of variables that may be selected, along with their associated p-values. In a second step, we used the variable selection procedure on these sequences of p-values for each α in a relatively dense grid of values (from 0.01 to 0.5, with a step of 0.01), with both thresholding and ForwardStop. In this second step, we observe the procedure's performance on a large number of values of α in a computationally affordable manner since it does not require recomputing the empirical p-values through simulation-calibration, as this was done in the first step.



  §.§.§ Control of FWER and FDR


Figure <ref> and <ref> show how the FWER and FDR of the thresholding and ForwardStop procedures vary with α in the same 8 example scenarios as in Figure <ref>: 10 active variables, either zero or maximal correlation between covariates, minimal or maximal SNR, linear or dense binary model. In the the supplementary material, Figures 1 and 2 of  show the variation of the FWER and FDR of both procedures with α in the 8 scenarios of Figure <ref> (sparse binary and Poisson models with the same other parameters) while Figures 3 and 4 display the variations of their sensitivity in those 16 scenarios.





Although theoretical results on controlling the FWER or FDR are valid only in the linear model, the shape of these curves varies little across model types. In scenarios where it is easiest to detect a signal, i.e., those with high SNR and especially among those, the uncorrelated scenarios, it is observed that the ForwardStop criterion is significantly less conservative than thresholding at an equal α, with higher FWER and FDR. This is explained by a substantial number of small p-values among the first values of the sequence (p_k), corresponding to easily detected variables, which, by averaging, pull down the quantity p_k^FS at higher k, making it significantly lower than p_k.


In all uncorrelated scenarios, as expected, FWER control is very close to α with the thresholding criterion (blue curves in Figure <ref> are close to or below the bisectors). Additionally, in these scenarios, we observe FDR control at a level close to or lower than α with the ForwardStop criterion (red curves in Figure <ref>). Due to the more conservative nature of thresholding compared to ForwardStop, thresholding also controls the FDR below α in uncorrelated scenarios, but ForwardStop does not control the FWER. The gain in sensitivity of ForwardStop compared to thresholding is modest, reaching a maximum of 0.127 across all 252 scenarios.


On the other hand, in scenarios with many active variables and strong correlation among covariates, neither the FWER nor the FDR is controlled by either halting criterion. The strong correlation between active and inactive variables often allows a nominally inactive variable to be selected by Lasso before the corresponding active variable. The selected variable can then have a statistically significant association with the response which is a manifestation of the actual association between the active variable and the response. This cannot be attributed to the nominally active variable, because it is not known and not detected by the Lasso. The status of these "false positives" carrying a statistical signal has been discussed in the literature <cit.>.



  §.§.§ Comparison with Lockhart et al.'s Covariance Test
 

To evaluate the sensitivity of the variable selection procedure via simulation-calibration, it is useful to compare it to another selection method aiming for the same objectives. This competitor is the covariance test (CovTest) by <cit.>, which, like the simulation-calibration test, measures the significance of variables entering the Lasso path by assigning each one a p-value.


Due to reasons of implementation of CovTest, we focused on the 63 linear model scenarios. We used CovTest with the equivalent of algorithm <ref>, with the thresholding halting criterion at α = 0.05, which is supposed to control the FWER at this level. Figure <ref> shows the FWER and sensitivity of both methods across all linear model scenarios.



We observe a failure to control the FWER in correlated scenarios by both CovTest and simulation-calibration (points ∘ and × in figure <ref>). In these scenarios, simulation-calibration tends to produce false positives more frequently than the CovTest when the signal-to-noise ratio (SNR) is low, but less frequently when it is high (see the impact of SNR in figure 6 of the supplementary material). However, the FWER is controlled at 5% or less in the 21 linear scenarios with no correlation, both by the CovTest and by simulation-calibration (points ∙), but with a clear difference between the two methods.


The FWER observed with the CovTest selection procedure falls well below its nominal level in certain scenarios, even reaching 0 — meaning no false positives observed in 500 simulations — in the scenario with 10 regressors, SNR = 0.3.


In contrast, the FWERs observed with the simulation-calibration procedure across the 21 uncorrelated scenarios all hover close to their nominal levels, ranging from 0.042 to 0.072. This implies that the number of simulated data sets on which the procedure produces at least one false positive ranges between 21 and 36 out of 500. Given the Bonferroni correction, this is consistent with the hypothesis that the true FWER is equal to α in each of these scenarios, i.e., that the number of data sets exhibiting at least one false positive in each scenario follows a Binomial(500, 0.05) distribution. Indeed, the p-value associated with the largest of the 21 observed FWERs is 0.0196 = 0.412 / 21.


This stronger conservatism of the CovTest procedure in uncorrelated scenarios results in notable gaps in sensitivity to the advantage of the simulation-calibration procedure. The sensitivity difference is positive or zero in each of these 21 scenarios, exceeds 0.1 in 6 of them (all with at least 5 active regressors and an SNR of at least 0.1), and reaches a maximum of 0.388 in the scenario with 10 regressors, SNR = 0.3. Therefore, the simulation-calibration test represents an improvement over the CovTest under the ideal condition of no correlation between covariates, with substantial gains in sensitivity allowed by a controlled increase in FWER that does not significantly exceed its nominal level.



§ APPLICATION TO PHARMACOVIGILANCE DATA
 

We illustrated the sequential variable selection procedure by using it on data from the national pharmacovigilance database (BNPV). We used the same data preprocessing as described in <cit.>, producing a database of n = 452914 spontaneous notifications of adverse drug reactions from January 1, 2000, to December 29, 2017, with 6617 distinct adverse events (coded at the Preferred Term level of the Medical Dictionary for Regulatory Activities, MedDRA) and p = 1692 distinct drugs (coded at the 5th level of the Anatomical Therapeutic Chemical hierarchy) reported at least 10 times. We focused on a binary outcome, the adverse event "Drug-Induced Liver Injury" (DILI). It is one of the most frequent adverse events with 25187 occurrences, accounting for 5.56% of all spontaneous notifications. We used a logistic regression model on the drug exposures, which are binary covariates. To reduce computation time, we used only N = 50 simulations in the computation of empirical p-values.


When a variable is selected with a positive coefficient estimate, we considered it to be a pharmacovigilance signal. To assess the performance of our method, as in <cit.>, we used the reference set DILIrank of known pharmacovigilance signals related to drug-induced liver injuries <cit.>. It includes 203 negative controls (drugs known not to be associated with DILI) and 133 positives (drugs known to be associated with DILI).


Given initial results, it appeared necessary to preprocess the data before performing the Lasso. Figure 7 in the supplementary material shows that the p-values estimated by simulation-calibration are zero for the first 14 variables on the Lasso path, then distributed around 0.1. This sudden change coincides with a peculiarity in the correlation structure of the exposure matrix X: the 14th exposure selected, trimethoprim (ATC code J01EA01), is exceptionally correlated (ρ = 0.9998) with another exposure, sulfamethoxazole (J01EC01), a drug with which trimethoprim is almost always co-prescribed.





We removed from the matrix X the 6 covariables that had a correlation of at least 0.9 with any of the variables on the Lasso path (computed before this removal), and re-used the Lasso and the sequence of simulation-calibration tests on these purged data (p = 1686). The Lasso yields identical results since the removed variables were not selected by it. However, p-values estimated by simulation-calibration (Figure <ref>) are lower starting from the 15th variable, and no longer show a regime change at that point. They do show a comparable regime change further along the Lasso path, after the p-value associated with the 86th selected variable (P01BDXX, diaminopyrimidines). This variable is also correlated with another variable in the database (P01BD01, pyrimethamine), but at a lower level: ρ = 0.307, similar to the other correlations involving selected variables.


Figure <ref> shows that using ForwardStop instead of thresholding, in addition to resulting in less conservative selection at the same α, has a stabilizing effect: while the empirical p-values exhibit fluctuations due to the limited number of simulations from which they are estimated, the quantity of interest p^FS on which ForwardStop is based presents a much smoother profile, being calculated by averaging from the empirical p-values.



Table <ref> shows the performances of the variable selection procedure — with the thresholding or ForwardStop halting criterion, at the same level α = 0.2 — calculated from the status of drug exposures listed in the DILIrank reference set <cit.>. Only exposures where Lasso estimates a positive association with DILI are considered signals. The status of the exposures is also indicated in figures <ref> and <ref> allowing visualization of performances at other levels of α. Both approaches are markedly more conservative than the information criteria (including ForwardStop, which itself is less conservative than thresholding). This allows them to achieve a lower false positive rate at the expense of lower sensitivity. Daunorubicin (L01DB02) is the sole false positive according to DILIrank.



§ DISCUSSION

We have proposed a test for the significance of variables entering the Lasso path, which <e use sequentially to select a model. It tests the null hypothesis H_0 (A) which states that a known set of variables—A—includes all the active variables, and it focuses on λ_A, the largest value of λ at which a variable not belonging to A is selected. Rejecting the test means selecting this variable.


We built the test to circumvent the difficulty of using λ_A directly as a test statistic. Just as a p-value is the probability under the null hypothesis that the test statistic exceeds its observed value, the statistic we consider, p_A, is a conditional probability that λ_A exceeds its observed value: conditional on the correlation structure relating the response Y to the already selected variables X_A, given by θ_A. This is why we call it a conditional p-value. We estimate it by the Monte Carlo method provided by algorithm <ref>: simulation-calibration. It involves simulating response vectors under the null hypothesis and then calibrating them on the condition on θ_A. This distributes the response vectors following the conjunction of the null hypothesis and this condition, which allows using them to estimate the conditional p-value.


This can be seen as a generalization of the permutation selection by <cit.>, which simulates a population of response vectors with the same distribution as the observed y, and considers the population of the λ at which a variable enters the Lasso path computed on these simulated vectors. In permutation selection, the focus in on the λ at which the very first variable enters the Lasso path, i.e., in our notation, λ_∅. Permutation is also a way to impose conserving θ_∅ = (β_0, σ̂), the empirical mean and standard deviation of y. Permutation selection resembles the test of H_0 (∅) by simulation-calibration in these aspects. Since permutation selection retains the median of λ_∅, using it to decide only on the selection or not of the first variable on the Lasso path would essentially amount to conducting this test at α = 0.5 (the only difference being in the method of calibrated simulation, permutation, or post-simulation calibration of arbitrary vectors).


We have proven in sections <ref> to <ref> the validity of our method in the case of the linear model. The conditional p-value p_A (Y) under the null hypothesis follows a distribution which stochastically dominates the uniform distribution over [0,1] (lemma <ref>), and its estimator p_A (Y) is unbiased, consistent with the number N of simulations chosen to perform (lemma <ref> and its consequences), and itself dominates the uniform distribution over the discrete set of values it can take (lemma <ref>). This means that rejecting the null hypothesis based on a threshold α on p_A (Y) controls the type I error rate at this level α, with a small residual term. From a variable selection perspective, this type I error means the selection (necessarily erroneous) of an additional variable while all active variables are already selected, i.e., belong to A. Furthermore, we have proven (section <ref> and supplementary matereial) that even when some active variables have not been selected, the risk of selecting an inactive variable by simulation-calibration is controlled at this same level. This latter result requires orthogonality between active and inactive covariates. The control of the error under H_0 (A) does not depend on any assumptions about the correlation structure.


In discrete generalized linear models (binary and Poisson), the theoretical framework is different since there is no conditional p-value that follows or dominates the continuous uniform distribution over [0,1]. Nevertheless, the production of an empirical p-value by simulation-calibration adapts well to these models. Theoretically, p_A (Y) is seen as the estimator of a probability whose definition approaches that of the conditional p-value (section <ref>). Practically, the "calibration" part of the algorithm becomes more complex: it is iterative and stochastic (algorithm <ref>) while linear calibration was a simple affine function. Although it is not guaranteed by a theorem, the simulation study (section <ref>) shows that the distribution of empirical p-values is indistinguishable from the uniform distribution in the majority of non-linear simulation scenarios and deviates slightly in others, without any scenario of failure to control the type I error at usual levels of α.


Simulation-calibration has the disadvantage of requiring a computation time that can be significant under certain circumstances. It involves, for each of the N simulated vectors, adjusting the unpenalized model restricted to A in order to perform the calibration, then applying the Lasso to the entire set of covariates to the calibrated vector. In nonlinear models, calibration is itself iterative and requires repeatedly readjusting the unpenalized restricted model. In the application to the BNPV data, the conditions for a very long computation time were met: large data (452914 by 1686), a binary model, and especially an A set of selected variables that reaches a significant size (over 100). It was therefore necessary to repeatedly adjust relatively large-dimensional logistic models without penalty. By contrast, the simulation study was carried out with a maximum of 10 active variables, so the number of selected variables did not exceed this order of magnitude, which limited the computation time.


Our test is comparable, in its objectives, to the covariance test (CovTest) of <cit.>, which is also a test of significance for a variable preselected by the Lasso. We focus on the distribution of λ_A alone and not, like Lockhart et al., on the evolution of the Lasso between two consecutive values of λ where a variable enters the Lasso. Therefore, we avoid a situation of low power of the covariance test: when two active variables enter the Lasso at values of λ close from each other, the CovTest's test statistic is the difference between the same quantity computed on Lasso results which differ only a little; this can lead artificially to selecting none of the two active variables even if both are significantly associated with the outcome. This could explain simulation-calibration's better sensitivity in the absence of correlation between covariates that we observe in section <ref>. In the presence of correlation, however, the performance of simulation-calibration is not always better than that of the selection by CovTest, with lower FWER at equal sensitivity in high SNR scenarios, but higher FWER at equal sensitivity in low SNR scenarios (Figure <ref> and, in the supplementary material, Figure 6).


Using the procedure on data from the French national pharmacovigilance database has highlighted an undesirable phenomenon that could help explain the poorer performance in some cases of correlation. We observe in Figure <ref> a change of pattern in the sequence of estimated p-values where, once a certain variable j_0 belongs to A, p_A (y) is always relatively high. The detail of the Lasso performed in the simulation-calibrations provides an explanation for this phenomenon. We observe that when it occurs, there is a variable j_1 ∉ A correlated with j_0 and possibly with other variables in A such that, among the population of λ_A (y^(l)) obtained by simulation-calibration (see algorithm <ref>), a significant proportion are high values associated with an early selection of j_1 by the Lasso applied to y^(l). Indeed, calibrating the y^(l) on one or more associations with variables correlated with j_1 can lead to an association with j_1 itself, which is captured by the Lasso. These high λ_A (y^(l)) lead to the estimation of a high p_A (y).


This phenomenon reflects the influence of j_A on the distribution of λ_A under the null hypothesis. It might be more relevant, from the perspective of the test's power, to compare the observed λ (y), not as we have done with the distribution of λ_A (Y) conditional on θ_A (Y), but with its distribution under a double conditioning by θ_A (Y), j_A (Y), or more generally taking into account j_A (Y). In practice, this would involve introducing a weighting in the calculation by averaging p_A (Y).


This phenomenon results, when j_0 ∈ A, in a loss of power in the selection of all variables, including those that are not correlated with j_0 or j_1. This loss of power in certain cases of correlation could offset the tendency of the test by simulation-calibration, observed in simulations without correlation, to be more powerful but less conservative than the CovTest which would explain why in some correlated scenarios, variable selection by simulation-calibration is less conservative than that by CovTest while it has practically the same sensitivity.


Thesimulation-calibration procedure's lower FWER compared to that by CovTest in high-SNR correlated scenarios could be explained by the fact that due to the high SNR, it is more common for all active variables to be selected before inactive variables on the Lasso path. When this is the case, H_0 (A) is verified (thus the risk of false positive is controlled) in all tests of the iterative procedure that are likely to produce a false positive. For the risk control to fail, it is necessary both for the active variables to be mixed with inactive variables on the Lasso path – for there to be an iteration where an unverified H_0 (A) is tested which can nevertheless lead to the selection of an inactive variable –, and for active variables to be correlated with inactive variables – so that the extended theorem of control of the selection error (section <ref>) does not apply.



