Abstract
Background
Glaucoma is a progressive optic neuropathy in which longitudinal visual field (VF) prediction is essential for timely clinical decision‐making and preventing irreversible VF loss. However, accurate prediction remains challenging due to the noisy, high‐dimensional nature of VF data and their complex longitudinal patterns. Capturing these patterns may require increasingly sophisticated prediction models, but greater model complexity can reduce interpretability. Therefore, prediction models need to capture temporal patterns while balancing interpretability and computational feasibility. This study is aimed at evaluating and comparing the performance of four methods, including personalized models (first‐order autoregressive [AR(1)] and ordinary linear regression (OLR)), a population‐level long short‐term memory (LSTM), and a naive baseline (BL) forecast for glaucoma VF prediction.
Materials and Methods
Data were simulated for patients with early‐stage glaucoma using a linear mixed effects model (LMM) to account for noise in VF measurements and the impact of individual effects on glaucoma progression. Future VFs were predicted from examinations 6–11 based on previous examinations. Prediction accuracy was assessed using root mean squared error (RMSE) and pointwise mean absolute error (PMAE).
Results
Evaluation of the four approaches showed that OLR and naive BL forecast achieved higher errors than AR(1) and LSTM models. LSTM performed best when the number of examinations was small (in the initial examinations, indeed) in slow, moderate, and fast progression rate scenarios. However, with a larger number of examinations, prediction errors increased for LSTM compared with AR(1), resulting in higher accuracy for AR(1), particularly under moderate and fast progression rate scenarios.
Conclusion
The performance of the models shows dependence on the number of examinations in the present study. These findings highlight the importance of temporal information in longitudinal VF prediction.
Keywords: autoregressive model, glaucoma, LSTM, OLR, visual field prediction
1. Introduction
Glaucoma is recognized as a leading cause of irreversible blindness worldwide [1–3]. The disease is characterized by retinal ganglion cell (RGC) death and functional visual field (VF) defects [4]. Standard automated perimetry (SAP) is employed as a primary functional test for the detection and monitoring of glaucoma [5, 6]. Early and reliable detection of VF progression in glaucoma is critical since VF deterioration is permanent [7].
Quantitative assessments for monitoring glaucoma progression are generally based on either global VF indices or pointwise VF measurements. Global indices provide summary measures of VF function but may overlook localized defects. In contrast, pointwise analyses are more sensitive to localized functional changes; however, this higher sensitivity is accompanied by lower specificity [8, 9]. Among pointwise measurements, total deviation values (TDVs) are particularly suitable for longitudinal modeling because they quantify localized sensitivity loss after adjustment for age‐related decline, enabling standardized comparisons across patients and over time.
A wide range of statistical and machine learning approaches has been developed to predict future VF using either global indices or pointwise measurements [10–23]. Among conventional statistical methods, pointwise linear regression (PLR), which fits ordinary linear regression (OLR) independently at each VF location, has been widely used for longitudinal VF analysis because of its simplicity and clinical interpretability [13, 14, 16, 17, 19, 20, 24]. However, OLR does not account for the temporal dependence among repeated VF measurements within each eye. Time series models such as the first‐order autoregressive (AR(1)) model explicitly accounts for this temporal dependence and, therefore, represents an appropriate alternative for longitudinal VF prediction. Despite this theoretical advantage, its application in glaucoma VF prediction has been limited, with only a few previous studies evaluating the AR(1) model for longitudinal VF measurements [24, 25].
More recently, deep learning techniques, particularly recurrent neural networks (RNNs) and long short‐term memory (LSTM) networks have gained attention because of their ability to learn complex temporal patterns directly from sequential VF data [10, 12] .Previous studies have demonstrated that LSTM‐based and hybrid deep learning architectures outperformed conventional regression approaches in VF prediction [10–12, 21, 23]. Machine learning methods have also shown effectiveness in other aspects of VF analysis, including glaucoma detection based on TDVs and pattern deviation values (PDVs), further highlighting their ability to capture complex VF patterns [15]. However, the successful implementation of these approaches requires large, high‐quality datasets and substantial computational resources. Moreover, comparisons among prediction models have predominantly relied on real‐world datasets, in which reliable assessment of VF progression is challenging due to the substantial variability of SAP measurements and the absence of a universally accepted gold standard for defining true progression [26–32]. These limitations make the assessment and comparison of predictive models difficult [33]. Therefore, differences in predictive performance may reflect variations in dataset characteristics rather than genuine methodological superiority.
Computer simulation offers a controlled and reproducible framework in which progression rates, follow‐up schedules, baseline (BL) disease severity, and measurement variability can be specified a priori, and large‐scale longitudinal VF datasets can be generated. These capabilities facilitate objective evaluation and improvement of progression detection and prediction methods and reduce the limitations associated with real‐world data [9, 33, 34]. Early simulation studies [18] generated longitudinal pointwise threshold sensitivities using a simple linear regression model. Although this approach enabled controlled simulation of localized VF deterioration, it did not explicitly model the dependence among repeated measurements within the same eye or the heterogeneity in BL status and progression rates across eyes. More recently, Sabouri et al. [35] introduced a simulation framework based on linear mixed‐effects models (LMMs) for generating longitudinal VF series.
This LMM‐based simulation accounted for between‐eye heterogeneity and dependence among repeated measurements within each eye through random intercepts and random slopes [36, 37]. Although this framework was developed for global MD, pointwise VF series share a similar longitudinal repeated‐measures structure. LMMs provide a controlled longitudinal data‐generating framework that can incorporate systematic progression, between‐eye heterogeneity, within‐eye dependence, and measurement variability. Compared with simple linear regression‐based simulation [18], this approach, therefore, provides a more realistic representation of the longitudinal structure of pointwise VF data.
Despite these advances, several important gaps remain. Previous simulation studies have focused either on global VF indices or on pointwise simulation using simple linear regression. To our knowledge, no existing framework has applied LMM to simulate longitudinal pointwise VFs while accounting for within‐eye dependence and between‐eye heterogeneity. Moreover, comparative evaluations of OLR, AR(1), and deep learning approaches have mainly relied on real‐world datasets, where the true progression trajectory is unknown. Finally, it is unclear whether the increased complexity of deep learning methods translates into superior predictive performance or whether simpler statistical approaches may perform comparably in specific clinical settings.
Accordingly, this study pursued three primary objectives: (1) to develop an LMM‐based framework for simulating longitudinal pointwise TDVs in early‐stage glaucoma, (2) to compare the predictive performance of conventional statistical models (OLR and AR(1)) and a deep learning model (LSTM), and (3) to examine the influence of disease progression rate and the number of examinations on the models′ performance. In addition, a naive BL forecast, assuming the subsequent observation is identical to the most recent measurement, was included as a benchmark for evaluating the more sophisticated prediction models. Based on these findings, we provide recommendations regarding the appropriate application of OLR, AR(1), and LSTM methods for longitudinal VF prediction in glaucoma.
2. Materials and Methods
2.1. Simulation of Clinical Data
We simulated threshold sensitivity data for 36,000 eyes with a focal inferior nasal defect at BL examination. The simulation period was 5 years with 6‐month intervals, comprising 11 examinations per eye and a total of 396,000 simulated VFs. BL age was set at 40–50 years. Simulation parameters were defined based on the opinion of a glaucoma specialist, who assessed reliable BL VFs with false negative rates < 12% and false positive rates < 4%, as well as fixation losses < 15%. Pointwise VF data were derived from SAP performed with the Humphrey Field Analyzer (Carl Zeiss Meditec), using the 24‐2 test pattern and the SITA Standard strategy. For all analyses and simulations, 52 test locations were used after excluding the two blind‐spot locations from the 24‐2 grid. A BL VF pattern from a glaucoma patient with a focal inferior nasal defect was selected by the glaucoma specialist as the reference pattern for simulation, and age‐matched normative VF data were used to determine the expected interval of threshold sensitivities at each test location. These intervals were used to initialize the BL threshold sensitivity values in the simulated VF data. For the simulation, at locations without glaucomatous damage in the selected BL VF pattern, the lower bound of each interval corresponded to the amount of sensitivity reduction that led to a change in the location from normal to a total deviation probability level of 5%, and the upper bound represented the expected normal sensitivity at that location. For the three locations related to the inferior nasal step defect (locations 27, 28, and 35; Figure 1), sensitivity levels that led to changes in total deviation probability level of 2% and 1% were also determined.
Figure 1.

Numbered visual field map (excluding the two blind spots).
Six progression patterns (VF deterioration over time) and three progression rates were defined for the simulation of clinical data in different scenarios, yielding 18 separately simulated datasets, each containing 2000 eyes (Table 1 and Figure 2). Threshold sensitivity values at each of the 52 test locations were generated using an LMM with random intercept and slope parameters. The presence of random intercept and slope in the LMM causes each eye to have its own unique trajectory at BL and over time [37–39]. The equation of this model for threshold sensitivity values over time is written as follows:
| (1) |
Table 1.
Simulation scenarios based on progression rates and progression patterns.
| Progression pattern | Progression rate | Number of deteriorating locations over time |
|---|---|---|
| Superior nasal defect (small defect) | −0.25 dB/year (slow) | Deterioration in three new locations with no change in the baseline defect (Figure 2a) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) | ||
| Inferior nasal defect (small defect) | −0.25 dB/year (slow) | Deterioration in three new locations with no change in the baseline defect (Figure 2b) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) | ||
| Deepening inferior nasal defect and early arcuate damage (medium defect) | −0.25 dB/year (slow) | Deterioration in eight locations: a deepening of the defect at baseline along with five new locations (Figure 2c) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) | ||
| Deepening inferior nasal defect and early arcuate damage (medium defect) | −0.25 dB/year (slow) | Deterioration in eight new locations with no change in the baseline defect (Figure 2d) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) | ||
| Inferior nasal defect and advanced inferior arcuate damage (large defect) | −0.25 dB/year (slow) | Deterioration in 16 locations: a deepening of the defect at baseline along with 13 new locations (Figure 2e) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) | ||
| Inferior nasal defect and inferior arcuate damage with central threatening, and early superior nasal defect (large defect) | −0.25 dB/year (slow) | Deterioration in 19 new locations with no change in the baseline defect (Figure 2f) |
| 0.75 dB/year (moderate) | ||
| −1.5 dB/year (fast) |
Figure 2.

Simulated progression patterns (visual field deterioration over time). Six progression patterns ((a)–(f)) were applied on a baseline visual field with a focal inferior nasal defect in early‐stage glaucoma with different progression rates.
where T S i j k indicates the threshold sensitivity (dB) for ith eye, jth visit, and kth location. β 0k and β 1k are the fixed intercept and slope parameters, respectively; b 0i k and b 1i k are the corresponding random intercept and random slope effects; t i j denotes time; and ε i j k represents the error term or noise.
The fixed intercept parameter (β 0k ) is the population mean BL threshold sensitivity for each location, whereas the random intercept (b0ik) represents eye‐specific deviations from this population mean. To determine the fixed and random intercepts, each VF test location was assigned a BL threshold sensitivity interval defined by a glaucoma specialist. Likewise, the fixed slope (β 1k ) is the population progression rate, and the random slope parameter (b 1i k ) captures eye‐specific variation around this population progression rate. To determine the fixed and random slopes for deteriorating locations, we defined three intervals for the progression rate based on previous studies [35, 40]: slow (0 to −0.5 dB/year), moderate (−0.5 to −1.0 dB/year), and fast (−1.0 to −2.0 dB/year), with corresponding population means of −0.25, −0.75, and −1.50 dB/year, respectively. Locations assumed to remain stable over time, regardless of whether they were healthy or damaged at BL, were assigned an age‐related decline interval of −0.15 to −0.05 dB/year with a population mean of −0.1 dB/year. The random intercept and random slope parameters were generated from a continuous uniform distribution within the predefined intervals, whereas the fixed intercept and slope parameters were set to the corresponding interval means. The error term, or noise component, represents the variability of repeated measurements for each eye. To add noise to the true threshold sensitivity values, a normal distribution was used based on the study by Rabiolo et al. [18]. At each test location, the noisy threshold sensitivity value was obtained by adding a random Monte Carlo value generated from this distribution to reproduce the inherent variability in VF measurements. The mean of the distribution at each test location was considered the simulated true threshold sensitivity, and the standard deviation (SD) was calculated using the function proposed by Gardiner and Crabb [41] as ln (SD) = −0.081 ×sensitivity (dB) +3.27. An eccentricity‐weighting factor, proposed by Spry et al. [9], was added to the SD to increase the variability of peripheral test locations. It is defined as follows:
| (2) |
where i and j denote the coordinates of each location in the 10 × 10 grids. Noisy threshold sensitivity values simulated outside the range of 0–36 dB were truncated to 0 and 36, respectively. Finally, TDVs and PDVs were generated for each simulated examination using the “VisualFields” package in RStudio. Table 2 summarizes the simulation setup.
Table 2.
Parameters of the linear mixed effects model used to simulate threshold sensitivity in patients with early‐stage glaucoma.
| Parameter | Simulation setup |
|---|---|
| Fixed intercept (β 0k ) | Mean of the expected interval for threshold sensitivity defined by glaucoma specialist at each of the 52 locations. |
| Random intercept (b 0i k ) | Continuous uniform distribution with parameters −1/2 interval distance of threshold sensitivity defined by glaucoma specialist and +1/2 interval distance of threshold sensitivity defined by glaucoma specialist. |
| Fixed slope (β 1k ) | Deteriorating locations: −0.25 dB/year (slow), −0.75 dB/year (moderate), and −1.5 dB/year (fast) |
| Stable locations (healthy or damaged at baseline): −0.1 dB/year (age‐related decline) | |
| Random slope (b 1i k ) | Deteriorating locations: continuous uniform distribution with parameters −0.25 dB/year and +0.25 dB/year (slow and moderate); continuous uniform distribution with parameters −0.50 dB/year and +0.50 dB/year (fast) |
| Stable locations (healthy or damaged at baseline): continuous uniform distribution with parameters −0.05 dB/year and +0.05 dB/year (age‐related decline) | |
| Error term (ε i j k ) | Normal distribution with mean equal to the simulated true threshold sensitivity value and SD equal to exp(−0.081 × sensitivity (dB) +3.27)) |
2.2. Methods to Predict the VF
For each progression rate, datasets from the six predefined progression patterns were combined into a single dataset, resulting in three simulated datasets representing slow, moderate, and fast progression rates (Figure 3). Four methods were employed to predict future VF examinations (6–11) for each of the three simulated datasets, including a naive BL forecast, OLR, AR(1), and LSTM‐based RNN. For each method, separate models were fitted to predict the successive examinations using all preceding VF examinations as inputs. For example, the model predicting the sixth examination used the BL and the second to fifth examinations as inputs, whereas the model predicting the seventh examination used the BL and the second to sixth examinations. The same sequential prediction framework was applied to all subsequent examinations.
Figure 3.

Flowchart of the data preparation process. Eighteen simulated datasets (each containing 2000 eyes) were generated from combinations of six progression patterns and three progression rates. These datasets were then aggregated according to progression rate (slow, moderate, and fast), resulting in three final datasets, each containing 12,000 eyes.
For all analyses, TDVs were used as the primary data representation because they capture age‐adjusted, pointwise measurements of VF sensitivity across the 52 test locations, allowing the models to learn spatial patterns of disease progression. Although threshold sensitivity values were simulated during data generation to derive longitudinal TDVs, they were not used directly for model training or prediction. Similarly, global indices such as MD and pattern standard deviation (PSD) were not used because they summarize the VF into a single value and, therefore, do not preserve the location‐specific information required for predicting future VF changes. Because TDVs were referenced to age‐matched normative values and the same TDV representation and normalization procedure was applied across all models, age‐related adjustment did not favor any prediction method. Consequently, differences in performance reflect the models′ ability to capture temporal and spatial dependencies. Before implementing the models, the simulated data were normalized to the range of [0, 1] using min‐max scaling. Each dataset (12,000 eyes) was randomly divided into training (10,200 eyes; 85%) and testing (1800 eyes; 15%) sets, with 10% of the training data reserved for validation. Following previous studies [10, 12], the LSTM model was trained using the training dataset and evaluated using the test dataset. Its prediction performance was then compared with that of naive BL forecast, OLR, and AR(1) in the test dataset. All analyses were performed in RStudio using the “forecast” package for AR(1) and the “tensorflow” and “keras” packages for the LSTM model.
2.3. Naïve BL Forecast
For each eye and each VF location, the predicted TDV at the next examination was assumed to be equal to the TDV observed at the immediately preceding examination:
| (3) |
This parameter‐free forecasting approach serves as a reference benchmark, reflecting the assumption that no measurable change occurs between consecutive examinations.
2.4. OLR
For each VF location and each eye, six independent OLR models (Model 1 to Model 6) were fitted separately (3,744,000 models for each scenario). In each model, TDVs were regressed against time according to
| (4) |
where Y t is the TDV at time t, T denotes the time at each examination from the BL, β 0 and β 1 are the model regression coefficient estimated by the least square (LS) method, and e t is the random error term.
2.5. AR(1) Model
Similar to OLR, six separate AR(1) models (Model 1 to Model 6) were fitted independently for each eye and each VF location (3,744,000 models for each scenario). In each model, TDVs from all examinations preceding the target visit were used to predict the TDV at the target visit.
Because TDV series exhibit a systematic temporal trend due to disease progression, each series was first detrended [25]. Accordingly, the observed series for each eye and VF location was modeled as
| (5) |
Here, Y t denotes the TDV at time t (in years from BL), β 0 and β 1 are the intercept and slope estimated by OLS, T denotes the time at each examination from the BL, and e t is the residual error term. The detrended series was obtained as the residuals from the fitted linear model:
| (6) |
The AR(1) model was then fitted to the detrended series:
| (7) |
In this equation, Ø is the autoregressive parameter (−1 < Ø < 1), reflecting the temporal correlation between successive examinations, and Ɛ t is a white noise error term with zero mean and constant variance. ∅ ∅ was estimated by maximum likelihood estimation (MLE) using all available lag‐1 pairs from the detrended series. The one‐step‐ahead prediction on the detrended scale was obtained as
| (8) |
where t t a r g e t is the target examination (sixth to eleventh examination), is the detrended residual at t t a r g e t , and is the detrended residual at the immediately preceding examination. The final TDV prediction was calculated by restoring the estimated linear trend:
| (9) |
As follow‐up examinations were conducted at a 6‐month interval over a 5‐year period, the target examination (sixth to eleventh examination) corresponded to a time point ranging from 2.5 to 5 years from BL, which was substituted into the estimated linear trend in Equation ((9).
2.6. RNN
RNNs based on LSTM cells consist of four components: the forget gate, the input gate, the output gate, and the memory cell, which are defined as follows:
| (10) |
| (11) |
| (12) |
| (13) |
Here, w f , w i , w o , and w C denote the weights, and b f , b i , b o , and b C represent the corresponding biases. The sigmoid and tanh functions are nonlinear activation functions applied within the gates. The previous hidden state, h (t−1), together with the memory cell C (t), allows the network to retain, discard, and update temporal information through elementwise multiplication (⊗). The hidden state is then computed as
| (14) |
The operation of the three gates and the memory cell governs the selective retention, discarding, and updating of temporal information, thereby facilitating the learning of long‐term dependencies [10].
2.7. Proposed Method
In our proposed method, a moving “sliding window” of length five was constructed independently for each eye, such that in each window, five consecutive VF examinations were used as model inputs and the subsequent examination was used as model output. As more examinations became available, the number of training windows increased from one (using the first six examinations) to six (using all eleven examinations). For model evaluation, however, only the final (most recent) window for each test eye was used to predict the subsequent VF examination. The few previous studies [10, 12] that used LSTM considered only the first five examinations as inputs for predicting the sixth examination, without applying a sliding window structure due to the limited number of examinations.
To prevent information leakage between training and testing sets, data splitting was performed at the eye level before window construction. The dataset (12,000 eyes) was randomly divided into a training set (10,200 eyes; 85%) and a testing set (1800 eyes; 15%), after which sliding windows were generated independently for each eye. Regardless of the training strategy, the primary objective of all models was to generate the predicted value for the subsequent VF examination. Thus, the main distinction between the models concerned how the model parameters were estimated. Unlike OLR and AR(1), which were fitted separately for each eye and each VF location, only six LSTM models (with the number of training windows increasing from one in Model 1 to six in Model 6) were trained on pooled data from all eyes and locations. Accordingly, OLR and AR(1) represented personalized prediction models with eye‐specific parameter estimation, whereas LSTM represented a population‐level prediction model with a single set of learned weights applied identically to all test eyes.
Each LSTM model received 104 normalized input features comprising 52 TDVs and 52 PDVs. The network architecture was selected empirically for each prediction model. The architecture consisted of two stacked LSTM layers followed by a fully connected output layer with 52 neurons, each corresponding to one VF test location. The first LSTM layer contained 64 hidden units, and the second layer varied according to the progression scenario and prediction model (slow: 50–140 units; moderate: 30–120 units; and fast: 60–170 units). The final output consisted of the predicted TDVs at all 52 test locations.
Models were trained using the Adam optimizer (learning rate = 0.0001) with mean squared error (MSE) as the loss function, a batch size of 512, and a maximum of 100 epochs. To reduce overfitting, early stopping was applied based on validation loss, with a patience of four epochs and a minimum improvement threshold of 0.001, restoring the weights corresponding to the best validation performance. L2 weight regularization was applied to the output layer, and in some architectures, a dropout layer (rate = 0.1) was inserted between the final LSTM and output layers. Training continued until the early stopping criterion was satisfied. The complete training hyperparameters are summarized in Table 3.
Table 3.
Summary of model training hyperparameters.
| Parameter | Optimizer | Learning rate | Loss function | Batch size | Maximum epochs | Validation split | Early stopping monitor | Dropout |
|---|---|---|---|---|---|---|---|---|
| Value/description | Adam | 0.0001 | MSE | 512 | 100 | 10% of training data | Validation loss | 0.1 |
2.8. Model Evaluation Metrics
To evaluate the prediction performance of the models, two accuracy metrics were used, including the root mean squared error (RMSE) and the pointwise mean absolute error (PMAE) of TDV [10–12, 42]. The RMSE for each eye was computed using the following formula:
| (15) |
where n denotes each test location in the VF examination. The PMAE across all eyes for each VF test location was computed using the following formula:
| (16) |
where i represents the eye, n represents the test location of VF examination, and TDVi,n represents TDV of eye i, test location n. RMSE and PMAE metrics were calculated for naive BL forecast, LSTM, AR(1), and OLR. To compare the accuracy metrics among four methods, repeated measures one‐way analysis of variance was used, followed by paired t‐tests with Bonferroni correction for pairwise comparisons. Cohen′s d or effect size was calculated for pairwise comparisons. Model evaluation was performed using SPSS software Version 23.0.
3. Results
Tables 4, 5, and 6 summarize the overall RMSE and PMAE prediction errors for the slow, moderate, and fast progression rate scenarios, respectively. Multiple comparisons among the four models were significant across all progression rate scenarios. In the slow progression rate scenario, pairwise comparisons showed that the RMSE and PMAE for LSTM at the sixth examination were 3.33 ± 0.37 dB and 2.61 ± 0.69 dB, respectively, which were significantly lower than those of AR(1) and OLR (Cohen′s d > 0.8). As the number of examinations increased, the prediction errors for LSTM gradually increased, whereas those for AR(1) decreased, so that, at the eleventh examination, AR(1) performed slightly better than LSTM in RMSE (3.49 ± 0.44 dB vs. 3.56 ± 0.40 dB; Cohen′s d = 0.20) and moderately better than LSTM in PMAE (2.66 ± 0.82 dB vs. 2.79 ± 0.71 dB; Cohen′s d = 0.63). OLR and naive BL forecast consistently exhibited higher prediction errors than both LSTM and AR(1) in all examinations, with large effect sizes (Cohen′s d > 0.8). Compared with OLR, the naive BL forecast also yielded higher prediction errors from the seventh examination onward, with Cohen′s d ranging from 0.24 to 1.55 for RMSE and from 1.62 to 3.48 for PMAE (Table 4)
Table 4.
Comparison among OLR, LSTM, and AR(1) models in slow progression rate scenario.
| 6th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.43 | 2.24 | 0.90 | 1.20 | 0.24 | 2.16 | ||||
| PMAE (dB) |
|
|
|
|
2.76 | 2.09 | 0.98 | 2.66 | 1.62 | 1.99 | ||||
| 7th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.32 | 1.94 | 0.63 | 1.42 | 0.38 | 2.15 | ||||
| PMAE (dB) |
|
|
|
|
3.26 | 2.04 | 0.69 | 3.21 | 1.87 | 2.14 | ||||
| 8th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.20 | 1.51 | 0.13 | 1.61 | 0.88 | 1.99 | ||||
| PMAE (dB) |
|
|
|
|
2.82 | 1.48 | 0.18 | 2.79 | 2.47 | 1.80 | ||||
| 9th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.21 | 1.21 | 0.04 | 1.75 | 1.11 | 1.84 | ||||
| PMAE (dB) |
|
|
|
|
2.58 | 1.06 | 0.28 | 3.14 | 3.48 | 1.74 | ||||
| 10th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.11 | 0.99 | 0.07 | 1.91 | 1.38 | 1.85 | ||||
| PMAE (dB) |
|
|
|
|
3.00 | 0.95 | 0.39 | 3.06 | 2.99 | 1.77 | ||||
| 11th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.01 | 0.76 | 0.20 | 1.89 | 1.55 | 1.81 | ||||
| PMAE (dB) |
|
|
|
|
3.09 | 0.67 | 0.63 | 3.18 | 3.06 | 1.65 | ||||
Abbreviations: AR(1), first‐order autoregressive; BL, naive baseline forecast; LSTM, long short‐term memory; OLR, ordinary linear regression; PMAE, pointwise mean absolute error; RMSE, root mean squared error; SD, standard deviation.
Table 5.
Comparison among OLR, LSTM, and AR(1) models in moderate progression rate scenario.
| 6th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.41 | 2.00 | 0.61 | 1.15 | 0.21 | 1.88 | ||||
| PMAE (dB) |
|
|
|
|
3.21 | 1.80 | 0.41 | 2.87 | 1.98 | 1.63 | ||||
| 7th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.21 | 1.41 | 0.21 | 1.38 | 0.46 | 1.70 | ||||
| PMAE (dB) |
|
|
|
|
3.35 | 1.31 | 0.07 | 3.17 | 1.90 | 1.55 | ||||
| 8th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.12 | 0.94 | 0.13 | 1.58 | 0.95 | 1.51 | ||||
| PMAE (dB) |
|
|
|
|
2.99 | 0.76 | 0.43 | 2.78 | 2.18 | 1.19 | ||||
| 9th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.13 | 0.65 | 0.40 | 1.71 | 1.15 | 1.42 | ||||
| PMAE (dB) |
|
|
|
|
2.57 | 0.37 | 0.62 | 2.90 | 2.98 | 1.00 | ||||
| 10th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
0.95 | 0.95 | 0.47 | 1.77 | 1.39 | 1.49 | ||||
| PMAE (dB) |
|
|
|
|
3.16 | 0.25 | 0.72 | 3.26 | 2.94 | 1.19 | ||||
| 11th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
0.95 | 0.06 | 0.98 | 1.94 | 1.57 | 1.15 | ||||
| PMAE (dB) |
|
|
|
|
3.05 | 0.13 | 0.63 | 2.99 | 2.75 | 0.57 | ||||
Abbreviations: AR(1), first‐order autoregressive; BL, naive baseline forecast; LSTM, long short‐term memory; OLR, ordinary linear regression; PMAE, pointwise mean absolute error; RMSE, root mean squared error; SD, standard deviation.
Table 6.
Comparison among OLR, LSTM, and AR(1) models in fast progression rate scenario.
| 6th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.29 | 1.65 | 0.46 | 1.08 | 0.23 | 1.57 | ||||
| PMAE (dB) |
|
|
|
|
2.47 | 1.39 | 0.17 | 2.67 | 1.17 | 1.35 | ||||
| 7th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.17 | 1.21 | 0.14 | 1.31 | 0.48 | 1.52 | ||||
| PMAE (dB) |
|
|
|
|
2.90 | 1.07 | 0.19 | 2.80 | 1.76 | 1.31 | ||||
| 8th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
1.06 | 0.67 | 0.31 | 1.54 | 0.92 | 1.25 | ||||
| PMAE (dB) |
|
|
|
|
2.58 | 0.47 | 0.58 | 2.71 | 2.21 | 0.95 | ||||
| 9th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
0.93 | 0.05 | 0.85 | 1.60 | 1.16 | 0.91 | ||||
| PMAE (dB) |
|
|
|
|
3.12 | 0.12 | 0.93 | 3.24 | 2.90 | 0.54 | ||||
| 10th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
0.88 | 0.27 | 1.12 | 1.67 | 1.35 | 0.76 | ||||
| PMAE (dB) |
|
|
|
|
3.13 | 0.25 | 0.86 | 2.92 | 2.59 | 0.45 | ||||
| 11th examination | OLR | AR(1) | LSTM | BL | Cohen′s d | |||||||||
| OLR vs. AR(1) | OLR vs. LSTM | AR(1) vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||||||
| RMSE (dB) |
|
|
|
|
0.89 | 0.63 | 1.47) | 1.80 | 1.50 | 0.53 | ||||
| PMAE (dB) |
|
|
|
|
2.93 | 0.52 | 0.97) | 3.08 | 3.01 | 0.17 | ||||
Abbreviations: AR(1), first‐order autoregressive; BL naive baseline forecast; LSTM, long short‐term memory; OLR, ordinary linear regression; PMAE, pointwise mean absolute error; RMSE, root mean squared error; SD, standard deviation.
Similarly, in the moderate progression rate scenario, at the sixth examination, the lowest RMSE (3.51 ± 0.38 dB) and PMAE (2.75 ± 0.75 dB) were observed for LSTM. However, the effect size between LSTM and AR(1) was medium for RMSE (Cohen′s d = 0.61) and small for PMAE (Cohen′s d = 0.41). As the number of examinations increased, LSTM prediction errors increased. From the eighth examination onward, AR(1) performed better than LSTM in PMAE (Cohen′s d = 0.43–0.72), and from the ninth examination onward, it also performed better than LSTM in RMSE (Cohen′s d ranging from 0.40 to 0.98). At the eleventh examination, AR(1) achieved the lowest RMSE and PMAE, which were 3.63 ± 0.47 dB and 2.76 ± 0.85 dB, respectively. Additionally, the differences between OLR and LSTM decreased as the number of examinations increased, with negligible effect sizes observed for both RMSE and PMAE at the eleventh examination. Across all examinations, AR(1) and LSTM outperformed the naive BL forecast (Cohen′s d > 0.8), and OLR surpassed it from the seventh examination onward (Table 5).
In the fast progression rate scenario at the sixth examination, the RMSE for LSTM and AR(1) were 3.74 ± 0.44 dB and 3.96 ± 0.56 dB, respectively, corresponding to a small effect size (Cohen′s d = 0.46). At the seventh examination, the difference between LSTM and AR(1) was negligible (Cohen′s d = 0.14). In subsequent examinations with Cohen′s d increasing from 0.31 to 1.47, the difference between the two models increased, and AR(1) performed better than LSTM. Additionally, from the eighth examination onward, AR(1) exhibited a lower PMAE than LSTM, with Cohen′s d increasing from 0.58 at the eighth examination to 0.97 at the eleventh examination. As the number of examinations increased, the prediction errors for OLR decreased compared with those for LSTM, so that, at the eleventh examination, OLR achieved RMSE and PMAE of 4.30 ± 0.62 dB and 3.25 ± 1.10 dB, respectively. The differences between OLR and LSTM corresponded to medium effect sizes (Cohen′s d = 0.63 for RMSE and 0.52 for PMAE). Across all examinations, the RMSE and PMAE for AR(1) were lower than those of OLR and naive BL forecast, with consistently large effect sizes (Cohen′s d > 0.8). OLR also outperformed the naive BL forecast from the seventh examination onward, with an increasingly larger difference over subsequent examinations .In contrast, LSTM outperformed the naive BL forecast at earlier examinations, but its advantage gradually diminished as the number of examinations increased (Cohen′s d for RMSE decreased from 1.57 at the sixth examination to 0.53 at the eleventh examination) (Table 6).
Table 7 describes the overall RMSE across the VF sectors for the sixth examination in the slow progression rate scenario. As in the study by Kim et al. [10], the prediction error was compared among the six sectors of the 24‐2 VF proposed by Garway‐Heath et al. [43], as well as central and peripheral zones. The LSTM prediction errors were lower than those of AR(1) and OLR across all sections. Effect sizes for LSTM and AR(1) comparisons were small to medium in most sections (Cohen′s d = 0.26–0.52). However, a larger effect in the peripheral zone (Cohen′s d = 0.86) and a negligible effect in the nasal sector (Cohen′s d = 0.07) were observed. In contrast, effect sizes for OLR and LSTM comparisons (Cohen′s d = 0.63–2.06) and OLR and AR(1) comparisons (Cohen′s d = 0.50–1.29) were medium to large across all sections. The naive BL forecast showed higher RMSE than both LSTM and AR(1) across all VF sectors and zones, whereas the differences between OLR and naive BL forecast were small.
Table 7.
Comparison among OLR, LSTM, and AR(1) models in the six Garway‐Heath sectors and peripheral and central zones of the 24–2 visual field test for the sixth examination in the slow progression rate scenario.
| Prediction error (RMSE, dB), mean ± SD | (Cohen′s d) | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| OLR | LSTM | AR(1) | BL | AR(1) vs. LSTM | AR(1) vs. OLR | OLR vs. LSTM | BL vs. AR(1) | BL vs. OLR | BL vs. LSTM | |||||
| Superotemporal |
|
|
|
|
0.33 | 0.69 | 0.98 | 0.57 | 0.12 | 0.92 | ||||
| Superonasal |
|
|
|
|
0.41 | 0.92 | 1.25 | 0.78 | 0.14 | 1.19 | ||||
| Nasal |
|
|
|
|
0.07 | 0.60 | 0.63 | 0.50 | 0.10 | 0.58 | ||||
| Inferonasal |
|
|
|
|
0.35 | 0.60 | 0.91 | 0.52 | 0.06 | 0.89 | ||||
| Inferotemporal |
|
|
|
|
0.52 | 0.67 | 1.16 | 0.55 | 0.12 | 1.08 | ||||
| Temporal |
|
|
|
|
0.28 | 0.50 | 0.74 | 0.40 | 0.10 | 1.44 | ||||
| Central |
|
|
|
|
0.26 | 0.86 | 1.03 | 0.75 | 0.13 | 0.99 | ||||
| Peripheral |
|
|
|
|
0.86 | 1.29 | 2.06 | 1.09 | 0.22 | 2.00 | ||||
Abbreviations: AR(1), first‐order autoregressive; BL naive baseline forecast; LSTM long short‐term memory; OLR ordinary linear regression; RMSE, root mean squared error; SD, standard deviation.
Figure 4 illustrates the overall RMSE ± SD and PMAE ± SD of OLR, LSTM, AR(1), and naive BL forecast across all scenarios for examinations 6–11, whereas Figure 5 shows the number of eyes categorized based on the RMSE prediction error for the sixth examination in the slow progression rate scenario. For LSTM, most eyes exhibited prediction errors of 2–3 dB or 3–4 dB. Specifically, 347 eyes had prediction errors of 2–3 dB, and 1364 eyes had prediction errors of 3–4 dB. For OLR, 299 eyes had prediction errors of 3–4 dB, 1187 eyes had prediction errors of 4–5 dB, and 304 eyes had prediction errors of 5–6 dB. For AR(1), 102 eyes had prediction errors of 2–3 dB, 1242 eyes had prediction errors of 3–4 dB, and 441 eyes had prediction errors of 4–5 dB. For naive BL forecast, 390 eyes had prediction errors of 3–4 dB, 1171 eyes had prediction errors of 4–5 dB, and 227 eyes had prediction errors of 5–6 dB.
Figure 4.

The overall RMSE and PMAE values obtained by predicting the TD values for the 6th–11th examinations. The plot shows the RMSE and PMAE values obtained by predicting the TD values. Error bars indicate standard deviations, as shown in Tables 4, 5 and 6. RMSE, root mean squared error; PMAE, pointwise mean absolute error; TD, total deviation.
Figure 5.

Number of eyes categorized according to the prediction error of RMSE for the sixth examination in the slow progression rate scenario. RMSE, root mean squared error.
Figure 6 indicates the PMAE across 52 VF locations for the sixth examination in the slow progression rate scenario. LSTM had significantly lower prediction errors compared with AR(1) (red dots), OLR (blue dots), and naive BL forecast (green dots) in most locations.
Figure 6.

Pointwise mean absolute error (PMAE) of predicted total deviation value (TDV) for the sixth examination in the slow progression rate scenario. The darker colors indicate higher error. The darker color indicates higher error. Red dots indicate significant differences between AR(1) and LSTM; blue dots indicate significant differences between LSTM and OLR; green dots indicate significant differences between BL and LSTM (paired t‐test). OLR, ordinary linear regression; LSTM, long short‐term memory; AR(1), first‐order auto regressive; BL, naive baseline forecast.
4. Discussion
VF prediction is important in patients with glaucoma because it allows clinicians to assess future VF function and helps them decide on optimal treatment and follow‐up schedules. It is challenging to choose a model that is both interpretable and capable of predicting complex patterns incorporating noise in VF data.
Although deep learning methods such as RNNs have shown high accuracy for predicting complex patterns in VF data, their ambiguity in interpretation (often referred to as a black box) and their need for large training datasets make them difficult to implement in clinical settings. Conversely, conventional statistical methods, such as OLR, are more interpretable but are unable to account for temporal dependence among repeated VF measurements. AR(1) models address this limitation by explicitly modeling serial correlation while retaining the interpretability and simplicity of a statistical approach. Accordingly, the present study compared personalized statistical models (OLR and AR(1)) with a population‐level deep learning model (LSTM) to predict longitudinal VF measurements in glaucoma. For this purpose, different progression rate and pattern scenarios were designed, and the LMM approach was used to realistically simulate pointwise VF measurements over time. The performance of the models was then evaluated under slow, moderate, and fast progression scenarios with varying numbers of examinations.
Across all progression scenarios, AR(1) consistently achieved significantly lower prediction errors than OLR. This result is expected because, unlike OLR, AR(1) explicitly accounts for temporal dependence between consecutive VF measurements. Although most VF locations were stable or exhibited slow progression, AR(1) outperformed the naive BL forecast, particularly as the number of examinations increased. These results indicate that the superior performance of AR(1) cannot be explained solely by the inherent stability of the VF series; rather, modeling temporal dependence captures meaningful predictive information beyond persistence alone.
Several studies [13, 14, 19, 20] have reported that OLR demonstrates satisfactory performance for VF prediction in glaucoma. More recently, however, increasingly sophisticated statistical and machine learning methods have been proposed to improve predictive accuracy [10–12, 21, 23–25]. For example, Childress et al. [24] evaluated the OLR, AR(1), and first‐order irregular autoregressive [IAR(1)] models and found that both autoregressive models showed better predictive performance than OLR. In another related study, Naderi et al. [25] examined the predictive performance of AR(1), OLR, exponential, and polynomial models in patients with glaucoma. AR(1) showed superior pointwise predictive accuracy, particularly in advanced glaucoma, where it had significantly lower RMSE and PMAE than the other models. Consistent with these findings, AR(1) also significantly outperformed OLR across all VF progression rate scenarios in the present study.
Murata et al. [11] compared variational Bayesian linear regression (VBLR) with OLR for predicting TDVs by varying the number of input VF data from 2 to 10. VBLR exhibited more accurate predictions than OLR (p < 0.001), with an overall RMSE of 4.5 ± 2.4 dB when predicting the sixth VF examination. Although the LSTM model in the present study yielded lower RMSE values across all simulated progression scenarios, direct comparison between the two studies should be interpreted cautiously because of substantial differences in the study design, simulation framework, and evaluation methodology. Consequently, our findings do not establish the superiority of LSTM over the VBLR approach.
Similarly, Park et al. [12] demonstrated that LSTM significantly outperformed OLR when predicting the sixth VF examination and attributed this improvement to the incorporation of reliability indices (false positive rates, false negative rates, and fixation losses) as additional model inputs. Kim et al. [10] likewise reported that both bidirectional gated recurrent unit (bi‐GRU) and LSTM had superior predictive performance compared with OLR when predicting the sixth VF examination. Our findings are consistent with those reported by Park et al. [12] and Kim et al. [10]. The present study showed that the population‐level LSTM model achieved significantly better predictive performance than both personalized statistical models (OLR and AR(1)) at the sixth examination across all progression rate scenarios. This is likely because only five previous examinations were available for estimating the parameters of the personalized models. Under these conditions, OLR and AR(1) must estimate regression and autoregressive parameters, respectively, from a small number of longitudinal observations, which may lead to greater uncertainty in parameter estimation. In contrast, the LSTM model benefits from population‐level learning by leveraging information from all training eyes to estimate a common set of model parameters, which may allow it to generate more robust predictions when only a limited number of previous examinations are available.
Contrary to expectations, increasing the number of examinations did not improve the performance of LSTM compared with AR(1) in the present study. Rather, AR(1) obtained superior predictive accuracy, particularly under moderate and fast progression scenarios. These results should be interpreted in the context of the distinction between personalized and population‐level modeling frameworks. Because AR(1) was estimated separately for each eye, additional examinations progressively provided more information for estimating eye‐specific temporal dynamics and autoregressive parameters. This enhanced the model′s ability to characterize each eye′s longitudinal progression trajectory. In contrast, the LSTM model relied on a single set of population‐level weights that remained unchanged during prediction. Consequently, the improved performance of AR(1) with longer follow‐up likely reflects the advantages of personalized model fitting rather than an inherent superiority of autoregressive modeling over deep learning.
Although the predictive performance of the AR(1) model improved as the number of examinations increased, the number of examinations was still small, with the longest follow‐up in the present study comprising only 11 examinations over a 5‐year period. Hence, the superior performance of the AR(1) model should be interpreted in the context of short‐term prediction. This observation is in agreement with previous studies indicating that AR(1) models are particularly suitable for short‐term forecasting [44, 45].
More recently, Abbasi et al. [21] compared a Hybrid‐VF‐Net architecture integrating convolutional neural networks (CNNs) and RNNs with standalone CNN and RNN models for VF prediction without incorporating reliability indices. The hybrid architecture had the lowest prediction error and outperformed both individual models. Furthermore, prediction accuracy declined for less reliable VF examinations, emphasizing the importance of reliability indices in VF prediction. These findings suggest that including complementary clinical information beyond VF measurements alone may further enhance predictive performance.
Finally, the present study has several limitations. First, although the proposed simulation framework accounted for between‐eye heterogeneity, within‐eye dependence, clinically defined progression scenarios, and realistic measurement variability, its reliance on a linear mixed‐effects model resulted in predominantly linear progression trajectories. Consequently, the simulated datasets may not fully capture the nonlinear progression patterns observed in clinical glaucoma populations. Nevertheless, by incorporating random effects, the framework accounted for heterogeneity in disease progression trajectories and BL VFs across eyes, generating simulated data that more closely reflected the complexity and variability of real‐world clinical data than those generated using simple linear regression [18]. Second, because the simulated datasets were generated from an inherently linear process, the better performance of the AR(1) model compared with the LSTM may have been partly attributable to a bias in favor of linear models. Third, although simulated data have the advantage of controlling progression rates, measurement variability, and follow‐up schedules while facilitating standardized model comparisons, they rely on predefined parametric assumptions and thus cannot fully reproduce the complex variability, noise characteristics, and irregular missing patterns observed in real‐world datasets. Fourth, the LSTM input representation in this study treats VF test locations as flat feature vectors, effectively ignoring the anatomically structured spatial relationships between them. While spatiotemporal architectures [21] could potentially capture these dependencies to improve accuracy, our primary objective was focusing on the temporal dynamics of progression using a standard LSTM model. Fifth, reliability indices—including false‐positive rates, false‐negative rates, and fixation losses—were not incorporated because these variables could not be realistically simulated within the current framework. Sixth, the autoregressive coefficient was estimated from relatively short longitudinal series and may, therefore, be subject to Hurwicz small‐sample bias, which generally leads to underestimation of the autoregressive coefficient. However, a previous simulation study [46] suggests that this estimation bias has a relatively little effect on point forecasts although it may have greater implications for prediction‐interval coverage and multistep forecasting. Given that the AR(1) models in the present study were primarily used for one‐step‐ahead point prediction, the practical impact of this bias on predictive accuracy is likely to be limited.
Accordingly, future research should evaluate the performance of these prediction models using real world clinical VF data as a complementary approach to the present simulation study. Using anatomically informed spatiotemporal architectures may further improve predictive accuracy by modeling spatial and temporal dependencies. In addition, future studies should integrate reliability indices together with VF measurements into deep learning models to enable a more clinically representative comparison with personalized statistical approaches. Because the naive BL forecast used in the present study relies solely on the most recent measurement, it may be affected by variability in a single VF measurement. A simple running‐average forecast could therefore be considered as an additional BL in future studies, as it may reduce the influence of variability in a single recent measurement and provide a reference for assessing the robustness of prediction model comparisons. Finally, simulating data with irregular follow‐up intervals and applying the IAR(1) model, which is more appropriate for irregularly spaced longitudinal data, could be a valuable direction for future research.
5. Conclusion
This study developed a simulation framework for comparing pointwise VF prediction models over time and showed that predictive performance depended on both the modeling approach and the amount of longitudinal data available. The population‐level LSTM achieved the highest accuracy with limited number of examinations, whereas the personalized AR(1) model became more accurate as additional examinations were conducted. Given that the length of the VF series was still relatively short as the number of examinations increased, these findings support the suitability of AR(1) for short‐term prediction. Overall, when a sufficient number of examinations is available, the results of the present study emphasize the value of personalized temporal modeling and highlight the potential importance of accounting for individual progression trajectories in longitudinal VF prediction.
Author Contributions
S.N.: data generation, data analysis, and writing of the draft. E.M.: programming advisor. M.M.: clinical advisor and draft approver. M.D.: statistical advisor and drafted the manuscript. S.P.: research idea, study design, and improvement of the manuscript.
Funding
This study was supported by the Research Council of Shiraz University of Medical Sciences, 29800.
Disclosure
This paper was extracted from Soheila Naderi′s Ph.D. thesis under the supervision of Prof. Saeedeh Pourahmad.
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting Information
Additional supporting information can be found online in the Supporting Information section.
Supporting information
Supporting Information 1 The codes for generating the simulated data are provided in File S1.
Supporting Information 2 The codes used for data analysis are provided in File S2.
Supporting Information 3.
Acknowledgments
This work was supported by the Research Council of Shiraz University of Medical Sciences (SUMS), Shiraz, Iran (grant no. 29800), and was approved by the Ethics Committee of SUMS (IR.SUMS.MED.REC.1403.393). The authors acknowledge the use of ChatGPT (GPT‐3.5) for grammar and language editing assistance during the preparation of this manuscript. The content and scientific integrity of the manuscript remain the sole responsibility of the authors. The authors would like to thank the Center for Development of Clinical Research of Nemazee Hospital and Dr. Nasrin Shokrpour for editorial assistance.
Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
References
- 1. Gupta P. and Yadav K. S., Applications of Microneedles in Delivering Drugs for Various Ocular Diseases, Life Sciences. (2019) 237, 116907, 10.1016/j.lfs.2019.116907, 31606378. [DOI] [PubMed] [Google Scholar]
- 2. Resnikoff S., Pascolini D., Etya′Ale D., Kocur I., Pararajasegaram R., Pokharel G. P., and Mariotti S. P., Global Data on Visual Impairment in the Year 2002, Bulletin Of The World Health Organization. (2004) 82, no. 11, 844–851, 10.1076/opep.11.2.67.28158, 15640920. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 3. Tham Y.-C., Li X., Wong T. Y., Quigley H. A., Aung T., and Cheng C.-Y., Global Prevalence of Glaucoma and Projections of Glaucoma Burden Through 2040: A Systematic Review and Meta-Analysis, Ophthalmology. (2014) 121, no. 11, 2081–2090, 10.1016/j.ophtha.2014.05.013, 24974815. [DOI] [PubMed] [Google Scholar]
- 4. Buckingham B. P., Inman D. M., Lambert W., Oglesby E., Calkins D. J., Steele M. R., Vetter M. L., Marsh-Armstrong N., and Horner P. J., Progressive Ganglion Cell Degeneration Precedes Neuronal Loss in a Mouse Model of Glaucoma, Journal of Neuroscience. (2008) 28, no. 11, 2735–2744, 10.1523/JNEUROSCI.4443-07.2008, 18337403. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 5. Alencar L. M. and Medeiros F. A., The Role of Standard Automated Perimetry and Newer Functional Methods for Glaucoma Diagnosis and Follow-Up, Indian journal of ophthalmology. (2011) 59, no. 7, S53–S58, 10.4103/0301-4738.73694. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 6. Swaminathan S. S., Berchuck S. I., Jammal A. A., Rao J. S., and Medeiros F. A., Rates of Glaucoma Progression Derived From Linear Mixed Models Using Varied Random Effect Distributions, Translational Vision Science & Technology. (2022) 11, no. 2, 16–16, 10.1167/tvst.11.2.16, 35138343. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 7. Asano S., Murata H., Matsuura M., Fujino Y., and Asaoka R., Early Detection of Glaucomatous Visual Field Progression Using Pointwise Linear Regression With Binomial Test in the Central 10 Degrees, American Journal of Ophthalmology. (2019) 199, 140–149, 10.1016/j.ajo.2018.11.004, 30465746. [DOI] [PubMed] [Google Scholar]
- 8. Chauhan B. C., Drance S. M., and Douglas G. R., The Use of Visual Field Indices in Detecting Changes in the Visual Field in Glaucoma, Investigative Ophthalmology & Visual Science. (1990) 31, no. 3, 512–520, 2318590. [PubMed] [Google Scholar]
- 9. Spry P. G., Bates A. B., Johnson C. A., and Chauhan B. C., Simulation of Longitudinal Threshold Visual Field Data, Investigative Ophthalmology & Visual Science. (2000) 41, no. 8, 2192–2200, 10892862. [PubMed] [Google Scholar]
- 10. Kim H., Lee J., Moon S., Kim S., Kim T., Jin S. W., Kim J. L., Shin J., Lee S. U., Jang G., Hu Y., and Park J. R., Visual Field Prediction Using a Deep Bidirectional Gated Recurrent Unit Network Model, Scientific Reports. (2023) 13, no. 1, 11154, 10.1038/s41598-023-37360-1, 37429862. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 11. Murata H., Araie M., and Asaoka R., A New Approach to Measure Visual Field Progression in Glaucoma Patients Using Variational Bayes Linear Regression, Investigative Ophthalmology & Visual Science. (2014) 55, no. 12, 8386–8392, 10.1167/iovs.14-14625, 25414192. [DOI] [PubMed] [Google Scholar]
- 12. Park K., Kim J., and Lee J., Visual Field Prediction Using Recurrent Neural Network, Scientific Reports. (2019) 9, no. 1, 10.1038/s41598-019-44852-6, 31182763. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 13. Bengtsson B., Patella V. M., and Heijl A., Prediction of Glaucomatous Visual Field Loss by Extrapolation of Linear Trends, Archives of Ophthalmology. (2009) 127, no. 12, 1610–1615, 10.1001/archophthalmol.2009.297, 20008716. [DOI] [PubMed] [Google Scholar]
- 14. McNaught A., Hitchings R. A., Crabb D., and Fitzke F., Modelling Series of Visual Fields to Detect Progression in Normal-Tension Glaucoma, Graefe′s Archive For Clinical And Experimental Ophthalmology. (1995) 233, no. 12, 750–755, 10.1007/BF00184085, 8626082. [DOI] [PubMed] [Google Scholar]
- 15. Moradi M., Hashemabad S. K., Vu D. M., Soneru A. R., Fujita A., Wang M., Elze T., Eslami M., and Zebardast N., Pyglaucometrics: A Stacked Weight-Based Machine Learning Approach For Glaucoma Detection Using Visual Field Data, Medicina. (2025) 61, no. 3, 10.3390/medicina61030541, 40142352. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 16. Nouri-Mahdavi K., Hoffman D., Ralli M., and Caprioli J., Comparison of Methods to Predict Visual Field Progression in Glaucoma, Archives of Ophthalmology. (2007) 125, no. 9, 1176–1181, 10.1001/archopht.125.9.1176. [DOI] [PubMed] [Google Scholar]
- 17. O′Leary N., Chauhan B. C., and Artes P. H., Visual Field Progression in Glaucoma: Estimating the Overall Significance of Deterioration With Permutation Analyses of Pointwise Linear Regression (PoPLR), Investigative Ophthalmology & Visual Science. (2012) 53, no. 11, 6776–6784, 10.1167/iovs.12-10049, 22952123. [DOI] [PubMed] [Google Scholar]
- 18. Rabiolo A., Morales E., Mohamed L., Capistrano V., Kim J. H., Afifi A., Yu F., Coleman A. L., Nouri-Mahdavi K., and Caprioli J., Comparison of Methods to Detect and Measure Glaucomatous Visual Field Progression, Translational Vision Science & Technology. (2019) 8, no. 5, 2–2, 10.1167/tvst.8.5.2, 31555493. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 19. Kwon Y. H., Kim C.-S., Zimmerman M. B., Alward W. L., and Hayreh S. S., Rate of Visual Field Loss and Long-Term Visual Outcome in Primary Open-Angle Glaucoma, American Journal Of Ophthalmology. (2001) 132, no. 1, 47–56, 10.1016/S0002-9394(01)00912-6. [DOI] [PubMed] [Google Scholar]
- 20. Pereira M. L. M., Kim C.-S., Zimmerman M. B., Alward W. L., Hayreh S. S., and Kwon Y. H., Rate and Pattern of Visual Field Decline in Primary Open-Angle Glaucoma, Ophthalmology. (2002) 109, no. 12, 2232–2240, 10.1016/S0161-6420(02)01248-4, 12466164. [DOI] [PubMed] [Google Scholar]
- 21. Abbasi A., Gowrisankaran S., Lin W.-C., Song X., Antony B. J., Wollstein G., Schuman J. S., and Ishikawa H., A Hybrid Deep Learning–Based Approach for Visual Field Test Forecasting, Ophthalmology Science. (2025) 5, no. 5, 100803, 10.1016/j.xops.2025.100803, 40520474. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 22. Pathak M., Demirel S., and Gardiner S. K., Nonlinear Trend Analysis of Longitudinal Pointwise Visual Field Sensitivity in Suspected and Early Glaucoma, Translational Vision Science & Technology. (2015) 4, no. 1, 8–8, 10.1167/tvst.4.1.8, 25694844. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 23. Wen J. C., Lee C. S., Keane P. A., Xiao S., Rokem A. S., Chen P. P., Wu Y., and Lee A. Y., Forecasting Future Humphrey Visual Fields Using Deep Learning, PloS One. (2019) 14, no. 4, e0214875, 10.1371/journal.pone.0214875, 30951547. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 24. Childress C., Time Series Analysis of Longitudinally Collected Standard Autoperimetry Data in Glaucoma Patients, 2023, College thesis, Muray State University. [Google Scholar]
- 25. Naderi S., Arman A., Shahparast S., Vermeer K. A., Lemij H. G., and Pourahmad S., Future Visual Field Prediction in Glaucoma: An Application of First-Order Autoregressive Approach, International Ophthalmology. (2026) 46, no. 1, 10.1007/s10792-026-04211-7, 42635842. [DOI] [PubMed] [Google Scholar]
- 26. Bryan S. R., Vermeer K. A., Eilers P. H., Lemij H. G., and Lesaffre E. M., Robust and Censored Modeling and Prediction of Progression in Glaucomatous Visual Fields, Investigative Ophthalmology & Visual Science. (2013) 54, no. 10, 6694–6700, 10.1167/iovs.12-11185, 24030462. [DOI] [PubMed] [Google Scholar]
- 27. Chauhan B. C., Garway-Heath D. F., Goñi F. J., Rossetti L., Bengtsson B., Viswanathan A. C., and Heijl A., Practical Recommendations for Measuring Rates of Visual Field Change in Glaucoma, British Journal of Ophthalmology. (2008) 92, no. 4, 569–573, 10.1136/bjo.2007.135012, 18211935. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 28. Gardiner S. K., Mansberger S. L., and Demirel S., Detection of Functional Change Using Cluster Trend Analysis in Glaucoma, Investigative Ophthalmology & Visual Science. (2017) 58, no. 6, BIO180–BIO190, 10.1167/iovs.17-21562, 28715580. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 29. Medeiros F. A., Weinreb R. N., Moore G., Liebmann J. M., Girkin C. A., and Zangwill L. M., Integrating Event- and Trend-Based Analyses to Improve Detection of Glaucomatous Visual Field Progression, Ophthalmology. (2012) 119, no. 3, 458–467, 10.1016/j.ophtha.2011.10.003, 22266109. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 30. Vesti E., Johnson C. A., and Chauhan B. C., Comparison of Different Methods for Detecting Glaucomatous Visual Field Progression, Investigative Ophthalmology & Visual Science. (2003) 44, no. 9, 3873–3879, 10.1167/iovs.02-1171, 12939303. [DOI] [PubMed] [Google Scholar]
- 31. Vianna J. R. and Chauhan B. C., How To Detect Progression In Glaucoma, Progress In Brain Research. (2015) 221, 135–158, 10.1016/bs.pbr.2015.04.011. [DOI] [PubMed] [Google Scholar]
- 32. Zhu H., Russell R. A., Saunders L. J., Ceccon S., Garway-Heath D. F., and Crabb D. P., Detecting Changes in Retinal Function: Analysis With Non-Stationary Weibull Error Regression and Spatial Enhancement (ANSWERS), PloS One. (2014) 9, no. 1, e85654, 10.1371/journal.pone.0085654, 24465636. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 33. Li Y., Eizenman M., Shi R. B., Buys Y. M., Trope G. E., and Wong W., A Data-Driven Model for Simulating Longitudinal Visual Field Tests in Glaucoma, Translational Vision Science & Technology. (2023) 12, no. 6, 27–27, 10.1167/tvst.12.6.27, 37382576. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 34. Russell R. A., Garway-Heath D. F., and Crabb D. P., New Insights Into Measurement Variability in Glaucomatous Visual Fields From Computer Modelling, PloS One. (2013) 8, no. 12, e83595, 10.1371/journal.pone.0083595, 24386230. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 35. Sabouri S., Haem E., Masoumpour M., Vermeer K. A., Lemij H. G., Yousefi S., and Pourahmad S., Frequency of Visual Fields Needed to Detect Glaucoma Progression: A Computer Simulation Using Linear Mixed Effects Model, Journal of Glaucoma. (2023) 32, no. 5, 355–360, 10.1097/IJG.0000000000002155, 37054400. [DOI] [PubMed] [Google Scholar]
- 36. Hamilton J. D., Time Series Analysis, 2020, Princeton University Press, 10.2307/j.ctv14jx6sm. [DOI] [Google Scholar]
- 37. Everitt B. and Hothorn T., An Introduction to Applied Multivariate Analysis With R, 2011, Springer Science & Business Media, 10.1007/978-1-4419-9650-3. [DOI] [Google Scholar]
- 38. Fitzmaurice G. M., Laird N. M., and Ware J. H., Applied Longitudinal Analysis, 2011, John Wiley & Sons, 10.1002/9781119513469. [DOI] [Google Scholar]
- 39. Schielzeth H., Dingemanse N. J., Nakagawa S., Westneat D. F., Allegue H., Teplitsky C., Réale D., Dochtermann N. A., Garamszegi L. Z., and Araya-Ajoy Y. G., Robustness of Linear Mixed-Effects Models to Violations of Distributional Assumptions, Methods In Ecology And Evolution. (2020) 11, no. 9, 1141–1152, 10.1111/2041-210X.13434. [DOI] [Google Scholar]
- 40. Jammal A. A., Thompson A. C., Mariottoni E. B., Urata C. N., Estrela T., Berchuck S. I., Tseng H. C., Asrani S., and Medeiros F. A., Rates of Glaucomatous Structural and Functional Change From a Large Clinical Population: The Duke Glaucoma Registry Study, American Journal of Ophthalmology. (2021) 222, 238–247, 10.1016/j.ajo.2020.05.019, 32450065. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 41. Gardiner S. K. and Crabb D. P., Examination of Different Pointwise Linear Regression Methods for Determining Visual Field Progression, Investigative Ophthalmology & Visual Science. (2002) 43, no. 5, 1400–1407, 11980853. [PubMed] [Google Scholar]
- 42. Murata H., Zangwill L. M., Fujino Y., Matsuura M., Miki A., Hirasawa K., Tanito M., Mizoue S., Mori K., Suzuki K., Yamashita T., Kashiwagi K., Shoji N., and Asaoka R., Validating Variational Bayes Linear Regression Method With Multi-Central Datasets, Investigative Ophthalmology & Visual Science. (2018) 59, no. 5, 1897–1904, 10.1167/iovs.17-22907, 29677350. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 43. Garway-Heath D. F., Poinoosawmy D., Fitzke F. W., and Hitchings R. A., Mapping The Visual Field To The Optic Disc In Normal Tension Glaucoma Eyes11the Authors Have No Proprietary Interest In The Development Or Marketing Of Any Product Or Instrument Mentioned In This Article, Ophthalmology. (2000) 107, no. 10, 1809–1815, 10.1016/S0161-6420(00)00284-0. [DOI] [PubMed] [Google Scholar]
- 44. Iaousse M., Jouilil Y., Bouincha M., and Mentagui D., Comparative Simulation Study of Classical and Machine Learning Techniques for Forecasting Time Series Data, International Journal of Online and Biomedical Engineering. (2023) 19, no. 8, 56–65, 10.3991/ijoe.v19i08.39853. [DOI] [Google Scholar]
- 45. Siami-Namini S. and Namin A. S., Forecasting Economics and Financial Time Series: ARIMA vs. LSTM, arXiv preprint 16 Mar. 2018 10.48550/arXiv.1803.06386. [DOI]
- 46. Ledolter J., Estimation Bias in the First-Order Autoregressive Model and Its Impact on Predictions and Prediction Intervals, Communications in Statistics - Simulation and Computation. (2009) 38, no. 4, 771–787, 10.1080/03610910802645354. [DOI] [Google Scholar]
Associated Data
This section collects any data citations, data availability statements, or supplementary materials included in this article.
Supplementary Materials
Supporting Information 1 The codes for generating the simulated data are provided in File S1.
Supporting Information 2 The codes used for data analysis are provided in File S2.
Supporting Information 3.
Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
