Beta Regression with Autoregressive Errors for Interrupted Time Series Analysis of Proportion and Rate Outcomes: A Simulation Study

Technical Summary: Beta Regression with Autoregressive Errors for Interrupted Time Series Analysis

Problem Statement
Interrupted time series analysis (ITSA) is a standard design for evaluating interventions in health and policy research. However, when outcomes are proportions or rates bounded within the unit interval (0, 1), standard analytical approaches face two distinct forms of model misspecification:

1. Bounded Support: Ordinary Least Squares (OLS) assumes an unbounded, constant-variance conditional mean. Applying OLS to bounded data can yield fitted values outside the admissible range (0, 1) and ignores the intrinsic mean-variance relationship where variance shrinks near the boundaries. While Beta regression addresses this by modeling the conditional mean on a logit scale and a precision parameter, it typically assumes independent observations.
2. Autocorrelated Errors: Time series data frequently exhibit serial dependence. Ignoring autoregressive (AR) error structures leads to biased standard errors, inflated Type I error rates, and unreliable confidence intervals. The standard remedy for continuous outcomes—using OLS with Heteroskedasticity- and Autocorrelation-Consistent (HAC) standard errors (e.g., Newey–West)—is often applied to bounded outcomes via quasi-binomial Generalized Linear Models (GLMs). However, this two-stage approach (fitting a mean model and correcting variance post-estimation) fails to jointly model the bounded conditional density and the autoregressive error structure.

No existing estimator simultaneously addresses the bounded nature of the outcome and the autoregressive error structure within a single likelihood framework.

Methodology
The study introduces betark, a Stata implementation of a joint conditional maximum likelihood estimator for Beta regression with AR($k$) errors. The methodology relies on the recursive substitution approach (Rocha and Cribari-Neto; extended by Ferreira et al.), which yields a closed-form conditional Beta likelihood for an autoregressive process of arbitrary order $k$.

• Model Specification: The conditional mean $\mu_t$ is linked to covariates via a logit function, $g(\mu_t) = x'_t\beta$. The disturbance term $\xi_t = g(\mu_t) - x'_t\beta$ follows an AR($k$) process: $\xi_t = \sum_{i=1}^k \rho_i [g(y_{t-i}) - x'_{t-i}\beta] + r_t$. Crucially, the autoregressive feedback is constructed from realized, link-transformed lagged outcomes rather than latent errors, respecting the bounded support of $y_t$.
• Estimation: The model is estimated by maximizing the joint conditional log-likelihood over the mean coefficients ($\beta$), precision parameters ($\delta$), and autoregressive coefficients ($\rho$) simultaneously. This joint estimation ensures that the resulting covariance matrix inherently reflects the autocorrelation structure, eliminating the need for separate sandwich or two-stage corrections.
• Comparator: The study evaluates betark against the current standard alternative: a quasi-binomial GLM with a logit link, paired with Newey–West HAC standard errors (implemented in Stata as glm with vce(hac nwest #)). This comparator estimates the mean model first and applies a nonparametric variance correction post-estimation.

Simulation Design
A Monte Carlo simulation study was conducted to evaluate finite-sample properties across:

• Design: Single-group ITSA with a single intervention point.
• Parameters: Three series lengths ($T \in \{100, 200, 400\}$), three AR orders ($k \in \{1, 2, 3\}$), and three autocorrelation scenarios (mild positive, oscillatory, and high persistent).
• Outcomes: Bias, standard error ratio, 95% confidence interval coverage, Type I error, and statistical power.
• Sensitivity: Analyses examined AR order misspecification (off by one lag) and variations in starting mean and pre-intervention trends.

Key Results

1. Bias: Both betark and the GLM+HAC estimator were essentially unbiased across all scenarios. Differences in inferential performance were driven entirely by variance estimation, not point estimation.
2. Type I Error and Coverage:
   • Under mild positive and oscillatory autocorrelation, both methods performed reasonably, though the GLM+HAC estimator tended to be overly conservative (wider intervals) in oscillatory scenarios.
   • Under high persistent autocorrelation, the GLM+HAC estimator suffered severe inflation in Type I error (exceeding 60% in short series with AR(3) processes) and poor coverage.
   • betark provided better-calibrated inference in most scenarios. While betark also exhibited elevated Type I error under highly persistent AR(3) processes, the inflation was consistently and substantially lower than the GLM+HAC comparator (e.g., 43.0% vs. 65.4% for AR(3) at $T=100$).
3. Standard Error Calibration: betark's standard error ratios were closer to the ideal value of 1.0. The GLM+HAC estimator systematically underestimated standard errors (ratios as low as 0.22) under high persistence, leading to the observed power inflation and Type I error excess.
4. Robustness: Misspecifying the AR order by a single lag had only a modest effect on betark's performance (1–2 percentage point difference in coverage), suggesting the method is robust to moderate lag-order uncertainty. Sensitivity analyses confirmed that findings held across varying starting means and pre-intervention trends.

Significance and Claims
The paper claims to provide the first systematic evaluation of the finite-sample inferential properties of a joint conditional maximum likelihood Beta-AR($k$) model for ITSA. The primary contribution is the demonstration that betark offers better-calibrated inference than the quasi-binomial GLM with Newey–West HAC errors, particularly under conditions of high persistence and higher-order autoregression where the consequences of estimator choice are most severe.

The authors emphasize that the "power advantage" of the GLM+HAC estimator in high-persistence scenarios is an artifact of anticonservative standard errors (inflated Type I error) rather than genuine statistical efficiency. While betark represents a significant improvement, the study modestly acknowledges that even betark does not fully resolve the inferential challenge under the most demanding conditions (AR(3) high persistence) at moderate sample sizes ($T=400$), where Type I error remains elevated. The paper concludes that researchers analyzing bounded outcomes with suspected high-order serial dependence should characterize the autocorrelation structure and exercise caution, noting that betark is preferable but not a complete solution for all extreme scenarios.