{"id":"7b6618ed-3fa4-4208-a739-d587fb03f30b","arxiv_id":"2411.14185","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A new cAIC formula that assumes multivariate normal random effects performs robustly across continuous and discrete observation distributions in nonlinear mixed-effects model simulation.","lead":"Two formulas for computing the conditional Akaike information criterion in nonlinear mixed-effects models are proposed, and simulations show that one of them works for Gaussian, gamma, negative binomial, and Tweedie data. The paper offers practitioners a numerical recipe for model selection in mixed-effects settings where cluster sizes are small.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The claim that the neglected O(1/nt) and second-order terms in Eq. (6) are negligible rests on unpublished simulations; for the paper's own small-cluster regime this is unverified and load-bearing.","rationale":"The reader's weakest-assumption analysis identifies exactly the same load-bearing point: the derivation of Eq. (6) drops an O(1/nt) term and unspecified second-order Taylor terms, with support coming from an unpublished numerical study and an informal order argument. My independent reading of Appendix C confirms that the omitted curvature difference is not bounded in the paper, and the O(1/nt) expression in Eq. (24) is acknowledged to be hard to evaluate, including at boundary parameter values. The paper does provide real evidence in its favor: a complete Gaussian simulation code, reproducible-looking tables, and the observation that Eq. (6) reduces to the known Liang et al. (2008) formula in the linear mixed-model case. Those strengths mean the concern does not refute the method; it means the central claim is not yet fully established. A conditional acceptance, with the concrete test above as a condition, matches the evidence. I therefore keep the reader's CONDITIONAL verdict unchanged rather than moving to reject or accept.","tokens_in":18160,"tokens_out":9059,"duration_ms":92296,"concrete_test":"In the Sec. 3.1 Gaussian setting (T=50, nt=3), compute the full O(1/nt) term from Eq. (24) using the true Sigma_u, the fitted COV(theta_hat) from sdreport, and partial Psi_hat / partial theta from the joint Hessian; additionally estimate the dropped second-order curvature difference by Monte Carlo over y*. Recompute the relative bias and, more decisively, the cAIC rankings across the four delta models. If including these terms changes any model ranking relative to the true cAI, or shifts the bias correction by more than 5% of BCtrue, the paper's assertion that Eq. (6) is a distributionally robust cAIC criterion is not supported as stated.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central recommendation of Method 2 (Eq. (6)) depends on two unproved cancellations. First, the O(1/nt) term in Eq. (6), which the paper itself states is -2 trace{Sigma_u^{-1} (partial Psi_hat / partial theta^T) COV(theta_hat) (partial Psi_hat^T / partial theta)}, is not shown to be small for the cluster sizes (3 or 5) emphasized throughout the paper. The paper argues the term is negligible because partial Psi_hat / partial theta is 'typically small', but no bound or simulation evidence is provided; the conclusion explicitly relies on 'unpublished simulation results'. Second, Appendix C drops the 'difference between the second order terms' between Eq. (21) and Eq. (22), justified only by an 'unpublished numerical study'. If either neglected contribution is material, the bias correction in Eq. (6) is itself biased, and the claimed distributional robustness of cAIC model rankings is unsupported. This is not an internal inconsistency, but it is the weakest load-bearing assumption in the argument for the headline claim.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes two estimators for the conditional Akaike information criterion (cAIC) in nonlinear mixed-effects models. Method 1 (Eq. 4) is a Taylor-expansion-based bias correction for continuous data and gives explicit derivative formulas for fixed-effects estimators and empirical Bayes random-effect predictors with respect to the data (Lemma 1). Method 2 (Eq. 6) assumes the random effects are modeled as multivariate normal and uses a Hessian-based trace correction. The paper reports simulation studies for Gaussian, gamma, negative binomial, and Tweedie observations with small fixed cluster sizes (nt = 3 or 5), evaluating the relative bias of the bias-correction term. Based on these simulations, the authors recommend Method 2 as a distributionally robust cAIC criterion for model selection.","tokens_in":18405,"tokens_out":9040,"duration_ms":84822,"significance":"If the central claim holds, the paper fills a practical gap by offering a cAIC formula applicable to nonlinear mixed-effects models with non-Gaussian, skewed, discrete, or zero-inflated responses, going beyond existing cAIC work that often assumes Gaussian or known-dispersion settings. The paper provides analytical derivative expressions, leverages TMB for automatic differentiation, and ships R/TMB code for the Gaussian example, which are concrete strengths. However, the central formulas rest on two approximations that are justified only by references to unpublished numerical studies, and the paper's own discussion acknowledges that the O(1/nt) term is not negligible in the small-cluster regime it emphasizes. These issues make the theoretical claim less secure than the simulations alone suggest.","major_comments":[{"comment":"The O(1/nt) term in Eq. (6) is dropped on the basis of informal reasoning, but the paper's stated regime is small, bounded cluster sizes (nt = 3 or 5). The paper itself states that this term 'approaches zero only as cluster sizes increase' (p.5), and the subsequent two arguments—that differences between models are dominated by O(T) terms and that ∂Ψhat/∂θ^T is 'typically small'—do not constitute a bound. The only direct support is the assertion of 'unpublished simulation results' in the Conclusion (p.9). Since a non-negligible O(1/nt) term would make the bias correction in Eq. (6) itself biased, this is load-bearing for the headline recommendation of Method 2. The authors should either provide a rigorous condition under which the term is negligible, or present a reproducible numerical study that quantifies its magnitude across a grid of nt values and model configurations, or show that the simulation conclusions are insensitive to its inclusion.","section":"§2.2, Eq. (6), and p.5"},{"comment":"Both derivations neglect the difference between second-order Taylor terms. In Appendix B the text says, 'If we neglect the difference between the second order terms of (12) and (16) (our unpublished numerical study also indicates that the difference is small),' and in Appendix C the same step is introduced with 'Here we neglect the difference between the second order terms.' These are unverifiable approximations at the core of both Eq. (4) and Eq. (6). Because the paper is a methodological contribution, the justification cannot be an unpublished numerical study. Please include the numerical study in the supplementary material with sufficient detail to reproduce it, or provide a formal argument that the neglected difference is of lower order than the retained terms.","section":"Appendix B (p.11) and Appendix C (p.13)"}],"minor_comments":[{"comment":"The relative-bias estimates are reported without Monte Carlo standard errors or confidence intervals. Given that many values are close to zero (e.g., -0.002 in Table 1), it is unclear how much of the deviation is Monte Carlo noise; adding standard errors would substantially improve interpretability.","section":"§3, Tables 1-4"},{"comment":"The simulations for Eq. (4) and Eq. (6) use different numbers of outer and inner iterations (e.g., nout=1500/ninner=20000 for Eq. (6) versus nout=500/ninner=1000 for Eq. (4) in the Gaussian example). If the 'true BC' is recomputed with these different Monte Carlo sizes for each method, the comparison of the two methods mixes formula error with Monte Carlo error. Please state whether the same true-BC estimate was used for both methods and, if not, report the Monte Carlo uncertainty.","section":"§3.1 and §3.2"},{"comment":"There is a grammatical error: 'According to the the conditional inference' should read 'According to the conditional inference.'","section":"Appendix C, p.12"},{"comment":"The discussion of the O(1/nt) term first gives an expression involving the true random-effect covariance Σ_u and then notes that the MVN assumption for p(u) is unnecessary and that Σ_u^{-1} can be replaced by the model precision matrix; stating the final usable expression once in equation form would avoid ambiguity about which covariance is intended.","section":"§2.2, p.5"}],"recommendation":"major_revision","confidential_remarks":"The paper's core methodological claim is plausible and the simulations are encouraging, but the reliance on 'unpublished simulation results' at two load-bearing points is unusual for a methods paper and should be addressed in revision. If the authors can supply the unpublished numerical studies or a rigorous order argument, the contribution would be publishable. The scope matches stat.ME; the practical relevance for fisheries and ecology applications is a strength."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Nan, quick take on arXiv:2411.14185. This is a genuinely useful note. The field has cAIC for linear and generalized linear mixed models, but nonlinear mixed-effects models with non-Gaussian, discrete, or zero-inflated observations lacked a workable criterion. The paper supplies two: Eq. (4) via Taylor expansion with derivatives of estimators w.r.t. data (Lemma 1), and Eq. (6), a simple trace correction that only needs the joint and RE-prior Hessians. Method 2 is the real contribution. It extends Yu–Yau and Liang et al. to NMMs with arbitrary links and small clusters, and the simulation evidence is consistent: relative biases in the bias-correction are small across Gaussian, gamma, negative binomial, and Tweedie cases, with cluster sizes 3 and 5. That is not a trivial result.\n\nWhat I like: the paper is honest about what is approximate. It explicitly flags the O(1/nt) term in Eq. (6) and says why it might be small. It also ships working R/TMB code for the Gaussian example, which is a real artifact. The simulations are structured sensibly, and the 'true' BC is computed by simulation from the generating model, which is standard in this literature.\n\nNow the soft spots, proportionate to how soft they are. The load-bearing concern is real: the drop of the O(1/nt) term and of the second-order differences in Appendix C is justified by an 'unpublished numerical study' and an informal order argument. For the small-cluster regime the paper itself emphasizes, that's not a proof. The paper would be stronger with either a bound on the term or the unpublished simulations included as a supplement. That said, the published simulations do empirically support the claim across many configurations, so the concern is not that the method fails; it's that the main theoretical justification for the drop is not in the paper. Second, the RB tables have no Monte Carlo standard errors, so we don't know the noise in those numbers. Minor. Third, code is only for the Gaussian case; the non-Gaussian and Tweedie simulations are not reproducible from the supplement. Also minor, but worth asking for.\n\nThe citation pattern looks fine. The Zheng–Cadigan self-citations are to published papers and are used for intermediate steps; that's appropriate.\n\nBottom line: this is a useful note, worth a serious referee. The recommendation is to accept after revision, with the authors asked to supply the unpublished simulation results or a theoretical argument for the dropped terms, plus code for at least one non-Gaussian case.","headline":"Practical cAIC estimator for nonlinear mixed models that works in simulations; the dropped O(1/nt) term deserves a real justification before the paper is accepted.","tokens_in":18903,"tokens_out":2656,"would_cite":true,"duration_ms":24683,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper derives two conditional-AIC estimators for nonlinear mixed-effects models and recommends Eq. (6) as the distributionally reliable one for any observation type.","keywords":["conditional Akaike information","nonlinear mixed-effects models","model selection","bias correction","empirical Bayes","Tweedie distribution","random effects","automatic differentiation"],"falsifier":"Use the paper's own Monte Carlo protocol (outer draws of true random effects, inner draws of prediction data) with tiny cluster sizes such as $n_t=2$, a skewed observation distribution, and a random-effect variance on the boundary of its parameter space. If the relative bias of Eq. (6)'s bias correction exceeds the paper's typical benchmark of roughly 0.035, or if adding the explicit $O(1/n_t)$ term from the Supplementary Material changes which model wins, the central recommendation fails.","tokens_in":17977,"feed_emoji":"📊","tokens_out":12452,"duration_ms":104062,"temperature":0.7,"pith_summary":"The paper tries to give practitioners a usable way to compute the conditional Akaike information criterion for nonlinear mixed-effects models, where cluster sizes are often small and observations may be skewed, discrete, or zero-inflated. It proposes two formulas: Method 1, a derivative-based correction for continuous data, and Method 2, a correction that assumes only that the modeled random-effect distribution is multivariate normal and is claimed to work for any observation type. The paper's central recommendation is Method 2, Eq. (6), whose penalty self-adjusts between counting only fixed parameters and counting fixed plus random effects. If the claim is right, one formula covers Gaussian, gamma, negative binomial, and Tweedie data with flexible links and small clusters, removing a major obstacle to model selection in nonlinear mixed models.","feed_headline":"One formula fixes conditional AIC for nonlinear mixed models","feed_subtitle":"Recommended criterion holds bias low for normal, gamma, negative binomial, and Tweedie data at small cluster sizes.","key_machinery":"The load-bearing identity is Theorem 2's Eq. (6), which expresses the optimism of the fitted conditional log-likelihood as $2p_c + 2q - 2\\operatorname{trace}\\{\\ddot{\\ell}_j^{-1}\\ddot{\\ell}_r\\}$. Here $\\ddot{\\ell}_j$ is the second derivative of the joint log-likelihood with respect to the random effects and $\\ddot{\\ell}_r$ is the second derivative of the random-effect log-density, so the trace term is the effective number of random-effect parameters. When the data are uninformative about the random effects this trace is close to $q$ and the penalty shrinks toward $2p_c$; when the data pin the random effects down it is close to $0$ and the penalty approaches $2(p_c+q)$. The derivation reaches this trace by a second-order Taylor expansion of the conditional log-likelihood around the true parameters followed by an integration-by-parts identity for the multivariate normal random-effect model.","core_discovery":"Working from the conditional Akaike information, the paper derives two estimators. Method 1 (Eq. 4) uses derivatives of the maximum marginal likelihood estimator $\\hat{\\theta}$ and the empirical Bayes predictor $\\hat{\\Psi}$ with respect to the data; the paper shows it is near-exact for Gaussian data with the identity link but degrades for skewed continuous data and cannot handle discrete or zero-inflated observations. Method 2 (Eq. 6) assumes the modeled random-effect distribution is multivariate normal and obtains the bias correction $2p_c + 2q - 2\\operatorname{trace}\\{\\ddot{\\ell}_j^{-1}\\ddot{\\ell}_r\\} + O(1/n_t)$. In the paper's simulations this formula keeps relative bias small across normal, gamma, negative binomial, and Tweedie observations with log links and cluster sizes of 3 or 5, which is the stated basis for recommending Method 2.","pith_inferences":["As an editorial extension, including the explicit $O(1/n_t)$ expression derived in the Supplementary Material when cluster sizes are very small could be a sharp test of Eq. (6); the authors note their unpublished simulations saw no improvement, but the published evidence for that is absent.","Because Eq. (6) assumes only that the modeled random effects are multivariate normal, the same trace identity should transfer to latent Gaussian process and spline models, where that assumption is standard.","As an editorial extension, the derivative identities behind Method 1 are reusable for sensitivity analysis and delta-method standard errors, not only for cAIC."],"forward_implications":["One formula, Eq. (6), applies to Gaussian, gamma, negative binomial, and Tweedie observations with flexible links and small bounded cluster sizes, so practitioners no longer need a bespoke cAIC derivation for each setting.","The criterion's effective penalty moves between $2p_c$ and $2(p_c+q)$ depending on how much information the data carry about the random effects, which lets it mimic both marginal and conditional AIC in the right limits.","The needed derivatives can be obtained numerically through automatic differentiation, so the formula is automatable for the kinds of nonlinear mixed models used in fisheries, ecology, pharmacology, and public health.","Method 1 remains a reasonable cAIC only for Gaussian data with identity links; the paper's simulations show it should not be trusted for skewed, discrete, or zero-inflated data."],"supporting_citations":[{"why":"Defines the conditional Akaike information and establishes the cAIC framework that both estimators target.","marker":"Vaida and Blanchard (2005)"},{"why":"Provides the derivative-based cAIC formula that Eq. (4) generalizes and the exact Gaussian identity-link benchmark.","marker":"Liang et al. (2008)"},{"why":"Supplies the delta-variance approximations used to justify the first-order expansions and covariance identities in the derivations.","marker":"Zheng and Cadigan (2021)"},{"why":"Gives the conditional-variance result used to derive the $O(1/n_t)$ expression for Eq. (6).","marker":"Zheng and Cadigan (2023)"},{"why":"Defines the empirical Bayes predictor that enters the conditional likelihood and the bias correction.","marker":"Kass and Steffey (1989)"},{"why":"Provides the $o(1)$ approximation-order baseline for conventional AIC against which the paper measures its $O(1/n_t)$ error.","marker":"Cavanaugh and Neath (2019)"},{"why":"Supplies the automatic differentiation engine used to compute the derivatives needed by Lemma 1 and Eq. (6).","marker":"Kristensen et al. (2016)"}],"fun_headline_variants":["Distribution-robust cAI for nonlinear mixed-effects models","Recommended cAI criterion holds bias low for non-Gaussian data","Method 2: robust cAI for skewed and discrete data","Simulation-backed cAI: robust for skewed and count data"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The claim stands on the assumption that the Taylor-series remainder and the $O(1/n_t)$ term dropped in Eq. (6) are small in practice; the paper supports this with an informal order argument in Appendix C and an unpublished numerical study, not with a proof or a published check.","fun_headline_variants_meta":{"raw":{"variants":["Distribution-robust cAI for nonlinear mixed-effects models","Recommended cAI criterion holds bias low for non-Gaussian data","Method 2: robust cAI for skewed and discrete data","Simulation-backed cAI: robust for skewed and count data"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001115,"raw_usage":{"total_tokens":4603,"prompt_tokens":863,"completion_tokens":3740,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":479,"completion_tokens_details":{"reasoning_tokens":3670}},"tokens_in":479,"tokens_out":3740,"duration_ms":27027,"temperature":1.0,"reasoning_tokens":3670,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T15:25:38.957534+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Use the paper's own Monte Carlo protocol (outer draws of true random effects, inner draws of prediction data) with tiny cluster sizes such as $n_t=2$, a skewed observation distribution, and a random-effect variance on the boundary of its parameter space. If the relative bias of Eq. (6)'s bias correction exceeds the paper's typical benchmark of roughly 0.035, or if adding the explicit $O(1/n_t)$ term from the Supplementary Material changes which model wins, the central recommendation fails.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the conditional Akaike information and establishes the cAIC framework that both estimators target."},{"cited_title":"Wu, and G","cited_arxiv_id":null,"evidence_quote":"Provides the derivative-based cAIC formula that Eq. (4) generalizes and the exact Gaussian identity-link benchmark."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the delta-variance approximations used to justify the first-order expansions and covariance identities in the derivations."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the conditional-variance result used to derive the $O(1/n_t)$ expression for Eq. (6)."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the empirical Bayes predictor that enters the conditional likelihood and the bias correction."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the $o(1)$ approximation-order baseline for conventional AIC against which the paper measures its $O(1/n_t)$ error."},{"cited_title":"Nielsen, C","cited_arxiv_id":null,"evidence_quote":"Supplies the automatic differentiation engine used to compute the derivatives needed by Lemma 1 and Eq. (6)."}],"review_version":1}