A Bayesian Hierarchical Spatial Longitudinal Model Improves Estimation of Local Macular Rates of Change in Glaucomatous Eyes

Abstract

Purpose: Demonstrate that a novel Bayesian hierarchical spatial longitudinal (HSL) model improves estimation of local macular ganglion cell complex (GCC) rates of change compared to simple linear regression (SLR) and a conditional autoregressive (CAR) model.

Methods: We analyzed GCC thickness measurements within 49 macular superpixels in 111 eyes (111 patients) with four or more macular optical coherence tomography scans and two or more years of follow-up. We compared superpixel-patient-specific estimates and their posterior variances derived from the latest version of a recently developed Bayesian HSL model, CAR, and SLR. We performed a simulation study to compare the accuracy of intercept and slope estimates in individual superpixels.

Results: HSL identified a significantly higher proportion of significant negative slopes in 13/49 superpixels and a significantly lower proportion of significant positive slopes in 21/49 superpixels than SLR. In the simulation study, the median (tenth, ninetieth percentile) ratio of mean squared error of SLR [CAR] over HSL for intercepts and slopes were 1.91 (1.23, 2.75) [1.51 (1.05, 2.20)] and 3.25 (1.40, 10.14) [2.36 (1.17, 5.56)], respectively.

Conclusions: A novel Bayesian HSL model improves estimation accuracy of patient-specific local GCC rates of change. The proposed model is more than twice as efficient as SLR for estimating superpixel-patient slopes and identifies a higher proportion of deteriorating superpixels than SLR while minimizing false-positive detection rates.

Translational Relevance: The proposed HSL model can be used to model macular structural measurements to detect individual glaucoma progression earlier and more efficiently in clinical and research settings.

Introduction

Detection of change over time is often important for the proper treatment of chronic diseases and is crucial to the management of glaucoma. Glaucoma is particularly suited to applying statistical and machine learning models to detect disease and to identify its progression, because most disease-related outcome measures are quantifiable.^1–6 Relevant outcomes include the optic disc size, neuroretinal rim, retinal nerve fiber layer, or macular thickness measures along with a host of quantitative variables related to visual fields. 

Global, regional, or local structural and functional measures have frequently been used as metrics to detect differences of a specified amount from baseline, commonly called event analyses, or to estimate rates of change, also known as trend analyses.^7–9 Longitudinal analyses of local structural or functional measures frequently rely on repeated simple linear regression (SLR) of such measures from single eyes against time to estimate eye-specific rates of change at macular superpixels, optic disc or retinal nerve fiber layer sectors, or visual field test locations or clusters.^9–13 Although SLR is easy to do and seems intuitive, there are multiple issues with this approach, which although occasionally acknowledged, have mostly been ignored.^14–16 Simple linear regression models ignore the spatial relationships among local structural (sectors, superpixels) or functional (visual field clusters or locations) measures. Another shortcoming of SLR is that valuable population information from the cohort is not used to refine estimated rates of change. Similarly, there is no formal way to specify the correlation of baseline measurements (the intercepts) with slopes (rates of change). This is an important issue because baseline measurements influence the magnitude of rates of change in an optical coherence tomography (OCT) sector or superpixel.⁹ Accounting for baseline structural measurements can lead to a reduction of the estimated variability for rates of change (slopes). A few prior studies have proposed linear mixed or hierarchical models to address some of the above shortcomings.^14,15,17,18 For example, Montesano et al.¹⁴ applied a hierarchical linear mixed model to 24-2 visual field data to estimate global and local rates of change in individual eyes. To account for spatial correlation, Betz-Stablein et al.¹⁹ developed a model using conditional autoregressive (CAR) priors on intercepts and slopes for individual-level 24-2 visual field data. However, in either case, no attempt was made to include the cohort's population data to help estimate individual eye parameters. The cohort's data has information on the distribution of possible slopes, and using this information as a prior results in more accurate slope estimates. Such models are challenging to specify correctly and run and can require a significant investment of time and computer CPU power. 

Our team has developed several versions of a novel Bayesian hierarchical spatial longitudinal (HSL) model to improve estimation of local macular thickness rates of change in a prospective cohort of glaucoma patients.^4,20 This ongoing project aims to provide a longitudinal framework to estimate global and most importantly local rates of change more precisely across the macula within individual eyes while at the same time overcoming the inadequacies of the SLR approach. Each successive model improves on the previous version with additional modeling features. The current model adds spatially correlated patient-specific random effects and spatially correlated visit effects. The goal of the current work is to compare the performance of the latest HSL model to that of SLR and a CAR model. We hypothesized that the HSL model would provide more accurate estimates of local macular rates of change thanks to reduced variance for those rates of change compared to SLR and more appropriate shrinkage toward the population averages by location than CAR. 

Methods

One hundred eleven eyes (111 patients) from the Advanced Glaucoma Progression Study (AGPS) were included in this study. The AGPS is an ongoing longitudinal study at the University of California Los Angeles. The university's Institutional Review Board approval was obtained for this study. The study adhered to the tenets of the Declaration of Helsinki and conformed to Health Insurance Portability and Accountability Act policies. All patients provided written informed consent at the time of enrollment in the study. Inclusion criteria for enrolled eyes were as follow: (a) clinical diagnosis of primary open-angle glaucoma, pseudoexfoliative glaucoma, pigmentary glaucoma, or primary angle-closure glaucoma and (b) evidence of either central damage on 24-2 visual field (VF), defined as two or more points within the central 10° with P < 0.05 on the pattern deviation plot or VF mean deviation (MD) worse than −6 dB. Exclusion criteria consisted of baseline age less than 40 years or greater than 80 years, best-corrected visual acuity <20/50, refractive error exceeding 8 diopters of sphere or 3 diopters of cylinder, or significant retinal or neurological disease affecting OCT measurements. Eyes with four or more OCT scans and two or more years of follow-up were included. We analyzed observations up to 4.25 years after baseline. Data from visits less than 0.2 years after a previous visit were omitted.²⁰ 

Macular OCT Imaging

Macular volume scans were obtained with a Spectralis spectral-domain OCT (Heidelberg Engineering, Heidelberg, Germany). The Posterior Pole Algorithm of the Spectralis OCT acquires 30° × 25° volume scans of the macula (61 B-scans spaced approximately 120 µm apart) centered on the fovea and repeated nine to 11 times to reduce speckle noise. Proprietary software of Spectralis OCT, the Glaucoma Module Premium Edition, was used to automatically segment individual retinal layers before data export. Images were reviewed for segmentation errors and image artifacts. Segmentation errors were manually rectified with the OCT device's built-in software. After segmentation, the individual layer thickness measurements are provided as 8 × 8 arrays of 3° × 3° superpixels for the central 24° × 24° region centered on the fovea. Because of substantial measurement noise in the peripheral macular regions, we only included the ganglion cell complex (GCC) thickness within 7 × 7 arrays of superpixels from Spectralis OCT macular volume scans after excluding the most inferior row and nasal column of superpixels. The GCC thickness was calculated by adding the thickness measurements of the retinal nerve fiber layer, ganglion cell layer, and inner plexiform layer. 

Data Management and Statistical Modeling

Our methods have been described previously.²⁰ We identified potential outliers as observations with very large increases or decreases between consecutive measurements and removed observations that produced the greatest reduction in the sum of absolute consecutive differences resulting in removing approximately 0.5% of observations as outliers. We removed all observations for a single person in a single superpixel if we identified two or more outliers in that superpixel.²⁰ We then fit a Bayesian normal hierarchical random effects model using the JAGS package in R (R2jags).^21,22 The current version of our Bayesian HSL model includes (a) a macula-wide global intercept and slope, (b) superpixel-specific intercept and slope deviations from the global intercept and slope (superpixel-level random effects), (c) macula-wide patient-specific intercept and slope effects (global patient random effects), (d) patient-superpixel specific intercept and slope random effects, and (e) macula-wide visit effects. Letting y_ijk denote a single observation of GCC thickness (µm) for patient i at time t_ij in superpixel k, the model is  

\begin{eqnarray*} {y}_{ijk} &=& \ {\mu }_0 + {\mu }_1{t}_{ij} + {\alpha }_{0k} + {\alpha }_{1k}{t}_{ij} + \ {\beta }_{0ik} + {\beta }_{1ik}{t}_{ij}\\ &&+\; {\phi }_{0{\rm{k}}} P{C}_{0i}\ + {\phi }_{1k}P{C}_{1i}{t}_{ij}\ + {\phi }_{2k}V{E}_{ij} + {\epsilon }_{ijk} \end{eqnarray*}

 

\begin{eqnarray*} {\epsilon }_{ijk}\sim N\left( {0,\ \sigma _k^2{\rm{exp}}\left( {P{C}_{2i}} \right)\ } \right), \end{eqnarray*}

where µ₀ is the global macula-wide intercept; µ₁ is the global macula-wide slope; α_0k is the population average intercept in superpixel k; α_1k is the population average slope in superpixel k; β_0ik is the patient-superpixel intercept for the ith patient; β_1ik is the patient-superpixel interaction slope for the ith patient; ϕ_0k is the spatial effects eigenvector component for the patient-specific intercepts; ϕ_1k is the spatial effects eigenvector component for the patient-specific slopes; ϕ_2k is the spatial visit effects eigenvector component; PC_0i is the macula-wide patient-specific intercept random effect; PC_1i is the macula-wide patient-specific slope random effect; PC_2i is the macula-wide patient-specific log residual variance random effect; VE_ij is the macula-wide visit random effect; the population component of the patient-superpixel log residual variance is \(\log \sigma _k^2\). Our model uses novel, semi-parametric spatial effects for the intercept ϕ_0kPC_0i, slope ϕ_1kPC_1i, residual variance \(\sigma _k^2{\rm{exp}}( {P{C}_{2i}} )\), and visit effects ϕ_2kVE_ij. In each component, terms with a k subscript ϕ_0k, ϕ_1k, \(\sigma _k^2\), and ϕ_2k define the spatial pattern of variation across the macula whereas the random effects PC_0i, PC_1i, exp(PC_2i), and VE_ij are random patient effects, or visit effects, that indicate how much patient i's observations deviate from the population averages. Each of these four terms are factor analytic models with one factor each where the ϕ_0k, ϕ_1k, \(\sigma _k^2\), and ϕ_2k are called loadings and the random effects are also called factors.²³ Superpixel-specific population intercepts and slopes are the sum of the macula-wide global intercept and slope plus the superpixel-specific intercept and slope deviations. Similarly, a superpixel-patient slope (intercept) is the sum of (a), (b), (c), and (d): the sum of the macula-wide global slope (intercept), superpixel-specific slope (intercept), the macula-wide patient-specific intercept slope (intercept), and the patient-superpixel specific slope (intercept) random effect. Patient-superpixel log residual variances are modeled as the sum of a superpixel component and a patient-specific component. The outlier removal algorithm and full model and priors are given in the Supplemental material. 

In Bayesian models, inference is made by summarizing the posterior distributions of the parameters of interest. The posterior distribution is obtained by combining the likelihood (information in the data) with the prior (prior knowledge about the unknown parameters) and quantifies the uncertainty in the parameters after observing the data. The posterior mean (i.e., mean of the posterior distribution of a parameter of interest) is often taken as the Bayesian point estimate; similarly, the posterior standard deviation (SD), that is, the SD of the posterior distribution, summarizes the uncertainty in the posterior. Posterior means and SDs of superpixel-patient-specific intercepts and slopes were calculated from the HSL model and compared to superpixel-patient-specific intercepts and slopes obtained from a Bayesian SLR model with a flat prior analyzing longitudinal data from each superpixel-patient separately and a Bayesian CAR model with intrinsic CAR priors analyzing data from each patient separately. 

The Bayesian SLR model produces posterior means and posterior Credible Intervals (CrI) that are identical to the classical point estimate and confidence intervals (CI) from classical least squares regression on the same data. The Bayesian SLR posterior SDs are slightly larger than classical SLR standard errors (SEs) by a factor of \({( {\frac{{{n}_i - 4}}{{{n}_i - 2}}} )}^{\{ {1/2} \}}\), where n_i is the number of observations for patient i because SEs are not standard deviations and SEs do not account for the classical estimates being t distributed with n_i − 2 degrees of freedom. Using the posterior SD from this Bayesian SLR model puts SLR, CAR, and HSL on an equal footing. 

For the CAR model inspired by Betz-Stablein et al.,¹⁹ we assumed GCC thicknesses in neighboring superpixels would be more similar than in more distant superpixels. We defined the neighborhood structure, or adjacency matrix, such that superpixels sharing an edge are weighted as 1. For each patient i, we model spatial dependence in intercepts and slopes  

\begin{eqnarray*} {y}_{jk} = \ {\alpha }_0 + {\alpha }_1{t}_j + \ {\beta }_{0k} + {\beta }_{1k}{t}_j + {\epsilon }_{jk} \end{eqnarray*}

 

\begin{eqnarray*} {\epsilon }_{jk}\sim\ N\left( {0,\ {\sigma }^2} \right), \end{eqnarray*}

where α₀ is the overall eye intercept; α₁ is the overall eye slope; β_0k is the superpixel-specific spatial deviation from the overall eye intercept; β_1k is the superpixel-specific spatial deviation from the overall eye slope; and σ² is the residual variance specific to patient i. We fit the CAR models using the R package NIMBLE v0.13.^24,25 The full CAR model and priors are given in the Supplemental material. 

We identified superpixel-patient-specific rates of change as worsening or improving (significantly negative or positive) when the upper or lower limit of the symmetric 95% CrI was less than or greater than zero, respectively. We compared the proportion of significant negative or significant positive slopes at each superpixel from the HSL and SLR models using an α = 0.05 level classical McNemar's test with the null hypothesis being that the two proportions were equal. 

We summarize differences (SLR − HSL and CAR − HSL) of posterior means and of SDs for individual patient-superpixel intercepts and slopes. We compared differences of posterior means averaged within each of 49 superpixels and also averaged across all 5419 patient-superpixel profiles. We similarly compared the SLR/HSL and CAR/HSL ratios of posterior SDs for superpixel averaged intercepts and slopes, and patient-superpixel intercepts and slopes. Ratios of posterior variances (i.e., ratios of squared SDs) are measures of the improved estimation efficiency of the better model. Individual patient-superpixel intercept and slope posterior means, posterior SDs, and 95% CrI lengths from HSL, CAR, and SLR were plotted against each other to compare the inferences from the three approaches. A patient-superpixel profile is the set of measurements over time for a single patient in a single superpixel. We omitted 18 patient-superpixel profiles in the data cleaning step because there were two or more outliers in a single patient-superpixel profile and omitted two more patient-superpixel profiles with constant GCC at all time points, as SLR, unlike the HSL model or CAR, is unable to provide an appropriate inference. Thus 5419 is slightly smaller than 111 × 49 = 5439 patient-superpixel profiles. 

Simulation Study

We ran a simulation study to provide comparisons between the accuracy of HSL, CAR, and SLR estimates for rates of change (slopes); this comparison is only possible in a simulation study where true individual eye slopes are known. In the initial step, a true slope and intercept for each superpixel-patient combination (49 superpixels in each of 111 patient-eyes) was set by sampling each from a normal distribution with mean equal to the sum of the superpixel-specific population and macula-wide patient-specific parameters and variance equal to the superpixel intercept and slope random effects variance. All superpixel-specific and patient-specific parameters were set as the posterior means from our model fit to the actual data; log residual variance was set as the posterior mean of the patient plus superpixel log components; visit times were copied from the data set and visit effect parameters were set to the posterior mean of those parameters. In the next step, we generated 100 data sets with random residual errors and visit effects. 

We recorded bias, estimator variance, and root mean squared error (RMSE) from the simulation for all patient-superpixel intercepts and slopes. We took the posterior means from the HSL and CAR models as the estimates, and the least squares estimate was the SLR estimate. Bias was calculated as the average estimate over simulated datasets minus truth, which is a measure of whether the model estimates the true value on average. Estimator variance is the variance over simulations of the estimate around its average estimate over simulated data sets, a measure of precision. Estimator variance is estimator SD squared. RMSE is the square root of the sum of squared bias and estimator variance  

\begin{eqnarray*} {\rm{RMSE}} &=& \sqrt {\frac{{\mathop \sum \nolimits_{m = 1}^{100} {{\left( {{\rm{estimat}}{{\rm{e}}}_m - {\rm{trut}}{{\rm{h}}}_m} \right)}}^2}}{{100}}}\\ &=& \ \sqrt {{\rm{bia}}{{\rm{s}}}^2 + {\rm{estimator\ variance}}} \ , \end{eqnarray*}

where m indexes data sets. The RMSE is the gold standard for assessing the accuracy of a model for estimating the true parameter value. We calculated average over simulations of 95% Credible Interval Length (CrIL), which is the length of the interval with 2.5% posterior probability to the left and to the right of the interval for HSL and CAR and the usual 95% CI for SLR. Finally, 95% CrI Coverage is the proportion of the time where the 95% CrI contains the truth. For bias, estimator SD, RMSE, and 95% CrIL, smaller is better. For 95% CrI Coverage, values close to 95% are preferred. For the estimator SD, RMSE, and 95% CrIL, we report the SLR/HSL and CAR/HSL ratios; ratios greater than one represent better performance of HSL. We report superpixel summaries by averaging reported measures over all patients at each superpixel location and global summaries by averaging each component over all patients and superpixel locations. 

Results

A total of 39,625 GCC superpixel measurements in 49 superpixels from 815 visits were included in the analysis. Supplementary Table S1 presents clinical and demographic characteristics of the study sample. Mean (SD) baseline 10-2 visual field mean deviation was −8.9 (5.9) dB. Mean (SD) follow-up time was 3.6 (0.4) years with an average (SD) of 7.3 (1.1) OCT scans per eye. 

Table 1.

 

Averages Over all Patient-Superpixel Intercepts and Slopes of RMSE, CrIL, CrI Coverage Probability

 

Simulation Results

Table 1 presents the average RMSE, 95% CrIL, and 95% CrI coverage probability for SLR, CAR, and HSL averaged over all 5439 patient-superpixel for intercepts and for slopes. Compared to SLR, the HSL model had substantially better RMSE (intercept SLR/HSL mean ratio: 1.38; slope SLR/HSL mean ratio: 2.05) and 95% CrIL (intercept SLR/HSL mean ratio: 1.70; slope SLR/HSL mean ratio: 2.25), and similar 95% CrI coverage probability. Compared to CAR, the HSL model had better RMSE (intercept CAR/HSL mean ratio: 1.25; slope CAR/HSL mean ratio: 1.65), similar 95% CrIL, and much better 95% CrI coverage probability (HSL vs. CAR coverage probability for intercepts: 0.949 vs. 0.882; for slopes: 0.947 vs. 0.810). 

Patient-Superpixel Intercepts and Slopes

Supplementary Table S2 presents summaries of the grand mean and tenth to ninetieth percentile range of posterior means, posterior SDs, estimator SDs, 95% CrILs, RMSEs, and 95% CrI coverage probabilities for intercepts and slopes across all 5439 patient-superpixel combinations for HSL, CAR, SLR, and the differences between the means (SLR − HSL and CAR − HSL) and ratios (SLR/HSL and CAR/HSL) of the SDs, CrILs, and RMSEs. 

Table 2.

 

Summary of the Grand Mean of Posterior Means and Posterior SD for Intercepts and Slopes for all 5419 Patient-Superpixel Combinations for SLR, CAR, and HSL Models for the Data From the Patient Cohort

 

Intercepts

For the 5439 individual patient-superpixel intercepts, SLR, CAR, and HSL models demonstrated similar average, tenth-, and ninetieth-percentile posterior means. The SLR model had larger (worse) estimator SDs (ratio range 1.12–2.31) and larger (worse) 95% CrILs (ratio range 1.23–3.09) than HSL for all patient-superpixels. Although CAR had similar posterior SDs and 95% CrILs to HSL, HSL had higher (better) 95% CrI coverage probabilities for 80.2% (4364 out of 5439) of patient-superpixel intercepts. Compared to SLR and CAR, the HSL model had lower (better) RMSE values for 95.1% (5174/5439) and 92.7% (5028/5439) of patient-superpixel intercepts, respectively. 

Slopes

Averaged across all 5439 estimated patient-superpixel slopes, SLR, CAR, and HSL had similar average slope estimates. However, the tenth to ninetieth percentile range of SLR slope posterior means (−1.36, 0.48 µm/year) was much wider than HSL (−1.17, 0.19 µm/year) and CAR (−1.01, 0.14 µm/year). The SLR model had larger (worse) estimator SDs than HSL (SLR/HSL ratio range 1.17–11.84) and 95% CrILs (SLR/HSL ratio range 1.32–6.44) for all patient-superpixel slopes. The CAR model had larger estimator SDs than HSL (CAR/HSL ratio range 0.71–7.57) with comparable 95% CrILs (CAR/HSL ratio range 0.55–2.75); however, HSL had better 95% CrI coverage probabilities for 4965/5439 (92.1%) patient-superpixel slopes. The average 95% CrI coverage probabilities for patient-superpixel slopes were 0.951, 0.810, and 0.949 for SLR, CAR, and HSL, respectively. Compared to SLR and CAR, the HSL model had lower RMSE values for 5211/5439 (95.8%) and 5097/5439 (93.7%) patient-superpixel slopes, respectively. Figure 1 plots histograms of mean squared error (MSE), squared bias, and estimator variance for patient-superpixel slopes for all patients and superpixels for the SLR (blue), CAR (green), and HSL (red) models. The HSL model had noticeably smallest MSE because of having much smaller estimator variance compared to CAR and SLR. Equivalent histograms are plotted for the intercepts in Supplementary Figure S1, where the HSL model also demonstrates noticeably smaller MSE compared to CAR and SLR. 

Figure 1.

 

Histograms of MSE (upper left), squared bias (upper right), and estimator variance (lower left) for patient-superpixel slopes comparing SLR, CAR, and HSL models pooled across all superpixels from the simulation study. Mean squared error is squared bias plus estimator variance; the x-axis values are directly comparable in all three plots. Lower MSE indicates better model performance. Counts of large values omitted due to truncating the x-axis at 1: MSE (HSL = 43; CAR = 241; SLR = 500), squared bias (HSL = 31; CAR = 120; SLR = 0), and estimator variance (HSL = 0; CAR = 8; SLR = 492).

Superpixel-Averages of Intercepts and of Slopes

For the intercepts, there was significant bias averaged across patients for HSL (range −0.054 to 0.082 µm) and CAR (range −0.486 to 0.290 µm) for 28/49 and 45/49 superpixels, respectively, while SLR had no bias in any superpixels. While the magnitude of the observed bias for HSL was not clinically relevant both in comparison to the size of the intercepts (posterior mean range 52.9–101.6 µm) and compared to the posterior SDs, the magnitude of bias for CAR was more comparable to the posterior SDs. In contrast, SLR had substantially larger (worse) estimator SD (ratio range 1.29–1.74), RMSE (ratio range 1.19–1.58), and 95% CrILs (ratio range 1.47–1.94) than HSL for all 49 superpixels. The CAR model had larger estimator SD (ratio range 0.97–1.40) and RMSE (ratio range 1.15–1.36) than HSL across superpixels, but similar 95% CrILs (ratio range 0.76–1.33). There were no important differences (fraction coverage − 0.95) in 95% CrI coverage probability for HSL and SLR, but noticeable differences in parafoveal and nasal superpixels for CAR (HSL range −0.014 to 0.008; CAR range −0.169 to 0.012; SLR range −0.004 to 0.006). Superpixel averages of estimator SD, RMSE, 95% CrIL, and 95% CrI coverage probability for intercepts and slopes are shown in Supplementary Figures S2, S3, S4 and S5, respectively. 

For the slopes, there was significant bias for HSL (range −0.040 to 0.038 µm/year) and CAR (range −0.155 to 0.252 µm/year) for 27/49 and 45/49 superpixels, respectively. The magnitude of the bias for HSL was generally modest in comparison to the slopes (range −0.982 to 0.025µm/year) and was small in comparison to posterior SDs; however, the magnitude of the bias for CAR was on par with the slope estimates. For the 49 superpixel slopes, SLR had larger (worse) estimator SD (ratio range 1.55–4.56), RMSE (ratio range 1.33–3.08), and 95% CrILs (ratio range 1.63–3.14) than HSL across superpixels. CAR had larger estimator SD (ratio range 0.94–3.05) and RMSE (ratio range 1.27–2.28) than HSL across superpixels, but similar 95% CrILs (ratio range 0.73–1.69). There were no substantive differences (coverage − 0.95) in 95% CrI coverage probability for HSL and SLR, but noticeable differences in parafoveal and nasal superpixels for CAR (HSL range −0.033 to 0.024; CAR range −0.273 to −0.006; SLR range −0.004 to 0.005). Figure 2 displays heat maps of the intercept and slope average RMSE ratios (CAR/HSL and SLR/HSL) across superpixels. In all superpixels, HSL outperforms CAR and SLR with all ratios >1. Figure 3 gives the proportion of significant negative and positive slopes detected by HSL, CAR, and SLR when the true slope is negative. The HSL model detected a higher proportion of significant negative slopes than SLR in 38/49 superpixels, with notably larger differences in the central superpixels. It also detected a lower proportion of significant positive slopes than SLR in all superpixels. In contrast, CAR detected a higher proportion of significant negative slopes in 36/49 superpixels, but also a higher proportion of significant positive slopes in all 49 superpixels than HSL when the true slopes were negative. The proportions of significant negative and positive slopes detected by HSL, CAR, and SLR when the true slope is positive are shown in Supplementary Figure S6. The HSL model detected a lower proportion of significant negative slopes in 42/49 superpixels than SLR when the true slopes were positive and a lower proportion of significant positive slopes in 43/49 superpixels than SLR. In addition, HSL detected a lower proportion of significant negative slopes in all 49 superpixels and a lower proportion of significant positive in 39/49 superpixels than CAR. A detailed breakdown of the proportion of significant slopes detected by superpixel is shown in Supplementary Table S3. 

Figure 2.

 

Heat map of the average RMSE ratio for SLR over HSL or CAR over HSL for intercepts (top) and slopes (bottom) by superpixel location from the simulation study. HSL outperforms SLR and CAR in all superpixels as ratios greater than one favor HSL over the alternative model.

Figure 3.

 

Heat map of the percentage of significant negative slopes (top) and significant positive slopes (bottom) detected by SLR, CAR, and HSL models when the true slope is negative in the simulation study.

Analysis of the AGPS Data

Table 2 presents summaries of the grand mean and tenth to ninetieth percentile range of posterior means and SDs for intercepts and slopes across all 5419 patient-superpixel combinations for HSL, CAR, and SLR and the differences between the means (CAR − HSL and SLR − HSL) and ratios of the SDs. Globally, HSL, CAR, and SLR posterior means were similar on average for both the intercepts (CAR − HSL mean difference [tenth, ninetieth percentile]: 0.01 [−1.36, 1.39 µm]; SLR − HSL difference: 0.01 [−1.53, 1.55 µm]) and slopes (CAR − HSL mean difference [tenth, ninetieth percentile]: −0.02 [−0.76, 0.71 µm/year]; SLR − HSL difference: −0.02 [−0.85, 0.82 µm/year]); however, substantial differences existed across patients and superpixels. The striking difference between HSL and SLR was the systematically higher posterior SDs for SLR compared to HSL particularly for the slopes (median 1.69, tenth, ninetieth percentile = 0.93, 3.28). The mean and tenth and ninetieth percentiles of posterior means and SDs of slopes by superpixel for HSL, CAR, and SLR are shown in Supplementary Table S4. Although the posterior means of the slopes are similar, the ranges are noticeably smaller and the posterior SDs are also uniformly smaller across all superpixels for HSL and CAR. Across 5419 patient-superpixel curves, the HSL posterior SDs were smaller than those of SLR for 75.2% of intercepts and 87.1% of slopes but larger than those of CAR for 63.4% of intercepts and 54.6% of slopes. Figure 4 displays the median and tenth and ninetieth percentile of posterior SD ratios of CAR/HSL and SLR/HSL. Across all 49 superpixels, HSL has smaller posterior SDs than SLR; HSL has smaller posterior SDs than CAR in temporal superpixels and larger posterior SDs in nasal superpixels. Figure 5 shows scatter plots of the posterior means from SLR against HSL. There is noticeable shrinkage towards the population mean in the HSL estimates in the peripheral superior and temporal superpixels. The scatter plots of the posterior means from CAR against HSL are shown in Supplementary Figure S7, where there is noticeable shrinkage in the HSL estimates in the temporal regions as well. We provide the scatter plots of posterior means from SLR against HSL and from CAR against HSL on the same axes for all 49 superpixels in Supplementary Figures S8 and S9, respectively. 

Figure 4.

 

Heat map of the median (tenth, ninetieth percentile) slope SD ratio as a function of superpixel location for the patient cohort data. The ratio was defined as either CAR over HSL or SLR over HSL. Ratios >1 indicate better performance for HSL. On average, HSL outperforms SLR at all superpixel locations.

Figure 5.

 

Scatter plots of slope posterior means from simple linear regression model (SLR, y-axis) against those from the hierarchical spatial longitudinal model (HSL, x-axis) in each superpixel for the patient cohort data (right eye format). The SLR posterior means are much more variable than the HSL means. There is noticeable shrinkage towards the population mean in the HSL estimates in the peripheral superior and temporal regions. Each plot is square with its own axes with the x- and y- axes having the same range. The red dashed line represents the x = y diagonal.

The HSL model identified a higher proportion of significant negative slopes compared to SLR and a lower proportion compared to CAR (HSL = 17.6%; CAR = 26.6%; SLR = 15.6%); it detected a lower proportion of significant positive slopes compared to both CAR and SLR (HSL = 1.2%; CAR = 6.9%; SLR = 4.6%). Supplementary Figure S10 presents the McNemar's test results comparing the proportion of significant negative slopes between HSL and SLR or HSL and CAR, where a higher proportion of significant negative slopes was identified in 13/49 superpixels for HSL compared to SLR, 4/49 superpixels for SLR compared to HSL, and 31/49 superpixels for CAR compared to HSL. Supplementary Figure S11 shows the McNemar's test results for comparing the proportion of significant positive slopes, where a lower proportion of significant positive slopes was identified in 21/49 superpixels for HSL compared to SLR and 36/49 for HSL compared to CAR. 

Discussion

We examined the ability of SLR, CAR, and HSL to accurately estimate rates of GCC thinning within macular superpixels in a cohort of eyes with central or moderate to advanced glaucoma damage at baseline. Our novel HSL model resulted in lower posterior SD for both intercepts (SLR/HSL median ratio: 1.35) and slopes (SLR/HSL median ratio: 1.69, Table 1) indicating a marked improvement in the certainty in estimated intercepts and rates of change for HSL over SLR. The simulation study showed a significantly higher performance by HSL compared to SLR in terms of detecting actual change; on simulated models, HSL detected 21% of slopes as significantly negative while SLR detected only 13% when the true slopes were negative. A smaller and still significantly better performance was observed with the cohort data. Although CAR also reduced the posterior SD, CAR estimates were more biased by superpixel location and offered substantially reduced 95% CrI coverage probability (mean intercept probability = 0.882; slope probability = 0.810). Based on the simulation study findings, CAR under-reported posterior SDs or overestimated coverage, implying that in the AGPS data analysis, the more frequent declarations of significance are overly optimistic. Therefore, with the advantage of the reduction in noise and appropriate coverage, HSL could detect changes in GCC more efficiently and earlier with relevant clinical and research implications for earlier detection of glaucoma progression. 

Simple linear regression is still frequently used to estimate global or local rates of change of structural or functional measures in the field of glaucoma. We show that SLR has numerous weaknesses in this context. Data from the patient cohort are not used to help draw inferences about individual eyes, sectors or superpixels. Spatial correlations are ignored, and correlation of baseline thickness and slopes are not accounted for. Visit effects cannot be accommodated in SLR. These limitations in SLR lead to substantially larger uncertainty in estimating individual patient-superpixel rates of change. In contrast, our HSL model accommodates all these features of the data and thus reduces uncertainty in estimating rates of change. 

Our Bayesian hierarchical spatial longitudinal model addresses these shortcomings in SLR by modeling the structure of the data with patients nested in a cohort, spatial correlations across the macula, and visit effects. The HSL model provides intercept and slope estimates with substantially higher accuracy by incorporating information from the cohort and other superpixels from the same person and reduces the posterior variances for both intercepts and slopes. We implemented this updated version of the HSL model developed in our research laboratory in this study. The current version allows random variation across superpixels of individual intercepts and slopes from a global patient estimate. Similarly, intercepts and slopes within individual superpixels are allowed to randomly vary from the population intercept and slope. Population level intercepts and slopes vary across the macula and log residual variances are modeled with superpixel and individual patient components. We recently identified spatial distributions of patient-intercepts, patient-slopes, patient-residual variances, and residuals across superpixels through residual analysis of a prior model that did not include spatial effects except for superpixel level parameters.²⁶ The current model in this study improves on the model in Mohammadzadeh et al.²⁶ by including such spatial effects for patient-superpixel intercepts, slopes, and residual variances. Moreover, we now include visit effects, which model correlations in residuals across superpixels for a single patient-visit. These spatial effects are modeled as factors in a factor analysis with one factor per dimension. 

Visit effects can be estimated with our HSL model but are impossible to implement in SLR. Visit effects model the correlation in regression residuals across superpixels in a single visit. Presence of visit effects means that all GCC thickness measurements from a single patient-visit tend to be randomly above or below the patient's actual GCC thickness. Estimating visit effects requires multivariate modeling of all superpixel thickness measures in one model. We believe this is the first study to identify visit effects in structural OCT data. Visit effects in visual field data are well known,^27–29 though modeling of visual field visit effects has been rarely done.³⁰ 

Our findings indicate a more than twofold reduction in the variance of slope estimates with the HSL model compared to SLR. This significant reduction in posterior variance shows that the HSL model is much more efficient than SLR in using the available data and it allows HSL to detect significant rates of change earlier as compared to SLR. For example, HSL identified a higher proportion of significant negative slopes (17.6% vs. 15.6%) and lower proportion of significant positive slopes (1.2% vs. 4.6%) as compared to SLR. Although these numbers may not seem clinically impressive, the overall superior performance of the HSL is clinically relevant as it not only increases negative hit rates, an indication of higher sensitivity for identifying actual decreasing thicknesses, but it also results in marked reduction of significant positive slopes, an indication of potentially higher specificity as GCC thickness is not expected to increase over time. These findings were confirmed in the simulation study. 

There are previous studies that used hierarchical linear models for detection of longitudinal change in visual fields; our study is unique as it addresses detection of longitudinal changes in macular structural measurements.^{14,16,19,30,31} Montesano et al.¹⁴ developed a Bayesian model for visual field data accounting for the within-eye hierarchical structure, data censoring, and the heteroskedastic variance as a function of the mean threshold sensitivity; they found that time to detect progression was shorter for Bayesian models compared to SLR with permutation analysis of point-wise linear regression.¹⁴ 

In the article by Betz-Stablein et al.,¹⁹ a conditional autoregressive prior was used to model spatial correlation across the visual field intercepts and slopes within a single person but data from multiple patients were not included in a single model. Intercepts and slopes were not modeled as correlated, and conditional on neighboring superpixels, distant superpixels were considered independent. Our CAR model, although inspired by Betz-Stablein et al.,¹⁹ is novel; like in Betz-Stablein et al., our CAR model was fit to data from each patient separately. Although the variances in slope estimates were on par with HSL, there was substantially more bias and reduced 95% CrI coverage probability, translating to a higher proportion of significant slopes identified even when the true slopes were of the opposite sign; thus the CAR model had higher error rates. In contrast, our model also has random patient-superpixel intercepts and slopes, but then allows for finding global patterns of associations in both intercepts and slopes, as we demonstrated in our previous work.²⁶ It also accommodates residual variance varying by both superpixel and patient and allows for visit effects as well. 

Based on the simulation study, HSL has substantially lower RMSE than SLR and CAR for both individual patient intercepts and slopes and for superpixel aggregated effects compared to SLR. The SLR model does not have population or superpixel parameters that summarize the population of patients globally or in a superpixel, so any superpixel level inference from SLR is ad hoc; by taking a Bayesian approach with our SLR models, we were able to create an appropriate, if simplistic, inference for SLR results aggregated across patients within a superpixel. For patient-superpixel intercepts and slopes, HSL had lower RMSE in an overwhelming proportion of intercepts (95%) and slopes (96%). Estimates from HSL had thus lower noise as compared to SLR. 

Our simulation results demonstrated a small bias in some estimated HSL intercepts and slopes. This phenomenon of shrinkage toward the prior mean is a well-known feature of Bayesian estimators. By borrowing information from the population of patients, the Bayesian models shrink all estimates toward population averages. Shrinkage depends on the play between uncertainty in the estimate without population information and the variability in the population. Unusually high or low estimates without shrinkage are typically due to noise in addition to possibly higher or lower underlying true values; these values are shrunken by greater absolute amounts in hierarchical Bayesian models than estimates nearer to the center of the distribution. Posterior means from a hierarchical model are more stable than classical non-Bayesian estimates such as from SLR. This shrinkage mitigates against erroneous high/low slope estimates and, hence, helps prevent making aggressive therapeutic decisions when relying on uncertain and possibly erroneous estimates, pending additional data. The magnitude of the superpixel-averaged slope bias was at most 10% of the average slope across all patients and superpixels (−0.040 µm/year vs. −0.39 µm/year). Despite this bias, the HSL model was still able to identify a higher proportion of worsening slopes while at the same time minimizing the significant positive slopes, which is desirable within both clinical and research frameworks. 

Our model, as with any model, assumes data generated by the model. Because we have normality of residuals built into the model, it does not necessarily model well a distribution with long tails. The obvious outliers are quite extreme, and we tried to remove all the obvious outliers. We believe these are due to measure error. Undoubtedly there are still outliers in the data set. A future extension of our methods may allow for diagnosing outliers via the model itself, or we may develop a novel and improved outlier removal algorithm. 

Optical coherence tomography data are somewhat noisy. If clinicians use raw data with outliers, they will not be using the best data possible, and the conclusions could be flawed. Undoubtedly OCT image quality will improve with time and outliers will lessen or disappear. In the meantime, our outlier removal algorithm does not require model fitting; clinicians could use it in clinical practice before trying to fit any model including SLR. 

In clinical practice, the HSL model can be fit to longitudinal data from the entire cohort. The estimated rates of change for the cohort of patients provide information on progression rates at the patient-superpixel level in the same cohort. Once posterior samples are obtained, observations at future time points for any patient can be predicted using the posterior predictive distribution. For a new patient from the same population, intercepts and slopes can be predicted using the posterior samples without refitting the model. Periodically, as more data is collected over time, the HSL model can be refitted to update model parameters and better guide predictions. For clinicians at a new institution, it would be best to refit the model to data from the new institution and then continue along the same path. 

The implications of our findings go beyond macular structural measures. The proposed framework can be applied to other structural measures such as retinal nerve fiber layer or neuroretinal rim measurements with modifications for differing geometry and scale of measurements given the fact that all the limitations of SLR apply to those structural measures as well. We are also developing a similar hierarchical model for visual field measurements. The proposed framework is a useful starting point for analysis of visual fields; however, modifications will be required given properties of visual fields such as increased variability with worsening thresholds (heteroscedasticity), censoring, and possible loss of information once threshold sensitivity at individual test locations drops to 15 to 19 dB or below.^14,32 

In conclusion, we present a novel Bayesian HSL model that improves estimation accuracy of local GCC rates of change. In a simulation study, SLR and CAR have median MSE ratios over the proposed model of 3.3 and 2.4, respectively, for estimating superpixel-patient slopes; in both the simulation study and in the patient cohort data, HSL identifies a higher proportion of deteriorating superpixels when compared to SLR while minimizing positive detection rates. This efficiency is found by more fully utilizing already available information from measurements on a cohort of glaucoma patients and jointly analyzing measurements on all superpixels. Our findings have important implications for improved detection of glaucoma deterioration both clinically and in the research setting. 

Acknowledgments

Supported by an NIH R01 grant (R01-EY029792), an unrestricted Departmental Grant from Research to Prevent Blindness, and an unrestricted grant from Heidelberg Engineering. 

Presented as a poster at the Annual Meeting of the Association for Research in Vision and Ophthalmology, May 2022, Denver CO. 

Disclosure: E. Su, None; V. Mohammadzadeh, None; M. Mohammadi, None; L. Shi, None; S.K. Law, None; A.L. Coleman, Alcon Foundation and Laboratories Thea SAS (2022-2023) (C); J. Caprioli, None; R.E. Weiss, None; K. Nouri-Mahdavi, Topcon Healthcare (R), Heidelberg Engineering (F)