A Real Data-Driven, Robust Survival Analysis on Patients who Underwent Deep Brain Stimulation for Parkinson’s Disease by Utilizing Parametric, Non-Parametric, Semi-Parametric Approaches

1. Introduction

Parkinson’s Disease (PD) is a devastating neurodegenerative disorder that affects the brain neuron system and medical professionals and research scientists have been studying PD for a long period of time. However, the cause of PD still remains undetected. If the cause of PD is unknown it is called idiopathic PD and is the most common type of PD. Also, patients who get PD before 50 years of age are known as young-onset Parkinson’s disease. Unintended/Uncontrollable movements such as shaking stiffness, and loss of balance are some of the major symptoms of PD, among others. Symptoms of PD start slowly, progress over time, and finally, get worse over time. Patients in the later stage of PD are unable to walk or talk, [1,2]. Moreover, PD is a well-known symptom of the devastating dementia disease.

Over 10 million people around the globe suffer from PD and more than 50,000 people are diagnosed with PD each year only in the United States. However, these numbers do not represent thousands of undetected cases of PD around the world, [3]. In our previous studies, we have utilized a real data-driven machine learning modeling approach to detect PD patients with a very high degree of accuracy, [4] and invented a real data-driven analytical model to identify the causes of PD through Unified Parkinson’s Disease Rating Scale (UPDRS), [5].

Research scientists have been working through the last few decades on therapies and drugs to slow down or cure this devastating disease. Lepadopa and Dopamine agonist are the most common medications among others to slow down this devastating disease, [6]. Also, Deep Brain Stimulation (DBS) is a surgical therapy that has been used as a treatment to improve the movement functions and the quality of life of patients with movement disorders. Since movement disorders are one of the major symptoms of PD, DBS is widely used to treat PD. DBS is one of the most powerful innovative treatments that showed promising results for the movement symptoms of advanced PD, such as bradykinesia, essential tremor, and other motor complications of PD. Also, it reduces other medication usages of an individual patient with PD, [7,8,9].

In the literature we found that several studies have confirmed that the survival times of PD patients increase with DBS therapy compared to the PD patients who are taking other treatments, [10]. For instance, we found several articles that have been trying to analyze the long-term survival outcome of PD patients who underwent DBS, and all the articles uniformly conclude the significant overall improvement of the patient’s motor symptoms. Also, medication like levodopa equivalent daily dose was reduced over 55% for 5 years of follow-up after DBS, [11,12]. Furthermore, we found that a greater number of publications shows over 98% survival rate after 5 years of DBS treatment and also it shows a very high survival rate even after 3 years of DBS treatment, [13,14,15]. In another study the estimated long-term (10 years) survival probability was identified as 51%, [2]. However, in the literature, most researchers utilized the most widely used nonparametric Kaplan-Meier approach to perform the analysis and made conclusions about the survival of the patients who underwent DBS for PD, but it is neither a powerful nor a robust method. Furthermore, researchers utilized the semiparametric Cox-Proportional Hazards approach without addressing its underlying important assumptions and made conclusions about the survival of the patients and these results may be misleading because of the inappropriateness/invalidness of the developed models with unsatisfied underlining assumptions.

In this study, we analyze the survival of the patients who underwent DBS for PD more sophistically and investigate the sustainability of the long-term survival of PD by utilizing well-defined robust statistical methodologies such as Nonparametric, Semiparametric, and Parametric. Thus, the rest of this paper is organized as follows. Section 2 and Section 3 consist of a descriptive analysis of the data and a brief introduction to sophisticated statistical methodologies, respectively. Finally, Section 4 and Section 5 discuss the results and the conclusions of the study highlighting the usefulness of this study and concluding remarks.

2. Materials and Methods

2.1. The Real Data

The real data consists of information on individuals who underwent the DBS treatment for PD starting from 1999 to 2006 with at least 10 years of follow-up at Pennsylvania Hospital. We had the opportunity to get access to this data through the contact of Dr.Ashwin Ramayya (MD), [2]. The data consists of a total of 320 individuals who underwent DBS for PD. However, we considered only complete observations (deceased patients) for our study of survival analysis and it consists of only 149 individual patients after removing the censored observations and the missing observations without using any data imputation method to replace those missing survival times of the individual patients.

Age, Age at Implant, Gender, Survival time (in months), Initial implant side, and Number of revisions are the six individual covariates that we have used to develop the present study. The final data frame consists of individual patients in the range of age from 40 years to 93 years old and the mean age of the patients who underwent DBS surgery is 64.41 years (± 8.59 standard deviation). Also, 73% out of the total individual patients are Males. Furthermore, 69% of the patients underwent bilateral DBS for PD. The following data diagram in Figure 1 and Table 1 illustrate the exploratory data analysis for some continuous covariates.

2.2. Methodology

In this study we utilize the well-known non-parametric Kaplan-Meier (KM) approach, robust Parametric approach, and widely used semi-parametric multivariate Cox-PH method to analyze the long-term survival outcomes of the individuals who underwent DBS for PD. In the following section, each of the proposed models and their usefulness is briefly discussed.

2.2.1. Kaplan-Meier Survival Analysis

Kaplan-Meier (KM) estimation of the survival time is a fully non-parametric approach for analyzing the time of occurring the event of interest (time until death) and it was first introduced in 1958, [16]. KM is a well-known approach in the field of survival analysis and reliability analysis. However, KM is not as powerful nor robust as the Parametric or Bayesian survival analysis since it does not make any distributional assumption about underlying survival time. KM involves estimating the survival probabilities of patients at a certain point of time as follows,

$$S_t={\displaystyle \frac{Numberofpatientslivingatthestart\left(t_0\right)-Numberofpatientsdiesattimet}{Numberofpatientslivingatthestart\left(t_0\right)}}.$$

Suppose that there are n individual patients at the start of the time ($t_0$) and their respective observed event times (died at time $t_i$) are $t_1<t_2<t_3<...<t_k$ with $d_i$ units fail at time $t_i$. Then the KM estimator of the survival function $S_{KM}\left(t\right)$ is given by following equation,

$$S_{KM}\left(t\right)=\prod _{i:t_i<t}{\displaystyle \frac{(n_i-d_i)}{n_i}}=\prod _{i:t_i<t}(1-{\displaystyle \frac{d_i}{n_i}}),$$

(1)

where $t>0$, $n_i$ is the number of patients at risk at time $t_i$ of the study, $d_i$ is the number of patients who died at time $t_i$. We proceed and calculate survival probabilities using KM estimators given by the Equation (1) under the following assumptions, [17,18,19],

• Patients who have been censored have the same survival prospects as those who continue the follow-up.
• Survival probabilities are the same for individuals who were recruited early and late in the study.
• Event happens at the time specified.

2.2.2. Parametric Survival Analysis

In this section, we briefly discuss the more powerful and robust parametric survival analysis and its usefulness.

Let T be a non-negative random variable, which measures the survival time of an individual patient, and its probabilistic behavior is characterized by the well-defined probability density function $f\left(t\right)$. Then the probability that an event occurs prior to time t can be defined as the cumulative distribution function and is given below by Equation (3.1.2),

$$F\left(t\right)=P(T\le t)={\int }_0^{t}f\left(t\right)dt,t>0,$$

(2)

where $f\left(t\right)$ is the probability density function of the survival time.

The survival of an individual patient after time t is given by the Survival Function and it can be derived from the above CDF. The generalized analytical form of the survival function can be derived as follows in Equation (3).

$$S\left(t\right)=1-F\left(t\right)=P(T>t)={\int }_t^{\infty }f\left(t\right)dt,t>0.$$

(3)

One can utilize Equation (3.1.2) and Equation (3) to derive the hazard function to measure the chances of instantaneous death of a given patient at any time t as follows in Equation (4),

$$h\left(t\right)=\underset{\delta t\to 0}{lim}{\displaystyle \frac{P(t<T\le t+\delta t|T>t)}{\delta t}}={\displaystyle \frac{f\left(t\right)}{S\left(t\right)}}.$$

(4)

2.2.3. Multivariate Cox-Proportional Hazards Model

In the context of survival/reliability analysis, the hazard function plays a major role. It explains the instantaneous risk of an event at time t. That is, the probability of the desired event occurring any time after time t is the hazard function under the assumption that individual observation has not experienced the desired event until time t.

Let T be a random variable that represents the survival time of a given patient. Then the probability that an event occurs prior to or at time t is defined by the cumulative distribution function (CDF) and is given below,

$$F\left(t\right)=P(T\le t)={\int }_0^{t}f\left(t\right)dt,t>0,$$

(5)

where the probability density function can be derived as follows,

$$f\left(t\right)={\displaystyle \frac{dF\left(t\right)}{dt}}.$$

(6)

The probability of a given patient surviving after time t (Survival function) can be defined by the following Equation (3.1.2),

$$S\left(t\right)=1-F\left(t\right)=P(T>t)={\int }_t^{\infty }f\left(t\right)dt,t>0.$$

(7)

Utilizing Equation (6) and Equation (3.1.2), the hazard function for an individual patient at time t can be derived as follows,

$$h\left(t\right)=\underset{\delta t\to 0}{lim}{\displaystyle \frac{P(t<T\le t+\delta t|T>t)}{\delta t}}={\displaystyle \frac{f\left(t\right)}{S\left(t\right)}}.$$

(8)

To estimate the hazard function the Cox Proportional Hazards Model (Cox-PH) was introduced by D. R Cox in 1972, [20], which specifies the hazard function directly rather than deriving it as the ratio of probability density function and survival function. The analytical form of the Cox-PH model is given by Equation (9),

$$h_i\left(t\right)=h_0\left(t\right)exp\left(\sum _{j=1}^p{\beta }_j{X}_j+\sum _{j=1}^p{\alpha }_{jK}{X}_j\ast X_K\right),$$

(9)

where $h_0\left(t\right)$ is the baseline hazard function and it depends on only the survival time. ${\beta }_j$ is the coefficient parameter for $j^{th}$ covariate and ${\alpha }_{jk}$ is the coefficient parameter for the interaction between $j^{th}$ and $k^{th}$ covariates. In the Cox-PH model, the baseline hazard function $h_0\left(t\right)$ is unspecified and it is treated as a non-parametric function of time t. However, the Cox-PH model assumes the time independence of the covariates, linearity in the covariates, additivity, and proportional hazards. Thus, the Cox-PH model is considered as a semi-parametric survival model, [21]. To estimate the parameters of the Cox-PH model we cannot use the full likelihood approach since the baseline hazard function is unspecified. Thus, the partial likelihood approach proposed by Cox in 1979 is utilized to estimate the model coefficients of the Cox-PH model, [22].

In many situations, researchers are interested in the factor exp(coef) of the model rather than the estimating entire $h_i\left(t\right)$ which determines how an individual patient varies as a proportion to the baseline hazard function. Equation (9) can be rearranged to find the hazard ratio for a given individual patient as given below,

$$HR={\displaystyle \frac{h_i\left(t\right)}{h_0\left(t\right)}}=exp\left(\sum _{j=1}^p{\beta }_j{X}_{ij}\right).$$

(10)

Suppose that we have two types of treatment groups, namely treatment group (T) and placebo group (P). Then we can extract the following important and useful information about the two treatment groups using the HR,

• HR $>1({\beta }_i>0$). The treatment group has a higher hazard over the placebo group. (i.e., the placebo group is favored.)
• HR $\approx 1$. The treatment group and placebo group have no significant difference.
• HR $<1({\beta }_i<0$). The treatment group has a lower hazard over the placebo group. (i.e., the treatment group is favored.)

3. Results

Prior to model the survival times, we assess whether the median survival time of the Male patients who underwent DBS for PD is statistically significantly different from that of female patients by utilizing hypothesis testing such as Wilcoxon rank sum test and log-rank test, [23,24,25]. In the following section, we have defined our primary null hypothesis and alternative hypothesis,

$H_0$: Median Survival Time of Male patients is not different from that of Female patients.

Vs

$H_1$: Median Survival Time of Male patients is different from that of Female patients.

Figure 2. Visualization of Log-Rank Test Results

Figure 2. Visualization of Log-Rank Test Results

The Wilcoxon rank sum test gave a p-value of 0.029 and for the log-rank test p-value is 0.013. Both Wilcoxon and log-rank tests uniformly confined that p-values are less than the 5% significant level. Thus, we could reject the null hypothesis and concluded that the median survival time of male patients is different from that of females. Since the median survival time of males is significantly different from that of females we consider males and females as two different entities and continue our analysis separately for male patients and female patients.

3.1. Univariate Survival Analysis

In this subsection, we summarize the results of parametric and nonparametric univariate survival analysis for PD patients who had DBS.

3.1.1. Kaplan Meier (KM)

KM survival probabilities have been estimated for both male patients and female patients separately because initially, we identified that their survival times are statistically significantly different. The KM estimation results for the Male patients and Female patients are summarized in Table 2 and graphically illustrated by Figure 3 and Figure 4 for Females and Males, respectively.

When carefully investigating the above KM plots for the two genders it is clear that the median survival time of female patients is higher than the median survival time of the Male patients after performing DBS for PD.

3.1.2. Parametric Survival Analysis

In the case of Parametric survival analysis, it is essential to identify the probability distribution functions that characterize the probabilistic behavior of the survival time of the Male patients and Female patients separately since their median survival times are statistically significantly different. Thus, we identified the best-fitted probability distribution functions for both male and female patients that characterize the probabilistic behavior of their survival times and utilized well-defined goodness of fit tests, Kolmogorov-Smirnov, Anderson-Darling, and Chi-Squared to validate the resulted probability distribution functions. The goodness of fit test results are summarized in Table 3, given below.

The goodness of fit test results given by Table 3 is uniformly confirmed that the probabilistic behavior of the survival time of the Male patients are characterized by the 3-parameter Weibull distribution and its analytical form is given by the mathematical function 11, given below. Furthermore, parameters of the Weibull distribution were estimated by using the most commonly applied Maximum Likelihood Estimation (MLE) method, [26].

$$f_M\left(x\right|\theta ,\eta ,\tau )={\displaystyle \frac{\eta }{\theta }}{(x-\tau )}^{\eta -1}exp\{-{\displaystyle \frac{1}{\theta }}{(x-\tau )}^{\eta }\},\eta >0,\theta >0,x>\tau ,$$

(11)

where $\theta$, $\tau$ and $\eta$ are Shape, Scale, and Location parameters, respectively. MLE estimates for the model parameters are $\theta$= 4.9593, $\eta$ = 186.21, $\tau$= -76.262.

Once the parameters are estimated the shape of the Weibull distribution that characterizes the probabilistic behavior of the survival times of the male patients is graphically illustrated by Figure 5, given below.

From Figure 5 one can compute the probabilities of the survival times of the patients who underwent DBS for PD between time $t_i$ and $t_{i+1}$. As an illustration, we can compute the probability of survival between 25 months to 50 months for a PD patient who underwent DBS as $P(25\le t\le 50)$$\approx 0.00475-0.002\approx 0.00275$. In Figure 5 we can clearly see the probability of a patient who will survive between 25 months to 50 months is approximately $0.3\%$.

Once we have the analytical form that characterizes the probabilistic behavior of the survival time of the male patients, the Cumulative Distribution Function, (CDF) for male patients ($F_M\left(t\right)$) can be derived as follows in Equation (12) and its graphical illustration for male patients is by Figure 6.

$$F_M\left(t\right|\theta ,\eta ,\tau )=P(T\le t)={\int }_0^t{f}_M\left(x\right|\theta ,\eta ,\tau )dx,t>0,$$

$$F_M\left(t\right|\theta ,\eta ,\tau )=1-e^{-{\left({\displaystyle \frac{t-\tau }{\theta }}\right)}^{\eta }}.$$

(12)

Identifying the CDF for a given random variable is a crucial achievement in probabilistic analysis because the CDF determines the probability or likelihood that a random variable will take at a specific value or less than that value. Thus, one can estimate the probability that an individual patient who underwent DBS for PD will survive until time t. Suppose a PD patient who underwent DBS wants to know what is the probability that he live until 140 months. That is $P(t\le 140)\approx 0.88$ as sketched in Figure 6.

Furthermore, the patient can estimate his/her survival probability beyond 140 months by utilizing the survival function. Since we have already estimated the CDF the derivation of the Survival Function for the male patients and its analytical form is given by Equation (13) and its behavior by Figure 7.

$$S_M\left(t\right|\theta ,\eta ,\tau )=1-F_M\left(t\right|\theta ,\eta ,\tau )=P(T>t)={\int }_t^{\infty }{f}_M\left(x\right|\theta ,\eta ,\tau )dx,t>0,$$

$$S_M\left(t\right|\theta ,\eta ,\tau )=e^{-{\left({\displaystyle \frac{t-\tau }{\theta }}\right)}^{\eta }}.$$

(13)

Thus, as shown in Figure 7 an individual patient surviving beyond 140 months can be estimated using Equation (13) and that is approximately $0.12$.

According to the summary of the goodness of fit test results given in Table 3, the survival times of the Female patients are characterized by the Lognormal (3P) distribution and its analytical form is given by the following Equation (14).

$$f_F\left(x\right|\mu ,\sigma ,\gamma )={\displaystyle \frac{1}{(x-\gamma )\sigma \sqrt{2\pi }}}exp\{-{\displaystyle \frac{{[ln(x-\gamma )-\mu ]}^2}{2{\sigma }^2}}\},\gamma <x,-\infty <\mu <\infty ,\sigma >0,$$

(14)

where $\mu$, $\gamma$, and $\sigma$ are Scale, Location, and Shape parameters of the 3P-Lognormal distribution, respectively. We have utilized the Maximum Likelihood Estimation (MLE) approach to estimate those three parameters of the Lognormal distribution as illustrated in the literature, [27,28]. Thus, the estimated parameters are $\mu$ = 7.138, $\sigma$ = 0.03655, $\gamma$ = -1149.2. The graphical visualization of the shape of the Lognormal distribution is given by the following Figure 14. For example, one can interpret the estimated probability that a PD patient will survive between 50 months and 100 months by using the probability distribution of the survival time of the female patients as shown in Figure 14. That is, $P(50\le t\le 100)=0.0085-0.0032=0.0053$.

Once we derived the analytical form that characterizes the probabilistic behavior of the survival times of the female patients, the Cumulative Distribution Function (CDF) for female patients ($F_F\left(t\right)$) can be defined as the following Equation (15), given below, and its graphical illustration for female patients is given by Figure 9. We can use this CDF to estimate the probability of survival of the female patients until time t$P(t\le T)$. For example, suppose that we want to estimate the probability of surviving until 140 months ($P(t\le 140)$) using the CDF of the female patients given in Equation (15). Then the estimated probability is $P(t\le 140)\approx 0.75$. The graphical visualization of the estimated value is shown in Figure 9.

$$F_F\left(t\right|\mu ,\sigma ,\gamma )=\varphi \left({\displaystyle \frac{ln(t-\gamma )-\mu }{\sigma }}\right).$$

(15)

Finally, the estimated CDF can be used to derive the Survival Function for the female patients and its analytical form is given by Equation (16) and its graphical visualization of the behavior by Figure 10.

$$S_F\left(t\right|\mu ,\sigma ,\gamma )=1-F_F\left(t\right|\mu ,\sigma ,\gamma )=1-\varphi \left({\displaystyle \frac{ln(t-\gamma )-\mu }{\sigma }}\right).$$

(16)

Concerning the graphical illustration of the survival functions of the male and female patients, it is noticeable that the female patients have significantly higher median survival time than the median survival time of the male patients. This outcome is confirmed by the univariate results obtained from the KM survival analysis. However, the univariate survival analysis was performed under the assumption that all the patients died under the same cause regardless of their other information.

3.2. Multivariate Survival Analysis

3.2.1. Cox-PH Model

The Multivariate Cox-Ph model developing procedure was started with four individual risk factors that are identified by the medical professionals and believed to significantly contribute to the survival of a patient who underwent DBS, [Reference for Risk factors]. Thus, we initialize the model development procedure with four individual risk factors for both genders, but the performance of the models was not impressive. Thus, we utilized all individual risk factors and their two-way interactions along with the step-wise elimination method. Equation (17) expresses the analytical form of the final cox-proportional hazard model for male patients and Equation (18) gives the final analytical form for the female patients.

$$\widehat{h_M\left(t\right)}=h_0\left(t\right)exp(0.03074{\mathbf{X}}_1-2.90540{\mathbf{X}}_{2\left(1\right)}-3.09802\sqrt{X_{\mathbf{3}}}+2.63631{\mathbf{X}}_{2\left(1\right)}:\sqrt{X_{\mathbf{3}}}),$$

(17)

$$\widehat{h_F\left(t\right)}=h_0\left(t\right)exp(0.17810{\mathbf{X}}_1+2.45162{\mathbf{X}}_{2\left(1\right)}-2.07761{\mathbf{X}}_3),$$

(18)

where $h_M\left(t\right)$ is the hazard function for male patients at time t,

$h_F\left(t\right)$ is the hazard function for female patients at time t,

$h_0\left(t\right)$ is the arbitrary baseline hazard function,

$X_{i\left(l\right)}$ is the $i^{th}$ covariate with level l (for a categorical risk factor) as described in Section 2.

The above analytical models can be rearranged into a linear relationship between risk factors and logarithmic of the relative hazard ratios as follows for males and females, respectively.

$$log\left({\displaystyle \frac{\widehat{h_M\left(t\right)}}{h_0\left(t\right)}}\right)=0.03074{\mathbf{X}}_1-2.90540{\mathbf{X}}_{2\left(1\right)}-3.09802\sqrt{{\mathit{X}}_{\mathbf{3}}}+2.63631{\mathbf{X}}_{2\left(1\right)}:\sqrt{{\mathit{X}}_{\mathbf{3}}}$$

(19)

$$log\left({\displaystyle \frac{\widehat{h_F\left(t\right)}}{h_0\left(t\right)}}\right)=0.17810{\mathbf{X}}_1+2.45162{\mathbf{X}}_{2\left(1\right)}-2.07761{\mathbf{X}}_3$$

(20)

The Table 4, given below, consists of the information about the developed Cox-PH models for male patients and female patients including the estimates of model coefficients (parameters), hazard ratios (exp(coef)), standard errors of coefficients and 95% lower and upper confidence intervals for the hazard ratios. Also, Table 5, given below, provides the global statistical significance values for the overall significance of the proposed Cox-PH models by using three hypothesis tests: Likelihood ratio test, Wald test, and Score test and all the p-values ($<<0.05$) uniformly confirm that the proposed Cox-PH models are statistically significant at 5% level of significance. These three tests are asymptotically equivalent. However, for small samples, the likelihood ratio test slightly outperforms the Wald test and both outperform the Score test.

Since the hazard ratio (exp(coef)) of the Cox-PH model measures the relative risk it is possible to interpret the model coefficients in a meaningful way. If the model coefficient is positive ($\widehat{\beta }>0$), i.e., hazard ratio, (exp(coef)) > 1 then the event of occurrence is high with a bad prognostic factor. However, if the coefficient is negative ($\widehat{\beta }<$ 0), i.e., hazard ratio, (exp(coef)) < 1 then the event of occurrence is low with a good prognostic factor. In the case of proposed models, both male and female patients consist of the coefficient of $X_1$ (${\widehat{\beta }}_{MX1}=0.050>0$ and ${\widehat{\beta }}_{FX1}=0.114>0$) and are positive. It means the age at implant of both male and female patients has a higher risk for death with a bad prognostic. In the case of coefficients of $X_2$, the behavior of the two models is different from each other. For the female patients $X_2$ (${\widehat{\beta }}_{FX2}=1.014>0$) has positive coefficient and has higher risk for death with a bad prognostic while male patients $X_2$ (${\widehat{\beta }}_{MX2}=-2.74<0$) has negative coefficient and the event of occurrence is low risk with a good prognostic factor.

3.2.2. Validation of the Assumptions

In this section, we assess the five major assumptions of a Cox-PH model that need to be satisfied to perform a high-quality and robust survival analysis.

(1)

Proportional hazards assumption

We have performed the statistical test and also the graphical visualization based on Scaled Schoenfeld Residuals to validate the proportional hazards assumption on the Cox-PH model. Table 4, given below, consists of p-values for the test based on scaled Schoenfeld residuals. The P-values for all covariates are greater than $0.05$ at the 5% level of significance confirming the fact that our proposed models for male and female patients satisfy the proportional hazards assumption.

Next, we visualize the graphical approaches of the scaled Schoenfeld residuals given in Figure 11 and Figure 12 for Male patients and Female patients, respectively. Both figures indicate a zero slope pattern by confirming the proportional hazard assumption.

(2)

Linear functional form of Continuous variables

The linearity Assumption is one of the major underlying assumptions of the Cox-PH model that should be satisfied by the final model to achieve more robust results without misinterpretation. The graphical illustration of the Martingale residuals for the continuous covariates can be used to identify a violation of this assumption. Figure 13 and Figure 14 uniformly confirm that both Cox-PH models given by Equations (17) and 18 have satisfied the linearity assumption as given below for the continuous risk factors.

(3)

Check for Extreme Observations

Under the next assumption, we check for the extreme values in our final Cox-PH model by utilizing Devience residual plots as given by Figure 2 and Figure 16. Both figures don’t consist of any extreme outline observations. Thus, we can conclude that the proposed models have met the underlying assumption of no outliers of extreme observations.

(4)

Overall Fit of the Cox-PH Model

Finally, we assessed the overall fit of the proposed Cox-PH models using the cumulative hazard of the Cox-Snell residuals plot. Figure 17 and Figure 18, given below, graphically visualize the Cox-Snell residuals plots for both genders and they follow a straight line with a zero intercept and a unit slope by confirming the high quality of the proposed models.

Thus, the proposed models for Male patients and Female patients have satisfied the underlying assumptions of the Cox-PH model and its high quality.

4. Discussion and Conclusions

In the scope of this paper, we have utilized nonparametric, semi-parametric, and parametric survival analysis to perform systematic and deep univariate and multivariate survival analysis for the patients who underwent DBS for PD. Within the first part of this study, the univariate survival analysis was performed through widely used nonparametric KM estimation method and more powerful and robust parametric approaches under the assumption that all patients deceased due to the same cause of death. Thus, we were able to compute and compare the survival probabilities using nonparametric KM and parametric methods for an individual patient. Since the survival times of the male and female patients indicated a statistically significant difference as shown in Section 3, we performed survival analysis for male and female patients independently. In the sense of parametric survival analysis, 3-Parameter Lognormal Distribution probabilistically characterizes the survival times of the Female patients while the probabilistic behavior of the survival times of the Male patients is characterized by 3-Parameter Weibull Distribution.

Figure 19 and Figure 20, given below, graphically illustrate the comparison of the parametric and nonparametric survival curves for female patients and male patients, respectively. According to the graphical illustration, it is clear that KM method underestimates the survival probabilities of the female patients while it overestimates the survival probabilities for the male patients. At this point of our study, we rely on survival probabilities estimated by the robust and more powerful parametric approach because of its robustness as explained previously compared to the widely used KM method. As shown in Figure 19 the median survival time of the female patients is approximately 110 months while in Figure 20 the median survival time for male patients is approximately 95 months. Thus, we can conclude that the female patients who underwent DBS for PD have higher survival compared to that of male patients. However, in the case of univariate survival results, they are based on the assumption that all the patients’ underlined cause of death is the same without considering personalized risk factors or other relevant information.

Then, the effect of covariates was taken into account for both female and male patients separately and developed their own Cox-PH models. In the developing process, we considered the three individual covariates: Age at Implant, Initial Side of Implant, and Number of Revision. Furthermore, the two-way interactions of the individual covariates are taken into account. However, it is important to note that in the literature, researchers do not consider interaction effects when developing a Cox-PH model. In Section 3.2, we have developed the independent Cox-PH models for male patients and female patients. The final model for Male patients given by Equation (17) consists of three individual statistically significant covariates: Age at Implant, Initial Side of Implant, the square root of the Number of Revision and a statistically significant two-way interaction term between Initial Side of Implant and the square root of the Number of Revision. Moreover, the final Cox-PH model for female patients consists of only three individual statistically significant covariates: Age at Implant, Initial Side of Implant and Number of Revision as given by Equation (18). In the model 17, the individual covariate, Initial side has $HR<1$. However, when it comes to the interaction between the Initial Side of Implant and the Number of Revision it has $HR>1$. This is a perfect example to illustrate the importance of adding interaction terms. Furthermore, in Section 3.2.2 we have addressed all underlying assumptions of the Cox-PH model and they uniformly confirmed the robustness of our developed model and its high quality.

Once the Cox-PH model is finalized by addressing all underlying assumptions, the following important results are available concerning the long-term survival outcomes of the patients who underwent DBS for PD.

• For the covariate, Initial Side of Implant, the HR value (2.76) for female patients is greater than 1, which means the risk of death is higher with a bad prognosis. That is, Initial Right Side Implant has a higher risk of death than Initial Left Side Implant. For the male patients, the HR value (0.064) is less than 1, which means the risk of death is lower with a good prognosis.
• For the interaction between Initial Side of Implant and the square root of Number of Revisions, the HR values are greater than 1, which means the risk of death is higher with a bad prognosis. That is, Number of Revisons that a patient has done has a stronger effect on hazard ratio with the Right side of the Initial Implant compared to the Left side of the initial implant.
• Based on the hazard ratios, the most contributing individual covariates are Age at Implant, Initial Side of Implant and Number of Revision for male patients, and Initial Side of Implant, Age at Implant, Number of revision for female patients, respectively. Furthermore, the ranking of the individuals and interaction terms for both models are given by Table 7.