The preprint arXiv:2307.00929 (Gomtsyan, 2023) proposes a three-stage variable-selection procedure for negative binomial GLARMA regression — a count time series model whose conditional mean carries autoregressive-moving-average serial dependence — combining (1) Newton-Raphson estimation of the serial-dependence parameters, (2) ℓ1\ell_1-penalized quadratic approximation of the negative binomial log-likelihood for sparse regression coefficients, and (3) maximum-likelihood re-estimation on the selected support, all wrapped in stability selection. The procedure is supported only by simulations and one short RNA-Seq application. The preprint leaves this problem open: establish (or disprove) asymptotic guarantees — consistency, convergence rates, and selection/sign consistency — for the actual plug-in three-stage estimator in the high-dimensional dependent-count regime, and derive a principled, data-driven rule for the stability-selection threshold and variant with error control adapted to serially dependent counts. No work establishing such theory is known; the closest prior results cover only the ARMA-part estimator in a specific case, or the independent Poisson case without serial dependence or overdispersion.

Background

A time series of counts consists of nonnegative-integer observations, such as daily symptom counts or RNA-Seq read counts of genes tracked over successive samples. Two standard model families exist: parameter-driven models, in which serial dependence enters through a latent stochastic process and likelihood evaluation requires computationally expensive Monte Carlo methods, and observation-driven models, in which the conditional mean at time tt depends on past observations and past residuals through explicit recursions and hence is cheap to evaluate. The GLARMA (generalized linear autoregressive moving average) model of Davis, Dunsmuir and Streett, implemented in the glarma package, is the standard observation-driven choice: the linear predictor is Wtγ+Xt⊤βW_t\gamma + X_t^\top\beta, where WtW_t is a serial latent state built from past outcomes and past residuals via ARMA-type recursions with parameter γ\gamma, XtX_t are covariates with coefficient vector β\beta, and the observations follow a Poisson or negative binomial (overdispersed Poisson) distribution whose mean is a link function of the linear predictor; the negative binomial version adds a dispersion parameter α\alpha. When many covariates are recorded (pp comparable to or larger than the series length nn), one wants to select the sparse set of relevant covariates. Lévy-Leduc, Ouadah and Sansonnet (arXiv:1907.07085) and Gomtsyan, Lévy-Leduc, Ouadah, Sansonnet and Blein (Statistics, 2022) introduced a two-stage variable-selection method for Poisson GLARMA models — first estimating the serial-dependence parameter γ\gamma by least squares on working residuals, then applying an ℓ1\ell_1 (lasso-type) penalty to select β\beta — but established consistency only for the ARMA-part estimator γ^\hat\gamma, and only in a specific case, with no selection-consistency theory for β\beta. Gomtsyan, Lévy-Leduc, Ouadah and Sansonnet (arXiv:2303.14020) proved sign consistency of the ℓ1\ell_1 quadratic-approximation strategy for independent Poisson observations, a setting with no serial dependence, no negative-binomial overdispersion, and no plug-in of estimated nuisance parameters. The preprint (Gomtsyan, arXiv:2307.00929) extends this line to negative binomial GLARMA regression: a three-stage estimator that (i) estimates γ\gamma by Newton-Raphson maximization of the conditional negative binomial log-likelihood over γ\gamma, starting from initial values fitted by a GLM that ignores the serial-dependence part (the ARMA-type working-residual state), (ii) estimates β\beta by maximizing an ℓ1\ell_1-penalized quadratic (Taylor) approximation of the negative binomial log-likelihood, and (iii) re-estimates the dispersion and coefficients by maximum likelihood of a negative binomial GLM restricted to the selected support. The estimator is wrapped in stability selection (Meinshausen and Bühlmann, 2010): the transformed response Y∈Rp+1\mathcal{Y}\in\mathbb{R}^{p+1} of the quadratic approximation is split into random subsamples of size (p+1)/2(p+1)/2 (half its length — a subsample of the p+1p+1 transformed coordinates, not of the nn time points), and covariates selected in at least a threshold fraction of subsamples are retained. Stability selection comes with false-selection error bounds only under independence assumptions, which serially dependent count data violate. The preprint evaluates the method solely by simulation (true/false positive rates; Fourier-basis covariates, p=100p=100, ARMA orders q∈{1,2}q\in\{1,2\}, series lengths n∈{150,250,500,1000}n\in\{150,250,500,1000\}) and one RNA-Seq application with n=15n=15 time points, choosing the stability-selection threshold empirically (e.g. threshold 0.7 picked by simulation), and proves no asymptotic property of any stage.

Problem Statement

For the negative binomial GLARMA regression model and the three-stage estimation procedure of Gomtsyan (arXiv:2307.00929) — stage 1: the serial-dependence (ARMA) coefficients γ^\hat\gamma, obtained by Newton-Raphson maximization of the negative binomial log-likelihood over γ\gamma, starting from initial values fitted by a GLM that ignores the serial-dependence part; stage 2: β^\hat\beta by maximizing an ℓ1\ell_1-penalized quadratic (Taylor) approximation of the negative binomial log-likelihood; stage 3: the dispersion α^\hat\alpha and the coefficients re-estimated by negative binomial GLM maximum likelihood restricted to the selected support; overall wrapped in stability selection applied to the transformed response Y∈Rp+1\mathcal{Y}\in\mathbb{R}^{p+1} of the quadratic approximation, with subsamples of size (p+1)/2(p+1)/2 (the preprint also describes an iterated variant in which these stages are repeated until the ARMA-coefficient iterates stabilise; results for either the basic or the iterated variant are of interest) — establish the statistical theory that the preprint proves nowhere, or prove it unattainable. Specifically, in an explicitly stated asymptotic framework where the number of covariates pp grows with the series length nn and the true coefficient vector β⋆\beta^\star is sparse: (a) prove or refute consistency, convergence rates, and selection (sign) consistency of the recovered support of β⋆\beta^\star for the actual plug-in three-stage estimator — with the stage-1 error in γ^\hat\gamma and the stage-3 re-estimation propagated through the analysis, under regularity conditions that account for the serial dependence of the count observations and the negative binomial overdispersion — including what stability selection adds or costs under serially dependent data; and (b) derive from such analysis (or independently) a principled, data-driven criterion for the two tuning choices the preprint makes ad hoc — the selection-frequency threshold, and the choice between the ss_min variant (the penalty parameter λ\lambda of the ℓ1\ell_1 step fixed at the smallest value of the regularization grid of the standard R lasso implementation glmnet) and the ss_cv variant (λ\lambda chosen by cross-validation) — with false-selection error control adapted to dependent counts, replacing the current empirical threshold selection. Part (a) is the central target; part (b) is an additional target. The empirical aspects the preprint leaves uncharacterised (scalability to p≫100p\gg100, q>2q>2 or q=∞q=\infty recursions; selection uncertainty for n≈15n\approx15; formal comparison with full-likelihood or MCMC alternatives) are recorded as separate acknowledged open aspects and are not part of this problem's completion boundary.

References

  1. Marina Gomtsyan (2023), 'Variable selection in a specific regression time series of counts', arXiv:2307.00929, DOI 10.48550/arXiv.2307.00929, https://arxiv.org/abs/2307.00929
  2. Céline Lévy-Leduc, Sarah Ouadah, Laure Sansonnet (2019), 'Variable selection in sparse high-dimensional GLARMA models', arXiv:1907.07085, DOI 10.48550/arXiv.1907.07085, https://arxiv.org/abs/1907.07085
  3. Marina Gomtsyan, Céline Lévy-Leduc, Sarah Ouadah, Laure Sansonnet, Thomas Blein (2022), 'Variable selection in sparse GLARMA models', Statistics 56(4), 755-784, DOI 10.1080/02331888.2022.2090943, https://doi.org/10.1080/02331888.2022.2090943
  4. Marina Gomtsyan, Céline Lévy-Leduc, Sarah Ouadah, Laure Sansonnet (2023), 'Sign-consistent estimation in a sparse Poisson model', arXiv:2303.14020, DOI 10.48550/arXiv.2303.14020, https://arxiv.org/abs/2303.14020
  5. Nicolai Meinshausen, Peter Bühlmann (2010), 'Stability Selection', Journal of the Royal Statistical Society Series B: Statistical Methodology 72(4), 417-473, DOI 10.1111/j.1467-9868.2010.00740.x, https://doi.org/10.1111/j.1467-9868.2010.00740.x
  6. William T. M. Dunsmuir, David J. Scott (2015), 'The glarma Package for Observation-Driven Time Series Regression of Counts', Journal of Statistical Software 67(7), DOI 10.18637/jss.v067.i07, https://doi.org/10.18637/jss.v067.i07
  7. Guilherme Pumi, Ana Julia Alves Câmara (2026), 'Bayes Estimation of GLARMA Models With Applications', arXiv:2607.03586, DOI 10.48550/arXiv.2607.03586, https://arxiv.org/abs/2607.03586