Pith. sign in

REVIEW 3 major objections 7 minor 33 references

Capturing heterogeneous time-variation in covariate effects in non-proportional hazard regression models

T0 review · 3 major / 7 minor · reviewed 2026-08-10 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read The paper introduces a hazard-regression term $f_g(t)\cdot x$ that gives each subgroup its own non-linear, time-varying covariate effect, and shows it fits better than nested alternatives while shrinking when the effect is absent.

desk verdict A solid, incremental transfer of functional random coefficients to survival analysis; the machinery is standard, the simulation supports the main claims, but the exchangeability assumption for fixed subgroups is the real soft spot and is not stress-tested. read the letter →

arxiv 2501.13525 v1 pith:RYT3W5Q5 submitted 2025-01-23 stat.ME math.STstat.TH

Niklas Hagemann , Thomas Kneib , Kathrin Möllenhoff This is my paper · ORCID
classification stat.MEmath.STstat.TH MSC 62N0162N0262G08
keywords survivalanalysishazardregressionpiecewiseexponentialadditivemixedmodelstime-varyingcovariateeffectsfunctionalrandomfactorsmoothsnon-proportionalhazardspenalizedsplines
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Survival analysts often need covariate effects that change over time and differ across subgroups, sometimes with the time pattern itself differing across subgroups. The paper proposes to capture that third layer—heterogeneously time-varying effects—by adding a functional random coefficient to a piecewise exponential additive mixed model. The term $f_g(t)\cdot x_{ik}$ lets each group have its own non-linear time curve for a covariate's effect, built as a tensor-product interaction of a random group effect and a penalized spline in time. Simulations show this term fits better when such effects exist and is penalized toward the simpler nested model when they do not; the brain-tumor example finds that fraction genome altered affects survival with direction and timing that differ by diagnosis.

What carries the argument

The central object is the functional random coefficient, implemented as $s(g, t, \text{by}=x, \text{bs}=\text{fs})$: a functional random effect (factor smooth) used as a varying coefficient. It is an anisotropic tensor-product interaction of a group-index random effect (indicator basis with identity penalty) and a penalized B-spline in time, built with the numerically stable reparameterization of tensor-product smooths. This object carries the argument because it is the only term in the predictor that lets a covariate's effect have both a group-specific level and a group-specific non-linear time curve, while its penalty structure automatically shrinks the term toward simpler nested alternatives.

What would settle it

Simulate data in which the four group effects are fixed, widely separated curves with no shared variance structure, fit both the functional random coefficient and a fixed group-by-time interaction, and check whether the random-effect version visibly shrinks the outer curves toward a common mean; if it does, the exchangeability assumption fails in exactly the cases the method is meant to expose.

Watch

Extended reading notes

Core claim

The central claim is that a piecewise exponential additive mixed model (PAMM) can include covariate effects of the form $f_g(t)\,x_{ik}$, where $g$ indexes a subgroup and $t$ is time, using a functional random effect—a tensor product of an i.i.d. random intercept in $g$ and a P-spline in $t$. Because the effect is a separate curve per group, it can capture not only subgroup-specific time variation but subgroup-specific changes in that time variation, something no single existing survival framework did. The paper argues that the penalized spline construction keeps the term from overfitting when the true effect is simpler, shrinking toward the main-effect-plus-time model, and demonstrates the point with four simulation scenarios and a diagnosis-specific analysis of fraction genome altered in gliomas. On the case study, modeling FGA with this term improves fit over nested alternatives and reveals effects that would cancel out if the time variation were pooled across diagnoses.

Load-bearing premise

The method treats the subgroup labels as if they are random draws from one common distribution with a shared spread and shared smoothness, so each subgroup's curve is pulled toward the average rather than estimated as its own fixed curve.

Editorial extensions

If this is right

  • If correct, researchers can fit non-proportional hazard models where a covariate's effect curve is allowed to differ per subgroup without specifying which subgroup has which curve; the data choose.
  • The penalization means the same software can serve as a model-selection device: when heterogeneous time-variation is absent, the term shrinks toward the simpler model, so choosing interval cut-points or testing for such effects becomes less delicate.
  • The brain-tumor analysis suggests that pooling diagnosis-specific time-varying effects can hide real signals—FGA's effect declines for glioblastoma but rises for other diagnoses—so a single average time-varying coefficient can be misleading.
  • Because the term lives inside the Poisson-likelihood reparameterization of PAMMs, it inherits existing GAM software and can be applied to large survival datasets.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Inference: the random-effects assumption behind the group term is the part most in need of stress-testing; fixed subgroups with systematically different curves may be over-shrunk, so a fixed-effects version or a prior that allows curve similarity is a natural comparison.
  • Inference: the construction could transfer to settings where the grouping variable is not categorical—for example a continuous moderator—by replacing the random intercept in $g$ with a smooth function of that moderator.
  • Inference: one testable extension is to use the fitted group curves to build a clustering or classification rule, since the method yields one full hazard-effect trajectory per subgroup that could feed downstream analyses.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 7 minor

Summary. The paper proposes modeling subgroup-specific time-varying covariate effects in survival analysis via functional random coefficients of the form f_g(t) * x_ik, embedded in piecewise exponential additive mixed models (PAMMs). The functional random effect is implemented as a factor-smooth tensor product interaction using the mgcv specification s(g, t, by = x, bs = "fs"), with penalized P-splines providing regularization. The authors argue that this is the first unified framework for heterogeneously time-varying effects in hazard regression. They support the proposal with a simulation study comparing four nested models under four data-generating scenarios, and with a brain tumor case study examining the diagnosis-specific time-varying effect of fraction genome altered (FGA) on survival.

Significance. The methodological core is sound and practically valuable: the paper connects a recently proposed functional random coefficient construction to the well-established PAMM/Poisson-likelihood framework, and the implementation relies on existing, widely used software (mgcv, pammtools), which lowers the barrier for adoption. The simulation study is extensive (1000 repetitions, three sample sizes, four scenarios) and the authors appropriately acknowledge the one case of slight overfitting. The case study is clinically relevant and produces interpretable diagnosis-specific effect curves. If the exchangeability assumption of the group random effects is accepted, the contribution is a useful extension of PAMMs. However, the principal evidence for the central claim is a simulation whose data-generating process is the proposed model itself, and the case study uses five fixed diagnosis groups, so the paper does not yet demonstrate that the method is reliable when the exchangeability assumption is violated.

major comments (3)
  1. [Section 4.2, Eq. (5)] The simulation DGP in scenario (I) generates the group-specific time-varying effect f_g(t) from the same functional random coefficient family that the proposed model assumes. Therefore, scenario (I) only demonstrates that the estimator can recover curves generated under its own exchangeability assumption; it does not test the method's ability to capture heterogeneous time-variation for fixed subgroups, which are an explicit target in the Introduction (Section 1, final paragraph) and in the case study. I recommend adding a fixed-effects competitor such as s(t, by = diagnosis, bs = "ps") or a tensor-product smooth with diagnosis as a fixed factor, and adding a simulation scenario in which the true group curves are deterministic and systematically different across groups to assess shrinkage bias.
  2. [Section 5, Figure 3 and Table 3] With only G = 5 diagnosis groups, the single common variance component and shared smoothness penalty in Eq. (4) imply potentially strong shrinkage of all five curves toward a common mean, and the paper does not report the estimated variance component, the effective degrees of freedom, or pointwise confidence intervals for the five curves. Without a comparison to a fixed-effects specification, the conclusion that the FGA effect 'strongly varies' between diagnoses and that this variation is not a shrinkage artifact is not yet supported. Please report the variance estimate and the effective degrees of freedom, add pointwise confidence intervals to Figure 3, and fit a fixed-effects interaction model to the case study as a sensitivity check.
  3. [Section 4.3 and Abstract] The abstract states that 'using a penalized basis prevents overfitting in case of absence of such effects,' but the simulation results in scenario (III) show that the proposed model fits slightly better than the true nested model, which the authors acknowledge as slight overfitting. The suggested remedy of visual inspection is subjective and not quantified. This weakens the overfitting-prevention claim as stated. Please provide a quantitative measure, for example the proportion of simulation runs in which the functional random coefficient is estimated to have non-negligible time-variation under a null DGP, or an explicit discussion of how the abstract should be read in light of the scenario (III) finding.
minor comments (7)
  1. [Section 3] Typo: 'the reader is refereed to' should be 'the reader is referred to.'
  2. [Section 4.1] Typo: 'implemented in in the R package mgcv' contains a duplicated 'in.'
  3. [Section 6] The conclusion contains 'the simulation study outlines the superior fit of your approach'; 'your' should be 'our.'
  4. [Section 5, Table 3] The model labeled 'Heterogeneity and time-variation' is not fully specified in Section 5; please clarify its predictor structure, particularly whether it is the same as model (ii) in Section 4.2.
  5. [Section 5, Figure 3] Adding pointwise confidence intervals to the estimated FGA curves would help the reader judge the statistical evidence for the between-diagnosis differences, which is currently only visually implied.
  6. [Section 2.2 and 4.3] The claim that penalization 'mostly solves' the choice of the number of intervals is not accompanied by a sensitivity analysis with respect to the number of intervals or the number of inner knots; a brief simulation or a discussion would make this statement more precise.
  7. [Section 1 and 3] The statement that this is the first proposal of heterogeneously time-varying covariate effects in hazard regression models should be tempered, since Hagemann et al. (2024) already propose the same effect type in conditional logit models; the novelty is the application to hazard regression, not the effect type itself.

Circularity Check

0 steps flagged · score 2.0 of 10

No load-bearing circularity; the method is an explicit model specification, and the self-referential simulation DGP is a standard estimation check rather than a derivation of the target result.

full rationale

The paper's construction is explicit: Eq. (4) gives the random-effect basis f(g)=gamma_{1g}, and the functional random coefficient is the tensor-product interaction of this basis with a P-spline in t, implemented as s(g,t,by=x,bs='fs'). No equation defines the claimed effect in terms of the estimated output; the likelihood and REML estimation are imported from external PAMM/mgcv machinery (Bender et al. 2018; Wood 2011). The self-citations to Kneib et al. (2019) and Hagemann et al. (2024) provide the FRE construction, but this is standard tensor-product smoothing (also cited to Wood 2017) and is not invoked as a uniqueness theorem or to exclude alternatives. The simulation in scenario (I) is self-referential in design: the DGP is exactly model (i), f(x2,t,g)=f_g(t)*x2, so the better fit of model (i) over nested competitors is expected by construction; this is a standard check that the smoother can recover its own target, not a derivation of the method from data. Scenarios (II)-(IV) and the real-data case study provide non-tautological evidence on shrinkage and practical performance. The exchangeability assumption for the group index is a modeling choice, not a circular step; the paper does not claim to derive it. Overall, the central claim retains independent content, and no enumerated circularity pattern is present beyond a minor self-referential simulation design.

Assumptions & free parameters 1 free parameters · 6 assumptions · 0 invented entities

The method does not introduce new physical or mathematical entities; it relies on established PAMM, Poisson-likelihood, penalized-spline, and random-effects theory. The only hand-chosen modeling inputs are spline basis dimensions and penalty orders, which are penalized and therefore not central free parameters. The main conceptual assumption is exchangeability of the grouping levels.

free parameters (1)
  • Basis dimension / inner knots for P-splines = 9 inner knots (case study); defaults for simulations
    The authors choose cubic P-splines with first-order penalty and 9 inner knots in the case study; simulations rely on default or unspecified mgcv dimensions. The penalization is intended to reduce sensitivity to this choice, but the paper does not show a systematic sensitivity analysis.
assumptions (6)
  • domain assumption The hazard is piecewise constant within each time interval (Section 2.1, Eq. 2).
    This is the core PAMM assumption; the approximation quality depends on the interval widths, which are chosen to be all unique survival times in PAMM practice.
  • domain assumption Censoring is non-informative and independent of the survival time (Section 2.1).
    Required for the likelihood factorization used in the data restructuring and Poisson-likelihood inference.
  • standard math The survival likelihood is proportional to a Poisson regression likelihood under the piecewise exponential model (Section 2.3, Eq. 3).
    Proven result from Bender et al. (2018); the paper relies on it to estimate PAMMs with mgcv.
  • domain assumption Grouping levels are exchangeable random effects with a common variance and shared smoothness penalty (Section 3, Eq. 4).
    The functional random effect uses i.i.d. random intercepts for groups and a shared P-spline penalty; if groups are fixed categories, the shrinkage behavior changes and may distort the estimated curves.
  • domain assumption The time-varying effects are smooth and can be represented by P-splines with a first-order penalty (Sections 3 and 5).
    Smoothness and penalty order are modeling choices, not derived from the data; they affect the estimated functional random coefficients.
  • domain assumption REML smoothing parameter estimation and Wood (2012) approximate p-values are valid for the fitted PAMMs (Sections 4.1 and 5).
    Inference relies on standard asymptotic approximations implemented in mgcv; these are widely used but not exact in small samples or with strong random-effect shrinkage.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Capturing heterogeneous time-variation in covariate effects in non-proportional hazard regression models." pith.science (2026). https://pith.science/paper/RYT3W5Q5

@misc{pith2026250113525,
  author       = {Pith},
  title        = {Pith review of: Capturing heterogeneous time-variation in covariate effects in non-proportional hazard regression models},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RYT3W5Q5}},
  note         = {Machine review of arXiv:2501.13525}
}
read the original abstract

A central focus in survival analysis is examining how covariates influence survival time. These covariate effects are often found to be either time-varying, heterogeneous - such as being specific to patients, treatments, or subgroups - or exhibit both characteristics simultaneously. While the standard model, the Cox proportional hazards model, allows neither time-varying nor heterogeneous effects, several extensions to the Cox model as well as alternative modeling frameworks have been introduced. However, no unified framework for incorporating heterogeneously time-varying effects of covariates has been proposed yet. Such effects occur when a covariate influences survival not only in a heterogeneous and time-varying manner, but when the time-variation is also heterogeneous. We propose to model such effects by introducing heterogeneously time-varying coefficients to piecewise exponential additive mixed models. We deploy functional random effects, also known as factor smooths, to model such coefficients as the interaction effect of heterogeneity and time-variation. Our approach allows for non-linear time-effects due to being based on penalized splines and uses an efficient random effects basis to model the heterogeneity. Using a penalized basis prevents overfitting in case of absence of such effects. In addition, the penalization mostly solves the problem of choosing the number of intervals which is usually present in unregularized piecewise exponential approaches. We demonstrate the superiority of our approach in comparison to competitors by means of a simulation study. Finally, the practical application and relevance are outlined by presenting a brain tumor case study.

Figures

Figures reproduced from arXiv: 2501.13525 by the authors.

Figure 1
Figure 1. Effect of x2 in the DGPs of the four scenarios. The first subplot (from the left) shows scenario (I), i.e. heterogeneous time-variation, the second one scenario (II), i.e. the combination of heterogeneity and time-variation (without interaction), the third one scenario (III), i.e. heterogeneity only, and the last one scenario (IV), i.e. time variation only [PITH_FULL_IMAGE:figures/full_fig_p008_1.png] view at source ↗
Figure 2
Figure 2. Results of the simulation study with n = 400 in terms of the three fit measures. For each of the scenarios (I) - (IV) there is one block of consisting of four boxplots, one for each model (i) - (iv). 9 [PITH_FULL_IMAGE:figures/full_fig_p009_2.png] view at source ↗
Figure 3
Figure 3. Estimated Smooth effects for the log-baseline function and the functional random coefficient [PITH_FULL_IMAGE:figures/full_fig_p012_3.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

33 extracted references · 32 canonical work pages

  1. [1]

    Andersen, P. K. and Gill, R. D. (1982). Cox's regression model for counting processes: a large sample study. The annals of statistics , pages 1100--1120

    work page 1982

  2. [2]

    Bender, A., Groll, A., and Scheipl, F. (2018). A generalized additive model approach to time-to-event analysis. Statistical Modelling , 18(3-4):299--321

    work page 2018

  3. [3]

    and Scheipl, F

    Bender, A. and Scheipl, F. (2018). pammtools: Piece-wise exponential additive mixed modeling tools. arXiv:1806.01042 [stat]

    arXiv 2018

  4. [4]

    Bender, R., Augustin, T., and Blettner, M. (2005). Generating survival times to simulate cox proportional hazards models. Statistics in Medicine , 24(11):1713--1723

    work page 2005

  5. [5]

    P., Malta, T

    Ceccarelli, M., Barthel, F. P., Malta, T. M., et al. (2016). Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma. Cell , 164(3):550--563

    work page 2016

  6. [6]

    Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) , 34(2):187--202

    work page 1972

  7. [7]

    and Rodriguez-Bravo, V

    Dhital, B. and Rodriguez-Bravo, V. (2023). Mechanisms of chromosomal instability (cin) tolerance in aggressive tumors: surviving the genomic chaos. Chromosome Research , 31(2):15

    work page 2023

  8. [8]

    Eilers, P. H. C. and Marx, B. D. (1996). Flexible smoothing with B -splines and penalties . Statistical Science , 11(2):89--121

    work page 1996

Show all 33 references
  1. [9]

    Fahrmeir, L., Kneib, T., Lang, S., and Marx, B. D. (2022). Regression: Models, Methods and Applications . Springer, Berlin, Heidelberg

    2022

  2. [10]

    Friedman, M. (1982). Piecewise exponential models for survival data with covariates. The Annals of Statistics , 10(1):101--113

    1982

  3. [11]

    Graf, E., Schmoor, C., Sauerbrei, W., and Schumacher, M. (1999). Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine , 18(17-18):2529--2545

    1999

  4. [12]

    Gray, R. J. (1992). Flexible methods for analyzing survival data using splines, with applications to breast cancer prognosis. Journal of the American Statistical Association , 87(420):942--951

    1992

  5. [13]

    Hagemann, N., Guhl, D., Kneib, T., M\"ollenhoff, K., and Steiner, W. (2024). Dynamic heterogeneity in discrete choice experiments. Preprint available at SSRN: 4957076

    2024

  6. [14]

    Hennerfeind, A., Brezger, A., and Fahrmeir, L. (2006). Geoadditive survival models. Journal of the American Statistical Association , 101(475):1065--1075

    2006

  7. [15]

    Hess, K. R. (1994). Assessing time-by-covariate interactions in proportional hazards regression models using cubic spline functions. Statistics in Medicine , 13(10):1045--1062

    1994

  8. [16]

    Hofner, B., Hothorn, T., and Kneib, T. (2013). Variable selection and model choice in structured survival models. Computational Statistics , 28:1079--1101

    2013

  9. [17]

    Hofner, B., Kneib, T., Hartl, W., and K \"u chenhoff, H. (2011). Building cox-type structured hazard regression models with time-varying effects. Statistical Modelling , 11(1):3--24

    2011

  10. [18]

    Jachno, K., Heritier, S., and Wolfe, R. (2019). Are non-constant rates and non-proportional treatment effects accounted for in the design and analysis of randomised controlled trials? a review of current practice. BMC Medical Research Methodology , 19(1):103

    2019

  11. [19]

    and Fahrmeir, L

    Kneib, T. and Fahrmeir, L. (2007). A mixed model approach for geoadditive hazard regression. Scandinavian Journal of Statistics , 34(1):207--228

    2007

  12. [20]

    Kneib, T., Klein, N., Lang, S., and Umlauf, N. (2019). Modular regression - a lego system for building structured additive distributional regression models with tensor product interactions. TEST , 28:1--39

    2019

  13. [21]

    Li, H., Han, D., Hou, Y., Chen, H., and Chen, Z. (2015). Statistical inference methods for two crossing survival curves: a comparison of methods. PLoS One , 10(1):e0116774

    2015

  14. [22]

    R., Nakao, K., Zuraek, M

    Mehta, K. R., Nakao, K., Zuraek, M. B., et al. (2005). Fractional genomic alteration detected by array-based comparative genomic hybridization independently predicts survival after hepatic resection for metastatic colorectal cancer. Clinical Cancer Research , 11(5):1791--1797

    2005

  15. [23]

    Murphy, S. A. and Sen, P. K. (1991). Time-dependent coefficients in a cox-type regression model. Stochastic Processes and their Applications , 39(1):153--180

    1991

  16. [24]

    and Palmgren, J

    Ripatti, S. and Palmgren, J. (2000). Estimation of multivariate frailty models using penalized partial likelihood. Biometrics , 56(4):1016--1022

    2000

  17. [25]

    M., Grambsch, P

    Therneau, T. M., Grambsch, P. M., and Pankratz, V. S. (2003). Penalized survival models and frailty. Journal of Computational and Graphical Statistics , 12(1):156--175

    2003

  18. [26]

    W., Manton, K

    Vaupel, J. W., Manton, K. G., and Stallard, E. (1979). The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography , 16(3):439--454

    1979

  19. [27]

    Wood, S. N. (2003). Thin plate regression splines. Journal of the Royal Statistical Society Series B: Statistical Methodology , 65(1):95--114

    2003

  20. [28]

    Wood, S. N. (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 73(1):3--36

    2011

  21. [29]

    Wood, S. N. (2012). On p-values for smooth components of an extended generalized additive model. Biometrika , 100(1):221--228

    2012

  22. [30]

    Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, Second Edition . Chapman & Hall/CRC Texts in Statistical Science. CRC Press

    2017

  23. [31]

    N., Pya, N., and S\"afken, B

    Wood, S. N., Pya, N., and S\"afken, B. (2016). Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association , 111(516):1548--1563

    2016

  24. [32]

    N., Scheipl, F., and Faraway, J

    Wood, S. N., Scheipl, F., and Faraway, J. J. (2013). Straightforward intermediate rank tensor product smoothing in mixed models. Statistics and Computing , 23(3):341--360

    2013

  25. [33]

    Zucker, D. M. and Karr, A. F. (1990). Nonparametric survival analysis with time-dependent covariate effects: a penalized partial likelihood approach. The Annals of Statistics , 18(1):329--353

    1990

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.