Bias Analysis and Correction in Weighted-L1 Estimators for the First-Order Bifurcating Autoregressive Model

1. Introduction

The bifurcating autoregressive (BAR) model is widely employed to analyze binary tree-structured data, as illustrated in Figure 1, which frequently arise in various applications, most notably in cell lineage studies [1,2,3,4]. The BAR(1) model was initially proposed by [5] to model cell lineage data. This model extends the first-order autoregressive AR(1) model by accounting for a unique aspect: observations on two sister cells, which share the same mother cell, are correlated. In practical settings, the BAR(1) model is useful for modeling the single-cell reproduction progress [c.f. [3]. A number of studies have addressed the estimation and inference challenges associated with the BAR model [see, e.g., [5,6,7,8,9,10,11,12,13], among others].

Maximum likelihood (ML) estimators for the model coefficients and the correlation between errors for the BAR(1) model were introduced by [5] under the assumption of normally distributed errors. Building on this foundation, [8] examined the asymptotic properties of the ML estimators for the BAR(p) model, also assuming the normality of errors. In contrast, [10] explored least squares (LS) estimation for BAR(p) models without imposing distributional assumptions on the errors. In their study, they developed LS estimators for the model coefficients and introduced consistent estimators for the correlation between model errors. They further derived the asymptotic distribution of the LS estimators. However, the study did not address the finite-sample properties of these estimators. While LS estimators of BAR model coefficients are asymptotically unbiased, [13] demonstrated that the finite sample bias of these estimators can be substantial, potentially leading to inaccurate inferences. This outcome is consistent with the well-known fact that LS estimators for the AR(1) model often exhibit significant bias, especially when the autoregressive parameter is near its limitations ($\pm 1$) (e.g., [14]).

To address this issue, [13] introduced two bias correction techniques for the LS estimator of the autoregressive parameter in the BAR(1) model: one based on asymptotic linear bias functions and another using a bootstrapping approach. Their results indicated that both methods effectively reduce the bias, with the bootstrap bias-correction method proving more efficient than the asymptotic linear bias function method, particularly when the parameter is close to the upper boundary (i.e., near $+1$).

A related challenge in analyzing cell lineage data involves handling outliers. Several approaches have been proposed to develop robust estimators that mitigate the impact of outliers within the BAR model context, including quasi-maximum likelihood estimation methods (cf., [15] and [16]). In this direction, [12] developed robust estimators for the first-order BAR model based on weighted absolute deviations. They demonstrated that, when the weights are fixed, the resulting estimator is an $L_1$-type estimator, offering robustness against outliers in the response space and showing greater efficiency than the least squares (LS) estimator in the presence of heavy-tailed error distributions. In contrast, when the weights are random and dependent on the points in the factor space, the weighted $L_1$ (${\mathrm{WL}}_1$) estimator shows robustness against outliers in the factor space. Additionally, they derived the asymptotic distribution of the ${\mathrm{WL}}_1$ estimator and analyzed its large sample properties. However, they did not examine the bias of ${\mathrm{WL}}_1$ estimators under finite sample conditions (i.e., small or moderate sample sizes). The present study seeks to address this gap in the literature.

In this paper, we concentrate on examining the bias of the weighted $L_1$ (${\mathrm{WL}}_1$) estimators for BAR(1) models. Our findings indicate that the ${\mathrm{WL}}_1$ estimators of the BAR(1) autoregressive parameter (${\varphi }_1$, the parameter of interest) can show considerable bias, particularly when working with small to moderate sample sizes. To address this issue, we employ bootstrapping techniques to develop bias-corrected ${\mathrm{WL}}_1$ estimators. We assess both the bias and mean squared error (MSE) of these estimators. It is important to note that bias is typically influenced by the true, unknown parameter values, and that bias correction may not always be advisable. This is because bias reduction can lead to increased variance, potentially resulting in a higher MSE than that of the uncorrected estimates [17].

The layout of this paper is as follows. In the following section, we provide an overview of the BAR(1) model and the ${\mathrm{WL}}_1$ estimation of its parameters. Section 3 empirically examines the bias of the ${\mathrm{WL}}_1$ estimators for the autoregressive parameter in the BAR(1) model. In Section 4, we propose bias-corrected ${\mathrm{WL}}_1$ estimators. Section 5 presents a summary of the empirical findings, based on Monte Carlo simulations and a real data application, comparing the performance of the bias-corrected and uncorrected ${\mathrm{WL}}_1$ estimators. Finally, the paper concludes in Section 6 with a discussion of the key results and directions for future research.

2. The BAR(1) Model and ${\mathrm{WL}}_1$ Estimation

Consider the random variables $X_1,X_2,\cdots ,X_n$, which represent observations on a perfect binary tree with g generations. The initial observation, $X_1$, corresponds to generation 0, while the observations $X_{2^i},X_{2^i+1},\cdots ,X_{2^{i+1}-1}$ correspond to the $2^i$ observations in generation i, where $i=1,2,\cdots ,g$. The total sample size is given by $n=2^{g+1}-1$.

The first-order bifurcating autoregressive (BAR(1)) model is formulated as:

$$X_t={\mathbf{Y}}_t^{\prime }\varphi +{\epsilon }_t,\mathrm{for}\mathrm{all}t\ge 2,$$

(1)

where $X_t$ represents the observed value of a quantitative characteristic at time t, and ${\mathbf{Y}}_t^{\prime }=[1,X_{\lfloor t/2\rfloor }]$, with $X_{\lfloor t/2\rfloor }$ denoting the "mother" of $X_t$ for all $t\ge 2$, and $\lfloor u\rfloor$ representing the greatest integer less than or equal to u. The vector ${\varphi }^{\prime }=[{\varphi }_0,{\varphi }_1]$ contains the unknown model coefficients: the intercept, ${\varphi }_0$, and the autoregressive parameter, ${\varphi }_1$ (also referred to as the inherited effect or maternal correlation), where ${\varphi }_1\in (-1,1)$. This range for ${\varphi }_1$ ensures the stationary of the process.

The pairs $({\epsilon }_{2t},{\epsilon }_{2t+1})$ are supposed to be independently and identically distributed (i.i.d.) according to a joint distribution F. Each pair $({\epsilon }_{2t},{\epsilon }_{2t+1})$ has a zero mean vector and a variance-covariance matrix given by:

$${\Sigma }_{{\epsilon }_t}=\left(\begin{array}{cc}1& \theta \\ \theta & 1\end{array}\right){\sigma }^2,$$

where $\theta$ represents the linear correlation between ${\epsilon }_{2t}$ and ${\epsilon }_{2t+1}$, commonly referred to as the environmental effect or sister-sister correlation given their shared mother. Additionally, ${\sigma }^2$ represents the variance of the errors. It is further assumed that the pairs $({\epsilon }_{2t},{\epsilon }_{2t+1})$ and $({\epsilon }_{2s},{\epsilon }_{2s+1})$ are independent for $s\ne t$.

The underlying rationale for this correlation structure is rooted in the assumption that sister cells, particularly in the early stages of their development, are exposed to a common environment. Consequently, two distinct forms of correlation are expected: (i) environmental correlation between sisters, and (ii) maternal correlation deriving from inherited effects from the mother. In contrast, more distantly related cells, such as cousins, experience less environmental overlap, making it reasonable to assume that their environmental effects are independent.

[12] introduced the weighted L₁ estimators of the intercept and the autoregressive coefficient, ${\varphi }_0$ and ${\varphi }_1$, respectively, for the first-order BAR model by minimizing the following dispersion function

$$D\left(\varphi \right)=\sum _{2\le t\le n}w\left(X_{[t/2]}\right)\left|X_t-{\mathbf{Y}}_t^{\prime }\varphi \right|,$$

(2)

where $w(·)$ is a weight function chosen such that the product $w\left(X_{[t/2]}\right)X_{[t/2]}$ remains bounded, thereby providing robustness against outliers in the design points (i.e., outlying values of $X_{[t/2]}$). When $w(·)=1$, the estimator reduces to the standard ${\mathrm{L}}_1$-norm estimator. The dispersion function $D\left(\varphi \right)$ is non-negative, piecewise linear, and convex, provided that $w(·)>0$. These properties are sufficient to guarantee that $D\left(\varphi \right)$ attains a minimum, yielding the estimate ${\widehat{\varphi }}_n$.

2.1. Weighted and Non-Weighted Estimates for BAR(1) Model

In this section, we introduce two types of estimators for the BAR(1) model based on the weight function, $w(·)$, in (2). These estimators were studied by [12] in the context of bifurcating autoregressive (BAR) models.

2.1.1. Least Absolute Deviation (L₁) Estimator

The first estimator is the Least Absolute Deviation (L₁) estimator, which is a non-weighted version. This estimator is derived by setting the weight function in (2) to 1 for all x (i.e., $w\left(X_{[t/2]}\right)=1$). The L₁ estimator is known to be sensitive to bad leverage points in a multivariate autoregressive context [18]. [12] demonstrated that this sensitivity persists in the BAR(1) model when confronted with bad leverage points caused by additive outliers.

2.1.2. Weighted Least Absolute Deviation (WL₁) Estimators

The second type is the Weighted Least Absolute Deviation (${\mathrm{WL}}_1$) estimator, which incorporates a weighting scheme defined by:

$$w\left(X_{[t/2]}\right)=:w_t=min\left(1,{\left({\displaystyle \frac{k}{{\left(X_{[t/2]}-\right)}^{{\mu }^{\prime }}{\Sigma }^{-1}\left(X_{[t/2]}-\mu \right)}}\right)}^{\delta /2}\right),$$

(3)

The weight is determined based on the Mahalanobis distance between $X_{\lfloor t/2\rfloor }$ and a robust measure of location and dispersion, $(\mu ,\Sigma )$, utilizing the minimum covariance determinant (MCD) estimates as outlined in [19]. The constant k serves as a cutoff point, set at the 95th percentile of the Chi-square distribution with one degree of freedom, while the parameter $\delta$ is fixed at 2. The weighting scheme in (3) is similar to that used by [20].

This weighted estimator provides an advantage over the non-weighted L₁ estimator by being less sensitive to bad leverage points in both the response and factor spaces. For more detailed discussions, refer to [21] in the context of linear regression, [22] for autoregressive time series models, and [12] for bifurcating autoregressive models.

[12] derived the joint limiting distribution of the ${\mathrm{WL}}_1$ estimators of the BAR(1) model coefficients. Under some regularity conditions, they showed that

$$\sqrt{n}({\widehat{\varphi }}^{{\mathrm{WL}}_1}-\varphi )\stackrel{d}{\to }N(0,{\tau }^2{\mathbf{C}}^{-1}\Sigma {\mathbf{C}}^{-1}),$$

where $\tau ={\displaystyle \frac{1}{2}}f\left(0\right)$, $\mathbf{C}=E\left(w\left(X_{[t/2]}\right){\mathbf{Y}}_1{\mathbf{Y}}_1^{\prime }\right)$, $\Sigma =E\left(w\left(X_{[t/2]}\right){\mathbf{Y}}_1\Phi {\mathbf{Y}}_1^{\prime }w\left(X_{[t/2]}\right)\right)$, $\Phi =1+{\rho }_{sgn\left({ϵ}_{2t}\right),sgn\left({ϵ}_{2t+1}\right)}$, and ${\rho }_{sgn\left({ϵ}_{2t}\right),sgn\left({ϵ}_{2t+1}\right)}$ denotes the correlation coefficient between the signs of the pairs of errors $\left({ϵ}_{2t},{ϵ}_{2t+1}\right)$. The method of moments estimators for $\mathbf{C}$ and $\Sigma$ are given by

$$\widehat{\mathbf{C}}=\left(\begin{array}{cc}{\displaystyle \frac{1}{n-1}}\sum _{t=2}^{n}w\left(X_{[t/2]}\right)& {\displaystyle \frac{1}{n-1}}\sum _{t=2}^{n}w\left(X_{[t/2]}\right)X_{[t/2]}\\ {\displaystyle \frac{1}{n-1}}\sum _{t=2}^{n}w\left(X_{[t/2]}\right)X_{[t/2]}& {\displaystyle \frac{1}{n-1}}\sum _{t=2}^{n}w\left(X_{[t/2]}\right)X_{[t/2]}^2\end{array}\right),$$

$$\widehat{\Sigma }=\left(\begin{array}{cc}{\displaystyle \frac{1}{n-1}}\sum _{t=2}^n{w}^2\left(X_{[t/2]}\right)& {\displaystyle \frac{1}{n-1}}\sum _{t=2}^n{w}^2\left(X_{[t/2]}\right)X_{[t/2]}\\ {\displaystyle \frac{1}{n-1}}\sum _{t=2}^n{w}^2\left(X_{[t/2]}\right)X_{[t/2]}& {\displaystyle \frac{1}{n-1}}\sum _{t=2}^n{w}^2\left(X_{[t/2]}\right)X_{[t/2]}^2\end{array}\right),$$

and

$$\widehat{\Phi }=1+{\widehat{\rho }}_{sgn\left({ϵ}_{2t}\right),sgn\left({ϵ}_{2t+1}\right)},$$

where

$${\widehat{\rho }}_{sgn\left({ϵ}_{2t}\right),sgn\left({ϵ}_{2t+1}\right)}={\displaystyle \frac{1}{m}}\sum _{t=1}^{m}sgn\left({\widehat{ϵ}}_{2t}\right),sgn\left({\widehat{ϵ}}_{2t+1}\right),$$

with $m=(n-1)/2$ and ${\widehat{ϵ}}_t=(X_t-{\mathbf{Y}}_t^{\prime }\widehat{\varphi })$. [12] suggested to follow the estimation procedure proposed by [23] and [24] to estimate the parameter $\tau$:

$$\widehat{\tau }={\displaystyle \frac{1}{2}}\widehat{f}\left(0\right)={\displaystyle \frac{\left({\widehat{ϵ}}_{(\frac{n}{2}+\sqrt{n})}-{\widehat{ϵ}}_{(\frac{n}{2}-\sqrt{n})}\right)}{4}}.$$

3. Bias in the ${\mathrm{WL}}_1$ Estimators for the BAR(1) Model

In this section, we perform an empirical investigation using Monte Carlo simulations to assess the finite-sample bias of the ${\mathrm{WL}}_1$ estimator for the autoregressive coefficient ${\varphi }_1$, examining how it varies with the model parameters ${\varphi }_1$ and $\theta$. The simulations follow a specific setup. Perfect binary trees of sizes $n=31,63,127$, and 255 (corresponding to different numbers of generations, $g=4,5,6$, respectively) are generated for each combination of ${\varphi }_1=\pm (0.10,0.35,0.60,0.85)$ and $\theta =\pm (0.30,0.60,0.90)$, capturing a range of maternal and environmental effect correlation levels. The model intercept, ${\varphi }_0$, is fixed at 10 in all simulations. All generated trees are assumed to be stationary, with the initial observation $X_1$ randomly drawn from a larger simulated perfect binary tree of size 127.

Outliers are introduced by contaminating the errors through the use of two bivariate normal distributions. These distributions combine a standard bivariate normal distribution with contamination rates of 5% and 10%, drawn from an additional bivariate normal distribution with a mean of zero and a variance of 25. For each simulation scenario, $N=10,000$ binary trees are generated, and the ${\mathrm{WL}}_1$ estimator ${\widehat{\varphi }}_1^{\mathrm{WL}}$ is computed for each tree. Subsequently, we assess the empirical bias, absolute relative bias (ARB), variance, and mean squared error (MSE) of ${\widehat{\varphi }}_1^{\mathrm{WL}}$ as follows:

$$\widehat{\mathrm{BIAS}}\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1}\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1,\left(j\right)}-{\varphi }_1\right),$$

(4)

$$\widehat{\mathrm{ARB}}\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1}\right)=|\widehat{\mathrm{BIAS}}\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1}\right)/{\varphi }_1|,$$

(5)

where ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1,\left(j\right)}$ is the ${\mathrm{WL}}_1$ estimator from iteration $j=1,2,\cdots ,N$.

Figure 2 presents the empirical bias and ARB ($\times 100\%$) of ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$. The results indicate that the bias of ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$ can be significant across various combinations of ${\varphi }_1$ and $\theta$, particularly for smaller sample sizes. Specifically, it is evident that the volume of the bias is concerning when $\theta$ is positive and ${\varphi }_1$ exceeds $-0.5$. In general, ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$ tends to underestimate ${\varphi }_1$ for values of ${\varphi }_1$ greater than $-0.5$, with the degree of underestimation increasing as ${\varphi }_1$ approaches $+1$. Furthermore, the bias appears to exhibit a linear relationship with ${\varphi }_1$.

The extent of the bias is further emphasized by the ARB results, which show that for small sample sizes ($n=31$) and ${\varphi }_1$ near zero (in the range $-0.3$ to $+0.3$), the bias of ${\widehat{\varphi }}^{{\mathrm{WL}}_1}$ can range from 5% to 60% of the true value of ${\varphi }_1$ as $\theta$ varies from $-1$ to $+1$. For the smallest sample size ($n=31$), the relative bias remains considerably large for values of ${\varphi }_1$ greater than $-0.5$, provided $\theta$ does not approach $-1$. While the bias remains considerable in most cases for $n=63$ and $n=127$, it becomes negligible for larger sample sizes ($n=255$), except when ${\varphi }_1$ is near zero (in the range $-0.3$ to $+0.3$) and $\theta$ is positive.

4. Bias-Corrected ${\mathrm{WL}}_1$ Estimators for the BAR(1) Model

In this section, we outline three bootstrap methods designed to correct the bias of the ${\mathrm{WL}}_1$ estimators, ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$, for the autoregressive coefficient ${\varphi }_1$ in the BAR(1) model. Specifically, we introduce the single bootstrap bias-corrected version of ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$, as well as a fast double bootstrap bias-corrected version.

The model-based bootstrap technique has been widely used to correct bias in the least squares (LS) estimation of the AR(1) model, with its properties thoroughly explored in the literature (see [25,26,27,28], among others). This approach involves generating a large number of bootstrap samples based on the estimated model parameters and resampled residuals. These samples are then used to compute multiple replicates of the estimator, which allows for the estimation of the sampling distribution, variance, and bias of the estimator. Recently, [13] extended this method to correct bias in the LS estimation of the BAR(1) model coefficients. In this section, we focus on describing two model-based bootstrap bias correction methods specifically designed for the BAR(1) model.

We begin by rewriting the BAR(1) model, originally presented in Eq. (1), as follows:

$$\begin{array}{cc}\hfill X_{2t}=& {\varphi }_0+{\varphi }_1{X}_t+{\epsilon }_{2t},\hfill \\ \hfill X_{2t+1}=& {\varphi }_0+{\varphi }_1{X}_t+{\epsilon }_{2t+1},\hfill \end{array}$$

for $t\ge 1$. Here, $X_{2t}$ and $X_{2t+1}$ represent the observations for the two daughter nodes branching from $X_t$.

The single bootstrap bias correction procedure for the autoregressive coefficient in the BAR(1) model is described in Algorithm 1. It is important to note that the bias in the corrected estimate, ${\widehat{\varphi }}_1^{\mathrm{SBC}}$, is of the order $O\left(n^{-2}\right)$, whereas the original estimator, ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$, shows bias of the order $O\left(n^{-1}\right)$.

The double bootstrap technique involves performing the single bootstrap procedure twice. This method has been shown to enhance the accuracy of bootstrap bias correction in the AR(1) model [e.g., [29,30,31,32], and similarly for the BAR(1) model [see, [13]. However, a notable drawback of the standard double bootstrap method is its high computational cost. To address this issue, a modified double bootstrap bias correction algorithm, referred to as the fast double bootstrap, was proposed by [33]. This algorithm offers greater computational efficiency compared to the traditional double bootstrap, as it requires only one bootstrap resample in phase 2 (i.e., $B_2=1$) for each bootstrap sample generated in phase 1, in contrast to the conventional double bootstrap, which generates $B_2$ resamples within each phase 1 bootstrap sample. Consequently, while the computational complexity of the standard double bootstrap algorithm is of the order $O\left(B_1{B}_2\right)$, the fast double bootstrap reduces this to $O\left(B_1\right)$. Despite this significant reduction in computational cost, [13] observed that the bias correction performance of both algorithms remains comparable. The fast double bootstrap approach for correcting the bias in the ${\mathrm{WL}}_1$ estimator of ${\varphi }_1$ under the BAR(1) model is outlined in Algorithm 2.

Algorithm 1: Single Bootstrap Bias-Corrected ${\mathrm{WL}}_1$ Estimation for ${\varphi }_1$ (SBC).

Input Observed tree: ${\mathbf{X}}^n=(X_1,X_2,\cdots ,X_n)$ ${\mathrm{WL}}_1$ estimates of the BAR(1) model coefficients: ${\widehat{\varphi }}_0^{{\mathrm{WL}}_1}$ and ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$ Centered residuals: ${\tilde{\epsilon }}_t={\widehat{\epsilon }}_t-{(n-1)}^{-1}{\sum }_{i=2}^n{\widehat{\epsilon }}_i$ for $t\ge 2$ Number of bootstrap resamples: B 1: for each $b\leftarrow 1$ to B do 2:     Set $X_{1,b}^{*}=X^0$, the final observation from an initial binary tree § of size $n^0=31\S \S$, where $X_1$ is the first observation in the original tree 3:     for each $j\leftarrow 1$ to $m=(n-1)/2$ do 4:         Sample with replacement a pair $({\tilde{\epsilon }}_{2j,b}^{*},{\tilde{\epsilon }}_{2j+1,b}^{*})$ from the set of pairs $\{({\tilde{\epsilon }}_{2t},{\tilde{\epsilon }}_{2t+1});t\ge 1\}$ 5:         Compute $$X_{2j,b}^{*}={\widehat{\varphi }}_0^{{\mathrm{WL}}_1}+{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}{X}_{j,b}^{*}+{\tilde{\epsilon }}_{2j,b}^{*},$$ $$X_{2j+1,b}^{*}={\widehat{\varphi }}_0^{{\mathrm{WL}}_1}+{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}{X}_{j,b}^{*}+{\tilde{\epsilon }}_{2j+1,b}^{*}$$ 6:     end for 7:     Construct the bootstrap tree ${\mathbf{X}}_b^{*n}=(X_{1,b}^{*},X_{2,b}^{*},\cdots ,X_{n,b}^{*})$ 8:     Compute the ${\mathrm{WL}}_1$ estimate ${\widehat{\varphi }}_{1_b}^{*}$ from ${\mathbf{X}}_b^{*n}$ 9: end for 10: Estimate the bias of ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$ as ${\widehat{\beta }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}={\displaystyle \frac{1}{B}}\sum _{b=1}^B({\widehat{\varphi }}_{1_b}^{*}-{\widehat{\varphi }}_1^{{\mathrm{WL}}_1})$ Output:  The single bootstrap bias-corrected ${\mathrm{WL}}_1$ estimate of ${\varphi }_1$ is given by: $${\widehat{\varphi }}_1^{\mathrm{SBC}}={\widehat{\varphi }}_1^{{\mathrm{WL}}_1}-{\widehat{\beta }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}$$

§In each bootstrap iteration, the initial tree is generated from a BAR(1) model with coefficients set to the observed ${\mathrm{WL}}_1$ estimates, ${\widehat{\varphi }}_0^{{\mathrm{WL}}_1}$ and ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$, while the errors are sampled in pairs from the centered residuals. It is important to note that using $X_1$ from the observed tree as the initial observation in all bootstrap trees could introduce artificial correlation among the bootstrap samples. This issue may become more significant in higher-order BAR(p) models, where $2^p-1$ initial observations are required. In such cases, our method utilizes the final $2^p-1$ observations from the initial tree to avoid this potential problem.

§§Any appropriate initial tree size, $n^0$, can be chosen. The suggested size of $n^0=31$ provides a reasonable compromise between computational efficiency and result stability.

Algorithm 2: Fast Double Bootstrap Bias-Corrected ${\mathrm{WL}}_1$ Estimation for ${\varphi }_1$ (FDBC).

Input: Observed tree: ${\mathbf{X}}^n=(X_1,X_2,\dots ,X_n)$ ${\mathrm{WL}}_1$ estimates of the BAR(1) model coefficients: ${\widehat{\varphi }}_0^{{\mathrm{WL}}_1}$ and ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$ Centered residuals: ${\tilde{\epsilon }}_t={\widehat{\epsilon }}_t-{(n-1)}^{-1}{\sum }_{i=2}^n{\widehat{\epsilon }}_i$, for $t\ge 2$ Number of phase 1 bootstrap resamples: $B_1$ Number of phase 2 bootstrap resamples: $B_2$ 1: for each $b\leftarrow 1$ to $B_1$ do 2:     Set $X_{1,b}^{*}=X^0$, where $X^0$ is the last observation in an initial binary tree of size $n^0=31$, with $X_1$ as the first observation in ${\mathbf{X}}^n$ 3:     for each $j\leftarrow 1$ to $m=(n-1)/2$ do 4:         Sample with replacement a pair $({\tilde{\epsilon }}_{2j,b}^{*},{\tilde{\epsilon }}_{2j+1,b}^{*})$ from the set $\{({\tilde{\epsilon }}_{2t},{\tilde{\epsilon }}_{2t+1})\mid t\ge 1\}$ 5:         Compute: $$X_{2j,b}^{*}={\widehat{\varphi }}_0^{{\mathrm{WL}}_1}+{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}{X}_{j,b}^{*}+{\tilde{\epsilon }}_{2j,b}^{*},$$ $$X_{2j+1,b}^{*}={\widehat{\varphi }}_0^{{\mathrm{WL}}_1}+{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}{X}_{j,b}^{*}+{\tilde{\epsilon }}_{2j+1,b}^{*}$$ 6:     end for 7:     Construct the bootstrap tree ${\mathbf{X}}_b^{*n}=(X_{1,b}^{*},X_{2,b}^{*},\dots ,X_{n,b}^{*})$ 8:     Compute the first-phase bootstrap ${\mathrm{WL}}_1$ estimates ${\widehat{\varphi }}_{0_b}^{*}$ and ${\widehat{\varphi }}_{1_b}^{*}$ from ${\mathbf{X}}_b^{*n}$ 9:     Obtain the first-phase bootstrap residuals: $${\widehat{\epsilon }}_{2t,b}^{*}=X_{2t,b}^{*}-({\widehat{\varphi }}_{0_b}^{*}+{\widehat{\varphi }}_{1_b}^{*}{X}_{t,b}^{*}),$$ $${\widehat{\epsilon }}_{2t+1,b}^{*}=X_{2t+1,b}^{*}-({\widehat{\varphi }}_{0_b}^{*}+{\widehat{\varphi }}_{1_b}^{*}{X}_{t,b}^{*}),t\ge 1$$ 10:     Compute centered bootstrap residuals: $${\tilde{\tilde{\epsilon }}}_{t,b}^{*}={\widehat{\epsilon }}_{t,b}^{*}-{(n-1)}^{-1}\sum _{i=2}^n{\widehat{\epsilon }}_{i,b}^{*},t\ge 2$$ 11:     for each $k\leftarrow 1$ to $B_2=1$ do 12:         Apply steps 2-6 to ${\mathbf{X}}_b^{*n}$ using ${\widehat{\varphi }}_{0_b}^{*}$, ${\widehat{\varphi }}_{1_b}^{*}$, and the centered residuals ${\tilde{\tilde{\epsilon }}}_{t,b}^{*}$ 13:         Construct the second-phase bootstrap tree ${\mathbf{X}}_{k,b}^{**n}=(X_{1,kb}^{**},X_{2,kb}^{**},\dots ,X_{n,kb}^{**})$ 14:         Compute the second-phase bootstrap ${\mathrm{WL}}_1$ estimate ${\widehat{\varphi }}_{1_{kb}}^{**}$ from ${\mathbf{X}}_{k,b}^{**n}$ 15:     end for 16: end for 17: Estimate the single bootstrap bias of ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1}$: $${\widehat{\beta }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}={\displaystyle \frac{1}{B_1}}\sum _{b=1}^{B_1}({\widehat{\varphi }}_{1_b}^{*}-{\widehat{\varphi }}_1^{{\mathrm{WL}}_1})$$ 18: Compute the double bootstrap bias adjustment factor: $${\widehat{\gamma }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}={\widehat{\beta }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}-{\displaystyle \frac{1}{B_1{B}_2}}\sum _{b=1}^{B_1}\sum _{k=1}^{B_2}({\widehat{\varphi }}_{1_{kb}}^{**}-{\widehat{\varphi }}_{1_b}^{*})$$ Output:  The fast double bootstrap bias-corrected ${\mathrm{WL}}_1$ estimate for ${\varphi }_1$: $${\widehat{\varphi }}_1^{\mathrm{FDBC}}={\widehat{\varphi }}_1^{{\mathrm{WL}}_1}-{\widehat{\beta }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}-{\widehat{\gamma }}_{{\widehat{\varphi }}_1^{{\mathrm{WL}}_1}}$$

5. Empirical Results

In this section, we present the findings of an empirical investigation aimed at (i) evaluating the effectiveness of the proposed bias correction methods for the ${\mathrm{WL}}_1$ estimator of the autoregressive coefficient in the BAR(1) model, and (ii) comparing the relative performance of these bias correction techniques. The study involved extensive simulations. All computations were performed using R version 4.1.3 [34] on a machine with an Intel Xeon E5-2699Av4 processor, featuring a clock speed of 2.4 GHz, 44 cores, and 512 GB of memory. The simulations utilized 42 of the 44 available cores through the `parallel` package, while bifurcating trees were generated using the `bifurcatingr` package [35].

The "core" process follows the model described in , adhering to the same variance-covariance matrix. The "observed" process is defined as:

$$X_t^{*}=X_t+V_t,t\ge 1.$$

(6)

In this definition, when $V_t=0$ for all t, the "observed" process aligns with the "core" process, and any deviations are attributed to the distribution of the errors. As noted by [36], these deviations are referred to as Innovation Outliers (IO) model. In this, the distribution of errors is typically characterized by heavy tails. On the other hand, when $V_t\ne 0$, the process corresponds to the Additive Outliers (AO) model [36]. In this case, $V_t$ follows a mixed distribution given by $(1-\gamma ){\delta }_0(·)+\gamma G(·)$, where $\gamma$ represents the contamination proportion, ${\delta }_0$ denotes a point mass at zero, and G is the contaminating distribution function. For our simulations, G was selected as a $N({\mu }_2,{\sigma }_c^2)$ distribution. Therefore, the distribution of $V_t$ can be expressed as:

$$BCN(\gamma ,{\mu }_1,{\Sigma }_{{\epsilon }_t},{\mu }_2,{\Sigma }_{{\sigma }_c^2})=(1-\gamma )BN({\mu }_1,{\Sigma }_{{\epsilon }_t})+\gamma BN({\mu }_2,{\Sigma }_{{\sigma }_c^2}),$$

(7)

where $BCN(·)$ denotes a Bivariate Contaminated Normal distribution with a specified mean vector and variance-covariance matrix.

The variable $V_t$ is generated as a sequence of independent and identically distributed random variables. In the AO model, outliers are introduced through the incorporation of external influences into the core process, resulting in outliers that arise from factors beyond the error distribution.

In the context of our study, an outlier introduced into the design is referred to as a leverage point. We differentiate between two types of leverage points: “good” and “bad.” As demonstrated by [19], a leverage point is classified as “good” if it results in a small residual despite being an outlier in the design space. In contrast, it is considered “bad” if it results in a large residual. According to [19], the IO model tends to produce “good” leverage points with relatively minimal impact on estimates. In contrast, the AO model generates “bad” leverage points that can significantly affect estimates, even those that are robust [12].

5.1. Simulation Setup

The simulations contain a range of scenarios defined by varying sample sizes n, maternal correlation levels ${\varphi }_1$, and error correlations $\theta$, as outlined in Section 3. For each scenario, we generate $N=10,000$ trees and compute the ${\mathrm{WL}}_1$ estimator for the coefficient ${\varphi }_1$ from each tree. Specifically, we evaluate the ${\mathrm{WL}}_1$ estimator (${\widehat{\varphi }}^{{\mathrm{WL}}_1}$), the single bootstrap bias-corrected estimator (${\widehat{\varphi }}_1^{\mathrm{SBC}}$), and the fast double bootstrap bias-corrected estimator (${\widehat{\varphi }}_1^{\mathrm{FDBC}}$).

We excluded the case of $n=225$ from this simulation study, as the bias in ${\widehat{\varphi }}^{{\mathrm{WL}}_1}$ for this larger sample size was found to be negligible, as detailed in Section 3. This exclusion was implemented to facilitate computation time.

For both the ${\widehat{\varphi }}_1^{\mathrm{SBC}}$ and ${\widehat{\varphi }}_1^{\mathrm{FDBC}}$ estimators, the number of bootstrap replicates was set to $B=999$. The empirical bias and absolute relative bias (ARB) for each estimator were computed as outlined in Eqs. (4) and (5). Furthermore, the root mean squared error (RMSE) and mean absolute deviation (MAD) were calculated as follows:

$$\mathrm{RMSE}\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1}\right)=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{j=1}^N{({\widehat{\varphi }}_1^{{\mathrm{WL}}_1,\left(j\right)}-{\varphi }_1)}^2},$$

$$\mathrm{MAD}\left({\widehat{\varphi }}_1^{{\mathrm{WL}}_1}\right)={\displaystyle \frac{1}{N}}\sum _{j=1}^N\left|{\widehat{\varphi }}_1^{{\mathrm{WL}}_1,\left(j\right)}-{\varphi }_1\right|.$$

where ${\widehat{\varphi }}_1^{{\mathrm{WL}}_1,\left(j\right)}$ is the ${\mathrm{WL}}_1$ estimator from iteration $j=1,2,\cdots ,N$.

5.1.1. Innovation Outliers (IO) Model Setup

Under the Innovation Outlier (IO) model, where $V_t=0$, the “observed” process simplifies to the “core” process as specified in Eq. (1). The error distributions examined in this study include the bivariate standard normal distribution with correlation $\theta$, as well as four variants of the bivariate contaminated normal distributions. The contaminated distributions are characterized by different levels of contamination, specifically (1) Contaminating variance of 25 with a contamination rate of 5%, (2) Contaminating variance of 25 with a contamination rate of 10%, (3) Contaminating variance of 49 with a contamination rate of 5%, (4) Contaminating variance of 49 with a contamination rate of 10%. These distributions are denoted as $\mathrm{BCN}(5\%,0,{\Sigma }_1,0,{\Sigma }_{25})$, $\mathrm{BCN}(10\%,0,{\Sigma }_1,0,{\Sigma }_{25})$, $\mathrm{BCN}(5\%,0,{\Sigma }_1,0,{\Sigma }_{49})$, $\mathrm{BCN}(10\%,0,{\Sigma }_1,0,{\Sigma }_{49})$, respectively. The chosen parameter settings for these bivariate contaminated normal distributions are intended to represent both “mild” and “extreme” contamination scenarios.

5.1.2. Additive Outliers (AO) Model Setup

Under the Additive Outlier (AO) model, the errors $\left({\epsilon }_{2t},{\epsilon }_{2t+1}\right)$ follow an innovation process modeled by a bivariate normal distribution with zero mean, unit variance, and $\theta =0$. The distribution of $V_t$ is given by:

$$V_t\sim (1-\gamma ){\delta }_0(·)+\gamma G(·),$$

where $\gamma$ represents the proportion of contamination, ${\delta }_0$ is a point mass at zero, and G is a contaminating distribution function. The “observed” process is defined as:

$$X_t^{*}=X_t+V_t={\mathbf{Y}}_t^{\prime }\varphi +{\epsilon }_t+V_t,\mathrm{for}\mathrm{all}t\ge 1,$$

where ${\mathbf{Y}}_t^{\prime }\varphi$ represents the model component, ${\epsilon }_t$ is the error term, and $V_t$ represents the additive outlier.

The distributions considered for $V_t$ include the bivariate standard normal distribution with correlation $\theta$, as well as four bivariate contaminated normal distributions. These contaminated distributions have contaminating variances of 25 and 49, with contamination rates of 5% and 10%. Specifically, the distributions are (1) $\mathrm{BCN}(5\%,0,{\Sigma }_1,0,{\Sigma }_{25})$, (2) $\mathrm{BCN}(10\%,0,{\Sigma }_1,0,{\Sigma }_{25})$, (3) $\mathrm{BCN}(5\%,0,{\Sigma }_1,0,{\Sigma }_{49})$, (4) $\mathrm{BCN}(10\%,0,{\Sigma }_1,0,{\Sigma }_{49})$. In addition, It is assumed that $V_t$ is independent of ${\epsilon }_t$.

5.2. Simulation Results

In this section, a results summary of the simulation study is presented in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. Our analysis focuses on the absolute relative bias and root mean squared error (RMSE) for the ${\mathrm{WL}}_1$ estimator of ${\varphi }_1$, including both the single bootstrap and fast double bootstrap bias-corrected estimators. The key findings are outlined as follows:

• Generally, all proposed bias-correcting estimators substantially reduce the bias of the ${\mathrm{WL}}_1$ estimator across nearly all combinations of ${\varphi }_1$ and $\theta$, even with small sample sizes.
• Both the single bootstrap (${\widehat{\varphi }}_1^{\mathrm{SBC}}$) and fast double bootstrap (${\widehat{\varphi }}_1^{\mathrm{FDBC}}$) bias-corrected estimators effectively reduce bias. However, the fast double bootstrap estimator demonstrates superior performance compared to the single bootstrap estimator, particularly for small sample sizes (e.g., $n=31$). Despite this reduction, some residual bias persists when both $\theta$ and ${\varphi }_1$ are positive, with the bias increasing as $\theta$ and ${\varphi }_1$ approach the positive boundary.
• For larger sample sizes, $n=63$ and $n=127$, the fast double bootstrap estimator generally outperforms the single bootstrap estimator, though their performances are sometimes comparable. It is important to note that the fast double bootstrap method requires more computational time compared to the single bootstrap method.
• Both bootstrap methods lead to significant reductions in RMSE and mean absolute deviation (MAD) compared to the ${\mathrm{WL}}_1$ estimator across all sample sizes and positive values of $\theta$.
• The bias reduction achieved by both bootstrap estimators is particularly notable for ${\varphi }_1$ values near the boundaries (-1, +1), especially for small sample sizes (e.g., $n=31$). Additionally, the fast double bootstrap method outperforms the single bootstrap method in mitigating bias in the ${\mathrm{WL}}_1$ estimators.

Overall, these findings indicate that both the single and fast double bootstrap bias-corrected estimators are suitable for reducing bias in the ${\mathrm{WL}}_1$ estimation for the BAR(1) model under most combinations of ${\varphi }_1$ and $\theta$. Although the fast double bootstrap method is more effective in reducing bias, it requires more computational resources. Consequently, the single bootstrap approach may be preferred for its computational efficiency.

6. Conclusions

In this study, we revisited the weighted L₁ (WL₁) estimator within the framework of first-order bifurcating autoregressive BAR(1) models. While the WL₁ estimator is known for its robustness to outliers, it can show substantial bias for small and moderate sample sizes, akin to the least squares estimator. The extent and direction of this bias are influenced by the autoregressive coefficient (${\varphi }_1$) and the correlation between model errors ($\theta$).

To address this issue, we introduced bias-corrected versions of the WL₁ estimator using the bootstrap method. Specifically, we developed the single bootstrap and fast double bootstrap bias-corrected WL₁ estimators. Our extensive empirical analysis demonstrated that these bootstrap bias-corrected estimators effectively reduce the bias inherent in the WL₁ estimator. The observed reduction in bias corresponded with decreased root mean squared error (RMSE), thereby enhancing the overall efficiency of the corrected estimators. Although the fast double bootstrap method is somewhat more computationally demanding than the single bootstrap, it consistently outperformed the latter in most of the scenarios considered in our simulations.

While our research primarily focused on BAR(1) models, future investigations will extend to higher-order BAR models (i.e., BAR($p>1$)). Additionally, we plan to explore how bias and its correction impact related inferences, such as confidence intervals and hypothesis tests, to further assess the practical implications of our findings.