Variable Selection in Doubly Truncated Regression

We consider variable selection for model (1) based on the observed data points. The loss function (3) motivates a penalized approach for variable selection. Instead of restricting to fixed dimension of $X_{i}$, we allow the dimension to go to infinity with $n$, that is, $p$ depends on $n$, denoted by $p_{n}$. Then all the quantities that are functions of the covariates depend on $n$. For the simplicity of notation, we shall suppress the subscript $n$ when there is no confusion. Meanwhile, $n/\tilde{n}$ is assumed to converge to a positive constant less than 1 as $\tilde{n}$ goes to infinity.

Because of the $L_{1}$ form of the loss function (3), the LASSO-type penalties are naturally considered. Specifically, the objective function with usual LASSO penalty is given by

where $\lambda_{n}$ is a (data-dependent) tuning parameter. There exist two main concerns if we want to use the above LASSO penalty. The first one is about the computation. One can see that $l_{n}(\beta)$ is not convex in $\beta$, so the optimization is no longer convex. However, we note that $l_{n}(\beta)$ can be represented as

The non-convexity comes from the indicator function, while $|e_{i}(\beta)-e_{j}(\beta)|$ is convex in $\beta$. Thus, one can develop an iterative algorithm where the parameter in the indicator function is kept fixed in each iteration step and the optimization is implemented for the rest part of the objective function. Some more details of the computation are discussed in Section 3.2 and Appendix A.1.

The second one is that it is well known that the LASSO estimator does not have the unbiasedness property and there are scenarios where the LASSO selection is not consistent. To overcome these drawbacks, Zou (2006) proposed adaptive LASSO penalty, in which proper weights are introduced for penalizing different coefficients in the lasso penalty. The weights used by Zou (2006) are defined based on root-$n$-consistent estimator of $\beta$. He showed that with a proper choice of the tuning parameter $\lambda_{n}$, the adaptive LASSO possesses the oracle properties when the dimension of $\beta$ is fixed.

Trying to keep the oracle property, we also consider the adaptive LASSO approach for variable selection. Specifically, the objective function with adaptive LASSO is given by

where $\hat{w}_{j}$ is the adaptive weight of $\beta_{j}$. Originally, Zou (2006) proposed to use $\hat{w}_{j}=1/|\tilde{\beta}_{j}|^{\gamma}$, where $\tilde{\beta}_{j}$ is a root-$n$-consistent estimator of $\beta_{j}$ and $\gamma$ is a predetermined positive constant. In our case, apparently, even when $p_{n}$ diverges with $n$, we can still get $\hat{\beta}$ by minimizing $l_{n}(\beta)$ as long as $p_{n}$ does not grow too fast. Based on $\hat{\beta}$, we define $\hat{w}_{j}=1/|\hat{\beta}_{j}|^{\gamma}$ in (4), where $\hat{\beta}_{j}$ stands for the $j$-th component of $\hat{\beta}$. The proposed adaptive LASSO estimator, conceptually, is defined to be the minimizer of $l_{n}^{\mbox{\tiny AL}}(\beta)$.

As we mention in Section 3.1, $l_{n}(\beta)$ is not convex, neither is $l_{n}^{\mbox{\tiny AL}}(\beta)$, so we turn to an iterative procedure to get the proposed estimator. Specifically, we modify $l_{n}^{\mbox{\tiny AL}}(\beta)$ to

Let $\hat{\beta}^{(0)}$ be an initial estimate; for example, one may use $\hat{\beta}$ as the start. Then, in the $k$-th iteration step, let $\hat{\beta}^{(k)}=\mbox{argmin}_{\beta}l_{n}^{\mbox{\tiny AL}}(\beta,\hat{\beta}^{(k-1)})$ for $k\geqslant 1$. If $\hat{\beta}^{(k)}$ converges to a limit as the number of $k\to\infty$, then we use the limit as the proposed estimator.

The main reason that we use the above iterative algorithm is because in the $k$th iteration, the objective function can be written as

where ${\cal S}_{k}=\{(i,j)|L_{j}(\hat{\beta}^{(k-1)})<e_{i}(\hat{\beta}^{(k-1)})<R_{j}(\hat{\beta}^{(k-1)}),L_{i}(\hat{\beta}^{(k-1)})<e_{j}(\hat{\beta}^{(k-1)})<R_{i}(\hat{\beta}^{(k-1)})\}$. Thus, the optimization problem becomes a linear programming problem which can be solved efficiently. Moreover, we find that after some algebraic operations, the optimization can be transformed into a least absolute deviation (LAD) problem which can be solved efficiently through standard LAD procedure, such as the quantreg package in R. More details on the transformation are provided in the Appendix A.1.

Next we establish the oracle property of the proposed adaptive lasso approach under some regularity conditions. The oracle property requires the consistency of $\hat{\beta}$. Although Ying et al. (2020) proved root-$n$-consistency of $\hat{\beta}$ under suitable conditions, their result is only confined to fixed dimension. When $p_{n}$ diverges with $n$, the consistency of $\hat{\beta}$ needs further investigation. Let $\beta_{0}$ denote the true value of $\beta$. The following proposition, proved in the Appendix A.2, gives out the consistency of $\hat{\beta}$ with the diverging dimension.

We then show the existence of an adaptive lasso estimator and it is consistent with the convergence rate of $O_{p}((p_{n}/n)^{1/2})$. The following theorem, proved in the Appendix A.2, states the results.

To present the following oracle properties, we need some more notation. Let $\xi_{ij}(\beta)=I\{L_{j}(\beta)<e_{i}(\beta)<R_{j}(\beta),L_{i}(\beta)<e_{j}(\beta)<R_{i}(\beta)\}(X_{i}-X_{j})\mbox{sgn}\{e_{i}(\beta)-e_{j}(\beta)\}$, $h_{ij}(\beta)=\left|\left[(e_{i}(\beta)-e_{j}(\beta))\wedge(R_{j}-Y_{j})\wedge(Y_{i}-L_{i})\right]\vee(L_{j}-Y_{j})\vee(Y_{i}-R_{i})\right|$, and $\bar{l}(\beta)=\mathsf{E}(h_{ij}(\beta))$, where the expectation is taken with respect to the distribution of the observed data. We define two $p_{n}\times p_{n}$ matrices $V=\mathsf{E}(\xi_{ij}(\beta_{0})\xi_{ik}^{\top}(\beta_{0}))$ and $A=\partial^{2}\bar{l}(\beta)/\partial\beta\partial\beta^{\top}|_{\beta=\beta_{0}}$. We further partition $\beta_{0}$ into $\beta_{0,1}$ and $\beta_{0,2}$, where $\beta_{0,1}$ stands for the nonzero component of $\beta_{0}$ and $\beta_{0,2}$ for the zero component. The dimension of $\beta_{0,1}$ is denoted by $p_{n,1}$. Let $\hat{\beta}^{\mbox{\tiny AL}}_{1}$ and $\hat{\beta}^{\mbox{\tiny AL}}_{2}$ be the adaptive lasso estimator for $\beta_{0,1}$ and $\beta_{0,2}$, respectively. Moreover, let $V_{11}$ denote the $p_{n,1}\times p_{n,1}$ components of $V$ corresponding to $\beta_{0,1}$ and $A_{11}$ denote the $p_{n,1}\times p_{n,1}$ components of $A$.

We give out the proof of Theorem 2 in the Appendix A.2. By choosing a suitable $u_{n}$, it can be seen from Theorem 2 that the estimator of each individual component of $\beta_{0,1}$ is asymptotically normal. The results can be used to make inference on individual coefficients. The estimation of the variance is discussed in Section 3.4.

Typical methods of selecting tuning parameter include information criteria methods such as AIC, BIC, and data-driven procedures such as $K$-fold cross-validation. Because of the computational issue, we consider information criteria approaches here. It is well known that BIC usually leads to the consistency of selection for fixed dimension scenarios. However, as we allow the dimension of the parameter to diverge, the standard BIC does not guarantee the selection consistency. Thus, we consider the following modified BIC:

where $\hat{\beta}^{\mbox{\tiny AL}}_{\lambda}$ represents the adaptive lasso estimate obtained at the tuning parameter with level $\lambda$, $C_{n}$ is a sequence going to infinity with $(\log n/n)C_{n}\to 0$, and $\mbox{df}_{\lambda}$ stands for the size of the non-zero component of the estimate at level $\lambda$. We find that $C_{n}=\log\log p_{n}$ works well in this paper. The tuning parameter we choose is given by

The proposed modified BIC has the consistency of selection. Let ${\cal M}^{\ast}$ and ${\cal M}_{\hat{\lambda}}$ be the true model and the selected model, respectively. We have the following result on selection consistency.

Finally, we also need to give out the variance estimation of the proposed estimator for the nonzero component of $\beta_{0}$. Note that the asymptotic variance-covariance matrix of $\hat{\beta}^{\mbox{\tiny AL}}_{1}$ involves the density function of $\tilde{\varepsilon}$ and is complicated to estimate directly. We turn to use a random weighting method to obtain the variance estimate. To be specific, we defined the perturbed penalized estimator, denoted by $\hat{\beta}^{\ast}$, to be the minimizer of the following perturbed version of the objective function $l^{\mbox{\tiny AL}}_{n}(\beta)$:

where $W_{i}$, $i=1,\ldots,n$ are i.i.d. non-negative and bounded random variables with mean 0.5 and variance 1, independent of the observed data. Let $\hat{\beta}^{\ast}_{1}$ be component of $\hat{\beta}^{\ast}$ corresponding to $\hat{\beta}^{\mbox{\tiny AL}}_{1}$, i.e., the non-zero part of $\hat{\beta}^{\mbox{\tiny AL}}$. It can be shown that under suitable conditions, given the observed data, the conditional asymptotic distribution of $\hat{\beta}^{\ast}_{1}$ is the same as the asymptotic distribution of $\hat{\beta}^{\mbox{\tiny AL}}_{1}$. Thus, by generating repeatedly the sequence of $W_{i}$, $i=1,\ldots,n$, we can obtain a large number of replications of $\hat{\beta}^{\ast}_{1}$ and then approximate the variance of $\hat{\beta}^{\mbox{\tiny AL}}_{1}$.

We conduct some simulation studies to illustrate the finite sample performance of the proposed approach. We set $p_{n}=[7\tilde{n}^{1/5}]$, where $\tilde{n}$ is the size of the full sample, as defined in Section 2. The size of the non-zero components is set to be $p_{n,1}=[p_{n}/3]$. We consider two scenarios with $\tilde{n}=300$ and $\tilde{n}=500$, respectively. In the first scenario, $p_{n}=21$ and $p_{n,1}=7$. We set $\beta_{1,0}=(3.12,2.20,-0.86,0.92,-2.49,1.95,-1.32)^{\top}$. For the covariates $\tilde{X}=(\tilde{X}_{1},\ldots,\tilde{X}_{21})$, we independently generate $\tilde{X}_{1}$ from Binomial(0.25) (the Binomial distribution with success probability 0.25), $\tilde{X}_{2}$ from Binomial(0.8), $\tilde{X}_{8}$ from Uniform(0,2) (the uniform distribution from 0 to 2), $\tilde{X}_{9}$ from Binomial(0.5), and $\tilde{X}_{10}$ from Uniform(-2,0). All the rest covariats are generated from the standard normal distribution, with the correlation between $\tilde{X}_{i}$ and $\tilde{X}_{j}$ being $0.3^{|i-j|}$. In the second scenario, $p_{n}=24$ and $p_{n,1}=8$. We set $\beta_{1,0}=(3.12,2.20,-0.86,0.92,-2.49,1.95,-1.32,-2.13)^{\top}$. We generate $\tilde{X}_{1}$ from Binomial(0.25), $\tilde{X}_{2}$ from Binomial(0.8), $\tilde{X}_{9}$ from Uniform(0,2), $\tilde{X}_{10}$ from Binomial(0.5), $\tilde{X}_{11}$ from Uniform(-2,0), and the rest from the standard normal distribution, with the correlation between $\tilde{X}_{i}$ and $\tilde{X}_{j}$ being $0.3^{|i-j|}$. For the error distribution in the linear model, we consider the standard normal ($N(0,1)$) distribution and the extreme value distribution with location parameter 0 and scale parameter 1 ($EV(0,1)$). The left truncation variable $\tilde{L}$ is generated from Uniform($a$,$b^{\top}X$) and the right one $\tilde{R}=\tilde{L}+c$, where $a$, $b$, and $c$ are constants to yield truncation percentage of 30% and 40%. The percentages of left and right truncation are almost the same.

Under each scenario, 1000 replications are carried out. For variable selection, we define several criteria to assess the performance. The first one is the model error (ME), being defined as ME$=(\tilde{\beta}-\beta_{0})^{\top}\mathsf{E}(X_{i}X_{i}^{\top})(\tilde{\beta}-\beta_{0})$, where $\tilde{\beta}$ is an estimator of $\beta_{0}$ derived from a specific approach. The median and median absolute deviation (MAD) of the 1000 MEs are recorded to evaluate the prediction performance of different procedures. The second one is the average number of total estimated zero coefficients (TN), that is, the average number of zero estimates obtained in the 1000 replications, including correctly estimated zero number (CN) and incorrectly estimated zero number (IN). For CN, the closer it is to $p_{n}-p_{n,1}$ the better a procedure performs, while for IN, the closer it is to 0 the better the procedure performs. The third one is the ratio of the correctly identified models (RCM), that is, the proportion of the replications that can identify the true model correctly among the 1000 replications. The closer it is to 100%, the better a procedure performs. We consider three procedures, including the proposed selection procedure (Proposed), the variable selection procedure with adaptive lasso penalty without taking truncation into consideration (Naive), and the oracle procedure where the correct subset of covariates is used to fit the model (Oracle). The results are summarized in the following tables.

[Insert Table 1 here]

[Insert Table 2 here]

From the tables, we see that the proposed method outperforms the naive approach significantly in terms of the ME. This is expectable since the naive approach yields biased estimators for the regression coefficients. For the TN and RCM, the proposed method also has better results than the naive approach, while the extend of the improvement is not very large. The TN and RCM are criteria for variable selection performance. Although selecting important variables seems to be less affected by the double truncation, correcting the bias induced by double truncation still helps in increasing the precision of the variable selection procedure.

We analyze an astronomical data from the 46420 quasar catalog produced by the Sloan Digital Sky Survey (SDSS) team. In this data, the main purpose is to predict the red shift of a quasar by using some related features. Thus, the linear model is used and the dependent variable is the red shift of the quasars. Compared with the complete SDSS catalog, the data set we use here omits some technical columns and therefore contains 23 covariates. Some more details about the data set can be found in https://astrostatistics.psu.edu/datasets.

It might not be that straightforward to see why there exists double truncation in the red shift data. This mainly comes from the nature of the observation feature of the astronomical telescope. Because of the aperture of the astronomical telescope, the red shift of a quasar one can observe is limited from the above. For the current data, we consider the maximum value of the red shift (5.4135) as the right truncation bound of the dependent variable. That is, the red shift larger this value is right truncated. For the low red shift, the measurement might be unreliable. The values should be omitted if the value of the red shift is lower its estimated error. The estimated error is then treated as the left truncation bound of the red shift observation. Therefore, the red shift value of the quasar is considered to be doubly truncated. Among the 23 covariates, some are the measurement errors of the brightness which are determined by the SDSS team from knowledge of the observing conditions, detector background, and other technical considerations. These covariates are not included in the model. We use 13 main covariates (shown in Table 3) to fit the linear model. Due to the computation capacity limitation, we do not use the full data to fit the model. A random sample with sample size 929 is drawn to illustrate the proposed method. The sampled data is split into the training part and testing part according to the proportion of 7:3. We use the proposed variable selection procedure on the training set to get a list of ”important” covariates and then used the selected model to do prediction on the testing set. The results of the variable selection are summarized in Table 4.

[Insert Table 3 here]

[Insert Table 4 here]

The proposed variable selection procedure gives out 5 ”important” covariates with non-zero coefficients. We use the selected covariates to do the prediction on the testing set. In Figure 1, we draw the line plot to compare the predicted values and the true response values. We see that the predicted values (shown by the dotted line) are quite close to the true values (shown by the solid line). The selected covariates give out reasonable predictions.

[Insert Figure 1 here]