{"id":"bd966d10-d888-4b5e-bd28-e6dc5f1fbf23","arxiv_id":"1908.07477","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"high","formal_verification":"none","parameter_count":3,"one_line_summary":"The paper presents a ridge-regularized EM algorithm and a component-based regularized EM algorithm for GLMMs with an AR(1) random effect, with simulation results only for the ridge version.","lead":"The paper proposes two regularized EM algorithms for fitting generalized linear mixed models on panel data with an autoregressive time effect. It demonstrates the L2-penalized version on simulated Poisson data, while the component-based high-dimensional version is only sketched.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"EM M-step for the AR(1) parameter ρ is never derived or verified; if that update is incorrect, the central estimates are invalid even when the linearisation is exact.","rationale":"The reader identifies the linearisation approximation as the weakest assumption. I agree that approximation is unvalidated, but the more fundamental and specific gap is the EM update for the AR(1) parameter ρ, because it is the paper's novelty and is asserted without derivation. Even if the linearisation were exact (e.g., Gaussian identity link), a misspecified M-step would break the algorithm. This concern is concrete and testable, and the proposed Gaussian test isolates it from the linearisation issue. I therefore partially agree with the reader: their weakest_assumption is about the linearisation, while I point to the under-specified EM M-step; both support the same conditional verdict. The manuscript should be accepted only after the M-step is derived and validated, along with the broader simulations the reader requested.","tokens_in":4399,"tokens_out":7571,"duration_ms":85746,"concrete_test":"Run Algorithm 1 on data generated from a Gaussian linear mixed model with identity link, response y=Xβ+U1ξ1+U2ξ2+e, ξ2 AR(1), so the linearisation in Section 3 is exact by construction. Compare the algorithm's estimates of (β,σ1^2,σ2^2,ρ) against a reference fit (e.g., nlme::lme with corAR1 or a direct maximum likelihood/REML routine) over multiple replications. If the proposed estimates have material bias or fail to track the true ρ, the EM M-step is incorrect and the central claim fails; if they match the reference, the concern about the M-step does not land and the remaining risk is the linearisation approximation.","verdict_should_be":"UNCHANGED","load_bearing_attack":"In Algorithm 1 (step 2.c), the update θ[t+1]=argmax_θ Qpen(θ,θ[t]) with θ=(β,σ1^2,σ2^2,ρ) is asserted but not derived. The complete-data log-likelihood of the linearised model M depends on ρ through the inverse covariance matrix D^{-1} of the AR(1) random effect ξ2, which is tridiagonal with stationary edge corrections. The paper provides neither the conditional moments E[ξ2|z,θ[t]] and E[ξ2ξ2^T|z,θ[t]] nor the resulting M-step equations for σ2^2 and ρ. The M-step for ρ is especially delicate: it is not the sample autocorrelation of the conditional mean, because it must account for conditional variance and for the stationarity constraint |ρ|<1. Without this derivation, the Figure 1 results cannot be reproduced or checked, and the central claim that ρ is estimated reliably is unsupported. Note also that step 2.b requires the hat matrix S[t]_λ, which is never defined; this is another under-specified component of the same algorithm. The figure alone cannot establish that the implemented updates actually maximize Qpen.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a regularised EM algorithm for generalised linear mixed models with a random intercept per individual and a stationary AR(1) time-specific random effect, in a balanced panel-data setting. The algorithm alternates a Schall-type linearisation step with an L2-penalised EM step whose ridge parameter is chosen by GCV; a supervised component-based extension for high-dimensional predictors is sketched in Section 3.2. The numerical section reports convergence trajectories and MSE curves for a Poisson panel simulation with N=10 and T varying from 10 to 100, together with boxplots of estimated autocorrelations.","tokens_in":4625,"tokens_out":4640,"duration_ms":46022,"significance":"If the algorithm performs as claimed, it would fill a genuine gap, since existing regularised GLMM methods such as glmmLasso do not accommodate autocorrelated random effects. The paper is clearly written and the model formulation is natural for panel data. However, the significance is currently limited by the absence of derivations for the M-step, the lack of any comparison to simpler baselines, and the absence of numerical results for the high-dimensional component-based method. The paper does not ship code or machine-checked proofs, so the reported simulation cannot be independently checked from the text alone.","major_comments":[{"comment":"The M-step is stated only as the argmax of Qpen, with no derivation of the update equations. Since the complete-data log-likelihood of the linearised model M depends on rho through the inverse of the AR(1) covariance matrix of xi_2, which is tridiagonal with stationary edge corrections, the conditional moments E[xi_2|z,theta^(t)] and E[xi_2 xi_2^T|z,theta^(t)], and the resulting updates for sigma_2^2 and rho, are nontrivial and must be given. In particular, the rho update is not simply the sample autocorrelation of the conditional means, because the conditional variance and the stationarity constraint |rho|<1 must be accounted for. Without these equations, the central claim that beta, sigma_2^2, and rho are reliably estimated cannot be verified or reproduced.","section":"Section 3.1, Algorithm 1 (step 2.c)"},{"comment":"The hat-matrix S[t]_lambda is said to satisfy zhat[t]=S[t]_lambda z[t], but no explicit expression is provided. In a penalized linear mixed model this matrix depends on the current variance components and on lambda, and the GCV criterion also requires its trace; without this definition the algorithm is under-specified. The subsequent statement that lambda is chosen by GCV is therefore not implementable from the text.","section":"Section 3.1, Algorithm 1 (step 2.b)"},{"comment":"The only numerical evidence is a single Poisson setting with N=10 and T varying from 10 to 100. There is no comparison with an unregularised EM, with Schall's algorithm, or with any existing mixed-model package, no Monte Carlo standard errors, and no reported bias of the point estimates. The caption's claim that 'MSEs of beta, sigma_2^2 and rho decrease towards zero' is only indicative, and the 'good estimation' claim is not supported. Moreover, the high-dimensional component-based method in Section 3.2 is not evaluated at all.","section":"Section 4, Figure 1"},{"comment":"The working-variable approximation g(y_i) approx z_i = g(mu_i)+(y_i-mu_i)g'(mu_i) is used without checking its accuracy for the proposed model. With Poisson counts of small mean, strong autocorrelation, or a non-canonical link, this linearisation can introduce a substantial bias that the EM loop does not correct for. The paper neither reports a diagnostic (e.g., a simulation with small means) nor gives conditions under which the approximation error is negligible for the AR(1) random-effect structure.","section":"Section 3, Linearisation step"}],"minor_comments":[{"comment":"The notation involving Gamma[t]^{-1} is ambiguous; it should be written as a quadratic form (z - S z)^T (Gamma[t])^{-1} (z - S z) rather than as a norm with a subscripted inverse matrix.","section":"Section 3.1, GCV formula"},{"comment":"The definitions U1 = I_N tensor 1_T and U2 = 1_N tensor I_T are terse; a sentence explaining that U1 selects the individual-specific effect and U2 the time-specific effect would help the reader.","section":"Section 2, design matrices"},{"comment":"The x-axis of the boxplot of estimated autocorrelations is a long unbroken sequence of numbers; it should be reformatted with labeled tick marks at meaningful values.","section":"Section 4, Figure 1(d)"},{"comment":"The tuning parameters s and l are said to be chosen by cross-validation, but no cross-validation procedure is described; at least a brief explanation of the CV scheme is needed.","section":"Section 3.2"},{"comment":"The paper ends abruptly after the numerical section; a short conclusion discussing limitations and next steps would be helpful.","section":"General"}],"recommendation":"major_revision","confidential_remarks":"The paper is a very condensed research announcement. The authors need to provide the full derivation of the EM updates (especially for rho) and more extensive simulations, including a baseline comparison. The absence of any evaluation of the high-dimensional method is a concern for a journal that expects a complete contribution. The novelty relative to Eliot et al. (2011) is the AR(1) structure, and the paper should clarify exactly what is new beyond that reference."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"First thing you should know: this is a plausible extension of Eliot et al.'s ridge EM to GLMMs with an AR(1) time-specific random effect, but the paper is too under-specified to be usable as it stands. The central M-step for the AR(1) parameter ρ is never derived, and the only simulation has no comparison baselines.\n\nWhat's new: the combination of regularized fixed effects with an autoregressive random effect in a GLMM is not in the cited literature, and the component-based extension for p >> n is a natural direction. The model setup is clear, and the convergence trajectories in Figure 1 show their implementation does stabilize.\n\nThe biggest problem is Algorithm 1. Step 2.c writes θ[t+1] = argmax Qpen, but the conditional moments and update equations for σ2^2 and ρ are not given. For ρ, the M-step is not just the autocorrelation of the conditional mean; it has to account for the posterior variance and the |ρ|<1 constraint. Without that derivation, the simulation results can't be checked or reproduced. Step 2.b also references a hat matrix S_t_lambda that is never defined, so the GCV step is ambiguous.\n\nBeyond that, the evidence is thin. One Poisson simulation with N=10, no baseline methods, no varying N, and no small-count check. The high-dimensional variant has zero numerical results. The linearisation approximation is standard Schall, so I'm not worried in principle, but it should be tested in the settings where the method is actually meant to work.\n\nWho this is for: researchers working on panel GLMMs with serial correlation and many covariates. They would want the full derivations and code. Right now it reads like an extended abstract, not a complete methods paper.\n\nRecommendation: send to peer review, but with the clear expectation of major revision. The problem is worthwhile and the approach is credible, but the missing derivations and thin evidence are too substantial to accept as is.","headline":"Sensible combination of ridge EM and AR(1) random effects, but the key M-step is asserted rather than derived and the evidence is one unbenchmarked simulation.","tokens_in":5151,"tokens_out":2198,"would_cite":false,"duration_ms":23476,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that an L2-penalised EM algorithm, built on a first-order linearisation of the link function, reliably estimates the regression coefficients, variance components, and autoregressive parameter of a GLMM with an AR(1)…","keywords":["Generalized linear mixed models","EM algorithm","Ridge regularization","Autoregressive random effect","Panel data","Generalized cross-validation","Component-based regularization","Poisson regression"],"falsifier":"Run the algorithm on Poisson GLMM data with an AR(1) random effect and small expected counts (mean near 0.5), then check whether the mean squared errors of $\\beta$, $\\sigma_2^2$, and $\\rho$ still decrease to zero as $T$ grows; if they level off or increase, the linearisation step is not adequate.","tokens_in":4214,"feed_emoji":"📈","tokens_out":7663,"duration_ms":70181,"temperature":0.7,"pith_summary":"Generalised linear mixed models for panel data usually assume random effects with independent levels, which ignores the serial correlation of latent time-specific effects. This paper proposes an L2-penalised EM algorithm that fits a GLMM whose time-specific random effect follows a stationary AR(1) process, and a component-based regularised variant for high-dimensional predictors. A simulation in the Poisson case indicates that the algorithm converges in roughly a hundred iterations and that the mean squared errors of the regression coefficients, the temporal variance, and the autoregressive parameter decrease as the number of time points grows. The paper therefore aims to extend regularised likelihood estimation from linear mixed models to GLMMs with autocorrelated random effects.","feed_headline":"Ridge-penalised EM recovers AR(1) random effect in GLMMs","feed_subtitle":"Poisson panel-data simulations show errors for beta, variances, and the autoregressive parameter shrink as T grows.","key_machinery":"The engine of the method is Schall's linearisation step: each observation $y_i$ is replaced by the working variable $z_i = g(\\mu_i) + (y_i - \\mu_i) g'(\\mu_i)$, so the GLMM is approximated by a linear mixed model $z = X\\beta + U_1\\xi_1 + U_2\\xi_2 + e$ with heteroskedastic error variance $\\Gamma = \\operatorname{Diag}([g'(\\mu_i)]^2 \\operatorname{Var}(Y_i|\\xi))$. Estimation then proceeds by an EM algorithm for the penalised complete log-likelihood $L_{\\text{pen}} = L - (\\lambda/2)\\beta^T\\beta$, with $\\lambda$ selected at each iteration by a generalised cross-validation criterion modified for heteroskedastic errors; the M-step updates $\\xi$ by its conditional expectation, and the AR(1) covariance structure of $\\xi_2$ carries the new parameter $\\rho$ into the estimation.","core_discovery":"The central claim is that a regularised expectation-maximisation algorithm, built on a first-order linearisation of the link function, can estimate the fixed regression coefficients and the variance components of a GLMM when the second random effect is a stationary AR(1) process shared by all individuals. The paper derives the penalised complete log-likelihood, chooses the ridge parameter by a generalised cross-validation criterion adapted to heteroskedastic errors, and alternates a linearisation step with an EM step until the parameter updates are stable. Simulation evidence in the canonical Poisson case shows convergence trajectories and decreasing mean squared errors for $\\beta$, $\\sigma_2^2$, and $\\rho$ as $T$ grows, while the MSE for $\\sigma_1^2$ stays flat because $N$ is fixed. For the high-dimensional case, the paper describes a supervised component-based variant.","pith_inferences":["The linearisation assumption would be the first point to test: at small expected counts, the algorithm's bias could be compared against a simulation-based gold standard that does not rely on the working normal approximation.","The same EM scheme should transfer to other autocorrelation structures, such as ARMA(p,q) random effects, by replacing the AR(1) covariance matrix in the M-step; nothing in the derivation appears specific to order one.","The component-based high-dimensional algorithm is described without numerical evidence, so its practical tuning, including the trade-off parameter $s$ and the structural-relevance exponent $l$, remains an open empirical question."],"forward_implications":["For panel data with many redundant covariates, the L2-penalised EM returns estimates of $\\beta$, $\\sigma_1^2$, $\\sigma_2^2$, and $\\rho$ in a single iterative scheme.","The simulation results imply that increasing the time dimension $T$ sharpens the estimates of the regression coefficients, the temporal variance, and the autoregressive parameter, while the individual variance estimate is limited by the number of individuals $N$.","The ridge parameter does not need to be fixed in advance: it is re-estimated at each iteration via a generalised cross-validation criterion adapted to heteroskedastic errors.","In the high-dimensional setting, the same penalised EM objective can be replaced by a component-based criterion that balances model fit against structural relevance of the components."],"supporting_citations":[{"why":"Supplies the linearisation step that converts the GLMM into a working linear mixed model.","marker":"Schall (1991)"},{"why":"Supplies the L2-penalised EM algorithm for longitudinal data that the paper extends to an AR(1) random effect.","marker":"Eliot et al. (2011)"},{"why":"Justifies using EM for penalised likelihood estimation.","marker":"Green (1990)"},{"why":"Provides the GCV criterion for choosing the ridge parameter.","marker":"Golub et al. (1979)"},{"why":"Provides the modified GCV criterion for heteroskedastic errors, needed because the working model has non-constant error variance.","marker":"Andrews (1991)"},{"why":"Supplies the supervised component-based regularisation approach adapted for the high-dimensional case.","marker":"Bry et al. (2013)"},{"why":"Extends the component-based method to GLMMs and serves as the baseline for the component-based variant.","marker":"Chauvet et al. (2016)"}],"fun_headline_variants":["Regularised EM for GLMMs with AR(1) random effect","Ridge-penalised EM estimates AR(1) effect in GLMMs","Component-based regularised EM for high-dimensional GLMMs","L2-penalised EM for panel GLMMs with AR(1) effect","Autoregressive GLMMs fitted by regularised EM"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The whole algorithm rests on the first-order linearisation that treats $g(y_i)$ as approximately normal with working variable $z_i = g(\\mu_i) + (y_i - \\mu_i) g'(\\mu_i)$; if that approximation is bad, for example with small counts or high dispersion, the EM estimates can be biased.","fun_headline_variants_meta":{"raw":{"variants":["Regularised EM for GLMMs with AR(1) random effect","Ridge-penalised EM estimates AR(1) effect in GLMMs","Component-based regularised EM for high-dimensional GLMMs","L2-penalised EM for panel GLMMs with AR(1) effect","Autoregressive GLMMs fitted by regularised EM"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000213,"raw_usage":{"total_tokens":1368,"prompt_tokens":840,"completion_tokens":528,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":456,"completion_tokens_details":{"reasoning_tokens":433}},"tokens_in":456,"tokens_out":528,"duration_ms":5118,"temperature":1.0,"reasoning_tokens":433,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T14:15:58.388716+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the algorithm on Poisson GLMM data with an AR(1) random effect and small expected counts (mean near 0.5), then check whether the mean squared errors of $\\beta$, $\\sigma_2^2$, and $\\rho$ still decrease to zero as $T$ grows; if they level off or increase, the linearisation step is not adequate.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the linearisation step that converts the GLMM into a working linear mixed model."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Justifies using EM for penalised likelihood estimation."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the modified GCV criterion for heteroskedastic errors, needed because the working model has non-constant error variance."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the supervised component-based regularisation approach adapted for the high-dimensional case."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Extends the component-based method to GLMMs and serves as the baseline for the component-based variant."}],"review_version":1}