Pith. sign in

REVIEW 4 major objections 5 minor 22 references

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

T0 review · 4 major / 5 minor · reviewed 2026-07-10 · grok-4.5

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read A joint beta-AR likelihood gives better-calibrated ITSA inference for proportions than GLM with HAC standard errors, especially under high persistence.

desk verdict Useful Stata package and honest Monte Carlo for beta-AR ITSA; main claim is real under the matching DGP, but partly measures correct specification rather than robustness. read the letter →

arxiv 2607.07914 v1 pith:3GAYFUCE submitted 2026-07-08 stat.ME

Ariel Linden This is my paper
classification stat.ME MSC 62J1262M1062P10
keywords interruptedtimeseriesanalysisproportionoutcomesbetaregressionautoregressiveerrorsjointconditionalmaximumlikelihoodNewey–WeststandardTypeIerrorsimulationstudy
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Interrupted time series studies of rates and proportions often ignore either the bounded support of the outcome or serial correlation, or both. This paper introduces betark, a joint conditional maximum-likelihood estimator that puts a beta density for the proportion and an AR(k) error structure into one closed-form likelihood, so the reported standard errors already reflect autocorrelation. In a large Monte Carlo study of single-group ITSA designs, both betark and the common quasi-binomial GLM with Newey–West HAC standard errors are essentially unbiased for the slope-change effect. Betark’s confidence intervals and Type I error rates are better calibrated in most cells, with the biggest gap under highly persistent autocorrelation, where GLM+HAC Type I error exceeds 60 percent for short series. Modest lag-order misspecification and changes in starting mean or pre-trend cost little; the remaining weakness is that even betark stays anticonservative for high-persistence AR(3) series of length up to 400.

What carries the argument

The beta-AR(k) recursive substitution of Rocha and Cribari-Neto (and its later extension): the conditional mean on the logit scale is written as the usual ITSA linear predictor plus a sum of AR terms built from realized link-transformed lagged outcomes, yielding a closed-form conditional beta likelihood that is maximized jointly for the mean, precision, and AR coefficients.

What would settle it

Re-run the same Monte Carlo grid but generate the series from a different serial structure (for example ARMA or long-memory) or from a non-beta conditional distribution that still has mean on (0,1); if betark’s coverage and Type I error advantage over GLM+HAC disappears or reverses, the calibration claim does not generalize beyond the correctly specified case.

Watch

Extended reading notes

Core claim

In single-group interrupted time series analysis of beta-distributed proportions, a joint conditional maximum-likelihood beta-AR(k) estimator produces better-calibrated inference for the slope-change parameter than a quasi-binomial GLM with Newey–West HAC standard errors across most of the simulated designs, with the largest advantage under highly persistent autocorrelation where the HAC Type I error exceeds 60 percent at T=100. Both estimators are essentially unbiased; the difference is entirely in variance calibration.

Load-bearing premise

The simulation data are generated from exactly the same beta-AR(k) recursive model that betark estimates, so the method is correctly specified by construction while the comparator is only a quasi-likelihood mean model plus a nonparametric variance fix.

Editorial extensions

If this is right

  • Under mild or oscillatory autocorrelation both estimators are usable, so series length matters more than the choice of method.
  • Under high-persistence serial dependence betark is clearly preferable at every length examined, and the gap widens with AR order.
  • Analysts need not identify the exact lag order: misspecifying by one lag costs only one to two percentage points of coverage.
  • Even with betark, nominal p-values and intervals for high-persistence AR(3) series of length up to 400 should be treated with caution.
  • Future work should check multiple-group, panel, and longer-series designs under the same joint likelihood.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The practical recommendation is condition-dependent: characterize persistence first, then choose the estimator; high-persistence health-rate series are exactly where the joint likelihood pays off.
  • Because the DGP matches betark, the reported advantage is a lower bound on what a correctly specified parametric model can achieve relative to nonparametric HAC; real-data gains will be smaller when the AR form is only approximate.
  • A size-adjusted power comparison (thresholding each method to its own empirical null) would likely erase the raw power advantage of GLM+HAC and leave betark competitive on both size and power.
  • Extending the same recursive construction to panel or multi-intervention ITSA would let the method inherit the multiple-group logic already used for continuous Prais–Winsten estimators.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper introduces betark, a Stata joint conditional MLE for beta regression with AR(k) errors via recursive substitution (Rocha & Cribari-Neto; Ferreira et al.), and evaluates its finite-sample performance in single-group ITSA of proportion/rate outcomes against a quasi-binomial GLM with Newey–West HAC SEs. Across AR(1)–AR(3) scenarios, T ∈ {100,200,400}, and four slope-change effect sizes (2000 replications per cell), both estimators are essentially unbiased; betark yields better coverage, SE ratios, and Type I error in most cells, with the largest gains under high-persistence AR, where GLM+HAC Type I error exceeds 60% at T=100. Lag-order misspecification by ±1 and sensitivity to μ0 and pre-trend have modest effects. The authors candidly report that betark’s own Type I error remains elevated under AR(3) high persistence even at T=400.

Significance. If the comparative calibration results hold under the designs that matter in practice, this is a useful methods contribution: it supplies a joint likelihood estimator that addresses both bounded support and AR errors—two misspecifications that applied ITSA of rates/proportions often treats separately—and documents when the default fractional-GLM+HAC pipeline is badly anticonservative. Strengths include a transparent Monte Carlo design (108 primary cells, 2000 replications, stationarity-checked AR scenarios), external performance metrics (bias, SE ratio, coverage, Type I, power), honest reporting of residual Type I inflation for betark under AR(3) high persistence (Table 6), careful interpretation of raw power as contaminated by SE underestimation, and a usable Stata implementation. The work is incremental relative to the recursive beta-AR literature but fills a clear applied gap for ITSA.

major comments (4)
  1. [§2.5.1, eqs. (11)–(12); §2.5.4] §2.5.1 (eqs. 11–12) and §2.5.4: The Monte Carlo DGP is exactly the recursive beta-AR(k) conditional model whose joint likelihood betark maximizes (§2.2, eq. 6), while the comparator is a quasi-likelihood mean model plus nonparametric HAC (§2.3). The largest reported advantages (e.g., Type I 30.8% vs 60.2% under AR(2) high-persistent at T=100; SE ratios ~0.55–0.77 vs ~0.22–0.23 under AR(3) high-persistent, Tables 6 and 8) therefore partly measure the benefit of correct parametric specification rather than robustness when the true density or dependence form is not beta-AR(k). This is load-bearing for the Abstract claim and the §5.5 recommendation. Either add at least one alternative DGP family (e.g., logit-scale AR with non-beta errors, additive AR on the proportion scale truncated/reflected, or a different link/precision structure) or reframe the central comparative claim as conditional o
  2. [§2.5.6; Appendix B; Abstract] §2.5.6 and Appendix B: The only misspecification check varies lag order by ±1 under mild-positive AR and a single 50% effect. That does not support the Abstract’s broader robustness language (“Misspecifying the AR order by one lag … had only modest effects”) as evidence that betark remains well calibrated under realistic process misspecification. Expand the misspecification design (at least under high-persistent Scenario 3, and ideally under a non-beta or non-recursive dependence structure) or narrow the Abstract/§5.4 claims to lag-order robustness under mild AR only.
  3. [§5.5; Table 6] §5.5 and Table 6: The practical recommendation that “under high persistent autocorrelation … betark is clearly preferable … at every series length examined” is only partly supported. Preferability over GLM+HAC is clear in the reported cells, but betark’s own Type I error under AR(3) high persistent remains 43%/34%/23% at T=100/200/400—materially anticonservative even at the largest T. The recommendation should state more sharply that neither estimator delivers reliable nominal inference under that condition at the lengths studied, and that “preferable” does not mean “adequately calibrated.” Align §5.5 with the more cautious language already in §5.2 and the Conclusion.
  4. [§2.3; §2.5.4] §2.3: Matching the Newey–West bandwidth exactly to the true AR order k is an idealized comparator (acknowledged as such). In applied work the bandwidth is chosen by rule of thumb or data-driven criteria and is often misspecified relative to a structured AR(k). At least one sensitivity with automatic lag selection (or a fixed rule independent of true k) is needed before concluding that GLM+HAC’s failure under high persistence is not partly an artifact of the bandwidth protocol used here; otherwise state explicitly that the comparison is to an oracle-bandwidth HAC.
minor comments (5)
  1. [Figure 1] Figure 1 caption says the black line is the “betark-fitted trajectory” but the text states fitted values shown are from glm+HAC and that betark fits are indistinguishable. Align caption and body.
  2. [Table 1] Table 1 lists performance measures with a dangling “[16]” after the header; integrate the Burton et al. citation cleanly.
  3. [§3.1] §3.1 non-monotonic power for betark under high persistence is well explained later; a one-sentence forward pointer in §3.1 would help readers who stop at the power tables.
  4. [§1 / §5] Several related Linden preprints/papers on Prais–Winsten/Newey–West ITSA are cited heavily; a short paragraph situating betark relative to those continuous-outcome results would help readers who know that line of work.
  5. [§2.2; §2.5.4] Clarify whether precision ϕ is estimated freely in every betark fit (as the joint likelihood suggests) or held fixed; the simulation description emphasizes mean/AR structure more than precision estimation.

Circularity Check

0 steps flagged · score 1.0 of 10

No derivation circularity; only disclosed correct-specification DGP that structurally favors the parametric joint MLE over the quasi-likelihood+HAC comparator.

full rationale

This is a Monte Carlo methods paper, not a first-principles derivation that equates a claimed prediction to its inputs. The recursive beta-AR(k) likelihood (eqs. 5–6) is taken from external sources (Rocha & Cribari-Neto 2009; Ferreira et al. 2015) and implemented as joint CML; the simulation then generates data from exactly that process (eqs. 11–12, §2.5.1) and reports empirical bias, SE ratio, coverage, Type I, and power against a known true slope-change parameter. Those Monte Carlo frequencies are external benchmarks, not tautologies: coverage and Type I are not forced to 95%/5% by construction, and the paper itself documents residual anticonservatism for betark under high-persistent AR(3) even at T=400. The design therefore evaluates correct-specification finite-sample behavior (standard practice) and discloses the match; it does not rename a fitted constant as a prediction, import a uniqueness theorem from self-citation, or reduce Eq. X to Eq. Y by definition. Self-citations ([3],[5],[6],[11],[15]) supply prior ITSA software and continuous-outcome simulation templates; none is load-bearing for the comparative calibration claim. Misspecification is limited to lag order ±1 under mild AR (App. B), which is a scope limitation, not circularity. Score 1 only for the mild, disclosed design favoritism that makes the largest reported gains partly a correct-vs-quasi-likelihood contrast rather than pure robustness; no steps meet the enumerated circular patterns.

Assumptions & free parameters 5 free parameters · 5 assumptions · 1 invented entities

The central comparative claim rests on standard beta regression and a published recursive AR substitution, plus simulation design choices (starting mean, dispersion, AR coefficient vectors, series lengths, effect sizes) that define the Monte Carlo world. No new physical entities are postulated; betark is software for an existing likelihood. The load-bearing modeling axioms are the beta conditional density, the recursive link-scale AR feedback, and stationarity of the chosen AR processes.

free parameters (5)
  • cv (dispersion via coefficient of variation)
    Fixed at 0.05 throughout the primary design (ϕ≈3599 at μ0=0.10); chosen from preliminary bias checks (Appendix A), not estimated from external data, and not fully crossed with AR order and T.
  • starting mean μ0
    Fixed at 0.10 in the primary grid; sensitivity later varies {0.05,0.10,0.30,0.50}. Primary ranking of methods depends on this design choice for base-rate outcomes.
  • AR scenario coefficient vectors (Scenarios 1–3 for k=1,2,3)
    Hand-chosen mild, oscillatory, and high-persistent ρ vectors following prior ITSA simulation designs; they define persistence and thus drive the main performance gaps.
  • post-intervention percentage effect sizes (0%, 25%, 50%, 100%)
    Design anchors for solving logit-scale slopes; define power and non-null calibration cells rather than being estimated from data.
  • series lengths T ∈ {100,200,400}
    Design grid that bounds the asymptotic claims; AR(3) high-persistence calibration is incomplete even at the largest T.
assumptions (5)
  • domain assumption Conditional on F_{t-1}, y_t ~ Beta(μ_t, ϕ_t) with logit mean and (here) constant precision (eqs. 1–3).
    Standard Ferrari–Cribari-Neto beta regression; required for the joint likelihood and for the DGP.
  • domain assumption AR(k) dependence enters via recursive substitution on the link scale using realized g(y_{t-i}) (eqs. 4–5; Rocha & Cribari-Neto; Ferreira et al.).
    Core structural assumption of betark; also used as the simulation truth.
  • standard math All simulated AR processes are stationary (companion-matrix eigenvalues inside the unit circle).
    Stated as numerically confirmed before simulation (§2.5.2); excludes near-unit-root or nonstationary regimes.
  • ad hoc to paper Single-group ITSA mean structure with midpoint intervention, δ_step=0, primary δ_pre=0, treatment effect = post–pre slope change (eqs. 10–11).
    Design restriction that defines the estimand and excludes multi-group, multi-intervention, and seasonal structures.
  • domain assumption Newey–West HAC with Bartlett weights and bandwidth L set to true AR order k is a fair idealized comparator.
    §2.3; gives the semiparametric alternative its best lag order while still not modeling beta precision or parametric AR.
invented entities (1)
  • betark (Stata joint CML beta-AR(k) estimator) independent evidence
    purpose: Package the recursive conditional beta likelihood and report joint asymptotic SEs for ITSA of proportions.
    Software implementation of prior theory rather than a new scientific object; independent handle is the SSC module and finite-sample behavior under the stated DGP.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Beta Regression with Autoregressive Errors for Interrupted Time Series Analysis of Proportion and Rate Outcomes: A Simulation Study." pith.science (2026). https://pith.science/paper/3GAYFUCE

@misc{pith2026260707914,
  author       = {Pith},
  title        = {Pith review of: Beta Regression with Autoregressive Errors for Interrupted Time Series Analysis of Proportion and Rate Outcomes: A Simulation Study},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/3GAYFUCE}},
  note         = {Machine review of arXiv:2607.07914}
}
read the original abstract

Interrupted time series analyses (ITSA) of proportion and rate outcomes are frequently estimated using ordinary least squares regression despite the bounded nature of these outcomes. When methods appropriate for bounded outcomes are used, the standard approach is a quasi-likelihood generalized linear model (GLM) with heteroskedasticity- and autocorrelation-consistent (HAC) standard errors. However, no existing estimator jointly models the beta-distributed conditional density and autoregressive (AR) error structure. We introduce betark, a Stata implementation of a joint conditional maximum likelihood estimator for beta regression with AR(k) errors based on a recursive substitution that yields a closed-form conditional beta likelihood with autoregressive dependence of arbitrary order. Unlike two-stage approaches, betark jointly estimates the mean, precision, and AR(k) coefficients in a single likelihood, so reported standard errors directly account for autocorrelation without separate correction. A Monte Carlo study compared betark with a quasi-binomial GLM using Newey-West HAC standard errors across AR(1)-AR(3) processes, three series lengths, and four effect sizes in a single-group ITSA design. Both methods were essentially unbiased, but betark produced better-calibrated inference than GLM+HAC in most scenarios, with the largest gains under highly persistent autocorrelation, where GLM+HAC Type I error exceeded 60% for short series. Misspecifying the AR order by one lag and varying the starting mean and pre-intervention trend had only modest effects on performance. However, betark's own Type I error remained elevated under highly persistent AR(3) processes even for the longest series examined.

Figures

Figures reproduced from arXiv: 2607.07914 by the authors.

Figure 1
Figure 1. Proportion of patients above the glycemic threshold (observed, gray circles) with [PITH_FULL_IMAGE:figures/full_fig_p028_1.png] view at source ↗

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

22 extracted references · 22 canonical work pages

  1. [1]

    Campbell and Julian C

    Donald T. Campbell and Julian C. Stanley.Experimental and Quasi-Experimental Designs for Research. Rand McNally, Chicago, 1966

    work page 1966

  2. [2]

    Shadish, Thomas D

    William R. Shadish, Thomas D. Cook, and Donald T. Campbell.Experimental and Quasi- Experimental Designs for Generalized Causal Inference. Houghton Mifflin, Boston, 2002

    work page 2002

  3. [3]

    Conducting interrupted time-series analysis for single- and multiple-group com- parisons.Stata Journal, 15(2):480–500, 2015

    Ariel Linden. Conducting interrupted time-series analysis for single- and multiple-group com- parisons.Stata Journal, 15(2):480–500, 2015

    work page 2015

  4. [4]

    Beta regression for modelling rates and proportions

    Silvia Ferrari and Francisco Cribari-Neto. Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7):799–815, 2004

    work page 2004

  5. [5]

    Ariel Linden. Adjustment for autocorrelation in multiple-group (controlled) interrupted time series analysis and its effect on power: a simulation study of the Newey-West and Prais-Winsten methods, 2026. Preprint, Research Square

    work page 2026

  6. [6]

    Multiple-group (controlled) interrupted time series analysis with higher-order au- toregressive errors: A simulation study comparing Newey–West and Prais–Winsten methods

    Ariel Linden. Multiple-group (controlled) interrupted time series analysis with higher-order au- toregressive errors: A simulation study comparing Newey–West and Prais–Winsten methods. Journal of Statistical Computation and Simulation, 2026

    work page 2026

  7. [7]

    Turner, Andrew B

    Simon L. Turner, Andrew B. Forbes, Amalia Karahalios, Monica Taljaard, and Joanne E. McKenzie. Evaluation of statistical methods used in the analysis of interrupted time series studies: a simulation study.BMC Medical Research Methodology, 21:181, 2021

    work page 2021

  8. [8]

    Christian Bottomley, Mildred Ooko, Antonio Gasparrini, and Ruth H. Keogh. In praise of Prais–Winsten: an evaluation of methods used to account for autocorrelation in interrupted time series.Statistics in Medicine, 42(8):1277–1288, 2023

    work page 2023

Show all 22 references
  1. [9]

    Rocha and Francisco Cribari-Neto

    Andr´ ea V. Rocha and Francisco Cribari-Neto. Beta autoregressive moving average models. TEST, 18(3):529–545, 2009

    2009

  2. [10]

    Figueroa-Z´ u˜ niga, and M´ ario de Castro

    Guillermo Ferreira, Jorge I. Figueroa-Z´ u˜ niga, and M´ ario de Castro. Partially linear beta regression model with autoregressive errors.TEST, 24(4):752–775, 2015

    2015

  3. [11]

    BETARK: Stata module for computing beta regression with autoregressive- corrected errors for proportion outcomes, by joint conditional maximum likelihood, 2026

    Ariel Linden. BETARK: Stata module for computing beta regression with autoregressive- corrected errors for proportion outcomes, by joint conditional maximum likelihood, 2026. Sta- tistical Software Components s459773, Boston College Department of Economics

    2026

  4. [12]

    Papke and Jeffrey M

    Leslie E. Papke and Jeffrey M. Wooldridge. Econometric methods for fractional response vari- ables with an application to 401(k) plan participation rates.Journal of Applied Econometrics, 11(6):619–632, 1996. 16

    1996

  5. [13]

    Newey and Kenneth D

    Whitney K. Newey and Kenneth D. West. A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix.Econometrica, 55:703–708, 1987

    1987

  6. [14]

    Newey and Kenneth D

    Whitney K. Newey and Kenneth D. West. Automatic lag selection in covariance matrix estimation.Review of Economic Studies, 61(4):631–653, 1994

    1994

  7. [15]

    Extending Prais–Winsten regression to panel data with higher-order autoregres- sive errors: A simulation study, 2026

    Ariel Linden. Extending Prais–Winsten regression to panel data with higher-order autoregres- sive errors: A simulation study, 2026. Preprint, arXiv:2606.12596

    2026 arXiv

  8. [16]

    Altman, Patrick Royston, and Roger L

    Andrea Burton, Douglas G. Altman, Patrick Royston, and Roger L. Holder. The design of simulation studies in medical statistics.Statistics in Medicine, 25:4279–4292, 2006

    2006

  9. [17]

    Medicare disease management in policy context.Health Care Financing Review, 29:1–11, 2008

    Ariel Linden and Julia Adler-Milstein. Medicare disease management in policy context.Health Care Financing Review, 29:1–11, 2008

    2008

  10. [18]

    Evaluation methods in disease management: determining program effectiveness, October 2003

    Ariel Linden, John Adams, and Nancy Roberts. Evaluation methods in disease management: determining program effectiveness, October 2003. Position Paper for the Disease Management Association of America (DMAA)

    2003

  11. [19]

    Kullgren, Erin Krupka, Abigail Schachter, et al

    Jeffrey T. Kullgren, Erin Krupka, Abigail Schachter, et al. Precommitting to choose wisely about low-value services: a stepped wedge cluster randomised trial.BMJ Quality and Safety, 27:355–364, 2018

    2018

  12. [20]

    Biuso, Susan Butterworth, and Ariel Linden

    Thomas J. Biuso, Susan Butterworth, and Ariel Linden. A conceptual framework for tar- geting prediabetes with lifestyle, clinical and behavioral management interventions.Disease Management, 10(1):6–15, 2007

    2007

  13. [21]

    A user’s guide to the disease management literature: rec- ommendations for reporting and assessing program outcomes.American Journal of Managed Care, 11:113–120, 2005

    Ariel Linden and Nancy Roberts. A user’s guide to the disease management literature: rec- ommendations for reporting and assessing program outcomes.American Journal of Managed Care, 11:113–120, 2005. 17 Tables Table 1: Simulation design and inputs. Parameter Input values Regre...

    2005

  14. [22]

    Because the data-generating process is identical across estimators for a given AR order, the fitted trajectories shown here are fromglm+HAC;betarkfitted values are indistinguishable and are not shown separately. 28

Pith tools

Reviewed July 10, 2026 · model on record in the stance chip above.