**Purpose**:
We developed a statistical model to improve the detection of glaucomatous visual field (VF) progression as defined by the consensus of expert clinicians.

**Methods**:
We developed new methodology in the Bayesian setting to properly model the progression status of a patient (as determined by a group of expert clinicians) as a function of changes in spatially correlated sensitivities at each VF location jointly. We used a spatial probit regression model that jointly incorporates all highly correlated VF changes in a single framework while accounting for structural similarities between neighboring VF regions.

**Results**:
Our method had improved model fit and predictive ability compared to competing models as indicated by the deviance information criterion (198.15 vs. 201.29–213.38), a posterior predictive model selection metric (130.08 vs. 142.08–155.59), area under the receiver operating characteristic curve (0.80 vs. 0.59–0.72; all *P* values < 0.018), and optimal sensitivity (0.92 vs. 0.28–0.82). Simulation study results suggest that estimation (reduction of mean squared errors) and inference (correct coverage of 95% credible intervals) for the model parameters are improved when spatial modeling is incorporated.

**Conclusions**:
We developed a statistical model for the detection of VF progression defined by clinician expert consensus that accounts for spatially correlated changes in visual sensitivity over time, and showed that it outperformed competing models in a number of areas.

**Translational Relevance**:
This model may easily be incorporated into routine clinical practice and be useful for detecting glaucomatous VF progression defined by clinician expert consensus.

^{1–3}Both methods are potentially inefficient and neglect important spatial information contained in the data. We hypothesized that jointly considering changes in sensitivities at each VF location over time when determining progression and identifying regions across the VF where clinicians place increased emphasis when judging progression would preserve important spatial relationships, increasing the explanatory and predictive capabilities of the model. Such an approach increases the amount of available information used in determining progression and more closely resembles the decision-making process used by clinicians who focus on changes across the entire VF and patterns of change as opposed to individual locations separately. The purpose of this study was to develop, implement, and validate a statistical model for detecting glaucoma progression, as defined by the consensus of expert clinicians, and compare its performance to existing approaches.

^{4,5}In the first dataset, each of the VFs were obtained using the SITA standard strategy and VF progression status was determined based on the independent reviews of two expert clinicians. In the case of disagreement, a third clinician was consulted. Reviewers had access to an event-based determination of progression, the Glaucoma Progression Analysis output (GPA; Carl Zeiss Meditec, Inc.), and the print versions of the VF tests over time. Clinicians were instructed to define each series of VF tests as definitely progressing or not.

^{4}Table 1 displays the summary information from this modeling dataset.

^{5}Progression status was once again determined based on review of the VF series and GPA output for each patient, this time by a different set of five expert clinicians. The intraobserver and interobserver agreements were good to excellent and moderate, respectively.

^{5}Table 1 displays the summary information from this validation dataset.

*Y*|

_{i}*,*

**β***∼ Bernoulli{*

**θ***p*(

_{i}*,*

**β***)},*

**θ***i*= 1,…,

*n*and where

*Y*is the progression status of patient/eye

_{i}*i*,

*n*is the total number of unique patient/eye combinations included in the study (see Table 1),

*p*(

_{i}*,*

**β***) is the probability that patient/eye*

**θ***i*is diagnosed as progressing, Φ

^{−1}(·) is the inverse cumulative distribution function of the standard normal distribution,

*d*is the length of time in years that patient/eye

_{i}*i*has been followed in the study,

*= (*

**β***β*

_{0},

*β*

_{1})

*,*

^{T}*β*

_{1}is the parameter describing the association between follow-up time and the probability of being diagnosed as progressing by the clinician expert consensus,

*β*

_{0}is the intercept parameter,

*= {*

**θ***θ*(

*s*

_{1}),…,

*θ*(

*s*)}

_{m}*, and*

^{T}*m*= 52 is the number of VF locations (after removing the two locations within the physiologic blind spot; see Fig. 1). We note that progression often occurs at different rates in different eyes of the same individual and that the clinicians assessed progression for each eye individually. Therefore, we chose to model each eye independently.

*z*(

_{i}*) term is defined as the estimated slope from a simple linear regression analysis of VF sensitivities at VF location*

**s**_{j}*for person/eye*

**s**_{j}*i*. Increasingly negative values of this metric indicate deterioration in vision at a given location while values near zero indicate little change over time. The sum, , represents the total impact of changes in the VF over time at each of the individual VF locations on the probability of patient/eye

*i*being diagnosed as progressing by the clinician expert consensus.

*θ*(

*),*

**s**_{j}*j*= 1,…,

*m*. These spatially varying regression parameters describe the association between a change in VF sensitivities over time at each location and the probability of interest. Including a separate parameter for each VF location allows for increased modeling flexibility and for the possibility that deterioration in the VF at different locations is more (or less) informative with respect to being diagnosed as progressing by the clinician expert consensus. A priori we assume a constant association between sensitivity changes at each location and the probability of interest along with location-specific deviations such that

*θ*(

*) =*

**s**_{j}*θ*

_{0}+

*η*(

*), where*

**s**_{j}*θ*

_{0}represents the constant mean of the spatial process and

*η*(

*) is the slope deviation specific to VF location*

**s**_{j}*. The constant prior mean reflects our initial beliefs that changes in sensitivities over time at any location should have a similar impact on progression diagnosis by the clinician expert consensus, with the possibility that certain locations are more (or less) meaningful to clinicians being represented by the spatial deviations. The*

**s**_{j}*η*(

*) parameters allow for location-specific slopes for each VF location, potentially leading to a model that better describes the process used by clinicians to diagnose glaucoma progression. We allow these parameters to be spatially correlated as detailed in the Prior Specification Subsection so that estimation of their values will depend on the parameter values in surrounding VF locations if necessary. Equation 3 describes how these parameters are related to each other a priori in more detail. Therefore, the probit regression model in Equation 1 can be rewritten as: where represents a measure of the total level of deterioration across the entire VF over time and*

**s**_{j}*θ*

_{0}describes the relationship between this measure and the probability of interests.

*η*(

*). Providing spatial structure in this context helped to control the effect of multicollinearity which often can inflate standard errors and change the sign of parameter estimates. We included a separate predictor at each of the 52 VF locations for a patient/eye with the understanding that sensitivities in similar VF locations are likely positively correlated. Conditional on the spatial proximity matrix,*

**s**_{j}*W*, which describes the spatial relationship between each of the 52 VF locations, the

*η*(

*) parameters were jointly assigned a prior distribution that allowed neighboring parameters to share information spatially, if appropriate. We specified an intrinsic conditional autoregressive (ICAR) prior model such that and where*

**s**_{j}*= {*

**η***η*(

**s**_{1}),…,

*η*(

*)}*

**s**_{m}*,*

^{T}*(−*

**η***) = {*

**s**_{i}*η*(

**s**_{1}),…,

*η*(

**s**_{i}_{−1}),

*η*(

**s**_{i}_{+1}),…,

*η*(

*)}*

**s**_{m}*, is the variance parameter for the spatial process,*

^{T}*w*is the row

_{ij}*i*, column

*j*entry of

*W*describing the proximity between VF locations

*and*

**s**_{i}*, and*

**s**_{j}*w*

_{i}_{+}represents the row

*i*sum of

*W*.

^{6}

*W*is typically specified as

*w*= 1 if

_{ij}*and*

**s**_{i}*share a common border and*

**s**_{j}*w*= 0 otherwise. Locations are not considered to be neighbors of themselves, resulting in

_{ij}*w*= 0 for all

_{ii}*i*. However, the spatial structure across the VF is more complex due to each location's structural relationship with the optic disc. We, therefore, used a similar proximity matrix to that previously introduced.

^{4}In Figure 1, we display the optic disc regions (Fig. 1B) and their spatial association with each location of the VF (Fig. 1A).

^{7}Based on this structural mapping, we defined

*W*(

*δ*) as: with

*w*(

_{ij}*δ*) = 0 for all

*i*. The

*δ*parameter describes the strength of spatial association between neighboring VF locations in adjacent optic disc regions. A number of past studies also have investigated the use of a partially unknown proximity matrix while working in different settings.

^{8–10}

*β*

_{0},

*β*

_{1},

*θ*

_{0}∼ N(0,10

^{10}). The variance parameter for the process was given a vague yet proper uniform prior distribution, such that . The neighborhood proximity matrix parameter was given a uniform prior distribution such that δ ∼ Unif(0,1).

^{11}Summaries of the MC standard errors for presented posterior mean estimates are displayed in Table 2 and Figure 2. In Supplementary Figure S1, we present trace plots for all nonspatial parameters and randomly selected spatial deviations (

*η*(

*)). We noted that these plots are representative of the remaining parameters not displayed. All analyses were done using R Statistical Software and the model fitting details are presented in the Supplementary Material.*

**s**_{i}^{12}

- Model 1: The newly developed probit regression model for spatially correlated predictors with spatially varying VF parameters modeled using the prior distribution and partially unknown spatial proximity matrix.
- Model 2: A probit regression model with separate parameters for VF location that are modeled using exchangeable and independent normally distributed prior distributions.
- Model 3: A probit regression model including only a single spatially aggregated predictor without considering location.
- Models P1 to P4: A probit regression model using the first (Model P1), second (Model P2), third (Model P3), and fourth (Model P4) smallest
*P*values from simple pointwise linear regressions (PLRs), fit separately at each location on the VF, as a single predictor. - Models PS1 to Model PS4: A probit regression model using the first (Model PS1), second (Model PS2), third (Model PS3), and fourth (Model PS4) smallest estimated slopes that are associated with a
*P*value less than or equal to 0.01 from simple PLRs, fit separately at each location in the VF, as a single predictor.

*θ*

_{0}) across all VF locations. We expected this model to outperform Model 1 if there was no spatial variability across the VF in parameter estimates. This would suggest that clinicians weigh each VF location equally when determining progression status.

^{13}Simple PLR analyses are done at each of the available 52 locations on the VF. From each of these separately run analyses, a frequentist

*P*value is obtained indicating if the regression slope is significantly less than zero. In Model P1, we calculated the minimum

*P*value over the entire VF for a patient and entered it as a predictor in a probit regression model with the progression status, as defined by the consensus of expert clinicians, as the outcome. We then used the second (Model P2), third (Model P3), and fourth (Model P4) smallest

*P*values as predictors separately in their own models and compared results. Models PS1 to PS4 are similar except that they use the slopes from the PLR as the predictor and restrict to only those slopes that have a

*P*value less than or equal to 0.01. Patients who do not have a VF location that meets this

*P*value restriction are given a covariate value of zero to indicate this. These models account for statistical significance through the

*P*value restriction as well as the rate of progression through use of the slope values and, therefore. more closely represent progression methods used in practice.

*D̄*, is used to assess model fit while the effective number of parameters,

*p*, is used to describe the model complexity. The DIC then is defined as DIC =

_{D}*D̄*+

*p*, with smaller values being preferred.

_{D}^{14}

*D*is a posterior predictive loss model comparison metric that considers the balance between model fit (

_{k}*G*) and complexity (

_{k}*P*).

^{15}We work with the deviance of the Bernoulli distribution (with continuity correction) and set

*k*= 10

^{10}to calculate

*D*. As with DIC, smaller values of

_{k}*D*suggest an improved model fit.

_{k}*D*is most useful in comparing the predictive performance of a model while DIC is preferred if the model is used mainly for explanatory purposes.

_{k}*P*values. The optimal sensitivity and specificity values were determined using Youden's Index.

^{16}

*δ*, selected prior distribution, selected neighborhood definition, and choice of the prior distribution. These details and results also are presented in the Supplementary Material. The findings suggested that the observed inference is robust to each of these modeling choices. We also presented an exploratory analysis of the robustness of our findings based on the number of available VF tests in Supplementary Table S5. The results appear to be robust to the number of available VF tests overall.

*p̂*

_{i}_{0}(

*,*

**β***θ*

_{0},

*)] using the posterior predictive distribution. Recall that this dataset was not used in building the predictive model. We similarly predicted these same probabilities using the competing models.*

**η***i*= 1,…,

*n*

_{0}, where

*n*

_{0}is number of patient/eyes in the validation dataset (

*n*

_{0}= 100). Values of

*d̂*

_{i}_{0}near zero indicate accurate prediction while large values indicate the opposite. We explored the distribution of these diagnostic measures for each model to determine if any outliers exist. We also calculated and compared and for each model to identify which model provides the best predictions overall. Smaller values for these summary measures indicate a model is predicting well.

*β*

_{1}does not include zero and is positive. This indicates that a patient/eye followed for a longer period of time is more likely to be diagnosed as progressing by the group of clinicians. The 95% CI for

*θ*

_{0}does not include zero and is negative. This indicates that a decrease in any of the VF sensitivities over time (negative slope) leads to an increased probability of being diagnosed as progressing by the clinician expert consensus, regardless of spatial location. In Figure 2A, we display the posterior means for the location-specific slope parameters,

*θ*

_{0}/

*σ*+

_{z}*η*(

*), for Model 1. These parameters represent the increase in an individual's latent variable ( ) with a 1 dB decrease in sensitivity per year at the respective VF location for a patient/eye. A significant increase in significantly increases the probability of progression diagnosis. Therefore, more negative estimates of*

**s**_{j}*θ*

_{0}/

*σ*+

_{z}*η*(

*), as shown in the blue and green regions of Figure 2A, indicate that decreases in VF sensitivities over time at these locations have a greater impact on the clinicians' progression determination. Parameter estimates near zero indicate that decreases in these areas have less impact on the progression determination.*

**s**_{j}*s*

_{19},

*s*

_{20},

*s*

_{21}, and

*s*

_{22}of the inferior-temporal optic disc region are the only locations to have 95% CIs which fail to include zero [

*θ*

_{0}/

*σ*+

_{z}*η*(

**s**_{19}): (−43.89, −1.29);

*θ*

_{0}/

*σ*+

_{z}*η*[

**s**_{20}]: (−37.26, −4.64);

*θ*

_{0}/

*σ*+

_{z}*η*(

**s**_{21}): (−38.25, −4.48);

*θ*

_{0}/

*σ*+

_{z}*η*(

**s**_{22}): (−35.37, −1.19)].

*p*values, indicating poor overall fit to the data. The predictive measures suggest that Model 2 has improved predictive ability when compared to the other models. However, given the DIC results, Model 2 is likely overfitting the data and resulting in improved predictions as a result. In terms of model fit, Model 1 has the smallest DIC when compared to the competing models. Due to the increased flexibility of Model 1, we now have improved predictive performance when compared to the other models (excluding Model 2) as suggested by

_{D}*G*,

_{k}*P*, and

*D*(

_{k}*k*= 10

^{10}). The AUC value for Model 1 is significantly larger than the values from all other models (excluding Model 2) with the

*P*vaues from the paired statistical tests ranging from 0.000 to 0.018. Model 1 also has the highest optimal sensitivity among all models. The competing models are outperformed by Model 1 in terms of model fit and prediction, indicating that accounting for spatial variability across the VF is an important feature of the analysis.

*d̂*

_{i}_{0}, indicating poor prediction for those patient/eyes. Models PS2 to PS4 predict slightly better than Models P1 to P4 overall but are still outperformed by Model 1. These results provide further evidence that the newly introduced model is not suffering from overfitting and suggests that it may be clinically useful.

^{17,18}

^{1,2}Following glaucomatous VF over time using these global indices leads to a reduction of the data to a single measure that ignores the spatial nature of VF loss of glaucomatous progression.

^{19}The pitfall with PLR is that it is based on the erroneous assumption of progression independence between tests locations and only considers test points that show deepening and not enlargement; this assumption is contradicted by available evidence.

^{20}Other studies, such as the Advanced Glaucoma Intervention Study

^{21}and the Collaborative Initial Glaucoma Treatment Study,

^{22,23}have used scoring systems to condense VF data into a single numeric value. The scoring techniques so far have been used only in clinical trials and not in the clinical practice because they display larger long-term variability and fail to provide information with regard to the spatial characteristics of the defects.

^{3}Event-based methods, which define progression based on a prespecified level of deterioration, also have been considered.

^{24}Spatial filters have been offered as a method to reduce variability in VF sensitivity measurements before determining progression status.

^{25}Spatial disease mapping techniques also have been used to model the raw VF sensitivities over time with the goal of identifying significant progression.

^{4}We proposed a new analytical method that uses the sensitivities of all VF test locations and takes into account their spatial correlation. This model outperformed competing models in terms of model fit and predictive ability in the modeling and validation datasets.

*P*value are evaluated and a decision about progression is made based on these 52 values. Our method also attempts to summarize the information contained in these 52 values. However, instead of offering a new deterministic method to make the progression decision, our method uses the introduced statistical model to optimize the decision-making process using all of the information simultaneously. This data-driven method is shown to outperform competing methods with similar forms. We also note that without the use of the spatial prior distribution (Model 2), the method would be essentially ignoring the spatial location of vision deterioration. So while the forms of the model are similar conditional on the spatial random effects, unconditionally Model 1 allows for sharing of information across the VF during estimation. As shown in Figure 2, this results in stable and interpretable spatial patterns that are driven by the observed data as well as our prior understanding of the structure-function correlation between the VF and optic nerve head.

^{26}it is worth noting that a positive attribute of the proposed model is its flexibility as it may be applied easily to other glaucoma VF cluster map schemes. Studies comparing different cluster map schemes may be needed in the future to determine the scheme with the strongest model fit and highest predictive ability to detect glaucoma progression.

*. 1996 ; 80: 40–48.*

*Br J Ophthalmol**. 2008 ; 145: 343–353.*

*Am J Ophthalmol**. 2011 ; 52: 4765–4773.*

*Invest Ophthalmol Vis Sci**. 2013 ; 54: 1544–1553.*

*Invest Ophthalmol Vis Sci**. 2011 ; 118: 60–65.*

*Ophthalmology**. 1974 ; 36: 192–236.*

*J R Stat Soc Series B Stat Methodol**. 2000 ; 107: 1809–1815.*

*Ophthalmology**. 2005 ; 37: 265–285.*

*Geogr Anal**. 2007 ; 14: 433–452.*

*Environ Ecol Stat**. 2011 ; 15: 457–469.*

*Biostatistics**. Suffolk, UK: Chapman & Hall; 1996: 45–57.*

*Markov Chain Monte Carlo in Practice**R:*. Vienna, Austria: R Foundation for Statistical Computing, 2014. Available at: http://www.R-project.org/.

*A Language and Environment for Statistical Computing**. 1996 ; 37: 1419–1428.*

*Invest Ophthalmol Vis Sci**. 2002 ; 64: 583–639.*

*J R Stat Soc Series B Stat Methodol**. 1998 ; 85: 1–11.*

*Biometrika**. 1950 ; 3: 32–35.*

*Cancer**. 2003; 87: 726–730.*

*Br J Ophthalmol**. 2011 ; 21: 573–579.*

*Eur J Ophthalmol**. 1990 ; 31: 512–520.*

*Invest Ophthalmol Vis Sci**. 2004 ; 138: 1029–1036.*

*Am J Ophthalmol**. 1994 ; 10: 144–155.*

*Ophthalmology**. 2003 ; 44: 2613–2620.*

*Invest Ophthalmol Vis Sci**. 1999 ; 106: 653–62.*

*Ophthalmology**. 1997 ; 75: 184–188.*

*Acta Ophthalmol Scand**. 2004 ; 44: 839–848.*

*Vision Res**. 2005 ; 46: 3712–3717.*

*Invest Ophthalmol Vis Sci*