Determinants of Intraocular Pressure of Glaucoma Patients: A Case Study at Menelik II Referral Hospital, Addis Ababa, Ethiopia

: The main theme of the paper is the well-known problem of glaucoma which is the main cause of blindness worldwide and is also considered a major public health issue. It is usually associated with intraocular pressure above the normal range. The normal range is considered to be 10-21mmHg. Elevated intraocular pressure is a major risk factor for the development and/or progression of glaucoma, and intraocular pressure reduction is a well-known treatment strategy for slowing the progression of the disease. The objective of this article is to identify factors/covariates which affect intraocular pressure on glaucoma patients taking into consideration various demographic, socio-economic, and clinical factors. A retrospective longitudinal cohort study was conducted; the study was based on data from all glaucoma patients who visit at least 3 times repeatedly six waves from January 2016 to December 2018 at Menelik II Referral Hospital Eye Clinic. Profile plots, univariate and multivariate linear mixed effect models were used to explore the major risk factors for the progression of intraocular pressure of a patient. The predictor variables gender (p-value=0.0218), occupation (p-value=0.0025), blood pressure (p-value, 0.0263), diabetes (p-value=0.0139), ocular problem (p-value=0.0290) and type of treatment (p-value=0.0176) found statistically significant effects on intraocular pressure of glaucoma patient. The interaction effects, i.e. time with age (p-value<.0001), time with ocular problem (p-value=0.0002), time with cataract surgery (p-value=0.0002), time with duration of treatment (p-value=0.0014) and time with type of treatment (p-value=0.0262) were found statistically significant on intraocular pressure of glaucoma patient.


INTRODUCTION
Glaucoma refers to a group of eye conditions that lead to damage to the Optic nerve head with progressive loss of retinal ganglion cells and their axons. The two major categories of glaucoma are open-angle glaucoma (OAG) and closed-angle glaucoma. The angleǁ‖ in both cases refers to the drainage angle inside the eye that controls the outflow of the watery fluid (aqueous) that is continually being produced inside the eye [1]. Intraocular pressure is the fluid pressure inside the eye and tonometry is the method eye care professionals use to determine this. The tonometer is calibrated to measure pressure in millimeters of mercury (mmHg) [2]. Glaucoma is usually associated with intraocular pressure (IOP) above the normal range. The normal range is considered to be 10-21mmHg. We can define glaucoma as IOP above 21mmHg in an eye.
IOP is an important aspect in the evaluation of patients at risk from glaucoma. The increased pressure can damage the optic nerve and will loss of vision. Elevated intraocular pressure (IOP) is a major risk factor for the development and/or progression of glaucoma, and IOP reduction is a well-known treatment strategy for slowing the progression of the disease [2].
Glaucoma is one of the most common chronic eye diseases [2]. The prevalence increases with advancing age [3][4][5]. Even though glaucoma can occur at all levels 2. DATA AND METHODOLOGY

Data and Description of the Study Area
This article was based on data from all glaucoma patients who visit at least 3 times from January 2016 to December 2018 at Menelik II Referral Hospital Eye Clinic, Addis Ababa, Ethiopia.

Research Design
The research design is a cohort retrospective study on longitudinal data. The longitudinal data were collected repeatedly between the months of January 2016 to December 2018 from patients 'chart which contains epidemiological and clinical information of glaucoma patients. Each repeated measure was conducted within the one-month interval.

Eligible Criteria
Subjects who dropout (withdraw), intermittent missing values, and loss to follow up within the study period can be included in the study, even if, the repeated measurements are three-time and more but one time or two times patients will be excluded from the study.2.4. Study variables

Response Variable
The intraocular pressure (IOP) of glaucoma patients was used as the response variable which is a continuous variable.

Explanatory Variables
The following covariates were used for the study and coded as shown in Table 1.

Longitudinal Data
Longitudinal data means data sets in which the dependent variable is measured at several points in time for each unit of analysis [12]. Longitudinal data is a special case of repeated measurements. The observations are not independent and are characterized as having between-subject and within-subject variation, time-dependent covariates, and missing data [13].

Linear Mixed-Effects Model
The linear mixed-effects model assumes that the observations follow a linear regression where some of the regression parameters are fixed or the same for all subjects, while other parameters are random or specific to each subject [12,14]. Laird and Ware (1982) [15], describe a two-stage model concept in which the random-effects make up the second stage of the model. Meanwhile, population parameters, individual effects, and within-person variation make up the first stage of the model [15]. This study tried to incorporate a combination of fixed and random effects as predictor variables on the Intraocular pressure. The introduction of random effects affords several non-exclusive benefits. Biological datasets are often highly structured, containing clusters of non-independent observations units that are hierarchical, and the linear mixed-effects model allows us to explicitly model the non-independent in such data. The researchers want to select the linear mixed-effects model.
The general form of the linear mixed-effects model after combining the two stages is approximately normal as follows [12][13][14][15]. A two-stage LMM given by: Stage 1: linear regression model for each subject separately measure the response Y ij for an i th subject at the time ! !" , ! = 1, … , !: ! = 1, … , ! ! response vector Y i for i th subject.
matrix of known covariates and ! ! is q dimensional vector of unknown subject specific regression coefficients and ! ! is a vector of residual component! !" , j=1,2,…,! !
this model describes the observed variability with in subjects.
! !" is the ! ! dimensional identity matrix, Stage 2: Describes between subject variability, that is, explains variability, in the subject specific regression coefficients using known covariates.
! is a ! dimensional vector of unknown regression parameter !"# ! !~!! 0, ! , where D is a general (qxq) covariance matrix with (i,j) element ! !"! ! !" . Now, combining the two levels/ stage models (equation 3.1 and 3.2), we have: Where ! ! is ( ! ! !! ! ) covariance matrix which depends only on i through its dimension i.e. the set of unknown parameters in ! ! will not depend on i.
In Equation 3.3, ! ! represents a vector of continuous response for the i th subject ! ! is an ! ! ! ! design matrix, with known values of the covariates, ! ! , ! ! , … , ! ! , for each of the ! ! observations collected for the i th subject: The ! in Equation (3.3) is a vector of ! unknown regression coefficients or fixed-effect parameters associated with the ! covariates used in the construction of ! ! matrix: 3) is a design matrix that represents the known values of the q covariates, ! ! , ! ! , … , ! ! , for the i th subject.
The columns in the ! ! matrix represents observed values for the ! predictor variables for the i th subject, which have effects on the continuous response variable that vary randomly across subjects. In many cases, predictors with effects that vary randomly across subjects are represented in both the ! ! matrix and the ! ! matrix. In an LMM in which only the intercepts are assumed to vary randomly from subject to subject, the ! ! matrix would simply be a column of 1's.
The ! ! vector for the i th subject in equation (3.3) represents a vector of q random effect associated with the q covariates in the ! ! matrix: Finally, the ! ! a vector in equation (3.3) is a vector of ! ! residuals, with each element in i denoting the residual associated with an observed response at the occasion ! !" for the i th subject.

Assumptions of Linear Mixed Effects Model
The two basic distributional assumptions for the general linear mixed-effects model are The within-group errors are independent and identically normally distributed, with mean zero and variance ! ! I ! ! and they are independent of the random effects.
The random effects are normally distributed, with mean zero and covariance matrix D (not depending on the group), and are independent for different groups.

Methods of Parameter Estimation for LMM
In mixed-model, both maximum likelihood (ML) and restricted maximum likelihood (REML) are the two common methods for parameter estimation in this study.
The likelihood function of the covariance matrix of R and G in the case of ML and REML is given as follows.
and p is the rank of X.

Model Selection
It is an important task to select the best model, i.e., a model that is parsimonious in terms of the number of parameters used, and at the same time is best at predicting (or explaining variation in) the dependent variable [12]. We also use analytic tools, such as the hypothesis tests and the information criteria [12].
There is also the problem of multiple testing that comes with fitting and refitting the model. The issue is made more complicated in the case of longitudinal data where selecting the best model means not only selecting the best mean structure but also the most optimal variance-covariance structure for model selection criteria; i.e. AIC, BIC, and likelihood ratio test is used. The variance-covariance structure with the lower AIC value can be appropriate to select the model. Finally, we can be taken into account the individual Wald test. However, maybe unreliable, especially for small samples. The likelihood ratio test may be one form of assessing goodness of fit for nested models under the normality assumption.

Model Building Process
The process of building an LMM for a given set of longitudinal data was an iterative one that requires a series of model-fitting steps and investigations, and the selection of appropriate mean and covariance structures for the observed data. The model building typically involves a balance of statistical and subject matter considerations; there was no single strategy that applies to every application [12].
In this study, model building starts from a single covariate analysis approach to "screen" out potentially significant variables for consideration in the multi covariate model to identify the importance of each predictor. All variables that are significant at the 25% level are taken into multi covariate model.

Model Adequacy Checking (Diagnostics)
After fitting an LMM, it was important to carry out model diagnostics to check whether distributional assumptions for the residuals are satisfied and whether the fit of the model is sensitive to unusual observations.

Missing Data Handling
There are missing completely at random (MCAR), which refers to amusingness if the missing values are independent of both unobserved and observed data. There is also missing at random (MAR) and occurs if conditional on observed data, the amusingness is independent of the unobserved measurements. The other mechanism of amusingness is referred to as missing not at random (MNAR). In this case, the missing data process is neither MCAR nor MAR but is non-random.
For all methods, we were using the multiple imputation (MI) method. MI has emerged as one of the more common modern options in missing data handling. Essentially, MI uses a variety of advanced techniques, and also accurately predicting missing values is possible because repeated measures are often highly correlated to each other. Multiple imputation (MI) inference involves three distinct phases: the first phase is the missing data are filled in m times to generate m complete data sets, next the m complete data sets are analyzed by using standard procedures, and then finally, the results from the m complete data sets are combined for the inference. These parallel data sets can then be analyzed via standard methods and results combined to produce estimates and confidence intervals that are often more robust than other methods of dealing with missing values.

Statistical Data Analysis
We explore the data using basic descriptive statistics and a profile plot of the mean IOP throughout our study. Univariate and multivariable analyses are done. Graphical analysis is conducted to investigate the shape of the IOP trend and to assess, in an exploratory fashion, the importance of the considered explanatory variables.

RESULTS
The response variable, intraocular pressure of glaucoma patients is continuous. In this article linear mixed model is used to see the progression effect magnitude between the proposed independent variables and the response variable.

Missing Data Handling
The problem of dealing with missing values is a challenge in any analysis of data. There were some missing values in the collected data and treated by applying multiple imputations (MI). Essentially, MI uses a variety of advanced techniques and also accurately predicts missing values because repeated measures are often highly correlated to each other.

Data Descriptions and Descriptive Statistics
Medical cards of glaucoma patients have been reviewed from January 2016 to December 2018 at Menelik II Referral Hospital. The longitudinal response variable was measured for at least 3 visits and at most 6 visits. There were a total of 513 visits from 1042 medical cards. The descriptive statistics showed that 218 (42.1%) of them were females and 295(56.9%) of them were males out of 513 patients. Among this 312(60.2%) were from an urban area and the rest 201(38.8%) were from a rural area. From a total of 513 cases of patients, 205(39.6%) of them are a civil servant, 208(40.2%) of them are merchant, 50(9.2%) of them are farmer, 31(6.0%) of them are a daily worker and the rest 19(3.7%) of them are another work category. And also 129(24.9%) patients had family patient history and the rest 384(74.1) patients hadn't family patient history.

Summary Statistics for Continuous Variables
The patients of average intraocular pressure is 31.61 with a standard deviation of 12.17, and the average age in the year of glaucoma patients included in this study is 55.23 with a standard deviation of 13.96, the minimum age in year 32 and maximum age in year 79. From

Exploratory Data Analysis
The first step to analyze longitudinal data is to explore the pattern of data given over time through graphical display and summary statistics.

The Mean Plot of IOP for each Categorical Covariate
Graphical inspection of the sample means vectors is an important tool for understanding possible relationships among means over time. The following graphs were plots of IOP for each categorical covariate.
The Figure 2 shows that, the mean plot of IOP with the gender of glaucoma patients. As shown in the figure the IOP of males is high. As follow-up time increases, both males and females of glaucoma patients' IOP decrease.
The Figure 3 shows that, the mean plot of IOP with the place of residence of glaucoma patients. The figure exhibits that the initial IOP of glaucoma patients who come from the urban area is high. As the follow-up time increases patients' IOP who come from both areas decreases.
Below Figure 4 shows that, the mean plot of IOP with different types of occupation of glaucoma patients. As shown in the figure at the initial IOP of farmers are high while civil servants are low. As the follow-up time increases, the farmers and daily workers of glaucoma patients' IOP increase but in the case of the merchant of glaucoma patients, IOP decreases.
Below Figure 5 shows that, the mean plot of IOP with smoking status. The figure shows that at the initial smokers are high IOP. As follow-up time increases the average numbers of IOP for both smokers and non-smoker of glaucoma patients decrease.
Below Figure 6 shows that, the mean plot of IOP with alcoholic glaucoma patients. As shown in the figure the initial IOP of alcohol users is high. As follow-up time increases the average numbers of         intraocular pressure for both alcohol users and non-user of glaucoma patients decrease.
The above Figure 7 shows that the mean plot of IOP having blood pressure of glaucoma patients is high. As follow-up time increases the average numbers of IOP for both having blood pressure and not having blood pressure of glaucoma patients decrease.
The above Figure 8 shows that the mean plot of IOP with diabetes history of glaucoma patients. At all-time IOP of glaucoma patients having diabetes is high. As follow-up time increases the average numbers of intraocular pressure for both having diabetes and not having diabetes of glaucoma patients decreases.
The above Figure 9 shows that, the mean plot of IOP with different types of ocular problems of glaucoma patients. As shown in the figure at the initial IOP of ocular infection is high while the corneal scar is low. As the follow-up time increases cataract problem of glaucoma patients IOP increase.
The above Figure 10 shows that the mean plot of IOP of glaucoma patients having cataract surgery. As shown in the figure at the initial IOP of not having cataract surgery of glaucoma patients is high while having cataract surgery of glaucoma patients is low. As the follow-up time increases there exists a high difference between having and not having cataract surgery of glaucoma patients.
The above Figure 11 shows that, the mean plot of glaucoma patients' IOP with a duration of treatment. As shown in the figure at the initial IOP of the short   duration of treatment and medium duration of treatment are equal. As the follow-up time increases long duration of treatment of glaucoma patients IOP decrease.
The above Figure 12 shows that, the mean plot of glaucoma patients IOP with different types of treatment. As shown in the figure at the initial timolol with pilocarpin and timolol with Diamox and pilocarpin are equal while with pilocarpin is low. The follow-up time increases all type of treatments of glaucoma patients IOP decrease.

Model Selection
The primary goal of model selection is to choose the simplest model that provides the best fit to the observed data. There are several choices concerning which fixed and random effects should be included in a linear mixed model (LMM). The issue is made more complicated in the case of longitudinal data where selecting the best model means not only selecting the best mean structure but also the most optimal variance-covariance structure for model selection criteria; like AIC, BIC, and likelihood ratio test will be used. The variance-covariance structure with the lower AIC value can be appropriate to select the model.

Model Fitting for Fixed and Random Effects
In this article, model building starts from a single covariate analysis approach of first doing a single covariate analysis to screen out potentially significant variables for consideration in the multi covariate model to identify the importance of each predictor. All variables that are significant at 25% level, the modest level of significance from one explanatory single covariate regression model are taken into multi covariate model. The purely statistical method is to use an automatic process (stepwise regression), which can be forward the variables are added successively (the most significant at each step) until no variable adds significant information. To select statistically significant covariates for the independent mixed effect models with outcome variables IOP stepwise method was implemented. The models based on only fixed effects were selected with constant random effects at first and then after, the significance of random effects was also checked.

Covariance Structure Selection
For an analysis to be valid, the covariance among repeated measures must be modeled properly. Table 3 shows the value of fit statistics for different types of covariance structures to select the best model for intraocular pressure (IOP).
According to the above table, we can choose the model with the smallest AIC value. From the above correlation structure type Unstructured type had the smallest AIC=8756.6 values and the model had a better fit and greater parsimony than any other correlation structure type.

Univariate Linear Mixed Regression Model Analysis
By using linear mixed effect model regression analysis, we include potentially relevant variables in the model. The place of residence and smoking status of a patient has no statistical significant because their p-value greater than 0.25 (level of significance); but the rest variables has statistically significant variables.

Interaction Effect of Linear Mixed Effect Model
In the multivariable model, we start by including all covariates from univariate analysis. By considering the type three test result of the multivariate linear mixed model we select out possible interaction effects, by comparing the respective p-value with the level of significance alpha 0.25. We first generate all possible interactions of covariates with time then gender with time. The Age in the year of a patient over time, blood pressure of a patient with time, diabetes of a patient with time, the ocular problem of a patient with time, cataract surgery with time, duration of treatment with time, and type of treatment with time have statistical significance effect on intraocular pressure of a patient. The remaining interaction covariates with that of the place of residence, occupation, alcohol intake, and   Table 4 shows the significant factors that affect intraocular pressure of glaucoma patients at α=0.05 level of significance.

Interpretation of the Result
The above table indicates that the fixed effect intercept was 31.76 (SE=6.23) represents an estimate of the mean intraocular pressure at the time (in a month) = 0 excluding all covariates in the model. For a unit change of time the rate of change of the average number of IOP is decreased by 4.87 (p-value <.0001) keeping the other variables constant. And also it indicates that females were 3.31 times lower in their IOP than males through the follow-up time.
The average intraocular pressure of farmers were 6.03 times (P-value= 0.0001) higher than other occupation through the follow-up time keeping the effects of the other variables constant.
The average IOP of a patient who has no blood pressure is 3.34 (P-value=0.0263) times lower than a patient who has blood pressure. And also a patient who has no diabetes is 3.78 (P-value=0.0139) times lower than a patient who has diabetes keeping the other variables constant.
When we compare the average IOP of a patient (3.73) with a patient who has no ocular problem and who has other ocular problems then we find that no ocular problem is 3.73(P-value=0.0151) times lower than other ocular problem keeping the effects of the other variables constant.
The average IOP of a patient having treatment type timolol and pilocarpin is 4.74 (P-value=0.0149) times higher than a patient having treatment type timolol. And also a patient having treatment type timolol, Diamox and pilocarpin are 5.96 (P-value=0.0081) times higher than a patient having treatment type timolol keeping the other variables constant. For a unit change of time, the rate of change of average intraocular pressure of a patient difference between treatment type timolol with diamox and timolol is -1.99 (P-value=0.0026) keeping the effects of the other variables constant.

Variability of Error and Random Effect in the Fitted Model
Variability analyses of both fixed and random effect is also another important aspect. High variability indicates less accuracy or high error in the prediction of the association of outcome with respective risk factors. The subject-specific random intercept variance is found to be 16.2796 (SE =6.4753) (p=<0.0001) with 95% CL (3.5, 29.05). In addition to that the subject-specific random slope variance is estimated to be 0.01286 (SE=0.0350) (p=<0.0001) and 95% CI (0.01015, 0.02188), the estimated variance of the random error is 0.8233 (SE=0.03709) p= (<0.0001) and with 95% CI (0.7386, 0.8830).

Model Diagnosis
As we observed from the below plots of residuals versus predicted values; even if there are some outliers they may not be considered as real outliers because often points identified as outliers are simply a reflection of a skewed distribution. And also to check linearity assumption we plot observed with predicted. Figure 13 suggested that an equal number of predicted values lies above the diagonal line and below the diagonal line it was indicated that the variability of the error in intraocular pressure (IOP) was almost constant i.e. the error does not deviate from each other and this indicates that linearity assumption is satisfied.
Furthermore, according to the probability plots those were shown above. We observed that the normality assumption was supported through the upward nearly straight line of normal plots. In Figure 14 both P-P and Q-Q plot shows that the response variables are normally distributed because the points are scattered nearly on the line as well as the residual against fitted value plot didn't show any systematic pattern thus, it meets the assumption of the distribution of error terms.

Discussion of the Results
The purpose of this study was to determine factors that affect the progression of intraocular pressure of glaucoma patients by using a linear mixed effect model. Intraocular pressure is known as an important factor for glaucoma pathogenesis [16]. Therefore, it is important to identify various factors that may influence IOP. However, there is still debate about these factors, especially age. From the result of the linear mixed model, age in the year of a patient appears to be statistically non-significant for the progression of intraocular pressure of glaucoma patients. Naturally, a patient's IOP increases with increasing age but this study shows that age in the year of a patient has a negative effect on a patient's intraocular pressure. But, the age of a patient over time follow has a positive effect on the IOP of glaucoma patients. Several studies have shown that IOP increased with age in study populations, whereas others suggested no relationship with age. In a Barbados Eye Study, among the black participants, the mean IOP increased ~ 1mmHg for every increase of 10 years of age [17]. However, a Japanese study showed that IOP correlated inversely with age in men and women [18].
According to our study, gender had been found a significant effect on the progression of IOP of glaucoma patients. Our findings show that men's had greater IOP compared with women's. A pattern was also reported in the Gutenberg Health studies [19]. In contrast, the Barbados Eye study reported men had lower IOP [20] while the Framingham Eye study reported no association between sex and IOP [21].
In addition to age and gender, many systemic factors showed an effect with IOP. This study also showed that blood pressure and diabetes have a significant positive correlation with IOP. It has been postulated that increased blood pressure leads to elevated ciliary artery pressure, increasing the ultrafiltration of the aqueous humor and thus increasing IOP [22].
Occupation of a patient is found to be another significant variable for the progression of intraocular pressure of glaucoma patients. This study found that farmers have greater intraocular pressure than patients who have other work.  According to the study patients' duration of treatment is not an important variable for the progression of intraocular pressure of glaucoma patients. But the study also shows that patient's duration of treatment over time is a significant variable for the progression of intraocular pressure. If the patient's duration is too long and gets sufficient care with necessary treatments, intraocular pressure becomes decreases as a number of treatments follow-up time increases.
Regarding treatment type, the finding indicated that the patient with treatment type has a significant effect on the progression of intraocular pressure of glaucoma patients. The results of this study show that timolol is more effective for lowering IOP than another treatment type. The results also show that a patient's treatment type over time has a significant effect on the progression of intraocular pressure of glaucoma patients.
This study found that th e ocular problem of a patient is found to be another significant variable for the progression of intraocular pressure of glaucoma patients. A patient with the ocular problem has greater IOP than Patients with no ocular problem. The result also shows that a patient's ocular problem over time has a significant effect on the progression of IOP.
This study also found that the cataract surgery of a patient over time has a significant effect in lowering intraocular pressure. The lowering effect is also well established by Surgeons, he aware that cataract surgery in glaucoma patients significantly reduces intraocular pressure [23].