Samuel Berchuck is a Postdoctoral Associate in Duke University’s Department of Statistical Science and Forge-Duke’s Center for Actionable Health Data Science.
Joshua L. Warren is an Assistant Professor of Biostatistics at Yale University.
Analyzing Visual Field Data
In Part I of this series on statistic in glaucoma, we detailed the use of visual fields for understanding functional vision loss in glaucoma patients. Before discussing a new method for modeling visual field data that accounts for the anatomy of the eye, we discussed how visual field data is typically analyzed by introducing a common diagnostic metric, point-wise linear regression (PLR). PLR is a trend-based diagnostic that uses slope p-values from the location specific linear regressions to discriminate progression status. The motivation for PLR is straightforward, assuming that large negative slopes at numerous visual field locations is indicative of progression. This is characteristic of a large class of methods for analyzing visual field data that attempt to discriminate progression based on changes in the DLS across time. This technique is simple, intuitive, and effective; however, it is often limited due to the naivete of modeling assumptions, including the independence of visual field locations.
Ocular Anatomy in the Neighborhood Structure of the Visual Field
To properly account for the spatial dependencies on the visual field, Berchuck et al. 2018 introduce a neighborhood model that incorporates anatomical information through a dissimilarity metric. Details of the method can be found in Berchuck et al. 2018, but we provide a quick introduction. The key development is the specification of the neighborhood structure through a new definition of adjacency weights. Typically in areal data, the adjacency for two locations \(i\) and \(j\) is defined as \(w_{ij} = 1(i \sim j)\), where \(i \sim j\) is the event that locations \(i\) and \(j\) are neighbors. As discussed in Part I, this assumption is not sufficient due to the complex anatomy of the eye. To account for this additional structure, a more general adjacency is introduced that is a function of a dissimilarity metric, \(w_{ij}(\alpha_t) = 1(i \sim j)\exp\{-z_{ij}\alpha_t\}\). Here, \(z_{ij}\) is a dissimilarity metric that represents the absolute difference between the Garway-Heath angles of locations \(i\) and \(j\).
The parameter \(\alpha_t\) dictates the importance of the dissimilarity metric at each visual field exam \(t\). When \(\alpha_t\) becomes large, the model reduces to an independent process, and as \(\alpha_t\) goes to zero, the process becomes the standard spatial model for areal data. Based on the specification of the adjacency weights, \(\alpha_t\) has a useful interpretation with respect to deterioration of visual ability. In particular, \(\alpha_t\) changing over exams indicates that the neighborhood structure on the visual field is changing, which in turn implies damage to the underlying retinal ganglion cell structure. This observation motivates a diagnostic of progression that quantifies variability in \(\alpha_t\) across time. We choose the coefficient of variation (CV) and demonstrate that is a highly significant predictor of progression, and furthermore, independent of trend-based methods such as PLR.
Navigating the womblR Package
To make the method available to clinicians, the R package womblR was developed. The package provides a suite of functions that walk a user through the full process of analyzing a series of visual fields from beginning to end. The user interface was modeled after other impactful R packages for Bayesian spatial analysis, including spBayes and CARBayes. The package name combines Hadley’s naming convention for R packages (i.e., ending a package with the letter R) with the name of the author of the seminal paper on boundary detection, originally referred to areal wombling (Womble 1951).
We will now walk through the process of analyzing visual field data, estimating the \(\alpha_t\) parameters, and assessing progression status. The main function in womblR is the Spatiotemporal Boundary Detection with Dissimilarity Metric model function (STBDwDM). Inference for the method is obtained through Markov chain Monte Carlo (MCMC), which is a computationally intensive method that iterates between updating individual model parameters until enough posterior samples have been collected post-convergence for making accurate posterior inference. Because of the iterative nature of MCMC, the majority of computation is performed within a for loop, so the package is built on C++ through the packages Rcpp and RcppArmadillo. Because of the increased complexity of writing in C++, the pre- and post-processing of the model are done in R with the for loop implemented in C++. The MCMC method employed in womblR is a Metropolis-Hastings within Gibbs algorithm.
Just as a quick aside, with the more recent advent of probabilistic programming, this model could have been implemented using the Hamiltonian Monte Carlo methods used in software like Stan or PyMC3. These programs do not require the derivation of full conditionals, and push the MCMC algorithm to the background. There is undoubtedly a huge market for this type of software, and it is clearly playing a significant role in the popularization of Bayesian modeling. At the same time, implementing MCMC samplers using Rcpp with traditional MCMC algorithms can be instructive, and for those with experience, nearly as quick of a coding experience.
We now begin by formatting the visual field data for analysis. According to the manual, the observed data Y must first be ordered spatially and then temporally. Furthermore, we will remove all locations that correspond to the natural blind spot (which, in the Humphrey Field Analyzer-II, correspond to locations 26 and 35).
###Load package
library(womblR)
###Format data
blind_spot