Spatial and Temporal Characteristics of Visual Field Progression in Glaucoma Assessed by Parallel Factor Analysis.

PURPOSE: To explore spatial and temporal characteristics of glaucomatous visual field (VF) progression through multi-way decomposition of data.
METHODS: Six serial VF exams with intervals of 6.0 ± 1.0 months in 121 pre-perimetric glaucoma eyes and 80 perimetric glaucoma eyes were arranged into a three-dimensional cube. The data were decomposed using parallel factor analysis.
RESULTS: Three tri-linear components (i.e., spatial scores, temporal loadings, and subject-specific loadings) were identified. Component 1 clearly showed differences between superior and inferior hemispheres, linear trends over time, and wide variability in perimetric glaucoma. Findings were compatible with well-known characteristics of glaucomatous VF defects. Component 2 showed nasal and central areas in contrast with superior, inferior, and temporal peripheral locations, whereas component 3 showed a contrast between nasal and temporal hemispheres. Both components 2 and 3 failed to show clear temporal trends.
CONCLUSIONS: Identification of spatio-temporal patterns shows new possibilities for a multi-way decomposition methodology for earlier diagnosis and prediction of glaucomatous VF progression.

Glaucoma is a disease characterized by loss of retinal ganglion cells and the retinal nerve fiber layer, which eventually leads to deterioration of visual function [1]. Functional deterioration is currently measured by standard automated perimetry. An accurate and reliable method is needed to measure deterioration of the visual field (VF) to monitor the progression of glaucoma and evaluate the efficacy of therapies. Predicting future VF changes is currently one of the most challenging areas in managing glaucoma. To predict VF accurately, identifying the temporal and spatial characteristics of VF changes is essential. In structural aspects, VF data can be interpreted as multi-dimensional data, containing multiple VF locations for each test performed in every single patient. Past research on VF has focused on contemporary cross-sectional analysis [1234] on regression methods applied separately in different locations (pointwise approach) [56789], or on regression analysis based on global indices [101112].

Cross-sectional approaches cannot capture any temporal variability, (i.e., progression of VF defects) because the data is a collection of single tests on each subject. VF locations reflect specific locations in individual human retinas, and each location should be correlated with the others (possible correlations are 52 × 52 = 2,074 correlations) because of the anatomical patterns of retinal nerve fibers and ganglion cells [131415]. Because models based on point-wise approaches separately analyze discrete VF locations, these models have inherent difficulty incorporating such correlations among the locations. Some recently developed regression models capture spatial characteristics by analyzing clusters rather than individual points, either following point-wise regression [13] or prior to regression [1617]. Also notable is the application of a Monte-Carlo simulation in analysis of spatial correlations in retinal sensitivity changes [18]. Using pattern deviation plot values from a baseline to two follow-up tests, probabilities were calculated for negative changes in one of 52 VF locations (in a 24-2 pattern, excluding two blind-spot locations). The results were compared with Monte-Carlo simulation results, providing solid evidence of spatial correlations in VF defect patterns in a glaucoma population. Location-wise aggregate results cannot reflect any subject-specific characteristics, however, and the large number of parameters (2,704 probabilities for each given level of defect) makes it difficult to integrate this methodology into a practical diagnostic framework.

Hence, this study introduces a multi-dimensional approach implementing a multi-modal decomposition method known as parallel factor analysis/canonical polyadic decomposition (PARAFAC/CANDECOMP) to incorporate intra-subject, inter-subject, and temporal variability into a single framework for simultaneous characterization of temporal and spatial features of glaucomatous VF. Using this approach, we analyze the characteristics of longitudinal VF data in our cohort of glaucoma patients.

Materials and Methods

Subjects

Medical records of glaucoma patients who were followed with a Humphrey VF analyzer (Carl Zeiss Meditec, Dublin, CA, USA) at the glaucoma clinic in the department of ophthalmology at Asan Medical Center between January 2007 and August 2013 were retrospectively reviewed. Eyes with more than six serial regular reliable VF exams (false-positive errors <15%, false-negative errors <15%, and fixation loss <20%) at intervals of 6.0 ± 1.0 months (150 to 210 days) were enrolled. Initial exams were eliminated from analysis to minimize learning effects. All participants had best-corrected visual acuity of 20 / 40 or better, with a spherical refractive error between −6.0 and +4.0 diopters (D) and a cylinder correction within +3 D. All participants had a normal anterior chamber and an open-angle on slit-lamp and gonioscopic examinations. A glaucomatous VF defect was defined as a cluster of three points with a probability <5% on the pattern deviation map in at least one hemifield, including at least one point with a probability of <1%, a cluster of two points with a probability of <1% and a glaucoma hemifield test result outside normal limits, or a pattern standard deviation outside 95% of normal limits at baseline. Participants who had both glaucomatous optic disc changes (cup to disc ratio ≥0.7, diffuse or focal neural rim thinning, disc hemorrhage, or retinal nerve fiber layer defects) and a glaucomatous VF defect were defined as the perimetric glaucoma group (PG). Patients who had glaucomatous optic disc changes but a normal VF at baseline were assigned to the pre-perimetric glaucoma group (PPG). Exclusion criteria included the presence of any neuro-ophthalmic disorders, retinal diseases, a history of intraocular surgery other than uncomplicated cataract and glaucoma surgery at baseline, and/or a history of any intraocular surgery during the follow-up period. All included patients had a baseline VF mean deviation no worse than −9.0 decibels to control heterogeneity among subjects.

This study adhered to the Declaration of Helsinki principles and the study protocol was approved by the institutional review board of the Asan Medical Center, University of Ulsan, Seoul, Korea (2017-1311). Informed consent was waived due to the retrospective nature of the study.

Statistical analysis

The PARAFAC/CANDECOMP method was introduced in 1970 [1920] in psychometrics and has been widely adopted in chemometrics as well [21]. Recently, the method has gained some interest in space-time signal brain analysis, electroencephalography analysis [2223] and functional neuroimaging [2425]. The PARAFAC/CANDECOMP can be seen as a three-mode extension of principal component analysis and aims to reduce three-mode data table X by extracting the same number of components (F) for every mode [26]. Data table X is represented as the sum of each tri-linear component and residual E, usually presented as follows:

where a, b, c are elements of the loading matrices A, B, and C, respectively. Term e represents elements of residual matrix E.

The equation can be additionally presented in matrix form after dropping residuals (e) as:

where a, b, c are the f columns of the loading matrices A, B, and C, respectively.

An important advantage of PARAFAC/CANDECOMP in comparison to other techniques such as bilinear principal component analysis or independent component analysis is that the solution is unique without additional orthogonality or independence constraints as long as the right number of components is used and the signal-to-noise ratio is appropriate [212327]. In previous research, PARAFAC/CANDECOMP analysis was calculated with R ver. 3.0.2, function CP of R statistical package ThreeWay ver. 1.1.1 (http://cran.r-project.org/web/packages/ThreeWay/) [28]. Other statistics were analyzed with SPSS ver. 15.0 (SPSS Inc., Chicago, IL, USA).

Results

One hundred and twenty-one PPG and 80 PG eyes were included in the final analysis. The 121 PPG eyes were comprised of 63 male and 58 female eyes. In PG eyes, 43 and 37 were male and female eyes, respectively. Clinical characteristics of the participants assessed at baseline and at final follow-up examinations are described in Table 1. In PG subjects, measures of VF mean deviation were worse at both baseline and final visits, as expected (−0.32 ± 1.30 vs. −4.80 ± 3.24, p < 0.001). Visual acuity, intraocular pressure, and number of glaucoma medications showed no significant differences between the two groups.

Table 1

Patient characteristics

PPG = pre-perimetric glaucoma; PG = perimetric glaucoma; logMAR = logarithm of minimum angle of resolution; IOP = intraocular pressure.

*Mann-Whitney U test.

Fitting the model

A graphical representation of our model is shown in Fig. 1. The first mode represents specific locations on the VF map (spatial mode), and is labeled as “VF locations.” The second mode represents subjects (subject mode), and the third mode represents time (time mode). After centering raw data across the first and second modes, we ran PARAFAC/CANDECOMP analysis with one to six components (components 1 to 6) to determine the optimal number of components. The model's goodness-of-fit improved as the number of components increased, but a rapid decline in core consistency was shown with each additional components, especially after the fourth component (two components 100%, three components 90.16%, and four components 33.87%) [28]. Thus, a three-component-model with a fit of 51.76% (Fig. 2) was chosen for further analysis.

Fig. 1

Diagram representing the parallel factor analysis model. VF = visual field.

Fig. 2

Goodness of fit (%) with varying number of components (the number under each point denotes the number of components included in the analysis).

Spatial mode

The VF map with corresponding spatial scores of each component is presented in Fig. 3A–3D. Because raw data were centered prior to analysis, scores can be interpreted as the deviation from the m ean. For every V F test, the signs and scales of spatial scores of each component were analyzed after modifying the corresponding subject and time mode loadings. Hence, the signs of scores in the component can be inversed, and comparison of magnitude of each scale should be done only within the component.

Fig. 3

Spatial scores of three components. (A) Component 1, (B) component 2, and (C) component 3 are arranged to match visual field locations (both eyes presented as the right eye). (D) Reference map of visual field locations. SN = superonasal; ST = superotemporal; IN = inferonasal; IT = inferotemporal.

Component 1 clearly shows a contrast between superior and inferior hemispheres, suggesting that vertical hemispheric differences are the most prominent feature of glaucomatous VF defects in the study population. Scores are generally smaller in temporal areas, and larger along inferior and superior arcuate locations (superonasal SN01-04, SN10-13, superotemporal ST10-11; inferonasal IN10-13, IN20-22, and inferotemporal IT10-11, IT20). The spatial patterns of components 2 and 3 are less clear. In component 2, positive scores are mostly located in superior (superior rows 2, 3) and temporal (temporal columns 2, 3) peripheral areas, whereas nasal negative scores are mostly located in nasal (SN01-04, IN01-04, and IN11-13) and central areas (SN00 and ST00). Component 3 is slightly C-shaped, with a mainly horizontal contrast, having scores that are mostly nasally located except for two central (SM00 and IN00) negative signs, and the strongest positive scores in superior central locations (SN00 and ST00-01). Components 2 and 3 might represent other aspects of glaucomatous VF defects or might indicate normal variability.

Temporal mode

The temporal loadings of component 1 increased with the passage of time, showing a clear linear pattern. Components 2 and 3 did not show the same trend (Fig. 4A–4C).

Fig. 4

Temporal loadings of three components. (A) Component 1 showed clear linear pattern. (B) Component 2 and (C) component 3 did not show the same trend (time 1 = initial follow-up, time 6 = last follow-up).

Subject mode

The mean loadings differed significantly between groups in all three components, with differences being most prominent in component 1 (Table 2). Distribution in PPG subjects was more concentrated near zero, whereas PG subjects had wider and subnormal distributions in all three components, with outstanding variance in component 1 (Fig. 5A–5C). To further assess structural differences in component loadings between the groups, we calculated the correlation between components within subjects. Within PPG subjects, component 2 had highly significant correlations with component 1 (0.296, p < 0.01) and component 3 (−0.516, p < 0.01). In contrast, components 1 and 3 had no significant correlation. No pair of components with a high correlation was detected in PG subjects (Table 3).

Table 2

Distribution of subject loadings in PPG and PG

Values are presented as mean ± standard deviation (range).

PPG = pre-perimetric glaucoma; PG = perimetric glaucoma.

*t-test with equal variances not assumed; †p < 0.05 considered as statistically significant.

Fig. 5

Histogram of subject loadings by group. PG subjects had wider and subnormal distributions compared to PPG subjects in (A) component 1, (B) component 2, and (C) component 3. (A) Component 1 showed clear variability. PPG = pre-perimetric glaucoma; PG = perimetric glaucoma.

Table 3

Pearson correlation analysis of components in PPG and PG

PPG = pre-perimetric glaucoma; PG = perimetric glaucoma; C1 = component 1; C2 = component 2; C3 = component 3.

Discussion

Various methodologies have been developed for analyzing multi-dimensional data in neuroscience and are being extended to fields of clinical diagnosis [2930313233]. Applications in glaucoma have not been found in existing literature to date. Our PARAFAC/CANDECOMP analysis is a data-driven approach, drawing conclusions from the data itself rather than from pre-defined strict assumptions and models. Thus, this new approach enables simultaneous analysis of spatial, temporal and subject-specific characteristics, providing a more comprehensive and unbiased understanding of glaucomatous VF.

Turning to the data herein, component 1 had outstandingly high variability among three components, suggesting the importance of component 1 in explaining glaucomatous VF defects in the study population. On the spatial map, component 1 showed a clear vertical contrast, (i.e., differences between superior and inferior areas), which is a well-known feature of glaucomatous VF defects. In addition, component 1 showed a clear linear trend over time, suggesting increasing contrast between superior and inferior hemispheres in the study population over time. Such a clear temporal trend confers the possibility of developing a predictive scheme based on multi-mode data analysis, predicting not only the general rate of decline but also the patterns of progression.

The sharp contrast in distributions in component 1 loadings between PPG and PG shows the potential of multi-dimensional decomposition methodology for discerning glaucomatous VF defects from other types of VF defects. Patients with PPG in our study included those without apparent glaucomatous VF defects, but with clinical features of glaucomatous optic discs. Thus, it is not clear whether PPG results represent normal or early glaucomatous changes in retinal ganglion cells with non-significant VF defects [15]. Further prospective research including presumed normal controls and PPG patients might reveal novel characteristics of very early glaucomatous VF damage, thereby facilitating earlier diagnosis.

In PPG patients herein, component 2 is correlated with components 1 and 3. This relationship may be explained insofar as components 2 and 3 compensated for component 1 in patients with a pre-perimetric status. In PG patients, however, where the variability of loadings was much higher, the relationship becomes vague. There is a possibility that a different complementary set of components existed in the two groups. Further research with longer follow-up time and more components is necessary.

More detailed identification of glaucomatous VF patterns was not possible due to insufficient data and the limited number of components. With more data and sophisticated methodologies with more patterns, sub-patterns of glaucomatous VF defects might be identified along with pattern-specific progression rates.

Collecting large amounts of reliable VF data is difficult due to 1) the long timescale of VF progression, 2) limited follow-up time because of various clinical and personal factors, 3) variability in interval periods between subsequent tests, and 4) treatments and interventions affecting progression rates and patterns [34].

With the recent development of precise diagnostic modalities, intense research has emerged in assessing clinical implications of sophisticated measurements in the field of glaucoma—specifically in optical coherence tomography [3536]. In this context, we expect that multi-dimensional approaches will find wider use. Current methodologies include only VF retinal sensitivity data, but multi-way decompositions can be extended to include more dimensions, namely demographic data (age, sex, and other characteristics), physical examination data (intraocular pressure, cataract and angle grading, and other physical features), and other diagnostic modalities, such as optical coherence tomography, by combining spatial and structural information.

We expect that the future development of more suitable models and collection of longer-term follow-up data will overcome many of the existing difficulties described.