Survival of hemodialysis patients: modeling differences in risk of dialysis centers
Abstract
Objective. Dialysis is the most common renal replacement therapy for patients with end stage renal disease. This paper considers survival of dialysis patients, aiming to assess quality of renal replacement therapy at dialysis centers in Rio de Janeiro, Brazil, and to investigate differences in survival between health facilities.
Methods. A Cox proportional hazards model, allowing for timevarying covariates and prevalent data, was the basic method used to analyze the survival of 11 579 patients on hemodialysis in 67 health facilities in Rio de Janeiro State from January 1998 until August 2001, using data obtained from routine information systems. A frailty random effects model was applied to investigate differences in mortality between health centers not explained by measured characteristics.
Results. The individual variables associated with the outcome were age and underlying disease, with diabetes being the main isolated risk factor. Considering covariates of the health unit, two factors were associated with performance: bigger units had on average better survival times than smaller ones and units which offered cyclic peritoneal dialysis performed less well than those that did not. There were significant frailty effects among centers, with relative risks varying between 0.24 and 3.15, and an estimated variance of 0.43.
Conclusions. Routine assessment based on health registries of the outcome of any high technology medical treatment is extremely important in maintaining quality of care and in estimating the impact of changes in therapies, units, and patient profiles. The frailty model allowed estimation of variation in risk between centers not attributable to any measured covariates. This can be used to guide more specific investigation and changes in health policies related to renal transplant therapies.
 health services assessment
 renal replacement therapies
 survival frailty models
End stage renal disease (ESDR) is highly prevalent worldwide, with around 1 million patients undergoing renal replacement therapies. Dialysis is still the most common treatment, partly due to insufficient kidney donation to meet the demand for renal transplants. The mortality rate is one of the most important outcome measures and is greatly influenced by individual risk factors such as the type of disease causing ESDR, age, and comorbidities. The use of registry data to analyze patterns and outcome of care of ESDR in order to improve quality of dialysis facilities is stimulated through the International Federation of Renal Registries [1]. As there is no ‘best’ treatment suited for all patients, and improvement in medical technology is a continuous process, the development of methods to be used for quality assessment based on medical registries is of major concern.
The mortality of patients undergoing dialysis is a main primary outcome variable. It can vary according to treatment and individual patient characteristics, and between health facilities [2]. Some characteristics of health facilities, such as the size of the unit or the number of patients assisted, can affect the outcome [3,4]. Recommended procedures in hemodialysis, related to adequacy of water sources and filtering, infection control, and handling of hepatitis B serumpositive patients, should be routinely assessed. Socioeconomic and cultural profiles may also impact strongly on survival due to differences in coping with treatment and side effects. Other known factors related to outcome are nutritional status [5], dialysis dose [6], late nephrology referral [7], and the use of auxiliary medicine.
In Brazil, more than 95% of renal replacement therapies are funded and controlled by the public health system (SUS). Payment for services delivered to patients is based on a computerized information system that collects data on all procedures on every patient on a monthly basis. Most ESDR patients undergo hemodialysis, with <10% undergoing peritoneal dialysis. This choice is possibly related to local medical culture, as there is no evidence of impact on survival of either treatment [8]. The patient is allocated to a specific dialysis center (DC) for a variety of reasons. Proximity to residence is more important where there is only one DC in the same municipality or health district, but irrelevant in larger cities with more than one DC. Place of previous treatment of underlying disease can influence the choice of DC. Physicians’ suggestions and individual choices are taken into account, as long as the desired center can accept the patient. As all DCs are at least partially public funded, they should all have roughly the same resources, including solutions, medicines, and equipment. The use of auxiliary medicines and surgical interventions related to ESRD is standardized [9]. Data on renal replacement therapies available through this information system have not yet been used for evaluation of quality of care delivered. Recently, the Brazilian Health Ministry released a computer system specifically designed to be used by local level health managers to analyze basic data on high cost ambulatory procedures. However, for a more careful analysis of the outcome, more sophisticated tools must be used [10]. Differences between treatment centers certainly play a role in dialysis adequacy [6,11]. Variation often cannot be fully explained by differences in registered clinical procedures, observed patient profile, or other specific and known conditions. Hence, a shared frailty or random effect approach to investigating betweencenter variation after allowance for measurable risk factors is proposed. The frailty model used in this paper is an extension of the Cox proportional hazards model, which allows for independent random effects, one for each DC, as a way of modeling unobserved sources of variation in survival of patients sharing the same hemodialysis facilities. The use of random effects modeling in epidemiology has increased in recent years [12], although the introduction of such modeling into eventtime analysis has proceeded more slowly. This approach to investigation of outcome related to the use of health services can help uncover aspects not yet studied and help formulate health policy guidelines.
Methods
Data
Data on patients using all DCs in Rio de Janeiro State were obtained from the High Complexity Procedures Information System (APAC), on CDROM. This followup system began in January 1998, recovering information on all patients using hemodialysis at that time, including the date treatment began. There is no information on patients that died before this date. The period covered by this study includes all patients on hemodialysis until August 2001. Only hemodialysis survival was studied, as alternative treatments are rarely used in Brazil. The data set includes 11 579 patients, receiving hemodialysis in 67 DCs. During the followup period, 2550 patients changed their place of treatment. Some patients were lost to followup before the closing date of the study. The reasons for this additional Censoring of some survival times were moving to another state or abandoning treatment (10%), kidney transplant (9%), and change to peritoneal dialysis (2.5%).
The outcome is the time in months since the beginning of hemodialysis until death or the end of followup. Individual variables included in final model were: (1) age of patient at onset of disease, as a continuous variable in models; categorized in three groups (≤42, between 43 and 56, ≥57) in Kaplan–Meier plots; and (2) underlying disease, registered in the system using the 10th revision of the International Classification of Diseases, and categorized as hypertension, diabetes, primary kidney diseases, birth defects, or other causes. Hypertension was taken as baseline.
Gender was not included in the study described here, as in the preliminary analysis of the data it was found not to affect survival. There is no comorbidity, socioeconomic, or disease severity information in the data set.
DC characteristics available were:
• number of dialysis machines, grouped in three classes: (1) from three to 10; (2) between 11 and 19; and (3) ≥20, up to 50 dialysis machines in the DC; and
• other dialysis modalities offered by the units: 23 DCs offer both ambulatory peritoneal dialysis (APD) and cyclic peritoneal dialysis (CPD), 12 offer APD only, five offer CPD only, and 27 offer no alternative treatment. Note that this variable measures other possibilities of treatment offered by the unit, but not the actual individual treatment. If the unit offers peritoneal dialysis, patients can change from peritoneal dialysis to hemodialysis, and viceversa. Alternative treatments in the same unit can then be a reason for increasing censoring, as the followup stops when patients change treatment. For reference, change from peritoneal dialysis to hemodialysis is much more common than change in the other direction [8].
Other variables investigated, but not included in the final model as they showed no relation with survival time, were: total number of rooms, type of machine and type of water filters used, hierarchical level of the DC in the state health system, and presence of separated rooms for hepatitis Bpositive patients. Variables related to the patient profile within DCs, such as proportions of diabetic patients and proportions of patients in a queue for kidney transplantation, showed no significant contributions to the model. There were no detectable interactions between age and the other covariates.
Method
The basic method, Cox proportional hazards model [13], is defined as a semiparametric model because a parametric form is assumed for the covariate effects, whilst the baseline hazard function is treated nonparametrically. It is also called a proportional hazards model because the ratio of hazard rates of any two individuals is proportional and constant over time. The assumption of constant relative risk over time was checked using Schoenfeld’s residual analysis [14].
The data set comprises a combination of two types of data: ‘incident’, for patients beginning treatment after January 1998, and ‘prevalent’ for patients undergoing treatment at that time. A further complication is that some patients changed DC during their course of treatment. Fortunately, both of these analytical difficulties can be overcome using a Cox proportional hazards analysis with proper regard to risk set definition, using procedures available in the public domain statistical package R [15].
Let α(t/x) be the hazard rate at time t for an individual with covariate vector x. The Cox model α(t/z) = α_{0}(t) exp(βx), where α_{0}(t) is an arbitrary baseline hazard rate, β is a parameter vector, and, thus, the logarithm of α(t/x)/α_{0}(t) is in the spirit of the usual linear models formulation for the effects of covariates. The model compares the time until death or censoring of each subject with an ‘at risk’ population under observation with the same survival time. The function exp(βx) for the failure is then compared with corresponding values for patients in the risk set. This makes it easy to deal with timedependent covariates such as DC in use: we simply allocate the patients in the risk set to the DC in use at the time of interest. In effect, this divides the total time until death or followup for each patient into a number of subintervals, corresponding to periods at risk in the different DCs, ending with induced censoring when the patient changes DC, except of course for the very last interval which may be terminated by death. The same basic idea of risk set adjustment is used to handle the prevalent part of the data set. Each patient in this cohort is included in risk sets only for failure times that are longer than the period between the beginning of their treatment and January 1998. This then properly accounts for these patients not being at risk of ‘observed’ failure during the preJanuary 1998 period. Note that failure to adjust for the prevalent cohort data would lead to an upward bias in survival rates, as highrisk patients with early failure would be underrepresented in the data. As mentioned, this adjustment is incorporated in the R routine.
One possibility for analyzing the effect of the health facility is to include it as a dummy variable, estimating the effect of each one compared with some chosen baseline unit. This approach, although feasible for a small number of groups, does not allow estimation of the parameters of unit level covariates, as they are the same for all patients in each health unit. Hence, either the effect of covariates can be estimated or a single centerspecific effect parameter can be estimated, but not both simultaneously. An alternative approach, which allows for unit level covariates in addition to unit level effects, is the socalled random effects model. These models provide a flexible and powerful tool for the analysis of grouped data, which can be useful in describing heterogeneity and/or correlation among individuals. In the survival analysis context, random effects are assumed to act multiplicatively on the baseline hazard. This random effect can be thought of as a ‘frailty’, increasing group or individual susceptibility to death when it is large, and decreasing susceptibility when it is small.
Formalizing the frailty model, usually assumed to be an extension of the Cox proportional hazards model, let z be a random variable representing an unknown random effect, related to dialysis center, with unit mean and variance ε. Large values of ε reflect a greater degree of heterogeneity among DCs. For an individual with regressor variables vector (x) and frailty effect (z), the form often assumed for the hazard rate is α_{0}(t/z) = zα_{0}(t) exp(βx). The standard assumption is that individuals within each cluster share a common frailty effect, z, in our case the clusters being the DCs.
Individuals or groups with frailty values >1 tend to experience the event at a faster rate than under the basic Cox model, indicating poor performance. Conversely, when frailty values are <1, survival times tend to be longer. Note that when the variance of z approaches zero, the model reduces to the basic Cox model. The most common model proposed in the literature for the frailty effect is the oneparameter gamma distribution [16]. One attractive feature of the gamma distribution is that it is mathematically tractable. Computationally, frailties can be viewed as an unobserved covariate, estimated through the estimation–maximization algorithm [17]. An alternative approach, which is less computationally demanding, is to treat the frailty model as a penalized Cox model [18]. This algorithm is available in the public domain statistical package R, which was used for the analysis.
Results
Kaplan–Meier plots illustrating the effect of each variable are shown in Figure 1, with group definitions and number of observations given in the legend. Logrank tests for survival differences were all highly significant (P < 10^{−6}), due in a large part to the size of the data set. The main known risk factor for survival in ESRD patients is diabetes, as can be seen from the first plot. Age is important too, with a poorer prognosis for older patients. The DC variables show a protective effect of larger units, although we of course recall that this variable is to be considered as a proxy for the likely expertise available in the DC; simply adding machinery/equipment would not in itself affect prognosis. Hemodialysis survival is lower in the centers that offer only cyclic peritoneal dialysis as an alternative treatment. However, only five of the 67 DCs fall into this category, with 1139 patients (out of 11 579 in total), and so this result should be interpreted with caution.
Table 1 shows the relative risk for each covariate in survival models with and without frailty. Age at onset of disease is highly significant: each additional year leads to an increase in risk of 3%. As expected, diabetes is the most serious underlying disease, increasing risk by >50% in comparison with hypertension. Birth defect causes for kidney failure present the best prognosis, as in this group the chances for kidney transplant are greater. Primary kidney disease does not affect survival. Under the label of ‘other causes’, a number of different pathologies, from cancer to autoimmune diseases, which are all too rare to be analyzed individually, were grouped. Parameter estimates of individual covariates are quite similar for both models.
Covariates  No frailty  Frailty  
Relative  Relative  
risk  Lower 0.95  Upper 0.95  risk  Lower 0.95  Upper 0.95  
Individual  
Age at onset  1.034  1.031  1.036  1.035  1.033  1.038 
Diabetes  1.536  1.404  1.681  1.554  1.417  1.705 
Primary kidney disease  0.922  0.840  1.012  0.930  0.844  1.025 
Birth defect  0.616  0.477  0.797  0.646  0.499  0.836 
Other causes of ESRD  1.005  0.902  1.120  1.086  0.959  1.229 
Dialysis centers  
Number of machines  
11–19  0.912  0.837  0.994  0.954  0.648  1.403 
≥20  0.595  0.540  0.656  0.598  0.373  0.960 
Cyclic PD  2.371  2.047  2.747  2.359  1.153  4.828 
Ambulatory PD  1.333  1.196  1.486  1.541  0.940  2.526 
Both types of PD  0.832  0.747  0.926  0.863  0.542  1.374 
Marginal likelihood  −26 199.30  −25 870.32  
Frailty variance  0.425 

PD, peritoneal dialysis.
With respect to DC covariates, the unit size showed a strong protective effect, with a consistent trend: large units (10–19 machines) present a slightly smaller risk, while those with >20 machines show a reduction of 40%. Providing cyclic peritoneal dialysis as an alternative treatment showed the highest risk. Patients receiving hemodialysis in one of the five centers where cyclic dialysis was the only alternative treatment had a relative risk that was 2.3 times greater than those using DCs where no other treatment was offered. Offering only ambulatory peritoneal dialysis has also been found to be a risk factor, although only in the model without frailty effect, while offering both types of hemodialysis presents a protective effect in the same model.
When DC frailties are included in the model there is a significant increase in the marginal likelihood, estimated after integrating out the frailty effect. For just one degree of freedom in the model, the likelihood increased from −26 199.30 to −25 870.32. The estimated size of the DC covariate coefficients increased under the frailty effect model, as expected theoretically [19], but only by a small amount. At the same time the inclusion of the random effect broadens the confidence intervals, and three DC level covariates lost significance.
The estimated frailty variance (0.43) is significant and indicates the presence of a large variability among DCs. The relative risks of units estimated through the frailty effect varied from 0.235 and 3.146, with four DCs with less than half the relative risk and two units with more than double the risk (Figure 2). Based on the confidence interval of relative risk of each unit frailty, the DCs were classified into three groups: frailty significantly above one, frailty significantly below one, and the middle risk group. A Kaplan–Meier plot of these groups (Figure 3) displays risk profiles compatible with estimated frailties. The survival curves of four DCs with the highest risk, together with the mean survival in the significantly highrisk group show that the model is adequately capturing highrisk units.
Discussion
Assessment of outcome in any high technology medical treatment, such as renal replacement therapy, is extremely important in order to maintain quality of care. Furthermore, it allows estimation of the impact of new therapies and technology introduced in uncontrolled, normaluse conditions. For this task, the data sets have to be the routine ones, with all concomitant problems. Exploratory data analysis suggested that the quality of the information set was adequate for the purpose of this study. Moreover, data are updated on a monthly basis and we made available for analysis within 3 months, thus allowing regular and uptodate assessment.
Many variables included in the system, however, did not seem to be associated with patient survival. For instance, all variables related to the DC, with respect to equipment, physical installations, and filters, showed no relation to survival. Although this kind of information is useful at the time of the implementation of routines, collecting the same data on unit structure every month may be inefficient. On the other hand, the system should be sufficiently flexible to include new variables as soon as they are implemented in any service, allowing follow up of any subsequent changes in patient survival.
The individual variables included in the model behaved as expected. Other important factors not available, however, are the nutritional status at the beginning of dialysis [5] and the socioeconomic status of the patient. The latter could help to estimate problems in coping with the renal replacement therapy.
One important methodological issue tackled in this paper is the simultaneous use of prevalent and incident data, which is normal in observational studies. The technique to deal with it was easily implemented using R, as well as the use of a timedependent covariate.
The protective effect of larger DCs found in other studies [3,20] should be investigated in detail, as after controlling for the main covariates (age and underlying disease), larger units still had a reduced risk of approximately 40%. There was no evidence of patient profiles at DCs affecting survival.
It was not possible to compare our results related to the offer of dialysis modalities, because most DCs worldwide offer both kinds of dialysis. As only hemodialysis patients were included in the study, and followup starts with the actual beginning of hemodialysis, not with ESRD diagnosis, patients undergoing CPD may have moved from this therapy to hemodialysis, and thus had prolonged ESRD. This longer time of disease could affect observed survival during hemodialysis. However, this would also happen with ambulatory peritoneal dialysis and the offer of this modality in the DC did not affect survival as strongly, since this covariate lost significance in the frailty model. Moreover, CPD is not a first choice for treatment, and is used only when patients cannot undergo hemodialysis (usually three times a week) or cannot be trained to use peritoneal dialysis at home. CPD may also be a temporary resource when arterial fistula access is lost. Centers offering only CPD, and not APD, may be technologically impaired. A longer followup time is needed to evaluate this finding, as peritoneal dialysis is still little used in Rio: in the last month of the followup period, only 6% of patients were undergoing PD. However, these findings have to be seen in the context of assessing survival, quality of life, and global costs of both treatments. It should be noted that the choice of treatment modality is strongly influenced by the provider’s perspective, such as whether they choose high cost investments in expensive hemodialysis facilities as opposed to low fixed costs and high expenditures in peritoneal dialysis supplies [21].
The frailty model allowed identification of variation in risk not attributable to any measured covariates, and showed a large amount of variability between DCs with respect to performance. One possible explanation for this variability could be attributed to differences in the socioeconomic level of the population using different services, due to location. As there is no variable for measuring socioeconomic status, a possible way to incorporate this could be through the investigation of a geographical pattern in the frailty of DCs [22]. Preliminary analysis showed no visible pattern, although more investigation is needed.
Other possible explanations for significant variation in DC survival rates are in different patient mix or unregistered differences in care practices. The first one could be explored in a model still under development, including individual frailty besides the aggregated level one. The variations in care have to be investigated locally in order to uncover possible differences and to propose new variables to capture these differences. This is one important advantage of the frailty model used here: detecting unmeasured variation leading to differences in risk among units. The other possible way to model the influence of each DC in risk, through including a dummy variable for it in the model, does not allow us to include in the model unit level covariates, as they are the same for all patients in the DC. So if DC relative risk was estimated in this way we could not at the same time determine the effect of variables such as unit size or other treatments offered.
The focus of this study is specifically hemodialysis patients, because peritoneal dialysis therapies, though beginning to be used more widely throughout the state, are still little used in the region. However, as the information system contains all data on the trajectory of each patient, including changes of residence, unit, therapies, and related surgical interventions, it will be possible to devise a comprehensive regular assessment of high cost ambulatory therapies. Studying this kind of longitudinal data, in spite of methodological problems, would allow better understanding of the entire process.
Although there is no doubt that controlled clinical trials represent the gold standard with which to assess efficacy of any treatment under everyday conditions, in the prolonged treatment of a chronic disease, the use of registry data is critical in comparing different centers, patients, and places. The frailty model proposed in this paper can be useful for this purpose. Future developments are the incorporation of nested frailties, treating simultaneously the random effects of services and patients, and time varying frailty, to assess variation in performance over time.
In the Brazilian context, more investigation should be undertaken in longitudinal analysis of the entire history of patients undergoing renal replacement therapies, as well as devising ways to regularly assess all complex ambulatory procedures, with data available in this information system. Finally, local investigation of highrisk DCs should be undertaken in order to disclose factors leading to the differences in risk detected in this work.
Acknowledgments
The revision of this manuscript has benefited from comments from the editor and two anonymous referees. The research in this paper was funded by the Brazilian research agencies CNPq and CAPES, through grants to M.S.C. (no. 200509/004) and S.S. (no. BEX1139/967), and by the Portuguese agency Fundação para a Ciência e Tecnologia, providing a master studentship to I.P.S.C.
Footnotes

Address reprint requests to Dr Marilia Sá Carvalho, Departamento de Epidemiologia, Escola Nacional de Saúde Pública, Rua Leopoldo Bulhões, 1480/8° Andar, Rio de Janeiro, RJ 210041210, Brazil. Email: carvalho{at}procc.fiocruz.br
 International Journal for Quality in Health Care 15(3) © International Society for Quality in Health Care and Oxford University Press 2003; all rights reserved