Estimating the True Distribution of Visual Field Progression Rates in Glaucoma

Abstract

Purpose.: Bayesian methods allow the distribution of glaucomatous visual field progression rates in the population to constrain individual progression rate estimates. As the true population distribution is never known, it must itself be estimated from a finite number of noisy individual estimates of rate. We used simulations to investigate the relationship between the true distribution of progression rates and that estimated from noisy empirical data.

Methods.: We generated series of visual fields (3–10 per patient) using different variabilities (SD of 0.5–2.0 dB) for the summary index mean deviation (MD) to determine the relationship between the distribution of empirical estimates of progression rates determined by linear regression and the true underlying distribution of progression rates.

Results.: Estimating the underlying distribution from empirical data produced biases that broadened the distribution and made it more symmetrical, particularly for a short series of variable visual field estimates. Decreasing cohort size increased the variability in distribution parameter estimates, but produced no bias. Variability in distribution tails produced a 3.5-fold variation in the proportion of rapid progressors for the smallest cohort (200), falling to 1.8- and 1.3-fold for cohort sizes of 800 and 3200, respectively.

Conclusions.: The underlying distribution of glaucomatous visual field progression rates for the population is likely to be narrower, and less symmetric, than that predicted from empirical data. Therefore, care should be exercised when inferring the benefits of Bayesian estimators, particularly where prior information is itself derived from a small sample of noisy empirical estimates.

Introduction

The overall rate of decrease in visual field sensitivity in glaucoma, in decibels per year, can be estimated by performing a linear regression of the perimetric summary index mean deviation (MD) over time.^1–3 Such estimates are critical in establishing the efficacy of treatment and in predicting the likelihood that a patient will suffer profound vision loss within their lifetime. These estimates are noisy however, particularly when only a few visual field tests are available.⁴ Given this noise, Bayesian methods can help constrain rate estimates to sensible values and so potentially increase the use of these estimates.^5–8 

Conventional Bayesian methods require the prior distribution to be specified, which quantifies the frequency with which various progression rates are expected to occur in the population.⁵ Information from this prior is then combined with a patient's visual field data in order to produce an estimate of progression rate. Strictly, the prior distribution should be the distribution of the true underlying rates of progression. In practice, the underlying rate of progression in an individual, and hence also in the population, is never fully known. Because of this, any realisable prior distribution is subject to at least two potential sources of error. Firstly, the distribution is for the estimates of rate (r) rather than the true underlying rates (R), and estimates of rate are necessarily noisy. Secondly, the distribution is not that of the population, but rather a sample of the population. This sample may be taken from a previous study,⁵ or from the data itself in the case of Empirical Bayes methods.⁷ 

The shape of the prior distribution in glaucoma is typically asymmetric and with a relatively long tail for negative progression rates,^1,9 indicating that some patients have very rapid rates of visual field decline yet almost none show a rapid improvement. The shape of the distribution can be approximated by a modified hyperbolic secant,⁵ which has the properties that the width of either tail can be specified, along with the mode of the distribution.¹⁰ Various risk factors alter the rate of glaucoma progression,^6,9,11–14 and it has been proposed that systematic shifts in the parameters of a modified hyperbolic secant might be used to incorporate these risk factors into Bayesian prior distributions.⁵ Recent work has shown that the distribution of progression rates appears to shift as the risk factor of age changes.¹⁵ In order to know whether such shifts are significant it is necessary to determine the variability in the form of the prior that result when a population distribution is estimated from a finite sample of participants. The size of this sample may be substantially smaller than the nominal cohort size if an analysis on subgroups is performed. Furthermore, estimates of rate r used in previous distributions^1,9,15 are subject to random variation, which is assumed to be Gaussian when conventional linear regression is used. The influence of this random noise on an asymmetric distribution needs to be determined, although it might be expected that when rate estimates are highly variable this symmetrical Gaussian noise may act to make the distribution of progression rates appear artifactually symmetrical. 

In this paper, we use Monte Carlo simulations to explore how well the distribution of clinical progression rate estimates r predict the true underlying distribution of rates R, the prior distribution. We also use a statistical bootstrapping procedure to determine the magnitude of errors introduced when attempting to generate a prior distribution from a finite sample of clinical data. 

Methods

Reference Distribution

We estimated the distribution of visual field progression rates (dB/y) in glaucoma from published data of 583 patients with POAG and pseudoexfoliation examined in a hospital setting.¹⁶ Patients had to have a repeatable visual field defect at diagnosis, and be followed for more than 5 years with more than five visual field examinations performed over this period.¹⁶ We fitted these data with a modified hyperbolic secant¹⁰ using an maximum-likelihood procedure. Using the nomenclature of King-Smith et al.,¹⁰ the critical parameters describing the modified hyperbolic secant are B, C, and t, which describe the right-hand tail slope, the left-hand tail slope, and the mode, respectively. The best fit can be seen in Figure 1, which for the purposes of this paper is referred to as the reference distribution. During all iterations of the fitting procedure both here, and for the least-squares method described below, distributions were scaled to have an area of unity under the curve (evaluated between the limits of −10 and +10 dB/y) using the vertical scaling parameter A.¹⁰ Because of this, A was not a free parameter and distributions were fully described by parameters B, C, and t alone. 

Figure 1

 

Distribution of visual field progression rates for patients with POAG and pseudoexfoliation glaucoma, derived from data reported by Heijl et al.¹⁶ Data are fitted with a modified hyperbolic secant with parameters B, C, and t (using the nomenclature described by King-Smith et al.¹⁰) of 4.36, 1.33, and −0.39, respectively. Coefficient of determination (R²) for the fit is 0.98. Histogram bins for the modified hyperbolic secant were identical to those for the empirical data (0.2 dB/y wide bins), but are plotted as reduced bin widths to allow the data and the fit to be simultaneously visualized.

It should be noted that there is no theoretical reason for using a modified hyperbolic secant, and that other skewed distributions might well model the data similarly well. For the current paper, of primary importance is that such mathematic models provide a parametric description of the data, allowing shifts in the distribution to the readily quantified by determining the magnitude of shifts in the model's parameters. Mathematic models are also continuous functions, thereby allowing frequencies (e.g., the proportion of fast versus slow progressors) to be calculated in a way that is not constrained by the coarseness with which the empirical data is binned. 

Using the Distribution of Empirically Determined Rates to Infer the Distribution of the True Underlying Rates

As described previously, the true rate of progression R is generally unknown and so the distribution of R is estimated by examining the distribution of r, the empirical estimate of R determined by linear regression of a series of visual fields.⁵ We examined the relationship between these two distributions using Monte-Carlo simulations. The true distribution of R was set to be the same as the reference distribution. We then simulated longitudinal visual field data for a large number (200,000) of simulated patients using a method described previously⁵: in brief, true progression rates were selected with a frequency as described by the true distribution of R, and MD values generated at each timepoint by applying Gaussian jitter of a particular SD to the MD that was expected under noise-free conditions. This SD was fixed for each set of 200,000 simulated patients. Rates of progression r for these data were then determined by linear regression, and these rates then placed in 0.1-dB/y wide histogram bins. The resulting distribution was fitted with a modified hyperbolic secant for binned rates of −10.0 to +10.0 dB/y, using an iterative least squares method. We examined the influence of the number of visual fields (3–10 visual fields, at yearly intervals) and of MD variability (low variability [SD = 0.5 dB], moderate [1.0 dB], and high [2.0 dB], consistent with previous definitions^4,17). 

Bootstrapping to Determine Variability in Distribution Parameters

We determined the variability of the parameters defining glaucoma progression rate distributions using a bootstrapping procedure. The height of the reference distribution was appropriately scaled and rounded counts at each 0.1 dB/y midpoint used to generate a simulated cohort of patients. From this initial cohort, 10,000 new cohorts were generated via sampling-with-replacement with each then fitted with a modified hyperbolic secant using an iterative least squares method. We then calculated the 2.5% and 97.5% limits from the resulting distribution of the parameters B, C, and t defining the modified hyperbolic secant. We investigated the influence of doubling the cohort size, from a starting size of 100 through to 12,800. 

Results

Figure 2 shows how well distributions of empirical estimates r predict the true underlying, and generally unknown, distribution of rates R, for three different levels of variability. Variability in empirical estimates (quantified by σ_MD, the SD of the MD values) produced a significant bias in the distribution parameters (left hand panels) which was greater for higher levels of variability and when the number of visual fields in a series was small. The right hand panels give an indication of what change in shape of the modified hyperbolic secant results from such parameter biases: large amounts of variability and/or a small number of visual fields tends to broaden the distribution and also make it more symmetrical. For observers showing moderate (1.0 dB) or low (0.5 dB) levels of variability, distributions well approximated the true underlying distribution when series of approximately six or more visual fields were available. 

Figure 2

 

The influence of the number of visual fields measured on a patient (in yearly intervals) on the distribution of progression rates, for several MD variability levels. Distributions are from modified hyperbolic secant fits to the data from 200,000 simulated visual field series. Left hand panels give fit parameters for three different levels of MD variability (MD SD [σ_MD] of 0.5, 1.0, and 2.0 dB for downward-pointing triangles, circles, and upward-pointing triangles, respectively). Right hand panels show distributions for 3, 6, and 10 visual fields (thin curves of decreasing half width) as a function of MD variability, with the true underlying distribution given by the thick lines. Frequencies for the modified hyperbolic secant were calculated in 0.1-dB/y wide bins and are represented by the value of the function at the bin midpoint: these points are linked by lines, for ease of visualization.

Figure 3 shows the variability of estimates in the three parameters defining the modified hyperbolic secant as a function of the size of the cohort. For all parameters, variability decreased as cohort size increased. However, in contrast to Figure 2, there was almost no bias in the estimates, indicated by the median being almost indistinguishable from the expected values from the reference distribution (horizontal dashed lines) even for small cohorts. 

Figure 3

 

Parameters for a modified hyperbolic secant fitted to bootstrapped data, for cohort sizes ranging from 100 to 12,800. Symbols give median parameters values, with error bars showing the 2.5% and 97.5% limits. Horizontal dashed lines give the parameter value of the reference distribution used to seed each bootstrapped cohort.

Figure 4 shows the changes in shape of the modified hyperbolic secant expected from such parameter variability. The right hand panels show distributions created from parameters B (right hand tail) and C (left hand tail) at the 2.5% and 97.5% (thick and think lines, respectively) limits given in Figure 3, compared with the reference distribution. Parameters at the 2.5% limit produced distributions that underestimated the peak height of the distribution and overestimated the width of the tails, with parameters at 97.5% producing the opposite effect. As cohort size increased, the magnitude of these discrepancies decreased. Of particular interest are differences in the proportion of people showing rapid progression (≤−2 dB/y¹⁷). With a cohort of 200, there was a 3.5-fold variation in the proportion of rapid progressors between the 2.5% and 97.5% distributions, with this falling to 1.8- and 1.3-fold for cohort sizes of 800 and 3200, respectively. Variations were larger for the proportion showing small amounts of visual field improvements (≥0.5 dB/y). Figure 4 also highlights the variation in individual histogram frequencies expected simply through sampling variability (left hand panels, thin error bars). For a small cohort size (n = 200) the width of 95% limits were larger than the nominal histogram frequency, and there was substantial variation in the proportion of patients that could appear in the distribution tails. 

Figure 4

 

Results of the bootstrapped data, showing the variability in the frequency data and in the modified hyperbolic secants, as a function of cohort size. Gray histogram bars give the median frequency for the bootstrap. For the frequency data (left panels), vertical error bar show 2.5% to 97.5% limits. For the modified hyperbolic secants (right panels), thick lines show functions generated from the 2.5% limits for parameters B & C (which change the slope of the right and left tails, respectively), with the thin lines showing functions generated from the 97.5% limits of these parameters (see Fig. 3 for parameter values). As such, the functions indicate how the tails of the distribution change at these limits, rather than a fit to a particular dataset. Parameter t was fixed to be the same as the reference distribution from which the bootstrap data was generated. The representation of the modified hyperbolic secants is as described in Figure 2.

Discussion

Our results highlight two sources of variability that affect our ability to quantify the true distribution of visual field progression rates in glaucoma: that the rates analyzed are empirical estimates of the true underlying rate, and that studies necessarily only provide a sample of the population. Estimating the true, underlying distribution of rate in glaucoma from empirical data produces biases, particularly when attempting to do so from only a limited number of visual fields and/or when data is variable (Fig. 2). These biases both broaden the estimated distribution, and tend to make it more symmetrical. When sample sizes are small, estimates of distribution are poor and may result in greater than a 3-fold change in the estimated proportion of fast (≤−2 dB/y) progression rates (Fig. 4). Parametrically quantifying empirical distributions allows the variability in fitted parameters to be determined (Fig. 3) and so whether changes in distributions, such as may occur with changes in risk factors¹⁵ or glaucoma type,¹ are significant given the sample size used. Such limits could also be used to prospectively estimate the sample size required to find a distributional change of a particular magnitude. 

Our study examined two sources of variability independently, and indeed this is one of the strengths of performing simulations. This assumption of independence is reasonable in terms of bias: sampling produces uncertainly about a distributions parameters but does not introduce any bias (Fig. 3), whereas the use of empirical rates does produce biased estimates (Fig. 2). Real datasets will, however, be influenced by both sources of variability we analyses. For example, our results in Figure 2 were obtained from a very large dataset (200,000 series of visual fields) and so show the influence of estimating the true underlying distribution of rates from empirical data when there are essentially no limitations on study sample size. However, empirical studies reporting distributional data commonly have cohort sizes markedly smaller than this (e.g., 2300,¹⁵ 583,¹⁸ and 583 patients¹⁶), with even smaller sizes resulting when distributions of subgroups are analysed.^1,15 Such studies will therefore be subject to both uncertainties (but not bias) due to sampling (Figs. 3, 4) as well as biases from using empirical estimates (Fig. 2). 

Study Limitations

Like all simulations, our study has limitations. We used a fixed variability for MD, although there is evidence that MD variability scales with the magnitude of MD itself.¹⁹ Incorporating this relationship in previous Monte-Carlo investigations did not alter any of the main finding found using fixed levels of MD variability, however.⁵ Potentially of greater importance is that our study examines what happens when the variability of the entire cohort alters (Fig. 2). We believe that there are several reasons why such overall changes in variability might exist. Different studies may use different perimetric tests whose variability characteristics therefore differ (for example, SITA Standard versus SITA Fast²⁰). Some study cohorts may represent a highly selected and trained group in a clinical trial, with strict guidelines of acceptable perimetric performance,²¹ whereas other cohorts may represent less selected people from within a clinical population.¹⁶ In addition, robust regression methods may be used to estimate progression rates in some studies¹⁵ which would be expected to reduce variability compared with studies using ordinary least squares methods. Differences in MD variability between individuals still remain in all of the aforementioned situations, although there is evidence that it is a shift in the average MD variability in the group that is most influential on population-based statistical measures.⁵ 

We also assumed that the retest distribution of MD values for an observer was symmetrical, and so we modelled it by a normal distribution as done previously.^4,5 The index MD represents an average, albeit a weighted one, of the total deviation values (i.e., the difference at each test location between the patient's sensitivity and the median age-expected sensitivity). Even though the retest distribution of total deviations can be clearly skewed,²² one might anticipate this skew not to manifest in MD values because of the central limits theorem, which states that averages should be approximately normally distributed regardless of how skewed the data being averaged is. Indeed, data from Artes et al.¹⁹ suggests there is little skew in MD distributions (e.g., for fields with an MD of −15 dB, 90% probability limits for variability are approximately −2 dB to +1.6 dB). The cause of the small residual skew is not clear, although it may be due to the underlying assumptions of the central limits theorem not being fully met for perimetric data, or it may represent an artifact of the longitudinal method used to infer these variability data.¹⁹ Either way, the magnitude of any residual skew is small and so is unlikely to materially affect the findings of this study. 

A further potential limitation is that we performed our fitting using a least squares minimisation procedure. Such procedures assume that the data variance is independent of the magnitude of the measure (homoscedasticity), that errors are normally distributed, and that data points are independent.²³ None of these assumptions are perfectly met: variability increases with histogram percentage (Fig. 4), errors are asymmetrically distributed (Fig. 4), and the percentage in one histogram bin necessarily influences the percentages in others. Despite these limitations, our fitting procedure is capable of producing unbiased results that converge on the true underlying parameters values (Fig. 3). 

Summary

Our study shows how changes in cohort size influences our ability to parametrically quantify empirically determined progression rate distributions. Such empirical distributions are an estimate of the true, underlying distribution of glaucomatous visual field progression rates in the population, and hence are used to approximate the prior distribution used in Bayesian methods for improving visual field progression rate estimates in individuals. Our results show that the true, underlying distribution of glaucomatous visual field progression rates for the population is likely to be narrower, and less symmetric, than that predicted from empirical data, however. Care should therefore be exercised when inferring the benefits of Bayesian estimators. Our study provides a way for authors of Bayesian methods to critically appraise how well their prior is specified with respect to the number of visual fields, the variance of MD, and the number of patients used to quantify the distribution. We would encourage authors to report this critical appraisal when describing the results of their Bayesian methods. It is perhaps fortunate that poorly specified priors (i.e., ones derived from a small sample of noisy empirical estimates of rate) are wider than the true underlying distribution, and so Bayesian analyses using such priors should not artifacturally show an improved ability to estimate progression slopes based on this factor alone. However, the improvement expected from Bayesian methods may be overstated if, when developing the method, the prior is tightly matched to the empirical data analyzed and a similarly close matching cannot be achieved when the method is employed in a wider clinical setting.²⁴ 

Acknowledgments

Supported by grants from Australian Research Council Future Fellowship FT120100407 (Canberra, ACT, Australia). 

Disclosure: A.J. Anderson, None