REVIEW 4 major objections 5 minor 8 references
Regularising Generalised Linear Mixed Models with an autoregressive random effect
T0 review · 4 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read 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)…
desk verdict 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. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
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.
What would settle it
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.
Extended reading notes
Core claim
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.
Load-bearing premise
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.
Editorial extensions
If this is right
- 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.
Reading between the lines
- 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.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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.
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 (4)
- [Section 3.1, Algorithm 1 (step 2.c)] 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 3.1, Algorithm 1 (step 2.b)] 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 4, Figure 1] 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 3, Linearisation step] 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.
minor comments (5)
- [Section 3.1, GCV formula] 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 2, design matrices] 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 4, Figure 1(d)] 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 3.2] 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.
- [General] The paper ends abruptly after the numerical section; a short conclusion discussing limitations and next steps would be helpful.
Circularity Check
No significant circularity: linearisation, penalised EM, GCV and simulation evaluation are self-contained; the under-specified M-step for ρ is an incompleteness, not a circular reduction.
full rationale
The derivation chain is: order-1 linearisation of the link gives a working variable z and a linearised model M; an L2-penalised complete log-likelihood Lpen is defined; an EM Q-function is formed as Eξ|z[Lpen(θ;z[t],ξ)|θ[t],λ[t]] and the iterate is θ[t+1]=argmaxθ Qpen; the ridge parameter is chosen by GCV. None of these steps uses the paper's conclusions as an input. The simulation in Figure 1 is an evaluation of the algorithm on data generated from the same model, which is standard practice and not a prediction derived from the fitted parameters. The self-citations (Bry and Verron 2015 for the SR criterion, Chauvet et al. 2016 for SCGLR, Bry et al. 2013) introduce design choices and an extension baseline; they are not invoked as an external theorem that forces the paper's claimed estimates. The main weakness is that Algorithm 1's M-step for σ2_2 and ρ is only asserted via θ[t+1]=argmax Qpen, and the hat-matrix S[t]_λ is not defined; this is an under-specification/reproducibility risk, not a case where a fitted input is renamed a prediction or where an equation reduces to its own input. Likewise, the high-dimensional component-based method defines a new objective Lreg as a weighted combination of L and a structural-relevance criterion; this is a proposed construction, not a hidden restatement of the data. Therefore no circular step can be exhibited from the text.
Assumptions & free parameters
free parameters (3)
- Ridge penalty lambda =
Chosen by GCV at each iteration
- SR weight s =
Tuned by cross-validation
- SR power l =
Tuned by cross-validation
assumptions (4)
- domain assumption Conditional distribution of Y|xi belongs to the exponential family
- domain assumption Random effects xi1 ~ N(0, sigma1^2 I_N) and xi2 follows a stationary AR(1) process, independent of xi1
- domain assumption Balanced panel data with N individuals observed at T time-points
- domain assumption First-order linearisation of the link function provides an approximately normal working variable z
Cite this review
Pith. "Pith review of Regularising Generalised Linear Mixed Models with an autoregressive random effect." pith.science (2026). https://pith.science/paper/MYZHWL42
@misc{pith2026190807477,
author = {Pith},
title = {Pith review of: Regularising Generalised Linear Mixed Models with an autoregressive random effect},
year = {2026},
howpublished = {\url{https://pith.science/paper/MYZHWL42}},
note = {Machine review of arXiv:1908.07477}
}
read the original abstract
We address regularised versions of the Expectation-Maximisation (EM) algorithm for Generalised Linear Mixed Models (GLMM) in the context of panel data (measured on several individuals at different time-points). A random response y is modelled by a GLMM, using a set X of explanatory variables and two random effects. The first one introduces the dependence within individuals on which data is repeatedly collected while the second one embodies the serially correlated time-specific effect shared by all the individuals. Variables in X are assumed many and redundant, so that regression demands regularisation. In this context, we first propose a L2-penalised EM algorithm, and then a supervised component-based regularised EM algorithm as an alternative.
Figures
Reference graph
Works this paper leans on
-
[1]
Andrews, D.W. (1991). Asymptotic optimality of generalized CL, cross- validation, and generalized cross-validation in regression with het- eroskedastic errors. Journal of Econometrics , 47, 359 –
work page 1991
-
[5]
Graphical diagnoses relative to the L2-penalised EM algorithm. (a),(b): 40 trajectories of the L2-convergence criterion for parameters β and σ2 1 (A similar behaviour is observed for parameters σ2 2 and ρ). About a hundred iterations is necessary to achieve convergence. (c): MSEs of parameters β,σ 2 1,σ 2 2 and ρ on simulated data where N = 10 and T ∈ { 1...
work page 1979
-
[60]
Bry, X. and Verron, T. (2015). THEME: THEmatic model exploration through multiple co-structure maximization. Journal of Chemomet- rics, 29, 637 –
work page 2015
-
[223]
Green, P.J. (1990). On use of the EM for penalized likelihood estimation. Journal of the Royal Statistical Society, Series B , 52, 443 –
work page 1990
-
[377]
Bry, X., Trottier, C., Verron, T., and Mortier, F. (2013). Supervised com- ponent generalized linear regression using a pls-extension of the fisher scoring algorithm. Journal of Multivariate Analysis , 119, 47 –
work page 2013
-
[452]
Schall, R. (1991). Estimation in generalized linear models with random effects. Biometrika, 78, 719 –
work page 1991
-
[647]
Chauvet, J., Trottier, C., Bry, X., and Mortier, F. (2016). Extension to mixed models of the Supervised Component-based Generalised Linear Regression. COMPSTAT: Proceedings in Computational Statistics . Eliot, M., Ferguson, J., Reilly, M.P., and Foulkes, A.S. (2011). Ridge Re- gression for Longitudinal Biomarker Data. The International Journal of Biostati...
work page 2016
-
[727]
Schelldorfer, J., Meier, L., and B¨ uhlmann, P. (2014). Glmmlasso: an algo- rithm for high-dimensional generalized linear mixed models using l1- penalization. Journal of Computational and Graphical Statistics , 23, 460 – 477
work page 2014
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.