**Purpose**:
Develop a hierarchical longitudinal regression model for estimating local rates of change of macular ganglion cell complex (GCC) measurements with optical coherence tomography (OCT).

**Methods**:
We enrolled 112 eyes with four or more macular OCT images and ≥2 years of follow-up. GCC thickness measurements within central 6 × 6 superpixels were extracted from macular volume scans. We fit data from each superpixel separately with several hierarchical Bayesian random-effects models. Models were compared with the Watanabe–Akaike information criterion. For our preferred model, we estimated population and individual slopes and intercepts (baseline thickness) and their correlation.

**Results**:
Mean (SD) follow-up time and median (interquartile range) baseline 24-2 visual field mean deviation were 3.6 (0.4) years and −6.8 (−12.2 to −4.3) dB, respectively. The random intercepts and slopes model with random residual variance was the preferred model. While more individual and population negative slopes were observed in the paracentral and papillomacular superpixels, superpixels in the superotemporal and inferior regions displayed the highest correlation between baseline thickness and rates of change (*r* = –0.43 to –0.50 for the top five correlations).

**Conclusions**:
A Bayesian linear hierarchical model with random intercepts/slopes and random variances is an optimal initial model for estimating GCC slopes at population and individual levels. This novel model is an efficient method for estimating macular rates of change and probability of glaucoma progression locally.

**Translational Relevance**:
The proposed Bayesian hierarchical model can be applied to various macular outcomes from different OCT devices and to superpixels of variable sizes to estimate local rates of change and progression probability.

^{1}

^{,}

^{2}Detection of glaucoma progression is essential for preventing irreversible visual loss. Earlier detection of glaucoma worsening would allow clinicians to titrate treatment in a timely fashion. Rates of progression vary among patients. Identifying eyes that are progressing rapidly is important for the appropriate allocation of available resources. In some glaucoma eyes, progressive structural changes can be detected earlier than functional changes.

^{3}

^{–}

^{8}

^{9}

^{–}

^{14}There is now evidence supporting the utility of macular OCT imaging for detection of glaucoma progression, especially in later stages of glaucoma.

^{2}

^{,}

^{10}

^{,}

^{15}

^{–}

^{17}Most prior studies of macular thickness changes over time have used global outcome measures or thickness measurements in hemiretinal regions. These approaches may fail to identify local or regional progressive thinning of the inner macular layers. Hence, estimating macular thickness changes in smaller areas such as 3° × 3° superpixels is of interest.

^{4}

^{,}

^{18}

^{–}

^{23}Miraftabi et al.

^{24}reported that repeated thickness measurements across the macula at the level of 3° × 3° superpixels demonstrated a low and uniform within-session variability. Nouri-Mahdavi et al.

^{25}recently showed that between-session variability of local inner retinal thickness measurements was very low and that most of the total between-session variability could be explained by within-session variability. While measurement of change in localized areas of the macula is essential for detecting local changes, there is a significant correlation between thickness measurements and thickness changes in physically or physiologically adjacent regions of the macula.

^{20}A hierarchical approach, once fully developed, can address the varying behavior of superpixels across the macular region while at the same time taking into account any potential correlations in rates of change or their variance among superpixels of varying sizes.

^{26}

^{–}

^{30}Most such studies addressed functional measurements in glaucoma.

^{31}

^{–}

^{40}Random-effects models can more accurately estimate rates of structural change if appropriately implemented as compared to fixed-effects models that do not take population information into account.

^{41}

^{–}

^{45}Random-effects models are known to be more efficient and accurate in modeling subject-specific effects in longitudinal data, allowing population information to help estimate individual trends.

^{46}Random-effects longitudinal models, therefore, would be superior models for analyzing macular thickness measurements over time, both to accommodate within-subject correlations over time and to take advantage of such correlations to more efficiently identify progression.

*P*< 0.05 on the pattern deviation plot or VF mean deviation (MD) worse than −6 dB. Exclusion criteria were baseline age less than 40 years or greater than 80 years, best-corrected visual acuity <20/50, refractive error exceeding 8 diopters (D) of sphere or 3 D of cylinder, or any significant retinal or neurologic disease potentially affecting OCT measurements.

*y*is the GCC thickness observed on subject

_{ijk}*i*at time

*t*at subject

_{ij}*i*'s

*j*th visit, in superpixel

*k*, where the time of the first visit for all subjects is

*t*= 0. A profile is the time-ordered sequence of observations (

_{ij}*t*,

_{ij}*y*) from

_{ijk}*j*= 1 to

*n*, where

_{ik}*n*is the number of observations on superpixel

_{ik}*k*for subject

*i*. Profiles were plotted by connecting consecutive observations with line segments. Plotting data helps in identifying appropriate classes of models for fitting data. We plotted profiles for all subjects and superpixels to look for outliers, check for the time trend, and understand variability within a person and between people. Profile plots suggested that outliers were present.

*y*−

_{ijk}*y*

_{i(j − 1)k}| and consecutive visit absolute centered slopes, | (

*y*−

_{ijk}*y*

_{i(j − 1)k})/(

*t*−

_{ij}*t*

_{i(j − 1)}) + 0.5|. Slopes were centered around –0.5 µm/year as that was the mean of the pooled set of slopes across all subjects, superpixels, and pairs of consecutive visits. We plotted histograms of the differences and slopes and examined the largest values. Large absolute differences and large absolute centered slopes are caused by unusual GCC measurements that are high or low compared to preceding or subsequent visit measurements. We flagged absolute slopes greater than 24 µm/year with differences greater than 5 µm between sequential visits; this identified the

*ij*th and

*i*(

*j*− 1)th observations as candidates for removal. We calculated the sum of absolute visit differences in each superpixel for each subject,

*n*observations for subject

_{ik}*i*in superpixel

*k*. Using empirical residuals in profile plots makes it easier to understand time trends as compared to plotting raw data profiles.

^{45}We calculated averages and standard errors of the empirical residuals

*e*in 6-month windows and plotted the averages ±2 standard errors against time to understand the population time trend for each superpixel.

_{ijk}*k*separately with a set of hierarchical random-effects models. Our initially preferred model was linear in time, with a linear population time trend α

_{0k}+ α

_{1k}t

_{ij}and random intercepts β

_{0ik}and random slopes β

_{1ik}(RIAS),

*k*,

- α
_{0k}is the population intercept, the average intercept at time t_{ij}= 0; - α
_{1k}is the population average slope; - β
_{0ik}is the*i*th subject's random intercept: the unknown difference between subject*i*'s intercept and the population intercept α_{0k}; - β
_{1ik}is the*i*th subject's random slope: the unknown difference between subject*i*'s slope and α_{1k}; and - ε
_{ijk}is a random error at time t_{ij}.

_{1k}= β

_{1ik}= 0; a random intercept and fixed slope (RIFS) model, where β

_{1ik}= 0 but α

_{1k}is free to vary; and a random intercept and slope and fixed quadratic (RIASFQ) model, which adds a term \({{\rm{\alpha }}_{2{\rm{k}}}}t_{ij}^2\) to the right-hand side of Equation (1). This is a total of four models for the population and individual time trends: RI, RIFS, RIAS, and RIASFQ. We modeled residuals ε

_{ijk}as normal with unknown variance \(\sigma _k^2\),

_{ik}= β

_{i0k}in the RI models or β

_{ik}= (β

_{i0k}, β

_{i1k})′ in the RIAS models, as normal \({\beta _{ik}}|D\;\sim\;N( {0,{D_k}} )\), where

*D*is scalar for the RI models and

_{k}*D*is a 2 × 2 matrix for the RIAS models:

_{k}_{ik}were inverse gamma with unknown mean σ

_{mk}and unknown standard deviation σ

_{sk}. For the hyperparameters α

_{0k}, α

_{1k}, σ

^{2},

*D*, ρ

_{k}_{k}, and (σ

_{mk},σ

_{sk}), we used convenient prior distributions. Parameters of these distributions were specified so that they worked well and were appropriate for all superpixels. While this approach uses the cohort data to specify these priors, data from each superpixel represent only 1/36th of the entire superpixel data, so very little of each superpixel data contributed to the prior construction. Full model specifications and exact parameter values of the priors are given in the online supplement. We plotted summary measures of our analyses over the 6 × 6 grid of superpixels.

*P*values defined as the posterior probability that the population slope is greater than zero. This is a measure of how well the sign of the population slope is determined and corresponds asymptotically to a one-sided frequentist

*P*value. This number can be doubled to construct a two-sided

*P*value but then loses its Bayesian interpretation. Bayesian

*P*values that are near zero or 1 indicate a slope whose sign is well determined; we used

*P*< 0.025 and

*P*> 0.975 as traditional cutoffs to identify statistical significance. We discuss this in more detail in the online supplement.

_{1k}+ β

_{1ik}, and we calculated the posterior probability

*p*=

_{ik}*P*(α

_{1k}+ β

_{1ik}> 0|

*Y*) that the slopes were positive. Small values of

*p*< 0.025 identify subjects whose GCC thickness is decreasing at a significant rate over time. We report the fraction of eyes that are decreasing in thickness, determined as the fraction of eyes whose slopes α

_{ik}_{1k}+ β

_{1ik}are significantly negative or significantly positive.

^{47}

^{,}

^{48}We also explored an alternative criterion called the widely applicable Bayesian information criterion (WBIC) and compared the results to WAIC for model selection. WBIC led to similar results as WAIC, and results for WBIC are provided in the online supplement.

*Y*represents all the data for the

_{k}*k*th superpixel.

^{45}

*P*value of the population slopes, the posterior mean of the correlation between slopes and intercepts and Bayesian

*P*value of the correlation, and the proportion of individual subject slopes that were significantly negative. The largest population slopes were observed in superpixels 4.3 to 4.7, 5.6, 5.7, 6.5 and 6.6. Superpixels whose posterior mean population slope was below –0.2 microns per year (meaning decreasing more rapidly than 0.2 microns per year) had slopes that were significantly negative (one-sided Bayesian

*P*< 0.01 and mostly

*P*< 0.001). The fastest thinning (most negative) population superpixels were mostly located in the paracentral and nasal (papillomacular) region. The parafoveal superpixel population slopes were most likely to be significantly less than zero. Point estimates of the correlations between random intercepts and random slopes were negative in 34 of 36 superpixels and were significantly negative in 20 of the 36 superpixels; negative correlations stronger than −0.25 were always significant and correlations weaker than −0.22 were never significant. Superpixels in the superotemporal and inferior regions tended to demonstrate the highest correlation between the baseline thickness (the intercept from the hierarchical model) and their rates of change. The top five coefficients for the correlation between random slopes and baseline thickness at the level of superpixels ranged from –0.43 to –0.50.

*P*< 0.025 are given in the last plot of Figure 3. Superpixels 4.3 to 4.6, 5.6 to 5.7, and 6.5 to 6.6 had at least 30% of individual eyes that were identified as significantly decreasing (Bayesian

*P*< 0.025).

^{4}

^{,}

^{20}and thus we used GCC as the macular outcome of interest for this study.

*P*< 0.025) using 4 years of data from 112 subjects. Our Bayesian hierarchical random-effects model is helpful for identifying individual eyes that are thinning at a significant rate because it uses population-level information to help estimate slopes in individual eyes, increasing our chances of identifying individual eyes as decreasing compared to a simple linear regression model analyzing a single individual's observations. The Bayesian approach provides a justified and informative prior for subject-level slopes; the prior, in turn, yields additional information about a given slope above and beyond the information in an individual's observations.

^{16}

^{,}

^{35}

^{,}

^{38}

^{,}

^{40}Our current work represents an initial effort toward developing an optimized hierarchical model for fitting longitudinal structural macular data. Superpixels or even clusters of superpixels across the macula are correlated, with the correlation being a function of physical or physiologic distance between the superpixels.

^{20}We analyzed data for each superpixel separately. This made the modeling easier, at a cost of not allowing us to borrow strength from neighboring superpixels to better estimate trends at each superpixel and preventing us from making joint statements about individual trends across superpixels. Our results are valid despite the lack of a spatial component: posterior means of parameters estimate the corresponding underlying parameters, and standard deviations of the parameters are appropriate summary measures of the uncertainty. For example, posterior estimates, standard deviations, and 95% confidence interval of subject-specific slope estimates could be used as inputs to treatment decisions. Combined spatial-longitudinal models come with high computational costs and have yet to be developed for macular thickness measurements.

^{51}

_{ijk}might demonstrate autocorrelation over time, although we did not see that in our analysis. Recent research has suggested that the superpixels in which measurements behave similarly are not just neighbors but may belong to differently shaped nonconvex clusters.

^{20}Similarly, random intercepts, random slopes, and random errors may have complex spatial associations not easily modeled by local spatial autocorrelation models. Our future work will consider spatial location and physiologic relationships of superpixels (clusters) in expanded Bayesian hierarchical models.

*have been reached before deciding on a change in treatment. Rather, one might use a more generous cutoff of, for example, 0.1 instead of 0.025, for flagging a clinically significant thinning of the GCC in an individual eye.*

_{ik}**V. Mohammadzadeh**, None;

**E. Su**, None;

**S. Heydar Zadeh**, None;

**S.K. Law**, None;

**A.L. Coleman**, None;

**J. Caprioli**, None;

**R.E. Weiss**, None;

**K. Nouri-Mahdavi**, Heidelberg Engineering

*Am J Ophthalmol*. 1989; 107(5): 453–464. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2017; 175: 37–44. [CrossRef] [PubMed]

*Ophthalmol Glaucoma*. 2019; 2(1): 36–46. [CrossRef] [PubMed]

*Ophthalmology*. 2020; 127(7): 888–900. [CrossRef] [PubMed]

*Ophthalmology*. 2005; 112(3): 366–375. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 1999; 40(10): 2242–2250. [PubMed]

*Arch Ophthalmol*. 1991; 109(1): 77–83. [CrossRef] [PubMed]

*Prog Retin Eye Res*. 2005; 24(3): 333–354. [CrossRef] [PubMed]

*Surv Ophthalmol*. 2020; 65(6): 597–638. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2013; 54(3): 1941–1949. [CrossRef] [PubMed]

*J Glaucoma*. 2014; 23(4): 195–198. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2018; 7(1): 20. [CrossRef] [PubMed]

*J Glaucoma*. 2017; 26(3): 208–215. [CrossRef] [PubMed]

*J Ophthalmol*. 2018; 2018: 6581846. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2012; 53(7): 3817–3826. [CrossRef] [PubMed]

*Ophthalmology*. 2017; 124(12S): S57–S65. [CrossRef] [PubMed]

*Ophthalmology*. 2012; 119(2): 308–313. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2016; 57(11): 4815–4823. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2017; 101(8): 1052–1058. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2020; 9(7):50.

*Prog Retin Eye Res*. 2013; 32: 1–21. [CrossRef] [PubMed]

*JAMA Ophthalmol*. 2015; 133(12): 1438–1444. [CrossRef] [PubMed]

*Arch Ophthalmol*. 2011; 129(12): 1529–1536. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2016; 5(4): 5. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2019; 207: 18–36. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2021; 222: 238–247. [CrossRef] [PubMed]

*Sci Rep*. 2020; 10(1): 2808. [CrossRef] [PubMed]

*Ophthalmol Glaucoma*. 2019; 2(2): 72–77. [CrossRef] [PubMed]

*J Glaucoma*. 2017; 26(8): 721–725. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2016; 166: 29–36. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2018; 59(5): 1897–1904. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2015; 56(3): 1603–1608. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2014; 55(12): 8386–8392. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2013; 54(3): 2198–2206. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2012; 53(6): 2760–2769. [CrossRef] [PubMed]

*J Glaucoma*. 2012; 21(3): 147–154. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2012; 53(4): 2199–2207. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2012; 153(6): 1197–1205.e1191. [CrossRef] [PubMed]

*Ophthalmology*. 2012; 119(3): 458–467. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2011; 52(8): 5794–5803. [CrossRef] [PubMed]

*Commun Stat Simul Comp*. 2019: 1–25.

*Stat Med*. 2017; 36(11): 1735–1753. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2010; 149(2): 187–188.e181. [CrossRef] [PubMed]

*Applied longitudinal analysis*. New York, NY: John Wiley & Sons; 2012.

*Modeling Longitudinal Data*. Springer Science & Business Media

*;*2005.

*Ann Stat*. 1976; 1: 384–395.

*Stat Comput*. 2013; 24(6): 997–1016. [CrossRef]

*J Machine Learn Res*. 2010; 11(12): 3571–3594.

*R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing; 2013, https://www.R-project.org/.

*J Stat Soft*. 2016; 71(9): 1–25. [CrossRef]

*JAMA Ophthalmol*. 2016; 134(5): 496–502. [CrossRef] [PubMed]