knitr::opts_chunk$set(echo=TRUE,eval=TRUE,comment="#>")
Table of content
This package offers an implementation of the $SIR$ (Sliced Inverse Regression) method, along with a thresholded version of $SIR$ that allows variable selection. The purpose of this vignettes is to explain the function implemented in the package, and how they are used. The first section is dedicated to the theoretical explaination of $SIR$, and its thresholded version. Reading this section is therefore optional, but can be useful to better understand how the package functions work.
Consider the following semi-parametric model: $$Y = g(\beta'X,\epsilon)$$ where:
It is common to present the model as $Y\perp X\mid \beta'X$ which means that $Y$ is independent of $X$ conditional on $\beta'X$. Thus, it is possible to replace $X\in\mathbb{R}^p$ with the index $\beta'X\in\mathbb{R}$ without loss of information in the regression of $Y$ on $X$. Thus, we do dimension reduction of the explanatory part $X$ from $p$-dimension to $1$-dimension, without specifying either the relationship between the index $\beta'X$ and $Y$ or the distribution of the error term $\epsilon$.
In the framework of this model, $g$ being unknown, we cannot identify $\beta$, but we can identify its direction. We then define the $EDR$ (Effective Dimension Reduction) space $E$ as the linear subspace of $\mathbb{R}^p$ generated by $\beta$. Any vector belonging to $E$ is then an $EDR$ direction.
The $SIR$ method is used to obtain the estimate of an $EDR$ direction. In theory, the $SIR$ method requires a condition on the distribution of the explanatory variable $X$ called the linearity condition: $\forall b \in\mathbb{R}^p, \mathbb{E}[b'X\mid\beta'X]$ is linear in $X'\beta$. In practice, it is impossible to verify this assumption because $\beta$ is unknown and unidentifiable. However, this assumption is valid when the distribution of $X$ is elliptical, which is the case for the multidimensional normal distribution.
Let's now consider a monotonic transformation $T$ (which in the case of the SIR method will correspond to a slicing function). Under the previous model and the linearity condition, (Duan & Li, 1991) showed that the centered inverse regression curve verifies : $$\mathbb{E}[X\mid T(Y)] - \mu \in \mbox{Span}(\Sigma\beta)$$ Thus, the subspace generated by this curve, ${\mathbb{E}[X \mid T(Y)] - \mathbb{E}[X] : Y \in \mathcal{Y}}$ (where $\mathcal{Y}$ is the support for the response variable $Y$), belongs to the $EDR$ space. Therefore, the centered inverse regression curve can clearly be used to find the $EDR$ space.
A direct consequence of this result, is that the variance covariance matrix of this curve, $\Gamma = \mathbb{V}(\mathbb{E}[X \mid T(Y)])$, is degenerate in any direction $\Sigma$-orthogonal to $\beta$. Therefore, the eigenvector associated with the largest nonzero eigenvalue of the matrix of interest $\Sigma^{-1}\Gamma$ is an $EDR$ direction, so it generates the $EDR$ space $E$.
In the slicing step of the $SIR$ method, the support of $Y$ is sliced into $H$ slices ${s_1,\dots,s_H}$. With a function $T$ corresponding to this slicing, the matrix $\Gamma$ can then be written as $$\Gamma=\sum_{h=1}^Hp_h(m_h-\mu)(m_h-\mu)'$$ where $p_h=P(Y\in s_h)$ is the theoretical proportion of $Y_i$ falling in the $s_h$ slice and $m_h=\mathbb{E}[X\mid Y\in s_h]$ is the theoretical mean of $X_i$ associated with $Y_i$ falling in the $s_h$ slice. This matrix $\Gamma$ is therefore easy to estimate.
\
Let us consider a sample ${(X_i,Y_i), i=1,\dots,n}$ generated from the previous model. In practice, it is necessary to replace the theoretical moments $\mu, \Sigma, p_h$ and $m_h$ by the empirical moments. Then, an estimator of a basis of the $EDR$ space $E$ is obtained via the eigenvector $\hat{b}$ associated to the largest eigenvalue of the estimator $\widehat{\Sigma}_n^{-1}\widehat{\Gamma}_n$ of $\Sigma^{-1}\Gamma$ where :
$\widehat{\Sigma}n= \frac{1}{n} \sum{i=1}^n (X_i - \bar{X}_n)(X_i -\bar{X}_n)'$ is the empirical variance covariance matrix of $X$
$\widehat{\Gamma}n= \sum{h=1}^H\hat p_{h,n} (\hat m_{h,n} - \bar{X}n)(\hat m{h,n} -\bar{X}_n)'$ is the empirical variance covariance matrix of the empirical means of $X$ per slices $s_h$
with :
$\bar{X}n= \frac{1}{n}\sum{i=1}^nX_i,$ the empirical mean of the sample
$\hat p_{h,n} = \frac{1}{n}\sum_{i=1}^{n} \mathbb{I} {[Y_i\in {s_h} ]} = \frac{\hat{n}{h,n}}{n}$ the empirical proportion of $Y_i$ falling in the slice $s_h$ (where $n{h,n}$ is the number of $Y$ falling in the slice $s_h$)
$\hat m_{h,n}=\frac{1}{\hat{n}{h,n} }\sum{i\in {s_h}}X_i$ the empirical mean of $X_i$ such that $Y_i$ falls into the slice $h$
It is possible to perform variable selection in SIR regression. The method presented here is a computational method using soft or hard thresholding of the matrix of interest. Let us define the soft thresholding ($S$ as in soft) and hard thresholding ($H$ as in hard) operators that apply to any matrix $M$ of dimension $p\times p$. For any real $\lambda\geq 0$ and $(i,j)\in {1,\dots,p}\times {1,\dots,p}$ :
$$ S_\lambda(M){i,j}= sign(M{i,j}) \times \begin{cases} \lvert M_{i,j} \rvert -\lambda & \mbox{if } \lvert M_{i,j} \rvert -\lambda >0, \ 0 & \text{otherwise.} \end{cases} $$
$$ H_\lambda(M){i,j}= \begin{cases} M{i,j} & \mbox{if } \lvert M_{i,j} \rvert -\lambda >0, \ 0 & \text{otherwise.} \end{cases} $$
with $sign$ the function that returns the sign of the considered element ($1$ if positive, $0$ if null and $-1$ if negative).
We then introduce two new methods:
For each approach, an eigenvector associated to the largest eigenvalue of $S_\lambda(\widehat{\Sigma}n^{-1}\widehat{\Gamma}_n)$ or $H\lambda(\widehat{\Sigma}_n^{-1}\widehat{\Gamma}_n)$ is computed, for a given value of $\lambda$.
The choice of the thresholding hyper-parameter $\lambda$ must provide a balance between
Then, once the $p^{\star}$ useful variables have been identified, a new $SIR$ model can be estimated on the basis of these $p^{\star}$ selected variables. Thus, the estimation of $\beta$ will be easier because the estimation is done in a space of reduced dimension.
SIR
functionApply a single-index $SIR$ on $(X,Y)$ with $H$ slices. This function allows to obtain an estimate of a basis of the $EDR$ (Effective Dimension Reduction) space via the eigenvector $\hat{b}$ associated with the largest nonzero eigenvalue of the matrix of interest $\widehat{\Sigma}_n^{-1}\widehat{\Gamma}_n$. Thus, $\hat{b}$ is an $EDR$ direction.
SIR_bootstrap
functionApply a single-index $SIR$ on $B$ bootstraped samples of $(X,Y)$ with $H$ slices.
SIR_threshold
functionApply a single-index $SIR$ on $(X,Y)$ with $H$ slices, with a parameter $\lambda$ which apply a soft/hard thresholding to the interest matrix $\widehat{\Sigma}_n^{-1}\widehat{\Gamma}_n$.
SIR_threshold_opt
functionApply a single-index SIR on $(X,Y)$ with $H$ slices, with a soft/hard thresholding of the interest matrix $\widehat{\Sigma}n^{-1}\widehat{\Gamma}_n$ by an optimal parameter $\lambda{opt}$. The $\lambda_{opt}$ is found automatically among a vector of n_lambda
$\lambda$, starting from 0 to the maximum value of $\widehat{\Sigma}n^{-1}\widehat{\Gamma}_n$. For each feature of $X$, the number of $\lambda$ associated with a selection of this feature is stored (in a vector of size $p$). This vector is sorted in a decreasing way. Then, thanks to strucchange::breakpoints
, a breakpoint is found in this sorted vector. The coefficients of the variables at the left of the breakpoint, tend to be automatically toggled to 0 due to the thresholding operation based on $\lambda{opt}$, and so should be removed (useless variables). Finally, $\lambda_{opt}$ corresponds to the first $\lambda$ such that the associated $\hat{b}$ provides the same number of zeros as the breakpoint's value.
For example, for $X \in \mathbb{R}^{10}$ and n_lambda=100
, this sorted vector can look like this :
| X10 | X3 | X8 | X5 | X7 | X9 | X4 | X6 | X2 | X1 | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | | 2 | 3 | 3 | 4 | 4 | 4 | 6 | 10 | 95 | 100 |
Here, the breakpoint would be 8.
SIR_threshold_bootstrap
functionApply a single-index optimally soft/hard thresholded $SIR$ with $H$ slices on n_replications
bootstraped replications of $(X,Y)$. The optimal number of selected variables is the number of selected variables that came back most often among the replications performed. From this, we can get the corresponding $\hat{b}$ and $\lambda_{opt}$ that produce the same number of selected variables in the result of SIR_threshold_opt
.
This method allow to validate the robustness of the model.
In this section :
Let's create $Y = (X\beta)^3 + \epsilon$ with:
set.seed(4) n <- 200 # Sample size p <- 30 # Number of variables in X p_star <- 10 # Number of relevant variables in X X <- mvtnorm::rmvnorm(n,sigma=diag(p)) # X ~ N(0,I_p) dimnames(X) <- list(1:n, paste("X", 1:p, sep = "")) # Rename columns of X eps <- rnorm(n, sd = 20) # Error beta <- matrix(c(rep(1,p_star),rep(0,p-p_star)),ncol=1) # Beta = Heaviside function rownames(beta) <- colnames(X) # Rename rows of X Y <- (X %*% beta)**3 + eps # The model
We can now plot Y versus the reconstructed index $X\beta$
plot(X %*% beta, Y, xlab = "true index")
library(SIRthresholded)
The SIR
function is applied to these data.
res_SIR = SIR(Y = Y, X = X, H = 10,graph = FALSE)
The object res_SIR
is an object of class SIR
and contains many numerical results (see ?SIR for more details), summarized in the S3 methods print
and summary
.
summary(res_SIR)
You can now access the estimation of the EDR direction with res_SIR$b
. Thus, it is possible to compute the correlation between $\beta$ and $\hat{b}$ :
cor(c(beta),c(res_SIR$b))
Here are the two graphical outputs obtained with the method plot
of the objects of class SIR
.
plot(res_SIR,choice="estim_ind") plot(res_SIR,choice="eigvals")
The first plot shows the reconstructed index by plotting $Y$ versus $X\hat{b}$ with $\hat{b}$ the estimated $EDR$ direction.
The second one shows the ten largest eigenvalues of the interest matrix. The $EDR$ direction is the eigenvector $\hat{b}$ associated with the largest eigenvalue. This plot can be useful to check if the largest eigenvalues are very close or not. In the case where they are close, it could be interesting to try to choose the second or third largest eigenvector to reconstruct the index. You can access the interest matrix with res_SIR$M1
, and then compute the eigenvectors and eigenvalues with the function eigen()
.
By looking at $\hat{b}$, we notice that the coefficient associated to some variables are close to $0$. We assume that some variables are useless and that it would be interesting to do variable selection.
An other approach is to directly use the optimally thresholded version of $SIR$ which will provide an estimate of the $EDR$ direction, and select the relevant variables. Let's start with $ST-SIR$ (soft thresholded SIR) with 100 $\lambda$ values to test :
res_STSIR = SIR_threshold_opt(Y=Y, X=X, H=10, n_lambda=100, thresholding="soft", graph=FALSE)
The object res_STSIR
is an object of class SIR_threshold_opt
and contains many numerical results (see ?SIR_threshold_opt for more details), summarized in the S3 methods print
and summary
.
summary(res_STSIR)
Here are some graphical outputs obtained with the method plot
of the objects of class SIR_threshold_opt
.
plot(res_STSIR,choice="cos2_selec")
This plot shows the evolution of $cos^2(\hat{b}{ST-SIR},\hat{b}{SIR})$ and the percentage of variables selected by the model, according to the value of $\lambda$. This graph is very informative as it shows the compromise between variable selection and distortion of the estimated direction. The $\hat{\lambda}_{opt}$ must provide a balance between these two.
plot(res_STSIR,choice="opt_lambda")
This plot shows how the optimal lambda has been chosen. For each feature of $X$, the percentage of $\lambda$ associated with a selection of this feature is stored in a vector, which is sorted and plotted here. Thanks to strucchange::breakpoints
, a breakpoint is found in this sorted vector, which allows to deduce the $\lambda_{opt}$ as explained in the Overview.
plot(res_STSIR,choice="regul_path")
This plot shows the regularization path of $\hat{b}$ according to the value of $\lambda$. Here, we can see that according to $\lambda$, some variables rapidly converge to 0, whereas the 10 relevant variables are more difficult to toggle to 0.
Then, let's plot plot Y versus the estimated first index $X\hat{b}$.
plot(res_STSIR,choice="estim_ind")
We can also see some graphical outputs of SIR_threshold_opt
with hard thresholding :
res_HTSIR = SIR_threshold_opt(Y=Y, X=X, H=10, n_lambda=100, thresholding="hard", graph=FALSE) plot(res_HTSIR,choice="cos2_selec") plot(res_HTSIR,choice="regul_path")
The hard thresholding method results in a more abrupt variable selection according to the $\lambda$s, wich explains that the curves are straighter with bigger jumps between values.
As the soft thresholding deforms the matrix of interest with the $\lambda$ parameter, the estimation of $\hat{b}$ is biased. That is why, at the end of the S3 summary method, it is proposed to create a new SIR model based only on the relevant variables, selected by the method. Let's copy paste the command line from the result of SIR_threshold_opt
and see the results:
res_SIR_pstar = SIR(Y=res_STSIR$Y, X=res_STSIR$X_reduced, H=res_STSIR$H) summary(res_SIR_pstar)
Let's plot Y versus the estimated first index $X\hat{b}$
plot(res_SIR_pstar,choice="estim_ind")
We directly notice that this estimated index is very close to the theoretical index plotted earlier.
As done previously, let's compute the correlation between $\beta$ and $\hat{b}{ST-SIR-p^\star}$ and compare it with the results of classic SIR. As the real $\beta$ vector has a length of 30 and $\hat{b}{ST-SIR-p^\star}$ has a length of 10, we need to create a new vector of size 30 to be able to compare them.
b_extended <- matrix(rep(0,p),nrow=1) # Create the empty vector colnames(b_extended) <- colnames(X) # Assign the colnames of X # Assign the values of b_extended according to the colnames b_extended[which(colnames(b_extended) %in% colnames(res_SIR_pstar$b))] = res_SIR_pstar$b
Let's compute the correlation :
cor(c(beta),c(b_extended))
We notice that the estimation of $\beta$ is better than the classic $SIR$ method (95% for $SIR$ and 99% for $ST-SIR-p^\star$). The 20 useless variables have been eliminated, and the estimation of the $EDR$ direction is very good.
Wether using hard or soft thresholding, it is always relevent to create a new $SIR$ model considering only the relevant variables. It allows to improve the estimation of the $EDR$ direction.
If you want to experiment the thresholding operation by yourself, you can call SIR_threshold
and give your own $\lambda$ to threshold the interest matrix :
res_SIR_thresh = SIR_threshold(Y, X, H = 10, lambda = 0.04, thresholding = "hard") summary(res_SIR_thresh)
To test the robustness of the model, it is relevant to perform bootstrapping. This method is even more relevant when the sample size is small. Lets call SIR_threshold_bootstrap
on the data, with 10 replications (only 10 to simplify vignettes building), a bootstrapped sample size multiplied by 2 (k=2) and the hard thresholding method :
res_SIR_boot = SIR_threshold_bootstrap(Y,X,H=10,n_lambda=100,thresholding="hard", n_replications=10,k=2,graph = FALSE)
The object res_SIR_boot
is an object of class SIR_threshold_bootstrap
and contains many numerical results (see ?SIR_threshold_bootstrap for more details), summarized in the S3 methods print
and summary
.
summary(res_SIR_boot)
Here are some graphical outputs obtained with the method plot
of the objects of class SIR_threshold_bootstrap
.
plot(res_SIR_boot,choice="size")
This barplot shows the distribution of the number of variables keeped by the model.
plot(res_SIR_boot,choice="selec_var")
This barplot shows the percentage of selection by the model for every variables.
plot(res_SIR_boot,choice="lambdas_replic")
This boxplot shows the distribution of the $\lambda_{opt}$ choosen over the replications. Moreover, the $\lambda_{opt}$ that has been selected (see how in the overview) is indicated with a purple horizontal line.
plot(res_SIR_boot,choice="coefs_b")
This plot shows the distribution of every coefficients of the $\hat{b}$ computed over the replications. The final coefficients of $\hat{b}$ that have been computed (see how in the overview) are indicated with a purple dot.
As done with SIR_threshold_opt
, we can also create a new SIR model considering only the relevant variables.
In this section :
For this example, we chose the wine dataset from the UCI machine learning repository (the official link is not working anymore, please find another version of the dataset on that link https://gist.githubusercontent.com/tijptjik/9408623/raw/b237fa5848349a14a14e5d4107dc7897c21951f5/wine.csv). As explained on the page description of the dataset, these data are the results of a chemical analysis of wines grown in the same region in Italy. The analysis determined the quantities of 12 constituents found in the 178 samples of wine. Here, we'll try to find the relevant variables in order to perform regression on the response variable Nonflavanoids according to the 12 other constituents. First, let's load the dataset.
wine <- read.csv("https://gist.githubusercontent.com/Clement-W/26d99a28ab89929b6321f70a04535451/raw/169a36568d1fc7c5a7c508e8d3e720d5040744cd/wine.csv", header = FALSE) colnames(wine) <- c('Type', 'Alcohol', 'Malic', 'Ash', 'Alcalinity', 'Magnesium', 'Phenols', 'Flavanoids', 'Nonflavanoids', 'Proanthocyanins', 'Color', 'Hue', 'Dilution', 'Proline') # Extract the response variable Y <- wine$Nonflavanoids # Remove the response variable (Nonflavanoids) and the class information (type of cultivars) X <- wine[, -which(names(wine) %in% c("Type","Nonflavanoids"))] head(cbind(Y,X),3) print(dim(X))
The X data contains 178 samples of 12 variables. As suggested in the description page of this dataset, we'll standardise the variables with the scale()
method. Here, scaling the data is mandatory to ensure that the SIR method works properly.
X = scale(X) Y = scale(Y)
Now, we can use SIR_threshold_opt
to do variable selection with an optimal thresholding parameter. The sample size is relatively small. To improve the estimations of $p_h$ and $m_h$ in the SIR method (more details in the first part of this vignette), we reduce the number of slices to 5, so there is more sample in each slices.
res1 = SIR_threshold_opt(Y=Y, X=X, H=5, n_lambda=300, thresholding="soft", graph=FALSE) summary(res1)
In this first model, 2 variables are selected : Ash and Flavanoids. Let's explore the graphs to better understand what is going on.
plot(res1,choice="estim_ind")
We first note that the model is very noisy. However, we still notice a linear tendency between the estimated index and Y.
Let's now focus on variable selection.
plot(res1,choice="cos2_selec") plot(res1,choice="regul_path")
These two graphs helps us to understand the link between variable selection (the coefficients of $\hat{b}$) and $\lambda$. Here, the $\lambda_{opt}$ allows to select 2 variables. In the first plot, we notice that if $\lambda$ were a bit less than $\lambda_{opt}$, we could have selected 2 or 3 more variables, while improving the estimation of $\hat{b}$. In fact, in the second plot, we notice that the coefficients of $\hat{b}$ associated to the variables Dillution, Magnesium and Hue are set to 0 by $\lambda_{opt}$.
As the data is noisy, the sample size is small, and the $\lambda_{opt}$ obtained with SIR_threshold_opt
is questionable, it is relevant to perform SIR_threshold_opt on bootstrapped replications of the dataset. Let's perform 500 replications, with a bootstrapped sample size multiplied by 2 :
# To lighten the build of the vignette, the result of this command was saved into a RData file. #res2 = SIR_threshold_bootstrap(Y=Y, X=X, H=5, n_lambda=200, thresholding="hard", n_replications = 500 , graph=FALSE) load("../R/sysdata.rda") # load res2 summary(res2)
Let's explore some of the graphical outputs :
plot(res2,choice="size")
Over the 500 replications, 50% of the models selected 2 variables. Then, approximately 45% of the models selected 3,4,5 or 6 variables. Let's take a look at the percentage of variable selection.
plot(res2,choice="selec_var")
As expected, a large majority of the models have selected the variables Ash and Flavanoids. But it is interesting to see that 4 other variables were regularly selected (between 30 and 40% of the replications) : Magnesium, Phenols, Hue and Dilution. We notice that these variables were closed to be selected in SIR_threshold_opt
.
plot(res2,choice="coefs_b")
Looking at the boxplot of the value of $\hat{b}$ over the 500 replications, we also find the same results. The coefficients of $\hat{b}$ associated to Magnesium, Phenols, Hue and Dilution are greater than 0 in median. It would then be interesting to create a SIR model, that contains these variable.
Let's look at the regularization path plotted earlier, and find the $\lambda$ which select the 4 potentially relevant variables. By looking at the top axis of the regularization path, we find the first lambda which select 6 variables, which is approximately 0.07 :
res3 = SIR_threshold(Y=Y, X=X, H=5, lambda = 0.071, thresholding="hard", graph=FALSE) summary(res3)
Now, the 6 variables are selected with this estimation of $\hat{b}.$. Let's estimate a new SIR model only on these relevant variables.
Let's copy paste the command from the previous output of the S3 summary method, and plot the estimated index versus Y :
res4 = SIR(Y=res3$Y, X=res3$X_reduced, H=res3$H, graph = FALSE) summary(res4)
plot(res4,choice="estim_ind")
We can note that the estimated index seems to be a bit less noisy, the linear tendency is easier to see.
Visually, the link between Y and X seems to be linear. To estimate the quality of our model, and to evaluate the relevance of the selected variable, we can build a linear model on the data.
Let's start with a model containing every variables :
summary(lm(Y~X))
The R-squared is equal to 44.33%, so aproximately 44% of the variability is explained by this model. We also notice that the coefficients associated to the variables Ash, Magnesium, Flavanoids, Hue and Dilution are significants.
Let's compare with the R-squared obtained with the variables selected by the previous SIR model :
The first model res1 which selected Flavanoids and Ash :
r
summary(lm(Y~Flavanoids+Ash,data=wine))$r.squared
The model with a custom thresholding parameter res3 which selected Flavanoids, Ash, Magnesium, Dilution, Phenols, Hue :
r
summary(lm(Y~Ash+Magnesium+Phenols+Flavanoids+Hue+Dilution,data=wine))$r.squared
The best compromise seems to be the SIR model that keep 6 variables. The R-squared is almost identical (-1%) and half of the variables have been kept.
The package SIRthresholded helped us to find the relevant variables, in order to create the best non-parametric regression on the response variable Nonflavanoids. This example shown that only using SIR_threshold_opt
can be limiting on real data, especially when the data is noisy with a small sample size. The method SIR_threshold_bootstrap
helped us to counter this issue, and retrieve important information about the potential relevant variables. Finally, the best model is res4
which is a SIR model only based on the 6 relevant variables : Flavanoids, Ash, Magnesium, Dilution, Phenols, Hue.
Duan, N., & Li, K.-C. (1991). Slicing Regression : A Link-Free Regression Method. The Annals of Statistics, 19(2), 505‑530.
Dua, D. and Graff, C. (2019). UCI Machine Learning Repository. Irvine, CA: University of California, School of Information and Computer Science.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.