LDR a Package for Likelihood-Based Suﬃcient Dimension Reduction

We introduce a software package running under Matlab that implements several recently proposed likelihood-based methods for suﬃcient dimension reduction. Current capabilities include estimation of reduced subspaces with a ﬁxed dimension d , as well as estimation of d by use of likelihood-ratio testing, permutation testing and information criteria. The methods are suitable for preprocessing data for both regression and classiﬁ-cation. Implementations of related estimators are also available. Although the software is more oriented to command-line operation, a graphical user interface is also provided for prototype computations.


Introduction
Since the introduction of sliced inverse regression (SIR; Li 1991) and sliced average variance estimation (SAVE; Cook and Weisberg 1991) there has been considerable interest in dimension reduction methods for the regression of a real response Y on a random vector X ∈ R p of predictors. A common goal of SIR, SAVE and many other dimension reduction methods is to estimate the central subspace S Y |X (Cook 1994(Cook , 1998, which is defined as the intersection of all subspaces S ⊆ R p with the property that Y X|P S X. Informally, these methods pursue the estimation of the fewest linear combinations of the predictors that contain all the regression information on the response. SIR uses a sample version of the first conditional moment E(X|Y ) to construct such an estimator, while SAVE uses sample first and second E(XX T |Y ) conditional moments. Although SIR and SAVE have found wide-spread use in applications, both have well known limitations. In particular, the subspace S sir estimated by SIR is typically a proper subset of S Y |X when the response surface is symmetric about the origin. SAVE was developed in response to this limitation and provides exhaustive estimation of S Y |X under mild conditions (Li and Wang 2007;Shao, Cook and Weisberg 2007), but its ability to detect linear trends is generally inferior to SIR's. For these reasons, SIR and SAVE have been used as complementary methods, with satisfactory results often being obtained by informally combining their estimated directions (see, for example, Bura and Pfeiffer 2003;Cook 2003;Ye and Weiss 2003;Li and Li 2004;Zhu, Othaki and Li, 2005;Pardoe, Yin and Cook 2007). Another method called directional regression (DR) was recently proposed by Li and Wang (2007) and it was shown to provide exhaustive estimation of S Y |X under mild conditions. In a series of recent articles, Cook (2007) and Forzani (2008, 2009a-b) took a substantial step forward in the development of dimension reduction methods based on the first two conditional moments. The idea behind this new methodology is to estimate S Y |X via maximum likelihood. Likelihood-based methods provide exhaustive estimation of S Y |X under the same conditions as SIR, DR and SAVE. Unlike those methods, however, a likelihood-based objective function is employed to acquire the reduced dimensions, giving these methods √ n consistency and asymptotical efficiency, which is not a claimed attribute for any previous method. The dimension d of S Y |X can be estimated using likelihood ratio testing or an information criterion like AIC or BIC, and conditional independence hypotheses involving the predictors can be tested straightforwardly. Furthermore, it was demonstrated both theoretically and with simulations that these methods have good robustness properties, even under substantial deviations of the error distribution from its nominal specification.
The goal of this article is to introduce software tools for these recently proposed likelihoodbased methods for sufficient dimension reduction. We have called the package LDR (standing for Likelihood-based Dimension Reduction) and it runs under Matlab, which allows use of computational tools developed outside of the statistics community. In particular, likelihood maximization under some of the models discussed below requires optimization on a Grassman manifold G (d,p) , which is the set of all d dimensional subspaces of R p (Edelman, Arias and Smith 1999). Optimization on Grassmann manifolds is fairly common in some areas such as computer science and signal processing, but is relatively rare in statistics. The software is oriented to command-line operation and it is intended to provide an easy-to-use interface similar to Weisberg's DR package for R (http://cran.us.r-project.org/web/packages/dr). Fixed-dimension as well as dimension selection tasks are supported, both for continuous and for discrete responses. Functions to compute other estimators such as SIR, SAVE and DR are also provided, although methods for dimension estimation are not yet available for these methods.
The paper is organized as follows. In Section 2 we briefly review some theory concerning likelihood-based sufficient dimension reduction. In Section 3, we describe the main features of the software and then give some examples of how to use it in Section 4. Finally, we briefly describe how to use the graphical user interface in Section 4.3. The package and the complete software documentation is available at http://liliana.forzani.googlepages.com/ldrpackage.

Likelihood-based sufficient dimension reduction
The likelihood-based approach to dimension reduction uses model-based inverse regression of X on Y to gain reductive information for the forward regression of Y on X. We assume that the data consist of n independent observations on (X, Y ), with X ∈ R p , Y ∈ R, n > p and that X|Y ∼ N (µ y , ∆ y ). Then we can write where ε ∼ N (0, ∆ y ), ∆ y > 0. The goal is to find the maximum likelihood estimator for the central subspace, S Y |X . Let d = dim(S Y |X ) and let α denote a p × d semi-orthogonal matrix whose columns form a basis for S Y |X , so Y |X ∼ Y |α T X. Likelihood-based methods are designed to estimate S Y |X = span(α) under different structures for the mean µ y and the covariance ∆ y functions in model (1). These methods are summarized in the sections that follow. In preparation, let M = span{µ y − µ|y ∈ S Y } ⊆ R p , where S Y denotes the sample space of Y , and letâ stand for the maximum likelihood estimator of a parameter a.

Principal Component, PC.
In this model, which was introduced by Cook (2007), it is assumed that the errors are isotonic, ∆ y = σ 2 I p . As a consequence of this structure, S Y |X = M and therefore we have the representation µ y = µ + αν y , where ν y ∈ R d is unknown for all y ∈ S y . Under this setting, the maximum likelihood estimator of S Y |X is the span of the first d eigenvectors of the sample covariance matrixΣ of X.

Isotonic Principal Fitted Components, IPFC.
As in the PC model, it is assumed that the errors are isotonic. Consequently S Y |X = M and we can again represent the conditional means as µ y = µ + αν y . However, following Cook (2007), the coordinate vectors ν y are now modeled as where f y ∈ R r is a known vector-valued function of y with linearly independent elements and β ∈ R d×r , d ≤ min(r, p), is an unrestricted rank d matrix. For this model the maximum likelihood estimator of S Y |X is the first d eigenvectors ofΣ fit , whereΣ fit is the sample covariance matrix of the fitted vectors from the multivariate linear regression of X y on f y , including an intercept. When Y is discrete or slicing is used to approximate the coordinates of ν y with step functions,Σ whereX y the average predictor vector in slice y,X the overall average, and h is the number of slices.

Structured Principal Fitted Components, SPFC.
This model was introduced by Cook and Forzani (2009-a). The coordinate vectors ν y are again modeled as in (2), but now a linear structure is used to model ∆ = m i=1 δ i G i , where m ≤ p(p + 1)/2, G 1 , . . . , G m are known real symmetric p × p linearly independent matrices and the elements of δ = (δ 1 , . . . , δ m ) T are functionally independent. It is required also that ∆ −1 have the same linear structure as ∆: where e i ∈ R p contains a 1 in the i-th position and zeros elsewhere. This basic structure can be modified straightforwardly to allow for a diagonal ∆ with sets of equal diagonal elements, and for a non-diagonal ∆ with equal off-diagonal entries and equal diagonal entries. In the latter case, there are only two matrices G 1 = I p and G 2 = ee T , where e ∈ R p has all elements equal to 1.
In this case the central subspace S Y |X = ∆ −1 M and then µ y = µ + ∆αβf y . The maximum likelihood estimator for S Y |X is the span of the first d eigenvectors of ∆ −1Σ where ∆ indicates the MLE of ∆. The MLE ∆ is relatively complicated and an iterative algorithm is evidently required for its computation 2.4. Extended Principal Fitted Components, EPFC.
In this version the coordinate vectors ν y are again modeled like (2) so that µ y = µ + αβf y , but now ∆ is modeled as ∆ = αΩα T + α 0 Ω 0 α T 0 , where Ω ∈ R d×d and Ω 0 ∈ R p−d,p−d are positive definite matrices and still S Y |X = span(α). Then the maximum likelihood estimator for S Y |X maximizes over S(α) ∈ G (d,p) the log likelihood function whereΣ res =Σ −Σ fit withΣ fit defined in the IPFC model. See Cook (2007) for further discussion.

Principal Fitted Components, PFC.
In this model, which was studied by Cook and Forzani (2009-a) the coordinate vectors ν y are again modeled as in (2), but ∆ is now an unstructured positive definite matrix that is functionally independent of y. The central subspace is again S Y |X = ∆ −1 M and the MLE of res , withΣ res defined as in the EPFC model.

Covariance Reduction, CORE.
Covariance reduction addresses a rather different dimension reduction issue. Consider the problem of characterizing the covariance matrices ∆ y , y = 1, . . . , h, of a random vector X observed in each of h normal populations. There is no interest in the population means µ y . The methodology reviewed here, which was proposed by Cook and Forzani (2008), is based on reducing the sample covariance matrices to an informational core that is sufficient to characterize the variance heterogeniety across the h populations. It may be a useful alternative to the spectral modeling in many applications.
For this problem the central subspace, S Y |X = span(α) is redefined to be the smallest subspace satisfying ∆ y = ∆ + P T α(∆y) (∆ y − ∆)P α(∆y) , where P α(∆y) = α(α T ∆ y α) −1 α T ∆ y the projection onto S Y |X using the inner product defined by ∆ y . We assume that the data consist of n y independent observations of X y , y = 1, . . . , h. Then the maximum likelihood estimator of S Y |X is the subspace span(α) that maximizes over G (d,p) the log likelihood function where c is a constant depending only on p, n y and∆ y , y = 1, . . . , h denote the sample covariance matrix from population y computed with divisor n y , and∆ = 1 n h y=1 n y∆y .

Likelihood Acquired Directions, LAD.
In this model, which was proposed by Cook and Forzani (2009-b), we consider both a general conditional covariance ∆ y and a general mean µ y . It requires a discrete response Y . When the response is continuous it is typical to follow Li (1991) and replace it with a discrete version constructed by partitioning its range into h slices. The central subspace S Y |X = span(α) is the smallest subspace that satisfies the conditions (i) ∆ y = ∆ + P T α(∆y) (∆ y − ∆)P α(∆y) and (ii) span(α) ⊂ ∆ −1 span(µ y − µ). Condition (i) is the same as that encountered in the CORE model, while condition (ii) incorporates the conditional means µ y . Assume as before that the data consist of n y independent observations on X y , y = 1, . . . , h. Then the MLE for S Y |X maximizes over span(α) ∈ G (d,p) the log likelihood function where∆ y is as defined in the CORE model.

Envelope Models for Multivariate Linear Regression, MLM
The dimension reduction methodology described in this section is different from the previous methods because it involves a multivariate response Y ∈ R r and a multivariate predictor X ∈ R p , assumed to be non-stochastic, that are related through the standard normal multivariate linear model: where β 0 ∈ R r , β ∈ R r×p and ε ∼ N r (0, Σ). The number of responses r is often large in modern applications of this model. In such situations the response vector may contain irrelevant linear combinations whose conditional distribution does not change with X. It may also contain redundant information characterized by Σ having relatively large eigenvectors. Developed recently by Cook, Li and Chiaromonte (2009), envelope models allow for the possibility that such irrelevant or redundant information may be present in Y, and can yield maximum likelihood estimators of β with substantially smaller variation that the usual maximum likelihood estimators.
The redundant and irrelevant information in (5) can be characterized by using the reducing subspaces of Σ. Let (Γ, Γ 0 ) ∈ R r×r be an orthogonal matrix, where the columns of Γ ∈ R r×d form a basis for the smallest reducing subspace of Σ that contains span(β). Then the parameters in model (5) have the structure β = Γη, and Σ = ΓΩΓ T + Γ 0 Ω 0 Γ T 0 , where η ∈ R d×r is an unconstrained coordinate matrix, and Ω ∈ R d×d and Ω 0 ∈ R r−d×r−d are positive definite matrices. The theory developed by Cook, Li and Chiaromonte (2009) shows that substantial gains in efficiency are possible if the linear combinations Γ T 0 Y contain both the irrelevant and redundant information in Y.
Fitting with this parameter structure requires estimation of β 0 , span(Γ), η, Ω and Ω 0 . The partially maximized log likelihood for estimating span(Γ) is similar to the objective function for extended principal fitted components described in Section 2.4.

Software overview
The software is organized to allow users to run different types of processing by calling just a few interface functions. Several auxiliary tools are then internally called to perform computations according to the specified options on those interface functions. As the main example, a single function called ldr manages all methods for maximum likelihood estimation of the central subspace. Arguments in this function specify the model to use, the response type -continuous or discrete -and the protocol for the dimension d. This protocol allows d to be specified or estimated using an information criterion, a likelihood-ratio test or in some case a permutation test. Separate functions mlm are used for the envelope methodology discussed in Section 2.8. In Section 4 we describe how to use the software in detail.
The package is designed to run mainly from the command-line, but we have also built a simple graphical user interface to make usage more intuitive and analysis more interactive. Although the graphical interface does not currently provide all the features available from the command-based interface, it is powerfull enough to run prototype computations using standard values for some parameters of the models. A brief tutorial for using this interface is given in Section 4.3.
Our implementation emphasizes easy-reading code over performance. It is mainly based upon built-in functions in Matlab and has been tested for compatibility starting up from version 6.5 Release 13. Nevertheless, some functions require the Statistics Toolbox too. The LDR package relies on Lippert's sg-min toolkit for gradient-based optimization over Stiefel-Grassmann manifolds (Lippert and Edelman 2000). This toolkit is freely available at http://www-math.mit.edu/∼lippert/sgmin.html. Some minor modifications of sg-min were necessary for this application and this modified version of the sg-min is available with our package. The modifications are described in the documentation that comes with the code.

Available methods
The LDR package offers an implementation of the models described in Section 2. SIR, SAVE and DR are incorporated as separate Matlab functions with the same names. SIR and SAVE are available in Weisberg's DR package, and DR is available from its authors (Li and Wang 2007). We have included basic implementations here in order to provide a more self-contained tool for dimension reduction. We have also used a separate Matlab function PC for the PC model, since this model gives the same solution as the usual principal components and we do not have a novel way to estimate d.
For the rest of the models, which are the primary focus of the LDR package, we have incorporated methods to estimate the central subspace, including estimation of its dimension d. Available procedures for choosing d include information criteria such as AIC and BIC (Burnham and Anderson 2002), likelihood-ratio testing, and permutation testing (Cook and Yin 2001). Table 1 shows the available testing methods for each model. All of these functions run both for continuous and discrete responses and share the common interface function ldr.

Managing data, results and plots
Interface functions like ldr have workspace variables as arguments and do not allow managing data files directly. Despite several procedures for reading data being already available in Matlab, they usually do not provide support for data files with headers. We have included functions loadDATA and getDATA to allow for this. The former returns a data matrix, the text in the header and the labels of the variables if all of them exist. On the other hand, getDATA returns the response vector and the predictors directly, but does not return the header. In both of these procedures, data must be numeric and organized with rows as cases and columns as variables, with columns separated by white space. If another delimiter is used to separate columns, you can specify it as an optional input.
The header can be just a line of labels, one for each column and with the same delimiter, or it can have more text describing the data. In either case, the header and variable names should be separated by a line feed. For illustration, assume that datafile.txt consists of the following semi-colon delimited data: This is the header for seven variables var1 var2 var3 var4 var5 var6 var7 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; The command loadDATA can be used to read this data file: [data,header,labels] = loadDATA('datafile.txt',';'); The semi-colon at the end of the last line is the Matlab command to supress printing. To read the same file with white space used as the delimiter use the command [data,header,labels] = loadDATA('datafile.txt'); If the response Y is in the first column of matrix data and all other columns are predictors X, then we can create these data arrays in our workspace with the commands Y = data(:,1); X = data(:,2:end); We can get the same result by typing:  If we want just the first five columns as predictors we should type: [Y,X] = getDATA('datafile.txt',';',7,1:5); The graphical user interface also allows for using any column in a data array as the response vector and any number of the remaining columns as the predictors. Nevertheless, for full compatibility with the software, the response should be the first column of the data matrix with the predictors occupying the remaining columns.
While saving results in text files is a straightforward operation in Matlab, we have included a function saveastxt to make it even easier for those who are not familiar with the language.
As it is usual in Matlab, help for this or any other function in the LDR package is available with a command of the form help function-name. In addition, we have also added some tools for plotting results to supplement those in Matlab. The function plotDR allows for an easier labeling of results and automatical selection of the most suitable plot according to the dimension of the reduced subspace and the type of response. A brief overview of these capabilities is given in Section 4. More detailed description is available through the software documentation.
4. Using the software 4.1. Methods for estimating the central subspace Except for envelope models, a single function called ldr provides a unified interface for all dimension reduction methods based on maximum likelihood estimation of the central subspace.
Other methods such as SIR, SAVE, DR and PC are called from separate interface functions. Their usage is very similar and is discussed briefly in this section.
To start using the software, set the current directory to the package's location in your disk and type setpaths. This command adds all the directories in the package to the Matlab path, so that all functions there become available for computation.
Five arguments are mandatory when calling the ldr function and they must be provided in the order that follows: i) the response vector Y ; ii) the matrix of predictors X; iii) an identifier for the model to be used, such as 'PFC' or 'epfc'; iv) the type of response, which may be 'cont' or 'disc' depending on whether Y is continuous or discrete; and v) the dimension d of the central subspace or an identifier such as 'BIC' for an estimation method along with any necessary argument. Acronyms previously introduced in Section 2 have been used to identify the available models. Identifiers for estimation methods are 'AIC' for Akaike's information criterion, 'BIC' for Bayes information criterion, 'LRT' for likelihood-ratio test and 'PERM' for permutation test. Matching is case-insensitive, so you can type either 'LRT' or 'lrt'. The identifiers 'LRT' and 'PERM' should be followed by a test level, typically .01 or .05.
Outputs are the same for all models and they are: i) the projection of the predictors onto the estimated central subspace; ii) a matrix whose columns are the estimated basis vectors for the central subspace; iii) the maximum value of the log likelihood; and the iv) the estimated dimension in case a criterion for d was used.
Depending on the model, additional and optional arguments can be given. As an example, to slice a continuous response into ten slices for processing data under the LAD model, an additional string 'nslices' and the value 10 should be added to the arguments. In this case, a call to the function ldr should look like: [WX,W,L,d] = ldr(Y,X,'LAD','cont','aic','nslices',10); where WX is the projection of the predictors onto the estimated central subspace, W is a matrix whose columns are the estimated basis vectors for the central subspace, L the maximum value of the log likelihood and d is the dimension of the central subspace estimated by using Akaike's information criterion.
Methods such as EPFC, LAD and CORE require a starting basis for the central subspace prior to iteration. We determine the starting basis internally by searching over several estimators like SIR, SAVE, DR, PLS and PC and selecting the one that gives the largest value of the likelihood. Nevertheless, in some cases an external starting basis may be useful. A starting basis can be supplied by adding the text string 'initval' followed by a matrix whose column form a starting basis. Iterative estimation of the central subspace is carried out using the sg-min package. All the available parameters in this toolkit have been retained as optional inputs in LDR. In addition, you can set the maximum number of iterations to be done. By default, computation will stop after 10000 iterations if convergence has not been achieved. Another optional argument available for all methods is '-v', which enables a verbose mode to get progress information of the running process. This input should always be given as the last one. The previous call to LAD in verbose mode with an initial basis B ∈ R p×d for the central subspace and with no more than 1000 iterations would look like 'cont','aic','nslices',10,'initval',B,'maxiter',1000,'-v').
A list of the remaining optional arguments and their identifiers is given in Table 2.
To further illustrate how to use the package, consider a script to study the identification of hand-written digits {0, 1, . . . , 9}. The 44 subjects were asked to write 250 random digits. Each digit yields a 16-dimensional feature vector, consisting of 8 pairs of randomly sampled twodimensional locations on the digit. The 44 subjects were divided into two groups of size 30 and 14, in which the first formed the training set with sample size 7, 494 and the second formed the test set with sample size 3, 498. The data set is available from the UCI machine-learning repository at ftp://ftp.ics.uci.edu/pub/machine-learning-databases/pendigits. We focus on dimension reduction of the 16-dimensional feature vectors for the training set, which serves as a preparatory step for developing an efficient classifier. For clarity, we first consider the digits 0, 6 and 9. The reduced data set comprises 2,219 cases. Suppose data are stored in the text file digits.txt where the first column takes values 1, 2 or 3 indicating if the observation corresponds to a 0, 6 or 9 respectively. Remembering the file convention discussed in Section 3.2, we first load the response vector and the predictors by typing: setpaths; % adds all directories in the LDR package to Matlabs's path. data = load('digits.txt'); % reads data from text file. Y = data(:,1); % assigns the first column of data to the response. X = data(:,2:end); % assigns the remaining columns to the predictors.
Now we look for a subspace of dimension d = 2 using the LAD model: [WX,W] = ldr(Y,X,'LAD','disc',2); Here we used default options for optimization. To plot results, we can call the provided function plotDR as: plotDR(Y,WX,'disc','LAD'); Here, arguments 'disc' and 'LAD' allows for labeling data according to the response and for setting axes labels. The function plotDR behaves like standard Matlab plotting functions. For instance, a second call to plotDR will overwrite the first plot, unless the Matlab command for a new figure has been issued. To carry out a similar analysis but using the CORE model, we could just replace 'LAD' with 'CORE' in all the statements above. When using SIR, SAVE, DR or PC the dimension d should be given, as no testing method for d is currently available for them in the LDR package. If we choose d = 2, corresponding commands should look as: Processing data with continuous responses is completely analogous, but it often allows for additional arguments. For illustration, let's review a script to study the data from Figure 5 in Cook and Forzani (2009-b). Suppose you have already loaded your continuous data into workspace variables Y and X (the complete script and data are available in the directory data). To process under the LAD model, suppose 5 slices are suitable and d = 1. Then, use a statement such as ldr(Y,X,'LAD','cont',1,'nslices',5); If we were looking for a similar analysis but wanted instead to estimate d by using likelihoodratio testing with level 5% we could type: (Y,X,'LAD','cont','lrt','nslices',5,'alpha',0.05); Similarly, if we were considering estimating d with a permutation test using 500 permutations with level of 5 %, we would type: [WX,W,f,d] = ldr (Y,X,'LAD','cont','perm','nslices',5,'alpha',0.05,'npermute',500); These arguments are inappropriate for models like IPFC, SPFC, EPFC and PFC. For them, instead, the f y 's should be given as an n × r matrix, say FY, with rows f T y . In this package we provide an auxiliary function get_fy which returns f y as a polynomial of order r, where the order of polynomials r to be used for f y should be given. Calling for the PFC model, for example, should look like: Further examples of how to use the software are provided in the papers folder within LDR. They reproduce results and figures published when introducing the methods discussed in Section 2 and they can help in getting used to the main features available in the package while testing the methods.

Envelope models
The basic command for computing an estimated basis G of span(Γ) in the envelope model (2.8) is [GX,G,L,dhat] = mlm_fit(Y,X,dim), where dim stands for the dimension of the envelope. The essential output consists of the estimate G of Γ, the value L of the log likelihood at the MLE's and the estimated dimension dhat of the envelope. GX is an exploratory quantity and not used routinely in this application. The following options are avaiable for dim: If dim is set to an interger d, eg., mlm_fit(Y,X,2), then dhat = d and the envelope model is fitted with the specified dimension. If dim is set to 'aic' or 'bic' then the indicated information criterion AIC or BIC is used to estmimate d. These commands may take a long time, depending on the size of the regression. If dim is set to the pair 'lrt',.05, eg,. mlm_fit(Y,X,'lrt',.05), then likelihood ratio testing at level .05 is used to estimate d.
The following commands are available following the initial fit that produces G.
• [betaem,eta,Omega,Omega0,S1,S2] = mlm_empars(Y,X,G) returns the estimated parameters from the fit of the envelope model. betaem = G x eta is the estimated coefficient matrix Omega and Omega0 are the estimates of Ω and Ω 0 , S1 is the estimate of ΓΩΓ T and S2 is the estimate of Γ 0 Ω 0 Γ T 0 .
• mlm_emses(Y,X,G) returns the matrix of standard errors of the elements betaem from the fit of the envelope model.
• mlm_seratios(Y,X,G) returns the matrix of ratios of standard errors of the full model and envelope model estimates of β. The ratios are (full model se's)/(envelope model se's).
Other commands are documented in the folder mlm that comes with the Matlab code.

Using the graphical user interface
The graphical user interface (GUI) is expected to provide an easy-to-use tool to perform small prototype tasks using the functions available in the package. The GUI starts by typing demo in the command window, provided current directory is the LDR folder. Once it has been loaded, data analysis using the GUI typically follows a sequence of steps: 1. Load data by clicking on the LOAD DATA button and then follow the dialog box until finding the desired data file.
2. Specify the type of the response according to the data by choosing the right option in the popup menu. For continuous responses, you should also give the number of slices for the slicing procedure if you plan to apply methods such as LAD, CORE, SIR, SAVE or DR. Similarly, you can use this editable text box to type specific input values for each method. As another example, if you plan to apply models such us PFC, IPFC, SPFC or EPFC, you should type here the order of the polynomials you want to use in order to estimateΣ fit . Only polynomial bases can be specified when using the GUI.
3. Set the dimension reduction method to apply. Available methods are listed in a popup menu.
4. Select the dimension of the reduced subspace or choose a method to test for it. Currently, the GUI allows only d = 1, 2, 3. Testing criteria include all those listed in Table 1. Some testing criteria, in particular permutation tests, take a long computation time and can freeze the GUI for awhile.
5. Start the computation by clicking the RUN button. You will then be asked to select the response vector from the data file. Notice you can select only one vector as the response. By default, data corresponding to the first column in the data file will be highlighted.
Clicking on the RUN button again to select the predictors. A default set is highlighted, but you can change selection if you want to use a subset of them for computations. A third click on the button sets the selected predictors and starts processing data.
By default, processing data with the GUI only returns results through plots. Nevertheless, you can check the Show results option to allow the program to print results on the command window. Furthermore, you can save results as text files by clicking on the SAVE RESULTS button. When you do so, you are first asked to select what results you want to save. All analysis performed so far are listed and you can choose as many as you like. For each selection, you are allowed to save two files. The first one contains the response plus the transformed predictors, while the second one contains only the generating vectors for the reduced subspace. You can skip saving either of them by clicking on the CANCEL button in the respective save dialog.
After some data is loaded, the popups below the figure axes get filled with identifiers for each column in the data file. You can get scatter plots of pairs of variables by selecting one of them in the popup corresponding to the vertical axis and the other one in the popup menu related to the horizontal axis. For a sequential display of scatter plots for some given variable, you can instead set the related identifier in one of these popups and then click on the SCATTER PLOTS button. You can then move through the scatter plots by clicking the UP and DOWN buttons placed below the SCATTER PLOTS button.

Conclusion
This paper introduces new software for sufficient dimension reduction based on maximum likelihood estimation of the central subspace. Matlab implementation of these methods allows easy integration of complex optimization tools and provides a flexible environment for graphics capabilities and further model development. Scripts are given to reproduce all published results concerning the methods discussed above, in agreement with the reproducible research philosophy. We plan to continue adding methods to the package and to finish code verification for full compatibility with the OCTAVE programming language in order to provide tools suitable to run under GNU software.

Aknowledgements
We are indebted to Eduardo Tabacman for his patient testing of the code and for bringing us some tips about the sg-min package. We also thank Ross Lippert for allowing us to distribute a modified version of his code with our package. This work was supported in part by National Science Foundation grant DMS-0704098 awarded to RDC.