A Sequential Rejection Testing Method for High-Dimensional Regression with Correlated Variables

1.1 Relation to other work

Our proposed method is based on the multi sample splitting method from Meinshausen et al. [7] and the sequential rejection principle of Goeman and Solari [17]. It is a generalization and power improvement over the multi sample splitting technique for inference of single variables [7] and for hierarchically ordered clusters of variables [15]. The improvement in power is strict, and in analogy to the gain of power of Holm’s procedure [23] over the Bonferroni adjustment. Thus, even if the increased power might be only small for some datasets, one cannot do worse with the new procedure. The only price to pay is a slightly more complicated algorithm: we provide an implementation in the R-package hdi. We note that the mentioned improvement, based on the work by Goeman and Solari [17] and Goeman and Finos [18], is not entirely trivial to derive in the setting of the multiple sample splitting scheme. The modular set-up based on four building blocks presented here is addressing the issue how the improvement can be constructed.

1.2 Outline of the paper

In Section 2 we describe the four basic building blocks of the method and the assumptions that are sufficient to establish in Section 3 its strong FWER control. In Sections 4.1 and 4.2, respectively, we focus on the inference of two specific kinds of cluster collections: singletons and hierarchically ordered clusters. In Section 4.3 we show how logical relationships can be used to improve the power. Finally, we provide in Section 5 a comparison based on empirical results for error control and power, with a focus on minimal true detections, and we apply the new method to a real dataset.

2.1 Screening of variables

We consider variable screening where an estimator Sˆ(b)⊆{1,…,p}, based on data corresponding to Nin(b), is aiming at including all active variables S0. A prime example is the Lasso [19], while a detailed empirical comparison of five popular screening procedures can be found in Bühlmann and Mandozzi [20]. Assume that the screening procedure satisfies the following properties for any sample split b:

The sparsity property in (A1) implies that for each sample split b it holds that |Sˆ(b)|<|Nout(b)|, a condition which is necessary for applying classical tests as described in Section 2.2 below. The δ-screening property in (A2) ensures that all the relevant variables are retained with high probability (where δ>0 is typically small). We indicate in Section 3.1 that under some assumptions, the Lasso satisfies (A1) and (A2).

2.2 Testing and p-values

The idea is to perform a classical statistical test on the other half sample from Nout(b) in a low-dimensional problem with variables from Sˆ(b) only. For each sample split b, based on the second half of the sample corresponding to Nout(b), consider a testing procedure, e.g. the classical partial F-test (see also Section 3.1), that provides correct p-values pC,(b) for the null hypothesis H0,C∩Sˆ(b) for each screened set Sˆ(b), in the sense that for each nominal level α∈(0,1)

We note that the probability is with respect to the data generating random variables corresponding to the second half Nout(b), and the null-hypothesis is fixed with respect to Nout(b). Due to the screening property (A2), when δ→0, the null-hypothesis H0,C∩Sˆ(b) approximates the unconditional hypothesis H0,C which we aim to test for. If C∩Sˆ(b)=∅ define pC,(b)=1. This provides a (correct) p-value pC,(b) for each cluster C∈c and each sample split b∈{1…B}.

2.3 Multiplicity adjustment

Consider for each sample split b and each cluster C∈c a multiplicity adjustment procedure mC(b):2c→[1,∞] that for each collection r of rejected clusters provides a multiplicity adjustment mC(b)(r)≥1 and satisfies the following properties:

where we define 1/∞=0. Such a family of multiplicity adjustments for b=1,…,B are often naturally induced from a global multiplicity adjustment procedure mC.

2.4 Aggregation of p-values

Consider a collection of screened sets of variables Sˆ(b), a cluster C∈c, a collection of p-values pC,(b) for the null-hypothesis H0,C∩S(b) (which approximates H0,C, see comment after (A3)) and a collection of multiplicity adjustments mC(b)≥1 (we drop here the dependence on r).

The goal is to aggregate the p-values pC,(1),…,pC,(B) to a single p-value which is adjusted for multiplicity. An aggregation procedure is a monotone increasing function aggr:[0,1]B→[0,1]. Assume it satisfies the following property:

2.5 The procedure

Our procedure is based on the four building blocks above. First, we proceed with screening of the variables based on the first half sample from Nin(b) (Section 2.1), e.g., in Section 5.1 we use the Lasso with regularization parameter chosen by 10-fold cross-validation (see also Section 3.1). Then, we construct the p-values based on the second half sample from Nout(b) by using the partial F-test (Section 2.2 and see also Section 3.1). This leads to a (correct) p-value pC,(b) for each cluster C∈c and each sample split b∈{1…B}.

The multiplicity adjustment is done sequentially (Section 2.3). Based on a chosen significance level α∈(0,1) and for a collection of currently rejected sets r, define the successor of r as

Start from “no rejections” r0=∅, define ri+1=ri∪n(ri) and r∞=limi→∞ri (although r∞ is never constructed due to finite-ness of all possible subset of the variables). Concrete choices of mC(1)(r),…,mC(B)(r) are discussed in Section 4.

Finally, we aggregate the p-values as indicated in Section 2.4. Concrete aggregation methods are described in Proposition 1 in Section 3.1.

3.1 Screening, testing and aggregation: their properties

We discuss here some choices for screening, testing and aggregation which we use in the implementation in the R-package hdi. The issue of sequential multiplicity adjustment is treated separately in Section 4.

For variable screening, we use the Lasso with regularization parameter chosen by 10-fold cross-validation. Theoretical justification of the sparsity and screening property (A1) and (A2) can be derived by assuming a compatibility or restricted eigenvalue condition on the fixed design matrix X and a beta-min assumption requiring that minj∈S0|βj0|≫|S0|log(p)/n is sufficiently large: we refer to Bühlmann and van de Geer ([13], Ch. 2.7 and Ch. 6) for the details.

For construction of the p-values (in the low-dimensional setting, due to variable screening in the first half of the sample) we use the partial F-test. Then, assuming fixed design X and Gaussian errors, condition (A3) holds.

For aggregation of the p-values, ensuring that (A6) holds, we have the following result for two slightly different methods.

A proof, which was basically given in Meinshausen et al. [7], can be found in the Appendix.

4.1 Inference of single variables

This first example is paradigmatic for the advantages of the modular approach: a simple improvement of the multiplicity adjustment procedure allows for a better power, basically in the same way as in a low-dimensional setting in Goeman and Solari [17].

Concretely, we consider the problem of inferring single variables, i.e., the collection of clusters c={{i};i=1,…,p}. The method proposed in Meinshausen et al. [7] corresponds to the method of Theorem 1 with the aggregation procedures of Proposition 1 and the following Bonferroni-based [21, 22] multiplicity adjustment procedure:

As the multiplicity adjustments are independent from the (previously) rejected collection of sets, the monotonicity property (A4) is trivially satisfied, while the single-step property (A5) follows from

The power of the method can be improved taking instead of (2) the following Bonferroni-Holm-based [23] multiplicity adjustment procedure:

The monotonicity property (A4) is still satisfied since

for r⊆s, whereas

proves the single step property (A5).

4.2 Inference of hierarchically ordered clusters of variables

When dealing with the challenge of inferring hierarchically ordered clusters of variables, e.g. from the tree-structured output of a hierarchical clustering algorithm, one considers a collection of clusters c={Ci}i where for any two clusters Ci,Ci′∈c, either one cluster is a subset of the other, or they have an empty intersection. The method proposed in Mandozzi and Bühlmann ([15], Section 2), which is based on the procedure of Meinshausen [24], corresponds to the one as in Theorem 1 with the aggregation methods of Proposition 1 and the following multiplicity adjustment:

Here, anc(C) denotes the ancestors in a hierarchically ordered cluster tree. To check the monotonicity property (A4), consider r⊆s. For C∈c with anc(C)⊆r it holds anc(C)⊆s and hence mC(b)(r)=mC(b)(s), while for C∈c with anc(C)/⊆r one has mC(b)(r)=∞≥mC(b)(s). The single step property (A5) follows from

where in the inequality we have used the fact that for two sets in the sum above, one cannot be a subset of the other and hence, by definition of the hierarchy c, they are disjoint.

4.3 Exploiting logical relationships: Shaffer improvements

Logical relationships between hypothesis can be exploited to improve the power of the sequential rejection procedure. A first example of such an improvement for hierarchically ordered clusters was given in Meinshausen [24], while in Goeman and Finos [18] the improvement is applied to the inheritance procedure. Since those improvements are based on the considerations of Shaffer [25] they are called “Shaffer improvements”. For the high-dimensional setting a possible Shaffer improvement consists of multiplying the multiplicity adjustment mC(b)(r) with the Shaffer factor

where a set u⊆c is called congruent if, by the logical implications, it can be a complete set of false hypothesis (e.g. for a collection c of hierarchically ordered hypothesis u⊆c is congruent if for each C∈u it holds anc(C)⊆u and at least one offspring leaf node of C is in u).

Note that multiplication with the Shaffer factor never decreases the power of the method since by the monotonicity property (A4), sC(b)(r)≤1. Moreover, if r is congruent, then sC(b)(r)=1 and since the collection f of all false hypothesis is congruent, the Shaffer improvement doesn’t affect the validity of eq. (8). Finally, for r⊆s,

and hence the Shaffer improvement doesn’t affect the validity of eq. (7) neither.

We want to apply this Shaffer improvement to the inheritance procedure described in Section 4.2.1. Following the same reasoning as in Goeman and Finos ([18], Section 6), with the weights wC(b)=|Sˆ(b)∩C| we get the Shaffer factor

where si(C)=chpa(C)∖{C} denotes the siblings of C, l⊂c denotes the collection of leaf nodes, uC(b)=∑D∈si(C)wD(b) and vC(b)=minD∈si(C)wD(b). If c is a binary tree the Shaffer factor becomes

Unlike as for the inheritance procedure in (5), the Shaffer factor (6) for the procedure in (4) is always 1. Nevertheless, a possibility how to exploit logical relationships to improve the power of the procedure (4) for binary trees, which provides a Shaffer improvement very similar to the one above, is illustrated in Mandozzi and Bühlmann [15].

5.1 Implementation of the methods and considered scenarios

In this section we compare the performance of the four methods illustrated in Sections 4.1 and 4.2 and refined in Section 4.3, i.e. single variable method with Bonferroni multiplicity adjustment (2), hierarchical method with Bonferroni-based adjustment (4) along with Shaffer improvement as in Mandozzi and Bühlmann [15], single variable method with Bonferroni-Holm multiplicity adjustment (3) and hierarchical method with inheritance procedure (5) along with Shaffer improvement (6). In the following we refer to the first two methods as the “non-sequential methods” (strictly seen, the hierarchical method with Bonferroni-based adjustment is actually sequential, but there previous rejections are not used to improve subsequent multiplicity corrections) and the latter two methods as the “sequential methods”.

We consider the same implementation of the methods and the same scenarios (with exactly the same sample splits) as in Mandozzi and Bühlmann [15], although here we use only standard hierarchical clustering for the hierarchical methods. Concretely, the following choices have been made for implementation:

1. construction of the clusters with standard hierarchical clustering (using the R-function hclust) with distance between two covariables equal to 1 minus the absolute correlation between the covariables, and using complete linkage;

2. screening with the Lasso [19] with regularization parameter chosen by 10-fold cross-validation;

3. B=50 sample splits (for each scenario exactly the same splits as in Mandozzi and Bühlmann [15]);

4. for aggregation, the p-values PhC in Proposition 1 are computed over a grid of γ-values between γmin=0.05 and 1 with grid-steps of size 0.025;

5. nominal significance level α=5%.

1. 42 scenarios based on 7 designs;

2. for each design we consider 6 settings by varying the number of variables p in the model and the signal to noise ratio defined by

   SNR=(β0)TXTXβ0n−1σ−2,

   namely for p=200 we use SNR=4 and SNR=8, for p=500 we use SNR=8 and SNR=16 and for p=1000 we use SNR=16 and SNR=32;

3. 3 designs based on synthetic data (“equi correlation”, “high correlation within small blocks” and “high correlation within large blocks”) and 4 designs based on semi-real data (“Riboflavin with normal correlation”, “Breast with normal correlation”, “Riboflavin with high correlation”, “Breast with high correlation”);

4. sparsity s0=6 for the two “Riboflavin”-designs and s0=10 for the other five designs.

5.2 Familywise error rate control (FWER)

For each of the 42 scenarios described in Section 5.1 we consider exactly the same 100 independent simulation runs as in Mandozzi and Bühlmann ([15], section 4.2.2) by varying only the synthetic noise term ε and count the number where at least one false selection is made. According to Theorem 1, we expect this number to be at most 100α=5 (α=0.05). The results for the Bonferroni-based methods can be seen in Mandozzi and Bühlmann ([15], Table 1): FWER control holds for 40 of the 42 scenarios and in 37 scenarios there is no false selection at all.

The results for the methods with sequential rejection are very similar, the only differences being that for the “high correlation within small blocks”-design with p=500 and SNR=8 the number of runs with at least a false selection increases (compared to Bonferroni-type methods) from 7 to 9 for the single variable method, and from 7 to 13 for the hierarchical method, respectively; for the same design with p=1000 and SNR=16 the number of runs with at least a false selection increases from 5 to 6 for both the single variable and hierarchical method. For all other scenarios, inclusively the “high correlation within large blocks”-design with p=200 and SNR=4, where the non-sequential hierarchical method slightly failed to control FWER (6 runs with at least a false detection), the sequential methods exhibit the same FWER control as their non-sequential counterparts.

Summarizing, FWER holds for all four methods in 39 out of 42 scenarios and the designs where it doesn’t fully hold are “high correlation within small blocks” and “high correlation within large blocks”, which is not surprising since each active predictor is highly correlated with one or more inactive predictors from S0c and hence it is rather difficult for our screening method (the Lasso) to guarantee that Sˆ⊇S0.

5.3 Power

For measuring the power we consider four different aspects: the one-dimensional statistics defined in Mandozzi and Bühlmann ([15], section 4.2.1) as “Performance 1” and “Performance 2” (see below), the number of minimal true detections (MTDs, i.e., smallest significant groups of variables of any cardinality, containing at least one active variable, see below) and singleton true detections (STDs, i.e., MTDs with cardinality 1). Concretely, a cluster is said to be a MTD if it satisfies all of the following:

1. C is a significant cluster, e.g., has p-value <5% (“Detection”);

2. There is no significant sub-cluster D⊂C (“Minimal”);

3. C∈f, i.e., there is at least one active variable in C (“True”);

For each of the 42 scenarios outlined in Section 5.1, we consider exactly the same 100 independent simulation runs obtained in Mandozzi and Bühlmann ([15], section 4.2.3-4) by varying the synthetic noise term ε and the synthetic regression vector β0. We then calculate the average Performance 1, Performance 2, number of MTDs and number STDs, over the 100 simulation runs. The results are shown in Table 1 (for the single variable methods each MTD is an STD and by definition Performance 1 and 2 are the same).

Considering both low and high SNR, the methods with sequential rejection improve the considered power measures in comparison to the analogous method without sequential rejection in 207 out of 252 cases, the absolute improvement being at least 0.05 for MTDs and STDs, and at least 0.5% for Performance 1 and Performance 2 in 133 cases out of 252 cases. For better interpretation of these results: an absolute improvement of 0.05 MTDs (resp. STDs) basically means that in one out of 20 runs one more MTD (resp. STD) could be detected. Averaging over all scenarios, the improvement given by the sequential rejection procedures lies between 0.04 and 0.06 for MTDs and STDs, and between 0.5% and 0.7% for Performance 1 and Performance 2. The biggest gain with sequential rejection can be found in the “high correlation within small blocks”-design with p=200 and low SNR: it consists of 0.48 more STDs, 0.33 more MTDs and an absolute increase of 3.8% of Performance 1 and 3.3% of Performance 2, respectively. This basically means that in half of the runs the method with sequential rejection could find one STD more and in one third of the runs it could find one MTD more. Other particularly favorable scenarios for an improvement with sequential rejection are the “equi correlation”-design and the “breast normal correlation”-design, both with p=200 and low SNR and the “high correlation within small blocks”-design with high SNR and p=200, resp. p=500.

In general, the improvement given by the sequential rejection procedures decreases with increasing number p of covariables and is substantial only when the power of the method without sequential rejection is intermediate. These empirical findings are not surprising, since looking at how the methods are defined and in particular at the eqs (2)–(5), we conclude that an improvement with the sequential rejection methods is only possible if the related non-sequential method provides at least an STD (and gets more likely the more STDs are provided by the non-sequential method). Moreover, an improvement with sequential rejection is more likely to happen when the number |Sˆ| of screened variables is small.

For a better illustration of what kind of an improvement is possible using sequential rejection, we show in Figures 1 and 2 the dendrograms (in gray) for a paradigmatic simulation run of the “equi correlation” – and the “high correlation within small blocks”-design, respectively, both with p=200 and SNR=4.

Figure 1:

Dendrograms for a paradigmatic simulation run of the “equi correlation”-design with p=200 and SNR=4. The active variables are labeled in black and the truly detected non-zero variables along the hierarchy are depicted in black.

Figure 1 illustrates that sequential rejection allows the detection of a further singleton, increasing the number of STDs from 6 to 7 and the number of MTDs 8 to 9. In Figure 2 sequential rejection allows to detect a singleton that could previously only be detected together with another non-relevant variable in a cluster of cardinality 2, increasing the number of true STDs from 4 to 5 (while the number of MTDs remains to be 6).

Figure 2:

Dendrograms for a paradigmatic simulation run of the “high correlation within small blocks”-design with p=200 and SNR=4. The active variables are labeled in black and the truly detected non-zero variables along the hierarchy are depicted in black.

Finally, we have performed a simulation with the same scenarios (and the same sample splits) as in Mandozzi and Bühlmann ([15], section 4.3), i.e. “small blocks”-designs and “large blocks”-designs with 8 different correlations

The full results are shown in Table 2 in the Appendix. While the methods with sequential rejection control the FWER in exactly the same scenarios where it is also controlled by the non-sequential methods, they increase the average number of MTDs from 5.51 to 5.62 for the single variable method, and from 8.11 to 8.18 for the hierarchical method, and the number of STDs for the hierarchical method from 5.44 to 5.55, with improvements for a single scenario up to 0.48 MTDs and 0.56 STDs (averaged over 100 runs).

The empirical results can be summarized as follows. The methods with sequential rejection essentially controls the FWER in the same way as the non-sequential methods. Regarding power, sequential rejection allows for improvements, to a similar extent for the single variable and the hierarchical procedures. As already noted in Mandozzi and Bühlmann [15], for the non-sequential methods, the hierarchical methods have similar STDs as the single variable methods but allow for substantially more MTDs. Thus, our proposed hierarchical method with the inheritance procedure can be seen as the best of the considered methods.

5.4 Real data application: Motif Regression

We consider here a problem of motif regression [26] from computational biology. We apply the four methods described above, plus the two hierarchical methods (with and without sequential rejection) using the recently proposed canonical correlation clustering of Bühlmann et al. [14], to a real dataset with n=287 and p=195, used in Meinshausen ([24], section 4.3) and Mandozzi and Bühlmann ([15], section 4.4). The sequential rejection methods detects exactly the same significant structures as non-sequential methods, namely a single variable and a cluster containing 165 variables (the latter can be detected only with the hierarchical method with canonical correlation clustering). This can barely be considered as surprising, as with only one STD by the non-sequential methods, further improvements by the sequential methods are rather unlikely (see Section 5.3 for more explanation and empirical evidence).