**Purpose**:
To compare the ability of linear mixed models with different random effect distributions to estimate rates of visual field loss in glaucoma patients.

**Methods**:
Eyes with five or more reliable standard automated perimetry (SAP) tests were identified from the Duke Glaucoma Registry. Mean deviation (MD) values from each visual field and associated timepoints were collected. These data were modeled using ordinary least square (OLS) regression and linear mixed models using the Gaussian, Student’s *t*, or log-gamma (LG) distributions as the prior distribution for random effects. Model fit was compared using the Watanabe–Akaike information criterion (WAIC). Simulated eyes of varying initial disease severity and rates of progression were created to assess the accuracy of each model in predicting the rate of change and likelihood of declaring progression.

**Results**:
A total of 52,900 visual fields from 6558 eyes of 3981 subjects were included. Mean follow-up period was 8.7 ± 4.0 years, with an average of 8.1 ± 3.7 visual fields per eye. The LG model produced the lowest WAIC, demonstrating optimal model fit. In simulations, the LG model declared progression earlier than OLS (*P* < 0.001) and had the greatest accuracy in predicted slopes (*P* < 0.001). The Gaussian model significantly underestimated rates of progression among fast and catastrophic progressors.

**Conclusions**:
Linear mixed models using the LG distribution outperformed conventional approaches for estimating rates of SAP MD loss in a population with glaucoma.

**Translational Relevance**:
Use of the LG distribution in models estimating rates of change among glaucoma patients may improve their accuracy in rapidly identifying progressors at high risk for vision loss.

^{1}With such low frequency of testing, OLS-derived rates of change would take more than 7 years to detect eyes progressing at a moderate rate of visual field loss.

^{2}

^{3}

^{–}

^{5}which allow data regarding the overall population to influence these estimates; the accuracy of estimates can be increased by “borrowing strength” from the population when fewer data points are available for a particular patient. Mixed model estimates include a fixed-effect component that represents the overall rate of a population and a random-effect component that reflects the degree of deviation of an individual eye from the population average. This process creates eye-specific intercepts and slopes. Although not yet incorporated in routine clinical practice, estimates of rates of change using linear mixed models have been widely applied in research settings.

^{3}

^{–}

^{8}

^{9}

^{–}

^{11}Prior work has demonstrated that the assumption of normally distributed random effects may cause biased estimations of parameters when heterogeneity is present in a population, as would be expected in the rates of progression of glaucoma patients.

^{12}Thus, fast progressors may not be properly identified due to shrinkage to the population mean in a Gaussian model.

*t*and log-gamma (LG), would allow for more accurate estimation of rates of change and detection of eyes exhibiting fast progression.

^{13}Institutional review board approval was obtained for this analysis, and a waiver of informed consent was provided due to the retrospective nature of this work. All methods adhered to the tenets of Declaration of Helsinki for research involving human participants.

^{13}In addition, eyes that underwent trabeculectomy or aqueous shunt surgery were identified using CPT codes. For those eyes, only visual fields obtained before surgery were included, given the likely abrupt postsurgical alteration in the rate of change of SAP MD.

*Y*represents the SAP MD value at time

_{it}*t*of eye

*i*, β

_{0}represents the fixed intercept for the overall population, β

_{1}represents the fixed slope for the overall population, and β

_{0}

*and β*

_{i}_{1}

*represent eye-specific random intercepts and slopes, respectively. In all models, the prior distributions for β*

_{i}_{0}, β

_{1}, and the error term (ε

_{it}) were normally distributed. However, the prior distributions of the random effects differed as noted below. A correlation term with an unstructured correlation matrix was included in the model to account for associations between intercept and slope values. Of note, random effects were placed at the eye level; a more complex model with the eye nested within the patient did not provide additional improvement in the model, thus the simpler model is described here. The correlation between intercept and slope was modeled using an unstructured covariance for the Gaussian and Student’s

*t*-distribution models, whereas a covariance structure previously described was used for the log-gamma model.

^{14}

*t*, and LG distributions were used to model the random effects of intercepts and slopes. The LG distribution is a left-skewed distribution that is sufficiently flexible to allow for more extreme negative values while maintaining a peak close to zero.

^{15}Recent work has suggested that the LG distribution may be a more appropriate distribution in estimating the intercepts and slopes of MD and VFI, given the inherent left skew of these data

^{16}; the majority of eyes have values of intercepts and slopes near zero, but a smaller proportion of eyes have more extreme values.

*t*distributions, the

*brms*package in R was utilized. This package computes estimates of the posterior distributions using Stan, which is a C++ probabilistic Bayesian programming interface using open-source Hamiltonian Monte Carlo (HMC) sampling (Stan Development Team). HMC sampling is thought to be superior to traditional Markov chain Monte Carlo sampling, as this method can achieve a more effective exploration of the posterior probability space without inducing high rates of autocorrelation.

^{17}For the LG distribution, the prior distribution was directly coded into Stan via the

*rstan*R package.

_{0}and β

_{1}), were calculated. Mean posterior estimated intercepts and slopes were calculated for each eye by adding the fixed and random effects of each draw and averaging these values for all draws corresponding to each eye. Eyes were defined as progressors if the one-sided Bayesian

*P*value was less than 0.05. OLS progressors were defined as those with a statistically significant negative rate of change (

*P*< 0.05, one-sided).

*t*, and LG models that were subsequently used to evaluate the simulated eyes. The remaining 20% of the observed dataset was used to create a distribution of residuals for use in the simulations as detailed below.

^{18}Non-progressors were defined as those with a slope of 0 dB/y. Slow, moderate, fast, and catastrophic progressors were defined as eyes with a slope between 0 and −0.5 dB/y, −0.5 and −1.0 dB/y, −1.0 and −2.0 dB/y, and −2.0 and −4.0 dB/y, respectively. These categories have been previously defined

^{10}

^{,}

^{13}and were chosen to simulate eyes with varying rates of disease progression.

^{1}

^{,}

^{19}timepoints with intervals of 6 or 12 months were chosen. At each timepoint, the “true” MD value was based on the simulated intercept and slope. For example, assuming a “true” intercept of −4 dB and a “true” slope of −1 dB/y, “true” MD values would be −4, −4.5, −5.5, −6.5, −7.5, −8.5, −9, and −9.5 dB at the simulated timepoints. As visual field data are affected by noise, we added a residual value to each “true” MD value, according to a previously described methodology.

^{2}

^{,}

^{20}

^{,}

^{21}

^{22}For each test in the sequence of visual fields, a residual value was randomly sampled from the distribution corresponding to the “true” MD. This noise component was then added to the “true” value. For example, for a “true” MD of −4 dB, the distribution of residuals corresponding to −4 dB would be randomly sampled and a residual of +0.5 dB might be selected. This sampling would result in a simulated MD value of −3.5 dB. Using randomly selected residuals for the example above, a simulated set of values might be −3.5, −3, −4.7, −7, −7.8, −9, −8.3, and −11 dB. This simulated eye with these data points mimicking “real-world” observations and their inherent variability was then evaluated by OLS and Bayesian models as described below.

*t*, and LG models (which had been trained on 80% of the dataset) to obtain estimates of the eye-specific intercepts and slopes. Performance of the models was assessed within each simulation setting using bias and by calculating the rates of declared glaucomatous progression. Bias was defined as the difference between the true and estimated posterior slope; negative values of bias indicate underestimation of slope, positive values indicate overestimation of slope, and values closer to zero reflect more optimal prediction. Bias values were pooled across the intercept groups and were compared at each timepoint and for each progressor group using the Kruskal–Wallis test with the Dunn test for pairwise comparisons with Šidák correction. Given multiple comparisons, Bonferroni correction was applied, and an alpha value of 0.05/24 = 0.0021 was used to determine statistical significance. The 95% credible intervals of bias for each model at each setting were also calculated. These intervals reflect that there is a 95% probability that the true value of bias lies within the calculated range. The 95% credible intervals that exclude zero indicate significant underestimation of the slope (e.g., a 95% credible interval of −2.5 to −0.5 indicates that the model significantly underestimates the slope).

*P*value cutoff used for declaring progression in these simulated eyes was set to only allow 2.5% of non-progressor eyes to be erroneously identified as progressors (i.e., a false-positive rate of 2.5%).

*P*value cutoffs were specified uniquely for each intercept range and timepoint. For each setting, the log-rank test was used to determine if the curves were significantly different. Finally, hazard ratios were calculated to determine if the Bayesian models differed in time to declaring progression, using Cox proportional hazards regression. Median time to declared progression was calculated for the different settings as the timepoint at which ≥50% of simulated eyes were declared as progressors.

*t*model included both extreme negative and positive values, whereas the range of slopes of the LG model captured extreme negative values without extreme positive slopes (Fig. 2; “eye-specific slopes” in Table 2). The LG model produced the lowest WAIC value, indicating that the LG model provided the optimal fit for the data compared with the Gaussian and Student's

*t*models (Table 2). When comparing the results of predictive modeling using a limited number of visual fields, Bayesian models consistently performed better compared with OLS, producing lower MSPE values for each predicted visual field MD value (Fig. 2). Overall mean MSPE values of the OLS, Gaussian, Student's

*t*, and LG predictions were 232.6 ± 91.3, 5.2 ± 0.3, 24.2 ± 9.6, and 7.9 ± 0.7, respectively, with significant differences noted between each Bayesian model and OLS (

*P*< 0.01 for pairwise comparisons) at each timepoint until five visual fields were utilized in the models. At this point, the Student's

*t*predictions were no longer significantly different compared with those of OLS (

*P*= 0.84), but MSPE from the Gaussian and LG models remained significantly lower than those of OLS (

*P*= 0.01 and

*P*= 0.02, respectively). When seven visual fields had been utilized, all Bayesian model predictions were non-significant compared with OLS. Of note, differences in the MSPE of the Gaussian, Student's

*t*, and LG predictions were not statistically significant at any timepoint.

*t*models identified a greater number of eyes with faster rates of MD loss among all eyes and progressors. For example, the Gaussian model identified only 8.0% of all progressors as having fast progression and only 0.5% as having catastrophic progression. The LG model identified almost 2 times more progressor eyes as having fast progression (15.2%) and over 5 times more as having catastrophic progression (2.7%) (Table 4).

*t*models in all settings, most notably among fast and catastrophic progressors (Fig. 3). Among fast and catastrophic progressors, mean bias values from the LG, Student's

*t*, and Gaussian models were −0.51 ± 0.49, −0.62 ± 0.51, and −1.20 ± 0.67 dB/y, respectively (

*P*= 0.008, Kruskal–Wallis). When evaluating 95% credible intervals of bias, Gaussian models persistently underestimated the true slope. Gaussian credible intervals excluded zero when using the first three visual fields among moderate progressors; when using the first three, four, or five visual fields among fast progressors; and when using the first three, four, five, six, or seven visual fields among catastrophic progressors. In contrast, all 95% credible intervals of the Student's

*t*and LG models contained zero, indicating that these models did not severely underestimate the slope.

*P*> 0.05, Cox hazard ratio), they were significantly quicker to identify progression compared to OLS among moderate, fast, and catastrophic progressors (

*P*< 0.001, Cox hazard ratio). Median times to progression in the moderate, fast, and catastrophic progressors were consistently lower among Bayesian models compared with OLS (Table 5). The average median time to progression was lower in the LG model (2.8 years) compared with the Student’s

*t*(3.0 years), Gaussian (3.2 years), and OLS (4 years) models.

^{2}previously demonstrated that with the use of OLS 80% of eyes progressing at −2 dB/y would be identified as progressors only after 2.1 years if three visual fields were performed per year (i.e., after six visual fields were completed). Although the benefit of Bayesian linear mixed models over OLS appeared to decline when seven tests were available (Fig. 2), obtaining visual fields at a sufficiently high frequency to procure such a large number of tests is often challenging in clinical practice. The reduction in time to progression using a minimal number of visual fields with Bayesian modeling may be of great value to the clinician. Median time to progression was lower among Bayesian models, especially for the LG model (Table 5).

^{16}previously demonstrated the value of the LG model, as it provided a better fit for SAP data derived from 203 patients in a prospective study compared with a Gaussian model. The authors also constructed a joint longitudinal model using functional SAP and structural optical coherence tomography data, which demonstrated a stronger correlation between functional and structural rates of change when the LG model was utilized. Our work confirms the better fit of the LG model in a much larger dataset.

^{3}

^{,}

^{4}We believe that this discrepancy is due to the greater number of tests that were available in the current dataset (mean of 8.1 visual fields per eye). When evaluating those eyes identified as progressors by OLS but not by the Gaussian, Student’s

*t*, and LG models, the average OLS rates of change were −0.30 ± 0.25 dB/y (interquartile range [IQR], −0.35 to −0.14), −0.31 ± 0.26 dB/y (IQR, −0.38 to −0.14), and −0.28 ± 0.21 dB/y (IQR, −0.34 to −0.14), respectively. These values are reflective of slow rates of change, which would not be as worrisome to the clinician and would be unlikely to lead to severe vision loss. In contrast, when evaluating eyes identified as progressors by the Gaussian, Student’s

*t*, and LG models but not by OLS, the average OLS rates of change were −0.64 ± 0.57 dB/y (IQR, −0.78 to −0.30), −0.62 ± 0.57 dB/y (IQR, −0.76 to −0.26), and −0.68 ± 0.58 dB/y (IQR, −0.86 to −0.31) respectively. OLS was unable to confirm progression among these eyes displaying a faster rate of change, which would be of greater concern and clinical importance. Although the Bayesian models may have identified fewer progressors, the clinical relevance of the progressors identified by these models appears to be greater.

^{3}The magnitudes of these rates are crucial to clinical care; whereas slow progressors may be carefully observed, fast progressors may need to be treated more aggressively in order to prevent vision loss. Therefore, accurate estimation of rates of change is essential to characterizing the nature of a patient's disease. The LG model was able to accurately identify fast progressors while still characterizing the majority of eyes as slow or non-progressors.

^{23}

^{–}

^{25}a linear approximation is likely a reasonable approximation within the limited timeframe used to make most clinical decisions. We also assumed a constant correlation between intercept and slope regardless of disease severity. In clinical practice, severe glaucoma patients are often aggressively monitored and treated, leading to a reduction in correlation between baseline disease (i.e., the intercept) and rate of change (the slope). In addition, the retrospective data collection does not provide insight into augmentation of medical therapy, which could affect rate of progression. It is possible that additional medical or laser therapies may have occurred between visual field tests. The censoring protocol described above only pertained to surgical glaucoma cases. Finally, other potential distributions exist to model random effects which were not evaluated in the present study. We empirically chose the Gaussian, Student’s

*t*, and LG models for comparison given clinical knowledge regarding the distribution of SAP MD in large populations. Further studies should investigate whether other distributions may provide advantages compared with the ones assessed in our work.

*t*distributions. The LG model is sufficiently flexible to accurately characterize non-progressors, slow progressors, and fast progressors. Although Gaussian and LG models are comparable in predicting future SAP MD values, Gaussian models tend to underestimate fast progressors. The LG model was optimal in predicting the rates of change with greatest accuracy while rapidly identifying progressors. These findings may have significant implications for estimation of rates of visual field progression in research and clinical practice.

**S.S. Swaminathan**, Sight Sciences (C), Ivantis (C), Heidelberg Engineering (S), Lumata Health (C);

**S.I. Berchuck**, None;

**A.A. Jammal**, None;

**J.S. Rao**, None;

**F.A. Medeiros**, Alcon Laboratories (C, L, S), Allergan (C, L), Aerie Pharmaceuticals (C), Galimedix (C), Stuart Therapeutics (C), Google (S), Genentech (S), Apple (L), Bausch & Lomb (F), Carl Zeiss Meditec (C, S), Heidelberg Engineering (L), nGoggle (P), Reichert (C, S), National Eye Institute, National Institutes of Health (S)

*HS&DR*2014; 2: 1–102.

*Ophthalmology*. 2017; 124: 786–792. [PubMed]

*Invest Ophthalmol Vis Sci*. 2011; 52: 5794–5803. [PubMed]

*Invest Ophthalmol Vis Sci*. 2012; 53: 2199–2207. [PubMed]

*Ophthalmology*. 2012; 119: 458–467. [PubMed]

*Am J Ophthalmol*. 2012; 153: 1197–1205.e1. [PubMed]

*Ophthalmology*. 2016; 123: 552–557. [PubMed]

*Ophthalmology*. 2021; 128: 48–57. [PubMed]

*Invest Ophthalmol Vis Sci*. 2014; 55: 4135–4143. [PubMed]

*Acta Ophthalmol*. 2013; 91: 406–412. [PubMed]

*Invest Ophthalmol Vis Sci*. 2016; 57: 2012–2020. [PubMed]

*J Am Stat Assoc*. 1996; 91: 217–221.

*Am J Ophthalmol*. 2020; 222: 238–247. [PubMed]

*Bayesian Anal*. 2018; 13: 253–310.

*Biometrics*. 2008; 64: 29–38. [PubMed]

*Biom J*. 2015; 57: 766–776. [PubMed]

*Genet Sel Evol*. 2019; 51: 73. [PubMed]

*Clinical decisions in glaucoma*. St. Louis, MO: Mosby; 1993.

*Br J Ophthalmol*. 2013; 97: 843–847. [PubMed]

*JAMA Ophthalmol*. 2018; 136: 329–335. [PubMed]

*Br J Ophthalmol*, https://doi.org/10.1136/bjophthalmol-2020-318104.

*PLoS One*. 2013; 8: e83595. [PubMed]

*Invest Ophthalmol Vis Sci*. 2006; 47: 5356–5362. [PubMed]

*Invest Ophthalmol Vis Sci*. 2004; 45: 1823–1829. [PubMed]

*Prog Retin Eye Res*. 2007; 26: 688–710. [PubMed]