How Useful Is Population Data for Informing Visual Field Progression Rate Estimation?

Abstract

Purpose.: Bayesian estimators allow the frequency of visual field progression rates in the population (the prior distribution) to constrain rate estimates for individuals. We examined the benefits of a prior distribution accounting for one of progression's major risk factors—whether intraocular pressure is treated—to gauge the maximum benefit expected from developing priors for other glaucoma risk factors.

Methods.: Our prior distribution was derived from published data from either treated (matched–prior condition) or untreated (unmatched–prior condition) glaucoma patients. We simulated MD values (6-monthly) with true underlying progression rates drawn from the same distribution as the prior for the matched–prior condition. We estimated rates through linear regression, and determined the likelihood of obtaining this estimate as a function of a range of true underlying progression rates (the likelihood function). The maximum likelihood estimate of rate was the most likely value of the posterior distribution (the product of the prior distribution and likelihood function).

Results.: For short (4) visual field series, the matched–prior condition, unmatched–prior condition, and linear regression gave median errors (estimated minus true rate) of 0.02, 0.20, and 0.00 dB/y, respectively. Positive predictive values for determining rapidly progressing (<−1 dB/y) rates were 0.46, 0.42, and 0.38, with negative predictive values of 0.93, 0.94, and 0.95. For more extended series the magnitude of the differences between techniques decreased, although the order was unchanged.

Conclusions.: Performance shifts in Bayesian estimators of visual field progression are modest even when prior distributions do not reflect large risk factors, such as IOP treatment.

Introduction

How best to determine ongoing visual field loss in glaucoma patients remains a significant challenge. One commonly performed technique involves performing a linear regression on the summary index, mean deviation (MD). ^1–3 In addition to determining whether there is significant progression, linear regression also allows the rate of progression to be estimated, and so predictions made about whether a patient is likely to suffer substantial visual impairment within their lifetime given their current level of treatment. 

Estimates of rates typically are subject to substantial amounts of noise (variability) when few visual fields exist, becoming increasingly better estimated as a more extended series of visual field tests become available for a patient. ⁴ Recent simulation studies have demonstrated the ability of linear regression (also known as ordinary least squares) to determine significant progression ⁵ as well as define the rate of progression ⁶ as a function of the number of visual field tests performed. Such simulations demonstrate that it is largely impossible to estimate progression rates to any satisfactory degree until several visual fields have been performed. Although the presence of significant progression may be determined with multiple testing with the first two years, ⁵ it has been suggested that the rate cannot be predicted usefully before five years given reasonable testing protocols. ⁶ It is unlikely that clinicians wait for such a period before estimating the rate of visual field progression, however, and so progression estimation is made despite only a limited series of visual fields being available. Therefore, being able to improve progression rate estimators, especially when the number of visual field tests performed is low, would be of benefit. Additionally, clinical practitioners use more than just visual field information to determine the status of glaucoma patients over time, and some preliminary investigations are now incorporating additional clinical measurements and evaluations into their modelling procedures. ^7,8  

Conventional linear regression makes no assumptions regarding the true rate of progression underlying the data, despite the fact that some rates of progression are a priori more common than others. ^3,4,9 Being able to use such a priori information may be useful in constraining estimates of progression rate and, therefore, improving their utility, however. For example, a rate estimate that suggests very rapid progression after a few visual fields is more likely the result of the inherent variability of MD indices ¹⁰ rather than a truly rapidly progressing field defect. Bayesian techniques provide a formal framework in which population-based a priori information and empirical data from a patient can be combined to estimate a particular variable, for example, progression rate. They have the general property that the estimate of the variable is influenced most heavily by the population data when empirical data from the patient are scant, with the population data having decreasing influence as more empirical data accumulates. They have been used successfully in visual science to estimate visual sensitivities, ¹¹ and appear in modified forms to estimate visual sensitivities in certain clinical perimeters. ^12,13  

A priori information is quantified by the use of a prior distribution, which for glaucoma would quantify the frequency with which various progression rates are expected to occur. Risk factors, such as patient age, ⁹ the presence of bilateral field loss, ¹⁴ thin central corneal thickness, ¹⁵ or the presence of optic nerve head change ^8,16 or hemorrhages, ¹⁷ exist that modify the likely rate of progression. Therefore, prior distributions could be tailored to particular patients, based on the presence or absence of risk factors for glaucoma progression. Determining appropriate prior distributions likely requires a significant investment in time and money if obtained via population studies, however. As increasingly more specific risk factors are assessed, the proportion of the glaucomatous population having these factors will decrease and so obtaining appropriately large samples from which to determine distributions becomes increasingly difficult. Before obtaining such distributional data, it would be useful to estimate the influence of priors incorporating major risk factors for glaucoma progression to gauge the maximum benefit expected from developing priors for other, lesser glaucoma risk factors. 

In our study, we examined the benefits of having a prior distribution accounting for one of progression's major risk factors—whether IOP is treated—in a simple maximum-likelihood Bayesian estimator for progression rate. We then assess the statistical properties of our estimator, in comparison with conventional linear regression, using a simulation technique in which the true rates of visual field progression are known and the statistical properties of the simulated data are described fully. In a clinical setting, the utility of Bayesian estimators of progression rate is likely to be determined in large part by the clinician having a clear understanding of the processes underlying the Bayesian estimation technique, without which the specific advantages and disadvantages of the technique are difficult to appreciate. Our current study also aims to present the general processes underlying our Bayesian estimation technique in a way that may be understood by a nonspecialist reader, and so aid this understanding. 

Methods

We defined the true rate of visual field progression of an individual as R, with the estimate of this rate, derived from a series of visual field examinations, being r. In the current study, estimates of R obtained through linear regression (or ordinary least squares) and our maximum-likelihood Bayesian (hereafter referred to as Bayesian) method are denoted r _OLS and r _MLB, respectively. 

Prior Distribution

Our Bayesian technique requires information about the population distribution of the parameter R. We took published estimates of progression rates for MD (in dB/y) on the Humphrey Field Analyzer (HFA) to form prior distributions reflecting either treated (Canadian Glaucoma Study, ⁹ combining progressing and nonprogressing patients data) or untreated (Early Manifest Glaucoma Trial, ³ high-tension patients only) primary open angle glaucoma (Fig. 1, upper panels). We fitted these data with a modified hyperbolic secant, ¹⁸ and from this continuous function created a discrete probability mass function ¹⁹ normalized to give an area under the function of unity (Fig. 1, lower panel). Figure 1 shows that most observers with treated open angle glaucoma can be expected to have progression at a slow rate, although a small proportion might be expected to show signs of rapid progression (<−1 dB/y). As the distribution of progression rates for our simulated patients was identical to our prior distribution for treated glaucoma, the priors for treated versus untreated patients are referred to as matched and unmatched priors, respectively. 

Figure 1

 

Distribution of progression rates for treated (Canadian Glaucoma Study, upper) and untreated (Early Manifest Glaucoma Study, middle) primary open angle glaucoma, and the derived prior distributions (lower: solid line = treated, dashed line = untreated). Using the nomenclature described by King-Smith et al. ¹⁸ parameters B, C, and t for the modified hyperbolic secants describing the prior distributions are 3.55, 1.65, and 0.0076, and 6.31, 0.660, and 0.103 (treated and untreated conditions, respectively).

Likelihood Function P(r|R)

As r determined by linear regression is an imprecise estimate of the true progression rate, a given value for r is compatible with a range of underlying progression rates R. This can be quantified by generating the likelihood function giving the probability of obtaining the particular value r given a true rate of progression R (i.e., the conditional probability P(r|R)). For linear regression, the shape of this likelihood function is Gaussian with a mean of r and a variance:  

 

where σ² _MD is the variance of the errors about the linear regression line in the y-direction, x_i is the position of the ith value in the x-direction, and x̄ is the average position along the x-direction.  

Bayesian Estimation Technique

Our Bayesian technique is designed to determine the most likely estimate of a patient's true progression rate R, given an empirical estimate of progression from linear regression r _OLS combined with knowledge about the likely values R can take in the population (as given in the prior distribution). In the absence of any patient information (i.e., prior to visual field testing), the most likely value of progression rate is given by the peak of the prior distribution (= −0.01 dB/y, Fig. 2A). In the left panels, linear regression of MD values for three sequential visual fields in a particular patient provides an estimate of their progression rate r _OLS (Fig. 2B). The Equation can be used to determine the likelihood of obtaining this estimate given a true underlying progression rate R (Fig. 2C). According to Bayes theorem, ²⁰ the prior distribution and the likelihood function then are multiplied together to form a posterior distribution (ignoring the presence of a normalizing constant), which gives the probability of R given our finding of r (the conditional probability P(R|r), Fig. 2D). The most likely value for R is at the peak of the posterior distribution (Fig. 2D, solid vertical line). The 95% credible intervals around this maximum likelihood estimate are given by determining the R values corresponding to the 2.5% tails of the distribution (Fig. 2D, dashed vertical lines) via integration – the probability that R lies between these limits is 0.95. Intuitively, if the likelihood function, which is related to measurement uncertainty, is broader than the prior distribution, it will add very little information and so the most likely estimate of rate will be determined almost exclusively by the population-derived information present in the prior distribution. 

Figure 2

 

Demonstration of our Bayesian model for a series of 3 visual fields (left) and 8 visual fields. (A, E) Show the prior distribution, being the distribution of visual field progression rates in the population. (B, F) Show the visual field series for a patient, along with an ordinary least squares linear regression of these data to give the estimated rate of progression r _OLS. (C, G) Give the likelihood functions, which quantify the likelihood of obtaining the value r _OLS given a true underlying rate of progression R. (D, H) Give the posterior distribution, which quantifies the probability of R given the particular value r _OLS and the prior distribution: it is formed by multiplying together the prior distribution and the likelihood function, according to Bayes theorem. The most likely value for R in our Bayesian model (= r _MLB) is given by the peak of the posterior distribution (solid vertical line). The 95% credible intervals around this estimate are produced by integrating under the posterior distribution to find the 2.5% and 97.5% tails (vertical dashed lines).

If a further five visual fields were obtained on the same patient (Fig. 2, right panels), a new estimate of progression rate r _OLS can be obtained and the likelihood function again derived. As predicted by the Equation, as the number of visual field estimates increases the width of the likelihood function will tend to decrease (Figs. 2C versus 2G) and so its influence on determining the peak of the posterior distribution increases. In the limiting case, an infinite series of visual fields would result in an infinitely narrow likelihood function and an infinitely narrow posterior distribution, both peaking at the value R = r _OLS. It should be noted, however, that this decrease in the likelihood function width is only what is expected on average: had the first three visual fields in Figure 2 happened to fall almost exactly on a straight line by chance, the likelihood function after three fields may have been narrower than that obtained with a more extended series. 

Simulation Details -- Visual Field Series

For each simulated patient, we created MD values for 13 visual fields performed at 6-month intervals. We selected the true rate of progression R for each series by a pseudorandom number generator (Matlab “random” function, Matlab R2010b for Macintosh; MathWorks, Natick, MA) with a frequency as given by the matched-prior distribution. The MD estimate for a given visual field then was generated by taking the height of a line with slope of R (in dB/y) and intercept of zero, and jittering it by a random value selected from a normal distribution with a particular standard deviation. For our main simulation, this standard deviation was 1.0 dB and is the same as the moderate variability condition assumed in previous simulation studies. ⁶ An assumption of linear regression is that variability is constant for all measures in the series (homoscedasticity), and use of a fixed standard deviation jitter allows us to compare the statistical performance of linear regression and our Bayesian estimator in the absence of assumption violations for either technique. Assuming a constant variability also allows us to ignore the absolute magnitude of our MD measures and so analyze rates in isolation, meaning that we do not need to consider the distribution of MD values in addition to the distribution of progression rates. MD variability does, however, increase as a function of the defect depth in clinical data, ¹⁰ although the increase is smaller than that seen in individual sensitivity measurements. ²¹ Therefore, we performed an additional simulation where this increase was modelled in the simulated patient data: variability linearly increased from a standard deviation of 0.5 dB at an MD of 0 dB through to 1.2 dB at −10 dB, and remaining constant outside of these limits, approximating the trend seen in full threshold fields on the HFA. ¹⁰ We selected baseline MD values (i.e., those for the first visual field in each series) using the histogram frequencies reported in the Canadian Glaucoma Study data, applied to the midpoint MD values of each histogram bin, ²² giving baseline MD values ranging from +1 to −13 dB. In contrast to our simple model, MD variability should reduce at very high values of MD as an increasing proportion of the field shows absolute field defects (sensitivity = 0 dB). The frequency of such fields in our modelled data set is rare, however. 

For each series of visual fields, we generated estimates r_OLS for 10 subsets of the series (the first three visual fields, through to the first 12 visual fields) and also for the complete visual field series. For each simulation we generated 30,000 visual field series. In our simulations, all distributions and likelihood functions were represented discretely in 0.01 dB/y wide bins between extremes of −10 to +10 dB/y. 

Results

Figure 3 shows the distribution of errors (r minus R) in the estimates returned from linear regression (upper panel), and the matched-prior Bayesian (middle panel) and unmatched-prior Bayesian (lower panel) techniques, as a function of the number of fields in the series (two visual fields per year). In comparison with the matched-prior Bayesian estimator, linear regression produces larger errors for short series of fields; this is highlighted in the upper right-hand panel, which gives the difference between each method's 25%, 75%, and 95% confidence intervals (97.5% − 2.5%) limits, with positive values denoting smaller limits for the Bayesian estimator. By comparison, the difference between the two Bayesian estimators (right-hand lower panel) is small with each producing similar confidence limits (dashed line). The matched–prior condition gives smaller 25% error limits, but poorer 75% limits, consistent with the difference in bias, or median error, between the two techniques (circles, left panels). For series of nine visual fields or greater, performance of all techniques largely is indistinguishable. This suggests that a potential benefit of the Bayesian estimators is that a more reliable assessment can be achieved when the number of visual fields is limited. 

Figure 3

 

Errors in estimating the true progression rate R, for our linear regression (upper) and matched–prior and unmatched–prior Bayesian method (middle and lower, respectively) techniques, assuming a visual field test every six months. Dashed outer lines represent the 2.5 and 97.5% limits of the error data, with triangles giving the 25 and 75% limits, and circles the median. Upper right shows the difference between the matched–prior Bayesian and linear regression 25% limits, 75% limits, and 95% confidence intervals (97.5% limit − 2.5% limit; up triangles, down triangles, and dashed line, respectively). Positive values mean that the Bayesian estimator had smaller limits. The vertical and oblique crosses give the difference between 75% limits when MD variability was 2.0 and 0.5 dB, respectively. Lower right shows a similar analysis of the difference between the matched– and unmatched–prior Bayesian methods. Rerunning the simulation (matched–prior condition, MD variability = 1.0 dB) produced a median change across all limits of 0.003, 0.000, and 0.000 dB/y (r_OLS , matched-prior r_MLB and unmatched-prior r_MLB , respectively) and no change greater than 0.04, 0.04, and 0.02 dB/y.

Although the bias in our Bayesian estimators is comparatively small overall, Figure 4 shows it is related to the underlying slope R and can be very large for very small field series. Specifically, slopes from the matched-prior Bayesian method are reduced systematically for rapidly declining and rapidly improving MD values, as such rapid changes are very unlikely given the prior distribution. The bias reduces as the number of fields in the series increases, becoming effectively zero for large series (lower right panel). Figure 4 also shows that the width of the 95% confidence intervals for errors also is a function of R for the matched-prior Bayesian estimator, with errors sometimes being greater than those for linear regression for infrequently encountered rates of progression. However, for the most common rates of progression (R between −1 and 1 dB/y = 84% of results), errors are consistently narrower for the Bayesian estimator. 

Figure 4

 

Error in slope estimate r as a function of the true slope R, for linear regression (open symbols) and for the matched–prior Bayesian technique (closed symbols). Circles give the median error, with triangles showing the 2.5% and 97.5% limits. On the x-axis, data points are plotted at the midpoint of the 1 dB/y wide bins used to group the data. Expected frequencies in each bin, based on the distribution shown in Figure 1, are shown in the lower right. However, for this simulation only, R values were distributed evenly across the analyzed range to allow approximately even numbers of data in each bin.

Figure 5 shows the positive (PPV) and negative (NPV) predictive values for each technique for determining moderate (<−0.5 dB/y) and rapid (<−1.0 dB/y) progression among the simulated cohort. Both Bayesian methods produced an improvement in positive predictive values compared to linear regression despite little decline in negative predictive values, particularly for rapid progression with short series of fields. For rapid progression, use of an unmatched-prior changed positive predictive values by less than 0.05. Rerunning of the matched–prior condition with MD variability increasing with MD produced similar shaped curves, except with higher PPV (dotted lines) and NPV (not shown: median increase = 0.006, maximum increase = 0.02 compared to the matched-prior, fixed variability condition) curves. The average value for σ_MD was 0.845 for this simulation, being slightly lower than the fixed value of 1.0 used in our main simulation. Reducing MD variability to a fixed value of 0.845 in our matched–prior condition reduced the discrepancy between simulations where MD variability either was fixed or variable (maximum PPV discrepancy of 0.013 and 0.022, and maximum NPV discrepancy of 0.005 and 0.003, for the <−0.5 dB/y and <−1.0 dB/y conditions respectively), indicating that the change in predictive indices seen when MD variability increased with MD largely is due to the reduced average value for σ_MD under this condition. The influence of alterations in MD variability can be seen in the error limits given in Figure 3 (upper right panel), where reducing variability to 0.5 dB marked reduced errors (oblique crosses) while raising it to 2.0 dB markedly increased errors (vertical crosses), relative to simple linear regression. Such changes further highlight how measurement variability has a large effect in determining whether using prior-knowledge helps to estimate progression rates—when measurements are most noisy, the benefits of prior knowledge are greatest. 

Figure 5

 

PPV and NNV for two progression rates. The matched–prior Bayesian, unmatched–prior Bayesian and linear regression are shown by the closed symbols, solids lines, and unfilled symbols, respectively. PPV results for when a matched-prior was used, but when the variability of the MD increased with the magnitude of the MD, are shown by the dotted lines.

Discussion

Our results showed that failure to consider a major risk factor for visual field progression—whether or not IOP is treated—did not alter dramatically the performance of a Bayesian estimator of visual field progression. Therefore, it would be anticipated that priors incorporating lesser risk factors will have even less influence on the performance of Bayesian estimates of visual field progression, and so undertaking population-based studies specifically aimed at quantifying the shape of these priors may not be justified. Despite this, use of a prior reflecting some knowledge that very rapid progression rates are relatively uncommon can improve the estimation of the rate of visual field progression compared to linear regression, especially when the number of visual fields in a series is low. Although our Bayesian method can introduce a variable bias into the estimation of progression rate r, this bias tends to be smallest for the most common range of progression rates (Fig. 4) and rarely results in errors that are greater than those expected from conventional linear regression. This bias does have the potential to underestimate very rapid rates of progression, although our observation that NPV are barely altered by our Bayesian technique (Fig. 5) suggests that this effect is of limited magnitude. There is evidence that even experienced clinicians can over-call visual field progression in some circumstances, ^23,24 and so the improvement in PPV seen with our Bayesian technique, in comparison with linear regression (Fig. 5), may be of use in reducing such false alarms. Our Bayesian technique is conceptually simple and so may allow its advantages and limitations to be more readily grasped by nonspecialists, in comparison with comparatively more complex Bayesian models designed to detect significant progression ⁷ or to incorporate structural data into progression rate estimates. ⁸  

A prior distribution should be best matched to the presumed population distribution wherever practical, and so how best to describe these priors mathematically must be considered. The existence of several large population-based studies in glaucoma ^25–27 provides an opportunity to examine the general shape these distributions may take. Our modified hyperbolic secant is able to capture the asymmetric tails in glaucoma progression distributions seen in some of these studies although there is no theoretical underpinning for this distribution's use, and other functions may well provide better fits. Showing that one distribution is a statistically better fit than another often is very difficult, however. For example, Spry & Johnson failed to show significant or meaningful differences between various candidate models to describe age changes in perimetric data, despite data from over 560 participants, ²⁸ The existence of three parameters in our modified hyperbolic secant that alter distinct aspects of the distribution—in particular, the slope of the fall-off for each tail and the mode of the distribution ¹⁸ —may have practical advantages when incorporating information into prior distributions. For example, Heijl et al. found that untreated pseudoexfoliative glaucoma could produce devastatingly rapid progression in excess of −10 dB/y, yet untreated normal tension glaucoma never produced progression rates over −5 dB/y, and so adjusting the slope of one tail of the prior distribution could incorporate such information relating to glaucoma type. ³ In addition, some literature relating to risk factors for progression do not report distributions but rather changes in average progression rates. ^14,29 In the absence of any other information, these changes might be modelled by adjusting the parameter that determines the mode of the modified hyperbolic secant. For most clinical situations, however, the prior distribution or the precise risk factors likely will not be known with confidence, and so our simulation results provide some reassurance that the performance characteristics of Bayesian estimators are robust to modest mismatches between the prior and the true distribution of progression rates in the population. Our results showed that the largest changes in performance occur when going from having no prior (ordinary least squares) to some sensibly shaped, if only approximate, prior distribution. 

Estimating Variability in MD

It is likely that some individuals will have MD variability significantly greater or smaller than the average value for the population. While the value σ_MD , derived from the linear regression (see the Equation), provides an estimate of MD variability in an individual, this estimate is itself very noisy when the number of visual fields available is small. It would be possible to modify our Bayesian estimator such that a population-based value for σ_MD is used until there is evidence an individual's value for σ_MD is significantly different from this population value. This value would depend upon the type of test being performed: MD variability for SITA-Fast, ³⁰ for example, appears somewhat larger than typical values for full-threshold fields. ¹⁰ Rather than a dichotomous selection of σ_MD (population versus individual) based on significance testing, the value σ_MD itself could be estimated using a Bayesian technique that is updated as each new visual field is collected, with the most likely value taken as the value for σ_MD in the above estimator, or σ_MD and r _MLB simultaneously optimized. Similar two-dimensional Bayesian parameters estimations already have been described for psychophysical estimation of sensitivity. ³¹ It would need to be demonstrated that the additional complexity of such a technique provides a useful improvement, however. This may not be the case, as reliably estimating individual values for σ_MD likely requires longer series of visual fields than for reliably estimating progression rates. These series almost certainly would be longer than the limit beyond which Bayesian estimators cease to show any advantage over linear regression. An analogous situation regarding the difficulty in estimating variability is seen when determining frequency-of-seeing curves, where the central position of the curve (the threshold) can be estimated in fewer trials than the variability (slope) of the curve. ³¹  

General Applicability of Bayesian Methods

Although the Bayesian method described in our report was developed for linear regression of the summary index MD, the general principals could be applied to other measures of glaucomatous progression. Most clearly related would be point-wise estimates of progression rates through linear regression, where prior information of the likely rates of progression at a particular point could be used to constrain these estimates. Prior distributions giving the likelihood of progression at a point given the status of surrounding points (e.g., level of sensitivity loss, estimated progression rate) also could be developed. Furthermore, priors relating to structural data from imaging devices might also provide useful constraints on the likely distribution of progression rates. ⁸ Nonlinear regression methods also could be accommodated, ³² provided the nature of the likelihood function is known. In addition, it is not required that visual fields be measured at regular intervals. Recommendations already exist for how nonuniform spacing of field tests may improve progression detection, ^5,33,34 some of which themselves are based on the application of Bayes' theorem, ³⁵ and our method can accommodate such spacings. 

While our study shows that Bayesian estimators of visual field progression rates demonstrated several advantages over simple linear regression in many situations, even when the prior distribution is not an exact match to the true patient distribution, it is important to remember a Bayesian estimate of rate will not be ideal for all individuals in all situations. Importantly for glaucoma, the typically small number of very rapidly progressing patients will have their progression rates underestimated when few visual fields are available, as it is more likely that such rapid progression reflects variability in the data rather than a true underlying progression rate (Fig. 4). Discrepancies between rate estimates can be denoted appropriately by having simple linear regression and Bayesian progression rates always available simultaneously to the clinician, possibly augmented by a qualitative comment—akin to the plain language analysis seen in the Glaucoma Hemified Test ³⁶ —describing the nature of any disagreement. Once detected, any discrepancy then can be acted upon if required, for example, by collecting addition visual fields at more closely spaced intervals if other information (e.g., intraocular pressure control, optic nerve head changes) suggests that true rapid progression is occurring. Clinicians and basic scientists can make such informed decisions only when they have a high level of familiarity with how a Bayesian estimation method works, as well as with the technique's advantages and disadvantages. 

Acknowledgments

The authors thank Paul Artes and Neil O'Leary for informative discussions relevant to this study. 

Footnotes

 Supported by Australian Research Council Future Fellowship FT120100407 (AJA).

Footnotes

 Disclosure: A.J. Anderson, None; C.A. Johnson, Haag-Streit (C)