Predicting Visual Function from the Measurements of Retinal Nerve Fiber Layer Structure

Abstract

Purpose.: To develop and validate a method of predicting visual function from retinal nerve fiber layer (RNFL) structure in glaucoma.

Methods.: RNFL thickness (RNFLT) measurements from scanning laser polarimetry (SLP) and visual field (VF) sensitivity from standard automated perimetry were made available for 535 eyes from three centers. In a training dataset, structure–function relationships were characterized by using linear regression and a type of neural network: radial basis function customized under a Bayesian framework (BRBF). These two models were used in a test dataset to (1) predict sensitivity at individual VF locations from RNFLT measurements and (2) predict the spatial relationship between VF locations and positions at a peripapillary RNFLT measurement annulus. Predicted spatial relationships were compared with a published anatomic structure–function map.

Results.: Compared with linear regression, BRBF yielded a nearly twofold improvement (P < 0.001; paired t-test) in performance of predicting VF sensitivity in the test dataset (mean absolute prediction error of 2.9 dB [SD 3.7] versus 4.9 dB [SD 4.0]). The predicted spatial structure–function relationship showed better agreement (P < 0.001; paired t-test) with anatomic prior knowledge when the BRBF was compared with the linear regression (median absolute angular difference of 15° vs. 62°).

Conclusions.: The BRBF generates clinically useful relationships that relate topographical maps of RNFL measurement to VF locations and allows the VF sensitivity to be predicted from structural measurements. This method may allow clinicians to evaluate structural and functional measures in the same domain. It could also be generalized to use other structural measures.

Assessing the way in which structural and functional measures in glaucoma interact is clinically important. Visual loss is assumed to follow from, and correlate to, structural loss caused by the disease process. It would be clinically useful to know the magnitude and location of structural loss that will result in visually important functional loss. However, current clinical devices for measuring structural and functional deficits are far from accurate and have imperfect precision. Structural measures of optic nerve head (ONH) topography or retinal nerve fiber layer thickness (RNFLT) can be made using optical imaging techniques, but these provide only surrogate measures of the biological variable of real interest—namely, the retinal ganglion cell (RGC) count and function. At the same time, standard automated perimetry (SAP), the clinical cornerstone of functional testing in glaucoma, is subject to considerable measurement variability and is also a poor surrogate for RGC count and function. ¹ Despite their limitations, these techniques are currently central to the diagnosis and management of glaucoma. It would, therefore, be beneficial if structure and function measurements were directly linked in some way, allowing clinicians to corroborate decisions by considering the measurements in tandem. 

Several studies have been conducted in an attempt to quantify the structure–function relationship using clinical measurements. These typically proceed by taking one summary value to represent function (for example, mean deviation [MD] of the visual field from SAP) and one number to represent the structural data (for example, average neuroretinal rim area or mean RNFLT), then assessing the curvilinear (e.g., log-linear) or monotonic association between the two variables via R ², Pearson, or Spearman coefficients. ^{2 –5} This approach has two major flaws: The use of summary data loses spatial information and may reduce power, and these association measures and regression models assume a particular shape of the relationship. Furthermore, these analyses fail to take account of spatial associations in the data, an integral attribute of glaucomatous loss. 

Although some investigators have attempted to overcome these limitations by analyzing smaller regions from the structural measures (for example, sectoral values of rim area) and groups or individual points in the visual field (VF), ^{6 –12} we propose that much more could be gained by analyzing the data in its high-resolution form. For example, in scanning laser polarimetry (SLP), RNFLT estimates are yielded over an image space of several thousand pixels. These are high-dimensional data and any method linking structural measures to the 50 or so individual points in the VF ideally should take this into account. Moreover, individual points from both structure (pixel or sector values) and function (areas of VF or individual locations) are more likely to interact as groups rather than single independent measurements. 

In recent years, a special class of artificial neural networks (ANNs) using radial basis functions (RBFs) have received considerable attention. ^{13 –16} This class of networks differs from the ANNs that use a multilayer perceptron (MLP) approach, in that the nonlinearity of the model is embedded in a hidden layer of the network. This hidden layer consists of basis functions. The key element of the RBF is that predictions can be made in a multidimensional space that consists of mini-distributions of possible predictions, analogous to a kernel or window when interpolating data. RBF networks are particularly useful when mapping two- or three-dimensional images where interpolation is a prerequisite, with some image features being preserved while others are not necessarily mapped exactly. In this way, they seem most suitable to the problem of mapping the individual points from the structural and functional measures. Moreover, the RBF method developed for our purposes is formulated under a Bayesian probabilistic framework to tackle the considerable variability in measurements as well as to form an automatic learning process. This customization leads to what we denote a Bayesian RBF (BRBF) network. 

The purpose of this study is to develop a general statistical method for linking clinically used measures of structure and function. We demonstrate the use of the method on measurements of the RNFLT derived from SLP with variable corneal compensation (GDxVCC; Carl Zeiss Meditec AG, Jena, Germany) and individual VF locations from SAP with the purpose of providing predictions of one measure from another and to generate a topographic map of the spatial relationship between the imaged structure data and individual points in the VF. It is imperative that any new method should be developed and tested on more than one data set. For this work, we used three large data sets collected at three independent clinical centers. The method, which aimed to make predictions, is compared to the correlation approach for mapping structural and functional measures typified by Gardiner et al. ⁹ In addition, we illustrate how the method can be extended to high-resolution measurements. 

Methods

Subjects

The study sample was derived from three independently acquired populations from Moorfields Eye Hospital, London (MEH), Rotterdam Eye Hospital, The Netherlands (REH) and the Blue Mountains Eye Study, Australia (BMES). 

MEH Data

Thirty-four healthy subjects (34 eyes), 43 glaucoma patients (43 eyes), and 30 ocular hypertension subjects (30 eyes) were enrolled. Inclusion criteria for the healthy subjects were a normal VF, intraocular pressure (IOP) < 21 mm Hg, no previous history of ocular disease, and no family history of glaucoma in a first-degree relative. For the glaucomatous patients, inclusion criteria were history of raised IOP (>21 mm Hg), reproducible VF defects, and absence of other disorders that may cause VF loss. The VF defects were defined as a defect of two or more contiguous points in the pattern deviation probability maps with P < 1% loss or greater, three or more contiguous points with P < 5% loss or greater, or a 10 dB difference across the nasal horizontal midline at two or more adjacent points in the total deviation plot. Ocular hypertension subjects were recruited from clinic on the basis of having a raised IOP >21 mm Hg in two consecutive visits and having at least two normal VFs. Optic disc appearance was not used to categorize subjects. However, the optic disc was evaluated to exclude anatomic abnormalities such as coloboma or drusen. For all participants, one eye was randomly selected for study if both were eligible. All subjects had visual acuity of 20/40 or better, with ametropia <7 diopters, and had no other significant ocular abnormality or concomitant ophthalmic condition. 

REH Data

Forty six healthy subjects (46 eyes) and 76 glaucoma patients (76 eyes) took part in a study that has been described in detail in Reus and Lemij. ¹⁷ In short, healthy subjects had normal VF, optic discs without structural abnormalities, IOP < 21 mm Hg, no previous history of ocular disease, and no family history of glaucoma. The patients had a glaucomatous appearance of the optic disc and a corresponding nerve fiber bundle VF defect, as described by Keltner et al. ¹⁸ with SAP. All subjects had a visual acuity of 20/40 or better and had no other significant ocular abnormality. For all participants, one eye was randomly selected for the study. 

BMES Data

VFs and images were available for 1540 subjects from a large population-based study of visual impairment, common eye diseases, and other health conditions from an elderly community in Australia. A description of the BMES is given elsewhere. ¹⁹ Two hundred thirty healthy subjects (230 eyes) and 76 patients (76 eyes) diagnosed with glaucomatous optic neuropathy were selected from this population under strict measurement quality criteria (described later). Only one randomly selected eye per subject was used throughout. The criteria used for defining glaucomatous VF loss was an abnormal result of the glaucoma hemifield test (GHT) plus one or more of the following VF defect classifications (top row of test points in the 24-2 pattern were excluded to reduce the effect of lid artifact), which could not be explained on the basis of nonglaucomatous ocular, or neurologic, causes: (1) at least four contiguous points on the pattern deviation plot depressed at the P < 0.5% level; (2) at least two horizontal points in nasal step locations with a pattern deviation plot depressed at the P < 0.5% level; and (3) advanced glaucomatous field loss (hemispheric or severe generalized field loss, with residual temporal or central islands). Glaucoma was diagnosed if glaucomatous VF loss was present, combined with matching neural rim loss at the optic disc, and gonioscopy showed no evidence of angle closure, rubeosis, or secondary glaucoma, other than pseudoexfoliation syndrome. 

All three datasets were collected in accordance with the tenets of the Declaration of Helsinki from studies that had local research ethics committee approval, with all participants giving informed consent. Data were anonymized and transferred to a secure database held at City University London. 

Measurements and Data

This study made use of measurements of the RNFLT yielded from scanning laser polarimetry (SLP: GDxVCC; Carl Zeiss Meditec AG) and individual VF locations from SAP. The principles of SLP have been described by Greenfield. ²⁰ This instrumentation estimates the thickness of the peripapillary RNFLT by measuring the summed retardation in a polarized scanning laser beam reflected from the fundus. Retardation measurements at various points around the optic disc are used to construct a thickness map of the retinal nerve fiber layer (RNFL) in micrometers. SLP with variable corneal compensation (VCC) has been shown to improve the estimates of RNFLT compared with earlier versions of the technology. ²¹ All SLP images from all three centers had quality scores ≥8 and typical scan scores ≥80. Single SLP images were available from all the subjects from all the centers. The 64-sector thickness profile on the peripapillary annulus (provided by the GDx software) for each subject and the raw images (.mif files) were transferred to a secure database. 

The VF measurements were acquired in all cases with the Humphrey Field Analyzer (HFA) II (Carl Zeiss Meditec, Dublin, CA), with either full-threshold testing or the Swedish Interactive Threshold Algorithm (SITA) standard test program in the standard 24-2 test pattern. For the MEH and REH data, all VFs were considered reproducible as well as reliable. The VFs from MEH data were all tested with the SITA standard program. VF reliability indices applied were fixation losses ≤15% and false-positive and false-negative response rates ≤25%. In the REH data, 29 healthy subjects and 73 glaucomatous subjects were tested with the full-threshold program and the others with the SITA standard program. Reliability indices applied were fixation losses ≤25% and false-positive and false-negative response rates ≤20% for the full-threshold test and ≤7% for the SITA standard test. Higher false-negative response rates were accepted in eyes with advanced field loss (up to 33%). All VFs in BMES were tested with the full-threshold program. Stricter reliability criteria were applied: fixation loss and false-positive and false-negative responses ≤15%. In all eyes, the two VF test points above and below the blind spot were excluded, and the remaining 52 raw sensitivity values (in decibels) were transferred to a secure database for analysis. 

In this study, the MEH and REH datasets (229 eyes in total) were used for development of the models, and the BMES dataset (306 eyes in total) was used for independent testing to demonstrate the generalization of the method. The image and VF quality criteria were used only for the purpose of selecting reliable measurements and were not used in the modeling. Moreover, the proposed model does not need the initial labeling of the subject data as being from glaucomatous or healthy eyes, and under the Bayesian framework, no manual input parameters are required. 

Statistical Models

What follows is a description of the principal methods and models developed. Technical and mathematical details can be found elsewhere ²² and are reviewed in the 1. The BRBF model was developed and written in commercial software (MATLAB, ver. 7.4.0 R2007a; The MathWorks, Inc., Natick, MA). An executable version of this code is freely available from the authors. 

We sought to generate a model that will map an individual structural measurement (RNFLT values from a software generated peripapillary profile or an individual RNFLT value at a pixel in the image) to a sensitivity value (in decibels) at an individual VF location. Such a model will allow the prediction of a VF sensitivity value from RNFLT values and vice versa. In the remainder of this section, we compare the BRBF model to the linear model that has been used widely to assess the association between structural and functional measurements. 

Starting with a linear model, consider the case in which we attempt to predict an individual VF sensitivity value denoted ŷ _d (where d is one of the 52 locations in a 24-2 HFA VF) from a series of RNFLT values denoted x _i for i from 1 to m. This can be expressed by the following oversimplified but illustrative equation:  

 

where c_d is a constant offset. The symbol on ŷ _d indicates that it is a prediction rather than the real measured value denoted by y _d . In this example, the equation has 64 peripapillary thickness profile values (m = 64), each with its own coefficient that quantifies the contribution of each x value to the prediction. Thus, each y value can be predicted by a linear combination of x values. With some actual data we can find some real numbers for the w terms by least-squares regression. This calculation attempts to fit an equation that minimizes the difference between the predicted and measured values. It yields an individual equation for each y value that can be enumerated across all the points, to predict a complete VF from a given vector of x values. This classic multiple linear regression can be adapted to select only those x values that are statistically significant for the prediction of y values, by using techniques such as stepwise multiple regression ²³ with the forward-selection scheme. In the linear model described by equation 1, the w values (divided by their standard errors) quantify the amount of meaningful contribution made by x values to predict the y values. The largest absolute w term (with respect to the variability in estimating the term) would indicate the x value that affects the y _d value the most, in the sense that a change in this x value results in the largest change in y _d . Similarly, the next largest absolute w term would indicate the second most important term and so on. Equivalently, one could look at the relationship between the y _d value and each x value separately and simply calculate the correlation coefficient based on the raw data or on the ranks of the data (Spearman's ρ) and end up with a similar result. Loosely speaking, this is the approach of Gardiner et al., ⁹ who used neuroretinal rim area estimated from scanning laser ophthalmoscope measurements (Heidelberg Retina Tomograph [HRT]; Heidelberg Engineering, Heidelberg, Germany) as the surrogate structural measurement for glaucomatous damage. Particularly, in the linear model implemented for this study, the VF sensitivity was unlogged from the decibel value and the prediction from this model was converted back to the decibel scale when compared with the measured sensitivity. We will refer to this method as a classic linear model. 

This classic linear model approach makes several restrictive and incorrect assumptions about the data. First it assumes that each x value is independent of all the other x values, ²⁴ whereas, in reality, the x and y values are topographically and physiologically related and may interact as groups. Although one could try to demarcate these groups on the basis of a physiological relationship between the x values or of an anatomic map, it would be preferable for the groups to be learned from the data rather than imposing any relationship from incomplete prior knowledge. Second, this approach assumes that the relationship between y and x is either linear or becomes linear after some transform (typically logarithmic). In reality, the relationship between y and x may be more complex, with the nature of the association probably changing across the measurement range of y. Put simply, at different stages of disease, the apparent relationship between y and x could switch from linear, to noisy, to curvilinear and occasionally could be censored because of the imprecision and range of the measurement. This notion that the association may change at different levels of functional loss has also been asserted in published studies. ^25,26 The third difficulty with the classic linear model is that outlier points exert an overly strong influence and can yield a false association. 

The RBF attempts to model the relationship between y and x without the limiting assumptions associated with the classic linear model just described. As an illustrative example, consider one x value, say x ₁, appearing to be correlated to a y value, say y ₁. This apparent relationship might be explained by the very strong relation of x ₁ to x ₂, which in turn is very strongly associated with y ₁. Thus, the apparent x ₁ and y ₁ relationship may well be much weaker or not even significant. As it will be shown in the derived structure–function relationship (see Figs. 3 4–5), this covariance in the relationship between the y and x values is modeled better with the RBF approach. Furthermore, the central idea of the RBF is the basis functions, each of which performs very much like a dynamic window or kernel that moves across the data, both spatially and at various stages in disease severity, identifying groups of measurements that appear to behave in a similar pattern. The non-normalized Gaussian basis function used in this study has an activation field that has a center—that is, a particular input value at which it has a maximal output. The output tails off as the input moves away from this point. In this way, those hidden basis functions that have centers similar to the input x patterns will have stronger activation and will thus contribute more to the prediction of y. On the other hand, those basis functions with weak activation will be isolated and will not affect the prediction. Moreover, the RBF learns the parameters from the data and this is customized for our purpose by using a BRBF. The method makes predictions in multiple dimensions by extending the standard relevance vector machine ²⁷ where the output is restricted to one dimension. The VF sensitivities at different locations are implicitly correlated by sharing the same basis functions and parameters of weights (see α in the 1). In the BRBF model, the original VF sensitivity in decibels was used. 

In general, BRBF is not restricted by the number of inputs in the model or, more precisely, by the dimensionality of the data, so that it can be adapted to use the RNFLT values at each pixel, rather than only the 64-sector peripapillary thickness profile. In this study, examples of the SLP images used for predicting VF are shown in Figure 2. We used an annular region, centered on the ONH. The inner and outer diameters were 2.3 and 4.9 mm, respectively, so that the annulus was 26 pixels wide (compared with the eight-pixel calculation ring in the GDx software from which the 64-sector RNFLT peripapillary profile is computed). In this annulus, there are 16,512 pixels each with a retardation value. We hypothesized that the predicted spatial structure–function relationships would be strengthened by avoiding the data reduction to 64-sector RNFLT values. This hypothesis will be validated by the improved structure–function relationship derived from the BRBF by using the measurement in its higher dimensional form. 

Another challenge in the modeling process involves handling the large dimensionality of the SLP data. If the dimensionality M and the number of data points N of a dataset satisfy M ≥ N, the dimensionality of the dataset can be reduced from M to N − 1, with minimal loss of information using principal component analysis (PCA), ²⁸ because the N data points span at least one linear hyperplane in this M-dimensional space. Using this technique the 16,512 dimensional SLP image vectors were reduced to 228 dimensions for analysis and transformed back to the original SLP vectors for the purpose of visualization and evaluation of results. Using the classic linear model on a reduced SLP image vector with 228 elements and a 24-2 VF with 52 sensitivity values will still result in a prohibitive number of weights (11,856) to be fitted, which will cause significant overfitting to the noise in the data. 

Testing the Model

It is well established that if validation of a modeled relationship between dependent and independent variables is performed on the same dataset as the selection of input variable and weights, then the model estimates will be overly optimistic. ^29,30 It is also well known that in developing a model, the input variables and weights selected may vary across different samples. ³¹ Therefore, the models were developed on the MEH and REH data alone, leaving the BMES data as a test dataset. 

The predictive performance of the classic linear model and the BRBF was evaluated by point-by-point analysis of the predictions of the VF sensitivity in the 306 VFs from the BMES dataset. The predictive performance was summarized by the mean of the absolute prediction errors in the 52 points of the VF. 

Code was written in commercial software (MATLAB ; MathWorks) to display the results and present the output of the linear and BRBF models via a graphic user interface. One graphic output was a structure–function map in a format similar to that described by Gardiner et al. ⁹ The structure–function relationship was defined by the corresponding functional change given a subtle structural change, which is mathematically modeled by the derivative of the BRBF (1). This derivative describes the relationship between y (VF sensitivity) and x (RNFLT) at each VF location. The other output from this analysis consisted of point-by-point predictions of each subject's VF, as represented by the HFA grayscale (which was replicated for this purpose). These outputs were considered for (1) the classic linear model, (2) the BRBF I model (based on the 64 summary RNFLT values output from the GDx software), and (3) for the BRBF II model (based on the reduced 228 dimensional data derived from the PCA on the 16,512 individual pixel retardation values in the broad annulus centered on the ONH). 

Results

A summary of the measurements—HFA mean deviation (MD), HFA pattern standard deviation (PSD), and GDx nerve fiber indicator (NFI)—for each of the datasets is given in Table 1. 

Table 1.

 

A Summary of the Measurements of the Three Datasets (MEH, REH, BMES) Used in This Study

Data Set	Ocular Status	Sex (n)	Age (y)	MD (dB)	PSD (dB)	NFI
MEH	Healthy	19 M, 11 F	40.6 ± 15.6 (21, 75)	1.26 ± 0.76 (0.09, 3.30)	1.47 ± 0.28 (1.12, 2.27)	16 ± 8 (2, 30)
	OHT	13 M, 21 F	60.7 ± 11.1 (21, 75)	1.07 ± 0.82 (0.10, 3.06)	1.35 ± 0.27 (1, 2.1)	20 ± 10 (5, 43)
	Glaucoma	28 M, 15 F	60.8 ± 13.1 (31, 84)	−4.02 ± 2.55 (−12.00, 1.15)	4.93 ± 2.93 (1.49, 12.50)	42 ± 18 (11, 80)
REH	Healthy	23 M, 23 F	60.4 ± 12.1 (23, 77)	0.38 ± 0.99 (−1.55, 2.73)	1.63 ± 0.26 (1.13, 2.30)	21 ± 9 (2, 43)
	Glaucoma	47 M, 29 F	62.2 ± 10.1 (30, 82)	−9.52 ± 8.43 (−30.39, 1.25)	8.35 ± 4.32 (1.99, 15.92)	63 ± 21 (21, 98)
BMES	Healthy	99 M, 131 F	69.3 ± 6.5 (60, 87)	0.21 ± 1.07 (−1.26, 3.03)	1.53 ± 0.29 (1.12, 2.55)	19 ± 10 (2, 49)
	Glaucoma	25 M, 51 F	72.0 ± 6.3 (61, 78)	−7.94 ± 6.55 (−29.67, 1.56)	6.97 ± 3.76 (1.67, 15.56)	65 ± 24 (27, 98)

 The data are expressed as the mean ± SD (min, max). MD, PSD, and GDx NFI results are provided for the three datasets. OHT, ocular hypertension.

The mean (SD) absolute prediction error of VF sensitivities in the 306 eyes from the BMES data were 4.9 dB (4.0 dB) for the classic linear model. In comparison, both BRBF models yielded a nearly twofold improvement (P < 0.001; paired t-test) in performance (mean [SD]): BRBF I 2.9 dB (3.7 dB) and BRBF II 2.8 dB (3.8 dB). In the BRBF I and II models, the training process described in the 1 selected 49 and 73 basis functions, respectively, in the hidden layer. 

Figure 1 summarizes the predictive performance of the classic linear model and the BRBF I model across the range of VF sensitivity measurements. BRBF II demonstrated a predictive performance similar to that of BRBF I. Published test–retest differences (5th and 95th percentile limits) ³² across all locations from two VF tests from each of 49 individuals have been superimposed on Figure 1b. On inspection, these limits are similar to the 90% prediction limits when using the BRBF to predict the VF from a GDx RNFLT measurement of the same individual. Note that predictions at higher sensitivities (>30 dB) tend to be slightly lower than the actual values, whereas at lower sensitivities (<20 dB), the predictions tend to be higher. 

Figure 1.

 

Distributions of the error between the predicted and the measured sensitivity for each VF location in 306 eyes from the BMES data, stratified by VF sensitivity. Each error bar summarizes the predictive performance over a 2-dB range from 0 to >36 dB. Thin vertical lines: 90% prediction limits (5th and 95th percentile of error), each box indicates the interquartile range of the prediction error (25th and 75th percentile error) with the line in the box indicating the median error. The dotted line of unity indicates perfect prediction (no error). The predictive performances of (a) the classic linear model and (b) the BRBF I model are shown. (b, solid lines) Previously published (5th and 95th percentiles) test–retest limits ³² for VF data derived from the point-wise differences between two VFs tested over a short period.

Figure 2 gives some case examples of the predictions. In some cases, the classic linear model overestimates the defect severity of the VF (Fig. 2I, 2III) and in other cases predicts a less damaged VF (IV), when compared with the true measured VF. In Figure 2II the classic linear model matches the overall average sensitivity of the VF but fails to capture the spatial location of this loss. In each case, the BRBF better estimates the true VF, with spatial features of the measured defects generally retained. In Figure 2IV the BRBF model manages to predict the advanced defect severity. More examples are provided as Supplementary Material. 

Figure 2.

 

Model predictions for four cases from the BMES dataset. For each case (I–IV), the top row shows, from left to right, HFA VF grayscales for the measured VF, the VFs predicted from the classic linear regression, and the BRBF with 64-sector RNFLT (BRBF I), respectively. The predicted VF from BRBF II with 16,512-pixel annulus from the GDxVCC SLP images was in general similar to the one from BRBF I. The row of graphics (below) shows the SLP image annulus and the corresponding 64-sector RNFLT (black line) used to predict the VFs. The colored lines in the 64-sector RNFLT plot indicate the top five RBFs with the strongest activation (those contributing the most to the prediction). More examples are provided as Supplementary Material.

Figures 3, 4, and 5 provide topographical maps of the spatial relationship between the RNFLT positions on the peripapillary annulus and VF sensitivity at all the points in the VF, by means of classic linear regression, BRBF I, and BRBF II, respectively. The agreement between the derived structure–function relationship and the anatomic benchmark was also summarized by calculating the absolute angular difference between the direction with the strongest derived relationship and the anatomic benchmark. This difference was calculated as the median of the absolute angular difference across all locations in a spatial relationship map, such as Figures 3, 4, and 5. Both visual inspection on the spatial relationship and the quantified agreement with anatomic prior knowledge showed that the structure–function relationship derived by the classic linear regression model (Fig. 3) has little concordance (median absolute angular difference of 62°) with the anatomic benchmark. ³³ This correspondence improves (P < 0.001; paired t-test) with the BRBF I model (median absolute angular difference of 15°) in general, and improves further (P < 0.001; paired t-test) still with the BRBF II model (median absolute angular difference of 12°), especially at the points around central vision and blind spot. This result is consistent with the BRBF models learning to encode the structure–function relationship during the training. 

Figure 3.

 

A topographical map describing the relationship between the 64-sector RNFLT profile and individual VF locations, as described by classic linear regression. The map is in the shape of an HFA 24-2 VF test for a right eye. Each location is represented by a circular graph made up of 64 bars representing the correlation value between the VF location and the RNFLT at the corresponding angle. In this instance, the correlations are derived from scaled values of the weights given in equation 1. The length of the bar indicates the magnitude of the correlation: red bars, positive correlation; blue bars, negative correlation. The green bar with an asterisk indicates the location of expected strongest correlation on the basis of the anatomically derived map from Garway-Heath et al. ³³

Figure 4.

 

A topographical map describing the relationship between the RNFLT profile and individual VF points as described by BRBF I using the 64-sector RNFLT profile. The composition of the graph is the same as Figure 3, with the green bar with an asterisk indicating the location of the expected strongest correlation on the basis of the anatomically derived map from Garway-Heath et al. ³³ In this instance, the correlations are derived from scaled values in equation 13.

Figure 5.

 

A topographical map describing the relationship between the RNFLT profile and individual VF points as described by BRBF II, using the 16,512 pixel retardation values from the SLP image. The composition of the graph is the same as Figures 3 and 4 with the green bar with an asterisk indicates the location of expected strongest correlation on the basis of the anatomically derived map from Garway-Heath et al. ³³ In this instance, the correlation values are derived from scaled values in equation 13 and are summarized in 64 sectors for the purpose of comparison.

Discussion

The application of ANNs to both functional and structural measurements in glaucoma is not a new idea. ^{34 –40} Most of these applications have tended to use a conventional MLP ANN. Technically, there are several advantages of the RBF over MLP ANNs. For example, with the latter, there is an input, and the distributed pattern “lights up” all hidden units, to contribute to the prediction of the output, which makes them combine and interfere with each other. This typically yields a highly nonlinear training process with mathematical difficulties that result in a slow convergence of a training procedure. ¹⁶ Moreover, the complexity of hidden unit patterns causes difficulty when interpreting the result, because the connection among units and the hidden-unit output do not have any physical or realistic meaning: they are simply numbers that can produce a correct output. Hence, MLP ANNs are less suitable for the mapping of points in different measurement spaces, which requires a detailed understanding of the hidden layer output and other manipulation (e.g., derivatives) within the network. In this study, although an MLP provided a prediction accuracy comparable to the BRBF (data not shown), the spatial aspects of the structure–function relationship were poorly predicted. This deficiency probably resulted from complex interactions among the large number of weights in the MLP model. In contrast, RBF substitutes hidden layers in MLP with a set of basis functions, which leaves the solution of weights in a linear space. The Gaussian basis function forms a local representation in the hidden unit, each of which can be understood as a representative of similar input patterns (Fig. 2). With a given RNFLT input, only a few representatives will be activated and contribute to the VF prediction (Fig. 2). Furthermore, the RBF handled within a Bayesian probabilistic framework, ²⁷ as developed for our purposes, is unlike most other ANN approaches because it is independent of any subjective input parameter and thus requires no model validation on a test dataset, provided that the training dataset is representative and sufficiently large to enable modeling of the many different states of the VF. However, for this study, the trained model is still tested on a separate, independent dataset to illustrate the generalization of the model performance. One technical limitation of the current BRBF model, despite its good performance, is that it assumes that the variability in the VF measurements is largely Gaussian, which is not optimal, given that it is often skewed and heavily tailed. ³²  

VF locations have been related to sectors of the ONH using an anatomic map ³³ derived by overlaying an appropriately scaled VF grid on RNFL photographs and tracing nerve fiber bundles or defects from various VF locations to the ONH margin. Gardiner et al. ⁹ produced a topographical map of the relationship between sectors of the ONH and locations in the VF by considering the linear correlation between Heidelberg retina tomography neuroretinal rim area and VF sensitivity at each point in the VF. The classic linear model is akin to that derived by Gardiner et al. except that the SLP measurements from GDx imaging are used. The maps derived from the BRBF models (Figs. 4, 5) indicated a closer concordance with the anatomically derived landmark than did the map derived from the classic linear model (Fig. 3). In addition, the predictions on point-wise VF sensitivities in the validation dataset were, on average, better with the BRBF method than with the classic linear model. A plausible explanation for the improvements afforded by the BRBF technique is that it models the spatial and quantitative structure–function relationship more precisely. The technical and statistical advantages of the BRBF model over the classic linear model support this notion. For example, the BRBF considers that the VF points, and indeed RNFLT values acquired from different discreet areas, interact as groups rather than as independent measurements. The BRBF also accommodates the covariance and nonindependence of the measurements. It is less affected by outlying observations and makes no assumption about the linearity of relationships. Of course, one difficulty in generating any type of map, driven by data or anatomic observation, is the restriction in the sampling of VF points. These are, for example, probably not optimally placed for estimating RGC density nearer the fovea. 

The range and distribution of differences between the measured VF sensitivity values and those predicted from the RNFLT by the BRBF, at different levels of sensitivity, is shown in Figure 1b. This profile is remarkably similar to published limits for test–retest variability when two VFs are measured within a short time. ³² This similarity suggests that, on average, a VF predicted by the BRBF from RNFLT values has measurement noise equivalent to that found in a newly measured field. This finding is not as exciting as it may first appear, because it is well established that the measurement noise in VFs is already very high, prohibiting straightforward clinical diagnosis of glaucomatous defects and monitoring progression. Nevertheless, this finding illustrates that the range and scale of the average predictive performance of the BRBF model is much better than the classic linear model approach, which completely fails to predict the full range of VF values (Fig. 1a). 

Although the mean absolute prediction error is reduced with BRBF I and II, the standard deviations of the absolute prediction errors of both BRBF models are still large (3.7 and 3.8 dB). This can also be observed in more detail in the prediction limits shown in Figure 1b. Despite of the general similarity between the prediction limits and VF test–retest limits, predictions at the normal end of the range still tend to be lower than the actual VF measurements and at the damaged end, predictions tend to be higher than actual VF measurements. The reported retest differences, ³² on average, are very small at the normal or healthy end. At the damaged end, the retest values tend to be higher, but the median BRBF predictions are slightly higher than the retest values. This likely reflects the difficulty the prediction from RNFLT images have in identifying small, focal defects. Moreover, the floor effect in the VFs and SLP images ^41,42 and the atypical scan pattern in SLP images, which may be associated with glaucoma severity, ⁴³ may be additional causes of the overestimation at the lower end of the VF sensitivity. Furthermore, because the diagnosis of subjects in REH datasets includes structural criteria, the normal subjects in the training dataset may have supernormal structure and the glaucomatous subjects have greater than average structural damage. This potential bias on subject selection may distort the structure–function relationship. Therefore, the training process may be improved by including a range of glaucoma severity in subjects defined only by VF loss. Another factor that may confound the reported prediction accuracy is that the models were trained with VFs tested with both SITA and full-threshold programs, whereas they were evaluated on VFs tested only with full-threshold program. This may, in part, account for the observed tendency for the prediction to overestimate the sensitivities in the test dataset due to the higher sensitivities obtained with SITA compared with the full-threshold program. ⁴⁴ However, this effect would be small because the sensitivity difference with two programs is just 1.3 dB, on average. ⁴⁴  

Recent investigations attempting to elucidate the structure–function relationship in glaucoma generally have the goal of assessing the relative accuracy of structural and functional tests throughout the course of the disease. ^45,46 For example, loss of function without loss in structure that does not adhere to a particular structure–function model's prediction may be an indicator of a nonglaucomatous process, measurement imprecision, or some artifact in the image of the structure or function test. Our long-term goal is to use the BRBF technique to provide a relevant clinical tool that indicates concordance between the VF and the chosen surrogate measure for structural loss. For example, when a VF and structural measure derived from GDx, HRT, or optical coherence tomography (OCT) is available, a chart mapped in VF space will be provided indicating areas where the measurements are in concordance (within a certain range of accuracy and precision) and where they are not, it could provide clinically useful information about the reliability of the individual measurements or diagnostically useful information. 

It is an imperative that any new statistical method should be developed and tested on more than one dataset. ³⁰ We had access to three large, independent datasets, each collected at one of three clinical centers. The inclusion criteria for glaucomatous and healthy subjects were generally consistent across the three samples. However, as the purpose of this study was not to determine diagnostic performance, the precise definitions for glaucoma were less important. In fact, the mixture of data can be viewed as an advantage in the study design. However, further testing on different datasets, especially where realistic estimates of measurement precision have been performed (from test–retest measurement), is still required and is being undertaken in our laboratory and elsewhere. 

The BRBF method is not limited to one type of input of structural measurement or imaging device. It was shown to handle input of the GDx RNFLT peripapillary profile (64 values) as well as the PCA-reduced RNFLT values at 16,512 pixels in a wide peripapillary annulus. It could be used on neuroretinal rim area values from scanning laser ophthalmoscopy technology or RNFLT values derived from OCT technology, or any other surrogate measure of glaucomatous structural loss. Moreover, we have demonstrated, albeit qualitatively with the maps in Figure 5, that using the surrogate measures of RNFLT in their high-dimensional form provides mapping that is closer to the expected structure–function relationship. We speculate that the next generation of Fourier domain OCT instrumentation, now finding its way into clinics and providing volumes of data for RNFLT, will be particularly amenable to this method. 

In conclusion, we have introduced a new statistical method for describing the relationship between functional and structural measurements used in the clinical evaluation of glaucoma. Evidence from a dataset independent of those used to derive the model indicates that the BRBF method has advantages over standard statistical approaches for modeling these relationships, and estimates of functional deficits from structural measures yielded from this method are better than those derived from a classic linear regression approach. This method can provide a platform from which clinically useful tools can be derived for mapping and charting concordance between VF measurements and RNFLT measurements in glaucoma. 

Supplementary Materials

Supplementary Material - (PDF) 

Footnotes

 Supported by a Pfizer Investigator-Initiated research grant and in part by the Department of Health's NIHR Biomedical Research Centre for Ophthalmology at Moorfields Eye Hospital and UCL Institute of Ophthalmology. The views expressed in the publication are those of the authors and not necessarily those of the Department of Health.

Footnotes

 Disclosure: H. Zhu, None; D.P. Crabb, Pfizer Inc. (F); P.G. Schlottmann, None; H.G. Lemij, None; N.J. Reus, None; P.R. Healey, None; P. Mitchell, None; T. Ho, None; D.F. Garway-Heath, Carl Zeiss Meditec (F, C); Optovue, Inc. (F); Heidelberg Engineering (F), Pfizer Inc. (F)

Appendix

Our proposed method of linking the structural and functional measurements is formalized by a function , which predicts a functional measurement y from a structural one x : y = (x). This calculation sets a mathematical framework for the question: what would be the corresponding functional change (Δy) for the structural change Δx : Δy = (x + Δx) − (x)? Because we are interested in subtle structural change to model the slow progression of RNFL damage, we assume that Δx becomes very small and tends toward 0. The structure–function relationship is then defined with the general equation  

 

where is differentiable with regard to x, and ∇_x is the gradient of at x. Since we are sampling subjects rather than considering the whole population, the final term in this equation must be expressed as a statistical expectation of ∇_x or, for simplicity, the mean of ∇_x. 

We examined classic linear regression and BRBF as the choice of function . The latter was extended from the relevance vector machine (RVM),²⁷ where the output is originally only one-dimensional. The extension was similar to the model derived by Thayananthan et al.⁴⁷ In particular, VF sensitivity is assumed to be  

 

where y_dⁿ is the dth element in the measurement of the nth subject, xⁿ is the RNFLT measurement (64-point profile or SLP image) vector, w_d is a weight vector, ε_d is an additive 0-mean Gaussian noise N(0, δ_d²) with variance δ_d², and the radial basis function vector φ(xⁿ) is defined to be M+1 dimensional for M bases: φ(xⁿ) = [1, φ₁(xⁿ), φ₂(xⁿ),…, φ_M(xⁿ)]^T, where each element is a radial basis function with center θ_m and an isotropic covariance  

 

If all weight vectors w_d are organized into a matrix W columnwise, then, using Bayesian methodology, we assign a prior over W 

 

where Ψ is a diagonal matrix whose elements are α₁⁻¹, α₂⁻¹,…, α_M⁻¹ on diagonal, and 0 otherwise. Each α_m⁻¹ represents the average variance of the weights for the mth basis. 

According to Bayesian methodology, priors of hyperparameters are defined over α _m and β _d = δ _d ⁻²  

 

where α and β are vectors of α _m and β _d , respectively, and gamma(a_m | a, b) is a gamma distribution with parameters a and b. 

This framework forms the objective function, which is the probability of all parameters W, α, and β, given the observations Y  

 

where Y is a matrix, the columns of which are y ⁿ for n from 1 to N. 

In the first item in equation 7, it is straightforward to infer that w_d is independent from any other weight vectors, given Y, α, and β. Consequently  

 

where the covariance matrix and mean for the posterior distribution of w_d are  

 

where Φ(X) is a matrix: Φ(X) = [φ(x¹),…, φ(x^N)]^T. 

The latter item in equation 7 can be calculated as: p(α, β | Y) ∝ p(Y | α, β)p(α)p(β) where  

 

With no prior knowledge on α and β, these two parameters are assumed to be “uniformly” distributed so p(α) and p(β) have little impact on p(α, β | Y). When an approximation similar to that of RVM is used,²⁷ the hyperparameters are optimized by setting the derivative of equation 10 to 0  

 

where μ_wd^m is the mth element of the vector μ_wd, and Σ_wd^mm is the diagonal element at mth row and column in Σ_wd. 

The parameters are inferred by iterating between equations 9 and 11 until convergence. The radial basis centers θ _m are initialized to contain all x ⁿ . As with an RVM, many of α _m become infinite during the training process so the corresponding radial bases are removed accordingly. The radial basis parameter η _m is optimized by gradient descent. 

Given the inferred parameters, the distribution of predicted values, y ′_d given a test example x′, is computed as  

 

where the prediction is made by the mean of the distribution μ _{w d} ^Tφ(x′). 

Therefore, the structure–function relationship in equation 2 is implemented by  

 

where  

 

and  

 

is a vector −2η_mφ_m(xⁿ) (xⁿ − θ_m), according to equation 4.