{"id":"645c085b-418c-40c8-8178-89411d8240e5","arxiv_id":"2607.10432","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"EM-INLA fits Bayesian hierarchical quantile regression with closed-form EM steps and INLA's fast Gaussian approximations, plus a sandwich correction for calibrated intervals.","lead":"Researchers combined the EM algorithm with INLA to fit quantile regression on large hierarchical data without slow MCMC sampling. The method handles a 400,000-student Colombian national exam in minutes, estimating effects at the low and high ends of the score distribution.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"EM-INLA is not a standard EM for the stated marginal likelihood: the E-step conditions on point estimates of the random effects rather than integrating them jointly with v, and the M-step objective omits the random-effect priors that the INLA fit actually uses; no monotone-ascent target is establish","rationale":"I read the paper in good faith. The GIG moment derivations and closed-form variance updates are correct, and the validation against lqmm/QREM plus the sandwich coverage results are genuine evidence that the algorithm works in the tested settings. However, the central claim rests on the iterates solving a well-defined estimation problem. The paper's own Section 6 flags the lack of monotone ascent; my reading of Sections 2-3 and Appendix A shows the issue is more fundamental than INLA's numerical approximation. The E-step in Eq. (10) conditions on point estimates of alpha, but alpha is a latent random vector, so the standard EM monotonicity theorem does not apply; and Eq. (14) omits the random-effect prior terms that the INLA call in Algorithm 1 includes. Thus no written objective is provably increased. The stopping rule on sigma (Eq. 20) is not sufficient for full-vector convergence. The 25-replication HMC benchmark, partly run without convergence, cannot fully support the headline accuracy and speedup claims. These concerns do not disprove the method's practical value, but they mean the paper should remain CONDITIONAL until the convergence/objective gap is either resolved, or explicitly scoped with a diagnostic showing the approximation is benign in the claimed regimes.","tokens_in":32687,"tokens_out":15154,"duration_ms":159315,"concrete_test":"On a small dataset (n=300, J1=15, J2=5) generated under Scenario M1, instrument the implementation to evaluate the exact marginal observed-data log-likelihood under the ALD working model (integrating over v and alpha by numerical quadrature or a long HMC reference run) after every EM-INLA iteration and at the converged point. If the sequence is not monotonically increasing, or if a single Newton step from the converged point improves the marginal log-likelihood by more than 1%, EM-INLA is not converging to a maximizer of the stated objective; repeating this diagnostic across M1-M4 would determine whether the approximation is benign in the tested settings.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim that EM-INLA 'recovers the true parameters' and maximizes an empirical-Bayes ALD working likelihood requires the iterates to form a valid EM. Section 3.1's E-step computes moments for v_i conditional on point estimates of alpha (Eq. 10), while Section 3.4's variance update (Eq. 19) treats alpha as latent and uses its posterior moments. These two treatments are inconsistent: a standard EM would integrate over the joint posterior of (v, alpha), or treat alpha as parameters throughout. The M-step via INLA (Algorithm 1, lines 4-6) maximizes a penalized Gaussian likelihood that includes the alpha priors, but Theorem 3.1's Q-function (Eq. 14) omits those prior terms. Thus no single objective is demonstrably increased. The authors concede in Section 6 that 'strict monotone ascent ... is not mathematically guaranteed,' but the gap is stronger: the paper does not identify a well-defined objective that EM-INLA monotonically increases. The stopping rule (Eq. 20) monitors only sigma, so a flat sigma trajectory can stop away from a fixed point of the full vector (beta, alpha, sigma, sigma_k^2). The simulation evidence is encouraging, but it does not by itself establish the general accuracy claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes EM-INLA, an empirical-Bayes algorithm for hierarchical quantile regression based on the asymmetric Laplace (ALD) working likelihood. The ALD is written as a normal-exponential mixture, and the algorithm iterates between an E-step that computes closed-form Generalized Inverse Gaussian moments for the latent variables and an M-step that fits a weighted Gaussian regression with INLA. Closed-form updates are derived for the ALD scale and for the random-effect variances. The authors also propose an analytic cluster-robust sandwich correction for the fixed-effect intervals and compare it with a cluster bootstrap. The method is evaluated in four simulation scenarios against lqmm/qrLMM, QREM, and Hamiltonian Monte Carlo, and is applied to the 2023 Prueba Saber 11 data with 415,472 students nested in schools and municipalities. The central claims are: EM-INLA avoids MCMC, recovers parameters comparably to HMC, achieves 15x–50x speedups, produces calibrated intervals via the sandwich, and scales sublinearly with sample size.","tokens_in":33027,"tokens_out":4992,"duration_ms":60310,"significance":"If the central claims hold, this is a practically valuable contribution: it would bring empirical-Bayes hierarchical quantile regression to data sizes at which MCMC is infeasible, while retaining closed-form updates and a single-pass uncertainty correction. The manuscript has notable strengths: the derivations in Appendix A are explicit and the σ update reproduces the known ALD maximum-likelihood estimator in the baseline case; the simulations compare against independent external implementations (lqmm, qrLMM, QREM, HMC/Stan); the reproducibility statement indicates that code is available; and the interval-coverage problem is addressed directly, with the sandwich estimator checked against a cluster bootstrap. However, the load-bearing point that EM-INLA is a valid EM algorithm with a well-defined monotonically increasing objective is not established, and the authors themselves acknowledge that strict monotone ascent is not guaranteed. This gap must be resolved or the algorithm must be reframed as an approximate fixed-point procedure with appropriate convergence diagnostics.","major_comments":[{"comment":"The algorithm is presented as an EM algorithm, but no single objective is shown to be monotonically increased. The E-step in Eq. (10) conditions the latent variables v_i on point estimates of the location parameters, while the M-step for random-effect variances in Eq. (18)–(19) treats alpha as latent and uses posterior moments. Meanwhile, Theorem 3.1’s Q-function Eq. (14) omits the random-effect prior terms that the INLA fit in Algorithm 1, lines 4–6, actually maximizes. Thus there is an inconsistency between the Q-function used in the derivation and the criterion actually optimized by INLA. The Section 6 concession that 'strict monotone ascent ... is not mathematically guaranteed' is weaker than the actual problem: the paper does not identify a well-defined objective that EM-INLA increases. This is load-bearing for the 'EM' description and for the claim that the algorithm maximizes an e","section":"§3.1–3.4 and Algorithm 1, especially Eq. (10), Theorem 3.1 Eq. (14), and Prop. 3.2 Eq. (19)"},{"comment":"The stopping rule monitors only the relative change in σ. The paper states that once σ stabilizes, 'the location and variance updates reach a steady state,' but no argument is given for why flatness of σ implies convergence of β, α, and the random-effect variances. This matters because the objective-function gap above means that a flat σ trajectory could stop at a point that is not a fixed point of the full update vector. At minimum, the simulation section should report the full-vector changes over iterations, or the stopping rule should include the location parameters and variance components. Without this, the empirical convergence claims are incomplete.","section":"§3.5, Eq. (20)"},{"comment":"The analytic sandwich estimator is derived from a first-order expansion of the penalized criterion, and its applicability to two-level nested random effects is asserted rather than proved. The authors note that the derivation assumes Gaussian random effects, exogeneity, and a constant residual density, and they relegate a formal asymptotic treatment to future work. The simulation coverage in Table 5 is encouraging, but it is based on one design with J2 = 10 top-level clusters and does not exercise the assumptions (e.g., non-Gaussian random effects, clustered covariates, or a larger number of top-level units). Since the application and the paper’s uncertainty-quantification claims rely on Eq. (21), the manuscript should either provide a more formal justification or clearly state the conditions under which the intervals are expected to be calibrated, with targeted simulations for those con","section":"§3.5, Eq. (21), and §4.5"}],"minor_comments":[{"comment":"The introduction describes the method as an 'exact EM algorithm,' but Section 6 states that the M-step is approximate because it uses an INLA fit. Please replace 'exact' with 'closed-form EM-type' or similar throughout, so the terminology is consistent.","section":"Section 1"},{"comment":"The conditioning notation in Eq. (10) is inconsistent with later expressions: Eq. (10) conditions on current parameter estimates, while Prop. 3.2 conditions on '˜y, β̂'. Please unify the notation and clarify that the posterior moments in Eq. (19) are INLA approximations.","section":"Section 3.1, Eq. (10)"},{"comment":"Step 6 uses Var(α_j^(k) | ·) from INLA, but the variance update in Eq. (19) is derived as an EM conditional expectation. Since the INLA fit conditions on the current hyperparameters and uses a Laplace approximation, the variance used is not the exact EM conditional variance. This approximation should be stated explicitly next to Proposition 3.2.","section":"Algorithm 1, step 6"},{"comment":"The HMC results under M2 appear to be dominated by a small number of divergent fits (RMSE of σ2 between 16.8 and 22.3 against a true value of 4). Reporting the number of divergent transitions or other convergence diagnostics for the HMC runs would make the comparison more interpretable and would strengthen the claim that EM-INLA is 'comparable to HMC.'","section":"§4.4, Table 3"},{"comment":"The QREM variance collapse under M4 is mentioned qualitatively. Please report the frequency of boundary estimates (zero variance) across replications, so the reader can assess how typical the collapse is.","section":"§4.3, Table 2"},{"comment":"The code availability statement says replication code is available at GitHub but no URL is given. Please provide a link or repository identifier for journal review.","section":"Reproducibility"},{"comment":"There are minor typographical issues, e.g., 'Wepropose' in the abstract and 'understating' where 'understating uncertainty' is intended. A careful copyedit is recommended.","section":"Throughout"},{"comment":"Please verify the status and URLs of the references dated 2026, since some appear to be online-first or arXiv preprints; give the most citable version.","section":"References"}],"recommendation":"major_revision","confidential_remarks":"The stress-test concern about the missing monotone-ascent objective is, on close reading, valid: the paper does not define a single objective that EM-INLA increases, and the authors' own Section 6 statement is a concession of this point. However, I do not see this as an irreparable flaw; the simulations are extensive and the method is likely useful as an approximate algorithm. The authors should either supply a formal convergence argument for a well-defined objective or reposition EM-INLA as an approximate fixed-point method and add diagnostics on all parameters. The sandwich derivation also deserves more rigor or a narrower claim. If these issues are addressed, the paper could be suitable for publication."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a genuinely useful computational paper, not a theoretical one. The EM-INLA combination is new as far as I know, the closed-form GIG moments and variance updates are neat, the sandwich correction restores near-nominal coverage, and the scale-up to 415k students with 15k random effects is real. I'd send it to a methods journal without hesitation.\n\nWhat it does well: Appendix A derivations check out; the sigma update reduces to the known ALD MLE in the symmetric no-random-effects case, which is a good sanity check. The simulation design covers four error structures, compares against three credible alternatives (lqmm, qrLMM, QREM), and uses HMC as a Bayesian benchmark. The cluster bootstrap as an independent check on the sandwich is the right way to validate an uncertainty estimator. The Saber 11 application is a nice real-world demonstration, and the runtime scaling plots are informative.\n\nSoft spots: the main one is also the most interesting. This is not a standard EM. The E-step in Section 3.1 conditions on point estimates of the random effects, while the variance update in Section 3.4 treats them as latent with posterior moments. The M-step replaces the exact maximizer with an approximate INLA fit. So there is no single objective that the iterates provably increase; the stress-test note is correct that the gap is stronger than just \"monotonicity is not guaranteed\" — it's not clear what objective the algorithm is climbing. The authors are honest in Section 6 about the lack of a formal guarantee, and the simulations show the method works in practice, but the theory is absent. That's acceptable for a computational paper if the evidence is strong, and it mostly is.\n\nThe other soft spots are minor. The sandwich correction is a heuristic, with asymptotics deferred; the coverage results are good, but I'd want a more careful statement of when it can fail. The HMC benchmark under M2 is run with a fixed small budget — HMC does not converge and the reported speedups are partly against a non-converged algorithm. The paper is transparent about this, but the headline \"15x to 50x faster than HMC\" is a bit inflated. Finally, the code is promised on GitHub but no URL appears, so the reproducibility claim can't be checked today.\n\nWho is this for? Applied statisticians and method developers working on large-scale quantile regression. The paper deserves a serious referee: the novelty is real, the empirical work is solid, and the limitations are acknowledged. I'd accept it subject to a convergence discussion that is honest about the lack of a formal objective and a code link that actually works.","headline":"A practical EM+INLA hybrid for hierarchical quantile regression that scales to 400k observations; the EM theory is loose, but the closed-form updates and empirical validation make it worth refereeing.","tokens_in":33501,"tokens_out":2105,"would_cite":true,"duration_ms":27042,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62J05","62G08"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper proposes EM-INLA, an algorithm that makes Bayesian hierarchical quantile regression tractable for very large datasets by replacing MCMC with closed-form EM steps and INLA approximations.","keywords":["Bayesian quantile regression","Asymmetric Laplace Distribution","EM algorithm","Integrated Nested Laplace Approximations","hierarchical models","empirical Bayes","cluster-robust sandwich","large-scale inference"],"falsifier":"On a moderate simulated dataset where the exact M-step is computable, evaluate the penalized ALD log-likelihood after every EM-INLA iteration; any non-negligible decrease from one iteration to the next would falsify the implicit monotone-ascent premise of the algorithm.","tokens_in":61,"feed_emoji":"📊","tokens_out":4473,"duration_ms":108687,"temperature":0.7,"pith_summary":"This paper aims to show that Bayesian quantile regression for hierarchical data, traditionally a computationally heavy MCMC problem, can be solved by an EM algorithm whose M-step is a weighted Gaussian regression handled by INLA, with scale and variance updates in closed form. If correct, the method delivers point estimates and calibrated uncertainty for quantile-specific covariate effects in settings with hundreds of thousands of observations and thousands of random effects, where Hamiltonian Monte Carlo is impractical. The authors demonstrate accuracy comparable to HMC in simulations, speedups of 15x to over 50x, and a successful fit to 415,472 students nested in 13,801 schools and 1,116 municipalities, revealing tail-specific patterns in socioeconomic effects.","feed_headline":"Quantile regression without MCMC: 50x faster, 400k students","feed_subtitle":"EM-INLA matches Hamiltonian Monte Carlo accuracy at a fraction of the cost, fitting 400k students in minutes.","key_machinery":"The engine is the normal-exponential mixture representation of the Asymmetric Laplace Distribution: Y = mu + theta V + kappa sqrt(sigma V) Z, with V exponential and Z standard normal, which makes Y conditionally Gaussian given V. The E-step then reduces to closed-form moments of a Generalized Inverse Gaussian distribution with index 1/2, and the M-step becomes a weighted Gaussian regression on pseudo-responses, solved by INLA; the scale and random-effect variances update analytically. The cluster-robust sandwich formula is the companion device that converts the working-likelihood posterior spread into calibrated frequentist intervals.","core_discovery":"The central claim is that the Asymmetric Laplace working likelihood, via its normal-exponential mixture, turns each EM M-step into a weighted Gaussian regression that INLA can fit, so no MCMC and no smoothing of the check loss are needed. The E-step conditional distributions of latent variables are Generalized Inverse Gaussian with closed-form moments, the scale parameter has an exact closed-form update, and each random-effect variance is updated by posterior mean plus variance from INLA. A final INLA run conditional on the converged hyperparameters gives posterior marginals, and a one-pass cluster-robust sandwich correction restores near-nominal interval coverage that the naive empirical-Ba","pith_inferences":["Because the method inherits INLA's prior flexibility, the same EM-INLA loop should extend to spatial, temporal, or spatio-temporal random effects, giving scalable quantile versions of areal and longitudinal models.","The absence of a formal convergence guarantee for the approximate M-step suggests a targeted diagnostic: monitoring the penalized ALD objective directly, not just sigma, would catch drift when INLA's approximation is crude.","The empirical-Bayes conditioning on hyperparameters means practitioners should treat the sandwich or bootstrap intervals as the inferential output, not the raw posterior marginals; this is a design choice, not a bug.","Fitting each quantile independently invites quantile crossing; combining EM-INLA with non-crossing or rearrangement methods would likely preserve its speed while restoring coherent quantile functions."],"forward_implications":["EM-INLA removes the need for MCMC in hierarchical quantile regression, making quantile analyses of administrative-scale datasets feasible in minutes.","Simulations indicate point estimates match Hamiltonian Monte Carlo across normal, heavy-tailed, and heteroscedastic error scenarios, with 15x to 50x speedups.","The closed-form random-effect variance update stays strictly positive, avoiding the zero-variance boundary collapse seen in frequentist mixed-model quantile fits.","Naive intervals from the final INLA call undercover; the analytic cluster-robust sandwich correction brings coverage to near nominal in one extra pass, matching a cluster bootstrap benchmark.","On the 2023 national test data, quantile effects reveal that parental education and internet access matter most in the upper tail, while the private-school premium is largest at the bottom."],"fun_headline_variants":["Quantile regression without MCMC: EM-INLA scales to 400k students","EM-INLA: MCMC-free Bayesian quantile regression, 50x speedup","No MCMC, 50x faster: EM-INLA for large-scale quantile regression","EM-INLA: Bayesian quantile regression for big data, no MCMC","From 15x to 50x faster: EM-INLA replaces MCMC in quantile regression"],"cache_read_input_tokens":34816,"weakest_assumption_plain":"The M-step replaces the exact weighted-Gaussian maximizer with INLA's approximate fit, and the paper does not prove that the objective still ascends; if that approximation error accumulates across iterations, the algorithm could converge to a point far from the true marginal maximum.","fun_headline_variants_meta":{"raw":{"variants":["Quantile regression without MCMC: EM-INLA scales to 400k students","EM-INLA: MCMC-free Bayesian quantile regression, 50x speedup","No MCMC, 50x faster: EM-INLA for large-scale quantile regression","EM-INLA: Bayesian quantile regression for big data, no MCMC","From 15x to 50x faster: EM-INLA replaces MCMC in quantile regression"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000587,"raw_usage":{"total_tokens":2607,"prompt_tokens":770,"completion_tokens":1837,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":514,"completion_tokens_details":{"reasoning_tokens":1725}},"tokens_in":514,"tokens_out":1837,"duration_ms":12672,"temperature":1.0,"reasoning_tokens":1725,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-02T07:13:39.704102+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a moderate simulated dataset where the exact M-step is computable, evaluate the penalized ALD log-likelihood after every EM-INLA iteration; any non-negligible decrease from one iteration to the next would falsify the implicit monotone-ascent premise of the algorithm.","supporting_citations":[],"review_version":2}