A piece-wise additive model of survival data with non-linear rut

INTRODUCTION

In a ground breaking paper, Cox proposed a model for survival data in the proportional hazard frame work. Sir David Cox observed that if the proportional hazard assumption holds (or, is assumed to hold), then it is possible to estimate the effect of parameter(s) without any consideration of the hazard function.[1] Developing models in complex analysis may violate this assumption as the model may incorporate time-varying effects and covariates, thereby relaxing the proportional hazard assumption, and allowing the hazard ratio to depend on time t. A common approach to model time-varying effects is by piece-wise constant functions, as these are flexible enough to capture any shape of the baseline hazard or covariate effects,[2] for a frequentist and[3] and more recently[4] for a Bayesian approach.

The hazard model was further extended by Hennerfeind et al., Kneib and Fahrmeir,[5,6] to incorporate a flexible spatial generalization of the Cox model, a structured geo-additive predictor, including a spatial component for geographical effects and non-parametric terms for modeling unknown functional forms of the log-baseline hazard rate and non-linear effects of continuous covariates.

Generalized additive models (GAMs) are statistical models in which the conventional multiple linear regression is generalized to permit a much broader class of time-varying and non-linear functional form of continuous covariates and their effects, but still with additive relationships between response and predictor variables. GAMs, derived from the work of Hastie and Tibshirani,[7,8] provide flexible and effective means of moving out of the “linear rut”[9] in which a considerable amount of bio-statistical modeling is still located. Conventionally, non-linearity is handled through transformations and then estimation of linear models. Such an approach requires the researcher to have good knowledge of the correct functional form before the model is fitted when, in reality, the choice is very open and theory is usually vague. In contrast, the GAM represents an “adaptive” approach in which the data help guide the choice of appropriate functional form.[10]

In approaching these complexities in building a model, data are partitioned in bits of intervals and regularized estimation through penalized spline is carried out. The complexity of covariates included in models and the estimation method adopted by different authors such as Marano et al., Lang and Brezger[4,11] has, however, caught our attention in this study.

Motivated by this, we studied and compared the sensitivity of models under different censoring percentages, for several sample sizes within the framework of Weibull distributions whose shape and scale parameters are obtained by varying the variances while keeping the mean at 1; through the unit root function implemented in R.

The study further investigates: The performance of single additive models (SAMs) and the modified piece-wise additive extension or piece-wise additive models (PAMs) under various censoring percentages and sample sizes employing three levels Weibull baseline variance.

METHODOLOGY

The risk data used for this paper were simulated from a Weibull baseline hazard distribution which was used to generate survival times for sample sizes of 100, 500, and 1000, respectively. Various censoring levels or percentages of no censoring “0%,” low “about 25%,” moderate “about 50%,” and high “about 75%” were used.

Model Specification

The Cox hazard model

The baseline hazard rate is unspecified and assumes that covariates x=(x₁,…,x_p) act multiplicatively on the hazard rate through the exponential link function.[1]

Model Extension

Where, λ_i (t) is the hazard function for individual i at time t.

γ₁ x_i1…γ_p x_ip: Are time invariant covariates.

γ_j x_in (t): Time-dependent covariates

λ₀ (t) is the baseline hazard function for an individual whose covariates x₁…x_p all have values of 0.

An additive representation of model (1) is given by:

This is a reparameterization of the Cox model

Where, f₀ (t)=logλ₀ (t) is the functional form of the baseline hazard, which implies exp(f₀ (t)). Other aspects of the model include the functions f₁ (t) z₁…f_p (t) z_p which are the functional form of time-varying covariates z₁,…,z_p, the functions f_j (w₁)…f_q (w_q) are possibly non-linear effects of metrical covariates w₁,…,w_q and γ is the usual linear part of the predictor for some categorical covariates.[5,12]

Modification of Piece-wise Additive Models (PAMs)

The proposed model in bits of intervals is given as

with its various terms defined as:

The function f_h=logλ_h is the baseline effect for the h^th interval of PAM

The functions f₁ (t) z₁₁…f_p (t) z_ph are the functional form of time-varying covariates z₁₁,…,z_ph in the hth interval, the functions f_j (w₁₁),…, f_p (w_qh) are possibly nonlinear effects of metrical covariates w_1h,…,w_qh in the h^th interval, and f_spat (s_ih) is a structured spatial effect, where s, s = 1, , S is either a spatial index, with s = s_i if subject i in the h^th bit (interval) is from area s or it is an exact spatial coordinate s =(x_i, y_s), for example, for centroids of regions or if exact locations of individuals are known.

γ is the usual linear part of the predictor.

Where, for each subject i there is a product of H_i terms, H_i being the number of intervals in which the subject is followed. In the expression above, d_ih is the status of the i^th subject within the interval T_h (0 = alive or censored, 1 = failed); D_ih is the time spent in T_hby the subject. From expression (5), it may be seen that L_PAM is proportional to the product of Poisson likelihoods for D_ih with mean parameters:

. As a consequence, the expression of the Poisson regression model is:

Where, h(i) indicate the interval where t_i falls, that is, the interval where individual i died or was censored,

α_h = log(λ_h) are log-hazard parameters, and the term log(D_ih) is an offset.

Bayesian P-splines Approach for Modeling the Unknown Functions

A number of competing approaches are available for modeling and estimating non-linear function f_j of continuous covariates. These include smoothing splines,[8] local polynomials,[13] regression splines with adaptive knot selection,[14,15] and P-splines.[16,17] In this study, Bayesian version of penalized splines is employed following.[11] P-splines were introduced in a frequentist setting by Eilers and Marx, Marx and Eilers.[16,17] The basic idea of P-splines is to assume that an unknown smooth function f_j (x_j) of a continuous covariate x_j can be approximated by a linear combination of B-spline basis functions B_m.[18] That is, denoting the m-th basis function by B_jm, we then obtain

The basic functions B_m are B-splines of degree l defined over a grid of equally spaced knots

, where s is the number of the equally spaced knots.

The Bayesian P-splines method is based on a hierarchical model with non-informative priors for the regression coefficients (b) and a Gaussian random walk (RW) prior of order d for the coefficients of the hazard function (B-spline), conditional to a smoothing parameter τ² the general expression of the RW prior as suggested by Lang and Brezger, Kooperberg and Intrator[11,19] is the following:

The penalty matrix K_j is of the form K_j = D′D, where D is a first or second order difference matrix.

For an independent and identical random effect, the penalty matrix is the identity matrix, that is, K_j=I. The variance parameter

controls the tradeoff between flexibility and smoothing and an inverse gamma prior (the conjugate prior) is assumed, that is,

.

Gaussian Random Field (GRF) Priors

For georeferenced data, it is commonly assumed that v_i = v(s_i) arises from a Gaussian random field (GRF) {v(s),sϵS} such that v = (v₁,…,v_m) follows a multivariate Gaussian distribution as

, where τ² measures the amount of spatial variation across locations and the (i, j) element of R is modeled as R[i, j] = ρ(s_i, s_j). Here ρ(.,.) is a correlation function controlling the spatial dependence of v(s). In “survregbayes” package in R, the powered exponential correlation function

is used, where φ > 0 is a range parameter controlling the spatial decay over distance, vϵ (0,2] is a pre-specified shape parameter, and ‖s–s’‖ refers to the distance (e.g., Euclidean, great-circle) between s and s’. Therefore, the prior GRF (τ²,∅) is defined as

i = 1,…, m where P_ij is the (I, j) element of R^–1[20]

Hazard Model with Regularized Functions

Where, B₀ is the vector form of a B-spline of degree 0 defined for the follow-up period, and π^β and π^τ2 are generic prior densities for the regression coefficients and the smoothing parameter.[20,21]

Regularized piece-wise additive model (RPAM)

The time-dependent effects for each covariate are:

; j = 1,…,p. Thus, for each Z_j, its values are multiplied by a piece-wise constant function:

; in the parameters.

. This enables the effect of each Z_1j to vary in each interval T_h of the original partition of the follow-up:

for t∈T_h. On the other hand, the effects of w₁,…,w_q have been represented through the general expression of B-splines, with the vector

including the spline basis of w_i calculated for the i^th individual. Thus, in this case, the degree and the knots of each spline are fixed following conventional rules. π^β and π_τ2 are generic prior densities for the regression coefficients and the smoothing parameter.

Model Specification to Advance Simulation

Model 1:

Model 2:

Where, λ_PI is the hazard function when partitioning is ignored (PI) or SAM

Where, λ_PD is the hazard function when partitioning is done (PD) or piece wise additive model (PAM).

Test for Non-proportionality

To test the hypothesis that the proportional hazard assumption is valid, the following statement of hypothesis is made.

H₀: δ₁=δ₂=⋯=δ_p (Assumption is valid)

H₁: At least one of the δ_i’ s is not equal to zero (Assumption violated)

Decision rule: Reject H₀ if P ≤ α (level of significance)

Residual measures are used to investigate the departure from the proportional hazard assumption. Schoenfeld residuals are used to test the assumption of proportionality. Schoenfeld residuals are usually calculated at every failure of time under the proportional hazard assumption and usually not defined for censored observations. The overall significance test is called the global test of the model, sighted in Adeniyi and Akinrefon.[22]

Data Analysis

The data are simulated using the functional form of time-varying covariate by Bender et al.[23] given as

The functional form of the continuous covariates as in Brezger[24] is given as:

where x_i~U(–3,3).

For spatial frailty, we propose, S = pnorm(v) and v~mvrnorm(1,S); if S=pnorm(v) then S~mvrnorm, where, S is the covariance matrix for spatial correlation.

Coordinates for spatial correlations follow the uniform distribution. s₁ = runif(N,0,40) and s₂ = runif(N,0,100),[25] obtained the shape and scale parameters of the Weibull distribution from the following

And

for a convenience choice of mean 1 and variance 0.5. Using the uniroot function in R. Parameters were given to be approximately α=1.435523 and η=1.101321. We considered studying the impact of increasing and decreasing the variance of the Weibull distribution while keeping the mean at 1. The result is displayed in Table 1.

Table 1: Shape and scale parameters of the Weibull distributions

Data simulation and analysis were carried out in R using the coda package for spBayesSurv, version 3.6.2. Comparisons were done using deviance information criterion (DIC) (smaller is better) which places emphasis on the relative quality of model fitting and log pseudo-marginal likelihood (LPML) (larger is better) focuses on the predictive performance. Both criteria are readily computed from the MCMC output.