Basis-Adaptive Selection Algorithm in dr-package

Sufficient dimension reduction (SDR) turns out to be a useful dimension reduction tool in high-dimensional regression analysis. Weisberg (2002) developed the dr-package to implement the four most popular SDR methods. However, the package does not provide any clear guidelines as to which method should be used given a data. Since the four methods may provide dramatically different dimension reduction results, the selection in the dr-package is problematic for statistical practitioners. In this paper, a basis-adaptive selection algorithm is developed in order to relieve this issue. The basic idea is to select an SDR method that provides the highest correlation between the basis estimates obtained by the four classical SDR methods. A real data example and numerical studies confirm the practical usefulness of the developed algorithm.


Introduction
Sufficient dimension reduction (SDR) in the regression of y ∈ R 1 |X ∈ R p = (x 1 , . . ., x p ) T replaces the original p-dimensional predictors X with its lower-dimensional linearly transformed predictors M T X without any loss of information on y|X, which is equivalently expressed: where stands for independence, M is a p × q matrix and q ≤ p.
A space spanned by the columns of M to satisfy (1) is called a dimension reduction subspace.Hereafter, S(M) denotes the column subspace of a p × q matrix M, and M T X is called a sufficient predictor.The intersection of all possible dimension reduction subspaces is called the central subspace S y|X , if it exists.The main goal of SDR is to infer S y|X , which is done through the estimations of its true structural dimension d and orthonormal basis matrix.
According to Cook (1998a), for a non-singular transformation of X such that Z = A T X, the following relation holds: S y|X = AS y|Z .Considering a standardized predictor Z = Σ −1/2 {X − E(X)}, we have that S y|X = Σ −1/2 S y|Z , where Σ = cov(X) and Σ −1/2 Σ −1/2 = Σ −1 .Typically, SDR methods estimate S y|Z first, then back-transform it to S y|X .Hereafter kernel matrices to restore S y|Z for each method will be denoted as M • , and it will be assumed that M • is exhaustively informative to S y|Z so that S(M • ) = S y|Z .With Z-scale predictors, the classical but most popularly used SDR methodologies include: (i) sliced inverse regression (SIR; Li, 1991): M SIR = cov{E(Z|y)}.If y is categorical, a sample version of E(Z|y) is straightforward.If y is many-valued or continuous, y is categorized by dividing its range into h slices.
(ii) sliced average variance estimation (SAVE; Cook and Weisberg, 1991): The sample version of cov(Z|y) is constructed in the same way as SIR.
(iii) principal Hessian directions (pHd; Li, 1992;Cook, 1998b): where β is the ordinary least squares on the regres- sion of y|Z.The sample version of M pHd is constructed by replacing the population quantities with their usual moment estimators.
(iv) covariance method (covk Yin and Cook, 2002): , where w = {y − E(y)}/ var(y).As can be easily seen, E(Zw j ) in M covk is the covariance between w j and Z as well as the ordinary least squares coefficient of w j |Z for j = 1, . . ., k.
These four SDR methods can be implemented in the dr-package.The package provides the basis estimates and dimension test results.A detailed review on the dr-package is given in Weisberg (2002).The following question is a critical issue in using dr: which one among the four methods has to be used?Even with the same data, one can obtain dramatically different dimension reduction results from the four methods in dr, yet there is no clear and up-to-date guideline on method selection.It is known that SIR and covk work better when a linear trend exists in regression, while SAVE and pHd are effective under a nonlinear trend, particularly in the case of a quadratic relationship.However, in practice, these may not be useful guidelines for choosing an SDR method, because it is not easy to check the The R Journal Vol.10/2, December 2018 ISSN 2073-4859 existence of non/linear trends in regression.Even in the case that a linear trend exists, it is still not clear which of SIR and covk should be used, as they yield different dimension test results.
The purpose of the paper is to develop a basis-adaptive selection algorithm for selecting an SDR method.The basic idea is to select the one that gives the highest correlation between the basis estimates obtained by all possible pairs of the four SDR methods.To measure the correlation, a trace correlation r (Hooper, 1959) will be used.The algorithm is data-driven and can be enhanced to other SDR methods, although it is highlighted with the use of the dr-package in this paper.
The organization of this article is as follows.We develop the idea of a basis-adaptive selection algorithm, and discuss its selection criteria.Next, we present a real data example and simulation studies.Finally, we summarize our work.

Development of algorithm
We begin this section with soil evaporation data in Section 6 of Yin and Cook (2002).The data contains 46 observations of daily soil evaporation, daily air and soil temperature curves, daily humidity curves and wind speed.For illustration purposes, we consider a regression analysis of daily soil evaporation given the integrated area of and the range of the daily air and soil temperatures.We will revisit the data in a later section.
Before developing selection criteria among the four SDR methods, there are some aspects that must be considered first.A selection based on the criteria should be data-driven, unless it is purposely The R Journal Vol.10/2, December 2018 ISSN 2073-4859 pre-selected.In addition, the criteria should be generally extended to the other SDR methods, although the four methods in dr are highlighted in this paper.
In order to have some idea of how to select a method, consider a simulated example of y|X = (x 1 , . . ., x 10 ) T = x 1 + ε.The column of η = (1, 0, . . ., 0) T spans S y|Z , and the sufficient predictor of x 1 is well estimated not only by SIR but also by cov2.Therefore, it is expected that the estimates by SIR and cov2 will be more highly correlated than those of any other pairs of the four methods.Assuming that X and ε are independently normally-distributed with µ = 0 and σ 2 = 0.1 2 , the averages of the absolute correlations between the pairs of the estimates from the four methods as well as between x 1 and the estimate from each method are computed through 500 iterations.
sir.cov <-sir.save<-sir.phd<-save.cov<-save.phd<-phd.cov<-NULL sir.eta <-save.eta<-phd.eta<-cov.eta<-NULL set.seed( 1 round(apply(cbind(sir.eta, cov.eta, save.eta, phd.eta), 2, mean), 3) sir.eta cov.eta save.etaphd.eta 0.941 0.939 0.215 0.267 The highest average of the six pairwise absolute correlations is highest between SIR and covk; in addition, either of SIR or covk estimates x 1 well.If the regression model is changed to y|X = x 2 1 + ε, the pair of SAVE and pHd yield the highest averages in the absolute correlations among the six pairs, and both SAVE and pHd are good SDR methods for the regression.
Let us think about this in a reverse manner.Suppose that the estimates of sufficient predictors by a pair of SAVE and pHd are more highly correlated than those of any other pairs of the SDR methods.Then, it would be not bad reasoning to believe that SAVE or pHd should be preferable to the data, although we do not know what the true model is.That is, the pair of the two methods that gives the highest correlation would be not a bad choice for estimating S y|X , although it is not guaranteed to be the best.We will formalize this idea.The similarity of two matrices with the same column rank is equivalent to the distance between their column subspaces.Finally, a correlation between of ηT a Z and ηT b Z can be alternatively measured by a trace correlation r tr (Hooper, 1959) between S( ηa ) and S( ηb ): The trace correlation r tr varies between 0 and 1, and higher values indicate that S( ηa ) and S( ηb ) are closer.If the two subspaces coincide, then we have r tr = 1.The trace correlation is expected to be maximized under the true dimension of d, which is usually unknown.For a smaller choice of d, the true basis are underestimated, so each method will not normally estimate the same parts of S y|X .This implies that r tr will be smaller compared to it under the true dimension.In contrast, in the case of a larger selection of d, the true basis are overestimated.Then, there is redundancy in the estimates, and the redundancy will be random for each method.This indicates that smaller correlations should be expected.This discussion is confirmed through the following simulation example.
2. Compute r tr between the basis estimates from each pair for m = 1, . . ., d max .For d = 1 or 2, the cov2 is commonly used, and the covd is employed for d ≥ 3.
3. Choose a pair of the methods to provide the maximum r tr in Step 2. The pair chosen in this step will be called the initial pair.
4. Remove the initial pair and all of the other pairs not containing one of the methods of the initial pair for all values of d.After completing this step, the surviving pairs must contain one method of the initial pair.
5. Search a second pair that has the highest r tr among all of the remaining pairs.This pair will be called the final pair.
6.The common method in the initial and final pairs is the representative SDR method to the data.
Steps 4-5 are required to select one of the methods of the initial pair in Step 2. This approach of selecting SDR methods through the proposed algorithm is called basis-adaptive selection (BAS).The BAS is data-driven and not necessarily limited in the four SDR methods of SIR, SAVE, pHd and covk, as one can extend it to other SDR methods.

The bas1 function
The bas1 function runs the BAS algorithm for SIR, SAVE, pHd and covk.The function requires dr and has the following arguments: bas1(formula, data, nsir=5, nsave=4, k=4, plot=TRUE) The arguments of nsir, nsave and k determine the numbers of slices for SIR and SAVE as well as the moment for covk, respectively.The default values are 5, 4 and 4, in order.If plot=TRUE, the function bas1 returns a scatter plot of the trace correlations against the various choices of d, up to min(nsir-1, k) for the six pairs of the four methods in dr.
The values returned by the bas1 function are selection, sir, save, phd and covk.Their descriptions are as follows: • selection: a selected method among SIR, SAVE, pHd and covk by the BAS algorithm • sir: the SIR application object with the number of slices equal to nsir • save: the SAVE application object with the number of slices equal to nsave  2. From Table 2, it can be seen that the selection results are different from the numbers of slices in SIR.For example, if the numbers of slices in SIR was 3, BAS recommended covk.For i = 4 and 5, BAS chose pHd and SIR, respectively.Thus, additional work needs to be done in order to choose between covk, pHd and SIR in the soil evaporation data.
For this purpose, we considered how reasonably the methods recommended in Table 2 estimated the true dimension of the central subspace.The dimension estimation results are already summarized in Table 1, and the nominal level 5% was used.First, the two methods of cov2 and cov3 were inspected.According to Table 1, cov2 and cov3 determine that d ≥ 2 and d ≥ 3, respectively.This indicates that the order of polynomial in covk should be bigger than or equal to 4 for further dimension determination.However, in the case k = 4, covk yields d ≥ 4, so the dimension reduction is meaningless because p = 4. Therefore, we conclude that covk would not be recommended one for this data.Next we consider pHd, which infers that d = 0 according to Table 1.Therefore, the pHd should also be ruled out for the final choice.In Table 1, the SIR with 5 slices determines that d = 2, which is the most reasonable among the three recommendations.Thus, one may continue the regression analysis with the two-dimensional sufficient predictors from the SIR results.
Each model was iterated 500 times with n = 100 and n = 200.Throughout all numerical studies, five slices were commonly used for SIR and SAVE, and four was the maximum polynomial in covk.
As a summary, the selection percentages of the methods by BAS are reported in Table 3.According to Table 3, the BAS selects SIR and covk 99.8% of the time for Model 1 and SAVE and pHd 97.8% of the time for Model 2 with n = 100, respectively.These results are consistent with the discussion given in the previous section regarding the development of the algorithm.In Model 1, covk is preferred to SIR with n = 100, but the two are almost equally selected with n = 200.For Model 2, the pHd is recommended more than SAVE with n = 100 and n = 200.In Model 3, the pHd and the covk are two dominant methods according to BAS, although the covk is selected more frequently with n = 100.For Model 4, the covk is recommended most frequently, but the selection percentages of SIR rapidly grow with n = 200.set.seed(5);sel1 <-sel2 <-sel3 <-sel4 <-sel12 <-sel22 <-sel32 <-sel42 <-NULL ## starting loop for (i in 1:500){ ## model construction for n=100

First
, we list all six possible pairs of the four methods in the dr-package: (1) (SIR, SAVE); (2) (SIR, pHd); (3) (SIR, covk); (4) (SAVE, pHd); (5) (SAVE, covk); (6) (pHd, covk).Let ηa and ηb be orthonormal basis estimates of S y|Z by any pair among the six under d = m.If d = 1, the correlation between ηT a Z and ηT b Z can be simply computed using the usual Pearson correlation coefficient.However, if d ≥ 2, the correlation between ηT a Z and ηT b Z is not straightforward.Now it should be noted that the correlation between ηT a Z and ηT b Z depends on the similarity between ηa and ηb regardless of the value The R Journal Vol.10/2, December 2018 ISSN 2073-4859 of d.

Table 2 :
Method selection results by BAS in soil evaporate data

Table 3 :
Percentages of the selection of a method among SIR, SAVE, pHd, and covk by BAS