Using Unsupervised Learning with Independent Component Analysis to Identify Patterns of Glaucomatous Visual Field Defects

Abstract

purpose. Clustering by unsupervised learning with machine learning classifiers was shown to segment clusters of patterns in standard automated perimetry (SAP) for glaucoma in previous publications. In this study, unsupervised learning by independent component analysis decomposed SAP field patterns into axes, and the information represented by these axes was evaluated.

methods. SAP fields were used that were obtained with the Humphrey Visual Field Analyzer (Carl Zeiss Meditec, Dublin, CA) from 189 normal eyes and 156 eyes with glaucomatous optic neuropathy (GON) determined by masked review with stereoscopic optic disc photographs. The variational Bayesian independent component analysis mixture model (vB-ICA-mm) partitioned the SAP fields into the most informative number of clusters. Simultaneously, the model learned an optimal number of maximally independent axes for each cluster.

results. The most informative number of clusters in the SAP set was two. vB-ICA-mm placed 68.6% of the eyes with GON in a cluster labeled G and 98.4% of the eyes with normal optic discs in a cluster labeled N. Cluster G optimally contained six axes. Post hoc analysis of patterns generated at −1 SD and +2 SD from the cluster G mean on the six axes revealed defects similar to those identified by experts as indicative of glaucoma. SAP fields associated with an axis showed increasing severity, as they were located farther in the positive direction from the cluster G mean.

conclusions. vB-ICA-mm represented the SAP fields with patterns that were meaningful for glaucoma experts. This process also captured severity in the patterns uncovered. These findings should validate vB-ICA-mm as a data-mining technique for new and unfamiliar complex tests.

This study focuses on the problem faced by physicians in interpreting tests for the diagnosis and management of diseases. The discipline of glaucoma diagnosis is particularly rich with data-intense tests that are currently being applied or developed to determine the status of glaucoma in an eye. Diagnosing and staging glaucoma both necessitate the use of several tests, including structural and functional assessment of the optic nerve. One of these, visual field testing, has been used by experts for several decades as the prime indicator for both diagnosis of glaucoma and the detection of progression. This wide expertise can be used to evaluate the ability of machine learning classifiers to find information that can be valuable in clinical practice. The hypothesis proposed is that machine learning classifiers can learn useful information about visual field patterns. This information can be validated against the prevailing experience of experts. Such validation would give credence to the application of the same process to new and unfamiliar tests. 

Visual fields have been used to diagnose glaucoma since 1856. ¹ In 1889 Bjerrum began to uncover patterns of visual field defects, such as a comet-shaped arcuate scotoma or a nasal step scotoma, with quantitative perimetry. ² Consequently, for more than 110 years, generations of experts in glaucoma have accumulated knowledge to recognize patterns of visual field defects that indicate glaucoma. With the advent of standard automated perimetry (SAP) and the development of statistical field analysis packages such as Statpac (Carl Zeiss Meditec, Inc., Dublin, CA), the depth of defect within the patterns of loss and the relationship of adjacent test locations to each other could be quantified. We applied machine-learning data-mining techniques to uncover visual field patterns associated with glaucoma and to compare these patterns with the subjective qualitative and semiquantitative patterns evolved from experience by human experts. 

There are two ways that classifiers can learn from data concerning medical conditions. One way is to learn to diagnose disease, predict outcomes, and look for change. ^{3 4 5 6 7 8} Another way is by finding useful patterns in large groups of patients. In the present study, we applied a variational Bayesian independent component analysis mixture model (vB-ICA-mm) to subdivide the results of SAP field tests performed in a cohort of normal subjects and patients with glaucomatous optic neuropathy (GON), to find patterns of interest in the SAP fields, and to propose how this knowledge can improve medical care. 

A classifier may learn from training examples the knowledge it needs to make decisions. Specifically, a classifier may learn to diagnose glaucoma by distinguishing abnormal SAP fields from normal fields after being trained with a set of SAP fields labeled with the correct diagnosis. This process is called learning with a teacher or supervised learning. 

When the diagnosis is not supplied with each of the training samples, the objective of the learning algorithm no longer is classification or diagnosis. The goal can be to organize the input data into meaningful structures or groups of patterns. This process is called learning without a teacher or unsupervised learning. For our purposes, the goal can be to organize data, such as the set of visual fields, into differing clusters of patterns with similar members. Hence, another name given to the process is clustering, and the investigation of the clusters segmented from the data is cluster analysis.  

Machine learning classifiers impose fewer constraints on the data than do statistical classifiers, permitting better adaptation to the data and thus better organization. Unsupervised learning with these newer classifiers has the potential of grouping the patterns in the data in a manner that is more useful than that achieved by statistical methods. 

The decomposition of the data by different types of unsupervised learning yields different structures. Instead of clusters, in the present study we evaluated a different structure that relies on axes. Component analysis projects the data within each cluster located in multidimensional space onto axes that meaningfully represent the data. Although the axes do not produce clusters, representation of the data with axes still yields useful information about the patterns in the data. A further refinement, principal component analysis (PCA) projects d-dimensional data on a lower dimension subspace of s-orthogonal axes. Though the expectation is that the orthogonal axes are independent, in reality, the axes may not be independent. We chose independent component analysis (ICA), because it incorporates a measure of independence to produce axes that are maximally independent. ^{9 10} There may be data distributions in which components are nonlinearly related or clustered so that they are difficult to describe by one ICA model. The ICA mixture model ⁹ is a nonlinear ICA technique that extends the linear ICA method by learning multiple ICA models and weighting them in a probabilistic manner. The ICA mixture model settles on the optimal number of axis sets simultaneously with the generation of the axes (Figs. 1A 1B) . 

Representing all the data (e.g., normal and glaucomatous eyes together) with a single set of axes may produce suboptimal representation of the data. Multiple sets of axes are more likely to match the local conditions (which, in this study, turned out to be a set mostly of normal eyes and a set mostly of eyes with glaucoma). Instead of relying on a chosen number of clusters with fixed dimensionality, the ICA mixture model learns the dimensionality and number of classes. The exact computation of the marginal likelihood is computationally intractable; thus, we used variational Bayesian approximation techniques. ¹¹ The variational Bayesian framework helps to capture the number of axes in the local axis set and reduces the computational complexity (by bounding intractable integrals). The amalgamation of all these processes is the variational Bayesian ICA mixture model (vB-ICA-mm) ⁶ that we applied in this study to fields in healthy eyes and glaucomatous eyes to identify patterns of field loss associated with this diagnosis. 

Methods

Participant Selection and Testing

Participants.

We used the visual field set from the first 12 years of the ongoing National Eye Institute–sponsored longitudinal Diagnostic Innovations in Glaucoma Study (DIGS) of visual function in glaucoma. Normal participants in this study were recruited from the community, staff, and spouses or friends of patients. Patients with primary open-angle glaucoma were recruited from the Hamilton Glaucoma Center. Informed consent was obtained from all participants. The study was approved by the Institutional Review Board of the University of California at San Diego, and it adhered to the tenets of the Declaration of Helsinki. 

Exclusion criteria for both groups included unreliable visual fields (defined as a value more than 33% for fixation loss, false-negative errors, or false-positive errors), ¹² angle abnormalities on gonioscopy, diseases other than glaucoma that could affect the visual fields, and medications known to affect visual field sensitivity. Subjects with a best-corrected visual acuity worse than 20/40, spherical equivalent outside ±5.0 D, and cylinder correction greater than 3.0 D were excluded. Poor quality stereoscopic photographs of the optic nerve head also served as an exclusion for the glaucoma category. A family history of glaucoma was not an exclusion criterion. 

Inclusion criteria for the normal category required that subjects have normal findings in dilated eye examinations, open angles, and no evidence of visible GON. Normal optic discs had a cup-to-disc ratio asymmetry ≤0.2; intact rims without hemorrhages, notches, or excavation; and an absence of nerve fiber layer defects in the adjacent retina. Normal subjects had intraocular pressures (IOP) ≤22 mm Hg and no history of elevated IOP. If both of the eyes met the inclusion criteria, one of the eyes was selected at random, to ensure independence between eyes. 

GON as an Indicator of Glaucoma.

Because the goal was to analyze SAP visual fields, visual fields were not used to determine whether an eye was glaucomatous. The classification of an eye as glaucomatous or normal and the labeling of its visual field for the post hoc analysis of the results were based on the appearance of the optic disc. The designation of the optic disc as glaucomatous or normal was accomplished with masked evaluations by two independent graders of a stereoscopic disc photograph taken within 6 months of the visual field test. Inconsistencies between the graders’ evaluations were resolved by consensus or through adjudication by a third masked evaluator. Color simultaneous stereoscopic photographs were obtained (TRC-SS camera; Topcon Instrument Corp. of American, Paramus, NJ) after maximal pupillary dilation. Stereoscopic disc photographs were evaluated for all eyes, with the exception of a subset of normal subjects (95 eyes) in whom photography had not been performed early in the collection of the normal database. All normal subjects had no evidence of optic disc damage during dilated slit lamp indirect ophthalmoscopy with a hand-held 78-D lens. The final selection of eyes totaled 345, including 189 normal eyes (mean age, 50.0 ± 6.7 years [SD]) and 156 eyes with GON (mean age, 62.3 ± 12.4 years). 

Visual Field Testing.

All subjects had automated full-threshold standard visual field testing with the Humphrey Visual Field Analyzer (HFA; Carl Zeiss Meditec, Inc.) with program 24-2 or 30-2. Though most clinicians use the 24° 24-2 program to test for glaucoma, some of the subjects were in other studies that required field testing to be performed with the 30° 30-2 program. The visual field locations in the 30-2 fields that are not in 24-2 fields are at the edge of the field and do not contribute much to the accurate diagnosis of glaucoma. These locations were deleted from the 30-2 field data and displays, to make the SAP data consistent. 

Representation of Data

Input for the Classifier.

The absolute sensitivity (in decibels) of the 52 visual field locations (L) plus age formed a vector in 53-dimensional input space for each of the 345 SAP fields of normal and glaucomatous eyes. The 52 threshold values were extracted from the HFA by computer (Peridata ver. 6.2; Peridata Software GmbH, Hürth, Germany). Each feature vector was x = (L1, … , L54, age), excluding locations L18 and L31, because they fell in the normal blind spot. Age was included because both normal and glaucomatous SAP fields are affected by age, and age had been used in some studies that incorporated supervised learning. ^{3 4 5 7 8}  

Partitioning Data into Clusters and Adjusting Axes.

The unsupervised learning was performed with the vB-ICA-mm. This method is the core of the new direction in data exploration we propose. A detailed, mathematically rigorous description of this method is available in the ^Appendix. A general description follows. Starting with 345 subjects evenly distributed in c clusters and with random initialization of axes, the axes and the probability of each cluster were learned. The number of clusters, c, was increased from c = 1 until c = 5, seeking the value beyond which no further gain in information would be obtained. Subjects belonging to the same cluster defined the mean and templates for that cluster. For each SAP field, its cluster assignment was then recomputed according to its likelihood value given by ICA. ICA alone would have sought one set of axes. The mixture model of ICA created clusters and sought an optimal set of axes for each cluster. Bayesian learning was applied to overcome overfitting by maximum-likelihood estimation. The cluster assignments and axis adjustments were iterated until no further change occurred in the cluster assignment. ^{6 9} The possibility of ending up in a high local minimum of the error surface was minimized by selecting the model with the best marginal likelihood value from 100 different random initializations, thus ending up as close to the global minimum as possible. 

Post Hoc Analysis

Validation of Structure.

The process was validated by observing whether the structure obtained by vB-ICA-mm was appropriate for the data. As described in the Results section, the data were best represented by two clusters. Because the data contained visual fields from normal and glaucomatous eyes, we first evaluated the two clusters, created without knowledge of the diagnosis, for their proportion of normal and glaucomatous eyes. The second validation came from evaluation of the SAP patterns in each cluster represented by the axes uncovered by ICA. The representation by these axes of patterns similar to those identified by generations of glaucoma experts through more than a century of accumulated experience would validate the structure uncovered by unsupervised learning with machine learning classifiers. 

Placing Fields around an Axis.

Experience with cluster analysis methods that perform only clustering demonstrates that analysis of members inside the clusters gives us some understanding of the data organization achieved by unsupervised learning. ¹³ Instead of creating subclusters, as was reported with the variational Bayesian mixture of factor analysis mixture model, ¹³ vB-ICA-mm creates axes within the clusters. The axes were analyzed by generating patterns in SAP field input space at specific points along the axes created by the vB-ICA-mm. We did not want to ignore potential information in the SAP fields. Information was also sought by examining the SAP fields in addition to examining patterns generated at points on the axes. The SAP fields in a cluster generated by the vB-ICA-mm were thus organized into clouds around the axes created by the vB-ICA-mm, and these clouds were analyzed. 

To analyze the individual SAP fields in cluster c, each field was assigned to an axis. In 53-dimensional space, the angle was calculated at the cluster c centroid between the vector for any individual SAP field and the vector of each of the axes. The individual field was assigned to the axis with which the SAP field vector had the smallest angle. This created a cloud of points in space around each axis, with the axis running up the center of the cloud (Fig. 2) . One can imagine the shape of each cloud to approximate a “hypercone” (multidimensional manifestation of a cone in three-dimensional space) expanding away from the glaucoma mean. Each axis had a positive cloud and a negative cloud. Each field was assigned to one axis in one direction, and there was no overlap in the clouds. This allowed us to look at the fields falling within each cloud and, we hoped, to determine what was similar about those fields in one cloud and how they differed from fields falling within the other clouds. 

Visual Field Interpretation.

To determine whether the structure disclosed by the unsupervised learning was meaningful to the study of visual fields in glaucoma, we performed a post hoc classification of fields. Two visual field experts (PAS, CB) were masked to the mean patterns generated by the vB-ICA-mm along the positive and negative directions of each of the resultant axes. They were then asked to review the actual grayscale printouts from the SAP visual fields arranged along each axis according to the distance in standard deviation from the mean of cluster c, to discern whether a consistent pattern of field loss was shown in each cloud, and if so, to describe the pattern within each of the resultant groups of fields. Disagreements were resolved by consensus. 

Results

Post Hoc Analysis

After convergence with unsupervised learning, post hoc analysis found cluster 1 to contain 107 eyes with GON and 3 eyes that had normal optic discs. Cluster 2 contained 186 eyes with normal optic discs and 49 eyes with GON. Without knowing the diagnosis during the learning phase, vB-ICA-mm placed 68.6% of the eyes with GON in cluster 1 and 98.4% of the eyes with normal optic discs in cluster 2. For ease of reading, from this point on, cluster 1 will be called cluster G, and cluster 2 will be called cluster N. 

Clusters and Axes Created by the vB-ICA-mm

With a data set of 345 eyes containing SAP fields of 189 normal eyes or 156 eyes with GON, the model that performed best was that with two clusters of one and six axes. To determine the optimal number of axes within each global cluster, the contribution of each axis for reconstructing the input was plotted against the axis number of choices provided by the vB-ICA-mm. In cluster G, the contribution value dropped to near zero after 6 axes and was zero after 12 axes (Fig. 3) . Examination of axes 7 through 12 provided little additional information; these axes may have represented noise. Hence, the most informative number of axes in cluster G was six, which meant that most of the information in the visual field data could be described using six maximally independent axes. In cluster N, the contribution of the axis declined to zero after one axis; hence, only one axis was necessary to represent cluster N. 

Figure 1Bshows in a two-dimensional projection the distribution of clusters G and N collapsed from the 53 dimensions, with the vB-ICA-mm-derived axes superimposed. The two clusters overlapped, and each cluster appeared to have an ellipsoidal distribution. Each axis in cluster G passed through the centroid of that cluster in 53-dimensional space. We related the findings of the vB-ICA-mm to clinical practice by examining specific points along each axis, at −1 and +2 SD away from the centroid (Fig. 4) . The mean pattern generated at the centroid of clusters G is also displayed in Figure 4 . A similar examination of the axis in cluster N was not conducted, because it was the only axis. 

Pattern Display

Simulation of the Total Deviation Plot in Generated Patterns.

Patterns generated at a particular point in the 53-dimensional space are equivalent to absolute sensitivity patterns. Visual field experts do not evaluate the absolute sensitivities plotted in the visual field printout. Instead, they typically rely on the total deviation (TD) or pattern deviation (PD) plots supplied by the Statpac analysis (Carl Zeiss Meditec, Inc.), which take into account the deviation from age-matched healthy eyes (TD) and the effect on the fields of global factors such as cataract (PD). Thus, the numerical TD plot was simulated by subtracting these generated absolute sensitivity patterns from the mean absolute sensitivity pattern at the centroid of cluster N, the cluster holding nearly all the normal eyes. Then, the simulated numerical TD pattern was converted into a grayscale pattern. Another way to consider these derived plots is as 53-dimensional vectors originating at the centroid of the mainly normal cluster N and ending at specified points along any of the six axes in cluster G. Figure 1Cshows graphically how the generated patterns on the axes are made to simulate the TD plot of Statpac. 

Pattern Shift along Each Axis.

The centroid of cluster G from the vB-ICA-mm result could be considered the mean glaucoma visual field (Fig. 4) , and the centroid of cluster N could be considered the mean normal visual field. Although cluster N contained 49 eyes with GON along with 98.4% of the normal eyes, neither the two visual field experts nor supervised and unsupervised learning algorithms noted field abnormalities in most of these 49 eyes. ^{4 13} The mean glaucoma field was not quite uniformly depressed, with greater depression than average at the nasal step zone, a slightly greater depression than average in the superior hemifield, and the least depression just inferior to fixation (Fig. 4) . 

Each axis passed through the mean glaucoma visual field at the centroid of cluster G. The six axes were, by definition, maximally independent. Hence, the patterns represented by each axis were maximally different: the farther along the axis away from the glaucoma mean, the greater the deviation from the glaucoma mean. The axis direction was considered to be positive if the distance from the normal mean field always increased as the shift from the glaucoma mean increased along the axis. As the distance from the glaucoma mean increased in the negative direction, the distance of a point on the axis to the normal mean initially decreased until the minimum distance of the axis to the normal mean was reached (Fig. 1B) . Thereafter, the distance from a point on the axis to the normal mean increased as the distance from the glaucoma mean increased. 

Analysis of Generated Patterns and SAP Fields in Axis-Derived Clouds.

Extraction of Information from Vectors.

The information collected from each axis in cluster G was the generated field pattern on the axis at 0 SD, −1 SD, and +2 SD displayed as TD-like plots along with the actual printed Statpac grayscale plots for each of the SAP fields assigned to the axis (Fig. 4) . The generated plot at 0 SD, being the mean pattern at the centroid of cluster G, was common to all axes. There were no individual SAP fields placed at or near the centroid, with the closest fields being 0.4 to 0.5 SD from the centroid on the negative side of axes 1 to 5. (There were no fields on the negative side of axis 6.) 

On the negative side of the axis, the mean distance from the centroid for all axes was −0.75 SD. Hence, we rounded up and analyzed the generated pattern in this direction at −1 SD. On the positive side, the mean distance from the centroid was 2.50 SD. Thus, we rounded down to analyze the generated pattern on the axis at +2 SD. 

Table 1shows the average, minimum, and maximum distance in standard deviation in the positive and negative directions along each axis from the cluster G mean of the 110 SAP fields assigned to the axes. 

Visual Fields Associated with Each Axis.

We labeled the printed Statpac TD plot for each study eye in cluster G with the axis to which that visual field was assigned (see the Methods section: Placing Fields around an Axis) and with the distance of that field from the glaucoma mean in units of standard deviation. For each axis, the visual field experts evaluated the printed fields ranked from most minus to most plus and looked for the common pattern elements in the fields assigned to that axis. 

Axis 1: At −1 SD, the vB-ICA-mm-generated field had a mild generalized depression that was slightly deeper superiorly and nasally (Fig. 4) . The negative side of this axis had the least correlation between the generated fields and the SAP fields and the least in common among the SAP fields. The 16 SAP fields ranged from −1.5 to −0.5 SD. The common feature of the 16 SAP fields was the tendency to have a single small depression, but that depression varied in location among the fields (superior nasal step, superior paracentral focus, superior arcuate, enlarged blind spot, and temporal wedge). The generalized depression in the generated field represented an average of all the spot depressions in the individual SAP fields on the negative side. At +2 SD, the generated field had a greater generalized depression deepest in both superior and inferior nasal steps and a slightly greater than average depression in the inferior hemifield. The seven SAP fields on the positive side ranged from +0.8 to +6.4 SD. They had multiple deep defects involving three or four quadrants with no consistent patterns. The field defects became larger or deeper as the distance of the field increased from the glaucoma mean. The farthest SAP field at +6.4 SD was too advanced to show a pattern other than severe depression everywhere. 

Axis 2: The generated field at −1 SD had mild inferior hemifield depression with the greatest depression at the nasal steps (Fig. 4) . The six SAP fields ranged from −0.9 to −0.5 SD. They also had mild arcuate defects inferiorly or a focus of depression at a nasal step. The generated field at +2 SD revealed a strong superior hemifield defect augmented nasally. The nine SAP fields ranged from +1.0 to +4.5 SD. They also displayed mostly superior nasal arcuate and hemifield defects. The defects became more severe as the distance from glaucoma mean increased (Fig. 5) . 

Axis 3: The generated field at −1 SD showed a nasal hemifield depression greater along the nasal edge and a little exaggerated at the superior nasal step (Fig. 4) . The eight SAP fields ranged from −1.0 to −0.4 SD. The common pattern elements in TD plots of the fields in the negative direction were nasal step depressions and arcuate depressions at the nasal edge of the 24° SAP field. The generated field at +2 SD revealed superior hemifield with superior temporal arcuate (wedge) depression combined with a superior nasal step. The 10 SAP fields ranged from +0.6 to +4.9 SD. The fields in the positive direction also tended to combine a superior temporal wedge or a superior outer arcuate depression with a depression at the superior temporal wedge. The defects became more severe as the distance from the glaucoma mean increased. 

Axis 4: The generated field at −1 SD showed peripheral depression somewhat greater than central depression (Fig. 4) . The SAP fields in general had a peripheral ring of defect or focal defects scattered about the periphery. The generated field at +2 SD had a line of depression extending from the blind spot to the nasal step, passing just superior to fixation. To contrast this pattern with the positive pattern on axis five, we call this pattern the arrow shaft. The central depression was greater than the peripheral depression. The SAP fields on the positive side ranged from +1.2 to +2.8 SD. The depression pattern was horizontal linear from the blind spot to the nasal step, passing just superior to fixation. The defects became deeper or larger as the SAP field was located farther from the glaucoma mean. 

Axis 5: The generated field at −1 SD on the negative axis had generalized depression greater in the superior hemifield, centrally, and at the superior nasal step (Fig. 4) . The four SAP fields ranged from −0.6 to −0.4 SD. The patterns were not uniform, with a tendency for more depressions superiorly and at the superior nasal step. The generated field at +2 SD on the positive axis was weighted along the inferior (more) and superior (less) nasal arcuate edge. To contrast this pattern with that on the positive side of axis 4, we called this the arrowhead. The five SAP fields ranged from +1.3 to +6.1 SD. The defects had the arrowhead at the nasal periphery or were inferior nasal. Except for one field, there was a sense of increasing defect as the distance increased from the glaucoma mean. 

Axis 6: The generated field at −1 SD combined mildly greater defects at the superior edge, superior nasal, and superior temporal central (Fig. 4) . There were no SAP fields on the negative side of axis 6. It is not known whether a larger or different data set would have fields in the space defined by the negative side of axis 6 or whether no real fields would ever exist in this region. As with the other axes, the pattern on the negative sides of the axis tended to be the inverse of the pattern on the positive side. The generated pattern at +2 SD on the axis had a linear depression from inferior to fixation to the inferior nasal step. The 11 SAP fields ranged from +1.4 to 4.1 SD. The defects were linear from inferior to fixation to the inferior nasal step or were mostly inferior nasal quadrant. The defects tended to broaden or deepen as the distance increased from the glaucoma mean. 

In summary, the axis patterns generated by the vB-ICA-mm in general were in good agreement with the pattern types identified by experts in this set of fields. 

Discussion

Representation of Data Structure

vB-ICA-mm settled on two clusters. Without knowing the diagnosis during the learning phase, vB-ICA-mm placed 68.6% of the eyes with GON in cluster G (true-positive rate) and 98.4% of the eyes with normal optic discs in cluster N (true-negative rate). In a previous report, ⁴ mixture of Gaussian, a comparable classifier that had knowledge of the correct diagnosis during the learning phase (supervised learning), had a sensitivity of 67% when specificity was restricted to 100%. Consequently, the cluster structure obtained with unsupervised learning by vB-ICA-mm correlated well with that obtained by supervised learning. 

In cluster G, most of the information in the visual field data from patients with glaucoma could be described using six maximally independent axes. The need for six axes means that there were several patterns. In cluster N, one axis was sufficient to describe all the data, indicating that the normal-appearing fields were essentially uniform. 

Change along an Axis

Why should we be interested in generating axes to represent the data when we can generate clusters of patterns that are internally similar but differ from each other? ^{11 14} The generated field patterns and the analyzed SAP fields each had 52 dimensions; with the addition of age, they were in 53-dimensional space. We will discuss the shapes in three-dimensional space, because it is easier to grasp cognitively. Cluster G can itself be segmented into clouds. The elongated clouds that are created around axes within cluster G are likely to have a different type of similarity of its members than clusters close to spherical in shape. Clusters of glaucoma fields created by unsupervised learning with variational Bayesian mixture of factor analysis were somewhat spherical. ¹³ With that type of clustering, it is possible for one spherical cluster to have mostly mild defects and another cluster to have mostly advanced defects. 

Representing the data with axes can give us different insights. The generated patterns represented by the axes differed from each other. The patterns resembled classic glaucoma patterns, such as altitudinal hemifield depression, defects in the arcuate region, reduced sensitivity adjacent to fixation, and temporal wedge. Frequently, there was a combination of classic defects in an axis. The pattern on the negative side of an axis tended to be the inverse of the pattern on the positive side—for example, central versus peripheral, superior versus inferior, or temporal versus nasal. Another observation was for each axis in general to capture not only a unique pattern but also mild and severe forms of that pattern. 

The axes in the glaucoma cluster pass through the glaucoma mean, but an increase in the severity of a field pattern with respect to the normal mean should more closely match the progression of glaucoma away from normal. It is difficult to model the data in a way that captures severity along axes coming out of the normal mean. The severity along ICA axes in relation to the glaucoma mean serves as a surrogate for the sense of severity of a field from normal. Axes that are close to perpendicular to the vector from the normal mean to the glaucoma mean capture the sense of severity the least, and axes that are close to parallel to the vector from the normal mean to the glaucoma mean capture the sense of severity the most. A question arises of whether moving in the positive direction along an axis represents increasing severity of the same pattern, or whether moving in the positive direction shows increasing severity by adding new patterns. This issue is addressed in a separate report that detects worsening of disease by tracking the position of serial fields different from those used in this study. Those fields were tracked within the vB-ICA-mm space developed in this study. ¹⁵ That report describes increasing severity as worsening of the same pattern rather than the addition of new patterns. Another question is the effect of the age difference between the normal and glaucoma participants. In cluster analysis performed on the same participants, the age difference was found not to influence the patterns. ¹³ Because severity, which is incorporated in each axis, is not a factor that distinguishes the axes, the age difference between the normal and glaucoma participants also does not affect the representation of data created with vB-ICA-mm. The representation of patterns is based on glaucoma and not age. 

Validation of Method

vB-ICA-mm produces both clusters and representative axes within the clusters. Without knowing the goals of the analysis of the glaucoma and normal fields, and without the fields being labeled with an indication of glaucoma, the clustering accomplished with vB-ICA-mm separated the SAP fields into two clusters. Cluster N had 98.4% of the SAP fields from normal eyes, and cluster G had 68.6% of the fields from eyes with GON. If the goal had been classification instead of data organization, these values would have represented specificity and sensitivity, respectively. The cluster structure obtained with unsupervised learning by vB-ICA-mm correlated well with that obtained by supervised learning. 

The patterns generated along the ICA axes matched the descriptions of the SAP fields provided by the two visual field experts in this study. More important, the patterns generated along the ICA axes were similar to those discovered by generations of human experts with decades of experience analyzing visual fields in eyes with glaucoma. This finding validates the vB-ICA-mm as a method of separating out patterns of defects consistent with glaucoma. 

In summary, the vB-ICA-mm process can be relied on to analyze a complex test for glaucoma in a meaningful way. The data were separated with unsupervised learning by vB-ICA-mm into normal and glaucoma clusters as well as could be accomplished with machine learning classifiers trained with supervised learning on the diagnosis and by glaucoma experts. The patterns disclosed by the axes representing the glaucoma cluster correlated to patterns discovered by years of human experience with glaucomatous visual fields. The validation of vB-ICA-mm on tests we understand indicates that this approach to data mining should teach us the patterns of significance in tests we do not yet understand. Last, the type of data representation achieved by ICA axes captures the sense of severity, giving this type of analysis the potential of detecting increasing severity of disease in a quantitative manner. ¹⁵  

Appendix 1

The application of the variational Bayesian ICA mixture model in this project is adapted from the formal presentation in Chan. ⁶ vB-ICA-mm automatically determines the number of clusters and dimensions (axes) of each cluster. The input data (52 field locations plus age for each eye) are denoted by x, in 53-dimensional space, and vB-ICA-mm models the data density by clusters, p(x) = \({{\sum}_{c}}P(c)p(\mathbf{x}{\vert}c)\) , where P(c) is the probability mass of clustercandp(x|c) is the probability density ofxwithin clusterc. Within each cluster \(c\) , the data are modeled by the linear combination of independent sources,x ^c = A ^c s ^c +υ ^c +ς ^c , whereA ^c = (A ₁ ^c ,A ₂ ^c , ⋮ ,A _n ^c ) is the mixing matrix of the independent axes (A ₁ ^c , ⋮ ,A _n ^c),s ^c = (s ₁ ^c , ⋮ ,s _n ^c ) ^T are the activation coefficients along the axes, and \(\mathbf{{\upsilon}}^{\mathit{c}}\) is the centroid, \(\mathbf{{\varsigma}}^{\mathit{c}}\) the noise, and n the dimensionality (number of axes) of cluster c. The noise is modeled by Gaussian distribution with zero mean and covariance \(\mathbf{{\psi}}^{\mathbf{c}}\) . So the distribution of x in any cluster can be written asp(x|c) = p(x|A ^c ,υ ^c ,ψ ^c) = ∫N(x|A ^c s ^c + υ ^c,ψ ^c)p(s ^c)ds ^c, where \(N\) denotes Gaussian distribution, the activation coefficients \(\mathbf{s}^{c}\) are assumed to be independent, and the density of each source, \(s_{m}^{c}\) , is modeled by k mixtures of Gaussian, p(s_m ^c) = \({{\sum}_{k}}{\pi}_{mk}^{c}N(s_{m}^{c}{\vert}{\phi}_{mk}^{c},\ {\beta}_{mk}^{c}),\) where each \({\pi}_{mk}^{c}\) is a mixture weight and \(N\) denotes Gaussian distribution whose mean is \({\phi}_{mk}^{c}\) and variance is \({\beta}_{mk}^{c}\) . The prior for mixing matrix \(\mathbf{A}^{c}\) is also Gaussian with zero mean and covariance \(\mathbf{{\alpha}}\) : \(p(\mathbf{A}_{nm}^{c}{\vert}{\alpha}_{m}^{c}){=}N(\mathbf{A}_{nm}^{c}{\vert}0,{\alpha}_{m}^{c}).\) Often, maximum likelihood estimation overfits data, and Bayesian learning overcomes overfitting by introducing priors. The priors introduced for the parameters π, φ, β, α, υ, ψ, P(c) are D, N, Γ, Γ, N, Γ, D, respectively, where N, Γ, and D are Gaussian, Gamma, and Dirichlet distributions, respectively. 

To illustrate variational Bayesian learning, it is convenient to collect all the parameters and call them θ = {π, φ, β, α, υ, ψ, P(c)}. Given all priors for these parameters, denoted by p(θ), the marginal likelihood is given by p(x) = ∫p(x|θ)p(θ)d θ. This Bayesian approach automatically performs model selection, but exact Bayesian learning is rarely computationally tractable, because it is hard to obtain the posterior distribution of parameters p(θ|x). We need a simpler function, q(θ), to approximate the true posterior of parameters, and the log marginal likelihood is lower bounded by  

\[\mathrm{log}p(\mathbf{x}){\geq}\ {{\int}}q(\mathbf{{\theta}})\mathrm{log}p(\mathbf{\mathrm{x}}{\vert}\mathbf{{\theta}})d\mathbf{{\theta}}{+}\ {{\int}}q(\mathbf{{\theta}})\mathrm{log}\frac{p(\mathbf{{\theta}})}{q(\mathbf{{\theta}})}d\mathbf{{\theta}}.\]

 

We introduced separable distribution over the parameters, and the closed form of \(q(\mathbf{{\theta}})\) was obtained, but with different parameters. The learning of vB-ICA-mm was accomplished as follows. For each number of clusters, c (c = 1,2,… 5), each number of axes, m (m = up to 20), vB-ICA-mm did the following: 

Step 1: initialize the parameters of \(q(\mathbf{{\theta}})\) . 

Step 2: Derive the new axes by re-estimating the parameters of \(q(\mathbf{{\theta}})\) with \(q(\mathbf{{\theta}})\) fixed. 

Step 3: Based on the new axes, recalculate the cluster probability of each datum and update the parameters of \(q(\mathbf{{\theta}})\) . 

The vB-ICA-mm is then iterated between steps 2 and 3 until the axes and the cluster probability stop changing (convergence). 

Readers are referred to Chan et al. ⁶ for further details of the functional form for priors, close functional form of q(θ), and the learning rules. 

As a local method, vB-ICA-mm is likely to get stuck at local minima. Random initialization is a way to overcome the local minima problem by permitting selection of the best local minimum. Consequently, the vB-ICA-mm was set to repeat the initial randomization, as in step 1, 100 times, so that for each c and m, we had 100 models. All the steps were repeated while simultaneously varying the number of clusters and the number of axes. We chose the final model by comparing all the models based on their marginal likelihood values (the larger the value, the better) and the classification accuracy. 

As described in the results section, the best model we got for this analysis was with two clusters: one cluster with only one axis and the other cluster with six axes. vB-ICA-mm assigned a posterior probability, \(p(c{\vert}\mathbf{x})\ {\propto}\ p(\mathbf{x}{\vert}c)P(c)\) , to each SAP field, which enabled the program to determined the class of SAP with maximum posterior probability. Each cluster was then examined and labeled according to most of the GON and normal data points (eyes) within. The cluster with one axis that contained mostly normal eyes was called cluster N, and the cluster with six axes that contained mostly eyes with GON was called cluster G. Within cluster G, each data point was projected onto the six axes and was assigned into the cloud of the axis with which the vector for that data point had the smallest angle at the centroid of cluster G. On a desktop computer with 2GH Pentium 4 processor (Intel, Mountain View, CA) running MatLab (The MathWorks, Natick, MA) using the ICA code, training time was approximately 15 minutes for generating clusters and axes. 

 This article is based on a thesis that was prepared in partial fulfillment of the requirements for membership in the American Ophthalmological Society.

 Supported by National Eye Institute Grants EY13235 (MHG) and EY08208 (PAS) and funds from Research to Prevent Blindness (PAS) and the Howard Hughes Medical Institute (TS).

 Submitted for publication October 1, 2004; revised March 9, 2005; accepted June 28, 2005.

 Disclosure: M.H. Goldbaum, None; P.A. Sample, Carl Zeiss Meditec, Inc. (F); Z. Zhang, None; K. Chan, None; J. Hao, None; T.-W. Lee, None; C. Boden, None; C. Bowd, None; R. Bourne, None; L. Zangwill, None; T. Sejnowski, None; D. Spinak, None; R.N. Weinreb, Carl Zeiss Meditec, Inc. (F)

 The publication costs of this article were defrayed in part by page charge payment. This article must therefore be marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

 

Figure 1.

 

(A) Two-dimensional scatterplot of participants with (+) and without (○) GON. Units are in standard deviations away from the mean of cluster G. (B) Two-dimensional projection of 53-dimensional space. vB-ICA-mm converged to two clusters: named from post hoc analysis the glaucoma cluster (cluster G) and normal cluster (cluster N). The glaucoma cluster decomposed into six axes through the glaucoma mean. Axes 1 to 3 are shown; axes 4 to 6 are not depicted, to avoid cluttering. The positive direction of each axis moves farther from the normal mean. (C) Relationship of points on axes in cluster G to the mean of the normal cluster. Large arrows: vectors from the cluster N mean to a point on the negative side and a point on the positive side of axis 2. Grayscales at the tips of these vectors simulate the grayscale of TD plots. Increasing severity away from the normal mean most closely matches increasing deviation from the glaucoma mean on the positive side of the glaucoma mean.

Figure 2.

 

Distribution of SAP fields assigned to an axis. Each axis has a cloud around the positive axis and a cloud around the negative axis. The broken lines indicate the boundary surface for the cloud. Individual absolute sensitivity plots are distributed around the axis they are closest to, inside the cloud on either side of the glaucoma cluster mean. The vector that connects the normal cluster mean with an individual field simulates the TD plot (shown in Fig. 1C ).

Figure 3.

 

Contribution of axis versus axis number. The norm of A _m , on the y-axis represents the contribution of axis m in reconstructing the input. Note the small contribution of axes 7 through 12. There was no contribution after 12 axes.

Figure 4.

 

Generated patterns on axes. Patterns generated at cluster G (mostly GON) mean, −1 SD on each axis in cluster G, and +2 SD on each axis in cluster G, with the plots displayed as if they were grayscale images of TD plots of the left eye.

Table 1.

 

The Average, Minimum, and Maximum Distance in Standard Deviation in the Positive and Negative Directions along Each Axis from the Glaucoma Mean of the 110 SAP Fields Assigned to the Axes in the Glaucoma Cluster

	Axis											Total
		1	2	3	4	5	6
Negative side of glaucoma mean
n	18	7	8	15	4	0	52
Average	−0.93	−0.75	−0.72	−0.79	−0.56
Min; max	−1.5; −0.5	−1.0; −0.5	−1.0; −0.4	−1.3; −0.4	−0.6; −0.4
Positive side of glaucoma mean
n	7	9	10	16	5	11	58
Average	+1.76	+2.63	+2.50	+1.78	+3.86	+2.49
Min; max	+0.8; +6.4	+1.0; +4.5	+0.6; +4.9	+1.2; +2.9	+1.4; +6.1	+1.4; +4.1
Total	25	16	18	31	9	11	110

Figure 5.

 

Demonstration of severity along an axis. Representative grayscale of two SAP fields from axis 2 showing increasing defect as the distance from cluster G mean increases in the positive direction. Notice that the SAP field grayscale matches the pattern of the generated field at +2 SD in Figure 4and that the +4.2-SD pattern is deeper but not different in shape from the +1.0-SD pattern.