## Abstract

Glaucoma is an optic neuropathy accompanied by vision loss which can be mapped by visual field (VF) testing revealing characteristic patterns related to the retinal nerve fibre layer anatomy. While detailed knowledge about these patterns is important to understand the anatomic and genetic aspects of glaucoma, current classification schemes are typically predominantly derived qualitatively. Here, we classify glaucomatous vision loss quantitatively by statistically learning prototypical patterns on the convex hull of the data space. In contrast to component-based approaches, this method emphasizes distinct aspects of the data and provides patterns that are easier to interpret for clinicians. Based on 13 231 reliable Humphrey VFs from a large clinical glaucoma practice, we identify an optimal solution with 17 glaucomatous vision loss prototypes which fit well with previously described qualitative patterns from a large clinical study. We illustrate relations of our patterns to retinal structure by a previously developed mathematical model. In contrast to the qualitative clinical approaches, our results can serve as a framework to quantify the various subtypes of glaucomatous visual field loss.

## 1. Introduction

Glaucoma [1], one of the leading causes of blindness worldwide, is an optic neuropathy accompanied by both characteristic optic nerve damage and the loss of the spatial array of sensitivity, so-called *visual field* (*VF*) *loss*. Typically, the VF is assessed by standard automated perimetry in which patients respond to light dots presented at fixed locations in the field of vision by a machine (perimeter) which determines a sensitivity threshold for each location. Glaucomatous VF loss is not random, but rather follows specific patterns which are related to the retinal structure. Axons of retinal ganglion cells form bundles of retinal nerve fibres (RNFs) which project to specific areas of the optic disc. Defects in RNF trajectories and on the respective locations of their projections to the optic disc can be related to specific glaucomatous functional loss.

While there have been many suggestions for the classification of VF loss [2], classifications that investigate detailed *subtypes* are rare. Such subtypes, however, are important for studies investigating structure–function relationships, including optic nerve fibre layer anatomy or the genetic backgrounds of glaucoma. For example, research in glaucoma genomics has found three gene regions associated with vulnerability to paracentral scotomas in primary open angle glaucoma: p53 [3], the CAV1/CAV2 region [4] and the GUCY1A3/GUCY1B3 regions [5].

A detailed *qualitatively derived* classification scheme has been proposed in the course of the ocular hypertension treatment study (OHTS) which was a large-scale randomized clinical trial in which VFs of patients with raised intraocular pressure (ocular hypertension) were investigated. Ocular hypertension is known as a major risk factor for glaucoma [6]. In the course of performing the OHTS, a classification of VF abnormalities was reported [7]. The authors visually inspected abnormal VFs and qualitatively developed a classification scheme of 17 mutually exclusive categories, divided into two groups: first, *nerve fibre bundle abnormalities*, which are likely to be glaucomatous and are expected to be strongly represented in our glaucoma service VFs, and second, *non-nerve fibre bundle abnormalities*, which include probable artefact abnormalities as well as abnormalities owing to ophthalmic diseases other than glaucoma.

The OHTS patterns are defined over verbal descriptions and have been derived using ophthalmic background knowledge. A very different approach to identify patterns is *unsupervised statistical learning*, which means, computer algorithms are applied to large bodies of VF data in order to autonomously (hence, unsupervised) categorize and classify the data. These methods are typically less biased by prior knowledge concerning the nature of glaucomatous field loss. Herein, we summarize and categorize previously applied statistical learning approaches to VFs and compare them with the approach introduced in this work. Where appropriate, we apply the approaches to our large VF dataset which will be described in §2.1. To visualize general principles of the statistical learning approaches, we additionally illustrate some of the approaches by simple two-dimensional hypothetical examples for better understanding.

### 1.1. Two approaches: prototypes and components

Unsupervised statistical learning procedures to characterize subtypes or patterns of a dataset can be classified into two categories. One group of approaches attempts to obtain partitions and their respective *prototypes*, that is, representative patterns for each partition. A typical family of procedures for this purpose is cluster analysis, where the data samples are grouped into clusters within each of which the similarity of samples is higher than between separate clusters. The OHTS scheme described above [7] implicitly assumes that VFs can be partitioned into such groups as ‘arcuate’, etc., and the verbal descriptions of these classes have the character of prototypes.

Another group of approaches tries to unveil structure in the data, such as dimensionality or latent sources. These procedures explore *components* of the data. In the component approaches applied to VFs, these components are typically represented as axes in the data space.

Figure 1 summarizes and illustrates these two approaches. In the simple hypothetical example (*a*), the prototype approach is illustrated by the popular cluster analysis algorithm *k*-means [8], an extension of which has recently been applied to VFs [9]. In our example, *k*-means partitions the cloud of points into four clusters. The cluster means act as representative patterns (prototypes). In (*c*), we apply *k*-means to our VF dataset (see §2.1). We show the cluster centres of a four cluster solution. The colours denote total deviations (TDs; red, sensitivities below normal values; blue, above normal values). The two grey-shaded locations within each VF plot are the locations closest to the blind spot and were omitted. On the bottom left corner of each VF plot, we denote the cluster size that is the percentage of VF measurements in the respective cluster.

The simple two-dimensional example of the component approach (*b*) shows principal component analysis (PCA), which has recently been applied to VF data [10,11]. Formally, PCA applied to VFs yields the orthogonal axes in the 52-dimensional data space which capture the greatest variance. The axes follow the direction of the eigenvectors (EVs) of the covariance matrix, with the tails shifted to the mean. In contrast to the prototype approach, PCA does not aim to cluster the data space and therefore does not yield prototypical patterns for subgroups. In (*d*), we plot the first four EVs of the covariance matrix. Compared with the patterns in (*c*), they do not look like typical TD plots from glaucoma patients, although their appearance covers important structural aspects of the data, such as the global average (EV 1), upper/lower hemifield differences (EV 2), left/right or centre/surround patterns (EV 3 and 4).

Next, we investigate popular component approaches for their applicability to VFs. We represent a VF measurement as a 52-dimensional vector **x** (without loss of generality, we assume that our measurements are normalized to zero mean). We represent our sample of *N* VF measurements by an *N* × 52 matrix **X** the rows of which are the single VF measurements. Typically, we aim to reduce the dimensionality of our data space, that is, we search for *m* < 52 components to represent our data (in figure 1*d*, for instance, we show only the first four of the 52 components). Let **s** be the *m*-dimensional vector of the first *m* components (sorted decreasingly by variance). If we sort the EVs **a*** _{i}* of the covariance matrix of

**X**decreasingly according to their corresponding eigenvalues and represent

**a**

_{1}to

**a**

*as columns of a matrix*

_{m}**A**,

**s**is given by a linear combination of

**x**and

**A**.

PCA can be formulated as a *generative model* for the data. We can represent a VF measurement **x** as1.1where **n** is a noise vector with mean zero, possibly generated by microsaccades, patient positioning, lid anatomy, etc. That means, we represent our data as a linear combination of *m* components **s** and a mixing matrix **A** with additive noise **n**. The representation in (1.1) suggests the interpretation that a VF measurement is composed of latent variables **s**. However, in a linear combination of latent variables, we might want to distinguish between noise components common to all variables and noise components which are unique to each of the latent variables. Formally, this means that the covariance matrix of the noise **n** in (1.1) should be diagonal and of full rank. This model is called *factor analysis* (*FA*) [12]. By contrast, in PCA, the components lie in the linear space spanned by the samples **x**, so that the covariance matrix of **n** does not have these properties. That means, if we assume that glaucomatous VF loss can be caused by latent variables, for instance, particular structural damage of optic nerve fibres, or certain genetic predispositions, FA might be a more promising model than PCA.

It has been frequently discussed that the criteria for the decomposition of **x** into **A** and **s** in (1.1) hold for any *orthogonal transformation*, that is, given any orthogonal matrix **R**, we can set **A** ← **AR** and **s** ← **R**^{T}**s** (**R**^{T}: transpose of **R**) without violating the criteria (e.g. [13, ch. 14.7.1] or [14, ch. 6.3]). In other words, FA and the generative model of PCA are not unique.

Both PCA and (ordinary) FA assume Gaussian distributions of the latent variables. However, in reliability studies, a heteroscedasticity of sensitivity thresholds has been shown, as the dispersion increases with decreasing threshold. Furthermore, with decreasing thresholds, the skewness of the distributions tends to become negative (see, for instance, fig. 4 in [15]). Therefore, instead of PCA and FA, it might be beneficial to consider non-Gaussian procedures. A popular statistical learning method which is related to FA, but requires non-Gaussian distributions of the latent variables is *independent component analysis* (*ICA*) [14], which yields not only uncorrelated but also statistically independent components (see [14, ch. 2.3]). While for Gaussian-distributed components, uncorrelatedness and independence are identical, for all other distributions, statistical independence is a special case of uncorrelatedness, and for non-Gaussian variables, ICA provides unique solutions.

Figure 2*a* illustrates the comparison of PCA and ICA applied to a non-Gaussian set of random points which are uniformly distributed along two non-orthogonal axes. The orthogonality constraint of PCA does not result in a useful representation of the data, whereas the ICA axes follow the orientations of the two sets of random dots.

Note that the algorithms introduced above, no matter whether prototype or component approaches, have in common that they can be represented by the simple linear combination written in equation (1.1). Each VF measurement is assigned an individual vector **s**, which is a ‘component’ in PCA or ICA or a binary vector that contains exactly one value 1 in *k*-means. In addition, there is a mixing matrix **A** which is common to all VF measurements and which determines how the individual ‘components’ are mixed to model the data samples. For PCA, for instance, it contains the EVs of the covariance matrix, or for *k*-means cluster centres.

### 1.2. Hybrid approaches

Recently, attempts were made to combine prototypes and components. Chan and co-workers proposed a hybrid of prototype and component approaches and modelled VFs by clusters of linear mixtures of components. They aimed to determine both the number of clusters and the dimensionality within each cluster. In contrast to maximum-likelihood approaches, which overestimate dimensionality owing to overfitting, Bayesian model selection approaches avoid this problem by integrating out parameters (see, for instance, [16, ch. 5.3]). As a full Bayesian treatment was not tractable, the authors apply an analytical approximation called the *variational* Bayesian (vB) method (e.g. [16, ch. 21]). Their hybrid approach was applied to two component algorithms: FA and ICA. Figure 2*b* illustrates the concept of vB-ICA: ICA is applied to three clusters of random dots, separately, and for each cluster, two axes are learned.

The authors applied these algorithms to 189 normal and 156 glaucomatous VFs. Both algorithms favoured optimal solutions with two clusters which represented normal and glaucomatous fields, respectively. While for the vB-FA approach, no further dimensionality details are reported, the vB-ICA approach yielded a glaucomatous cluster with 12-dimensions and a normal cluster with six-dimensions. Sample *et al.* [17] recalculated vB-FA on the same dataset and propose five instead of two clusters. One of the clusters represented the normal VFs, whereas the other four clusters represented typical patterns of glaucomatous VF loss, namely early defects, deep superior hemifield defects, deep inferior together with less deep superior defects and deep defects in both hemifields, respectively. The internal structure of the five clusters is not further detailed by the authors, which means, although the methodology they apply combines prototype and component approaches, the final result is only a special case of a cluster analysis and therefore a prototype approach.

Subsequently, Goldbaum *et al.* [18] republished the vB-ICA results from Chan *et al.*, but argued that six of the 12 axes of the glaucomatous cluster contributed only weakly to the reconstruction of the input and should therefore be omitted. In a later study [19], they applied the same method to a larger dataset of VFs. Their best solution had three clusters, with cluster 1 (two axes) representing normal VFs, cluster 2 (two axes) representing mild defects and cluster 3 (five axes) representing moderate to severe defects. Figure 3 illustrates TD patterns along these five axes, adopted after data given in figure 3 and table 2 in [19]. The directions of the five axes were determined by the cluster centroid and each of the columns of the mixing matrix **A** in equation (1.1). The patterns are located at ±2 standard deviations (s.d.) from the centroid on each axis.

As summarized above, the influential vB-ICA algorithm was proposed in [20] to assess the inherent dimensionality of glaucomatous VF loss, but has been applied to identify prototypical patterns of field loss [18,19]. However, ICA, which is applied to each cluster separately, yields statistically independent *axes* but not patterns. The patterns proposed by Goldbaum *et al*. have been obtained by selecting two points on each axis, namely the points at −1 s.d. and +2 s.d. in [18] and ±2 s.d. in [19]. On the one hand, as illustrated in figure 3, these patterns cover aspects that are related to glaucomatous VF loss, such as the separation between upper and lower hemifield and typical arcuate shapes. On the other hand, some of the patterns, particularly for −2 s.d., contain positive TDs, that is supernormal VF locations, which would be atypical for glaucomatous fields. Moreover, while ICA provides five independent axes for the severe-to-moderate cluster, it is unclear how to interpret these axes in an ophthalmic context. In figure 3, there is no obvious ophthalmological relation or common principle related to the axes, and the authors only label the single, isolated fields shown at ±2 s.d. of each axis instead of describing the axis itself.

Here, we propose the approach *archetypal analysis* (AA) which originates in statistical physics and engineering [21] and is not based on axes within the data space. In short, as all approaches introduced above, AA can be formalized using equation (1.1). We require for all *N* data samples **x** that all elements of **s** are non-negative and that Furthermore, we write the mixing matrix as **A** = **X**^{T}**B**, with data matrix **X** as introduced above and matrix **B** of dimension *N* × *m*, all of its elements non-negative, and its columns summing up to 1. Given these conditions, AA minimizes ||**X** − **AS**||, where the columns of *m* × *N* matrix **S** are the respective vectors **s** of the *N* samples. Note that, in contrast to the component approaches, the mixing matrix of AA is itself a *linear combination* of the data.

The *m* columns **a**_{1} to **a*** _{m}* of the mixing matrix

**A**constitute the representative patterns and are called

*archetypes*of the data. They are convex combinations of the data samples, which means, in contrast to the respective vectors for PCA or ICA, they are patterns located at the edges of the data space.

For our purposes, a relevant property of archetypes is that they lie on the *convex hull* of the data space. While prototypes from cluster analysis methods like *k*-means also resemble real VFs, each representative pattern is the centre of the respective cluster and therefore averaged over all the features summarized in that cluster. Archetypes, however, are not located at the centres but at the ‘edges’ (convex hull) of the data space, as illustrated by the simple hypothetical example in figure 2*a*: the four archetypes (asterisks) lie on the ‘corners’ of the data space. In this example, the data points are distributed along lines that can be modelled by ICA axes, and the archetypes are located close to these axes. However, AA can also be applied to datasets that are not distributed along axes.

In this study, we statistically learn archetypes of VF loss by using our VF samples as input. Compared with patterns on ICA axes, we expect VF archetypes to be closer to real VFs and therefore easier to interpret. Compared with cluster analysis VF prototypes, we expect VF archetypes to better emphasize distinctive and representative features of glaucomatous VF loss.

## 2. Material and methods

This retrospective study was approved by the institutional review board of Massachusetts Eye and Ear Infirmary.

### 2.1. Visual fields

The sensitivity thresholds of the protocol Swedish interactive thresholding algorithm [22] Standard 24-2 of glaucoma patients or suspects from a large clinical glaucoma practice were transferred from three Humphrey field analysers (HFA-II, Carl Zeiss Meditec AG, Jena, Germany) to a computer using Peridata (Peridata Software GmbH, Hürth, Germany).

From an initial dataset of 31 591 VF measurements from 8077 patients, measured between October 2006 and April 2012, all VFs with a false-negative rate less than 20% and a fixation loss rate greater than 33% [23] were defined unreliable and excluded. Of the remaining 24 412 reliable measurements of 7300 patients (age quartiles in years: 25%: 51, 50%: 63, 75%: 74), only the most recent measurement of each of the eyes of each patient was included to avoid possible distortions owing to multiple measurements of the same patient. If there were only reliable measurements for one of the eyes of a patient, only this eye was included.

The remaining 13 231 VF measurements were used in this study. For each of these measurements, TD values, that is, the difference of measured sensitivity and the respective population norm at each VF location, were extracted via Peridata. Note that the Peridata norms are similar, but not identical to the STATPAC norms of the HFA. All eyes were converted to ‘right eye format’.

The TDs of the 13 231 measurements are provided as the electronic supplementary material. Of the 54 VF testing locations from the 24-2 protocol, the two locations closest to the blind spot were excluded, as for these locations no TD norms exist.

TD probability plots have been generated using the Peridata norms. While TD values are based on age-specific norms and implicitly contain the age of the patient, the TD distributions differ over age. As the TD probability values are quantiles of these distributions, TD probability plots are specific to the age of each respective patient. Archetypes are purely based on TDs and hence do not explicitly contain age information. In order to generate TD probability plots for the archetypes nevertheless, we used the TD probability values for the age of 63, which was the median age of our dataset.

Mean deviations (MDs) and pattern standard deviations (PSDs) were calculated according to Heijl *et al*. [24] as a sum of TDs weighted by location-specific variance of normal field measurements (MD) and the weighted standard deviation of the TDs (PSD), respectively. The normative variances were taken from [25].

### 2.2. Statistical learning and modelling

All statistical learning was performed using *R* [26] and the *Archetypes* R library [27]. The optimal number of archetypes *m* was determined by 10-fold cross-validation: the data samples were randomly divided into 10 partitions, where each of the 10 subsets was used once as the testing set, whereas the remaining nine subsets were used as the training set. The residual sums of squares (RSS), normalized by sample size, were calculated as a function of *m* for *m* = 2, … , 25, with a cumulative computation time of about one week on an Intel i7 CPU (3.20 GHz). This function was expected to be U-shaped, with decreasing RSS values for small *m*, but increasing reconstruction errors for large *m* owing to overfitting. Model selection was performed by the Bayesian information criterion (BIC). See the electronic supplementary material, methods for further details.

All archetypal analyses were calculated with 100 different random starting values. The model with the respectively lowest RSS was chosen. The relative weights of the archetypes given in some plots are calculated as follows: over all VFs in the dataset, the coefficients of each archetype were summed. These sum values were normalized to sum to 100%.

### 2.3. Retinal nerve fibre trajectories

RNF trajectories were calculated based on the model by Jansonius *et al.* [28,29]. We implemented the model in R. For matching VF locations to RNF trajectories, we used the means and 95% central ranges given in fig. 7 of Jansonius *et al*. [29]. In this work, no parametric distribution models are given for the probabilities of RNF trajectories at a specific VF location to project to a specific angle at the optic nerve head. However, the means and 95% ranges indicate substantially asymmetric distributions. We model these distributions parametrically by *skew normal distributions* [30] with the following density function2.1with a location parameter *ξ*, a scale parameter *ω* and a shape parameter *α*. The mean of the skew normal distribution is given by2.2Let *Φ*(*x*, *ξ*, *ω*, *α*) be the cumulative density function of (2.1) and *Φ*^{−1} its inverse. We obtained *ξ*, *ω*, and *α* by setting the lower bound of the 95% range to *Φ*^{−1} (0.025), the upper bound to *Φ*^{−1} (0.975), and the mean to (2.2) and solved the respective system of equations numerically for every VF location.

These location-specific estimates of the distribution of angles were then used to calculate a *weighted mixture* distribution of all VF locations for every archetype, where the mixture weights were chosen according to the TDs at each location. The TDs are given in dB, a logarithmic scale. It has been noted [31] that for structure–function relations the loss in axon thickness is linearly related to the loss in sensitivity on a linear, but not a log scale. Therefore, we ‘anti-logged’ the dB values for calculating the weights. In the nerve fibre trajectory figures, the resulting mixture distributions are shown by colouring the trajectories with colours between white (zero) and red (maximum of each mixture distribution).

## 3. Results and discussion

### 3.1. Descriptive statistics of the visual field measurements

Figure 4 shows MDs versus PSDs (see §2.1) of our VFs, together with the MD quartiles. Both MD and PSD are known to correlate with glaucomatous VF loss. The middle 50% of the MD distribution, illustrated by the area between the vertical dashed lines, is located between −6.5 dB and −1.1 dB, which corresponds to mild-to-moderate VF defects. This implies that a large number of patients monitored at the glaucoma clinic were glaucoma suspects or glaucoma patients without substantial VF loss.

### 3.2. Archetypes of glaucomatous visual field loss

Figure 5 shows the cross-validation results for numbers of archetypes *m* with 2 ≤ *m* ≤ 25. As expected, the median of the reconstruction error (black line) decreases for small values of *m*. The global minimum of the median is at *m* = 24. For *m* values from 25 on, we expect an increase of the reconstruction error owing to overfitting. As illustrated by the plus symbols, the smallest *m* with a positive BIC difference was *m* = 17, which implies an optimal number of 17 archetypes.

The optimal solution is shown in figure 6. The archetypes are sorted by their relative weights in the dataset (see §2.2). Note that, in contrast to most qualitative VF classification schemes, AA does not partition the data space, so that the relative weights are not the relative numbers of VFs in the respective class. Instead, each VF measurement can be modelled as a weighted sum of *all* archetypes, and archetypes may be correlated.

By far, the largest relative weight (41%) is assigned to archetype #1 with a slightly supernormal MD of 2.5 dB. This archetype represents normal VFs. The mean MD of the normal population is by definition zero. As described in the Introduction, AA places prototypes at the convex hull of the VF data space, which ‘exaggerates’ distinctive features. The distinctive feature of archetype #1 is its absence of VF impairments. This feature is exaggerated first by the slightly supernormal MD and second by slightly higher TD values in the periphery than in the centre, whereas glaucomatous VFs tend to be more impaired peripherally than centrally. The other 16 archetypes represent special types of impairments which will be discussed in §3.4.

Figure 7 illustrates how AA is used to decompose a VF of a 42-year-old glaucoma patient into the 17 archetypes. In this example, four of the archetypes have substantial coefficients, particularly archetype #10 contributes strongly (85.8%). This decomposition transforms the 52-dimensional VF measurement into the 17-dimensional archetypal space. Figure 7*c* shows how the original VF can be reconstructed from that 17-dimensional space by calculating a weighted sum of the archetypes and their respective coefficients. The R source code to decompose a VF measurement into the 17 archetypes is provided as the electronic supplementary material.

Figure 8 illustrates archetype decompositions of seven arbitrarily chosen glaucoma VFs which were not part of the dataset that was used for learning the archetypes. Only archetypes with coefficients greater than or equal to 5% are shown. In addition, MD and PSD (see §2.1) are specified, which are general characteristics of VFs and known to be related to glaucomatous vision loss [32]. While PSD represents localized sensitivity deviations typically observed in glaucoma patients, MD measures global loss of sensitivity, which may or may not be related to glaucoma. The coefficient of archetype #1 (the normal archetype), however, in relation to the sum of the coefficients of all other archetypes, can be regarded as a specific measure of how typical or common a VF is with respect to glaucomatous functional loss. The comparison of the archetype decompositions in figure 8*d,g*, for example, illustrates this relation to glaucoma in contrast to MD: measurement D is a predominantly upper hemifield defect which extends along the horizontal midline, which is more typical for glaucoma than the localized defect in the lower hemifield in G. Indeed, the coefficient of archetype #1 for the latter measurement is substantially higher (37.9%) than the corresponding coefficient in D (11.3%), whereas the MD values are nearly identical (−3.7 and −3.5 dB, respectively).

The sum of the coefficients of typical glaucomatous archetypes always exceeds the coefficient of archetype #1 in the glaucomatous VF loss examples in figure 8, so that our archetypal decomposition has diagnostic potential. However, the main benefit of this decomposition is not so much the distinction of glaucomatous from non-glaucomatous VFs but instead a functional signature of a VF reflected by the archetype coefficients. This functional signature might, for instance, correspond to structural signatures or specific gene defects or can be monitored for progression.

### 3.3. Comparison with ocular hypertension treatment study

We found that all the 16 abnormal archetypes matched definitions of OHTS classes [7]. Figures 9 and 10 show which archetypes match the OHTS definitions of nerve fibre abnormalities and non-nerve fibre abnormalities, respectively.

The close match between the archetypes and the OHTS classes is particularly noteworthy as the OHTS classes have been derived without statistical learning but are instead based on clinical experience and known structure–function relationships related to the optic nerve fibre layer. In contrast, our archetypes were generated by purely mathematical procedures that are agnostic to the RNF layer and its relationship to functional loss.

The comparison with the archetypes reveals essential differences to the OHTS classification, as the verbal, descriptive definitions of these classes are, in contrast to the intention of the authors, not mutually exclusive, not even between the two main categories nerve fibre related versus non-nerve fibre related. The visual inspection of the colour plot of archetype #13 (the archetype indices refer to figure 6), for instance, suggests deep lower hemifield glaucomatous VF loss together with weaker field loss in the upper hemifield, which is compatible with the OHTS category ‘altitudinal’ (figure 9, top). However, the MD of this archetype is −20.9 dB, and its TD probability plot (figure 9, top, right) illustrates ‘severe widespread VF loss (MD ≤ −20 dB)’, which fully matches the OTHS definition of ‘total loss’, which, in turn, is categorized as a ‘non-nerve fibre abnormality’. In addition, the TD probability plot of archetype #6 strongly suggests the total loss class (figure 10, top), but our colour plot shows a ring-shaped, arcuate pattern of particularly severe field loss, so that this archetype represents most likely an advanced glaucomatous loss, which makes it a nerve fibre abnormality rather than a non-nerve fibre abnormality.

### 3.4. Characteristics of the archetypes

In accordance with the OHTS classification, we find archetypes which suggest relationships to RNF trajectories and archetypes for which there is no indication of such a relationship. Among the latter, archetypes #12 and #15 (all indices refer to figure 6) indicate hemianopia, an impairment typically related to strokes. Archetype #7 might be related to cataract or other non-nerve fibre abnormalities like non-specific photoreceptor or choroidal disease. Archetype #17 does not indicate any obvious relation to an ophthalmic disorder. Its relative weight in the dataset is only 0.56%, and a closer inspection revealed that for only one single, possibly artefactual, VF measurement the coefficient for archetype #17 is higher than for all other archetypes. Therefore, this archetype should be considered an overfit to our specific dataset.

The remaining 12 archetypes indicate relations to RNF trajectories. Archetypes #2, #3, #8 and #14 predominantly affect the upper visual hemifield, archetypes #5, #9, #10, #13 and #16 the lower hemifield, archetypes #6 and #11 are annular (ring-shaped) and archetype #4 is restricted to the area temporal to the blind spot. While we consider archetype #4 as nerve fibre related here, as it is known that localized temporal defects can be glaucomatous, it needs to be noted that the occurrence of localized temporal VF loss in relation to nerve fibre defects is uncommon (for instance, it was only found in 3.3% of the measurements in OHTS [7]). Besides, the 24-2 protocol samples only a small region temporal to blind spot outside which the shape of the defect is unclear. Therefore, VFs with dominant coefficients of archetype #4 are not necessarily related to nerve fibre defects.

Figures 11–14 illustrate the nerve-fibre-related categories further by comparing each archetype with the locations of nerve fibre impairments predicted by the model by Jansonius and co-workers (see §2.3) and adds two illustrative VF measurements which have coefficients more than 75% for the respective archetype.

#### 3.4.1. Retinal nerve fibre trajectories

The RNF trajectory plots show mixtures of probability distributions weighted according to the TDs. Given the VF of the respective archetype, and given that the sensitivity loss at each location is indeed nerve fibre related, the model estimates which of the nerve fibres might be impaired and to what degree. Note that the model can infer only about those nerve fibre trajectories that are related to those VF locations that have been measured. That means, additional nerve fibres which project to locations where no sensitivities have been measured might also be damaged. In addition, the model necessarily neglects individual retinal parameters. Therefore, the illustrations should be regarded as a rough estimate of structural damage and its location.

Apart from that, the model has two singularities, namely at ±60° around the optic nerve head, where the nerve fibre trajectories are too unspecific to be covered. In some of the plots, these regions appear as a gap, which has been annotated accordingly in the respective plot.

VF locations get increasingly non-specific with respect to nerve fibre trajectories the greater their distance from the optic nerve head, and superior structural damage, which are related to lower visual hemifield locations, are generally more unspecific than inferior damage. Furthermore, although our VF measurements focus on locations outside the macula, our plots verify two glaucoma-related *macular findings* [33]: the nerve fibre trajectories of the upper central archetype (figure 12), which resemble a pattern previously termed ‘arrow shaft’ by Goldbaum *et al*. [19], is restricted to the so-called *macular vulnerability zone* (see fig. 9 in [33]). This zone is related to early central upper hemifield damage which is typically hard to detect by the coarse spatial resolution of the 24-2 protocol. In addition, all our trajectory plots, including the heavily impaired total loss archetype (figure 13), confirm the so-called central isle within the upper macula which is functionally reflected by a region of relative sensitivity preservation to the lower right of fixation, which has been shown with the macula-specific 10–2 measurement protocol (fig. 12*c* in [33]).

Hence, AA provides a statistical verification of findings for which our 24-2 measurements were inappropriately coarse.

#### 3.4.2. Upper versus lower hemifield archetypes

Figures 11 and 12 demonstrate that the specific upper hemifield archetypes have matching analogues in the lower hemifield. Particularly, upper and lower altitudinal, central and nasal ATs roughly match. Upper peripheral AT, which is relatively centralized between nasal and temporal hemifield, has a temporal tendency and therefore matches roughly lower peripheral temporal AT. Generally, upper hemifield archetypes are more localized or focused than their lower hemifield counterparts. In the lower, but not in the upper hemifield, AA prefers a separation of peripheral nasal and peripheral temporal patterns, which is in accordance with the OHTS definition of partial arcuate VF loss as a defect contiguous with either the blind spot or the nasal meridian.

The assumed symmetry of upper and lower hemifield is a necessary condition for several diagnostic tools, particularly the *glaucoma hemifield test* (GHT, see [34]), which divides the hemifields into symmetric sectors and compares summary indices of these sectors between the hemifields. While our findings show symmetries on a very general level by matching upper and lower archetypes, they challenge the strict symmetry assumption over exactly matching sectors.

Although the two central archetypes roughly correspond to each other, they are not symmetric, particularly not around fixation. This is compatible with nerve fibre trajectory asymmetries in the macula [33]: the VF locations closest to fixation of upper central might reflect the macular vulnerability zone, whereas the corresponding locations in lower central, and to some extent also in the other lower archetypes, might reflect the central isle of preservation, which is located in the lower temporal part of the VF close to fixation (see fig. 12 in [33]).

#### 3.4.3. Nasal and temporal archetypes

Archetype #5, which resembled, but did not fully match, the ‘nasal step’ definition of OHTS, was regarded by us as a lower hemifield correspondent to the clearly nasal pattern of archetype #3. Note that these predominantly nasal archetypes do not respect the horizontal meridian. A reason for this might be slight anatomical variants among patients which, particularly for nasal locations, lead to mappings in the opposite hemifield, as discussed, for example, in [35], but it might also signify simultaneous vision loss in both hemifields.

The *temporal* archetype does not respect the horizontal midline, with is in accordance with the nerve fibre trajectory model by Jansonius *et al*. [29] and compatible with the anatomic fact that the nerve fibre layer does not have a horizontal raphe in the nasal field. The OHTS definition of the temporal wedge class does not mention horizontal midline effects, although the example given in the study clearly respects the midline (fig. 1 in [7]).

#### 3.4.4. Annular archetypes

While OHTS classifies the total loss pattern as non-nerve fibre abnormality, our colour plot, in contrast to the TD probability plot, clearly indicates the ‘central isle’ of relative sensitivity preservation discussed above, which is why we combine total loss together with the tunnel archetype in a category of annular patterns.

The tunnel archetype resembles a pattern which is suggested to be artefactual in OHTS (peripheral rim). Note, however, that the OHTS peripheral rim (fig. 3 in [7]) is restricted to 30-2 VF locations outside the locations of our 24-2 measurements.

## Data accessibility

The dataset of the TD measurements, the numerical values of the 17 archetypes and the R codes to compute the nerve fibre model as well as to perform decompositions of VF measurements into archetypes are uploaded as the electronic supplementary material.

## Funding statement

T.E. and P.B. were supported by National Institutes of Health Grant No. R01 EY018664. L.R.P. was supported by the Harvard Glaucoma Center of Excellence, a Harvard Medical School Scholar award and the Research to Prevent Blindness Foundation in New York. L.Q.S. was supported by the American Glaucoma Society Young Clinician Scientist Award, 2012. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

- Received October 10, 2014.
- Accepted November 13, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.