{"id":"3d9846aa-f283-4c40-87ed-1d4aea6202d7","arxiv_id":"1908.06687","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":3.0,"correctness_risk":"low","formal_verification":"none","parameter_count":2,"one_line_summary":"Bayesian survival models fitted with different R packages yield hazard ratio estimates close to maximum likelihood, with two notable exceptions.","lead":"This paper compares R packages for Bayesian survival analysis on a real clinical trial dataset. It finds most give hazard ratio estimates close to standard Cox results, with a few exceptions.","discovery_kind":"review","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Broad 'close agreement' claim rests on a single ALLOZITHRO mock dataset; the methods most sensitive to implementation choices (dynsurv, Murray splines) could behave differently under other event rates, censoring, or non-PH hazards.","rationale":"The paper's central claim is a practical guidance statement about agreement between Bayesian implementations and MLE, and the presented analyses with MCMC diagnostics and sensitivity checks support that claim for the single dataset studied. The weakest assumption is indeed the generalizability from one bootstrap-resampled ALLOZITHRO dataset. The authors acknowledge this in the Limitations, but the abstract and conclusion state the result more broadly. The concern is load-bearing because flexible-baseline methods (dynsurv, Murray LRTP splines, and to a lesser extent rstanarm M-splines) already show sensitivity to implementation choices on this dataset, and those sensitivities are plausibly dataset-dependent. A simulation study across event rates, censoring patterns, sample sizes, and PH violations would directly test whether the agreement pattern persists. This does not move the reader's verdict: the paper remains a useful review-style comparison, but the general guidance should be conditional on the dataset characteristics or backed by such a simulation. No internal mathematical error or unsupported reproducibility claim beyond the missing mock dataset/code is identified, and the reported diagnostics give reasonable assurance for the given example.","tokens_in":25651,"tokens_out":8814,"duration_ms":100920,"concrete_test":"Conduct a small simulation study crossing sample size (N=200, 500), event rate (20%, 50%), true HR (1.0, 1.5, 2.0), and PH versus non-PH (e.g., Weibull with a time-varying treatment effect). For each scenario, generate, say, 200 datasets, fit the same R functions with their default settings, and record the posterior mean log HR for each package. Report the mean and worst-case absolute deviation from the true log HR. If any package shows a worst-case deviation greater than 0.1 log HR (about 10% on the HR scale) in a realistic scenario, the abstract's 'close agreement' claim should be qualified by dataset characteristics; if all packages remain within that band, the concern is resolved.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim that Bayesian PH models implemented in various R packages give posterior HR estimates close to the Cox MLE is supported by one bootstrap-resampled dataset derived from ALLOZITHRO (Section 3: 465 patients, 231 events, HR=1.44, PH test P=0.60). The Limitations state that 'results may be dependent on the dataset,' but the abstract and conclusion present the agreement without that caveat. This is a genuine limitation of the evidence, not an internal inconsistency. The dataset is moderately large and has a clear event rate; a smaller trial with 20% events, heavy censoring, or a mild PH violation could make the flexible-baseline methods—which already produced the largest discrepancies—deviate much more. The Discussion itself shows the dynsurv estimate moving from log HR 0.260 to 0.406 when the interval partition is changed, illustrating that implementation choices interact with data features. Without a simulation or second case study, the generality of the headline conclusion is under-supported. This is the weakest load-bearing premise of the paper's main guidance claim.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper compares Bayesian implementations of proportional hazards survival models available in R, applied to a bootstrap-resampled mock dataset derived from the ALLOZITHRO trial (465 patients, 231 events, Cox log-HR estimate 0.366). The authors review partial-Bayesian conjugate analyses and fully Bayesian parametric, piecewise exponential, and nonparametric baseline models implemented in SemiCompRisks, INLA, survHE, rstanarm, dynsurv, spBayesSurv, and the Murray et al. (2016) spline function. They report posterior mean log-HR estimates mostly in the range 0.34 to 0.36, with exceptions for dynsurv's bayesCox (0.260), Murray et al.'s LRTP splines (0.223), and rstanarm M-splines (0.289), and they provide MCMC diagnostics and a sensitivity analysis with homogenized MCMC parameters. The paper concludes that Bayesian PH models in R provide results in close agreement with maximum likelihood estimates and enable probabilistic statements not otherwise available.","tokens_in":25855,"tokens_out":4302,"duration_ms":41992,"significance":"If the central claim holds, the paper provides useful practical guidance for clinical trialists who want to perform post hoc Bayesian analyses of survival endpoints. The implementation is careful in several respects: convergence diagnostics are reported for every MCMC run, a sensitivity analysis with optimized sampling parameters is included, and the package-by-package parameterization is documented in detail. The paper also honestly identifies the dataset dependence as a limitation. However, the evidence base is a single mock dataset, two of the implementations produce estimates that are not in close agreement, and the mechanisms behind those discrepancies are not explained. The useful comparative material is therefore accompanied by a generalizability claim that is currently under-supported.","major_comments":[{"comment":"The central claim of close agreement between Bayesian implementations and the Cox MLE rests on a single bootstrap-resampled dataset derived from the ALLOZITHRO trial. The authors acknowledge in the Limitations that 'results may be dependent on the dataset,' but the abstract and conclusion state the agreement without that caveat. Since two of the implementations (dynsurv bayesCox and Murray et al. LRTP splines) already deviate substantially on this dataset, a simulation study varying event rates, censoring proportions, and the degree of PH violation, or a second real-data case study, is needed before the headline guidance can be considered generalizable. Without such evidence, the claim that Bayesian PH models 'providing results in close agreement with the maximum likelihood estimates' is too strong for the presented support.","section":"Section 3, Limitations"},{"comment":"There is a direct inconsistency between the reported estimand and the package output. Table 2 states that SemiCompRisks' BayesSurv_HReg() 'reports posterior median estimates (not mean),' yet Table 3, Figure 2, and the abstract present all results under the heading 'Mean logHR.' The SemiCompRisks Weibull and PEM entries are therefore likely posterior medians rather than means, which compromises the comparison of posterior means that is the paper's central quantity. The authors should either compute posterior means from the MCMC samples or clearly relabel these entries as medians throughout, and the abstract's 'posterior mean hazard ratio' wording would then need adjustment.","section":"Table 2 and Table 3"},{"comment":"The exceptions to the agreement claim are not explained mechanistically. The dynsurv default partition moves the log-HR estimate from 0.260 to 0.406 when the partition is changed (Discussion), and Table 2 notes that dynsurv's default MCMC output includes burn-in iterations; the authors attribute the deviation only to the larger number of intervals without testing whether the untrimmed MCMC output contributes. Similarly, the Murray et al. LRTP spline estimate of 0.223 with a 95% CrI covering zero is reported without an explanation of whether the discrepancy reflects the spline prior, the partition choice, or the time rescaling mentioned in the Discussion. Because the paper's purpose is to guide implementation choices, unexplained variability driven by implementation details in exactly these methods is a load-bearing gap.","section":"Section 3.2 and Discussion"}],"minor_comments":[{"comment":"The function name 'indepdtCoxph()' is a typo; the correct spBayesSurv function is 'indeptCoxph()', and Figure 2 uses the lower-case label 'indeptcoxph'. Please standardize.","section":"Section 3.2"},{"comment":"The caption of Figure 3 says the middle and right panels refer to an exponential model via survHE/Stan and an MPT prior, respectively, but the plotted middle panel is labeled 'SemiCompRisks PEM' and the right panel 'MPT - Weibull'. The caption does not match the figure content.","section":"Figure 3 caption"},{"comment":"The text says 'The pnorm() and pqnorm() functions in R enable...' but the second function should be qnorm(); pnorm and qnorm are the relevant functions for the posterior probability and quantile computations described.","section":"Section 2.1"},{"comment":"The SemiCompRisks prior is listed in Table 1 as 'β ∝ 1', which is an improper flat prior, while Section 3.2 describes it as 'uniform with mean zero'; this wording conflates an improper prior with a proper uniform distribution and should be corrected.","section":"Table 1 and Section 3.2"},{"comment":"The caption refers to 'survreg2(..., survmodel=\"PH\", ...)', but the function name in the package is survregbayes2(); please correct the caption.","section":"Figure 5 caption"},{"comment":"The name 'Kalbfleish' is a typo for Kalbfleisch, matching reference 21; please correct the spelling.","section":"Section 2.2.1"}],"recommendation":"major_revision","confidential_remarks":"The median-versus-mean issue in the SemiCompRisks entries should be verified against actual package output; if confirmed, it directly affects the abstract's wording. The single-dataset limitation is acknowledged by the authors, but a simulation component may be substantial; if the authors can add even a small simulation study, the paper would be considerably stronger. The paper also relies on a development version of rstanarm, which may date quickly, but that is a normal caveat for software-comparison papers."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"What you should know: this is a practical comparison paper, not a methods paper. The authors fit a large set of Bayesian proportional-hazards survival models in R to one mock dataset derived from the ALLOZITHRO trial, and show that most posterior hazard ratio estimates land near the Cox MLE. The genuinely useful part is the survey: Table 1 gives a clear, side-by-side view of baseline hazard priors and log-HR priors across packages; the MCMC diagnostics are reported; and a sensitivity analysis with homogenized/optimized MCMC parameters confirms the main results. That is real service for practitioners who need to pick a package and want to know what defaults do.\n\nThe soft spots are in proportion to how soft they are. The headline claim—'Bayesian models provide results in close agreement with maximum likelihood'—rests on a single bootstrap-resampled dataset, 465 patients and 231 events. The authors do acknowledge in the Limitations that results may depend on the dataset, but the abstract and conclusion do not carry that caveat. The two methods that disagree (dynsurv's bayesCox and Murray et al.'s LRTP spline model) are not explained mechanistically; the dynsurv discrepancy is shown to be driven by the time-interval partition, moving from log-HR 0.260 to 0.406 with a 5th-percentile partition. That is an implementation sensitivity, not a statistical paradox, but without simulation we have no idea how these methods behave under heavy censoring, low event counts, or non-proportional hazards. Also, no code or data is provided, even the mock dataset, which undercuts the paper's value as a reproducible reference. A minor inconsistency: Table 2 notes that SemiCompRisks reports posterior medians, while the text calls the same numbers posterior means; presumably the authors computed means from the MCMC samples themselves, but they should say so.\n\nWho is this for? Someone who wants to run a Bayesian PH survival analysis in R and needs a map of which packages exist and what their defaults do. It is not a methodological advance and does not settle any open statistical question. I would not cite it in my own work, but I would point a colleague to it. The evaluation is careful and honest enough that the paper deserves a serious referee, not a desk reject. My recommendation: send it to peer review, but ask the authors to soften the abstract's generality, release the mock data and analysis code, and if possible add one simulation or a second case study. That would turn a useful reference into a defensible one.","headline":"Useful practical map of Bayesian survival PH packages in R, but the 'close agreement' headline is under-supported by a single mock dataset and two unexplained outlier methods.","tokens_in":26358,"tokens_out":2043,"would_cite":false,"duration_ms":23499,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62N01","62P10"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper shows that Bayesian survival models in R, despite very different baseline-hazard priors and samplers, give hazard-ratio estimates close to the Cox maximum-likelihood value, plus probabilities the frequentist approach lacks.","keywords":["Bayesian survival analysis","proportional hazards model","clinical trials","hazard ratio","R software packages","Markov chain Monte Carlo","baseline hazard priors","post hoc analysis"],"falsifier":"Re-run the same set of R packages on trials with heavy censoring (around 70–80%), small event counts, or crossing survival curves that clearly violate the proportional-hazards assumption, and compare each posterior mean log hazard ratio with the Cox maximum-likelihood estimate. The paper's claim would be undercut if the inter-package spread widens materially beyond the range it reports (about 0.22 to 0.41, including the partition-shifted dynsurv estimate), or if a single default setting moves an estimate as much as dynsurv's time partition did (a shift of about 0.14 on the log scale).","tokens_in":25441,"feed_emoji":"📊","tokens_out":22476,"duration_ms":180345,"temperature":0.7,"pith_summary":"This paper aims to establish that a post hoc Bayesian analysis of a randomized trial with a right-censored survival endpoint is practical in standard R software and stable across implementations: the choice of package changes the treatment-effect estimate only slightly. Fitting proportional-hazards models to a mock dataset resampled from the ALLOZITHRO trial, across packages that differ sharply in how they model the baseline hazard (parametric exponential or Weibull, piecewise exponential, transformed Bernstein polynomials, mixtures of Polya trees, splines) and in how they sample (JAGS, Stan, integrated nested Laplace approximation), the authors obtain posterior mean hazard ratios that converge on the frequentist maximum-likelihood estimate of HR = 1.44. With a few documented exceptions traced to time-partition defaults, these estimates cluster between 1.39 and 1.43 on the hazard-ratio scale. The payoff is that Bayesian survival reanalysis can be offered as a safe complement to the standard Cox analysis in trials, yielding probabilistic statements such as $\\Pr(\\mathrm{HR} > 1.5)$ that the frequentist approach cannot produce.","feed_headline":"Most Bayesian survival models in R match the classical estimate","feed_subtitle":"Posterior hazard ratios stayed near the frequentist value of 1.44 across every baseline model family tested.","key_machinery":"The load-bearing object is the Cox proportional-hazards model $H(t \\mid Z_i) = \\exp(Z_i'\\beta)\\,H_0(t)$, with the treatment effect $\\beta = \\log(\\mathrm{HR})$ as the parameter of interest and the baseline hazard $H_0$ (equivalently $S_0$ or $h_0$) left to be specified. The comparison is carried by how each implementation parameterizes $H_0$: fully parametric (exponential or Weibull), piecewise exponential with a time partition, or nonparametric through transformed Bernstein polynomial priors, mixtures of Polya trees, or splines on the log-hazard — always paired with a normal prior on $\\beta$ whose variance ranges from 1 to $10^{10}$. For the partial Bayesian approach the machinery is instead a conjugate normal update, $\\beta \\mid y \\sim N\\big((\\mu\\sigma^2 + y\\sigma_0^2)/(\\sigma^2 + \\sigma_0^2),\\, \\sigma^2\\sigma_0^2/(\\sigma^2 + \\sigma_0^2)\\big)$, which requires no sampling. Estimation otherwise proceeds by MCMC (JAGS or Stan) or by integrated nested Laplace approximation. The mechanism producing the paper's conclusion is that the likelihood dominates the prior on $\\beta$: as long as that prior is flat or weakly informative, the posterior mean lands near the maximum-likelihood value no matter how the baseline hazard is modeled, and the one modeling choice that can genuinely move the estimate is the partition of the time axis in piecewise models.","core_discovery":"The paper's central claim, on its own terms, is that Bayesian proportional-hazards models for right-censored clinical-trial survival data can be implemented through various R packages, and these implementations — despite heterogeneity in the modeling of the baseline hazard, in the priors on the treatment effect, and in the estimation machinery — produce posterior mean estimates of the log hazard ratio that agree closely with the maximum-likelihood estimate of the semiparametric Cox model. On the ALLOZITHRO-derived example the frequentist estimate is $\\log(\\mathrm{HR}) = 0.366$, i.e. HR = 1.44 (95% CI 1.11 to 1.87). The parametric and piecewise-exponential implementations give posterior means between 0.345 and 0.362, the nonparametric-baseline implementations (transformed Bernstein polynomial or mixture of Polya trees priors) between 0.337 and 0.358, and the conjugate partial-Bayesian analysis between 0.332 and 0.351, corresponding to hazard ratios from 1.39 to 1.43. Three implementations fall outside that cluster: the dynsurv package's bayesCox function (0.260), rstanarm's M-spline model (0.289), and Murray et al.'s LRTP-spline model (0.223). The paper traces the dynsurv discrepancy to its default time partition built from observed failure times, and shows that re-parameterizing the partition to quantile-based or equally spaced intervals moves that estimate to about 0.40, toward the other methods. The authors conclude that Bayesian approaches can be used more routinely in survival trials, with the added value of probabilistic statements such as $\\Pr(\\mathrm{HR} > 1.5)$ that cannot be obtained from the frequentist analysis alone.","pith_inferences":["An implication the paper leaves implicit: its own Table 3 shows that the quantities that justify Bayesian analysis, posterior probabilities such as $\\Pr(\\mathrm{HR} > 1.5)$, vary much more across implementations (from about 0.08 to 0.38) than the hazard-ratio point estimates do, so the choice of package matters more for the probabilistic statements than for the effect estimate.","Because the posterior mean was insensitive to prior variance from 1 up to $10^{10}$, the agreement plausibly extends to other weakly informative priors; it says nothing about informative priors that conflict with the data, which would be the natural test of where the agreement breaks.","A testable extension is to run the same comparison on trials with non-proportional hazards or heavy censoring, since the paper's conclusion that implementation choice barely matters would be expected to fail exactly where the Cox model itself is a poor summary of the data."],"forward_implications":["A clinical team performing a post hoc Bayesian reanalysis of a randomized trial with a survival endpoint can use almost any of the reviewed R packages and expect a posterior mean hazard ratio within a few percent of the Cox maximum-likelihood estimate on the hazard-ratio scale.","Bayesian implementations add decision-relevant probabilistic statements — such as the posterior probability that the hazard ratio exceeds a threshold like 1.5 — that the frequentist analysis cannot supply, making their added value interpretive rather than a different point estimate.","In piecewise-exponential models the time partition is the parameterization that matters: default partitions built from observed failure times can shift the estimate (dynsurv moved from 0.260 to about 0.40), while quantile-based or equally spaced partitions bring it into line with the other methods.","The conjugate partial-Bayesian approach with skeptical and enthusiastic priors is a reliable and computation-free route to Bayesian estimates under proportional hazards, which supports its use as a simple default for trial reanalysis.","If the agreement generalizes, Bayesian analysis of survival endpoints can become a standard complement to the primary frequentist analysis of confirmatory trials, converting a hazard-ratio estimate into directly usable probability statements for clinicians and policymakers."],"supporting_citations":[{"why":"The ALLOZITHRO randomized trial report; supplies the trial whose resampled mock dataset grounds the entire comparison.","marker":"9"},{"why":"Supplies the normal prior on the log hazard ratio whose variance encodes the implicit sample size, the basis of the partial Bayesian approach.","marker":"10"},{"why":"Supplies the skeptical and enthusiastic priors and the conjugate normal posterior update that constitute the recommended non-fully-Bayesian reference method.","marker":"12"},{"why":"Supplies the integrated nested Laplace approximation engine used for the exponential, Weibull, and random-walk piecewise-exponential models.","marker":"16"},{"why":"Supplies the dynsurv bayesCox() function whose default failure-time partition produces the main outlying estimate and motivates the partition sensitivity analysis.","marker":"17"},{"why":"Supplies the spBayesSurv implementations — piecewise exponential, transformed Bernstein polynomial, and mixture of Polya trees baseline priors.","marker":"18"},{"why":"Supplies the survHE interface through which the exponential and Weibull models are fit with Hamiltonian Monte Carlo/Stan and INLA.","marker":"19"},{"why":"Supplies SemiCompRisks' BayesSurv_HReg() used for the parametric Weibull and piecewise exponential models.","marker":"20"},{"why":"Supplies the flexible LRTP-spline model whose estimate falls outside the agreement cluster and whose partition-choice guidance informs the sensitivity analysis.","marker":"29"},{"why":"Supplies the proportional-hazards diagnostic test whose non-significant result (P = 0.60) justifies the PH modeling setting.","marker":"31"}],"fun_headline_variants":["Bayesian survival models in R mostly match Cox, one lags","One R package breaks rank in Bayesian survival analysis","Bayesian HR close to frequentist across most R packages","R survival Bayesian models: consistent HR, except dynsurv","Bayesian survival in R: near Cox agreement, one outlier"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise, acknowledged by the authors in their Limitations paragraph ('results may be dependent on the dataset'), is that the single mock dataset resampled from the ALLOZITHRO trial represents typical right-censored trial data; if other censoring rates, event proportions, or departures from proportional hazards change how the packages behave, the close agreement between Bayesian and frequentist estimates may not transfer.","fun_headline_variants_meta":{"raw":{"variants":["Bayesian survival models in R mostly match Cox, one lags","One R package breaks rank in Bayesian survival analysis","Bayesian HR close to frequentist across most R packages","R survival Bayesian models: consistent HR, except dynsurv","Bayesian survival in R: near Cox agreement, one outlier"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000573,"raw_usage":{"total_tokens":2770,"prompt_tokens":1068,"completion_tokens":1702,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":684,"completion_tokens_details":{"reasoning_tokens":1619}},"tokens_in":684,"tokens_out":1702,"duration_ms":11941,"temperature":1.0,"reasoning_tokens":1619,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T12:36:20.287745+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Re-run the same set of R packages on trials with heavy censoring (around 70–80%), small event counts, or crossing survival curves that clearly violate the proportional-hazards assumption, and compare each posterior mean log hazard ratio with the Cox maximum-likelihood estimate. The paper's claim would be undercut if the inter-package spread widens materially beyond the range it reports (about 0.22 to 0.41, including the partition-shifted dynsurv estimate), or if a single default setting moves an estimate as much as dynsurv's time partition did (a shift of about 0.14 on the log scale).","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"The ALLOZITHRO randomized trial report; supplies the trial whose resampled mock dataset grounds the entire comparison."},{"cited_title":"Incorporating Bayesian ideas into health-care evaluation","cited_arxiv_id":null,"evidence_quote":"Supplies the normal prior on the log hazard ratio whose variance encodes the implicit sample size, the basis of the partial Bayesian approach."},{"cited_title":"& Parmar, M","cited_arxiv_id":null,"evidence_quote":"Supplies the skeptical and enthusiastic priors and the conjugate normal posterior update that constitute the recommended non-fully-Bayesian reference method."},{"cited_title":"& Yan, J","cited_arxiv_id":null,"evidence_quote":"Supplies the dynsurv bayesCox() function whose default failure-time partition produces the main outlying estimate and motivates the partition sensitivity analysis."},{"cited_title":"spBayesSurv: Fitting Bayesian Spatial Survival Models Using R","cited_arxiv_id":"1705.04584","evidence_quote":"Supplies the spBayesSurv implementations — piecewise exponential, transformed Bernstein polynomial, and mixture of Polya trees baseline priors."},{"cited_title":"survHE: Survival analysis in health economic evaluation","cited_arxiv_id":null,"evidence_quote":"Supplies the survHE interface through which the exponential and Weibull models are fit with Hamiltonian Monte Carlo/Stan and INLA."},{"cited_title":"H., Lee, C., Alvares, D","cited_arxiv_id":null,"evidence_quote":"Supplies SemiCompRisks' BayesSurv_HReg() used for the parametric Weibull and piecewise exponential models."},{"cited_title":"A., Hobbs, B","cited_arxiv_id":null,"evidence_quote":"Supplies the flexible LRTP-spline model whose estimate falls outside the agreement cluster and whose partition-choice guidance informs the sensitivity analysis."}],"review_version":1}