A Functional Regression Model of the Retinal Nerve Fiber Layer Thickness in Healthy Subjects

Abstract

Purpose: A new functional regression model is presented to explain the intersubject variability of the circumpapillary retinal nerve fiber layer (RNFL) thickness in healthy subjects.

Methods: To evaluate the functional regression approach we used data from 202 healthy volunteers, divided equally into training samples (TS) and validation samples (VS). Covariates included RNFL, fovea distance, fovea angle, optic disk ratio, orientation and area provided by Fourier-domain–optical coherence tomography, age, and refractive error. Root mean square errors (RMSE) were calculated for each of the 256 sectors and for the 12 clock-hour sectors in the TS and VS and were compared to the RMSE of the previous model and the standard deviation of the raw data.

Results: With the functional regression approach, we were able to explain on average 27.4% of the variation in the TS and 25.1% of the variation in the VS. The new model performed better compared to a multivariate linear regression model. It performed best in the superior-temporal and inferior-temporal clock-hour sectors where the percentage of RMSE reduction ranged between 26.3% and 44.1% for the TS and between 20.6% and 35.4% for the VS.

Conclusions: The new functional regression approach improves on the multivariate linear regression model and allows an even larger reduction of the amount of intersubject variability, while at the same time using a substantially smaller number of parameters to be estimated.

Translational Relevance: The demonstrated reduction of interindividual variation is expected to translate into an improved diagnostic separation between healthy and glaucomatous subjects, but this remains to be demonstrated in further studies.

Introduction

Glaucoma is a progressive degenerative optic neuropathy caused by damage of ganglion cell axons leading to irreversible loss of ganglion cells and their axons. This results in thinning of the circumpapillary retinal nerve fiber layer (RNFL), loss of visual field (VF), and without proper treatment, to permanent blindness. Current diagnosis methods consider imaging techniques such as optical coherence tomography (OCT), which, succinctly, measures RNFL thickness and compares it with a normative database.¹ However, these measurements currently present a high intersubject variance,² mainly due to individual anatomical factors. Because it is the most frequent cause of irreversible blindness in industrialized nations, with an estimate of about 112 million people affected by the disease,³ it is crucial to find more accurate methods that can improve diagnosis. 

Our recent efforts have been focused on developing and implementing a compensation method that takes into account individual anatomical parameters that may contribute and compensate for the intersubject variability.^4–6 Specifically, a multivariate linear regression approach including retinal vessel density (RVD)—a parameter reflecting the locations and thicknesses of all measurable circumpapillary retinal vessels—as well as optic disc shape descriptors and fovea parameters allowed a significant decrease in intersubject variability of RNFL.⁶ We reported a reduction of 18% on average, and up to 29% measured in 12 clock-hour sectors, of the coefficient of variation of RNFL thickness. This previous approach considered 256 independent models around the optic disc (OD) with values of RNFL thickness at a given sector as the dependent variable and all putative parameters, including the values of RVD at the respective sector, as independent factors. Although achieving an actual reduction of intersubject variability, this approach had, from a statistical perspective, some serious shortcomings. Model selection was performed for each of the 256 sectors individually, resulting in potential overfitting due to the large number of parameters to be estimated. Furthermore, the apparent correlation of RNFL between neighboring sectors was not taken into account at all. We propose here a functional regression approach that is meant to overcome these weaknesses. 

Functional regression analysis is becoming more and more important in many research fields. The basic ideas are presented by Ramsay and Silverman,⁷ whereas more recent developments in terms of application have been described by Morris.⁸ Here we modeled for each subject the entire RNFL curve, composed of 256 thickness values, using only 16 basis functions. The coefficients corresponding to these basis functions were then used as dependent variables in regression models with eight regressors: vessel thickness plus seven other explanatory variables that were used in the previous multivariate approach. Model selection was performed for each of the 16 linear models. The resulting number of coefficients to be estimated was much smaller than in the multivariate approach, while achieving a better fit of the data, as will be demonstrated by analyzing both the training data set as well as an independent validation sample. 

Methods

Subjects and Parameter Extraction

Available data from 202 healthy subjects were randomly split in a training data set (TS) and a validation data set (VS). The research followed the tenets of the Declaration of Helsinki. Informed consent was obtained from the subjects after explanation of the nature and possible consequences of the study. The demographics of the overall sample and description of the ophthalmic examinations performed for all subjects included are described elsewhere.⁶ Shortly, OCT examinations of the OD (cube, 200 × 200) and the macula (cube, 512 × 128) were acquired with Fourier-domain OCT (FD-OCT) (Cirrus; Carl Zeiss Meditec, Inc., Dublin, CA, USA). Scans with a quality index lower than 6 (ranging from 0 to 10) or movement artifacts within the measurement circle were excluded. The image data were exported and analyzed with the Cirrus Research Browser (software version 6.0.2.81), which included segmentation of the RNFL of the 3.4-mm peripapillary circle. 

The automated parameter extraction was performed as previously published when concerning OD and fovea descriptors. Retinal vessel segmentation was performed using the same methods previously described.⁶ Using the vessel trees generated in the OD-centered images, we considered all vessels within a band of diameter around the OD center extending from 3.28 to 3.64 mm to integrate a 256-sector vessel profile. This discrete profile was based on vessel thickness and position relative to the OD (angle). The thickness of each vessel was attributed to a specific sector (out of 256 sectors, each 1.4° wide) according to its angle. As opposed to the linear regression model previously presented, we considered only a discrete profile representative of retinal vessel location in a circumpapillary profile, that is, without considering a Gaussian convolution.⁶ 

Statistical Analysis

The functional regression approach that was implemented in R version 3.3.2, making use of the fda.usc package,⁹ will be only briefly outlined here. A more detailed description is given in the Supplementary Material. The first step consisted of landmark registration with respect to three characteristic points of the RNFL curves: the minimum near the fovea and two flanking maximum peaks. After aligning individual curves according to these landmark points, 16 functional principal components were computed and used as base functions to approximately model the aligned RNFL curves. The corresponding principal component analysis (PCA) loadings were then considered as dependent variables in a regression model including vessel thickness and seven additional subject-specific characteristics as explicatory variables. To be able to compare our results with those from Pereira et al.,⁶ the subject-specific variables included age, spherical refractive error (RE), OD area (ODA), OD orientation (ODO), OD ratio (ODR), as well as fovea angle (FA), which is the angle between the fovea and the OD centers and a horizontal line, and fovea distance (FD) from the OD center. The final step consisted of model selection based on the Akaike Information Criterion (AIC) resulting in a model with 49 parameters. 

The functional regression approach and the multivariate model of Pereira et al.⁶ were applied to the 101 subjects of TS. RNFL thickness profile, estimated with the functional regression approach for the aligned curves, was transformed back to the original coordinate system by linear interpolation, and residuals were calculated at the 256 original sectors. The different models were assessed according to the root mean square error (RMSE) at each sector. Moreover, for each individual the overall RMSE and the RMSE over 12 clock-hour sectors were calculated for both approaches, and differences were tested using pairwise t-tests. The validation data set underwent the same kind of analysis to compare multivariate linear regression and functional data analysis. The validation step, performed in a sample different from the one used to train the model, ensures that differences in performance were not due to overfitting the training data set. 

Results

After backward selection using AIC, 49 coefficients were estimated in our functional regression model. By multiplication of the covariate specific coefficients with the basis functions, we obtained temporal-superior-nasal-inferior-temporal (TSNIT) profiles for each subject characteristic (Fig. 1). These curves show the amount of expected change in RNFL thickness according to one unit change in the subject characteristic. 

Figure 1

 

Modified TSNIT profiles on the aligned coordinate system for the seven parameters: age, ODR, ODO, ODA, FD, FA, and refractive error.

Age did not show a strong effect throughout the complete RNFL profile, with −0.036 μm per year on average. An increase in ODR was associated with a thicker RNFL in the area around the peaks, leaving the peaks relatively unchanged. In addition, RNFL was reduced in the nasal and temporal parts of the measurement circle. Subjects with a larger ODO, meaning a superior pole rotated toward the macula, tended to have thicker RNFL nasal inferior and thinner RNFL nasal superior. RNFL thickness generally increased with bigger ODA, and did so especially in the areas around the two peaks. 

Regarding fovea parameters, larger FD was clearly associated with increased RNFL thickness in the temporal parts of the measurement circle and with reduced RNFL thickness in nasal inferior. A bigger FA (i.e., the fovea position being higher up in the macula) resulted in larger RNFL thickness around the superior peak and even much more so on the nasal side of the superior peak. On the nasal side of the inferior peak, subjects with larger FA showed a steeper decline of RNFL thickness. 

Regarding refractive error, a positive value (hyperopia) was associated with slightly higher RNFL thickness in the nasal area but lower RNFL thickness in the peaks and on the temporal side of the peaks. 

Compensated RNFL Model Reveals a Clinically Significant Reduction of RNFL Variation

Figure 2 provides for each sector the RMSE over all subjects in the TS (blue curve). We compared the RMSE of the multivariate model from Pereira et al.⁶ (red curve) with the RMSE of the simplest possible approach of considering the overall mean curve without any subject-specific information (null model, black curve). Note that in this last case the RMSE is identical to the standard deviation of the sample. 

Figure 2

 

RMSE over null model (black), multivariate model (red), and functional regression approach (blue) for the training data set.

RMSE averaged over the 12 clock-hour sectors is presented in Table 1 for all three approaches. Corresponding results for the coefficient of variation of RMSE are presented in Part B of the Supplementary Material. Numbering of the clock-hour sectors starts with sector 1, corresponding to the 9 o'clock sector in the temporal side of the eye. With the functional regression model, the average RMSE was reduced from 18.97 to 13.78 μm, meaning that we were able to explain, on average, 27.4% of the variation in the TS, significantly more than the multivariate linear regression model approach (P-value t-test < 0.0001), which reduced the RMSE from 18.97 to 15.42 μm, which translates to a reduction of 18.7%. The new model performed best in the superior-temporal and inferior-temporal areas (clock-hour sectors 2–4 and 10–11), where the percentage of reduction ranged between 26.3% and 44.1%, respectively. In these sectors, the reduction was also significantly higher than the reduction obtained the multivariate model (15.4% to 22.0%). Only in sector 8 did the new functional approach performed significantly worse than the multivariate linear model, with a reduction of 11.7% compared to 21.2%. 

Table 1

 

RMSE of the 12 Clock-Hour Sectors Over the 101 Subjects of the Training Data Set for the Three Different Models^a

 

Furthermore, we compared the width of the 5th- and 95th-percentile band of errors in the TS between the three models (Fig. 3). This comparison is especially relevant when identifying pathologic subjects. It is evident that the 5th- to 95th-percentile band of the functional regression model was smaller over all 256 sectors than the raw data percentile band. It was also smaller than the percentile band of the multivariate approach in the temporal-superior and temporal-inferior area of the circumpapillary retinal band while having comparable width in the nasal and temporal areas. 

Figure 3

 

Five to ninety-five percent percentile bands of residuals: black, null model; red, multivariate model; blue, functional regression model.

Comparing the Effect of RNFL Alignment, RVD, and Other Parameters

To investigate the amount of variance reduction due to each of the three components (the alignment, the vessel correction, and the influence of the subject's characteristics) we compared the results of several reduced models to the full model. Table 2 shows the RMSE and which percentage of the raw data variance could be explained (in percentages) for each model for each of the 12 clock-hour sectors separately and for the average over all 256 sectors. A reduction of RMSE by 14.8% occurs already because of the alignment, which especially leads to reduced variation in the clock-hour sectors 1 to 4 and 9 to 12, but increases variation in the clock-hour sectors 5 to 8. Quite a large proportion of the variance can be explained by adding retinal vessels to the model, with additional 9.4% reduction in the mean RMSE over all 256 sectors. The model combining alignment and vessels results in variance reduction in all of the clock-hour sectors (between −5.6% and −40.7%). By correcting for the other seven characteristics of subjects, we gain 4.2% compared to the alignment-only model and 3.2% when comparing the full model to the model including only alignment and vessel correction. 

Table 2

 

RMSE of the Null Model and Reduction in Percentage for the Other Models in Each Clock-Hour Sector and Averaged Over All 256 Sectors for the Training Data Set

 

Successful Validation of the Functional Regression Model in a Separate Sample

In the validation step, the functional regression and multivariate linear regression models obtained in the TS were applied to an independent data set of 101 healthy volunteers, the VS. All physiological parameters were calculated in a way similar to those in the TS. The resulting RMSEs from the three models are presented in Figure 4 for all 256 sectors and in Table 3 for the 12 clock-hour sectors. Results for the coefficient of variation are presented again in Part B of the Supplementary Material. 

Figure 4

 

RMSE over null model (black), multivariate model (red), and functional regression approach (blue) for the validation data set.

Table 3

 

RMSE in for the 12 Clock-Hour Sectors and Averaged Over All 256 Sectors for the Three Different Models in the Validation Data Set

 

In the VS the functional regression model reduced the average RMSE from 19.55 to 14.65 μm, a reduction of 25.1%. Similar to the TS, the model performed best in the superior and inferior areas (clock-hour sectors 2–3 and 10–11) where the percentage of reduction in the VS ranged between 20.6% and 35.4%. The multivariate model reduced the RMSE from 19.55 to 15.87 μm, a reduction by 18.8%, which is significantly less than the functional regression model (P-value t-test < 0.0001) 

In addition, the amount of variance reduction due to each of the three components (the alignment, the vessel correction, and the influence of the subject's characteristics) was quite comparable between the VS and the TS. Table 4 shows the RMSE for each of the clock-hour sectors and averaged over all 256 sectors (mean in micrometers), as well as the resulting reduction compared with the null model (in percentages) for each model. The reduction of RMSE due to alignment only was 13.5% (compared to 14.8% in the TS), again with reduced variation in sectors 1 to 4 and 9 to 12, but with increased variation in sectors 5 to 8. Taken together with vessel correction, 23.6% of the RMSE could be explained on average over all 256 sectors (which is similar to the reduction in TS with 24.2%). Only the reduction of the RMSE due to the subject's characteristics was smaller in the VS than in the TS: further 1.0% if applied to the aligned data but without correcting for the vessels and 1.5% if included in a model with alignment and vessel correction. 

Table 4

 

RMSE of the Null Model and Reduction in Percentage for the Other Models in Each Clock-Hour Sector and Averaged Over all 256 Sectors for the Validation Data Set

 

Discussion

In this study we present an improvement of the multivariate linear regression model previously published by our team.⁶ This new model is based on functional regression, a relatively new development in statistical methodology, which allows for functional modeling of the circumpapillary RNFL thickness profile. In this way, the RNFL thickness is treated as a continuous function and correlation between different sectors is properly accounted for, in contrast to the previous approach in which independent linear models for each sector were considered. A further advantage of the new method is that the number of parameters to be estimated is much smaller than in the previous method, which reduces the danger of overfitting the training data set. Each RNFL curve is no longer represented by 256 points but rather by 19 parameters corresponding to three landmark parameters and 16 basis functions. In terms of estimation, we are dealing with a total of 128 parameters, out of which 49 were selected with AIC to obtain one model representing the RNFL of the complete circumpapillary measurement circle, as compared to 2048 parameters in 256 separate models in the multivariate linear regression approach. The new method also overcomes the issue of having a different best-fitting model in each sector, and as seen in Figure 1, the influence of subject-specific covariates is modeled as a smooth function over all sectors in stark contrast to the previous model. 

Our particular functional regression approach includes an initial landmark registration. In the past, it has been proposed to align the RNFL measurement circle according to the disc-fovea line, aiming to partly reduce the interindividual variability. However, this approach equals a simple rotation of the RNFL measurement circle, which has been demonstrated to fail in reducing interindividual variability.^10–12 The difference in our present approach is that both RNFL peaks and the temporal RNFL minimum are used to align the RNFL measurement circle in a way that results in different amounts and directions of rotation around the circle for a given subject. This landmark registration procedure reduces, on average, the intersubject variability by about the same amount as the joint influence of vessels and other subject-specific covariates. There are some clock-hour sectors where our landmark registration slightly increases the RMSE, but only in those sectors where even the null model gives already relatively small errors, which was nasally between 1 and 5 o'clock. In those sectors where the null model is doing particularly badly (close to the peaks P1 and P2, see Supplementary Figure S1), landmark registration gives the largest reduction in RMSE. 

When comparing the TSNIT profiles of the individual parameters, it is of interest to note that the TSNIT profiles of the functional regression model (Fig. 1) are presented in the aligned coordinate system, whereas Pereira et al.⁶ gave results from the multivariate linear regression model in the original coordinate system. In general, all factors had fairly similar behavior along the TSNIT profile for both approaches. In accordance with our previously reported results, age did not show a strong effect throughout the complete RNFL profile, with a loss of 0.036 μm per year on average. In the literature, this loss has mostly been found to be between 0.2 and 0.5 μm per year.^13–21 The reason for the weaker correlation in our study might be that our samples are composed by mainly young subjects with a mean age around 31 years. Some authors have published evidence that age-dependent decline in RNFL may be steeper after an age of 40 to 50 years,^19,20 although this has not been confirmed in a recent longitudinal study.¹⁸ 

An increased ODR (more elongated OD) resulted in a broader RNFL peak and a reduced RNFL thickness nasally and temporally in the functional regression model. The former has not yet been described, to the best of our knowledge, while the latter has also been a finding in our multivariate linear model.⁶ The difference between the models might be explained by the fact that the peak location is variable, which is now reflected in the model by the alignment but was not so in the multivariate linear model. Our results confirm previous findings of our own group and others that bigger ODA is associated with thicker RNFL over all sectors but is especially strong in the areas around the two peaks.^6,22 

As in our previously published multivariate model, increasing FD still is associated with increasing RNFL in temporal areas.⁶ This might reflect the temporal shift in the RNFL profile in subjects in which the distance between fovea and OD center was increased, as previously reported.²³ FA presents a highly positive association with RNFL in superior through superior-nasal areas and a negative influence in nasal and inferior sectors, which may represent the rotation of the superior RNFL peak toward temporal sectors in cases in which the fovea is located more inferiorly to the OD. These findings essentially confirm those obtained with our old approach.⁶ 

Regarding refractive error, we found a negative association with RNFL in temporal areas, which confirms our previous findings.⁶ However, with the functional regression model the maximum of this negative association is around the peak areas of RNFL, while in the multivariate linear model there was a positive association at the superior peak. The difference between the two models likely is related to the landmark registration, which now corrects for different peak locations. In the old model, these differences had to be accounted for by other factors such as the vessels and also by the refractive error. It has been previously reported^6,24 that in myopic eyes both the retinal vessel arcades and the arcuate bundles of the RNFL tend to be displaced to the temporal side. 

The functional data analysis increases the reduction of intersubject variability compared to the multivariate model approach. We observed a reduction of 27.4% of RMSE (as opposed to 18.7% in the multivariate approach) in the training sample and a reduction of 25.1% (compared to the previous reduction of 18.8%) in the validation sample. In both samples, the difference between the multivariate linear approach and the new functional data analysis was statistically significant. The validation of the functional data model in an independent data set indicates that there was no serious problem of overfitting. However, despite theoretical concerns about possible overfitting, as discussed above, the previous multivariate model was performing quite well in the validation step. 

The improvement by the functional regression model as compared to our old model was mainly concentrated on the superior-temporal and the inferior-temporal sectors, which are of major clinical relevance. These sectors represent the locations with the best discrimination between glaucoma suspects and glaucoma patients.²⁵ Furthermore, the majority of glaucoma patients show RNFL defects at these locations.²⁶ Reducing the interindividual variation may thus be especially valuable in those sectors. 

What would be the relevance of reducing the intersubject variability of RNFL measurement? It has been stated by others that if interindividual variability can be reduced, it should be possible to improve the sensitivity and specificity of tests based upon OCT RNFL thickness.²⁷ This is to be expected since usually a factor with known association with RNFL thickness (for example, age) can be used to either correct for its impact on RNFL thickness or to construct an age-corrected normative database—two sides of the same coin. The effect of age on RNFL thickness as measured by OCT has been first described by Bowd et al. in 2002,¹³ where it explained less than 17% of the total interindividual variance. Since then, the effect of age on RNFL has been confirmed by numerous authors,^{14–17,20–21,28} but to the best of our knowledge nobody as yet has evaluated whether compensation for age reduces interindividual variation of RNFL thickness or improves diagnostic separation. Nevertheless, the normative databases of OCT devices take into account the effect of age. In our opinion, this discrepancy between clinical use and lack of scientific proof reflects the natural acceptance by human reasoning that correcting for age will improve the diagnostic performance of OCTs. It is obvious for us that an aged person with thus reduced but healthy RNFL must not be compared with normal limits of young people or even a mixed population without taking into account the effect of age on RNFL. Consequently, as yet nobody has tested this hypothesis. If we take into account several physiologic factors, as we do here, that are all associated with RNFL thickness, the same argument should hold: If these factors explain a considerable part of the interindividual variation (and they do), then this is of clinical relevance and should result in a better diagnostic separation. 

There are, of course, caveats and therefore limiting factors for the above-outlined reasoning: The factors we are using should not themselves be involved in the disease process: We think this is not the case, but we cannot completely exclude this possibility. An effect of overfitting might lead to an overestimation of the effect of the associations, but this is unlikely since we have confirmed our findings in a separate validation sample. In addition, we might not have considered all factors that have an influence on RNFL in the current model. Previous reports demonstrate that correcting for additional factors may further reduce intersubject variability of RNFL²⁹ and, in the case of our multivariate model, may also influence some of the described correlations.^24,30 In terms of age distribution, both data sets (training and validation) are composed of relatively young subjects (under 40 years old), which may mask or minimize some meaningful impacts that are currently not explicit in the analysis. This fact, however, does not invalidate any of the results presented, and it does not contradict any of the previous reports correlating age and RNFL thickness.^20,21,28,31 Moreover, it would be desirable to redefine OD descriptors. The fact that three descriptors (orientation, ratio, and area) are considered may bring some competitiveness between parameters, specifically between ratio and orientation, both dependent on the major axis defined by the OD contour. 

One limitation when applying the model to glaucoma patients with moderate to advanced damage might be that at least one of the major RNFL peaks could be missing. This may impact the landmark registration, and as yet it remains unclear how this would impact the performance of the model. However, the main purpose of applying a model such as ours will be to enhance early diagnosis of glaucoma by reducing the interindividual variance and thus be able to detect a smaller amount of RNFL reduction with an increased specificity. Although this remains to be demonstrated in future studies, the limitation discussed above will likely not be very relevant in those patients because degeneration has not yet been severe enough to eliminate RNFL peaks. 

This new update of our previously developed model confirms and further emphasizes the crucial importance of considering anatomical parameters, specific to each subject, as a strong influence on RNFL thickness distribution. The idea presented previously is now confirmed and improved and likely would be suitable for clinical routine. As with age in present normative databases, a more complex model such as ours can be used to calculate expected values of the RNFL. Based on this expected value and the variance of the model, a lower limit of normality may then be calculated. We would call such normal limits individual normal limits, since no two subjects would be equal in all parameters included in the model. 

To summarize, this manuscript presents a statistically more sophisticated method to explain, correct for, and thereby reduce interindividual variation in RNFL thickness. The new method uses functional regression analysis and significantly improves on the previously published multivariate linear model in terms of reducing interindividual variation of the RNFL by 25% (compared to 19%). The largest improvements are observed in the inferior-temporal and the superior-temporal sectors, with up to 35% reduction of variation. These sectors also represent the most important regions to detect early glaucomatous RNFL damage. The relevant reduction of interindividual variation is expected to translate into an improved diagnostic separation between healthy and glaucomatous subjects, but this will have to be demonstrated in future work. 

Acknowledgments

Supported by Grant LS11-046 from Vienna Science and Technology Fund. 

Disclosure: I. Pereira, None; E. Pablik, None; F. Schwarzhans, None; H. Resch, None; G. Fischer, None; C. Vass, None; F. Frommlet, None