# Enrichment Score: a better quantitative metric for evaluating the enrichment
capacity of molecular docking models

Ian Scott Knight
<EMAIL_ADDRESS>Slava Naprienko
<EMAIL_ADDRESS>John J. Irwin
<EMAIL_ADDRESS>
(October 2022)

###### Abstract

The standard quantitative metric for evaluating enrichment capacity known as
LogAUC depends on a cutoff parameter that controls what the minimum value of
the log-scaled x-axis is. Unless this parameter is chosen carefully for a
given ROC curve, one of the two following problems occurs: either (1) some
fraction of the first inter-decoy intervals of the ROC curve are simply thrown
away and do not contribute to the metric at all, or (2) the very first inter-
decoy interval contributes too much to the metric at the expense of all
following inter-decoy intervals. We fix this problem with LogAUC by showing a
simple way to choose the cutoff parameter based on the number of decoys which
forces the first inter-decoy interval to always have a stable, sensible
contribution to the total value. Moreover, we introduce a normalized version
of LogAUC known as enrichment score, which (1) enforces stability by selecting
the cutoff parameter in the manner described, (2) yields scores which are more
intuitively meaningful, and (3) allows reliably accurate comparison of the
enrichment capacities exhibited by different ROC curves, even those produced
using different numbers of decoys. Finally, we demonstrate the advantage of
enrichment score over unbalanced metrics using data from a real retrospective
docking study performed using the program DOCK 3.7 on the target receptor
TRYB1 included in the DUDE-Z benchmark.

##  1 Introduction

In the field of computational drug discovery, the technique of molecular
docking is widely used to predict how small molecules will interact with
larger molecules (typically receptor proteins)[1]. The more favorable a given
interaction is predicted to be, the higher the small molecule’s expected
binding affinity with the receptor. Generally, a higher predicted binding
affinity is taken to mean that that small molecule is more likely to actually
bind to that receptor. A small molecule that actually binds to a receptor is
called a ligand.

A given docking program constitutes a docking model that estimates the binding
affinity of candidate small molecules to a receptor of interest by performing
a constrained minimization of the model’s scoring function[1]. Many scoring
functions have been proposed[1][2][3][4][5][6][7], ranging from force field
calculations to knowledge-based methods and even machine learning methods.

The quality of a docking model is measured by how well the model can
distinguish ligands from non-ligands[1][19]. This is done by examining how the
model scores a dataset of small molecules consisting of (1) known ligands
(called actives) and (2) molecules known or expected not to bind to the
receptor (known as decoys). The practice of using a dataset of actives and
decoys to measure how well a docking model can distinguish between ligands and
non-ligands is known as retrospective docking[8]. This is the principal model
validation technique available for docking models. Once properly validated, a
docking model may be used to perform prospective docking, which means scoring
molecules of unknown activity[8]. A prospective docking screen usually
involves the scoring of many millions or even billions of molecules[19].

Decoys can be generated for a given receptor in a number of ways[19][11], but
there are also several established datasets of actives and decoys available as
benchmarks[8][9][10]. Decoy sets are often designed to be particularly
difficult to discriminate from actives for certain targets, such that if a
given docking model can successfully discriminate between them in the setting
of retrospective docking then this warrants a stronger belief in the model’s
usefulness as a tool for enrichment in the setting of prospective docking[19].

An important goal in the field of computational drug discovery is to find a
way to quantitatively measure the capacity of a given docking model to enrich
a set of molecules by reliably predicting mostly favorable interactions for
actives and mostly unfavorable interactions for decoys[14]. There are a number
of enrichment metrics that have been developed[14][16][18]. For example, the
metric known as enrichment factor (EF) captures the idea of enrichment by
equating it to the proportion of actives present in some top fraction of best
scoring molecules[15]. Quantitative metrics like EF also allow researchers to
compare the performance of different docking models on the same dataset of
actives and decoys[16].

The metric known as LogAUC[13] is one of the most popular metrics for
evaluating the quality of molecular docking models[9][10][19][21]. However, it
comes with a significant drawback: it depends on a cutoff parameter[13][17].
This cutoff parameter controls what the minimum value of the log-scaled x-axis
is, which is a mandatory consequence of the choice to use the logarithmic
scale. Unless this parameter is chosen carefully depending on the dataset, one
of the two following situations occurs: either (1) some fraction of the first
inter-decoy intervals of the ROC curve are simply thrown away and do not
contribute to the metric at all, or (2) the very first inter-decoy interval
contributes too much to the metric (even compared to the contribution of the
second inter-decoy interval), and this comes at the expense of all inter-decoy
intervals following the first interval.

We fix this problem with LogAUC by showing a simple way to choose the cutoff
parameter based on the number of decoys which forces the first inter-decoy
interval to always have a stable, sensible contribution to the total value.
Moreover, we introduce a normalized version of LogAUC known as enrichment
score, which (1) enforces stability by selecting the cutoff parameter in the
manner described, (2) is more intuitively meaningful, and (3) allows reliably
accurate comparison of the enrichment capacities exhibited by different ROC
curves.

##  2 Receiver Operating Characteristic

The statistical tool known as the receiver operating characteristic (AKA “ROC
curve”) is a plot of how the true positive rate and false positive rate of an
indexed family of binary classifiers change as the index is varied. The ROC
curve itself is a step function:

$\operatorname{ROC}(x)=\begin{cases}y_{1},&\text{if $k_{0}<x\leq k_{1}$,}\\\
y_{2},&\text{if $k_{1}<x\leq k_{2}$,}\\\ &\vdots\\\ y_{n},&\text{if
$k_{n-1}<x\leq k_{n}$,}\\\ \end{cases}$

where $y_{i}$ are true positive rates and $k_{i}$ are false positive rates.

Figure 1: An illustration of how to interpret the receiver operating
characteristic.[22]

To compute the ROC, we usually use the family of binary classifiers indexed on
the real line where each particular classifier’s index $i$ specifies a cutoff
value below which it predicts $\operatorname{True}$ (“active”) and otherwise
predicts $\operatorname{False}$ (“decoy”):

$\text{classifier}_{i}(x)=\begin{cases}\text{True},&\text{if $x<i$}\\\
\text{False},&\text{if $x\geq i$}\end{cases}$

Now that we have our family of binary classifiers, we determine the true
positive rate and false positive rate of each classifier in the family using
the following dataset collected for a given set of molecules: (1) each
molecule’s docking model score $s_{i}\in\mathbb{R}$ and (2) its true
active/decoy status $b_{i}\in\mathbb{B}$, where
$\mathbb{B}=\\{\operatorname{True},\operatorname{False}\\}$ is the set of
Booleans.

This description the family of binary classifiers typically used for computing
ROC assumes that the docking model being used indicates more favorable
predictions (i.e. those more likely to be “active”) by assigning lower scores.
For example, if the docking model score is an approximation of the energy of
the physical system formed by the scored molecule bound to the receptor, then
a more negative score will indicate a more negative energy, which is more
favorable. This is why the above classifier function predicts
$\operatorname{True}$ for input below its index rather than above. If the
docking model being used indicates more favorable predictions by assigning
higher scores rather than lower scores, then the classifier function should be
modified to instead predict $\operatorname{True}$ for input above its index.

##  3 Linear-log plot and LogAUC

Calculated using a given docking model’s predictions for a given set of
molecules, the ROC curve is often used to analyze the trade-off between true
positive rate and false positive rate in an indexed family of classifiers as
its index is varied. While this visual presentation is useful, it is desirable
to capture the enrichment capacity of a docking model in a single quantitative
metric. The metric known as “LogAUC” is most commonly used for this purpose.
In this section, we demonstrate the mathematical derivation of LogAUC and show
how to calculate it explicitly for an ROC curve.

First we define the literal logarithmic area under the curve for an arbitrary
function $f$ in the form of an explicit integral. Let $f\colon[0,1]\to[0,1]$
be an integrable function (any function you meet in practice is integrable).
We consider the plot of $f$ with the $x$-axis on a logarithmic scale and the
$y$-axis unchanged (on a linear scale). This kind of plot is known as a
linear-log plot. This logarithmic scaling of the $x$-axis is used in order to
preferentially weight the early behavior of the function in determining the
total value of LogAUC. In this way, LogAUC uses the well-understood
mathematical operation of the logarithm to formalize the idea that earlier
enrichment is significantly more important than later enrichment; in fact,
because the logarithm is the inverse operation of exponentiation, the choice
to use logarithmic scaling implies that earlier enrichment is, in particular,
exponentially more important than later enrichment.

Since $\ln(x)$ is undefined at $x=0$, we need to decide an offset $0<a<1$ from
zero to serve as the lowest value of $x$. This means we actually consider $f$
on $[a,1]$. The function $f(x)$ on $[a,1]$ on a logarithmic scale turns into
$f(e^{x})$ on $[\ln(a),0]$. (We remind that $\ln(a)<0$ for $0<a<1$.) The area
$A(f;a)$ of $f$ on a linear-log plot with the offset $a$ is given by

$\displaystyle A(f;a)$ $\displaystyle=\int_{\ln(a)}^{0}f(e^{x})\,dx$ (1)
$\displaystyle=\int_{a}^{1}f(y)\,\frac{dy}{y}.$

The maximum area is attained when $f=1$ on the entire interval $[0,1]$, giving
$A_{\operatorname{max}}=A(f=1;a)=-\ln(a)$.

The usual approach at this moment is to divide this area of function in the
linear-log plot by its maximum possible area to get the quantity known as
LogAUC:

$\displaystyle\operatorname{LogAUC}(f;a)$
$\displaystyle=\frac{A(f;a)}{A(1;a)}$ (2)
$\displaystyle=\frac{1}{-\ln(a)}\int_{a}^{1}f(x)\,\frac{dx}{x}.$ (3)

For example, for the random classifier $\operatorname{rand}(x)=x$, we get

$\displaystyle\operatorname{LogAUC}(\operatorname{rand};a)$
$\displaystyle=\frac{1}{A(1;a)}\int_{a}^{1}x\,\frac{dx}{x}$
$\displaystyle=\frac{1-a}{-\ln(a)}.$

According to the literature [13, 9, 10, 19, 21], there is the standard offset
value $a=$1.0\text{\times}{10}^{-3}$$ that is typically used by practitioners
who employ LogAUC to evaluate the enrichment capacity of their docking models.

If $a=$1.0\text{\times}{10}^{-3}$$, then

$\displaystyle\operatorname{LogAUC}(\operatorname{rand};a)$
$\displaystyle=\frac{1-0.001}{\ln(1)-\ln(0.001)}$ $\displaystyle\approx
0.144622,$

which matches the value specified in the literature for the LogAUC of the
random classifier[10].

But notice that the value $\operatorname{LogAUC}$ of the random classifier
depends on the parameter $0<a<1$. For example, the LogAUC of the random
classifier with the offset $a=10^{-4}$ $\approx 0.108563$. It shows how
unintuitive $\operatorname{LogAUC}$ is as a metric for enrichment capacity
because (1) its value for the random classifier depends on the offset
parameter, and (2) it does not provide an easy way to compare a given
classifier against the random classifier.

This makes it seem desirable to use a consistent value for $a$ across all
datasets, so as to allow comparison of their LogAUC measurements.

###  3.1 Dependence on the x-axis offset parameter $a$

As has been noted elsewhere[20], LogAUC measurements of different ROC curves
are comparable only if they use the same cutoff parameter $a$. However, it is
extremely important to point out that if $a$ is not set carefully for a given
dataset, then there are two problems that can arise:

(1) $a$ is too high, resulting in some fraction of the lowest ranked molecules
being simply discarded.

(2) $a$ is too low, resulting in the weight of the first term in the weighted
sum becoming arbitrarily large, thereby causing the first term to dominate the
score at the expense of the later terms.

In both cases, some information in the dataset goes unused towards calculating
LogAUC. In the former case, it is the information embodied in some number of
the lowest ranked Booleans, and in the latter case, it is the information
embodied in the higher ranked Booleans. In order to make use of all the
information available in a given dataset, $a$ must be set in a principled way
that avoids both of these problems.

So, by using the same $a$ for all datasets, one still runs the risk of losing
the ability to accurately compare different ROC curves against each other.
Thus, in order to ensure that we may accurately compare different ROC curves,
including those of different datasets, we need to set $a$ depending on the
dataset in a way that ensures that $a$ is neither too large nor too small for
that dataset.

To solve this problem, we propose setting the value of $a$ such that weight of
the first term in the weighted sum is always 1. This is achieved by setting
$a$ dynamically based on the number of decoys: $a=\frac{1}{en}$.

Suppose $f$ is a step function with values $y_{i}$ each on the interval
$[\frac{i}{n},\frac{i+1}{n}]$. Note that an ROC curve satisfies this
definition.

$f(x)=\begin{cases}y_{1},&\text{if $k_{0}<x\leq k_{1}$}\\\ y_{2},&\text{if
$k_{1}<x\leq k_{2}$}\\\ &\vdots\\\ y_{n},&\text{if $k_{n-1}<x\leq k_{n}$}\\\
\end{cases}.$

Then we can compute the integral (Equation 1) explicitly as a weighted sum of
the values $y_{i}$ of the step function:

$\displaystyle A(f;a)$
$\displaystyle=y_{1}\int_{a}^{\frac{1}{n}}\frac{dx}{x}+\sum_{i=2}^{n}y_{i}\int_{\frac{i-1}{n}}^{\frac{i}{n}}\frac{dx}{x}.$
(4)

Since an ROC curve is a step function, we use this weighted sum to calculate
$A(\operatorname{ROC};a)$.

Substituting our proposed value for $a$ into (Equation 4), we get:

$\displaystyle A(f;a=\frac{1}{en})$
$\displaystyle=y_{1}\int_{\frac{1}{en}}^{\frac{1}{n}}\frac{dx}{x}+\sum_{i=2}^{n}y_{i}\int_{\frac{i-1}{n}}^{\frac{i}{n}}\frac{dx}{x}$
$\displaystyle=y_{1}+\sum_{i=2}^{n}y_{i}\ln\left(\frac{i}{i-1}\right).$

Let $n=1000$. Then the terms in the sum look like

$y_{1}+0.69y_{2}+0.41y_{3}+0.29y_{4}+\dots+0.001y_{n}$

Exactly as we want: the first term contributes the most, and the respective
contributions of following terms diminish exponentially.

The last formula has obvious meaning: the first term has a weight of 1 by the
choice of $a=\frac{1}{en}$ and all succeeding terms have exponentially
diminishing weights, with the final term’s weight $\ln(\frac{n}{n-1})$ coming
arbitrarily close to zero with high enough $n$. Hence, the values $y_{i}$
matter less and less as $i$ increases. Moreover, setting $a$ in this manner
forces all the weights to be invariant to $a$, which would not be the case
otherwise due to the fluctuation of the first weight. Since the weight of the
first term is always 1, we never have to worry about $a$ being too high or too
low.

##  4 Enrichment Score

Taking into account the discussion above, we propose an alternative metric
called Enrichment Score given by

$\leavevmode\resizebox{368.57964pt}{}{$\begin{aligned}
\operatorname{EnrScore}(f)=\frac{A(f;a=\frac{1}{en})-A(x;a=\frac{1}{en})}{A(1;a=\frac{1}{en})-A(x;a=\frac{1}{en})}\end{aligned}$}.$
(5)

The enrichment score satisfies

  1. 1.

$\operatorname{EnrScore}(f=\operatorname{rand})=0$,

  2. 2.

$\operatorname{EnrScore}(f=1)=1$,

  3. 3.

$\operatorname{EnrScore}(f=0)<0$,

  4. 4.

The weight of the first term contributes $1$, and the weights of subsequent
terms diminish exponentially.

Let us elaborate on why the enrichment score provides an intuitive
understanding of the enrichment capacity of a given docking model.

First, the enrichment score of the random classifier $f=x$ is always zero. It
is intuitively appealing, as a random classifier merely shuffles data and
should therefore be expected to provide zero enrichment.

If the enrichment score of a classifier is positive, then that means that the
classifier performs better than the random classifier. If the score is
negative, it performs worse than the random classifier. The higher the score,
the better the classifier.

The best classifier $f=1$ gives the maximum possible score, $1$. The value of
the worst classifier $f=0$ varies depending on $a$, but it is always negative.

While the best and random classifier’s enrichment scores are both invariant,
it should be noted that the minimum value (i.e., the worst classifier’s
enrichment score) depends on $a$. However, it should also be the case that all
models with negative enrichment scores are equally useless as tools for
enrichment, so what the minimum score is should not matter in practice.

Let us compare the interpretability of the enrichment score to that of LogAUC
with arbitrary parameter $a$. If $\operatorname{LogAUC}(f;a)=0.12$, does it
perform better or worse than the random classifier? One cannot tell. To get
the answer, one needs to compute the LogAUC of the random classifier (which
depends on $a$) and compare values, which is neither intuitive nor convenient.

We remark that the enrichment score is closely related to the LogAUC. To wit,
the enrichment score is a normalization of LogAUC involving a shifting
(subtraction of random classifier LogAUC) and re-scaling (to enforce maximum
enrichment score of 1). The relation is expressed by the following formula:

$\begin{aligned}
\operatorname{EnrScore}(f)=\frac{\operatorname{LogAUC}(f;a=\frac{1}{en})-\operatorname{LogAUC}(\operatorname{rand};a=\frac{1}{en})}{1-\operatorname{LogAUC}(\operatorname{rand};a=\frac{1}{en})}\end{aligned}$

Therefore, one can still acquire an intuition of enrichment score by looking
at a linear-log ROC plot. We emphasize that the enrichment score is not the
area itself but the normalized LogAUC. (On that note, it should be remembered
that LogAUC itself is also not simply the log-scaled area under the curve, as
shown in (Equation 2).)

As one is more interested in a positive score (better than random classifier)
than a negative score (worse than random classifier), it makes more sense to
re-scale the score to induce the property that the maximum score
(corresponding to the performance of the best classifier) is 1. And the
careful choice of the offset parameter $a$ ensures that you don’t throw away
all of your data or inadequately weight the first term.

##  5 EnrichmentScore vs. unbalanced metrics: some examples

In this section, we show how a poor choice of the offset parameter $a$ leads
to one of the two scenarios previously described and demonstrate that using
enrichment score avoids these problems.

For the first case, let us examine what happens when $a$ is too large.
Consider an arbitrary ROC curve derived from a dataset with number of decoys
$n=10000$. Taking $a=1.0\times 10^{-3}$, which is the standard value we find
in the literature, we see clearly that the behavior of the ROC curve over the
range $[0,0.001)$ does not inform the enrichment score at all. This range
accounts for the first 10 inter-decoy intervals: $[0.0001,0.0002)$,
$[0.0002,0.0003)$, …, $[0.0009,0.001)$. So, using such a high value for $a$
means throwing away the information contained in those intervals. These
discarded inter-decoy intervals can have at most the value of the first kept
inter-decoy interval $y_{kept}$ (in this case, $y_{kept}=ROC(0.001)$).
Throwing these intervals away allows the resultant enrichment score to become
higher than it would be if these intervals were included in its calculation.
This should be sufficient to demonstrate why having too high a value for $a$
is worth avoiding.

Next, let us examine what may happen when $a$ is too small for a given
dataset. We consider a dataset formed by docking against one of the receptors
in the DUDE-Z dataset[10], a benchmark for molecular docking programs
providing actives and property-matched decoys for 43 target receptors.
Specifically, we docked against the target receptor labeled TRYB1, which is
the protease known as tryptase beta-1. We used DOCK 3.7[12] for our docking
program.

We performed two retrospective docking jobs on TRYB1, one with number of
decoys $n=100$ and the other with $n=50$. In both cases, we use all 38 ligands
included for TRYB1 in DUDE-Z as actives. As for decoys, for the dataset
defined by $n=100$, we randomly selected 100 of the 1911 decoys included for
TRYB1. In order to select the 50 decoys needed for the dataset defined by
$n=50$, we (1) docked on the dataset defined by $n=100$, (2) sorted the decoys
by docking score, and (3) deleted every other decoy in the sorted list (i.e.,
we removed the 2nd best scored decoy, 4th best scored decoy, …, 100th best
scored decoy). This procedure is intended to make the ROC curve for the
dataset defined by $n=50$ remain approximately the same shape as the ROC curve
for the dataset defined by $n=100$.

Using the value $a=$1.0\text{\times}{10}^{-3}$$ (previously mentioned to be
the standard $a$ according to the literature), we get the following linear-log
ROC plots for $n=100$ (Figure 2) and $n=50$ (Figure 3), respectively.

Figure 2: ROC curve for TRYB1 using $n=100$ and
$a=$1.0\text{\times}{10}^{-3}$$. Figure 3: ROC curve for TRYB1 using $n=50$
and $a=$1.0\text{\times}{10}^{-3}$$.

Notice that enrichment scores of the two plots vary, with a difference of
0.036.

Now, let us plot the same data using $a=\frac{1}{en}$ in the following Figure
4 and Figure 5.

Figure 4: ROC curve for TRYB1 using $n=100$ and
$a=$3.679\text{\times}{10}^{-3}$$. Figure 5: ROC curve for TRYB1 using $n=50$
and $a=$7.358\text{\times}{10}^{-3}$$.

$a=$3.679\text{\times}{10}^{-3}$$ was used for $n=100$, and
$a=$7.358\text{\times}{10}^{-3}$$ was used for $n=50$. Notice that the
enrichment scores for both are nearly identical, with a difference of only
0.004, which is considerably lower than the difference of 0.036 induced by
using $a=$1.0\text{\times}{10}^{-3}$$. This demonstrates that setting
$a=\frac{1}{en}$ avoids the variability in enrichment score witnessed for
different $n$ when using the apparently standard
$a=$1.0\text{\times}{10}^{-3}$$. Moreover, notice how far both enrichment
scores produced using $a=$1.0\text{\times}{10}^{-3}$$ are from those using
$a=\frac{1}{en}$. For $n=50$, the difference is a staggering 0.095. At this
point, it should be noted that the difference in enrichment score using too
low $a$ is determined by the value of the first interdecoy interval
$y_{1}=ROC(0)$, which the enrichment score comes arbitrarily close to with
small enough $a$. Therefore, in the limit, having too low $a$ is tantamount to
throwing away all inter-decoy intervals except the first one. The following
figures (fig. 6, fig. 7, fig. 8) demonstrate this for the dataset defined by
$n=50$ using smaller and smaller $a$.

Figure 6: ROC curve for TRYB1 using $n=50$ and
$a=$1.0\text{\times}{10}^{-4}$$. Figure 7: ROC curve for TRYB1 using $n=50$
and $a=$1.0\text{\times}{10}^{-5}$$. Figure 8: ROC curve for TRYB1 using
$n=50$ and $a=$1.0\text{\times}{10}^{-10}$$.

We see now that the direction of error in enrichment score induced by a poor
choice of $a$ is determined by the direction of error in choice of $a$: too
high $a$ causes enrichment score to become higher than it should be, and too
low $a$ causes enrichment score to become lower than it should be. In this
way, the respective causes for the problems of using too high $a$ or too low
$a$ complement each other.

##  6 Data and Software Availability

The software DOCK 3.7 is available through two licenses: (1) an academic and
non-profit license, which is free for academic researchers and those engaging
in non-profit research, and (2) a commercial license, which requires paying a
license fee. More information is available at the DOCK 3.7 website:
https://dock.compbio.ucsf.edu/DOCK3.7/. The DUDE-Z dataset is publicly
available at https://dudez.docking.org/.

##  7 Funding

J.J.I. was funded by NIGMS (Grant No. GM133836).

## References

  * [1] Torres, P. H., Sodero, A. C., Jofily, P., & Silva-Jr, F. P. (2019). Key topics in molecular docking for drug design. International journal of molecular sciences, 20(18), 4574.
  * [2] Halperin, I., Ma, B., Wolfson, H., & Nussinov, R. (2002). Principles of docking: An overview of search algorithms and a guide to scoring functions. Proteins: Structure, Function, and Bioinformatics, 47(4), 409-443.
  * [3] Warren, G. L., Andrews, C. W., Capelli, A. M., Clarke, B., LaLonde, J., Lambert, M. H., … & Head, M. S. (2006). A critical assessment of docking programs and scoring functions. Journal of medicinal chemistry, 49(20), 5912-5931.
  * [4] Jain, A. N. (2006). Scoring functions for protein-ligand docking. Current Protein and Peptide Science, 7(5), 407-420.
  * [5] Li, J., Fu, A., & Zhang, L. (2019). An overview of scoring functions used for protein–ligand interactions in molecular docking. Interdisciplinary Sciences: Computational Life Sciences, 11(2), 320-328.
  * [6] Shen, C., Ding, J., Wang, Z., Cao, D., Ding, X., & Hou, T. (2020). From machine learning to deep learning: Advances in scoring functions for protein–ligand docking. Wiley Interdisciplinary Reviews: Computational Molecular Science, 10(1), e1429.
  * [7] Li, H., Sze, K. H., Lu, G., & Ballester, P. J. (2021). Machine‐learning scoring functions for structure‐based virtual screening. Wiley Interdisciplinary Reviews: Computational Molecular Science, 11(1), e1478.
  * [8] Irwin, J. J. (2008). Community benchmarks for virtual screening. Journal of computer-aided molecular design, 22(3), 193-199.
  * [9] Mysinger, M. M., Carchia, M., Irwin, J. J., & Shoichet, B. K. (2012). Directory of useful decoys, enhanced (DUD-E): better ligands and decoys for better benchmarking. Journal of medicinal chemistry, 55(14), 6582-6594.
  * [10] Stein, R. M., Yang, Y., Balius, T. E., O’Meara, M. J., Lyu, J., Young, J., … & Irwin, J. J. (2021). Property-unmatched decoys in docking benchmarks. Journal of chemical information and modeling, 61(2), 699-714.
  * [11] Imrie, F., Bradley, A. R., & Deane, C. M. (2021). Generating property-matched decoy molecules using deep learning. Bioinformatics, 37(15), 2134-2141.
  * [12] Coleman, R. G., Carchia, M., Sterling, T., Irwin, J. J., & Shoichet, B. K. (2013). Ligand pose and orientational sampling in molecular docking. PloS one, 8(10), e75992.
  * [13] Mysinger, M. M., & Shoichet, B. K. (2010). Rapid context-dependent ligand desolvation in molecular docking. Journal of chemical information and modeling, 50(9), 1561-1573.
  * [14] Bender, A., & Glen, R. C. (2005). A discussion of measures of enrichment in virtual screening: comparing the information content of descriptors with increasing levels of sophistication. Journal of chemical information and modeling, 45(5), 1369-1375.
  * [15] Chen, H., Lyne, P. D., Giordanetto, F., Lovell, T., & Li, J. (2006). On evaluating molecular-docking methods for pose prediction and enrichment factors. Journal of chemical information and modeling, 46(1), 401-415.
  * [16] Truchon, J. F., & Bayly, C. I. (2007). Evaluating virtual screening methods: good and bad metrics for the “early recognition” problem. Journal of chemical information and modeling, 47(2), 488-508.
  * [17] Clark, R. D., & Webster-Clark, D. J. (2008). Managing bias in ROC curves. Journal of computer-aided molecular design, 22(3), 141-146.
  * [18] Sheridan, R. P., Singh, S. B., Fluder, E. M., & Kearsley, S. K. (2001). Protocols for bridging the peptide to nonpeptide gap in topological similarity searches. Journal of chemical information and computer sciences, 41(5), 1395-1406.
  * [19] Bender, B. J., Gahbauer, S., Luttens, A., Lyu, J., Webb, C. M., Stein, R. M., … & Shoichet, B. K. (2021). A practical guide to large-scale docking. Nature protocols, 16(10), 4799-4832.
  * [20] Thomas, T., Chalmers, D. K., & Yuriev, E. (2016). Homology Modeling and Docking Evaluation of Human Muscarinic Acetylcholine Receptors. In Muscarinic Receptor: From Structure to Animal Models (pp. 15-35). Humana Press, New York, NY.
  * [21] Ye, W. L., Shen, C., Xiong, G. L., Ding, J. J., Lu, A. P., Hou, T. J., & Cao, D. S. (2020). Improving docking-based virtual screening ability by integrating multiple energy auxiliary terms from molecular docking scoring. Journal of Chemical Information and Modeling, 60(9), 4216-4230.
  * [22] Thoma, M. (2018). File:Roc-draft-xkcd-style.svg. Wikimedia. Retrieved October 18, 2022, from https://commons.wikimedia.org/wiki/
File:Roc-draft-xkcd-style.svg#filelinks.


# Semiparametric Shape-restricted Estimators
for Nonparametric Regression

Kenta Takatsu  Department of Statistics and Data Science, Carnegie Mellon
University  Tianyu Zhang  Department of Statistics and Data Science, Carnegie
Mellon University  Arun Kumar Kuchibhotla  Department of Statistics and Data
Science, Carnegie Mellon University

###### Abstract

Estimating the conditional mean function that relates predictive covariates to
a response variable of interest is a fundamental task in statistics. In this
paper, we propose some general nonparametric regression approaches that are
widely applicable under very mild conditions. The method decomposes a function
with a Lipschitz continuous $k$-th derivative into a sum of a $(k-1)$-monotone
function and a parametric component. We implement well-established shape-
restricted estimation procedures (such as isotonic regression) to handle the
“nonparametric” components of the true regression function and combine them
with a simple sample-splitting procedure to estimate the parametric
components. The resulting estimators inherit several favorable properties from
the shape-restricted regression estimators. Notably, it is (practically)
tuning parameter-free, converges at the minimax rate, and exhibits a locally
adaptive rate when the true regression function is “simple”. Finally, a series
of numerical studies are presented, confirming these theoretical properties.

##  1 Introduction

This article considers the nonparametric regression problem of estimating the
conditional mean function based on observed covariates and response variables,
without assuming the truth taking any specific parametric forms. Specifically,
we consider $n$ independent and identically distributed (IID) observations
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}\sim P_{0}(X,Y)$, where the predictive covariate
variable vector $X_{i}$ takes values in some measurable space $\Omega$ (often
a subset of $\mathbb{R}^{d}$) and the response $Y_{i}$ is a real number. We
are interested in estimating the unknown conditional mean function
$f_{0}(X):=\operatorname{\mathbb{E}}[Y|X]$ that minimizes the mean-squared
prediction error of the response among all the square-integrable functions of
$X$. Defining the “error variables” as $\xi_{i}:=Y_{i}-f_{0}(X_{i})$, we can
write the relationship between the covariates and the response as

$\displaystyle Y_{i}=f_{0}(X_{i})+\xi_{i}\quad\text{for}\quad i=1,\dots,n.$
(1)

Many existing works assume $\xi_{i}$’s to be independent of $X_{i}$’s, which
is more justified if $\xi_{i}$ represents “external noise” contaminating the
samples. We note that no such independence assumptions are made in this paper.
The only property of the variables $\xi_{i}$’s implied directly by the
relationship (1) is that $\mathbb{E}[\xi_{i}|X_{i}]=0$.

Nonparametric regression methods usually assume the truth $f_{0}$ belong to
some function space. These function spaces often impose constraints on the
smoothness properties of the contained functions. One common class of
nonparametric functions considered in the literature is the Hölder class.
Formally, let $k=(k_{1},k_{2},\dots,k_{d})$ be a $d$-dimensional index set
where each $k_{i}$ is a non-negative integer and $|k|=\sum_{i=1}^{d}k_{i}$.
For each $f:\Omega\mapsto\mathbb{R}$ where
$x=(x_{1},x_{2},\dots,x_{d})\in\Omega\subseteq\mathbb{R}^{d}$, differentiable
up to the order $k\geq 1$, we define an operator $D^{k}$ as

$D^{k}f=\frac{\partial^{|k|}f(x)}{\partial^{k_{1}}x_{1}\dots\partial^{k_{d}}x_{d}},$

and $D^{0}f=f$. For $\beta,L>0$, the Hölder class $\Sigma(\beta,L)$ on
$\Omega$ consists of functions that satisfy the following condition:

$\displaystyle\Sigma(\beta,L)$
$\displaystyle:=\left\\{f:\Omega\mapsto\mathbb{R}\,|\,\|D^{m}f\|_{\infty}\leq
L\quad\text{for all}\quad m\in\mathbb{Z}_{+}^{d}\text{ with
}|m|=0,1,\dots\lfloor\beta\rfloor\quad\text{and}\right.$
$\displaystyle\qquad\left.|D^{k}f(y)-D^{k}f(x)|\leq L\|x-y\|^{\beta-|k|}\text{
with }|k|=\lfloor\beta\rfloor\quad\text{for all}\quad x,y\in\Omega\right\\}.$

When $\beta=1$, the Hölder class coincides with the class of $L$-Lipschitz
functions.

The performance of numerous classical methods have been analyzed for the
$\Sigma(\beta,L)$ class, including $k$-nearest neighbor, Nadaraya-Watson
estimator, local polynomial estimator, smoothing splines, series estimator,
tree-based methods, RKHS, random forest, and so on. These methods typically
involve a tuning parameter, such as a bandwidth for local smoothing methods or
the number of basis functions for the series estimator. The choice of the
parameters is crucial for both establishing the convergence rates and
achieving better practical performances. However, the optimal choice of such
tuning parameters is often unavailable in practice since it almost always
depends on the underlying distribution $P_{0}$. Therefore, data-adaptive model
selection procedures, such as cross-validation or Lepski’s method (Lepskii,,
1991), become necessary for estimating the optimal tuning parameters. While
these procedures can automate the model selection and achieve provably optimal
convergence rates, they are not entirely parameter free. For instance, both
methods require the candidate sets where the optimal parameter is contained,
which is often difficult to verify. Additionally, the performance of the
resulting estimator can still be sensitive to the specific choices made due to
the randomness of the selection process (Lei,, 2020).

In this paper, we propose a class of simple estimators based on shape-
restricted methods, which can be seen as tuning parameter free. The key
concept underlying our methodology can be summarized through the following
simple yet significant observations: For any $L$-Lipschitz function $f$, there
exists a corresponding function $g$ such that

$f(x)=g(x)-Lx\quad\text{where $g$ is a non-decreasing function}.$ (2)

Similarly, for any function $f\in\Sigma(2,L)$, there exists a function $g$
such that

$f(x)=g(x)-Lx^{2}\quad\text{where $g$ is a convex function}.$ (3)

We will discuss this decomposition with greater generality in Section 3. These
results can also be found in the optimization literature (Zlobec,, 2006).
Following (2) and (3), the original problem of nonparametric regression can be
decomposed into two subproblems: (1) estimating the shape-restricted
regression function $g$, and (2) estimating the parameter $L$. One may expect
the first problem to be more challenging as it involves the estimation of a
function under shape constraints, which is an infinite-dimensional object.
Therefore, the properties of the proposed class of estimators are expected to
resemble those of shape-restricted regression estimators. This is indeed the
case as we demonstrate in this article.

Given the decompositions listed in (2) and (3), the proposed estimator
naturally builds on existing nonparametric methods for estimating shape-
restricted regression functions. Shape-restricted methods differ significantly
from other nonparametric methods in several ways. Firstly, unlike popular
nonparametric estimators including kernel smoothing, local polynomials, or
orthogonal series estimators, shape-restricted estimators are often tuning
parameter free, eliminating the need for model selection procedures. Secondly,
despite their nonparametric nature, shape-restricted estimators can often be
computed efficiently. For instance, a least squares estimator over the
monotone function class can be solved in a near linear time. In addition to
their computational advantages, shape-restricted estimators often exhibit the
property of adaptive risk bounds. This implies that the estimator converges at
a faster rate than the worst-case, that is, the minimax rate over general
monotone functions, when the true function associated is “simple”. In the case
of estimating non-decreasing $f_{0}$, for instance, the least squares
estimator does not only achieve the minimax rate of $n^{-2/3}$ in terms of the
squared risk, but also achieves the parametric $n^{-1}$ rate (ignoring a
logarithmic term) if $f_{0}$ is non-decreasing piecewise constant. Remarkably,
the estimator achieves this faster convergence rate automatically without any
knowledge of the underlying function (and hence it is called adaptive).
Finally, shape-restricted regression can be extended to multivariate cases,
for instance, by assuming that the underlying function is additive in each
coordinate. When this assumption holds, the proposed estimator is no longer
confined to univariate problems. By directly building upon shape-restricted
regression estimators, our nonparametric estimators for the class
$\Sigma(\beta,L)$ also inherit the aforementioned properties. To summarize,
our proposed nonparametric estimators possess the following properties:

  1. 1.

Nonparametric consistency: The proposed estimator does not rely on parametric
assumptions and consistently estimates a regression function that belongs to a
nonparametric class of functions.

  2. 2.

Optimal and adaptive convergence rate: The proposed estimator demonstrates the
minimax optimal rate of convergence for the studied function classes. The
estimator is also adaptive in the sense that its risk automatically converges
at a faster rate if $f_{0}$ is “simple”.

  3. 3.

No tuning parameters: The proposed procedures essentially have no tuning
parameters, eliminating the need for the cumbersome model selection process.

  4. 4.

Efficient computation time: The proposed estimators can often be constructed
with an almost linear time expense.

  5. 5.

Support for multivariate covariates: The general estimation procedure can be
extended to multivariate covariates as long as a corresponding multivariate
shape-restricted estimator exists. When the additivity of a regression
function is assumed, the method can be extended with great generality.

##  2 Literature Review

The proposed method in this article is closely related to the extensive body
of literature on nonparametric shape-restricted regression estimation. The
problem is particularly well-studied in the univariate covariate setting,
where popular shape constraints include monotonicity and convexity. A common
approach is calculating the least squares estimator (LSE) using the observed
samples over the function class of interest (Guntuboyina and Sen,, 2018). When
the monotonicity constraint is considered, the LSE is commonly known as the
isotonic regression, which can be efficiently computed using the Pool Adjacent
Violators Algorithm (PAVA) (Ayer et al.,, 1955). The univariate regression
estimator under the convexity has also been studied in the literature
(Groeneboom et al.,, 2001; Guntuboyina and Sen,, 2015; Ghosal and Sen,, 2017).
Recently, there has been a surge in theoretical analysis associated with the
convergence rate of shape-restricted LSEs. The literature often emphasizes the
remarkable adaptive property of LSEs, where the estimator demonstrates a
faster convergence rate depending on the local structure around the true
regression function (Zhang,, 2002; Chatterjee et al.,, 2015; Guntuboyina and
Sen,, 2018; Bellec,, 2018). These results are typically derived under strong
assumptions regarding the distribution of $\xi$, such as sub-Gaussianity;
Zhang, (2002) is an exception as the errors here are only assumed to have a
finite second moment. Recent works have focused on relaxing such assumptions
and explored the behavior of LSEs in the presence of heavy-tailed errors (Han
and Wellner,, 2018, 2019; Kuchibhotla and Patra,, 2022). Multivariate
applications of shape-restricted methods have been investigated particularly
in recent years although a comprehensive understanding of theoretical
behaviors is still under development (Guntuboyina and Sen,, 2018). Han et al.,
(2019) studied the LSE estimator for the multivariate isotonic regression
while Han and Zhang, (2020) proved that an alternative method must be
considered to achieve minimax-optimal adaptivity for all dimensions. The
multivariate convex estimators have been studied in Seregin and Wellner,
(2010) and Seijo and Sen, (2011). Finally, additivity is a common structural
assumption in regression analysis. In this context, studying shape constraints
in conjunction with additivity helps maintain the theoretical properties of
shape-restricted methods. Some recent developments can be found in the works
of Mammen and Yu, (2007); Han and Wellner, (2018) and Chen and Samworth,
(2016).

As outlined earlier, the proposed estimator leverages a decomposition that
separates a nonparametric function into shape-restricted and parametric
components. The concept of decomposing nonparametric regression into
parametric components has been previously explored in the field of
semiparametric regression. In particular, the literature has investigated two-
stage estimation procedures that aim to improve the initial parametric
estimator through nonparametric methods. This approach has been adopted to the
context of density estimation (Hjort and Glad,, 1995; Hjort and Jones,, 1996),
regression estimation (Eguchi et al.,, 2003), conditional distribution
functions (Veraverbeke et al.,, 2014), and additive models (Lin et al.,,
2009). These methods often exhibit a faster rate of convergence when the
initial parametric estimator is correctly specified. However, this differs
from the notion of local adaptivity in the shape-restricted LSE literature. In
this case, adaptivity is not automatic and often relies on prior knowledge of
the true data-generating distribution to attain a faster convergence rate.

#### Notation.

For any positive integer $n\geq 1$, we denote by $[n]$ the index set
$\\{1,2,\dots,n\\}$. For any $j\in[d]$ with $d\geq 1$, we denote by $e_{j}$
the $d$-dimensional vector of zero’s with one at the $j$th position. For a
univariate function $f$, $f^{(k)}$ with a positive integer $k$ denotes the
$k$th derivative of $f$.

##  3 Decompositions of Nonparametric Regression Functions

As briefly mentioned in the introduction, the key concept behind the proposed
method is the general decomposition of a regression function within the Hölder
class into its shape-restricted and parametric component. We first discuss
this decomposition in detail.

A real-valued univariate function $g:\mathbb{R}\mapsto\mathbb{R}$ is
$k$-monotone if its $(k-1)$-th derivative is non-decreasing (Chatterjee et
al.,, 2015). Common examples are monotone $(k=1)$ and convex $(k=2)$
functions. Note that this definition differs from the one common in the
context of nonparametric density estimation (Balabdaoui and Wellner,, 2005;
Gao and Wellner,, 2009). We first provide the following decomposition for
univariate functions:

######  Proposition 1.

For any $k$-times differentiable function $f:\mathbb{R}\mapsto\mathbb{R}$ with
an $L$-Lipschitz $k$-th derivative, and any $\alpha\geq L$, there exists a
$(k+1)$-monotone function $g_{\alpha}:\mathbb{R}\mapsto\mathbb{R}$ such that
$f(x)=g_{\alpha}(x)-\\{\alpha/(k+1)!\\}x^{k+1}$.

###### Proof.

Let $f^{(k)}$ denote the $k$-th derivative of $f$. By the assumption that $f$
is $k$-times differentiable, it follows that $f(x)+\\{\alpha/(k+1)!\\}x^{k+1}$
is also $k$ times differentiable with the $k$-th derivative $f^{(k)}(x)+\alpha
x$. It thus remains to show $g_{\alpha}^{(k)}(x):=f^{(k)}(x)+\alpha x$ is non-
decreasing, which implies that $g_{\alpha}$ is $(k+1)$-monotone. This holds
because for any $y\geq x$,

$\displaystyle g^{(k)}_{\alpha}(y)-g^{(k)}_{\alpha}(x)$
$\displaystyle=\left(f^{(k)}(y)+\alpha y\right)-\left(f^{(k)}(x)+\alpha
x\right)\geq-L|y-x|+\alpha(y-x)\geq 0,$

which follows by the Lipschitz continuity of $f^{(k)}$ and $\alpha\geq L$.
This concludes the claim. ∎

While Proposition 1 is presented for a univariate function, analogous results
also hold for several multivariate functions, as demonstrated by the following
examples. To begin, we define coordinate-wise $k$-monotone functions. Given
the index set $k:=(k_{1},\dots,k_{d})$ of non-negative integers and
$\Omega\subseteq\mathbb{R}^{d}$, a function $g:\Omega\mapsto\mathbb{R}$ is
coordinate-wise $k$-monotone if for each $j\in[d]$ and for all $x\in\Omega$,
the $(k_{j}-1)$-th derivative of the univariate mapping $t\mapsto g(x+te_{j})$
is non-decreasing in $t\in\mathbb{R}$. We now extend Proposition 1 to the
context of multivariate covariates in the following result:

######  Proposition 2 (Coordinate-wise Lipschitz derivatives).

Suppose we have a function $f:\Omega\mapsto\mathbb{R}$ where
$\Omega\subseteq\mathbb{R}^{d}$. For each $j\in[d]$,
$x=(x_{1},\dots,x_{d})\in\Omega$, and an index set $k=(k_{1},\dots,k_{d})$ of
non-negative integers, we define a univariate mapping

$t\mapsto f_{j,x}(t)=\frac{\partial^{k_{j}}}{\partial
t^{k_{j}}}f(x+te_{j}),\quad\mbox{for}\quad t\in\mathbb{R}.$ (4)

For each $j\in[d]$, suppose $f_{j,x}$ is Lipschitz continuous with the
Lipschitz constant $L_{j}$ such that

$|f_{j,x}(t)-f_{j,x}(0)|\leq L_{j}|t|\quad\text{for all}\quad
x\in\Omega\quad\text{and}\quad t\in\mathbb{R}.$ (5)

Then for any sequence $(\alpha_{1},\ldots,\alpha_{d})$ such that
$\alpha_{j}\geq L_{j}$, there exists a coordinate-wise $(k+1)$-monotone
function $g_{\alpha}$ such that

$\displaystyle
f(x)=g_{\alpha}(x)-\sum_{i=1}^{d}\frac{\alpha_{i}}{(k_{i}+1)!}x_{i}^{k_{i}+1}.$

###### Proof.

For each $j\in[d]$, it follows that

$\displaystyle\frac{\partial^{k_{j}}}{\partial t^{k_{j}}}g_{\alpha}(x+te_{j})$
$\displaystyle=\frac{\partial^{k_{j}}}{\partial
t^{k_{j}}}\left(f(x+te_{j})+\sum_{i=1}^{d}\frac{\alpha_{i}}{(k+1)!}(x+te_{j})^{k+1}\right)$
$\displaystyle=\frac{\partial^{k_{j}}}{\partial
t^{k_{j}}}\left(f(x+te_{j})+\sum_{i\neq
j}\frac{\alpha_{i}}{(k_{i}+1)!}x_{i}^{k_{i}+1}+\frac{\alpha_{j}}{(k_{j}+1)!}(x_{j}+t)^{k_{j}+1}\right)$
$\displaystyle=f_{j,x}(t)+\alpha_{j}(x_{j}+t).$

It remains to show that the univariate function
$t\mapsto\frac{\partial^{k_{j}}}{\partial t^{k_{j}}}g_{\alpha}(x+te_{j})$ is
non-decreasing in $t\in\mathbb{R}$. From the derivation above,
$\frac{\partial^{k_{j}}}{\partial t^{k_{j}}}g_{\alpha}(x+te_{j})$ equals
$f_{j,x}(0)+\alpha_{j}x_{j}$ when $t=0$. Thus for $t>0$, we have

$\displaystyle
f_{j,x}(t)+\alpha_{j}(x_{j}+t)-(f_{j,x}(0)+\alpha_{j}x_{j})\geq-
L_{j}|t|+\alpha_{j}t\geq 0.$

Hence, this concludes that $g_{\alpha}$ is coordinate-wise $(k+1)$-monotone. ∎

As a concrete application of Proposition 2, we obtain the following result:

######  Remark 1 (Coordinate-wise Lipschitz functions).

Suppose $f:\Omega\mapsto\mathbb{R}$ for $\Omega\subseteq\mathbb{R}^{d}$ is a
coordinate-wise Lipschitz function with a Lipschitz constant
$L:=(L_{1},\dots,L_{d})$ where

$|f(x)-f(x+e_{j}h)|\leq|h|L_{j}\quad\text{for any}\quad h\in\mathbb{R}.$

There exits a coordinate-wise monotone function
$g_{\alpha}:\Omega\mapsto\mathbb{R}$ such that
$f(x)=g_{\alpha}(x)-\alpha^{\top}x$ for any
$\alpha=(\alpha_{1},\dots,\alpha_{d})$ satisfying $\alpha_{i}\geq L_{i}$ for
all $i=1,\dots,d$.

The above result follows as the direct application of Proposition 2 where
$k_{1}=\dots=k_{d}=0$.

We now introduce additional multivariate extensions of Proposition 1. However,
readers who are primarily interested in understanding the proposed methods can
skip to Section 4 as the remaining part of this section may not be crucial for
this purpose. First, a similar result to Proposition 2 holds when the gradient
of $f:\Omega\mapsto\mathbb{R}$ is Lipschitz continuous in $L_{2}$-norm.

######  Proposition 3 (Multivariate function with a Lipschitz continuous
gradient).

Suppose $f:\Omega\mapsto\mathbb{R}$ for $\Omega\subseteq\mathbb{R}^{d}$ has a
Lipschitz continuous gradient in $L_{2}$-norm such that

$\|D^{1}f(y)-D^{1}f(x)\|_{2}\leq L\|x-y\|_{2}$

for some $L>0$. Then there exits a multivariate convex function
$g_{\alpha}:\Omega\mapsto\mathbb{R}$ for $\Omega\subseteq\mathbb{R}^{d}$ such
that $f(x)=g_{\alpha}(x)-\alpha x^{\top}x$ for any $\alpha$ satisfying
$\alpha\geq L$.

###### Proof.

A multivariate function $g:\Omega\mapsto\mathbb{R}$ is convex if and only if

$(y-x)^{\top}\left(D^{1}g(y)-D^{1}g(x)\right)\geq 0$

for any $x,y\in\Omega$. By definition of $g_{\alpha}$ and $f$, it follows that

$\displaystyle(y-x)^{\top}\left(D^{1}g_{\alpha}(y)-D^{1}g_{\alpha}(x)\right)$
$\displaystyle=(y-x)^{\top}\
\left(D^{1}f(y)-D^{1}f(x)\right)+\alpha\|x-y\|_{2}^{2}$
$\displaystyle\geq-\|y-x\|_{2}\left\|D^{1}f(y)-D^{1}f(x)\right\|_{2}+\alpha\|x-y\|_{2}^{2}$
$\displaystyle\geq-L\|y-x\|_{2}^{2}+\alpha\|x-y\|_{2}^{2}\geq 0$

which follows by Hölder’s inequality and the assumed Lipschitz continuity.
Thus, the function $g_{\alpha}$ is convex. ∎

We note that the parameter $L$ is a $d$-dimensional vector in Proposition 2
while it is a scalar constant in Proposition 3. The final example involves the
generalized additive index model, which was initially studied by Chen and
Samworth, (2016). In this model, the multivariate function admits the additive
decomposition:
$f(x_{1},x_{2},\dots,x_{d})=\sum_{i=1}^{m}f_{i}(\beta_{i}^{\top}x)$ for some
$\beta_{i}\in\mathbb{R}^{d}$ and $1\leq i\leq m$. Since each $f_{i}$ is a
univariate function, we can extend Proposition 1 to the context of the
additive model in the following proposition:

######  Proposition 4 (Shape-restricted multi-index model).

Suppose $f:\Omega\mapsto\mathbb{R}$ where $\Omega\subseteq\mathbb{R}^{d}$ and
suppose there exists $\beta_{i}\in\mathbb{R}^{d}$ ($1\leq i\leq m$) such that
$f(x_{1},x_{2},\dots,x_{d})=\sum_{i=1}^{m}f_{i}(\beta_{i}^{\top}x)$ for all
$x\in\Omega$. If $f_{i}$ is $k_{i}$-times differentiable function with an
$L_{i}$-Lipschitz $k_{i}$-th derivative for all $1\leq i\leq m$, then for any
$\alpha_{i}\geq L_{i}$, there exists $(k_{i}+1)$-monotone funct

$\displaystyle
f(x_{1},\ldots,x_{d})=\sum_{i=1}^{m}g_{i,\alpha_{i}}(\beta_{i}^{\top}x)-\sum_{i=1}^{m}\\{\alpha_{i}/(k_{i}+1)!\\}(\beta_{i}^{\top}x)^{k_{i}+1}$

where $g_{i,\alpha_{i}}$ is a $(k_{i}+1)$-monotone function. In particular, if
$k_{i}=0$ for all $1\leq i\leq d$, then, for any $\alpha_{i}\geq L_{i}$, there
exists non-decreasing functions $g_{i,\alpha_{i}}(\cdot)$ such that

$f(x_{1},\ldots,x_{d})=\sum_{i=1}^{d}g_{i,\alpha_{i}}(\beta_{i}^{\top}x)-\ell^{\top}x,$

where $\ell:=\sum_{i=1}^{d}\alpha_{i}\beta_{i}$.

The “standard” additive model is a special case of Proposition 4 as follows:

######  Remark 2.

Taking $\beta_{i}=e_{i}$ for $1\leq i\leq m=d$ in Proposition 4 implies the
decomposition of the following additive model:

$\displaystyle\sum_{i=1}^{d}f_{i}(x_{i})=\sum_{i=1}^{d}g_{i,\alpha_{i}}(x_{i})-\sum_{i=1}^{d}\\{\alpha_{i}/(k_{i}+1)!\\}x_{i}^{k_{i}+1}$
(6)

where $g_{i,\alpha_{i}}$ is a $(k_{i}+1)$-monotone function.

We have demonstrated that functions in the nonparametric class can be
decomposed into the shape-restricted and parametric components. This
decomposition holds true for a range of applications, including $k$-times
differentiable functions for any integer $k$ or certain multivariate cases.
Based on this observation, we introduce a class of nonparametric regression
estimators that leverages shape-restricted regression estimators, as detailed
in the subsequent section.

######  Remark 3 (Application to the Sobolev space $W^{s,p}(\Omega)$).

Although all the results in this section are presented in terms of Lipschitz
continuous $k$-th derivatives, they can be extended to the Sobolev space using
the Sobolev embedding theorem (Adams and Fournier,, 2003; Kinnunen,, 2017).
Assume $\Omega$ is an open subset of $\mathbb{R}^{d}$. For a multi-index $k$,
the Sobolev space $W^{s,p}(\Omega)$ consists of functions $f$ in
$L^{p}(\Omega)$ such that for every multi-index $k$ with $|k|\leq s$ the weak
derivative $D^{k}f$ belongs to $L^{p}(\Omega)$. We thus have

$\displaystyle W^{s,p}(\Omega):=\left\\{f:\Omega\mapsto\mathbb{R}\,|\,f\in
L^{p}(\Omega)\quad\text{and}\quad D^{k}f\in L^{p}(\Omega)\quad\text{for
all}\quad|k|\leq s\right\\}.$

The Sobolev space does not require pointwise differentiability of the
functions it contains, but instead characterizes them based on the
integrability of their weak derivatives. The Sobolev embedding theorem
establishes that a function in the Sobolev space is a subset of a suitably
chosen Hölder space, which allows us to extend the results of Proposition 1
with the corresponding $k$. The theorem states

$W^{s,p}(\Omega)\subset\Sigma(s-d/p,L)\quad\text{when $ps>d$}.$

As a limiting case, we simply have $W^{1,\infty}(\mathbb{R})=\Sigma(1,L)$,
equivalently the space of $L$-Lipschitz functions.

##  4 Proposed Estimators

In this section, we present a general estimation procedure for nonparametric
regression functions that admit the decomposition described in Section 3.
While we focus on the estimation of $L$-Lipschitz functions for exposition,
the following method can be extended to the general case with minimal
modifications. To establish the underlying intuition, we begin by considering
the estimation of the univariate function. By applying Proposition 1 with
$k=0$, any $L$-Lipschitz function can be decomposed as
$f(x)=g_{\alpha}(x)-\alpha x$ for any $\alpha\geq L$ and $g_{\alpha}$ is a
non-decreasing function. If the true Lipschitz constant were known, we could
easily estimate $g_{L}$ by regressing $\\{y_{i}+Lx_{i}\\}_{i=1}^{n}$ on
$x_{1},\dots,x_{n}$ using an isotonic regression. However the exact value of
$L$ is rarely available in practice, and thus it needs to be estimated.

Although a simple LSE over the class of functions $g_{\alpha}(x)-\alpha x$ may
seem plausible, it is not immediately effective. This is due to the fact that
any finite data can be interpolated using a Lipschitz function by selecting an
adequately large constant $\alpha$. As a result, this class of functions can
always achieve zero squared error for any observed data, which is reminiscent
of overfitting. Hence, it is crucial to incorporate an estimation procedure
for $L$ beyond relying solely on simple least squares.

Algorithm 1 outlines a simple procedure based on sample-splitting. The
procedure can be summarized as follows: First, we split the data into two
independent sets. For each value of $\alpha$ in a pre-specified candidate set,
we construct the estimator $\widehat{g}_{\alpha}(x)$ based on the first set.
This can be any isotonic estimator including the LSE (Guntuboyina and Sen,,
2018), the monotone spline estimator (Meyer,, 2008), or the smoothed monotone
estimator (Mammen,, 1991; Dette et al.,, 2006), in the univariate case. In the
multivariate case, one can consider the block min-max estimator of Han and
Zhang, (2020) or the Bayesian monotone estimator (Wang and Ghosal,, 2023). The
value of $\widehat{\alpha}$ and thereby the final estimator
$\widehat{g}_{\widehat{\alpha}}(x)-\widehat{\alpha}^{\top}x$ is selected if it
minimizes the estimated risk on the second set.

Algorithm 1 Lipschitz Function Estimation using Monotone Estimators

Observation sequence
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}\subset\mathbb{R}^{d}\times\mathbb{R}$, set of
candidate vectors $\mathcal{L}\subset\mathbb{R}^{d}$.

An estimator of $f_{0}(x)=\mathbb{E}[Y|X=x]$ assuming $f_{0}$ is coordinate-
wise Lipschitz.

  1. 1.

Randomly split $[n]=\\{1,2,\ldots,n\\}$ into two disjoint subsets:
$\mathcal{I}_{1}$ and $\mathcal{I}_{2}$.

  2. 2.

For each $\alpha\in\mathcal{L}$, compute an isotonic regression estimator
$\widehat{g}_{\alpha}$ based on
$(X_{i}\,Y_{i}+\alpha^{\top}X_{i})_{i\in\mathcal{I}_{1}}$. Define the
resulting estimator for fixed $\alpha$ as

$\widehat{f}_{\alpha}(x)=\widehat{g}_{\alpha}(x)-\alpha^{\top}x.$

  3. 3.

Find the “optimal” vector $\widehat{\alpha}$ based on
$\\{(X_{i},Y_{i}):i\in\mathcal{I}_{2}\\}$ by:

$\widehat{\alpha}:=\operatorname*{arg\,min}_{\alpha\in\mathcal{L}}\sum_{i\in\mathcal{I}_{2}}(Y_{i}-\widehat{f}_{\alpha}(X_{i}))^{2}.$

  4. 4.

Return the final estimator
$\widehat{f}_{\widehat{\alpha}}(x):=\widehat{g}_{\widehat{\alpha}}(x)-\widehat{\alpha}^{\top}x$.

While the proposed algorithm can be interpreted as performing cross-validation
over the parameter $\alpha$, it also possesses unique and desirable
properties. Firstly, the core decomposition based on Proposition 1 holds for
any $\alpha\geq L$ without requiring $\alpha=L$. This means that a “good”
estimate of the regression function $f_{0}$ in terms of mean squared error can
be obtained without precisely identifying the true value of $L$. As we
demonstrate later, the estimator exhibits robust performance in terms of the
mean squared error when a sufficiently large $\alpha$ is chosen (see Figure
4). This stands in contrast to standard model-selection procedures based on
cross-validation. Second, the third step in Algorithm 1 can be implemented
using a numerical optimization program such as the Broyden-Fletcher-Goldfarb-
Shanno (BFGS) algorithm (Byrd and Nocedal,, 1989). Consequently, in practice,
the only restriction is that the candidate set $\mathcal{L}$ must contain
large enough values of $\alpha$ such that $\alpha\geq L$ holds for some
$\alpha\in\mathcal{L}$. We will also demonstrate that the proposed procedure
is robust against the randomness introduced by the data split, holding for a
fixed observation sequence (see Figure 5), which is a desirable yet uncommon
property for methods based on cross-validation.

The extension to functions with Lipschitz $k$-th derivatives is
straightforward by replacing $\alpha x$ with the appropriate polynomials as
provided by Proposition 1. In this case, $\widehat{g}_{\alpha}$ can be
obtained using a suitable estimator for $(k+1)$-monotone functions. For
example, when $k=1$, this corresponds to an estimator for convex functions.
The same approach can be applied to the multivariate examples presented in
Propositions 2 and 3 by adopting the coordinate-wise isotonic regression or
multivariate convex function and appropriate polynomial terms accordingly.

Next, in view of Proposition 4, we outline the procedure when the additivity
assumption is plausible. We focus on the estimator for the additive Lipschitz
functions. However, the same algorithm can be straightforwardly extended to
additive functions where each component has Lipschitz $k$-th derivatives. The
detailed procedure is outlined in Algorithm 2. While the basic concept behind
the estimation procedure remains the same, there are notable differences in
steps 2 and 3 compared to Algorithm 1. In this case, we perform a coordinate-
descent procedure to update $\widehat{g}^{(i)}$ for each $i=1,\dots d$
component until the empirical risk converges.

Algorithm 2 Additive Lipschitz Function Estimation using Monotone Estimators

Observation sequence
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}\subset\mathbb{R}^{d}\times\mathbb{R}$ where
$X_{i}=(X_{i1},\dots,X_{id})$, set of candidate vectors
$\mathcal{L}\subset\mathbb{R}^{d}$.

An estimator of $f_{0}(x)=\operatorname{\mathbb{E}}[Y|X=x]$ assuming that
$f_{0}$ is additive and each component is Lipschitz.

  1. 1.

Initialize estimators $\widehat{g}^{(i)}(x)$ for $i=1,2,\dots,d$. These can,
for instance, be constant functions.

  2. 2.

Randomly split $[n]=\\{1,2,\dots,n\\}$ into two disjoint subsets:
$\mathcal{I}_{1}$ and $\mathcal{I}_{2}$.

  3. 3.

For each $\ell\in\\{1,2,\dots,d\\}$ and $\alpha\in\mathcal{L}$, compute an
isotonic regression estimator $\widehat{g}^{(\ell)}_{\alpha}$ based on
$(X_{i\ell},\,Z_{i,\alpha}-\sum_{j\neq\ell}g^{(j)}(X_{ij}))_{i=1}^{n}$ where
$Z_{i,\alpha}:=Y_{i}+\alpha^{\top}X_{i}$. For instance, one may define the LSE
for each coordinate as

$\widehat{g}_{\alpha}^{(\ell)}:=\operatorname*{arg\,min}_{g\in\mathcal{G}}\sum_{i\in\mathcal{I}_{1}}\left\\{\left(Z_{i,\alpha}-\sum_{j\neq\ell}\widehat{g}^{(j)}(X_{ij})\right)-g(X_{i\ell})\right\\}^{2},$

iteratively obtaining an isotonic regression for each coordinate while keeping
estimators for the remaining coordinates fixed.

  4. 4.

Repeat step 3 until the empirical risk associated with $\mathcal{I}_{1}$ given
by

$\sum_{i\in\mathcal{I}_{1}}\left(Z_{i,\alpha}-\sum_{j=1}^{d}\widehat{g}_{\alpha}^{(j)}(X_{ij})\right)^{2}$

converges. Define the resulting estimator for fixed $\alpha$ as
$\widehat{f}_{\alpha}(x)=\sum_{j=1}^{d}\widehat{g}_{\alpha}^{(j)}(x)-\alpha^{\top}x$.

  5. 5.

Find the “optimal” vector $\widehat{\alpha}$ based on
$\\{(X_{i},Y_{i}):i\in\mathcal{I}_{2}\\}$ by:

$\widehat{\alpha}:=\operatorname*{arg\,min}_{\alpha\in\mathcal{L}}\sum_{i\in\mathcal{I}_{2}}(Y_{i}-\widehat{f}_{\alpha}(X_{i}))^{2}.$

  6. 6.

Return the final estimator
$\widehat{f}_{\widehat{\alpha}}(x):=\sum_{j=1}^{d}\widehat{g}_{\widehat{\alpha}}^{(j)}(x)-\widehat{\alpha}^{\top}x$.

##  5 Convergence Rates

The convergence rate of (shape-restricted) LSEs has been well-studied in the
literature. In this section, we present the convergence rates of the proposed
estimator when $\alpha$ is fixed (i.e., not random) and the shape-restricted
LSE $\widehat{g}_{\alpha}$ is used in Step 2 of Algorithm 1. Specifically, we
let $Z_{i,\alpha}:=Y_{i}+\alpha X_{i}$, and the shape-restricted least squares
estimator (Guntuboyina and Sen,, 2018) is given by

$\displaystyle\widehat{g}_{\alpha}=\operatorname*{arg\,min}_{g\in\mathcal{G}}\sum_{i=1}^{n}\left(Z_{i,\alpha}-g(X_{i})\right)^{2}$
(7)

where $\mathcal{G}$ is a collection of functions that satisfy constraints such
as monotonicity.

We now introduce the necessary notation and conditions. For the marginal
distribution $P_{X}$ of $X$ and $P_{X}$-square integrable function $g$, we
denote the squared $L_{2}(P_{X})$-norm as $\|g\|^{2}_{L_{2}(P_{X})}:=\int
g^{2}(x)\,dP_{X}$. Using this notation, the mean squared error of the
estimator integrated under the marginal distribution of $X$ is given by
$\|\widehat{f}_{n}-f_{0}\|^{2}_{L_{2}(P_{X})}=\int(\widehat{f}_{n}-f_{0})^{2}\,dP_{X}$,
which is the primary object of interest in this section. Furthermore, we
assume the univariate regression model (1) where $X_{1},\dots
X_{n}\in\mathcal{X}:=[0,1]$. We then introduce the following assumptions:

  1. (A1)

There exists a constant $C_{1}>0$ such that $\|f_{0}\|_{\infty}<C_{1}$ and
$\|\widehat{f}_{n}\|_{\infty}=O_{P}(1)$.

  2. (A2)

There exists a constant $C_{2}>0$ such that
$\operatorname{\mathbb{E}}[Y^{2}\mid X=x]\leq C_{2}$.

  1. (A2’)

There exists a constant $C_{3}>0$ such that
$\operatorname{\mathbb{E}}[\xi^{2}\mid X]\leq C_{3}$ $P_{X}$-almost everywhere
and $\operatorname{\mathbb{E}}[|\xi|^{3}\mid X]$ is uniformly bounded.

  1. (A3)

The Lebesgue density of $X$ satisfies $c_{4}<dP_{X}<C_{4}$ for some
$c_{4},C_{4}>0$.

  1. (A3’)

$P_{X}$ is the uniform distribution.

Broadly speaking, these assumptions can be grouped into three categories:
(A1), (A2)–(A2’), and (A3)–(A3’). We note that (A2’)$\implies$(A2) and (A3’)
$\implies$ (A3). (A1) is a common requirement that assumes the uniform
boundedness of $f_{0}$ and the boundedness in probability of
$\widehat{f}_{n}$. Assumptions (A2) and (A2’) impose the moments of $Y$ (and
$\xi$). Traditionally, more stringent structures are assumed for $\xi_{i}$’s,
such as sub-Gaussianity, but recent works have focused on analyzing scenarios
involving heavy-tailed $\xi_{i}$’s (Han and Wellner,, 2018; Kuchibhotla and
Patra,, 2022). Assumptions (A3) and (A3’) pertain to the structure of the
marginal distribution of $P_{X}$. Notably, the independence between
$\xi_{i}$’s and $X_{i}$’s is not required for the following result we state.

######  Theorem 1.

For any $L$-Lipschitz function $f_{0}$ and $\alpha\geq L$, assuming (A1), (A2)
and (A3),

$\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}=O_{P}(n^{-1/3}).$

Similarly for $f_{0}\in\Sigma(2,L)$ and $\alpha\geq L$, assuming (A1), (A2’)
and (A3’),

$\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}=O_{P}(n^{-2/5}\log n).$

###### Proof.

Under the regression model (1) with $L$-Lipschitz function $f_{0}$,
Proposition 1 with $k=0$ implies that

$\displaystyle Y_{i}=f_{0}(X_{i})+\xi_{i}=g_{\alpha}(X_{i})-\alpha
X_{i}+\xi_{i}$

for any $\alpha\geq L$ where $g_{\alpha}$ is a non-decreasing function. It
then follows that

$\displaystyle\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}$
$\displaystyle=\|(\widehat{g}_{\alpha}-\alpha x)-(g_{\alpha}-\alpha
x)\|_{L_{2}(P_{X})}=\|\widehat{g}_{\alpha}-g_{\alpha}\|_{L_{2}(P_{X})}.$

The last quantity has been investigated extensively in the literature on LSEs,
and similar results can be obtained under different assumptions. For instance,
Theorem 4.1 of Balabdaoui et al., (2019) (with $d=1$ and
$\alpha_{0}=\widehat{\alpha}=1$) states that under (A1), (A2) and (A3),

$\displaystyle\|\widehat{g}_{\alpha}-g_{\alpha}\|_{L_{2}(P_{X})}=O_{P}(n^{-1/3})$

which concludes the first claim.

For $f_{0}\in\Sigma(2,L)$, Proposition 1 with $k=1$ implies that

$\displaystyle Y_{i}=f_{0}(X_{i})+\xi_{i}=g_{\alpha}(X_{i})-\alpha
X_{i}^{2}+\xi_{i}$

where $g_{\alpha}$ is a convex function. Theorem 3.1 of Kuchibhotla and Patra,
(2022) (with $\alpha=1/2,s=2/3,\nu=1/3,\Phi=C\sqrt{\log n}$ and $A=(\log
n)^{1/4}$) as well as its discussion in Section 3.1 shows that assuming (A1),
(A2’) and (A3’)

$\|\widehat{g}_{\alpha}-g_{\alpha}\|_{L_{2}(P_{X})}=O_{P}(n^{-2/5}\log n).$
(8)

This proves the second claim. ∎

######  Remark 4 (On the logarithmic factor in (8)).

The logarithmic factor in the upper-bound (8) can be relaxed even for a heavy-
tailed $\xi_{i}$’s under the independence between $\xi_{i}$ and $X_{i}$. See,
for instance, Theorem 3 of Han and Wellner, (2019) for the detailed
conditions.

Next, we present the convergence rate of the estimator when it is expected to
be locally adaptive. The results are based on Han and Wellner, (2018), which
provides the adaptive rate of shape-restricted LSEs. The following assumptions
are introduced to the regression model (1):

  1. (A4)

$\xi_{1},\dots,\xi_{n}$ are IID and independent of $X_{1},\dots,X_{n}$.

  2. (A5)

$\xi$ has finite $L_{2,1}$-norm, defined as
$\|\xi\|_{2,1}:=\int_{0}^{\infty}P(|\xi|>t)^{1/2}\,dt<\infty$.

We conjecture that (A4) can be relaxed in view of the recent work by
Kuchibhotla and Patra, (2022). (A5) assumes finite $L_{2,1}$-moment of $\xi$,
which is only slightly stronger than assuming finite second moment of $\xi$
but not stronger than assuming $2+\delta$ moment for any $\delta>0$. We define
$\mathcal{M}_{m}$ as a collection of non-decreasing $m$-piecewise constant
functions and $\mathcal{C}_{m}$ as a collection of convex $m$-piecewise linear
functions. We now have the following result:

######  Theorem 2.

Suppose $f_{0}(x)=f_{m}(x)+Lx$ for all $x$, for some
$f_{m}\in\mathcal{M}_{m}$. Then, under (A1), (A3’)–(A5) and
$0<\alpha-L=O(\sqrt{\frac{m}{n}\log^{2}n})$,

$\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}=O_{P}\left(\sqrt{\frac{m}{n}\log^{2}n}\right).$

Similarly, if $f_{0}(x)=f_{m}(x)+Lx^{2}$ for some $f_{m}\in\mathcal{C}_{m}$,
then under (A1), (A3’)–(A5) and $0<\alpha-L=O(\sqrt{\frac{m}{n}\log^{2}n})$,

$\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}=O_{P}\left(\sqrt{\frac{m}{n}\log^{2}n}\right).$

###### Proof.

By the assumption, there exists $L$ and $f_{m}\in\mathcal{M}_{m}$ such that
$f_{0}(x)=f_{m}(x)-Lx$. It now follows by Proposition 1 that

$\displaystyle
P\left(\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}\geq\delta_{n}\right)$
$\displaystyle=P\left(\|(\widehat{g}-\alpha x)-(f_{m}-Lx+\alpha x-\alpha
x)\|_{L_{2}(P_{X})}\geq\delta_{n}\right)$ $\displaystyle\leq
P\left(\|\widehat{g}-f_{m}\|_{L_{2}(P_{X})}+\|(L-\alpha)x\|_{L_{2}(P_{X})}\geq\delta_{n}\right)$
$\displaystyle\leq
P\left(\|\widehat{g}-f_{m}\|_{L_{2}(P_{X})}\geq\delta_{n}/2\right)+P\left(|L-\alpha|\|x\|_{L_{2}(P_{X})}\geq\delta_{n}/2\right).$
(9)

Theorem 3 by Han and Wellner, (2018) states that under (A1) and (A3’)–(A5),

$\displaystyle P\left(\|\widehat{g}-f_{m}\|_{L_{2}(P_{X})}\geq
c\sqrt{\frac{m}{n}\log^{2}n}\right)=\varepsilon$ (10)

for any $\varepsilon\in(0,1)$ where the constant $c$ depends only on
$\varepsilon,\|\xi\|_{2,1},\|f_{m}\|_{\infty}$ and $\mathcal{G}$. Since
$\alpha-L=O(\sqrt{\frac{m}{n}\log^{2}n})$ and $\|x\|_{L_{2}(P_{X})}=1/9$, we
can select $\delta_{n}=\sqrt{\frac{m}{n}\log^{2}n}$ to make the right-hand
side of the above inequality arbitrarily small. This concludes the proof.

When $f_{0}$ has representation $f_{m}(x)-Lx^{2}$ for some
$f_{m}\in\mathcal{C}_{m}$, it follows analogously that

$\displaystyle
P\left(\|\widehat{f}_{n}-f_{0}\|_{L_{2}(P_{X})}\geq\delta_{n}\right)\leq
P\left(\|\widehat{g}-f_{m}\|_{L_{2}(P_{X})}\geq\delta_{n}/2\right)+P\left(|L-\alpha|\|x^{2}\|_{L_{2}(P_{X})}\geq\delta_{n}/2\right).$

The rest of the proof follows as an analogous application of Theorem 3 by Han
and Wellner, (2018). ∎

Contrary to Theorem 1 where the result holds for all $\alpha\geq L$, the
adaptive rate of Theorem 2 requires $\alpha-L=O(\sqrt{m\log^{2}n/n})$. This
indicates that the proposed method requires slightly more accurate estimate of
$L$ in order to obtain the adaptive rate.

######  Remark 5 (The rates under fixed design points).

In a nonparametric regression model, it is also common to consider a scenario
where the design for $X$ is fixed, meaning that $X_{1},\dots,X_{n}$ are
deterministic and are not random or IID. An illustrative example is the equi-
distant design, where $X_{i}$ is defined as $i/n$ for $i=1,\dots,n$. The upper
bound on the squared risk and the local adaptivity of the LSEs under the fixed
design have been shown by Zhang, (2002); Chatterjee et al., (2015); Bellec,
(2018); Guntuboyina and Sen, (2018) for sub-exponential or sub-Gaussian
$\xi$’s.

######  Remark 6 (The minimax optimal convergence rates for multivariate
$\mathcal{X}$).

The preliminary results in this section focus on the univariate (and convex
set) $\mathcal{X}\subset\mathbb{R}$. For the multivariate function class
corresponding to Proposition 2, it is known that the LSEs fail to adapt to
non-decreasing piecewise constant functions at the minimax rate in $L_{2}$ for
$d\geq 3$ (Han et al.,, 2019). Hence, an alternative approach such as the
block max-min estimator must be considered in order to recover minimax
optimality (Han and Zhang,, 2020). We conjecture that the minimax rate of the
block max-min estimator proved by Han and Zhang, (2020) can be extended to our
setting.

Finally, we discuss the convergence rate of shape-restricted additive models.
Remark 2 implies that an additive function can be decomposed as follows:

$\displaystyle\sum_{i=1}^{d}f_{i}(x_{i})$
$\displaystyle=\sum_{i=1}^{d}g_{i,\alpha_{i}}(x_{i})-\sum_{i=1}^{d}\\{\alpha_{i}/(k_{i}+1)!\\}x_{i}^{k_{i}+1}$
$\displaystyle=\mu^{*}+\sum_{i=1}^{d}g^{*}_{i,\alpha_{i}}(x_{i})-\sum_{i=1}^{d}\\{\alpha_{i}/(k_{i}+1)!\\}x_{i}^{k_{i}+1}.$

The constant offset $\mu^{*}$ is introduced to ensure that
$\int_{0}^{1}g^{*}_{i,\alpha_{i}}(x)\,dx=0$ for all $i=1,\dots,d$. By applying
the analogous arguments to the proofs of Theorems 1 and 2, we can characterize
the risk of the proposed estimator for a fixed
$\alpha=(\alpha_{1},\dots,\alpha_{d})$ using the LSEs of
$\mu^{*}+\sum_{i=1}^{d}g^{*}_{i,\alpha_{i}}(x_{i})$. Formally, we define the
LSE of the additive function for given $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ with
$X_{i}=(X_{i1},\dots,X_{id})$ as

$\displaystyle(\widehat{\mu},\widehat{g}_{1},\dots\widehat{g}_{d}):=\operatorname*{arg\,min}\sum_{i=1}^{n}\left\\{Y_{i}-\mu-\sum_{j=1}^{d}g_{j}(X_{ij})\right\\}^{2}$

where $\operatorname*{arg\,min}$ is taken over $\mu\in\mathbb{R}$ and
$g_{j}\in\mathcal{G}_{j}$ for $j=1,\dots,d$. We define an oracle estimator of
the $\ell$th component of the additive function as follows:

$\displaystyle\widetilde{g}_{\ell}:=\operatorname*{arg\,min}_{g\in\mathcal{G}}\sum_{i=1}^{n}\left\\{Y_{i}-\mu^{*}-\sum_{j\neq\ell}g^{*}_{j,\alpha_{j}}(X_{ij})-g(X_{i\ell})\right\\}^{2}.$

This estimator is not practically feasible as it assumes the knowledge of
$\mu^{*}$ and $g^{*}_{j,\alpha_{j}}$ for $j\neq\ell$. Under the assumption
that the covariate space $\Omega\subseteq\mathbb{R}^{d}$ forms a Cartesian
product set $\mathcal{X}_{1}\times\dots\times\mathcal{X}_{d}$ and
$\mathcal{X}_{i}\subseteq\mathbb{R}$ for all $i=1,\dots,d$, Lemma 3.1 of
Guntuboyina and Sen, (2018) states that $\widetilde{g}_{i}=\widehat{g}_{i}$
for $i=1,\dots,d$. With this general result, the risk behavior of the LSE
$\widehat{g}_{i}$ of $g^{*}_{i,\alpha_{i}}$ is characterized by the behavior
of the oracle estimator $\widetilde{g}_{i}$. Thus, the convergence rates
provided by Theorems 1 and 2 are applicable to the context of the additive
model. In the subsequent section, we demonstrate that our estimator for
additive models, as outlined by Algorithm 2, achieves the convergence rate
prescribed by Theorems 1 and 2 when each component of additive functions
satisfies the necessary shape constraint.

##  6 Numerical Studies

###  6.1 Univariate Nonparametric Regression

We conduct numerical studies to assess the finite-sample properties of the
proposed procedures. We first consider cases with univariate covariates where
$X$ follows a uniform distribution $\mbox{Unif}[0,1]$. The response variables
are generated according to $Y_{i}=f(X_{i})+\xi_{i}$ for $i=1,2,\dots,n$ with
several different regression functions $f$. The noise $\xi_{i}$’s are IID.
$N(0,0.1^{2})$ variables across all scenarios. The samples sizes vary from
$10^{2}$ to $10^{4}$. For each scenarios, we repeat the experiment 300 times.

For the proposed method, we split the index set $\\{1,2,\dots,n\\}$ into two
disjoint subsets $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ such that
$\\{(X_{i},Y_{i}):i\in\mathcal{I}_{2}\\}$ contains $\sqrt{n}$ observations. We
also compare our method with other existing nonparametric estimators:

  1. 1.

Kernel ridge regression (KRR): The kernel function is given by
$K(x,z):=1+\min(x,z)$, which corresponds to the first-order Sobolev space
(i.e., Example 12.16 of Wainwright, (2019)). We use $10$ grids for the
penalization tuning parameter of ridge regression. Due to the computationally
intensive nature of KRR, we do not explore finer choices of the tuning
parameters.

  2. 2.

Gradient boosting machines (GBM): The shrinkage parameter is set to 0.01, and
we choose the maximum depth of each tree from the set $\\{2,5\\}$. The total
number of trees is selected from $\\{100,1000,2000,4000,8000\\}$. The
remaining parameters are set to their default values according to the GBM
library in R.

  3. 3.

Random forest regression (RF): The number of trees is selected from
$\\{50,100,500,1000,5000\\}$. We use the estimator implemented by R package
randomForest (Liaw and Wiener,, 2002).

  4. 4.

Penalized sieve estimator with cosine basis: We employ 50 basis functions and
select the penalization tuning parameter from (approximately) 100 default
grids. The estimators are realized using R package Sieve (Zhang and Simon,,
2022).

For all procedures, the model selection is performed through cross-validation
based on the validation set $\\{(X_{i},Y_{i}):i\in\mathcal{I}_{2}\\}$ with a
sample size of $\sqrt{n}$. For sample sizes exceeding 2000, we omit the
results obtained from KRR and RF due to computational limitations. We consider
the following four “true” regression functions:

#### Scenario 1

We examine a Lipschitz function defined as:

$f_{1}(x):=(1-3x)I_{[0,1/3]}(x)+(-1+3x)I_{[1/3,2/3]}(x)+(3-3x)I_{[2/3,1]}(x).$
(11)

The proposed estimator, along with other nonparametric regression estimators,
is expected to converge at a rate of $n^{-2/3}$ in terms of their (integrated)
squared risk.

#### Scenario 2

We consider a case where our proposed estimators are expected to be locally
adaptive, converging at a parametric rate of $n^{-1}$ (up to a logarithmic
term). Define $M_{m}(x):=\sum_{i=1}^{m}iI_{[(i-1)/m,i/m]}(x)$. The proposed
estimator is anticipated to be adaptive to the following truth function:

$f_{2}(x;m,\beta):=M_{m}(x)+\beta x.$ (12)

For this specific scenario, we use $m=3$ and $\beta=1$.

#### Scenario 3

The next two scenarios focus on the application of the convex regression. The
following example corresponds to a smooth function, in the Hölder sense,
defined as:

$f_{3}(x;\gamma):=\sin(\gamma(2x-1)).$

For this particular scenario, we select $\gamma=4$. The proposed estimator is
anticipated to converge at a rate of $n^{-4/5}$ ignoring a logarithmic term.

#### Scenario 4

The final example illustrates another scenario where our proposed estimator is
expected to be locally adaptive. Specifically, the estimator is designed to
adapt to any function that can be decomposed as a sum of a convex
$m$-piecewise linear function and a quadratic term. We define $C_{m}(x)$ as a
convex $m$-piecewise linear function with $1/m$ equally sized segments over
$X$, and the slopes are determined as $(-1,0,1,\dots,m-2)$. Additionally, we
enforce the condition $C_{m}(0)=0$, thereby defining a unique convex
$m$-piecewise linear function. We generate observations from one such
regression defined as:

$f_{4}(x;m,\beta):=C_{m}(x)+\beta x^{2}.$

In this particular scenario, we consider the case with $m=3$ and $\beta=1$.

Figure 1: Estimated regression functions based on $500$ observations. The left
two panels illustrate the estimated regression functions derived from an
isotonic regression and a linear parameter, while the right two panels
illustrate the estimations based on a convex regression and a quadratic
parameter.

Figure 1 displays a single sample path of the estimated regression function
using a sample size of $n=500$. The red line represents the estimated
regression function, while the black line represents the true regression
function. The estimator appears to be consistent with the true curves,
including the regression functions that contain non-differentiable points
(Scenario 1) as well as discontinuities (Scenario 2).

(a) Average MSEs across 300 replications for different sample sizes. The solid
black line represents the expected convergence rate of the proposed method.

(b) Boxplot of MSEs across 300 replications for sample sizes $n=100,1000$, and
$10000$.

Figure 2: Performance comparison. Univariate cases.  | Scenario 1 | Scenario 2 | Scenario 3 | Scenario 4
---|---|---|---|---
Theoretical rate |  $-$0.667 |  $-$1 |  $-$0.8 |  $-$1
LSE + parametric |  $-$0.738 |  $-$1.201 |  $-$0.799 |  $-$0.895
Sieve |  $-$0.449 |  $-$0.276 |  $-$0.349 |  $-$0.710
GBM |  $-$0.568 |  $-$0.720 |  $-$0.590 |  $-$0.645
Table 1: Estimated slopes of Figure 2(a) based on linear regression using
observations with $n\geq 2000$. The displayed values correspond to the
estimated rate $\gamma$ for $n^{\gamma}$.

Next, we study the convergence rate of mean squared errors (MSE) for various
methods as the sample sizes vary. To estimate the MSEs, we generate $10^{5}$
new data points and calculate the prediction errors. Figure 2(a) displays the
average MSEs for various sample sizes, while Figure 2(b) displays box plots
representing the MSE distributions specifically for sample sizes of
$n=100,1000$, and $10000$. The solid black lines in Figure 2(a) represent the
advertised rates of convergence for the proposed method, namely $n^{-2/3}$ for
Scenario 1, $n^{-4/5}$ for Scenario 3, and $n^{-1}$ (ignoring a logarithmic
rate) for Scenarios 2 and 4. Additionally, Table 1 presents the estimated
slope based on linear regression using observations with sample sizes $n\geq
2000$.

We observe that in the small-sample regimes, the proposed method deviates from
the theoretical rate of convergence. However, as the sample size increases,
the method aligns more closely with the expected rate of convergence. In
Scenario 1, most of the nonparametric regression methods demonstrate
comparable performance, with the proposed method performing particularly well
for larger sample sizes. Similar conclusions can be drawn for Scenario 3, with
the performance of the proposed method being especially pronounced. In the two
scenarios where the proposed method is expected to be adaptive, exhibiting a
parametric convergence rate, it outperforms all other methods for sample sizes
larger than 1000. These findings highlight the adaptability and effectiveness
of the proposed method as the sample size increases.

###  6.2 Multivariate Regression with Additive Structure

In this section, we consider the multivariate extensions assuming an additive
structure for the underlying regression functions, as discussed in Remark 2
from Section 3. For this set of numerical studies, we also investigate the
Generalized Additive Model (GAM) with the dimension of basis set to $30$. We
first consider two examples with 2-dimensional covariates $x=(x_{1},x_{2})$:

  * •

Scenario 1 (2d): $f(x):=f_{1}(x_{1})-f_{1}(x_{2})$,

  * •

Scenario 2 (2d): $f(x):=f_{2}(x_{1};3,1)-f_{2}(x_{2};3,1)$,

where the component functions $f_{1},f_{2}$ are defined in (11) and (12).
Similar to the univariate cases, we anticipate that our proposed method will
converge essentially at a rate of $n^{-2/3}$ for Scenario 1 and $n^{-1}$ for
Scenario 2.

We also consider two 5-dimensional examples:

  * •

Scenario 3 (5d): $f(x):=f_{1}(x_{1})-f_{1}(x_{2})+x_{3}-x_{4}+1$,

  * •

Scenario 4 (5d):
$f(x):=f_{2}(x_{1};1,0)+f_{2}(1-x_{2};3,3)+f_{2}(x_{3};3,3)+f_{2}(1-x_{4};1,3)+f_{2}(x_{5};1,3)$.

We also expect our proposed method to converge at a rate of $n^{-2/3}$ for
Scenario 3 and $n^{-1}$ for Scenario 4.

(a) Average MSEs across 300 replications. The solid black line represents the
expected convergence rate of the proposed method.

(b) Boxplot of MSEs across 300 replications for sample sizes $n=100,1000$, and
$10000$.

Figure 3: Performance comparison. Additive multivaraite cases.

In the case of additive functions, Figure 3(a) displays the average MSEs while
Figure 3(b) displays box plots representing the distribution of observed MSEs.
Additionally, Table 2 presents the estimated slope based on linear regression
using observations with sample sizes $n\geq 2000$. Similar to the univariate
case, the theoretical rate of convergence aligns with the empirical behavior
for larger sample sizes. We observe that the generalized additive model
performs well for regression functions even with non-differentiable points
(Scenarios 1 and 3). However, it struggles to accurately estimate functions
with discontinuities (Scenarios 2 and 4). In contrast, our proposed method
performs well when each additive component is a non-decreasing piecewise
constant function plus a linear term.

| Scenario 1 (2d) | Scenario 2 (2d) | Scenario 3 (5d) | Scenario 4 (5d)
---|---|---|---|---
Theoretical rate |  $-$0.667 |  $-$1 |  $-$0.667 | $-1$
LSE + parametric |  $-$0.717 |  $-$0.885 |  $-$0.828 |  $-$0.814
Sieve |  $-$0.553 |  $-$0.528 |  $-$0.577 |  $-$0.571
GBM |  $-$0.068 |  $-$0.633 |  $-0$.107 |  $-$0.563
GAM |  $-$0.738 |  $-$0.089 |  $-$0.780 |  $-$0.104
Table 2: Estimated slopes of Figure 3(a) based on linear regression using
observations with $n\geq 2000$. The displayed values correspond to the
estimated rate $\gamma$ for $n^{\gamma}$.

###  6.3 Robustness to the Tuning Parameters

Next, we examine the robustness of the proposed method regarding two aspects:
(1) the specific values of the parameter $\alpha$ and (2) the randomness from
cross-validation. First, we investigate the MSEs of the proposed method when
the parameter $\alpha$ is pre-specified. Hence, the resulting estimator
$\widehat{g}_{\alpha}-\alpha x$ depends on the data only through the estimated
function $\widehat{g}_{\alpha}$. Figure 4 displays the MSEs (on a logarithmic
scale) as the value of $\alpha$ changes. The results are presented for sample
sizes of $n=500,1000$ and $5000$. To recall, Proposition 1 holds for any
$\alpha\geq L$ where $L$ is the true Lipschitz constant, implying that we
expect the MSEs to be robust once the value of $\alpha$ surpasses $L$. This
property is demonstrated in Figure 4.

We also investigate the robustness of the proposed method under the random
split of cross-validation. To study this, we generate the observations once
and then examine the MSEs of the resulting estimators for different cross-
validation splits. We repeat this process 300 times, each time using a
different split, and then examine the MSEs. The top plot of Figure 5 displays
the box plots of the MSEs while the bottom plot shows the box plots of the
values of $\alpha$ obtained from the cross-validation splits. We observe that
the MSEs across different splits are concentrated even for small sample sizes.
For example, in Scenario 1 with $n=500$, the majority of MSEs fall within the
range of $10^{-3.25}$ to $10^{-3}$, indicating very small variability between
cross-validation splits. Similar conclusions can be drawn for different
scenarios. The bottom plot of Figure 5 displays that the variability of
$\alpha$ between random splits is also small. This indicates the robustness of
the proposed procedure in the selection of $\alpha$.

Figure 4: The MSEs of the proposed estimator as the values of $\alpha$ change.
Different line types represent results for different sample sizes. In Scenario
1, the Lipschitz constant is 3, while in Scenario 2, the Lipschitz constant is
32.

Figure 5: The box plots of the MSEs and the estimated parameter across
different random splits during the cross-validation step. The initial
observations are fixed for each same size, and the variability of statistics
is induced solely by the random split of data.

###  6.4 Additional Results on Adaptive Rates

Lastly, we investigate the behavior of the proposed estimator in scenarios
where the method is expected to converge at the faster rate. In particular,
Scenario 2 represents non-decreasing $m$-piecewise constant functions with
linear terms and Scenario 3 represents $m$-piecewise convex affines with
quadratic terms. Theorem 2 implies that the rate of convergence of the
proposed estimator in these scenarios is expected to display a linear
relationship as the number of $m$ changes. To verify this property, we
generated $5000$ observations from Scenarios 2 and 4 over the values of $m$ in
$\\{1,2,3,4,5\\}$. The results are presented in Figure 6. As expected, the
average MSEs over 300 repetitions increase linearly as the value of $m$
increases.

Figure 6: The MSEs of the proposed estimator for two scenarios ($n=5000$) with
expected parametric rate of convergence: $m$-piecewise linear functions and
$m$-piecewise convex affines. The convergence rate is expected to exhibit a
linear relationship with respect to the number of segments $m$.

##  7 Conclusion

This article introduces a new class of nonparametric regression estimators
that leverages existing shape-restricted regression methods. The proposed
method takes advantage of the decompositions of nonparametric functions within
the Hölder class into shape-restricted and parametric components. We then
propose an estimation procedure based on sample-splitting, which practically
eliminates the turning parameter. Our proposed method inherits favorable
properties from shape-restricted regression estimators, including efficient
computation and optimal convergence rates. Moreover, our estimator
demonstrates adaptability to specific regression functions such that it
automatically converges at the parametric rate (up to a logarithmic term).

In our subsequent investigation, we establish more rigorous theoretical
properties, specifically the rate of convergence under dependent,
heteroscedastic and heavy-tail noises (Kuchibhotla and Patra,, 2022).
Moreover, we provide these theoretical results without assuming a fixed value
for $\alpha$, thereby formally proving the properties of the proposed
procedure while accounting for the randomness introduced by the sample-
splitting steps. In practice, the uncertainty quantification is crucial in
regression analysis and the construction of valid confidence intervals is
desired. While precisely characterizing the limiting distribution in shape-
restricted regression is challenging due to irregularity, we anticipate that
the recently developed HulC procedure (Kuchibhotla et al.,, 2021) can be
applicable in this context.

## Acknowledgements

The authors gratefully acknowledge support from NSF DMS-2210662.

## References

  * Adams and Fournier, (2003) Adams, R. A. and Fournier, J. J. (2003).  Sobolev spaces.  Elsevier.
  * Ayer et al., (1955) Ayer, M., Brunk, H. D., Ewing, G. M., Reid, W. T., and Silverman, E. (1955).  An empirical distribution function for sampling with incomplete information.  The annals of mathematical statistics, pages 641–647.
  * Balabdaoui et al., (2019) Balabdaoui, F., Durot, C., and Jankowski, H. (2019).  Least squares estimation in the monotone single index model.
  * Balabdaoui and Wellner, (2005) Balabdaoui, F. and Wellner, J. A. (2005).  Estimation of a k-monotone density, part 1: characterizations consistency, and minimax lower bounds.  arXiv preprint math/0509080.
  * Bellec, (2018) Bellec, P. C. (2018).  Sharp oracle inequalities for least squares estimators in shape restricted regression.  The Annals of Statistics, 46(2):745–780.
  * Byrd and Nocedal, (1989) Byrd, R. H. and Nocedal, J. (1989).  A tool for the analysis of quasi-newton methods with application to unconstrained minimization.  SIAM Journal on Numerical Analysis, 26(3):727–739.
  * Chatterjee et al., (2015) Chatterjee, S., Guntuboyina, A., and Sen, B. (2015).  On risk bounds in isotonic and other shape restricted regression problems.
  * Chen and Samworth, (2016) Chen, Y. and Samworth, R. J. (2016).  Generalized additive and index models with shape constraints.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 78(4):729–754.
  * Dette et al., (2006) Dette, H., Neumeyer, N., and Pilz, K. F. (2006).  A simple nonparametric estimator of a strictly monotone regression function.  Bernoulli, 12(3):469–490.
  * Eguchi et al., (2003) Eguchi, S., Yoon Kim, T., and Park, B. U. (2003).  Local likelihood method: a bridge over parametric and nonparametric regression.  Nonparametric Statistics, 15(6):665–683.
  * Gao and Wellner, (2009) Gao, F. and Wellner, J. A. (2009).  On the rate of convergence of the maximum likelihood estimator of ak-monotone density.  Science in China Series A: Mathematics, 52(7):1525–1538.
  * Ghosal and Sen, (2017) Ghosal, P. and Sen, B. (2017).  On univariate convex regression.  Sankhya A, 79(2):215–253.
  * Groeneboom et al., (2001) Groeneboom, P., Jongbloed, G., and Wellner, J. A. (2001).  Estimation of a convex function: characterizations and asymptotic theory.  The Annals of Statistics, 29(6):1653–1698.
  * Guntuboyina and Sen, (2015) Guntuboyina, A. and Sen, B. (2015).  Global risk bounds and adaptation in univariate convex regression.  Probability Theory and Related Fields, 163(1-2):379–411.
  * Guntuboyina and Sen, (2018) Guntuboyina, A. and Sen, B. (2018).  Nonparametric shape-restricted regression.  Statistical Science, 33(4):568–594.
  * Han et al., (2019) Han, Q., Wang, T., Chatterjee, S., and Samworth, R. J. (2019).  Isotonic regression in general dimensions.
  * Han and Wellner, (2018) Han, Q. and Wellner, J. A. (2018).  Robustness of shape-restricted regression estimators: An envelope perspective.  arXiv preprint arXiv:1805.02542.
  * Han and Wellner, (2019) Han, Q. and Wellner, J. A. (2019).  Convergence rates of least squares regression estimators with heavy-tailed errors.
  * Han and Zhang, (2020) Han, Q. and Zhang, C.-h. (2020).  Limit distribution theory for block estimators in multiple isotonic regression.  The Annals of Statistics, 48(6):3251–3282.
  * Hjort and Glad, (1995) Hjort, N. L. and Glad, I. K. (1995).  Nonparametric density estimation with a parametric start.  The Annals of Statistics, pages 882–904.
  * Hjort and Jones, (1996) Hjort, N. L. and Jones, M. C. (1996).  Locally parametric nonparametric density estimation.  The Annals of Statistics, pages 1619–1647.
  * Kinnunen, (2017) Kinnunen, J. (2017).  Sobolev spaces.  Department of Mathematics and Systems Analysis, Aalto University.
  * Kuchibhotla et al., (2021) Kuchibhotla, A. K., Balakrishnan, S., and Wasserman, L. (2021).  The hulc: Confidence regions from convex hulls.  arXiv preprint arXiv:2105.14577.
  * Kuchibhotla and Patra, (2022) Kuchibhotla, A. K. and Patra, R. K. (2022).  On least squares estimation under heteroscedastic and heavy-tailed errors.  The Annals of Statistics, 50(1):277–302.
  * Lei, (2020) Lei, J. (2020).  Cross-validation with confidence.  Journal of the American Statistical Association, 115(532):1978–1997.
  * Lepskii, (1991) Lepskii, O. (1991).  On a problem of adaptive estimation in gaussian white noise.  Theory of Probability & Its Applications, 35(3):454–466.
  * Liaw and Wiener, (2002) Liaw, A. and Wiener, M. (2002).  Classification and regression by randomforest.  R News, 2(3):18–22.
  * Lin et al., (2009) Lin, L., Cui, X., and Zhu, L. (2009).  An adaptive two-stage estimation method for additive models.  Scandinavian journal of statistics, 36(2):248–269.
  * Mammen, (1991) Mammen, E. (1991).  Estimating a smooth monotone regression function.  The Annals of Statistics, pages 724–740.
  * Mammen and Yu, (2007) Mammen, E. and Yu, K. (2007).  Additive isotone regression.  Lecture Notes-Monograph Series, pages 179–195.
  * Meyer, (2008) Meyer, M. C. (2008).  Inference using shape-restricted regression splines.  The Annals of Applied Statistics, pages 1013–1033.
  * Seijo and Sen, (2011) Seijo, E. and Sen, B. (2011).  Nonparametric least squares estimation of a multivariate convex regression function.
  * Seregin and Wellner, (2010) Seregin, A. and Wellner, J. A. (2010).  Nonparametric estimation of multivariate convex-transformed densities.  Annals of statistics, 38(6):3751.
  * Veraverbeke et al., (2014) Veraverbeke, N., Gijbels, I., and Omelka, M. (2014).  Preadjusted non-parametric estimation of a conditional distribution function.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 399–438.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional statistics: A non-asymptotic viewpoint, volume 48.  Cambridge university press.
  * Wang and Ghosal, (2023) Wang, K. and Ghosal, S. (2023).  Posterior contraction and testing for multivariate isotonic regression.  Electronic Journal of Statistics, 17(1):798–822.
  * Zhang, (2002) Zhang, C.-H. (2002).  Risk bounds in isotonic regression.  The Annals of Statistics, 30(2):528–555.
  * Zhang and Simon, (2022) Zhang, T. and Simon, N. (2022).  Regression in tensor product spaces by the method of sieves.  arXiv preprint arXiv:2206.02994.
  * Zlobec, (2006) Zlobec, S. (2006).  Characterization of convexifiable functions.  Optimization, 55(3):251–261.