Subdata selection algorithm for linear model discrimination

Abstract

A statistical method is likely to be sub-optimal if the assumed model does not reflect the structure of the data at hand. For this reason, it is important to perform model selection before statistical analysis. However, selecting an appropriate model from a large candidate pool is usually computationally infeasible when faced with a massive data set, and little work has been done to study data selection for model selection. In this work, we propose a subdata selection method based on leverage scores which enables us to conduct the selection task on a small subdata set. Compared with existing subsampling methods, our method not only improves the probability of selecting the best model but also enhances the estimation efficiency. We justify this both theoretically and numerically. Several examples are presented to illustrate the proposed method.

Appendices

Appendix

Technical details

Proof of Theorem 1

For any given subdata $X^*$, by applying the entropy decomposition in information theory (Sebastiani and Wynn 2000, Equation (2)), the joint entropy of $Y^*$ and the parameter $\varTheta$ can be decomposed as

$$\begin{aligned} \begin{aligned} \mathrm {Ent}(Y^*,\varTheta |X^*)&= \mathrm {Ent}(\varTheta |X^*)+ E_{\varTheta }\{\mathrm {Ent}(Y^*|\varTheta ,X^*)\} \\&= \mathrm {Ent}(\varTheta )+ \mathrm {Ent}(\varvec{\varepsilon }^*), \end{aligned} \end{aligned}$$

(A.1)

where $\varvec{\varepsilon }^*$ stands for the corresponding error term of $(Y^*,X^*)$ in model (1). The second equality holds by the model assumption since all the randomness of $Y^*$ comes from the error term conditional on $\varTheta ,X^*$ and $\varTheta$ is functionally independent of $X^*$. This implies that $\mathrm {Ent}(Y^*,\varTheta |X^*)$ is a constant up to the subdata size r.

Also note that $\mathrm {Ent}(Y^*,\varTheta |X^*)$ can also be decomposed as

$$\begin{aligned} \mathrm {Ent}(Y^*,\varTheta |X^*)=\mathrm {Ent}(Y^*|X^*)+E_{Y^*}\{\mathrm {Ent}(\varTheta |Y^*,X^*)\}. \end{aligned}$$

(A.2)

That is to say maximizing $\mathrm {Ent}(Y^*|X^*)$ indicates minimizing the overall expected deviance loss $E_{Y^*}\{\mathrm {Ent}(\varTheta |Y^*,X^*)\}$.

Now we turns to calculate $\mathrm {Ent}(Y^*|X^*)$. Without loss of generality, we assume that the first $p_{(t)}$ columns of X be the model matrix of the true model. Thus the prior of $\beta _{(t)}$ comes from $N(\beta _\mathrm{prior,(t)},\sigma _f^2I_{p_{(t)}})$, where $\beta _\mathrm{prior,(t)}$ corresponds to the first $p_{(t)}$ entries of $\beta _\mathrm{prior}$. Note that $Y^*=X_{(t)}^*\beta _{(t)}+\varvec{\varepsilon }$ and the prior of $\beta _{(t)}$ obeys $N(\beta _\mathrm{prior,(t)},\sigma _f^2I_{p_{(t)}})$. Thus the marginal distribution of $Y^*$ is normal with mean $X_{(t)}^*\beta _\mathrm{prior,(t)}$ and variance $\sigma ^2I_t+\sigma _f^{2}X_{(t)}^{*{\mathrm {T} }}X_{(t)}^*$ under model (1). The desired results come from the facts

$$\begin{aligned} \begin{aligned} \mathrm {Ent}(Y^*|X^*)&=\log \det (\sigma ^2I_r+\sigma _f^{2}X_{(t)}^{*}X_{(t)}^{*{\mathrm {T} }})+c_1\\&=\log \det (\sigma ^{-2}\sigma _f^{2}X_{(t)}^{*{\mathrm {T} }}X_{(t)}^{*}+I_t)+c_2,\\&=\log \det (X_{(t)}^{*{\mathrm {T} }}X_{(t)}^{*}+\sigma ^{2}\sigma _f^{-2}I_t)+c_3, \end{aligned} \end{aligned}$$

(A.3)

where $c_1,c_2,c_3$ are some constant up to the subdata size r. The second equality comes from the matrix determinant lemma, i.e., $\det (A+BC)=\det (A)\det (I+CA^{-1}B)$ for some matrices A, B, C with $A>0$. $\square$

Proof of Theorem 2

It is sufficient to show that $X^{*{\mathrm {T} }}X^*\le (\sum _{i=1}^n\delta _i h_{ii}) X^{{\mathrm {T} }}X$ in the sense of Loewner ordering. Let ${{x}}_i$ be the ith row of X. For any $a\in \mathbb {R}^{p}$, noting that $X^{\mathrm {T} }X$ is a full rank matrix, a can be represent as $a=(X^{\mathrm {T} }X)^{-1/2} b$ for some $b\in \mathbb {R}^{p}$. Then,

$$\begin{aligned} a^{\mathrm {T} }{{x}}_i {{x}}_i^{\mathrm {T} }a= & {} b^{\mathrm {T} }(X^{\mathrm {T} }X)^{-1/2} {{x}}_i^{\mathrm {T} }{{x}}_i(X^{\mathrm {T} }X)^{-1/2} b \nonumber \\\le & {} \mathrm {tr}\{(X^{\mathrm {T} }X)^{-1/2} {{x}}_i^{\mathrm {T} }{{x}}_i(X^{\mathrm {T} }X)^{-1/2}\}\Vert b\Vert _2^2 \nonumber \\= & {} h_{ii}\{b^{\mathrm {T} }(X^{\mathrm {T} }X)^{-1/2}(X^{\mathrm {T} }X)(X^{\mathrm {T} }X)^{-1/2} b\} \end{aligned}$$

(A.4)

$$\begin{aligned}= & {} h_{ii} a^{\mathrm {T} }(X^{\mathrm {T} }X) a, \end{aligned}$$

(A.5)

where $\mathrm {tr}(\cdot )$ is the trace operator and $(X^{\mathrm {T} }X)^{-1/2}(X^{\mathrm {T} }X)^{-1/2}=(X^{\mathrm {T} }X)^{-1}$. Therefore, ${{x}}_i {{x}}_i^{\mathrm {T} }\le h_{ii}X^{\mathrm {T} }X$ and the desired result comes from summing over the both side of the inequality. $\square$

For clarity, we begin with the proof of the following lemma since some results in the following lemma will be used in the proof of Theorem 3.

Lemma 1

Assume that $n^{-1}X^TX$ goes to a positive definite matrix. Let $\hat{\beta }_k^*$ be the MLE based on selected subdata set according to Algorithm 1 for the kth candidate model. As $r\rightarrow \infty ,n\rightarrow \infty$, the following result holds:

$$\begin{aligned} \mathrm{Var}(\hat{\beta }_k^*)=O\left( \frac{1}{n\sum _{i=1}^rh_{(ii)}}\right) . \end{aligned}$$

(A.6)

Proof of Lemma 1

According to Algorithm 1, $X^*=U_\varGamma \varSigma V^{\mathrm {T} }$. Then it is sufficient to show that

$$\begin{aligned} c\sum _{i=1}^{r}h_{(ii)} \le \lambda _{\min } (U_\varGamma ^TU_\varGamma ) \le \lambda _{\max } (U_\varGamma ^TU_\varGamma )\le \sum _{i=1}^{r}h_{(ii)}, \end{aligned}$$

(A.7)

for some constant c, where $\lambda _{\max }(A)$, $\lambda _{\min }(A)$ stand for the maximum and minimum eigenvalue of A, respectively. Since $U_\varGamma ^{\mathrm {T} }U_\varGamma$ is positive definite through Algorithm 1, therefore $\lambda _{\max } (U_\varGamma ^{\mathrm {T} }U_\varGamma )\le \mathrm {tr}(U_\varGamma ^{\mathrm {T} }U_\varGamma )=\sum _{i=1}^{r}h_{(ii)}$. By the definition of the condition number, it holds that

$$\begin{aligned} \lambda _{\min } (U_\varGamma ^{\mathrm {T} }U_\varGamma )=\lambda _{\max } (U_\varGamma ^{\mathrm {T} }U_\varGamma )/\kappa (U_\varGamma ^{\mathrm {T} }U_\varGamma )\ge \mathrm {tr}(U_\varGamma ^{\mathrm {T} }U_\varGamma )/(pT), \end{aligned}$$

where the last inequality comes from the fact that $\mathrm {tr}(U_\varGamma ^{\mathrm {T} }U_\varGamma )\le p\lambda _{\max }(U_\varGamma ^{\mathrm {T} }U_\varGamma )$.

From (A.7), it follows that

$$\begin{aligned} c\sum _{i=1}^{r}h_{(ii)}V\varSigma ^2 V^T \le X^{*T}X^* \le \sum _{i=1}^{r}h_{(ii)}V\varSigma ^2 V^T, \end{aligned}$$

(A.8)

and the desired results follows by noting $X^TX=V\varSigma ^2 V^T$ and $\mathrm{Var}(\hat{\beta }_k)=(P^TX^{*T}X^*P)^{-1}$ for the projection matrix P such that $X_{(k)}=XP$ where $X_{(k)}$ is the design matrix for model $S_k$. $\square$

Now, let us turn to proof Theorem 3.

Proof of Theorem 3

Denote $\mathcal {M}^C$ be the set of correct candidate models, and $\mathcal {M}^I=\mathcal {M}-\mathcal {M}^C$ be the set of incorrect candidate models. We first show that

$$\begin{aligned} \varDelta (k)=\liminf _r\min _{S_k\in \mathcal {M}^I }\Vert \mu ^*-H_{(k)}^*\mu ^*\Vert ^2/\log r\rightarrow \infty , \end{aligned}$$

(A.9)

where $\mu ^*$ stands for the mean of the selected data, $H_{(k)}^*=X_{(k)}^*(X_{(k)}^{*{\mathrm {T} }}X_{(k)}^*)^{-1}X_{(k)}^{*{\mathrm {T} }}.$

For any candidate model in $\mathcal {M}^{I}$, say $S_k$ as an example, let the model matrix for the closest correct model be $\tilde{X}_{(k)}^*:=(X_{(\check{c})}^*,X_{(k)}^*)$. Here $X^*_{(\check{c})}$ stands for the “complementary” design, which consists of the columns of $X_{(t)}^*$ that are not included in $X_{(k)}^*$. Denote the regression coefficient vector corresponding to $X_{(\check{c})}^*$ as $\beta _{(\check{c})}$, which is a subvector of $\beta _{(t)}$. Direct calculation yields

$$\begin{aligned} \Vert \mu ^*-H_{(k)}^*\mu ^*\Vert ^2= & {} \inf _{\alpha }\Vert X^*_{(\check{c})}\beta _{(\check{c})}-X^*_{(k)}\alpha \Vert ^2 \end{aligned}$$

(A.10)

$$\begin{aligned}= & {} \inf _{\alpha }\{(\beta _{(\check{c})}^{\mathrm {T} },\alpha ^{{\mathrm {T} }})(\tilde{X}_{(k)}^{*{\mathrm {T} }}\tilde{X}^*_{(k)})(\beta _{(\check{c})}^{\mathrm {T} },\alpha ^{{\mathrm {T} }})^{{\mathrm {T} }}\}. \end{aligned}$$

(A.11)

Utilizing the results in (A.8), we have

$$\begin{aligned} \left( c_2n\sum _{i=1}^{r}h_{(ii)}\right) \left( \frac{1}{n}X^{\mathrm {T} }X\right) \le X^{*{\mathrm {T} }}X^* \le \left( n\sum _{i=1}^{r}h_{(ii)}\right) \left( \frac{1}{n}X^{\mathrm {T} }X\right) , \end{aligned}$$

(A.12)

for some constant $c_2$. Note that the $\tilde{X}_{(k)}$ is a submatrix of X up to a column permutation. Thus $\lambda _{\min }(\tilde{X}_{(k)}^{*{\mathrm {T} }}\tilde{X}^*_{(k)})\ge \lambda _{\min }(\tilde{X}^{*{\mathrm {T} }}\tilde{X}^*)=O\left( n\sum _{i=1}^{r}h_{(ii)}\right)$, where $\lambda _{\min }(\cdot )$ stands for the smallest eigenvalue of a squared matrix. From (A.10), we have

$$\begin{aligned} \liminf _{n,r\rightarrow \infty }\min _{S_k\in \mathcal {M}^I }\Vert \mu ^*-H_{(k)}^*\mu ^*\Vert ^2/\log r\ge \liminf _{n,r\rightarrow \infty }\min _j(n\sum _{i=1}^{r}h_{(ii)})\Vert \beta _{(t)j}\Vert ^2/\log r\rightarrow \infty , \end{aligned}$$

which implies (A.9) holds.

For convenience, let ${\text {BIC}^*}(S_k)=r\log \left( r^{-1}\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(k)})^2\right) +(p_{(k)}+1)\log r.$ From (3.7) in Shao (1997), for any model $S_k$ in $\mathcal {M}^I$, we have

$$\begin{aligned} \sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(k)})^2-\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(t)})^2\ge \varDelta (k)\ge p\log r>0, \end{aligned}$$

(A.13)

which implies $\log (\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(k)})^2)-\log (\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(t)})^2)>r\log (1+p\log r/r)$ under the assumption $\hat{\sigma }^*\not \rightarrow 0.$ Therefore,

$$\begin{aligned} {\text {BIC}^*}(S_k)-{\text {BIC}^*}(S_t)\ge r\log (1+{p\log r}/{r})-(p-p_{(k)})\log r\rightarrow \infty .\nonumber \\ \end{aligned}$$

(A.14)

Similarly, for any model $S_{k^\prime }$ in $\mathcal {M}^C$ with $p_{(k^\prime )}>p_{(t)},$ where $p_{(t)}$ is the column dimension of $X_{(t)}$, it is straightforward to see that

$$\begin{aligned} r\log \left( r^{-1}\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(k^\prime )})^2\right) -r\log \left( r^{-1}\sum _{i=1}^r(y_i^*-{\hat{\mu }}_i^{*(t)})^2\right) \rightarrow \chi ^2_{p_{(k^\prime )}-p_{(t)}},\nonumber \\ \end{aligned}$$

(A.15)

according to the log likelihood ratio test (see van der Vaart 1998, Chapter 16). Therefore, it holds that

$$\begin{aligned} {\text {BIC}^*}(S_{k^\prime })-{\text {BIC}^*}(S_t)=O(\log r)\rightarrow \infty . \end{aligned}$$

(A.16)

Combining (A.14) and (A.16), we can get the desired result. $\square$

Proof of Theorem 4

This is the direct result from Lemma 1.

Fig. 5

Selection accuracies for different r based on UNIF (solid line with circles), LEVSS (dashed line with triangles), LEV (dotted line with squares) and IBOSS (dotdash line with plus signs) method for the six cases listed in Sect. 5.1. The models are selected based on the forward regression procedure via BIC

Fig. 6

Log MSPEs for different r based on UNIF (solid line with circles), LEVSS (dashed line with triangles), LEV (dotted line with squares) and IBOSS (dotdash line with plus signs) methods for the six cases listed in Sect. 5.1. The models are selected based on the forward regression via BIC

Additional simulation results on forward regression

Since the number of possible models increases exponentially with p, all-subset regression is only feasible for the cases that p is relatively small. Alternatively, a forward selection approach is usually adopted. More precisely, the forward regression starts from the null model, and iteratively adds one variable to the currently “best” model which yields the lowest value for the BIC at a time. This process is repeated until no more variables should be added into the currently “best” model. In this part, we adapt the forward regression to illustrate the proposed method. Of course, a backward elimination procedure and a step-wise regression procedure can also be adopted. Since the three methods have similar performance, we only report the results on forward regression.

Fig. 7

Selection accuracies for different r based on UNIF (solid line with circles), LEVSS (dashed line with triangles), LEV (dotted line with squares) and IBOSS (dotdash line with plus signs) method for the six cases listed in Sect. 5.1. The models are selected via Lasso

In accordance with Sect. 5.1, we also demonstrate our method as well as the uniform subsampling, leveraging score subsampling and IBOSS through the six cases. The selection accuracies for the six cases are presented in Figure 5. The log MSPEs are also provided in Figure 6 to evaluate the performance of prediction.

From Figures 5, and 6, one can see that the forward regression results based on BIC  are very similar to the results on the all-subset regressions on BIC.

Additional simulation results on Lasso

Now, we will exam our method’s performance on model selection according to Lasso (Tibshirani 1996).

To be aligned with the settings described in Sect. 5.1, we also demonstrate our methods as well as the other three methods (i.e., uniform subsampling, leverage score subsampling, and the IBOSS) on the six cases listed at the beginning of Sect. 5.1 and evaluate the selection performance through the selection accuracies and MSPEs. The Lasso method is conducted through the glmnet package (Simon et al. 2011) and the tuning parameters are selected through 10-fold cross-validation according to the cv.glmnet() function. As for the leverage score subsampling, the Lasso is conduct as in Leng and Leung (2011).

Results on the selection accuracies are presented in Figure 7. It can be seen that the selection results based on Lasso are very similar to the results on the subset selection based on BIC.

To see the benefits of the model selection, we also report the log MSPEs in Figure 8 with $n_\mathrm{test}=500.$ From Figure  8, we can clearly see that the IBOSS and our method (LEVSS) are uniformly performance better than the uniform subsampling method.

Fig. 8

Log MSPEs for different r based on UNIF (solid line with circles), LEVSS (dashed line with triangles), LEV (dotted line with squares) and IBOSS (dotdash line with plus signs) methods for the six cases listed in Sect. 5.1. The models are selected via Lasso