A Data-Driven Model for Simulating Longitudinal Visual Field Tests in Glaucoma

Abstract

Purpose: To develop a simulation model for glaucomatous longitudinal visual field (VF) tests with controlled progression rates.

Methods: Longitudinal VF tests of 1008 eyes from 755 patients with glaucoma were used to learn the statistical characteristics of VF progression. The learned statistics and known anatomic correlations between VF test points were used to automatically generate progression patterns for baseline fields of patients with glaucoma. VF sequences were constructed by adding spatially correlated noise templates to the generated progression patterns. The two one-sided test (TOST) procedure was used to analyze the equivalence between simulated data and data from patients with glaucoma. VF progression detection rates in the simulated VF data were compared to those in patients with glaucoma using mean deviation (MD), cluster, and pointwise trend analysis.

Results: VF indices (MD, pattern standard deviation), MD linear regression slopes, and progression detection rates for the simulated and patients’ data were practically equivalent (TOST P < 0.01). In patients with glaucoma, the detection rates in 7 years using MD, cluster, and pointwise trend analysis were 24.4%, 26.2%, and 38.4%, respectively. In the simulated data, the mean detection rates (95% confidence interval) for MD, cluster, and pointwise trend analysis were 24.7% (24.1%–25.2%), 24.9% (24.2%–25.5%), and 35.7% (34.9%–36.5%), respectively.

Conclusions: A novel simulation model generates glaucomatous VF sequences that are practically equivalent to longitudinal VFs from patients with glaucoma.

Translational Relevance: Simulated VF sequences with controlled progression rates can support the evaluation and optimization of methods to detect VF progression and can provide guidance for the interpretation of longitudinal VFs.

Introduction

Glaucoma is a group of optic neuropathies characterized by progressive degeneration of retinal ganglion cells and irreversible vision losses.¹ Visual field (VF) testing, or perimetry, is the standard clinical assessment of functional vision integrity in glaucoma. In VF testing, differential light sensitivity (DLS) at discrete locations in the central and peripheral regions of the retina is estimated based on the patient's responses to visual stimuli. Because of the subjective and probabilistic nature of determining psychophysical thresholds, VF measurements can significantly vary from test to test.^2–4 Timely detection of VF deterioration is essential for the management of patients with glaucoma as it impacts decisions regarding treatment and follow-up strategies.⁵ However, due to the high variability of VF test results, accurate detection of VF progression remains challenging.^5,6 

Various methods have been proposed to improve the detection of glaucomatous VF progression.^5–11 However, the lack of consensus on the gold standard to detect VF progression in glaucoma makes it challenging to evaluate the performance of these methods.^6,12 To this end, developing simulation models that generate longitudinal VF data with known progression rates will support the evaluation and optimization of methods to detect VF progression.¹² 

To generate longitudinal VF data for progression analysis, one needs to develop models that simulate both the variability of longitudinal VF tests and the manner that glaucomatous VF defects progress. To describe VF variability, Russell et al.⁴ derived empirical models based on a large-scale longitudinal VF data set. In their model, VF variability for a specific DLS is modeled by the distribution of residuals from pointwise linear regression. This nonparametric method has the benefit of not depending on prior assumptions regarding the statistics of the residuals. However, independent pointwise noise sampled from the empirical distributions was used in their simulation to generate the simulated VF tests.¹³ This makes the implicit assumption that noise at different VF test points is uncorrelated, resulting in the underestimation of the variability of VF global indices (e.g., mean deviation [MD]).^13,14 For this reason, Wu and Medeiros¹⁴ developed noise templates that account for spatially correlated VF noise. Their improved model had better agreement with the statistical characteristics of noise in longitudinal VF data of patients with glaucoma. 

To describe the deterioration of glaucomatous VF defects, the above studies^13,14 did not use explicit progression models. Instead, they fitted models to longitudinal VF sequences from patients to generate estimates of VF loss at different time points. Using such methods, it is difficult to differentiate the effects of noise and actual VF progression on the fitted model. In a separate effort, Gardiner and Crabb¹⁵ and Gardiner et al.¹⁶ developed a model with controlled progression rates to simulate longitudinal VF data, which enables quantitative assessment of methods to detect VF progression. However, this method has some limitations. First, their simulation is based on a limited number of manually defined baseline tests, which may not be sufficient to represent the wide range of VF defects in patients with glaucoma. Second, they assume that VF deterioration occurs in only one anatomically correlated region in the field. This assumption is restrictive and not consistent with clinical observations.^17–19 

In this study, we combined and expanded the methods proposed in previous works^12–16 to develop a data-driven model to simulate glaucomatous longitudinal VF data with controlled progression rates. Our model uses statistical characteristics learned from longitudinal VF data of patients with glaucoma to automatically generate patterns of VF progression in several regions of the retina. To improve the description of spatially correlated noise in longitudinal VF tests, we used a unified pool of noise templates that was generated from VF tests of patients with glaucoma to account for spatially correlated noise in VF tests. The resulting model can automatically generate longitudinal VF sequences with different progression patterns from any baseline VF. Thus, it enables the generation of a large longitudinal VF data set with known progression rates to support the quantitative evaluation and optimization of VF progression detection methods. 

Methods

Longitudinal Visual Field Cohorts

Longitudinal VF data sets from three studies were used to derive the statistical characteristics of glaucomatous VF tests. These data sets include the Rotterdam VF data set,^9,20 the Canadian Glaucoma Study (CGS) VF data set,²¹ and a longitudinal VF data set from the glaucoma clinic of Toronto Western Hospital (TWH). All VF tests in the three data sets were performed with the Humphrey Field Analyzer (Carl Zeiss Meditec, Dublin, CA, USA) with the 24-2 or 30-2 test patterns, using either Full Threshold (Rotterdam and CGS) or SITA Standard (TWH) algorithms. In this study, VFs that were tested with the 30-2 pattern were converted to the 24-2 pattern by keeping the common test locations. Patients with a diagnosis of primary open-angle glaucoma, normal-tension glaucoma, and chronic angle-closure glaucoma who had been followed for at least 5 years (including ≥10 reliable VF tests for each eye) were included in the study. Based on previous work,^11,14,22,23 the criteria for reliable VF tests were set to fixation losses and a false-negative rate <33% and a false-positive rate <20%. This study was approved by the research ethics board of the University Health Network, Toronto, Ontario, Canada, and was conducted in adherence with the Declaration of Helsinki. 

Simulation of Longitudinal VF Tests

In this section, we first introduce a data-driven model that automatically generates VF progression patterns from the baseline test. Subsequently, a noise template model is presented. Finally, the process of generating longitudinal VF data and methods to evaluate the simulation are introduced. Figure 1 shows the overall workflow of the proposed simulation model. 

Figure 1.

 

Overall workflow of the proposed model for simulating longitudinal visual field tests in glaucoma.

Data-Driven Model for VF Progression

Based on the method by Gardiner et al.,¹⁶ VF progression was simulated to occur in clusters of VF points with high anatomic correlation. The anatomic correlation was measured by Garway–Heath angles²⁴ (i.e., the angles of optic nerve fiber bundles underneath VF points entering the optic disc). In the simulation, a VF test point was first specified as the “progression center.” Then, the corresponding progression cluster could be uniquely determined by the progression center and VF points that have Garway–Heath angles lying within ±10° of the progression center. 

The number of progression centers/clusters for each simulated eye was stochastically determined from the empirically derived cumulative distribution of number of eyes containing ≥ n progression clusters in patients with glaucoma. As eyes with baseline scotomas are more likely to show VF progression than those without VF defects, the empirical cumulative distributions of progression clusters for eyes with and without baseline scotomas were computed separately. Similar to other studies,^25,26 a baseline scotoma was defined as three or more contiguous points in the baseline VF test with a pattern deviation probability of <5%. Since the cumulative distribution function is monotonic over the range of 0 to 1, a random number between 0 and 1 was used to automatically select the number of simulated progression centers for each baseline VF test. Therefore, the number of simulated progressions shares the same statistical characteristics as those observed in patients with glaucoma. If the selected number of progression centers equals 0, the simulated eye remains stable. Figure 2 shows the empirical cumulative distribution of progression clusters in VF sequences of patients with glaucoma with (Fig. 2a) and without (Fig. 2b) baseline scotoma. 

Figure 2.

 

Empirical cumulative distributions for the number of progression clusters in visual field sequences of patients with glaucoma with (a) and without (b) baseline scotomas.

For eyes with at least one simulated progression center/cluster, a method was developed to determine the locations of progression clusters. Our data show that (a) VF regions with baseline defects are more susceptible to progression, and (b) the probability of progression is higher in regions with more severe baseline defects. As such, the method was designed to ensure that the simulated progression regions are consistent with the morphologic pattern and severity of baseline defects. Specifically, the average DLSs of clusters centered on each test point in the baseline VF are first calculated. Then, VF test points are sorted in ascending order according to their corresponding cluster average DLS. Thus, the first progression center is selected from the top of the sorted list, and the corresponding cluster is used as a progression cluster in the simulation. (Clusters with a low mean DLS have limited capacity to show progression due to the flooring effect of VF measurements. Therefore, clusters with average DLS that are less than 10 decibels are not used as progression clusters.) All points in the corresponding cluster are removed from subsequent selections. This process is repeated until all required progression centers/clusters are found. In the simulation, VF points within progression clusters are assigned linear rates of deterioration. All other points in the VF are considered stable with a linear deterioration rate of −0.1 decibels/year (dB/y), representing the age-related progression.^15,27 The progression rate is a hyperparameter in our proposed model and can be set to any desired value. To facilitate comparisons between simulated VF data and data from patients with glaucoma, instead of assigning a uniform progression rate to all progression clusters, the progression rate for each cluster was randomly selected from the probability distribution of the mean progression rates of clusters at the same location in patients with glaucoma. 

Noise Template Model for VF Variability

We used a nonparametric method introduced by Wu and Medeiros¹⁴ to collectively describe the short- and long-term VF variability.²⁸ First, pointwise linear regression models were fitted to VF sequences from stable eyes in our data sets. The fitted values were treated as the “true sensitivity,” and the residuals were taken as noise or variability in the field. Residual values associated with the same true sensitivity were then pooled across locations and patients to construct empirical distributions of VF variability for different sensitivity values. Next, to construct noise templates, residuals at each location (in dB) were mapped to probability scores (0–1) based on their percentile ranks in the empirical distribution of VF variability for the true sensitivity. The probability scores of all locations in a VF collectively constitute the noise template for this field, representing the normalized magnitude and spatial patterns of noise. By allowing for spatial correlations between noise at different test points, noise templates can describe global visit effects²⁹ (e.g., fatigue) that affect the DLS measurements at several VF points during the test. Unlike the study by Wu and Medeiros,¹⁴ where patient-specific noise templates were used (i.e., creating different noise template pools for each patient), we pooled noise templates across all patients and used them interchangeably. Pooling the noise templates from different patients may offer several advantages for the simulation, which are explained in the Discussion section. As glaucomatous VF measurements in superior and inferior hemifields are commonly assumed independent,³⁰ we vertically flipped each obtained noise template around the horizontal meridian to augment the number of templates in the pool. Furthermore, we chose to use stable VF sequences to derive the residuals. As the primary source of variation in stable VF tests is measurement noise, the impact of nonlinear VF progression in estimating variability can be minimized by using stable eyes. In this study, a stable VF sequence was defined as a sequence with MD linear regression slope falling in the range of ±0.5 dB/y. 

Simulating Longitudinal VF Data

The average of the first two tests from the VF sequence of patients with glaucoma was used as baseline VFs for the simulation. VF tests at different follow-up time points were generated using the cluster progression model presented earlier. The generated VF tests represent the expected DLS values at each follow-up visit, referred to as the expected field. To simulate the noise in each follow-up VF test, a noise template was randomly sampled from the noise template pool. Then, a noise field was reconstructed based on the probability scores in the selected noise template and the distribution of residual errors for the DLS in the expected field. Last, the noise field was added to the expected field to construct the measured VF at each follow-up visit. 

Evaluating the Simulated Longitudinal VF Data

We evaluated the model by comparing the MDs, pattern standard deviations (PSDs), and MD linear regression slopes of simulated longitudinal VFs to those of patients with glaucoma. We also compared the standard deviation of MD linear regression residuals between the simulated data and the data from patients with glaucoma to demonstrate that using noise templates leads to more accurate modeling of VF variability. 

To evaluate the modeling of VF progression, we compared the percentage of eyes showing statistically significant MD deterioration (i.e., detection rates) in the simulated VF data at different time points to those in VF data of patients with glaucoma. Following the convention in previous studies,^10,11 statistically significant MD deterioration is defined as a negative MD linear regression slope with P < 0.05. With the same criteria, sensitivity analysis was used to compare the detection rates in three sets of simulated longitudinal VF data. These three simulated VF data sets were generated using (A) the proposed model (i.e., multiple progression centers with noise templates), (B) multiple progression centers with independent pointwise noise, and (C) a single progression cluster with noise templates. Furthermore, we compared the detection rates in simulated and patients’ VF data using two other trend-based methods: cluster and pointwise trend analysis. In these comparisons, we used the detection criteria proposed by Gardiner et al.¹¹ (i.e., for cluster trend analysis: ≥3 clusters with negative mean cluster slope with P < 0.117; for pointwise trend analysis: ≥9 points with negative slopes with P < 0.138). These criteria were derived with a test–retest data set so that each detection method could achieve 95% specificity in detecting VF progression.¹¹ 

In our experiments, we used 1008 baseline VFs and generated 10 different longitudinal sequences for each baseline field. Each simulated VF sequence consists of 15 VF tests (including the baseline) with 6-month intervals, representing a 7-year follow-up period. 

Statistical Analysis

The two one-sided test (TOST) procedure³¹ was used to analyze the equivalence between simulated VF data and data from patients with glaucoma. In the TOST, two composite null hypotheses (H₀₁ and H₀₂) are tested: H₀₁: µ₁ − µ₂ ≥ Δ_U and H₀₂:  µ₁ − µ₂ ≤ Δ_L, where µ₁ and µ₂ are the mean values of parameters that describe the simulated and patients’ data, respectively, and Δ_U and Δ_L are the upper and lower bounds reflecting the clinical tolerance for equivalence. If both null hypotheses can be rejected (P < 0.05), µ₁ and µ₂ are considered practically equivalent. In this study, the equivalence bounds were set to ±0.5 dB for MD and PSD, ±0.1 dB/y for MD slope, and ±4% for detection rate. The selection of equivalence bounds for MD and PSD (±0.5 dB) is based on the clinically acceptable levels for prediction errors.²² For MD slope, the equivalent bounds were set to the expected progression rate due to aging (−0.1 dB/y). For detection rate, the bounds were determined based on the change in detection rate when using the upper and lower bounds for MD slope as the detection criteria (with MD trend analysis). Moreover, we adopted one-sided Mann–Whitney U tests instead of one-sided t-tests to deal with the nonnormality of our data. Levene's test for equality of variance was used to analyze differences between the variance of MD linear regression residuals in simulated VF data and data from patients with glaucoma. 

Model Availability

Our simulation model is publicly available at https://github.com/lcapacitor/glaucomatous-longitudinal-vf-simulator. 

Results

Data Characteristics

A total of 20,096 VF tests from 1008 eyes of 755 patients with glaucoma were included in this study. The median (interquartile range [IQR]) MD and age in the baseline tests were −3.35 (−1.16 to −7.78) dB and 59.7 (50.9 to 67.5) years, respectively. (The IQR for MD is reported in the format of “Q3 to Q1.”) The median (IQR) follow-up length was 11.1 (9.6 to 16.1) years, with a median number of tests per eye of 19 and a median test interval of 0.5 years. Table 1 summarizes the data characteristics of our longitudinal VF data sets. 

Table 1.

 

Data Characteristics of the Longitudinal VF Data Sets

 

To derive the distributions of pointwise linear regression residuals and the noise templates, 17,273 VF tests from 864 stable eyes of patients with glaucoma were used. A pool of 34,546 noise templates was obtained after vertically flipping the original 17,273 templates around the horizontal meridian. 

Comparison of Simulated and Patients’ Longitudinal VF Data

The mean and standard deviation (indicated in parentheses) for MD, PSD, and MD linear regression slope for the VF data from patients with glaucoma were −5.0 (5.6) dB, 5.1 (3.7) dB, and −0.19 (0.41) dB/y. For the simulated VF data, the mean (standard deviation) MD, PSD, and MD linear regression slope were −5.1 (5.5) dB, 5.6 (3.3) dB, and −0.19 (0.32) dB/y, which were all practically equivalent to those from patients with glaucoma (TOST P < 0.001). Furthermore, the same comparisons were carried out in different severity groups: the mild (initial MD >−6 dB), moderate (initial MD ∈ [ − 12,   − 6] dB), and severe (initial MD < −12 dB) groups. The simulation results demonstrated good agreement with those from patients with glaucoma in the different groups (TOST P < 0.01). A detailed comparison can be found in Table 2. 

Table 2.

 

Comparisons of VF Indices Between Simulated and Real Data in Different Severity Groups

 

Figure 3 shows histograms of MD linear regression residuals for patients with glaucoma and two sets of simulated VF data. Similar to the results from Wu and Medeiros,¹⁴ there was no significant difference between the standard deviation of MD linear regression residuals for VF data of patients with glaucoma and simulated data when noise templates were used to model VF variability (Figs. 3a, 3b). By contrast, the standard deviation of MD residuals when using the independent noise model for each VF test point was significantly lower (Levene's test P < 0.0001, Fig. 3c). 

Figure 3.

 

Histograms of linear regression residuals for the mean deviation in (a) longitudinal VF data of patients with glaucoma, (b) simulated longitudinal VF data with spatially correlated noise templates, and (c) simulated longitudinal VF data with independent noise model to represent pointwise VF variability. The variance in (a) and (b) is similar, whereas the variance of (c) is smaller.

VF Progression Analysis Using Simulated Longitudinal VF Data

Figure 4 shows the progression detection rates over a period of 2 to 7 years from baseline when MD trend analysis was applied to VF sequences of patients with glaucoma and three sets of simulated VF data (see “Evaluating the Simulated Longitudinal VF Data”). After 3, 5, and 7 years from the baseline, the progression detection rates using MD trend analysis in VF data of patients with glaucoma were 8.5%, 16.9%, and 24.4%, respectively. The detection rates for the simulated data generated by our model (set A) were practically equivalent to those in patients with glaucoma (TOST P < 0.0001). Specifically, the mean (95% confidence interval [CI]) detections rates after 3, 5, and 7 years from the baseline were 7.7% (7.3% − 8.0%), 17.1% (16.4% − 17.8%), and 24.7% (24.1% − 25.2%), respectively. In comparison, the mean (95% CI) detection rates at the same time points for the simulated data with independent pointwise noise (set B) were 20.5% (19.7% − 21.3%), 26.9% (26.1% − 27.8%), and 40.7% (39.7% − 41.7%). When assuming the simulated VF progression occurs in only one cluster of points (set C), the mean (95% CI) detection rates were 2.2% (2.0% − 2.5%), 6.3% (6.0% − 6.7%), and 11.7% (11.3% − 12.2%), respectively. 

Figure 4.

 

Percentage of eyes showing VF progression (i.e., detection rates) over a period of 2 to 7 years from the baseline test. The mean deviation trend analysis was used to detect VF progression. The solid red curve indicates the detection rates in longitudinal VF data of patients with glaucoma. The green dashed curve shows the mean detection rates in simulated data (set A) generated by our model (using multiple progression clusters per eye and spatially correlated pointwise noise, denoted as Multiclusters + Correlated Noise). The orange dash-dotted curve represents the mean detection rates in simulated data (set B) generated by a model using multiple progression clusters and independent pointwise noise (denoted as Multiclusters + Independent Noise). The purple dotted curve indicates the mean detection rates in simulated data (set C) generated by a model using a single progression cluster per eye and spatially correlated pointwise noise (denoted as Single Cluster + Correlated Noise). The shaded areas of the curves represent the 95% CI of the detection rates.

Similarly, the detection rates when using cluster and pointwise trend analysis in the simulated VF data were practically equivalent to those of VF data from patients with glaucoma (TOST P < 0.01). Specifically, when using cluster trend analysis for data from patients with glaucoma, the detected rates were 9.3%, 18.4%, and 26.2% after 3, 5, and 7 years from the baseline, respectively. In comparison, the mean (95% CI) detection rates using cluster trend analysis with simulated VF data generated by our model were 10.2% (9.5% − 11.0%), 18.7% (18.0% − 19.5%), and 24.9% (24.2% − 25.5%) at the same time points. For pointwise trend analysis, the detected rates in patients with glaucoma were 17.3%, 28.1%, and 38.4% after 3, 5, and 7 years, while the detection rates with the simulated VF data were 18.9% (18.3% − 19.6%), 29.4% (28.6% − 30.2%), and 35.7% (34.9% − 36.5%) at the same time points, respectively. 

Example: Prediction of Follow-Up VF Tests

In Figure 5a, we present VF tests of a patient with glaucoma who was followed over a 7-year period. The baseline VF has an arcuate-shaped scotoma in the superior field, and the VF sequence manifests a moderate progression with MD linear regression slope of −0.5 dB/y. Figures 5b–e demonstrate four examples of simulated VF sequences with different progression rates using the patient's first VF test as the baseline. The examples show that when the simulated VF progression rate matches the patient's VF progression rate (−0.5 dB/y in MD slope), the model generates (predicts) future VF tests (years 1 to 7) that are similar to those measured in the patient's follow-up tests (compare Figs. 5a, 5c). When the simulated VF progression rates differ from the patient's VF progression rate, the differences between predicted and measured VF tests at follow-up visits are larger (compare Figs. 5a to Figs. 5b−e). 

Figure 5.

 

(a) VF of a patient with glaucoma over a 7-year follow-up period. (b–e) Simulated VF sequences for the same baseline VF test (the first test in the patient's VF sequence) with different progression rates. The patient data show moderate VF progression with MD linear regression slope of −0.5 dB/y. (b) A simulated stable VF sequence (MD slope of 0 dB/y). (c) A simulated VF sequence with the same progression rate (MD slope of −0.5 dB/y) as the patient's data (notice the similarity between the progression patterns in panels a and c). (d, e) Simulated VF sequences with higher progression rates than that of the patient's data (MD slopes of −1.0 dB/y and −1.5 dB/y, respectively).

The simulated VF sequences with various progression rates can be viewed as predictors of future VF tests for the patient (e.g., Figs. 5b−e). Thus, one can estimate the VF progression rates of patients by comparing differences (e.g., mean square error [MSE]) between patients’ VF tests at follow-up visits and simulated/predicted VFs. For each patient, the progression rate of the simulated VF sequence with the lowest difference is used as an estimate of the patient's VF progression rate. To demonstrate this approach, VF progression rate of the patient in Figure 5a was estimated. We first generated 100 longitudinal VF sequences for each of the four progression rates shown in Figures 5b−e. These sequences represent stable VFs (0 dB/y in MD slope), VFs with moderate progression (−0.5 dB/y in MD slope), VFs with fast progression (−1.0 dB/y in MD slope), and VFs with very fast progression (−1.5 dB/y in MD).³² For each progression rate, the average MSEs between the patient's VF tests at follow-up visits and the 100 simulated VF sequences were calculated. To construct the 95% CI for the MSE, the above process was repeated 10 times. Figure 6 shows the cumulative sum of average MSEs between the simulated VFs for the four different progression rates and the patient's VFs over 7 years. After the first year, the MSEs between the simulated and patient's VF data for all four progression rates are not significantly different. After 2 years, the cumulative MSEs between simulated VF sequences with progression rates ≤−0.5 dB/y and the patient's data are significantly lower (t-test P < 0.01) than the cumulative MSEs with simulated stable VFs (0 dB/y in MD slope). In other words, the data in Figure 6 suggest that after 2 years, the patient's VFs are deteriorating at a rate equal to or worse than −0.5 dB/y. After 3 years, the cumulative MSEs between the patient's data and simulated VFs with moderate and fast progression rates (−0.5 dB/y and −1.0 dB/y in MD slope) are significantly lower (t-test P < 0.01) than the MSEs with simulated VFs that are either stable (0 dB/y in MD slope) or progressing at a very fast rate (−1.5 dB/y in MD slope). After 4 years, Figure 6 shows that the cumulative MSEs between the simulated, moderately deteriorating VFs (i.e., −0.5 dB/y in MD slope) and the patient's data are significantly lower (t-test P < 0.01) than the cumulative MSEs with all the other three simulated sequences. Note that after 4 years, the estimated VF progression rate (−0.5 dB/y) is the same as the clinical determination (moderate progression) obtained by measuring the MD slope over the 7 years. Moreover, even before the data in Figure 6 converged to a single solution in year 4, one could get insights into the range of progression rates that are consistent with the measured VFs at each time point. 

Figure 6.

 

Graphs of the cumulative sums of the average MSE between a VF sequence of a patient with glaucoma (Fig. 5a) and four sets of 100 simulated VF sequences with different progressing rates. The progressing rates represent stable VF sequences (0 dB/y in MD slope), VF sequences with moderate progression (−0.5 dB/y in MD slope), VF sequences with fast progression (−1.0 dB/y in MD slope), and VF sequences with very fast progression (−1.5 dB/y in MD slope). The error bars represent the 95% CI for the MSE at each time point.

Discussion

In this study, we presented a novel data-driven simulation model for glaucomatous longitudinal VF tests. The model uses baseline VFs of patients with glaucoma and statistical features derived from longitudinal VF data to automatically generate VF sequences with controlled progression rates. Using multiple progression clusters to model VF progression and a unified pool of spatially correlated noise templates to model test variability resulted in simulated VF sequences that are practically equivalent to those from patients with glaucoma. Models that use only one progressing cluster to describe VF progression or models that use uncorrelated VF noise to describe VF test variability did not achieve such equivalence (Fig. 4). 

To evaluate the modeling of VF progression, we compared the detection rates in simulated VF data at different follow-up time points to those measured in data from patients with glaucoma. In these comparisons, we used three trend-based detection methods for VF progression, which analyze deteriorations occurring in the whole VF, the cluster of points, and a single VF point, respectively. As shown in the “VF Progression Analysis Using Simulated Longitudinal VF Data” section, the detection rates in simulated and patient data at different time points for all three trend-based analysis methods were practically equivalent. The ability to achieve practically equivalent detection performance for global (whole field), regional (clusters), and local detectors suggests that the simulated VF sequences can faithfully reproduce the global, regional, and local changes in longitudinal VFs of patients with glaucoma. In addition, when using different MD slope thresholds to define VF progression (i.e., progression is defined as statistically significant deterioration with MD <−0.25 dB/y, <−0.5 dB/y, and −1 dB/y), the detection rates obtained from our simulation results were also practically equivalent to those obtained from patients with glaucoma (TOST P < 0.01). 

To represent longitudinal VF test variability, our model used pooled noise templates from all patients rather than using patient-specific templates as suggested in previous studies.¹⁴ The rationale for this change is that (a) our method does not require a patient's entire VF history for simulations. Thus, it may be able to simulate the projected trajectory of any patient from a single baseline test. (b) It may be able to better capture the statistics of the global visit effects by pooling across all eyes. (c) Patient-specific templates require shuffling within the same eccentricity for data augmentation, which may not be entirely consistent with anatomic correlations between VF points.²⁴ Nevertheless, we compared the results of using unified versus patient-specific noise templates. Both methods demonstrated similar capabilities in characterizing the MD of patients with glaucoma. Specifically, the MD variabilities (measured by the standard deviation of MD linear regression residuals¹⁴) were 1.4 dB, 1.3 dB, and 1.3 dB for simulated data using patient-specific templates, simulated data using our unified pool, and data from patients with glaucoma, respectively. 

Because of the large number of patients in this study, VF indices estimated from either the simulated or patient data have very small standard errors of the mean (SEMs) compared with accepted errors in clinical practice. For example, SEMs for MD in patient and simulated data are 0.04 dB, while the corresponding acceptable errors in clinical practice are 1 dB (see “Statistical Analysis” section). The same observation also can be found for all VF indices. Due to the small SEMs in this study, the hypothesis that VF indices for the simulated and patient data are the same is rejected (at P < 0.05), even when the differences are not clinically significant. To overcome this problem, we used the equivalence test (i.e., the TOST procedure) to show that differences between data generated by our model and data from patients with glaucoma lie between clinically acceptable bounds for equivalence. The result of the TOST procedure demonstrated that all VF indices for the simulated and patient data are practically equivalent (P < 0.01). 

As shown in Figure 5, the novel simulation model can generate longitudinal VF tests with different progression rates for the same baseline test. This means that a large longitudinal VF data set can be generated from a relatively limited number of baseline tests. Obtaining VF sequences with fast progression rates from clinical records of patients with glaucoma is challenging since most longitudinal VF data from them are stable. For instance, only 13% of the VF sequences in our data sets have MD linear regression slope worse than −0.5 dB/y, and VF sequences with fast progressing rates are often interrupted by clinical interventions such as surgery or change in medication. However, data sets that include such sequences are crucial when developing data-demanding methods (e.g., machine learning algorithms) to detect VF progression. Therefore, our simulation model can be used to generate VF sequences to augment data from real patients. 

In the “Example: Prediction of Follow-up VF Tests” section, we described a method to estimate VF progression rates by comparing differences between the simulated and patient's VF tests at follow-up visits. Using this method, we estimated the progression rate of the VF sequence in Figure 5a. The result (Fig. 6) demonstrates that the range of the estimated VF progression rate is tightened as the number of follow-up visits increases and eventually converges to a value that is consistent with the clinical measurement of the progression rate. Specifically, after 2 years from baseline, the estimated progression rate is at least −0.5 dB/y (i.e., the estimated progression rate is in the range of −0.5 dB/y to −1.5 dB/y at this point). After 3 years, the estimated VF progression rate is narrowed down to the range of −0.5 dB/y to −1.0 dB/y. After 4 years and thereafter, the estimated VF progression rate converges to −0.5 dB/y, which is the same as the clinical determination (moderate progression) obtained by measuring the MD slope over the 7-year period. Therefore, as early as 2 years after the baseline VF test, the data in Figure 6 could be used to provide clinicians with insights regarding the range of progression rates that are consistent with the measured VFs from the patient. It should be noted that this example only illustrates the potential translational use of our simulation model. Quantitative and systematic evaluations are required to validate the performance and will be carried out in the follow-up publication. 

Our model has several limitations. We used a linear progression model to simulate DLS thresholds at different time points. Various studies suggested that VF progression can be better characterized by nonlinear models.^33–36 However, our model can be easily extended to incorporate nonlinear models of VF progression. Moreover, to achieve a large sample size, we mixed VF data sets tested using SITA Standard and Full Threshold algorithms, which may result in different noise characteristics due to postprocessing in the SITA Standard. Finally, our model assumes that the disease progresses without interruptions and does not explicitly account for clinical interventions, which can slow down the progression. 

In conclusion, this study presents a data-driven simulation model to automatically generate glaucomatous VF sequences with controlled progression rates. The model is based on statistical features learned from large-scale longitudinal VF data of patients with glaucoma. We have shown that VF indices and VF progression detection performance in the simulated and patients’ data are practically equivalent. Simulated VF sequences with controlled progression rates can support the evaluation and optimization of methods to detect VF progression and provide guidance for interpreting measured changes in longitudinal VFs of patients with glaucoma. 

Acknowledgments

Supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant 458039 and Vision Science Research Program, University Health Network, Toronto, Ontario, Canada (to YL and RBS). 

Disclosure: Y. Li, None; M. Eizenman, None; R.B. Shi, None; Y.M. Buys, None; G.E. Trope, None; W. Wong, None