Rank-based regression with bootstrapping and the least square regression for analysis of small sample interrupted time-series data: simulation studies and application to illegal organ transplants

Abstract

For interrupted time series with small samples and outliers, the standard least squares regression for estimating intervention effects can be biased when the assumption of normality is violated. This article proposes a rank-based regression method combined with bootstrapping to address outliers and non-normality. The proposed method minimizes the Jaeckel dispersion function based on ranked residuals, offering robustness against outliers, skewed error distributions, and correlated errors. Simulation studies and a real-world example compare the performance of the proposed method with the least squares method under various conditions. Both produce unbiased estimates for the intervention effects under normality. However, with non-normal errors, outliers, and correlated errors, the rank-based bootstrapping method outperforms least squares producing narrower confidence intervals, higher efficiency, lower type I error rates, and greater power. Application to the Istanbul declaration on illegal organ activities confirms the improved robustness and precision of the rank-based method with bootstrapping over least squares.

KEYWORDS:

• Interrupted time-series analysis
• Jaeckel dispersion function
• simulations
• power of the test
• mean square error (MSE)
• type I error
• bootstrapping

1. Introduction

In various fields, such as public health, clinical practice, economics, political sciences and road safety analysis, it is often necessary to assess the success of an intervention or a program [1–5]. Evaluating policy interventions involves analyzing interrupted time series data (ITS) – data collected before and after the interventions. In public health, medicine and road safety analysis, ITS analysis is commonly used to evaluate the effectiveness of policy interventions [1–4,6].

Various statistical approaches are used to measure these effects in ITS [2,3,5], but the two most frequently used approaches are either regression-based methods, such as standard regression and generalized least squares, or ARIMA (autoregressive, integrated, moving average) models [5]. Among the regression methods, segmented regression analysis is the most widely used [2]. For a detailed description of segmented regression analysis, see [2–4].

Studies often compared segmented regression models with other models. Wang et al. [3] compared the segmented regression model to both a log linear binomial regression model, and a full Bayesian regression model to estimate the effect of the mandatory helmet law on the number of head injuries. Turner et al. [7] compared six different commonly used statistical models on 180 published data series in the public health area to see how the choice of the statistical model impacted statistical significance and estimation. Models examined in the Turner et al. [7] study included ordinary least squares (OLS) without adjustment for autocorrelation, OLS with Newey-West standard errors, Prais-Winsten generalized least squares, restricted maximum likelihood and autoregressive integrated moving average (ARIMA).

Through simulation, Turner et al. [8] also investigated the effect of autocorrelation, one of the important features of interrupted time series data [9], on inferences for the intervention effect. Although the impact of autocorrelation on standard regression individual parameter tests is minimal for small sample ITS data [10,11], methods that do not account for the autocorrelation provide less-reliable standard errors of the estimates, specifically downwardly biased standard errors [10]. Turner et al. [8] showed that autocorrelation is difficult to detect in small samples, so even if a formula-based method is employed to adjust for autocorrelation (OLS with Newey-West standard errors, Prais-Winsten generalized least squares) it won’t be effective if there is not an accurate estimate of the true autocorrelation. Short time series are common in public health, ecological and environmental studies [12], and Turner et al. [13] also showed that public health ITS data often have small samples.

The accuracy of segmented regression depends on the absence of outliers in the data, and the normality of the error distribution [1–3,14,15]. When the error distribution deviates from normality or outliers are present in the data, confidence intervals and hypothesis testing for the intervention effects become inefficient and biased using a segmented regression model [14,16]. However, with large samples, normality assumptions are robust, and valid results can be obtained from the segmented regression model using the normal approximation [14]. But the normality assumption is not robust in small sample studies, and the normal approximation is not appropriate.

To address the influence of outliers, various robust regression techniques are used [17–20], including the least absolute deviation (LAD) method, Huber’s M-estimator, least median square, least trimmed square estimators (LTS), S-estimator, MM-estimators [21] and Rank-based robust regression [22–28]. Robust regression analysis in ITS is limited [7] and rank-based regression using the Jaeckel dispersion score function, which is highly robust to outliers and distributional anomalies [5,16,24–27,29–37], is currently in limited use in ITS [7].

In this article, we propose and explore a rank-based regression with bootstrapping to address the problems of small sample sizes, outliers in the data, non-normality in the error distribution and any potential autocorrelation in ITS data analysis. Through simulation, we investigated the effectiveness of our proposed statistical procedure on the estimation of the intervention effects, their confidence intervals and type I errors when there are deviations from the normality assumption, autocorrelated errors are present and outliers exist in the data. We also analyzed a real dataset on illegal organ transplants to examine the effect of the Istanbul Declaration on reducing such transplants using both the LS and RB with bootstrapping methods.

The article is organized as follows. Section 2 introduces the estimation and hypothesis procedures for intervention effects in interrupted time-series data using both the least square (LS) method and the method of rank-based (RB) regression with bootstrapping. In Section 3, comprehensive simulation studies are conducted under various conditions to compare the performances of the LS method and the RB with bootstrapping method. The simulations assessed the estimates of the intervention effects, the efficiency of the estimates using mean square error, type I error and power. In Section 4, the LS and RB methods are used for the analysis of the real dataset on illegal organ transplants to evaluate the effectiveness of the 2008 Istanbul Declaration in reducing illegal organ trafficking. Finally, Section 5 presents discussions, practical implications of the study and conclusions. All analyses were conducted using R software version 4.3.0, including the Rfit package [36] for estimating the intervention effects using rank-based regression.

2. Methods

2.1. Estimation of intervention effects

The segmented regression model for interrupted time-series data is defined as (1) $Y_t={\beta }_0+{\beta }_1{T}_t+{\beta }_2{X}_t+{\beta }_3{X}_t\left(T_t-\frac{n}{2}\right)+{\epsilon }_t,t=1,2,\dots ,n.$(1) In this model, $Y_t$ represents the outcome variable measured at equally spaced time points $t$. $T_t$ denotes the time elapsed since the start of the study. The variable $X_t$ is a dummy variable that takes a value of 0 during pre-intervention periods, and 1 otherwise, indicating the presence of the intervention. The term $X_t\left(T_t-\frac{n}{2}\right)$ represents the interaction between $X_t$ and $T_t$, with $T_t$ centered to prevent multicollinearity.

The validity of the model depends on some assumptions, such as (i) there must be time points before and after an intervention, with evenly spaced intervals between time points, (ii) the relationship between the outcome and time is linear within pre-intervention and post-intervention segments, (iii) the errors ${\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n$ are independent and identically normally distributed with mean zero and constant variance, (iv) there are no outliers or missing data points, (v) the intervention is expected to immediately change the level or slope, not gradually over time.

Model (1) can be decomposed into two parts: $Y_t={\beta }_0+{\beta }_1{T}_t+{\epsilon }_t$ for the pre-intervention outcomes, and $Y_t=({\beta }_0+{\beta }_2-{\beta }_3\frac{n}{2})+({\beta }_1+{\beta }_3)T_t+{\epsilon }_t$ for the outcomes after the intervention. Consequently, ${\beta }_0$ represents the initial level of the outcome variable and ${\beta }_1$ corresponds to the slope or trend of the outcome variable before the intervention. The term ${\beta }_2-{\beta }_3\frac{n}{2}$ denotes the change in the intercept between the pre-and post-intervention models, while ${\beta }_3$ represents the difference in slopes or trends between the pre- and post-intervention periods. Generally, ${\beta }_2$ does not possess a direct interpretation but contributes to adjusting the intercept in the post-intervention model. However, in the special case where ${\beta }_3=0$ but ${\beta }_2\ne 0$, ${\beta }_2$ represents a flat and steady increase or decrease in $Y_t$, in addition to the constant trend.

For convenience, we can express model (1) in matrix notation as follows: $Y={\beta }_{0}1+\mathrm{X\beta }+\epsilon ,$where $Y=(y_1,y_2,\dots ,y_n{)}^T$ is the $n\times 1$ vector of responses, $X$ is a $n\times 3$ design matrix, $\beta =({\beta }_1,{\beta }_2,{\beta }_3{)}^T$, and $\epsilon =({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_n{)}^T$ is a $n\times 1$ vector of error terms.

The least squares (LS) estimate of $\beta$ is obtained by minimizing the sum of squared errors defined as $L=\sum _{i=1}^n{\epsilon }_i^2=\sum _{i=1}^n(y_i-{\beta }_0-x_i^{\mathrm{\prime }}\beta {)}^2.$The least squares estimate of $\beta$ is given by $\hat{\beta }=(X^{\mathrm{\prime }}X{)}^{-1}{X}^{\mathrm{\prime }}y,$where ${\hat{\beta }}^T=({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3)$. This method performs well when the data contains no outliers, and the errors follow a normal distribution.

Alternatively, for the rank-based estimates of the parameter vector $\beta ,$ the Jaeckel dispersion score function [16] is minimized using Wilcoxon scores. The dispersion score function, denoted as $D(\beta )$, is defined as: (2) $D(\beta )=\sqrt{\{12\}{(n+1)}^{\frac{1}{2}}\sum _{i=1}^n\left[R(y_i-x_i^{\mathrm{\prime }}\beta )-\frac{n+1}{2}\right](y_i-x_i^{\mathrm{\prime }}\beta )},$(2) where $R(.)$ represents the rank of $y_i-x_i^{\mathrm{\prime }}\beta$ among the residuals, $y_1-x_1^{\mathrm{\prime }}\beta ,y_2-x_2^{\mathrm{\prime }}\beta ,\dots ,y_n-x_n^{\mathrm{\prime }}\beta$. The dispersion function $D(\beta )$ is a non-negative, piece-wise linear, continuous convex function of $\beta$. Minimizing this dispersion function provides an estimate of $\beta$, which generally does not have a closed-form expression. Therefore, iterative computer methods are typically required to obtain numerical solutions.

The rank-based estimation approach offers robustness to outliers and deviations from normality in the error distribution. By minimizing a linear function of the residuals with rank-based coefficients, the effects of outliers are incorporated in a linear rather than quadratic manner. This robustness makes the rank-based estimates more reliable compared to the least squares (LS) method.

It is important to note that the Jaeckel dispersion function [16] is invariant to the intercept term ${\beta }_0$. Consequently, ${\beta }_0$ cannot be simultaneously estimated with the other elements of $\beta .$Hattmansperger and McKean [32] proposed using the full-model residuals, denoted as ${\hat{\epsilon }}_i=y_i-x_i^{\mathrm{\prime }}\hat{\beta }$, for $i=1,2,\dots ,n,$ to estimate ${\beta }_0$. The estimator for ${\beta }_0$ is obtained as the median of the residuals, denoted as ${\hat{\beta }}_0=\text{median}\{{\epsilon }_i,i=1,2,\dots ,n\}$.

2.2. Inference

The parameters ${\beta }_2$ and ${\beta }_3$ in model (1) quantify the intervention effect. If there is no intervention effect, their values will be zero. Thus, the hypothesis of interest can be stated as $H_0:{\beta }_2={\beta }_3=0$ versus $H_a:$ at least one of ${\beta }_2$ and ${\beta }_3$ is not zero.

Under the null hypothesis, model (1) reduces to a simplified form: (3) $Y_t={\beta }_0+{\beta }_1{T}_t+{\epsilon }_t,t=1,2,\dots ,n.$(3)

Let ${\hat{\beta }}_F$ and ${\hat{\beta }}_R$ be the value of $\beta$ that minimizes the Jaeckel dispersion score function $D(\beta )$ in (2) under the full model (1) and the reduced model (3), respectively. Consequently, $D({\hat{\beta }}_F)$ and $D({\hat{\beta }}_R)$ represent the overall minimum of the Jaeckel dispersion measure under the full model and the reduced model, respectively. We define $D^{\ast }=D({\hat{\beta }}_R)-D({\hat{\beta }}_F)$, as the drop or reduction in Jaeckel dispersion when fitting the full model (intervention effect) rather than the reduced model (no intervention effect) under the null hypothesis.

The Jaeckel-Hettmansperger-McKean test statistic [16,32] is $\mathrm{HM}=\frac{D^{\ast }/2}{\hat{\tau }/2}\sim F(2,n-3-1),$where $\tau$ is the scale parameter and for the Wilcoxon score, $\tau =[12{]}^{-\frac{1}{2}}{\left[{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}^2(t)\mathrm{dt}\right]}^{-1}$. We use the consistent estimator of $\tau$ proposed by Koul et al. [30] to calculate the HM test. The estimation of $\tau$ can be based on a window- or kernel-type estimation of the probability density function $f(.)$ of errors [7,10].

The decision rule is to reject $H_0$ if $\mathrm{HM}\ge F_{2,n-4,\alpha }$, where $F_{2,n-4,\alpha }$ is the upper $\alpha$ percentile of an $F$ distribution with 2, and $n-4$ degrees of freedom. The estimate of $\beta$ is consistent and asymptotically normal [35] for large sample sizes and no autocorrelation: $\hat{\beta }\sim N(\beta ,{\tau }^2{(X^{T}X)}^{-1}).$

We can approximate a $(1-\alpha )$100% confidence interval for ${\beta }_i$ as ${\hat{\beta }}_i\pm t_{\frac{\alpha }{2},n-4}\mathrm{se}({\hat{\beta }}_i),i=2,3,$where $\mathrm{se}({\hat{\beta }}_i)=\hat{\tau }\sqrt{(X^{T}X{)}_{\mathrm{ii}}^{-1}),}$ and $(X^{T}X{)}_{\mathrm{ii}}^{-1}$ represents the $i\text{th}$ diagonal element of $(X^{T}X{)}^{-1}$.

Similarly, under the least squares (LS) method with large sample sizes and no autocorrelation, the estimate of $\beta$ is also consistent and normally distributed: $\hat{\beta }\sim N(\beta ,{\sigma }^2{(X^{\mathrm{\prime }}X)}^{-1}).$For the parameters, ${\beta }_2$ and ${\beta }_3$, the $(1-\alpha )100\mathrm{\%}$ confidence interval is defined as ${\hat{\beta }}_i\pm t_{\frac{\alpha }{2},n-4}\mathrm{se}({\hat{\beta }}_i),i=2,3.$

Here, $\mathrm{se}({\hat{\beta }}_i)=\hat{\sigma }\sqrt{(X^{T}X{)}_{\mathrm{ii}}^{-1})}$ where $(X^{T}X{)}_{\mathrm{ii}}^{-1}$ is the $i\text{th }$ diagonal element of $(X^{T}X{)}^{-1}$.

Moreover, the test procedure for testing $H_0:{\beta }_2={\beta }_3=0$ versus $H_a:$ at least one of ${\beta }_2$ and ${\beta }_3$ is not zero under the LS method involves the $F$-test statistic: $F=\frac{\frac{\mathrm{SSE}(R)-\mathrm{SSE}(F)}{\mathrm{dfE}(R)-\mathrm{dfE}(F)}}{\mathrm{MSE}(F)/\mathrm{dfE}(F)}\sim F(2,n-4).$

We reject $H_0$ if $F\ge F_{2,n-4,\alpha }$; otherwise, we do not reject $H_0$.

The confidence interval and hypothesis testing procedures defined above are based on the normality assumption or large sample sizes, and no autocorrelation. However, when only small sample sizes are available, the approximation may not work well. We recommend using bootstrapping to construct confidence intervals and carry out hypothesis testing for the intervention effect (parameters ${\beta }_2$ and ${\beta }_3$) for both rank-based and standard LS regression approaches when there are small samples, and deviations from normality. If outliers and correlated errors exist as well, we only recommend using the rank-based method with bootstrapping. A brief description of the general algorithm for the bootstrap confidence intervals for ${\beta }_2$ and ${\beta }_3$ and the bootstrap p-value for the linear hypothesis of $H_0:{\beta }_2={\beta }_3=0$ is presented below. For more details about bootstrapping, see [38–42].

To construct $(1-\alpha )100\mathrm{\%}$ bootstrap confidence intervals for ${\beta }_2$ and ${\beta }_3$, the following steps are followed:

1. Fit the full model $Y_t={\beta }_0+{\beta }_1{T}_t+{\beta }_2{X}_t+{\beta }_3{X}_t(T_t-\frac{n}{2})+{\epsilon }_t$, $t=1,2,\dots ,n$ using LS and RB methods. Then calculate the fitted response values, ${\hat{Y}}_i$ and the estimated residuals, ${\hat{\epsilon }}_i$, from the full model.

2. Obtain a bootstrap sample of size $n$ with replacement from the estimated residuals, ${\hat{\epsilon }}_i$. Denote the bootstrap sampled residuals as ${\epsilon }_i^{\ast }$.

3. Find $Y_i^{\ast }$ using $\hat{Y}+{\epsilon }_i^{\ast }$ and fit the full model again using the LS and RB methods. Extract the estimate of ${\beta }_2$ and ${\beta }_3$ from the full model results. Label them ${\hat{\beta }}_2^{\ast i}$ and ${\hat{\beta }}_3^{\ast i}$.

4. Repeat steps 1 through 3 $R$ times and construct a confidence interval by finding the associated percentiles in the bootstrap distributions of ${{\hat{\beta }}_2}^{\ast }$ and ${{\hat{\beta }}_3}^{\ast }$.

Similarly, to calculate the bootstrap p-value for the linear hypothesis of $H_0:{\beta }_2={\beta }_3=0$, the general algorithm for the bootstrap p-value is as follows:

1. Fit the full model in (1) and the null model in (3) using LS and RB methods. Then calculate the observed value of the test statistics: $\mathrm{HM}$ and $F$. Denote it as $F_0$.

2. Take a bootstrap sample with replacement from the residuals obtained under the null model. Denote the bootstrap sampled residuals from the null model as ${\epsilon }_i^{\ast }$, and define $Y_i^{\ast }=\hat{Y}+{\epsilon }_i^{\ast }$, where $\hat{Y}$ is the fitted value from the null model. Fit the full and null model using $Y_i^{\ast }$ and calculate the Jaeckel-Hettmansperger-McKean test statistic, and the $F$-statistics. Denote it as $F^{\ast }$. It implies that fitting, and testing are performed under the assumption that $H_0$ is true.

3. Repeat steps 1 and 2 R times. The bootstrap p-value is defined as $p^{\ast }=\frac{I(F^{\ast }\ge F_0)}{R}$.

3. Simulation study

This section begins with a simulation study designed to evaluate the performance of two methods, RB and LS regression with formula-based asymptotic SEs, for estimating intervention effects in time series data. This study examines their performance under different error distributions, with and without the presence of outliers, and different levels of autocorrelation.

To calculate bias, mean squared error (MSE) and relative mean squared error (RMSE) of the estimates for the intervention effect parameters, ${\beta }_2$ and ${\beta }_3$, the response variable $y_i,i=1,2,\dots ,n$ is generated using the model: $Y_t={\beta }_0+{\beta }_1{T}_t+{\beta }_2{X}_t+{\beta }_3{X}_t\left(T_t-\frac{n}{2}\right)+{\epsilon }_t$with the parameter vector $\beta =({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3{)}^T=(1,1,0,0{)}^T$ and design matrix:

Table

$X_0$	$T_t$	$X_t$	$X_t\ast \left(T_t-{\displaystyle \frac{n}{2}}\right)$
1	1	0	0
1	2	0	0
$⋮$	$⋮$	$⋮$	0
1	${\displaystyle \frac{n}{2}}$	0	0
1	${\displaystyle \frac{n}{2}}+1$	1	1
1	${\displaystyle \frac{n}{2}}+2$	1	2
$⋮$	$⋮$	$⋮$	$⋮$
1	$n$	1	$n$

To mimic the autocorrelation commonly seen in time series data, the errors ${\epsilon }_t$ follow the first-order autoregressive AR (1) process: ${\epsilon }_t=\rho {\epsilon }_{t-1}+u_t$where $u_t$ are independent and identically distributed with mean zero and finite variance ${\sigma }^2$; and $\rho$ is the lag-1 autocorrelation of the errors. The variance-covariance matrix of the error vector of $\mathit{\epsilon }=({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_t{)}^T\text{ is }V(\mathit{\epsilon })={\sigma }^2\mathrm{\Omega }=\frac{{\sigma }^2}{1-{\rho }^2}\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{c}1\\ \rho \end{array}& \begin{array}{c}\rho \\ 1\end{array}\end{array}& \cdots & \begin{array}{c}{\rho }^{n-1}\\ {\rho }^{n-2}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{\rho }^{n-1}& {\rho }^{n-2}\end{array}& \cdots & 1\end{array}\right].$

The white noise terms $u_t$ are generated from three different symmetric distributions: (i) a normal distribution with a mean of 0, and standard deviation of 1, (ii) a $t-$ distribution with 5 degrees of freedom and (iii) a double exponential distribution (Laplace distribution) with a location parameter of 0 and a scale parameter of 1. The distributions are scaled by $\sqrt{1-{\rho }^2}$ to achieve the specified variance at different levels of $\rho =(0,0.2,0.4,0.6,0.8)$.

To introduce outliers, approximately 13% of the sample are generated using the uniform distribution over (0,30), and the responses are replaced with the generated uniform data. The RB and LS regression methods are then applied to the generated response variable $y_i,i=1,2,\dots ,n$ to estimate the parameter vector $\beta$ and its standard error. The analytical form of the variance of the estimates obtained through both the LS and RB is ${\sigma }^2(X^{T}X{)}^{-1}{X}^T\mathrm{\Omega }X(X^{T}X{)}^{-1}$. The standard error of the estimates of ${\beta }_2$ and ${\beta }_3$ obtained through both bootstrapping and asymptotic formula under LS and RB, are then compared to the true value found through the analytical formula under a given known scenario.

The bias of the estimate is defined as the difference between the expected value of the estimator and its true value. $\text{Bias = }\frac{1}{R}\sum _{i=1}^R{\hat{\beta }}_{\mathrm{ji}}-{\beta }_j,j=2,3.$The MSE of the bias is defined as the average squared difference between the estimated values and the actual value of the parameter. $\text{MSE = }\frac{1}{R}\sum _{i=1}^R({\hat{\beta }}_{\mathrm{ij}}-{\beta }_j{)}^2,j=2,3.$

The relative mean squared error (RMSE) is the ratio of the mean squared error of the RB method to the LS method for the intervention effects ${\beta }_2$ and ${\beta }_3$. The gain in efficiency is defined as (100-RMSE) *100.

For this simulation study, two different sample sizes – a relatively small sample, $n=16$ and a large sample, $n=50$ are considered and the simulation is repeated 1000 times.

Tables 1 and 2 present the gain in efficiency for ${\beta }_2$ and ${\beta }_3$ of RB method with bootstrapping with respect to the LS method under $H_0:{\beta }_2={\beta }_3=0$. Instead of using formula-based SEs through the LS and RB methods, bootstrapping SEs are calculated. Comparisons are made for different error distributions, two different levels of sample sizes, different levels of autocorrelation, and with and without outliers in the data.

Table 1. Efficiency gain when no outliers exist in the data and errors follows AR (1).

		RB method				LS method
$n$	$\rho$	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	Efficiency gain for ${\beta }_2$ in %	Efficiency gain for ${\beta }_3$ in %
Errors ∼ N (0,1)
16	0	0.00	1.04	−0.01	0.23	0.01	1.00	−0.01	0.22	−8.62	−9.73
	0.2	−0.03	1.11	0.00	0.26	−0.03	1.08	0.00	0.25	−6.06	−5.89
	0.4	−0.04	1.21	0.00	0.29	−0.04	1.18	0.00	0.29	−4.92	−4.71
	0.6	−0.05	1.21	0.01	0.32	−0.04	1.17	0.01	0.31	−5.83	−5.52
	0.8	0.00	1.00	0.00	0.29	0.00	0.97	0.00	0.28	−5.63	−1.50
50	0	0.02	0.59	0.00	0.04	0.01	0.57	0.00	0.04	−5.91	−4.51
	0.2	0.00	0.68	0.00	0.05	0.00	0.67	0.00	0.05	−2.36	−1.86
	0.4	−0.01	0.83	0.00	0.06	−0.01	0.82	0.00	0.06	−3.51	−1.51
	0.6	0.01	0.92	0.00	0.07	0.01	0.91	0.00	0.08	−2.63	0.55
	0.8	−0.04	1.09	0.00	0.09	−0.04	1.08	0.00	0.09	−3.62	−1.97
Errors∼ Laplace (0,1)
16	0	0.03	1.35	−0.01	0.29	0.06	1.43	0.00	0.31	10.55	13.64
	0.2	−0.03	1.45	−0.01	0.35	−0.07	1.53	−0.01	0.36	10.00	6.91
	0.4	0.02	1.62	0.00	0.40	0.02	1.64	0.00	0.41	2.44	4.34
	0.6	0.10	1.71	0.02	0.45	0.10	1.68	0.02	0.45	−4.18	0.66
	0.8	0.00	1.40	−0.01	0.41	−0.01	1.35	−0.01	0.41	−7.12	−1.11
50	0	0.00	0.70	0.00	0.05	0.00	0.80	0.00	0.06	23.24	26.43
	0.2	−0.01	0.84	0.00	0.06	−0.01	0.95	0.00	0.06	21.47	19.10
	0.4	0.01	1.04	0.00	0.07	0.02	1.10	0.00	0.08	11.70	13.08
	0.6	−0.01	1.27	0.00	0.09	−0.01	1.30	0.00	0.10	4.01	7.35
	0.8	0.08	1.55	−0.01	0.12	0.09	1.54	−0.01	0.12	−1.59	0.37
Errors ∼ t (5)
16	0	0.02	1.26	−0.01	0.27	−0.02	1.32	−0.01	0.28	8.02	3.63
	0.2	−0.05	1.42	0.00	0.31	−0.06	1.44	0.00	0.31	2.61	3.17
	0.4	−0.02	1.58	0.00	0.35	−0.01	1.60	−0.01	0.36	3.16	2.61
	0.6	0.01	1.56	0.02	0.39	0.00	1.54	0.02	0.40	−2.84	0.85
	0.8	0.01	1.26	−0.01	0.38	0.01	1.22	−0.01	0.38	−7.76	0.07
50	0	0.03	0.67	0.00	0.05	0.02	0.71	0.00	0.05	11.62	17.58
	0.2	−0.02	0.80	0.00	0.06	−0.02	0.88	0.00	0.06	16.04	15.51
	0.4	−0.04	0.99	0.00	0.07	−0.03	1.04	0.00	0.07	9.75	9.85
	0.6	0.03	1.21	0.00	0.09	0.02	1.23	0.00	0.09	4.37	6.83
	0.8	0.02	1.37	0.00	0.12	0.02	1.34	0.00	0.12	−3.83	1.24

Note: SEs under both methods are estimated through bootstrapping.

Table 2. Efficiency gain when outliers exist in the data and errors follow AR (1).

		RB method				LS method
$n$	$\rho$	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	Efficiency gain for ${\beta }_2$ in %	Efficiency gain for ${\beta }_3$ in %
Errors ∼ N (0,1)
16	0	0.06	1.59	0.00	0.37	−0.01	3.76	−0.01	0.86	82.10	81.17
	0.2	0.00	1.74	0.00	0.45	−0.13	3.93	−0.05	0.92	80.53	76.39
	0.4	−0.01	1.82	−0.01	0.45	−0.13	3.91	−0.03	0.91	78.30	75.41
	0.6	0.05	1.73	0.01	0.43	−0.05	3.88	0.00	0.89	80.16	76.02
	0.8	−0.01	1.40	−0.01	0.36	−0.06	3.83	0.01	0.92	86.69	84.82
50	0	−0.02	0.77	0.02	0.05	0.08	3.33	0.01	0.28	94.69	96.27
	0.2	−0.04	0.87	0.01	0.06	−0.09	3.27	0.00	0.28	92.93	95.25
	0.4	−0.03	0.92	0.02	0.07	−0.03	3.38	0.01	0.28	92.55	94.29
	0.6	0.00	1.05	0.02	0.08	0.08	3.43	0.01	0.30	90.71	92.70
	0.8	−0.04	1.13	0.01	0.10	−0.07	3.36	0.00	0.29	88.74	89.24
Errors ∼ Laplace (0,1)
16	0	0.22	2.00	0.02	0.49	0.27	3.81	0.06	0.91	72.45	71.24
	0.2	0.19	2.14	−0.03	0.51	0.35	4.02	−0.02	0.95	71.70	70.90
	0.4	0.05	2.21	−0.03	0.52	−0.05	3.89	−0.01	0.92	67.86	68.17
	0.6	0.02	2.24	0.00	0.55	−0.13	3.95	−0.01	0.97	67.77	68.20
	0.8	0.01	1.83	0.00	0.49	0.06	3.80	0.00	0.94	76.67	73.00
50	0	−0.03	0.92	0.02	0.06	−0.10	3.29	0.01	0.27	92.10	94.62
	0.2	−0.06	1.01	0.02	0.07	0.00	3.29	−0.01	0.28	90.53	93.12
	0.4	−0.10	1.21	0.02	0.09	−0.12	3.39	0.00	0.29	87.26	91.10
	0.6	−0.05	1.38	0.01	0.10	0.00	3.40	−0.01	0.29	83.57	87.18
	0.8	−0.05	1.66	0.01	0.13	−0.09	3.65	0.00	0.30	79.31	80.31
Errors ∼ t (5)
16	0	0.04	2.01	0.00	0.49	−0.06	3.92	0.00	0.92	73.81	72.47
	0.2	0.04	2.03	0.01	0.47	−0.04	3.94	0.01	0.94	73.53	74.95
	0.4	−0.06	2.11	0.00	0.51	0.00	4.05	0.00	0.94	72.96	70.15
	0.6	0.00	2.04	−0.01	0.50	−0.04	3.95	−0.01	0.96	73.38	72.81
	0.8	0.04	1.70	−0.01	0.45	−0.06	3.78	0.06	0.97	79.76	78.22
50	0	−0.03	0.87	0.02	0.06	0.06	3.33	−0.01	0.29	93.11	95.13
	0.2	−0.06	1.03	0.01	0.07	−0.02	3.35	−0.01	0.27	90.60	93.41
	0.4	−0.06	1.14	0.02	0.08	−0.06	3.43	0.00	0.28	88.99	92.05
	0.6	−0.02	1.35	0.01	0.10	0.13	3.48	−0.01	0.29	85.00	89.08
	0.8	−0.08	1.53	0.02	0.12	−0.07	3.57	0.01	0.30	81.63	83.55

Note: SEs under both methods are estimated through bootstrapping.

In Table 1, the results indicate that when there are no outliers in the data, and the normality assumption of errors holds, the LS method demonstrates modest gains in efficiency relative to the RB method, but the relative efficiency for both parameters ${\beta }_2$ and ${\beta }_3$ of the LS method with respect to RB method decreases as autocorrelation and sample sizes increase. On the other hand, when errors are non-normal, the relative efficiency for both ${\beta }_2$ and ${\beta }_3$ of the RB method are positive relative to the LS method. The gain in efficiency for ${\beta }_3$ decreases with autocorrelation but increases as sample sizes increases. For ${\beta }_2$, the gain in efficiency is negative when autocorrelation is within a range of 0.6–0.9.

Table 2 shows the most significant gains in efficiency for the RB method. It reveals that whenever outliers are present in the data, and no matter how much autocorrelation exists, the RB method achieves remarkably high gains in efficiency ranging between 71% and 96% for ${\beta }_2$ and ${\beta }_3$. These substantial improvements are consistent across all error structures, including normality with higher efficiency gain than other long-tailed error distributions. The gain slightly decreases with autocorrelation but increases as the sample size increases.

Bias and MSE of the estimates, ${\hat{\beta }}_2$ and ${\hat{\beta }}_3$, under the null hypothesis $H_0:{\beta }_2={\beta }_3=0$ are presented in Tables 3 and 4. Table 3 indicates that in the absence of outliers, and at different levels of autocorrelation, and different sample sizes, the MSE is relatively similar for the two methods, consistent with the efficiency comparisons in Table 1. It is slightly lower for the LS method under normality, while the RB method exhibits slightly lower MSE under non-normality and lower levels of autocorrelation (0 to 0.4). For higher levels of autocorrelation (0.6 to 0.8) and non-normality, the MSE is slightly lower for the LS method when estimating ${\beta }_2$. Both methods produced almost the same values for the MSE when estimating ${\beta }_3$ under these conditions. But as the sample size increases, MSE in both methods for both ${\beta }_2$ and ${\beta }_3$ decreases, but the decrease in MSE for ${\beta }_3$ is significantly higher. The bias is also quite comparable for both methods and remains close to 0 across all error structures.

Table 3. Bias of the estimate and MSE when no outliers exist, and errors follow AR (1).

		RB method				LS method
		${\hat{\beta }}_2$		${\hat{\beta }}_3$		${\hat{\beta }}_2$		${\hat{\beta }}_3$
$n$	$\rho$	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
Errors ∼ N (0,1)
16	0	0.00	1.08	−0.01	0.05	0.01	0.99	−0.01	0.05
	0.2	−0.03	1.23	0.00	0.07	−0.03	1.16	0.00	0.06
	0.4	−0.04	1.45	0.00	0.09	−0.04	1.39	0.00	0.08
	0.6	−0.05	1.45	0.01	0.10	−0.04	1.37	0.01	0.10
	0.8	0.00	0.99	0.00	0.08	0.00	0.94	0.00	0.08
50	0	0.02	0.35	0.00	0.00	0.01	0.33	0.00	0.00
	0.2	0.00	0.46	0.00	0.00	0.00	0.45	0.00	0.00
	0.4	−0.01	0.69	0.00	0.00	−0.01	0.67	0.00	0.00
	0.6	0.01	0.85	0.00	0.01	0.01	0.83	0.00	0.01
	0.8	−0.04	1.20	0.00	0.01	−0.04	1.16	0.00	0.01
Errors ∼ Laplace (0,1)
16	0	0.03	1.83	−0.01	0.08	0.06	2.04	0.00	0.10
	0.2	−0.03	2.11	−0.01	0.12	−0.07	2.35	−0.01	0.13
	0.4	0.02	2.61	0.00	0.16	0.02	2.67	0.00	0.17
	0.6	0.10	2.94	0.02	0.20	0.10	2.82	0.02	0.20
	0.8	0.00	1.95	−0.01	0.17	−0.01	1.82	−0.01	0.17
50	0	0.00	0.49	0.00	0.00	0.00	0.63	0.00	0.00
	0.2	−0.01	0.71	0.00	0.00	−0.01	0.90	0.00	0.00
	0.4	0.01	1.08	0.00	0.01	0.02	1.22	0.00	0.01
	0.6	−0.01	1.62	0.00	0.01	−0.01	1.68	0.00	0.01
	0.8	0.08	2.41	−0.01	0.02	0.09	2.37	−0.01	0.02
Errors ∼ t (5)
16	0	0.02	1.60	−0.01	0.07	−0.02	1.74	−0.01	0.08
	0.2	−0.05	2.03	0.00	0.09	−0.06	2.08	0.00	0.10
	0.4	−0.02	2.49	0.00	0.13	−0.01	2.57	−0.01	0.13
	0.6	0.01	2.44	0.02	0.15	0.00	2.37	0.02	0.16
	0.8	0.01	1.60	−0.01	0.15	0.01	1.48	−0.01	0.15
50	0	0.03	0.45	0.00	0.00	0.02	0.51	0.00	0.00
	0.2	−0.02	0.65	0.00	0.00	−0.02	0.77	0.00	0.00
	0.4	−0.04	0.97	0.00	0.00	−0.03	1.08	0.00	0.01
	0.6	0.03	1.46	0.00	0.01	0.02	1.52	0.00	0.01
	0.8	0.02	1.86	0.00	0.01	0.02	1.79	0.00	0.01

Table 4. Bias of the estimate and MSE when outliers exist, and errors follow AR (1).

		RB method				LS method
		${\hat{\beta }}_2$		${\hat{\beta }}_3$		${\hat{\beta }}_2$		${\hat{\beta }}_3$
$n$	$\rho$	Bias	MSE	Bias	MSE	Bias	MSE	Bias	MSE
Errors ∼ N (0,1)
16	0	0.06	2.53	0.00	0.14	−0.01	14.12	−0.01	0.74
	0.2	0.00	3.01	0.00	0.20	−0.13	15.47	−0.05	0.86
	0.4	−0.01	3.31	−0.01	0.20	−0.13	15.28	−0.03	0.82
	0.6	0.05	2.98	0.01	0.19	−0.05	15.01	0.00	0.79
	0.8	−0.01	1.95	−0.01	0.13	−0.06	14.67	0.01	0.85
50	0	−0.02	0.59	0.02	0.00	0.08	11.07	0.01	0.08
	0.2	−0.04	0.76	0.01	0.00	−0.09	10.68	0.00	0.08
	0.4	−0.03	0.85	0.02	0.00	−0.03	11.39	0.01	0.08
	0.6	0.00	1.10	0.02	0.01	0.08	11.79	0.01	0.09
	0.8	−0.04	1.27	0.01	0.01	−0.07	11.28	0.00	0.09
Errors ∼ Laplace (0,1)
16	0	0.22	4.04	0.02	0.24	0.27	14.58	0.06	0.82
	0.2	0.19	4.61	−0.03	0.26	0.35	16.27	−0.02	0.91
	0.4	0.05	4.87	−0.03	0.27	−0.05	15.14	−0.01	0.84
	0.6	0.02	5.01	0.00	0.30	−0.13	15.57	−0.01	0.94
	0.8	0.01	3.36	0.00	0.24	0.06	14.42	0.00	0.88
50	0	−0.03	0.85	0.02	0.00	−0.10	10.79	0.01	0.07
	0.2	−0.06	1.03	0.02	0.01	0.00	10.82	−0.01	0.08
	0.4	−0.10	1.48	0.02	0.01	−0.12	11.52	0.00	0.08
	0.6	−0.05	1.90	0.01	0.01	0.00	11.53	−0.01	0.09
	0.8	−0.05	2.76	0.01	0.02	−0.09	13.31	0.00	0.09
Errors ∼ t (5)
16	0	0.04	4.03	0.00	0.24	−0.06	15.37	0.00	0.85
	0.2	0.04	4.11	0.01	0.22	−0.04	15.51	0.01	0.88
	0.4	−0.06	4.44	0.00	0.26	0.00	16.42	0.00	0.89
	0.6	0.00	4.15	−0.01	0.25	−0.04	15.57	−0.01	0.92
	0.8	0.04	2.89	−0.01	0.21	−0.06	14.25	0.06	0.95
50	0	−0.03	0.76	0.02	0.00	0.06	11.07	−0.01	0.09
	0.2	−0.06	1.06	0.01	0.01	−0.02	11.21	−0.01	0.07
	0.4	−0.06	1.30	0.02	0.01	−0.06	11.74	0.00	0.08
	0.6	−0.02	1.82	0.01	0.01	0.13	12.12	−0.01	0.09
	0.8	−0.08	2.34	0.02	0.01	−0.07	12.72	0.01	0.09

However, when outliers are present in the data, Table 4 shows the RB method performs exceedingly better than the LS method in terms of MSE for both the ${\beta }_2$ and ${\beta }_3$ parameters. The MSE is much lower for the RB method under all three error structures and all levels of autocorrelation. As expected, the MSE declines as sample sizes increase. The bias is again near 0 for both methods under all three error structures and is consistent with the efficiency comparisons in Table 2.

Tables 5 and 6 show the type I errors of the global HM and F test as well as the individual t-test through both methods for testing the hypothesis $H_0:{\beta }_2={\beta }_3=0$ and $H_0:{\beta }_i=0,i=2,3$, respectively. The type I error measures the rate at which intervention effects are detected when in fact, there are no such effects. In the absence of outliers, both methods have type I error rates almost equal to the nominal level $\alpha =5\mathrm{\%}$ when errors are independent. Type I error is slightly better for the LS method under normality, but slightly better for the RB method under non-normality when there are no outliers in the data. Interestingly, the type I error rates decreased for both methods under all three error structures if outliers were present in the data and no autocorrelation was present. However, the decrease is more evident for the RB method, resulting in rates as low as around 3%. When autocorrelation is present, with and without outliers, the type I error inflates in both methods and becomes more pronounced as the autocorrelation increases.

Table 5. Type I error of asymptotic test HM and F for $H_0:{\beta }_2={\beta }_3=0$ when error follows AR (1).

		Without outliers		With outliers
$n$	$\rho$	RB method	LS method	RB method	LS method
Errors ∼ N (0,1)
16	0	0.06	0.06	0.03	0.04
	0.2	0.10	0.10	0.04	0.05
	0.4	0.20	0.20	0.08	0.06
	0.6	0.33	0.33	0.13	0.08
	0.8	0.45	0.48	0.15	0.06
50	0	0.07	0.06	0.05	0.06
	0.2	0.13	0.13	0.11	0.05
	0.4	0.26	0.26	0.15	0.05
	0.6	0.46	0.47	0.30	0.07
	0.8	0.66	0.68	0.52	0.06
Errors ∼ Laplace (0,1)
16	0	0.05	0.05	0.03	0.04
	0.2	0.10	0.11	0.07	0.06
	0.4	0.20	0.20	0.09	0.06
	0.6	0.37	0.34	0.17	0.09
	0.8	0.47	0.46	0.18	0.08
50	0	0.05	0.05	0.04	0.04
	0.2	0.13	0.12	0.08	0.06
	0.4	0.25	0.25	0.17	0.07
	0.6	0.46	0.44	0.29	0.07
	0.8	0.69	0.69	0.52	0.08
Errors ∼ t (5)
16	0	0.06	0.05	0.04	0.05
	0.2	0.10	0.09	0.05	0.05
	0.4	0.21	0.22	0.09	0.07
	0.6	0.34	0.33	0.15	0.08
	0.8	0.48	0.48	0.17	0.08
50	0	0.06	0.05	0.06	0.06
	0.2	0.13	0.13	0.10	0.04
	0.4	0.28	0.27	0.17	0.05
	0.6	0.44	0.45	0.29	0.07
	0.8	0.66	0.67	0.53	0.09

Table 6. Type I error for individual t-test for testing $H_0:{\beta }_2=0$ and $H_0:{\beta }_3=0$ under AR(1).

		Without outliers				With outliers
		RB method		LS method		RB method		LS method
$n$	$\rho$	${\beta }_2$	${\beta }_3$	${\beta }_2$	${\beta }_3$	${\beta }_2$	${\beta }_3$	${\beta }_2$	${\beta }_3$
Errors ∼ N (0,1)
16	0	0.05	0.05	0.04	0.05	0.04	0.03	0.04	0.03
	0.2	0.09	0.08	0.08	0.10	0.04	0.06	0.05	0.05
	0.4	0.10	0.16	0.11	0.19	0.06	0.08	0.06	0.05
	0.6	0.18	0.26	0.18	0.29	0.07	0.11	0.06	0.06
	0.8	0.21	0.39	0.24	0.41	0.09	0.15	0.06	0.08
50	0	0.06	0.05	0.06	0.05	0.03	0.04	0.02	0.06
	0.2	0.10	0.10	0.10	0.11	0.06	0.06	0.03	0.07
	0.4	0.18	0.19	0.20	0.21	0.08	0.12	0.03	0.08
	0.6	0.29	0.34	0.30	0.34	0.13	0.16	0.04	0.07
	0.8	0.45	0.53	0.47	0.55	0.23	0.30	0.04	0.09
Errors ∼ Laplace (0,1)
16	0	0.06	0.05	0.06	0.04	0.04	0.05	0.05	0.05
	0.2	0.09	0.10	0.09	0.10	0.05	0.06	0.04	0.05
	0.4	0.13	0.17	0.15	0.18	0.08	0.10	0.06	0.07
	0.6	0.19	0.28	0.17	0.28	0.10	0.13	0.06	0.08
	0.8	0.25	0.40	0.24	0.44	0.11	0.17	0.06	0.09
50	0	0.05	0.05	0.05	0.05	0.04	0.05	0.03	0.07
	0.2	0.09	0.13	0.08	0.12	0.05	0.07	0.03	0.06
	0.4	0.17	0.19	0.18	0.20	0.09	0.13	0.03	0.07
	0.6	0.29	0.33	0.28	0.33	0.16	0.21	0.03	0.07
	0.8	0.43	0.54	0.43	0.54	0.26	0.34	0.05	0.07
Errors ∼ t (5)
16	0	0.04	0.06	0.05	0.05	0.04	0.05	0.04	0.04
	0.2	0.08	0.10	0.09	0.10	0.06	0.06	0.05	0.05
	0.4	0.14	0.17	0.13	0.19	0.06	0.08	0.05	0.06
	0.6	0.17	0.28	0.18	0.30	0.08	0.13	0.05	0.08
	0.8	0.23	0.39	0.23	0.41	0.07	0.16	0.05	0.08
50	0	0.05	0.06	0.06	0.05	0.04	0.04	0.03	0.07
	0.2	0.10	0.10	0.10	0.10	0.05	0.05	0.03	0.06
	0.4	0.19	0.21	0.19	0.21	0.07	0.11	0.02	0.07
	0.6	0.30	0.30	0.30	0.33	0.14	0.22	0.03	0.06
	0.8	0.44	0.52	0.44	0.52	0.26	0.35	0.04	0.11

When there is no autocorrelation, with and without outliers, the empirical type I error rates are equal to the nominal level of 0.05 for both the HM and F statistics for the general hypothesis $H_0:{\beta }_2={\beta }_3=0$. When outliers were present, and autocorrelation increased from 0 to 0.2, the empirical type I error remained around the nominal level of 0.05 for the F test under the LS method. For higher levels of autocorrelation (0.4 to 0.8), and outliers in the data across all error structures, type I error rates increased only slightly to 0.09 for the F test. In contrast, the performance of the HM test under the RB method is not satisfactory when autocorrelation increases and outliers are present. Type I error rates for the HM test got as high as 0.53 when the autocorrelation was 0.8 and the sample size was 50. Even with a small sample size ($n=16$), the type I error rate of the HM test at autocorrelation of 0.8 is 0.52. This shows that for small samples, the normal approximation of the HM test does not perform well compared to the F test under LS. Specifically, since the HM test involves the tau parameter, which is estimated through the error distribution, if the errors are not independent then the estimate of the distribution of tau might be affected. Our study has found the performance of the HM test to be unsatisfactory in the presence of autocorrelated errors. Similar results are found for the individual t-test for testing $H_0:{\beta }_i=0,i=2,3$, as can be seen in Table 6.

When no outliers are present, but autocorrelated errors exist, then neither the HM test nor the F test perform well. Both methods have similar and extremely inflated type I error rates at all levels of autocorrelation, reaching about 0.68 when the autocorrelation is 0.8.

Additional simulations were conducted to compare SEs calculated through bootstrapping, and the asymptotic formula under both RB and LS methods, to the analytical values of SE of ${\hat{\beta }}_2$ and ${\hat{\beta }}_3$ when errors are AR (1). The SEs of the estimates of ${\beta }_2$ and ${\beta }_3$ obtained through bootstrapping, the asymptotic formula, and the analytical formula under different autocorrelation levels, with and without outliers, and different sample sizes are shown in Tables 7 and 8. When errors are not autocorrelated ($\rho =0$), there are no outliers in the data, and under all three different error distributions, the formula-based asymptotic SEs and Bootstrap SEs for the estimates of ${\beta }_2$ and ${\beta }_3$ are approximately the same for both the LS and RB regression methods, and are equal to the value of analytical formula of SEs. It is also important to note that if the errors are normally distributed, the bootstrap SEs for RB regression are almost as good as the bootstrap SEs for LS regression, even under these ideal conditions.

Table 7. Comparison of SE estimated through asymptotic and bootstrap approaches when no outliers.

		RB method				LS method
		Bootstrap SE		Asymptotic SE		Bootstrap SE		Asymptotic SE		Analytical SE
$n$	$\rho$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$
Errors ∼ N (0,1)
16	0	1.05	0.23	1.10	0.24	1.00	0.22	1.01	0.22	1.01	0.22
	0.2	1.13	0.27	1.03	0.22	1.09	0.26	0.94	0.20	1.11	0.25
	0.4	1.18	0.29	0.92	0.20	1.16	0.28	0.84	0.18	1.17	0.29
	0.6	1.17	0.33	0.79	0.17	1.14	0.32	0.71	0.15	1.16	0.32
	0.8	0.93	0.29	0.56	0.12	0.90	0.29	0.50	0.11	0.95	0.30
50	0	0.59	0.04	0.59	0.04	0.58	0.04	0.56	0.04	0.57	0.04
	0.2	0.66	0.05	0.57	0.04	0.64	0.05	0.55	0.04	0.67	0.05
	0.4	0.78	0.06	0.56	0.04	0.77	0.06	0.54	0.04	0.79	0.06
	0.6	0.94	0.07	0.52	0.04	0.93	0.07	0.50	0.03	0.95	0.07
	0.8	1.09	0.09	0.43	0.03	1.07	0.09	0.41	0.03	1.09	0.09
Errors ∼ Laplace (0,1)
16	0	1.38	0.29	1.34	0.29	1.43	0.30	1.36	0.29	2.02	0.44
	0.2	1.45	0.34	1.28	0.28	1.53	0.35	1.28	0.28	2.21	0.51
	0.4	1.70	0.39	1.19	0.26	1.70	0.40	1.18	0.25	2.35	0.58
	0.6	1.64	0.44	1.02	0.22	1.61	0.44	1.00	0.22	2.32	0.64
	0.8	1.40	0.41	0.74	0.16	1.35	0.41	0.71	0.15	1.90	0.60
50	0	0.69	0.05	0.68	0.05	0.81	0.05	0.78	0.05	1.13	0.08
	0.2	0.85	0.06	0.69	0.05	0.93	0.07	0.78	0.05	1.34	0.09
	0.4	1.05	0.08	0.70	0.05	1.11	0.08	0.74	0.05	1.59	0.12
	0.6	1.29	0.10	0.68	0.05	1.32	0.10	0.70	0.05	1.90	0.15
	0.8	1.51	0.12	0.57	0.04	1.51	0.12	0.58	0.04	2.18	0.19
Errors ∼ t (5)
16	0	1.26	0.26	1.28	0.28	1.30	0.27	1.27	0.27	1.69	0.36
	0.2	1.37	0.31	1.18	0.25	1.45	0.33	1.16	0.25	1.84	0.42
	0.4	1.53	0.37	1.12	0.24	1.56	0.38	1.07	0.23	1.95	0.48
	0.6	1.43	0.39	0.95	0.20	1.43	0.39	0.91	0.20	1.93	0.53
	0.8	1.30	0.38	0.69	0.15	1.27	0.38	0.66	0.14	1.58	0.50
50	0	0.68	0.05	0.68	0.05	0.73	0.05	0.72	0.05	0.94	0.07
	0.2	0.78	0.06	0.66	0.05	0.85	0.06	0.70	0.05	1.11	0.08
	0.4	0.96	0.07	0.65	0.05	1.01	0.07	0.68	0.05	1.32	0.10
	0.6	1.18	0.09	0.63	0.04	1.22	0.09	0.63	0.04	1.58	0.12
	0.8	1.43	0.12	0.53	0.04	1.42	0.12	0.53	0.04	1.82	0.16

Table 8. Comparison of SE estimated through asymptotic and bootstrap approaches when outliers.

		RB method				LS method
		Bootstrap SE		Asymptotic SE		Bootstrap SE		Asymptotic SE		Analytical SE
$n$	$\rho$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$	${\hat{\beta }}_2$	${\hat{\beta }}_3$
Errors ∼ N (0,1)
16	0	1.79	0.38	1.87	0.40	3.84	0.88	3.81	0.82	1.01	0.22
	0.2	1.91	0.43	1.83	0.39	3.81	0.95	3.73	0.80	1.11	0.25
	0.4	1.63	0.40	1.70	0.37	3.64	0.86	3.69	0.80	1.17	0.29
	0.6	1.98	0.47	1.60	0.35	3.90	0.93	3.62	0.78	1.16	0.32
	0.8	1.31	0.41	1.44	0.31	3.78	0.94	3.65	0.79	0.95	0.30
50	0	0.77	0.05	0.82	0.06	3.26	0.29	3.74	0.26	0.57	0.04
	0.2	0.86	0.06	0.80	0.06	3.44	0.27	3.76	0.26	0.67	0.05
	0.4	0.90	0.07	0.77	0.05	3.26	0.28	3.74	0.26	0.79	0.06
	0.6	1.10	0.08	0.74	0.05	3.42	0.28	3.77	0.26	0.95	0.07
	0.8	1.16	0.10	0.62	0.04	3.33	0.30	3.72	0.26	1.09	0.09
Errors ∼ Laplace (0,1)
16	0	2.15	0.52	2.25	0.48	3.81	0.91	3.86	0.83	2.02	0.44
	0.2	2.11	0.52	2.13	0.46	3.92	0.95	3.87	0.84	2.21	0.51
	0.4	2.17	0.56	2.03	0.44	3.91	0.98	3.83	0.83	2.35	0.58
	0.6	2.21	0.56	1.84	0.40	4.11	0.93	3.72	0.80	2.32	0.64
	0.8	1.74	0.50	1.63	0.35	3.79	0.95	3.72	0.80	1.90	0.60
50	0	0.93	0.06	0.97	0.07	3.41	0.28	3.80	0.26	1.13	0.08
	0.2	1.04	0.07	0.97	0.07	3.39	0.28	3.81	0.26	1.34	0.09
	0.4	1.15	0.09	0.97	0.07	3.44	0.29	3.77	0.26	1.59	0.12
	0.6	1.45	0.11	0.93	0.06	3.46	0.30	3.79	0.26	1.90	0.15
	0.8	1.61	0.13	0.81	0.06	3.40	0.29	3.75	0.26	2.18	0.19
Errors ∼ t (5)
16	0	1.84	0.44	2.13	0.46	3.65	0.91	3.80	0.82	1.69	0.36
	0.2	2.04	0.49	2.03	0.44	3.97	0.93	3.82	0.82	1.84	0.42
	0.4	2.04	0.49	1.96	0.42	3.96	0.93	3.86	0.83	1.95	0.48
	0.6	2.01	0.50	1.76	0.38	3.90	0.94	3.70	0.80	1.93	0.53
	0.8	1.78	0.46	1.52	0.33	3.82	0.92	3.58	0.77	1.58	0.50
50	0	0.88	0.06	0.94	0.06	3.32	0.28	3.79	0.26	0.94	0.07
	0.2	0.95	0.07	0.92	0.06	3.45	0.29	3.81	0.26	1.11	0.08
	0.4	1.13	0.09	0.91	0.06	3.39	0.29	3.77	0.26	1.32	0.10
	0.6	1.31	0.10	0.86	0.06	3.62	0.29	3.76	0.26	1.58	0.12
	0.8	1.51	0.13	0.74	0.05	3.56	0.30	3.73	0.26	1.82	0.16

However, as the autocorrelation increases, the formula-based asymptotic SEs, which are not adjusted for autocorrelation under either approach, decreases. In contrast, the bootstrap SEs increase as autocorrelation increases. This indicates that the formula-based asymptotic SEs for both LS and RB are downwardly biased, as shown in Table 7. Bootstrap estimates of standard errors should be considered if there are any signs of autocorrelation in the data.

When there are outliers in the data, the formula-based asymptotic SEs under the LS method are almost twice the asymptotic SEs under the RB method. As before, asymptotic SEs under both methods remain the same as autocorrelation increases, which again produces downwardly biased SEs. However, the bootstrap SEs for ${\beta }_2$ and ${\beta }_3$ in the RB method again increase as autocorrelation increases, and these SEs are relatively smaller than the bootstrap SEs obtained through the LS method (Table 8). The SEs through bootstrapping under the RB method are almost equal to the true value of the SEs, which indicates bootstrap SEs captures autocorrelation in the data and SE estimates through RB are more reliable. But bootstrap SEs through the LS method are not close to the true value of the SEs, indicating SE estimates obtained through the LS method are less reliable.

The asymptotic SEs from LS regression are the best under normality with no outliers or autocorrelation. However, the bootstrapping methods are almost as good under these conditions and the difference might be considered negligible. Computing formula-based asymptotic SEs using either LS or RB regression is not appropriate in the presence of autocorrelation, and some type of bootstrapping is recommended. If no outliers are present, then bootstrap SEs from LS and RB regression will be very similar. If outliers are present, then much better results are obtained using RB with bootstrapping.

For completeness, power comparisons were also made between the LS and RB methods using asymptotic SEs. The power of a test is the probability of rejecting the null hypothesis when it is false. The power for the HM and F-test, as described in Section 2, is estimated through Monte Carlo simulations when there is no autocorrelation in the error structure. The simulations for power comparisons consider two error distributions: a standard normal distribution without outliers, and a non-normal distribution with outliers (t-distribution with 5 degrees of freedom). Approximately 13% of the sample observations are generated using the uniform distribution over (0,30) and added to the response variable as outliers. The response variable $Y_t$ is generated using model (1) with different effect sizes for ${\beta }_2$ and ${\beta }_3$. The null regression coefficients configuration $({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3)=$(1, 1, 0, 0), and effect sizes (0, 0.2, 0.4, … ,1.4) for ${\beta }_2$ and ${\beta }_3$  are considered.

Figure 1 shows the power comparisons in two situations: when errors follow a normal distribution without outliers (left panel) and when errors deviate from normality with the presence of outliers (right panel). In the ideal condition, the LS method outperforms the RB method as expected, but the difference is very small. However, when the error distribution deviates from normality and outliers are present, the RB method significantly outperforms the LS method. The LS method has very little power even at very large effect sizes when outliers are present.

Figure 1. Power with normally distributed errors and no outliers (left panel) and non-normality for errors and outliers existing (right panel).

If there is no autocorrelation present, the simulation study shows that the LS method outperforms the RB method only under normality for the error structure, and no outliers in the data. If there exist correlated and non-normal errors and/or outliers exist in the data, then the RB method outperforms the LS method. The difference is drastic if outliers are present in the data. If autocorrelated errors are present, and sample size is small, a bootstrap hypothesis technique would be appropriate because type I errors associated with using HM or F or the individual t-test for hypothesis testing of intervention effects can become drastically inflated for both methods in the presence of autocorrelation.

4. Application to illegal organ transplants data

In this section, we use a dataset collected and analyzed by Islam et al. [14,43] that consists of incidences of illegal organ transplants from 11 countries over the period of 2002–2015. This dataset includes data on organ transplantation-related crimes reported on the Internet from 11 countries, which were randomly selected from among the original participants of the Istanbul Declaration. The selected countries are Brazil, Columbia, Egypt, India, Iran, Mexico, Nigeria, Pakistan, Philippines, Thailand and Turkey. Additionally, the dataset contains each county’s population size, which was obtained from the World Bank database [44]. To adjust for population size, the number of illegal organ transplants are normalized by the population size of each country.

The 2008 Istanbul Declaration is the consensus of the participants to develop a legal and professional framework that governs organ donation and transplantation activities [45,46]. It also establishes a transparent regulatory oversight system that ensures donor and recipient safety, enforces standards and prohibits unethical practices [45,46]. The signee countries are encouraged to take actions to protect the poorest and most vulnerable groups from transplant tourism and the sale of tissues and organs. Additionally, the declaration seeks to reduce the prevalence of international trafficking in human tissues and organs [45–48] . The signee countries are obligated to introduce programs that control and reduce unethical practices related to organ trafficking within the Istanbul Declaration framework [49].

Despite these efforts, numerous reports indicate that organ trafficking, and transplant tourism continue to increase [43]. Therefore, it is of interest to determine whether the Declaration has had a significant impact. In this study, we estimate the intervention effects of the 2008 Istanbul Declaration on illegal organ transplants using both the RB regression and the LS regression approaches, as described in Section 2. The estimated intervention effects are then used to test $H_0:{\beta }_2={\beta }_3=0$ (i.e. the Declaration had no impact on organ trafficking) versus $H_a:$ at least one of ${\beta }_2$ and ${\beta }_3$ is not zero (i.e. the Declaration had a significant impact on organ trafficking).

Before estimating the effects of the Istanbul Declaration and conducting hypothesis testing, it is necessary to check the linearity condition of the RB and LS regression techniques described in Section 2. To do so, we present the time series plots for the number of organ transplants reported over the period 2002–2015 for each country in Figure 2. The data reveals a clear upward trend in organ trafficking for the period 2002 to 2015. However, the increase in organ trafficking is exponential over time in some countries, resulting in a non-linear trend. However, in other countries, the upward trend is linear. To achieve linearity, we took the logarithm of the rate of the incidences per million population. Figure 3 displays the nearly linear trend after the transformation of the data. In addition, the Durbin-Watson test at the 5% level of significance was used to check for serial correlation in each country’s data. The test results indicated no evidence of serial correlation in the data for any country. Based on results from the simulation (Table 5), this would imply that inflated type I errors should not be a major concern for this data set. We estimated ${\beta }_2$ and ${\beta }_3$ for the organ transplant data using both LS and RB methods and computed 95% confidence intervals using procedures based on normality, as well as bootstrap procedures as described in Section 2.

Figure 2. Time series plot for the rate of organ trafficking per million population.

Figure 3. Time series plot for the log-rate of organ trafficking per million population.

Tables 9 and 10 present LS estimates, standard errors and 95% confidence intervals for the impact of the Istanbul Declaration on illegal organ transplants. Specifically, Table 9 shows the results for LS estimates computed using the standard LS techniques and Table 10 shows the LS results using the bootstrap procedures described in Section 2. The small negative values of ${\beta }_3$ in Tables 9 and 10 show that while the number of illegal transplants has continued to increase over time, the Declaration has at least been successful in slowing down the rate of increase in illegal transplants in most countries. Except for Columbia, Nigeria and Turkey, almost all other countries experienced a statistically significant reduction in the rate of increase in organ transplants after the Declaration. Most of the values for ${\beta }_2$ are positive, but a direct interpretation is difficult. It is best just to think of ${\beta }_2$ as part of the adjustment to the intercept for the post-intervention model. For Turkey however, the coefficient for ${\beta }_3$ was almost 0, but ${\beta }_2$ was positive. This indicates that there was a modest flat increase in organ trafficking post-intervention in Turkey, which was in addition to the continuing steady trend.

Table 9. Results of standard ls regression analysis.

	Country	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	95% CI for ${\beta }_2$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	95% CI for ${\beta }_3$
Outliers	Brazil	0.35**	0.138	(0.04, 0.66)	−0.04	0.034	(−0.11, 0.04)
	Egypt	0.50**	0.162	(0.13, 0.86)	−0.05	0.040	(−0.14, 0.04)
	Iran	0.15	0.119	(−0.12,0.41)	−0.14**	0.029	(−0.20, −0.07)
	Nigeria	0.41**	0.150	(0.08,0.75)	0.03	0.037	(−0.05, 0.11)
	Turkey	0.32**	0.134	(0.02,0.62)	0.01	0.033	(−0.06, 0.08)
No outliers	Colombia	0.67**	0.215	(0.19,1.45)	0.04	0.053	(−0.08, 0.15)
	India	0.21**	0.088	(0.02,0.41)	−0.05	0.022	(−0.10, 0.01)
	Mexico	0.17	0.100	(−0.05,0.39)	−0.12**	0.025	(−0.17, −0.06)
	Pakistan	0.49**	0.136	(0.18,0.79)	−0.14**	0.033	(−0.21, −0.06)
	Philippines	0.56**	0.151	(0.22,0.89)	−0.09**	0.037	(−0.17, −0.01)
	Thailand	0.29	0.141	(−0.02,0.61)	−0.15**	0.035	(−0.23, −0.07)

Table 10. Bootstrap results of standard ls regression analysis.

	Country	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	95% CI for ${\beta }_2$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	95% CI for ${\beta }_3$
Outliers	Brazil	0.36**	0.139	(0.10, 0.63)	−0.04	0.035	(−0.10, 0.03)
	Egypt	0.49**	0.159	(0.16, 0.80)	−0.05	0.040	(−0.12, 0.03)
	Iran	0.15	0.123	(−0.11,0.38)	−0.14**	0.029	(−0.19, −0.08)
	Nigeria	0.42**	0.157	(0.12,0.71)	0.03	0.040	(−0.05, 0.10)
	Turkey	0.32**	0.130	(0.08,0.60)	0.01	0.034	(−0.06, 0.07)
No outliers	Colombia	0.67**	0.217	(0.25,1.11)	0.04	0.053	(−0.06, 0.14)
	India	0.21**	0.090	(0.04,0.39)	−0.05**	0.023	(−0.09, −0.01)
	Mexico	0.16	0.102	(−0.03,0.37)	−0.12**	0.025	(−0.17, −0.07)
	Pakistan	0.49**	0.136	(0.21,0.74)	−0.14**	0.034	(−0.20, −0.07)
	Philippines	0.56**	0.154	(0.25,0.85)	−0.09**	0.038	(−0.17, −0.02)
	Thailand	0.29**	0.138	(0.03,0.57)	−0.15**	0.034	(−0.21, −0.08)

The estimates obtained through asymptotic and Bootstrap methods under LS are consistent with each other, but the bootstrap method produced slightly narrower confidence intervals. Notably, the impact of the intervention with regards to the change in trend $({\beta }_3)$ in India is significant according to the bootstrap method, which differs from the LS method.

Similarly, Tables 11 and 12 present the findings from RB regression. Table 11 presents the results of RB regression using the asymptotic normality procedures and Table 12 presents the RB regression results using the bootstrap procedures. The estimate of differential changes in the slope of the linear model before and after the declaration, measured by ${\hat{\beta }}_3$, suggests that all countries except Turkey, Columbia and Nigeria experienced a reduction in the rate of increase in illegal organ transplants following the 2008 Istanbul Declaration. Iran, Mexico, Pakistan, Philippines and Thailand were the only countries that had a statistically significant decline in the rate of increase in illegal organ transplants. Consistent with LS regression, all countries had positive estimates of ${\beta }_2$. Table 9 shows that RB regression using bootstrap methods produced shorter confidence intervals than RB regression based on asymptotic normality. This resulted in identifying one more country, India, as having a statistically significant reduction in the rate of increase in illegal organ transplants. The biggest difference between Tables 11 and 12 is with regards to ${\beta }_2$. All estimated coefficients are still positive and similar in magnitude, but the smaller standard errors from bootstrapping resulted in nine statistically significant results instead of only five.

Table 11. Results of rank-based regression analysis.

	Country	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	95% CI for ${\beta }_2$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	95% CI for ${\beta }_3$
Outliers	Brazil	0.33**	0.123	(0.06, 0.61)	−0.04	0.030	(−0.11, 0.02)
	Egypt	0.48**	0.183	(0.07, 0.89)	−0.04	0.045	(−0.14, 0.06)
	Iran	0.16	0.123	(−0.11, 0.44)	−0.15**	0.030	(−0.22, −0.08)
	Nigeria	0.43**	0.148	(0.10, 0.76)	0.03	0.036	(−0.05,0.11)
	Turkey	0.34	0.158	(−0.02,0.69)	0.00	0.038	(−0.08, 0.09)
No outliers	Colombia	0.62	0.281	(−0.01, 1.24)	0.05	0.069	(−0.11, 0.20)
	India	0.20	0.123	(−0.08, 0.45)	−0.05	0.030	(−0.12, 0.02)
	Mexico	0.16	0.130	(−0.13, 0.44)	−0.12**	0.031	(−0.19, −0.05)
	Pakistan	0.50**	0.147	(0.17,0.83)	−0.12**	0.036	(−0.20, −0.04)
	Philippines	0.56**	0.127	(0.28, 0.85)	−0.09**	0.031	(−0.16, −0.02)
	Thailand	0.25	0.157	(−0.10, 0.60)	−0.15**	0.038	(−0.24, −0.06)

Table 12. Bootstrap results of rank-based regression analysis.

	Country	${\hat{\beta }}_2$	$\widehat{\mathrm{se}}({\hat{\beta }}_2)$	95% CI for ${\beta }_2$	${\hat{\beta }}_3$	$\widehat{\mathrm{se}}({\hat{\beta }}_3)$	95% CI for ${\beta }_3$
Outliers	Brazil	0.33**	0.122	(0.09, 0.58)	−0.05	0.029	(−0.10, 0.01)
	Egypt	0.48**	0.140	(0.21, 0.77)	−0.05	0.035	(−0.11, 0.02)
	Iran	0.16	0.100	(−0.03, 0.38)	−0.15**	0.023	(−0.20, −0.10)
	Nigeria	0.43**	0.134	(0.16, 0.72)	0.03	0.031	(−0.03,0.10)
	Turkey	0.37**	0.120	(0.10,0.56)	0.00	0.029	(−0.05, 0.06)
No outliers	Colombia	0.62**	0.195	(0.21, 1.00)	0.05	0.046	(−0.05, 0.14)
	India	0.19**	0.080	(0.04, 0.35)	−0.05**	0.020	(−0.09, −0.01)
	Mexico	0.15	0.090	(−0.02, 0.33)	−0.12**	0.022	(−0.17, −0.08)
	Pakistan	0.50**	0.118	(0.27,0.74)	−0.12**	0.028	(−0.18, −0.07)
	Philippines	0.56**	0.131	(0.28, 0.84)	−0.09**	0.031	(−0.16, −0.03)
	Thailand	0.25**	0.120	(0.02, 0.49)	−0.15**	0.030	(−0.21, −0.09)

To assist with the interpretation of the estimated parameters in context with organ trafficking, examples using the data for both Thailand and Turkey are discussed. Figure 4 shows the fitted model for Thailand both pre- and post-intervention. The trend is clearly still increasing, but at a lower rate post-intervention. This coincides with the estimate for ${\beta }_3=-0.15$. The estimate for ${\beta }_2=0.25$, but this should be interpreted as part of the overall adjustment to the intercept for the post-intervention model given below: ${\beta }_2-7{\beta }_3=0.25-7(-.15)=1.10$While ${\beta }_2$ will normally not be interpreted directly, there is an exception if ${\beta }_3=0$. In this case ${\beta }_2$ represents a flat and consistent increase in addition to the original and continuing trend. This is the case in Figure 4 for Turkey where the estimate for ${\beta }_2=0.37$ and the estimate for ${\beta }_3=0.$

Figure 4. Thailand and Turkey fitted model.

Table 13 shows the p-values obtained through the $\mathrm{HM}$, $F$ and Bootstrap tests for the general linear hypothesis $H_0:{\beta }_2={\beta }_3=0$. The results of the general linear hypothesis show that using RB regression, countries such as Columbia, India and Turkey experienced a non-significant change in illegal organ transplants, while Brazil and Turkey were non-significant countries through the standard LS regression method.

Table 13. Comparisons of HM and F-test for testing $H_0:{\beta }_2={\beta }_3=0$.

		Rank-based test			Standard LS test
	Country	Value of HM Test	p-value	Bootstrap p-value	Value of F-statistic	p-value	Bootstrap p-value
Outliers	Brazil	4.44**	0.04	0.042	3.50	0.07	0.069
	Egypt	4.36**	0.04	0.045	5.03**	0.03	0.028
	Iran	9.94**	0.00	0.004	11.27**	0.00	0.007
	Nigeria	4.50**	0.04	0.037	4.40**	0.04	0.038
	Turkey	2.89	0.10	0.133	2.98	0.10	0.105
No outliers	Colombia	4.10	0.05	0.065	5.36**	0.03	0.018
	India	3.00	0.10	0.109	4.71**	0.04	0.039
	Mexico	7.70**	0.00	0.012	12.35**	0.00	0.000
	Pakistan	10.03**	0.00	0.004	13.19**	0.00	0.003
	Philippines	9.84**	0.00	0.006	8.90**	0.01	0.005
	Thailand	8.25**	0.01	0.005	10.26**	0.00	0.001

LS and RB regression procedures, as outlined in Section 2, yielded mixed results with the organ trafficking data (Tables 10 and 12). In some cases, LS regression produced smaller confidence intervals but in other cases RB regression produced smaller confidence intervals. Based on the simulation results of Section 3, this would suggest that the error distributions were approximately normal, or close to it, and that there were no extreme outliers. Interestingly, the best results were obtained using RB regression with bootstrapping procedures. This produced the smallest confidence intervals for ${\beta }_2$ and ${\beta }_3$, and the most statistically significant results.

5. Discussions

Over the past two decades, interrupted time-series (ITS) designs have emerged as a significant approach in healthcare research and social sciences [5,7,8], though there is considerable variation in the statistical analysis of ITS designs. For large sample ITS studies, common methods include ordinary least squares with standard error adjustments, generalized least squares and ARIMA models. These models account for autocorrelation and provide better inferences with accurate calculation of standard errors, and reliable hypothesis tests. Normality assumptions are often robust in large sample studies, and valid results can be obtained through the normal approximation. But for small sample ITS studies, normal approximations are not appropriate when errors are non-normal or contain outliers. Confidence intervals and hypothesis tests become inefficient and biased using existing standard regression-based methods, even with autocorrelation corrections indicating that alternative methods are required. To address issues with small samples, particularly outliers and non-normality in the error distribution, this article introduces rank-based (RB) regression with bootstrapping to estimate the intervention effects and make inferences in ITS analysis.

Simulations compared LS regression with RB regression using relative efficiency and MSE. The performance of both methods depended on the error distribution and outliers. If no outliers are present, LS regression is slightly better if errors are normally distributed, while RB regression performs slightly better when the errors are non-normal. However, in the presence of outliers, RB regression significantly outperforms LS regression, regardless of the error distribution. Both methods produce low bias for trend and intervention effects under all error structures, and with and without outliers. Without outliers, LS is more efficient than RB, but with outliers, RB achieves substantially higher efficiency gains. These gains range from 71% to 96% across all error distributions, including normality as well as in the presence of autocorrelation.

Type I error rates were also examined for formula-based LS and RB regression. When no outliers are present, both methods have similar type I error rates that are near the nominal level of 0.05 across all error structures. When outliers are present in the data, the type I error rates decrease from the nominal rate for both methods under all three error structures. However, the tests derived for testing for the intervention effect under the RB method are more conservative than the tests obtained under the LS method when outliers exist, regardless of the underlying error distribution. This suggests the LS tests are more likely to incorrectly reject the null hypothesis when outliers are present, while the RB method provides more rigorous control of false positives.

Note that the standard formula-based methods of LS and RB regression do not account for any autocorrelation in the error structure. If autocorrelation is detected, standard errors can be severely underestimated and lead to inflated type I errors. To that end, additional simulations were conducted and included bootstrap estimation of the standard errors for both LS and RB regression, and under autocorrelation levels of 0, 0.2, 0.4, 0.6 and 0.8. As expected, the standard errors computed with the formula-based LS and RB regression models stayed relatively constant even as autocorrelation increased. However, the bootstrap estimation of the standard errors for both LS and RB correctly accounted for any autocorrelation present and increased as the autocorrelation increased. With no outliers present, both LS bootstrapping and RB bootstrapping produced similar results, but LS bootstrapping was slightly better under normality or autocorrelation levels of 0.6 and 0.8 while RB bootstrapping was better under non-normality and all levels of autocorrelation. When outliers were present, RB bootstrapping had much lower standard errors than LS bootstrapping, regardless of the error structure and all levels of autocorrelation.

The simulation results suggest the RB method with bootstrapping should be recommended for ITS data analysis on small samples, unless the error structure appears normally distributed with no outliers. The RB method provides more efficient estimation of intervention effects and inferences compared to the least squares (LS) method under non-normality, and especially if any outliers exist. In addition, the presence of autocorrelation is often underestimated in ITS studies, in particular for small sample sizes [8]. However, this is not a concern when using RB with bootstrapping because it will automatically adjust for any autocorrelation when calculating standard errors. The RB method with bootstrapping provides more robust and accurate estimates of intervention effects and statistical tests when assumptions for standard LS regression are violated, which is common with small ITS sample sizes.

The findings of the simulation are supported by analysis of a real dataset evaluating the effectiveness of the Istanbul Declaration on illegal organ transplants. The organ trafficking data was analyzed using the formula-based methods of LS and RB regression, as well as both LS and RB methods with bootstrapping. Results indicate that in most countries, while organ trafficking continues to rise following the Declaration, the rate of increase has decreased for most countries.

Residual analysis, presented in supplemental figures, showed potential outliers in the data for Brazil, Iran, Nigeria, Pakistan, Philippines and Thailand. The remaining countries had no outliers, and near or approximately normal residuals. However, regardless of the error structure in the organ trafficking data, RB with bootstrapping gave the best results for all eleven countries. It produced the smallest confidence intervals and the most significant results for the general linear hypothesis.

The simulation showed that under perfect normality and no outliers, standard formula-based LS regression is slightly better than RB with bootstrapping, but not by very much. Even if a researcher used RB with bootstrapping under ideal conditions, the difference between their results and the results from LS regression would be negligible (Table 11). How often do errors follow a normally distribution exactly? The organ trafficking example showed several countries with residual plots that appeared ‘approximately normal’, which would normally suggest using standard LS regression. However, even in these cases, RB with bootstrapping performed better than standard LS regression on real data. It always had the lowest standard errors, narrowest confidence intervals, and while it ultimately didn’t make much difference in terms of statistical significance, it did have two more statistically significant results. The real data demonstrates the robustness of the RB method and its superiority to LS regression in analysing interrupted time series data when sample sizes are small, and any type of deviation from normality exists in the error structure, or outliers are present.

Though the formula-based procedure for RB regression did not perform well on the real data in this study, simulations showed that it can outperform standard LS regression in the presence of outliers, or non-normality in the error distribution. However, for the real data in this study, the outliers were probably not extreme enough and the departures from normality not severe enough to show its effectiveness. Even so, if extreme outliers are present and departures from normality are severe, this study suggests through simulation that RB with bootstrapping will perform just as well or better than formula-based RB regression when sample sizes are small.

There are limitations to the applicability of our findings for RB with bootstrapping. The simulations assumed equal error variances, no time-varying variables and linearity. A transformation to achieve linearity was conducted on real data in this study, but including time-varying covariates and investigating procedures for unequal error variances could be valuable.

Another issue in ITS analysis is missing data. Dropping observations with missing observations can significantly impact estimation and inference with small samples. A potential solution is to impute missing pre-intervention data using prior observations, and impute missing post-intervention data using subsequent observations, before fitting the intervention model. This study also does not account for external factors influencing the outcome variable. However, ITS analysis uses full longitudinal data to model pre-intervention trends, inherently accounting for external effects [50]. If time-varying impacts like seasonality or population changes are known, they can be controlled by including them as covariates.

The analysis of real ITS data using the RB method with bootstrapping relies on computational implementation of bootstrapping procedure. In this study, percentile intervals were constructed to obtain confidence intervals for the trend and intervention effect parameters, assuming symmetry in the distribution of the estimates and this assumption is supported by the data. However, the distribution of the intervention effect may be asymmetric. The use of percentile intervals under an assumption of symmetry is a limitation, as it may result in inaccurate confidence intervals if the true distribution is skewed. Further research on the shape of the sampling distribution for the intervention effect obtained through RB with bootstrapping could inform whether percentile intervals are appropriate, or if bias-corrected and accelerated intervals should be preferred to account for potential asymmetry.

Further investigating the performance of RB with bootstrapping under changing error variances, missing data and additional covariates would be valuable insights into its robustness for real-world ITS analyses. Such extensions will be important in future work to establish guidelines for properly applying this method. While this simulation study has limitations in replicating more diverse scenarios of ITS studies, the current findings demonstrate the potential of RB regression with bootstrapping to improve the validity of ITS studies with small samples, outliers, non-normal errors and even hard to detect autocorrelation. Additional research on its properties under various assumptions will help translate these benefits into practice.

5.1. Practical implications of the study

Some practical implications of this study for practitioners analyzing interrupted time series data include:

1. Exploratory data analysis is crucial. Before fitting a model, practitioners should check for outliers, non-normal errors, autocorrelation and other assumption violations using plots and statistical tests. This will inform the choice of analysis methods.

2. If outliers or non-normal errors are present, which is common in real-world data, rank-based regression with bootstrapping is recommended over standard linear regression. It provides more unbiased, efficient and powerful analysis as shown in the simulations. In fact, this study showed using real data that if even modest departures from normality exist, RB with bootstrapping will outperform LS regression.

3. Bootstrapping should be used for inference when sample sizes are small, which is typical in ITS studies. Normal approximations can be unreliable, as shown in the simulations, and bootstrap confidence intervals and hypothesis tests outperformed asymptotic methods. In addition, RB with bootstrapping will automatically adjust for any autocorrelation that is present in the error structure. Current formula-based methods such as Newey-West and Prais-Winsten only account for an AR (1) autocorrelation structure, and there is no guarantee that the true autocorrelation structure is AR (1). Even if it is, the value for $\rho$ must be estimated as well, and simulation studies have shown that $\rho$ is often underestimated or even undetected with small samples [8]. Bootstrapping eliminates these difficulties.

4. The rank-based method with bootstrapping is easy to implement in available software like R. This makes adopting these robust techniques accessible even without statistical expertise.

5. This study provides guidance to practitioners on checking assumptions and choosing appropriate statistical models for reliable ITS analysis in fields like public health, medicine, policy evaluation and environmental studies.

6. Conclusion

This study demonstrates that rank-based regression with bootstrapping is advantageous for interrupted time series analysis when standard regression assumptions are violated, and sample sizes are small. The simulations and real data application show this method provides more precise effect estimates, rigorous statistical inference and higher power-especially with non-normal errors and outliers. As small samples and assumption deviations are common in ITS studies which evaluate policy impacts over time, rank-based modelling with bootstrapping has the potential to promote more valid analyses in fields like public health and policy research. The choice of regression methods should be determined with the aid of exploratory data analysis. By checking assumptions, and applying appropriate methods, practitioners can improve the reliability of quantitative evaluations assessing the effectiveness of interventions through interrupted time series data.

Disclosure statement

No potential conflict of interest was reported by the author(s).