Introduction

The Primary biliary cholangitis (PBC) or primary biliary cirrhosis is a autoimmune disease leading to bile release in the liver[1]. One of the major role of the liver is bile production, a secretion involved in cholesterol and other cell wastes processing. For that matter, bile is being further transported to the digestive apparatus parlty through bile ducts. PBC is being experienced as bile conduct damages causing liver tissue scarring upon bile release leading to liver malfunction or failure of such organ[2]. As liver is one of the body vital organs, PBC in its latest stages unfortunalty leads to patient death if he is not subject to liver transplant.

In order to better diagnose and cure PBC, it is of the outmost importance to understand the risk factors and metrics indicators of such disease as well as the potential effect of given drug therapies. Mayo clinic trial have been aiming at performing a randomized placebo clinic test for D-penicillamine[3], a drug meant for PBC treatment. They have gathered a data set of 424 PBC patients in a ten year interval. Although most interests have been around drug therapy investigation, this trial also assesses other auxilary variables such as clinical and biological metrics.

This study consists of a Survival Analysis based on Mayo Clinic trial data set. In fact, it aims at investigating the outcome of D-pencilamine as well as other clinical and bilogical metrics on the survival probability of patients with PBC. Such investigation has been conducted through Keplen-Meyer estimates of the survival function. Further analysis of the data also lead to Cox regression to fit the hazard function of the resulting model. Explory data anlysis was also performed in this study to get a better grasp of the data set other characteristics.

Exploratory Data analysis

The Mayo Clinic data set regroups a total of 424 patients with PBC. 312 of those patients did participate to the randomized placebo trial metting the eligibility critetria where as the remaining 122 did not participate but agreed to be followed up for basic measurments and survival and a total of 6 of them was lost to follow up making them only 106. Since this study focuses on the survival probaility and hazard functions, it is important that the dataset contains a time variable (FUDAYS) indicating the time to the event of interest (here being the death of the patient) as well as indication on censored data (STATUS), 0 for censored patients, 1 for censored patients due to liver transplant, 2 for death. Therefore, the data set also contains the following variables :

First data analysis showed a certain amount of missing data. In fact, the extra trail patients were not followed up for ascites, hepatomes, alkaline phosphatase, cholesterol, aspartate transaminsase and urine cooper. As our data is both composed of continuous and categorical variables, data exploration had to be perfomed accordingly and as such, histograms were used to display the categorical variables distributions where as boxplots of were realized for the continuous ones.

\label{fig:EDA-barplot} Barplot of selected variables: Drug treamtment,Sex of patients, presence of Ascites, presence of Hepatomes, presence of Edma, Stage of the disease, Status of disease, Age of patients

Barplot of selected variables: Drug treamtment,Sex of patients, presence of Ascites, presence of Hepatomes, presence of Edma, Stage of the disease, Status of disease, Age of patients

Data distribution across categories allows to apprehend how balanced the various fatures are. As depicted in figure , the trial comprised less male patients, less patients having ascites and no edema. The majority of patients enrolled in the trial are in late stages of the desease (3 or 4). On the other hand, it is seen that the number of patients in both the control and treamtent group is more less the same as well as for the prensence and absence of hepatomes and patient distribution across age groups. It is also observed that transplanted patients only account for 25 of the cases. Finally as pbc is responsible for liver failure and induced death, it is important to investigate such event of interest related to disease stage.

\label{fig:EDA-barplot5} Barplot of censored data repartition across disease stages

Barplot of censored data repartition across disease stages

For that purpose, the number of patients per status for each stage of the disease was ploted in Fig. . At early stages, most of the data is censored, but the further we move in disease stages the more the enrolled patients die (48 and 84 in stage 3 and 4 versus 0 and 23 in stage 1 and 2). With that being said, one should mention that liver failure endured by the patients in late stages of the disease have caused their death (being the event of interest) as natural death cannot be the only case here by looking back at the patients repartition per age group in (roughly a 100 per age group). It can also be seen that essentially most liver transplants were being exeperienced in later stages (3 to 4) indicating sever liver damages at that point.

Getting insights on data distribution is also a key aspect of exploratory analysis on bivariate data. As such, boxploting the different variables in Figure provides a five number visual summary on each distribution. Length of the boxplot is defined as the interquartile range which is the difference between first and third quartile. The larger it is, the sparser the data is. Hence, one can mention the high diversity in trial days as it encompass a sparse time frame periode with a median of 1730 as well as for the age of the enrolled patients (median : approximatly 51 years old). As for the other variables boxplots appear to be quite narrow, they seem to display many outliers. However, such outlier are much more numerous for bilirubin, cholesterol, urine cooper and alkaline phosphatase concentrations. Those outliers might come from higher concentrations or higher metric values experienced for uncesored data. Indentifying explanatory variables will help us understand the underlying reason for such observations.

\label{fig:EDA-box}Boxplots of continuous variables

Boxplots of continuous variables

Keplen Meyer

To asses the effect of D-penicillamine treatment therapy, survival as a function of drug treatment was estimated using the Keplan-Meier estimator as displayed in figure . Keplen-meyer survival prbability function [4] is given by the following equation : \[ S(t_{i}) = S(t_{i-1})(1-\frac{d_i}{n_i}).\]

where \(S_{ti}\) is the probability of being alive at \(ti\), \(d_{i}\) the number of events at ti and \(n_{i}\) number of patients alive at \(ti\). One could notice that the therapy does not seem to be effective as the survival probability is approximately similar between treatment and control group. More over, this is confirmed as the difference was not found to be statistically different (p=0.9281) following a Log-rank test.

logrank_test(Surv(FUDAYS, STATUS) ~ DRUG, data=mydat)
## 
##  Asymptotic K-Sample Logrank Test
## 
## data:  Surv(FUDAYS, STATUS) by DRUG (1, 2, extra)
## chi-squared = 0.14925, df = 2, p-value = 0.9281
\label{fig:KM-Drug}Keplen Meyer estimation of the survival function for standard and test therapies Log rank p-value = 0.9219

Keplen Meyer estimation of the survival function for standard and test therapies Log rank p-value = 0.9219

One can conduct the survival analysis further and investigate the remaing variables effect on the survival probability. Regarding the categorical variables, ascites and spiders presence, edema presence with repsect to diuritic use and disease stages were found to have signicative impact on such function as depicted in figure

\label{fig:KM-Cat}Keplen Meyer estimation of the survival function for the selected categorical varibales (log rank test p-values): Ascites (p-value = 0.00069), Spiders (p-value = 0.001814), Edema (p-value = 0.0001758), Stages = ( p-value = 3.165e-05)

Keplen Meyer estimation of the survival function for the selected categorical varibales (log rank test p-values): Ascites (p-value = 0.00069), Spiders (p-value = 0.001814), Edema (p-value = 0.0001758), Stages = ( p-value = 3.165e-05)

Hence, it is seen that the absence of ascites and spiders yields to a signficantly lower probability of survival. Same comment can be made for edema and stage levels, the higher they are, the lower the survival probability is. Therefore, those metrics can be considered as indicative of the disease seriousness and impact on patient.

Regarding categorical variables impact on the Keplen-Meyer estimate of the survival time, based on the Fig. min, Q-1, median, Q-3 and max values of each variable, one can perform data munging by turning them into categorical variables and assess their impact using this method.

\label{fig:KM-Cont}Keplen Meyer estimation of the survival function for the selected categorical varibales (log rank test p-values): Age (p-value = 0.006366), Bilirubin (p-value = 5.476e-05), Cholesterol (p-value = 0.01477), Protime = ( p-value = 0.001538), Albumin = (p-value = 2.592e-05)

Keplen Meyer estimation of the survival function for the selected categorical varibales (log rank test p-values): Age (p-value = 0.006366), Bilirubin (p-value = 5.476e-05), Cholesterol (p-value = 0.01477), Protime = ( p-value = 0.001538), Albumin = (p-value = 2.592e-05)

Hence, by looking at Fig. it can be seen that age, bilirubin and albumin concentration had signifcant impact on the Keplen-meyer estimates. The older the patient is,the higher his bilirubin concentration and the lower his albumine one are, the lower his survival probability will be. Significant effect has also been experienced for protime duration and cholesterol concentration.

Cox proportional hazards model

In order to move on with our survival analysis, Cox-proportional hazards model is an appropiate, commonly used tool to investigate the impact of multiple covariates on the survival time[5]. In fact, Keplen-Meyer estimates only permit univariate categorical investigation of such relationship. Therefore, Cox mutilivariate regression assumes linearity on the log hazard scale and is given by Eq. :

\[\begin{equation} \label{eq:coxph} h(t) = h_{0}(t) \exp\Bigg(\sum_{i=1}^{k} \beta_{i} x_{i}\Bigg) \end{equation}\]

\(h\) being the hazard function (i.e. the instantaneous failure rate), \(x_{i}\) the covariates of interest, \(\beta_{i}\) the coefficient for covariate \(x_{i}\), and \(h_{0}\) the baseline hazard function. This last parameter depicts the common shape of survival time distribution for all patients.

For Cox analysis purpuses, only Mayo Clinic trial patient were kept as the remaining 106 patients were displaying many missing entries for tested variables. More over, the 19 liver transplant censored patients were viewed as a random sample and discarded for the rest of the study [6].

Multivariate cox analysis was first conducted on the “full model” i.e containing all the covariates. Such model was globaly found to be statistically significant for the three following metrics : Likelihood ratio test, Wald test, Score (logrank) test with p-values the order of 10e-15. Such model has also revealed the statistical significance of particular covariates as depicted in (Table ).

However, the Cox multivariate global as for some of the covariates did not fulfill the model’s assumptions[7]. In fact, for the model to be valid, the covariates as well as the model itself have to be time indepedent. Such assumption is tested and results are displayed in (Table ). Covariates as well as the model are to be time independent if the p-values displayed in the hazard test are not statistically significant. Therefore, it can be seen that some covariates as well as the model itself display signficant p-values (i.e <0.05) meaning that the null hypothesis stating the time independency cannot be rejected.

With that being said, a different approach had to be adopted in order to fit an accurate Cox regression model that validates time independency assumptions. Therefore, forward selection method was used to recursively add previously selected covariates while verifying at each time that Cox assumption still holds for the fitted model. Covariate selection was based on Keplen-Meyer univariate significant variables as well as full model Cox regression significant features. Therefore, bilirubin, albumine, urine cooper, prothrombin time, edema, cholesterol, stages, age, ascites and spiders were selected for a foward selection procedure. To that end, covariate transformation such as logarithmic transformation or stratification for categorical variables were utilised to pass Cox model assumptions when features were found to be time dependent. Model overall performance was juged using Akaike Information Criterion (AIC) as it ensures both model accuracy and simplicity. The final model is depicted in (Table ). Edema had be to be stratified and no other variables previuously selected and mentioned above have improuved the model once added (based on AIC) All covariates are significant as well as the global model (p<2e-16) are signifcant for a level 0.001.

Such model fulfills model’s assumptions as none of the shown p-values bellow (Table ) are statistically significant. Time indepedence can also be seen in the Schoenfeld residual plots Fig. as the residuals do not display a specific pattern nor systematic departures from a horizontal line are seen as they are indicative of non-proportional hazards. In fact, proportional hazards assumes that estimates are time independent.

\label{fig:Schoenfeld} Scaled Schoenfeld residuals of covariates from the final model over time. For reference, the y=0 line is shown in red.

Scaled Schoenfeld residuals of covariates from the final model over time. For reference, the y=0 line is shown in red.

One can also assess model’s goodness of fit. By using R Cox package, one is able to plot de resiudal deviance of the model with respect to the observations. As a diagnostic plot, it helps confirm that since not particular pattern is experienced and that residuals are more or less symtrically distributed around zero, the model is of an acceptable fit.
\label{fig:Dev} Diagnostic plot : model residual deviance, the y=0 line is shown in red.

Diagnostic plot : model residual deviance, the y=0 line is shown in red.

At last, plotting the survival curves for the “final” model with respect to edema stratification (Fig. ) gives a better grasp on survival differences for patients with pbc.

\label{fig::survfit3} Survival curves of the Cox model for each strata (i.e. each edema category). 0=no edema, 0.5 = present without or resolved with diuritics, 1= present with diuritics

Survival curves of the Cox model for each strata (i.e. each edema category). 0=no edema, 0.5 = present without or resolved with diuritics, 1= present with diuritics

Another submodel without log tranformation of most covariates (all except bilirubin to pass Cox assumptions) was also tested. Chi-square difference test was performed to ensure that the non loged submodel was not statistically different from the final “log” model. However, it appeared that the loged transformed version of final model still had the lowest AIC.

Cox proportional hazards model - Discussions

This investigation aimed at analysing the survival of Primary Biliary Cirrhosis when treated or not with a D-penicillamine. Using Cox multivariate regression,no statistical evidence that these treatments differed in terms of survival time were found. Among the various registered risk factors in this dataset, only the age, bilirubin and albumin concentration as well as “dema seriousness were statistically significant and included in the regression model.In fact for Cox assumption validation,”dema was stratified and continuous variables were log transformed. More over, any other variables included in the model did not significantly improve the godness of fit but only increased model complexity.

According to the estimated coefficient for the sub model in Table , the final model from this project is displayed in Eq. :

\[\begin{equation} \label{eq::final} h(t) = h_{0}(t) \exp(1.590\log(\text{Age}))+\exp(−2.540\log(\text{Alb}))+\exp(0.721 \log(\text{Bili}))+\exp(3.190 \log(\text{Prot))} +\exp(0.385 \log(\text{Coop}) \end{equation}\]

where \(h\) and \(h_{0}\) are defined like in Eq. , Alb: Albuminen, Bili: Bilirubine, Prot : Prothrombin time, Coop : Urine cooper.

References

[1] Poupon and Raoul, “Primary biliary cirrhosis: A 2010 update,” J. Hepatol., vol. 52, no. 5, pp. 745–758, May 2010.

[2] “Primary biliary cholangitis symptoms and causes,” Mayo Clinic..

[3] P. Vlachos, “StatLib Datasets Archive.”.

[4] T. G. Clark, M. J. Bradburn, S. B. Love, and D. G. Altman, “Survival analysis part I: Basic concepts and first analyses,” British Journal of Cancer, vol. 89, no. 2, pp. 232–238, Jul. 2003.

[5] D. R. Cox, “Regression Models and Life-Tables,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 34, no. 2, pp. 187–220, 1972.

[6] T. Therneau, C. Crowson, and E. Atkinson, “Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model,” p. 27.

[7] Steve Simon, “The Proportional Hazard Assumption in Cox Regression,” The Analysis Factor. Aug-2018.