Penalized spline smoothing in multivariable survival models with varying coefficients
Introduction
Modeling of survival data is largely dominated by the proportional hazard (PH) model introduced by Cox (1972). Even though the PH model appeals by simple numerical fitting based on the partial likelihood, the PH assumption often restricts the model in applications since it means that covariate effects remain constant over survival time. This assumption has been under major investigation and numerous papers suggest extensions and testing procedures, see for instance O’Sullivan (1988), O’Quigley and Pessione (1989) , Hastie and Tibshirani (1990), Gray (1994), Hess (1994), Abrahamowicz et al. (1996). For a general overview of estimation and tests in proportional hazard models, we also refer to Lin and Wei (1991), Sasieni (1999) or Grambsch and Therneau (2000). Allowing covariate effects to be dynamic in time leads to a varying coefficient model as generally introduced by Hastie and Tibshirani (1993). Here, constant covariate effects are replaced by smooth but unknown functions. Smooth estimation can then be carried out using e.g. spline fitting, as in Hastie and Tibshirani (1993), see also Kooperberg et al. (1995) or by applying local techniques, see e.g. Fan et al. (1997) or Cai and Sun (2003).
Smooth estimation in survival models is usually based on the partial likelihood function. There are, however, two points of criticism which should be raised against the use of the partial likelihood in the context of smoothing. First, in the simple case that covariate effects are in fact constant over time, that is if the PH assumption holds, the cumulative (integrated) hazard function in the likelihood function factorizes to the cumulative baseline hazard multiplied by the covariate effects. If the baseline hazard is then estimated by the empirical survivor function, the resulting profile likelihood for the parameters is equivalent to the partial likelihood suggested by Cox. This justification of the partial likelihood is due to Breslow (1972) (see also Cox, 1975 or Wong, 1986). However, if covariate effects do vary with time, that is if the PH assumption is violated, such factorization of the cumulative hazard does not exist and, consequently, the partial likelihood does not have any justification as profile likelihood function. Secondly, in partial likelihood estimation the baseline hazard is treated as nuisance component and not explicitly estimated. In applications, however, knowledge about the baseline hazard can be of interest, in particular if smooth, non-parametric regression is pursued. For this reason, it seems worthwhile to work directly with the likelihood function. This approach is pursued in this paper in order to fit a smooth, non-proportional hazard model. The integrated hazard function in the likelihood is thereby approximated using numerical integration based on a trapezoid approximation. This in turn leads to a simple likelihood function which resembles a Poisson model.
As smoothing technique we employ penalized spline fitting (P-spline). The approach was originally introduced by O’Sullivan (1986), but the procedure finally achieved general recognition with the paper by Eilers and Marx (1996). A comprehensive overview about the current state of the art is found in Ruppert et al. (2003). P-spline smoothing in survival models has been studied in Cai et al. (2002) for baseline hazard smoothing. The underlying idea of P-spline smoothing is to fit a smooth curve by using a high-dimensional basis. But instead of simple parametric fitting a penalized version is pursued to provide a smooth fit. The approach resembles standard spline smoothing as discussed, e.g. in Wahba (1978), or in its generalized form in Green and Silverman (1994). The major difference is that for spline smoothing the dimension of the corresponding spline basis grows with the sample size. In contrast, for P-spline smoothing a finite-dimensional basis is used, where the dimension is chosen in a rich and generous manner. The approach is numerically very handy. It also has strong links to linear mixed models (see Wand, 2003) and to penalized quasi-likelihood (PQL) estimation in generalized linear mixed models (GLMM), as discussed in Breslow and Clayton (1993). The connection becomes obvious if the penalty is rewritten as a priori distribution on the coefficients of the basis. In fact, the smoothing parameter steering the amount of penalization is then playing the role of the a priori variance in the resulting GLMM. We utilize this link for smoothing parameter estimation. It will be demonstrated that the PQL approach is numerically simple but fails to estimate reasonable smoothing parameters in low-intensity hazard models. Alternatively, an EM-based procedure as suggested in Booth and Hobert (1999) could be used for the price of increased numerical effort. We suggest a hybrid approach based on the numerically attractive PQL estimates combined with an Akaike criterion.
Our data example concerns the success (survival) of newly founded companies using the database of the “Munich Founder Study”. During 1985 and 1986 there were business registrations at the local chambers of commerce in Munich and surrounding administrative districts. About 5 years later in 1990, a sample of the businesses was drawn with subsequent interviews of about founders. Details on the study can be found in Brüderl et al. (1992). The recorded variables include the duration of business as response variable (actual duration time if business went bankrupt, censored observation if business is still in the market). Moreover, risk factors like for example whether the company was started with the intention of providing the main income for the founder (yes/no) or whether the business was planned in advance for more than 6 months (yes/no) were collected. More details are given later in the paper. The focus of interest is how these initial risk factors influence the chances of success (survival) of a company and more importantly how this effect vary with the time being in the market. In particular, we focus on young entrepreneurs aged 30 years and younger, leading to a subsample of 369 firms. Corresponding Kaplan–Meier estimates are shown in Fig. 1. The left curve shows the empirical survivor curve, dependent on the effect of the purpose of the business (main income or not) while the right plot shows the effect of planning on the probability of success. There is clear non-proportionality visible which will be taken into account in our modeling exercise.
The paper is organized as follows. In Section 2, we first motivate the use of P-splines for fitting non-proportional hazard models. We demonstrate how integrals of the hazard function can be approximated by trapezoid integration, yielding a Poisson-type model. We provide some asymptotic consideration and discuss practical adjustments of the fitting algorithm. In Section 3, we derive the link to GLMMs and discuss the estimation of the smoothing parameter. The data example and simulations are found in Section 4. A discussion finalizes the paper. Technical details are provided in the appendix.
Section snippets
P-spline fitting
Let denote the survival time of the ith individual or observational units and let be the corresponding right censored time, . We observe and define the censoring indicator if and otherwise. With we denote the p-dimensional covariate vector for the ith individual, which for simplicity of presentation is assumed to be time constant. The hazard function is then modeled aswith as baseline hazard and as vector of
PQL estimation
Penalized spline smoothing has strong affinities to PQL estimation in GLMM as discussed in Breslow and Clayton (1993) (see also McCulloch and Searle, 2001). For normal response models, this link is illuminated in depth in Wand (2003) (see also Ruppert et al., 2003). For non-normal response models we achieve the link in the following way. We consider coefficients , , as independent normally distributed variables withwhere is the (generalized) inverse of . The
Simulation
We simulate survival data for individuals on a discrete time grid using a constant drop out probability of 3% for each time interval t to . The two binary covariates and are randomly chosen with and . As dynamic effects we include as constant baseline hazard and and . In Fig. 3, we show for one simulation the principle of the hybrid smoothing parameter selection. Smoothing parameter estimation is started
Discussion
We demonstrated the use of P-splines for fitting non-proportional hazard models. Numerical integration was pursued which led to Poisson data. Multivariate smoothing parameter selection was carried out by a hybrid procedure, utilizing the link between P-spline smoothing and GLMM. In particular, complicated grid searching was avoided and the routine is numerically simple. A data example demonstrated the new insight which could be gained by allowing hazard functions to be dynamic in time.
References (37)
- et al.
Time-dependent hazard ratiomodeling and hypothesis testing with application in lupus nephritis
J. Amer. Statist. Assoc.
(1996) - et al.
Asymptotic Techniques for use in Statistics
(1989) A Practical Guide to Splines
(1978)- et al.
Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm
J. Roy. Statist. Soc. Ser. B
(1999) Comment on “regression and life tables” by D. R. Cox
J. Roy. Statist. Soc. Ser. B
(1972)- et al.
Approximate inference in generalized linear mixed model
J. Amer. Statist. Assoc.
(1993) - et al.
Bias correction in generalised linear mixed models with a single component of dispersion
Biometrika
(1995) - et al.
Survival chances of newly founded business organizations
Amer. Sociol. Rev.
(1992) - et al.
Mixed model-based hazard estimation
J. Comput. Graphical Statist.
(2002) - et al.
Local linear estimation for time-dependent coefficients in Cox's regression models
Scand. J. Statist.
(2003)
Local roughness penalities for regression splines
Comput. Statist.
Partial likelihood
Biometrika
Regression models and life tables (with discussion)
J. Roy. Statist. Soc. Ser. B
Analysis of Survival Data
Flexible smoothing with B-splines and penalties
Statist. Sci.
Local likelihood and local partial likelihood in hazard regression
Ann. Statist.
Proportional hazards tests and diagnostics based on weighted residuals (corr: 95v82 p668)
Biometrika
Modeling Survival Data
Cited by (53)
Association of Bioprosthetic Aortic Valve Leaflet Calcification on Hemodynamic and Clinical Outcomes
2020, Journal of the American College of CardiologyCitation Excerpt :We built a series of nested bivariate logistic models, and multiple testing was adjusted by Bonferroni correction (11). In survival analyses, the pattern of the association between AVCd and the clinical endpoint was initially examined with AVCd modeled as a penalized spline (16). We then modeled AVCd as both a continuous and a categorical variable according to the optimal threshold determined by maximally selected rank statistics (17) as well as quartile of the AVCd distribution.
On the maximum penalized likelihood approach for proportional hazard models with right censored survival data
2014, Computational Statistics and Data AnalysisA unifying switching regime regression framework with applications in health economics
2024, Econometric ReviewsAssessing nonlinearities and heterogeneity in debt sustainability analysis: a panel spline approach
2023, Empirical Economics