\subsection*{#1}\hskip3mm}
\def\vskip3mm\par{\vskip3mm\par}
\defdf{df}
\defadf{adf}
\defdf_{{\rm fit}}{df_{{\rm fit}}}
\defdf_{{\rm res}}{df_{{\rm res}}}
\defdf_{\yhat}{df_{{\widehat y}}}
\def\sigma_b{\sigma_b}
\def\sigma_u{\sigma_u}
\def{\widehat\sigma}_{\varepsilon}{{\widehat\sigma}_{\varepsilon}}
\def{\widehat\sigma}_u{{\widehat\sigma}_u}
\def{\widehat\sigma}_{\varepsilon}{{\widehat\sigma}_{\varepsilon}}
\def{\widehat\sigma}_b{{\widehat\sigma}_b}
\def{\widehat\sigma}{{\widehat\sigma}}
\def\sigma^2_b{\sigma^2_b}
\def\sigma^2_{\varepsilon}{\sigma^2_{\varepsilon}}
\def{\widehat\sigma}^2_{\varepsilon,0}{{\widehat\sigma}^2_{\varepsilon,0}}
\def{\widehat\sigma}^2_{\varepsilon}{{\widehat\sigma}^2_{\varepsilon}}
\def{\widehat\sigma}^2_b{{\widehat\sigma}^2_b}
\def{\rm df(II}|{\rm I)}{{\rm df(II}|{\rm I)}}
\def{\widehat m}_{\lambda}{{\widehat m}_{\lambda}}
\def\Dsc{{\cal D}}
\defA_{\epsilon}{A_{\epsilon}}
\defB_{\epsilon}{B_{\epsilon}}
\defA_b{A_b}
\defB_b{B_b}
\def\tilde{\bX}_r{\tilde{{\bf X}}_r}
\def\mbox{\tt main}{\mbox{\tt main}}
\def\mathop{\rm arg\min}{\bbeta,\bb}{\mathop{\rm arg\min}{\boldsymbol{\beta},{\bf b}}}
\def\Bsc{{\cal B}}
\def\mbox{\tt log(amt)}{\mbox{\tt log(amt)}}

\def\mathbb{A}{\mathbb{A}}
\def\mathbb{B}{\mathbb{B}}
\def\mathbb{Z}{\mathbb{Z}}
\def\mathbb{W}{\mathbb{W}}

\def\widehat{\mathbb{A}}{\widehat{\mathbb{A}}}


\def\phantom{-}{\phantom{-}}
\def\phantom{1}{\phantom{1}}


\section{Introduction}
When evaluating the effectiveness of a treatment, policy, or intervention, the desired measure of effectiveness may be expensive to collect, not routinely available, or may take a long time to occur. In these cases, it is sometimes possible to identify a surrogate outcome that can more easily/quickly/cheaply capture the effect of interest. For example, when evaluating an intervention designed to delay dementia onset, the time required to observe enough dementia diagnoses is often very long and surrogates that have been considered for intervention evaluation include mild cognitive impairment, adiponectin levels, neuroimages, amyloid plaques and neurofibrillary tangles \cite{small2006diagnostic,kantarci2004quantitative,teixeira2013decreased,levey2006mild,mathis2007impact}. Similarly, in studies evaluating treatments to prevent diabetes, surrogate measures for diabetes onset have included changes in body weight, fasting plasma glucose and hemoglobin A1c \cite{caveney2011diabetes,diabetes200910,choi2011hemoglobin,Parast2017}. 

Theory and methods for evaluating the strength of surrogate markers have been well studied in the context of a single surrogate marker measured in the course of a randomized clinical study. In particular, robust nonparametric methods have been proposed in \cite{parast2019assessing, parast2016robust}. However, fully nonparametric methods are not available or not reliable when the numbers of markers is more than one or two. In these cases, an initial model may be used to reduce the dimensionality of the surrogate \cite{agniel2020evaluation,parastevaluating,parast2016robust}, but an inappropriate initial model can produce badly biased estimates of the utility of the surrogate. Additionally, as the dimension of the surrogate nears or exceeds the sample size, initial parametric models may not have adequate properties or might not even be possible. Furthermore, assessments of surrogate markers are not solely limited to studies of randomized treatment \cite{de2004measurement,obirikorang2012total}. In fact, using surrogate outcomes based on complex observational data is finding purchase in all corners of science \cite{wang2020methods}.

In this work, we propose an efficient nonparametric method for evaluating high-dimensional surrogate markers in studies where the treatment need not be randomized. Our approach draws on a connection between quantifying the utility of a surrogate marker and the most fundamental tools of causal inference -- namely, methods for estimating the average treatment effect. We take advantage of state-of-the-art methods for incorporating flexible machine learning and/or sparse high-dimensional models into estimation of treatment effects via the cross-fitting approach of \cite{Chernozhukov2017}. Cross-fitting allows us to derive asymptotic distributions and compute confidence intervals for key quantities while putting minimal restrictions on what types of estimators are used. Furthermore, our estimator and results may be used regardless of the dimension of the surrogates.

The structure of the rest of the paper is as follows. In Section 2, we lay out notation and the setting in which we are working, and we motivate the importance of evaluating surrogate markers in light of recent advances in estimating surrogates from high-dimensional data. In Section 3, we detail assumptions necessary for identifying and interpreting parameters of interest, as well as our approach to estimating the quality of a surrogate marker. We derive the asymptotic behavior of our estimator in Section 4 and discuss variance estimation and inference. We evaluate the performance of the proposed approach using a simulation study in Section 
. We give final remarks and draw connections to other methods for evaluating surrogate markers in Section 6.

\section{Evaluation of Surrogate Markers}

\subsection{Notation}
Let $A$ denote a binary treatment, and let the primary outcome of the study be $Y$. Let there be a vector ${\bf S}$ of potential surrogate information, and let ${\bf X}$ be a vector of pre-treatment covariates. The primary quantity of interest is the treatment effect on the outcome 
\[
\Delta = \mathbb{E}\{ Y^{(1)} - Y^{(0)} \}
\]
where $Y^{(a)}$ is the potential or counterfactual outcome that would have been observed if the observed value of treatment were $A = a$, possibly contrary to fact. Similarly let ${\bf S}^{(a)}$ be the potential/counterfactual value the vector of surrogates would take if $A = a$. We desire to approximate $\Delta$ via a difference in some function of the surrogates ${\bf S}$ and covariates ${\bf X}$. Let $\psi_a({\bf X}, {\bf S})$ be such a scalar function that is here indexed by treatment $a$ but need not depend on treatment. Let the data observed in the current study be ${\cal O} = ({\bf X}_i, A_i, {\bf S}_i, Y_i)_{i=1,...,n}$, $n$ iid realizations of ${\bf O} = ({\bf X}, A, {\bf S}, Y)$. And let ${\cal F} = (\check{\bX}_i, \check{A}_i, \check{\bS}_i)_{i=1,...,n^*}$ be $n^*$ iid realization of ${\bf O}$ in a future study. The function $\psi_a({\bf x}, {\bf s})$ is useful and is of particular interest because for a future study, one could then use $\psi_i = \psi_1(\check{\bX}_i, \check{\bS}_i)I\{\check{A}_i = 1\} + \psi_0(\check{\bX}_i, \check{\bS}_i)I\{\check{A}_i = 0\}$ as the outcome instead of $Y$. The true treatment effect $\Delta$ could be approximated by estimating the effect of $\check{A}_i$ on $\psi_i$ conditional on $\check{\bX}_i$ via regression, weighting, matching, or any other reasonable approach.

\subsection{Importance of quantifying surrogate strength}

Understanding the strength of a potential surrogate is important for many reasons. In practice, information about surrogate strength will inform decisions regarding whether to measure the surrogate in a future study (especially if it is costly or invasive to measure) and/or whether to use the surrogate to assess treatment effectiveness in a future study. In addition, a number of novel statistical methods for using surrogates in future studies can only be applied when the surrogate is strong. For example, Parast {\em et al.} \cite{parast2019using} propose a robust nonparametric procedure to test for a treatment effect using surrogate marker information measured prior to the end of the study in a time‐to‐event outcome setting, but they rely on the assumption that the surrogate is sufficiently strong. As another example, Price {\em et al.} \cite{price2018estimation} propose constructing a function $\psi_a(\cdot)$ of the surrogates that leads to some optimality  properties, but it can be shown that this method should only be used with a strong surrogate. 

Specifically, they construct an optimal $\psi_a(\cdot)$ in terms of minimizing the following mean squared error (MSE)
\begin{equation}
{\cal M}_\psi = \mathbb{E}\left[e({\bf X}) \left\{Y^{(1)} - \psi_1({\bf X}, {\bf S}^{(1)})\right\}^2 \right] + \mathbb{E}\left[\{1-e({\bf X})\} \left\{Y^{(0)} - \psi_0({\bf X}, {\bf S}^{(0)})\right]^2 \right], \qquad e({\bf x}) = P(A = 1 | {\bf X})\label{MSEPSI}
\end{equation}
under the constraint that it satisfies the so-called Prentice definition, i.e., 
\begin{align}\label{prentice}
    E(Y^{(1)} - Y^{(0)}) = 0 \text{ if and only if } E\left\{\psi_1({\bf X}, {\bf S}^{(1)}) - \psi_0({\bf X}, {\bf S}^{(0)})\right\} = 0.
\end{align} 
The optimal transformations in this case are shown to be
$
\psi_a({\bf x}, {\bf s}) = \mathbb{E}\{Y^{(a)} | {\bf X} = {\bf x}, {\bf S}^{(a)} = {\bf s}\}.
$
In addition to this appealing optimality property, this proposal also resolves the so-called ``surrogate paradox." The paradox states that the treatment could have a positive causal effect on a univariate surrogate $S$ which could have a positive correlation with the outcome, but the treatment could yet have a negative effect on the outcome. The Price surrogate resolves the surrogate paradox by definition. The treatment effect on the surrogate $\mathbb{E}(\check{\Delta})$ is by definition equivalent to the treatment effect on the outcome: 

    \mathbb{E}(\check{\Delta}) = \mathbb{E}\{\psi_1({\bf X}, {\bf S}^{(1)})\} - \mathbb{E}\{\psi_0({\bf X}, {\bf S}^{(0)})\} = \Delta.\label{psi-mu} 


However, the Price surrogate's resolution of the surrogate paradox is in some sense too good. The surrogate paradox is thus resolved for every potential set of surrogates, good and bad. Even if the surrogates ${\bf S}$ are completely unrelated to the outcome $Y$, $\mathbb{E}(\check{\Delta}) = \Delta$. In fact, the power to detect a treatment effect in a future study using $\check{\Delta}$ actually increases as the surrogate becomes weaker (explaining less of the treatment effect), if the Prentice definition does not hold. Consider the toy example depicted in Figure \ref{toy-eg} (see Appendix \ref{fig1-sim} for details of the data generation). 
\begin{figure}
\centering
\begin{tabular}{c}
\includegraphics[width =\textwidth]{toy-example_price.jpg}
\end{tabular}\vspace{0.2in}
\caption{Distribution of outcomes $Y$ in an initial study (Study 1) and surrogate functions (approximating the outcome) $\psi_1$ (surrogate in treatment group) and $\psi_0$ (surrogate in control group) in a future study (Study 2) based on simulated data. Vertical dotted lines correspond to the means of the treatment groups in the first study. See Appendix \ref{fig1-sim} for simulation details. Surrogate functions $\psi_0, \psi_1$ were estimated via a correctly specified linear regression in the first study and then applied to the second study population. When the surrogate was strong, the distribution of the surrogate in the second study was quite close to the distribution of the outcome in the first study. However, when the surrogate was weak, the distribution of the surrogate in the second study clustered around the group mean from the first study.}
\label{toy-eg}
\end{figure}
As the strength of the surrogate increases, $\psi_a(\cdot)$ becomes more like the true outcome $Y_i^{(a)}$. As the surrogate becomes weaker, $\psi_a(\cdot)$ becomes more like $\bar{ Y}_a$ the mean outcome in treatment group $A = a$ from the first study. This means that a weak surrogate ensures that the treatment effect in the second study will be identical to the treatment effect in the first study because the distribution of $\psi_a(\check{\bX}_i, \check{\bS}_i)$ will cluster closely around the estimate of $E(Y^{(a)})$ from the first study (see specifically, the fourth panel of Figure \ref{fig1-sim}). In this scenario, the second study is providing no new information about the treatment. Thus, the Price surrogate should only be used with a strong ${\bf S}$, and its use with a weak or possibly even moderately strong set of surrogates may not be advisable.

Therefore, for both decision-making purposes and to use statistical methods that take advantage of strong surrogates, methods are needed to rigorously quantify the strength of the proposed set of surrogates before using (functions of) the surrogates in practice. In the following section, we propose a model-free method for estimating the strength of a possibly high-dimensional set of surrogates in observational studies.














\subsection{Proposed method to evaluate surrogate strength}\label{sec:pte}
In this section we propose a general method for evaluating the usefulness of a set of surrogates. This can be used to evaluate the suitability of the Price surrogate but can also be used with other surrogate transformations, and it is appropriate for ${\bf S}$ of any dimension. To evaluate the surrogates' usefulness, we follow \cite{parast2019assessing, parast2016robust, agniel2020evaluation} in evaluating the \textit{residual treatment effect}
\begin{align}
    \Delta_{{\bf S}} = \mathbb{E}\left\{\psi_1({\bf X}, {\bf S}) - \psi_0({\bf X}, {\bf S})\right\},  \label{def:delta-s}
\end{align}
where $\psi_a({\bf x}, {\bf s}) = E(Y^{(a)} | {\bf X} = {\bf x}, {\bf S}^{(a)} = {\bf s})$, as defined above. This quantity corresponds to the treatment effect conditional on the distribution of ${\bf S}^{(0)}$ and ${\bf S}^{(1)}$ both being equal to the distribution of ${\bf S}$ (all conditional on ${\bf X}$). If ${\bf S}^{(0)}, {\bf S}^{(1)}$ are independent of $Y^{(0)}, Y^{(1)}$ conditional on ${\bf X}$, then $\Delta_{{\bf S}} = \Delta$. In contrast, if all of the treatment effect can be attributed to ${\bf S}$, then $\Delta_{{\bf S}}$ should be close to 0. We may use the residual treatment effect to quantify \textit{how much} of the treatment effect the surrogates account for via the \textit{proportion of treatment effect explained} (PTE)
\begin{align}
    R_{{\bf S}} = \frac{\Delta - \Delta_{{\bf S}}}{\Delta} = 1 - \frac{\Delta_{{\bf S}}}{\Delta}.
\end{align}

The quantity $\Delta_{{\bf S}}$ is defined such that it aims to quantify the effect of the treatment on the outcome that is \textit{not} captured by the effect of the treatment on the surrogate, by forcing the conditional distribution of the potential surrogate in both treatment groups to be the same. Based on this definition, $\Delta - \Delta_{{\bf S}}$ then parallels the natural indirect effect quantity \cite{joffe2009related,vanderweele2013surrogate} which is often used to assess the proportion mediated in a mediation analysis. However, there is an important distinction between assessing a variable (or variables) for surrogacy versus mediation. As described by \cite{vanderweele2013surrogate}, the aim of identifying a mediator is determining whether the effect of treatment operates through the mediator itself e.g., through some biological pathway. Often, a good surrogate marker is similarly conceptualized as a variable through which the treatment operates, but this is not necessarily required; a variable can be a good surrogate if it captures the treatment effect on the outcome, even if the treatment effect does not operate through the variable itself. 

\begin{remark} While $R_S$ is easy to understand and interpret, and it has analogues in the mediation literature, it is not guaranteed to lie between 0 and 1. We give conditions in Section \ref{pte-assump} that guarantee that $R_S \in (0,1)$. In situations when these conditions are not or are unlikely to be met, one may desire a secondary definition of the PTE that can be estimated in all situations similar to \cite{wang2020model}. Conceptually, here we would consider a good surrogate to be one whereby the surrogate transformation substantially improves the MSE for predicting $Y^{(a)}$ over a \textit{null} transformation that does not use ${\bf S}$. The MSE of the null transformation would be defined as:
\begin{align*}
    {\cal M}_0 &= \mathbb{E}\left(e({\bf X}) \left[Y^{(1)} - m_1\{{\bf X}\}\right]^2 \right) + \mathbb{E}\left(\{1-e({\bf X})\} \left[Y^{(0)} - m_0\{{\bf X}\}\right]^2 \right
\end{align*}
where $m_a({\bf x}) = \mathbb{E}(Y^{(a)} |{\bf X} = {\bf x}).$
The secondary MSE-based definition of PTE could then be defined as the proportion of the null MSE that is reduced by controlling for the surrogates:
\begin{align}
    R_M = \frac{{\cal M}_0 - {\cal M}_\psi}{{\cal M}_0} = 1 - {\cal M}_\psi/{\cal M}_0
\end{align}
where ${\cal M}_\psi$ is defined in \eqref{MSEPSI}.
\end{remark}

\section{Assumptions and Estimation}

\subsection{Identifying assumptions}\label{assump-ident}
We make the three typical assumptions of treatment effect estimation: consistency, positivity, and no unmeasured confounding. Specifically, we assume that the observed values of ${\bf S}$ and $Y$ when $A = a$ are identical to the counterfactuals ${\bf S}^{(a)}$ and $Y^{(a)}$ such that
${\bf S} = {\bf S}^{(1)} A + {\bf S}^{(0)}(1 - A)$ and $Y = Y^{(1)} A + Y^{(0)}(1 - A).$
We furthermore assume that ${\bf X}$ contains all confounders of the effects of $A$ on the surrogates and the outcome, such that the treatment $A$ is as good as randomized conditional on the covariates ${\bf X}$
\begin{align}
\left.\{{\bf S}^{(0)}, {\bf S}^{(1)}, Y^{(0)}, Y^{(1)}\} \independent A \right| {\bf X}.\label{psr}
\end{align}
In addition, we assume two forms of positivity, which ensures that individuals in the two study arms are not too different from one another. First, the usual positive probability of receiving either treatment for some $\epsilon_1 > 0$
\[
P\{\epsilon_1 < e({\bf X}) < 1-\epsilon_1\} = 1,
\]
and a related assumption that further conditions on the surrogates
\[
P\{\epsilon_2 < \pi({\bf X}, {\bf S}) < 1-\epsilon_2\} = 1,
\]
for some $\epsilon_2 > 0$ where $\pi({\bf S}, {\bf X}) = P(A = 1 | {\bf S}, {\bf X})$ is the \textit{surrogate score} -- akin to the propensity score \cite{rosenbaum1983central} -- which has been proposed for a different purpose in \cite{Athey2016}. Notably, because $\pi({\bf x}, {\bf s}) = \frac{P({\bf S} = {\bf s} | {\bf X} = {\bf x}, A = 1)e({\bf x})}{P({\bf S} = {\bf s} | {\bf X} = {\bf x})} = \frac{P({\bf S}^{(1)} = {\bf s} | {\bf X} = {\bf x})e({\bf x})}{P({\bf S}^{(1)} = {\bf s} | {\bf X} = {\bf x})e({\bf x}) + P({\bf S}^{(0)} = {\bf s} | {\bf X} = {\bf x})\{1-e({\bf x})\}}$, these two positivity conditions ensure that the conditional distribution of the counterfactual surrogates under treatment and control cannot be too different from one another, i.e., ensuring overlap.

We additionally give regularity conditions facilitating estimation in Appendix \ref{assump-reg}.

\subsection{Required assumptions to ensure interpretation of $R_{\bf S}$ as a proportion}\label{pte-assump}

The interpretation of $R_{\bf S}$ as the PTE depends on it actually being a proportion, lying between 0 and 1. We adapt conditions in \cite{wang2002measure} and \cite{agniel2020evaluation} which ensure that $0 \leq R_{\bf S} \leq 1$. These conditions are as follows (assuming without loss of generality that $\Delta > 0$):
\begin{align}
\int \psi_1({\bf x}, {\bf s})d\left\{F_{{\bf X}, {\bf S}^{(1)}}({\bf x}, {\bf s}) - F_{{\bf X}, {\bf S}}({\bf x}, {\bf s})\right\} \geq \int \psi_0({\bf x}, {\bf s})d\left\{F_{{\bf X}, {\bf S}^{(0)}}({\bf x}, {\bf s}) - F_{{\bf X}, {\bf S}}({\bf x}, {\bf s})\right\} \text{ and }\label{as1}
\end{align}
\begin{align}
\int \left\{\psi_1({\bf x}, {\bf s}) - \psi_0({\bf x}, {\bf s})\right\}dF_{{\bf X}, {\bf S}}({\bf x}, {\bf s}) \geq 0.\label{as2}
\end{align}
If $\Delta < 0$, then one would require \eqref{as1} and \eqref{as2} to hold with inequalities reversed. 

Because $F_{{\bf X}, {\bf S}}$ is a mixture distribution $F_{{\bf X}, {\bf S}}({\bf x}, {\bf s}) = e({\bf x}) F_{{\bf X}, {\bf S}^{(1)}}({\bf x}, {\bf s}) + \{1-e({\bf x})\}F_{{\bf X}, {\bf S}^{(0)}}({\bf x}, {\bf s})$, the first condition \eqref{as1} may be rewritten:
\begin{align}
\int \left\{e({\bf x})\psi_1({\bf x}, {\bf s}) + \{1-e({\bf x})\}\psi_0({\bf x}, {\bf s})\right\}d\left\{F_{{\bf X}, {\bf S}^{(1)}}({\bf x}, {\bf s}) - F_{{\bf X}, {\bf S}^{(0)}}({\bf x}, {\bf s})\right\} \geq 0.\label{as1star}
\end{align}
The condition \eqref{as2} ensures that $\Delta_S \geq 0$, i.e., that the residual treatment effect is in the same direction as the overall treatment effect $\Delta$. Condition \eqref{as1star}, which ensures that $\Delta \geq \Delta_S$, requires that, roughly speaking, a propensity-weighted mixture of the two conditional mean functions $e({\bf x})\psi_1({\bf x}, {\bf s}) + \{1-e({\bf x})\}\psi_0({\bf x}, {\bf s})$ is larger when ${\bf s}$ takes values from the distribution of the counterfactual surrogates under treatment than if it took values from the distribution under control.  These two conditions are analogues of conditions (A6) and (A7) in \cite{agniel2020evaluation}.

\subsection{Identification}

In this section, we show how $\Delta_{{\bf S}}$ may be characterized via functionals of the observed data distribution. Specifically, 
\begin{align}
    \Delta_{{\bf S}} &= \mathbb{E}\left\{\mu_1({\bf X}, {\bf S}) - \mu_0({\bf X}, {\bf S})\right\}\label{om}\\
    &= \mathbb{E}\left[\frac{AY}{\pi({\bf X}, {\bf S})} - \frac{(1-A)Y}{\{1-\pi({\bf X}, {\bf S})\}}\right]\label{ipw}\\
    &= \mathbb{E}\left[\frac{AY - \{A - \pi({\bf X}, {\bf S})\}\mu_1({\bf X}, {\bf S})}{\pi({\bf X}, {\bf S})} - \frac{(1-A)Y + (A - \pi({\bf X}, {\bf S}))\mu_0({\bf X}, {\bf S})}{1-\pi({\bf X}, {\bf S})}\right],\label{dr}
\end{align} 
where $\mu_a({\bf x}, {\bf s}) = \mathbb{E}(Y | {\bf X} = {\bf x}, A = a, {\bf S} = {\bf s})$. The result \eqref{om} follows from the definition \eqref{def:delta-s} and the fact that $\psi_a({\bf X}, {\bf S}^{(a)}) = \mu_a({\bf X}, {\bf S})$ because of \eqref{psr}, and the other results follow shortly thereafter, following familiar paths as arguments for identification of the average treatment effect in other contexts. 
These results show that the residual treatment effect may be identified without knowledge of the mean function for the outcome via \eqref{ipw}, and gives an augmented inverse probability weighting \cite{bang2005doubly} version of the estimand \eqref{dr} .$\Delta$ may also be identified using similar functionals with $\mu_a({\bf X}, {\bf S})$ replaced by $m_a({\bf X})$ and $\pi({\bf X}, {\bf S})$ replaced by $e({\bf X})$:
\begin{align}
    \Delta &=  \mathbb{E}\left[\frac{AY - \{A - e({\bf X})\}m_1({\bf X})}{e({\bf X})} - \frac{(1-A)Y + (A - e({\bf X}))m_0({\bf X})}{1-e({\bf X})}\right]\label{drdelta}
\end{align}
where $m_a({\bf x}) = E(Y | {\bf X} = {\bf x}, A = a)$. In Section \ref{est-sec}, we propose an estimation approach using the expressions in \eqref{dr} and \eqref{drdelta}.

\begin{remark}
If one wished to pursue the alternative PTE definition using ${\cal M}_\psi$ and ${\cal M}_0$, these quantities may also be identified from the observed data. Let $\eta^2_a({\bf x}, {\bf s}) = \mathbb{E}\{Y - \mu_a({\bf x}, {\bf s})\}^2$ be the conditional (on ${\bf X}$ and ${\bf S}$) variance of $Y$ in treatment arm $A = a$. Further define $\sigma_a^2({\bf x}) = \mathbb{E}\left\{\eta^2_a({\bf x}, {\bf S}) | {\bf X} = {\bf x}, A = a\right\}$ to be the expectation of $\eta^2_a({\bf x}, {\bf S})$ over the conditional distribution of ${\bf S}$ given $A = a, {\bf X} = {\bf x}$. Then using a similar argument as above, we have that 
\begin{align*}
    {\cal M}_\psi &= \mathbb{E}\left[e({\bf X})\sigma^2_1({\bf X}) + \{1-e({\bf X})\}\sigma^2_0({\bf X})\right]\\ &= \mathbb{E}\left(A\eta_1({\bf X}, {\bf S}) + (1-A)\eta_0({\bf X}, {\bf S}) \right) \\
    &= \mathbb{E}\left(A\eta_1({\bf X}, {\bf S}) - \{A - e({\bf X})\}\sigma^2_1({\bf X}) + (1-A)\eta_0({\bf X}, {\bf S}) + (A - e({\bf X}))\sigma^2_0({\bf X})\right)
\end{align*}
with similar expressions holding for ${\cal M}_0$ \textit{mutatis mutandis}. The first equality is analogous to an outcome-model-type identification, as in \eqref{om}, the second equality is an inverse-probability-type identification, as in \eqref{ipw}, and the third is a doubly-robust-type identification, as in \eqref{dr}. One would need to modify slightly the cross-fitting procedure outlined in the next section to incorporate $\sigma_a^2({\bf x})$, which would require two disjoint folds of data for estimation: one to estimate $\eta_a^2({\bf x}, {\bf s})$ and a second to estimate $\sigma_a^2({\bf x})$.
\end{remark} 

\subsection{Estimation \label{est-sec}}
We propose an estimation approach utilizing the augmented inverse probability weighting versions of the estimands in \eqref{dr} and \eqref{drdelta}. The traditional form of these estimators would take the form (for $\Delta_S$)
\begin{align*}
    {\widehat\Delta}_{\bf S}({\cal O}, {\widehat\pi}, {\widehat\mu}_0, {\widehat\mu}_1) = n^{-1} \sum_{i=1}^n \frac{A_iY_i - \{A_i - {\widehat\pi}({\bf X}_i, {\bf S}_i)\}{\widehat\mu}_1({\bf X}_i, {\bf S}_i)}{{\widehat\pi}({\bf X}_i, {\bf S}_i)} - n^{-1} \sum_{i=1}^n \frac{(1-A_i)Y_i + \{A_i - {\widehat\pi}({\bf X}_i, {\bf S}_i)\}{\widehat\mu}_0({\bf X}_i, {\bf S}_i)}{1-{\widehat\pi}({\bf X}_i, {\bf S}_i)}
\end{align*}
with ${\widehat\mu}_a({\bf x}, {\bf s}) = {\widehat E}\{Y | {\bf X} = {\bf x}, A = a, {\bf S} = {\bf s}\}$ an estimate of the conditional outcome function and ${\widehat\pi}({\bf x}, {\bf s}) = {\widehat P}(A = 1 | {\bf S} = {\bf s}, {\bf X} = {\bf x})$, and (for $\Delta$)
\begin{align*}
    {\widehat\Delta}({\cal O}, {\widehat e}, {\widehat m}_0, {\widehat m}_1) = n^{-1} \sum_{i=1}^n \frac{A_iY_i - \{A_i - {\widehat e}({\bf X}_i)\}{\widehat m}_1({\bf X}_i)}{{\widehat e}({\bf X}_i)} - n^{-1} \sum_{i=1}^n \frac{(1-A_i)Y_i + \{A_i - {\widehat e}({\bf X}_i)\}{\widehat m}_0({\bf X}_i)}{1-{\widehat e}({\bf X}_i)}
\end{align*}
with ${\widehat m}_a({\bf x}) = {\widehat E}\{Y | {\bf X} = {\bf x}, A = a\}$  and ${\widehat e}({\bf x}) = {\widehat P}(A = 1 | {\bf X} = {\bf x})$ estimates of the corresponding mean functions and propensity score.
We discuss the doubly robust property further in Section \ref{discussion}.

 Here, we aim to specify our approach as generally as possible, such that in practice, one could choose any reasonable approach to obtaining estimates for ${\widehat\mu}_a({\bf x}, {\bf s})$,  ${\widehat\pi}({\bf x}, {\bf s})$, ${\widehat m}_a({\bf x})$, and ${\widehat e}({\bf x})$ while maintaining the ability to make inference about the PTE $R_{\bf S}$. Certainly, if the dimensions of ${\bf X}$ and ${\bf S}$ are low, one could use simple parametric functions or, if very low (e.g., $X, S$ both scalar random variables), one could use nonparametric kernel smoothing. As we are specifically interested in settings where ${\bf S}$, and possibly even ${\bf X}$, may be high-dimensional, we do not restrict ourselves to settings where the above may be feasible, and we allow for machine learning and high-dimensional regression approaches to be used for estimating these functions. In our simulations, we specifically examine two versions of our proposed estimation approach: one that uses the Super Learner \cite{van2007super} and one that uses the relaxed lasso \cite{meinshausen2007relaxed} to estimate each of the quantities: ${\widehat\mu}_a({\bf x}, {\bf s})$, ${\widehat\pi}({\bf x}, {\bf s})$, ${\widehat m}_a({\bf x})$, and ${\widehat e}({\bf x})$. Details are described in Section \ref{sims-sec}.

Importantly, in these settings when the dimensions of ${\bf X}$ and ${\bf S}$ may be high and/or when flexible machine learning models are used to estimate the quantities $\{{\widehat\mu}_0(\cdot), {\widehat\mu}_1(\cdot), {\widehat\pi}(\cdot), {\widehat m}_0(\cdot), {\widehat m}_1(\cdot)$, ${\widehat e}(\cdot)\}$, additional steps that preclude overfitting and facilitate asymptotic analysis should be considered. We propose to use the \textit{cross-fitting} approach of \cite{Chernozhukov2017}, where (for a given positive integer $K > 1$), the observed data are randomly partitioned into $K$ sets $({\cal O}_k)_{k=1,...,K}, \bigcup_k {\cal O}_k = {\cal O}, \bigcap_k {\cal O}_k = \emptyset$, with the indices corresponding to ${\cal O}_k$ denoted $I_k$, such that each set ${\cal O}_k$ has $m = |I_k| = n/K$ observations in it (assuming $n$ is a multiple of $K$). The estimator for $\Delta_{\bf S}$ then takes the form
\begin{align}
    {\widehat\Delta}_{\bf S} = K^{-1}\sum_{k=1}^K {\widehat\Delta}_{{\bf S},k} = K^{-1}\sum_{k=1}^K {\widehat\Delta}_{{\bf S}}({\cal O}_k, {\widehat\pi}_{-k}, {\widehat\mu}_{0-k}, {\widehat\mu}_{1-k})
\end{align}
where
\begin{align*}
    {\widehat\Delta}_{{\bf S}, k} &= m^{-1} \sum_{i \in I_k} \frac{A_iY_i - \{A_i - {\widehat\pi}_{-k}({\bf X}_i, {\bf S}_i)\}{\widehat\mu}_{1-k}({\bf X}_i, {\bf S}_i)}{{\widehat\pi}_{-k}({\bf X}_i, {\bf S}_i)} - \\
    &\qquad m^{-1} \sum_{i \in I_k} \frac{(1-A_i)Y_i + \{A_i - {\widehat\pi}_{-k}({\bf X}_i, {\bf S}_i)\}{\widehat\mu}_{0-k}({\bf X}_i, {\bf S}_i)}{1-{\widehat\pi}_{-k}({\bf X}_i, {\bf S}_i)}
\end{align*}
and $({\widehat\pi}_{-k}, {\widehat\mu}_{0-k}, {\widehat\mu}_{1-k})$ are estimated on all of the data excluding ${\cal O}_k$: ${\cal O}_k^c = {\cal O} / {\cal O}_k$. 

The estimator for $\Delta$ uses a similar cross-fitting procedure:
\begin{align}
{\widehat\Delta} &= K^{-1}\sum_{k=1}^K {\widehat\Delta}_{k} = K^{-1}\sum_{k=1}^K {\widehat\Delta}({\cal O}_k, {\widehat e}_{-k}, {\widehat m}_{0-k}, {\widehat m}_{1-k})\\
{\widehat\Delta}_{k} &= m^{-1} \sum_{i \in I_k} \frac{A_iY_i - \{A_i - {\widehat e}_{-k}({\bf X}_i)\}{\widehat m}_{1-k}({\bf X}_i)}{{\widehat e}_{-k}({\bf X}_i)} - m^{-1} \sum_{i \in I_k} \frac{(1-A_i)Y_i + \{A_i - {\widehat e}_{-k}({\bf X}_i)\}{\widehat m}_{0-k}({\bf X}_i)}{1-{\widehat e}_{-k}({\bf X}_i)}
\end{align}
where ${\widehat m}_{a-k}({\bf x}) = {\widehat E}(Y | A = a, {\bf X} = {\bf x})$ is the outcome model conditional only on the pre-treatment covariates ${\bf X}$ and ${\widehat e}_{-k}({\bf x}) = {\widehat P}(A = 1 | {\bf X} = {\bf x})$ is the estimated propensity score, both estimated on ${\cal O}_k^c$. Finally, the PTE may be estimated as
\[
{\widehat R}_{{\bf S}} = \frac{{\widehat\Delta} - {\widehat\Delta}_{{\bf S}}}{{\widehat\Delta}} = 1 - \frac{{\widehat\Delta}_{{\bf S}}}{{\widehat\Delta}}.
\]




\section{Asymptotic Behavior, Inference, and Variance Estimation}\label{asymp}
In this section, we show that our approach (and specifically, the use of cross-fitting) yields a convenient asymptotic distribution for the PTE estimator (given certain regularity conditions we outline), and we demonstrate how this may be used for making inference on ${\widehat R}$. 

The key conditions are on the rate of convergence of the nuisance functions $\{{\widehat\mu}_0(\cdot), {\widehat\mu}_1(\cdot), {\widehat\pi}(\cdot), {\widehat m}_0(\cdot), {\widehat m}_1(\cdot)$, ${\widehat e}(\cdot)\}$. These rates of convergence are allowed to be quite slow. Specifically, each estimated nuisance function is allowed to converge to the true nuisance function at a rate as slow as $n^{-1/4}$. Furthermore, because of sample-splitting, the form of the asymptotic variance and its plug-in estimate are convenient to use. We now state the main result more formally
\begin{proposition}
Let
\begin{align}
    \phi_1({\bf O}; \Delta, \eta_1) = \frac{AY - \{A - e({\bf X})\}m_1({\bf X})}{e({\bf X})} - \frac{(1-A)Y + \{A - e({\bf X})\}m_0({\bf X})}{1-e({\bf X})} - \Delta,
\end{align}
and
\begin{align}
    \phi_2({\bf O}; \Delta_{\bf S}, \pi, \mu_0, \mu_1) = \frac{AY - \{A - \pi({\bf X}, {\bf S})\}\mu_1({\bf X}, {\bf S})}{\pi({\bf X}, {\bf S})} - \frac{(1-A)Y + \{A - \pi({\bf X}, {\bf S})\}\mu_0({\bf X}, {\bf S})}{1-\pi({\bf X}, {\bf S})} - \Delta_{\bf S}.
\end{align}
If, in addition to the assumptions in Section \ref{assump-ident} and Appendix \ref{assump-reg}, we also have, for a sequence of positive constants $\delta_n, n= 1, ..., \infty$ that 
\begin{align*}
    &E\left[\{{\widehat\mu}_{a-k}({\bf X}, {\bf S}) - \mu_a({\bf X}, {\bf S})\}^2\right] \times E\left[\{{\widehat\pi}_{-k}({\bf X}, {\bf S}) - \pi({\bf X}, {\bf S})\}^2\right] \leq \delta_n n^{-\half},\\
    &E\left[\{{\widehat\mu}_{a-k}({\bf X}, {\bf S}) - \mu_a({\bf X}, {\bf S})\}^2\right] + E\left[\{{\widehat\pi}_{-k}({\bf X}, {\bf S}) - \pi({\bf X}, {\bf S})\}^2\right] \leq \delta_n, \\
    &P\{\epsilon_2 < {\widehat\pi}_{-k}({\bf S}, {\bf X}) < 1-\epsilon_2\} = 1
\end{align*}
and 
\begin{align*}
    &E\left[\{{\widehat m}_{a-k}({\bf X}) - m_a({\bf X})\}^2\right] \times E\left[\{{\widehat e}_{-k}({\bf X}) - e({\bf X})\}^2\right] \leq \delta_n n^{-\half},\\
    &E\left[\{{\widehat m}_{a-k}({\bf X}) - m_a({\bf X})\}^2\right] + E\left[\{{\widehat e}_{-k}({\bf X}) - e({\bf X})\}^2\right] \leq \delta_n,\\
    &P\{\epsilon_1 < {\widehat e}_{-k}({\bf X}) < 1-\epsilon_1\} = 1,
\end{align*}
then
\begin{align*}
    &\sqrt{n}\sigma^{-1}({\widehat R}_{\bf S} - R_{{\bf S}}) \rightarrow N(0, 1), \quad \mbox{where}\\
    &\sigma^2 = \Delta^{-2}E\left\{\phi_1({\bf O}, \Delta, e, m_0, m_1)^2\right\} + \Delta_{{\bf S}}\Delta^{-4}E\left\{\phi_2({\bf O}, \Delta_{\bf S}, \pi, \mu_0, \mu_1)^2\right\} -\\
    &\qquad 2\Delta_{{\bf S}}\Delta^{-3}E\left\{\phi_1({\bf O}, \Delta, e, m_0, m_1)\phi_2({\bf O}, \Delta_{\bf S}, \pi, \mu_0, \mu_1)\right\}.
\end{align*}
\end{proposition}
The consequence of this result is that a straightforward $(1-\alpha)\%$ confidence interval for $R_{\bf S}$ can be obtained as 
\begin{align}\label{asymp-ci}
    {\cal C} = ({\widehat R}_{\bf S} - z_{1-\alpha/2}{\widehat\sigma}, {\widehat R}_{\bf S} + z_{1-\alpha/2}{\widehat\sigma}),
\end{align} with $z_q$ the $q$th quantile of the standard normal distribution and
\begin{align*}
    &{\widehat\sigma}^2 = n^{-1}\sum_{i=1}^n\left({\widehat\Delta}^{-2}{\widehat\phi}_{1i}^2 + {\widehat\Delta}_{{\bf S}}^2{\widehat\Delta}^{-4}{\widehat\phi}_{2i}^2 - 2{\widehat\Delta}_{{\bf S}}{\widehat\Delta}^{-3}{\widehat\phi}_{1i}{\widehat\phi}_{2i}\right).
\end{align*}
where 
\[{\widehat\phi}_{1i} = \phi_1({\bf O}_i, {\widehat\Delta}, {\widehat e}_{-k_i}, {\widehat m}_{0-k_i}, {\widehat m}_{1-k_i}), {\widehat\phi}_{2i}=\phi_2({\bf O}_i, {\widehat\Delta}_{\bf S}, {\widehat\pi}_{-k_i}, {\widehat\mu}_{0-k_i}, {\widehat\mu}_{1-k_i})\}^{\transpose}
\] and $k_i = I\{k : i \in I_k\}$. In addition, it follows that $P(R_{\bf S} \in {\cal C}) \rightarrow 1- \alpha$ because, under these assumptions, based on the arguments in \cite{Chernozhukov2017},
\begin{align}
    \sqrt{n}\Sigma^{-\frac{1}{2}}({\widehat\btheta} - \boldsymbol{\theta})  \rightarrow N({\bf 0}, I)
\end{align}
where ${\widehat\btheta} = ({\widehat\Delta}, {\widehat\Delta}_{\bf S})^{\transpose}, \boldsymbol{\theta} = (\Delta, \Delta_{\bf S})^{\transpose}, \Sigma = \mathbb{E}(\boldsymbol{\phi}\bphi^{\transpose})$ and this holds when replacing $\Sigma$ by ${\widehat\Sigma} = n^{-1}\sum_{i=1}^n\boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\transpose}$, ${\widehat\bphi}_i = (\phi_{1i}, \phi_{2i})^{\transpose}$.








\section{Simulations \label{sims-sec}}

\subsection{Simulation Overview}
We used two approaches to estimate the six functions $\pi({\bf X}, {\bf S}), \mu_1({\bf X}, {\bf S}), \mu_0({\bf X}, {\bf S}), \pi({\bf X}), m_1({\bf X}), m_0({\bf X})$, yielding two different versions of our proposed estimator. In one, we used the Super Learner \cite{van2007super}, which finds an optimal combination of a set of candidate models or learners, as implemented in the \texttt{SuperLearner} package \cite{polley2019package} in \texttt{R}. We denote this estimator ``DR-SL" in what follows. In the second version, which we call "DR-lasso" we use the relaxed lasso \cite{meinshausen2007relaxed}, as implemented in the \texttt{glmnet} package \cite{glmnet}, to estimate all needed functions. For all simulations, we truncated estimates of the propensity and surrogate scores so that $0.01 \leq {\widehat\pi}_{-k}({\bf X}_i, {\bf S}_i), {\widehat e}_{-k}({\bf X}_i) \leq 0.99$ for all $i$.

It is our understanding that there are no currently available methods to estimate the PTE of a high-dimensional surrogate. However, since there are available methods to measure the strength of high-dimensional mediators, we compare our proposed approach to these available methods. While the goals of mediation analysis and surrogate markers analysis are different and the assumptions necessary differ, this allows us to offer some reasonable comparison to methodology that is currently available, rather than no comparison at all. Thus, to fairly compare, we purposefully set up our simulations such that ${\bf S}$ is a mediator in that it lies on the causal pathway between $A$ and $Y$. While mediation methods are attempting to estimate a distinct quantity from the PTE we are interested in, in the following simulations, the estimation of the ``proportion mediated" (which is used in mediation) and the PTE that we are interested in are, in fact, the same. 

Specifically, in our simulations we compare our proposed approach to three methods for high-dimensional mediation analysis: HIMA (High-dimensional mediation analysis, \cite{zhang2016estimating}) as implemented in the \texttt{HIMA} package \cite{hima}; BAMA (Bayesian mediation analysis, \cite{song2018bayesian}) as implemented in the \texttt{bama} package \cite{bama}; and HILMA (High-dimensional linear mediation analysis, \cite{zhou2020estimation}) as implemented in the \texttt{freebird} package. All three of these methods propose $p+1$ linear models: $p$ models for the surrogates/mediators
\[
S_{ij} = \alpha_{0j} + \alpha_jA_i + \sum_{k=1}^q\boldsymbol{\gamma}^*_{jk}X_{ik} + e_{ij}, j = 1,...,p; i = 1, ..., n,
\]
and one outcome model
\[
Y_i = \beta_0 + \Delta_S A_i + \sum_{j=1}^p\beta_jS_{ij} + \sum_{k=1}^q\gamma_kX_{ik} + \epsilon_i,
\]
though the implementation of freebird does not allow for the inclusion of covariates $X_{ik}$.
Using these models, the overall treatment effect can be identified as $\Delta = \Delta_S + \sum_{j=1}^p\alpha_j\beta_j$, and as above the PTE (or proportion mediated) may be identified as $R = 1 - \Delta_S/\Delta.$ The HIMA procedure employs sure independence screening and penalized least squares with the minimax concave penalty \cite{zhang2010nearly} to set some values of $\alpha_j\beta_j$ to be exactly 0. Similarly, BAMA employs Bayesian shrinkage estimation to estimate all parameters (including some at exactly 0) under normality assumptions on the error terms. Finally, the freebird approach debiases a pilot scaled lasso estimator. 

We computed 95\% confidence intervals (CIs) for $R$ from the results given in Section \ref{asymp} for the two proposed approaches. We computed 95\% credible intervals for the BAMA estimator from the 2.5\% and 97.5\% quantiles of the posterior distribution of $1 - \frac{\Delta_S}{\Delta_S + \sum_{j=1}^p\alpha_j\beta_j}.$ CIs were computed for the HIMA and freebird implementations because they are not available. 

\subsection{Simulation Setup}
We constructed two sets of simulations for a total of 10 simulation settings to assess the performance of our proposed approach in both low- and high-dimensional settings. 

For the first set of simulations, the data-generating mechanism for ${\bf S}^{(a)}$ was linear in ${\bf X}$, the data-generating mechanism for $Y^{(a)}$ was linear in ${\bf X}$ and ${\bf S}^{(a)}$, and the propensity score was linear on the log-odds scale. Given this data-generating mechanism, we would expect that all methods (proposed and comparisons) should perform reasonably well. We set the dimension of ${\bf S}$ ($p$) and of ${\bf X}$ ($q$) to be 100. We let  $X_{ij} \sim N(0, 1), A_i \sim \text{Bernoulli}\{\pi({\bf X}_i)\}, \pi({\bf X}_i) = \text{logit}(\boldsymbol{\gamma}^{\transpose}{\bf X}_i)$, where $\boldsymbol{\gamma} \sim N(0, 1)$. The surrogates were generated as ${\bf S}^{(a)}_i = \boldsymbol{\alpha}_a + \boldsymbol{\beta}_a^{\transpose}{\bf X}_i + {\bf e}_i$, where $\alpha_{11} = 0.75, \alpha_{12} = 0.25, \alpha_{01} = \alpha_{02} = 0, \alpha_{1j} \sim U(0,1), \alpha_{0j} \sim U(-0.5, 0.5), j = 3, ... p$. And $(\beta_{11}, \beta_{12}, \beta_{13}, \beta_{14}, \beta_{15}) = (-1, -0.5, 0, 0.5, 1), (\beta_{01}, \beta_{02}, \beta_{03}, \beta_{04}, \beta_{05}) = (-2, -1.5, -1, -0.5, 0), \beta_{aj} = 0, a=0,1, j > 5$. The outcome counterfactuals were given by $Y_i^{(a)} = a\Delta_S + \sum_{j = 1}^{25} X_{ij} + S_{i1}^{(a)} + S_{i2}^{(a)} + \epsilon_i$. The errors $e_{ij}$ and $\epsilon_i$ were $N(0, \sigma^2)$. We let the sample size vary between $n = 100$ and $n = 500$ and the level of noise between $\sigma = 0.1$ and $\sigma = 0.5$, and we set $R = 0.5$. 

In the second set of simulations, the data-generating mechanism was less linear, including interactions between covariates in both the propensity score and the model for ${\bf S}$, though the outcome model was still a simple linear combination of ${\bf X}$ and ${\bf S}$. This set of simulations should mimic what may happen in practice since all models are typically subject to some amount of mis-specification. We expected the nonlinearity to induce bias in the competing methods (which require linear models to hold). There were two pre-treatment covariates $X_1 \sim \text{Uniform}(-2, 5)$ and $X_2 \sim \text{Bernoulli}(0.5)$. The probability of treatment was determined by $P(A = 1 | {\bf X}) = \text{logit}(-X1 + 2X_1X_2)$. The surrogates under treatment were given by ${\bf S}^{(1)} = 1.5 + (X_1 + X_2){\cal I}_{10} + {\bf e}$ and under control by ${\bf S}^{(0)} = 2 + X_2{\cal I}_{10} - X_1X_2 + {\bf e}$ where ${\cal I}_{10}$ is a $p \times p$ matrix of 0s except it has 1s on the first 10 diagonal elements and ${\bf e} \sim N(0, \Sigma)$ and $\Sigma = 0.4{\bf 1}\bone^{\transpose} + 0.6I_p$ is a covariance matrix with common covariance $0.4$. Finally, the outcome counterfactuals were given by $Y^{(a)} = a\Delta_S + X_1 + X_2 + \sum_{j = 1}^{15}S_j^{(a)} + \epsilon, \epsilon \sim N(0, 1)$. We set the sample size to be $n=500$, and we varied the number of surrogates ($p = 50, 100, 500$), and the PTE (surrogates were either weak ($R = 0.41$) or strong ($R = 0.85$)).


\subsection{Results for first data-generating mechanism}
As expected, this data-generating mechanism was well approximated using linear models. When the sample size was large ($n = 500$) and $\sigma = 0.5$, all methods performed well, with relatively small bias (see Figure \ref{r-pl1}). The median ${\widehat R}$ estimates were 0.46 (BAMA), 0.48 (freebird), 0.49 (HIMA), 0.48 (DR-SL), and 0.49 (DR-lasso), all of which compare favorably to the true $R$ value of 0.5. The distribution of estimates for the proposed approaches was a bit tighter than for the competing methods. The empirical 2.5\% and 97.5\% quantiles were 0.30 and 0.62 for DR-SL, 0.35 and 0.58 for DR-lasso, 0.15 and 0.67 for HIMA, 0.11 and 0.63 for BAMA, and -0.85 and 0.71 for freebird. Results were similar for the sample size with less noise ($\sigma = 0.1$), except for the freebird approach which began to estimate the PTE to be exactly 1 in almost all simulations: the empirical 2.5\% and 97.5\% quantiles of the ${\widehat R}$ distribution for freebird were exactly 1. At the lower sample size ($n = 100$), HIMA did not run, while BAMA estimates tended to be near 0 and freebird estimates were again clustered tightly around 1. Our proposed approaches had median ${\widehat R}$ values of 0.46 (DR-SL) and 0.51 (DR-lasso), though with more variability than when $n = 500$. The median absolute deviation ($|{\widehat R} - R|$) was smaller for both DR-SL and DR-lasso than any of the competing approaches for all settings.

CI coverage for the proposed estimators tended to be conservative in all settings: coverage for DR-SL was 100\% ($n = 100, \sigma = 0.1$ and $n = 100, \sigma = 0.5$), 98\% ($n = 500, \sigma = 0.1$), and 99\% ($n = 500, \sigma = 0.5$), and for DR-lasso it was 100\% in all settings.
On the other hand, BAMA CIs only had nominal coverage at lower sample sizes: 95\% ($n = 100, \sigma = 0.1$), 98\% ($n = 100, \sigma = 0.5$), 88\% ($n = 100, \sigma = 0.1$), and 84\% ($n = 500, \sigma = 0.5$). Even when BAMA CIs had nominal coverage, they were about four times larger than the CIs for the proposed estimators: BAMA CI half-lengths were 1.38 ($n = 100, \sigma = 0.1$) and 1.64 ($n = 100, \sigma = 0.5$), while the corresponding half-lengths were 0.47 and 0.57 for DR-SL and 0.57 and 1.02 for DR-lasso.

\begin{figure}
\centering
\begin{tabular}{c}
\includegraphics[width =\textwidth]{sim2-ridgeplot.png}
\end{tabular}\vspace{0.2in}
\caption{Distribution of estimates of ${\widehat R}$ in first data-generating mechanism. Lighter shaded regions represent the distribution of the proposed estimators (``DR-SL" and ``DR-lasso"); darker shaded regions represent the distribution of the comparison estimators (``BAMA", ``freebird", and ``HIMA"). The true value of $R$ is given as a vertical dotted line at 0.5.}
\label{r-pl1}
\end{figure}

\subsection{Results for second data-generating mechanism}
This data-generating mechanism was not as well approximated by linear models, in part because of the interactions between $X_1$ and $X_2$. Thus, the linear model-based estimates of HIMA, BAMA, and freebird all performed poorly (see Figure \ref{r-pl2}). HIMA and BAMA typically had ${\widehat R}$ values around 0 and lower, while freebird again estimated $R$ to be 1 or above (the 2.5\% quantile of the empricial distribution of freebird estimates was at 1 in all settings). On the other hand, the two proposed approaches gave more reasonable if imperfect estimates. When the dimensionality of the surrogates was much smaller than the sample size ($p = 50$), the DR-SL estimator gave nearly unbiased results, with median ${\widehat R}$ values of 0.88 (when $R = 0.85$) and 0.41 (when $R = 0.45$), as it was able to adapt to the nonlinearity in the data-generating mechanism. As the dimensionality of the surrogate grew, the Super Learner estimates of nuisance functions began to lean more heavily on linear models (like the lasso), and performance degraded: median ${\widehat R}$ values were 0.78 ($R = 0.85, p = 100$), 0.37 ($R = 0.41, p = 100$), 0.65 ($R = 0.85, p = 500$), 0.31 ($R = 0.41, p = 500$). The DR-lasso estimator resulted in performance between HIMA/BAMA and DR-SL, typically underestimating $R$ by more than DR-SL but less than HIMA/BAMA. CI coverage was strong for DR-SL when $R = 0.41$ (97\%, 96\%, 98\% for $p = 50, 100, 500$, respectively) but was much weaker when $R = 0.85$ (91\%, 74\%, 69\% for $p = 50, 100, 500$, respectively). This subpar coverage appeared to be due to the bias of the estimator and not due to the width of confidence intervals. When CIs were adjusted for the empirical bias of the DR-SL estimator (by subtracting the empirical bias from both sides of CIs), nominal coverage was observed (96\%, 97\%, 95\% coverage for $p = 50, 100, 500$, respectively).

\begin{figure}
\centering
\begin{tabular}{c}
\includegraphics[width =\textwidth]{sim1-ridgeplot.png}
\end{tabular}\vspace{0.2in}
\caption{Distribution of estimates of ${\widehat R}$ in first data-generating mechanism. Lighter shaded histograms represent distribution of the proposed estimators (``DR-SL" and ``DR-lasso"); darker shaded histograms represent the distribution of the comparison estimators (``BAMA", ``freebird", and ``HIMA"). The true value of $R$ is given as a vertical dotted line.}
\label{r-pl2}
\end{figure}




\section{Discussion \label{discussion}}

We have proposed a very general approach to evaluating surrogate markers which can be applied in randomized experiments or in observational studies and can be used regardless of the dimensionality of either surrogates or any covariates needed to control for confounding. Our approach is nonparametric in that we have defined the proportion of the treatment effect explained by the surrogates without reference to any models, and we have shown how machine learning approaches like the super learner may be used to very flexibly estimate nuisance functions. Our simulation results suggest that our the Super Learner estimator (the estimator we called DR-SL) outperforms competing methods even when the underlying data-generating mechanism is linear and still gives reasonable results even when high dimensional linear models are mis-specified.

We note that the proposed estimator of the PTE, ${\widehat R}_{\bf S}$, is ``doubly-robust" in the sense of \textit{rate double robustness} discussed in \cite{smucler2019unifying}, where the rates of two estimators (for us, the estimators of $\mu_a$ and $\pi$ as well as the estimators of $m_a$ and $e$) may compensate for one another -- i.e., a slightly slower rate of convergence for $\pi$ may be compensated for by a slightly faster rate of convergence for $\mu_a$ so long as their product is $o_p(n^{-\half})$. In high dimensions, the approach of \cite{smucler2019unifying}, based on $\ell_1$-penalized regression, could be used to additionally obtain the slightly different property of \textit{model double robustness} where one of the estimators may converge to a different value besides the true parameter. Our development here, which we chose to be quite general and applicable in both low- and high-dimensional settings, could easily be adapted to use the estimation technique of \cite{smucler2019unifying} in high dimensions to recover model double robustness.

Our approach here is intimately tied to previous approaches for evaluating surrogate markers. The approaches in \cite{agniel2020evaluation} and in \cite{parast2016robust} can be seen to be similar to a version of our approach where the average treatment effect among controls is estimated under further assumptions. In \cite{agniel2020evaluation}, they further require strict randomization of treatment,
and they assume that ${\bf S}$ is a realization of a smooth continuous function. \cite{parast2016robust} make similar assumptions but take $S$ to be a scalar surrogate. They estimate a version of \eqref{om} among controls: $\Delta_{{\bf S}} = \mathbb{E}\left\{\mu_1({\bf S}) - \mu_0({\bf S}) | A = 0\right\}$ using kernel smoothing and taking advantage of the fact that treatment is randomized. Their estimates have the form
\[
{\widehat\Delta}_{{\bf S}} = n_0^{-1}\sum_{A_i = 0}\left\{{\widehat\mu}_1({\bf S}_i) - Y_i\right\}
\]
where $n_0 = \sum_{i=1}^n I\{A_i = 0\}$ and ${\widehat\mu}_1(\cdot)$ is estimated via kernel smoothing (possibly after dimension reduction). Our approach here could easily be adapted to estimate a similar quantity by using methods for doubly robust estimation of the average treatment effect on the treated \cite{shu2018improved, moodie2018doubly} by, for example,
\[
{\widehat\Delta}_{{\bf S}} = n^{-1}\sum_{i=1}^n\left\{A_i\frac{1-{\widehat\pi}({\bf S}_i)}{{\widehat\pi}({\bf S}_i)} - (1-A_i)\right\}\left\{{\widehat\mu}_1({\bf S}_i) - Y_i\right\}
\]
where ${\widehat\pi}({\bf s}) = {\widehat P}(A_i | {\bf S}_i = {\bf s})$ may be estimated using a model similar to the one used for ${\widehat\mu}_1({\bf s})$. Furthermore, our cross-fitting approach could be used to facilitate the use of machine learning methods in the estimation of ${\widehat\mu}_1(\cdot), {\widehat\pi}(\cdot)$, to include covariates to control for confounding, or to simplify or strengthen asymptotic results -- e.g., results for the kernel estimator used in \cite{agniel2020evaluation} obtain rates of convergence of $n^{-\half}$ only under very limited technical conditions.

Our proposed robust approach to examine the strength of a high-dimensional surrogate will be useful in practice as there is increased attention on utilizing high-dimensional data that can be measured earlier or with less cost compared to primary outcomes, in an effort to make conclusions about treatments and interventions sooner. We include software to implement our proposed methods in the R package \texttt{crossurr} available at \texttt{github.com/denisagniel/crossurr}. 

\bibliographystyle{plain}
