The Relationship between Variability and Sensitivity in Large-Scale Longitudinal Visual Field Data

Abstract

Purpose.: Evaluation of progressive visual field (VF) damage is often based on pointwise sensitivity data from standard automated perimetry; however, frequency-of seeing and test-retest studies demonstrate that these measurements can be highly variable, especially in areas of damage. The aim of this study was to characterize VF variability by the level of sensitivity using a statistical method to quantify heteroscedasticity.

Methods.: A total of 14,887 Humphrey 24-2 SITA Standard VFs from 2736 patients (2736 eyes) attending Moorfields Eye Hospital from 1997 to 2009 were studied retrospectively. The VF series of each eye was analyzed using pointwise linear regression of sensitivity over time, with residuals (difference from fitted-value) from each regression pooled according to both observed and fitted sensitivities.

Results.: The median (interquartile range) patient age, follow-up, and series length was 64 (54–71) years, 5.5 (3.9–7.0) years, and 6 (5–7) VFs, respectively. The inferred variability as a function of fitted-sensitivity was in good agreement with previous estimates. Variability was also described as a function of measured sensitivity, which confirmed that variability increased rapidly as the observed sensitivity decreased.

Conclusions.: This study highlights a new approach for characterizing VF variability by the level of sensitivity. A considerable strength of the method is that inference is based on thousands of clinic patients rather than the tens of subjects in test-retest studies. The results can help distinguish real VF progression from measurement variability and will be used in models for glaucoma progression detection.

Introduction

Standard automated perimetry (SAP) remains the gold standard for assessing and monitoring visual field (VF) impairment in glaucoma. However, a significant weakness of SAP is that threshold measurements are highly variable, as shown by frequency-of-seeing (FOS) ^1–3 and test-retest studies. ^4–11 In glaucoma management, VF variability represents an enormous hindrance for clinicians interpreting VF test results and determining progression. 

Previous studies have explored the relationship between VF variability and measured-sensitivity. Notably, Henson et al. ² collected FOS data, using SAP, from four VF locations in 71 subjects with a range of VF damage. The authors concluded that response variability (SD) increased with decreasing sensitivity, and summarized the relationship using the function: $lo{g}_e\left(SD\right)=A·sensitivity\left(dB\right)+B,$ where A = −0.081 and $B=3.27$; $R^2=\mathrm{0.57.}$ According to this function, as sensitivity decreases there is an exponential increase in response variability. However, the inferences from this study are somewhat limited in their relevance to clinical SAP test strategies, where measured-thresholds are not based on FOS curves. Furthermore, there was a paucity of measurements with low sensitivity, and no data below 10 dB. Artes et al. ⁹ overcame these limitations using test-retest data, and studied the variability properties of SAP threshold estimates from the Full Threshold and Swedish Interactive Threshold Algorithms (SITA) strategies of the Humphrey Field Analyzer (Carl Zeiss Meditec, Dublin, CA). The authors examined one eye each of 49 glaucoma patients, with a large range of VF damage, four times with each test strategy. Variability was shown to increase as the “best available estimate of sensitivity” (mean sensitivity from test-retesting) decreased up until approximately 10 dB, after which variability reduced as sensitivity declined further. Remarkably, for glaucomatous locations with low differential light sensitivity, such as 10-dB sensitivity, the authors showed that variability was so large that the 5th and 95th percentiles of the retest values spanned 0 dB to 28 dB, respectively. 

The aim of this study was to characterize variability by level of sensitivity in SITA-Standard VFs using a statistical method to quantify heteroscedasticity in longitudinal data. The main strength of the current approach is that inference is based on thousands of patients tested in standard clinic conditions, rather than the tens of subjects reported in formal test-retest studies. 

Methods

Study Sample

This was a retrospective analysis of 14,887 VFs from patients attending the Glaucoma Service of Moorfields Eye Hospital, London between 1997 and 2009. All VFs were made anonymous before analysis. In total, 2736 eyes from 2736 patients were included (if more than one eye was eligible, one eye was randomly selected). The study adhered to the tenets of the Declaration of Helsinki and was approved by research governance committees of the participating institutions. All data were transferred to a secure computer database at City University London. 

VF Testing

All VFs were carried out with the Humphrey Field Analyzer (Carl Zeiss Meditec) using the 24-2 test pattern with a Goldmann size III target and the SITA Standard testing algorithm. Unreliable VFs were discarded according to the following criteria: 20% or more fixation losses or 15% or more false-positive errors. VFs were not discarded according to false-negative errors because an increased false-negative rate is strongly associated with field status. ¹² Each VF series examined had to be at least five long for inclusion into the present study, and the first VF was then discarded to reduce perimetric learning effects. ^13,14  

Analyses

The relationship between variability and measured-sensitivity was analyzed by examining the results from linear regression of pointwise sensitivity over time. Linear regression estimates the relationship between a predictor variable $X$ and a response variable $Y$, expressed as the following equation: 

 

where $a$ is the intercept, $b$ is the slope (gradient), and $E$ is the error term. ¹⁵ The error term ( $E$) represents the part of the response variable ( $Y$) that is not explained by the predictor variable ( $X$). For the fitted function, the error term is estimated from the residuals, which are the vertical deviations from the fitted line. The simplest and most common method for fitting a regression line is to minimize the sum of squares of these residuals; this approach is known as ordinary least squares linear regression (OLSLR). The method assumes that the predictor variable is error-free and that the error-term is normally distributed with mean zero and constant variance across the range of the measurement. Constant variance is known as homoscedasticity (“equal scatter”); however, if the residuals are dependent on the on the magnitude of the measurement, this is referred to as heteroscedasticity (“unequal scatter”), which indicates that the variance of the error term is not uniform across observations. Importantly, heteroscedasticity does not cause OLSLR coefficient estimates (i.e., the slope and intercept terms) to be biased; however, estimates of the variance of the coefficients are biased. ¹⁶ Thus, OLSLR in the presence of heteroscedasticity provides an unbiased estimate for the relationship between the predictor variable and the response variable, but standard errors and consequently inferences on statistical significance are affected. The residuals from linear regression are informative when investigating heteroscedasticity. A scatterplot of the squared or absolute residuals against $X,Y$ or Ŷ (the fitted-value of $Y$) may be used to assess for nonconstant variance of the error term. 

To investigate heteroscedasticity in VF measurements, residuals from pointwise linear regression of sensitivity (the $Y$ variable in dB) against time (the $X$ variable in years) were examined for each eye's series of VFs (2736 eyes) using OLSLR and Tobit linear regression (TLR). ¹⁷ We have previously discussed the usefulness of TLR for the analysis of VF sensitivity measurements. ¹⁸ If the outcome variable in a linear model is censored (as is the case for threshold sensitivity), the assumptions of the OLSLR model are not valid; however, TLR provides a valid alternative. TLR uses a latent dependent variable, which respects left- and/or right-censoring and predicts the response only within the range specified. Thus, in contrast to OLSLR, TLR reconciles that a VF threshold of 0 dB may be 0 dB, or, a value less than 0 dB (were the perimeter to have a greater dynamic range). 

Regression models (TLR and OLSLR) were fitted to each eye's series of VF sensitivity values (for each VF test point separately) against time (the patient's age at each VF test); for example, see Figure 1. The two locations adjacent to the blind spot were excluded, giving 52 thresholds for each VF. Next, the least squares residuals (Ê) were extracted: 

 

and grouped (binned) according to $Y$ (measured-sensitivity) in the range [0 to 36] dB. Residuals were also binned according to the fitted-sensitivity value (rounded to the nearest whole decibel), Ŷ, as this value estimates “true” sensitivity. These two binning methods ask subtly different questions, which are now outlined. Binning by measured-sensitivity establishes variability conditional on the measured-threshold, and asks the question, given a measured-threshold what is the underlying range of values for the “true” value? This approach is akin to the method used by Wall et al., ⁴ where retest thresholds were compared with baseline threshold; in this case, as the authors state, “the first test is not meant to be the true sensitivity, as it has its own variability.” On the other hand, binning by fitted-sensitivity investigates variability according to the estimated true value, and asks the question, given a “true” threshold value what is the range of measured-thresholds expected for any given test? This approach mirrors the one used in Artes et al., ⁹ where thresholds were compared with the mean of several retest thresholds. The difference in the two binning strategies is shown in Figure 1. In this figure, residuals (indicated by the dashed lines) are associated with an OLSLR fitted-sensitivity bin of 31 dB; however, if the residuals are stratified by measured-threshold, they are pooled into the following bins: 30, 31, 31, 30, 33, 31, 33, and 28 dB. All statistical analyses were carried out in the open-source programming language, R. ¹⁹  

Figure 1. 

 

The solid line represents OLSLR of a patient's VF measurements (at a given point in their VF) against their age. The dashed lines represent the differences between the measurements and their fitted-values (“residuals”). In this example, the TLR line is no different from the OLSLR line.

Results

The study included 2736 eyes of 2736 patients. The Table summarizes the characteristics of the study sample: 

Table. 

 

Characteristics of the Study Sample

Measurement	Median (Interquartile Range)
Baseline age	63.7 (54.0–71.2) y
Baseline pointwise sensitivity	28 (24–30) dB
Follow-up period	5.5 (3.9–7.0) y
Number of VF tests	6 (5–7)

In total, 142,272 (52 × 2736) regression lines were fitted and almost one million residuals extracted for each linear regression method. The residuals associated with the fitted- and measured-sensitivity levels were evaluated as histograms for both the TLR and OLSLR models (see Fig. 2). Examination of the spread of residuals at different sensitivity levels highlights significant differences in the distributions of residuals; for both linear models, distributions were relatively compact at high VF sensitivities (26–36 dB) but stretched substantially as sensitivity decreased to a level of 10 dB. Sensitivity values near 10 dB were associated with residuals spanning almost the entire dynamic range of the instrument. Below 10 dB, there was some disparity between the distributions of residuals for the TLR and OLSLR models (see Fig. 2E). The results suggest that any reduction in variability below 10 dB is explained by the limited dynamic range of SAP, as evidenced by the negative skew in the distributions of residuals in Figure 2. At low sensitivity, in the range 0 to 21 dB, the distributions of residuals for OLSLR are negatively skewed due to the limited dynamic range of SAP, and this is not seen to the same extent in the corresponding TLR distributions. The SDs of the residuals for each fitted- and measured-sensitivity level are illustrated in Figures 3A and 3B, respectively; the blue lines are derived by binning residuals from OLSLR whereas the green lines are derived by binning residuals from TLR. These are contrasted with variability estimates from Artes et al. ⁹ (red points) and Henson et al. ² (yellow points). 

Figure 2. 

 

Distributions of residuals grouped according to fitted-values values (top row) and measured-values (bottom row) for (A, E) 0-, (B, F) 10-, (C, G) 20-, and (D, H) 30-dB levels. Distributions are illustrated as “back-to-back histograms”: the upper histogram (light gray bars) indicates residuals derived from OLSLR, whereas the lower histogram (dark gray bars) indicates residuals derived from TLR.

Figure 3. 

 

( A) Plot of the SD of residuals for each fitted-sensitivity level. Given a “true” threshold value, this plot illustrates the range of measured-thresholds expected for any given test. The solid blue line is derived by binning residuals associated with fitted-sensitivity from OLSLR, whereas the solid green line is derived by binning residuals associated with fitted-sensitivity from TLR. (B) Plot of the SD of residuals for each measured-sensitivity level. Given a measured threshold, this plot illustrates the underlying range of values for the “true” value. The dashed blue line is derived by binning residuals associated with measured-sensitivity from OLSLR, whereas the dashed green line is derived by binning residuals associated with measured-sensitivity from TLR. In both Figures 3A and 3B, variability estimates are contrasted with those from Artes et al. ⁹ (red points) and Henson et al. ² (yellow points). The blue background is a smoothed color density representation of the scatterplot for all absolute residuals.

Discussion

This study is original in its investigation of the relationship between variability and sensitivity in visual fields. Previous studies have investigated this relationship with small numbers of patients, or people with normal vision, using FOS experiments or test-retest data in a research environment. In contrast, in our study, variability was analyzed using linear regression analysis of retrospective large-scale longitudinal data. The results confirm that a reduction in VF sensitivity is accompanied by a large increase in variability. ^{1,4,6–9,20,21} For example, 25 years ago, Heijl et al. ²¹ analyzed test-retest VF data and concluded that, in areas of moderate to advanced glaucomatous damage (8 to 18 dB loss), the 95% prediction interval associated with the VF sensitivity measurement spanned the entire dynamic range of the instrument. Very similar results have been shown for the SITA standard testing algorithm. ^4,9 The association between a decline in VF sensitivity and an increase in response variability could be caused by a loss of ganglion cells (due to glaucomatous damage), or relocation of the stimulus to the peripheral visual field where there are fewer ganglion cells. ²² Previous studies have suggested that a reduction in the number of stimulated ganglion cells, by a reduction in stimulus size ^23,24 or pathological damage, ^2,9 may lead to an increase in response variability. Therefore, these results highlight the considerable difficulties in reliably identifying VF damage in areas of moderate to advanced glaucomatous VF loss. 

There are notable differences between our study and previous investigations of VF variability. Here, linear regression was used to estimate VF variability whereas test-retest studies assume that VF measurements remain unchanged over the testing period and take an average of threshold values to approximate true sensitivity. Previous research has suggested that the standard linear model is appropriate for tracking visual field progression. ^25,26 Moreover, standard pointwise linear regression remains the most popular method to investigate rates of VF progression. ^27–30 A limitation of our study is that the method estimates two parameters: a slope and an intercept term, whereas test-retest studies estimate the latter parameter only, as they assume that the measurement does not change in the period of time over which the data are captured. Accordingly, a limitation of test-retest data is that it assumes the absence of a perimetric learning effect, which can persist over tens of VF tests in some patients. ^13,14,31 Nevertheless, our findings agree well with those in Artes et al., ⁹ but are possibly more robust because of the large amount of data analyzed. Furthermore, our sample consists of patients from general ophthalmic practice and not “perimetry athletes” who are very familiar with SAP, and consequently may record less noisy measurements. Because our analysis was carried out on retrospective data, the method could be easily applied to other perimetry and imaging devices used in glaucoma assessment, saving time and money compared with prospective test-retest or FOS studies. 

Both our approach and the test-retest method assume that the fitted-value or mean of the measurements is the “true” value. Using the mean to approximate true sensitivity makes statistical assumptions about VF measurements; in particular, the arithmetic mean is heavily influenced by outliers, and is not robust to skewed and/or truncated distributions. This is supported by examining estimates of VF variability in the test-retest studies by Artes et al. ⁹ and Wall et al. ⁴ : when VF thresholds are equal to 0 dB, the arithmetic mean is not robust to the truncated nature of VF measurements, and the derived variability (expressed as the SD of retest thresholds) is consequently reduced (see Fig. 3A). This truncation effect is very evident in the distributions shown in Figure 2 and in similar figures in other studies. ^4,9 Our results emphasize that a reduction in VF variability for thresholds below 10 dB can largely be attributed to this truncation. 

TLR attempts to take into account the limited dynamic range of SAP (the truncation effect). Unlike OLSLR, TLR resolves the fact that a VF threshold of 0 dB could be 0 dB, or a value less than 0 dB. Thus, fitted-values from TLR are less affected by the floor effect than fitted-values from OLSLR. This is evident in Figure 2E: residuals binned by a measured-threshold of 0 dB are symmetrically distributed for the TLR model but are negatively skewed for the corresponding distribution derived from the OLSLR model. As expected, positive skew is apparent in both distributions in Figure 2A, where residuals are binned by a fitted-value of 0 dB, because no measured-value can fall below 0 dB. However, Figure 3A demonstrates that the SD of the residuals (grouped according to the fitted-value) is larger for TLR than OLSLR. This is explained by the fewer number of negative residuals associated with TLR as it respects the left-censoring of VF measurements, whereas OLSLR does not respect left-censoring. 

In contrast to previous research, our study also investigated variability as a function of measured-sensitivity. Stratification by measured-sensitivity indicates the variability associated with a measured-threshold, rather than the variability associated with “true” VF sensitivity. The former is arguably more important for the clinical scenario when interpreting VF results from a single test because the range of possible “true” values associated with the observed measurement is more informative than knowledge of the variability associated with the underlying true VF sensitivity. The dashed dark green and dark blue lines in Figure 3B illustrate the variability associated with measured-sensitivity for the TLR and OLSLR models, respectively, and Figures 2E to 2H illustrate the distribution of residuals associated with measured-sensitivity. Both sets of figures suggest that variability increases rapidly as measured-threshold decreases. This finding has implications for clinical SAP testing. Wall et al. ⁴ showed that threshold estimates below 20 dB have little value for predicting the value at retest; however, their results were based on only one baseline measurement and a single retest value. Our results are similar to the findings of Wall et al. ⁴ and support their suggestion that for the purposes of detecting VF progression, examination of damaged test locations could be worthless. 

Significantly, our study questions the validity of pointwise OLSLR to detect VF progression. As outlined earlier, linear regression estimates the relationship between a predictor variable and a response variable. Linear regression assumes that the variance of the error is constant (homoscedastic), and thus the regression is measured with equal precision throughout the range of the response variable. Put differently, if the error is not constant (heteroscedastic), then the regression is estimated less accurately where the noise is large compared with where the noise is small. Depending on the nature of the heteroscedasticity, significance tests can be inflated or reduced, leading to incorrect conclusions regarding, for example, the significance of a linear regression slope. Because OLSLR is designed to minimize the squared errors between the measurements and the regression line, if the error is heteroscedastic it will put more weight on measurements with the largest error terms and the least signal. Clearly, the error term associated with VF measurements is heteroscedastic and not normally distributed. Statistical solutions to account for errors that are not independent and identically distributed are available. For example, use of so-called robust SEs in linear regression could be used. ^16,32 Furthermore, recent research from Caprioli et al. ³³ suggest that pointwise exponential regression (PER) provides a robust estimate of rates of VF decay that may predict future global indices more accurately than standard linear regression. An interesting extension to our research would be to use the same methodology to investigate the relationship between variability and sensitivity level with the PER model; the results from such an analysis would indicate if the PER model overcomes the heteroscedasticity and non-normal distribution of residuals that are seen with OLSLR and TLR. 

One application of our results is in VF simulations. At present, there is no gold standard for the assessment of glaucoma progression, so our results could be used to generate a nonparametric model of VF variability. Visual fields simulated with known progression characteristics could be generated and different glaucoma progression detection algorithms could be compared. Until now, most VF simulation models ^30,34–36 have been based on Henson et al.'s ² equation for VF variability, which may be less appropriate than the results shown here, as results from Henson et al. ² were based on FOS curves with no thresholds below 10 dB. Furthermore, our results are based on thousands of patients in general care and not research patients with a superior knowledge of perimetry testing. Simulating VFs using our results would also allow progression criteria to be established with known sensitivity and specificity characteristics, which would be important for application in clinical trials. Finally, as stated earlier, because our study was carried out on retrospective data, the analysis could be easily used in other visual field and imaging devices, saving time and money compared with prospective test-retest or FOS studies. 

Footnotes

 Preliminary work was presented at the annual meeting of the Association for Research in Vision in Ophthalmology, Fort Lauderdale, Florida, May 2011.

Footnotes

 Supported in part by the UK National Institute for Health Research Health Services Research program (project 10/2000/68). A proportion of funding was from the Department of Health's National Institute for Health Research Biomedical Research Centre at Moorfields Eye Hospital NHS Foundation Trust and University College London Institute of Ophthalmology (RAR, RM, DFG-H). The chair at University College London (DFG-H) is supported by funding from the International Glaucoma Association. The views expressed are those of the authors and not necessarily those of the NHS, the National Institute for Health Research, or the Department of Health.

Footnotes

 Disclosure: R.A. Russell, None; D.P. Crabb, None; R. Malik, None; D.F. Garway-Heath, None