Abstracts

Introduction

Renal failure, also known as renal insufficiency and kidney failure, is a serious medical and public health problem, characterized by the inability of the kidneys to excrete toxic substances from the body in a proper way ¹1. Cardozo MT, Vieira IDO, Campanella LDA. Alterações nutricionais em pacientes renais crônicos em programa de hemodiálise. Rev Bras Nutr Cli´n 2006; 21:284-9.. There are numerous causes of kidney failure. Some of them lead to a rapid decline in renal function, with values less than 1 to 2% of normal levels, named acute renal failure (ARF). In addition, a gradual and progressive loss of a large part of the functioning nephrons can be observed, called chronic renal failure (CRF). According to Marques et al. ²2. Marques AB, Pereira DC, Ribeiro RCHM. Motivos e frequência de internação dos pacientes com IRC em tratamento hemodiali´tico. Arq Bras Ciênc Saúde 2005; 12:67-72., symptoms of kidney failure include headache, weakness, anorexia, nausea, vomiting, cramps, diarrhea, oliguria (insufficient secretion of urine), edema, confusion, thirst, loss of smell and taste, somnolence, hypertension and tendency to hemorrhage resulting from renal failure, as well as skin pallor, xerosis (pathological skin dryness), dysmenorrhea (cramping before or during menstruation), amenorrhea (absence of menses), testicular atrophy, impotence, attention deficit disorder, muscle spasms, and coma.

Patients who, for some reason, lost their kidney function and irretrievably reached the terminal stage of the disease, now have three different methods for treatment: hemodialysis, peritoneal dialysis and kidney transplantation. According to the Brazilian Society of Nephrology, hemodialysis promotes the removal of toxic substances, water and minerals from the body by filtering the blood. It is generally held three times per week in sessions lasting an average of three up to four hours, with the aid of a machine, in specialized clinics.

De-Lima et al. ³3. De-Lima JJG, da-Fonseca JA, Godoy AD. Dialysis, time and death: comparisons of two consecutive decades among patients treated at the same Brazilian dialysis center. Braz J Med Biol Res 1999; 32:289-95. and Andrade et al. ⁴4. Andrade LGM, Gabriel DP, Martin LC, Cruz AP, Balbi AL, Caramori JT, et al. Sobrevida em hemodiálise no hospital das cli´nicas da Faculdade de Medicina de Botucatu, Unesp: comparação entre a primeira e a segunda metades da década de 90. J Bras Nefrol 2005; 27:1-7. report that despite various technological innovations incorporated in hemodialysis, Brazilian’s studies did not show any great improvement in survival of patients with CRF in recent decades. The lack of increased survival of CRF carriers maintained on kidney dialysis could be explained by the increase in the average age of people who have started this treatment in recent years. Interest in the effects of age is justified by the recent growth in the number of elderly who start hemodialysis. According to Sims et al. ⁵5. Sims RJA, Cassidy MJD, Masud T. The increasing number of older patients with renal disease. BMJ 2003; 327:463-4., from the approximately 1 million people kept in this regular treatment in the world, more than half are older than 65 years.

The aging population and increasing lifetime expectancy, resulting from a demographic change in recent decades in Brazil, have contributed to changes in the morbidity and mortality profile and increased prevalence of chronic diseases, including chronic kidney disease (CKD) ⁶6. Grassman A, Gioberge S, Moeller S, Brown G. ESRD patients in 2004: global overview of patient numbers, treatment modalities and associated trends. Nephrol Dial Transplant 2005; 20:2587-93.. According to Bommer ⁷7. Bommer J. Prevalence and socio-economic aspects of chronic kidney disease. Nephrol Dial Transplant 2002; 17:8-12., hypertension and diabetes are the main risk factors for CKD and are becoming more frequent in the general population, contributing for the increased incidence of CKD. Oliveira et al. ⁸8. Oliveira MB, Romão JE, Zatz R. End-stage renal disease in Brazil: epidemiology, prevention, and treatment. Kidney Int 2005; 68:82-6.⁾ reported that in Brazil and Latin America about 15% of dialysis patients are diabetics.

In addition to age, hypertension and diabetes, other complication factors in patients undergoing hemodialysis are hepatitis B (HBS) and C (HCV). Despite the control of hepatitis transmission, it remains a huge problem in chronic hemodialysis units. Leão et al. ⁹9. Leão JR, Pace FHDL, Chebli JMF. Infecção pelo vírus da hepatite C em pacientes em hemodiálise: prevalência e fatores de risco. Arq Gastroenterol 2010; 47:28-34. and Ferreira et al. ¹⁰10. Ferreira RC, Teles AS, Dias MA, Tavares VR, Silva AS, Gomes SA, et al. Hepatitis B virus infection profile in hemodialysis patients in Central Brazil: prevalence, risk factors, and genotypes. Mem Inst Oswaldo Cruz 2006; 101:689-92. detected a statistically significant association between infection with hepatitis C and B and the duration of hemodialysis treatment.

In Brazil, there are insufficient studies that evaluate, at a national level, the survival and associated factors of patients in both dialysis modalities ¹¹11. Sancho LG, Dain S. Análise de custo-efetividade em relação as terapias renais substitutivas: como pensar estudos em relação a essas intervenções no Brasil. Cad Saúde Pública 2008; 24:1279-90.. Thus, regional and local studies are essential to identify such factors and to evaluate the survival of patients with CRF. Therefore, a survey was conducted at Maringá Kidney Institute to predict survival and identify associated factors in the lifetime of patients with renal insufficiency from the metropolitan region of Maringá, Paraná State, Brazil.

Methodology

Data set description

Data were collected from the Maringá Kidney Institute, for patients classified with CRF enrolled in a hemodialysis program from 1978 up to 2010. Due to the lack of specialized hospitals in the region, this institute is responsible for the majority of patients with renal failure in the metropolitan region of Maringá, which is formed by 26 cities and has more than 764,000 inhabitants.

The data set has a total of 177 patients (80 women and 97 men), of whom 23 underwent kidney transplant and 22 were diabetic. Most patients treated at the institute came from the Brazilian public health system and all patients in the survey were undergoing hemodialysis. The total of failure times (death) is 119, and 58 observations were considered censored (32% of censure) if the patient did not continue in the program for any reason or if patients did not die until the end of study. Moreover, hepatitis B and C antibodies were detected in 73 and 14 patients, respectively. Finally, average age of onset of treatment is 56 years old.

The variables considered in this survey were: (i) t _i: observed time (in days); (ii) δ _i: failure indicator (censored or observed); (iii) x _i1: age (in years) at the beginning of treatment; (iv) x _i2: gender (0-female, 1-male); (v) x _i3: marital status of the patient (living together, married, separated, divorced, single, widowed or others); (vi) x _i4: skin color indicator (yellow, white, black or pardo); (vii) x _i5: kidney transplant indicator (0-false, 1-true); (viii) x _i6: antibodies to hepatitis C (0-false, 1-true); (ix) x _i7: antibodies to hepatitis B (0-false, 1-true); (x) x _i8: diabetic indicator (0-false, 1-true), where i = 1,…,177.

Regression models

In practical applications, patients’ lifetimes are affected by explanatory variables, such as age, tumor size, lymph node status, among others, which can affect the healing probability of an individual. Thus, they should be added to statistical models to obtain better estimates as well as individual interpretations for these variables. Several surveys using the class of location models have been proposed in literature. A disadvantage of the class of location models is that the variance is not modelled explicitly regarding explanatory variables, but implicitly, through its dependence on the location parameter. Hence, an important consideration in regression models is the assumption that all observations have the same dispersion parameter. The noncompliance with this assumption may affect the efficiency of the estimates of the parameters.

As an alternative to the class of location models, Rigby & Stasinopoulos ¹²12. Rigby RA, Stasinopoulos DM. Generalized additive models for location, scale and shape. J R Stat Soc Ser C Appl Stat 2005; 54:507-54. proposed the generalized additive models for location, scale and shape (GAMLSS), in which the systematic part of the model is expanded, allowing not only the location but all the parameters of the conditional distribution of the response variable to be modelled as parametric and/or nonparametric functions of a set of explanatory variables.

Let T~D(t;θ), where D represents the distribution of T and θ ^T denotes the vector of its parameter. Consider independent observations t _i ’s conditional on the parameters θ _i (for i = 1,…,n) having probability density function (pdf) f(t _i ;θ _i ). To exemplify, consider a distribution having two parameters, say μ and σ, and define the elements of the vector θ ^T using two appropriate link functions,

where g _k (.), k = 1,2, denote appropriate link functions, β _k = (β _ok ,…,β _mk k) ^T is a parameter vector of length (m _k +1), m _k denotes the number of explanatory variables related to the k ^th parameter and Χ _k is a known model matrix of order nx(m _k +1). The choice of parameters to be modelled by explanatory variables will be discussed later in this study.

Next, consider a sample of n independent observations t ₁ ,…,t _n . Let c _i denote the censoring time, y _i = min(t _i ,c _i ) and δ _i = I(t _i ≤c _i ), where δ _i = 1 if t _i is a time-to-event and δ _i = 0 if it is right censored. From n observations, explanatory variables and censoring indicators (t ₁ ,δ _1, x _k1 ),...(t _n ,δ _n, x _kn ), the total log-likelihood function under non-informative censoring for the parameter vector θ = (β ₁ ^T _, β ₂ ^T ) ^T takes the form according to Cordeiro et al. ¹³13. Cordeiro GM, Ortega EM, Ramires TG. A new generalized Weibull family of distributions: mathematical properties and applications. J Stat Distrib Appl 2015; 2:13. (Section 13).

where S(t _i ;θ _i ) is the survival function relative to the distribution D, F and C denote the sets of individuals for which t _i (lifetime or censoring) and the vector of parameters are defined in (1) by specifying appropriate link functions. The numerical maximization of (2) can be easily performed using the GAMLSS package ¹⁴14. Stasinopoulos DM, Rigby RA. Generalized additive models for location, scale and shape (GAMLSS) in R. J Stat Softw 2007; 23:1-46. in the R software (The R Foundation for Statistical Computing, Vienna, Austria; http://www.r-project.org). Computational aspects are given in Figure 1.

Probabilistic distributions

In this survey, we use three different distributions as possible models to represent the lifetime of the patients. The simplest one-parameter distribution that may be appropriate to represent the response variable is the exponential distribution with E(T) = µ and Var(T) = µ ² , where the log link function is used to model µ with explanatory variables. Secondly, the log-normal distribution (two parameters) is also selected to explain the lifetime of the patients, with ET=exp⁡(μ+σ22) and Var(T) = [exp(σ ² )-1]exp(2µ+ σ ² ). In this case, we use the identity and log link functions to model µ and σ as functions of the explanatory variables, respectively. The last distribution used in this study is the Weibull distribution. We highlight here that the parameterization of this last model differs from the more common one generally used in the literature and is available in Stasinopoulos & Rigby ¹⁴14. Stasinopoulos DM, Rigby RA. Generalized additive models for location, scale and shape (GAMLSS) in R. J Stat Softw 2007; 23:1-46.. The main reason to use this parameterization is in fact that the parameter µ is exactly the mean of the distribution, making the interpretation of the results easier and simpler. The Weibull density is given by:

where βμ,σ= μ1σ+1, μ > 0 and σ > 0. The mean of T is given by E(T) = μ and the variance, VarT= μ2 2σ+1 21σ+1-1 where T(.) denotes the gamma function. The log link function is used to model μ and σ as a function of explanatory variables.

Model selection

Here, we consider the model selection process in three sequential steps. The first step consists in choosing the best distribution to represent lifetime and cure proportions. In the second step, we present a method to select the explanatory variables to fit the parameters of the model. Finally, in the third step, we investigated the model assumptions through a residual analysis.

Selecting the best model

In the first stage, the Akaike information criterion (AIC) and Bayesian information criterion (BIC) values are used to assess different fitted distributions. The AIC and BIC criteria are defined by AIC = GD+2k and BIC = GD+log(n)k, respectively, where GD represents the global deviance given by GD = -2l(θ), l(θ) is the maximized total log-likelihood function, n denotes the sample size and k denotes the number of fitted parameters. The model presenting the smallest values for these criteria is then selected. These criteria can be used to compare models disregarding regression structure as well as the class of location and heteroscedastic models. In general, the selected distribution to represent failure and censored times also provides a good fitted model in the presence of explanatory variables (see Cordeiro et al. ¹³13. Cordeiro GM, Ortega EM, Ramires TG. A new generalized Weibull family of distributions: mathematical properties and applications. J Stat Distrib Appl 2015; 2:13. and Ramires et al. ¹⁵15. Ramires TG, Cordeiro GM, Kattan MW, Hens N, Ortega EMM. Predicting the cure rate of breast cancer using a new regression model with four regression structures. Stat Methods Med Res, 2017; 1-17.). Moreover, an important function in survival analysis is the hazard function, which is defined as the event rate at time t conditional on survival until time t or later, that is, T ≥ t. In medicine, the hazard function is commonly called the force of mortality of the survival function. The empirical scaled total time on test (TTT) transform can be used to identify the shape of the hazard function, and then, we can verify if the chosen model was able to capture such shape.

Selecting the explanatory variables

The selection of the terms for all parameters is performed using a stepwise AIC procedure. Many different strategies could be applied for selection of the terms used to model all the distribution parameters. Here, we adopt the strategy described by Nakamura et al. ¹⁶16. Nakamura LR, Rigby RA, Stasinopoulos DM, Leandro RA, Villegas C, Pescim RR. Modelling location, scale and shape parameters of the Birnbaum-Saunders generalized t distribution. J Data Sci 2017; 2:221-38.. Basically, the strategy consists of two steps. In the first step, we adopt a forward selection procedure to choose an appropriate model for μ, with σ fitted as constants. Afterwards, we repeat the same procedure to select the model for σ, using the models already obtained in the previous steps as constants. For the second step, we perform a backward selection procedure to choose an appropriate model for μ, where σ is fitted as constant. At the end of the steps described, the final model may contain different terms for μ and σ. Note that if unnecessary terms are selected by the forward method in the first step, they will be removed in the backward method in the second step, thus the combination of these two methods is called stepwise AIC procedure. Details of this method can be found in Figure 1.

Residual diagnostics

To study departures from the error assumption and the presence of outlying observations, we can use the diagnostic tools in the GAMLSS package. The first technique is based on the normalized (randomized) quantile residuals ¹⁷17. Dunn PK, Smyth GK. Randomized quantile residuals. J Comput Graph Stat 1996; 5:236-44., which are given by r _i = Ф ^-1 (μ _i ), where Ф ^-1 () is the quantile function (qf) of the standard normal variable and μ _i = F(t _i ˅θ _i ), where F(t _i ˅θ _i ) denotes de cumulative density function. For censored response variables, μ is defined as a random value from an uniform distribution on the interval [F(t _i ˅θ _i ),1].

Although quantile residuals are widely used in the literature, it is not possible to identify specifically failures to fit the mean, variance, skewness and kurtosis in the response variable. As an alternative, we can use the worm plots (WP) ¹⁸18. van Buuren S, Fredriks M. Worm plot: a simple diagnostic device for modelling growth reference curves. Stat Med 2001; 20:1259-77.. These plots were pioneered by identifying regions (intervals) of an explanatory variable in which the model does not adequately fit the data. This is a diagnostic tool for checking the residuals for different ranges of one or two explanatory variables. In this study, we adopt an approach proposed by van Buuren & Fredriks ¹⁸18. van Buuren S, Fredriks M. Worm plot: a simple diagnostic device for modelling growth reference curves. Stat Med 2001; 20:1259-77., which is to fit cubic models to each of the detrended Q-Q plots with the resulting constant, linear, quadratic and cubic coefficients indicating differences (misfits) between the empirical and model residual mean, variance, skewness and kurtosis, respectively, within the range in the WP (for further details, see Stasinopoulos et al. ¹⁹19. Stasinopoulos MD, Rigby RA, Heller GZ, Voudouris V, De Bastiani F. Flexible regression and smoothing: using GAMLSS in R. London: Chapman and Hall/CRC; 2017.).

Results and discussion

We start the analysis by fitting exponential, Weibull and log-normal models disregarding regression variables. Table 1 shows maximum likelihood estimates (MLEs) and corresponding standard errors (SEs) (in parentheses) of the model parameters and GD, AIC and BIC statistic values for the fitted models. Note that we are considering the structure presented in (1), that is, for positive parameters we use the logarithmic link function μ = exp(β ₀ ). We also provide in Figure 2a the plots of estimated and empirical survival functions. Thus, we conclude that the Weibull distribution provides a good fit to these data.

The TTT plot showed in Figure 2b reveals that the hazard function has decreasing shape. On the other hand, Figure 2c shows that the estimated hazard functions of exponential, Weibull and log-normal models are constant, decreasing and unimodal (increasing and then decreasing), respectively, indicating that the Weibull distribution is appropriate to fit these data.

After selecting the Weibull distribution as the best one to fit these data, we propose a Weibull regression model. Using the steps described in the previous sections to select additive terms for the parameters, we obtain

and

For the selected explanatory variables to explain the parameters μ and σ from the Weibull model, we present in Figure 3 the empirical survival curves (Kaplan-Meier) stratified for the categorical variables. In these panels, we can see a great difference between the levels for each variable (gender, kidney transplant indicator, antibodies to hepatitis B and C), agreeing with the selected covariates in the final model.

Table 2 provides estimates, SEs, and p-values obtained from the fitted Weibull model based on GAMLSS. We highlight that all selected parameters are significant at the 10% significance level. Results in Table 2 reveal that the explanatory variables - age at the beginning of treatment, gender, kidney transplant indicator and antibodies to hepatitis B and C - are significant factors to explain the lifetime mean μ of the patients, which is expected from a wide literature that associates these variables as high-risk factors for patients (e.g. Tangri et al. ²⁰20. Tangri N, Grams ME, Levey AS, Coresh J, Appel LJ, Astor BC, et al. Multinational assessment of accuracy of equations for predicting risk of kidney failure: a meta-analysis. JAMA 2016; 315:164-74. and Fabrizi et al. ²¹21. Fabrizi F, Donato FM, Messa P. Association between hepatitis B virus and chronic kidney disease: a systematic review and meta-analysis. Ann Hepatol 2017; 16:21-47.). Moreover, the age at the beginning of treatment, kidney transplant indicator and antibodies to hepatitis B and C are significant factors to explain the variability of the shape parameter σ of survival times. A limitation of this study is that few explanatory variables were considered because only those variables were collected by the hospital. Further research, including other variables, may reveal factors not yet discovered in the literature.

Partial effects of the explanatory variables in the mean μ are displayed in Figure 4. From the model for the mean μ, we can see that the older the patient when he starts treatment, the lower is their lifetime expectancy (panel (a)). More precisely, for each additional year at the beginning of treatment, the patient’s average lifetime expectancy decreases 0.036 years (Table 2). The lifetime expectancy is lower in female patients (panel (b)) and in patients who have not undergone kidney transplantation (panel (c)). Moreover, panels (d) and (e) indicate that patients who have antibodies to hepatitis C and B, respectively, have greater lifetime expectancy.

Regarding the shape parameter σ (Figure 5), as we can see in panel (a), the variable age at the beginning of treatment has a positive linear relationship with the variability of patients’ lifetime , i.e. the variability of lifetime increases when patients waited to start treatment. Furthermore, according to panels (b)-(d), the variability also increases when the patient has undergone kidney transplantation and when there are antibodies to hepatitis B and/or C. Here, it is clear that the variability of patient’s lifetime must be modelled as well, and not only the mean (as more traditional models such as classical linear regression or generalized linear model usually do), which could mislead us on the selection of false risk factors for kidney failure.

Figure 6 shows some residual plots that can help to verify the adequacy and the assumptions of the chosen fitted model given in (3) and (4). Panel (a) suggests that the normalized quantile residuals have approximately normal distribution. Panel (b) shows that there are few points off the line in the low end of the range. Finally, as we can see in panel (c), the WP also suggests that there is a slightly problem in the lower tail of the distribution of T. The correct would be the fitted cubic model (red curve) to be a line on the X-axis. Nevertheless, the Weibull regression model based on the GAMLSS framework provides a reasonable fit to these data.

Finally, using the results presented in Table 2 we can calculate the average lifetime for each level of the explanatory variable selected in the regression model. Averages are easily obtained using the estimated linear structure:

Consider the high-risk group (x ₁ = maxx ₁ = 88, x ₂ = 0, x ₅ = 0, x ₆ = 0 and x ₇ = 0) and the low-risk group (x ₁ = minx ₁ = 17, x ₂ = 1, x ₅ = 1, x ₆ = 1 and x ₇ = 1). The average of the lifetime for high-risk and low-risk groups is 310 and 21,568 days, respectively, which correspond to 0.84 and 59 years, respectively. The same idea can be used to estimate the average lifetime for any level of explanatory variables.

Concluding remarks

Renal insufficiency is a serious medical and public health problem. Thus, regional studies are needed to identify associated factors that affect patients’ lifetime. Based on the fit of the Weibull model, we verify that the late onset of renal insufficiency, female gender, absence of kidney transplant and no antibodies to hepatitis B and C can be identified as negative factors that affect the lifetime of patients from the metropolitan area of Maringá. Hence, patients having such characteristics require special attention. Furthermore, the age at the beginning of treatment, kidney transplant indicator and antibodies to hepatitis B and C are significant factors to explain the variability of the survival time.

References

Publication Dates

History