Introduction
Hemodialysis was first introduced in the 1960s as a practical treatment for kidney failure. Hemodialysis is a process in which blood is removed from the body of a patient with kidney failure and returned to the body after being filtered and cleaned through a dialyzer. The dialysis machine and the filter called artificial kidney can remove wastes from the blood or add necessary materials to it. By this method, the dialysis machine controls the acid-base balance, the amount of water and soluble substances in the body. Chronic kidney failure (CKD) is one of the most common diseases in Iran and the world. This disease is common in countries that have a high health status, due to increased life expectancy [2]. Hemodialysis improves the survival rate of patients with CKD, but reduces their quality of life [4]. Older people are more susceptible to CKD, which happens gradually  [5]. CKD is a general term used for various disorders that affect the structure and function of the kidney. The annual incidence of this disease in most countries is more than 200 cases per one million people [6]. Studies on examining the survival rate of patients with dialysis indicated better survival of patients in countries such as Japan and European countries [7]. A study on patients with end-stage renal disease has shown a premature aging phenotype in the immune system, which has recently been considered as an important factor for increasing the risk of various complications [9].
Methods
This is a retrospective cohort study. The study population includes all patients over 60 years of age with hemodialysis hospitalized in the hemodialysis department of Hashminejad Hospital in Tehran, Iran from 2004 to 2019. Of these, 395 eligible patients were selected using a census method. The inclusion criteria were age > 60 years and the completion of medical files. Their data were analyzed in RStudio versions 3.4.3 and 3.5.2 and SAS version 9.4 applications in two stages: longitudinal and survival [25]. The packages “Lattice”, “nlme”, “ggbiplot”, “survival”, “devtools”, “plyr”, and “scales” were used RStudio software. In the first step, a multivariate mixed model was used for all longitudinal biomarkers based on a random effects approach. However, in case of a large number of longitudinal biomarkers, a pairwise modeling approach was used, in which all possible pairs of bivariate-mixed models are fitted and combined in the final step. In the second step, the Cox proportional-hazards regression model was used corresponding to each longitudinal indicator with different constructs, including random effects prediction at time t, true unobserved value at time t, time-dependent slopes at time t, and cumulative effect.
Results
Mean, standard deviation (SD) and percentage were used to summarize the data. The demographic and clinical data of 395 elderly patients were examined. The follow-up period of the patients was 15 years. By using the two-stage model, the factors affecting the occurrence of death was investigated. The demographic characteristics of participants are provided in Table 1.


Figure 1 shows the graph of general changes in each of the study longitudinal biomarkers using Spline Smoothing for patients with hemodialysis in two groups of dead (occurrence of event) and not dead (censored) during hospitalization. As can be seen, the effect of time on longitudinal biomarkers is not the same and is different in both groups. The estimation of the coefficients of the longitudinal multivariate model for creatinine, calcium, phosphate and parathormone biomarkers showed that the first- and second-order forms of the time effect were simultaneously entered in the model (Table 2).


According to the results, with each year of age increase, the mean creatinine level decreases by 0.001 (P<0.001) 
Some longitudinal biomarkers had high correlations with each other in intercept and slope (Figure 2). These correlations should be accounted for in the longitudinal model using the multivariate normal distribution for all random effects. The level of correlation between the four longitudinal biomarkers is very important; the correlation in terms of intercept and slope between each pair of bivariate-mixed models were calculated.
Discussion
The purpose of this study was to use two-stage survival-longitudinal modeling to identify the factors affecting the survival of hemodialysis patients in order to help the therapists plan for better control of these patients. Previous studies have focused on the joint modeling for one longitudinal marker and time-to-event data. We used the two-stage modelling for survival and multiple longitudinal biomarkers. In the time-dependent model, we do not have measurement error; to solve this problem and speed up the calculations, we used a two-stage model. The difference between the two-stage model and the combined model is in the calculation of the likelihood function. 
Güler et al. introduced the extended JMLS model for multivariate longitudinal data and survival data based on the two-stage approach [25]. In 2016, Ossareh et al. conducted a study on the survival of hemodialysis patients and predicting mortality using a single-center analysis of time-dependent factors. They showed that the mortality rate was lower in men than in women. In patients over 75 years old, survival rate was 19.9% compared to 77.2% in patients aged <20 years. In addition to low life expectancy, many age-related factors such as malnutrition, inflammation, and heart failure, were suggested to contribute to decreased survival rate in older patients [7]. This is consistent with our results where older age, male gender and cardiovascular diseases were reported as the predictors of death in hemodialysis patients.
Data censoring is important to consider when having internal repeated measures in hemodialysis patients. Internal longitudinal measurements are taken when missing data occur. Censored patients lose longitudinal biomarkers after the occurrence of this event. For this reason, the proposed two-stage model is only guaranteed to provide valid results in this particular case, but it can generally be used for external longitudinal markers and time-to-event data where repeated measurements are taken before the follow-up.
5. Conclusion
The two-stage survival-longitudinal modelling enables us to study the complex relationship between all longitudinal biomarkers and survival in the hemodialysis program. 

Ethical Considerations
Compliance with ethical guidelines
The study was approved by the Research Committee of University of Social Health and Rehabilitation Sciences (Code: IR.USWR.REC.1398.100).

Funding
The study is taken from the master's thesis of Behnoush Asadollahi Dehkordi in the Department of Biostatistics and Epidemiology, Faculty of Rehabilitation Sciences, University of Social Health and Rehabilitation Sciences.

Authors' contributions
All authors contributed equally in preparing all parts of the research.

Conflicts of interest
The authors declared no conflict of interest.


References

1. Abbaszadeh A, Javanbakhtian R, Salehee S, Motvaseliyan M. [Comparative assessment of quality of life in hemodialysis and kidney transplant patients (Persian)]. Journal of Shahid Sadoughi University of Medical Sciences. 2010; 18(5):461-8. [Link]
2. Mahdavi-Mazdeh M, Heidary Rouchi A, Norouzi S, Aghighi M, Rajolani H, Ahrabi S. Renal replacement therapy in Iran. Urology Journal. 2007; 4(2):66-70. [Link]
3. Faucett CL, Thomas DC. Simultaneously modelling censored survival data and repeatedly measured covariates: A Gibbs sampling approach. Statistics in Medicine. 1996; 15(15):1663-85. [PMID]
4. Illescas DE, Reyes JR, Rosendo MC, y López LRL, Aguilera AA, Magali CNI, et al. Sun-089 decrease on quality of life and functional assessment in elder patients with end-stage chronic kidney disease under hemodialysis for two years. Kidney International Reports. 2019; 4(7):s192. [DOI:10.1016/j.ekir.2019.05.486]
5. Debone MC, Pedruncci EDSN, Candido MDCP, Marques S, Kusumota L. Nursing diagnosis in older adults with chronic kidney disease on hemodialysis. Revista Brasileira de Enfermagem. 2017; 70(4):800-5.  [PMID]
6. Levey AS, Coresh J. Chronic kidney disease. The lancet. 2012; 379(9811):165-80. [DOI:10.1016/S0140-6736(11)60178-5]
7. Ossareh S, Farrokhi F, Zebarjadi M. Survival of patients on hemodialysis and predictors of mortality: A single-centre analysis of time-dependent factors. Iranian Journal of Kidney Diseases. 2016; 10(6):369-80. [Link]
8. Chandrashekar A, Ramakrishnan S, Rangarajan D. Survival analysis of patients on maintenance hemodialysis. Indian Journal of Nephrology. 2014; 24(4):206-13.  [PMID] [PMCID]
9. Xiang F, Chen R, Cao X, Shen B, Chen X, Ding X , et al. Premature aging of circulating T cells predicts all-cause mortality in hemodialysis patients. BMC Nephrology. 2020; 21(1):271. [DOI:10.21203/rs.2.21570/v2]
10. Sepanlou SG, Barahimi H, Najafi I, Kamangar F, Poustchi H, Shakeri R, et al. Prevalence and determinants of chronic kidney disease in northeast of Iran: Results of the Golestan cohort study. PloS One. 2017; 12(5):e0176540. [PMID]
11. Meguid El Nahas A, Bello AK. Chronic kidney disease: The global challenge. The Lancet. 2005; 365(9456):331-40. [PMID]
12. Saran R, Robinson B, Abbott KC, Agodoa LYC, Bhave N, Bragg-Gresham J, et al. US renal data system 2017 annual data report: Epidemiology of kidney disease in the United States. American Journal of Kidney Diseases. 2018; 71(3 Suppl 1):A7. [PMID]
13. Hadian B, Anbari K, Heydari R. [Epidemiologic study of CDK and ESRD and their risk factors in dialysis patients of Lorestan province (Persian)]. Yafteh. 2015; 16(3):44-53.  [Link]
14. Otaghi M, Hassani P, Zagharie Tafreshi M, Nikbakht Nasrabadi A. [Challenge in acceptance of heamodialysis: A part of a grounded theory (Persian)]. Journal of Research Development in Nursing & Midwifery. 2013; 10(2):26-34. [Link]
15. Jones CA, McQuillan GM, Kusek JW, Eberhardt MS, Herman WH, Coresh J, et al. Serum creatinine levels in the US population: Third National Health and Nutrition Examination Survey. American Journal of Kidney Diseases. 1998; 32(6):992-9. [DOI:10.1016/S0272-6386(98)70074-5]
16. Parmar MS. Chronic renal disease. BMJ (Clinical Research ed.). 2002; 325(7355):85-90.  [PMID]
17. Brunner FP, Brynger H, Challah S, Fassbinder W, Geerlings W, Selwood NH, et al. Renal replacement therapy in patients with diabetic nephropathy, 1980-1985: Report from the European Dialysis and Transplant Association Registry. Nephrology, Dialysis, Transplantation: Official Publication of the European Dialysis and Transplant Association - European Renal Association. 1988; 3(5):585-95.  [PMID]
18. Haghighi AN, Broumand B, D’Amico M, Locatelli F, Ritz E. The epidemiology of end-stage renal disease in Iran in an international perspective. Nephrology, Dialysis, Transplantation: Official Publication of the European Dialysis and Transplant Association - European Renal Association. 2002; 17(1):28-32. [DOI:10.1093/ndt/17.1.28] [PMID]
19. Rizopoulos D. Joint models for longitudinal and time-to-event data: With applications in R. New York: Taylor & Francis; 2012. [DOI:10.1201/b12208]
20. Gilani N. Joint modeling of repeated ordinal measures and time to event data for CHD risk assessment. Biometrics & Biostatistics International Journal. 2019; 8(6):204-12. [DOI:10.15406/bbij.2019.08.00290]
21. Rizopoulos D, Ghosh P. A Bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event. Statistics in Medicine. 2011; 30(12):1366-80.  [PMID]
22. Henderson R, Diggle P, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics. 2000; 1(4):465-80. [DOI:10.1093/biostatistics/1.4.465] [PMID]
23. Wulfsohn MS, Tsiatis AA. A joint model for survival and longitudinal data measured with error. Biometrics. 1997; 53(1), 330–9.  [PMID]
24. Sattar A, Sinha SK. Joint modeling of longitudinal and survival data with a covariate subject to a limit of detection. Statistical Methods in Medical Research. 2019; 28(2):486-502. [DOI:10.1177/0962280217729573] [PMID]
25. Guler I, Faes C, Cadarso-Suárez C, Teixeira L, Rodrigues A, Mendonca D. Two-stage model for multivariate longitudinal and survival data with application to nephrology research. Biometrical Journal. Biometrische Zeitschrift. 2017; 59(6):1204-20.  [PMID]
26. Donnelly C, McFetridge LM, Marshall AH, Mitchell HJ. A two-stage approach to the joint analysis of longitudinal and survival data utilising the Coxian phase-type distribution. Statistical Methods in Medical Research. 2018; 27(12):3577-94. [PMID]
27. Brilleman SL, Moreno-Betancur M, Polkinghorne KR, McDonald SP, Crowther MJ, Thomson J, et al. Changes in body mass index and rates of death and transplant in hemodialysis patients: A latent class joint modeling approach. Epidemiology. 2019; 30(1):38-47.  [PMID]
28. Mauff K, Steyerberg E, Kardys I, Boersma E, Rizopoulos D. Joint models with multiple longitudinal outcomes and a time-to-event outcome: a corrected two-stage approach. Statistics and Computing. 2020; 30:999–1014. [DOI:10.1007/s11222-020-09927-9]
29. Murawska M, Rizopoulos D, Lesaffre E. A two-stage joint model for nonlinear longitudinal response and a time-to-event with application in transplantation studies. Journal of Probability and Statistics. 2012; 2012:1-19. [DOI:10.1155/2012/194194]