**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.

^{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.

*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.

^{9}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.

^{14}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.

^{19}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.

^{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.

*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.

^{20}

^{20}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.

^{20}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*denote a single observation of GCC thickness (µm) for patient

_{ijk}*i*at time

*t*in superpixel

_{ij}*k*, the model is

_{0}is the global macula-wide intercept; µ

_{1}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

*i*th patient; β

_{1ik}is the patient-superpixel interaction slope for the

*i*th 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*is the macula-wide patient-specific intercept random effect;

_{0i}*PC*is the macula-wide patient-specific slope random effect;

_{1i}*PC*is the macula-wide patient-specific log residual variance random effect;

_{2i}*VE*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 ϕ

_{ij}_{0k}

*PC*

_{0i}, slope ϕ

_{1k}

*PC*

_{1i}, residual variance \(\sigma _k^2{\rm{exp}}( {P{C}_{2i}} )\), and visit effects ϕ

_{2k}

*VE*. In each component, terms with a

_{ij}*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*are random patient effects, or visit effects, that indicate how much patient

_{ij}*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*.

^{23}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.

*n*is the number of observations for patient

_{i}*i*because SEs are not standard deviations and SEs do not account for the classical estimates being t distributed with

*n*− 2 degrees of freedom. Using the posterior SD from this Bayesian SLR model puts SLR, CAR, and HSL on an equal footing.

_{i}^{19}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

_{0}is the overall eye intercept; α

_{1}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 σ

^{2}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.

*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.

*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

*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.

^{26}The current model in this study improves on the model in Mohammadzadeh et al.

^{26}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.

^{27}

^{–}

^{29}though modeling of visual field visit effects has been rarely done.

^{30}

^{14}

^{,}

^{16}

^{,}

^{19}

^{,}

^{30}

^{,}

^{31}Montesano et al.

^{14}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.

^{14}

^{19}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.,

^{19}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.

^{26}It also accommodates residual variance varying by both superpixel and patient and allows for visit effects as well.

^{14}

^{,}

^{32}

**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)

*Transl Vis Sci Technol*. 2020; 9(2): 1–19. [CrossRef]

*Med Image Comput Comput Assist Interv*. 2013; 16(Pt 2): 444–451. [PubMed]

*PLoS One*. 2020; 15(7): 1–15. [CrossRef]

*Am J Ophthalmol*. 2022; 237: 71–82. [CrossRef] [PubMed]

*Arch Ophthalmol*. 2002; 120: 701–713; discussion 829–830. [CrossRef] [PubMed]

*Arch Ophthalmol*. 2002; 120: 714–720; discussion 829–830. [CrossRef] [PubMed]

*Eye*. 2013; 27: 803–808. [CrossRef] [PubMed]

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

*Transl Vis Sci Technol*. 2020; 9(7): 1–14. [CrossRef]

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

*Invest Ophthalmol Vis Sci*. 2012; 53: 6939–6946. [CrossRef] [PubMed]

*Eye*. 2011; 25: 269–277. [CrossRef] [PubMed]

*PLoS One*. 2014; 9(8): e105611. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2021; 10(12): 1–16. [CrossRef]

*Transl Vis Sci Technol*. 2022; 11(2): 16. [CrossRef] [PubMed]

*PLoS One*. 2014; 9(1): 1–11.

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

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

*Invest Ophthalmol Vis Sci*. 2013; 54: 1544–1553. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2021; 10(4): 1–13. [CrossRef]

*R Foundation for Statistical Computing, Vienna, Austria.*2023.

*An Easy Guide to Factor Analysis*. Oxfordshire, UK: Routledge; 2014.

*J Comput Graph Stat*. 2017; 26: 403–413. [CrossRef]

*Ophthalmol Sci*. 2022; 2(3): 100187. [CrossRef] [PubMed]

*Transl Vis Sci Technol*. 2021; 10(13): 1–17. [CrossRef]

*Br J Ophthalmol*. 2008; 92: 569–573. [CrossRef] [PubMed]

*Sci Rep*. 2021; 11(1): 1–9. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2015; 56: 4283–4289. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2000; 41: 417–421. [PubMed]

*Invest Ophthalmol Vis Sci*. 2012; 53: 5985–5990. [CrossRef] [PubMed]