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.
44,056 characters · 6 sections · 34 citation commands
Exact Computation of Maximum Rank Correlation Estimator
\onehalfspacing
In this paper we provide a computation algorithm to get a global solution for the maximum rank correlation (MRC) estimator using the mixed integer programming (MIP) approach. The new algorithm returns a global solution. The MRC estimator was first proposed by han1987non to estimate the generalized regression model:
where $D:\mathbb{R}\mapsto\mathbb{R}$ is non-degenerate monotonic and $F:\mathbb{R}^2\mapsto\mathbb{R}$ is strictly monotonic in each arguments. The object of interest is the linear index parameter $\beta$. The model is general enough to include a binary choice model, a censored regression model, and a proportional hazards model as its example. han1987non proposed to estimate $\beta$ by maximizing Kendall's rank correlation coefficient:
where $1\left\{ \cdot \right\}$ is an indicator function. He showed the consistency of the MRC estimator and sherman1993limiting proved the $\sqrt{n}$-consistency and the asymptotic normality later. The flexible model structure leads to various extensions of the MRC estimator: for example, a quantile index model (khan2001two), a generalized panel model (abrevaya2000rank), a rank estimation of a nonparametric function (chen2002rank), a functional coefficient model (shin2010local), a random censoring model (khan2007partial), and a partial linear model (abrevaya2011rank).
There exist various semiparametric estimators in the class of single-index models (see, for example, the recent work by ahn2018simple and the references therein). Compared to them, the MRC estimator has the following advantages. First, it does not require any bandwidth selection since it does not involve any nonparametric estimation components. Second, it can be applied to various models without much modification (see, e.g., the references above and \citet*{khan2019inference} for the multinomial models). Finally, it is point robust in the sense that it provides a nontrivial identified set that includes the true parameter value when sufficient conditions for point identification are not satisfied (see the discussion in khan2018discussion for details).
Implementing the MRC estimator casts some computational challenges in practice, where the grid search method is not feasible. First, the objective function in (ref) is not differentiable in $\beta$ and we cannot apply a gradient-based optimization algorithm. Second, the objective function is not concave. Therefore, any solution found by a numerical algorithm could not be a global solution but a local one. This difficulty is well described in chay1998estimation, where they apply Powell's conjugate directions method, the simplex method with multiple starting points, and the piece-wise grid search method repeatedly to achieve a better solution in the empirical application. Even after these repeated searches, we are not sure whether the current solution is the global optimum. Finally, the objective function is the second order U-process and requires $O(n^2)$ computations for a single evaluation. abrevaya1999computation shows that the computation order can be reduced to $O(n\log n)$ by adopting the binary search tree structure. However, the fundamental local solution issue still remains.
The contribution of this paper is twofold. First, we propose a new computation algorithm that assures the global solution of the MRC estimator. We achieve this goal by transforming all indicator functions into binary parameters to be estimated along with additional constraints. We show that the proposed mixed integer programming (MIP) problem is equivalent to the original optimization problem. Although MIP is still an NP(non-deterministic polynomial-time)-hard problem (see, e.g. wolsey1998integer and johnson1978densest for details), we use a modern mixed integer programming (MIP) solver and confirm that it is feasible to get the solution within a reasonable time budget. The additional advantage of the MIP approach is that it provides us with the gap between the objective function value at the current best solution and the bound of the possible global maximum at any time point of the computation procedure. By this MIP gap, we can measure the quality of the interim solution when the time limit prevents us from waiting for the convergence of the procedure. Second, we consider an application of the best subset rank prediction and analyze the prediction performance. Building on chen2018best, we derive a non-asymptotic bound of the tail probability of the predictive performance measure. Since the objective function is defined as a second-order U-process, we develop a new technique to derive the finite sample tail probability bound for higher order U-processes.
We review some related literature. The MIP procedure is recently adopted in various applications in econometrics and statistics. florios2008exact show that the maximum score estimator of manski1975maximum can be reformulated as an MIP structure. \citet*{bertsimas2016best} consider the best subset selection problem and show that the MIP algorithm outperforms other penalty based methods in terms of achieving sparse solutions with good predictive power. chen2018best, chen2018exact investigate the binary prediction problem with variable selection and the instrumental variable quantile regression in the MIP formulation. kitagawa2018should apply the MIP procedure when they estimate the personalized optimal welfare policy. Finally, lee2018factor develop a MIP computation algorithm to estimate a two-regime regression model when the regime is determined by multi-dimensional factors. To the best of our knowledge, however, this is the first paper in the literature to apply the MIP approach when the objective function is defined as a higher order U-process.
The remainder of the paper is organized as follows. In section 2, we propose the MIP computation algorithm for the maximum rank correlation estimator. We show that the proposed algorithm is equivalent to the original optimization problem and illustrate how it can achieve a feasible solution. In Section 3, we consider the best subset rank prediction problem and derive the non-asymptotic tail probability bound of the performance measure. In section 4, we show the better performance of the proposed MIP algorithm by applying it to the female labor participation data of mroz1987sensitivity. Additional numerical evidence is provided through Monte Carlo simulation studies in section 5. We provide some concluding remarks in section 6.
In this section we describe the computational challenges of the maximum rank correlation (MRC) estimator and propose a new algorithm to compute a global solution of it. We illustrate the advantage of the new algorithm by investigating a simple numerical example.
We first discuss the computational difficulties of the MRC estimator. Recall that MRC is defined as follows:
where $\mathcal{B}$ is the parameter space of $\beta$ and $1\left\{ \cdot \right\}$ is an indicator function. Note that the objective function is neither differentiable nor concave in $\beta$. Furthermore, it is defined as a second-order U-process, which requires $O(n^2)$ order of computations for each evaluation of a candidate parameter value.\footnote{abrevaya1999computation proposes a nice algorithm that reduces the computation order to $O(n\log n)$ by using the binary search tree. However, it still does not guarantee the global solution.} As a result, we cannot apply any gradient-based optimization algorithm. Researchers usually adopt a simplex-based algorithm such as the Nelder-Meade method in MRC applications. However, it is difficult to get the global solution even with multiple starting points since the objective function is not globally concave. A grid search algorithm would give more robust solutions but the curse of dimensionality makes it infeasible in most cases when the dimension of $x$ is larger than 2.
In this paper we propose an alternative computational algorithm that is based on the mixed integer programming (MIP) procedure. Let $x_{ij} := x_i - x_j$ be a pairwise difference of $x_i$ and $x_j$. Let $\ensuremath{\varepsilon}$ be a small positive number, e.g.\ $\ensuremath{\varepsilon}=10^{-6}$, to denote an effective zero. Consider the following mixed integer programming problem: for $i,j =1,\ldots, n$ and $i \neq j$
where $ M_{ij} = \max_{\beta \in \mathcal{B}} \left\vert x_{ij}' \beta \right\vert + \ensuremath{\varepsilon} $. \footnote{Note that $M_{ij}$ is not a user-chosen turning parameter as it is the empirical bound of $\vert x_{ij}'\beta\vert$ determined by the parameter space and the data.} Since the objective function in (ref) is the linear function of the binary variables $\{d_{ij}\}$, the formulation becomes a mixed integer linear programming problem. We check the equivalence between the original problem in (ref) and the MIP problem in (ref)--(ref). Consider that the MIP problem chooses $\ensuremath{\widehat}{d}_{ij}=1$ for some $i,j$. Then, the constraint (ref) implies that the estimate for $\ensuremath{\widehat}{\beta}$ should satisfy $0<x_{ij}'\ensuremath{\widehat}{\beta} \le M_{ij}$, which is equivalent to $x_i\ensuremath{\widehat}{\beta} > x_j\ensuremath{\widehat}{\beta}$ for a large enough $M_{ij}$. Similarly, $\ensuremath{\widehat}{d}_{ij}=0$ is equivalent to $x_i\ensuremath{\widehat}{\beta} \le x_j\ensuremath{\widehat}{\beta}$. In sum, the constraint forces ${d}_{ij}=1\{x_{ij}'\beta >0\} = 1\{x_i'\beta > x_j'\beta\}$ given any $\beta \in \mathcal{B}$. Therefore, we can compute the global solution $\ensuremath{\widehat}{\beta}$ for (ref) by solving the equivalent MIP problem in (ref)--(ref). These two optimization problems give us the same numerical results but the MIP procedure has a clear computational advantage over the original problem, which is illustrated below.
Modern numerical solvers such as CPLEX and Gurobi make it possible to solve a large scale MIP problem by adopting branch-and-bound type approaches. We provide a heuristic explanation of how a vanilla branch-and-bound algorithm reduces the computational burden followed by a numerical example. Consider a binary tree representation for all possible values of $\{d_{ij}\}$ (for example, see Figure (ref)). A bottom node of the tree represents a different possible solution for $\{d_{ij}\}$, and $\beta$ can be easily solved by the linear programming procedure since $d_{ij}$ is fixed there. However, we have $2^{n(n-1)}$ bottom nodes in total and the brute force approach is still infeasible with a standard sample size. The branch-and-bound approach help eliminate a partition of the final nodes systematically. Suppose that we are in a node located in the middle of the tree, where only a part of $\{d_{ij}\}$ is fixed. Let $U^*$ be the current best objective function value.\footnote{An initial solution can be achieved from the linear programming problem at any bottom node of the tree.} Now we solve a subproblem after relaxing all $\{d_{ij}\}$ that are not fixed by parent nodes into continuous variables on the interval $[0,1]$. This relaxed subproblem can be solved easily by the linear programming procedure since it does not contain integer parameters anymore. There are two cases where we can reduce the computational burden. First, the objective function value of the relaxed problem, say $U^*_R$, is less than or equal to $U^*$. Since the objective function value of the original subproblem is always worse than that of the relaxed subproblem, we cannot achieve a better result than $U^*$ by solving any bottom nodes below the current node. Thus, we can drop all of them from the computation list. Second, $U^*_R > U^*$ and the solution of the relaxed problem satisfies the binary restriction for $\{d_{ij}\}$. This solution coincides with that of the original subproblem. Then, we update $U^* = U^*_R$ and can drop all bottom nodes below it from the computation list. While moving on to a child node, we solve a relaxed subproblem repeatedly and drop a partition of the bottom nodes from our computation list.
We provide a simple numerical example to illustrate how the branch-and-bound algorithm works.
Finally, we have some remarks on the implementation of the MIP procedure in (ref)--(ref). First, $M_{ij}$ can be computed by solving a separate linear programming problem and saving those values. Alternatively, we can also set a big number $M_{max}$ for all $i,j$, which is large enough to cover the absolute bound of the index $ x_{ij}'\beta$. The second order U-process property requires $O(n^2)$ computations for getting each $M_{ij}$ and it is usually faster to impose a constant number $M_{max}$ for all $M_{ij}$ than to solve the linear programming problem for each $i.j$ in our simulation studies. Second, it is well-known in the numerical optimization literature that any strict inequality should be switched into a weak inequality with some numerical precision bound. Thus, we change the second constraint in (ref) into
Third, when there are many tied observations in the dependent variable, we can reduce computation cost substantially by vectorizing paired observations and dropping the tied pairs as we have observed in Example (ref).
In this section we consider the application of a rank prediction problem. The goal is to find a linear index model that gives the best rank prediction of $y$ given $x$. The dimension of $x$ is potentially large and the model selection turns out to be the selection of the best predictors, i.e. the best subset. We propose an $\ell_0$-constraint maximum rank correlation estimation procedure and show that the MIP method in (ref)--(ref) can be immediately extended to this estimation problem. Building on chen2018best, we also provide the non-asymptotic bound of the rank prediction error. This bound implies that the dimension of $x$ can grow exponentially fast if the best subset size grows slowly enough, e.g. at a polynomial rate.
Suppose that we have a training set of $\{(y_i,x'_i): i=1,\ldots,n\}$, where $y$ can be either discrete or continuous. We want to learn the rank prediction rule for $y$ as well as to select the best $s$ predictors among $x$'s. Let $x=(x_1,x_{-1}')$ be $(p+1)$ covariates and we know that $x_1$ should be included in the predictor set. Let $\Vert \cdot \Vert_0$ be the $\ell_0$-norm, i.e.\ $\Vert \beta \Vert_0$ is the number of non-zero elements of the vector $\beta$. For any $k\neq l$, we propose the following prediction rule:
where $R_{\beta}(x_k,x_l)=1$ implies that $y_k$ is predicted to be larger than $y_l$. When we are given the whole prediction set $\{x_l: l=1\ldots, n_p\}$, the rank of $y_k$ is predicted by $\sum_{l=1}^{n_p} R_\beta(x_k,x_l)$. Let $F$ be the joint distribution of $(Y,X)$ and $Q:=P\times P$ be the product measure of $P$. Then, we choose the prediction rule as a sample analogue of
Recall that we also want to select the best $s$ predictors out of the total $p$ covariates of $x_{-1}$. Therefore, the prediction rule composed of the best $s$ predictors can be achieved by solving the the following $\ell_0$-constraint optimization problem:
where $\mathcal{B}_s : = \{\beta \in R^p: \Vert \beta \Vert_0 \le s \}$ and
We evaluate the performance of the predictor by the following measure:
where $S^*_{s}:=\sup_{\beta \in \mathcal{B}_s} S(\beta)$ and $\ensuremath{\widehat}{\beta}$ is the solution of the constraint maximization problem defined in (ref)--(ref) above. Note that $U_n \ge 0$ by the definition of $S_n^*$ and that a good prediction rule results in a small value of $S_n$ with a high probability. In the next theorem, we provide a non-asymptotic bound of $U_n$. Let $a \vee b := \max\{a,b\}$ and $r_n:= s\ln(p\vee n) \vee 1$.
Theorem (ref) shows that the tail probability of $U_n$ decreases exponentially in $r_n$. The probability bound in (ref) is non-asymptotic and holds for every $n$ if two inequality conditions (ref)--(ref) hold. Compared to the non-asymptotic bound of the best subset selection in chen2018best, Theorem (ref) requires an additional condition (ref) to bound the second order degenerate U-process. However, focusing on the leading terms, we confirm that both conditions hold if
Note that Theorem (ref) implies that $E(U_n) = O(n^{-1/2}\sqrt{s\ln(p\vee n)}) = o(1)$ if $s\ln(p\vee n)=o(n)$. Therefore, the best subset rank prediction performs well even when $p$ grows exponentially provided that $s$ increases slowly, e.g.\ at a polynomial rate.
We finish this section by formulating the $\ell_0$-constraint optimization problem as an MIP problem. Let $x_{-1,ij}:=x_{-1,i} - x_{-1,j}$ as before. For $i,j=1,\ldots,n$, $i\neq j$, $h=1,\ldots,p$, we consider the following constraint MIP problem:
where $\underline{\beta}_h$ and $\overline{\beta}_h$ are the lower bound and the upper bound of $\beta_h$, respectively. The constraint MIP problem in (ref)--(ref) is equivalent to the original constraint optimization problem. The objective function in (ref) is numerically same with $S_n$ since $d_{ij}$ is identical to $R_{\beta}$ for each $\beta$. Furthermore, the constraint (ref) makes $\beta_h=0$ whenever $e_h=0$. Thus, the $\ell_0$-norm constraint $\Vert \beta \Vert_0 \le s$ is achieved by the constraints (ref), (ref) and (ref). Note that the objective function can be also written in the familiar rank correlation form:
which is equivalent to (ref).
In this section we illustrate the advantage of the MIP procedure in an empirical application. We revisit the female labor force participation application in mroz1987sensitivity and estimate the binary choice model using the generalized regression model in (ref). Specifically, the composite functions are defined as $F(x'\beta,\ensuremath{\varepsilon}):=x'\beta+\ensuremath{\varepsilon}$ and $D(A):=1\{A>0\}$ so that it becomes a semiparametric binary choice model:
where the distribution of $\ensuremath{\varepsilon}_i$ is not specified. The parameter of interest is $\beta$ and we estimate it using the maximum rank correlation estimator defined in (ref). The outcome variable, $y_i$, is 1 if she participated in the labor force and 0, otherwise. We choose the following seven variables from the data for the covariate $x_i$: the number of kids less than 6-year-old ($kidslt6$), the number of kids aged between 6 and 18 ($kidge6$), years of education ($educ$), family income minus her income ($nwifeinc$) in \$1,000, years of experience ($exper$), experience squared ($expersq$), and age ($age$). We randomly draw 100 observations out of 753 for this computation exercise. Table (ref) reports summary statistics of both samples and we confirm that there is not much difference in terms of the mean and the standard deviation of each variable. We normalize the coefficient of $kidslt6$ to be -1. Note that the grid search method is infeasible given the sample size and the number of regressors in this application.
Table (ref) summarizes the estimation results. First, we estimate the model using the mixed integer programming procedure (MIP) with a time budget of 600 seconds. To compare its performance with the existing methods, we also estimate it using the following five methods: the Nelder-Mead simplex method with an initial value from OLS (Nelder-Mead 1), the Nelder-Mead method with multiple initial values until the time budget of 600 seconds is reached (Nelder-Mead 2), the iterative grid search method (Iter-Grid), the simulated annealing method (SANN), and the Markov Chain Monte Carlo (MCMC) method in chernozhukov2003mcmc. The parameter space was set to be $\mathcal{B}=[-10,10]^6$. The random starting points of Nelder-Mead 2 was generated from the uniform distribution on $\mathcal{B}$. We use the 2,001 equi-spaced grid points for each parameter for Iter-Grid. The Nelder-Mead method has been adopted in the applications of the MRC estimator, where the grid search is infeasible (for example, see cavanagh1998rank, abrevaya2003pairwise, khan2009inference). A more sophisticated version of the iterative grid search method is introduced by wang2007note and it is adopted in \citet*{fan2020rank} for their simulation studies with multi-dimensional regressors.
The estimation result in Table (ref) reveals several advantages of MIP over the existing alternative algorithms. First, MIP achieves the best objective function value among the candidate estimation methods within a reasonable time budget. Second, some estimates of $\hat{\beta}$ by alternative algorithms are qualitatively different from the solution of MIP. The coefficient estimate of $kidsge6$ by Nelder-Mead 1 shows the opposite direction. The estimates of $educ$ by alternative algorithms show much higher effects than MIP. MCMC shows the closest result although it is still suboptimal. Third, the Nelder-Mead algorithm with the multiple staring point for 600 seconds does not improve the result. In fact, the objective function value of Nelder-Mead 2 is lower than that of Nelder-Mead 1 which uses only one starting point of the OLS estimate. Finally, Figure (ref) shows how difficult the optimization problem is. We plot the empirical objective function values over the convex combinations of two $\beta$ estimates of MIP and Nelder-Mead 1. We can confirm that the objective function is not concave and that there exist many local maxima even between these two estimates.
In sum, inference based on inferior local solutions could lead researchers to imprecise or incorrect conclusions in practice although the theoretical properties of the MRC estimator are robust.
In this section we investigate the performance of the proposed MIP algorithm for the MRC estimator via Monte Carlo simulation studies. We focus on the achieved objective function value and the computation time in this section. All simulations are carried out on a computer equipped with AMD Ryzen Threadripper 1950X 16-Core processor and 64 Gigabytes of RAM.
We consider the following two regression models: for $i=1,\ldots, n$,
where $x_i$ is a $k$-dimensional vector generated from $N(0,I_k)$, $\ensuremath{\varepsilon}_i$ is an error term generated from $N(0,0.25^2)$, and $\beta$ is a parameter of interest. The true parameter value is set to be $\beta_0=(1,\ldots,1)$. Recall that we do not know the true transformation function (binary or censored) of the data generating process when we estimate $\beta$ by the MRC estimator. For identification, we normalize the first coefficient of $\beta$ to be 1. We compare the performance of the MIP algorithm with that of the Nelder-Mead algorithm (NM), the iterative grid search (Iter-Grid), the simulated annealing (SANN), and the Markov Chain Monte Carlo method (MCMC). For all methods, the parameter space is set to be $\mathcal{B}=[-10,10]^{k-1}$. The time budget is set to be 600 seconds. The Nelder-Mead algorithm with repeated random starting points (Nelder-Mead 2 in the previous section) does not perform better (especially for large $k$) than Nelder-Mead 1 and is dropped in these simulation studies.
We first consider small-scale designs and check if the MIP algorithm achieves the global objective function. We set the sample size and the dimension of regressors to be $n=(50,100)$ and $k=(2,3)$, respectively for these small-scale designs. We next extend them into $n=(200,400)$ and $k=(10,20)$ and check the performance of the MIP algorithm in the limited time budget. Therefore, we consider 8 different designs in total in each regression model (binary/censored). We conduct 10 replications of each simulation design.
Figures (ref)--(ref) report the Loss/Tie/Win ratios of each alternative algorithm against MIP. In the graph, `Loss' means the objective function value of the algorithm is lower than that of MIP. `Tie' and `Win' are defined similarly. Overall, MIP outperforms the altenative methods. In the case of Binary Regression in Figure (ref), MIP always achieves an equal or better objective function value than the alternative methods in all designs except one draw in $n=400, k=10$. In the small-scale design ($n=50, 100$ and $k=2,3$), SANN performs similar to MIP but the performance of MIP dominates in the large-scale design ($n=200,400$, $k=10,20$). It is interesting that MIP performs better as the dimension of $k$ increases when $n=400$. As we confirm in Table (ref) below, MIP finds a more precise solution (lower MIP gap) in a substantially shorter time when $n=400$ and $k=20$ than when $n=400$ and $k=10$.
We observe similar patterns in Censored Regression in Figure (ref). The outperfomance of MIP is clearer in the large-scale designs ($n=200,400$ and $k=10,20$). The overall performance of MIP in Censored Regression is better than that in Binary Regression when $n=200$. However, when $n=400$, MIP finds a worse solution than its competitors about 10-20%. This is because the implied parameter space of $d_{ij}$ has a much bigger dimension in Censored Regression than Binary Regression as it has less tied pairs of $(y_i,y_j)$. Recall that $d_{ij}$ is multiplied by 0 in the objective function when $y_i$ and $y_j$ are tied and we do not need to estimate such a $d_{ij}$.
Tables (ref)--(ref) provide some summary statistics of the computation time and the MIP gap. We first discuss the result of Binary Regression in Table (ref). In small-scale designs, MIP requires about the same computation time as the alternative algorithms and it finds the global solution in less than a second except $n=100$ and $k=3$. In large-scale designs MIP is still able to find the global solution within the allocated time budget of 600 seconds, except when $n=400$ and $k=10$. In that design MIP hits the time limit of 600 seconds more often and the MIP gap does not achieve 0%, i.e.\ we are not sure whether the solution is global or not. However, the gap size is quite small and less than 1%. It is noteworthy that MIP performs much better in terms of the MIP gap when $k$ is bigger in large-scale designs. The computation time is even dramatically reduced when $n=400$.
We turn our attention to the result of Censored Regression in Table (ref). As we discussed above, Censored Regression requires more computation time than Binary Regression and it mostly reaches the time limit of 600 seconds when $n=100$ and $k=3$. In large-scale designs, we observe the MIP gaps larger than 1%, which could be the reason that the solutions of MIP are sometimes worse than those of the alternative algorithms. Other patters are quite similar to those in Binary Regression including that the performance of MIP becomes better when $k$ is higher for large-scale designs.
In sum, the performance of the proposed MIP algorithm for the MRC estimator is satisfactory. It always finds the global solution in small-scale designs where the existing methods fail to do quite often. Furthermore, it performs better than the alternative algorithms even in large-scale designs by spending a feasible amount of computation time. The MIP gap also provides useful guidance for the quality of a solution in hand when a researcher should stop searching for the global solution because of the time limit.
In this paper we propose a feasible computation method to get a global solution for the maximum rank correlation estimator of han1987non using the mixed integer programming (MIP). We show that the proposed MIP method outperforms the alternative methods in the empirical example of female labor-force participation. One advantage of the proposed MIP method is that it can be easily extended to many constraint rank maximization problems as illustrated in the best subset rank prediction problem, where we also prove that the non-asymptotic bound of the tail probability decays exponentially. This result sheds light on the research of the high-dimensional rank estimation models, which we leave for future research.