The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
94,890 characters
Multiple Testing of a Function's Monotonicity
\title{Multiple Testing of a Function's Monotonicity}
\author[1]{\fnm{Wei} \sur{Zhao}}
\author*[1]{\fnm{David M.} \sur{Kaplan}}\email{[email removed]}
\affil*[1]{\orgdiv{Department of Economics}, \orgname{University of Missouri}, \orgaddress{\street{615 Locust St}, \city{Columbia}, \postcode{65211}, \state{MO}, \country{USA}}}
\abstract{
Instead of having a single ``yes'' or ``no'' result from a test of the global null hypothesis that a function is increasing, we propose a multiple testing procedure of the function's increasingness at several points.
If the global null is rejected, then multiple testing provides more information about why.
If the global null is not rejected, then multiple testing can provide stronger evidence in favor of increasingness, by rejecting null hypotheses that the function is decreasing.
Our approach uses high-level assumptions that apply to a broad class of causal and descriptive statistical models.
By inverting the proposed multiple testing procedure that controls the familywise error rate, we also generate ``inner'' and ``outer'' confidence sets for the set of points at which the function is increasing.
With high asymptotic probability, the inner confidence set is contained within the true set, whereas the outer confidence set contains the true set.
We also improve power with stepdown and two-stage procedures.
Simulation and empirical examples illustrate the new methodology, and all code is provided.
}
\pacs[MSC Classification]{62J15}
\keywords{Familywise error rate,
Inner confidence set,
Multiple testing procedure,
Outer confidence set}
\maketitle
\vspace*{-2\baselineskip}
Date: August 31, 2025\\[\baselineskip]
This version of the article has been accepted for publication, after peer review, but is not the Version of Record and does not reflect post-acceptance improvements or any corrections.
The Version of Record is available online at \url{https://doi.org/10.1007/s11749-025-00986-6}.
Use of this Accepted Version is subject to the publisher's Accepted Manuscript terms of use: \url{https://www.springernature.com/gp/open-research/policies/accepted-manuscript-terms}
\clearpage
\section{Introduction}
\label{introduction}
Unlike a global test of the single null hypothesis that a whole function is monotonically increasing, we propose a multiple testing procedure (MTP) for a function's increasingness at multiple points.
These multiple points correspond to multiple null hypotheses, for which our procedure controls the familywise error rate.
Testing where a function is decreasing is equivalent to testing where the negative of that function is increasing, so we focus on increasingness without loss of generality.
A global hypothesis test with a single null can only provide a ``yes or no'' result.
That is, it can only say either ``yes, we have enough evidence to reject the null that the entire function is increasing'' or ``no, we cannot reject the null that the entire function is increasing.''
In contrast, our MTP can provide more information by testing increasingness at several specific points.
Notationally, we refer to these points as values of ``$X$,'' but $X$ is not restricted to be a regressor; for example, it can represent quantile index values in a quantile regression.
The MTP assesses at which values of $X$ the function is increasing.
For example, imagine a function that is generally increasing except for a small decreasing segment near $X=x$.
The global null hypothesis should be rejected because it is false, but this result arguably misses the big picture that the function is mostly increasing.
Instead, the MTP can reject increasingness near $X=x$ specifically, without rejecting increasingness elsewhere, to provide a more precise assessment.
Beyond non-rejection, the MTP can provide even stronger evidence in favor of increasingness at certain points.
We know that ``non-rejection of a null'' is relatively weak evidence in favor of the null, given the possibility of a type II error (failure to reject false null), whose rate is not controlled.
However, if we reverse the original nulls from ``increasing'' to ``decreasing,'' then rejection provides strong statistical evidence in favor of increasingness at those points.
The MTP may still make an error, but it would be a type I error (rejection of true null), whose rate is controlled at the desired level, in the familywise sense.
Assessing a function's increasingness is a general statistical task with applications across many fields.
We provide a few examples here.
In finance, the below-mentioned work by \citet{RomanoWolf2013} is motivated by theories about financial returns increasing in a variety of asset or portfolio characteristics.
In public health, ``gradients'' encompass the positive or negative relationship between health and various socioeconomic variables, like the health--income gradient; for example, \citet{RidleyEtAl2020} review evidence on the causal relationship between poverty and mental health.
Generally, these gradients can be either descriptive or causal; our methods can apply to either.
In economics, theories related to information costs and rational inattention have implications for monotonicity of certain conditional probabilities, some of which have been tested experimentally; for example, see Section 5.3 of \citet{DewanNeligh2020} and Experiment 1.2 of \citet{DeanNeligh2023}, which both test for the correct choice probability increasing in the incentive level.
In climatology, \citet{FriedrichEtAl2020} test for atmospheric ethane increasing over time during some but not all periods, for which our multiple testing approach could inductively find such periods, rather than running a global test on a pre-specified period.
In ecology, a review paper by \citet{ZhangEtAl2015eco} focuses on non-monotonicity in many ecological settings, where our multiple testing approach can help gauge evidence for such non-monotonicity, including not merely its existence but specifically where in a function it occurs.
Additionally, the MTP can be useful for helping assess identifying assumptions.
For example, \citet{ManskiPepper2000} discuss partial identification of average treatment effects given various monotonicity assumptions, including monotone treatment selection combined with monotone treatment response.
This combination has a testable implication that the conditional expectation function is a weakly increasing function; see their (19) and footnote 9.
For example, in the empirical analysis of the return to schooling, the MTP results would be useful for assessing the combination of monotone treatment selection and monotone treatment response, which jointly imply the expectation of log wage increases with years of schooling.
Instead of using a global hypothesis test and abandoning whole assumptions because they are violated somewhere, the MTP can be used to find subpopulations where the assumptions are met, and then the causal analysis can be restricted to such subpopulations.
More specifically, we could either optimistically use the subpopulations where we do not reject mean log wage increasing with education, or we could more conservatively use the subpopulations where we reject that mean log wage decreases with education in favor of increasingness.
These subpopulations respectively correspond to our outer and inner confidence sets described below.
Inverting the MTP, we also construct ``inner'' and ``outer'' confidence sets.
The object of interest is the true set of points at which the function is increasing.
The inner confidence set is contained within the true set with high asymptotic probability.
Similarly, the outer confidence set contains the true set with high asymptotic probability.
The inner and outer confidence set mechanics come from \citet{Kaplan2024}, who builds on other work to develop them in the context of sets of utility functions.
Theoretically, we provide a unified framework that applies to a wide array of models.
Specifically, we assume the relevant points on the function have estimators that are jointly asymptotically normal, with a consistently estimable covariance matrix.
For example, this applies to many estimators of functions that have a causal interpretation, such as an average structural function \citep[][\S8.1]{BlundellPowell2003}, average index function \citep[][\S8]{LewbelEtAl2012}, or structural quantile function \citep[][\S3.1]{ImbensNewey2009}, which could be estimated by instrumental variables or control function methods.
This also applies to estimators of descriptive models, such as conditional mean and quantile functions.
As noted earlier, the function does not even need to be a function of a regressor, such as looking at a quantile regression slope coefficient as a function of the quantile index.
Our use of the covariance matrix allows our MTP to have asymptotically exact familywise error rate in the ``least favorable'' case of a flat function, whereas the Bonferroni approach generally is conservative.
Further, we offer two ways to improve power.
First, we use a stepdown procedure, like that of \citet{Holm1979}.
Second, we propose two-stage procedure following the approach of \citet{RomanoShaikhWolf2014}, using the least favorable null within a first-stage confidence set.
These two approaches have complementary strengths, and we show simulation results for a data-generating process where their combined power improvement exceeds the sum of their individual improvements.
The main limitation of our MTP is the use of a finite number of $X$ points.
However, this still allows $X$ to be continuous, as long as the number of evaluations points is finite, for example by using deciles.
However, with continuous $X$, the optimal choice of evaluation points remains an open question.
Naturally, the $X$ variable can also be discrete, or even ordinal, in which case potentially all possible values can be used.
Although we are not aware of any existing multiple testing methodology to assess monotonicity, there are related literatures on both monotonicity testing and multiple testing.
In the special case of a conditional mean function, some related papers propose tests of ``regression monotonicity'' with a single null.
For example, \citet{GhosalEtAl2000} and \citet{Chetverikov2019} have the single null ``$H_0\colon m(\cdot)$ is an increasing function'' for the conditional mean function $m(x)=\operatorname{E}(Y\mid X=x)$.
\Citet{HallHeckman2000} design a test statistic for the single null ``$H_0\colon m(\cdot)$ is nondecreasing on the interval $\mathscr{I}$'' without requiring computation of a curve estimator.
\Citet{RomanoWolf2013} also test a single null, the negation of the alternative hypothesis that financial returns are strictly increasing in a given characteristic; see their (2.11).
Alternatively, \citet{KostyshakLuo2021} propose a ``partial monotonicity parameter'' that measures the proportion of the population for which an increase in $X$ is associated with an increase in $Y$, complementing our methods that focus on \emph{where} in the domain the function is increasing.
On the multiple testing side, the general strategy for our basic MTP is similar to that of \citet{WuKaplan2025a}, who instead consider stochastic monotonicity with ordinal or discrete outcomes.
Besides the difference in application to stochastic monotonicity instead of function monotonicity, which entails different equations and formulas (like our \cref{supp_eqn:Zhat-asy-dist} vs.\ their Lemma 1), they do not provide refinements to improve power like in our \cref{sec:power}, and their approach with continuous outcomes uses a different approach.
There is also a literature on testing moment inequalities in general, some of which we apply here; for example, see Section 4 of the survey by \citet{CanayShaikh2017}, as well as the work of \citet{RomanoShaikhWolf2014} that we adapt to our setting.
\paragraph{Paper structure}
\Cref{sec:MTP,sec:CS} describe new methodology and provide formal theoretical results.
\Cref{sec:emp-CEF,sec:sim} contain empirical and simulation results, respectively.
\Cref{sec:power} provides two additional methods to improve power.
In the supplemental appendix, \cref{supp_sec:app-ext} shows how our methodology can apply to instrumental variables quantile regression as well as functional coefficient models, and \cref{supp_sec:app-proofs} collects proofs.
\paragraph{Notation and abbreviations}
Random and non-random vectors are respectively typeset as, e.g., $\bm{X}$ and $\bm{x}$,
while random and non-random scalars are typeset as $X$ and $x$,
and random and non-random matrices as $\mkern3mu\underline{\mkern-3mu \bm{X}\mkern-3mu}\mkern3mu$ and $\mkern3mu\underline{\mkern-3mu \bm{x}\mkern-3mu}\mkern3mu$.
Acronyms used include those for
conditional expectation function (CEF),
confidence set (CS),
familywise error rate (FWER),
multiple testing procedure (MTP),
pointwise rejection probability (PRP),
and
two-stage least squares (2SLS).
\section{Multiple Testing}
\label{sec:MTP}
This section presents our MTP and establishes its strong control of asymptotic FWER.
\subsection{Setting}
Consider learning whether or not a function $m(\cdot)$ is increasing.
This $m(\cdot)$ is a function of a scalar $x$, but the underlying model may include interactions between $X$ and other observed and unobserved variables, whose realizations are fixed when defining $m(\cdot)$.
We consider testing at $h$ different values of $X$.
These points are not necessarily the full support of $X$; for example, if $X$ is continuous, we could test at deciles of $X$.
For notational simplicity, these points are written as $x\in \mathcal{X}= \{1,\ldots,h \}$, but more generally $x=j$ can be interpreted as $x=x_j$ for arbitrary $\mathcal{X}=\{x_1,\ldots,x_h\}$.
The MTP decides whether or not to reject each hypothesis under consideration.
Specifically, the MTP tests the following $h-1$ nulls $H_{0x}$ simultaneously.
Defining $d_x\equiv m(x)-m(x+1)$,
\begin{equation}
\label{eqn:H0x}
\begin{split}
& H_{01}\colon d_1 \le 0\quad \textrm{(the function $m(\cdot)$ is increasing from $x=1$ to $2$)}\\
& H_{02}\colon d_2 \le 0\quad \textrm{(the function $m(\cdot)$ is increasing from $x=2$ to $3$)}\\
& ~~\vdots\\
&H_{0,h-1}\colon d_{h-1} \le 0\quad \textrm{(the function $m(\cdot)$ is increasing from $x=h-1$ to $h$)} .
\end{split}
\end{equation}
The MTP makes $h-1$ decisions, with $2^{h-1}$ possible outcomes.
In \cref{sec:ex-CEF}, we introduce a running example.
\subsubsection{Example: conditional expectation function (CEF)}
\label{sec:ex-CEF}
Define the conditional expectation function (CEF) as $m(x)\equiv\operatorname{E}(Y\mid X=x)$, $x\in \mathcal{X}=\{1,2,\dots,h\}$.
The null hypotheses are
\begin{equation*}
\begin{aligned}
&H_{01}\colon d_1=m(1)-m(2)=\operatorname{E}(Y\mid X=1)-\operatorname{E}(Y\mid X=2)\le0,\\
&H_{02}\colon d_2=\operatorname{E}(Y\mid X=2)-\operatorname{E}(Y\mid X=3)\le0,\\
& ~~\vdots\\
&H_{0,h-1}\colon d_{h-1}=\operatorname{E}(Y\mid X=h-1)-\operatorname{E}(Y\mid X=h)\le0.
\end{aligned}
\end{equation*}
Again, more generally, we could use any $x_1<\cdots<x_h$ and write the nulls as
\begin{equation*}
\begin{split}
H_{0j}\colon d_j&\le0
\textrm{ for }j=1,\ldots,h-1
,\\%\quad
&d_j \equiv m(x_j)-m(x_{j+1})
= \operatorname{E}(Y\mid X=x_j) - \operatorname{E}(Y\mid X=x_{j+1})
.
\end{split}
\end{equation*}
The MTP tests these $h-1$ nulls simultaneously.
The results help us learn whether the CEF is increasing or not at each point.
\subsection{Benefits of multiple testing for monotonicity}
\label{sec:MTP-benefits}
There are two main benefits of multiple testing in this context.
First, simply knowing that the global hypothesis test rejected or not does not have as much information as the MTP results.
The global hypothesis test only says whether the whole function $m(\cdot)$ is increasing or not.
The MTP further tells us in which region it is increasing or not.
For example, the global hypothesis test result ``reject $H_0$'' means that the function $m(\cdot)$ is not increasing everywhere in the domain, but that does not provide any information about where the increasing trend is broken.
In contrast, even the simple MTP result ``reject $H_{01}$ but not $H_{02}$'' not only tells us that the function $m(\cdot)$ is not increasing overall, but specifically says there is enough evidence to reject that it increases from $m(1)$ to $m(2)$.
Second, the MTP also can provide strong evidence in favor of $m(\cdot)$ increasing at certain points, by switching the direction of the null hypothesis.
In the global hypothesis test, non-rejection of $H_0\colon m(\cdot)\textrm{ is increasing everywhere}$ is relatively weak evidence in favor of increasingness, and in practice it is rare to have the strength of evidence required to reject the reversed null ``$H_0^*\colon m(\cdot)\textrm{ is not increasing everywhere}$.
Instead, with multiple testing of the reversed nulls
\begin{equation*}
H_{0x}^*\colon d_x\ge0\textrm{ for }x\in\mathcal{X} ,
\end{equation*}
the MTP can provide strong evidence in favor of $m(\cdot)$ increasing at certain points, even if there is not strong evidence at every single point.
\subsection{Assumptions}
\label{sec: assumption}
\Cref{a:asy-normal} is maintained throughout.
\begin{assumption}
\label{a:asy-normal}
The estimator $\hat{\bm{m}}\equiv(\hat{m}(1),\ldots,\hat{m}(h))'$ is asymptotically normal: given true value $\bm{m}\equiv(m(1),\ldots,m(h))'$,
\begin{equation*}
\sqrt{n}(\hat{\bm{m}}-\bm{m})
\xrightarrow{d} \mathrm{N}(\bm{0},\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a) ,
\end{equation*}
where positive definite matrix $\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a$ can be estimated consistently, $\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^a\xrightarrow{p}\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a$.
\end{assumption}
\Cref{a:asy-normal} has a couple notational simplicities that immediately generalize.
First, the specific $\sqrt{n}$ convergence rate does not affect any proof and can be replaced by $n^r$ for some $r>0$.
For example, this allows nonparametric estimators, where $\bm{m}$ is a finite-dimensional functional of some underlying population function, but its estimator $\hat{\bm{m}}$ need not have a $\sqrt{n}$ rate, like if $\bm{m}$ is the evaluation functional (points on the underlying function); for more discussion and technical results, see Section 3.4 of \citet{Chen2007}.
Second, recall that $m(j)$ can be interpreted as $m(x_j)$, in which case \Cref{a:asy-normal} refers to asymptotic normality of $\hat{\bm{m}}\equiv(\hat{m}(x_1),\ldots,\hat{m}(x_h))'$ for arbitrary $x_1<\cdots<x_h$.
\Cref{a:asy-normal} can allow for time series data.
In standard cases, the key is to use a covariance matrix estimator that is robust to dependence, such as the class of heteroskedasticity and autocorrelation consistent (HAC) estimators proposed and studied by \citet{NeweyWest1987} and \citet{Andrews1991HAC} and implemented in the R package \lstinline{sandwich} \citep{R.sandwich}.
For example, to estimate the CEF $\operatorname{E}(Y\mid X=x)$ with discrete or ordinal $X$, we can regress $Y$ on indicator variables for each $X$ value and no constant, and given assumptions on finite moments, stationarity, dependence, and the kernel function and lag length, Theorem 1(a) of \citet{Andrews1991HAC} establishes consistency of this class of covariance matrix estimators.
Conversely, of course, if there is nonstationarity or dependence is too strong, then our \Cref{a:asy-normal} does not hold.
For example, with a first-order autoregressive AR(1) process $y_{t}=\phi y_{t-1}+\varepsilon_{t}$, where $\varepsilon_{t}$ is white noise with variance $\sigma^2$, the OLS estimator $\hat{\phi}$ is not asymptotically normal if there is a unit root, $\phi=1$.
An interesting special case is when we wish to learn if the mean of a time series is increasing over time, i.e., monotonicity of $\operatorname{E}(Y_t)$ in $t$.
We can partition the time series sample $\{Y_t\}_{t=1}^{T}$ into $h$ blocks of consecutive observations and compute the sample mean of each block; equivalently, we can regress $Y_t$ on the $h$ block indicator variables and no constant, the corresponding estimated coefficients of which are $\hat{m}(j)$ for $j=1,\ldots,h$.
Asymptotically, if we consider a fixed number of observations per block and thus growing number of blocks $h$, then generally \Cref{a:asy-normal} does not hold; for example, the block averages $\hat{m}(j)$ are not normally distributed, unless the $Y_t$ are Gaussian already.
Alternatively, if we have a fixed number of blocks $h$ and growing number of observations per block, then \Cref{a:asy-normal} may indeed hold.
In practice, \citet{IbragimovMueller2010} emphasize that we should think about which asymptotic approximation seems more appropriate for our particular sample, given the number of observations per block, for example.
They further consider when blocks are large enough that they are approximately independent, so that cross-block correlations are negligible.
This also depends on the empirical setting; for example, a climate analysis may have significant dependence across decades, whereas stock market returns may have negligible serial correlation at even a daily sampling frequency.
Instead of \Cref{a:asy-normal}, we could make the high-level assumption that the asymptotic distribution of $n^r(\hat{\bm{m}}-\bm{m})$ exists and can be consistently estimated by a bootstrap or subsampling.
In that case, the critical value can be computed as in \cref{sec:cv-bootstrap}, which can also provide some finite-sample improvement even under \Cref{a:asy-normal}.
Thus, our approach can be applied in certain non-normal cases, too.
\Cref{a:asy-normal} is a relatively high-level assumption that holds in a wide variety of settings.
Below, we continue our example from \cref{sec:ex-CEF} to show how to establish this high-level assumption from lower-level assumptions.
\subsubsection{Example: CEF (continued)}
Continuing from \cref{sec:ex-CEF}, in the CEF example of $m(x)=\operatorname{E}(Y\mid X=x)$ for $x\in\{1,\ldots,h\}$, the following low-level assumptions are sufficient for \Cref{a:asy-normal}.
\begin{assumpC}
\label{a:CEF}
The subsample sizes are $n_x\equiv\sum_{i=1}^{n}\operatorname{\mathds{1}}\{X_i=x\}$, and they are fixed, with iid sampling from the respective conditional distributions.
For example, there are $n_1$ iid draws from the conditional distribution of $Y$ given $X=1$, $n_2$ iid draws from the conditional distribution of $Y$ given $X=2$, and generally $n_x$ iid draws from the conditional distribution of $Y$ given $X=x$.
Also, $n_1/n_x\to\gamma_x\in(0,\infty)$ for all $x \in \mathcal{X}$.
\end{assumpC}
Under \Cref{a:CEF}, using the central limit theorem and continuous mapping theorem \citep[e.g.,][Prop.\ 2.27 and Thm.\ 2.3]{vanderVaart1998},
the asymptotic distribution of $\sqrt{n_1} (\hat{\bm{m}}-\bm{m})$ is
\begin{equation}
\begin{aligned}
\label{eqn:CEF-mhat-asy-dist}
\sqrt{n_1} (\hat{\bm{m}}-\bm{m})
\xrightarrow{d} \mathrm{N}(\bm{0},\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a)
\end{aligned}
\end{equation}
as in \Cref{a:asy-normal}, where $\bm{0}\equiv(0,\dots,0)'$, $V^a_{(xj)}= \gamma_x v_x \operatorname{\mathds{1}}\{x=j\}$ represents the row $x$, column $j$ element of diagonal matrix $\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a$ for $x,j\in\mathcal{X}$, and $v_x\equiv\operatorname{Var}(Y\mid X=x)$.
\subsection{Multiple testing procedure}
\label{sec:plain-MTP}
Below, we describe our MTP and state its asymptotic FWER control.
Our MTP controls the ``overall type I error rate'' known as the familywise error rate (FWER).
This is one way to quantify the false positive rate that we want to control when looking across a family of hypotheses.
There are several alternatives to FWER, like $k$-FWER and false discovery proportion \citep{LehmannRomano2005fwer}, as well as the false discovery rate \citep{BenjaminiHochberg1995}, but such are left to future work.
One benefit of FWER control is that it justifies inversion of our MTP into confidence sets as in \cref{sec:CS}.
\Cref{def:FWER,def:FWER-control} follow \citet[\S9.1, p.\ 407]{LehmannRomano2022text}, jointly defining our MTP's desired property of strong control of asymptotic FWER.
\begin{definition}[familywise error rate]
\label{def:FWER}
For a family of null hypotheses $H_{0x}$ indexed by $x$, let $\mathcal{T}\equiv\{x:H_{0x}\textrm{ is true}\}$ be the set of indices of true hypotheses.
The FWER is the probability of rejecting any true null hypothesis.
Mathematically,
\begin{equation}
\label{eqn:FWER}
\textrm{FWER} \equiv
\operatorname{P}(\text{reject any }H_{0x}\text{ with }x \in \mathcal{T}) .
\end{equation}
\end{definition}
\begin{definition}[strong control of FWER]
\label{def:FWER-control}
Strong control of FWER requires $\textrm{FWER}\le\alpha$ for any combination of true and false $H_{0x}$.
Strong control of asymptotic FWER instead requires $\textrm{FWER}\le\alpha+o(1)$.
\end{definition}
\begin{method}
\label{meth:plain-MTP}
Given the $H_{0x}\colon d_x\le0$ in \cref{eqn:H0x}, the MTP rejects $H_{0x}$ when $\hat{t}_x>c_\alpha$, where the $t$-statistics are $\hat{t}_x=\hat{d}_x/\hat{s}_x$, the estimated standard errors $\hat{s}_x$ are the square roots of the diagonal elements of $\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^b/n$, where $\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^b$ is a consistent estimator for the $\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^b$ defined in \cref{supp_eqn:dhat-asy-dist}, and $c_\alpha$ is the critical value.
Specifically, $c_\alpha$ is the $(1-\alpha)$-quantile of the asymptotic distribution of $\hat{Z}^*$, which is the maximum of $h-1$ correlated random variables $\hat{Z}_x$:
\begin{equation*}
\hat{Z}^*\equiv\max_{x\in\{1,2,\dots,h-1\}} \hat{Z}_x
,\quad
\hat{Z}_x \equiv \frac{\hat{d}_x-d_x }{\hat{s}_x} \textrm{ for }x\in\{1,2,\dots,h-1\}
,
\end{equation*}
and the asymptotic joint normal distribution of the $\hat{Z}_x$ is in \cref{res:asy-Zhat}.
Details of simulating $c_\alpha$ are in \cref{sec:cv-sim,sec:cv-bootstrap}.
To test $H_{0x}^*\colon d_x\ge0$, simply run the above MTP for $H_{0x}\colon d_x\le0$ after replacing $Y$ with $-Y$ in the data.
\end{method}
\begin{lemma}
\label{res:asy-Zhat}
Under \Cref{a:asy-normal}, the asymptotic distribution of $\hat{\bm{Z}}\equiv ( \hat{Z}_1, \dots, \hat{Z}_{h-1} )$ is multivariate normal with mean zero and covariance matrix $\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu$ defined in \cref{supp_eqn:Zhat-asy-dist}: $\hat{\bm{Z}} \xrightarrow{d} \mathrm{N} (\bm{0},\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu)$.
\end{lemma}
\begin{theorem}
\label{res:plain-MTP}
Under \Cref{a:asy-normal}, \cref{meth:plain-MTP} has strong control of asymptotic FWER.
\end{theorem}
Although the null hypothesis is very different, \cref{meth:plain-MTP} can be compared to the global monotonicity test in Section 3.1 of \citet{RomanoWolf2013}.
(Their Section 3.3 method further provides a two-step procedure analogous to our \cref{meth:RSW}.)
The consider a more specific setting where $m(j)$ is the expected financial return of asset or portfolio $j$, and $\hat{m}(j)$ is a time series average of returns, but their approach readily generalizes to our setting.
Translating to our setting and notation, their (2.11) states their null and alternative hypotheses $H_0\colon \max_x d_x \ge 0$ vs.\ $H_1\colon \max_x d_x < 0$; that is, $H_1$ is that the function is strictly increasing, and $H_0$ is the negation of $H_1$.
Their Section 3.1 test rejects $H_0$ in favor of $H_1$ when $\max_x \hat{t}_x<z_\alpha$, the $\alpha$-quantile of the standard normal distribution, like $z_{0.05}\approx -1.64$.
This is more closely related to our MTP for the reversed nulls $H_{0x}^*\colon d_x\ge0$ from \cref{sec:MTP-benefits}, which rejects $H_{0x}^*$ when $\hat{t}_x < -c_\alpha$.
In the trivial case of $h=2$, $-c_\alpha=z_\alpha$, so both methods reject when $\hat{t}_1 < -c_\alpha$.
With $h>3$, there are a few different cases.
If $\max_x\hat{t}_x < -c_\alpha$, then the \citet{RomanoWolf2013} test rejects in favor of a strictly increasing $m(\cdot)$, and our method rejects each null in favor of $m(\cdot)$ strictly increasing between every pair of consecutive points.
If $-c_\alpha < \max_x\hat{t}_x < z_\alpha$, then their global test concludes that the function is strictly increasing everywhere, but there is at least one point where our MTP cannot reject in favor of strictly increasing.
In the most extreme case, $-c_\alpha < \min_x\hat{t}_x < \max_x\hat{t}_x < z_\alpha$, so none of the $H_{0x}^*$ are rejected by the MTP in favor of strict increasingness, yet the global test can reject the global null in favor of $m(\cdot)$ strictly increasing everywhere, which is equivalent to all $H_{0x}^*$ being false.
This illustrates that for the global null, the global test of \citet{RomanoWolf2013} has a power advantage over multiple testing.
Conversely, if $\min_x\hat{t}_x < -c_\alpha < z_\alpha < \max_x\hat{t}_x$, then their global test cannot reject in favor of increasingness, but our MTP can reject in favor of increasingness at certain $x$ values.
In one extreme example, imagine $\hat{t}_1=0$ but $\hat{t}_x < -c_\alpha$ for all $x>1$: then our MTP rejects in favor of $m(\cdot)$ strictly increasing over all $x\ge2$, but the global test merely fails to reject the global null of not-everywhere-increasing.
\footnote{Similar conclusions should apply to the analogous comparison of the global test of \citet{DavidsonDuclos2013} taking restricted first-order stochastic dominance as the alternative hypothesis with the MTP of \citet{GoldmanKaplan2018c}, where instead of $H_{0x}\colon m(x)-m(x+1)\le0$ over discrete $x$ they have $H_{0x}\colon F_1(x)-F_2(x)\le0$ over $x\in{\mathbb R}$ for CDFs $F_1(\cdot)$ and $F_2(\cdot)$, among other variations.}
So, the global test and multiple testing have complementary strengths in power.
The above begs the question: can the global test be combined with the closure method to make a higher-powered MTP?
The answer is no, for the following reason.
In order to reject $H_{0j}^*$ for a particular $j$, the closure method requires the global test to jointly reject every possible combination of $H_{0x}^*$ hypotheses including the particular $H_{0j}^*$.
This includes the joint test of all $H_{0x}^*$.
However, this joint test fails to reject when $\max_x\hat{t}_x > z_\alpha$, i.e., when even a single $t$-statistic is higher than the $z_\alpha$ critical value.
Thus, a necessary condition for \emph{any} $H_{0x}^*$ to be rejected by this closure method is for $\max_x\hat{t}_x < z_\alpha$, so the MTP reduces to ``reject all $H_{0x}^*$'' when $\max_x\hat{t}_x<z_\alpha$, and ``do not reject any $H_{0x}^*$'' otherwise.
Note that this inability to leverage the closure method is due to the increasing-everywhere being the alternative hypothesis; if it were the null hypothesis, and the joint test statistic is the minimum of the individual $t$-statistics, then applying the closure method simply reproduces our MTP in \cref{meth:plain-MTP}.
\subsubsection{Critical value: multivariate normal simulation}
\label{sec:cv-sim}
The critical value $c_\alpha$ must be simulated because the analytic formulation of the asymptotic distribution of $\hat{Z}^*$ is intractable.
\Cref{meth:plain-MTP} says $c_\alpha$ is the $(1-\alpha)$-quantile of the asymptotic distribution of $\hat{Z}^*$, which is a distribution of the maximum of $h-1$ correlated normal random variables with mean zero and covariance matrix $\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu$.
\Citet{NadarajahKotz2008} derive the exact distribution of the max of two Gaussian random variables, but it is very complex and in practice usually there are more than two.
To simulate $c_\alpha$, we can approximate the asymptotic distribution of $\hat{Z}^*$ using the following method.
\begin{method}[simulated normal critical value]
Run the following steps.
\begin{steps}
\item For $b=1,\ldots,B$, draw $\bm{Z}^{(b)}\stackrel{\mathit{iid}}{\sim}\mathrm{N}(\bm{0},\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu})$, where the covariance matrix estimator is defined just after \cref{supp_eqn:Zhat-asy-dist}.
\item Given each $\bm{Z}^{(b)}=(Z^{(b)}_1,\ldots,Z^{(b)}_{h-1})$, compute the maximum $Z^{(b)}_\text{max}\equiv\max_{x\in\{1,\ldots,h-1\}}Z^{(b)}_x$.
\item The simulated critical value $\hat{\hat{c}}_{\alpha}$ is the $(1-\alpha)$-quantile among the $B$ values of $Z^{(b)}_\text{max}$.
\end{steps}
\end{method}
The first ``hat'' on $\hat{\hat{c}}_{\alpha}$ represents the use of estimated covariance matrix $\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}$, whose estimation error converges in probability to zero as the sample size increases.
The second ``hat'' is because it is simulated.
The simulation error can be made arbitrarily small by taking a large enough number of simulation replications.
To choose the number of replications $B$ in practice, we have the following suggestion.
Specifically, we solve for the smallest $B$ such that there is a high probability $1-\epsilon$ that FWER does not exceed $\alpha+\delta$, abstracting away from potential estimation error in $\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}$ and asymptotic approximation error in the normal distribution.
For example, if $\epsilon=0.02$ and $\delta=0.01$, then we want a $98\%$ probability that FWER is at most one percentage point larger than the desired $\alpha$.
Let $c_{\alpha+\delta}$ be the critical value that achieves exactly $\alpha+\delta$ FWER under the least favorable null.
The simulated critical value will be at least $c_{\alpha+\delta}$ when no more than $(1-\alpha)R$ of the $Z^{(b)}_\text{max}$ are below $c_{\alpha+\delta}$.
Because the $Z^{(b)}_\text{max}$ are iid across $b=1,\ldots,B$, and the true probability of $Z^{(b)}_\text{max}\le c_{\alpha+\delta}$ is $1-\alpha+\delta$ under the least favorable null by definition of $c_{\alpha+\delta}$, then $\sum_{b=1}^{B}\operatorname{\mathds{1}}\{Z^{(b)}_\text{max}\le c_{\alpha+\delta}\} \sim \textrm{Binomial}(B,1-\alpha+\delta)$, a binomial distribution that depends only on the user-chosen values $B$, $\alpha$, and $\delta$.
Thus, to ensure a probability of at least $1-\epsilon$ that the FWER distortion due to simulation error does not exceed $\delta$, we can solve for the smallest $B$ that satisfies
\begin{equation*}
\operatorname{P}\bigl( \textrm{Binomial}(B,1-\alpha-\delta) \le (1-\alpha)B \bigr)
\ge 1-\epsilon .
\end{equation*}
For example, using $\alpha=0.05$, $\delta=0.01$, and $\epsilon=0.02$ yields $R=2180$; decreasing to $\epsilon=0.01$ yields $R=2820$.
Our simulations with $\alpha=0.05$ use $B=10^5$, which is near the solution $B=10{,}800$ when $\delta=0.005$ and $\epsilon=0.01$.
\subsubsection{Critical value: nonparametric bootstrap}
\label{sec:cv-bootstrap}
Alternatively, the critical value can be simulated using a bootstrap or subsampling method.
To give a concrete example, here we describe a Studentized nonparametric bootstrap appropriate for iid data.
We use this method in the simulation in \cref{sec:sim}.
\begin{method}[bootstrap critical value]
Run the following steps.
\begin{steps}
\item Let $\hat{d}_x$, $\hat{s}_x$, and $\hat{t}_x=\hat{d}_x/\hat{s}_x$ ($x=1,\ldots,h-1$) denote the estimates, standard errors, and $t$-statistics from the original sample of observations $\bm{W}_i$ ($i=1,\ldots,n$).
\item\label{step:cv-bs-resample} From the original sample, resample $\bm{W}^*_i$ ($i=1,\ldots,n$) with replacement.
\item\label{step:cv-bs-t} Using the bootstrap sample, compute estimates, standard errors, and $t$-statistics $\hat{m}^*_x$ ($x=1,\ldots,h$), $\hat{d}^*_x=\hat{m}^*_x-\hat{m}^*_{x+1}$, $\hat{s}^*_x$ (estimated standard error of $\hat{d}^*_x$), and $\hat{t}^*_x=(\hat{d}^*_x-\hat{d}_x)/\hat{s}^*_x$ ($x=1,\ldots,h-1$), as well as the maximum $t$-statistic $T^*=\max_x\hat{t}^*_x$. (Note $\hat{t}^*_x$ is centered at $\hat{d}_x$, not zero, where in the ``bootstrap world'' $\hat{d}_x$ is the true population value; $\hat{Z}^*_x$ is perhaps better notation than $\hat{t}^*_x$ but would partially conflict with our earlier $\hat{Z}$ notation.)
\item Repeat \cref{step:cv-bs-resample,step:cv-bs-t} $B$ times, to get $T^*_b$ for $b=1,\ldots,B$. (In \cref{sec:sim}, we use $B=1000$.)
\item The bootstrap critical value $\hat{c}_\alpha$ is the $(1-\alpha)$-quantile of the $T^*_b$ values.
\item As in \cref{meth:plain-MTP}, reject $H_{0x}$ when $\hat{t}_x>\hat{c}_\alpha$, where $\hat{t}_x$ is from the original sample.
\end{steps}
\end{method}
\subsubsection{Example: CEF (continued)}
In the CEF example, we have asymptotic normality of $\sqrt{n} (\hat{\bm{m}}-\bm{m})$ in \cref{eqn:CEF-mhat-asy-dist}, which in turn implies asymptotic normality of $\hat{\bm{Z}}$.
By sampling vectors from a mean-zero multivariate normal distribution with the appropriate estimated covariance matrix $\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}$, the simulated critical value is the $(1-\alpha)$-quantile among the simulated maxima.
\section{Confidence Sets}
\label{sec:CS}
Each confidence set (CS) that we propose is for a set-valued parameter, specifically the set of points $x$ at which $m(x)$ is increasing:
\begin{equation}
\label{eqn:true-S}
\mathcal{S}
\equiv \{ x : d_x \le 0 \}
= \{ x : m(x)-m(x+1) \le 0 \} .
\end{equation}
\subsection{Outer confidence set}
\label{outer confidence set}
An outer CS $\hat{\mathcal{S}}_o$ is the more common type of CS, containing the true set $\mathcal{S}$ with high asymptotic probability.
Mathematically,
\begin{equation}
\label{eqn:outer-CS}
\operatorname{P}(\hat{\mathcal{S}}_o\supseteq\mathcal{S})\ge1-\alpha+o(1) .
\end{equation}
Faced with uncertainty, an outer CS must err toward being too big and including some $x$ that are not in fact in $\mathcal{S}$.
Put differently, the outer CS must have strong evidence in order to exclude a certain $x$ from $\hat{\mathcal{S}}_o$, otherwise it will violate \cref{eqn:outer-CS}.
\Cref{meth:outer-CS} constructs an outer CS by inverting the MTP in \cref{meth:plain-MTP}, similar to Method 3 of \citet{Kaplan2024} and following the same logic as the (outer) confidence set for the identified set developed by \citet[Lem.~2.1]{RomanoShaikh2010}.
Its asymptotic validity is then given by \cref{res:outer-CS}.
\begin{method}[outer CS]
\label{meth:outer-CS}
Run the MTP in \cref{meth:plain-MTP} with $H_{0x}\colon d_x \le 0$ and let the outer CS be $\hat{\mathcal{S}}_o \equiv \{x : \textrm{MTP does not reject}\ H_{0x} \}$.
\end{method}
\Cref{meth:outer-CS} says that the outer CS collects all $x$ for which the MTP does not reject the corresponding null $H_{0x}$.
Thus, the true set will be a subset of the outer CS as long as none of the true null hypotheses is rejected.
As $1-\alpha$ decreases and thus $\alpha$ increases, each $H_{0x}$ becomes more likely to be rejected, so the outer CS becomes smaller.
\begin{corollary}
\label{res:outer-CS}
Given the conclusion of \cref{res:plain-MTP} and true set $\mathcal{S}$ in \cref{eqn:true-S}, the outer CS $\hat{\mathcal{S}}_o$ in \cref{meth:outer-CS} is asymptotically valid in the sense of \cref{eqn:outer-CS}.
\end{corollary}
\subsection{Inner confidence set}
\label{inner confidence set}
An inner CS $\hat{\mathcal{S}}_i$ is contained within the true set $\mathcal{S}$ with high asymptotic probability.
This idea seems to originate in (1) of \citet{ArmstrongShen2023}, and we generally follow the approach of \citet[][\S\S3,5.3]{Kaplan2024}.
Mathematically,
\begin{equation}
\label{eqn:inner-CS}
\operatorname{P}(\hat{\mathcal{S}}_i\subseteq\mathcal{S})
\ge 1 - \alpha + o(1) .
\end{equation}
Faced with uncertainty, the inner CS does the opposite of the outer CS, erring toward being too small and possibly omitting some $x$ that actually are in $\mathcal{S}$.
That is, the inner CS must have strong evidence in order to \emph{include} an $x$ value in $\hat{\mathcal{S}}_i$, otherwise it will violate \cref{eqn:inner-CS}.
Thus, the inner CS can be viewed as a more conservative ``estimate'' of $\mathcal{S}$ in the sense that it only includes $x$ values where there is strong evidence that $m(x)$ is increasing.
\Cref{meth:inner-CS} constructs an inner CS by inverting the MTP in \cref{meth:plain-MTP} with the null hypothesis inequality directions reversed, again following Method 3 of \citet{Kaplan2024}; its asymptotic validity is then given by \cref{res:inner-CS}.
\begin{method}[inner CS]
\label{meth:inner-CS}
Run the MTP in \cref{meth:plain-MTP} with $H_{0x}^*\colon d_x \ge 0$ (i.e., replacing $Y$ with $-Y$) and let the inner CS be $\hat{\mathcal{S}}_i \equiv \{x : \textrm{MTP rejects }H_{0x}^* \}$.
\end{method}
\Cref{meth:inner-CS} says that the inner CS collects all $x$ for which the MTP rejects the corresponding reversed null hypothesis $H_{0x}^*\colon d_x \ge 0$ in favor of increasingness at point $x$, so the inner CS equals the complement of the outer CS for $\{x:m(x)\le m(x+1)\}$, the set of points where $m(x)$ is decreasing.
Intuitively, strong evidence in favor of increasingness is the same as strong evidence against decreasingness, so the set of $x$ with strong evidence in favor of increasingness is the complement of the set of $x$ lacking strong evidence against decreasingness; that is, the inner CS for increasing points is the complement of the outer CS for decreasing points.
This symmetry is reassuring.
Regardless of whether we frame our inquiry in terms of increasingness or decreasingness, the confidence sets essentially partition the $x$ points into three groups: one with strong evidence in favor of increasingness, one with strong evidence in favor of decreasingness, and one without strong evidence in either direction.
As with the outer CS, the coverage of the inner CS is closely linked to the MTP.
Because the inner CS collects $x$ for which the MTP rejects $H_{0x}^*\colon d_x \ge 0$, the probability of any false rejection is the probability of incorrectly including an $x$ in the inner CS.
Thus, asymptotically, strong control of FWER implies correct coverage probability.
Also, as $1-\alpha$ decreases and thus $\alpha$ increases, each $H_{0x}^*$ becomes more likely to be rejected, so the inner CS becomes larger.
That is, as $\alpha$ increases, both the inner CS and outer CS approach the point estimate $\hat{\mathcal{S}}=\{x:\hat{d}_x\le0\}$, but the inner CS approaches from the inside whereas the outer CS approaches from the outside, with $\hat{\mathcal{S}}_i\subseteq\hat{\mathcal{S}}\subseteq\hat{\mathcal{S}}_o$.
\begin{corollary}
\label{res:inner-CS}
Given the conclusion of \cref{res:plain-MTP} and true set $\mathcal{S}$ in \cref{eqn:true-S}, the inner CS $\hat{\mathcal{S}}_i$ in \cref{meth:inner-CS} is asymptotically valid in the sense of \cref{eqn:inner-CS}.
\end{corollary}
\section{Empirical illustration: earnings and education}
\label{sec:emp-CEF}
Economists and others have long studied the relationship between earnings and years of education.
Here we consider not a causal relationship but a descriptive one, the conditional expectation function (CEF).
Further, as noted in \cref{introduction}, this CEF can also be used to assess the assumptions of \citet{ManskiPepper2000} for learning about the causal effect of education on earnings.
\subsection{Setup}
\label{sec:emp-CEF-setup}
In this example, $m(x)=\operatorname{E}(Y\mid X=x)$, and interest is in where log weekly income ($Y$) is increasing in years of education ($X$).
We use the dataset \lstinline{census2000} in package \lstinline{wooldridge} \citep{R.wooldridge}, originally provided by \citet{Wooldridge2010}.
It contains $\num[round-mode=none,group-digits=integer]{29501}$ observations.
In the dataset, our $Y$ variable is \lstinline{lweekinc} and our $X$ variable is \lstinline{educ}.
In \citeauthor{ManskiPepper2000}'s (\citeyear{ManskiPepper2000}) empirical example, they compute a uniform $95\%$ confidence band for the CEF $m(\cdot)$.
Then, they claim that the monotone treatment selection and monotone treatment response assumptions are consistent with the empirical evidence because the band includes a weakly increasing CEF.
However, such a band could include an increasing function even if the entire estimated CEF is decreasing.
The MTP can provide a stronger and more specific assessment.
\subsection{Results}
\begin{table}[htbp]
\centering
\caption{\label{tab:emp-CEF}CEF example, $\alpha=0.05$}
\sisetup{round-precision=3}
\begin{tabular}[c]{
r
S[table-format=1.2,round-precision=2]
S[table-format=5.0,round-precision=0]
S[table-format=1.2,round-precision=2]
S[table-format=-2.2,round-precision=2]
cc}
\toprule
\multicolumn{1}{c}{$x$} & \multicolumn{1}{c}{$\hat{m}(x)$}
& \multicolumn{1}{c}{$n_x$} & \multicolumn{1}{c}{$\hat{\hat{c}}_\alpha$}
& \multicolumn{1}{c}{$\hat{t}_x$} & \multicolumn{1}{c}{$H_{0x}\colon m(x)\le m(x+1)$}
& \multicolumn{1}{c}{$H_{0x}^*\colon m(x)\ge m(x+1)$} \\
\midrule
9 & 6.2586 & 374 & 2.3952 & -1.1968 & Not reject & Not reject \\
10 & 6.3104 & 621 & 2.3952 & 1.6616 & Not reject & Not reject \\
11 & 6.2475 & 601 & 2.3952 & -9.7814 & Not reject & Reject \\
12 & 6.4991 & 12433 & 2.3952 & -12.0561 & Not reject & Reject \\
13 & 6.6327 & 5424 & 2.3952 & -1.3792 & Not reject & Not reject \\
14 & 6.6552 & 2625 & 2.3952 & -17.7031 & Not reject & Reject \\
16 & 6.9397 & 7423 & &&&\\
\bottomrule
\end{tabular}
\end{table}
\Cref{tab:emp-CEF} includes the following columns.
The $x$ is the years of education.
The $\hat{m}(x)$ column shows the corresponding CEF point estimates, and $n_x$ is the number of observations with $X_i=x$.
The $\hat{\hat{c}}_\alpha$ is the simulated critical value; again, there is just one, repeated across rows for convenience.
The $\hat{t}_x$ column shows the $t$-statistics defined in \cref{meth:plain-MTP}, with more negative values indicating stronger evidence of $m(\cdot)$ decreasing from $x$ to $x+1$, and more positive values indicating stronger evidence of $m(\cdot)$ decreasing.
The $H_{0x}\colon m(x)\le m(x+1)$ column shows the MTP results for the null hypotheses that mean income is increasing in education from $x$ to $x+1$.
Similarly, the $H_{0x}^* \colon m(x)\ge m(x+1)$ column shows the MTP results for the reversed null hypotheses that mean income is decreasing in education from $x$ to $x+1$.
The testing is done separately for each column, so FWER is controlled within each column but not across both columns simultaneously.
\Cref{tab:emp-CEF} shows the following results.
From the $H_{0x}\colon m(x)\le m(x+1)$ column, we see the MTP does not anywhere reject that mean log weekly income is increasing with education.
(The estimated function actually decreases from $x=10$ to $x=11$, but not enough for even a pointwise $t$-test to reject at a $5\%$ level.)
However, we know that ``not rejecting'' is relatively weak evidence.
Switching the direction of null hypotheses produces stronger evidence in favor of increasingness, seen in the last column of the table.
Specifically, the MTP rejected $H_{0x}^*$ in favor of increasingness for $x\in\{11,12,14\}$.
That is, the MTP tells us there is strong evidence in favor of average log weekly income increasing with education for most values above $x\ge11$.
For example, rejecting $H_{0,12}^*$ says that the mean log weekly income is statistically significantly higher for individuals with $x=13$ years of education than individuals with $x=12$ years of education; even accounting for the fact that we are making several such comparisons simultaneously, there is a statistically significant increase associated with even just one year of college.
The one exception among $x\ge11$ is that $H_{013}^*$ is not rejected, even though the relevant subsample sizes $n_{13}$ and $n_{14}$ are relatively large.
The point estimates $\hat{m}(13)=6.63$ and $\hat{m}(14)=6.66$ are too close to reject $H_{013}^*\colon m(13)\ge m(14)$.
\Cref{tab:emp-CEF} also has enough information to compute the global test of \citet{RomanoWolf2013} as well as a Bonferroni MTP.
Given the $\hat{t}_x$ values, the MTP results will be identical for any critical value between $\hat{t}_{10}=1.66$ and $\lvert\hat{t}_{11}\rvert=9.78$, so the Bonferroni MTP gets the same results.
Similar to the simulation results in \cref{tab:sim-plain}, the Bonferroni critical value is very similar to ours in this dataset, differing by only $0.0012$.
For the test of \citet{RomanoWolf2013}, following the discussion after our \cref{res:plain-MTP} and using $\alpha=0.05$, the test rejects in favor of the alternative that $m(\cdot)$ is increasing when $\max_x\hat{t}_x<z_\alpha=-1.64$.
In \cref{tab:emp-CEF}, $\max_x\hat{t}_x=1.66>0$, so their test does not come close to rejecting.
This is reasonable: $\hat{m}(10)>\hat{m}(11)$, so we could not possibly claim that there is strong evidence in favor of $m(\cdot)$ increasing at every $x$.
However, it also misses our more nuanced MTP results that show strong evidence of the CEF increasing at certain $x$ values that include a large proportion of individuals in the sample.
\Cref{tab:emp-CEF} can also be used to construct inner and outer confidence sets for the true set of points $\mathcal{S} \equiv \{x\colon m(x)\le m(x+1)\}$ where mean log weekly income is increasing in education.
Following \cref{meth:outer-CS,meth:inner-CS}, the outer CS $\hat{\mathcal{S}}_o$ collects all $x$ for which the MTP does not reject the null $H_{0x}$ of increasingness, and the inner CS $\hat{\mathcal{S}}_i$ collects all $x$ for which the MTP rejects the reversed null $H_{0x}^*$ of decreasingness.
\Cref{tab:emp-CEF} shows that the MTP does not reject any $H_{0x}$ and rejects $H_{0x}^*$ for $x\in\{11,12,14\}$, so the confidence sets are
\begin{equation*}
\hat{\mathcal{S}}_o= \{9,10,11,12,13,14\}
,\quad
\hat{\mathcal{S}}_i=\{11,12,14\}
.
\end{equation*}
There is a high probability of sampling a dataset in which the true $\mathcal{S}$ is ``between'' the inner and outer CSs, which here would mean $\{11,12,14\}\subseteq\mathcal{S}\subseteq\{9,10,11,12,13,14\}$, suggesting the CEF increases at least at $11$, $12$, and $14$ years, and possibly at additional values, but there is not strong enough evidence to say either way.
For the monotone treatment selection and monotone treatment response assumptions assessment in \citet{ManskiPepper2000}, the results from \cref{tab:emp-CEF} indicate strong evidence that the testable implication of an increasing CEF is satisfied for points $x=11,12,14$.
For the rest of the points, there is not enough evidence to say for sure that it is satisfied, but also not enough evidence to reject it, which would then imply rejecting one or both of their assumptions.
For an empirical analysis relying on monotone treatment response and selection, this might suggest that we use the full dataset ($\hat{\mathcal{S}}_o$) for a main analysis, but restrict to $x\in\hat{\mathcal{S}}_i=\{11,12,14\}$ for a robustness check.
\section{Simulation with \texorpdfstring{\cref{meth:plain-MTP}}{Method \ref{meth:plain-MTP}} MTP}
\label{sec:sim}
In this section, the running CEF example is used to illustrate the finite-sample properties of our proposed MTP with FWER level $\alpha=0.05$.
The critical value is computed two ways: first, simulating from the asymptotic normal distribution as in \cref{sec:cv-sim} with $10^5$ draws, and second, using a nonparametric bootstrap as in \cref{sec:cv-bootstrap} with $10^3$ draws.
We also compare with a Bonferroni critical value.
In our setting, FWER is highest under the least favorable null where the function $m(\cdot)$ is flat, so we use a DGP with $h=4$ and $m(1)=m(2)=m(3)=m(4)$.
If the MTP controls FWER in this case, then the FWER will be even lower with strict inequalities like $m(1)<m(2)<m(3)<m(4)$.
Our MTP code is in R \citep{R.core} and uses the multivariate random normal function from the MASS package \citep{R.MASS}.
Results used R version 4.5.1.
Additional simulations are in \cref{sec:power-sim} to illustrate our alternative MTPs that further improve power.
\subsection{Setup}
\label{sec:sim-setup}
We generate the data as follows.
First, we take $n_x$ observations for each value $X_i=x$ ($x=1,2,3,4$), so there are $\sum_{x=1}^{4}n_x=4n_x$ total observations.
Second, we generate the $Y_i$ given its corresponding $X_i=x$.
The conditional distribution is $Y_i\sim\mathrm{N}(1,\sigma_x)$, where $\sigma_x^2\equiv\operatorname{Var}(Y\mid X=x)$ and either $\sigma_x=1$ or $\sigma_x=x$.
Normality is not the most challenging DGP but is sufficient to show patterns as the sample size increases and between our two critical value methods.
The conditional means are all $m(x)=\operatorname{E}(Y\mid X=x)=1$, which maximizes FWER.
That is, if FWER is controlled here, then it would be even lower with different $m(x)$.
Given the DGP, we compute the MTP properties as follows.
First, for each of $\num[round-mode=none,group-digits=integer]{100000}$ simulated datasets, we run the MTP in \cref{meth:plain-MTP} with $\alpha=0.05$.
Second, we compute the FWER.
Because all the nulls are true, the FWER is the proportion of simulated datasets in which the MTP rejects at least one $H_{0x}$.
Third, we also report the minimum, median, and maximum critical values $\hat{\hat{c}}_\alpha$, among simulation replications for a particular DGP and method, along with the Bonferroni critical value $\Phi^{-1}(1-\alpha/3)=2.128$.
\subsection{Simulation results}
\Cref{tab:sim-plain} shows the following patterns related to critical values.
First, as $n_x$ increases, the ``Normal'' (\cref{sec:cv-sim}) simulated critical value becomes more ``precise'' in the sense of lower variance across replications because the covariance matrix estimator's variance decreases, but even with $n_x=10$ the range is relatively precise: $[2.120,2.138]$.
The ``Bootstrap'' critical value also becomes more precise, but not as much as the Normal.
Second, also as $n_x$ increases, the median Normal critical value does not change, whereas the median bootstrap critical value dereases.
This reflects a finite-sample advantage of the bootstrap to capture non-normal sampling distributions.
In this case, with $n_x=10$ the marginal $t$-statistic distributions have somewhat thicker-than-normal tails and thus warrant a somewhat larger critical value, which is reflected by the median Bootstrap critical value but not the median Normal.
Third, for this DGP, the Bonferroni critical value $\Phi^{-1}(1-\alpha/3)=2.128$ is nearly asymptotically exact.
The Normal critical values are very close to this, even with the small $n_x=10$, and the corresponding FWER is either the same or within $0.001$.
That is, even in this case that heavily favors Bonferroni, our MTP adapts to get nearly the same results.
(This MTP also introduces the framework used for our refined two-stage procedure in \cref{sec:RSW}.)
Fourth, the MTP critical values are all well above the one-sided $t$-test critical value $1.64$, reflecting the fact that naively running individual $\alpha=0.05$ $t$-tests on each $H_{0x}$ would fail to control the FWER at $\alpha=0.05$.
\begin{table}[htbp]
\centering
\caption{\label{tab:sim-plain}Simulation examples, $\alpha=0.05$.}
\sisetup{round-precision=3}
\begin{tabular}[c]{c
S[table-format=5.0,round-precision=0,round-mode=places]
l c
S[table-format=1.3,round-precision=3]
S[table-format=1.3,round-precision=3]
S[table-format=1.3,round-precision=3]
S[table-format=1.3,round-precision=3]
}
\toprule
\multicolumn{1}{c}{$\sigma_x$} & \multicolumn{1}{c}{$n_x$} &
\multicolumn{1}{c}{Method} &
\multicolumn{1}{c}{$\hat{\hat{c}}_\alpha$ $[\min,\text{med},\max]$} &
\multicolumn{1}{c}{\textrm{FWER}} \\
\midrule
$1$ & 10 & Normal & [2.120,2.128,2.138] & 0.069500 \\
$1$ & 10 & Bootstrap & [1.982,2.381,3.511] & 0.039900 \\
$1$ & 10 & Bonferroni & 2.128 & 0.069400 \\
\midrule
$x$ & 10 & Normal & [2.121,2.128,2.138] & 0.071900 \\
$x$ & 10 & Bootstrap & [1.938,2.423,4.169] & 0.041200 \\
$x$ & 10 & Bonferroni & 2.128 & 0.071900 \\
\midrule
$1$ & 100 & Normal & [2.124,2.127,2.135] & 0.052100 \\
$1$ & 100 & Bootstrap & [1.925,2.146,2.385] & 0.049600 \\
$1$ & 100 & Bonferroni & 2.128 & 0.051900 \\
\midrule
$x$ & 100 & Normal & [2.124,2.127,2.135] & 0.051000 \\
$x$ & 100 & Bootstrap & [1.915,2.150,2.385] & 0.049000 \\
$x$ & 100 & Bonferroni & 2.128 & 0.050900 \\
\midrule
$1$ & 1000 & Normal & [2.125,2.128,2.131] & 0.052700 \\
$1$ & 1000 & Bonferroni & 2.128 & 0.052600 \\
\midrule
$x$ & 1000 & Normal & [2.125,2.128,2.131] & 0.054500 \\
$x$ & 1000 & Bonferroni & 2.128 & 0.054400 \\
\midrule
$1$ & 10000 & Normal & [2.126,2.128,2.131] & 0.048400 \\
$1$ & 10000 & Bonferroni & 2.128 & 0.048300 \\
\midrule
$x$ & 10000 & Normal & [2.126,2.128,2.131] & 0.049200 \\
$x$ & 10000 & Bonferroni & 2.128 & 0.049300 \\
\bottomrule
\end{tabular}
\end{table}
\Cref{tab:sim-plain} also shows patterns in the error rates.
First, with the Normal critical value, the FWER is near the nominal $\alpha=0.05$ at larger sample sizes, and even with $n_x=10$ the distortion is modest, with $7\%$ FWER.
Second, the Bonferroni error rates are nearly identical to the Normal rates in every case.
Third, compared to both Normal and Bonferroni, the Bootstrap critical value improves FWER at the smallest $n_x=10$, although even by $n_x=100$ this advantage has essentially disappeared.
(This is why we did not spend several more hours to run Bootstrap on the larger $n_x$ simulations.)
Of course, the specific $n_x$ where the advantage disappears depends on the underlying distribution, and more generally depends on the model and estimator, but generally this shows that it may be worth using bootstrap with smaller samples, but not worth the extra computation time on large samples.
As noted above, the FWER improvement seems due to the larger median Bootstrap critical value that captures some of the non-normality of the finite-sample $t$-statistic distributions.
\section{Procedures to Improve Power}
\label{sec:power}
We develop two ways to improve power without sacrificing strong control of FWER.
\Cref{sec:stepdown} describes a stepdown procedure that improves power when there are multiple false null hypotheses and at least one rejection by the original MTP in \cref{meth:plain-MTP}.
Essentially, the rejected hypotheses can be ignored when recomputing a smaller critical value that increases power against the remaining false hypotheses.
\Cref{sec:RSW} describes a two-stage procedure following the strategy of \citet{RomanoShaikhWolf2014}.
Complementing the stepdown, this two-stage procedure improves power when there are null hypotheses that are true and not binding, meaning the true $d_x$ is strictly below zero.
Roughly speaking, the procedure ``estimates'' how far from binding they are (how far below zero are those $d_x$) and lowers the critical value to account for the corresponding $t$-statistics being less likely to cause false rejections.
However, there can be DGPs for which the two-stage power is worse than the plain MTP in \cref{meth:plain-MTP}, like if all true null hypotheses are indeed binding equalities.
In contrast, by construction, the stepdown's power is always at least as high as the plain MTP's power.
\subsection{Stepdown procedure}
\label{sec:stepdown}
The idea of a stepdown procedure is from \citet{Holm1979}; see also Procedure 9.1.1 of \citet[\S9.1]{LehmannRomano2022text}.
Generally, a stepdown procedure starts by deciding whether or not to reject the hypothesis that has the largest test statistic.
If the corresponding null hypothesis is not rejected, then none of the other nulls are rejected, and the procedure stops.
If the null is rejected, then the critical value is adjusted, and the procedure moves to the next-largest test statistic.
Following this logic, \cref{meth:stepdown} describes our stepdown procedure.
\begin{method}[stepdown]
\label{meth:stepdown}
Run the following steps.
\begin{steps}
\item For iteration $i=0$, run \cref{meth:plain-MTP} with $H_{0x} \colon d_x \le 0$.
Set the iteration counter to $i=1$.
\item\label{meth:stepdown-K} Compute the set $\hat{K}^{(i)}\equiv\{x:H_{0x}\textrm{ not yet rejected}\}$, corresponding to null hypotheses not rejected in any iteration before $i$.
\item\label{meth:stepdown-cv} Compute critical value $c_{\alpha}^{\hat{K}^{(i)}}$ as the $(1-\alpha)$-quantile of the asymptotic distribution of $\max_{x \in \hat{K}^{(i)}} \hat{Z}_x$.
\item Reject any additional $H_{0x}$ for which $\hat{t}_x>c_{\alpha}^{\hat{K}^{(i)}}$.
If there are no additional rejections, or if all $H_{0x}$ have now been rejected, then stop.
Otherwise, increment $i$ by one and return to \Cref{meth:stepdown-K}.
\end{steps}
\end{method}
\Cref{prop: just stepdown} establishes the validity of the stepdown procedure in \cref{meth:stepdown}.
\begin{proposition}
\label{prop: just stepdown}
\Cref{meth:stepdown} tests $H_{0x} \colon d_x\le 0$ across $x\in\{1,2,\ldots,h-1\}$ with strong control of asymptotic \textrm{FWER}\ at level $\alpha$.
\end{proposition}
The upside of the stepdown procedure is that it weakly improves power for any DGP, but it has other limitations.
First, the critical values in \cref{meth:stepdown} are simulated based on the least favorable null with all not-yet-rejected $d_x=0$ in each step, which may be conservative.
Second, if the original MTP does not have at least one rejection, then the stepdown procedure makes no difference.
The method in \cref{sec:RSW} complements the stepdown procedure by addressing both of these limitations.
\subsection{Two-stage procedure}
\label{sec:RSW}
Our two-stage procedure follows \citet{RomanoShaikhWolf2014}.
First, we construct a $1-\beta$ confidence region (CR) for the true vector of differences $d_x=m(x)-m(x+1)$ using a small $\beta<\alpha$.
Second, we run the MTP at level $\alpha-\beta$ with the critical value calibrated to the least favorable null within the CR.
\Cref{sec:RSW-1,sec:RSW-2} detail these two stages.
Our original MTP in \cref{meth:plain-MTP} can be seen as the special case with $\beta=0$, in which case the CR is the entire parameter space and thus always includes the overall least favorable null with all $d_x=0$.
\subsubsection{First stage: confidence region}
\label{sec:RSW-1}
The first stage is to construct a CR that jointly covers the true differences $d_x$ with high asymptotic probability, and then to find the least favorable null within this CR.
\Cref{prop:CR} establishes the asymptotic validity of the CR in \cref{meth:CR}, which implicitly inverts an MTP like \cref{meth:plain-MTP} but reversing the direction of the $H_{0x}$ inequalities, implicitly similar to the CR in (4) of \citet{RomanoShaikhWolf2014}.
\begin{method}[confidence region]
\label{meth:CR}
The CR for the true vector $\bm{d}=(d_1,\ldots,d_{h-1})'$ is
\begin{equation*}
\widehat{\mathrm{CR}} \equiv
\bigl\{ \bm{\delta} : \delta_x \le \hat{d}_x - c_\beta \hat{s}_x \; \forall x\in\{1,\ldots,h-1\} \bigr\}
,
\end{equation*}
where $c_\beta$ is the $\beta$-quantile of the asymptotic distribution of $\min_{x\in\{1,\ldots,h-1\}} \hat{Z}_x$, $d_x\equiv m(x)-m(x+1)$ is the true difference, $\hat{d}_x=\hat{m}(x)-\hat{m}(x+1)$ is the estimated difference, $\hat{s}_x$ is the estimated asymptotic standard error of $\hat{d}_x$, and the asymptotic joint distribution of the $\hat{Z}_x\equiv(\hat{d}_x-d_x)/\hat{s}_x$ is in \cref{supp_eqn:Zhat-asy-dist}.
\end{method}
\begin{proposition}
\label{prop:CR}
Under \Cref{a:asy-normal}, the CR from \cref{meth:CR} covers the true point $\bm{d}$ with asymptotic probability $1-\beta$: $\operatorname{P}(\bm{d}\in\widehat{\mathrm{CR}})=1-\beta+o(1)$.
\end{proposition}
The least favorable null within the CR is characterized as follows.
The CR is an orthant whose corner point has components $\hat{d}_x - c_\beta\hat{s}_x$ for $x=1,\ldots,h-1$.
If the CR indeed contains the true $\bm{d}$, then the least favorable null is the least-negative point inside the intersection of this CR and the null hypothesis region with each $d_x\le0$.
That is, the hypothetical $\bm{d}$ within the CR that would generate the highest FWER is the point
\begin{equation}
\label{eqn:d-hat-star}
\min(\bm{0}, \hat{\bm{d}}^*)
\equiv
\bigl(
\min\{0,\hat{d}_1^*\} ,
\dots,
\min\{0,\hat{d}_{h-1}^*\}
\bigr)'
,\quad
\hat{d}_x^* \equiv \hat{d}_x-c_\beta \hat{s}_x .
\end{equation}
\subsubsection{Second stage: MTP}
\label{sec:RSW-2}
The second stage is to run an MTP like \cref{meth:plain-MTP} but with two modifications.
First, the FWER level is adjusted from $\alpha$ to $\alpha-\beta$ to account for possible error in the first stage.
Second, instead of calibrating the critical value to the overall least favorable null $\bm{d}=\bm{0}$, we calibrate the critical value to the new least favorable null in \cref{eqn:d-hat-star}.
\Cref{meth:RSW} describes this modified MTP, whose asymptotic validity is given in \cref{res:RSW}.
\begin{method}
\label{meth:RSW}
After running \cref{meth:CR}, the MTP rejects $H_{0x} \colon d_x \le 0$ when $\hat{t}_x > \hat{c}$, where the $t$-statistics are $\hat{t}_x=\hat{d}_x/\hat{s}_x$ with estimated standard error $\hat{s}_x$ defined below, and critical value $\hat{c}$ is defined as the $(1-\alpha+\beta)$-quantile of the $\max$ over $h-1$ jointly normal random variables with mean vector $\min(\bm{0}, \hat{\bm{d}}{}^*)/ \hat{\bm{s}}$ and covariance matrix $\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}$ defined below, where
\begin{equation*}
\frac{\min(\bm{0}, \hat{\bm{d}}{}^*)}{\hat{\bm{s}}}
\equiv \biggl(\frac{\min\{0, \hat{d}_1^*\}}{\hat{s}_1},
\dots,
\frac{\min\{0, \hat{d}_{h-1}^*\}}{\hat{s}_{h-1}} \biggr) .
\end{equation*}
That is, $\hat{c}$ is the $(1-\alpha+\beta)$-quantile of $\max \{ \mathrm{N}( \min(\bm{0}, \hat{\bm{d}}{}^*)/\hat{\bm{s}},\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}) \}$, with $\hat{s}_x= \sqrt{\hat{V}^b_{(xx)}/n}$, where $\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^b=\mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu' \hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^a \mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu$ as in \cref{supp_eqn:dhat-asy-dist}, with $\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^a$ the consistent estimator of $\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a$ from \cref{a:asy-normal}.
From \cref{supp_eqn:Zhat-asy-dist} and the text before and after it that defines $\mkern3mu\underline{\mkern-3mu \bm{A}\mkern-3mu}\mkern3mu$ and $\mkern3mu\underline{\mkern-3mu \bm{\hat{A}}\mkern-3mu}\mkern3mu$, $\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu \equiv \mkern3mu\underline{\mkern-3mu \bm{A}\mkern-3mu}\mkern3mu \mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^b \mkern3mu\underline{\mkern-3mu \bm{A}\mkern-3mu}\mkern3mu'=\mkern3mu\underline{\mkern-3mu \bm{A}\mkern-3mu}\mkern3mu \mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu' \mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu^a \mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu \mkern3mu\underline{\mkern-3mu \bm{A}\mkern-3mu}\mkern3mu'$, so $\hat{\mkern3mu\underline{\mkern-3mu \bm{\Sigma}\mkern-3mu}\mkern3mu}= \mkern3mu\underline{\mkern-3mu \bm{\hat{A}}\mkern-3mu}\mkern3mu\mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu'\hat{\mkern3mu\underline{\mkern-3mu \bm{V}\mkern-3mu}\mkern3mu}{}^a \mkern3mu\underline{\mkern-3mu \bm{G}\mkern-3mu}\mkern3mu \mkern3mu\underline{\mkern-3mu \bm{\hat{A}}\mkern-3mu}\mkern3mu'$.
\end{method}
In practice, $\hat{c}$ must be simulated because there is no closed-form expression for $\hat{c}$.
Similar to \cref{sec:cv-sim}, we take random draws of the Gaussian vector with the appropriate mean and covariance, and then we take the $(1-\alpha+\beta)$-quantile of the maxima of the vectors to get the simulated $\hat{c}$, denoted $\hat{\hat{c}}$.
\begin{theorem}
\label{res:RSW}
\Cref{meth:RSW} has strong control of asymptotic FWER at level $\alpha$.
\end{theorem}
\subsection{Simulation with refined MTPs}
\label{sec:power-sim}
We provide a small simulation to illustrate the potential power improvement of the stepdown and two-stage procedures.
We continue the CEF example but with different $m(x)$ values.
We compare the FWER and power of five different methods: the plain MTP (\cref{meth:plain-MTP}), Bonferroni, stepdown (\cref{meth:stepdown}), two-stage (\cref{meth:RSW}), and combined stepdown and two-stage.
As in \cref{sec:sim}, we used version 4.5.1 of R \citep{R.core} and the MASS package \citep{R.MASS}.
\subsubsection{Setup}
\label{sec:power-sim-setup}
The DGP is designed to have some features advantageous for the stepdown procedure and other features advantageous for the two-stage procedure.
We have $h=23$ values, $x\in\{1,\ldots,23\}$.
The stepdown helps when there are some initial rejections, so a smaller critical value can be used in the next iteration.
Thus, the stepdown benefits from some $H_{0x}$ being clearly violated, which in our context means $m(\cdot)$ decreasing steeply at certain points.
In our DGP, $m(\cdot)$ decreases from $x=12$ to $x=23$, with especially steep decreases from $x=13$ to $x=15$.
Complementing these stepdown benefits, the two-stage helps most when there are some $H_{0x}$ that are clearly satisfied, which means $m(\cdot)$ increasing.
In our DGP, $m(\cdot)$ increases from $x=2$ to $x=12$.
Because the stepdown and two-stage procedures complement each other, the combined power improvement is bigger than the individual improvements, rejecting more of the $16\le x\le 23$ points where $m(\cdot)$ is slowly decreasing and thus $H_{0x}$ is moderately violated.
We also have $m(1)=m(2)=0$, so $H_{01}$ is the most susceptible to false rejections.
Altogether, we set
\begin{equation}
\label{eqn:power-sim-m}
m(x) = \begin{cases}
0 & \text{if }\phantom{1}1\le x\le 2 ,\\[-0pt]
10(x-2) & \text{if }\phantom{1}3\le x\le 12 ,\\[-0pt]
100-10(x-12) & \text{if }13\le x\le 15 ,\\[-0pt]
70-0.42(x-15) & \text{if }16\le x\le 23 .
\end{cases}
\end{equation}
The data generation steps are the same as in \cref{sec:sim-setup}.
First, we take $n_x$ observations for each possible $X_i=x$ value, for $n_x\in\{50,100,200\}$ in turn.
Second, we generate the $Y_i$ given its corresponding $X_i=x$.
The conditional distribution is $Y_i\mid X_i\sim\mathrm{N}(m(X_i),1)$, with $m(\cdot)$ from \cref{eqn:power-sim-m}.
We simulate $1000$ datasets and compute the properties of the five methods.
For each simulated dataset, we run the plain MTP (\cref{meth:plain-MTP}), Bonferroni, stepdown (\cref{meth:stepdown}), two-stage (\cref{meth:RSW}), and combined stepdown and two-stage.
Each is run with nominal FWER level $\alpha=0.05$, and with $10^5$ draws to compute the critical value.
The two-stage procedure uses a first-stage $99\%$ CR by setting $\beta=0.01$.
Like before, we compute FWER as the proportion of simulated datasets in which any true null hypothesis is rejected.
Given \cref{eqn:power-sim-m}, $H_{0x}\colon m(x)\le m(x+1)$ is true for $1\le x\le 12$.
We also compute the power.
For each false $H_{0x}$, the simulated pointwise power is the proportion of simulated datasets in which that $H_{0x}$ is rejected.
To summarize power, for each method, we sum its pointwise power across all $13\le x\le 23$.
This can also be interpreted as the average number of false $H_{0x}$ rejected.
In the far right column of \cref{tab:sim2}, this is also expressed as a percentage of the total $11$ false $H_{0x}$.
For example, if on average $5.5$ false $H_{0x}$ are rejected, then this translates to $(5.5/11)\times100\%=50\%$.
\subsubsection{Simulation results}
\begin{table}[htb]
\centering
\caption{\label{tab:sim2}Simulation results, $\alpha=0.05$, $\beta=0.01$.}
\sisetup{round-precision=3}
\begin{tabular}[c]{
rl
S[table-format=1.3,round-precision=3]
S[table-format=2.2,round-precision=2]
S[table-format=2.1,round-precision=1]
}
\toprule
$n_x$ & \multicolumn{1}{c}{Method} & \multicolumn{1}{c}{\textrm{FWER}} &
\multicolumn{1}{c}{Power (sum)} & \multicolumn{1}{c}{Power (\% of $11$)} \\
\midrule
50 & Bonferroni & 0.00400 & 4.96800 & 45.1636 \\
50 & Plain MTP & 0.00400 & 4.99800 & 45.4364 \\
50 & Stepdown & 0.00500 & 5.20100 & 47.2818 \\
50 & Two-stage & 0.00500 & 5.31800 & 48.3455 \\
50 & Combined & 0.01100 & 5.84900 & 53.1727 \\
\midrule
100 & Bonferroni & 0.00400 & 7.46100 & 67.8273 \\
100 & Plain MTP & 0.00400 & 7.47100 & 67.9182 \\
100 & Stepdown & 0.00600 & 7.88100 & 71.6455 \\
100 & Two-stage & 0.00600 & 7.86200 & 71.4727 \\
100 & Combined & 0.02400 & 9.11200 & 82.8364 \\
\midrule
200 & Bonferroni & 0.00100 & 10.31700 & 93.7909 \\
200 & Plain MTP & 0.00100 & 10.32600 & 93.8727 \\
200 & Stepdown & 0.00100 & 10.53400 & 95.7636 \\
200 & Two-stage & 0.00100 & 10.46600 & 95.1455 \\
200 & Combined & 0.02600 & 10.87600 & 98.8727 \\
\bottomrule
\end{tabular}
\end{table}
\Cref{tab:sim2} shows the results, with the following patterns.
First, all five methods control FWER below the desired $\alpha=0.05$ level.
The plain MTP and Bonferroni have the smallest FWER, and the combined method the largest, but still well below $\alpha$.
Unlike in \cref{tab:sim-plain}, where the DGP was the least favorable null that leads to asymptotically exact FWER, the DGP here is far from least favorable in order to study power rather than FWER control, so the FWER will not approach $\alpha$ even with larger samples.
Second, the power improvements can be seen.
The qualitative patterns are the same for each $n_x$, so we focus on the $n_x=100$ results in more detail.
Individually, the stepdown and two-stage procedures both improve the aggregate power by around $0.4$; that is, compared to the plain MTP rejecting $7.5$ false $H_{0x}$ on average, they each reject $7.9$.
The combined improvement is even larger because the stepdown further reduces the critical value after the additional rejections from the two-stage method.
However, recall that the DGP was designed to have features that showcase the stepdown and two-stage procedures, so the improvements will not always be so large.
Finally, we can see that as $n_x$ increases, power increases toward $100\%$.
\section{Conclusion}
We have provided new multiple testing methods and corresponding confidence sets to provide richer results about a function's monotonicity than testing a single global null hypothesis.
Our assumptions cover a wide range of descriptive and causal statistical models.
This work can extend in several directions, including an asymptotically increasing number of evaluation points that leverages \citet{ChernozhukovEtAl2019many}, and the development of Bayesian credible sets.
For continuous $X$, another extension is to combine an existing global monotonicity test with the closure method \citep[e.g.,][\S9.2]{LehmannRomano2022text} as follows.
First, partition the support $\mathcal{X}$ into intervals $\mathcal{X}_j$ for $j=1,\ldots,h$.
Second, run the global test on every possible combination of such intervals.
Third, reject that the function is increasing over interval $\mathcal{X}_j$ if increasingness is rejected at level $\alpha$ for every combination of intervals that includes $\mathcal{X}_j$.
As shown by \citet[\S9.2.1]{LehmannRomano2022text}, this controls the familywise error rate at level $\alpha$.
\backmatter
\bmhead{Supplementary information}
The supplementary appendix has extensions to IVQR and functional coefficient models as well as full proofs of all our theoretical results.
Also provided is R code implementing our new methods and replicating all simulation and empirical results.
\bmhead{Acknowledgments}
Thanks to the following for their helpful feedback: Alyssa Carlson, Zack Miller, Shawn Ni, journal reviewers and editors, and participants in the annual meetings of the Missouri Valley Economic Association (2021) and Midwest Econometrics Group (2023).
This work is based on a chapter of the first author's PhD dissertation.
\section*{Declarations}
\bmhead{Conflict of interest}
We (the authors) declare no conflict of interest.
\bmhead{Data availability}
All data is available publicly and loaded automatically through the provided R code.
\bmhead{Code availability}
All code is available on the second author's website.
\footnote{\url{https://kaplandm.github.io/}}
\begin{thebibliography}{38}
\ifx \bisbn \undefined \def \bisbn #1{ISBN #1}\fi
\ifx \binits \undefined \def \binits#1{#1}\fi
\ifx \bauthor \undefined \def \bauthor#1{#1}\fi
\ifx \batitle \undefined \def \batitle#1{#1}\fi
\ifx \bjtitle \undefined \def \bjtitle#1{#1}\fi
\ifx \bvolume \undefined \def \bvolume#1{\textbf{#1}}\fi
\ifx \byear \undefined \def \byear#1{#1}\fi
\ifx \bissue \undefined \def \bissue#1{#1}\fi
\ifx \bfpage \undefined \def \bfpage#1{#1}\fi
\ifx \blpage \undefined \def \blpage #1{#1}\fi
\ifx \burl \undefined \def \burl#1{\textsf{#1}}\fi
\ifx \doiurl \undefined \def \doiurl#1{\url{https://doi.org/#1}}\fi
\ifx \textit{et al.} \undefined \fi
\ifx \binstitute \undefined \def \binstitute#1{#1}\fi
\ifx \binstitutionaled \undefined \def \binstitutionaled#1{#1}\fi
\ifx \bctitle \undefined \def \bctitle#1{#1}\fi
\ifx \beditor \undefined \def \beditor#1{#1}\fi
\ifx \bpublisher \undefined \def \bpublisher#1{#1}\fi
\ifx \bbtitle \undefined \def \bbtitle#1{#1}\fi
\ifx \bedition \undefined \def \bedition#1{#1}\fi
\ifx \bseriesno \undefined \def \bseriesno#1{#1}\fi
\ifx \blocation \undefined \def \blocation#1{#1}\fi
\ifx \bsertitle \undefined \def \bsertitle#1{#1}\fi
\ifx \bsnm \undefined \def \bsnm#1{#1}\fi
\ifx \bsuffix \undefined \def \bsuffix#1{#1}\fi
\ifx \bparticle \undefined \def \bparticle#1{#1}\fi
\ifx \barticle \undefined \def \barticle#1{#1}\fi
\bibcommenthead
\ifx \bconfdate \undefined \def \bconfdate #1{#1}\fi
\ifx \botherref \undefined \def \botherref #1{#1}\fi
\ifx \url \undefined \def \url#1{\textsf{#1}}\fi
\ifx \bchapter \undefined \def \bchapter#1{#1}\fi
\ifx \bbook \undefined \def \bbook#1{#1}\fi
\ifx \bcomment \undefined \def \bcomment#1{#1}\fi
\ifx \oauthor \undefined \def \oauthor#1{#1}\fi
\ifx \citeauthoryear \undefined \def \citeauthoryear#1{#1}\fi
\ifx \undefined \fi
\ifx \bconflocation \undefined \def \bconflocation#1{#1}\fi
\ifx \arxivurl \undefined \def \arxivurl#1{\textsf{#1}}\fi
\csname PreBibitemsHook\endcsname
\bibitem[\protect\citeauthoryear{Andrews}{1991}]{Andrews1991HAC}
\begin{barticle}
\bauthor{\bsnm{Andrews}, \binits{D.W.K.}}:
\batitle{Heteroskedasticity and autocorrelation consistent covariance matrix estimation}.
\bjtitle{Econometrica}
\bvolume{59}(\bissue{3}),
\bfpage{817}--\blpage{858}
(\byear{1991})
\doiurl{10.2307/2938229}
\end{barticle}
\bibitem[\protect\citeauthoryear{Armstrong and Shen}{2023}]{ArmstrongShen2023}
\begin{barticle}
\bauthor{\bsnm{Armstrong}, \binits{T.B.}},
\bauthor{\bsnm{Shen}, \binits{S.}}:
\batitle{Inference on optimal treatment assignments}.
\bjtitle{The Japanese Economic Review}
\bvolume{74}(\bissue{4}),
\bfpage{471}--\blpage{500}
(\byear{2023})
\doiurl{10.1007/s42973-023-00138-1}
\end{barticle}
\bibitem[\protect\citeauthoryear{Benjamini and Hochberg}{1995}]{BenjaminiHochberg1995}
\begin{barticle}
\bauthor{\bsnm{Benjamini}, \binits{Y.}},
\bauthor{\bsnm{Hochberg}, \binits{Y.}}:
\batitle{Controlling the false discovery rate: a practical and powerful approach to multiple testing}.
\bjtitle{Journal of the Royal Statistical Society: Series B}
\bvolume{57}(\bissue{1}),
\bfpage{289}--\blpage{300}
(\byear{1995})
\doiurl{10.1111/j.2517-6161.1995.tb02031.x}
\end{barticle}
\bibitem[\protect\citeauthoryear{Blundell and Powell}{2003}]{BlundellPowell2003}
\begin{bchapter}
\bauthor{\bsnm{Blundell}, \binits{R.}},
\bauthor{\bsnm{Powell}, \binits{J.L.}}:
\bctitle{Endogeneity in nonparametric and semiparametric regression models}.
In: \beditor{\bsnm{Dewatripont}, \binits{M.}},
\beditor{\bsnm{Hansen}, \binits{L.P.}},
\beditor{\bsnm{Turnovsky}, \binits{S.J.}} (eds.)
\bbtitle{Advances in Economics and Econometrics: Theory and Applications, Eighth World Congress}.
\bsertitle{Econometric Society Monographs},
vol. \bseriesno{2},
pp. \bfpage{312}--\blpage{357}.
\bpublisher{Cambridge University Press},
\blocation{Cambridge}
(\byear{2003}).
\bcomment{Chap. 8}.
\doiurl{10.1017/CBO9780511610257}
\end{bchapter}
\bibitem[\protect\citeauthoryear{Chernozhukov et~al.}{2019}]{ChernozhukovEtAl2019many}
\begin{barticle}
\bauthor{\bsnm{Chernozhukov}, \binits{V.}},
\bauthor{\bsnm{Chetverikov}, \binits{D.}},
\bauthor{\bsnm{Kato}, \binits{K.}}:
\batitle{Inference on causal and structural parameters using many moment inequalities}.
\bjtitle{Review of Economic Studies}
\bvolume{86}(\bissue{5}),
\bfpage{1867}--\blpage{1900}
(\byear{2019})
\doiurl{10.1093/restud/rdy065}
\end{barticle}
\bibitem[\protect\citeauthoryear{Chen}{2007}]{Chen2007}
\begin{bchapter}
\bauthor{\bsnm{Chen}, \binits{X.}}:
\bctitle{Large sample sieve estimation of semi-nonparametric models}.
In: \beditor{\bsnm{Heckman}, \binits{J.J.}},
\beditor{\bsnm{Leamer}, \binits{E.E.}} (eds.)
\bbtitle{Handbook of Econometrics}
vol. \bseriesno{6B},
pp. \bfpage{5549}--\blpage{5632}.
\bpublisher{Elsevier},
\blocation{Amsterdam}
(\byear{2007}).
\bcomment{Chap. 76}.
\doiurl{10.1016/S1573-4412(07)06076-X}
\end{bchapter}
\bibitem[\protect\citeauthoryear{Chetverikov}{2019}]{Chetverikov2019}
\begin{barticle}
\bauthor{\bsnm{Chetverikov}, \binits{D.}}:
\batitle{Testing regression monotonicity in econometric models}.
\bjtitle{Econometric Theory}
\bvolume{35}(\bissue{4}),
\bfpage{729}--\blpage{776}
(\byear{2019})
\doiurl{10.1017/S0266466618000282}
\end{barticle}
\bibitem[\protect\citeauthoryear{Canay and Shaikh}{2017}]{CanayShaikh2017}
\begin{bchapter}
\bauthor{\bsnm{Canay}, \binits{I.A.}},
\bauthor{\bsnm{Shaikh}, \binits{A.M.}}:
\bctitle{Practical and theoretical advances in inference for partially identified models}.
In: \beditor{\bsnm{Pakes}, \binits{A.}},
\beditor{\bsnm{Honor{\'e}}, \binits{B.}},
\beditor{\bsnm{Samuelson}, \binits{L.}},
\beditor{\bsnm{Piazzesi}, \binits{M.}} (eds.)
\bbtitle{Advances in Economics and Econometrics: Eleventh World Congress}.
\bsertitle{Econometric Society Monographs},
vol. \bseriesno{2},
pp. \bfpage{271}--\blpage{306}.
\bpublisher{Cambridge University Press},
\blocation{Cambridge}
(\byear{2017}).
\bcomment{Chap. 9}.
\doiurl{10.1017/9781108227223.009} .
\burl{https://doi.org/10.1017/9781108227223.009}
\end{bchapter}
\bibitem[\protect\citeauthoryear{Davidson and Duclos}{2013}]{DavidsonDuclos2013}
\begin{barticle}
\bauthor{\bsnm{Davidson}, \binits{R.}},
\bauthor{\bsnm{Duclos}, \binits{J.-Y.}}:
\batitle{Testing for restricted stochastic dominance}.
\bjtitle{Econometric Reviews}
\bvolume{32}(\bissue{1}),
\bfpage{84}--\blpage{125}
(\byear{2013})
\end{barticle}
\bibitem[\protect\citeauthoryear{Dewan and Neligh}{2020}]{DewanNeligh2020}
\begin{barticle}
\bauthor{\bsnm{Dewan}, \binits{A.}},
\bauthor{\bsnm{Neligh}, \binits{N.}}:
\batitle{Estimating information cost functions in models of rational inattention}.
\bjtitle{Journal of Economic Theory}
\bvolume{187},
\bfpage{105011}
(\byear{2020})
\end{barticle}
\bibitem[\protect\citeauthoryear{Dean and Neligh}{2023}]{DeanNeligh2023}
\begin{barticle}
\bauthor{\bsnm{Dean}, \binits{M.}},
\bauthor{\bsnm{Neligh}, \binits{N.}}:
\batitle{Experimental tests of rational inattention}.
\bjtitle{Journal of Political Economy}
\bvolume{131}(\bissue{12}),
\bfpage{3415}--\blpage{3461}
(\byear{2023})
\end{barticle}
\bibitem[\protect\citeauthoryear{Friedrich et~al.}{2020}]{FriedrichEtAl2020}
\begin{barticle}
\bauthor{\bsnm{Friedrich}, \binits{M.}},
\bauthor{\bsnm{Beutner}, \binits{E.}},
\bauthor{\bsnm{Reuvers}, \binits{H.}},
\bauthor{\bsnm{Smeekes}, \binits{S.}},
\bauthor{\bsnm{Urbain}, \binits{J.-P.}},
\bauthor{\bsnm{Bader}, \binits{W.}},
\bauthor{\bsnm{Franco}, \binits{B.}},
\bauthor{\bsnm{Lejeune}, \binits{B.}},
\bauthor{\bsnm{Mahieu}, \binits{E.}}:
\batitle{A statistical analysis of time trends in atmospheric ethane}.
\bjtitle{Climatic Change}
\bvolume{162},
\bfpage{105}--\blpage{125}
(\byear{2020})
\end{barticle}
\bibitem[\protect\citeauthoryear{Goldman and Kaplan}{2018}]{GoldmanKaplan2018c}
\begin{barticle}
\bauthor{\bsnm{Goldman}, \binits{M.}},
\bauthor{\bsnm{Kaplan}, \binits{D.M.}}:
\batitle{Comparing distributions by multiple testing across quantiles or {CDF} values}.
\bjtitle{Journal of Econometrics}
\bvolume{206}(\bissue{1}),
\bfpage{143}--\blpage{166}
(\byear{2018})
\end{barticle}
\bibitem[\protect\citeauthoryear{Ghosal et~al.}{2000}]{GhosalEtAl2000}
\begin{barticle}
\bauthor{\bsnm{Ghosal}, \binits{S.}},
\bauthor{\bsnm{Sen}, \binits{A.}},
\bauthor{\bsnm{Vaart}, \binits{A.W.}}:
\batitle{Testing monotonicity of regression}.
\bjtitle{Annals of Statistics}
\bvolume{28}(\bissue{4}),
\bfpage{1054}--\blpage{1082}
(\byear{2000})
\doiurl{10.1214/aos/1015956707}
\end{barticle}
\bibitem[\protect\citeauthoryear{Hall and Heckman}{2000}]{HallHeckman2000}
\begin{barticle}
\bauthor{\bsnm{Hall}, \binits{P.}},
\bauthor{\bsnm{Heckman}, \binits{N.E.}}:
\batitle{Testing for monotonicity of a regression mean by calibrating for linear functions}.
\bjtitle{Annals of Statistics}
\bvolume{28}(\bissue{1}),
\bfpage{20}--\blpage{39}
(\byear{2000})
\doiurl{10.1214/aos/1016120363}
\end{barticle}
\bibitem[\protect\citeauthoryear{Holm}{1979}]{Holm1979}
\begin{barticle}
\bauthor{\bsnm{Holm}, \binits{S.}}:
\batitle{A simple sequentially rejective multiple test procedure}.
\bjtitle{Scandinavian Journal of Statistics}
\bvolume{6}(\bissue{2}),
\bfpage{65}--\blpage{70}
(\byear{1979})
\end{barticle}
\bibitem[\protect\citeauthoryear{Ibragimov and M{\"u}ller}{2010}]{IbragimovMueller2010}
\begin{barticle}
\bauthor{\bsnm{Ibragimov}, \binits{R.}},
\bauthor{\bsnm{M{\"u}ller}, \binits{U.K.}}:
\batitle{$t$-statistic based correlation and heterogeneity robust inference}.
\bjtitle{Journal of Business \& Economic Statistics}
\bvolume{28}(\bissue{4}),
\bfpage{453}--\blpage{468}
(\byear{2010})
\end{barticle}
\bibitem[\protect\citeauthoryear{Imbens and Newey}{2009}]{ImbensNewey2009}
\begin{barticle}
\bauthor{\bsnm{Imbens}, \binits{G.W.}},
\bauthor{\bsnm{Newey}, \binits{W.K.}}:
\batitle{Identification and estimation of triangular simultaneous equations models without additivity}.
\bjtitle{Econometrica}
\bvolume{77}(\bissue{5}),
\bfpage{1481}--\blpage{1512}
(\byear{2009})
\doiurl{10.3982/ECTA7108}
\end{barticle}
\bibitem[\protect\citeauthoryear{Kaplan}{2024}]{Kaplan2024}
\begin{barticle}
\bauthor{\bsnm{Kaplan}, \binits{D.M.}}:
\batitle{Inference on consensus ranking of distributions}.
\bjtitle{Journal of Business \& Economic Statistics}
\bvolume{42}(\bissue{3}),
\bfpage{839}--\blpage{850}
(\byear{2024})
\doiurl{10.1080/07350015.2023.2252040}
\end{barticle}
\bibitem[\protect\citeauthoryear{Kostyshak and Luo}{2021}]{KostyshakLuo2021}
\begin{botherref}
\oauthor{\bsnm{Kostyshak}, \binits{S.}},
\oauthor{\bsnm{Luo}, \binits{Y.}}:
The Partial Monotonicity Parameter: A Generalization of Monotonicity.
Working paper available at \url{https://people.clas.ufl.edu/skostyshak/files/pmp.pdf}
(2021)
\end{botherref}
\bibitem[\protect\citeauthoryear{Lewbel et~al.}{2012}]{LewbelEtAl2012}
\begin{barticle}
\bauthor{\bsnm{Lewbel}, \binits{A.}},
\bauthor{\bsnm{Dong}, \binits{Y.}},
\bauthor{\bsnm{Yang}, \binits{T.T.}}:
\batitle{Comparing features of convenient estimators for binary choice models with endogenous regressors}.
\bjtitle{The Canadian Journal of Economics / Revue canadienne d'Economique}
\bvolume{45}(\bissue{3}),
\bfpage{809}--\blpage{829}
(\byear{2012})
\doiurl{10.1111/j.1540-5982.2012.01733.x}
\end{barticle}
\bibitem[\protect\citeauthoryear{Lehmann and Romano}{2005}]{LehmannRomano2005fwer}
\begin{barticle}
\bauthor{\bsnm{Lehmann}, \binits{E.L.}},
\bauthor{\bsnm{Romano}, \binits{J.P.}}:
\batitle{Generalizations of the familywise error rate}.
\bjtitle{Annals of Statistics}
\bvolume{33}(\bissue{3}),
\bfpage{1138}--\blpage{1154}
(\byear{2005})
\doiurl{10.1214/009053605000000084}
\end{barticle}
\bibitem[\protect\citeauthoryear{Lehmann and Romano}{2022}]{LehmannRomano2022text}
\begin{bbook}
\bauthor{\bsnm{Lehmann}, \binits{E.L.}},
\bauthor{\bsnm{Romano}, \binits{J.P.}}:
\bbtitle{Testing Statistical Hypotheses},
\bedition{4th} edn.
\bsertitle{Springer Texts in Statistics}.
\bpublisher{Springer},
\blocation{Cham, Switzerland}
(\byear{2022}).
\doiurl{10.1007/978-3-030-70578-7}
\end{bbook}
\bibitem[\protect\citeauthoryear{Manski and Pepper}{2000}]{ManskiPepper2000}
\begin{barticle}
\bauthor{\bsnm{Manski}, \binits{C.F.}},
\bauthor{\bsnm{Pepper}, \binits{J.V.}}:
\batitle{Monotone instrumental variables: With an application to the returns to schooling}.
\bjtitle{Econometrica}
\bvolume{68}(\bissue{4}),
\bfpage{997}--\blpage{1010}
(\byear{2000})
\doiurl{10.1111/1468-0262.t01-1-00144a}
\end{barticle}
\bibitem[\protect\citeauthoryear{Nadarajah and Kotz}{2008}]{NadarajahKotz2008}
\begin{barticle}
\bauthor{\bsnm{Nadarajah}, \binits{S.}},
\bauthor{\bsnm{Kotz}, \binits{S.}}:
\batitle{Exact distribution of the max/min of two {Gaussian} random variables}.
\bjtitle{IEEE Transactions on Very Large Scale Integration (VLSI) Systems}
\bvolume{16}(\bissue{2}),
\bfpage{210}--\blpage{212}
(\byear{2008})
\doiurl{10.1109/TVLSI.2007.912191}
\end{barticle}
\bibitem[\protect\citeauthoryear{Newey and West}{1987}]{NeweyWest1987}
\begin{barticle}
\bauthor{\bsnm{Newey}, \binits{W.K.}},
\bauthor{\bsnm{West}, \binits{K.D.}}:
\batitle{A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix}.
\bjtitle{Econometrica}
\bvolume{55}(\bissue{3}),
\bfpage{703}--\blpage{708}
(\byear{1987})
\end{barticle}
\bibitem[\protect\citeauthoryear{{R Core Team}}{2024}]{R.core}
\begin{bbook}
\bauthor{\bsnm{{R Core Team}}}:
\bbtitle{R: A Language and Environment for Statistical Computing}.
\bpublisher{R Foundation for Statistical Computing},
\blocation{Vienna, Austria}
(\byear{2024}).
\bcomment{R Foundation for Statistical Computing}.
\burl{https://www.R-project.org}
\end{bbook}
\bibitem[\protect\citeauthoryear{Ridley et~al.}{2020}]{RidleyEtAl2020}
\begin{barticle}
\bauthor{\bsnm{Ridley}, \binits{M.}},
\bauthor{\bsnm{Rao}, \binits{G.}},
\bauthor{\bsnm{Schilbach}, \binits{F.}},
\bauthor{\bsnm{Patel}, \binits{V.}}:
\batitle{Poverty, depression, and anxiety: Causal evidence and mechanisms}.
\bjtitle{Science}
\bvolume{370}(\bissue{6522}),
\bfpage{0214}
(\byear{2020})
\end{barticle}
\bibitem[\protect\citeauthoryear{Romano and Shaikh}{2010}]{RomanoShaikh2010}
\begin{barticle}
\bauthor{\bsnm{Romano}, \binits{J.P.}},
\bauthor{\bsnm{Shaikh}, \binits{A.M.}}:
\batitle{Inference for the identified set in partially identified econometric models}.
\bjtitle{Econometrica}
\bvolume{78}(\bissue{1}),
\bfpage{169}--\blpage{211}
(\byear{2010})
\doiurl{10.3982/ECTA6706}
\end{barticle}
\bibitem[\protect\citeauthoryear{Romano et~al.}{2014}]{RomanoShaikhWolf2014}
\begin{barticle}
\bauthor{\bsnm{Romano}, \binits{J.P.}},
\bauthor{\bsnm{Shaikh}, \binits{A.M.}},
\bauthor{\bsnm{Wolf}, \binits{M.}}:
\batitle{A practical two-step method for testing moment inequalities}.
\bjtitle{Econometrica}
\bvolume{82}(\bissue{5}),
\bfpage{1979}--\blpage{2002}
(\byear{2014})
\doiurl{10.3982/ECTA11011}
\end{barticle}
\bibitem[\protect\citeauthoryear{Romano and Wolf}{2013}]{RomanoWolf2013}
\begin{barticle}
\bauthor{\bsnm{Romano}, \binits{J.P.}},
\bauthor{\bsnm{Wolf}, \binits{M.}}:
\batitle{Testing for monotonicity in expected asset returns}.
\bjtitle{Journal of Empirical Finance}
\bvolume{23},
\bfpage{93}--\blpage{116}
(\byear{2013})
\doiurl{10.1016/j.jempfin.2013.05.001}
\end{barticle}
\bibitem[\protect\citeauthoryear{Shea}{2018}]{R.wooldridge}
\begin{botherref}
\oauthor{\bsnm{Shea}, \binits{J.M.}}:
Wooldridge: 111 Data Sets from ``Introductory Econometrics: A Modern Approach, 6e'' by Jeffrey M.\ Wooldridge.
(2018).
R package version 1.3.1.
\url{https://CRAN.R-project.org/package=wooldridge}
\end{botherref}
\bibitem[\protect\citeauthoryear{van~der Vaart}{1998}]{vanderVaart1998}
\begin{bbook}
\bauthor{\bsnm{Vaart}, \binits{A.W.}}:
\bbtitle{Asymptotic Statistics}.
\bpublisher{Cambridge University Press},
\blocation{Cambridge}
(\byear{1998}).
\doiurl{10.1017/CBO9780511802256}
\end{bbook}
\bibitem[\protect\citeauthoryear{Venables and Ripley}{2002}]{R.MASS}
\begin{bbook}
\bauthor{\bsnm{Venables}, \binits{W.N.}},
\bauthor{\bsnm{Ripley}, \binits{B.D.}}:
\bbtitle{Modern Applied Statistics with {S}},
\bedition{4}th edn.
\bpublisher{Springer},
\blocation{New York}
(\byear{2002}).
\burl{https://www.stats.ox.ac.uk/pub/MASS4/}
\end{bbook}
\bibitem[\protect\citeauthoryear{Wu and Kaplan}{2025}]{WuKaplan2025a}
\begin{botherref}
\oauthor{\bsnm{Wu}, \binits{Q.}},
\oauthor{\bsnm{Kaplan}, \binits{D.M.}}:
Multiple Testing of Stochastic Monotonicity.
Working paper available at \url{https://kaplandm.github.io/}
(2025)
\end{botherref}
\bibitem[\protect\citeauthoryear{Wooldridge}{2010}]{Wooldridge2010}
\begin{bbook}
\bauthor{\bsnm{Wooldridge}, \binits{J.M.}}:
\bbtitle{Econometric Analysis of Cross Section and Panel Data},
\bedition{2nd} edn.
\bpublisher{MIT Press},
\blocation{Cambridge, MA}
(\byear{2010}).
\burl{https://www.worldcat.org/oclc/831625495}
\end{bbook}
\bibitem[\protect\citeauthoryear{Zeileis}{2004}]{R.sandwich}
\begin{barticle}
\bauthor{\bsnm{Zeileis}, \binits{A.}}:
\batitle{Econometric computing with {HC} and {HAC} covariance matrix estimators}.
\bjtitle{Journal of Statistical Software}
\bvolume{11}(\bissue{10}),
\bfpage{1}--\blpage{17}
(\byear{2004})
\end{barticle}
\bibitem[\protect\citeauthoryear{Zhang et~al.}{2015}]{ZhangEtAl2015eco}
\begin{barticle}
\bauthor{\bsnm{Zhang}, \binits{Z.}},
\bauthor{\bsnm{Yan}, \binits{C.}},
\bauthor{\bsnm{Krebs}, \binits{C.J.}},
\bauthor{\bsnm{Stenseth}, \binits{N.C.}}:
\batitle{Ecological non-monotonicity and its effects on complexity and stability of populations, communities and ecosystems}.
\bjtitle{Ecological Modelling}
\bvolume{312},
\bfpage{374}--\blpage{384}
(\byear{2015})
\end{barticle}
\end{thebibliography}
\clearpage