Original Article

Modeling Factors Associated with Dialysis Adequacy Using Longitudinal Data Analysis: Generalized Estimating Equation Versus Quadratic Inference Function

Background

Chronic kidney disease (CKD) is an irreversible and progressive deterioration of kidney function. The main treatment is kidney transplantation, but patients must undergo dialysis until finding a transplanted kidney. Each year, more than 60 000 people worldwide die as a result of kidney disease. The incidence of chronic kidney failure is 242 cases per million people globally and is increasing by 8% annually.¹ The annual increase of this disease is reported to be 10%-12% in Iran.² The high prevalence of this disease continued to cause high costs in the health care system worldwide.³ In a meta-analysis on the global burden of CKD, the highest CKD prevalence (14.44%) was found in the United States and Canada, followed by 12.10%, 11.86%, and 11.68% in Chile, Europe, and Iran, respectively, while the lowest rate of 6.76% was observed in India and Bangladesh.⁴

Biochemical markers are routinely measured to monitor the nutritional and health status of end-stage renal disease (ESRD) patients on dialysis, which is one of the alternative treatments for ESRD patients.⁵ The dialysis adequacy (DA) index is determined based on the difference between pre- and post-dialysis urea measurements compared to pre-dialysis measurements at each dialysis session. The measurement of DA in these patients is highly important. In the case of severe inadequate dialysis, when chronic kidney failure reaches the end stage, it is impossible for the patient to continue living without replacement therapy. However, kidney transplantation is impossible in all patients.⁶ Therefore, the adequacy of hemodialysis must always be reviewed in hemodialysis patients. Various methods are employed to determine this adequacy, which is influenced by the ratio of the urea before and after dialysis. The method proposed by Daugirdas is most widely used in clinical practice.⁷ Previous reports from the Islamic Republic of Iran have shown that DA is low compared to developing countries.⁸ Several factors such as demographic and clinical characteristics may affect DA. However, none of the published studies considered the longitudinal structure of DA data and its correlation structure in their analysis.

To monitor the health status of hemodialysis patients, the DA index and biochemical markers such as creatinine and albumin, urea, hemoglobin, hematocrit, blood lipids, and blood pressure are measured regularly. Changes in these markers during dialysis and related factors are of particular importance in predicting these changes in the nutritional assessment of patients and monitoring their health.⁹ These data are collected almost exclusively in a longitudinal setting. One of the appropriate models for the analysis of this type of data is the marginal mean model, whose best-known model is the first-order generalized estimating equation (GEE1).^10,11 In this model, the maximum likelihood estimation (MLE) considers the correlation structure a confounding parameter. Fortunately, one of the good features of the GEE model is that the coefficient estimates are consistent and unbiased even when the chosen correlation structure is misspecified.^10-13 The second-order GEE (GEE2) estimates the parameters based on the correlation structure using the MLE method, which can promise better estimation efficiency than GEE1. A recently developed method called the quadratic inference function (QIF), introduced by Qu et al, uses the second moments in the MLE approach.^14,15 The latter is more of a new methodology that provides a direct method for the goodness-of-fit criteria and consistency of parameter estimates and may be more efficient, and its inferences are likely to be more robust to outliers.¹⁴ However, this new method has rarely been utilized in statistical analysis in medical research. The objectives of this study were to demonstrate the use of GEE and QIF in modeling marginal means for longitudinal data analysis, determine factors associated with changes in DA, and compare the relative efficiency (RE) of the three methods, GEE1, GEE2, and QIF.

Methods

Design and subjects

This was a longitudinal study of a historical cohort in which longitudinal measurements of DA were recorded at regular intervals in each hemodialysis patient. The study sample consisted of 153 ESRD patients on dialysis between 2018 and 2019 at Shahid Beheshti Hospital in Babol, Iran. Adult patients who were over 18 years and under haemodialysis were included in the study.

Data and measurements

The demographic and clinical characteristics of 153 ESRD patients were recorded from patient records using a checklist. These included age, gender, and clinical history of underlying diseases such as diabetes, hypertension (HTN), lipids, cardiovascular disease, and smoking, age at starting dialysis (year), duration of dialysis history (year), duration of each dialysis session (hour), number of dialyses per week, and weight (kg). The DA index as a dependent variable was calculated longitudinally based on the pre- and post-dialysis urea measurements and the urea clearance parameter using the Daugirdas method.⁷ These measurements were routinely recorded at the end of each month during the follow-up period. During the two-year follow-up period, these measurements were taken from each patient’s records at equal intervals every three months or eight times.

Statistical analysis

In statistical analysis, first, the package MICE was used for the multiple imputations of missing data. This package performs multiple imputations using Gibbs sampling with full conditional specification for each variable with missing elements. Then, the data normality of DA was tested using the Kolmogorov-Smirnov test, and all data were satisfied with the normality assumption. GEE1, GEE2, and QIF models were applied to estimate the regression coefficients. In the current study, the association between DA as a function of time, age, gender, age at dialysis initiation, history of cardiovascular disease, HTN, cholesterol level, smoking, diabetes, dialysis duration, dialysis history, and weight as independent variables and DA as a quantitative dependent variable was estimated, and the marginal regression model was developed as follows:

The coefficients of the models were estimated and tested with different correlation patterns using GEE1, GEE2, and QIF models. In the current study, various correlation structures (interchangeable, autoregressive, and unstructured) were employed in the analysis of correlated data. The correlation structures must be determined before estimating the correlation parameters.¹⁶ In the interchangeable correlation structure, both responses of a person have the same correlation. In this case, the main diameter of the correlation matrix is equal to one, and the other elements of the correlation matrix are equal in size.¹⁶ In the autoregressive structure (first-order autoregressive), when the repeated observations within the individual depend on the previous observation, autoregressive correlation is applied to estimate the parameters.^15,16 This correlation structure is mostly used because only one correlation parameter is estimated in this case. In this structure, the correlation value between two measurements of the same person equals the base correlation value multiplied by the absolute value of the time interval of these two measurements. In the unstructured method, there are fewer constraints on correlation parameters. In other words, there are different correlations for both observations within the same individual. This type of structure is used when none of the systematic structures are applicable.^16-20 The choice of correlation structures depends on the opinion of the researcher, and the correlation structure that best reflects the relationship between correlations is not always defined clearly. For large samples, the standard error estimates of the parameters rely more on the correlation structure than on parameter estimates.²⁰

In fitting the GEE model, the Akaike information criterion is not appropriate for selecting models based on likelihood. The GEE model is based on the pseudo-likelihood; for the selection of the best performing correlation structure and comparison of the models, a statistic called the quasi-likelihood information criterion (QIC) criterion was employed, which contains the pseudo-likelihood information. The model fit is better and more appropriate when the value is lower.^20-23 Thus, the fit of GEE1 and GEE2 was assessed with the statistical fit index QIC, QICu, and the QIF was also compared with the statistical fit index Akaike information criterion (AIC) and Bayes information criterion (BIC). In addition, the efficiency of each of these models was compared using the mean square error index of the coefficients and the estimated RE. The mice software was utilized for missing data imputation. Additionally, the library of QIF, GEE, and geepack in R version 4.0.3 software and SPSS software were used for data analysis, and statistical tests were considered at a significance level of 0.05.

Results

Demographic and clinical characteristics

Of the 153 ESRD patients who underwent hemodialysis, 78 (51.3%) and 74 (48.7%) cases were men and women, respectively. The missing data of DA at 8-time points of follow-up were ranged from 0% to 13%. The mean (SD) of DA and the rate of dialysis inadequacy (KT/V < 1.2) were 1.16 ± 0.24 (KT/V) and 40.5%, respectively. The DA index (KT/V) was more desirable in women than in men () and in the age group < 60 years compared with ≥ 65 years (). The clinical characteristics of the study participants are presented in Table 1. HTN, cardiovascular disease, smoking, high cholesterol, and diabetes were present in 138 (90.2%), 72 (47.1%), 51 (33.3%), 73 (47.7%), and 16 (10.5%) subjects, respectively. The mean age (SD) and dialysis experience were 58.34 ± 14.19 and 4.07 ± 3.9 years, respectively.

Close

Figure 1. The estimated marginal means of dialysis adequacy according to gender (KT/V)

Close

Figure 2. The estimated marginal means of dialysis adequacy according to age group (KT/V)

Factors associated with dialysis adequacy

Table 2 provides the estimated coefficients (95% CI) of the regression parameters in the GEE1, GEE2, and QIF models. The correlation structure that produced the lowest AIC, BIC, and QIC fitting criteria was presented among various tested correlation structures. The three models produced rather similar estimates of the regression coefficients of factors associated with DA. For example, in the GEE1 model, age at starting dialysis (P = 0.003), HTN (P = 0.023), diabetes (P = 0.007), and weight (P = 0.001) had a negative significant effect, while female gender (P = 0.002) and duration of dialysis (P = 0.001) exerted a positive significant impact on DA. The coefficient of age, high cholesterol, smoking, history of cardiovascular disease, and number of dialyses per week had no significant influence on DA. In fitting GEE1 and GEE2, the QIC and QICu were lower for the unstructured correlation than for the other correlation structures. However, in fitting QIF, the best correlation structure was exchangeable with lower AIC and BIC. Thus, it was chosen as the appropriate structure.

Comparison of the relative efficiency quadratic inference function vs. generalized estimating equation

The RE of three presented models (i.e., GEE1, GEE2, and QIF) was calculated using the formula presented by Qu et al.¹⁴ The results revealed that the GEE2 model was more efficient than GEE1 (RE = 1.113), and the QIF model was more efficient than GEE1 and GEE2 (RE = 1.175 and RE = 1.056, respectively).

Discussion

The results of the GEE1 and GEE2 methods for the DA response showed that the estimates of the regression coefficients were quite similar for both methods, and there was no significant difference between the two methods, which is similar to the statistical significance of the parameter estimates. As for the RE of the coefficient estimates, the GEE2 method performed better than the GEE1 method. It is also known that the GEE1 model is used to estimate the marginal mean parameters, and the correlation parameters are not estimated in this model. This version of GEE may lead to inefficient estimates of the regression parameters because the MLE method considers correlation parameters as nuisance parameters.^23-25 However, the GEE2 model estimates not only the marginal mean parameters but also the dependence correlation parameters. This may enhance the reliability of the estimates of regression coefficients. However, in a study by Zayeri et al,²⁶ both GEE1 and GEE2 methods were applied to the continuous microleakage rate data. Further, the correlation between replicates was not significant, and the estimates of both GEE1 and GEE2 methods were almost similar.

The results of the ongoing study suggested that the QIF method was more efficient than GEE1 and GEE2. However, the estimates of the regression coefficients were similar to GEE1 and GEE2. The results of the current study are consistent with the findings of other studies.^23,25 In a simulation study conducted by Qu et al using Poisson data, it was found that the QIF method is more efficient than the GEE method, and when the correlation structure was correctly identified, the GEE and QIF methods had the same efficiency.²⁵ In the longitudinal study by Odueyungbo et al, the comparison between GEE and QIF demonstrated that the QIF method had better RE than GEE.²³ In most studies, the QIF method has been proposed as a method that eliminates the shortcomings of the GEE and improves the efficiency of model parameters when the correlation is not correctly determined. A few studies also discussed that the QIF method may not provide adequate results.²⁷

In the present study, clinically interesting findings indicated that the majority of patients had inadequacy of dialysis ( < 1.2 KT/V). When dialysis treatment was first used to replace the work of the kidneys, no one knew how much dialysis treatment was needed to keep patients healthy. Despite the remarkable improvement in DA in hemodialysis patients in developed countries, the majority of patients in developing countries, including Iran, still have a KT/V < 1.2. According to the 2006 NKF/DOQI criteria, the target value for KT/V is 1.4, and KT/V above 1.2 is considered an acceptable value for DA, which was used in the present study. Many factors contribute to the failure to achieve the minimum dialysis requirements, including frequent episodes of low blood pressure, poor access to health resources, incorrect nutritional status, blood flow rate, type of filtration, and also variations in sampling, and/or information errors.^27-31

In the current study, the quality of dialysis had a statistically significant relationship with patient gender (i.e., women were more satisfied with the adequacy of dialysis than men), which is consistent with the findings of a study by Taziki and Kashi in other dialysis centers in northern Iran.³² This finding may be explained by lower muscle mass, less physical activity, and diet compliance in women. The higher level of DA in women, compared with men, is also consistent with the results of other studies.^33-38

The findings of the ongoing study indicated that a history of HTN, diabetes, weight gain, dialysis history, and duration of dialysis were significantly related to DA. In another study on Iranian dialysis patients, a significant relationship was found between weight and DA (i.e., patients with higher weight had lower adequacy).³ In another study, the mean DA was low, and gender, smoking, dialysis history, number of dialysis sessions per week, and duration of dialysis also had a significant relationship with DA.⁴ In studies conducted between 1999 and 2004 in dialysis centers in the south, east, center, west, and north of Iran, 86%, 90%, 80%, 100%, and 58% of patients had a KT/V ≤ 1.2 as inadequate dialysis, respectively.^32-34 However, in the present study, dialysis inadequacy was 41%, indicating that the level of DA in Iranian patients has improved compared to two previous decades. Nevertheless, the level of DA in the current study was lower than in developed countries. Higher levels of DA have been reported in this regard. For example, Del Pozo et al³⁵ and Maduell et al in Spain reported a value of 1.3 KT/V and 0.98 KT/V, respectively.³⁶ In contrast, in a recent study by Saeedi et al in Iran, the majority of patients (74.8%) had no DA; the mean DA was extremely low, and gender, dialysis history, number of dialysis sessions per week, dialysis duration, and smoking history were significantly related to the KT/V value.³⁷ The difference in DA may be explained by the power of the rotation dialysis machine, duration of dialysis, the experience of the staff that operating the dialysis machine, and some patients’ characteristics.

Moreover, in the present study, the average dialysis time was 3.58 hours per session, which was lower than the average dialysis time in Taiwan (4.53 hours) and in all of Europe (5.4 hours), but it was almost as long as the average dialysis time in the United States (3.68 hours) and Germany (3.7 hours).³⁸ Although time is an influential and independent factor in the quality of dialysis, it is influenced by the number of patients, the alarm time of the dialysis machine, the efficiency of the machine, the qualification of the staff, and the rotation of the machine per minute. Increasing dialysis duration may improve the quality of dialysis to some extent. Some hemodialysis center staff are aware of this issue but may not pay attention to it due to the high number of patients on each shift and inadequate equipment and resources.

This study had limitations. Although the study setting was a central dialysis facility covering more than half a million people in northern Iran, this study may limit the generalizability of the results. A future prospective multicenter study of ESRD patients may overcome this limitation. However, this study has several advantages in terms of longitudinal analysis to include the correlational structure of longitudinal data of DA as a quantitative dependent variable, whereas the published study almost exclusively used a cross-sectional data analysis of DA with a conventional regression model.

Conclusion

The results of the present study indicated that the majority of hemodialysis patients had an inadequate dialysis value (KT/V ≤ 1.2). History of HTN, diabetes, dialysis duration, and weight had a negative impact on DA. Based on the findings, GEE1, GEE2, and QIF models yielded quite similar estimates of the regression coefficients, while the QIF model led to higher efficiency in estimating the coefficients compared with GEE1 and GEE2.

Acknowledgments

The authors would like to thank the Research and Technology Vice-chancellor of Babol University of Medical Sciences and the esteemed staff of the Dialysis Department of Shahid Beheshti Hospital, Babol, Iran for their sincere cooperation.

Author’s Contribution

Conceptualization: Karimollah Hajian-Tilaki, Khadigeh Gholian, Roghayeh Akbari.

Data curation: Khadigeh Gholian, Roghayeh Akbari.

Formal analysis: Karimollah Hajian-Tilaki, Khadigeh Gholian.

Investigation: Karimollah Hajian-Tilaki, Khadigeh Gholian, Roghayeh Akbari.

Methodology: Karimollah Hajian-Tilaki, Khadigeh Gholian.

Project administration: Khadigeh Gholian.

Resource: Karimollah Hajian-Tilaki.

Software: Khadigeh Golian.

Supervision: Karimollah Hajian-Tilaki.

Validation: Karimollah Hajian-Tilaki, Khadigeh Gholian.

Visualiztion: Karimollah Hajian-Tilaki, Khadigeh Gholian, Roghayeh Akbari.

Writing–original draft: Khadigeh Gholian.

Writing–review & editing: Karimollah Hajian-Tilaki, Khadigeh Gholian, Roghayeh Akbari.

Competing Interests

The authors declare that they have no conflict of interests.

Ethical Approval

This study was approved by the Ethics Committee of Babol University of Medical Sciences (Ethics code: IR.MUBABOL.REC.1400.110). Research involving human participants or human data was performed in accordance with the Declaration of Helsinki. All participants of the study were informed about the study, and they gave written informed consent to be included in the study.

Funding

Not applicable.

References

1. End-stage renal disease attributed to acute tubular necrosis in the United States, 2001-2010. Am J Nephrol 2015; 41(1):1-6. doi: 10.1159/000369832 [Crossref]
2. Prevalence of chronic kidney diseases and its determinants among Iranian adults: results of the first phase of Shahedieh cohort study. BMC Nephrol 2022; 23(1):203. doi: 10.1186/s12882-022-02832-5 [Crossref]
3. The cost of care for people with chronic kidney disease. Can J Kidney Health Dis 2019; 6:2054358119835521. doi: 10.1177/2054358119835521 [Crossref]
4. Global prevalence of chronic kidney disease–a systematic review and meta-analysis. PLoS One 2016; 11(7):e0158765. doi: 10.1371/journal.pone.0158765 [Crossref]
5. Long-term survival of living donor renal transplants: a single center study. Indian J Nephrol 2010; 20(4):179-84. doi: 10.4103/0971-4065.73439 [Crossref]
6. Renal replacement therapy in patients with acute renal failure: a systematic review. JAMA 2008; 299(7):793-805. doi: 10.1001/jama.299.7.793 [Crossref]
7. Theoretical basis for and improvement of Daugirdas’ second generation formula for single-pool Kt/V. Int J Artif Organs 2015; 38(12):632-7. doi: 10.5301/ijao.5000456 [Crossref]
8. Hojjat M. Hemodialysis adequacy in patients with chronic renal failure. Iran J Crit Care Nurs 2009;2(2):61-6. [Persian].
9. Prevalence of chronic kidney disease in population-based studies: systematic review. BMC Public Health 2008; 8:117. doi: 10.1186/1471-2458-8-117 [Crossref]
10. Longitudinal data analysis for discrete and continuous outcomes. Biometrics 1986; 42(1):121-30. doi: 10.2307/2531248 [Crossref]
11. Longitudinal data analysis using generalized linear models. Biometrika 1986; 73(1):13-22. doi: 10.1093/biomet/73.1.13 [Crossref]
12. Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of Longitudinal Data. Oxford University Press; 2002.
13. Fitzmaurice GM, Laird NM, Ware JH. Applied Longitudinal Analysis. John Wiley & Sons; 2012.
14. Improving generalized estimating equations using quadratic inference functions. Biometrika 2000; 87(4):823-36. doi: 10.1093/biomet/87.4.823 [Crossref]
15. Correlated binary regression using a quadratic exponential model. Biometrika 1990; 77(3):642-8. doi: 10.1093/biomet/77.3.642 [Crossref]
16. Structured correlation in models for clustered data. Stat Med 2006; 25(14):2450-68. doi: 10.1002/sim.2368 [Crossref]
17. Working‐correlation‐structure identification in generalized estimating equations. Stat Med 2009; 28(4):642-58. doi: 10.1002/sim.3489 [Crossref]
18. Working correlation structure misspecification, estimation and covariate design: implications for generalized estimating equations performance. Biometrika 2003; 90(1):29-41. doi: 10.1093/biomet/90.1.29 [Crossref]
19. Working covariance model selection for generalized estimating equations. Stat Med 2011; 30(26):3117-24. doi: 10.1002/sim.4300 [Crossref]
20. Akaike’s information criterion in generalized estimating equations. Biometrics 2001; 57(1):120-5. doi: 10.1111/j.0006-341X.2001.00120.x [Crossref]
21. Goodness‐of‐fit tests for GEE with correlated binary data. Scand J Stat 2002; 29(1):101-10. doi: 10.1111/1467-9469.00091 [Crossref]
22. Modeling strategies in longitudinal data analysis: covariate, variance function and correlation structure selection. Comput Stat Data Anal 2010; 54(12):3359-70. doi: 10.1016/j.csda.2009.11.006 [Crossref]
23. Comparison of generalized estimating equations and quadratic inference functions using data from the National Longitudinal Survey of Children and Youth (NLSCY) database. BMC Med Res Methodol 2008; 8(1):28. doi: 10.1186/1471-2288-8-28 [Crossref]
24. Hedeker D, Gibbons RD. Longitudinal Data Analysis. John Wiley & Sons; 2006.‎.
25. Dose of hemodialysis and survival: differences by race and sex. JAMA 1998; 280(20):1764-8. doi: 10.1001/jama.280.20.1764 [Crossref]
26. First and second order generalized estimating equations and their application in analyzing longitudinal microleakage data. J Iran Dent Assoc 2013; 25(1):6-15.
27. A note on improving quadratic inference functions using a linear shrinkage approach. Stat Probab Lett 2011; 81(3):438-45. doi: 10.1016/j.spl.2010.12.010 [Crossref]
28. Body size, dose of hemodialysis, and mortality. Am J Kidney Dis 2000; 35(1):80-8. doi: 10.1016/s0272-6386(00)70305-2 [Crossref]
29. Prescription of twice-weekly hemodialysis in the USA. Am J Nephrol 1999; 19(6):625-33. doi: 10.1159/000013533 [Crossref]
30. Quality of current management form of vascular access in hemodialysis patients. Rev Cubana Med 2012; 51(2):106-16.
31. Comparison of GEE1 and GEE2 estimation applied to clustered logistic regression. J Stat Comput Simul 2009; 79(4):361-78. doi: 10.1080/00949650701786085 [Crossref]
32. Taziki OB, Kashi Z. Determination of dialysis sufficiency in the patients referring to dialysis center of fatemeh zahrah hospital of sari in 2000. J Mazandaran Univ Med Sci 2002;13(41):40-6. [Persian].
33. Borzo R, Amin R. Assessment of velocity of blood flow effect on quality of dialysis in haemodialysis patients. J Shahrekord Univ Med Sci 2006;8(2):60-6. [Persian].
34. Mozaffari N, Mohammadi M, Dadkhah B, Mahdavi A. Assessment of quality of dialysis in Ardabil hemodialysis patients. J Ardabil Univ Med Sci 2005;4(14):52-7. [Persian].
35. [Accumulated experience in the analysis of quality indicators in a haemodialysis unit]. Nefrologia 2009; 29(1):42-52. doi: 10.3265/Nefrologia.2009.29.1.42.1.en.full.pdf [Crossref]
36. Maduell F, Vera M, Serra N, Collado S, Carrera M, Fernández A, et al. [Kt as control and follow-up of the dose at a hemodialysis unit]. Nefrologia 2008;28(1):43-7. [Spanish].
37. Saeedi M, Zareie F, Javaheri F. The assessment of dialysis adequacy and its related factors among hemodialysis patients. Nursing Development in Health 2016;6(3-4):23-31. [Persian].
38. Comparisons of Kt/V evaluated using an online method and calculated from urea measurements in patients on intermittent hemodialysis. Hemodial Int 2006; 10 Suppl 2:S5-9. doi: 10.1111/j.1542-4758.2006.00121.x [Crossref]