Smooth Sigmoid Surrogate (SSS): An Alternative to Greedy Search in Decision Trees

Abstract

Greedy search (GS) or exhaustive search plays a crucial role in decision trees and their various extensions. We introduce an alternative splitting method called smooth sigmoid surrogate (SSS) in which the indicator threshold function used in GS is approximated by a smooth sigmoid function. This approach allows for parametric smoothing or regularization of the erratic and discrete GS process, making it more effective in identifying the true cutoff point, particularly in the presence of weak signals, as well as less prone to the inherent end-cut preference problem. Additionally, SSS provides a convenient means of evaluating the best split by referencing a parametric nonlinear model. Moreover, in many variants of recursive partitioning, SSS can be reformulated as a one-dimensional smooth optimization problem, rendering it computationally more efficient than GS. Extensive simulation studies and real data examples are provided to evaluate and demonstrate its effectiveness.

1. Introduction

Recursive partitioning or tree-structured methods [1,2], together with their various extensions such as multivariate adaptive regression splines MARS; [3], bagging [4], boosting [5], and random forests [6], have formed an important class of learning and modeling tools for pattern recognition, data mining, and statistical analysis. Greedy search (GS, also called exhaustive search or brute-force search) is the dominating method for identifying the best split or knot at each step of recursive partitioning. In this discrete process, all permissible splits are compared via some criterion and the best one is selected.

Consider data that consist of i.i.d. observations $\{(y_i,{\mathbf{x}}_i):i=1,\dots ,n\},$ where $y_i$ is the i-th observed value on the response or target variable Y and ${\mathbf{x}}_i\in {\mathbb{R}}^p$ is the associated input vector on predictor variables $(X_1,\dots ,X_p).$ The response Y can be either continuous, as in regression trees, or categorical, as in classification trees, while the predictor variables can be of mixed types. A split of data is induced by a binary question such as “is $X_j$ greater than c?” for a continuous $X_j$ or “does $X_j$ belong to $C_0$?” for a categorical $X_j,$ where c is a cutoff point (or cutpoint for short) and $C_0$ is a subset of categories or levels that $X_j$ can assume. A greedy search for the best split can be viewed as a two-step process. In the first step, the best way to split each predictor, i.e., the best cutoff point or the best subset, is sought; then, in the second step, a comparison is made across all p predictors.

Consider the first step, namely, how to determine the best cutoff point for a given predictor $X_j.$ There are two general approaches of splitting data in the literature: either by minimizing within-node impurity, or by maximizing between-node difference. The conventional approach, as in CART ([1], pp. 230–232), splits data by minimizing the within-node impurity. Let $𝚤\left(t\right)$ denote the impurity measure for node t and let $\Delta 𝚤=𝚤\left(t\right)-\left\{p\left(t_L\right)𝚤\left(t_L\right)+p\left(t_R\right)𝚤\left(t_R\right)\right\}$ denote the reduction in impurity due to a split of the data, where $t_L$ and $t_R$ denote the left and right child nodes, while $p\left(t_L\right)$ is the proportion of observations falling into $t_L$ and analogously for $p\left(t_R\right).$ For splits based on variable X, the best cutoff point $c_j^{★}$ achieves the greatest reduction in the sense that $c_j^{★}=arg{max}_c\Delta 𝚤.$ The other splitting approach is to maximize the between-node heterogeneity, as examined by LeBlanc and Crowley [7]. This approach is somehow more appealing from the statistical perspective and allows for the extension of tree methods to other data types, as any appropriately defined two-sample statistic may be used to bisect the data. In many scenarios, these two approaches lead to equivalent splitting criteria.

Although GS is the standard method for splitting data, it has a few deficiencies. First of all, the splitting statistic in GS is a random process that shows erratic patterns with substantial variations, as mentioned by [8] (Figure 1 on p. 8) in a recent rpart documentation, for example. As a result, it often occurs that local spikes overpower the overall peak associated with the true cutpoint, especially when the signal strength is relatively weak. Second, some commonly used splitting statistics in GS tend to favor unbalanced cuts that leave a child node with very small sample sizes, a problem referred to as end-of-cut preference (see Section 11.8 of CART [1]). This would eventually lead to tree structures that are hard to interpret. Figure 1 provides an illustration of the problems with simulated data. The goodness-of-split measures shows an overall pattern climaxing at the true cutpoint of 0.5, yet with high variation and many high local spikes; in particular, one prevailing local peak at the end is erroneously selected by GS. Third, variable selection bias has been considered inherent with GS. This refers to the phenomenon of a predictor with more values or levels being more likely to be selected as the splitting variable than a predictor with fewer values or levels. Loh and colleagues have pioneered excellent efforts in addressing this issue with a series of papers, as exemplified in Loh [9]. They also view GS as a two-step process, yet in a different way. Their main idea is to first select the most important variable and then evaluate splits based on the selected variable only. Another suggestion regarding this issue is to derive the distribution of maximally selected test statistics. This approach is exemplified by Miller and Siegmund [10], Shih and Tsai [11], and Hothorn and Zeileis [12]. Finally, searching over all permissible cutoff points for each predictor can be time-consuming. One common way for speeding up GS is to apply an updating formula to compute the splitting statistic for two consecutive cutoff points on a predictor. Examples can be seen in (Section 2.2, [13]), LeBlanc and Crowley [7], and Su, Tsai, and Wang [14] (Appendix A). Other strategies include splitting on percentiles, sub-sampling, and online incremental learning (see, e.g., [15]), which nevertheless reduces the search space and consequently provides no guarantee of success in finding the optimal cutoff point. Despite all of these efforts, greedy search can be time-consuming, especially when dealing with continuous variables or categorical variables that have many levels in large datasets.

In this article, we propose an alternative splitting method to greedy search, termed smooth sigmoid surrogate (SSS). The main idea is to replace the step functions involved in splitting with smooth sigmoid functions. The best cutoff point for each predictor can be estimated as a parameter in a nonlinear model. Very often, the estimation can be appropriately reformulated into a one-dimensional optimization problem that can be quickly solved. We demonstrate that this new method provides superior performance to GS along with correction or amelioration of its common deficiencies. First of all, SSS facilitates a parametric smoothing or regularization of the erratic splitting statistics in GS, yielding improved stability in the objective function to be optimized. As a result, SSS is more capable of identifying weak signals. As shown in Figure 1, smoothing preserves the overall parabola pattern while flattening the local spikes and reducing variations. A cutpoint close to the true one of 0.5 is successfully recovered for a range of smoothing parameter values. The smoothing effect also leads to substantial amelioration of the end-cut preference problem in practice. Furthermore, because the search for the best cutoff is cast in a nonlinear regression framework, conventional statistical inference can be exploited to facilitate convenient comparisons across different predictors for finding the best split of data. We show that this alternative splitting method helps to alleviate the variable selection bias problem. In addition, SSS potentially reduces computational cost without sacrificing performance. Lastly, SSS can be flexibly extended to many other recursive partitioning methods designated for different analytic purposes. In most scenarios, the optimization problem remains one-dimensional with appropriate reformulation.

We note that replacing step functions with smooth sigmoid functions is not an entirely new idea. Notably, this tactic has impacted the development of artificial neural networks (ANN). One well-known initial ANN proposal by Rosenblatt [16] consisted of single-layer networks called perceptrons, which contain threshold activation functions. In later developments such as multilayer perceptrons (MLPs) the step function is replaced by the differentiable logistic function, which has advanced ANNs tremendously. In the recursive partitioning context, the hierarchical mixture of experts (HME) architecture studied by Jordan and Jacobs [17] is essentially a mixture model where mixing probabilities are hierarchically modeled via multivariate sigmoid or softmax functions. Along similar lines, Ciampi, Couturier, and Li [18] proposed a ‘soft tree’ procedure at each node in which an observation goes to the left or right child node with a certain probability that is estimated via soft threshold functions or logistic models. To the best of our knowledge, however, our current proposal is the first attempt to use the sigmoid function as an alternative to GS for partitioning data. Therneau [8] also noted the potential of nonparametric smoothing in stabilizing the erratic pattern shown on the splitting statistic in GS; however, he applied nonparametric smoothing splines directly to the computed goodness-of-split measures, which could only add to the computational cost of GS. Comparatively, SSS facilitates a parametric way of smoothing that is embedded in the process that generates the splitting statistic.

The remainder of this paper is organized in the following manner: Section 2 introduces the SSS method for regression trees, where the response variable is continuous, and compares it to GS using various numerical studies in many aspects; classification trees with binary outcomes are treated in Section 3; Section 4 outlines the usage of SSS in several other recursive partitioning-based methods; finally, Section 4 ends the article with a brief discussion.

2. SSS for Regression Trees

In regression trees, the node impurity $𝚤\left(t\right)$ is commonly measured by the mean squared error $𝚤\left(t\right)={\sum }_{i\in t}{\left(y_i-{\overline{y}}_t\right)}^2/n_t.$ In this case, minimizing the within-node impurity is equivalent to maximizing a two-sample F test statistic that compares the two child nodes. For the time being, we consider continuous predictors only. The scenario for categorical predictors as well as binary variables will be discussed later. For simplicity, we denote $X_j$ as $X.$ The splitting statistic can be cast into the following linear model:

$$y_i={\beta }_0+{\beta }_1\delta (x_i;c)+{\epsilon }_i$$

(1)

where $\delta (x_i;c)=I\{x_i\ge c\}\triangleq {\delta }_i$ is the indicator function corresponding to the split and ${\epsilon }_i\stackrel{\mathrm{iid}}{\sim }N(0,{\sigma }^2).$ The two-sample F test is the same as the squared t test for testing $H_0:{\beta }_1=0$ vs. $H_a:{\beta }_1\ne 0.$ The indicator function in (1) can be written as $\delta (x;c)=I\{x-c\ge 0\}=H(x-c)$, where $H\left(x\right)=I\{x\ge 0\}$ denotes the Heaviside (unit) step function. Often used in integration, the value of $H\left(x\right)$ at $x=0$ does not really matter, whereas one natural choice is $H\left(0\right)=1/2$ for the sake of rotational symmetry.

Given that this search must be carried out over all variables, this can lead to a significant computational strain in a big data context whenever many features are present, with each containing several potential splits; it also suffers from other problems, as discussed earlier. This motivates us to consider a new way of splitting data. Our idea is to approximate the threshold indicator function with a smooth sigmoid surrogate (SSS) function.

2.1. Sigmoid Functions

A sigmoid function is a mathematical function with an ‘S’ shape. Very often, the sigmoid function refers to the special case of the logistic function (also called the expit function):

$$\pi \left(x\right)=\pi (x;a)={\{1+exp(-ax)\}}^{-1}={\displaystyle \frac{exp\left(ax\right)}{1+exp\left(ax\right)}}.$$

(2)

While there are various types of sigmoid functions, we primarily choose the logistic function due to its widespread use in statistics, such as in logistic regression, and the extensive research on its properties, e.g., $d\pi (x;a)/dx=a\pi (x;a)\{1-\pi (x;a\left)\right\}.$ Readers are referred to Balakrishnan [19], Johnson, Kotz, and Balakrishnan [20] and Rosenblatt [16] (Section 23) for a book-length account of the logistic function and the logistic distribution as well as their properties and applications.

Figure 2 plots the logistic function with different values of the shape or scale parameter $a=1,2,\dots ,100$. Note that $\pi (x-c;-a)=1-\pi (x-c;a)$, hence; $\pi (x-c;a)$ with $a<0$ approximates $I\{x<c\}.$ Either $\pi (x-c;a)$ or $\pi (x-c;-a)$ can be used to represent a split in the recursive partitioning setting. Therefore, without loss of generality (WLOG), our consideration is restricted to $a>0$. As a increases, $\pi \left(x\right)$ provides a better approximation of the indicator function $I\{x\ge 0\}$. Similarly, $\pi (x-c;a)$ approximates $I\{x\ge c\}.$ On the other hand, $\pi \left(x\right)$ becomes linear for small a. Because tree methods model data with threshold effects, we are primarily interested in seeking an approximation to the indicator function, where a relatively large a is desirable; hence, leaving a as a free parameter for estimation is not advisable. To fix a, we first note that the spread or scale of the predictor X on which the best cutoff is sought would affect the choice of a values. This concern leads to the specification of a as $a:=a/s_x$, where $s_x$ denotes the sample standard deviation of the predictor. This is equivalent to working with the standardized predictor. To make this more clear, standardization is a necessary first step in SSS splitting. In terms of approximation, note that ${lim}_{a\to \infty }\pi (x;a)=H\left(x\right)$ for any $x\ne 0$ with exponential rate, i.e., $\left|\pi (x;a)-H\left(x\right)\right|=O\{exp(-a)\}.$ Thus, a reasonably large constant a would provide a good approximation. Another possibility is to have a adaptive of the sample size n such that $a\equiv a_n$ satisfies ${lim}_{n\to \infty }{a}_n=\infty .$ Many rates can be used for $a_n$ because of the exponential convergence rate. The choice of a will be further explored and discussed in the subsequent sections.

2.2. Estimating the Best Cutoff Point

To continue, a natural approach is to replace $\delta (x_i;c)$ in (1) with a sigmoid function. The model then becomes

$$y_i={\beta }_0+{\beta }_{1}s(c;x_i)+{\epsilon }_i,$$

(3)

where $s(c;x_i)=\pi (x_i-c;a).$ Note that we have suppressed the parameter a in the notation, as it is fixed a priori. By choosing a large a that corresponds to a sharper ‘S’ shape and estimating the c corresponding to the best fit, we are essentially finding a sharp step function-type change in the data. Model (3) is a nonlinear parametric model involving four parameters $\{{\beta }_0,{\beta }_1,c,{\sigma }^2\}.$ Its estimation and related inference can be accomplished via nonlinear least squares [21]. In fact, the related optimization is a separable least squares problem (Section 9.4.1, [22]).

Estimating Model (3) involves multivariate optimization. However, we only need to estimate c for splitting purposes. This motivates us to approximate the splitting statistic directly and reduce the problem to a one-dimensional optimization with decision variable c only. We first note that the splitting statistic $\Delta 𝚤$ used in CART can be treated as an objective function for c and rewritten as follows:

$$\begin{array}{ccc}\hfill \Delta 𝚤\left(c\right)& =& \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2-\left\{\sum _{i\in t_L}{\left(y_i-{\overline{y}}_L\right)}^2+\sum _{i\in t_R}{\left(y_i-{\overline{y}}_R\right)}^2\right\}\hfill \\ & \widehat{=}& {\displaystyle \frac{1}{n_L}}·{\left(\sum _{i\in t_L}{y}_i\right)}^2+{\displaystyle \frac{1}{n_R}}·{\left(\sum _{i\in t_R}{y}_i\right)}^2\hfill \\ & =& {\displaystyle \frac{1}{n_L}}·{\left(\sum _{i\in t_L}{y}_i\right)}^2+{\displaystyle \frac{1}{n-n_L}}·{\left(\sum _{i=1}^n{y}_i-\sum _{i\in t_L}{y}_i\right)}^2,\hfill \end{array}$$

(4)

where $\{n_L,{\overline{y}}_L\}$ denote the sample size and the average response in the left child node, respectively; similarly, $\{n_R,{\overline{y}}_R\}$ for the right child node. Throughout this article, we use the notation $\widehat{=}$ to indicate ‘equivalence or correspondence up to some irrelevant constant’. In particular, we have ignored irrelevant terms that do not involve c in rewriting $\Delta 𝚤\left(c\right)$ in (4).

To approximate $\Delta 𝚤\left(c\right)$ in (4), it is only necessary to approximate $n_L$ and ${\sum }_{i\in t_L}{y}_i$ by replacing $\delta (x_i;c)$ with $s(c;x_i)$ such that

$$n_L=\sum _{i=1}^n\delta (x_i;c)\approx \sum _{i=1}^{n}s(c;x_i)$$

and

$$\sum _{i\in t_L}{y}_i=\sum _{i=1}^n{y}_i\delta (x_i;c)\approx \sum _{i=1}^n{y}_{i}s(c;x_i).$$

The approximated objective function, denoted as $\tilde{\Delta }𝚤\left(c\right)$, becomes

$$\tilde{\Delta }𝚤\left(c\right)={\displaystyle \frac{{\mathbf{s}}^T\left(\mathbf{y}{\mathbf{y}}^T\right)\mathbf{s}}{{\mathbf{j}}^T\mathbf{s}}}+{\displaystyle \frac{{(\mathbf{j}-\mathbf{s})}^T\left(\mathbf{y}{\mathbf{y}}^T\right)(\mathbf{j}-\mathbf{s})}{{\mathbf{j}}^T(\mathbf{j}-\mathbf{s})}},$$

(5)

where $\mathbf{y}=\left(y_i\right)$, $\mathbf{s}=\left\{s_i\right\}$ with $s_i=\left(s(c;x_i)\right),$ and $\mathbf{j}=\left(1\right)$ are n-dimensional vectors.

Suppose further that, WLOG, the response has been centered $\mathbf{y}:=({\mathbf{I}}_n-\mathbf{j}{\mathbf{j}}^T/n)\mathbf{y}$ such that ${\sum }_{i=1}^n{y}_i=0.$ It follows that ${\sum }_{i\in t_L}{y}_i=-{\sum }_{i\in t_R}{y}_i$; hence, $\Delta 𝚤\left(c\right)$ can be further simplified as ${\left({\sum }_{i\in t_L}{y}_i\right)}^2/\left\{n_L(n-n_L)\right\}$ up to some irrelevant constant. Its approximation $\tilde{\Delta }𝚤\left(c\right)$ reduces to

$$\tilde{\Delta }𝚤\left(c\right)={\displaystyle \frac{{\mathbf{s}}^T\left(\mathbf{y}{\mathbf{y}}^T\right)\mathbf{s}}{{\mathbf{s}}^T\left(\mathbf{j}{\mathbf{j}}^T\right)(\mathbf{j}-\mathbf{s})}}.$$

(6)

The objective Function (5) or (6) involves a scale parameter a. Whereas leaving a for estimation is important in HME [17] and soft trees [18], it is not a good idea for our purpose. If a is treated as a free parameter, a small a estimate that corresponds to a nearly linear relationship may provide a better fit, but may not be well suited for identifying the best cutoff point in ordinary tree procedures. This is because tree models account for the relationships between response and predictors via hard threshold functions. Second, if the best cutoff point is estimated for each predictor with a different scale parameter a, comparison across predictors becomes groundless. Specific recommendations on the choice of a will be provided in the subsequent sections. As discussed earlier, in order to fix a it is crucial to standardize the predictor $x_{ij}:=(x_{ij}-{\overline{x}}_j)/{\widehat{\sigma }}_j$, where $({\overline{x}}_j,{\widehat{\sigma }}_j)$ denote the sample mean and standard deviation of variable $X_j$, respectively. Alternatively, the range of the data may be normalized into the unit interval $[0,1]$ by setting $x_{ij}:=(x_{ij}-x_{\left(1\right)j})/(x_{\left(n\right)j}-x_{\left(1\right)j}),$ where $x_{\left(1\right)j}$ and $x_{\left(n\right)j}$ are the minimum and maximum of $X_j$, respectively.

With fixed a, the best cutoff point $\widehat{c}$ can be obtained by maximizing $\tilde{\Delta }𝚤\left(c\right)$ in (5) or (6) with respect to c, then transformed it back to the original data scale for interpretability. This is a one-dimensional optimization problem. It can be conveniently solved by many standard algorithms. In this paper, we use the method in [23] as implemented in the R [24] function optimize. Brent’s derivative-free method is very popular because it combines the bisection method, the secant method, and inverse quadratic interpolation by taking advantage of their best aspects.

One issue concerns non-concavity. The objective function in (5) or (6), although smooth, is not concave, as demonstrated in Figure 1. As a result, the algorithm may become stuck at a local optimum. In this case, global optimization procedures (see, e.g.,  [25,26]) can be helpful. Again thanks to the one-dimensional nature of the problem, one simple solution is to divide the search range into several intervals, find the maximum within each, and then identify the highest maximum among them. We have included this step as an option in our implementation, which may occasionally improve performance yet slow down the splitting procedure. This option is not used in any of the results reported throughout this paper. As shown in our numerical studies later, Brent’s method seems quite practically effective in locating the true cutpoint.

An immediate advantage of SSS over GS is computational efficiency. The following proposition provides an asymptotic quantification of the computational complexity involved in both GS and SSS splitting.

Proposition 1. 

Consider a typical dataset of size n in the regression tree setting. Both GS and SSS are used to find the best cutoff point of a continuous predictor X which has $O\left(n\right)$ distinct values. In terms of computation complexity, GS is at best $O\{ln(n\left)n\right\}$ with the updating scheme and $O\left(n^2\right)$ without the updating scheme, while SSS is $O\left(mn\right)$, where m is the number of iterations in Brent’s method.

The proof of Proposition 1 is provided in Appendix B. It is a common wrong impression that LS splitting with updating is only $O\left(n\right).$ This is because updating entails sorting the response values according to the X values. It turns out that this step would dominate the algorithm asymptotically in terms of complexity. Comparatively, SSS depends on the number of iterations in Brent’s method, denoted as $m.$ Although the number of iterations is affected by the convergence criterion and the desired accuracy, m is generally small, as Brent’s method has guaranteed convergence at a superlinear rate. Based on our limited numerical experience, m rarely reaches over 12 even for large n. In other words, the $O\left(mn\right)$ rate for SSS essentially amounts to the linear rate $O\left(n\right).$ Therefore, SSS can be more computationally efficient in least squares splitting. Additionally, SSS can be even more advantageous in other recursive partitioning scenarios where an updating formula is not readily available, especially when evaluating each individual split point involves iterative model estimation procedures.

To speed up GS, subsampling or splitting on percentiles are two common strategies. Clearly, subsampling can be combined with SSS as well. Nevertheless, reduced samples could weaken performance and leave weaker signals undetectable. For percentile splitting, the number g of percentiles used is an issue. In a single tree analysis, the investigator is keenly interested in the ‘best’ cutoff point for the sake of interpretability and relevance to the scientific or business question at hand. Using percentiles reduces the fineness level to 1% at best in identification of the best cutpoint, which is undesirable, especially when dealing with big data. Moreover, computing all g percentiles of X is not computationally thrifty. In terms of complexity, this step is only $O\left(gn\right)$, and g could be much larger than $m,$ i.e., the number of iterative steps in SSS.

2.3. Simulation Studies on Single Splits

To evaluate SSS and investigate several relevant issues, we simulated data from the following model:

$$\mathrm{Model} \mathrm{A}:y=1+\beta ·I\{x\le c_0\}+\epsilon \mathrm{with}\epsilon \stackrel{IID}{\sim }\mathcal{N}(0,1),$$

(7)

where the predictor X is generated from a uniform [0, 1] distribution. Two sample sizes $n\in \{50,500\}$ were investigated, as recursive partitioning always involves both small sample and large sample problems. We considered two true cutoff points $c_0\in \{0.50,0.80\},$ corresponding to the balanced and unbalanced scenarios, and three values for $\beta \in \{0.0,0.2,1.0\}$, corresponding to null, weak, and strong signals, respectively. These considerations form ten combinations in total, as the choice of $c_0$ no longer matters in the null case of ${\beta }_1=0.$ For each model configuration, 200 simulation runs were used. For each simulated dataset, both GS and SSS were used to identify the best cutoff point $\widehat{c}.$ For SSS, a spread of choices were considered for $a\in \{1,2,\dots ,100\}.$ With the identified cutpoints from all simulation runs, we obtained the empirical density curves of $\widehat{c}$, as presented in Figure 3.

Figure 3a,b corresponds to the null model case. The results of GS show extraordinarily high peaks at both ends. This phenomenon can be explained by the ‘end-cut preference’ (ECP) problem (see Section 11.8 of CART; [1]), which is known for the fact that the least squares criterion tends to favor unbalanced splits, i.e., splits in which one child node has only a very small proportion of observations. Comparatively, the smoothing effect of SSS substantially ameliorates the ECP problem. Although SSS is intended to approximate GS, weighting introduced by the SSS smoothing process mitigates the ECP problem markedly by changing the summations, e.g., ${\sum }_{i\in t_L}{y}_i$, into weighted sums ${\sum }_{i=1}^n{s}_i{y}_i.$ A smaller a is associated with more distributed weights, and seems more beneficial for this purpose. It is also of note that $\mathrm{logistic}\left(0\right)=1/2$ in SSS regardless of the choice of a, as opposed to 0/1 indicator values in GS. These features make cutting at the extreme ends less likely for SSS. To remedy ECP for least squares splitting in regression trees, Torgo [27] suggested setting a minimum number of observations allowable for either child node, which can be implemented via the minbucket argument in the R function rpart [28]. A similar strategy can be made available in SSS as well. Because Brent’s method asks for a range over which the optimum is sought, it is convenient to define the search interval as $(q_{\gamma },q_{1-\gamma })$, where $q_{\gamma }$ and $q_{1-\gamma }$ denote the sample $\gamma$-th and $(1-\gamma )$-th quantiles of $X,$ respectively. In our implementation, we have included $\gamma =0.02$ as an option.

Figure 3c–f corresponds to the weak signal scenarios, where it can be seen that SSS outperforms GS by a great deal in finding the true cutpoint with all choices of a. The failure of GS in Model B can be partly attributable to the ECP problem. Breiman et al. [1] showed that ECP in LS splitting is a joint result of the Kolmogorov inequality and the law of the iterated logarithm in theory, as further confirmed in a recent article of [29] for other splitting criteria. Ishwaran [29] argued that although ECP is unfavorable for single tree analysis, it is a beneficial feature for random forests, as it maximizes the sample size for a tree to recover from a bad split based on a pure noise predictor. However, we note that the ECP effect can undesirably block out weaker signals, as illustrated here. In addition to the adverse effect of ECP, the splitting statistics computed with GS can be rather turbulent, with considerable variation, especially when the signal grows weak. An illustration from one simulation run is provided in Figure 1. In this scenario, the overall weak signal is likely to be surmounted by some local spikes owing to the large variation. Smoothing can be very helpful here, as in general smoothing preserves the global pattern more than the local ones. Because SSS facilitates a parametric mode of smoothing, in which a plays the role of a smoothing parameter, it is reasonable to expect that SSS helps to remove irregular local fluctuations in order to better isolate the weak signal with a more stable objective function, especially with a smaller a. A relatively smaller a is even preferable in recovering weak signals, as a smaller a helps to smooth out more local erratic patterns. Figure 3g–j illustrates the strong signal scenarios, in which both GS and SSS perform exceptionally well.

We also computed the mean squared error MSE = ${\sum }_{i=1}^{200}{({\widehat{c}}_i-c_0)}^2/200$ as a performance measure for each non-null model configuration. The results are presented in Figure 4. Figure 4a–d depicts the weak signal scenarios and reinforces our earlier conclusion that SSS is more advantageous than GS in these contexts. The figure also reveals a pattern where MSE decreases as a decreases for a certain range of a values, say, $a\ge 5.$ However, excessively small a a values result in poor approximations, leading to diminished performance. This indicates a tradeoff between approximation and smoothness in terms of the choice of a. Figure 4e–h represents the strong signal scenarios, where both GS and SSS perform very well. Notably, SSS outperforms GS in the small sample scenario with $n=50$ across all choices of a considered. In contrast, with larger samples $n=500$ the results are mixed, and the difference between SSS and GS becomes negligible for reasonably large a values (e.g., $a\ge 10$). Additionally, both GS and SSS demonstrate better performance with the balanced cutpoint $c_0=0.5$ compared to the unbalanced cutpoint $c_0=0.8$, as illustrated by the MSE scales in Figure 4.

To further investigate the sensitivity with respect to a, Figure 5 plots $\widehat{c}$ versus a in each simulation run for 100 simulation runs with $\beta =1.$ It can be seen that the lines are essentially flat, except for very small a values. This indicates that the performance of SSS remains fairly stable with respect to a, particularly for reasonably large values. Similar patterns are observed in the weaker signal case ($\beta =0.2$), although these results are not presented here. This suggests that while the scale parameter a functions as a smoothing parameter and further research on its fine-tuning may be valuable, setting a as a constant seems to be a sensible approach.

Additionally, the averaged best cutoff point $\widehat{c}$ at each a overlaps well with the true cutoff point $c_0$, except in Panel (b), where the estimation is biased downwards. GS exhibits a similar downward bias; see panel (i) in Figure 3. Figure 5b highlights the increased difficulty in estimating an unbalanced cutoff point, especially with a small sample size. On the other hand, the notable difference between the mean $\widehat{c}$ values and $c_0$ in Panel (b) can be partially attributed to the nonrobustness of means. If we consider the median $\widehat{c}$ values (represented by the blue line), they align much more closely with the true value $c_0$. As the sample size increases, the bias becomes negligible, as shown in panel (d).

Another alternative for specifying a is to make a adaptive to the sample size $n.$ To investigate this, we reran the experiments by considering sample sizes while varying in {50, 100, 150, …, 1000}. Two strategies for setting a values were tried: constant a in $\{10,30,50\}$, and adaptive choices of a in $\{n,\sqrt{n},ln\left(n\right)\}.$ A total of 200 simulation runs was taken for each scenario and the MSE values were recorded. Because the scale for MSE drops dramatically as n increases, we computed the relative difference in MSE for the selected cutpoints from SSS as opposed to GS; thus, a negative value indicates improvement upon GS. Figure 6 plots the results versus n for each choice of a. It can be seen that SSS constantly outperforms GS with weaker signals $(\beta =0.2)$ and that a smaller a seems preferable. In the stronger signal $(\beta =1.0)$ scenario, the constant choices $a\in \{30,50\}$ work very well, while the adaptive choice $a=\sqrt{n}$ outperforms the other two.

In summary, a larger a yields a better approximation of the threshold effect, while a smaller a offers increased smoothness. With a larger a, the performance of SSS is relatively more stable, as shown in Figure 3, Figure 4 and Figure 5. Conversely, a smaller a can be beneficial for detecting weak signals and addressing the end-cut preference (ECP) problem. Our empirical results demonstrate that the performance of SSS remains stable as a varies, as well as that SSS outperforms GS in most scenarios with a reasonable choice of a. While some optimization of a could be explored, we believe that a reasonable choice of a will suffice for practical applications. Based on our findings, we tentatively recommend fixing a in the range of [20, 60], while the adaptive choice of $a=\sqrt{n}$ also seems advisable.

To investigate computational efficiency, we generated data from Model A in (7) with $\beta =1$ and $n\in \{20,30,\dots ,100,200,300,\dots ,10,000\}.$ Because $x\sim \mathrm{unif}(0,1)$, the number of distinct values that it has increases with $n.$ We considered different constant choices of $a\in \{1,2,\dots ,100\}$ and an adaptive choice $a=\sqrt{n}.$ For greedy search, we used the R function rpart [28] by setting its options rpart.control such that it produces one single split only and all other features, such as competitive and surrogate splits, are suppressed. The implementation of SSS was accomplished solely in R with the optimize function. Figure 7a plots the CPU time (in seconds) that SSS and GS spent on ten simulation runs. It can be seen that SSS consistently outperforms GS. Figure 7b plots the averaged number m of iterative steps in Brent’s method for SSS. It can be seen that m ranges from 7 to 10 in this example. In addition, a smaller a, corresponding to a smoother objective function, leads to a smaller number of iterative steps, as expected. In practice, similar experiments can be performed on specific datasets to determine the approximate threshold at which SSS becomes more advantageous than GS.

2.4. Variable Selection Bias

The variable selection bias problem is deemed inherent in recursive partitioning. By optimally selecting a goodness-of-split measure across all permissible splits from all variables, greedy search gives preference to variables that have more distinct levels or values. As a remedy, Loh [9] (GUIDE) proposed first selecting the most important splitting variable and then finding its best cutoff point via greedy search. In addition to its negligible variable selection bias, the method is advantageous in reducing the computational cost. A similar idea was followed in [30] (MOB). However, determination of the most important splitting variable is an equally if not more difficult problem. A predictor that is important in its threshold effect may not be so in other senses. In GUIDE, residuals from linear models are first grouped. Every predictor is treated or otherwise made categorical to form two-way contingency tables with residual-based grouping. The most important variable is then selected based on Pearson ${\chi }^2$ tests. In MOB, the selection is based on a stability assessment of the estimated coefficients in a linear model. One potential problem with both methods is that the variable selected via linear models may not be the one with the most important threshold effect for the partitioning purpose. Furthermore, these methods lose power in detecting the true threshold signal on continuous or nominal predictors, leading to a different type of bias in variable selection.

In other approaches, the best cutoff point is associated with a maximally selected random variable. Akin to change-point problems (see, e.g., [31,32]), asymptotic inference on maximally selected statistics resorts to standard results in stochastic processes (see, e.g., [10]), while their exact distributions can be determined with permutation-based methods (see, e.g., [12], 2008). Shih and Tsai [11] evaluated these approaches in partitioning data and demonstrated their potential in reducing variable selection bias via simulations. In these approaches, the p-value associated with the best cutpoint for each predictor is first computed and possibly transformed. These p-values are then compared across predictors to determine the best split of the data. Different methods might be used to obtain the p-value for the best cutpoint of each predictor depending on the type of the predictor. Nevertheless, the associated distributions for maximally selected statistics may not be easily available, and the computation could be intensive when either the sample size n or the number of distinct levels K become large. For example, the permutation-based methods are clearly not a good choice for this purpose, as the exact p-value is often computationally infeasible while the approximate p-value does not have the necessary precision in order to compare with p-values obtained from the Brownian bridge for example. Moreover, most of the above-mentioned works on the distribution of maximally selected statistics are applicable or feasible only to continuous or at least ordinal predictors, and not to nominal variables in the strict sense.

The SSS method can help to get around the variable selection bias problem by providing a simple alternative way of evaluating splits. The significance assessment of the best cutoff point on a continuous predictor can be made via its connection to a nonlinear parametric model (3). While other testing procedures are available, it is convenient to use the likelihood ratio test for $H_0:{\beta }_2=0$ in Model (3). Computing the LRT entails estimation of Model (3). This is a separable nonlinear least squares problem that can be efficiently solved by the variable projection method of Golub and Pereyra [33]. The main idea is that $\mathit{\beta }$ can be readily estimated given $c.$ We define matrix $\mathbf{S}=\left[\mathbf{j},\mathbf{s}\right]\in {\mathbb{R}}^{n\times 2}.$ The nonlinear least square criterion is ${\Vert \mathbf{y}-\mathbf{S}\mathit{\beta }\Vert }_2^2$ with $\mathit{\beta }={({\beta }_0,{\beta }_1)}^T.$ For a given c, the estimate of $\mathit{\beta }$ can be obtained via ordinary LS: $\widehat{\mathit{\beta }}={\left({\mathbf{S}}^T\mathbf{S}\right)}^{-1}{\mathbf{S}}^T\mathbf{y}.$ Plugging it into the LS criterion leads to an objective function for c only:

$$\Vert \mathbf{y}-\mathbf{S}\widehat{\mathit{\beta }}{\Vert }^2=\Vert \mathbf{y}-\mathbf{S}{\left({\mathbf{S}}^T\mathbf{S}\right)}^{-1}{\mathbf{S}}^T{\mathbf{y}\Vert }^2={\Vert {\mathbf{P}}_{{\mathcal{C}}^{\perp }\left(\mathbf{S}\right)}\mathbf{y}\Vert }^2$$

(8)

where $\mathcal{C}\left(\mathbf{S}\right)$ denotes the column subspace of matrix $\mathbf{S}$, ${\mathcal{C}}^{\perp }\left(\mathbf{S}\right)$ is the subspace perpendicular to $\mathcal{C}\left(\mathbf{S}\right)$, matrix ${\mathbf{P}}_{{\mathcal{C}}^{\perp }\left(\mathbf{S}\right)}={\mathbf{I}}_n-\mathbf{S}{\left({\mathbf{S}}^T\mathbf{S}\right)}^{-1}{\mathbf{S}}^T$ denotes the orthogonal projector to ${\mathcal{C}}^{\perp }\left(\mathbf{S}\right)$, and $\Vert ·\Vert$ denotes the Euclidean norm unless otherwise specified. Estimates of c can be obtained by minimizing (8), which is equivalent to maximizing $\tilde{\Delta }𝚤\left(c\right)$ in (5). A popular implementation is to combine the Gauss–Newton method with variable projection, as considered by Kaufman [34] and summarized in Osborne [35]. The algorithm is iterative with an approximately quadratic convergence rate. The corresponding pseudocode for estimating Model (3) is provided in Appendix A.

In our approach, the best cutoff $\widehat{c}$ is efficiently obtained by optimizing (8) or the approximated splitting statistic directly. Given $\widehat{c},$ Model (3) becomes linear and $\widehat{\mathit{\beta }}$ can be obtained via OLS. The LRT statistic (see, e.g., Section 12.4 of [21], 2003) is provided by

$$\mathrm{LRT}=nlog{\displaystyle \frac{{\mathbf{y}}^T({\mathbf{I}}_n-{\mathbf{jj}}^T/n)\mathbf{y}}{\Vert \mathbf{y}-\mathbf{S}\left(\widehat{c}\right)\widehat{\mathit{\beta }}{\Vert }^2}},$$

where $\mathbf{S}\left(\widehat{c}\right)$ denotes the $\mathbf{S}$ matrix evaluated at $\widehat{c}.$ This test compares Model (3) with the null model $y_i={\beta }_0+{\epsilon }_i$. Because a is fixed, the resulting LRT statistic follows a ${\chi }^2\left(2\right)$ distribution under $H_0.$ Alternatively, an F test can be used (see Section 5 of Seber and Wild [21]). This result suggests a conjecture that there may be a connection between the ${\chi }^2\left(2\right)$ distribution and the asymptotic distribution derived from the Brownian bridge [10], which could be further investigated in future research. We define the associated logworth

$$\mathrm{logworth}=-{log}_{10}Pr\left\{{\chi }^2\left(2\right)\ge \mathrm{LRT}\right\}$$

as the minus base-10 logarithm of the resulting p-value when referring to a ${\chi }^2\left(2\right)$ distribution. The logworth provides a coalescent measure of the strength for the best cutoff point of each predictor regardless of the variable type or the methods used to evaluate its best cutpoint. The best split of the data is then the one with the maximum logworth among all predictors.

We next discuss how to handle other types of predictors in combination with SSS. For a 0–1 binary variable, the LRT for $H_0:{\beta }_1=0$ in Model (1) with $c=0.5$ can be used to evaluate the single split. The corresponding p-value, together with the associated logworth, can be computed by referring to the ${\chi }^2\left(1\right)$ distribution. For an ordinal variable with K distinct values, we process with SSS by treating as if it is continuous. However, when $3\le K\le 10$ is small, it is possible to resort to the maximally selected statistics approaches. Several approaches are available for computing the associated p-value. A method that combines the change-point approach of Hawkins [36] and Genz [37], as implemented in the R package maxstat [38], is particularly appealing for the small K scenario in terms of computational feasibility and efficiency.

For a nominal variable X with K distinct levels, we first ‘ordinalize’ its levels by sorting its group means and then proceed by treating it as if it were continuous. This standard strategy in GS greatly reduces the total number of splits to be evaluated, and is theoretically justified (see CART, Section 9.4, [1]); however, both its steps introduce overoptimism or bias in terms of the association between Y and X. The ordinalization step may make an irrelevant X artificially correlated with Y, while the maximal selection of splitting criteria over the ordinalized levels of X induces additional overoptimism. Clearly, SSS can be used in the latter step to find the best cutpoint and help avoid bias owing to maximal selection. However, the overoptimism stemming from the ordinalization step remains an issue. In general, how to splitting nominal variables in regression trees without bias is not satisfactorily addressed in the extant literature, and further research is warranted.

To demonstrate this, we simulate data from a null model where the response Y has nothing to do with any of the nine either predictors $\{X_1,X_2,\dots ,X_9\}$. For the above discussed reasons, we have restricted ourselves to ordinal or continuous predictors only. The numbers of distinct values for $\{X_1,X_2,\dots ,X_9\}$ were set as {2, 3, 4, 5, 10, 20, 50, 100, 500}, respectively. Two sample sizes $n=50$ and $n=500$ were considered. Figure 8 plots the frequency of each variable being selected for splitting data by SSS (combined with maximally selected statistics) and GS. For SSS, we tried different choices of a in {10, 30, and 50}. It can be seen that greedy search suffers severely from selection bias, while SSS substantially alleviates the problem. Different choices of a for SSS yield very similar results, with a slight preference to smaller a values.

For further exploration, we consider another trial where data were generated from

$$\mathrm{Model} \mathrm{B}:y=1+\beta ·I\{x_1\le 0.5\}+\epsilon \mathrm{with}\epsilon \stackrel{IID}{\sim }\mathcal{N}(0,1),$$

(9)

where $X_1$ is binary 0–1 valued and we include a second variable $X_2\sim \mathcal{N}(0,1).$ In this model, Y is only related to the binary predictor $X_1$; nevertheless, the unrelated predictor $X_2$ has distinct values of $K\approx O\left(n\right).$ We then allowed the coefficient $\beta$ to vary from $\{0.01,0.02,\dots ,1\}.$ We recorded the percentage of correct selection of the true split $X_1\le 0.5$ made by GS and SSS. For SSS, we also tried different choices of $a\in$ {10, 30, 50, 100}. Figure 9 plots the results against $\beta$ for two sample sizes, $n=50$ and $n=500.$ SSS clearly outperforms GS in picking up the correct split. For weak signals, it is more likely for GS to have spurious splits introduced into the tree structure. As expected, both perform better as the signal grows stronger or the sample size grows larger. For SSS, a smaller a tends to yield a higher percentage of correct selection, yet with seemingly unimportant differences.

In the above simulation, we included only one binary predictor in order to better isolate the specific challenges of variable selection bias. As one reviewer noted, the effects of multiple predictors on the response can complicate the delineation of a variable’s impact. While future research with more extensive simulations is necessary to address this issue, we would like to emphasize that the approach used in recursive partitioning is somewhat analogous to coordinate descent in optimization, where each variable is handled separately. This allows us to focus on the individual impact of each predictor at each split.

2.5. Tree Performance

Although our primary focus is on one single split of data, it is interesting to see how SSS compares with GS in whole-tree performance as well. To this end, we implemented SSS using the user-written splitting functions [8] in rpart. While this may not be a good idea with regard to computing speed, it allows us to utilize many useful features available in rpart [28], such as surrogate splits for missing value treatment.

2.5.1. Simulation Studies

We first consider simulation studies. Data were generated from the following models:

$$\begin{array}{cc}\hfill \mathrm{Model} \mathrm{C}:& y=0.5+0.5·x_1+0.5·I\{x_2\le 0.5\}+\epsilon \end{array}$$

(10)

$$\begin{array}{cc}\hfill \mathrm{Model} \mathrm{D}:& y=2+2·x_1·I\{x_2\le 0.5\}+\epsilon \end{array}$$

(11)

with $\epsilon \stackrel{IID}{\sim }\mathcal{N}(0,1)$. Models C and D are expected to have four and three terminal nodes and correspond to weaker and stronger signals, respectively. Five predictors were generated, with $X_1$ being binary 0–1 and $X_2$, …, $X_5$ from $\mathrm{unif}(0,1),$ where only $X_1$ and $X_2$ are related to Y in either an additive or interaction form.

There are a few approaches available for determining the best tree size. To facilitate the comparison, a unified ‘test sample’ approach (see, e.g., Section 11.4 of CART, [1]) was taken for both SSS and GS. Specifically, we first constructed a large initial tree based a training sample ${\mathcal{L}}_1$ of size n and pruned it to obtain a nest sequence of subtrees. Then, we generated an independent validation sample ${\mathcal{L}}_2$ of the same size n and applied each subtree in the sequence to predict the responses in ${\mathcal{L}}_2$. The subtree yielding the minimum prediction mean squared error (PMSE) was identified as the ‘best’. Two sample sizes $n=100$ and $n=500$ were considered in our experiments. For SSS, the scale parameter was fixed at $a=50.$

To compare the predictive accuracy, another independent testing sample ${\mathcal{L}}_3$ of size $n^{\prime }=5000$ was generated from the same model. In every simulation run, the best tree structure we identified was applied to ${\mathcal{L}}_3$ for prediction. The resulting PMSE $={\sum }_{i\in {\mathcal{L}}_3}{(y_i-{\widehat{y}}_i)}^2/n^{\prime }$ was recorded. A total of 200 simulation runs was carried out, and the averaged PMSE is reported in

Table 1. Comparison of SSS vs. GS for whole tree performance in regression problems. For simulated data, the table reports the final tree size and PMSE, both averaged over 200 simulation runs. For benchmark datasets, the performance measures are based on one single tree procedure. The computing time refers to the duration (in seconds) taken to construct the initial tree using the entire dataset. For simulated data, the reported time is the average from 20 simulation runs.

	Data			SSS				GS
	Set	$\mathit{n}$	$\mathit{p}$	Size	PMSE	Time		Size	PMSE	Time
Simulated	Model C	100	5	2.295	1.155	0.0715		1.800	1.163	0.0350
		500	5	4.080	1.036	0.0790		3.850	1.045	0.2365
	Model D	100	5	3.385	1.121	0.0675		3.380	1.136	0.0285
		500	5	3.135	1.060	0.0725		3.190	1.058	0.2315
Benchmark	airfoil	1503	5	35	19.795	0.20		36	19.227	1.14
	auto	398	7	9	23.940	0.12		9	24.606	0.22
	concrete	1030	8	59	58.326	0.27		42	71.689	1.30
	energy: heating	768	10	35	0.978	0.11		33	0.968	0.40
	energy: cooling	768	10	33	5.709	0.17		34	5.790	0.34
	wiki4HE	907	48	9	0.401	1.57		9	0.434	4.22

. Note that ${\mathcal{L}}_3$ remained the same across all simulation runs. Table 1 also reports the ‘best’ tree size averaged over the simulation runs. It can be seen that SSS consistently outperforms GS in detecting weak signals. In the weak signal scenario (Model C), SSS leads to a final tree size that is closer to the true one of 4 and yields a smaller averaged PMSE. SSS is also highly competitive for Model F with stronger signals.

2.5.2. Real Data Examples

For real-world example illustrations, we first considered the 1987 baseball salary data, which have been widely analyzed in the literature. The final CART tree, selected with the 0-SE rule, has eight terminal nodes, as plotted in Figure 10a. Loh [9] commented that the two splits based on team86 and team87 (highlighted with dashed lines) are hard to interpret and may be attributable to the selection bias of greedy search. The 1-SE rule yields a final tree of four terminal nodes, which is a complete full tree of depth 3. Figure 10b plots the final SSS tree with $a=50,$ which has five terminal nodes. This tree structure was determined following the tree methodology based on the goodness-of-split [7], which combines split-complexity pruning with bootstrap-based bias correction. Of note is that the two categorical variables concerning baseball team membership no longer show up in the final SSS tree. The yrs variable is deemed important in the GUIDE tree, which can be confirmed via variable importance ranking in random forests as well. It does not show up in the CART tree, but does in the SSS tree. We also conducted a 10-fold cross-validation to compare their predictive performance. The mean squared prediction error measures for the SSS tree, 0-SE CART tree, and 1-SE CART tree were 0.1867, 0.1194, and 0.1915, respectively. Thus, the final SSS tree is highly competitive in terms of predictive accuracy.

For further illustration, we selected five more benchmark datasets from the Regression category at UCI Machine Learning Repository (http://archive.ics.uci.edu/ml/, accessed on 20 September 2024). These were Airfoil Self-Noise (airfoil), Auto MPG (auto), Concrete Compressive Strength (concrete), Energy Efficiency (energy), and Teaching Uses of Wikipedia (wiki4HE). As an “off-the-shelf” tool, tree analysis requires little data preparation. Of brief note, there are two responses in energy, namely, heating load and cooling load, which were analyzed separately. In wik4HE, we summed up the item scores for the last category (“experience with Wikipedia usage”) and used it as the response variable. All missing values were handled using the surrogate splits approach.

The following paradigm was taken to analyze each dataset. We randomly partitioned the dataset into three parts, the training sample ${\mathcal{L}}_1$, validation sample ${\mathcal{L}}_2$, and testing sample ${\mathcal{L}}_3$, with approximately equal proportions 1:1:1. A large tree was grown and pruned using ${\mathcal{L}}_1$ and the best subtree was selected according to the minimum PMSE based on ${\mathcal{L}}_2$. The best tree size and the PMSE measure recomputed from ${\mathcal{L}}_3$ are reported in Table 1. It can be seen that SSS generally compares favorably to GS in predictive performance in most cases, although not always. This can be attributed to its superior ability to detect weak signals and its immunity to issues such as end-cut preference and variable selection bias. Despite the potential advantages of SSS over GS in various scenarios, its superiority in prediction accuracy ultimately depends on the specific dataset being analyzed.

2.5.3. Computing Time Comparison

To compare computing times, we recorded the duration (in seconds) required to construct the initial large tree using the entire dataset. To ensure a fair comparison, we implemented both SSS and GS with the user-defined splitting function from the R package rpart [8]. The maximum tree depth was set to 10 and the minimum node size was 20. All programs were executed on a ThinkPad workstation equipped with a 2.40 GHz CPU and 36 GB of memory. The results are presented in Table 1. For simulated data, the reported times are the average from 20 simulation runs.

Table 1 shows that SSS outperforms GS across all benchmark datasets. However, for simulated data, GS is faster than SSS when $n=100$, although it quickly becomes slower as the sample size increases to $n=500.$ It is noteworthy that while the increase in sample size has a minimal impact on the computing time of SSS, it significantly influences the performance of GS. This intriguing result suggests that adaptive use of SSS and GS in practical implementations is advisable depending on factors such as sample size n and the number m of distinct values of a covariate. When n is small or m is limited, GS tends to be faster than SSS. These scenarios are commonly encountered as the tree continues to grow. Thus, a threshold might be established for n and m to determine when to use GS instead of SSS, applying GS when either n or m falls below the threshold.

3. SSS for Classification Trees

SSS can be extended to classification trees in a similar manner. Consider data $\{(y_i,{\mathbf{x}}_i):i=1,\dots ,n\}$ with binary responses $y_i\in \{0,1\}$. Commonly-used node impurities for classification trees include the misclassification error rate, Gini index, and entropy. The Gini index corresponds to treating binary Y as if continuous and applying OLS directly with the linear regression model. With the entropy measure, minimizing the within-node impurity is equivalent to maximizing the LRT for testing $H_0:{\beta }_1=0$ in a logistic model:

$$\mathrm{logit}\left({\pi }_i\right)=log{\displaystyle \frac{{\pi }_i}{1-{\pi }_i}}={\beta }_0+{\beta }_1\delta (x_i;c)$$

(12)

where ${\pi }_i=Pr(y_i=1|{\mathbf{x}}_i).$

Among the three common measures, the misclassification error rate is not recommended due to its serious defects, as discussed in Section 4.1 of CART [1]. The Gini index and entropy behave similarly. They both work well with strong signals and generalize well to multiclass responses, but suffer from the end-cut preference (ECP) problem. To remedy this, Morgan and Messenger [39] introduced a ‘delta splitting rule’ by replacing the squared Euclidean distance with the absolute deviation in the Gini index; however, this is slightly harder to compute. Breiman et al. [1] considered a ‘twoing rule’ in CART, which was shown to be equivalent to the delta splitting rule.

Just as earlier, replacing the threshold term $\delta (x_i;c)$ with an SSS term yields a nonlinear logistic model:

$$\mathrm{logit}\left({\pi }_i\right)={\beta }_0+{\beta }_{1}s(c;x_i).$$

(13)

To avoid optimization with multiple parameters, it is preferable to first form an objective function by approximating a legitimate splitting statistic and then optimizing it with respect to the cutpoint c only. This strategy, along a similar line to Neyman [40], shows considerable yield in terms of computational simplicity. For a node t, let $n_1={\sum }_{i\in t}{y}_i$ denote the number of observations with $y=1$, $n_L$ the number of observations in the left child node, and $n_{L1}$ the number of observations with $y=1$ in the left child node. The following proposition expresses the reduction in impurity $\Delta 𝚤$ in classification trees in terms of $n_L$ and $n_{L1}$ only.

Proposition 2. 

Let ${\pi }_t$ denote the proportion of observations with $y=1$ in node t. Let $\Delta 𝚤$ be the reduction in node impurity when node t is split into $t_L$ and $t_R$. With entropy impurity measure $𝚤\left(t\right)=-{\pi }_{t}log{\pi }_t-(1-{\pi }_t)log(1-{\pi }_t)$, we have

$$\begin{array}{ccc}\hfill \Delta 𝚤& \widehat{=}& n_{L1}log{n}_{L1}+(n_L-n_{L1})log(n_L-n_{L1})+(n_1-n_{L1})log(n_1-n_{L1})\hfill \\ & & +\{(n-n_1)-(n_L-n_{L1})\}log\{(n-n_1)-(n_L-n_{L1})\}\hfill \\ & & -n_{L}log{n}_L-(n-n_L)log(n-n_L)\hfill \end{array}$$

(14)

up to some constant. With Gini index $𝚤\left(t\right)={\pi }_t(1-{\pi }_t)$, we have

$$\Delta 𝚤\widehat{=}-{\displaystyle \frac{n_{L1}(n_L-n_{L1})}{n_L}}-{\displaystyle \frac{(n_1-n_{L1})\{(n-n_L)-(n_1-n_{L1})\}}{n-n_L}}.$$

(15)

Proposition 2 provides the grounds for an efficient way of estimating the best cutoff point c via SSS. To do so, we need to approximate only $n_L={\sum }_{i=1}^n\delta (x_i;c)$ with ${\tilde{n}}_L={\sum }_{i=1}^n{s}_i$ and $n_{L1}={\sum }_{i=1}^n{y}_i\delta (x_i;c)$ with ${\tilde{n}}_{L1}={\sum }_{i=1}^n{y}_i{s}_i$ in either (14) or (15). The approximated $\Delta 𝚤$ becomes a smooth objective function of c alone; we denote it as $\tilde{\Delta }𝚤\left(c\right).$ The best cutoff point $\widehat{c}$ can be found as $\widehat{c}=arg{max}_c\tilde{\Delta }𝚤\left(c\right).$ This is again a one-dimensional optimization problem, and can be efficiently solved with the Brent [23] method.

SSS is flexible with the general idea of approximation and optimization. It is obvious that any other two-sample test statistic, such as Pearson ${\chi }^2$, can be used in this approach as well. An analogous strategy for dealing with variable selection bias, as in Section 2.4, can be used in classification trees. Estimation of the nonlinear logistic model (13) can be carried out by combining the iteratively re-weighted least squares (IRWLS) method with the Gauss–Newton plus variable projection method discussed earlier. The algorithm is described in Appendix A. The LRT for testing $H_0:{\beta }_1=0$ with the nonlinear logistic model (13) can be shown to have the following form:

$$\mathrm{LRT}=2\sum _{d\in \{L,R\}}\sum _{k=0,1}{\tilde{n}}_{dk}log{\displaystyle \frac{{\tilde{n}}_{dk}}{E_{dk}}},\mathrm{with}{E}_{dk}={\displaystyle \frac{n_k{\tilde{n}}_d}{n}},$$

which is similar to its counterpart in the linear logistic model except that the counts $n_{dk}$ involving the split are now replaced with their approximated surrogates ${\tilde{n}}_{dk}.$ The statistical significance of LRT refers to a ${\chi }^2\left(2\right)$ distribution except for the scenario with binary predictors where LRT $\stackrel{H_0}{\sim }{\chi }^2\left(1\right).$ Alternatively, Pearson ${\chi }^2$ of form ${\sum }_{d\in \{L,R\}}{\sum }_{k=0,1}{({\tilde{n}}_{dk}-E_{dk})}^2/E_{dk}$ can be used instead. Comparing across predictors, each with its best cutoff point, the split that yields the smallest p-value or the largest logworth is selected to partition the data. As before, both ordinal variables (with small $3\le K\le 10$ distinct levels) and nominal variables can be treated separately, resulting in a hybrid approach. In particular, Boulesteix [41,42] provided explicit formulas for computing the exact p-values associated with the maximally selected ${\chi }^2$ variables in the classification scenario for both nominal and ordinal variables, as implemented in the R package exactmaxsel. However, the combinatorial approach of [42] becomes computationally infeasible when both n and K are large, a scenario that warrants further research.

3.1. Simulation Studies

To illustrate, we simulated data from a logistic model

$$\mathrm{Model}\mathrm{E}:Pr\{y=1\}=\mathrm{logistic}\left(-0.5+I\{x\le 0.5\}\right),$$

(16)

with $x\sim \mathrm{Unif}[0,1].$ Two sample sizes ($n=50$ and $n=500$) and 500 simulation runs were considered. We also experimented with different choices of a {5, 10, 15, …, 100 } with SSS. Figure 11a,b provides the empirical density curves of $\widehat{c}$ identified by the entropy-based GS and SSS for $n=50$ and $n=500$, respectively. It can be seen that GS and SSS have very similar performances for the large ($n=500$) samples. When $n=50$, SSS substantially outperforms GS. The distribution of GS shows two bumps at either end. This can be explained again by the end-cut performance problem. To see this, Figure 11c presents a smooth scatterplot of the best cutoff point found by GS versus that found by SSS with $a=30$ in the smaller samples scenario with $n=50$. Because there are many overlapping points, we have smoothed the scatterplot in such a way that the frequencies of overlapped points are represented in the electromagnetic spectrum. The zigzag pattern clearly indicates that GS selects more cuts that are near to either end of the data range. This is because the signal strength of Model E in (16) is not strong enough to overpower the inherent end-cut tendency, while SSS largely avoids this issue. To this end, a smaller a seemingly provides even better performance for SSS. Figure 11d plots the MSE with respect to the true cutoff point $c_0=0.5$ versus the a value. SSS has very comparable performance with GS in the case of larger samples, but triumphs with smaller samples for all choices of $a.$

To investigate the variable selection bias, we simulated data from a null logistic model with similar settings to those in Section 2.4. Specifically, each dataset contained nine features $\{X_1,X_2,\dots ,X_9\}$, with the possible number of distinct values being equal to {2, 3, 4, 5, 10, 20, 50, 100, 500}, respectively. Three choices of a, namely, {10, 30, 50}, were experimented with for SSS. The bar plots in Figure 12 provide a comparison of the selection frequencies from SSS and GS, where the hybrid method with SSS largely eliminates the selection bias of GS.

3.2. Example: Credit Card Default Data

To illustrate the use of SSS in classification trees, we considered the credit card default dataset [43] available at the UCI Machine Learning Repository. The primary objective of the study was to predict the default status of credit card holders based on 23 predictor variables related to their financial information, such as business financial statements, customer transactions, and repayment records. The dataset contains 30,000 observations, presenting a relatively large-scale problem. More details about the dataset can be found in [44].

As before, we randomly partitioned the data into three sets, the training sample ${\mathcal{L}}_1$, validation sample ${\mathcal{L}}_2$, and test sample ${\mathcal{L}}_3$, each containing one-third of the entire data. A large initial tree was grown and pruned using entropy as the impurity measure with ${\mathcal{L}}_1$. The final tree structure was determined through validation on ${\mathcal{L}}_2$ based on the misclassification error. Finally, the resulting tree structure was evaluated using the test sample ${\mathcal{L}}_3$. Two performance metrics were computed: the misclassification error rate, and the area (AUC) under the receiver operating characteristic (ROC) curve.

Figure 13 displays the final tree structures for both GS and SSS. The GS tree contains eight terminal nodes, while the SSS tree has sixteen terminal nodes. The splits in the top levels are essentially the same for both models. Based on the test sample, the misclassification error rates are 0.1792 for the SSS tree and 0.1798 for the GS tree, indicating that they are quite similar. These results align closely with those reported in [44]. Figure 13 also shows the ROC curves. The SSS tree achieves an AUC value of 0.760, significantly surpassing the AUC value of 0.692 for the GS tree. This indicates that the SSS tree is more effective at distinguishing between credible and non-credible cardholders.

4. Discussion

In this paper, we have explored an intuitive alternative to greedy search in decision trees by approximating step functions with smooth sigmoid functions. We demonstrate that the proposed SSS method can be reformulated as a one-dimensional optimization in various recursive partitioning procedures, resulting in improved computational efficiency. SSS facilitates a parametric smoothing or regularization of the erratic goodness-of-split measures computed in greedy search. The smoothing effect provides SSS with remarkable immunity to the end-cut preference problem in GS. By smoothing out local peaks that stem from large variation and end-cut preference, SSS outperforms GS substantially in identifying weak signals. Moreover, its natural connection with nonlinear regression renders a way of alleviating the variable selection bias problem. To this end, our ongoing research efforts are being directed towards investigating how valid statistical inferences can be made based on the best split on nominal predictors and fully addressing overoptimism owing to ‘ordinalization’.

We have confined our attention to the use of SSS in settings where the splitting statistic is of explicit form. This allows us to formulate identification of the best cutoff point as a one-dimensional smooth optimization problem. To solve the optimization, we have chosen the Brent [23] method for its easy accessibility and fast and stable performance. Other algorithms can be used instead to provide faster convergence rates or carry out global optimization. SSS can also be useful in scenarios without a closed form for the splitting criterion or for updating (see Chen et al. [45] and Su et al. [46] for examples where evaluation of each cutpoint involves iterations). The approximation idea of SSS can be incorporated into the associated models, which can then be estimated with the aid of separable least squares, as discussed in Appendix A.

Our discussion has mainly focused on searching for the best cutoff point of a predictor and one single split of data. We have demonstrated that SSS can be a superior alternative to GS, with several appealing advantages. Our results suggest that implementation of many recursive partitioning methods and their extensions can benefit from the SSS approach. Tree-based methods have been extended to accommodate various types of data, and their applications have expanded to serve a wider range of analytical purposes. A key idea is to formulate the search for the optimal cutoff point as a one-dimensional optimization problem whenever possible. For example, incorporating SSS into ensemble learning algorithms may enhance predictive accuracy. The use of SSS in random forests [6,47] may lead to better variable importance ranking. Future research to see how the benefits of SSS carry over to various recursive partitioning settings and impact their final results is warranted. Additionally, SSS can be especially beneficial for large-scale or high-dimensional data. Tree models process variables individually and perform automatic variable selection, making them well suited for high-dimensional analysis even in cases where $p\gg n.$ However, increased dimensionality can lead to greater computational burden. In these scenarios, the efficiency of SSS makes its incorporation especially advantageous, helping to alleviate the challenges associated with high dimensionality. Implementing SSS in a lower-level programming language is another promising future research avenue regarding computational efficiency concerns.