Using stan_jm for parametric proportional hazards regression only

This is an awesome idea! Just to make sure I understand correctly (and for future readers) let me try and formulate the idea in my own words (or mathematical notation):

  • W.l.o.g. set the origin of time to 0.
  • Suppose the event (either 0 or 1) for individual i has occurred at time t_i.
  • Suppose there are n_i unique event times before t_i and denote them by t_1 < t_2 < \dots <t_{n_i}.
  • For individual i we construct n_i+1 observations of the form (t_1, 0), ( t_2 - t_1, 0), (t_3-t_2, 0),…, (t_{n_{i}} - t_{n_{i-1}},0) and (t_i - t_{n_{i}},1). The first value in the tuples is the length of the interval and the second the number of events from individual i in the corresponding interval.
  • We model these induced (or transformed) observations for individual i through n_{i} + 1 independent Poisson regressions with rates \lambda(t) = \lambda_0(t) e^{x_i' \beta} and additional offset accounting for different interval lengths (aka “exposure times” in the Poisson regression language). Here t is any reference point in the corresponding interval (e.g. the start or the end of the interval) and the functional form of the baseline hazard rate \lambda_0 can be modelled through a smoothing spline. [Here x_i is a vector containing covariates for individual i and \beta a latent parameter column vector].