{"id":"c8fa8bb8-63ed-4142-b306-ee79f0e46a0f","arxiv_id":"2506.20060","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"The hdbayes package provides a unified Stan-based R implementation of six historical-data priors for GLMs, with consistent syntax, marginal likelihood functionality, and two worked data analyses.","lead":"This paper introduces hdbayes, an R package that implements six Bayesian priors for incorporating historical data into generalized linear models, all sharing a consistent glm-like syntax and built on Stan. A generalist reader might care because it removes a software barrier that has kept these borrowing methods from being used and compared in real clinical and observational analyses.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Multi-historical NPP in Eq. (6) uses per-dataset normalizing constants rather than the joint normalizing constant; for H>1 the omitted joint constant depends on a0, so the claimed NPP implementation may not be a normalized power prior.","rationale":"The central claim is that hdbayes implements the normalized power prior for GLMs with historical data, including multiple historical datasets (Eq. (6), Table 1). The most load-bearing condition for that claim is that the posterior for a0 is the NPP posterior. Algorithm 1 estimates Z_h(a0h) separately for each historical dataset and uses these per-dataset values in Stan. For H>1, the product of per-dataset normalizing constants is not the normalizing constant of the product of powered historical likelihoods; the missing joint integral depends on a0. This is not a numerical approximation issue that more bridge samples would fix; it is a structural mismatch between the stated model and the implemented model. The reader's weakest assumption concerned bridge-sampling and LOESS accuracy, which is a legitimate but secondary concern; the issue identified here is already present with exact per-dataset Z_h values. I would not reject the paper: the single-dataset NPP, NAPP, PP, BHM, RMAP, CP, and LEAP implementations are not implicated, and the multi-dataset NPP could be corrected or restricted to H=1 in the documentation. Thus the reader's CONDITIONAL verdict remains appropriate, with the condition now including a check of the multi-historical NPP definition.","tokens_in":20770,"tokens_out":14497,"duration_ms":165150,"concrete_test":"Take one historical dataset, split it into two exchangeable halves so the combined historical likelihood equals the original, and fit glm.npp twice: once with H=2 on the halves and once with H=1 on the combined data. If the posterior of the total discount s=a01+a02, or of beta, differs materially between the two fits, the per-dataset product normalization is not a valid NPP for H>1. An analytic companion is to verify numerically whether integral pi0(theta) L(theta)^{a01} L(theta)^{a02} dtheta equals Z_1(a01) Z_2(a02) for the split data; equality is required for Eq. (6) to match the joint NPP.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Equation (6) defines the multi-historical NPP as a product of per-dataset normalized power priors: each factor is L_h(theta)^{a0h} pi0(theta)^{1/H}/Z_h(a0h). For H>1, this product is not a normalized density in theta: the prior mass conditional on a0 is C(a0) = integral pi0(theta) prod_h L_h(theta)^{a0h} dtheta / prod_h Z_h(a0h), which is generally not 1 and varies with a0. The NPP for the combined historical data requires Z_joint(a0) = integral pi0(theta) prod_h L_h(theta)^{a0h} dtheta in the denominator, not prod_h Z_h(a0h). Algorithm 1 and glm.npp use per-dataset bridge-sampling estimates, passed to Stan after interpolation, so the a0 posterior is weighted by the omitted a0-dependent ratio Z_joint(a0)/prod_h Z_h(a0h). Concretely, with two identical historical datasets, the historical likelihood is L(theta)^{a01+a02}; the posterior of the total discount s=a01+a02 should depend only on s, but the product-normalized implementation gives different weight to splits such as (0.5,0.5) versus (0.1,0.9). Thus for H>1, glm.npp does not implement the stated NPP. For H=1, the product reduces to the joint constant, so the single-dataset NPP is not affected. The reader's bridge-sampling concern is secondary: even with exact per-dataset Z_h estimates, this structural mismatch remains.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper introduces hdbayes, an R package that implements several historical-data priors for Bayesian generalized linear models: the power prior, normalized power prior, normalized asymptotic power prior, Bayesian hierarchical model, robust meta-analytic predictive prior, commensurate prior, and latent exchangeability prior. The package is built on Stan and cmdstanr with a consistent formula/family interface, and includes functionality for marginal-likelihood-based model selection. The manuscript reviews each prior, illustrates the package on an HIV/AIDS clinical trial dataset, and shows how time-to-event data can be handled through Poisson GLM representation. The paper also proposes numerical approaches for normalizing constants and mixture-weight updating, and provides code examples and output summaries.","tokens_in":21126,"tokens_out":8067,"duration_ms":78886,"significance":"If the implementations are correct, hdbayes is a valuable contribution: it is a CRAN-released package that unifies several historical-data priors for GLMs under a single, user-friendly, HMC-based interface, and it provides the first publicly available implementation of the latent exchangeability prior. The consistent syntax and inclusion of marginal likelihood tools are genuine strengths. However, the manuscript contains a load-bearing mathematical error in the definition of the multi-historical normalized power prior (Eq. 6), and it overstates the package's GLM coverage in Section 5 and Table 1. The numerical approximation chain for normalizing constants is also not validated beyond a single example, which is insufficient for a software paper claiming to implement these priors. These issues affect the reliability of the package for multi-historical settings and require correction before the paper can be considered for publication.","major_comments":[{"comment":"For H>1, the product of per-dataset normalized power priors in Eq. (6) is not a properly normalized density over (β, φ) conditional on a0. The integral of the product is C(a0) = ∫ π0(β,φ) ∏_{h=1}^H L_h(β,φ|D0h)^{a0h} dβ dφ / ∏_{h=1}^H Z_h(a0h), which generally depends on a0 and is not equal to 1. Consequently, the posterior over a0 under Eq. (6) is weighted by the omitted factor Z_joint(a0)/∏_h Z_h(a0h), where Z_joint(a0) is the joint normalizing constant. As a concrete counterexample, for two identical historical datasets the likelihood depends only on the total discount s = a01 + a02, but the product-normalized form gives different prior weight to splits such as (0.5,0.5) versus (0.1,0.9). Thus glm.npp for H>1 does not implement the normalized power prior as stated. Algorithm 1 must either be revised to use the joint normalizing constant or the package/paper should restrict the NPP to a single historical dataset and clearly disclose this limitation.","section":"Section 2.3, Eq. (6) and Algorithm 1"},{"comment":"The paper claims in Section 5 that 'the full suite of GLMs and link functions available in the stats package is available in hdbayes' and Table 1 marks 'All models in stats::glm' for hdbayes. However, Section 2.2 states that binomial models with more than one trial are not supported by glm.pp. Since binomial models with a trial count greater than one are a standard part of stats::glm, these statements are contradictory and the 'full suite' claim is false as written. The paper should either document the actual supported family/link combinations precisely, or add support for binomial trials; at minimum, the overclaim in Section 5 should be removed.","section":"Section 2.2 vs Section 5 and Table 1"},{"comment":"The accuracy of the normalizing constant estimation chain for the NPP—bridge sampling on a grid of power-prior samples, LOESS smoothing, and linear interpolation inside Stan—is not validated beyond the single HIV example. The paper notes in Section 2.3 that for the normal linear model with a conjugate normal-gamma initial prior the normalizing constant is known, and lm.npp avoids estimation, but no comparison between glm.npp's bridge-sampling estimate and the analytic value is reported. No simulation study or independent check is provided to show that the estimated Z_h(a0h), and hence the posterior of a0 and β, is reliable. Since the NPP posterior depends critically on Z_h, the paper should include at least one validation study, such as recovery of known normalizing constants in the normal linear model or a comparison with an alternative implementation, to support the correctness of the approximation.","section":"Sections 2.3 and 4.1"}],"minor_comments":[{"comment":"The sentence 'The implementation allows for binomial models with the number of trials exceeding 1, which is a limitation of the implementation in the hdbayes implementation glm.pp' is garbled. It appears that 'allows' should be 'does not allow' or the sentence should be restructured to state clearly that glm.pp does not support binomial models with more than one trial while other implementations do.","section":"Section 2.2, paragraph on binomial models"},{"comment":"The expression for the asymptotic covariance matrix is typeset in a confusing way: 'a−1 0h h I(ˆβ0h, ˆϕ0h|D0h) i−1' should be written as a0h^{-1} [I(ˆβ0h, ˆϕ0h|D0h)]^{-1} to indicate the inverse of the Fisher information multiplied by the reciprocal of the discounting parameter.","section":"Section 2.4, Eq. (7)"},{"comment":"The text says the bridge sampling estimate uses 'samples from the prior,' but steps 2–3 sample from the power prior (4) with the given discounting parameter. Please rephrase to 'samples from the power prior posterior' or similar for clarity.","section":"Algorithm 1, line 4"},{"comment":"In Eq. (14), the denominator uses Z_V(D0), but the vague prior normalizing constant for the current data should be Z_V(D). Similarly, Algorithm 2 line 2 writes p_V(θ|D0) where the current-data posterior is intended. Also, the vague prior π_V is not explicitly defined in Section 2.6; please state that it is the same initial prior as in (5) or provide the formula.","section":"Section 2.6, Eq. (14) and Algorithm 2"},{"comment":"The sentence 'In hdbayes, we assume that the first specified data set in the list is the current data sets' contains a typo: 'data sets' should be 'data set'.","section":"Section 4.1, code description"},{"comment":"The description of the two-step marginal likelihood computation has awkward phrasing: 'after taking MCMC samples of the prior' and 'estimating the normalizing constant of the posterior' should be clarified to indicate that the MCMC samples are from the prior and posterior densities respectively, with bridge sampling applied to each.","section":"Section 3, Eq. (18) and surrounding text"}],"recommendation":"major_revision","confidential_remarks":"The multi-historical NPP issue in Eq. (6) is serious and should be verified against the package source code before acceptance; if the package indeed implements this product-normalized form, the H>1 functionality is mathematically incorrect. The overclaim in Section 5 about full GLM coverage is also an editorial problem that should be fixed. The lack of validation for the normalizing constant approximation is concerning but can be addressed with additional experiments in a revision."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: this is a solid software paper with a real package that fills a practical gap, but Eq. (6) has a load-bearing error for multiple historical datasets. The multi-historical NPP defines the prior as a product of per-dataset normalized power priors, each divided by its own Z_h(a0h). That product is not normalized: the joint normalizing constant for the combined historical likelihood, Z_joint(a0) = ∫ π0(β,ϕ) ∏_h L_h(β,ϕ)^{a0h} dβ dϕ, generally differs from ∏_h Z_h(a0h), and the ratio depends on a0. So the a0 posterior is weighted by an omitted, a0-dependent factor. For H=1 this reduces correctly; for H>1 it's a different prior. Since the package explicitly advertises support for multiple historical datasets, this needs to be addressed—at minimum acknowledged, better corrected.\n\nWhat is good: the package is on CRAN, the syntax is consistent and genuinely close to glm(), and it's the first public R implementation of LEAP. Bringing six priors under one Stan-based front-end with HMC is a useful contribution. The time-to-event Poisson-trick example is a nice touch. The paper also honestly discloses the binomial limitation in glm.pp in Section 2.2, even though the Discussion then overstates the 'full suite' claim.\n\nOther soft spots are minor in comparison: the NPP normalizing constant estimation uses a grid of PP samples, bridge sampling, LOESS smoothing, and linear interpolation inside Stan. That approximation chain is only validated on one dataset, and the paper acknowledges the NAPP approximation risk. The RMAP mixture-weight computation via bridge sampling is standard practice.\n\nOn balance, the central promise—usable, coherent implementations of these priors for GLMs—holds for the PP, BHM, RMAP, CP, LEAP, and NAPP, but not for the multi-historical NPP as written. The paper deserves peer review; a referee should push for either a corrected joint-constant implementation or a clear statement that the package only implements the NPP for H=1, with multi-historical support withdrawn until fixed.","headline":"A genuinely useful R package for historical-data priors, but the multi-historical NPP formula is structurally not a normalized power prior and should be fixed before publication.","tokens_in":21637,"tokens_out":2395,"would_cite":true,"duration_ms":23605,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62J12"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper introduces hdbayes, an R package that implements seven historical-data borrowing priors for generalized linear models behind one consistent glm-style syntax, with Stan-based sampling and marginal-likelihood tools.","keywords":["historical data","power prior","normalized power prior","commensurate prior","latent exchangeability prior","meta-analytic predictive prior","generalized linear models","R package"],"falsifier":"Fit the normalized power prior to a small normal linear model where the normalizing constant Z(a0) is known in closed form, using glm.npp with the package's grid estimates, and compare the resulting posterior for the discounting parameter and coefficients with the exact posterior; any systematic discrepancy that grows as the sample shrinks would show the approximation chain is unreliable.","tokens_in":20581,"feed_emoji":"📦","tokens_out":7429,"duration_ms":76742,"temperature":0.7,"pith_summary":"This paper claims that the main barrier to using historical-data priors in Bayesian regression is software fragmentation, and that a single R package can remove it. hdbayes implements seven such priors—the power prior, normalized power prior, normalized asymptotic power prior, Bayesian hierarchical model, robust meta-analytic predictive prior, commensurate prior, and latent exchangeability prior—for generalized linear models, all with essentially the same glm-style syntax and all sampling done by Stan's Hamiltonian Monte Carlo. The authors demonstrate the package on two clinical trial datasets, including a logistic regression comparison across all priors and a proportional-hazards analysis via the Poisson likelihood representation. If the package works as claimed, applied statisticians can try many borrowing strategies in one workflow and compare them fairly.","feed_headline":"Seven historical-data priors, one glm-style R package","feed_subtitle":"Switch among power, meta-analytic, commensurate, and LEAP priors without changing syntax; Stan does the sampling.","key_machinery":"The load-bearing object is the wrapper function glm.prior(formula, family, data.list, prior.args, ...), which mimics R's glm() and turns any of the seven priors into a Stan program. The computational workhorses are Hamiltonian Monte Carlo/NUTS sampling from Stan and bridge-sampling estimation of normalizing constants, with the grid-plus-LOESS-plus-linear-interpolation approximation used for the normalized power prior. The normalized asymptotic power prior uses the asymptotic normal approximation to the power prior under a log transformation of the dispersion parameter, giving a closed-form prior, while the latent exchangeability prior uses a finite mixture model to discount historical observations individually.","core_discovery":"The central claim is that a cohesive implementation is possible: every one of these priors can be expressed as a small wrapper around a common formula, family, and data-list interface, with sensible defaults, and inference can be carried out in Stan without hand-tuning. The package's key methodological move is to handle the unnormalized priors by estimating normalizing constants with bridge sampling: for the normalized power prior, it samples the power-prior posterior on a grid of discounting values, smooths the estimated log-normalizing constants with LOESS, and feeds the smoothed curve into Stan as a lookup grid with linear interpolation; for the robust MAP prior, it estimates the marginal likelihoods of the informative and vague components and reweights their posterior samples. The paper also introduces a normalized asymptotic power prior that avoids normalizing-constant estimation entirely by using the asymptotic multivariate normal form of the power prior on log-dispersion.","pith_inferences":["The paper does not test the accuracy of the bridge-sampling-plus-smoothing normalizing constant estimates beyond one dataset, so users should check stability across grid choices and MCMC seeds before relying on normalized power prior or robust MAP results.","The same wrapper architecture could be extended natively to longitudinal or time-to-event models, and to propensity-score integrated priors, as the discussion notes.","The consistent interface could serve as a benchmark for future software: any new historical-data prior becomes adoptable if it fits the same formula, family, and data-list contract.","The empirical comparison on the AIDS data suggests that borrowing strength varies by method, so reporting several priors rather than one may be the honest default in clinical analyses."],"forward_implications":["Users can fit all seven priors with one consistent interface, so method comparisons no longer require translating between different packages and sampling schemes.","Marginal likelihoods for unnormalized priors become computable, enabling Bayes-factor model and link selection that was previously awkward.","Time-to-event analyses can be conducted through the Poisson representation of piecewise-constant-hazard proportional hazards models, extending the package beyond GLM responses.","Stan's NUTS sampling removes the need to tune Metropolis or slice samplers, reducing a common obstacle to using these priors in high dimensions.","The package can stack multiple historical data sets for power priors and the latent exchangeability prior, allowing a single analysis to pool several sources."],"supporting_citations":[{"why":"Defines the power prior, the foundation for most priors implemented in the package.","marker":"Ibrahim and Chen (2000)"},{"why":"Shows that treating the discounting parameter as random requires estimating the normalizing constant, motivating the normalized power prior.","marker":"Duan et al., 2006"},{"why":"Provides the two-step normalized-power-prior approach that the package simplifies.","marker":"Carvalho and Ibrahim, 2021"},{"why":"Defines commensurate priors and the spike-and-slab prior on commensurability parameters.","marker":"Hobbs et al., 2012"},{"why":"Introduces the meta-analytic predictive prior and its robust mixture version.","marker":"Schmidli et al., 2014"},{"why":"Defines the latent exchangeability prior implemented as glm.leap.","marker":"Alt et al., 2024a"},{"why":"Bridge sampling is used to estimate normalizing constants for the normalized power prior, robust MAP prior, and marginal likelihoods.","marker":"Gronau et al., 2020"},{"why":"Stan is the probabilistic programming backbone providing Hamiltonian Monte Carlo and NUTS sampling.","marker":"Carpenter et al., 2017"}],"fun_headline_variants":["Six priors, one R package: hdbayes unifies historical data for GLMs","hdbayes: All your historical-data priors in one GLM package","Bridge sampling meets historical priors: the hdbayes R package","From power prior to LEAP: hdbayes standardizes historical data","Six historical-data priors, one R call: hdbayes"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the package's numerical estimates of the constants that normalize the priors—produced by sampling on a grid of discount values, smoothing, and interpolating—remain accurate enough in real problems, since the paper demonstrates this only on a single dataset.","fun_headline_variants_meta":{"raw":{"variants":["Six priors, one R package: hdbayes unifies historical data for GLMs","hdbayes: All your historical-data priors in one GLM package","Bridge sampling meets historical priors: the hdbayes R package","From power prior to LEAP: hdbayes standardizes historical data","Six historical-data priors, one R call: hdbayes"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000304,"raw_usage":{"total_tokens":1693,"prompt_tokens":838,"completion_tokens":855,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":454,"completion_tokens_details":{"reasoning_tokens":751}},"tokens_in":454,"tokens_out":855,"duration_ms":9161,"temperature":1.0,"reasoning_tokens":751,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T22:55:59.196546+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Fit the normalized power prior to a small normal linear model where the normalizing constant Z(a0) is known in closed form, using glm.npp with the package's grid estimates, and compare the resulting posterior for the discounting parameter and coefficients with the exact posterior; any systematic discrepancy that grows as the sample shrinks would show the approximation chain is unreliable.","supporting_citations":[],"review_version":1}