**Purpose.**:
We generated a variational Bayes model to predict visual field (VF) progression in glaucoma patients.

**Methods.**:
This retrospective study included VF series from 911 eyes of 547 glaucoma patients as test data, and VF series from 5049 eyes of 2858 glaucoma patients as training data. Using training data, variational Bayes linear regression (VBLR) was created to predict VF progression. The performance of VBLR was compared against ordinary least-squares linear regression (OLSLR) by predicting VFs in the test dataset. The total deviation (TD) values of test patients' 11th VFs were predicted using TD values from their second to 10th VFs (VF2-10), the root mean squared error (RMSE) associated with each approach then was calculated. Similarly, mean TD (mTD) of test patients' 11th VFs was predicted using VBLR and OLSLR, and the absolute prediction errors compared.

**Results.**:
The RMSE resulting from VBLR averaged 3.9 ± 2.1 (SD) and 4.9 ± 2.6 dB for prediction based on the second to 10th VFs (VF2-10) and the second to fourth VFs (VF2-4), respectively. The RMSE resulting from OLSLR was 4.1 ± 2.0 (VF2-10) and 19.9 ± 12.0 (VF2-4) dB. The absolute prediction error (SD) for mTD using VBLR was 1.2 ± 1.3 (VF2-10) and 1.9 ± 2.0 (VF2-4) dB, while the prediction error resulting from OLSLR was 1.2 ± 1.3 (VF2-10) and 6.2 ± 6.6 (VF2-4) dB.

**Conclusions.**:
The VBLR more accurately predicts future VF progression in glaucoma patients compared to conventional OLSLR, especially in short VF series.

^{1}Visual loss impacts on the quality of life of glaucoma patients,

^{2}and, as glaucomatous VF damage is irreversible, it is very important to accurately measure and predict progression of VF defects so timely treatment decisions can be made.

^{3}In the Humphrey Guided Progression Analysis (GPA) software, VF progression is measured using a trend analysis, in which simple ordinary least-squares linear regression (OLSLR) is done on the mean deviation (MD) or the visual field index (VFI) over time.

^{4}The software PROGRESSOR (Medisoft Ltd., Leeds, UK) also uses OLSLR to analyze VF progression in a pointwise manner.

^{5}One concern with applying OLSLR is that the method relies on the hypothesis that VF progression occurs linearly. Indeed, previous reports suggest that nonlinear methods, such as the exponential model or an episodic model, can better estimate VF progression,

^{6–8}although other reports give inconsistent conclusions.

^{9,10}Thus, the most appropriate model for prediction of VF progression still is controversial and the choice may vary throughout the disease process.

^{11}

^{12}and spatial information,

^{13}and it is rational that VF progression models can be improved significantly by using information about the distinctive patterns associated with different types of VF defects. Glaucoma specialists know that glaucomatous VF progression often follows characteristic patterns, such as the nasal step defect and the Bjerrum scotoma.

^{14}Moreover, patients with different types of VF damage often have distinct risk factors

^{15}that can be observed at presentation.

^{16}However, to date, almost all methods that measure pointwise VF progression have failed to take into account the correlation among neighboring test points.

^{17–19}

^{20}Many real-world complex data can be represented, using ML methods, in a low-dimensional subspace corresponding to just tens of categories to which the data belong.

^{21}Accordingly, VF data from glaucoma patients could be clustered into groups based on different patterns of VF damage using ML methods; moreover, a VF progression prediction model could be generated that uses this clustering information. Bayesian methods are used commonly to incorporate prior information into a statistical model to improve its accuracy.

^{12,22,23}Recently, Zhu et al.

^{24}reported a novel approach to measure VF progression that modeled spatial correlations as well as nonstationary variability. In this study, we constructed a different type of Bayesian model to predict VF damage that uses the spatial and temporal patterns associated with glaucomatous VF progression as described in a training data set. We then investigated the usefulness of this novel “variational Bayes linear regression” (VBLR) approach over conventional OLSLR to predict VF progression in an unrelated test data set.

**Table 1**

**Table 1**

Demographics | Value |

Follow-up time, y, mean ± SD | 4.4 ± 1.8 |

Age, y, mean ± SD | 58.7 ± 14.2 |

MD of initial VF, dB, mean ± SD | −7.0 ± 6.7 |

**Table 2**

**Table 2**

Demographics | Value |

Follow-up time, y, mean ± SD | 5.5 ± 1.2 |

Age, y, mean ± SD | 55.4 ± 13.6 |

Sex, male:female | 253:294 |

MD of initial VF, dB, mean ± SD | −7.8 ± 6.2 |

Type of glaucoma, POAG, NTG, ACG, PE, steroid, others | 323, 458, 28, 16, 36, 50 |

*n*

^{th}VF in their series;

*D*is the dimension of the vector

_{t}**t**

*and is equal to 52 in this study. Let*

_{n}**n**

*be the set of indices of data obtained from the*

_{m}*m*

^{th}eye,

**T**

*denotes the set*

_{m}**w**

*is the parameter vector of the*

_{m}*m*

^{th}eye (where the first half and latter half of this vector include the intercept and slope coefficients of all 52 test VF points, respectively). Next, let

*x*denote the interval from the first VF test of the

_{n}*n*

^{th}data,

**Φ**(

*x*) denotes a matrix defined as

_{n}**Φ**(

*x*) =

_{n}*D*then is the dimension of vector

_{w}**w**

*(equal to 104 in this study). Then,*

_{m}*λ*represents the scalar of the magnitude of reliability of VFs obtained from the

_{m}*m*

^{th}eye. We assumed the data,

**t**

*, were independently drawn from a Gaussian distribution with mean vector*

_{n}**Φ**(

*x*)

_{n}^{T}

**w**

*and inverse of covariance matrix*

_{m}**L**

*is a 52 by 52 matrix. Consequently, the likelihood is given by:*

_{m}**L**

*is not constrained to the diagonal.*

_{m}**w**

*,*

_{m}*λ*, and

_{m}**L**

*were random variables that followed a Gaussian mixture distribution, a Gamma mixture distribution, and a Wishart mixture distribution, respectively. The Gaussian mixture model*

_{m}^{25}is a type of soft clustering method and can represent the presence of subgroups in a multimodal population. Unlike hard clustering methods, such as K-means methods, the Gaussian mixture model allows the derivation of a probability density distribution, and hence, it is suitable to be integrated in a probabilistic statistical model; for instance, the certainty or reliability of the estimated VF sensitivity can be inferred. The likelihoods were given by: and where

*z*∈{0,1},

_{mk}*ζ*∈{0,1},

_{mh}*γ*∈{0,1},

_{mg}*η*

_{h}∈[0,1],

*π*

_{k}∈[0,1],

_{g}∈[0,1],

*K*,

*H*, and

*G*are the number of components in each mixture distribution,

*W*(·|·) denotes the Wishart distribution, and

*G*(·|·) denotes the gamma distribution. If we marginalize Equation 2 over the

*ζ*weighted by Equation 3, the Wishart mixture distribution is given. Similarly, is given by Equations 4 and 5, and is given by Equations 6 and 7. In this model, the slopes and intercepts of eyes were clustered during the training phase because

_{mh}*p*(

**w**

*|*

_{m}*λ*) was a Gaussian mixture model. The full joint density function of the data

_{m}**T**≜ {

**T**

_{1},…,

**T**

*} and all the latent variables were given by: where*

_{M}*M*is the number of eyes. We estimated the hyperparameters

**ν**,

**W**,

**μ**,

**Λ**,

**a**,

**b**,

**η**,

**π**, and

**θ**, where

**ν**≜ {

*ν*

_{1},…,

*ν*},

_{H}**W**≜ {

**W**

_{1},…,

**W**

*},*

_{H}**μ**≜ {

**μ**

_{1},…,

**μ**

*},*

_{K}**Λ**≜ {

**Λ**

_{1},…,

**Λ**

*},*

_{K}**a**≜ {

*a*

_{1},…,

*a*},

_{g}**b**≜ {

*b*

_{1},…,

*b*},

_{g}**η**≜ {

*η*

_{1},…,

*η*},

_{h}**π**≜ {

*π*

_{1},…,

*π*}, and

_{k}**θ**≜ {

*θ*

_{1},…,

*θ*}, with the training data, using the E-M algorithm by maximizing the log-likelihood

_{g}^{26,27}for the E-M algorithm, variational approximation in the E-step of E-M algorithm, and variational Bayes method in the prediction phase. See the Supplementary Material for the detailed updating equations.

**t**

_{N}_{+1}be the (

*N*+ 1)

^{th}VF of a patient at time

*x*

_{N}_{+1}where

*N*is the number of VF tests already observed. Since we have not observed

**t**

_{N}_{+1}yet,

**t**

_{N}_{+1}was introduced as a latent variable. Hence, the complete likelihood was:

*p*(

**t**

_{N}_{+1}|

**T**) is calculated approximately using the variational Bayes method.

*D*= 52 dimensional column vector filled with ones. Predictive accuracy for mTD was compared using absolute errors.

_{t}^{28}

*P*< 0.05, Table 3) Furthermore, the RMSE associated with the VBLR approach using only the initial two VFs (VF2-3) was smaller than that with seven VFs (VF2-8) using the OLSLR conventional method. Similar differences in absolute errors were observed for each of the 52 points in the VF (Fig. 2).

**Figure 1**

**Figure 1**

**Figure 2**

**Figure 2**

**Table 3**

**Table 3**

OLSLR | VBLR | PValue | |

VF2-2 | – | 5.3 (±2.8) | |

VF2-3 | – | 5.0 (±2.7) | |

VF2-4 | 19.9 ± (12.0) | 4.9 (±2.6) | <1e-10 |

VF2-5 | 12.4 ± (7.1) | 4.7 (±2.5) | <1e-10 |

VF2-6 | 8.9 ± (4.6) | 4.5 (±2.4) | <1e-10 |

VF2-7 | 6.8 ± (3.4) | 4.4 (±2.3) | <1e-10 |

VF2-8 | 5.6 ± (2.8) | 4.2 ± (2.2) | <1e-10 |

VF2-9 | 4.7 ± (2.2) | 4.1 ± (2.1) | <1e-10 |

VF2-10 | 4.1 ± (2.0) | 3.9 ± (2.1) | 2.7e-7 |

*P*< 0.05). The RMSE with VBLR using only the initial two VFs (VF2-3) was smaller than that using six VFs (VF2-7) and OLSLR.

**Figure 3**

**Figure 3**

**Table 4**

**Table 4**

OLSLR | VBLR | PValue | |

VF2-2 | – | 2.2 ± (2.2) | |

VF2-3 | – | 2.0 ± (2.0) | |

VF2-4 | 6.2 ± (6.6) | 1.9 ± (2.0) | <1e-10 |

VF2-5 | 4.0 ± (4.1) | 1.8 ± (1.9) | <1e-10 |

VF2-6 | 2.9 ± (3.1) | 1.7 ± (1.7) | <1e-10 |

VF2-7 | 2.1 ± (2.2) | 1.5 ± (1.7) | <1e-10 |

VF2-8 | 1.8 ± (1.9) | 1.5 ± (1.6) | 3.0e-7 |

VF2-9 | 1.4 ± (1.5) | 1.4 ± (1.5) | 0.17 |

VF2-10 | 1.3 ± (1.4) | 1.2 ± (1.3) | 0.13 |

^{29}while other studies suggested that even longer series are necessary.

^{30,31}Under the assumption that the VF is measured every six months, it would take, on average, roughly 3 to 4 years to obtain sufficiently accurate results; however, many patients could suffer detrimental VF loss in this period of time. Other reports recommend that six VF tests should be obtained in two years,

^{32,33}but it is often unrealistic to adhere to this testing frequency in all patients

^{34}due to limited clinical capacity, especially in regions where ophthalmic facilities are limited. This highlights the clinical usefulness of the VBLR approach: with two or even a single VF test, the VBLR approach can predict future VFs as accurately as OLSLR with a series of six VF tests. In other words, with the VBLR approach, it is possible to predict what a patient's VF will look like in 3 to 4 years' time using only their baseline VF, while OLSLR would require a further 6 to 8 VF tests to make a prediction with equivalent accuracy to VBLR. Consequently, with the VBLR approach, we could identify a patient with rapid progression much earlier than we used to using OLSLR. Therefore, clinical decisions could be made much earlier helping to save patients' vision related quality of life. In addition, earlier identification of patients with slow progression could spare the expense of unnecessary treatments.

^{35}; however, the learning effect can persist over many VFs.

^{36}Hence, we also compared the prediction accuracy of OLSLR and VBLR, excluding the first two or three VFs. The VBLR continued to outperform OLSLR, with a magnitude similar to the results shown in Figures 1 through 3. We also compared the prediction accuracy of OLSLR and VBLR in POAG and NTG groups separately, since previous studies have reported that the pattern of VF deterioration differs between these disease groups.

^{37–39}Our results (not shown) suggested there was no significant difference in prediction accuracy between the groups.

^{40–42}One of the limitations of the current study is that this type of information, along with IOP, was not incorporated into the current VBLR approach. However, the model could be extended to incorporate such information by introducing a mixture of expert models—this should be carried out in a future study to attempt to further improve prediction accuracy.

^{43}In the current analysis, VF series often belong to different clusters at different time points, which is inevitable because the clustering depends on spatial and temporal patterns. The purpose of the current study was to improve VF predictions, not to cluster VFs; however, it would be interesting to analyze the relationship between VF patterns and predictive ability; future studies are merited to shed light on this issue. Furthermore, we did not investigate how varying parameters K, H, and G in the VBLR model influenced the performance of VBLR; this was simply to avoid over-fitting. A future study could investigate this question, however, a further dataset would be require to infer the optimum values without over-fitting.

^{6}recently reported that an exponential model might be more suitable than OLSLR to measure pointwise VF sensitivity decay. We initially applied exponential regression to our data, but prediction accuracy was not as good as expected, compared to that with OLSLR (results not shown). There is a possibility that the exponential model fits better in a longer series of VFs since follow-up was 8 years on average in the report of Caprioli et al.

^{6}Russell and Crabb

^{44}proposed Tobit regression as an alternative to the exponential model, since Tobit regression respects the floor effect of VF sensitivity data; however, the Tobit model is not suitable for the variational Bayes framework due to its nonlinear optimization and the resulting problem of local optima. The exponential model could be applied using a local variational method,

^{45}but this would require applying further approximation. On the other hand, Bryan et al.

^{10}recently concluded that, despite violation of its assumptions, the classical uncensored OLSLR model provides the best fit for VF data and performs best at predicting future VFs when compared to the Tobit and exponential regression models. Moreover, they stated that more advanced regression models exploring the temporal–spatial relationships of glaucomatous progression are required to reduce prediction errors to clinically meaningful levels.

**H. Murata**, P;

**M. Araie**, None;

**R. Asaoka**, P

*. 2006; 90: 262–267. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 2003; 14: 100–105. [CrossRef] [PubMed]*

*Curr Opin Ophthalmol**. 2003; 87: 726–730. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 2012; 53: 6776–6784. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 1996; 80: 40–48. [CrossRef] [PubMed]*

*Br J Ophthalmol**. 2011; 527: 4765–4773. [CrossRef]*

*Invest Ophthalmol Vis Sci**. 1986; 101: 1–6. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 2013; 54: 5505–5513. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 1995; 233: 750–755. [CrossRef]*

*Graefe's Arch Clin Exp Ophthalmol**. 2013; 54: 6694–6700. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 2012; 53: 2390–2394. [CrossRef] [PubMed]*

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

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

*Invest Ophthalmol Vis Sci**. 1969; 8: 84–91. [PubMed]*

*Invest Ophthalmol**. 2011; 118: 1782–1789. [CrossRef] [PubMed]*

*Ophthalmology**. 2006; 50: 38–43. [CrossRef] [PubMed]*

*Jpn J Ophthalmol**. 1993; 116: 684–691. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 2001; 10: 121–128. [CrossRef] [PubMed]*

*J Glaucoma**. 1993; 100: 69–75. [CrossRef] [PubMed]*

*Ophthalmology**. 2002; 43: 162–169. [PubMed]*

*Invest Ophthalmol Vis Sci**. 2013; 35: 2765–2781. [CrossRef]*

*IEEE Trans Pattern Anal Mach Intel**. 2012; 153: 1197–1205. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 2012; 53: 2199–2207. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**. 2014; 9: e85654. [CrossRef] [PubMed]*

*PLoS One**Finite Mixture Models*. New York, NY: John Wiley & Sons; 2004.

*. 2008; 95: 012015. [CrossRef]*

*J Phys Conf Ser**. 1998; 11: 271–282. [CrossRef] [PubMed]*

*Neural Netw**. 1985; 173 (suppl): 19–21.*

*Acta Ophthalmol**. 1982; 60: 267–274. [CrossRef]*

*Acta Ophthalmol**. 2000; 41: 2192–200. [PubMed]*

*Invest Ophthalmol Vis Sci**. 2010; 94: 1404–1405. [CrossRef] [PubMed]*

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

*Br J Ophthalmol**. 2013; 3: e002067.*

*BMJ Open**. 1996; 114: 19–22. [CrossRef] [PubMed]*

*Arch Ophthalmol**. 2008; 85: 1043–1048. [CrossRef] [PubMed]*

*Optom Vis Sci**. 1995; 233: 610–616. [CrossRef]*

*Graefe's Arch Clin Exp Ophthalmol**. 1980; 24: 621–664. [CrossRef] [PubMed]*

*Surv Ophthalmol**. 1984; 97: 730–737. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 1989; 108: 636–642. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 1986; 102: 402–404. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 1987; 104: 577–580. [CrossRef] [PubMed]*

*Am J Ophthalmol**. 1991; 3: 79–87. [CrossRef]*

*Neural Comput**. 2011; 52: 9539–9540. [CrossRef] [PubMed]*

*Invest Ophthalmol Vis Sci**Pattern Recognition and Machine Learning*. Vol. 1. New York, NY: Springer; 2006.