Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,280 characters · 13 sections · 104 citation commands
Best Subset Binary Prediction
Prediction of binary outcomes is an important topic in economics and various scientific fields. Let $Y\in \{0,1\}$ be the binary outcome of interest and $ W$ a vector of covariates for predicting $Y$. Assume that the researcher has a training sample of $n$ independent identically distributed (i.i.d.) observations $\left( Y_{i},W_{i}\right) _{i=1}^{n}$ of $(Y,W)$. For $w\in \mathcal{W}$, let
where $\mathcal{W}$ is the support of $W$, $\theta $ is a vector of parameters, and $1\left\{ \cdot \right\} $ is an indicator function that takes value 1 if its argument is true and 0 otherwise.
One reasonable prediction rule is to choose $\theta $ such that it maximizes the probability of making the correct prediction $P\left[ Y=b_{\theta }(W) \right] $. However, this is infeasible in practice since the joint distribution of $(Y,W)$ is unknown. A natural sample analog is to maximize the sample average score which equals the proportion of correct predictions under the prediction rule ((ref)) in the training sample. This maximization problem is equivalent to the maximum score estimation in binary response models and is pioneered by Manski:1975, Manski:1985. Thus, we call the corresponding prediction rule the maximum score prediction rule. See Manski:Thompson:1989, Jiang:10, and Elliott2013 for prediction in the maximum score approach.
This prediction problem has the same structure as the binary classification problem, which is extensively studied in the statistics and machine learning literature. For example, see the classic work of DGL:1996 among many others. In this literature, the empirical risk minimization (ERM) classifier over the class of binary classifiers specified by ((ref)) is defined as a minimizer of the empirical predictive risk, which is taken to be one minus the objective function of the maximum score prediction problem. In other words, the ERM classification rule is identical to the maximum score prediction rule.
In this paper, we address the covariate selection issue in the framework of predicting the binary outcome $Y_{i}$ using the class of linear threshold-crossing prediction rules $b_{\theta }(W_{i})$ defined by ((ref)). We study the best subset selection procedure by which the covariates are chosen among a collection of candidate explanatory variables by maximizing the empirical score subject to a constraint on the maximal number of selected variables. In other words, we investigate theoretical and numerical properties of the $\ell _{0}$-norm constrained maximum score prediction rules.\footnote{ Here, the $\ell _{0}$-norm of a real vector refers to the number of non-zero components of the vector.}
To the best of our knowledge, Greenshtein2006 and Jiang:10 are the only existing papers in the literature that explicitly considered the same prediction problem as ours. Greenshtein2006 considered a general loss function that includes maximum score prediction as a special case in the i.i.d.\thinspace setup. Greenshtein2006 focused on a high dimensional case and established conditions under which the excess risk converges to zero as $n\rightarrow \infty $. Jiang:10 focused on the prediction of time series data and obtained an upper bound for the excess risk. Neither Greenshtein2006 nor Jiang:10 provided any numerical results for the best subset maximum score prediction rule. In contrast, we focus on cross-sectional applications and emphasize computational aspects.
The main contributions of this paper are twofold: first, we show that the best subset maximum score prediction rule is minimax rate-optimal and second, we demonstrate that it can be implemented via mixed integer optimization. The first contribution is theoretical and builds on the literature of empirical risk minimization tsybakov2004, massart2006. Specifically, we obtain non-asymptotic upper and lower risk bounds when the dimension of potential covariates is possibly much larger than the sample size $n$. Our upper and lower risk bounds are minimax rate-optimal when the maximal number of selected variables is fixed and does not increase with $n$. The existing results of finite-sample upper and lower risk bounds for the binary prediction problem focus on the case where there is no variable selection and the set of covariates is fixed and low-dimensional. Our risk bound results extend to the setup under the $\ell _{0}$-norm constraint when the set of potential covariates is high-dimensional.\footnote{ Raskutti2011 developed minimax rate results for high-dimensional linear mean regression models. We have used in the derivation of our lower risk bound a technical lemma of their paper Raskutti2011, which is based on the approximation theory literature. Nonetheless, our results are not directly obtainable from Raskutti2011, who considered the least squares objective function.}
The second contribution is computational. We face two kinds of computational challenges. One challenge comes from the nature of the objective function and the other is from the best subset selection. The maximum score objective function is a piecewise constant function whose range set contains only finitely many points. Hence the maximum of the score maximization problem is always attained yet the maximizer is generally not unique. It is known that computing the maximum score estimates regardless of the presence of the $ \ell _{0}$ constraint is NP (non-deterministic polynomial-time)-hard Johnson:78. See Manski:Thompson:1986 and PINKSE1993 for first generation algorithms for maximum score estimation.
Our computation algorithm is based on the method of mixed integer optimization (MIO). Florios:Skouras:08 provided compelling numerical evidence that the MIO approach is superior to the first-generation approaches. Kitagawa:Tetenov:2015 used an MIO formulation that is different from Florios:Skouras:08 to solve maximum score type problems. The objective of interest in Kitagawa:Tetenov:2015 is to develop treatment choice rules by maximizing an empirical welfare criterion, which resembles the maximum score objective function. They derived minimax optimality and used the MIO formulation to implement their algorithm. Neither Florios:Skouras:08 nor Kitagawa:Tetenov:2015 were concerned with the variable selection problem.
These second generation approaches are driven by developments in MIO solvers and also by availability of a much faster computer compared to the period when the first generation algorithms were proposed. Florios:Skouras:08 reported that they obtained the exact maximum score estimates using horowitz1993's data in 10.5 hours. In this application, the sample size was $n=842$ and there were 4 parameters to estimate.
It is well known that use of a good and tighter parameter space can strengthen the performance of a global optimization procedure including the MIO approach. In this paper, we propose a data driven approach to refine the parameter space. Using a state-of-the-art MIO solver as well as a tailor-made heuristic to choose the parameter space, it took us less than 5 minutes to obtain the exact maximum score estimates using the same dataset with the same number of parameters to estimate.\footnote{ This numerical result can be found in Online Appendix (ref) of the paper.} This is a dramatic improvement at the factor of more than 100 relative to the numerical performance reported in Florios:Skouras:08 . In other words, we demonstrate that hardware improvements combined with the advances in MIO solvers and also with a carefully chosen parameter space have made the maximum score approach empirically much more relevant now than ten years ago.
The second numerical challenge is concerned with constrained optimization with the $\ell _{0}$-norm constraint. It is well known that the $\ell _{0}$ -norm constraint renders the variable selection problem NP-hard even in the regression setup where the objective function is convex and smooth Natarajan:95,bertsimas2016. Recently bertsimas2016 proposed a novel MIO approach to the best subset variable selection problem when least squares and least absolute deviation risks are concerned. They demonstrated that the MIO approach can efficiently deliver a provably optimal solution to the resulting $\ell _{0}$-norm constrained risk minimization problem for a variety of datasets with practical problem size. Our implementation of the best subset maximum score prediction rules combines insights from bertsimas2016, Florios:Skouras:08, and Kitagawa:Tetenov:2015. We present two alternative MIO solution methods that complement each other.
In practical applications, it is useful to consider an approximate solution by adopting an early termination rule. In our empirical application, by setting an explicitly pre-specified optimization error, we were able to obtain approximate maximum score estimates with horowitz1993's data in around 10 minutes when both an intercept term and one specific random covariate were always selected, and there were 9 additional auxiliary covariates that were subject to the constraint where at most 5 of them could be selected. This suggests that fast developments in computing environments will enable us to solve an empirical problem at a practically relevant scale in very near future. We provide additional numerical evidence in Monte Carlo experiments in a high-dimensional setup when the number of potential covariates is larger than the sample size.
The remainder of this paper is organized as follows. In Section (ref), we describe our prediction rule. Section (ref) establishes theoretical properties of the proposed prediction rule. In Section (ref), we present computation algorithms using the MIO approach, and in Section (ref), we conduct a simulation study on the performance of our prediction rule in both low and high dimensional variable selection problems. In Section (ref), we illustrate usefulness of our prediction rule in the empirical application of work-trip mode choice using horowitz1993's data. We then conclude the paper in Section (ref). Proofs of all theoretical results and supplementary material of this paper are collated in online appendices.
In this section, we describe our proposal of the best subset maximum score prediction rule. Following MD:99 and Danilov2004, we distinguish between focus covariates that are always included in the prediction rule and auxiliary covariates of which we are less certain. We thus decompose the covariate vector $W$ as $W=(X,Z)$, where $X$ is a $\left( k+1\right) $-dimensional vector of focus covariates and $Z$ is a $p$-dimensional vector of auxiliary covariates.
Noting that $b_{a\theta }(w)=$ $b_{\theta }(w)$ for any positive real scalar $a$, we adopt the same scale normalization method as in Horowitz:92 and Jiang:10 by restricting the magnitude of the coefficient of one of the focus covariates to be unity. Specifically, write $X=(X_{0}, \widetilde{X})$ where $X_{0}$ is a scalar variable and $\widetilde{X}$ is the remaining $k$-dimensional subvector of focus covariates. The parameter vector $\theta $ in ((ref)) is decomposed accordingly as $\theta =(\alpha ,\beta ,\gamma )$, where $\alpha \in \{-1,1\}$ and $\left( \beta ,\gamma \right) \in \Theta $, which is a subset of $\mathbb{R}^{k+p}$. In this notation, the binary prediction rule has the following form:
We consider a parsimonious variable selection method by which the constituted prediction rule does not include more than a pre-specified number of auxiliary covariates. For any $p$ dimensional real vector $c$, let $\left\Vert c\right\Vert _{0}\equiv\sum\nolimits_{j=1}^{p}1\{c_{j}\neq0\}$ be the $\ell_{0}$-norm of $c$. We carry out the $\ell_{0}$-norm constrained covariate selection procedure by solving the constrained maximization problem
where the objective function $S_{n}$ is defined as
and the $\ell_0$-norm constrained parameter space is given as
for a given positive integer $q$.
As we discussed in the introduction, solving for the exact maximizer for ( (ref)) is desirable yet can be computationally challenging. It is hence practically useful to consider an approximate solution, which is constructed below, to the maximization problem ((ref)).
For any $\varepsilon \geq 0$, let $(\widehat{\alpha },\widehat{\beta }, \widehat{\gamma })\in \{-1,1\}\times \Theta _{q}$ be an approximate maximizer with $\varepsilon $ tolerance level such that
We refer to the prediction rule defined by $1\{\widehat{\alpha }x_{0}+ \widetilde{x}^{\prime }\widehat{\beta }+z^{\prime }\widehat{\gamma }\geq 0\}$ as the approximate best subset maximum score binary prediction rule. \footnote{ The dependence of $(\widehat{\alpha },\widehat{\beta },\widehat{\gamma })$ on $\varepsilon $ is suppressed for simplicity of notation.} The value of $ \varepsilon $ can be specified for early termination of the solution algorithm to the problem ((ref)). In Section (ref), we will present an algorithm that allows for computing an approximate solution to ((ref)) within a definite approximation error bound specified by the tolerance level $\varepsilon $. In what follows, we use PRESCIENCE as shorthand for the approximate best subset maximum score binary prediction rule.\footnote{ It comes from the aPpRoximate bEst S(C)ubset maxImum scorE biNary prediCtion rulE.}
In this section, we study the theoretical properties of PRESCIENCE. Let $F$ denote the joint distribution of $\left( Y,W\right) $. For $\left( \alpha ,\beta ,\gamma \right) \in \{-1,1\}\times \Theta $, let
Note that $S(\alpha ,\beta ,\gamma )$ depends on the joint distribution $F$. Given a cardinality bound $q$, let
That is, $S_{q}^{\ast }$ is the supremum of $S(\alpha ,\beta ,\gamma )$ given the $\ell _{0}$-norm constraint.
Following the literature on empirical risk minimization (see, e.g., DGL:1996, Lugosi:2002, tsybakov2004, massart2006, Greenshtein2006 and Jiang:10 among many others), we assess the predictive performance of PRESCIENCE by bounding the difference
where $(\widehat{\alpha },\widehat{\beta },\widehat{\gamma })$ is defined by ((ref)). The difference $U_{n}$ is non-negative by the definition of $S_{q}^{\ast }$. Hence, a good prediction rule will result in a small value of $U_{n}$ with a high probability and also on average.
Throughout this section, we assume that $\varepsilon =0$ for simplicity. Before presenting our theoretical results, we first introduce some notation. For any two real numbers $a$ and $b$, let $a\vee b\equiv \max \{a,b\}$. Let $ s\equiv k+q$ and
Theorem (ref) establishes that the tail probability of $ U_{n}$ decays exponentially in $r_{n}$. Moreover, this result is non-asymptotic: inequality ((ref)) is valid for every sample size $n$ for which condition ((ref)) holds. By comparing the leading terms on both sides of inequality ((ref)), we can see that condition ((ref)) is satisfied, for instance, if
Hence, condition ((ref)) is satisfied easily when $p\vee n$ takes a relatively large value compared to $s$. If $(p\vee n)$ diverges to infinity, then $s$ can diverge at a sufficiently slow rate.
Theorem (ref) implies that
provided that
holds. This allows the case that
In other words, the predictive performance of PRESCIENCE remains good even when the number of potentially relevant covariates ($p$) grows exponentially, provided that the number of selected covariates ($s$) can only grow at a polynomial rate. GreenshteinRitov2004 and Greenshtein2006 consider the case where $p$ grows at a polynomial rate. For this case, condition ((ref)) implies that $s=o\left( n/\ln n\right) $, which coincides with the optimal sparsity rate established by GreenshteinRitov2004 and Greenshtein2006 under which a sequence of predictor selection procedures subject to the sparsity constraint can be shown to be persistent.
The result ((ref)) is derived under the i.i.d.\thinspace setup but does not hinge on other regularity conditions on the underlying data generating distribution $F$. This rate result can be improved under additional assumptions on the distribution $F$. In this section, we consider a condition that is called the margin condition in the literature under which we may obtain a sharper result on the upper bound of $E\left[ U_{n}\right] $. As before, the derived bound will be non-asymptotic.
It is necessary to introduce additional notation. Let
That is, $\mathcal{B}_{q}$ is the class of all prediction rules in (ref) with the $\ell _{0}$-norm constraint. For $w\in \mathcal{W}$, let
For any measurable function $f:\mathcal{W\mapsto }\mathbb{R}$, let $ \left\Vert f\right\Vert _{1}=E\left[ \left\vert f(W)\right\vert \right] $ denote the $L_{1}$-norm of $f$. The functions $\eta $ and $b^{\ast }$ as well as the $L_{1}$-norm $\left\Vert \cdot \right\Vert _{1}$ depend on the data generating distribution $F$. For any indicator function $b:\mathcal{W} \mapsto \left\{ 0,1\right\} $, let
We now state the following regularity condition.
Condition (ref) is termed as the margin condition in the literature (see, e.g., mammen1999, tsybakov2004 and massart2006). For any binary predictor $b$,
so that $\widetilde{S}(b)$ is maximized at $b=b^{\ast }$. Hence, Condition (ref) implies that the functional $ \widetilde{S}\left( \cdot \right) $ has a well-separated maximum. Suppose that there exist universal positive constants $C$ and $\alpha $ such that
for all $t>0$. Then by modifying the proof of Proposition 1 of tsybakov2004 slightly, we can show that (ref) holds with $\vartheta =(1+\alpha )/\alpha $. See tsybakov2004 for further discussions on the margin condition.
Recall that it is not necessary to assume (ref) to establish the risk consistency, as shown in Theorem (ref) . We show below that we can obtain a tighter upper bound on $E\left[ U_{n} \right] $ under (ref). Let
The next theorem, which is an application of massart2006, establishes a finite-sample bound on $E\left[ U_{n}\right] $ under the margin condition.
For $p\vee n$ sufficiently large, we have that $\rho _{n}\leq 5s\ln (p\vee n) $; thus, inequality ((ref)) can hold under condition ((ref)) in large samples, provided that $h$ is fixed or does not go to zero too rapidly.
The first term on the right-hand side of inequality (ref) represents the bias term. Equation ((ref)) implies that there is no bias term, namely $\widetilde{S}\left( b^{\ast }\right) =\sup\nolimits_{b\in \mathcal{B}_{q}}\widetilde{S}\left( b\right) $ if $b^{\ast }\in \mathcal{B}_{q}$. Therefore, Theorem (ref) implies that
provided that $b^{\ast }\in \mathcal{B}_{q}$.\footnote{ For the case with $\varepsilon >0$, it is also straightforward to modify Theorem (ref) such that the rate result ((ref) ) continues to hold provided that
} The rate of convergence in (ref) doubles that in (ref) when $h$ is fixed and $\vartheta =1$. We notice that, if $b^{\ast }\notin \mathcal{B}_{q}$, the upper bound derived in Theorem (ref) would asymptotically reduce to the non-zero bias term $\widetilde{S}\left( b^{\ast }\right) -\sup\nolimits_{b\in \mathcal{B}_{q}}\widetilde{S}\left( b\right) $ and hence the margin condition alone does not suffice for deducing there is improved rate of convergence. Nevertheless, the rate result (ref) still holds regardless of the validity of the presumption that $b^{\ast }\in \mathcal{B}_{q}$.
We now remark on the condition that $b^{\ast }\in \mathcal{B}_{q}$ in the context of the binary response model specified below. Suppose that the outcome $Y$ is generated from a latent variable threshold crossing model Manski:1975, Manski:1985:
where $\theta ^{\ast }$ denotes the true data generating parameter vector and $\xi $ is an unobserved latent variable whose distribution satisfies that
Let $\theta _{0}\equiv \arg \sup\nolimits_{\theta \in \{-1,1\}\times \Theta _{q}}$ $S(\theta )$. For simplicity, assume that $\theta _{0}\in \{-1,1\}\times \Theta _{q}$ so that the maximum is attained.
Manski:1988 showed that, for the binary response model specified by ((ref)) and ((ref)), the true parameter value $\theta ^{\ast }$ is identified relative to another value $\theta $ if and only if the event that $W^{\prime }\theta $ and $W^{\prime }\theta ^{\ast }$ have different sign occurs with positive probability. Therefore, Proposition (ref) implies that $b^{\ast }\in \mathcal{B}_{q}$ if and only if the “pseudo-true” value $\theta _{0}$ is observationally equivalent to $\theta ^{\ast }$. In particular, this implies that $\theta ^{\ast }=\theta _{0}$ if $\theta ^{\ast }$ is point-identified.
It would be interesting to study the role of the bias when $b^{\ast }\notin \mathcal{B}_{q}$ using the framework of sieve estimation Chen:2007. As pointed by Elliott2013, what matters is how well we can approximate the value of optimum $\sup\nolimits_{\theta \in \{-1,1\}\times \Theta _{q}}S(\theta )$, not the optimizer $\arg \sup\nolimits_{\theta \in \{-1,1\}\times \Theta _{q}}S(\theta )$. However, it would be much more demanding to develop non-asymptotic theory when the bias is present in our framework. We leave this as a topic for future research.
In this section, we derive a minimax lower bound under the margin condition. In particular, we focus on the case that $s=k+q$ is low-dimensional in that $ s$ does not grow with sample size $n$ and also consider a sufficient condition for the margin condition.
Condition (ref) is satisfied with $ \vartheta =1$ whenever
massart2006 introduced ((ref)) as an easily interpretable margin condition requiring that the conditional probability $ \eta (w)$ should be bounded away from $1/2$. Condition ((ref)) holds under certain regularity assumptions on the binary response model as indicated in the following proposition.
Conditions (i) and (ii) in Proposition (ref) assume that $\left\vert W^{\prime }\theta ^{\ast }\right\vert $ is bounded away from zero and the density of $\xi $ conditional on $W=w$ is bounded away from zero in a neighborhood of zero. While the latter condition is mild, the former is non-trivial. Condition (i) can hold easily when all components of $W$ are discrete, which is not uncommon in microeconometric applications of binary response models (see e.g., komarova2013 and magnac2008). In the presence of continuous covariates, this condition becomes more restrictive.
For any real vector $u$, let $\left\Vert u\right\Vert _{E}=\sqrt{u^{\prime}u} $ denote the Euclidean norm of $u$. To state a minimax lower bound, we first define the following class of distributions.
The first two conditions in the definition of $\mathcal{P}(h,\mathcal{B} _{q}) $ have already been introduced before. The new condition (iii) imposes that the Euclidean norm $\left\Vert \theta -\widetilde{\theta }\right\Vert _{E}$ is equivalent to the $L_{1}$-norm $\left\Vert b_{\theta }-b_{ \widetilde{\theta }}\right\Vert _{1}$ for two values $\theta $ and $ \widetilde{\theta }$ that differ only in the components corresponding to the auxiliary covariate coefficients. This condition is concerned with restrictions on the distribution of the covariate vector $W$. The following proposition gives sufficient conditions for verifying this norm equivalence condition.
For any subset $J\subset \{1,...,p\}$, let $Z_{J}$ denote the $\left\vert J\right\vert $-dimensional subvector of $Z\equiv (Z^{(1)},\ldots ,Z^{(p)})^{\prime }$ formed by keeping only those elements $Z^{(j)}$ with $ j\in J$. Let $\mathcal{I}_{q}\equiv \cup _{\left( \beta ,\gamma \right) \in \Theta _{q}}\text{Supp}(\widetilde{X}\beta +Z^{\prime }\gamma )$, where $ \text{Supp}(V)$ denotes the support of the random variable $V$.
Condition (a) in Proposition (ref) is mild. The first part of condition (b) holds with $L_{2}=\overline{L}\sqrt{ 2q}$ if $\max_{j\in \{1,...,p\}}\left\vert Z^{(j)}\right\vert \leq \overline{ L}$ with probability 1 for some universal positive constant $\overline{L}$. The second part of condition (b) is related to the sparse eigenvalue assumption used in the high dimensional regression literature (see, e.g. Raskutti2011). For example, suppose that $Z$ is a random vector with mean zero and the covariance matrix $\Sigma $ whose $(i,j)$ component is $ \Sigma _{i,j}=\rho ^{|i-j|}$ for some constant $\rho >0$. Then the smallest eigenvalue of $\Sigma $ is bounded away from zero where the lower bound is independent of the dimension $p$ (van2009). Thus, in this case, $E\left( Z_{J}Z_{J}^{\prime }\right) $ is bounded below by that same lower bound.
We now state the result on the minimax lower bound for the predictive performance of PRESCIENCE.
For any estimator $\widehat{b}$ taking value in $\mathcal{B}_{q}$, Theorem (ref) implies that, as long as $\mathcal{P}(h,\mathcal{B} _{q})$ is non-empty, there is some distribution $F$ under which the average predictive risk $E\left[ S_{q}^{\ast }-\widetilde{S}(\widehat{b})\right] $ cannot be smaller than the lower bound term stated in (ref). Comparing the upper and lower bounds given by ((ref)) and ((ref)), we can deduce conditions under which these two bounds coincides in terms of rate of convergence such that the PRESCIENCE approach is rate-optimal in the minimax sense. Suppose that $ (s,c_{l},c_{u})$ are fixed and does not increase or decrease with $n$. Then the risk lower bound is of order
Comparing (ref) to (ref) evaluated at $\vartheta =1$, we see that, if $h$ is also a universal constant and $p$ grows at a polynomial or exponential rate in $n$, then the upper and lower bound results induce the same convergence rate and hence the PRESCIENCE approach is minimax rate-optimal. On the other hand, when Condition (ref) holds with $\vartheta >1$, the rate given by ((ref)) is slower than that given by ((ref)) such that the convergence rate implied by the risk lower bound need not be attained and therefore the PRESCIENCE approach may not be rate-optimal.
We now present algorithms for solving the maximization problem ((ref)). It is straightforward to see that solving ((ref)) is the same as solving
In what follows, we focus on solving the sub-problem
because the other case corresponding to $\alpha =-1$ can be solved by replacing the value of $X_{0i}$ with that of $-X_{0i}$ and then applying the same solution method as developed for the case ((ref)).
By ((ref)) and noting that $Y_{i}\in\{0,1\}$, solving the problem ((ref)) amounts to solving
We assume that the parameter space $\Theta$ is bounded and takes the polyhedral form:
for some real constant matrices $A_{1}$ and $A_{2}$ and some real constant vector $B$. Let
denote the smallest cube containing all values of $\gamma$ in the pair $ \left( \beta,\gamma\right) $ confined by $\Theta$. Writing $\gamma =(\gamma_{1},...,\gamma_{p})$, we have that, if $\left( \beta,\gamma\right) \in\Theta$, then $-\infty<\underline{\gamma}_{j}\leq\gamma_{j}\leq \overline{ \gamma}_{j}<\infty$ for $j\in\{1,...,p\}$. Let
Our implementation builds on the method of mixed integer optimization (in particular, bertsimas2016, Florios:Skouras:08, and Kitagawa:Tetenov:2015) and present two alternative solution methods that complement each other. The values $\left( M_{i}\right) _{i=1}^{n}$ can be computed by formulating the maximization problem in ((ref)) as linear programming problems, which can be easily and efficiently solved by modern numerical software. Hence these values can be computed and stored beforehand as inputs to the algorithms that are used to solve the MIO problems described below.
Our first solution method is based on an equivalent reformulation of the maximization problem ((ref)) as the following constrained mixed integer optimization (MIO) problem:
where $\delta $ is a given small and positive real scalar (e.g. $\delta =10^{-6}$ as in our numerical study).
We now explain the equivalence between ((ref)) and ((ref)). Given $(\beta ,\gamma )$, the inequality constraints ( (ref)) and the dichotomization constraints ((ref)) enforce that $d_{i}=1\{X_{0i}+\widetilde{X}_{i}^{\prime }\beta +Z_{i}^{\prime }\gamma \geq 0\}$ for $i\in \{1,...,n\}$. Therefore, maximizing the objective function in ((ref)) for $(\beta ,\gamma )\in \Theta $ subject to the constraints ((ref)) and ((ref)) is equivalent to solving the problem ((ref)) using all covariates. This part of formulation is similar to the MIO formulation used by Kitagawa:Tetenov:2015 for solving the maximum score type estimation problems without the variable selection constraint.
Following bertsimas2016, we implement the best subset variable selection feature through the additional constraints ((ref)), ((ref)) and ((ref)). The on-off constraints ((ref)) and ((ref)) ensure that, whenever $e_{j}=0$, the auxiliary covariate $Z_{j}$ is excluded in the resulting PRESCIENCE. Finally, the cardinality constraint $\left\Vert \gamma\right\Vert _{0}\leq q$ is enacted through the constraint ((ref)), which restricts the maximal number of the binary controls $e_{j}$ that can take value unity.
Modern numerical optimization solvers such as CPLEX, Gurobi, MOPS, Mosek and Xpress-MP can be used to effectively solve the MIO formulations of the PRESCIENCE problem. Most of the solution algorithms employed by the MIO solvers can be viewed as complex and advanced refinements of the well-known branch-and-bound method for solving MIO problems.\footnote{ See Online Appendix (ref) for further details of the branch-and-bound method.} Along the branch-and-bound solution process, we can keep track of two important values: the best upper and lower bounds on the objective value of the MIO problem ((ref)). The best lower bound corresponds to the objective function evaluated at the incumbent solution, which is the best feasible solution discovered so far. The best upper bound can be deduced by taking the maximum of the optimal objective values of all the linear programming relaxation formulations of the branching MIO sub-problems that have been solved so far. Let $MIO\_gap$ denote the difference between these two bounds. Note that the incumbent solution becomes optimal when the $MIO\_gap$ value reduces to zero.
We can use the $MIO\_gap$ value to solve for the $\varepsilon $-level PRESCIENCE introduced in Section (ref). To see this, consider an early termination rule by which the solution algorithm is terminated whenever $MIO\_gap\leq \varepsilon $ where $\varepsilon $ is a given tolerance level. Let $(\widehat{\beta },\widehat{\gamma })$ be the incumbent solution upon termination of the MIO solver. Because $(\widehat{ \beta },\widehat{\gamma })$ is in the feasible solution set of the problem ( (ref)), by constraints ((ref)) and ((ref)), we have that $\left\Vert \widehat{\gamma } \right\Vert _{0}\leq q$ . Moreover, by constraints ((ref)) and ((ref)), we have that $\widehat{d}_{i}=1\{X_{0i}+\widetilde{ X}_{i}^{\prime }\widehat{\beta }+Z_{i}^{\prime }\widehat{\gamma }\geq 0\}$ for $i\in \{1,...,n\}$ so that $S_{n}(1,\widehat{\beta },\widehat{\gamma })$ is equal to the objective function in ((ref)) evaluated at $( \widehat{d}_{1},...,\widehat{d}_{n})$. Since ((ref)) and ((ref)) are equivalent maximization problems, it thus follows from the construction of $MIO\_gap$ value that
Given the termination condition, we can therefore see that
which yields an approximately optimal solution with the optimization tolerance level $\varepsilon $\ for the problem ((ref)).
The PRESCIENCE can also be solved by incorporating the constraints ((ref)), ((ref)) and ((ref)) in the MIO formulation of Florios:Skouras:08 for solving the maximum score estimation problem. We now present this alternative MIO formulation below.
Consider the constrained maximization problem:
The problem ((ref)) without the constraint $\left\Vert \gamma \right\Vert _{0}\leq q$ reduces to the type of maximum score estimation problem studied by Florios:Skouras:08. The objective function in ((ref)) coincides with that in ((ref)) with probability 1 as long as the sum $X_{0}+\widetilde{X}^{\prime }\beta +Z^{\prime }\gamma $ is continuously distributed. This condition holds provided that the distribution of $X_{0}$ conditional on $(\widetilde{X},Z)$ is continuous. With such a continuous covariate, we can also solve ((ref)) by solving the following MIO formulation of ((ref)):
Florios:Skouras:08 showed that maximizing the objective function in ( (ref)) for $(\beta ,\gamma )\in \Theta $ subject to the constraints ((ref)) and ((ref)) is equivalent to solving the problem ((ref)) using all covariates. This can be seen from the fact that the objective function of the MIO problem ((ref)) is strictly increasing in $d_{i}$ so that, given $(\beta ,\gamma )$, it is optimal to set $d_{i}=1\{\left( 2Y_{i}-1\right) (X_{0i}+\widetilde{X} _{i}^{\prime }\beta +Z_{i}^{\prime }\gamma )\geq 0\}$ under the constraints ( (ref)) and ((ref)). Along similar arguments to those discussed for the problem ((ref)), it is also straightforward to verify that the variable selection constraint $ \left\Vert \gamma \right\Vert _{0}\leq q$ is imposed through the constraints ((ref)), ((ref)) and ((ref)). Therefore, the maximization problems ((ref)) and ((ref)) are equivalent.
The MIO formulations ((ref)) and ((ref)) depend on the specification of the parameter space $\Theta$. It is well known that use of a good and tighter parameter space can strengthen the performance of a global optimization procedure. Given an initial specification of $\Theta $, we propose below a data driven approach to refine the parameter space.
Recall that $W_{i}=(X_{0i},\widetilde{X}_{i},Z_{i})$ is the entire covariate vector. Let $\widetilde{W}_{i}$ denote the vector $(\widetilde{X}_{i},Z_{i})$ . For $i\in \{1,...,n\}$, let $\widehat{P}_{i}$ be an estimate of ${P} _{i}\equiv P(Y_{i}=1|W_{i})$. Define the following sets recursively:
and, for $m\in \{2,...,k+p\}$,
where, for $j\in \{1,...,k+p\}$, the quantities $\widehat{l}_{j}$ and $ \widehat{u}_{j}$ are defined respectively by
If the binary outcome $Y_{i}$ is generated from the model specified by ((ref)) and ((ref)) and the conditional probability $P(Y=1|W)$ is nonparametrically estimated, the interval $[ \widehat{l}_{j},\widehat{u}_{j}]$ is a nonparametric estimate of the identified set for the $j$th component of the parameter vector $t=(\beta ,\gamma )$. In this case, the sign-matching constraints ((ref)) can be regarded as the empirical counterparts of the inequalities stated in the set
which contains those $t$ values that are observationally equivalent to the true data generating parameter value komarova2013, chen2015. Our procedures for computing $\widehat{l}_{j}$ and $\widehat{u} _{j}$ are modified versions of horowitz1998's linear programming formulations of the identified bounds on the parameter components. The formulations ((ref)) and ((ref)) differ from those of Horowitz in that we further tighten the domain of $t$ in these optimization problems by exploiting the information of the upper and lower bound values that have been solved so far.
When the covariate vector is of high dimension, nonparametric estimation of $ P(Y=1|W)$ would suffer from the curse of dimensionality problem. In this paper, we consider estimating this conditional probability by parametric methods such as the logit or probit approach. In Monte Carlo experiments and an empirical application, we estimate $({P}_{i})_{i=1}^{n}$ by the fitted choice probabilities from the logit regression of $Y$ on all the covariates.
Noting that the parametric model for estimating $P(Y=1|W)$ may be misspecified, for $\tau \geq 1$, we construct a conservative space $\widehat{ \Theta }\left( \tau \right) $, which is a $\tau $-enlargement of the space $ \prod\nolimits_{j=1}^{k+p}\left[ \widehat{l}_{j},\widehat{u}_{j}\right] $ as given below:
We can solve the MIO problems ((ref)) and ((ref)) with the refined parameter space $\widehat{\Theta } \left( \tau \right) $ in place of the original space $\Theta $.\footnote{ Both the covering cube $\mathbf{C}$ and the quantities $\left( M_{i}\right) _{i=1}^{n}$ depend on the input parameter space. Hence, for the warm-start formulations of ((ref)) and ((ref)), these objects are also computed under the refined space $\widehat{\Theta } \left( \tau \right) $.} Using the terminology used in bertsimas2016, we shall refer to these refined MIO representations as the warm-start MIO formulations of the PRESCIENCE problem. The value of $\tau $ is treated as a tuning parameter for solving the warm-start MIO problems. The original formulations ((ref)) and ((ref)) based on the space $\Theta $ are referred to as the cold-start MIO formulations.
Computation of the refined space $\widehat{\Theta }\left( \tau \right) $ requires solving $2\left( k+p\right) $ simple linear programming problems. This task can be done very efficiently even when $p$ is relatively large. On the other hand, the space $\widehat{\Theta }\left( \tau \right) $ is not always constructible since the problems ((ref)) and ((ref)) may not admit any feasible solution. This may occur due to the misspecification issue of using parametric choice probability estimates. Alternatively, it can also occur when the postulated binary response model specified by ((ref)) and ((ref)) itself is misspecified. As illustrated by Monte Carlo simulations and a real data application in Online Appendices (ref) and (ref), when the refined space $\widehat{ \Theta }\left( \tau \right) $ is available, solving the warm-start MIO formulations can be computationally far more efficient than solving their corresponding cold-start versions.
We conclude this subsection by commenting that our warm-start approach does not work well when $p$ is greater than $n$. In this case, irrespective of the knowledge of the true choice probabilities, the dimension of the vector of unknown coefficients is larger than the number of inequalities given by the constraints ((ref)) such that these constraints may become ineffective for tightening the original parameter bounds. It is a topic for future research how to devise a good wart-start option for the high-dimensional setup.
In this section, we study the performance of the PRESCIENCE method in Monte Carlo experiments. Throughout this paper, we used the MATLAB implementation of the Gurobi Optimizer to solve the MIO problems. Moreover, all numerical computations were done on a desktop PC (Windows 7) equipped with 32 GB RAM and a CPU processor (Intel i7-5930K) of 3.5 GHz.\footnote{ The MATLAB codes for implementing the PRESCIENCE approach are available from the authors via the website \url{https://github.com/LeyuChen/Best-Subset-Binary-Prediction}. This implementation requires the Gurobi Optimizer, which is freely available for academic purposes.}
Let $V=(V_{1},...,V_{p+1})$ be a multivariate normal random vector with mean zero and covariance matrix $\Sigma $ with its element $\Sigma _{i,j}=\left( 0.25\right) ^{\left\vert i-j\right\vert }$. The binary outcome is generated according to the following setup:
where $\theta ^{\ast }$ denotes the value of the true data generating parameter vector, $W=(X,Z)$ is a $(p+2)$ dimensional covariate vector with the focus covariates $X=(X_{0},\widetilde{X})=(V_{1},1)$ and the auxiliary covariates $Z=(V_{2},...,V_{p+1})$, and $\xi $ is a $N(0,1)$ random variate independent of $V$. We set $\theta _{1}^{\ast }=1$, $\theta _{2}^{\ast }=0$, and $\theta _{j}^{\ast }=0$ for $j\in \{4,...,p+2\}.$ The coefficient $ \theta _{3}^{\ast }$ is chosen to be non-zero such that, among the $p$ auxiliary covariates, only the variable $Z_{1}$ is relevant in the data generating processes (DGP).
We consider the following two specifications for $\theta _{3}^{\ast }$ and $ \sigma (W)$:
As before, the parameter vector $\theta $ in ((ref)) is decomposed as $\theta =(\alpha ,\beta ,\gamma )$ where $\alpha $, $\beta $ and $\gamma $ are coefficients associated with $X_{0}$, $\widetilde{X}$ and $Z$, respectively. The parameter space for the PRESCIENCE approach is specified to be
over which we compute PRESCIENCE via solving its corresponding MIO problem.
There were $100$ simulation repetitions in each Monte Carlo experiment. For each simulation repetition, we generated a training sample of $n$ observations for estimating the coefficients $\theta $ and a validation sample of $5000$ observations for evaluating the out-of-sample predictive performance. The training sample size $n$ was set to be 100 for DGP(i) and 50 for DGP(ii). For each DGP setup, we performed simulations with both the low and high dimensional covariate configurations. For the low dimensional case, we set $p=10$ for both DGP(i) and (ii). For the high dimensional case, we set $p=200$ for DGP(i) and $p=60$ for DGP(ii).
We considered the following class of prediction methods:
where PRESCIENCE$\left( q\right) $ denotes the PRESCIENCE approach with a cardinality bound $q$ imposed on the auxiliary covariates, PRE_CV denotes the PRESCIENCE approach using a data driven value of $q\in \{1,2,3\}$ via the 5-fold cross validation procedure, and logit_lasso and probit_lasso respectively denote the $\ell _{1}$-penalized logit and probit maximum likelihood estimation (MLE) approaches friedman2010. Throughout this simulation study, we employed the cold-start MIO formulation ((ref)) to solve the PRESCIENCE problems. For the simulation experiment with $p<n$, we computed the exact solution to each PRESCIENCE problem. For the high dimensional case with $p>n$ , we solved for the PRESCIENCE solution with the tolerance level $ \varepsilon $ specified according to the rule
Note that this early termination rule is compatible with the order of magnitude stated in the condition ((ref)) for the convergence rate result (ref). For the logit_lasso and probit_lasso approaches, we used the MATLAB function lassoglm to implement these two penalized MLE approaches for which we calibrated the lasso penalty parameter value over a sequence of 100 values via the 10-fold cross validation procedure. We used the default setup of lassoglm for constructing this tuning sequence among which we made the following three choices, $\left\{ \lambda _{\min },\lambda _{1se},\lambda _{2se}\right\} $, of the penalty parameter value. To describe these, let $\{I_{j}:j=1,\ldots ,K\}$ be the partition of data, where $K=10$, and let $\hat{L}(I_{j},\lambda )$ denote the minus log likelihood function evaluated using data in $I_{j}$ but estimating the model using data in $\bigcup\nolimits_{i\in \{1,...,K\}\backslash \{j\}}I_{i}$ with a given penalty parameter value $ \lambda $. The value $\lambda _{\min }$ refers to the $\lambda $ value that minimized the mean cross validated deviances ($K^{-1}\sum_{j=1}^{K}\hat{L} (I_{j},\lambda )$), whereas $\lambda _{1se}$ and $\lambda _{2se}$ respectively denote the largest penalty parameter values whose corresponding mean cross validated deviances still fall within the one- and two-standard errors of $K^{-1}\sum_{j=1}^{K}\hat{L}(I_{j},\lambda _{\min })$.\footnote{ Here, standard errors are computed over the 10 cross-validation folds. Choice of the lasso tuning parameter based on $\lambda _{1se}$ is also known as the "one-standard-error" rule, which is commonly employed in the statistical learning literature Hastie2009.}
For each $m\in $ $\mathcal{M}$, let $\widehat{\theta }(m)$ denote the coefficients computed under the prediction method $m$. Let $in\_Score$ denote the average of the in-sample objective values $S_{n}(\widehat{\theta } (m))$ over all the simulation repetitions. In each simulation repetition, we approximated the out-of-sample objective value $S(\widehat{\theta }(m))$ using the generated validation sample. Let $out\_Score$ denote the average of $S(\widehat{\theta }(m))$ over all the simulation repetitions. It is straightforward to see that the theoretically best prediction rule $b^{\ast }(w)$ in this simulation design takes the form $b^{\ast }(w)=1\left\{ w^{\prime }\theta ^{\ast }\geq 0\right\} $. Hence, we also assess the predictive performance of a given prediction method $m\in \mathcal{M}$ by its relative score, which is ratio of the score evaluated at $\widehat{ \theta }(m)$ over that evaluated at $\theta ^{\ast }$. Let $in\_RS$ and $ out\_RS$ respectively denote the average of in-sample relative scores $S_{n}( \widehat{\theta }(m))/$ $S_{n}(\theta ^{\ast })$ and that of out-of-sample relative scores $S(\widehat{\theta }(m))/S(\theta ^{\ast })$ over all the simulation repetitions.
We also examine the variable selection performance of the prediction method. We say that a variable $Z_{j}$ is effectively selected under the prediction method $m$ if and only if the magnitude of $\widehat{\theta }_{j+2}(m)$ is larger than a small tolerance level (e.g. $10^{-6}$ as used in our numerical study) which is distinct from zero in numerical computation. Let $Corr\_sel$ be the proportion of the auxiliary covariate $Z_{1}$ being effectively selected. Let $Orac\_sel$ be the proportion of obtaining an oracle variable selection outcome where, among all the auxiliary covariates, $Z_{1}$ was the only one that was effectively selected. Let $Num\_irrel$ denote the average number of effectively selected auxiliary covariates whose true DGP coefficients are zero.
We now present the simulation results under the setup of DGP(i). First, we report the computational performance of our MIO solution algorithm to the PRESCIENCE problems. Table (ref) gives the summary statistics of the MIO computation time in CPU seconds across simulation repetitions. From this table, we can see that the MIO problems for the PRESCIENCE computation were solved very efficiently in the DGP(i) simulations where the number of the auxiliary covariates could be the double of the sample size yet the maximum computation time was only around 5 minutes. It is also interesting to note that the PRESCIENCE computation time was not monotone in $q$. This feature might be due to the branching strategy heuristics of the MIO branch-and-bound solution algorithms.
We next turn to the statistical performance of the binary prediction method. In Tables (ref) and (ref), we compare the aforementioned predictive and variable selection performance measures for the various prediction approaches given in ((ref)). As shown in these two tables, regardless of $p$, the in-sample fit in terms of $in\_Score$ and $in\_RS$ for the PRESCIENCE($q$) method increased with $q$. This finding is expected because the in-sample objective function ((ref)) is monotone in $q$ by design for the PRESCIENCE approach. Nonetheless, both tables indicate that $out\_Score$ and $out\_RS$ also declined as $q$ increased, thus resounding with the known issue that in-sample overfitting may result in poor out-of-sample performance. When the true number of effective auxiliary covariates is unknown, one can choose the $q$ value that maximizes the mean cross validated score. From Tables (ref) and (ref), we find that the PRE_CV approach indeed balanced well the in-sample and out-of-sample predictive performances. Moreover, its predictive performance measures were also comparable to those given by the logit_lasso and probit_lasso approaches.
We now discuss the variable selection results. Table (ref) indicates that all the prediction approaches had high $Corr\_sel$ rates and hence were capable of effectively selecting the relevant covariate $Z_{1}$. However, the good performance in the $Corr\_sel$ criterion may arise at the risk of overfitting. Therefore, we also have to take into account the performance in excluding irrelevant auxiliary covariates. The simulation design implies that the case with $q=1$ is the most parsimonious PRESCIENCE setup that correctly specifies the number of effective auxiliary covariates in the DGP. Therefore, it is not surprising that PRESCIENCE(1) performed the best in terms of $Num\_irrel$. We note that the PRE_CV approach also performed very well in excluding the irrelevant variables. In fact, for the PRESCIENCE based approaches, only PRESCIENCE(1) and PRE_CV could yield a non-zero probability of inducing an oracle variable selection outcome. For the penalized MLE approaches, the logit_lasso and probit_lasso coupled with the larger penalty parameter value $\lambda _{2se}$ also performed well in terms of $Orac\_Sel$ and $Num\_irrel$, albeit at the cost of a slight reduction of the out-of-sample predictive performances.
We also observe a similar pattern in the results of the DGP(i) setup with $ p=200$. From Table (ref), we find that the PRE_CV approach also balanced very well the requirement for including the relevant but excluding the irrelevant variables. It is also noted that the logit_lasso and probit_lasso approaches in this high dimensional simulation setup tended to select more irrelevant variables than the PRESCIENCE approach, hence suffering from a larger extent of overfitting.
To save space, we present details of the simulation results under the setup of DGP(ii) in Online Appendix (ref). The results are similar to those under DGP(i).
We illustrate usefulness of PRESCIENCE in the empirical application of work-trip mode choice. We used the transportation mode dataset analyzed by horowitz1993. This dataset has also been well studied for illustration of econometric methods developed for binary response models (e.g., see Florios:Skouras:08, benoit2012, and the references therein). The previous literature focused on estimating slope coefficients in the binary response model; however, in this section, we are mainly interested in the numerical performance of alternative MIO algorithms and the result of covariate selection.
The data consist of 842 observations sampled randomly from the Washington, D.C., area transportation study. Each record in the dataset contains the following information for a single work trip of the traveler: the chosen transportation mode, the number of cars owned by the traveler's household ($ CARS$), the transit out-of-vehicle travel time minus automobile out-of-vehicle travel time in minutes ($DOVTT$), the transit in-vehicle travel time minus automobile in-vehicle travel time in minutes ($DIVTT$) and the transit fare minus automobile travel cost in dollars ($DCOST$).
The dependent variable $Y$ is the traveler's chosen mode of transportation such that $Y=1$ if the choice is automobile and 0 otherwise. Following Florios:Skouras:08, we standardized each of explanatory variables to have mean zero and unit variance. Following horowitz1993 and Florios:Skouras:08, we specified the coefficient of $DCOST$ to be unity and did not estimate that parameter. We set the focus covariates $ X=(X_{0},\widetilde{X})$ to be $\left( DCOST,1\right) $, where the constant term was included to capture the regression intercept and the parameter $ \alpha $ was set to be unity. The resulting PRESCIENCE problem hence reduced to the maximization problem ((ref)). We implemented the two MIO formulations developed in Section (ref) for solving this problem. To compare their computational performance, we report the CPU time (in seconds) and the number of branch-and-bound nodes that the MIO solver had explored to reach the optimal solution. The former depends on both the computing hardware and software configurations whereas the latter only depends on the solution algorithms employed by the MIO solver.
For the auxiliary covariates, we set
The covariate specification ((ref)) allows us to approximate a smooth function of the variables $(CARS,DOVTT,DIVTT)$ by its quadratic expansion.\footnote{ In Online Appendix (ref), we report empirical results using $ Z=\left( CARS,DOVTT,DIVTT\right) $.} We are interested in the data driven selection of these expansion terms through the PRESCIENCE procedure. In this setup, we have that $k=1$ and $p=9$. We specified the parameter space $ \Theta $ to be $\left[ -10,10\right] ^{10}$. We computed the PRESCIENCE for each $q\in \{1,3,5\}$. Since there are 10 unknown parameters in this setup, we solved for the PRESCIENCE solutions with a non-zero tolerance level which was also specified according to the rule ((ref)). For $ n=842$, this amounts to setting the MIO tolerance level to be about 4.4%. To further reduce the computational cost, we adopted the warm-start strategy in the resulting MIO formulations. We set $\tau =1.5$ and constructed $( \widehat{P}_{i})_{i=1}^{n}$ using the fitted choice probabilities from the logit regression of $Y$ on all the covariates to derive the refined space $ \widehat{\Theta }\left( \tau \right) $ from the initial parameter space $ \Theta $.
Table 4 presents the refined parameter bounds derived from $\widehat{\Theta } \left( \tau \right) $ in this setup. The results of Table (ref) indicate that the size (measured in the volume of a ($k+p$) dimensional cube) of $\widehat{\Theta }\left( \tau \right) $ is only about 0.99% of that of $\Theta $; thus, there is a considerable reduction in the parameter search space by using the refined space $\widehat{\Theta }\left( \tau \right) $ in place of the original space $\Theta $.
We now present in Table (ref) the estimation results for the setup with the covariate specification ((ref)). From Table (ref), we can see that the two MIO formulations ((ref)) and ((ref)) yield the same set of selected variables across all the three cases of $q$. The parameter estimates computed from both formulations are also qualitatively similar in general. The variable $CARS$ remains to be selected in all these cases and its parameter estimate is also of the largest magnitude among all parameter estimates of the quadratic expansion variables. Moreover, there is very little loss in the goodness of fit from adopting only $CARS$ as the auxiliary covariate.
We now remark on the computational performance. For the cases of $q\in \{1,5\}$, formulation ((ref)) clearly outperformed formulation ((ref)) in both the CPU time and the number of branch-and-bounds used in the computation; however, both MIO formulations performed quite well for these two variable selection cases. By contrast, for the case of $q=3$, both approaches incurred more computational cost. In particular, it took around 18.3 hours and noticeably more branch-and-bound nodes to solve the formulation ((ref)) in the $q=3$ scenario. On the whole, these results suggest that the two MIO formulations ( (ref)) and ((ref)) are valuable complements for implementing the PRESCIENCE procedures.
By construction, the maximized average score of the PRESCIENCE approach increases with the specified value of $q$. In this empirical application, our results presented so far indicate that the in-sample predictive performance of the parsimonious predictive model using an intercept term and only the two variables $\left( DCOST,CARS\right) $ seems comparable to that of a more complex model using a richer set of covariates. We now investigate this issue further via the method of cross validation (CV).
We conducted the 5-fold CV analysis to assess the out-of-sample predictive performance of the PRESCIENCE methods for the selection of the variables specified by ((ref)). As in Table (ref), we considered the cases of $q\in \{1,3,5\}$ for which we implemented the corresponding PRESCIENCE procedures with a non-zero tolerance level. Because the training sample in each CV fold contains around 80% of the original observations, we set the MIO tolerance level for early termination to be about 4.9% by the rule ((ref)).
Table (ref) summarizes the 5-fold CV results, which are based on the averages over the performance results computed in each CV fold. From Table (ref), we can see that, for both MIO formulations, the in-sample maximized objective values were already very similar across the three cases of $q$ though they did strictly increase with $q$. Moreover, irrespective of the MIO formulations, the parsimonious case of $q=1$ had the best out-of-sample performance.
We now inspect the variables selected in each of the 5 CV folds. Figure (ref) summarizes the parameter estimates computed in each CV fold for the 9 auxiliary covariates specified by ((ref) ). From this figure, we also note that $CARS$ was selected across all $q$ cases in all CV folds. Its parameter estimate was also of a relatively large magnitude when compared to those of other selected variables. These cross-validation results further strengthen the finding that $CARS$ may be the most important predictive variable for the work-trip mode choice.
In this paper, we consider the variable selection problem for predicting binary outcomes. We study the best subset selection procedure by which the covariates are chosen by maximizing the maximum score objective function subject to a constraint on the maximal number of selected covariates. We establish non-asymptotic upper and lower risk bounds for the resulting best subset maximum score binary prediction rule when the dimension of potential covariates is possibly much larger than the sample size. The derived upper and lower risk bounds are minimax rate-optimal when the maximal number of selected variables is fixed and does not increase with the sample size. For implementation, we show that the variable selection problem of this paper can be equivalently formulated as a mixed integer optimization problem, which enables computation of the exact or an approximate solution with a definite approximation error bound.
The present paper takes the maximum score approach for the binary prediction problem. There is a large body of the literature that studies maximum score estimation in various other aspects since the seminal work by Manski:1975, Manski:1985. In the context of semiparametric binary response models, advances of the maximum score approach have been made in terms of point identification Manski:1988, partial identification Manski2002, komarova2013, blevins2015, chen2015, asymptotic distribution kim1990,Seo2016, panel data Manski87, CMvS:95, abrevaya2000, time series Moon2004, Guerre2006, deJong2011, dynamic network formation Graham:2016, nonparametrically generated regressors Chen-et-al:14, and so on. The numerical approach employed in this paper can be adapted to these contexts.
Efficient computation is particularly demanding when it is necessary to obtain maximum score estimates repeatedly many times. This difficulty naturally arises, for example, in the context of resampling delgado2001, abrevaya2005, Lee:Pun:06, Patra2015 and change-point problems Lee2008, Lee-et-al:11. It would be an interesting future research topic to investigate numerical performance of our method (e.g. the warm-start procedure in Section (ref)) in these computation-intensive problems.
The maximum score approach has produced many offsprings: smoothed maximum score estimation Horowitz:92, horowitz2002, multinomial choice estimation Matzkin:93, Fox2007, integrated maximum score estimation Chen2010, Bayesian method benoit2012, alternative estimation based on local nonlinear least squares blevins-khan-2013-dfbr, blevins-khan-2013, Khan2013, estimation using local polynomial smoothing Chen:Zhang:2015, and non-Bayesian Laplace type estimation JPW15, JPW14 among many others. Some of these alternative estimation methods are equipped with algorithms that are easier to compute than the maximum score estimation. It is an interesting open question how to accommodate a variable selection problem in these alternative methods. One starting point can be the work by Jiang:10 who considered empirical risk minimization under the $\ell _{0}$ constraint with Gibbs posterior as an alternative to maximum score prediction. It might be also interesting to consider penalized estimation with the $\ell _{0}$, $\ell _{1}$ and/or $\ell _{2}$ penalty. The maximum score approach is also closely related to the maximum utility estimation framework LieliNieto-Barthaburu, LieliWhite, Elliott2013, LieliSpringborn for binary decision under model uncertainty. It would be interesting to generalize our results to this framework. These are also topics for future research.