Cox proportional hazards with a flexible baseline

cox() is the flexible-parametric proportional hazards model brms fits: the baseline hazard is an M-spline in time, h0(t) = sum_j s_j M_j(t), and the cumulative baseline hazard is the I-spline integral of the same basis over the same weights, which makes it monotone by construction. The weights s form a simplex, which is what identifies the baseline against the intercept. They are estimated as sbhaz_raw, their Kbhaz - 1 free coordinates in the chart a mo() simplex uses (see "Monotonic effects" in frmtmb::frm()); cox_baseline() returns the simplex itself. The default basis is brms's: df = 5 cubic M-splines with an intercept, internal knots on response quantiles and boundary knots just outside the observed range.

Usage

cox(
  link = "log",
  df = 5,
  degree = 3,
  intercept = TRUE,
  knots = NULL,
  Boundary.knots = NULL,
  bhaz = NULL
)

Arguments

link: The link for the hazard ratio mu: "log" (the default), "identity", "softplus" or "squareplus", the set brms 2.23.0 allows. Quoted or not.
df: The number of M-spline basis functions in the baseline hazard (brms's bhaz(df = 5) default). With knots, it is length(knots) + degree + intercept, and need not be given.
degree: The spline degree of that basis (cubic by default).
intercept: Keep the basis function that is non-zero at the lower boundary knot.
knots: The internal knots of the basis. The default, NULL, puts df - degree - intercept of them on equally spaced quantiles of the observed times, as brms and splines2::mSpline() do. Given, they must be distinct and lie strictly inside the boundary knots.
Boundary.knots: The two boundary knots. The default, NULL, is brms's: the range of the observed times widened by 1/50 of its width at each end, with the lower knot not below zero. Given, they must cover every observed time.
bhaz: Not an argument of brms 2.23.0's cox(). Given, it is refused with the spellings to use instead.

Value

A frmtmb_family object, ready for frmtmb::frm().

Details

The hazard is h0(t) exp(eta), so a coefficient is a log hazard ratio, exactly as in survival::coxph(). Right, left, and interval censoring come through the ordinary cens() addition term: an event contributes the density and a censored observation the survivor function, which is what this family's log-density and log-CDF are. Random effects are the point: time | cens(c) ~ x + (1 | g) is a frailty model, and the Laplace approximation integrates the frailties out.

The baseline is semiparametric only in spirit: it has df parameters, not one per event time, so coefficients agree with coxph() closely rather than exactly. A survival response has no mean, so fitted() and frm_linpred(type = "response") are refused; frm_linpred(type = "link") gives the log hazard ratio. simulate() is not available.

Maximum likelihood often puts one or more baseline weights ON the simplex boundary, at exactly zero: an unpenalized flexible baseline does that. A weight of zero is a finite point of the chart, which the fit reaches with a finite Hessian, and a zero weight the likelihood rises from is a saddle the fit is moved off (frmtmb_control(mo_search = )). brms keeps the weights interior with its dirichlet(1) prior on sbhaz; frmtmb.sample's frm_sample() samples the baseline under the same prior, and set_prior("dirichlet(...)", class = "sbhaz") puts another one on it. Lower df until the baseline is one the data supports, and read cox_baseline() to see which weights are at zero.

Where it lives

cox() was in frmtmb itself before this package existed. brms has it in its core, so family = "cox" and brmsfamily("cox") resolve to this function once this package is loaded, and refuse by name, naming this package, when it is not.

The baseline simplex reaches frmtmb through frmtmb_family(simplex_pars = c(sbhaz = "sbhaz_raw")), and its basis through frmtmb_family(aterm_data = ). brms's bhaz() addition term reaches it through frmtmb_register_aterm(option = ). Those seams are all the family needs from core.

brms's bhaz() addition term

brms 2.23.0 sets the baseline basis on the response side: time | cens(c) + bhaz(df = 6) ~ x with family = cox(). This package takes that spelling as well. bhaz(df =, degree =, intercept =, knots =, Boundary.knots =) sets the argument of the same name of cox(), so these two models are the same model, with the same basis to rounding as brms's standata():

frm(bf(time | cens(c) + bhaz(df = 6) ~ x), family = cox(), data = d)
frm(bf(time | cens(c) ~ x), family = cox(df = 6), data = d)

The arguments of bhaz() are evaluated in the formula's environment. An argument that both bhaz() and cox() set must have the same value in both, or the fit stops and names the two. df and knots must agree too, in either place: length(knots) + degree + intercept basis functions. brms ignores df when knots is given; frmtmb refuses the pair.

Some differences from brms are refusals. bhaz(gr = g) gives each level of g a baseline of its own (a stratified Cox model), and this family has one baseline: fit each stratum as a model of its own, or let g shift the one baseline by adding it to the formula. Boundary.knots must cover every observed time. brms takes boundary knots inside the data and extrapolates the basis there, where the cumulative baseline hazard can be negative, so that the survivor function exceeds one. An argument name that brms drops without a message, such as bhaz(int = FALSE) or bhaz(Bound = c(0, 30)), is refused with the name that was probably meant. A shortened name that brms reads as an argument (bhaz(d = 6) for df, kn = for knots, deg = for degree) is read the same way here. warn.outside is taken and has no effect, as in brms.

brms 2.23.0's cox() has no bhaz = argument (brms before 2.22.0 took the basis that way). cox(bhaz = ) stops and names the two spellings above.

Examples

set.seed(2)
n <- 200
dd <- data.frame(x = rnorm(n))
dd$time <- rexp(n, exp(-0.5 + 0.7 * dd$x))
dd$out <- rbinom(n, 1, 0.3)        # 1 = right censored
cox_fit <- frm(bf(time | cens(out) ~ x), family = cox(), data = dd)
fixef_by_dpar(cox_fit)$mu          # log hazard ratios
#> (Intercept)           x 
#>   2.0927616   0.7388227 
cox_baseline(cox_fit)              # the baseline hazard weights
#>           s1           s2           s3           s4           s5 
#> 2.071234e-02 1.506379e-01 2.563727e-01 1.556952e-12 5.722771e-01