Nonlinear Trend Analysis of Longitudinal Pointwise Visual Field Sensitivity in Suspected and Early Glaucoma

Abstract

Purpose: : We have shown previously that a nonlinear exponential model fits longitudinal series of mean deviation (MD) better than a linear model. This study extends that work to investigate the mode (linear versus nonlinear) of change for pointwise sensitivities.

Methods: : Data from 475 eyes of 244 clinically managed participants were analyzed. Sensitivity estimates at each test location were fitted using two-level linear and nonlinear mixed effects models. Sensitivity on the last test date was forecast using a model fit from the earlier test dates in the series. The means of the absolute prediction errors were compared to assess accuracy, and the root means square (RMS) of the prediction errors were compared to assess precision.

Results: : Overall, the exponential model provided a significantly better fit (P < 0.05) to the data at the majority of test locations (69%). The exponential model fitted the data significantly better at 85% of locations in the upper hemifield and 58% of locations in the lower hemifield. The rate of visual field (VF) deterioration in the upper hemifield was more rapid (mean, −0.21 dB/y; range, −0.28 to −0.13) than in the lower hemifield (mean, −0.14 dB/y; range, −0.2 to −0.09).

Conclusions: : An exponential model may more accurately track pointwise VF change, at locations damaged by glaucoma. This was more noticeable in the upper hemifield where the VF changed more rapidly. However, linear and exponential models were similar in their ability to forecast future VF status.

Translational Relevance: : The VF progression appears to accelerate in early glaucoma patients.

Introduction

Early detection and accurate monitoring of glaucomatous visual field (VF) loss is important for providing better patient care. Modeling pointwise VF change, as measured by Standard Automated Perimetry (SAP), is clinically important because it can identify localized patterns of progression, which may be masked when using global indices, such as mean deviation (MD). Determining and identifying a pattern of VF change helps clinicians to choose the best strategies when managing patients with glaucoma. However, SAP exhibits high test–retest variability that confounds early detection of deterioration in visual function.^1,2 In clinical practice and research, attempts have been made to measure and predict the rate of VF change and various techniques have been developed to measure such change over time. 

Trend analyses performed on longitudinally collected VF sensitivities remain one of the most frequently used methods for monitoring VF decay.^3–5 Studies done in the past have examined various statistical techniques to identify models that can describe VF decay and predict future VF test results. McNaught et al.² analyzed data using curve fitting software. They showed that complex polynomial models provided the best fit to VF data, but were less accurate when used for forecasting. They recommended less complex linear models for fitting and prediction, and argued against using curvilinear models. Caprioli et al.⁶ explored VF progression in glaucoma using linear, quadratic, and nonlinear exponential models. They concluded that glaucomatous VFs progressed nonlinearly and an exponential decay model provided the best fit and better prediction for VF data. Kummet et al.⁷ systematically evaluated criteria that can be applied to pointwise regression to assist in deciding if clinically useful progression has occurred. Most recently, Bryan et al.⁸ examined fit and predictive ability of various linear and nonlinear models. None of the models they examined were found to be superior for fitting and prediction. However, they recommended using an uncensored linear model as the best compromise between fitting and forecasting. 

One of the major shortcomings of the previous studies examining pointwise VF sensitivities is that these studies ignored the hierarchical structure of ophthalmic data and/or the potential for temporal correlation that can occur when measurements are made repeatedly from the same eyes over time. Also, an approach to fit individual eye data cannot be generalized for a larger population. More complex models that take into account the hierarchical data structure and potential temporal correlation between multiple within eye measurements are needed to find a unified approach for producing the best model fits and predictions. A better approach could be the use of multilevel mixed-effects models. Multilevel mixed-effects models account for grouping within the data (within an eye and between fellow eyes of the same individual) as well as temporal autocorrelation of within group errors. Consequently they generate more valid P values for assessing the significance of change in VF series. 

In this study, we investigated multilevel linear and nonlinear (exponential) mixed effects models to describe VF change over time. Figure 1 presents a schematic diagram showing two approaches examined in this study. A linear progression model (blue line in Fig. 1) assumes a constant rate of VF change over the entire series of tests. An exponential progression model (red line), in which VF sensitivity declines exponentially with time, assumes instead that the rate of change worsens over time, being relatively slow when VF damage is minor. but accelerating as the disease progresses. Note that this differs from the “pointwise exponential regression” approach of Caprioli et al.,⁶ which assumes that the rate slows over time. This is an important issue, not only because it can inform us about the glaucomatous disease process, but also because the two models produce very different predictions of the likely VF status several years in to the future, as can be seen from Figure 1. 

Figure 1. 

 

Schematic figure showing linear versus exponential progression.

In our recent paper,⁹ an analysis of VF MD demonstrated that a nonlinear model fit seemed better for MD. The overall goal of the current study is to validate and extend those findings, assessing whether pointwise VF sensitivity data are better described by a linear or nonlinear (exponential) model by examining sensitivity data at each of the 52 nonblind spot test locations in the 24-2 VF. In addition, we quantify and compare the ability of the linear and exponential models for predicting future VF results. 

Methods

Data

The cohort used in this study has been described previously.⁹ Data from 475 eyes of 244 clinically managed participants with early glaucoma or with high-risk ocular hypertension from the ongoing Portland Progression Project at Devers Eye Institute in Portland, OR, USA were used. The study protocol was approved by the Legacy Health Institutional Review Board. This study complies with the provisions of the Declaration of Helsinki. Consent was obtained from all participants after they were informed about the risks and benefits of participation. 

Initially, participants were tested annually with various functional and structural tests. In 2009, testing was switched to 6-month intervals for all participants. At baseline, participants either had early glaucoma or had ocular hypertension (untreated intraocular pressure repeatedly >22 mm Hg) plus one or more risk factors for developing glaucoma as determined by their eye care provider. Risk factors included age >70,^10,11 systemic hypertension,¹² peripheral vasospasm,¹³ migraine,¹⁴ self-reported family history of glaucoma,¹⁵ and/or previously documented glaucomatous optic neuropathy or suspicious optic nerve head appearance (cup-disc ratio asymmetry >0.2, neuroretinal rim notching or narrowing, disc hemorrhage^16,17), African ancestry,¹¹ and diet-controlled diabetes.¹⁸ Participants having visual acuity worse than 20/40 in either eye or who had worse than mild cataract or media change at baseline were excluded. Other exclusion criteria included any other disease or use of any medications likely to affect the VF, or having undergone intraocular surgery (except for uncomplicated cataract surgery). 

The SAP VF testing was performed using a Humphrey Field Analyzer II (Carl Zeiss Meditec, Inc., Dublin, CA).²⁰ The 24-2 test pattern and the Swedish Interactive Threshold Algorithm (SITA) were used to collect all VF data.²¹ An optimal lens correction was placed before the tested eye and an eye patch was used to occlude the fellow eye. All participants had previous experience with VF testing before entering the study and most had performed multiple previous tests. The VF tests with >33% fixation losses or false negatives, or >15% false positives, were considered unreliable and excluded. Eyes with baseline MD ≤ −6 dB also were excluded from the analyses. Participants with six or more observations per eye meeting the reliability criteria were included in the analyses. Distributions of the sensitivity estimates (left) and MDs (right) at baseline are illustrated in Figure 2. Likewise, the average of baseline sensitivity data and corresponding standard errors are shown in Figure 3. 

Figure 2. 

 

Distribution of baseline sensitivities (left) and MDs (right).

Figure 3. 

 

Average of baseline sensitivities (left) and standard errors between eyes (right) at nonblind spot locations in a 24-2 map.

Models

To find the best fitting and predictive model for pointwise sensitivity data, linear (Sens = μ + β₁.age + β₂.t + ε) and nonlinear (exponential, Sens = μ − e ^(β₁^{+ β}₂^{.age) t} + ε) models with two levels of grouping (participant, eye within participant) were fitted to the sensitivity data. Within group errors (ε) of the fitted models were assessed using the continuous autoregressive (CAR1 [Φ]) method to determine whether serial correlation was present in the data. The CAR1 (Φ) mode serial correlation suggests that the correlation between error terms of longitudinal data decreases exponentially with the duration between them. Follow-up time (year) and the age at time of testing (a risk factor for glaucomatous progression¹¹) were used as covariates. The following four mixed effects models were constructed and compared at each VF location, as detailed in our previous study⁹; M1, linear change over time, uncorrelated residuals; M2, linear change over time, auto-correlated residuals; M3, nonlinear (exponential) change over time, uncorrelated residuals; and M4, nonlinear change over time, auto-correlated residuals. 

The models allow for each eye to have a different disease status at baseline and to have a different rate of VF change over time. Individual differences at baseline and in the rate of VF change over time were incorporated into the models by introducing random effects that characterize the baseline and rate of change in each of the four models. 

Assessment of Model Fits and Forecasts

Before fitting the models, the last test data were removed from the series at each location and the remaining series were fitted. To assess model fits, nested models (M1 vs. M2 and M3 vs. M4) were compared using analysis of variance (ANOVA). Non-nested models (models M1 vs. M3 and models M2 vs. M4) were compared based on Akaike Information Criterion (AIC); AIC = (2θ − 2log[L]) is a commonly used penalized model selection criteria for comparison and selection of models, where θ and L are the total number of parameters to be estimated and the likelihood function for the model respectively. The term −2log(L) the model fit to the data and 2θ acts to penalize models with a greater number of free parameters. So, AIC not only measures the goodness of fit, but also includes a penalty. The penalty discourages a tendency to over fit the data that would occur if we were to simply select a more complex model. Given the set of candidate models, the one with the smallest AIC is the preferred model. At each VF location, Models M1 and M2 were compared to identity the best linear model, and Models M3 and M4 were compared to identify the best exponential model. The best linear and exponential models at each VF location then were compared to find the best fitted model at the given VF location. 

The ability of the models to predict the last observation in each VF series was judged based on the accuracy and precision of the predicted value. To assess predictive ability, sensitivities measured on the last test date were forecast using a model fit from the earlier test dates in the series. Prediction at the last observation within a VF series was made using the fitted linear and exponential models at each VF location and the errors (observed-predicted at last test date) were calculated. The Wilcoxon sum rank test was used to compare the mean absolute prediction errors (MAPE) to assess each model's accuracy. The root mean square errors (RMSE) were compared to assess precision. 

Results

We analyzed 5602 VF tests from 475 eyes of 244 clinically managed participants. The mean age at baseline was 60.59 (± SD = 10.95) years. The mean number of VF tests for each eye was 14 and the mean follow-up duration was 8.36 (± 3.4) years (Table 1). 

Table 1. 

 

Characteristic of the Study Population

 

To choose the best linear and exponential model, an ANOVA comparison was first made among nested models. For the linear models (M1 and M2), model M1 fitted the data significantly better (P < 0.05) than model M2 at two VF locations, and model M2 fitted the data significantly better (P < 0.05) than model M1 at four VF locations. Models M1 and M2 were statistically equivalent (P > 0.05) at 46 (88%) VF locations. Similarly, between the two exponential models (M3 and M4), model M3 performed significantly better (P < 0.05) than model M4 at six VF locations (P < 0.05), and model M4 was significantly better (P < 0.05) than model M3 at three VF locations. In addition, models M3 and M4 were not significantly different (P > 0.05) at 43 (83%) VF locations. 

To assess goodness of fit of the linear versus exponential models, AICs were compared for the best linear model against the best exponential model at each VF location. Figure 4 (left) plots the difference between the AICs of the best linear and the best exponential model at each location displayed as a 24-2 pattern for a right eye. The exponential models had smaller (better) AICs than their linear counterparts at 37 (71%) of the 52 test locations. As seen in Figure 4, the exponential model outperformed the linear model predominantly in the upper hemifield. At 85% of locations in the upper hemifield, the exponential model performed significantly better (P < 0.05) than the linear model. In the lower hemifield, the exponential model provided a significantly better (P < 0.05) fit to the data than the linear model at 15 (58%) of the VF locations. 

Figure 4. 

 

Absolute difference of AIC for the fitted linear (blue, locations where linear model has smaller AICs) and exponential (red, locations where nonlinear has smaller AICs) models at each location in the 24-2 VF. The size of each circle corresponds to the magnitude of the absolute value of the AIC difference, with larger circles indicating bigger differences (left); estimates of the mean rate of VF change in dB/y at each test location (right).

In Figure 4 (right), the mean linear rates of change at each of the test locations are shown. For these data, the observed rate of change was more negative (deterioration) in the upper hemifield. The overall linear rate of change was −0.21 (range, −0.28 to −0.13) dB/y in the upper hemifield, compared to −0.14 dB/y (range, −0.2 to −0.09) in the lower hemifield. It is apparent from Figure 4 that the exponential model performed better predominantly at locations with more rapid deterioration. 

Figures 5 and 6 compare the ability of the linear and exponential models to predict sensitivity at the last test date. The MAPE from the linear model (Fig. 5, left) were statistically equivalent to MAPE from the exponential model (Fig. 5, right) at 51 of 52 test locations. At only one location was the MAPE from the linear model significantly smaller then the MAPE from the exponential model. Table 2 summarizes overall predictive ability of the linear and exponential models across the VF locations. The mean MAPE from the linear model was not significantly different than the mean MAPE from the exponential model (P = 0.26). 

Figure 5. 

 

Mean of the absolute prediction errors at each test location from the fitted linear (left) and exponential models (right).

Figure 6. 

 

The RMS of the prediction errors at each test location from the fitted linear (left) and exponential models (right).

Table 2. 

 

Predictive Ability of the Linear Versus Exponential Models

 

Discussion

Analysis of longitudinally collected perimetry data remains one of the most commonly used approaches for detecting and monitoring change in visual function, especially in patients with glaucoma. Having a reliable assessment of whether visual function is deteriorating is of paramount importance when making clinical management decisions. If the rate of VF change is to be estimated accurately and monitored at individual VF test locations, it is important to know the mode of VF progression over time. Our study demonstrated that the nature of VF change over time depends on the disease stage. The VF locations that display some glaucomatous damage tend to change over time in a manner best fitted by an exponential model. The ability to predict future performance tended to be similar for the exponential and linear models. 

The exponential model discussed in our study is different from the exponential model described by Caprioli et al.⁶ The model we describe has individualized offsets (baseline value) and rate of change. Moreover, when extrapolated back in time, the proposed exponential model used here asymptotes to a flat region that can be thought of as representing the predisease state. 

In our data, we determined that there was an apparent serial correlation in the errors/residuals from the fitted models if traditional linear regression was used. When implementing a generalized least square regression, we observed highly correlated errors at each VF location (Fig. 7). The mean of the autocorrelation parameters (Φ) across the VF was 0.71 (range, 0.57–0.85). Each estimated Φ had narrow 95% confidence intervals, indicating it was significantly different from zero. Interestingly, the higher values of Φ were observed at VF locations where the exponential model fitted the data better, which is not unexpected. If a linear fit is imposed upon data that are accelerating downwards exponentially, the residuals will tend to be positive in the center of the series, and negative at early and late visits in the series. This causes the appearance of significant temporal autocorrelation, which actually is the result of the linear model being a poor description of the mode of change occurring within the series. Given the highly correlated errors that result when linear models are applied to longitudinal series of VF data, studies that have used traditional linear regression may not be optimal, as suggested by Pathak et al.⁹ 

Figure 7. 

 

Autocorrelation parameter estimates at each VF location.

It was clear that the exponential model is preferable to the linear model mainly at more rapidly changing locations. It seems unlikely that the conclusions are driven entirely by eyes changing from “normal and stable” to “glaucomatous and progressing,” since such locations would have a slower rate of change. The lack of a significant difference between models at more stable locations could be due to there being little difference between fits at these locations. The differences in the fit and predictive ability of the linear and exponential models for progressing versus stable eyes can be seen in Figure 8. For a rapidly-changing eye (left panel), the predicted path from the fitted exponential model was in close agreement with the observed data. However, for a stable eye (right panel), where there was little change over time, the difference was minimal. As seen in the left panel in Figure 8, when there is clear deterioration occurring, the prediction of what sensitivity will be in two years' time are several dB worse when using the exponential model (red line). The standard assumption of linear change could lead clinicians to underestimate the likely severity of damage in the near future, with negative consequences for patients. 

Figure 8. 

 

Prediction curves: linear (blue), exponential (red).

Our results suggested that sensitivities actually decline exponentially in early glaucoma, even though VF variability and the high prevalence of stable VF locations in this cohort make linear fits adequate in many cases. The rate of VF change will tend to worsen over time, even in closely followed individuals. This implies that a patient will reach severe vision loss sooner than would be predicted if a constant rate of change (linear model) is assumed. More aggressive treatment may be necessary relatively early in the disease, especially for younger patients, to prevent visual disability. 

The results presented here may in part be a consequence of the logarithmic scaling used when sensitivities are reported in dB. It has been reported that the structure–function relation is linear when sensitivities are expressed on the same scale.^22,23 An exponential decline in dB sensitivity, according to our model, corresponds to a constant rate of change on a linear scale. Therefore, if we assume that change in dB is linear, we are implicitly assuming that the proportion of remaining retinal ganglion cells that die per year is approximately constant throughout the series (e.g., 10% of the remaining cells per year). By contrast, the exponential model used in this study implicitly assumes that the number of retinal ganglion cells that die per year remains approximately constant. Distinguishing between these two possibilities may provide important information about the glaucomatous disease process. 

One of the limitations of this study was that VF locations within an eye were treated as if they were independent. Even though we took account of potential temporal correlation between the multiple measures made within an eye, we did not account for any potential spatial correlations between VF locations. Modelling change in the VF that takes account of correlation across time and space would add an additional layer of complexity. 

In summary, we have shown that VF sensitivity change over time appears to be nonlinear, most noticeably in areas that were progressing more rapidly. These locations were predominantly in the upper hemifield in our dataset. A consequence of this finding is that extrapolating results from linear models of VF change will overestimate the rate of change occurring early in a series, but perhaps more importantly, will underestimate the current and future rates of change, leading to underestimates of the likely future level of damage. 

Acknowledgments

Supported by National Institutes of Health (NIH, Bethesda, Maryland, USA) Grants NIH R01-EY020922 (SKG) and NIH R01-EY19674 (SD), and The Legacy Good Samaritan Foundation, Portland, Oregon, USA. The authors alone are responsible for the content and writing of the paper. 

Disclosure: M. Pathak, None; S. Demirel, None; S.K. Gardiner, None