Faculty of Economics and Business - Amsterdam School of Economics
Master thesis in Econometrics
Nonparametric Estimation for Cancer
Survival Rates
Manuela Puente Beccar
Student number: 10604790
Supervisor
:
K. J. van Garderen
Second reader
:
J. C. M. van Ophem
Intership made at
:
Pharmerit International
Company supervisor :
I. Majer
1
Index
1 Introduction ... 3 2 Literature review ... 4 3 Data description ... 7 4 Theoretical background ... 8 4.1 Survival Analysis ... 84.1.1 Discrete time: non-parametric approach ... 9
4.1.2 Continuous time: parametric approach ... 10
4.2 Smoothing Techniques ... 11
4.2.1 Kernel regression ... 11
4.2.2 Local polynomial regression ... 12
4.2.3 Splines ... 12
4.3 Confidence intervals ... 13
4.3.1 Standard confidence interval ... 14
4.3.2 Bootstrap ... 14
4.3.3 Profile Likelihood ... 14
5 Estimation results ... 14
5.1 Smoothing methods ... 14
5.2 Confidence Intervals ... 19
5.3 Analysis of the three methods ... 20
6 Parametric method ... 21
6.1 A mixture of distributions ... 21
7 Simulation of survival data ... 24
7.1 Appropriateness of the simulations ... 24
2
8.1 Closeness of the smoothed curves ... 26
8.2 Width of the confidence interval ... 26
8.3 Robustness check ... 27
8.4 Bias concern ... 29
9 Selection of the smoothing technique ... 32
10 Comparison to the General Population death rates ... 32
10.1 Comparison of the hazard rate for the real data set ... 32
10.2 Simulation procedure for the confidence intervals ... 33
10.3 Comparison of the hazard rate for simulated data sets: Approach 1 ... 34
10.4 Comparison of the hazard rate for simulated data sets: Approach 2 ... 35
11 An economic implication ... 37
12 Conclusions ... 38
3
1 Introduction
The development of the medical science increases the survival probabilities of cancer patients year upon year. For patients diagnosed with invasive melanoma the death rates after some years may even be comparable to those of the general population. This is what we investigate in this thesis. The research question is whether melanoma patients have similar death rates as the general population that has not been ill, or if they still face a statistically significant higher risk.
To answer this question we need to analyze the hazard rate of melanoma patients which we can obtain with survival models. Traditional distribution functions in survival analysis fail to replicate the survival function of melanoma patients by underestimating it. This fact has been acknowledged by many authors and parametric models for a group of “cured” patients have been developed in order to overcome this problem. The most popular approach is a mixture model of Weibull distributions (Farewell (1982), Rosenberg (1995), Peng & Dear (2000), Rutherford et al. (2013)). The idea behind the mixture model is the existence of two different groups of patients, a high risk group whose proportion is reduced in time and a low risk group.
However, as in any parametric approach, if the data does not follow the assumed distribution, inference based on these models can be erroneous. Besides, it has been argued that supporting the distributional assumptions is difficult from a biological perspective as will be explained later on Farewell (1986).
Pharmerit International, a research organization in health economics that was interested in the research question asked me to answer it and advised me to compare the hazard rate of melanoma patients obtained from a Kaplan-Meier survival analysis to the tabulated death rates of the general population in The Netherlands as part of my internship program. They proposed to use different smoothing techniques such as the traditional Kernel method, cubic splines, P-splines and maybe another method. I began my work already knowing the characteristics of the data and knowing that the covariates present in the data were not very informative to estimate the survival probabilities since given the absence of a control group, the patients that received more treatment probably needed it the most.
For these reasons the present work is based on approximations of the cumulative hazard and the hazard rate obtained from the non-parametric survival function estimated with the Kaplan-Meier method due to its widespread use in survival analysis. The smoothing techniques are needed given the high variability of the hazard functions: they enable their analysis and comparison. In this work I make use of three different techniques, Kernel regression, local polynomial regression and penalized splines, in order to find the best curve to describe death rates of melanoma patients. I did not use P-splines for reasons explained later on.
There is abundant literature devoted to the specific properties and performance of the individual methods (Tanner & Wong (1983), Lambert et al. (2011), Rutherford et al. (2013)) but not much information on the performance of one versus the others. In addition, most studies use Proportional Hazard models or Accelerated Failure Time models, while in my case the focus will be
4
primarily on the shape of the hazard function and ask if they are essentially different. Using for instance a common baseline hazard is therefore fundamentally inappropriate.
I compare the different smoothing techniques through simulations from a parametric mixture model in a Monte Carlo analysis. If the assumed data generating process is valid then the different techniques are comparable in many simulated samples. After choosing the best smoothing technique I compare the lower bound of the confidence interval of the estimated hazard rate with the death rates of the general population.
2 Literature review
A common approach for survival analysis is a parametric specification of the density function of survival times and from there also specifications of the survival function and the hazard rate. However, as mentioned before the standard distributions fail to replicate the survival probabilities of cancer patients. Farewell (1982) proposed a mixture of distributions arguing that in the presence of individuals that respond to some treatment while others do not, the standard parametric distributions are not able to capture the characteristics of the survival function since analyzing the follow-up time for the long term survivals as censored observations of a single probability distribution is not appropriate. Instead he proposed the estimation of a mixture of distributions each for a different group of patients. He worked with a mixture of Weibull distributions and estimated the parameters with maximum likelihood. He argued that in some cases the interpretation of the parameters may be appealing from both a biological and a statistical point of view, however he advises that the use of such models should be restricted to problems in which there is strong evidence for the possibility that the individuals come from two separate populations, because inference from those models will be made about the two populations whether they exist or not.
An alternative is to use semiparametric and nonparametric models. There is abundant literature in the use of these models in survival analysis and specifically in studies with cancer patients, however a large majority is based on the proportional Cox model since most of the literature is interested on estimating the effect of covariates when there is a treatment and a control group, so a relative measure is all that is needed. This model is semiparametric because it does not require the specification of a baseline hazard: by comparing two groups all that matters is the relation between the hazard of one group with respect to the other and that will be a function of the covariates only. For example Cai and Sun (2003) focus on the estimation of the Cox proportional hazards model and make use of local linear regression in order to estimate the coefficient functions ( ) in the model
( | ) ( ) (∑ ( )
)
where ( ) is the baseline hazard and represents the covariates. They present the asymptotic consistency and normality of the local linear estimation as well as a consistent variance estimator
5
and test their estimation method with two simulated samples and one real data set from the study of gastric cancer patients.
A smaller number of articles are interested in the shape of the hazard rate and the authors make use of different smoothing techniques as well as of different ways to estimate the hazard rate. Tanner and Wong (1983) make use of Kernel smoothing and investigate the properties of the failure rate from:
̂( ) ∑ ( ) ( ( ))
where K is the symmetric nonnegative Kernel function and is a dummy for censored observations. The aim of this paper is to derive small and large sample expressions for the mean and the variance of the hazard estimator along with other properties in the case of censored data. Articles that are more related to my investigation are for example those by Rosenberg, Royston and Parmar and Rutherford et al. so I will present a more detailed summary of their work.
Rosenberg (1995) models the hazard rate and the survivorship function from censored survival data with cubic B-splines. B-splines have minimum support among all sets of basis functions for the space of cubic splines and thereby ease computations. He finds better results in terms of mean squared error compared to a non-parametric estimate proposed by Whittemore and Keller (1986). The non-parametric estimator was often less biased than the spline procedure for both the hazard and survivor functions, thus, the performance gain of the spline was frequently achieved by trading off greater bias for reduced variability. However, the spline procedure was both less biased and less variable in the right tail, which would indicate that for my purpose this procedure is more appropriate.
Rosenberg makes use of three methods for the calculation of confidence intervals: the delta method, profile likelihood and bootstrap. The three methodologies gave similar results, but the author advices the use of the bootstrap for three reasons: it incorporates nonnegativity constraints; in case of small to moderate samples, it reflects possible skewness and one set of bootstrap replications can be used to calculate confidence intervals simultaneously for all functionals of interest, whereas a different profile must be calculated for each functional of interest in order to estimate confidence intervals based on the profile likelihood method. Nonetheless, the error bars calculated by the profile method are nearly as informative as the bootstrap but require 16% as much computation time.
Royston and Parmar (2002) begin by acknowledging that the hazard rate is highly erratic but it is of potential medical interest. To estimate it informatively (that is smoothly) either a smoothing technique or a parametric model may be appropriate. Another motivation to find an estimation of the hazard rate is dealing with non-proportional hazards, which disables the semiparametric approach of the Cox model. Their aim is to develop an approach which is flexible, clinically interpretable, robust, transparent and reasonably easily implemented in standard software and
6
they propose to do so by modelling the logarithm of the baseline hazard function as a cubic spline function of log time.
In their article, they point out two disadvantages of cubic splines, one is that they are constrained to be linear beyond certain extreme observations and the second one is that they are not globally monotone. However, the linearity constraint imposes monotonicity in the tail regions where the observed data are sparse, whereas in regions where data are dense, monotonicity is effectively imposed by the data themselves. As a result, the use of cubic splines appears to guarantee monotonicity in data sets of any reasonable size, which they consider to be at least 50 uncensored observations.
Regarding the placement of the internal knots they explain that an optimized knot selection does not appear critical for a good fit and may even be undesirable in that the fitted curve may follow small-scale features of the data too closely.
In the paper written by Rutherford et al. (2013) interest lies on approximating complex shapes of the hazard rate for the purpose of obtaining absolute measures of risk from time-to-event data. They state that even in the estimation of a Cox model finding a fit for the baseline hazard has proved to be advantageous in being able to estimate absolute measures of risk and that the standard parametric models such as exponential, Weibull or gamma are not flexible enough to capture multiple turning points that might be present in some complex hazard functions.
They take the logarithm of the Weibull survival function to work in the log-cumulative hazards scale:
[ ( )] ( ) ( )
and model ( ) using cubic-spline functions varying its complexity with the number of knots which also define the degrees of freedom (which vary from 1 to 10 in their estimations) and compare the results also with a mixture of Weibull distributions.
They perform 1.000 simulations from the mixture model with different scenarios by simulating the true hazard from varying parameters for a mixture of Weibull distributions, and they take sample sizes of 300, 3.000 and 30.000. In order to compare the estimations they measure the absolute area difference between the fitted and true functions using numerical integration.
They find that the mixture Weibull approach performs similarly to a flexible parametric model with 3 degrees of freedom in terms of model fit, highlighting that it is not as flexible as the cubic-splines models with more degrees of freedom.
Another important result is that the use of restricted cubic splines does not rely heavily on the ‘correct’ selection of the number of knots. With sufficient number of knots, the spline functions are flexible enough to capture complex hazard functions closely. Using too few degrees of freedom leads to a poor fit, however, care must be taken not to specify too many degrees of freedom as over-fitting can become a problem, particularly for modestly sized datasets.
7
Another important article is that written by Eilers and Marx (1996). They propose a difference penalty on coefficients of adjacent B-splines instead of the usual approach to penalize the curvature of the function, and they call this new approach P-splines. As I mentioned in the introduction I was asked to work with this approach. In the article the authors assert that the computations are relatively inexpensive, easily incorporated into standard software, P-splines conserve moments of the data (means, variances) and have polynomial curve fits as limits.
Several comments of other authors included in the same paper question whether P-splines are really better. As an example, Chong Gu stated that P-splines can be as useful as other variants of penalized regression splines but not really advantageous. According to him the mechanical handling of the difference penalty is very interesting computationally, but that does not concern the end users. M. C. Jones stated that everyone can understand the usual derivative penalty, however, with P-splines the intuition is unfortunately taken away: he questions “Which is more interpretable: a roughness penalty on a curve or on a series of coefficients?” This is the reason why I decided to disregard this approach in my investigation.
So far we can see that there are many different approaches to estimate the hazard rate. The investigation I perform is based on the hazard rate derived from the Kaplan-Meier estimation of the survival function because this is the standard procedure in survival analysis and it is thus easily understood by anyone in need of a smooth version of the hazard rate. From the literature review I also rescue the following:
Cubic splines seem an appropriate approach, but there is no basis to support these ones over Kernel smoothing or local linear regression. The placement of the knots is not very important.
Mixture models are not flexible enough to capture characteristics of the hazard rate.
The article by Rosenberg supports the bootstrap procedure to estimate confidence intervals but taking into account that the profile likelihood method may also be interesting due to reduced computation time.
In order to simulate survival times it seems that the most appropriate approach is a mixture of Weibull distributions.
The use of P-splines has no added value compared to the standard cubic splines.
3 Data description
The data comes from the Nederlandse Kankerregistratie and it comprises survival data on 783 melanoma patients (stage IIIC and IV), in the Netherlands, diagnosed in the period 2003-2011. This means that the data presents right censoring refered as Type I censoring in the literature.
The death rates of the general population come from the Human Mortality Database. They are calculated as a ratio of the death count for a given age-time interval divided by an estimate of the exposure-to-risk in the same interval. We used the last data set available (2009).
8
4 Theoretical background
4.1 Survival Analysis
Survival analysis, most commonly known as duration analysis in econometrics, is the statistical analysis of data where the variable of interest is time until an event occurs. In biostatistics the event can be recovery, relapse, and more usually, death.
Most of the times, the data that needs to be analyzed is censored so we have some information about individual survival time but we do not know the exact survival time. If we are interested in survival time until death and we observe a group of individuals, then the survival time of all the individuals alive at the point in time in which our observations stop will be right censored. This is known as generalized Type I censoring: individuals enter the study at different times and the terminal point of the study is predetermined. This type of censoring is considered independent from the survival process but must be taken into account in the analysis.
One way of analyzing censored survival data is to assume two random variables: the survival time and the censoring time and we can only observe the minimum of the two. This way, the censoring dummy will be:
{
This way the analysis is performed with the paired data ( ) and by the assumption of independence, the distribution of does not need to be modelled since the analysis conditional on will not differ. So from now on, when we refer to survival time, we will be referring to .
Two important concepts in survival analysis are the survival and the hazard functions, directly related to the density function of the survival times. The survival function gives the probability that survival time equals or exceeds a moment in time, . All survival functions are non-increasing, they start at 1 and as time tends to infinity they equal 0.
( ) [ ]
( ) [ ] ( ) where ( ) is the cumulative distribution function of the survival times.
The hazard function is the instantaneous probability of the occurrence of the event conditional on survival up to time t. The hazard rate is always non-negative and need not be bounded. For continuous duration it is defined as
( )
[ ]
9
( ) ( ( ))
( ) ∫ ( ) ( ( )) where ( ) is the cumulative hazard function.
4.1.1 Discrete time: non-parametric approach
In discrete time, for we have:
( ) [ ] ( ) ( )
( )
( ) [ ] ( ) ( ) ∑
The most common approach for a non-parametric estimation of survival functions is the Kaplan-Meier method. To construct the Kaplan-Kaplan-Meier estimator we order the data in the following way: Let be the observed discrete failure times and the number of events (deaths) in the interval . We also need the number of censored observations in those intervals , and , the number of individuals at risk (alive) up to time .
The Kaplan-Meier estimator of the survival function will be given by: ( ) ( )
Note that decreases not only due to the values but also due to , so in this way our estimated already takes into account the censored observations.
The variance for the survival function from which standard pointwise confidence intervals can be derived is given by Greenwood (1926):
̂ [ ̂( )] ̂( ) ∑
( )
The Kaplan–Meier estimator can be shown to be the nonparametric Maximum Likelihood Estimator (Kalbfleisch & Prentice (1980)).
In the Kaplan-Meier estimation there is a problem of bias if the largest time point is censored, since then the value of ( ) beyond this point is undetermined because we do not know when this last survivor would have died if the survivor had not been censored. Efron (1967) and Gill (1980)
10
suggests estimating ( ) by 0 or ( ) respectively beyond the largest study time. Both estimators have the same large-sample properties and converge to the true survival function for large samples.
Another discrete time approach is known as the Nelson-Aalen estimator, which focuses on estimating the cumulative hazard instead of the survival function as follows:
( ) ∑
The Nelson-Aalen estimator and the Kaplan-Meier estimator are asymptotically equivalent.
4.1.2 Continuous time: parametric approach
Commonly used distributions used in survival analysis include the exponential, the Weibull and the log-normal.
With the Weibull distribution we can have monotonically increasing, decreasing or constant hazard functions (when the Weibull reduces to the exponential) while the lognormal distribution also allows a hazard rate that first increases and then decreases. The survival functions are the following:
Exponential : ( )
Weibull : ( )
Log-normal : ( )
The estimation of the parameters is usually achieved by maximum likelihood where the likelihood contribution of each observation will vary if the observation is censored. The contribution to the likelihood of actual failure observations is ( ) while for right-censored observations we know only that the individual survived, so the contribution is ( ), the survival function. The log-likelihood function is then:
( ) ∑
( ) ( ) ( )
where is a variable equal to 1 if the observation is not censored, the censoring dummy previously introduced. Given the properties stated above regarding the relations between the survival function, the density function and the hazard rate, the log-likelihood can also be written as:
( ) ∑
( ) ( )
The standard parametric functions are not always able to describe the clinical reality which is why also mixtures of these distributions are used to find a better parametric fit.
11
4.2 Smoothing Techniques
There are many different smoothing techniques and non-parametric estimation methods. Most of them, including Kernel regression and local polynomial estimation, are local averaging methods which consider only the neighborhood of each smoothing point, while others, like splines, can be seen as a global approach.
4.2.1 Kernel regression
The Kernel regression method is an averaging method where declining weights are assigned to observations more distant to the estimation point. The weights are determined by Kernel functions that satisfy some restrictions. For instance, they must be continuous and integrate to unity. The Nadaraya-Watson Kernel regression estimates the regression function ( ) in the model ( ) where ( ) is the expectation of and is an error term according to the following approximation:
̂( ) ∑ (
)
∑ ( )
where here is the bandwidth to consider the different close to and ( ) is the Kernel function. The most common Kernel functions are the normal and the Epanechnikov:
Normal : ( ) ( )
Epanechnikov : ( ) (| | )
where 1 is 0 when | | , so the Epanechnikov Kernel is restricted to [ ] while the Normal Kernel is unrestricted. Both Kernel functions give progressively heavier weights to points close to . The estimation results are very similar for any choice of Kernel function, but this one also affects the bandwidth and different bandwidths do have an important effect on the final smoothed curve.
A small bandwidth prevents bias, but taking into account a small number of observations increases the variance. To determine the optimal bandwidth the approach is to find the one that minimizes the mean of squared errors of the density function (estimated by Kernel smoothing of the histogram of ( )) at .
[( ̂( ) ( )) ]
So, for a measure of performance at every value of , we must minimize the Integrated Mean Square Error
12
If we assume that the density functions is Gaussian, this method is known as the Plug-in methodology for finding the optimal bandwidth which is less computationally burdensome than other methods and it usually provides a good starting point Cameron & Trivedi (2005).
Since the optimal bandwidth varies with the Kernel, the choice of the Kernel function will be important because of its influence on the bandwidth choice.
Another important consideration in the use of Kernel smoothing is that at the beginning and ending of the points there will be less or zero observations at one side of the smoothing point. To take this into account, the use of asymmetric Kernel functions has been proposed (Gasser & Müller (1979)) where the weights assigned to observations at one side of the smoothing point differ from those assigned to the other side.
4.2.2 Local polynomial regression
The Kernel regression estimator assumes that ( ) is a constant. If we allow the estimation of ( ) to be a polynomial then we are dealing with local polynomial regression. If the degree of the polynomial is one we have local linear regression.
To visualize this approach we can approximate ( ) as follows with Taylor series: ( ) ( ) ( )( ) ( )( ) ( )( )
( )
( ) ( ) ( ) ( ) ( )
So we must find the polynomial of degree that best approximates ( ) in a neighborhood of , and then assign less weight to approximations on points further away from usually with Kernel functions, i.e. we choose the that minimize the following function:
∑ ( ) ( )
( )
With local polynomial regression the determination of the bandwidth and the Kernel function is equally important and the same method as before can be applied.
Here we can see that if we choose , the solution to the minimization problem is the Nadaraya-Watson estimator and it is easy to grasp that even if the Kernel regression approach was appealing, the approximation of ( ) by ( ) disregarding the rest of the terms in the Taylor expansion can be quite poor for some non-linear functions.
4.2.3 Splines
With the previous methods the irregularity of the smoothed functions was controlled by the bandwidth of the Kernel functions. Another approach is to find a function that fits the data penalizing the irregularity:
13
( ) ∑( ( )) ∫( ( ))
The first term in the right hand side is the MSE previously explained. The second term is a restriction on the second derivative of the function, i.e. the curvature. This means that for different values of we can determine how smooth we impose our function to be. As we can see, we are willing to let the MSE increase as long as it is compensated with an increase in smoothness. The solutions for ( ) are piecewise cubic polynomials estimated in-between a selected number of knots. These polynomials have continuous first and second derivatives which coincide at the knots.
To find the value of , one approach is again minimize the squared error but in this case the expected prediction error:
( ) [ ( )]
where are new data. An estimator of PSE can be found by the leave-one-out Cross-validation function:
( ) ∑ ( ( ))
where ( ) is the smoothing spline estimator fitted from all data less the observation. However, in the case of repeated values of , the leave-one-out procedure will not function properly. For this reason it is advised to use the Generalized Cross Validation function given by Carew et al. (2003):
( ) ∑ ( ( ̂ ) ( ( ))
where ( ) is the trace of the linear transformation that maps the data to their fitted values, also known as hat matrix.
Cross-validation can also be used to determine the optimal bandwidth in Kernel regression but it is computationally burdensome and does not outperform the plug-in estimate. It is however the most common approach for penalized splines.
4.3 Confidence intervals
For the purpose of this work we are interested in finding the confidence interval for each point of our estimated curves. This is analogous to finding the pointwise confidence bands. Several methods have been developed for the estimation of confidence intervals. Some of them rely heavily on distributional assumptions while others are completely non-parametric.
14
4.3.1 Standard confidence interval
This way of estimating the confidence interval is the standard procedure and it is based on asymptotic theory, since it is supported by the Central Limit Theorem. The confidence interval is estimated by ̂ ̂ √ where ̂ is the estimated parameter, ̂ the estimated standard deviation and the sample size.
4.3.2 Bootstrap
In situations in which the standard error is unknown or the approximation to normality seems too extreme, it is preferable to make use of bootstrap confidence intervals. The bootstrap procedure consists of treating the data set as the entire population. Then, new samples are obtained by sampling from the original data with replacement. By resampling many times we obtain a large number of possible draws from our data. If the number of bootstrap samples is very large, we could treat the set of drawn observations at each evaluation point as all the possible values that our function could take at such a point. Efron (1981) showed that resampling the pairs ( ) is equivalent to resampling from the observed survival times and the unobserved survival times introduced before ( ) .
If, after the resampling, we order the different resampled pairs and take the 5th and 95th percentile, we have an approximation of the values that are inside a confidence level 90% of the times and, as in a one-sided test, which values are above the lower bound 95% of the times. This method for calculating confidence intervals is completely non-parametric since it does not rely on any parametric assumptions.
4.3.3 Profile Likelihood
The 95% profile confidence interval on a given point x is defined by the smallest and largest values of x such that the profile deviance ( ) ( ) exceeds the 0.95 percentile of a chi-square distribution with 1 degree of freedom where is the parameter that maximizes the likelihood function.
5 Estimation results
5.1 Smoothing methods
To begin our non-parametric analysis we first need an estimation of the survival function. To do so we made use of the Kaplan-Meier method, which is the most popular tool in survival analysis. As previously mentioned the Nelson-Aalen estimator is an alternative for the Kaplan-Meier estimator and even if they are asymptotically equivalent, it is stated by Klein and Moescherberg (1950) that the Nelson-Aalen estimator has better small sample properties and Colosimo et al. (2010) find a better performance of the Nelson-Aalen estimator for percentile estimation.
The definition of small sample differs in every case, so we performed the estimation of the survival function by both methods finding dismissible differences. We chose to continue the analysis with the Kaplan-Meier estimator given its widespread use.
15
For the estimation we used the software R1 where many packages are freely available. A Kaplan-Meier estimation was performed with the data by creating a survival function with the survfit function2 of the Survival package.
This is the survival curve with the 95% confidence interval in dotted lines, calculated from the Greenwood formula. We can observe that the curve levels-off at a survival probability close to 20% beyond a survival time of 60 months.
With this Survival function the cumulative hazard was estimated as ( ( )) and the hazard rate was approximated from ( ) .
1
R Core Team (2014). R: A language and environment for statistical computing. R Foundation for statistical Computing, Vienna, Austria.
2
16
The smoothness of the cumulative hazard is similar to the survival function, but the hazard rate is so variable that is hard to picture its real behavior. This is why we need to make use of smoothing techniques.
The first smoothing technique applied was the Nadaraya-Watson estimator with Epanechnikov Kernel and plug-in estimate for the bandwidth choice with the R functions ksmooth and dpik3 of the KernSmooth package. The resulting smoothed functions are:
The optimal bandwidth for the hazard rate provides a smoothed hazard difficult to work with, because it replicates the irregularities of the data.
3 Matt Wand (2014). KernSmooth: Functions for Kernel smoothing for Wand & Jones (1995). R package
17
The inclusion of different weights for observations close to the ending points provided an increasing curve in the right tail. From what is known about the survival characteristics of melanoma patients, this result is not logical, so we dropped the idea of differentiating the smoothing technique for points close to the end of the observations.
The second smoothing technique we applied was local linear regression with Gaussian Kernel and plug-in methodology to estimate the bandwidth with the function locpoly of the same Kernsmooth package.
With this method we observe that the cumulative hazard shows a stable trend and the optimal bandwidth selection gives a much smoother hazard rate than the Kernel method but we still observe a high peak at the end of the hazard rate curve, where a relatively high estimated hazard rate in one point is reflected even with the selected level of smoothness.
The third technique we used was cubic splines with a Generalized Cross-Validation4 procedure to define the smoothing parameter. The command used is smooth.spline, a built-in function in R.
4
The Generalized Cross Validation is used instead of ordinary one to maintain consistency throughout this work since in some smoothings performed later on there are repeated points in the x-axis and the ordinary Cross Validation method is not recommended in such cases.
18
Observing the plots, this methodology seems the most appropriate to find a smooth approximation of the hazard rate and the cumulative hazard specially because our aim is to compare the right tail of the hazard rate and with this methodology we observe a steady behavior of the hazard rate from month 70 onwards.
As we can see in the graphs, the cumulative hazard functions have a less erratic behavior and hence their estimation is more precise than that of the hazard rate. This provides an alternative for the estimation of the hazard rates by taking the first derivative of the estimated cumulative hazard in the case of local polynomial and splines (the first derivative of the linear functions and the cubic functions respectively). This approach is similar to the one taken by Rutherford et
al.(2013) where they make use of cubic splines to estimate the logarithm of the cumulative hazard
19
Performing the estimations in this manner we must consider that in the spline smoothing technique the resulting degrees of freedom corresponding to the optimal in case of the hazard rate is 9.6, while for the cumulative hazard it is 86. If we take the derivative of the smoothed hazard with 86 degrees of freedom, the resulting hazard rate is too rough, so the estimation is performed with the 9.6 degrees of freedom selected for the hazard rate.
In case of the local polynomial the bandwidth does not change, since it is determined by the values in the x-axis.
5.2 Confidence Intervals
The method used in the computation of confidence intervals is the nonparametric bootstrap. We chose this method because it is suitable for both parametric and non-parametric methods and it allows the construction of confidence bands. We estimated point-wise confidence intervals instead of simultaneous coverage confidence bands since we are interested in the acceptable deviation of each point in the curve: we do not need the confidence intervals to cover the corresponding true values simultaneously.
Our interest lies on the lower bound of the confidence interval, since our aim is to compare our estimations with hazard rates of the general population, which are certainly lower than those of cancer patients.
Taking 1.000 bootstrap samples of the survival data, re-estimating the hazard rate (for the Kernel method) and the cumulative hazard (for the other two) for each replication, applying the different smoothing techniques and taking the 5th and the 95th percentile we estimated the confidence interval for the hazard rate obtaining the following results:
20
Initially I also considered the estimation of confidence intervals through different methods, but by estimating the hazard rate as the derivative of the cumulative hazard, other confidence intervals such as the standard procedure based on the estimated variance or from profile likelihood would be for the cumulative hazard.
5.3 Analysis of the three methods
The shape of the estimated confidence intervals provides interesting information for the analysis of the three methods.5
First, in the case of Kernel smoothing the lower bound is actually above the smoothed curve of the estimated hazard rate. This can occur by definition 5% of the times even if the true curve is indeed inside the confidence interval.
Second, the upper bounds are noticeable high, especially the points corresponding to the lower survival times in the Kernel smoothing and in the right tail for Kernel and local polynomial. They are not very relevant since our concern is whether the general population death rates lie between the real curve and the lower bound, but it highlights two characteristics of our data: the high variability in the estimated hazard rate for the low survival times and the presence of a reduced number of events for large survival times. If we think of the way the bootstrap procedure works, we can imagine that if two high values of the hazard rate are randomly chosen while a low value between them is not, we can obtain a smoothed curve that overestimates the hazard rate. In the same manner, for large survival times, if the event that occurs approximately in month 90 is chosen while the previous one is not, we observe a considerable change in the hazard rate that will be even more marked in the case of the local polynomial for the considerable change in the slope of the cumulative hazard.
5
For the Kernel smoothing the bandwidth was doubled to reduce the roughness of hazard curve. For the other two methods the hazard rate presented is the approximated with the derivative of the cumulative hazard.
21
These problems are not present for the spline smoothing technique, since this method does not rely on the presence of observations close to the evaluation points; this characteristic is derived from the fact that splines are a global approximation of the points instead of a local one.
With any of the smoothing techniques we are able to compare the hazard rate of melanoma patients with that of the general population: all of them allow us to have an estimate of the hazard rate at any point in the interval of the observation period since without them we could only compare the existing observations. However, under the previous considerations, it seems that the spline smoothing technique is the most appropriate method.
Our comparison of the smoothing techniques per se is mainly visual. The other considerations were possible by analyzing the bootstrapped data sets. With the estimation of a parametric model we would be able to simulate different samples and hence delve our analysis and validate our conclusions.
6 Parametric method
To simulate survival times with similar characteristics to those of melanoma patients, we need to find a parametric method that fits our data. Using the traditional distribution functions in survival analysis, we find that even the best approximation fails to replicate the survival function by underestimating it.
As we can see, the best approximation is given by the lognormal distribution, but still the fit compared to the non-parametric Kaplan-Meier estimation is poor.
6.1 A mixture of distributions
To find a better parametric fit to our data we applied a mixture of distributions. A mixture of Weibull distributions is a common approach in the literature (Farewell (1982), Rosenberg (1995),
22
Peng & Dear (2000), Rutherford et al. (2013) but we also applied a mixture of exponential distributions in case the estimation of less parameters would be beneficial.
The parameters were estimated by maximizing the following log-likelihood function: ( ) ∑
( ) ( )
where ( ) ( ( )) and ( ) ( ) are derived from a mixture of Weibull distributions, with weight , shape parameter and scale parameter (Rutherford et al. (2013)).
( ) ( ) ( ) ( )
( )
( ) ( ) ( )
( ) ( ) ( )
( )
The only difference for the mixture of exponential distributions is that and equal 1.
To proceed with the maximum likelihood estimation I need to find the initial values for the optimizing iterations. For this reason I first applied the kmeans methodology which clusters the data in to k groups: each observation is assigned to the group with the nearest mean. The function is kmeans, a built-in function in R.
With the initial values I performed the maximum likelihood estimation and found that the mixture of Weibull distributions had a lower AIC than that of the mixture of exponentials (5099.071 vs 5099.98). The following is the parametric fit:
23
Parameters Weibull 1 Weibull 2
̂ 1.10 0.94
̂ 0.058 0.006
̂ 0.76
First of all it is important to note that the parametric fit is inside the confidence interval of the Kaplan-Meier survival curve. This means that it can be reasonably used as an approximation since it is not statistically different. The parametric hazard rate is an even smoother version of the spline-smoothed Kaplan-Meier. It may be not flexible enough to provide hazard rates with multiple turning points (Rutherford et al. (2013)) and for small survival times the parametric hazard predicts lower hazard rates than those actually observed, so even if the estimation of the survival function is pretty accurate, the estimation of the hazard rate is not so much.
Analyzing the estimated parameters and hence our parametric model should help us understand our data. A shape parameter smaller than one indicates that the failure rate decreases over time, while a scale parameter larger than one indicates the opposite. Under the “cure model” framework we can interpret that the resulting survival function arises by the combination of two distinct groups of people. Group one, which comprises 76% of the patients, with a lower survival curve due to a shape parameter larger than 1, consists of patients that get increasingly worse until they die. Group two would represent the patients that get cured since their failure rates decrease. In the case of melanoma patients this analysis is questionable, since stating that a group of patients is cured is pushing the reach of cancer treatments too far. This is the reason why, even if the model is statistically appropriate, the interpretation is not too realistic. As stated by Farewell (1986) the use of mixture models should be restricted to problems in which there is strong scientific evidence for the existence of two populations.
It is for these two reasons (not enough flexibility and a not realistic interpretation) that the parametric model will be used exclusively to simulate survival data with similar characteristics to that of melanoma patients and it will not be an alternative methodology to compare with the smoothing techniques of the nonparametric hazard.
24
7 Simulation of survival data
Now that we have parameters that allow us the replication of our survival data, we can also simulate new data with the same characteristics and even study the behavior of the right tail of the hazard rate had we more or less observations in our data set.
Bender et al. (2005) propose simulating survival times generating numbers from a uniform distribution. They explain that if is a random variable with distribution function , then ( ) follows a uniform distribution and so does , hence:
( ) ( ( ))
And therefore ( ). When using a mixture of distributions it is no longer possible to find a closed form expression for the survival time , so the roots of ( ) are found by numerical methods (Crowther & Lambert (2013)).
After obtaining the simulated survival times with the function uniroot, a built-in function in R, I applied administrative censoring at 10 years of follow up as in our original data. The term administrative censoring refers to the right-censoring that occurs when the study observation period ends.
7.1 Appropriateness of the simulations
First it is important to be sure how appropriate is this parametric model from which to simulate. The following table shows the average minimum and maximum values reached by 10 rounds of 100 simulated samples at the 25th, 50th and 75th percentiles of the survival function.
Real data n=200 n=780 n=2000
Percentile 75 5.1 3.8 - 7.6 4.6 - 6.5 5.0 - 6.4
Percentile 50 13.8 10.1 - 17.9 11.7 - 15.6 12.3 - 14.9
Percentile 25 33.9 23.4 - 66.8 27.5 - 44.1 29.7 - 40.8
Average of minimum and maximum survival times at survival probabilities of 0.25, 0.50 and 0.75 obtained by simulations of sample size 200, 780 and 2.000 along with the corresponding values in the real data set.
The true values of the percentiles are in every case inside the range of values covered by the simulations. In this case our comparison is in the horizontal axis which might seem strange but this is the standard procedure in survival analysis: the comparison is between survival times and not survival probabilities.
Some survival curves are plotted here for the different sample sizes to get an idea of the behavior of the simulated data along with the graph corresponding to the preceding table.
25
n=200 n=780 n=2000
It is clear from the plots that with our model we can simulate survival times that approximate the original survival curve and the only difference we obtain by varying the sample size is the overall variation of the estimated survival function.
8 Analysis from the simulation
In order to further assess the performance of each smoothing technique, making use of the possibility of simulating data sets with our validated parametric model, I created two indicators: one that measures the proximity between the smoothed and unsmoothed points and another one that measures the width of the confidence interval. Finally, a third procedure was developed in order to assess the robustness of the different methods.
26
8.1 Closeness of the smoothed curves
The first indicator measures the distance between the real simulated values in the right tail of the distribution and the smoothed values to have an idea of which smoothing technique is closer to the real curve. It is important to note that this measure is not trivial since the smoothed curves are constructed to provide a smoothed approximation of the estimated hazard rate from the Kaplan-Meier curve and not to replicate it. The indicator for each simulated data set is constructed as the average of (the absolute or squared) distance between the last points of the estimated hazard rate (since the 6th year of survival) and the corresponding smoothed values. The table shows the average over 1.000 simulations.
n=200 n=780 n=2.000
Indicator 1 Abs. value Sq. value Abs. value Sq. value Abs. value Sq. value
Kernel 0.0394 0.0219 0.0257 0.0044 0.0206 0.0017
Local poly. 0.0312 0.0330 0.0199 0.0057 0.0162 0.0021
Splines 0.0331 0.0333 0.0205 0.0057 0.0166 0.0022
Average of the absolute value or squared difference between the last points of the estimated hazard rate from the simulations and the corresponding smoothed curves. For sample sizes of 200, 780 and 2.000 the distance is calculated between the last 10, 25 and 80 points respectively since those are approximately the number of observations present in the real data set.
When we observe the mean of the absolute value of the distance (Abs. value) the local polynomial estimated as the derivative of the smoothed cumulative hazard seems the best technique for each sample size, closely followed by the derivative of the spline-smoothed cumulative hazard6. However, the mean squared error (Sq. value) points to the Kernel smoother as the best technique since it is smaller for each sample size. This is because the mean squared error gives a larger penalty to observations that are more distant, and by estimating the hazard rate from the derivative of the cumulative hazard for both local polynomial and splines these “outliers” do not affect the smoothed hazard curves as much as they affect the Kernel smoothed hazard.
Since our aim is to compare the right tail of the hazard rate with the death rates of the general population and we observe that the hazard rate levels off after approximately six years from diagnoses, taking into account the effect of the large penalty for “outliers” does not seem more appropriate. So from this indicator the conclusion would be that the smoothing technique that resembles more the real data is the local polynomial closely followed by the splines.
As the sample size increases we observe that the average distance is smaller for all the techniques, because with a larger number of observations the estimation is more precise.
8.2 Width of the confidence interval
The second indicator measures the distance between the smoothed curve and the estimated confidence interval in order to determine which method provides narrower confidence intervals. For each simulated data set the indicator is the average of (the absolute or squared) distance
6
The values for the local polynomial smoothing and the splines smoothing applied directly to the hazard rate were larger than 1 for all sample sizes, so these techniques were dismissed and every time we refer to local polynomial or splines we will be referring to the derivative of the smoothed cumulative hazard.
27
between the smoothed values and the lower bound. The overall indicator is again an average over 1.000 simulations.
n=200 n=780 n=2.000
Indicator 2 Abs. value Sq. value (10-2) Abs. value Sq. Value (10-2) Abs. value Sq. Value (10-2)
Kernel 0.0217 0.4302 0.0165 0.0853 0.0113 0.0210
Local poly. 0.0102 0.5359 0.0019 0.0004 0.0013 0.0002
Splines 0.0043 0.0022 0.0018 0.0003 0.0012 0.0001
Average of the difference in absolute value or squared between the smoothed curve and the bootstrap estimated lower bound at those points for survival times that exceed 70 months.
What we observe is that, as we would expect, with increasing sample size the distance to the confidence interval diminishes. However, for any sample size, the lower bound of the estimated confidence bands is closer to the smoothed hazard function when we use the splines technique. With a large sample size the local polynomial technique provides also a quite narrow confidence interval.
8.3 Robustness check
As a robustness check I decided to perform a little perturbation in the data to check for the magnitude of variation of the smoothed curves. I added two observations corresponding to deaths in two consecutive days at a survival time of 110 months. The idea behind this procedure is that the addition of two observations should not alter the conclusions of any analysis.
I added two uncensored survival times of 110 and 110.03 months. The plots of the estimated hazard rate and the smoothing methods of one random draw are the following:
28
From the graphs it is evident that the local polynomial and the splines are less affected by the perturbation. This is because they are approximated from the smoothed cumulative hazard which is not as affected as the hazard rate by individual events.
These are random draws of the sample size of 200 and 2.000:
n=200
29
From the plots we can imagine that the effect of two observations in the case of n=200 will be large while for n=2.000 it will have barely any effect.
I calculated the average of the distances between the smoothed hazard curves after and before the perturbation for each simulation. Here I present the average of 1.000 simulations:
n=200 n=780 n=2.000
Perturbation Abs. value Sq. value (10-2) Abs. value Sq. value (10-2) Abs. value Sq. value (10-2)
Kernel 0.27359 9.94977 0.02673 0.12977 0.00463 0.00522
Local poly. 0.00270 0.00112 0.00091 0.00011 0.00037 0.00002
Splines 0.00317 0.00130 0.00081 0.00007 0.00032 0.00001
Average of the distance in absolute value or squared between the smoothed curves after and before the perturbation.
For n=200 the local polynomial smoothing technique is less affected by the inclusion of two observations with both measures of error. However, with larger sample size the spline smoothing technique seems more robust. Another interesting but expected result is that with increasing sample size the effect of the perturbation is less important.
8.4 Bias concern
A final consideration is the possible presence of bias in our estimates. It is known in the literature (Stute (1994), Khan & Shaw (2013)) that the Kaplan-Meier estimator may present bias if the largest survival time is censored, which is our case. In order to assess if our Kaplan-Meier estimations were biased I simulated 5.000 data sets of different sample sizes and clustered the number of patients at risk in each month in order to get an approximation of the survival probabilities and the hazard rate at each month to compare it with the hazard rate estimated from the parametric model.
I took the mean and the median cumulative hazard in each month and plotted them along with the parametric estimation of the cumulative hazard evaluated at each month. This is the result for the sample size of 780:
30
In this case we observe no bias, in fact, the ratio between the Kaplan-Meier and the parametric estimations at each month had an average of 0.007 for 1.000 simulations and of 0.001 for 5.000 simulations, so we can expect that it approximates 0 while the number of simulations tends to infinity.
However, following the same procedure with the hazard rate I got the following:
31
For the other two sample sizes under study I found that the cumulative hazard is still unbiased while the bias of the hazard rate diminishes when the sample size increases.
This means that if we apply the Kernel smoothing technique we first need to correct for the bias. For the other two techniques no correction is needed since they are derived from the cumulative hazard.
32
9 Selection of the smoothing technique
Before the estimation of the parametric model and the analysis of simulated data, we were leaning towards the penalized splines technique. This was based basically on the comparative stability of the smoothed hazard rate and the property of an almost linear right tail.
With the analysis performed making use of simulated data sets we can conclude that penalized splines are the best smoothing technique for our purpose for the following reasons:
Spline smoothing is as good as smoothing with local polynomials in terms of proximity of the smoothed curve to the real data if we argue that the absolute value of the errors is a better measure than the mean squared error for our concern.
Penalized splines present the narrowest confidence interval.
They are more robust to outliers for medium and large sample sizes.
There is no need for bias correction.
For all these reasons in order to perform an analysis with survival data that concerns the comparison of the hazard rate I recommend estimating the hazard rate non-parametrically taking the derivative of the smoothed cumulative hazard with penalized splines. The only consideration needed is that in order to obtain a sufficiently smoothed hazard, the degrees of freedom used must be those that maximize the generalized cross validation score at smoothing the hazard. However, it is important to notice that when the sample size is large, the performance of each smoothing technique is very similar.
10 Comparison to the General Population death rates
We could say that the risk of dying for melanoma patients is comparable to that for the general population if we observe that the lower bound of the estimated confidence interval of the hazard rate is lower than the tabulated death rates of the general population.
10.1 Comparison of the hazard rate for the real data set
The average age of diagnoses in the database is approximately 61. For this reason I chose for the comparison the death rates of people between 61 and 71 years of age. This is the plot for the real data set where we observe that the general population death rate is indeed at some points above the lower bound which would lead to the conclusion that death rates of melanoma patients after 6 years of diagnoses and treatment are not statistically different from those of general population.
33
10.2 Simulation procedure for the confidence intervals
With the same arguments as before, the conclusion cannot be based in the analysis of a single data set, which is why I resort again to simulations.
In order to simulate censored survival data, two survival distributions are required, one for the uncensored survival times that would be observed if the follow-up had been sufficiently long to reach the event and another representing the censoring mechanism.
As explained before, if the censoring is independent from the survival times then censoring times do not need to be modelled because by estimating the survival function or the hazard rate properly taking into account the censored observations, these functions will be unaltered. However, in the estimation of confidence intervals, the presence of more or less observations is important, so when estimating proper confidence intervals the simulation of censoring times is needed.
The empirical survival distribution from a real data set provides a reasonable choice for the survival distribution (Burton et al. (2006)). Random non-informative right censoring with a specified proportion of censored observations can be generated assuming a particular distribution for the censoring times. The survival times incorporating both events and censored observations are calculated then by combining the uncensored survival times and the censoring times. If the uncensored survival time for a case is less than or equal to the censored time, then the event is considered to be observed and the survival time equals the uncensored survival time, otherwise the event is considered censored and the survival time equals the censored time.
The distribution I chose for the censoring times was the uniform distribution ( ). Its support is defined by those two parameters that are its minimum and maximum values. For censoring times we need the minimum value to be zero and the maximum value can be found by iteration to approximate the censoring probability in the real sample. I found the value of by iterations knowing that the true value was between 120 and 240 and then taking the mid-point until the
34
average censoring probability of 1.000 simulations differed less than 0.001 to that of the real sample. The values for and the number of iterations required were:
n=200 n=780 n=2.000
198.8 202.5 202.5
Iterations 5 4 4
Final value of and the number of iterations needed for each sample size
Having all the parameters needed I performed 1.000 simulations with this methodology estimating for each simulated data the corresponding confidence interval with the bootstrap methodology.
10.3 Comparison of the hazard rate for simulated data sets: Approach 1
To determine if the hazard rates of melanoma patients cannot be statistically different from those of the general population, I constructed a variable to count in how many simulations the general population rate is inside the confidence interval at least for one moment in time. This measure is chosen since given the non-perfect linearity of both the hazard rate and the lower bound of the confidence interval, it is possible that some points are above the lower bound, some subsequent points are below and the following points above again. For the comparison I consider enough if at least one point is above the lower bound.Sample size Simulations above the lower bound Not modelling censoring Modelling censoring
200 86.4% 94.5%
780 9.9% 28.6%
2.000 0% 1.2%
Percentage of simulations for which at least one point of the lower bound was above the regression line of the general population death rates
First we observe that when the censoring mechanism is not modelled the presence of a larger number of observations provides narrower confidence intervals and thus the acceptance of the hypothesis a smaller number of times.
However, even modelling the censoring mechanism for a sample size similar to our data set 29% of the simulations comparable to the death rates of general population is not enough to assert as a general conclusion that the death rates are comparable. We can only say that the death rates of melanoma patients are comparable to those of the general population when we simulate samples of 200 observations, since then the uncertainty around the observations gives sufficiently wide confidence intervals.
There are two important points to discuss in this first methodology. The first one is that the confidence interval of each data set is estimated by bootstrapping the values of that simulation only, and it may be that the coverage probability is not exactly 90%. Even if the coverage probability is correct, a second point to discuss is that the counter is activated regardless of there being one point or 50 points of the general population death rates above the lower bound. If there is only one point, considering that we are dealing with pointwise confidence intervals, it is likely
35
that the conclusion is erroneous since we expect 5% of the lower bounds above the true value, so between months 70 and 120 we would expect 2 or 3 erroneous estimations.
For these reasons I propose a second approach to answer the research question.
10.4 Comparison of the hazard rate for simulated data sets: Approach 2
In order to check if the estimated confidence intervals had the correct coverage probability I took 100 bootstrap samples from each of the 1.000 simulated data sets to determine the percentage of times that the estimated values at each point in the grid were inside the confidence interval and also above the lower bound. As we saw before, the hazard rate for small survival times differs considerably when it is the spline smooth version or the parametric estimation, so the average over the whole grid of the coverage probability is much lower than the average over the final points in the grid in which we are interested, i.e. from month 70 onwards.Nominal value Coverage
Inside Conf. Interval
All the curve
90% 70.53%
In the tail 92.23%
Above lower bound All the curve In the tail 95% 83.77% 97.61%
Coverage probability of the confidence intervals and the lower bound for all the curve and from month 70 onwards
For the confidence interval of our real data set, what we observe is that the coverage probability along the curve is considerably smaller than the nominal value estimated with a sample size of 780 observations while in the right tail the coverage probability is almost correct. However the lower bound in the tail (what we are interested in) is lower than it should be since we obtain in average almost 98% of the hazard rate points above the estimated lower bound.
Here I present the graph of the confidence interval of the real data set along with that estimated by the simulations (with the correct coverage probability) and in the left the parametric hazard rate with the correct confidence interval estimated for each sample size.
36
Hazard rate real data Hazard rate and confidence intervals
If I add the regression line of the general population death rates what we observe is that only the hazard rate of a sample size of 200 is comparable.
This second approach has a coverage probability that is the same as the nominal value and the 95% probability statement is correct. The disadvantage is that under this approach I rely completely in the parametric model, while this thesis consists in the use of smoothing techniques to the non-parametric estimation of the hazard with the argument that the parametric model is not flexible enough to capture the true hazard and it is not sustainable from a biological perspective.
