**Purpose**:
To quantify visual field (VF) variability as a function of threshold sensitivity and location, and to compare weighted pointwise linear regression (PLR) with unweighted PLR and pointwise exponential regression (PER) for data fit and prediction ability.

**Methods**:
Two datasets were used for this retrospective study. The first was used to characterize and estimate VF variability, and included a total of 4,747 eyes of 3,095 glaucoma patients with six or more VFs and 3 years or more of follow-up. After performing PER for each series, standard deviation of residuals was quantified for each decibel of sensitivity as a measure of variability. A separate dataset was used to test and compare unweighted PLR, weighted PLR, and PER for data fit and prediction, and included 261 eyes of 176 primary open-angle glaucoma patients with 10 or more VFs and 6 years or more of follow-up.

**Results**:
The degree of variability changed as a function of threshold sensitivity with a zenith and a nadir at 33 and 11 dB, respectively. Variability decreased with eccentricity and was higher in the central 10° (*P* < 0.001). Differences among the methods for data fit were negligible. PER was the best model to predict future sensitivity values in the mid term and long term.

**Conclusions**:
VF variability increases with the severity of glaucoma damage and decreases with eccentricity. Weighted linear regression neither improves model fit nor prediction. PER exhibited the best prediction ability, which is likely related to the nonlinear nature of long-term glaucomatous perimetric decay.

**Translational Relevance**:
This study suggests that taking into account heteroscedasticity has no advantage in VF modeling.

^{1}The timely identification of clinically significant rates of perimetric progression should prompt consideration of treatment escalation to preserve visual function.

^{2–5}Ordinary least square regression (OLSR) of global indices, VF clusters, or single locations are statistical tools commonly used to assess and quantify glaucoma progression, but they incorrectly assume data homoscedasticity. Heteroscedasticity does not invalidate results obtained with OLSR, but in its presence, this model may not be the best linear unbiased estimator, whereas other linear and nonlinear models may be preferable.

^{6}We have previously shown that pointwise exponential regression (PER) can better predict future changes than pointwise OLSR, and the vulnerability of OLSR to data heteroscedasticity could represent one plausible explanation.

^{7–9}

^{10}

^{11}

*U*test.

^{2}Three different regression models (exponential, linear, logistic) were employed to estimate variability. For each VF series, pointwise regressions were calculated on the threshold sensitivities (dB) over time, excluding the two locations corresponding to the blind spot. For the exponential model, the linear trend (negative or positive) was first determined with simple linear regression, and one of two models (decreasing or increasing) was used. The exponential model was based on a logarithm-transformed linear model and was mathematically defined by the following equations:

*y*is the observed threshold sensitivity (dB),

*x*the time (years), α the intercept, β the slope of the regression line, ε the random error, and Y equal to normal age-matched and location-matched sensitivity + 2 standard deviations (SDs).

^{12,13}

^{2}we fitted an exponential function, which is expressed by the following logarithm-transformed linear model:

*k*

_{1},

*k*

_{2}, and

*k*

_{i}_{–1}are the sensitivity values corresponding to the first, second, and

*i*th − 1 knots, respectively;

_{1}, β

_{2}, and β

*are the slopes to be estimated from the data for the first, second, and*

_{i}*i*th piece of regression, respectively.

*y*

_{[}

_{i}_{]}is the

*i*th known value of observed threshold sensitivity (dB), and

*x*

_{[}

_{i}_{]}is the

*i*th known value of the follow-up (years). The unweighted linear (and exponential) models assume that the variance of

*y*(and ln[

*y*]) are constant for all the observations. The weighted linear regression model has the same specifications of unweighted OLSR, except for the assumption of constant variance. Instead of assuming equal variance, we assume that the

*i*th observation,

*y*

_{[}

_{i}_{]}, has variance

*i*th observation.

^{6}For the exponential model, the linear trend (negative or positive) was first determined with simple linear regression, and one of two models (decreasing or increasing) was used. The exponential model was based on a logarithm-transformed linear model and was mathematically defined by the following equations:

*Y*[

*i*] is equal to normal age-matched and location-matched sensitivity + 2 SD.

^{12,13}

^{12,13}The exponential model asymptotically tends toward the floor or ceiling, respectively, and does not require censoring.

_{1–5}, VF

_{1–6}, VF

_{1–7}, VF

_{1–8}, and VF

_{1–9}). To predict VF number 15, we sequentially performed a regression of first five to the first 13 VFs adding two VFs at every iteration (VF

_{1–5}, VF

_{1–7}, VF

_{1–9}, VF

_{1–11}, and VF

_{1–13}). To predict VF number 20, we sequentially performed a regression of first five to the first 17 VFs adding three VFs at every iteration (VF

_{1–5}, VF

_{1–8}, VF

_{1–11}, VF

_{1–14}, and VF

_{1–17}).

*P*< 0.0001), the intermediate segment (14–32 dB) had a negative slope with a significant variability reduction (

*P*< 0.0001), while the last segment (32–35 dB) had a positive slope with a variability increase (

*P*< 0.0001).

*P*= 0.40), while it slightly increased with the logistic model (

*P*= 0.004).

*P*< 0.001); conversely, no difference in variability was found outside the central 10° (

*P*= 0.61).

*P*= 0.17) or MD rate of change (

*P*= 0.29).

*P*< 0.0001) for VF

_{1–5}, VF

_{1–6}, and VF

_{1–7}, while unweighted PLR had the smallest RMSE (

*P*< 0.0001) for VF

_{1–8}, and VF

_{1–9}. Weighted linear regression had a higher RMSE (

*P*< 0.0001) than unweighted PLR for every comparison.

*P*< 0.0001) for VF

_{1–5}, VF

_{1–7}, VF

_{1–9}, and VF

_{1–11}, while unweighted PLR had the smallest RMSE (

*P*< 0.0001) for VF

_{1–13}. Once again, weighted linear regression had a lower prediction ability (

*P*< 0.0001) than unweighted PLR at every comparison, except for VF

_{1-11}(

*P*= 0.99).

*P*< 0.0001) for VF

_{1–5}, VF

_{1–8}, VF

_{1–11}, and VF

_{1–14}, while unweighted PLR had the smallest RMSE (

*P*< 0.0001 versus PER and

*P*= 0.007 versus weighted PLR) for VF

_{1–17}. Weighted linear regression had a lower prediction ability (

*P*< 0.05) than unweighted PLR for every comparison, except compared with unweighted PLR for VF

_{1–11}(

*P*= 0.35).

^{14–17}Beyond the dichotomy between progression and nonprogression, the quantification of the rate of progression has added importance.

^{18}Even healthy patients experience a physiologic decay due to senescence or cataract, but at slow rates.

^{13,19}Also, glaucoma patients may progress at different speeds, and it is important to distinguish patients with slow rates of decay from those with high rates of decay, because the latter may require more intense follow-up and aggressive therapeutic intervention.

^{18}The identification and quantification of VF progression, however, is hampered by the intrinsic variability of the perimetric examination. Several factors have been associated with VF fluctuation, including patient and technician experience, patient motivation, fatigue, uncorrected refractive error, ethnicity, cognitive level, percentage of false-positive and false-negative test responses, time of the day, and season.

^{10,20–24}In addition, the stage of glaucoma is linked to variability.

^{25}

^{26}Previous studies have investigated the relationship between threshold sensitivity and variability with various methods. Old studies defined the variability as the fluctuation of the variance of the threshold sensitivity values in glaucoma patients clinically judged as stable.

^{27–29}Zulauf and colleagues

^{27}evaluated 29 stable glaucoma patients, and found a mean fluctuation of 4.25 dB

^{2}over the entire VF with a significant association with threshold sensitivity. Werner et al.

^{28}measured variability in 67 glaucoma patients, and confirmed a correlation between variability and threshold sensitivity. Boeglin et al.

^{29}quantified the relationship between variability and sensitivity, and found that variability was lowest at 32 dB, then peaked at 10 dB, and decreased below 10 dB. This methodology, however, has limitations because it is impossible to know if a series is truly stable and changes considered fluctuation may instead be actual VF deterioration.

^{3,30,31}In a large cohort of patients, Chauhan et al.

^{30}found a correlation between the slope of FOS curve, which is a measure of variability, and threshold sensitivity in healthy subjects, glaucoma suspects, and glaucoma patients. Henson et al.

^{3}studied the FOS curves in a heterogeneous cohort of patients, which included 23 healthy subjects and 25 POAG subjects, and confirmed that variability is inversely related to the threshold sensitivity. These findings were further confirmed by Spry and colleagues.

^{31}Nevertheless, FOS studies have several limitations because they are conducted on selected “research” subjects whose variability could be lower than the general glaucoma population, are costly, time-consuming, and suffer from a small sample size.

^{4,32–34}With this method, a group of patients is repeatedly tested in a short period of time with the assumption that glaucoma does not progress over a short time period and all the differences in sensitivity values represent variability.

^{4}Artes and colleagues

^{4}confirmed that variability increases with the severity of glaucoma damage until 10 dB and then decreases. Most of the limitations of FOS studies are also present in test–retest studies. In addition, the latter approach dismisses the possibility of a learning effect, which can be prolonged in some patients.

^{35}

^{2}proposed a method based on linear regression of retrospective large-scale longitudinal data from a clinical practice. This approach has several advantages, including the possibility to study very large cohorts of patients tested in a clinical setting, which are highly representative of the general glaucoma population. In our study, we applied a similar method to that described by Russell et al.,

^{2}and we obtained similar results both in terms of shape and absolute values of the variability curve.

^{3}reported that the variability increased with the reduction in threshold sensitivity as follows: ln(SD) = −0.081 ⋅ Sensitivity (dB) + 3.27, and this formula has been widely used to generate noise in the field of computer-simulated longitudinal VF data.

^{36–38}Gardiner

^{5}estimated the variability on a large cohort of patients, and found the following similar relationship: ln(SD) = −0.070 ⋅ Sensitivity (dB) + 2.70. Neither study, however, quantified the variability for low-sensitivity values, and assumed a similar trend for locations less than 10 and less than 15 dB, in the former and latter study, respectively. We found an analogous relationship because the variability increased with the reduction of threshold sensitivity but with a less steep slope, as expressed by the formula, ln(SD) = −0.027 ⋅ Sensitivity (dB) + 1.79. Previous studies assumed a constant increase of variability below 10 dB, but we demonstrated that it restart decreasing below 10 dB. When we ignored data below 10 dB, variability increased according to the following formula:

^{39}The previous formula by Henson et al.

^{3}has been widely employed to add variability to the simulated sequences, but it is limited by a small sample size, with no measured values below 10 dB. In addition to providing a new descriptive mathematic relationship, we also report the exact amount of variability for each threshold sensitivity value (Supplementary Table S1), and these data can be used to simulate VF sequences more realistically. Because residuals do not follow a normal distribution in the low sensitivity range due to floor effect, the computation of the noise in a Gaussian form, which is based on the SD values, may be inaccurate. In this regard, recent studies

^{40,41}have used the actual distribution of the residuals, empirically calculated on large cohort of patients with the same methodology employed in the present study. These methods, however, are difficult to replicate because raw data has not been shared in detail. In this article (see datasets in Supplementary Files S1–S3), we share with the readers the residuals, fitted, and observed values estimated with the three different models (exponential, linear, and logistic) for every eye, location, and time point. These values may be used to generate a decibel-by-decibel distribution of the residuals to replicate the aforementioned methods.

^{42}examined patients with glaucoma, glaucoma suspects, and healthy subjects, and found that fluctuation was not related to the test location in the glaucoma suspects and healthy subjects, but was slightly increased in glaucoma patients in the upper hemifield, which is the most frequently affected area. Werner et al.

^{28}evaluated 67 glaucoma patients, and found that fluctuation increased with eccentricity, but the effect of test location disappeared with the correction for differences in sensitivity. These results were replicated by Boeglin and colleagues

^{29}in 93 clinically stable eyes of 67 glaucoma patients. In a group of healthy subjects, Heijl et al.

^{12}found contradictory results with increased fluctuation in the peripheral locations. In a subsequent study, Heijl and colleagues

^{33}evaluated a small cohort of glaucoma patients with a test–retest strategy, and found that an increase of variability with eccentricity is seen in patients with mild-to-moderate damage, but the difference was no longer detected for locations with worse glaucomatous damage, defined as total deviation values worse than −10 dB. In a recent study, Gardiner

^{5}tested the impact of eccentricity on variability through linear regression of large-scale longitudinal data and found an increased variability at the peripheral test locations only for sensitivity values above 28 dB. Surprisingly, we observed an increase of fluctuation related to eccentricity, as the variability was higher within 10° from fixation compared with more than 20°. Although these results are similar with those by Gardiner,

^{5}the explanation of these results is not straightforward. Differences in study design, sample size, test strategy, and degrees of VF test can justify discrepancies in the study by Heijl et al.

^{12}According to the hill of vision, central locations have normally higher threshold sensitivity values compared with peripheral locations, and the same sensitivity value may represent a more significant level of damage in the central area because sensitivity at central locations is normally higher than at peripheral locations. The results did not change when we repeated the analysis on the total deviation map, which accounts for deviation from normal age-matched values for each test location.

^{26}Violation of the homoscedasticity assumption did not bias the results of standard OLSR, but this model should not be considered the best linear unbiased estimator.

^{6}

^{6}The first method relies on the application of some function to the dependent variable to reduce or minimize the unequal variance. Log transformation is a popular approach that does not preserve the linear relationship between the variables, but leads to an exponential relationship. Exponential regression works by compressing larger values more than smaller one, and, therefore, it can compensate for heteroscedasticity when the variance increases with the mean. In the specific setting of VF modeling, however, exponential regression is not able to account for heteroscedasticity because the variability decreases as the threshold sensitivity value increases. Other statistical methods that preserve the linear relationship are available (e.g., Box-Cox or signed modulus power transformations), but they have not been applied to VF modeling.

^{6}Robust regression refers to a broad group of models resistant to the presence of multiple outliers, some of which can also deal with heteroscedasticity.

^{6,43}Comparative studies did not find any advantage of robust regressions over classical OLSR in terms of both model fit and prediction, despite the theoretic advantages of the former models.

^{44,45}Weighted linear regression is a third method to deal with heteroscedastic data, and unlike the OLSR, treats each observation differently as it assigns more or less weight to those measures with smaller or higher variance, respectively.

^{46}Weighted linear regression requires the precise and correct estimation of the weights in a large sample, which has been used by a relatively small number of studies.

^{46}Our calculated weights were estimated with a very large glaucomatous population.

^{44,45}may suggest that data heteroscedasticity is not an important primary issue in VF modeling. It may also suggest that ignoring data heteroscedasticity does not significantly affect the goodness-of-fit and the prediction ability of most of the currently used methods to identify and measure perimetric progression, which are based on simple linear regression.

^{7}all methods performed similarly for the data fit, and differences were statistically but not clinically significant. In contrast, prediction ability was considerably different among the methods, and this provides further evidence for the concept that goodness of fit and prediction ability do not coincide.

^{7,47}PER was the best model to predict future sensitivity values in the medium and long term, and models tended to converge with the increase of VF number to make a prediction. We have previously demonstrated that VF loss across the entire perimetric range is best described by a nonlinear (i.e., logistic) rather than a linear function.

^{48}

^{49,50}

**A. Rabiolo**, None;

**E. Morales**, None;

**A.A. Afifi**, None;

**F. Yu**, None;

**K. Nouri-Mahdavi**, None;

**J. Caprioli**, None

*Ophthalmology*. 2014; 121: 733–740.

*Invest Ophthalmol Vis Sci*. 2012; 53: 5985–5990.

*Invest Ophthalmol Vis Sci*. 2000; 41: 417–421.

*Invest Ophthalmol Vis Sci*. 2002; 43: 2654–2659.

*Invest Ophthalmol Vis Sci*. 2018; 59: 3667–3674.

*Heteroskedasticity in Regression: Detection and Correction*. Thousand Oaks: SAGE Publications, Inc.; 2013.

*Invest Ophthalmol Vis Sci*. 2014; 55: 7881–7887.

*Invest Ophthalmol Vis Sci*. 2012; 53: 118.

*Jpn J Ophthalmol*. 2014; 58: 504–514.

*Ophthalmology*. 2017; 124: 1612–1620.

*Arch Ophthalmol*. 1987; 105: 1544–1549.

*Optom Vis Sci*. 2001; 78: 436–441.

*Lancet*. 2015; 385: 1295–1304.

*Ophthalmology*. 1999; 106: 653–662.

*Ophthalmology*. 1999; 106: 2144–2153.

*Control Clin Trials*. 1994; 15: 299–325.

*Am J Ophthalmol*. 2008; 145: 191–192.

*JAMA Ophthalmol*. 2014; 132: 1296–1302.

*Invest Ophthalmol Vis Sci*. 2012; 53: 7010–7017.

*Acta Ophthalmol (Copenh)*. 1994; 72: 189–194.

*Invest Ophthalmol Vis Sci*. 2000; 41: 2006–2013.

*JAMA Ophthalmol*. 2018; 136: 329–335.

*JAMA Ophthalmol*. 2017; 135: 734–739.

*Ophthalmology*. 2003; 110: 1895–1902.

*Surv Ophthalmol*. 2002; 47: 158–173.

*Perimetry Update 1990/91*. Amsterdam: Kugler Publications; 1991: 183–188.

*Perimetry Update 1990/91*. Amsterdam: Kugler Publications; 1991: 175–181.

*Am J Ophthalmol*. 1992; 113: 396–400.

*Invest Ophthalmol Vis Sci*. 1993; 34: 3534–3540.

*Invest Ophthalmol Vis Sci*. 2001; 42: 1404–1410.

*Invest Ophthalmol Vis Sci*. 1999; 40: 648–656.

*Am J Ophthalmol*. 1989; 108: 130–135.

*Am J Ophthalmol*. 1990; 109: 109–111.

*Optom Vis Sci*. 2008; 85: 1043–1048.

*Invest Ophthalmol Vis Sci*. 2002; 43: 1400–1407.

*Br J Ophthalmol*. 2002; 86: 560–564.

*Invest Ophthalmol Vis Sci*. 2011; 52: 3237–3245.

*Invest Ophthalmol Vis Sci*. 2000; 41: 2192–2200.

*PLoS One*. 2013; 8: e83595.

*Transl Vis Sci Technol*. 2018; 7 (3): 22.

*Arch Ophthalmol*. 1984; 102: 876–879.

*Comp Stat Data Analysis*. 2016; 104: 209–222.

*Invest Ophthalmol Vis Sci*. 2013; 54: 6694–6700.

*Invest Ophthalmol Vis Sci*. 2015; 56: 2334–2339.

*Transformation and Weighting in Regression*. London: Chapman & Hall, Ltd.; 1988.

*Graefes Arch Clin Exp Ophthalmol*. 1995; 233: 750–755.

*JAMA Ophthalmol*. 2016; 134 (5): 496–502.

*Ophthalmology*. 1991; 98: 1529–1532.

*Ger J Ophthalmol*. 1995; 4: 175–181.