Original Article

Modeling Time to Blindness of Glaucoma Patients: A Case Study at Jimma University Medical Center

Background

Glaucoma is one of the most common causes of blindness and patterns of visual field loss due to retinal ganglion cell degeneration.¹ Risk factors for the two main types of glaucoma, (i.e., open-angle and closed-angle glaucoma disease) with different patterns of disease occurrence include increasing age, high intraocular pressure (IOP), family history of glaucoma, race and ethnicity, diabetes, as well as female gender.^2-13Glaucoma imposes a substantial burden on society in terms of higher medical costs, lost productivity, patient morbidity, and the number of ophthalmic consultations¹⁴and results in different psychological problems, including depression, inferiority complex, anxiety, and repudiation due to the feeling of low self-esteem.¹⁵

Glaucoma is a major public health problem that accounts for 8% of the world’s irreversible blindness and is the second largest cause of blindness following cataracts.¹⁶ An estimated number of people with glaucoma is 60.5 million worldwide, about 13.9% of whom had become blind due to glaucoma.¹⁶ Moreover, the number of people with glaucoma is expected to escalate to 111.8 million by 2040.^1,17 In Africa, glaucoma is responsible for 15% of the world’s blindness.¹ In Ethiopia, an estimated 62 thousand people suffer from irreversible sight loss caused by glaucoma, as the fifth most common cause of blindness.^18,19

There is a dearth of published studies on glaucoma in Ethiopia, particularly in the study area with problems in routine health services. Therefore, this study aimed to provide a modeling of time to blindness of the left and right eyes of glaucoma patients considering the different factors that affect the time to the blindness in glaucoma patients. The time to the blindness of the left and right eyes in glaucoma patients can be predicted and statistically estimated using the bivariate survival model.²⁰ In this study, various techniques were considered for bivariate survival analysis, as a statistical tool for the analysis of time to the blindness of glaucoma patients’ data. Bivariate survival data is a term used to describe the data that measure the time to a given event of interest. In this study, the event of interest was the time of the blindness of the left and right eyes of glaucoma patients.

This study aimed to determine the dependence between the time of the blindness of the left and right eyes using the Copula model popularized by Clayton.²¹This model is an important tool for bivariate survival data and estimates the dependence between variables. In this study, the copula model has been used to couple the marginal survival functions of two cluster observations and form a joint survival function. Parametric Archimedean copula models were also considered to estimate the dependence between the time of the blindness of left and right eyes of glaucoma patients.

The advantage of the parametric method over the semi-parametric method is that, having a baseline distribution for parametric method was good for the sake of simplicity and completeness of further statistical analysis like, quantiles and so on, which account for the popularity of parametric distributions.^22,23

Methods

Study area, design, and period

This study has been conducted at the Ophthalmology Department of Jimma University Medical Center (JUMC), which is one of the oldest hospitals in Ethiopia and the only teaching and referral hospital in the southwestern part of the country. Currently, JUMC with an 800-bed capacity provides many health care services for at least 15 million people in its catchment area. This institution-based retrospective cohort study was performed from January 1, 2016, to August 30, 2020.

Data collection

The data card (related to both left and right eyes of the glaucoma patients) extracted from the ophthalmic patients included such information as socio-demographic and clinical information collected from January 1, 2016, to August 30, 2020, at JUMC, Southwest Ethiopia.The starting point was when the patients started follow-up or were diagnosed at the hospital, and the occurrence of blindness was the endpoint. Admission records were obtained from a total of 315 glaucoma patients.

Ethical clearance

The study protocol was approved by the Research Ethics Review Board of Jimma University, College of Natural Sciences, and the medical director of the Hospital. Due to the retrospective nature of the study, the requirement for obtaining written informed consent from the study participants was waived by the Ethics Review Board, and data were kept anonymous and confidential.

Study population and variables

The study population included all glaucoma patients registered at the Ophthalmology Department of JUMC, Southwest Ethiopia, with regular follow-ups from January 1, 2016, to August 30, 2020. The response variable of the study was the time of the blindness of left and right eyes of glaucoma patients, measured in months between the patient’s admission to the hospital and the occurrence of left and right eye blindness. Right censoring was applied in this study. The exclusion criteria included patients who lost to follow-up during the study period, those who withdrew from the study, and the subjects who did not experience the event before the termination of the study. The events in this study included the occurrence of blindness in the left and right eyes. The independent variables included age, gender (female, male), place of residence (rural, urban), diabetes disease (no, yes), type of medication (Timoglue, Diamox, and Timolol), duration of treatment (short, medium, and long), stage of glaucoma (early, moderate, advanced, and absolute), cup-disc ratio (≤ 0.7, > 0.7), and IOP (normal and not normal).

Inclusion and exclusion criteria

All glaucoma patients ≥ 40 years who were fully informed about study variables were included in the study, whereas the patients with insufficient information were excluded from the study.

Statistical analysis

The endpoints can be strongly correlated in this study, implying an appropriate statistical model that expresses the dependence between the times of two events. However, such standard tools as Cox regression, are not suitable for the simultaneous analysis of the two event times. This study provides an advanced statistical model that incorporates the dependence between the two endpoints in terms of copulas.^24-26 Therefore, in this particular study, the copula model was applied for the bivariate survival analysis. A copula can be used to link two event times by specifying their dependence structure. Copulas have primary and direct applications in the simulation of dependent variables. The copula function is used to model the bivariate survival data.²⁷

Let T₁,T₂ be the two bivariate event times, with the marginal survival function $S_j\left(t_j\right)=P_r\left(T_j>t_j\right),j=1,2$ and the joint survival function $S\left(t_1,t_2\right)=P_r\left(T_1>t_1,T_2>t_2\right)$ , assuming that there are n independent subjects in a study. Where T₁,T₂ are subject to the right censoring for the subject i = 1,2,…,n, we observe $D_i=\left\{\left(Y_{ij},{\sigma }_{ij},Z_{ij}\right);Y_{ij}=\min \left(T_{ij},c_{ij}\right),{\sigma }_{ij}=I\left(T_{ij}<c_{ij}\right),j=1,2\right\}$ , where c_ij is the censoring time, T_ijσ_ij are the censoring indicator, and Z_ij is the vector of covariates. By Sklar’s theorem,²⁸ there exists a unique function c_η that connects two marginal survival functions to the joint survival function so long as the marginal survival function S_j is continuous.

Here, the function c_η is called copula, and its parameter η measures the dependence between the two margins. It allows the dependence to be modeled separately from the marginal distribution. A commonly used family of copula functions is the Archimedean copulas which include Clayton copula, Gumbel-Hougaard copula, Frank copula, and Joe copula.^29,30

Model selection, diagnostics, and goodness of fit test

In this study, the Akaike information criterion (AIC) was used to compare different candidate models, and the model with a relatively small value of AIC has been considered to be a better fit.^31,32 Model diagnostics is one of the fundamental steps to confirm the model assumptions after obtaining a final model using the model selection method. In this study, a diagnostic test was applied for bivariate survival data with families of copula model which is based on Kendall’s τ, Q-Q plot2, and Scatter plot to verify the appropriateness of the distributional assumptions and the adequacy of the model assumed.²⁷ Overall, the goodness-of-fit test for the fitted models was assessed using the Pseudo Maximum Likelihood Estimator.³³

Results

Baseline information and descriptive statistics

In total, 211 (66.9%) out of 315 glaucoma patients were blind, whereas 104 (33.1%) patients were censored. Glaucoma patients were in the age range of 40-84 years with a median survival time of 12 months. The minimum and maximum follow-up times in this study were determined to be one month and 60 months, respectively. During the study period, 107 (33.9%) and 208 (66.1%) patients were female and male, respectively, about 190 (60.3%) and 125 (39.7%) of whom were living in rural and urban areas. In total, 174 (55.2%) and 111 (35.3%) patients had abnormal IOP and were diabetic, respectively. Most of the patients were treated with Timoglue medication, and 88 (27.9%) out of 128 (40.6%) patients were blind. The remaining 97 (30.8%) and 90 (28.6%) patients were treated with Diamox and timolol medication, respectively (Table 1). Moreover, 97 (30.7%) and 98 (31.1%) patients had a medium and long treatment duration, and the other 120 (38.2%) patients had a short treatment duration. A total of 36 (11.4%), 66 (21%), 145 (46.0%), and 68 (21.6%) patients had early, moderate, advanced, and absolute glaucoma, respectively. Similarly, 141 (44.8%) patients had a normal or ≤ 0.7 cup-disk ratios, and the remaining 174 (55.2%) patients did not have a normal or greater than 0.7 cup-disk ratio. In total, 65 (20.6%) female and 146 (46.3%) male patients were blind, respectively. Similarly, 132 (41.9%) and 79 (25.0%) blind patients were from rural and urban areas, respectively. Eventually, 77 (24.4%), 60 (19.0%), and 74 (23.5%) patients were blind and had short, medium, and long treatment duration, respectively (Table 1).

Multivariable analysis

Univariable analysis of all parameters regarding glaucoma patients’ time to blindness was analyzed before choosing variables for the model. The multivariable semi-parametric model included the variables that were significant at a cut of point 0.25 in univariable analysis. Accordingly, at a 0.05 level of significance, gender, place of residence, diabetic disease, IOP, type of medication, duration of treatment, cup-disk ratio, and the stage of glaucoma were significant and affected the time it took for glaucoma sufferers to go blind (Table 2).

Based on the results of the semi-parametric model, the level of significance was determined at 0.05. The effects of each predictor on the response variable were determined considering the Exp (β), which is the hazard ratio and can be interpreted as the predicted change in the hazard for a unit increase in the predictor. After controlling other variables as constant, the estimated hazard ratio of the blindness of male to female patients was HR = 1.005. The result indicates that the blindness risk of male patients was 1.005 times higher than that in female patients (P = 0.01). Similarly, the estimated hazard ratios of blindness for patients who took Timoglue and timolol medication to those who took Diamox medication were HR = 0.887 and HR = 0.976, respectively, indicating that the blindness risk of the patients who took Timoglue and timolol medication was estimated to be 0.887 (P = 0.003) and 0.976 (P = 0.008) times less than that in patients who took Diamox medication, respectively. Eventually, the estimated hazard ratios of blindness for patients who had early, moderate, and advanced glaucoma to patients who had absolute glaucoma were obtained at HR = 0.582, HR = 0.485, and HR = 0.887, respectively. This indicated that the blindness risk in patients who had early, moderate, and advanced glaucoma was 0.582 (P = 0.002), 0.485 (P = 0.001), and 0.887 (P = 0.003) times less than that in patients who had absolute glaucoma, respectively (Table 2).

Parametric model selection and comparison

In this particular study, we considered the Archimedean copula families with a marginal distribution, such as Weibull, Gompertz, and log-logistic. The AIC values of Clayton, Gumbel, Frank, and Joe copula families with Weibull margin were 982.214, 1042.566, 3627.856, and 1090.264, respectively. The Kendall’s tau value of Clayton with Weibull margin was 0.83, this indicates there is high dependence between the time to blindness of left and right eyes of glaucoma patients. Accordingly, the Clayton copula model with Weibull distribution was the minimum among other AIC values of the Archimedean copula models, indicating it to be the best fit for the glaucoma patients’ data sets.

A parametric model for various copula family

The Clayton copula model with Weibull margin of the predictor variables, such as place of residence, age, gender, type of medication, diabetic disease, stage of glaucoma, IOP, duration of treatment, and the cup-disk ratio, was significantly associated with the time of the blindness of the left and right eyes of glaucoma patients at 0.05 level of significance (Table 3). From this analysis, the estimated hazard ratio of blindness for patients’ age was exp (0.01) = 1.01, implying that for a unit increased in age, the hazard ratio of the blindness of patients was significantly increased by 1.01 (P = 0.002), keeping all other variables constant. The estimated hazard ratio of the blindness in patients who lived in urban areas to patients who lived in the rural areas was exp (-0.23) = 0.79, indicating that the blindness risk in patients who lived in urban areas was 0.79 times (P = 0.002) less than that in patients who lived in rural areas (Table 3).

Similarly, the estimated hazard ratio of blindness in male to female patients was exp (0.02) = 1.02. Based on the results, the blindness risk of male patients was 1.02 (P = 0.002) times higher than that in female patients. The estimated hazard ratios of blindness for patients who took Timoglue and timolol medication to patients who took Diamox medication were exp (-0.06) = 0.94 and exp (-0.04) = 0.96, respectively, indicating that the risk of the blindness in patients who took Timoglue and timolol medication was estimated to be 0.94 (P = 0.002) and 0.96 (P = 0.002) times less than that in patients who took Diamox medication, respectively (Table 3).

The results indicated that the estimated hazard ratio of blindness in patients who had diabetic disease relative to patients who had no diabetic disease was exp (0.04) = 1.04, indicating that the blindness risk in patients who had the diabetic disease was 1.04 (P = 0.002) times higher than that in patients who had no diabetic disease. The estimated hazard ratio of blindness for patients who had early, moderate, and advanced glaucoma relative to patients who had absolute glaucoma were exp (-0.33) = 0.72, exp (-0.51) = 0.60, and exp (-0.04) = 0.96, respectively, indicating that blindness risk in patients who had early, moderate, and advanced glaucoma were 0.72 (P = 0.002),0.60 (P = 0.002), and 0.96 (P = 0.002) times less than the blindness risk in patients who had absolute glaucoma, respectively (Table 3).

In line with the above interpretation, the estimated hazard ratio of blindness in patients who had abnormal IOP to patients who had normal IOP was exp (0.03) = 1.03, indicating that the blindness risk in patients who had abnormal IOP were 1.03 (P = 0.002) times higher than the blindness risk in patients who had normal IOP. The estimated hazard ratios of blindness in patients who had a short and medium duration of treatment to patients who had a long treatment duration were exp (-2.48) = 0.084 and exp (-1.26) = 0.28, respectively, indicating that the blindness risk in patients who had a short and medium treatment duration was 0.084 (P = 0.002) and 0.28 (P = 0.002) times less than that in patients who had a long treatment duration, respectively. Eventually, the estimated hazard ratio of blindness in patients who had a cup-disk ratio ≤ 0.7 to patients who had a cup-disk ratio > 0.7 was exp (-0.02) = 0.98, indicating that blindness risk in patients who had a cup-disk ratio ≤ 0.7 was 0.98 (P = 0.002) times less than that in patients who had cup-disk ratio > 0.7 (Table 3). Therefore, a high dependence existed between the time of the blindness of the left and right eyes in glaucoma patients ().

Close

Figure 1. Scatter plot: testing the degree of dependence of Clayton copula with Weibull margin, based on Kendall’s tau value.

Model diagnostics

Once the models are fitted, all the necessary model assumptions should be verified. Standard types of plots are often applied to check and validate the assumptions behind the copula model.

The overall goodness-of-fit test of the copula model can be assessed using the scatter plot. shows the lack of a systematic pattern, indicating that the model fits the data.

Close

Figure 2. Scatter plot of fitted Clayton copula model.

Discussion

A bivariate survival model such as copula model were used to estimate the dependence between the two responses (left and right eyes) of glaucoma patients, by considering the explanatory variables. The copula framework allows the dependency structure between the responses to be isolated from their marginal distributions.²⁵The method consists of introducing copulas as an alternative to the correlation coefficient commonly used as a measure of dependence. An algorithm based on the marginal distributions of random variables is applied to construct the Archimedean copulas. An alternative dependence measure is a copula that overcomes the limitations of correlation as a measure of dependence.^25,26,34 A copula has proved to be useful in a variety of modeling situations and is a relatively new concept that has been applied in survival data analysis.^35,36

Copulas are functions that join or couple multivariate distribution functions to their one-dimensional marginal distribution function. Advantages of using copulas in modeling include (i) allowance to model both linear and non-linear dependence, (ii) arbitrary choice of a marginal distribution, and (iii) capability of modeling extreme endpoints. This study aimed to describe the bivariate survival models, such as the copula model, which can be employed as an alternative for any multivariate data set and estimates the dependence between correlated endpoints. Implementation of Archimedean copula models, based on the copula approach, has been illustrated through the analysis of glaucoma patients’ data.²⁸

In this study, Archimedean copula family models, such as Clayton, Gumbel, Frank, and Joe were explored for modeling time to the blindness of the left and right eyes in the glaucoma patients. The assumption of the model was checked using the measure of dependence Kendall’s tau, normal Q-Q plot2, and scatter plot. All plots of the original data indicated that there was not any deviation from the model assumption and no need for transformation. The analyses of the bivariate survival data, such as parametric bivariate survival analysis, were used to investigate the determinant factors of the time to the blindness of the left and right eyes. The baseline parameter of the Clayton copula model with Weibull margin was statistically significant in the model, indicating the correlation between the bivariate responses. This finding was consistent with those obtained by Sun et al.³⁷Based on the results, the statistical significance of the estimated parameter is a piece of evidence indicating the Clayton copula model to be a better fit than other models.

The prognostic factors considered in this study included age, gender, place of residence, type of medication, diabetes disease, stage of glaucoma, treatment duration, IOP, and cup-disk ratio. All the above predictor variables were found to be the determinant factors for the time to the blindness of the left and right eyes in glaucoma patients, using the univariable analysis. Therefore, these covariates were used in the multivariable analysis to compare the parametric Archimedean copula models. Analysis using the best model (i.e. the Clayton copula model) with the Weibull margin showed that age was an important socio-demographic factor for the time to the blindness in glaucoma patients, implying that the blindness risk increases with age. Similarly, the results of another study conducted by Rossetti et al³⁸ showed that higher age was a significant risk factor for blindness in glaucoma patients. Moreover, time to the blindness of glaucoma patients was significantly associated with the type of medication and duration of treatment. Moreover, the risk of blindness was reduced in the glaucoma patients who took Timoglue and timolol medication and those who had a short and medium treatment duration. Consistently, the results of another study³⁹ showed that the use of various class of glaucoma medication and short and medium treatment duration reduced the hazard of death or blindness.

In this study, the risk of blindness in glaucoma patients was associated with IOP and diabetes. The risk of blindness in glaucoma patients who had abnormal IOP or IOP greater than 21 mmHg was higher compared to patients who had normal IOP, and the patients with diabetes had a high risk of blindness. However, these variables were insignificant in the study conducted by de Voogd et al,⁴⁰ indicating that the results obtained in this study were more reliable than those obtained by de Voogd et al. In addition, the predictor variable of the stage of glaucoma was significantly associated with the time to the blindness of glaucoma patients in this study. The glaucoma patients who had early, moderate, and advanced glaucoma had a lower risk of blindness, compared to those who had absolute glaucoma. This result contradicts the findings of Caprioli & Coleman⁴¹and Drance et al.⁴² Regarding the cup-disk ratio, the risk of blindness in patients who had a cup-disk ratio greater than 0.7 was higher compared to patients who had a cup-disk ratio less than or equal to 0.7. This finding was consistent with those obtained by Gardiner et al³¹ which showed an increased incidence of blindness in patients with a cup-disc ratio larger than 0.7.

The dependence of Kendall’s tau value of Clayton copula model with Weibull margin was estimated to be τ = 0.83, showing an extreme dependence with time to the blindness of the left and right eyes in glaucoma patients. The study results were in line with those obtained by Kendall.⁴³ Eventually, the results obtained in this study were confirmed by Sun et al,³⁷ whose results showed the significance of the copula model for the bivariate data set. Therefore, it can be concluded that the Clayton copula model is the best fit for the data compared to other Archimedean copula family models.

Strengths and limitations of the study

The application of an advanced statistical model that simultaneously incorporates time to two events (i.e. time to the blindness of the left and right eyes in glaucoma patients) is a strength of this study. Regarding the limitations of the present study, one can refer to the fact that this study has been conducted based on secondary data which might be incomplete and biased. Moreover, there was poor data recording on different charts of patients ’information. As for the Archimedean copula family models, some outputs were missing due to the lack of computer computing capacity and the unavailability of enough software packages for the copula model.

Conclusion

Overall, regarding the performance of Archimedean copula family models in terms of model parsimony and goodness of fit test, the Clayton copula model with Weibull margin was a better fit based on its AIC value. Place of residence, age, gender, type of medication, diabetes, stage of glaucoma, duration of treatment, IOP, and cup-disk ratio were significantly associated with and affected by the time to the blindness of the left and right eyes in glaucoma patients. Awareness should be given to the community to reduce the burden of glaucoma.

Acknowledgments

The authors would like to acknowledge Jimma University, the College of Natural Sciences, and the Research and Postgraduate Coordinating Office for their financial support and the permission to carry out this study.

Authors’ contribution

MGG contributed to the study concept, designing, and drafting of the manuscript. RHB, DKD, and KSW KTD participated in data collection and drafted the manuscript, and GMA reviewed the final manuscript.

Availability of data

The data sets analyzed in this study are available from the corresponding author on a reasonable request.

Conflict of interest

The authors declare that they have no competing interests regarding the publication of the present study.

Funding

This study was financially supported by the College of Natural Sciences, Jimma University, Jimma, Ethiopia. The supporting bodies had no role in data collection, analysis, and preparation of the manuscript or in the decision to publish.

References

1. Tham YC, Li X, Wong TY, Quigley HA, Aung T, Cheng CY. Global prevalence of glaucoma and projections of glaucoma burden through 2040: a systematic review and meta-analysis. Ophthalmology 2014; 121(11):2081-90. doi: 10.1016/j.ophtha.2014.05.013 [Crossref]
2. Durowade KA, Salaudeen AG, Akande TM, Musa OI, Olokoba LB, Ibrahim T. Prevalence and risk factors of glaucoma among adults in rural and urban communities of Ilorin West Local Government Area, North-Central Nigeria. Int J Clin Med Res 2016; 3(1):6-12.
3. Wolfs RC, Klaver CC, Ramrattan RS, van Duijn CM, Hofman A, de Jong PT. Genetic risk of primary open-angle glaucoma Population-based familial aggregation study. Arch Ophthalmol 1998; 116(12):1640-5. doi: 10.1001/archopht.116.12.1640 [Crossref]
4. Tielsch JM, Katz J, Sommer A, Quigley HA, Javitt JC. Family history and risk of primary open angle glaucoma The Baltimore Eye Survey. Arch Ophthalmol 1994; 112(1):69-73. doi: 10.1001/archopht.1994.01090130079022 [Crossref]
5. Wolfs RC, Grobbee DE, Hofman A, de Jong PT. Risk of acute angle-closure glaucoma after diagnostic mydriasis in nonselected subjects: the Rotterdam Study. Invest Ophthalmol Vis Sci 1997; 38(12):2683-7.
6. Rudnicka AR, Mt-Isa S, Owen CG, Cook DG, Ashby D. Variations in primary open-angle glaucoma prevalence by age, gender, and race: a Bayesian meta-analysis. Invest Ophthalmol Vis Sci 2006; 47(10):4254-61. doi: 10.1167/iovs.06-0299 [Crossref]
7. Jiang X, Varma R, Wu S, Torres M, Azen SP, Francis BA. Baseline risk factors that predict the development of open-angle glaucoma in a population: the Los Angeles Latino Eye Study. Ophthalmology 2012; 119(11):2245-53. doi: 10.1016/j.ophtha.2012.05.030 [Crossref]
8. Varma R, Ying-Lai M, Francis BA, Nguyen BB, Deneen J, Wilson MR. Prevalence of open-angle glaucoma and ocular hypertension in Latinos: the Los Angeles Latino Eye Study. Ophthalmology 2004; 111(8):1439-48. doi: 10.1016/j.ophtha.2004.01.025 [Crossref]
9. Seah SK, Foster PJ, Chew PT, Jap A, Oen F, Fam HB. Incidence of acute primary angle-closure glaucoma in Singapore An island-wide survey. Arch Ophthalmol 1997; 115(11):1436-40. doi: 10.1001/archopht.1997.01100160606014 [Crossref]
10. Lavanya R, Wong TY, Friedman DS, Aung HT, Alfred T, Gao H. Determinants of angle closure in older Singaporeans. Arch Ophthalmol 2008; 126(5):686-91. doi: 10.1001/archopht.126.5.686 [Crossref]
11. Zhao D, Cho J, Kim MH, Friedman DS, Guallar E. Diabetes, fasting glucose, and the risk of glaucoma: a meta-analysis. Ophthalmology 2015; 122(1):72-8. doi: 10.1016/j.ophtha.2014.07.051 [Crossref]
12. Zhou M, Wang W, Huang W, Zhang X. Diabetes mellitus as a risk factor for open-angle glaucoma: a systematic review and meta-analysis. PLoS One 2014; 9(8):e102972. doi: 10.1371/journal.pone.0102972 [Crossref]
13. Cheng JW, Cheng SW, Ma XY, Cai JP, Li Y, Wei RL. The prevalence of primary glaucoma in mainland China: a systematic review and meta-analysis. J Glaucoma 2013; 22(4):301-6. doi: 10.1097/IJG.0b013e31824083ca [Crossref]
14. Rein DB, Zhang P, Wirth KE, Lee PP, Hoerger TJ, McCall N. The economic burden of major adult visual disorders in the United States. Arch Ophthalmol 2006; 124(12):1754-60. doi: 10.1001/archopht.124.12.1754 [Crossref]
15. Varma R, Lee PP, Goldberg I, Kotak S. An assessment of the health and economic burdens of glaucoma. Am J Ophthalmol 2011; 152(4):515-22. doi: 10.1016/j.ajo.2011.06.004 [Crossref]
16. Bourne R, Steinmetz JD, Flaxman S, Briant PS, Taylor HR, Resnikoff S. Trends in prevalence of blindness and distance and near vision impairment over 30 years: an analysis for the Global Burden of Disease Study. Lancet Glob Health 2021; 9(2):e130-e43. doi: 10.1016/s2214-109x(20)30425-3 [Crossref]
17. Quigley HA, Broman AT. The number of people with glaucoma worldwide in 2010 and 2020. Br J Ophthalmol 2006; 90(3):262-7. doi: 10.1136/bjo.2005.081224 [Crossref]
18. Berhane Y, Worku A, Bejiga A, Adamu L, Alemayehu W, Bedri A. Prevalence and causes of blindness and low vision in Ethiopia. Ethiop J Health Dev 2007; 21(3):204-10. doi: 10.4314/ejhd.v21i3.10050 [Crossref]
19. Giorgis AT. Raising public awareness of glaucoma in Ethiopia. Community Eye Health 2012; 25(79-80):46.
20. Yashin AI, Vaupel JW, Iachine IA. Correlated individual frailty: an advantageous approach to survival analysis of bivariate data. Math Popul Stud 1995; 5(2):145-59. doi: 10.1080/08898489509525394 [Crossref]
21. Clayton DG. A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika 1978; 65(1):141-51. doi: 10.1093/biomet/65.1.141 [Crossref]
22. Ugarte MD. Book Review: Kleinbaum DG, Klein M 2005: Survival Analysis. A Self-Learning Approach. 2nd ed. New York: Springer; 2007.
23. Klein JP. Semiparametric estimation of random effects using the Cox model based on the EM algorithm. Biometrics 1992; 48(3):795-806.
24. Nelson RB. An Introduction to Copulas. 2nd ed. New York: Springer; 2006.
25. Genest C, Rivest LP. Statistical inference procedures for bivariate Archimedean copulas. J Am Stat Assoc 1993; 88(423):1034-43. doi: 10.1080/01621459.1993.10476372 [Crossref]
26. Joe H. Parametric families of multivariate distributions with given margins. J Multivar Anal 1993; 46(2):262-82. doi: 10.1006/jmva.1993.1061 [Crossref]
27. Attaf Z, Ghorbel A, Boujelbène Y. Dependence between non-energy commodity sectors using time-varying extreme value copula methods. Int J Econ Financ Manag 2015; 3(2):64-75. doi: 10.12691/ijefm-3-2-3 [Crossref]
28. Joe H. Dependence Modeling with Copulas. CRC Press; 2014.
29. Nelsen RB. An Introduction to Copulas. Springer Science & Business Media; 2007.
30. Akaike H. A new look at the statistical model identification. IEEE Trans Automat Contr 1974; 19(6):716-23. doi: 10.1109/tac.1974.1100705 [Crossref]
31. Gardiner SK, Johnson CA, Demirel S. Factors predicting the rate of functional progression in early and suspected glaucoma. Invest Ophthalmol Vis Sci 2012; 53(7):3598-604. doi: 10.1167/iovs.11-9065 [Crossref]
32. Okhrin O, Trimborn S, Waltz M. gofCopula: goodness-of-fit tests for copulae. R J 2021; 13(1):467-98.
33. Roberts DJ, Zewotir T. Copula geoadditive modelling of anaemia and malaria in young children in Kenya, Malawi, Tanzania and Uganda. J Health Popul Nutr 2020; 39(1):8. doi: 10.1186/s41043-020-00217-8 [Crossref]
34. Gross AJ, Lam CF. Paired observations from a survival distribution. Biometrics 1981; 37(3):505-11. doi: 10.2307/2530563 [Crossref]
35. Marshall AW, Olkin I. Families of multivariate distributions. J Am Stat Assoc 1988; 83(403):834-41. doi: 10.1080/01621459.1988.10478671 [Crossref]
36. Kumar P, Shoukri MM. Copula based prediction models: an application to an aortic regurgitation study. BMC Med Res Methodol 2007; 7:21. doi: 10.1186/1471-2288-7-21 [Crossref]
37. Sun T, Ding Y. CopulaCenR: copula based regression models for bivariate censored data in R. R J 2020; 12(1):266-82.
38. Rossetti L, Digiuni M, Montesano G, Centofanti M, Fea AM, Iester M. Blindness and glaucoma: a multicenter data review from 7 academic eye clinics. PLoS One 2015; 10(8):e0136632. doi: 10.1371/journal.pone.0136632 [Crossref]
39. French DD, Margo CE. Glaucoma medications and mortality: a retrospective cohort study. Ann Epidemiol 2010; 20(12):917-23. doi: 10.1016/j.annepidem.2010.08.009 [Crossref]
40. de Voogd S, Ikram MK, Wolfs RC, Jansonius NM, Witteman JC, Hofman A. Is diabetes mellitus a risk factor for open-angle glaucoma? The Rotterdam Study. Ophthalmology 2006; 113(10):1827-31. doi: 10.1016/j.ophtha.2006.03.063 [Crossref]
41. Caprioli J, Coleman AL. Intraocular pressure fluctuation a risk factor for visual field progression at low intraocular pressures in the advanced glaucoma intervention study. Ophthalmology 2008; 115(7):1123-9. doi: 10.1016/j.ophtha.2007.10.031 [Crossref]
42. Drance S, Anderson DR, Schulzer M. Risk factors for progression of visual field abnormalities in normal-tension glaucoma. Am J Ophthalmol 2001; 131(6):699-708. doi: 10.1016/s0002-9394(01)00964-3 [Crossref]
43. Kendall MG. A new measure of rank correlation. Biometrika 1938; 30(1-2):81-93. doi: 10.2307/2332226 [Crossref]