Nonparametric assessment of regimen response curve estimators

Our data consists of $n$ i.i.d. samples of $O_{i}=\{S_{1i},A_{1i},\dots,S_{Ti},A_{Ti},Y_{i}\}^{n}_{i=1}\sim\mathbb{P}_{0}\in\mathcal{M}$ where $\mathcal{M}$ is a nonparametric model. We denote $S_{t}$ as the vector of patient characteristics at time $t$ and $V\subseteq S_{1}$ as a subset of the baseline characteristics. The observed treatment is denoted $A_{t}\in\{0,1\}$ at time $t=1,\cdots,T$, and the outcome of interest is denoted $Y$. We use the overbar notation to denote the history of a variable. For example, $\bar{S}_{t}=(S_{1},\dots,S_{t})$ and $\bar{A}_{t}=(A_{1},\dots,A_{t})$ are the covariates and past treatment up to time $t$. We let $\bar{H}_{t}=(\bar{A}_{t-1},\bar{S}_{t})$ be all the information up to time $t$. We denote $\mathbb{P}_{n}$ as the empirical distribution and define $\mathbb{P}f=\int f(o)\,\mathrm{d}\mathbb{P}(o)$ for a given function $f(o)$ and distribution $\mathbb{P}$.

A dynamic regime consists of a set of decision rules that assigns the treatment at each stage based on the information at and before that stage. We define the deterministic dynamic treatment rule $d^{\theta}=I(\theta_{t}^{T}S_{t}>0)$ to be the decision rule at time $t$. $\pi_{0}^{t}\equiv P(A_{t}=d_{t}^{\theta}|\bar{H}_{t})$ is the true propensity score at time $t$ and $\pi_{n}^{t}$ denote the corresponding estimator. And, $Y^{\theta}$ denotes the potential outcome under the regime $\theta$. Suppose that we have the full data $X=\{(Y^{\theta}:d^{\theta}\in\mathcal{D}),V\}\sim P_{X}\in\mathcal{M}^{F}$ where $\mathcal{M}^{F}$ is a nonparametric full data model. We can then define the regimen response curve $m_{0}(\theta,V):=E_{P_{X}}(Y^{\theta}\mid V)$.

A dynamic marginal structural model can be used to approximate $m_{0}(\theta,V)$ by imposing a working model. We quantify the quality of the imposed working model using a cross-validated risk function. Let $B_{n}\in\{0,1\}^{n}$ denote the random partition of $\{1,\dots,n\}$ into two folds. A realization of $B_{n}=\{B_{n,1},B_{n,2},\cdots,B_{n,n}\}$ defines a particular split of the sample of $n$ observations into a training set $\{i\in\{1,\cdots,n\}:B_{n,i}=0\}$ and a validation set $\{i\in\{1,\cdots,n\}:B_{n,i}=1\}$. Denote by $\mathbb{P}_{n,b}^{0}$ and $\mathbb{P}_{n,b}^{1}$ the empirical distribution of the training sample and the validation sample. The type of cross-validation procedure is determined by the specific distribution of $B_{n}$. For example, in Monte Carlo cross-validation, where the sample is divided into training and validation sets with probability $p$, the distribution of $B_{n}$ puts mass $1/\binom{n}{np}$ on each binary realization of $B_{n}$. Another example is $B$-fold cross-validation, in which the distribution of $B_{n}$ puts mass $1/B$ on each possible realization. Throughout, we assume that this randomness is independent of all data.

Consider an estimation rule $m_{n}(\theta,V)\equiv m_{n}(\theta,V\mid\mathbb{P})$ that is fit using the distribution $\mathbb{P}$; we refer to $m_{n}$ as a marginal structural model (MSM) to be consistent with the literature, although it is the fit from a working model rather than a model itself. Let $m_{n,b}(\theta,V)\equiv m_{n}(\theta,V\mid\mathbb{P}_{n,b}^{0})$ denote the MSM $m_{n}(\theta,V)$ evaluated at the empirical distribution $\mathbb{P}_{n,b}^{0}$, i.e. the training data. We denote $m^{*}(\theta,V)$ as the pointwise probability limit of $m_{n,b}(\theta,V)$. Note, $m^{*}$ does not necessarily equal the true regimen-response curve $m_{0}$.

Define the full data conditional risk based on training observations as

where $\mathbb{E}_{B_{n}}$ denotes an expectation over possible partitions and $F$ is a user-specified distribution function over the parameter $\theta$. The averaging over the partitions in (1) prevents the estimand from depending on any randomness in the selection of training data, thus making $\Psi_{0}(m_{n})$ a data-adaptive estimand (Hubbard et al., 2016). We define the full non-data-adaptive for a given regimen-response curve $m$ as $\Psi_{0,\text{nda}}(m)\equiv\mathbb{P}_{0}\left\{\int\{Y^{\theta}-m(\theta,V)\}^{2}\,\mathrm{d}F(\theta)\right\}.$ The full data conditional risk, $\Psi_{0}(m_{n})$ can be expressed in term of non-data-adaptive risk, $\Psi_{0,\text{nda}}(m)$, through the relation, $\Psi_{0}(m_{n})=\mathbb{E}_{B_{n}}\left[\Psi_{0,\text{nda}}(m_{n,b})\right]$. In Section 3.1, we develop efficiency theory for $\Psi_{0,\text{nda}}(m)$ of a given regimen-response curve, $m$. Then, we will use those results to construct a nonparametric asymptotically linear estimators for $\Psi_{0}(m_{n})$.

For a given marginal structural model $m_{n}$, the definition of the risk $\Psi_{0}(m_{n})$ involves the potential outcome $Y^{\theta}$ which is not observed. However, under the following causal assumptions, it is identifiable using the observed data. Recall that $\pi_{0}^{t}\equiv P(A_{t}=d_{t}^{\theta}|\bar{H}_{t})$ is the true propensity score at time $t$.

Assumption 1a states that the potential outcome under an observed regime equals the observed value. This assumption can be violated if the treatment assignment of a unit impacts the outcome of others (Hernán and Robins, 2023; Kennedy, 2019). Assumption 1b states that there are no unmeasured confounders. This assumption holds in standard randomized studies. However, in observational studies it requires a rich enough set of characteristics $S_{1},\dots,S_{T}$. Assumption 1c states that regardless of a subject’s medical history, that person will have a nonzero probability of following any given treatment regime.

Applying the plug-in approach to Lemma 1, we have the inverse probability weight (IPW) estimator for the full data conditional risk of a working model $m_{n}(\theta,V)$.

where $\pi^{t}_{n,b}$ is the estimate of the propensity score applied to the training sample for the $b^{th}$ sample split. The consistency and convergence rate of $\Psi_{n,m}^{\mathrm{IPW}}(\pi_{n})$ depends entirely on the consistency and convergence rate of the estimators $\pi^{t}_{n}$ of $\pi^{t}_{0}$ for each time point $t=1,\dots,T$. It is common to estimate the propensity score using parametric models as it has the convergence rate of $O_{p}(n^{-1/2})$. Alternatively, data-adaptive methods can be used to increase the flexibility of the model. However, due to a slow convergent rate of nonparametric methods, the IPW estimator could be biased. We can see that through the following decomposition.

For a given marginal structural model $m(\theta,V)$, define the loss as

Then, we have

Assuming $\pi_{n}$ is a consistent estimator of $\pi_{0}$ and Donsker class condition, using standard empirical process theory, we can show that term (III) will converge to 0 at rate $n^{-1/2}$ (van der Vaart, 1998). Hence, $\Psi_{n,m}^{\mathrm{IPW}}$ is asymptotically linear only if term (II) also converges to 0 at the rate $n^{-1/2}$. Since that rate is not obtainable using nonparametric methods, term (II) will dominate the right hand side, causing $\Psi_{n,m}^{\mathrm{IPW}}$ to be biased. The same issue arises when $m_{n}(\theta,V)$ is a cross-fitted estimator of the underlying marginal structural model.

Our proposed approach for refining the inverse probability weighting estimator of $\Psi_{0}(m_{n})$ utilizes the Highly Adaptive Lasso model (HAL). The detail of this process will be discussed in Section 4.1. To provide the context, in this section, we will give a concise introduction to the HAL model and highlight its key characteristics. For an in-depth exploration of HAL, we refer our interested readers to Benkeser and van der Laan (2016); Ertefaie et al. (2023)

The highly adaptive lasso is a nonparametric regression function that estimates infinite dimensional functions by constructing a linear combination of indicator functions that minimize a loss function by constraining the $L_{1}$-norm complexity of the model. (Bibaut and van der Laan, 2019) shows that highly adaptive lasso estimators have a fast convergent rate of $n^{-1/3}$ (up to a log factor).

Let $\mathbb{D}[0,\tau]$ be the Banach space of $d$-variate real-valued cadlag functions on a cube $[0,\tau]\in\mathbb{R}^{d}$. For a function $f\in\mathbb{D}[0,\tau]$, the sectional variation norm is defined as

where $f_{r}(u)=f\{u_{1}I(1\in r),\dots,u_{d}I(d\in r)\}$ is the $r^{th}$ section of $f$ and the summation is over all subset of $\{1,\dots,d\}$. The term $\int_{0_{r}}^{\tau_{r}}\left|df_{r}(u)\right|$ is the r-specific variation norm.

For any $f\in\mathbb{D}[0,\tau]$ with $\|f\|_{\nu}<\infty$, Gill et al. (1995) shows that $f$ can be represented as a linear combination of indicator functions.

The latter integral in Equation 2 can be approximated using a discrete measure $F_{n}$ putting mass on the observed values of $X_{r}=\left(X_{j}:j\in r\right)$ denoted as $x_{r,i}$ for the $i$ the subject. We let $\beta_{r,i}$ denote the mass placed on each observation $x_{r,i}$ and $\phi_{r,i}=I\left(x_{i,r}\leq c_{r}\right)$. Then the approximated $f$ is

In this formulation, $\left|\beta_{0}\right|+\sum_{r\subset\{1,\ldots,d\}}\sum_{i=1}^{n}\left|\beta_{r,i}\right|$ is an approximation of the sectional variation norm of $f$. Consequently, the minimum loss based highly adaptive lasso estimator $\beta_{n}$ is defined as

where $L(.)$ is a given loss function. This is a constrained minimization that corresponds to lasso linear regression with up to $n\times\left(2^{d}-1\right)$ parameters (i.e., $\left.\beta_{r,i}\right)$. Each value of $\lambda$ will give a unique HAL estimator.

To construct the MR estimator, we need to find the efficient influence function of $\Psi_{0,\text{nda}}(m)$. Influence functions are important because a regular and asymptotically linear estimator can be expressed as the empirical mean of the influence function plus some negligible term that converge to 0 at the rate of $n^{-1/2}$ or faster (van der Laan and Robins, 2003; Tsiatis, 2006). The efficient influence function is the influence function whose variance corresponds to the efficiency bound. The influence function is unique in the nonparametric setting and equals the efficient influence function. The multiply robust estimator is defined using the efficient influence function, and we will show that it is asymptotically efficient when all the nuisance parameters converge to their true value at a rate of $n^{-1/4}$. This slower-than-parametric rate allows us to estimate the nuisance parameters with flexible nonparametric methods. Although the MR estimator is not efficient when the propensity score model is misspecified, it is still consistent when the outcome regression is correctly specified, or vice versa.

The following result provides the nonparametric efficient influence function for the inner term in the definition of the cross-validated risk $\Psi_{0}(m_{n})$.

The structure of the influence function in Theorem 1 involves a weighted average over $\theta$ of an inverse probability weighting term and an augmentation term for each of the time points. The function $Q^{t}_{0,m}$ can be represented as a recursive sequential regression formula, which is useful in practice.

The form of the efficient influence function for $\Psi_{0,\text{nda}}$ may be used to construct an estimator for $\Psi_{0}(m_{n})$ that is asymptotically efficient for certain choices of nuisance parameter estimators (Bickel et al., 1993). The estimator construction is more delicate than usual since the efficient influence function is for $\Psi_{0,\text{nda}}(m)$, with $m$ a given regimen-response curve, while the estimand is $\Psi_{0}(m_{n})=\mathbb{E}_{B_{n}}\{\Psi_{0,\text{nda}}(m_{n}(\cdot;\mathbb{P}^{0}_{n,b}))\}$. We first develop estimators as if $m$ were given, then we discuss how to adapt these estimators to instead target our estimand $\Psi_{0}(m_{n})$.

Consider estimation of $\Psi_{0,\text{nda}}(m)$ using the efficient influence function as a guide. Two common estimation approaches are to use one-step estimators, which add the the estimated efficient influence function to an initial estimator, and estimating-equation estimators, which treat the efficient influence function as an estimating function (Tsiatis, 2006; Kennedy, 2022). Since the efficient influence function in Theorem 1 only depends on $\Psi_{0,\text{nda}}(m)$ through its subtraction, each of these estimation strategies is equivalent. Thus we initially consider the estimator

where $\pi_{n}$ and $Q_{n,m}$ are estimators for $\pi_{0}$ and $Q_{0,m}$ respectively. Assuming these nuisance parameter estimators satisfy certain conditions, we can show that the above estimator is asymptotically efficient, i.e. has the efficient influence function of $\Psi_{0,\text{nda}}(m)$, found in Theorem 1, as its influence function.

Now consider estimation of our estimand $\Psi_{0}(m_{n})=\mathbb{E}_{B_{n}}[\Psi_{0,\text{nda}}\{m_{n}(\cdot\mid\mathbb{P}^{0}_{n,b})\}]$; we construct our estimator by adapting the estimator in (3) for $\Psi_{0,\text{nda}}(m)$. A plugin-like approach readily produces the estimator

in which an outer average over folds is included, each fixed regimen-response curve $m$ is replaced with $m_{n,b}$ which is estimated in the training data $\mathbb{P}^{0}_{n,b}$ of the fold, and each nuisance parameter is replaced with estimators within the training data.

The estimator in (4) is almost a suitable estimator for the cross-validated risk $\Psi_{0}(m_{n})$, however it still has one major deficiency which harms its performance. It averages over all observations, including those in the training set $\mathbb{P}_{n,b}^{0}$ as well as the validation set $\mathbb{P}_{n,b}^{1}$, which biases the estimator downwards. A corrected estimator, which we study throughout this paper, is

where $\pi_{n,b}$ and $Q_{n,b,m_{n,b}}$ are estimates of the nuisance parameters $\pi_{0}$ and $Q_{n,m_{0}}$ calculated with the training sample for the $b^{\text{th}}$ sample split.

The next result establishes that $\Psi_{n}^{\mathrm{MR}}$ is an asymptotic efficient estimator for the estimand $\Psi_{0}(m_{n})$. Before stating it, we make two assumptions.

Assumption 2 states that the estimators $\pi_{n}^{t}$ and $Q^{t}_{n,m_{n,b}}$ need to converge to its true value at a rate of $o_{p}(n^{-1/4})$ for any $\theta$. This assumption holds when these nuisance parameters are estimated using the highly adaptive lasso. Moreover, if either $\pi_{0}^{t}$ or $Q^{t}_{n,m_{n,b}}$ converge to its true value at the rate of $O_{p}(n^{-1/2})$, we only require consistency for the other nuisance parameter. The rate of $O_{p}(n^{-1/2})$ can be achieved using parametric methods. Assumption 3 is also a relatively mild assumption. It only requires the fitted marginal structure model, $m_{n,b}$, to have a limit $m^{*}$. And $m_{n,b}$ converges to $m^{*}$ at the rate of $o_{p}(n^{-1/4})$, enabling the utilization of both parametric and nonparametric methods to estimate the marginal structure model.

Theorem 2 shows that as the fitted working marginal structural model $m_{n,b}^{0}$ converges to its limiting value $m^{*}$ at the rate of $o_{p}(n^{-1/4})$, the multiply robust estimator is asymptotic linear with respect to the limit of $m_{n}$. The centering term in Theorem 2 is the data-adaptive risk $\Psi_{0}(m_{n})$ that depends on the estimator $m_{n}$, the true propensity score and outcome regression.

When two working models are under considerations, the results of Theorem 2 allow us to compare their estimators, $m_{n,1}$ and $m_{n,2}$, by applying the asymptotic linearity result once for each estimator. We may ascertain the quality of the estimators by comparing their risks $\Psi_{0}(m_{n,1})$ and $\Psi_{0}(m_{n,2})$ by testing the informal null hypothesis that these risks are equal.

In Theorem 2, we study the difference $\Psi^{\mathrm{MR}}_{n}(m_{n},\pi_{n},Q_{n,m_{n}})-\Psi_{0}(m_{n})$. This is because $\Psi^{\mathrm{MR}}_{n}(m_{n},\pi_{n},Q_{n,m_{n}})-\Psi_{0}(m^{*})$ is asymptotically linear under much more restrictive assumptions. In the following result, we study the asymptotics of $\Psi^{MR}_{n}(m_{n},\pi_{n},Q_{n,m_{n}})-\Psi_{0}(m^{*})$ to replace the data-adaptive target parameter $\Psi_{0}(m_{n})$ with the fixed target parameter $\Psi_{0}(m^{*})$ in the conclusions of Theorem 2.

The result of Theorem 3 can be used to compare two different modeling approaches known to be consistent with $\|m_{n,1}-m_{0}\|_{2}=o_{p}(n^{-1/4})$ and $\|m_{n,2}-m_{0}\|_{2}=o_{p}(n^{-1/4})$. Despite having the same asymptotic variance, the approximated variance based on $\varphi(\pi_{0},Q_{0,m_{n,1}};m_{n,1})$ and $\varphi(\pi_{0},Q_{0,m_{n,2}};m_{n,2})$ may differ. As an illustration, one could consider comparing the super learner with a library that includes HAL alongside other parametric and nonparametric methods, against a HAL-only model.

Although doubly robust procedures facilitate the use of data adaptive techniques for modeling nuisance parameters, the resultant estimator can be irregular with large bias and a slow rate of convergence when either of the nuisance parameters is inconsistently estimated (van der Laan, 2014; Benkeser et al., 2017). To mitigate this issue, we propose a sieve based IPW estimator called UIPW where the propensity scores are estimated using an undersmoothed highly adaptive lasso. Define the estimator as

where $\pi^{t}_{n,b,\lambda}$ is the highly adaptive lasso estimate of $\pi^{t}_{0}$ obtained using the b^th sample split for a given tuning parameter $\lambda$. Smaller values of $\lambda$ imply a weaker penalty, thereby undersmoothing the fit.

Our theoretical results show that for specific choices $\lambda$, the function $\Psi_{m_{n},n}^{\mathrm{UIPW}}(\pi_{n,\lambda})$ can solve the efficient influence function presented in Theorem 1. This implication indicates that the proposed undersmoothed highly adaptive lasso estimator achieves nonparametric efficiency. This accomplishment relies on the satisfaction of two key assumptions:

Assumption 4 typically holds in practical applications because the set of cadlag functions with finite sectional variation norms is extensive.Assumption 5, on the other hand, posits that the generated basis functions in $\pi^{t}_{b,\lambda}$ are sufficient to approximate $f_{t}$ within an $O_{p}(n^{-1/4})$ neighborhood of $f_{t}$. Without undersmoothing, the basis functions can effectively approximate $\pi^{t}_{b,\lambda}$ when $\pi^{t}_{0}$ is as complex as $f_{t}$. However, when the complexity of $f_{t}$ surpasses that of $\pi^{t}_{b,\lambda}$, more undersmoothing is necessary.

As discussed in Section 2.4, without undersmoothing, the estimator $\Psi_{m_{n},n}^{\mathrm{UIPW}}$ may fail to be asymptotically linear. An IPW estimator will be efficient if it solves the efficient influence function. Theorem 4 shows that the latter is achieved when $\mathbb{P}_{n}\sum_{t=1}^{T}D^{t}_{CAR}(\pi_{n,\lambda},Q_{0,m^{*}})=o_{p}(n^{-1/2})$. The theoretical rate might not provide practical guidance due to the unknown constant. However, given the finite sample at hand, one can achieve the best performance by selecting a tuning parameter that minimizes the $\mathbb{P}_{n}\sum_{t=1}^{T}D^{t}_{CAR}(\pi_{n,\lambda},Q_{0,m^{*}})$. Specifically, we define

where $Q_{b,m_{b}}$ is a cross-validated highly adaptive lasso estimate of $Q_{0,m^{*}}$ with the $L_{1}$-norm bound based on the global cross-validation selector.

The criterion (7) relies on knowledge about the efficient influence function. However, in some cases, deriving this function can be intricate (e.g., longitudinal settings with many decision points) or may not have a closed-form solution (e.g., bivariate survival probability in bivariate right-censored data) (van der Laan, 1996). To this end, we propose an alternative criterion, drawing from the approach in Ertefaie et al. (2023), that does not require the efficient influence function. Define

in which $\lVert\beta_{t,n,\lambda,b}\rVert_{L_{1}}=\lvert\beta_{t,n,\lambda,0}\rvert+\sum_{r\subset\{1,\ldots,d\}}\sum_{j=1}^{n}\lvert\beta_{t,n,\lambda,r,i}\rvert$ is the $L_{1}$-norm of the coefficients $\beta_{t,n,\lambda,r,i}$ in the highly adaptive estimator $\pi^{t}_{n,\lambda}$ for a given $\lambda$, and $\tilde{\Omega}_{r,i}(\phi,\pi^{t}_{n,\lambda,v})=\phi_{r,i}(\bar{H}_{t})\{A_{j}-\pi^{t}_{n,\lambda,v}(1\mid\bar{H}_{t})\}\{\pi^{t}_{n,\lambda,v}(1\mid\bar{H}_{t}))\}^{-1}$.

The main drawback of undersmoothing is the decreased convergence rate of the highly adaptive lasso fit. Our asymptotic linearity proof of $\Psi_{m_{n},n}^{\mathrm{UIPW}}(\pi_{n,\lambda})$ requires $\int\|\pi^{t}_{n}-\pi^{t}_{0}\|_{2}dF(\theta)=o_{p}(n^{-1/4})$ which is slower than the cross-validated based highly adaptive lasso fit of $O_{p}(n^{-1/3})$. This ensures that we can tolerate a certain degree of undersmoothing without compromising our desired asymptotic properties. Let $J$ be a number of features included in fit (i.e., $J=|\mathcal{J}_{n}|$). van der Laan (2023) showed that the convergence rate of a highly adaptive lasso is $(J/n)^{1/2}$. Hence to ensure the required $o_{p}(n^{-1/4})$, we must choose a $J$ such that $J<n^{1/2}$ and let $\lambda_{n,J}$ be the corresponding $L_{1}$ norm. We then propose our score based criterion as

In practice, the score-based criterion tends to undersmooth too much (i.e., pick a very small $\tilde{\lambda}_{n,t}$), and the max operator helps us to mitigate this issue.

In our simulation study, we evaluate and compare the performance exhibited by our proposed estimators, MR and UIPW to the conventional IPW estimator. We consider a two-stage treatment in our simulation (T = 2). The baseline covariates $S_{1,t=1},S_{2,t=1}$ are generated from a normal distribution with the mean equals 0 and the standard deviation equals 0.5. The time-dependent covariate $S_{1,t},S_{2,t}$ are function of the preceding $S_{1,t-1},S_{2,t-1}$. In particular,

The outcome $Y=S_{1,T}-S_{2,T}+A_{T}S_{1,T}+N(0,0.1^{2})$ is a function of $S_{1},S_{2}$ in the last stage. The decision rule is given as $d^{\theta}_{t}=I\{S_{1,t}<\theta\}$ which only depends on $S_{1,t}$. We generate $\theta$ from a normal distribution $f(\theta)=N(0,0.1^{2})$. The treatment $A_{t}$ is generated from a Bernoulli distribution with $P(A_{t}=1|S_{1,t},S_{2,t})=\text{expit}(0.1+0.2S_{1,t})$. Additional simulations where $A_{t}$ is randomized can be found in Section 13 of the Supplementary Material.

We consider two variants of the UIPW estimator, differing in how the tuning parameter $\lambda$ is selected. UIPW-Dcar chooses $\lambda$ based on the criterion 7, while UIPW-Score chooses $\lambda$ based on the criterion 9. The nuisance parameters, $Q^{t}_{n,m_{n}}$ and $\pi^{t}_{n}$, is estimated using the highly adaptive lasso. We impose the parametric working model $m=E[Y|\theta]=\beta_{0}+\beta_{1}\theta$. We perform 360 simulations on the sample sizes 250, 500, 850, and 1000. In Section 13 of the Supplementary Material, we present additional simulations when the working model $m$ is estimated using other regression models, both parametric and nonparametric.

We showcase the effectiveness of our proposed estimators via coverage plots, as seen in the right panel of Figure 1. Notably, the coverage of IPW consistently falls below 77% at a sample size of 250 and diminishes further with larger sample sizes. Conversely, the MR and UIPW-Dcar estimators maintain stable coverage around 90% at small sample sizes and achieving 93% and 95% at sample sizes of 850 and 1000, respectively. Regarding UIPW-Score, its initial coverage hovers around 85% at sample sizes of 250 and 500, but rapidly climbing to approximately 92% at a sample size of 1000.

The left panel of Figure 1 provides insight into the scaled bias, calculated by multiplying the bias by the square root of the sample size. The MR estimator displays the smallest scaled bias. Additionally, as the sample size increases, our proposed estimators exhibit a declining trend in scaled bias, indicating convergence towards the true value at a rate of $o_{p}(n^{-1/2})$.

Conversely, the scaled bias of the IPW estimator shows an upward trend with increasing sample size, signifying a slower convergence rate than $o_{p}(n^{-1/2})$. This slower convergence is attributed to the bias term in the IPW estimator, which does not diminish rapidly enough when estimating the propensity score model using a nonparametric method.