**Purpose.**:
In order to reduce noise and account for spatial correlation, we applied disease mapping techniques to visual field (VF) data. We compared our calculated rates of progression to other established techniques.

**Methods.**:
Conditional autoregressive (CAR) priors, weighted to account for physiologic correlations, were employed to describe spatial and spatiotemporal correlation over the VF. Our model is extended to account for several physiologic features, such as the nerve fibers serving adjacent loci on the VF not mapping to the adjacent optic disc regions, the presence of the blind spot, and large measurement fluctuation. The models were applied to VFs from 194 eyes and fitted within a Bayesian framework using Metropolis-Hastings algorithms.

**Results.**:
Our method (SPROG for Spatial PROGgression) showed progression in 42% of eyes. Using a clinical reference, our method had the best receiver operating characteristics compared with the point-wise linear regression methods. Because our model intrinsically accounts for the large variation of VF data, by adjusting for spatial correlation, the effects of outliers are minimized, and spurious trends are avoided.

**Conclusions.**:
By using CAR priors, we have modeled the spatial correlation in the eye. Combining this with physiologic information, we are able to provide a novel method for VF analysis. Model diagnostics, sensitivity, and specificity show our model to be apparently superior to current point-wise linear regression methods. ( http://www.anzctr.org.au number, ACTRN12608000274370.)

^{ 1 }The most common method for monitoring changes in the VF is through Standard Automated Perimetry (SAP).

^{ 2,3 }Progression is judged by looking at a series of SAP VF outputs over time. Due to the large amounts of variation in the VF it is difficult to detect true change from change due to measurement error and random variation.

^{ 4–6 }Many analytical methods exist to help determine progression. Currently, however, there is no consensus on the gold standard for analyzing VF data. This often leads to investigators developing their own methods to suit their particular study.

^{ 7 }

^{ 8 }reviewed literature prior to April 2009 and found that 47% of studies using quantitative methods calculated rates of progression, while the remainder dichotomized the outcome. The main methods can be divided into three categories: scoring methods, glaucoma change probability analysis, and point-wise linear regression (PLR) analysis. Other methods investigated include the VF index, neural networks, and spatial filters. A comparison of progression detection methods using a meta-analysis

^{ 7 }indicates that methods predicting the highest rates of progression included PROGRESSOR (Medisoft Ltd., Leeds, UK) (a form of PLR), the Threshold Noiseless Trend program, and a qualitative clinical method. The majority of the PLR methods using indices (e.g., mean deviation and VF index methods), and the PLR on threshold values with tighter significance criteria gave lower estimates of progression. However, due to the lack of a gold standard it is also hard to compare methods to determine whether one is better than another. Vesti et al.

^{ 9(p3878)}states that “The ideal method for analyzing VF change should be sensitive, detect progression with few examinations, maintain high specificity, and be resistant to fluctuation,” and it has been noted that a method which identifies the rate of progression is of more benefit in determining the future of vision loss, treatment, and thus quality of life in glaucoma patients.

^{ 10 }

^{ 11 }While this method succeeded in reducing variation, it was also less sensitive,

^{ 12 }and did not take into account the physiology of the eye, in that the VF locations along the temporal horizontal midline do not map in adjoining locations on the optic disc. Another spatial filter was created to explore spatial correlations in VFs in order to reduce noise. This was based upon the covariances between different points in the eye.

^{ 13 }While this method seemed to increase sensitivity and specificity in simulated data,

^{ 13 }it performed similarly to PLR when tested on longitudinal data.

^{ 14 }

^{ 15 }In the following sections, we present our methodology, results from fitting our model to data from the Lions Eye Institute trial registry and Vein Pulsation Study in Glaucoma (Lions Eye Institute, University of Western Australia), our discussion, and conclusions.

**Figure 1.**

**Figure 1.**

^{ 17 }Appropriate adjustments were, therefore, made to the observed field thresholds of participants using both techniques. All subjects had undergone at least two VF tests before participating in the trial, thus, minimizing the learning effect. VF output from left eyes were inverted allowing the same modeling as right eyes.

^{ 18,19 }for the underlying mean process to describe the dependence between regions sharing a common boundary. For the most part a similar principle applies to the VF measurements, since adjacent grid cells in the VF usually correspond to proximate areas on the optic disc. However, there is more complexity across the (horizontal) midline. Specifically, at the temporal side of the eye, pairs of adjacent cells with one above the midline and one below that correspond to highly separated pathways across the retina and into the optic disc, such that they can be treated as essentially independent

^{ 16 }(see dotted line in Fig. 1 separating the superior temporal and inferior temporal sections on the VF). We did consider splitting the temporal sector (Fig. 1) into superior and inferior halves but the degree of separation is not so clear compared with VF nasal fixation.

^{ 20,21 }Accordingly, we did not introduce any distinct bias in to the temporal sector as defined by Garway-Heath.

^{ 16 }Another complication in the retina is the blind spot, which if measured correctly will record values of 0 dB, and, therefore, does not undergo vision loss in the same manner as the rest of the eye.

*y*denote the measurement taken at VF grid location

_{ik}*i*at time

*t*, where the grid cells are numbered by row from top left to bottom right. These are assumed to follow a normal distribution, but to be censored to the left at zero, to account for the idea that a threshold of below zero could be obtained if a higher light intensity, thus, a lower decibel was tested on the Humphrey VF. That is,

_{k}*y*= max( $yik*$,0) where We implement our model within a Bayesian statistical framework. This requires that model parameters are assigned prior distributions that summarize any available prior information. The variance parameter

_{ik}*σ*

_{Ε}

^{2}describes the magnitude of variation due to measurement error. It is assigned an inverse gamma as follows: An informative prior was allocated to

*σ*

_{Ε}

^{2}, as we had preexisting information about the measurement error within the VF. A representative subset of 411 full threshold VFs from 34 eyes were examined. The difference in retested loci within each field was used to estimate variation. A mean variance of 9.4 dB (variance 36 dB) was calculated for VF readings at individual loci (

*y*). Values for the prior were converted to create an informative inverse gamma prior centered around 9.4 dB with a scale of 36 dB.

_{ik}*μ*

_{ik}. The spatial and temporal dependence in the data are generated by the manner in which we model

*μ*. Specifically, we model the spatiotemporal variation in the VF by The overall eye mean is represented by

_{ik}*α*, and the overall temporal trend described by

*β*. CAR prior distributions are employed to account for spatial correlated heterogeneity,

^{ 22 }as described below in Equations 4 to 6. The parameter

*δ*adjusts the mean for spatially correlated heterogeneity at each individual locus, while

_{i}*η*adjusts the trend.

_{i}*δ*parameter the prior is given by where According to our adjacency matrix described earlier,

_{i}*w*= 1 if

_{ij}*i*,

*j*are adjacent and of the same sector,

*w*= 0.3 if

_{ij}*i*,

*j*are adjacent and not of the same sector, and

*w*= 0 if they are nonadjacent. We employ a vague prior for

_{ij}*κ*

_{δ}(

*κ*∼ InverseGamma(1/2,

_{δ}*e*/2)) where

*e*is a small constant, in our case 0.01.

^{19}Hence, the conditional posterior for

*κ*is given by A similar type of CAR prior is applied to the space–time effects,

_{δ}*δ*, where

_{i}*η*and

_{i}*σ*

_{η}^{2}are defined as for

*δ*and

_{i}*σ*. That is A normal prior was employed for

_{δ}^{2}*α*(

*α*∼ Normal (24, 18)) where 24 is the mean threshold score and 18 is two times the SD across all loci of all eyes at time point zero. The doubling of the SD was included to reduce the influence of prior information. Time (

*t*from Equation 3) is centered at zero with a SD of 0.5. We then assign a Cauchy prior to our trend parameter (or global slope),

_{k}*β*. When little is known of the trend an uninformative prior Cauchy (0, 2.5) on a centered dataset with a SD of 0.5 can be used.

^{23}However, as we know more about the scale we assign

*β*∼ Cauchy (0, 5) instead of a scale 22.5 (2.5 times the SD of

*y*).

_{ik}^{ 24 }using Markov Chain Monte Carlo (MCMC) methods. MCMC works by drawing simulations of model parameters from a Markov chain whose stationary distribution matches the required posterior distribution.

^{ 25 }The Metropolis-Hastings (MH) algorithm is used to sample values from the Markov chain. A candidate value (

*θ*) is generated from the proposal distribution. The acceptance probability of

*θ*is calculated, and

*θ*is either accepted or rejected. If

*θ*is accepted it replaces the existing value and, thus, the chain moves, if rejected the existing value remains.

^{ 25 }Gibbs sampling is also used and is a special case of the MH algorithm where all candidate values are accepted. Random walk MH, where candidate values are sampled dependent on the current value of the chain, was used for sampling

*α*and

*β*.

- Set iteration counter
*z*= 0 and initialize*δ*^{(0)},*α*^{(0)},*η*^{(0)},*β*^{(0)}, $ \sigma \delta ( 0 ) $, $ \sigma \eta ( 0 ) $, $ \sigma \Epsilon ( 0 ) $ as specified in the following steps; - Set
*z*=*z*+ 1; - For
*i*in 1-52, sample $ \delta i ( z ) $,conditional on $ \delta 1 : i \u2212 1 ( z ) $, $ \delta i + 1 : 52 ( z \u2212 1 ) $,*α*^{(z−1)},*η*^{(z−1)},*β*^{(z−1)}, $ \sigma \delta ( z \u2212 1 ) $, $ \sigma \eta ( z \u2212 1 ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ using Equation 4; - Sample
*α*^{(z)}conditional on*δ*^{(z)},*η*^{(z−1)},*β*^{(z−1)}, $ \sigma \delta ( z \u2212 1 ) $, $ \sigma \eta ( z \u2212 1 ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ using random walk MH; - For
*i*in 1–52, sample $ \eta i ( z ) $ conditional on $ \eta 1 : i \u2212 1 ( z ) $, $ \eta i + 1 : 52 ( z \u2212 1 ) $,*δ*^{(z)},*α*^{(z)},*β*^{(z−1)}, $ \sigma \delta ( z \u2212 1 ) $, $ \sigma \eta ( z \u2212 1 ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ using Equation 8; - Sample
*β*^{(z)}conditional on*δ*^{(z)},*α*^{(z)},*η*^{(z)}, $ \sigma \delta ( z \u2212 1 ) $, $ \sigma \eta ( z \u2212 1 ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ using random walk MH; - Sample $ \kappa \delta ( z ) $ using Gibbs sampling (Equation 7) conditional on
*δ*^{(z)},*α*^{(z)},*η*^{(z)},*β*^{(z)}, $ \sigma \eta ( z \u2212 1 ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ and, hence, calculate $ \sigma \delta ( z ) $; - Sample $ \kappa \eta ( z ) $ using Gibbs sampling (Equation 7) conditional on
*δ*^{(z)},*α*^{(z)},*η*^{(z)},*β*^{(z)}, $ \sigma \delta ( z ) $, $ \sigma \Epsilon ( z \u2212 1 ) $ and, hence, calculate $ \sigma \eta ( z ) $; - Sample $ \sigma \Epsilon ( z ) $ conditional on
*δ*^{(z)},*α*^{(z)},*η*^{(z)},*β*^{(z)}, $ \sigma \delta ( z ) $, $ \sigma \eta ( z ) $ using random walk MH; and - Repeat steps 2 through 9 until required number of iterations is completed.

*α*and

*δ*

_{ i }s, and

*β*and

*η*

_{ i }s, is handled by replacing

*α*with

*α*+ ∑

*δ*/

_{i}*n*and recentering so that ∑

*δ*= 0 (i.e.,

_{i}*δ*←

_{i}*δ*– ∑

_{i}*δ*/

_{i}*n*).

^{ 26 }Parameters ∑

*η*

_{ i }and

*β*are handled in the same way. Models were run for 100,000 iterations, with burn in period of 12,000, and a thinning factor of 40. Parameter

*α*was initialized at 23, the mean sensitivity of all eyes over all time points. Trend

*β*was initialized at zero. SDs

*σ*

_{ δ }and

*σ*

_{ η }were initialized at 5 and 1, respectively, based on our experience with the variability of the VF measurements. Parameter

*η*is initialized at 0,

*δ*

_{i}at

*y*at time point 0. Convergence was assessed using the Geweke diagnostic.

_{i}^{ 27 }The deviance information criteria (DIC) was calculated and used to test the goodness of fit of our models.

^{ 25 }Like other information criteria it penalizes models, which provide a poor fit and/or are overly complex with superfluous parameters. For each eye, significant progression was defined as a one sided Bayesian

*P*value of less than or equal to 0.05 for the overall eye trend,

*β*. We term our method SPROG, for Spatial PROGression, and, henceforth, refer to our model by this name.

^{ 24 }was used to run point-wise linear regression for each of the 52 nonblind points on the 24-2 field. No outliers were removed. Progression was evaluated at the final time point based on the values of the slope and their significance (

*P*values). Progression was judged using several criteria as defined as in Table 1. Our results may differ slightly from Fitzke et al., Viswanathan et al., and Birch et al.

^{ 28–30 }as PLR was used instead of the PROGRESSOR software.

**Table 1.**

**Table 1.**

Method | Progression Criteria |

P1 | P ≤ 0.001 for ≥1 point^{31} |

P2 | P ≤ 0.001 for ≥2 points^{31} |

P3 | P ≤ 0.001 for ≥3 points^{31} |

P4 | P ≤ 0.001 for ≥4 points^{31} |

V1 | P ≤ 0.01 and slope ≤ −1 dB per year for ≥2 points over 3 consecutive tests^{9} |

V2 | P ≤ 0.01 and slope ≤ −1 dB per year for ≥3 points over 3 consecutive tests^{9} |

F1 | P ≤ 0.1 and slope ≤ −1 dB for inner points, ≤−2 dB for points beyond 15^{°} for ≥2 points^{28} |

F2 | P ≤ 0.1 and slope ≤ −1 dB for inner points, ≤−2 dB for points beyond 15^{°} for ≥3 points^{28} |

VS | P ≤ 0.05 and slope ≤ −1 dB per year for ≥1 nonedge* point, ≤−2 dB for edge points, for ≥1 point^{29} |

B1 | P ≤ 0.05 and slope ≤ −1 dB per year for ≥2 nonedge* points^{30} |

^{ 32 }No specific guidelines or criteria were given to the clinician.

*β*), at an average rate of −0.78 (range, −1.76 to −0.27) dB per year. Figure 2 shows the fit of our model by sector, for an individual eye compared with the observed thresholds. The local slopes (

*β*+

*η*) for an individual can be observed here through the predicted line. We can see that the model reduced the effect of overall variation in the eyes by allowing the predicted points to be influenced by their neighboring points. This also means the result is not overly affected by the occurrence of outliers. However, the model is still flexible enough to allow a certain point within a sector to differ if the mean of a certain point in the eye is consistently lower (Fig. 2, temporal sector). The spatial nature of the VF can be seen in Figure 3, where the local slopes (

_{i}*β*+

*η*

_{ i }) are presented by sector. This eye is progressing at an average rate of −1.48 dB per year. As expected, neighboring points are more similar than points farther apart.

**Figure 2.**

**Figure 2.**

**Figure 3.**

**Figure 3.**

*P*value of 0.1, and including a slope criteria for our model. We chose the value of −0.5 dB per year as generally any loss less than this is considered clinically significant progression. With the changed criteria SPROG still performed best.

**Figure 4.**

**Figure 4.**

**Table 2.**

**Table 2.**

Method | Sensitivity | Specificity |

P4 | 0.00 | 0.74 |

P3 | 0.04 | 0.99 |

P2 | 0.06 | 0.99 |

V2 | 0.12 | 1.00 |

P1 | 0.18 | 0.89 |

V1 | 0.18 | 0.99 |

SPROG (P = 0.05) | 0.64 | 0.66 |

SPROG (P = 0.1, slope < −0.5) | 0.68 | 0.73 |

B1 | 0.74 | 0.48 |

SPROG (P = 0.1) | 0.76 | 0.55 |

F2 | 0.78 | 0.49 |

F1 | 0.86 | 0.36 |

VS | 0.88 | 0.29 |

*P*value, as well as introducing a limit on the slope.

*P*value thresholds from the standard 0.05 or 0.001. We have not included ROC curves for methods with requirements on the numerical value of the slopes because these are not directly comparable. In particular, when the specificity is zero we will obtain a sensitivity of less than one because the slope conditions will not always be attained. The area under the curve (AUC) values, for the same methods, are presented in the caption of Figure 5. SPROG shows the largest AUC, showing that this method has the best overall ROC properties.

**Figure 5.**

**Figure 5.**

^{ 33 }These maps can be compared with the observed data seen in Figure 2. SPROG shows a global progression rate (

*β*) of −1.48 dB per year with local slopes ranging from −1.15 to −1.75. PLR produces a reasonably comparable estimate of the average progression rate (taken as the median of 52 local slopes, so as to control for extreme values) at −2.13, but gives much greater local variation (range, −7.87 to 4.77). One consequence of the erratic behavior of the PLR method is that it indicates very significant vision improvement in some regions of the eye (where the slopes are positive), an implausible finding that is not seen when using the SPROG methodology.

**Figure 6.**

**Figure 6.**

^{ 11–13 }PLR was then used for classification of the eye. These methods can be regarded as an approximation to the statistically more principled SPROG method, since all in affect involve averaging over local neighborhoods. However, our model has the additional advantage that these averages are calculated in conjunction with the other model terms (Equation 3) including measurement error, a temporal trend, and a spatial temporal trend. While previous spatial methods are able to account for spatial correlation and, thus, minimize measurement error, their benefits are still confined by the use of PLR for progression classification.

*y*) over time for each location (

_{ik}*i*). PLR methods require subjective cutoffs, which differ greatly in terms of what significance level is being used, the number of points that must be seen progressing at each test, and for the number of consecutive tests these points must be identified as progressing. This is largely due to the problem of multiple testing, 52 tests for each eye, and leads to either high false positive rates or low detection of true progression. A solution to multiple testing is to use the Bonferroni correction, which maintains that in multiple test scenarios an adjustment to the significance level should be made so that the null hypothesis has the same chance of being rejected, whether testing an individual case or over multiple cases.

^{ 34 }This is done by dividing the significance level (Table 1) by the number of hypotheses being tested (

*n*= 52). Without the correction, when significance is evaluated at the 5% level over 52 points, we would expect two to three points to show significant progression by chance; at 10%, we would expect approximately five points to be significant by chance.

^{ 7 }While this does give an idea of high- and low-ranking methods, it does not solve the problem of a lack of a gold standard.

^{ 7 }Other methods have to juggle between high sensitivity and high specificity by adjusting the various cutoffs. We propose our method reduces both false positives and false negatives by its novel approach to account for spatial correlation, thus, reducing the effects of high variation, which are largely responsible for false predictions. When comparing our model with the clinical reference, we see that our model is the closest of the conservative methods, as well as having the best overall ROC properties. We argue that it is better to slightly overestimate progression than underestimate it, as doing so could lead to a delay in treatment and potential visual loss in a patient.

^{ 7 }this could lead to a long delay before an eye is shown to be progressing.

^{ 4 }Measurement error is also known to be heterogeneous between fields dependent on glaucoma severity. When choosing how to model the magnitude of measurement error (

*σ*

_{Ε}), we looked at including this principle based on our repeated measurements data. We investigated the relationship between variance and distance from fixation, and variance, and mean. We found the power of the models to be low, therefore, in keeping to the principle of parsimony, we chose a constant value for measurement error. However, we acknowledge that this is an approach that may benefit from further research.

^{ 35 }However, the method of modeling spatial correlation is particularly important when considering defects that straddle the sector boundaries defined by Garway-Heath et al.

^{ 16 }While the CAR model allows for correlation between the sectors, we have placed less weight on correlations between sectors than within. The weighting 0.3 for between sectors was chosen as it gave the best combination of sensitivity and specificity over all the eyes. For eyes with defects straddling the sectors it may be of benefit to use a more even weighting scheme. Investigating modeling spatial correlation via a distance based method may prove worthwhile.

*. 2009; 360: 1113–1124. [CrossRef] [PubMed]*

*New Engl J Med**. 2010; 74: 1184–1190. [CrossRef] [PubMed]*

*Neurology**. 2005; 24: 333–354. [CrossRef]*

*Prog Ret Eye Res**. 1987; 105: 1544–1549. [CrossRef] [PubMed]*

*Arch Ophthalmol**. 1999; 40: 648–656. [PubMed]*

*Invest Ophthalmol Vis Sci**. New York: Oxford University Press; 2011.*

*Visual Fields: Examination and Interpretation**. 2012; 90: 10–19. [CrossRef] [PubMed]*

*Acta Ophthalmol**. 2012; 90: 101–108. [CrossRef] [PubMed]*

*Acta Ophthalmol**. 2003; 44: 3873–3879. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 2010; 24: S1–S7. [CrossRef] [PubMed]*

*Eye**. 1995; 79: 207–212. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 2002; 120: 173–180. [CrossRef] [PubMed]*

*Arch Ophthalmol**. 2004; 44: 839–848. [CrossRef] [PubMed]*

*Vis Res**. 2007; 48: 251–257. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 2010; 47: 2904–2910. [CrossRef]*

*Invest Ophthalmol Vis Sci**. 2000; 107: 1809–1815. [CrossRef] [PubMed]*

*Ophthalmology**. 2002; 42: 2654–2659.*

*Invest Ophthalmol Vis Sci**. Boca Raton, FL: Chapman & Hall CRC Press; 2004.*

*Hierarchical Modeling and Analysis for Spatial Data**. 1991; 43: 1–59. [CrossRef]*

*Ann I Stat Math**. 2009; 50: 3249–3256. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 2010; 51: 5657–5666. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. Chichester, West Sussex, England: Chapman & Hall CRC Press; 2003.*

*Disease Mapping with WinBUGS and MLwiN**. 2008; 2: 1360–1383. [CrossRef]*

*Annals of Applied Statistics**. Vienna, Austria: R Foundation for Statistical Computing; 2005.*

*R: A Language and Environment for Statistical Computing**. Boca Raton, FL: Chapman & Hall CRC Press; 2006.*

*Markov Chain Monter Carlo: Stochastic Simulation for Bayesian Inference**. 1995; 82: 733–746.*

*Biometrika**. Oxford, UK: Claredon Press; 1992: 169–193.*

*Bayesian Statistics 4**. 1996; 80: 40–48. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 1997; 81: 1037–1042. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 1995; 102: 1227–1235. [CrossRef] [PubMed]*

*Ophthalmology**. 1996; 37: 1419–1428. [PubMed]*

*Invest Ophthalmol Vis Sci**. 2008; 115: 1557–1565. [CrossRef] [PubMed]*

*Ophthalmology**. 1989; 108: 130–135. [CrossRef] [PubMed]*

*Am J Ophthalmol**. Thousand Oaks, CA: Sage Publications, Inc.; 2007: 104–107.*

*Encyclopedia of Measurement and Statistics**. 2005; 14: 35–59. [CrossRef] [PubMed]*

*Stat Methods Med Res*