**Purpose**:
The purpose of this study was to estimate the distribution of the true rates of progression (RoP) of visual field (VF) loss.

**Methods**:
We analyzed the progression of mean deviation over time in series of ≥ 10 tests from 3352 eyes (one per patient) from 5 glaucoma clinics, using a novel Bayesian hierarchical Linear Mixed Model (LMM); this modeled the random-effect distribution of RoPs as the sum of 2 independent processes following, respectively, a negative exponential distribution (the “true” distribution of RoPs) and a Gaussian distribution (the “noise”), resulting in a skewed *exGaussian* distribution. The *exGaussian*-LMM was compared to a standard Gaussian-LMM using the Watanabe-Akaike Information Criterion (WAIC). The random-effect distributions were compared to the empirical cumulative distribution function (eCDF) of linear regression RoPs using a Kolmogorov-Smirnov test.

**Results**:
The WAIC indicated a better fit with the *exGaussian*-LMM (estimate [standard error]: 192174.4 [721.2]) than with the Gaussian-LMM (192595 [697.4], with a difference of 157.2 [22.6]). There was a significant difference between the eCDF and the Gaussian-LMM distribution (*P* < 0.0001), but not with the *exGaussian*-LMM distribution (*P* = 0.108). The estimated mean (95% credible intervals, CIs) “true” RoP (−0.377, 95% CI = −0.396 to −0.359 dB/year) was more negative than the observed mean RoP (−0.283, 95% CI = −0.299 to −0.268 dB/year), indicating a bias likely due to learning in standard LMMs.

**Conclusions**:
The distribution of “true” RoPs can be estimated with an *exGaussian*-LMM, improving model accuracy.

**Translational Relevance**:
We used these results to develop a fast and accurate analytical approximation for sample-size calculations in clinical trials using standard LMMs, which was integrated in a freely available web application.

^{1}

^{–}

^{4}Because of the progressive and irreversible nature of the disease, VF sensitivity is only expected to worsen, resulting in a negative RoP, or remain unchanged over time. However, when the RoP is calculated through linear regression, the distribution of RoPs in clinical populations

^{5}

^{–}

^{11}is negatively skewed, but many patients show a positive RoP. Finding methods to estimate the underlying distribution of “true” RoPs would help the interpretation of the results of population-based studies and, importantly, of RCTs quantifying the effect of different treatment strategies.

^{9}

^{,}

^{10}proposed a hypergeometric-secant distribution for the RoPs in a clinical glaucoma cohort. Zhang et al.

^{11}and Swaminathan et al.

^{12}used a log-Gamma distribution for a similar purpose in their modeling of VF progression, where the log-Gamma distribution was used to model the distribution of the random effects in a hierarchical linear mixed effect model (LMM). Although these distributions fit the data by capturing some essential features (mainly, the skewness), they fail to give a description for a plausible underlying process of VF progression. In fact, a purely negative distribution of true RoPs is expected to arise to a negatively skewed distribution with positive values when measured in the presence of noise (see an example in Fig. 1). Additional factors, such as patients’ learning,

^{13}can contribute to the presence of observed positive RoP slopes and contaminate the true distribution. This relationship between true distribution and its noisy estimate from data has been already highlighted with simulations by Andrew Anderson.

^{9}

^{2}

^{,}

^{4}

^{14}

^{–}

^{16}All patient data were anonymized upon data extraction and transferred to a single secure database at City, University of London. Subsequent analyses of the data were approved by a research ethics committee of City, University of London. The study adhered to the Declaration of Helsinki and the General Data Protection Regulation of the European Union. The VFs included in this analysis were 24-2 tests performed with an Humphrey Field Analyzer (HFA), Goldmann III stimulus size and the Swedish Interactive Testing Algorithm (SITA Standard or SITA Fast). The whole database was composed of 576,615 VFs from 71,361 people, obtained between April 2000 and March 2015. VFs with a percentage of false positive errors ≥ 15% were excluded. No exclusion criteria were applied on fixation losses or false negative errors, following evidence from the literature.

^{17}

^{,}

^{18}Additional metrics, such as gaze tracking data, were not available in this dataset and could not be used to assess reliability. We included all patients with at least 10 VF tests performed over at least 4 years in one or both eyes and an MD worse than −2 decibels (dB) in at least 2 (not necessarily consecutive) VFs

^{19}

^{–}

^{21}in the same eye. This was used as a surrogate for the diagnosis of glaucoma, in the absence of a definitive label. Subjects with this level of damage and length of follow-up were likely to be either strong glaucoma suspects or persons with glaucomatous optic neuropathy. However, VF loss from other causes, such as vascular or neurological issues, could not be definitively excluded. VFs performed after any ocular surgery other than cataract were also excluded. Finally, only one eye from each patient was selected, at random if both were eligible, leaving 44,371 VFs from 3352 eyes. Because only one eye was included, the term eye and subject will be used interchangeably. This set of eyes has been used previously for other analyses.

^{4}

^{,}

^{22}

*exGaussian*). The formula for the PDF and the cumulative density function (CDF) of this distribution are reported in the Appendix. In simpler terms, each observed RoP can be thought of as the sum of a random draw from two independent distributions, an exponential (the “true” rate) and a Gaussian (the noise). Therefore, the mean of the resulting

*exGaussian*distribution is equal to the sum of the means of the two underlying exponential and Gaussian distributions. Similarly, the

*exGaussian*variance is equal to the sum of the two variances. Note that the PDF of the distribution of the sum of two random variables is obtained as the convolution of their PDFs. The

*exGaussian*PDF is the convolution of a Gaussian and exponential PDFs.

*exGaussian*(for reasons that will be explained), it is useful to understand whether it can be a practical approximation.

*exGaussian*distribution if the error is assumed Gaussian. The second component is intersubject variability, that is, not all patients have a “true” MD of 0 dB before developing glaucoma. If the distribution of “true” healthy MD values is also Gaussian, adding this element of noise would also result in an

*exGaussian*distribution. Finally, a third component is the actual time interval between the development of glaucoma and the first VF test. This is obviously unlikely to be the same for all subjects, as assumed in the ideal example. This interval can be described by a strictly positive random variable, such as one following a

*Gamma*distribution. The product of a Gamma and an exponential distribution is not an exponential distribution,

^{23}but would still result in a single-tailed negative distribution which can be approximated with an exponential distribution with mean equal to the product of the means of the original exponential and Gamma distributions. The exact PDF for the product distribution can be expressed analytically, but its parameters are difficult to reliably estimate, especially in the presence of additional Gaussian noise. A comparison between the exponential approximation and the exact PDF in a simulated example is provided in the Appendix.

*exGaussian*distribution for the intercept appears to be a reasonable approximation. For our main analyses, we chose a simplified model that estimated the Gaussian noise only from the SE of the intercept (see the next paragraph). An alternative model, that included the effect of intersubject variability, was also tested, but did not provide any improvement, because the estimated intersubject variability converged toward zero during fitting.

^{24}) within the R environment (R Foundation for Statistical Computing, Vienna, Austria). The LMM modeled the trend of MD over time, with two levels in the hierarchy, the population level and the subject level, because only one eye per subject was included. At the subject level, the model estimated an intercept and slope parameter for each subject. Like in standard Gaussian-LMMs, the response variable was assumed to be normally distributed and homoscedastic around the predicted value. This is a common assumption, although not necessarily accurately reflective of the data.

^{25}The choice of a Gaussian random effect distribution for the standard LMM does not have a strong theoretical justification, other than its relative simplicity, generalizability, and widespread use in the literature.

^{1}

^{,}

^{2}

^{,}

^{4}

^{,}

^{12}

^{,}

^{22}

^{,}

^{26}

^{–}

^{29}

*exGaussian*distribution. Instead of defining the

*exGaussian*PDF in the model, Bayesian computation allowed us to replicate the process underlying the generation of the observed data. This was achieved by sampling values from the exponential and Gaussian distribution separately and modeling the observed slope as the sum of the two values.

_{e}, i.e. the MD noise at the subject level). More details on the implementation of the Bayesian model are reported in the Appendix. One important aspect to note is that each subject in the dataset had a variable duration of follow-up time and number of tests. This introduced variation in the SE of the slope across subjects. Deriving the parameter sigma from the SE of the residuals allowed us to account for this, by calculating the expected SE of the intercept and slope for each subject (see Appendix). The data were also modeled with a standard Bayesian hierarchical LMM, using a Gaussian distribution for both intercepts and slopes. The two models were compared using the Wantanabe Akaike Information Criterion (WAIC) as implemented in the

*loo*package for R.

^{30}Note that, because of our specification of the

*exGaussian*-LMM, the number of estimated parameters is the same for both models, because the

*exGaussian*-LMM does not require an estimation of the between-subject level variance (see Appendix). This is, in fact, given by the sum of the variance of the exponential distribution (which is simply 1/λ

^{2}) and the variance of the Gaussian noise distribution (derived from the residual SE). However, the WAIC is influenced by the choice of prior distributions, which were partially constrained in our

*exGaussian*-LMM to improve stability with smaller datasets and shorter VF series (see Appendix). The model comparison was therefore based on a version of the

*exGaussian*-LMM for which the prior distributions were practically uninformative, so that the prior information was similar to that provided for the standard Gaussian-LMM (see Appendix). An additional comparison was performed by evaluating the empirical CDF of the RoPs and intercepts against the CDF of the estimated Gaussian and

*exGaussian*distributions with a Kolmogorov-Smirnov test.

^{31}

*exGaussian*-LMM on VF series that were progressively trimmed at the beginning, to observe whether a positive offset in the mean would reduce by eliminating earlier tests in the series, where the effect of learning would be the greatest. We kept the minimum series length at 10 VFs, as in the original analysis. The sample size was therefore reduced at each trimming step.

*P*value, derived from the

*P*direction.

^{27}

^{,}

^{28}

^{,}

^{32}Details are reported in the Appendix.

*exGaussian*-LMM were used to perform computer simulations of hypothetical treatment effects in RCTs. One advantage of our modeling is that it allows for a clear definition of a proportional treatment effect by scaling the parameter λ of the exponential distribution. This is different from previous attempts, where a change in slopes was somewhat poorly defined and achieved either by completely halting the progression of a proportion of patients

^{2}or by modifying observed RoP slopes by an additive factor to achieve the desired proportional change in the average RoP.

^{4}Another important difference is that previous attempts used observed regression slopes to model the “true” RoP, with simulated noise added to this “true” trend.

^{2}

^{,}

^{4}However, empirically calculated slopes will also be affected by noise, as explained previously. Instead, with our method, we simulated the “true” RoP slope for each eye as a draw from the “denoised” exponential distribution of slopes.

- 1. Perform a simple linear regression of MD over time for each eye on the real VF series.
- 2. From the linear regressions, calculate the SE of the residuals (σ
_{e}), an unbiased estimate of the individual variability. - 3. Select a sample size (
*N*) and a treatment effect (E). - 4. Randomly sample, without replacement,
*N*subjects for the placebo arm and*N*for the treatment arm. - 5. Generate “true” RoP values from an exponential distribution with rate λ for the placebo arm and λ/(1-E) for the treatment arm.
- 6. For each eye, generate a synthetic linear MD series for a time vector
*t*using the simulated RoP for the slope and the real baseline MD as the intercept. - 7. For each eye, add Gaussian noise using the subject specific variability calculated at point (2).
- 8. Fit an LMM and calculate the
*P*value of the interaction term between treatment and time (treatment effect) – Method 1 (see later). - 9. Fit simple linear regressions to the simulated data and perform an independent sample
*t*-test on the “observed” slopes comparing the two arms – Method 2 (see later). - 10. Repeat 1000 times from (4) to (9) for increasing
*N*and E.

*t*indicated the time of 16 tests over 2 years (with retest sessions), replicating the testing scheme of the United Kingdom Glaucoma Treatment Study (UKGTS) trial.

^{33}In reality, the variability of the MD varies with the level of damage.

^{25}However, global metrics, such as the MD, are predominantly influenced by individual performance.

^{4}

^{,}

^{34}

^{,}

^{35}Moreover, considering that we used the observed baseline MD as the intercept for our simulated series and that the simulated trial extended for only 2 years, it is reasonable to assume that the variability measured from each subject over at least 10 VFs (median of 11 years of follow-up) would be a realistic representation of the variability exhibited by that subject over such a hypothetical trial. Note that empirical approaches based on standardizing and permuting the observed residuals have the disadvantage of providing a biased estimate of variability (the SD of the residuals is smaller than the SE of the residuals).

*P*value for the difference in RoP between the two arms with an LMM using the

*lme4*

^{36}and

*lmerTest*

^{37}packages for R. The LMM had the MD as the response variable, with time (continuous) and treatment arm (factor) as fixed effect predictors. An interaction term between the treatment arm and time modeled the average difference in the RoP between the two arms. Random effects for both the intercept and the slope at the subject level were modelled as a bivariate Gaussian distribution, which included the correlation between the two parameters. Residuals were assumed conditionally independent. The

*lmerTest*package offers various methods to calculate

*P*values for LMMs, which differ in the way they calculate the degrees of freedom. A common choice is Satterthwaite's method, which provides an approximate t-distribution for the parameters.

^{37}

*P*value obtained for LMMs with Satterthwaite's method is well approximated by that calculated from a two-sample

*t*-test on the empirical RoPs obtained with simple linear regression. We confirmed this by performing a two-sample

*t*-test for all our simulations. We further exploited this property to obtain analytical approximations of the LMM power curves by calculating the expected variance and mean of the empirical exGaussian RoP distributions, based on the testing schedule, the average residual SE and the simulated exponential distribution for the “true” RoP. These values for the expected mean and variance were then used to analytically calculate the power of a

*t*-test for two samples with unequal variances (using the

*pwr*package for R).

*P*value for the treatment effect with each method was < 0.05. Confidence intervals for the power was calculated as 1.96 × SE

_{p}, where SE

_{p}is the standard error for the probability of a binomial process, calculated as \(S{E_p} = \sqrt {{p_{p < 0.05}}*( {1 - {p_{p < 0.05}}} )/N} \), where

*p*

_{p < 0.05}is the proportion of significant

*P*values in the simulations and

*N*is the total number of simulations. Naturally, the analytical approximation has the advantage of allowing quick adjustments of the parameters to accommodate for various trial designs, expected levels of variability and neuroprotective effects. We developed an interactive web application in the Shiny environment

^{38}to allow users to experiment with the effect of these parameters on the estimated power curves. The app is available at https://giovannimontesano.shinyapps.io/Sample_size/.

*exGaussian*-LMM (estimate and SE 192174.4 [SE = 721.2]) than with the Gaussian-LMM (192595 [SE = 697.4327], difference 157.2 [SE = 22.6]). This was confirmed by the Kolmogorov-Smirnov (KS) test for the slopes, which showed no significant difference between the empirical CDF of the RoPs and the estimated

*exGaussian*distribution (

*P*= 0.108), but a highly significant difference with the estimated Gaussian distribution (

*P*< 0.0001). A table comparing the WAIC for the two models for increasing series length (from 4 to 10 tests) is provided as Supplementary Material. For the distribution of the intercepts, the

*exGaussian*approximation was, as expected, not as good as for the observed RoPs (KS =

*P*< 0.0001). The empirical distribution showed a faster decay in the negative tail compared to the prediction. This is in agreement with theoretical expectations (see Appendix). Nevertheless, a model with a Gaussian distribution for the random effects on the intercepts provided a significantly worse fit (WAIC: 192473.3 [SE = 710.6]) compared to the full

*exGaussian*-LMM (difference: 149.4 [SE = 22.4]). The random effect estimates of intercepts and slopes obtained with the two LMMs are reported as Supplementary Material.

*t*-test on the empirical slopes in our simulated RCTs, for different neuroprotective effects and increasing sample sizes. The untreated “true” RoP was −0.38 dB/year (see the Table). The average residual SE of the MD was 1.97 dB. There was excellent agreement between the LMM and the

*t*-test, with the largest absolute difference in power being 0.1%. Disagreement in determining a significant effect was observed only in 16 of 40,000 simulations (0.04%).

*t*-test. The analytically predicted values were outside the 95% credible intervals (CIs) of the simulated values only in 4 simulations, one when simulating no effect, where the simulations resulted in a power lower than the expected 5% false positive rate, and 3 when simulating the 20% neuroprotective effect. These are highlighted in Figure 6. We also report a comparison with the power curves estimated with the residual SE estimated from test-retest data (see Supplementary Material), which are more likely to be reflective of test noise in controlled environments, such as those of RCTs.

^{9}

^{,}

^{10}

^{,}

^{12}particularly their skewness. This has important consequences. In terms of clinical interpretation, it shows that noise and learning can effectively mask the true rate of VF progression. However, this method allows the estimation of a plausible underlying distribution of “true” rates in actual clinical populations. This is important for the derivation of accurate population parameters. For example, Figure 4 shows how the average rate of progression would be biased by the effect of learning in a naïve calculation. Our interpretation of the parameters is supported by specific results in our analysis. In particular, progressively taking out the initial tests in the series showed that the mean of the Gaussian component of the estimated

*exGaussian*distribution (see Fig. 4) is likely modeling the average effect of learning. This effect came very close to 0 dB/year at the seventh test. This is in agreement with previous investigations,

^{13}which estimated that the effect of learning can extend to the sixth or seventh VF. The effect of learning on population estimates is also clearly demonstrated in our analysis by severity group (see Fig. 5). We show that the effect of learning, and the consequent bias on the average rate of progression, is stronger in patients with a more advanced baseline MD. One explanation for this result is the overestimation of baseline damage in series affected by a strong learning effect or regression to the mean, because patients who performed the worst in their first tests were more likely to be classified as having moderate or advanced baseline loss. Despite this, the “true” estimated RoP was significantly faster in patients with moderate or advanced baseline loss compared to patients with early damage, as expected. However, the difference was not significant between patients with moderate and advanced damage. This was likely due to inaccurate classification of baseline damage (again due to learning affecting the initial tests) and the perimetric floor effect in truly advanced patients, which can positively bias the RoP of global metrics.

^{22}

*exGaussian*distribution also offers straightforward physiological interpretations in the context of the progression of the disease, because it avoids the need to assume that the true rate of progression could be positive. Finally, the

*exGaussian*distribution effectively captures the skewness in the data and can be used in hierarchical models instead of the commonly used Gaussian distribution, reducing the effect of shrinkage towards the mean (see Supplementary Material). Such a shrinkage can bias the trajectory of fast progressing eyes toward the general trend of the population, especially in those for which a smaller number of VF tests is available. This is supported by the improvement in WAIC compared to a Gaussian-LMM and the lack of significant departure of the empirical CDF of the RoPs from the estimated

*exGaussian*, as shown by the KS test (

*P*= 0.108). Importantly, such a model requires the estimation of the same number of parameters as a Gaussian-LMM, as described in the Methods section. This is also another important difference from previous attempts using different skewed distributions, which instead required the skewness to be modelled with an additional parameter.

^{11}

^{2}to simply adding a positive effect to the observed RoPs to change the average rate by a predefined amount.

^{4}Modeling the underlying “true” rate with a negative distribution, such as a negative exponential, offers a univocal definition of percentage change in RoPs by simply scaling, with a multiplicative factor, the “true” distribution. This not only replicates the complex distributional change in the observed RoPs that would be observed with treatment, but also affects the individual rates in a manner that replicates the underlying physiological process. Moreover, in the absence of learning (Gaussian noise mean = 0), a proportional change in the “true” exponential distribution also results in the same proportional change of the resulting

*exGaussian*distribution. It is interesting to note that this is not true in the presence of learning. For example, in an RCT, the learning effect is likely to be the same in both arms. Therefore, the average difference in RoP, calculated, for example, with standard LMMs, would only depend on the average difference between the “true” RoPs (exponential). However, learning can positively bias the average RoP in both arms, and the same linear difference can result in various calculated “proportional” effects. This creates inconsistencies in trial design and in the interpretation and generalizability of the results, especially when naïve patients are recruited. For example, the treatment arm of the United Kingdom Glaucoma Treatment Study (UKGTS), in which true progression was likely to be slower, showed a positive median RoP (+0.12 dB/year), indicating a strong learning effect.

^{39}Modeling an

*exGaussian*distribution eliminates such an ambiguity because the linear difference in the average RoPs always corresponds to the same proportional change in the “true” distribution of rates.

^{4}Moreover, the

*exGaussian*distribution makes it straightforward to calculate the expected variance and mean of the empirical distribution of RoPs based on the expected test variability, length of the follow-up and planned testing schedule. This allows a simple analytical approximation of the results of standard LMMs with a two-sample

*t*-test, as explained in the Methods section and demonstrated in Figure 6, and makes it possible to translate our results into a user-friendly web application (https://giovannimontesano.shinyapps.io/Sample_size/). One important aspect to highlight is that the average residual SE estimated with linear regression from our sample was 1.97 dB, larger than the values reported in previous literature (closer to 1 dB).

^{34}

^{,}

^{40}This discrepancy is possibly explained by the fact that our sample is composed of real-life test series over a long period of time (≥ 10 years). On the one hand, therapeutic intervention might have introduced deviations from a simple linear trend over such a long follow-up. On the other hand, suboptimal testing conditions might have created larger test fluctuations, which would be lower in more controlled testing conditions, such as those of test-retest variability studies or prospectively planned data collections. Such controlled testing conditions are more likely to reflect the perimetric noise that would be found in RCTs. We provide, as Supplementary Material, an in-depth analysis of the combined data from two perimetric test-retest datasets.

^{22}

^{,}

^{39}

^{,}

^{41}

^{,}

^{42}The average test-retest SD from these datasets of MD was 0.94 dB, much closer to the literature, and was larger for patients with more advanced damage, as expected. Power calculations using these estimates are much more conservative and more likely to be appropriate for trial scenarios; these have been integrated in our web application and reported in Figure 6.

*exGaussian*to effectively describe the observed data, the VF series need to be relatively homogenous in terms of number of tests and duration of follow-up. This is because series with different characteristics are, by definition, not sampled from the same distribution, because they would have different expected variance for the noise component (the SE of the slope). However, the

*exGaussian*-LMM does account for such a heterogeneity in the estimation process, because the model calculates the expected SE of the slope based on the specific test series for each eye. This means that, in fact, each random slope is drawn from its own distribution, while simultaneously allowing the estimation of population parameters of interest. Naturally, there are cases where an exponential distribution might not reflect the underlying distribution of “true” rates. However, the exponential distribution seemed to provide a good fit for the observed data and allowed us to estimate an

*exGaussian*-LMM with the same number of parameters as a simple Gaussian-LMM. Moreover, when testing a more general Gamma distribution in preliminary fitting experiments, we found that the shape parameter converged to a value slightly smaller than 1, making the exponential distribution a good candidate for our implementation (a Gamma distribution with a shape parameter of 1 is a simple exponential). Another example is provided by a recruitment strategy for a hypothetical trial which selects participants using cutoffs on the RoP measured in the clinic. We provide a description of the changes that such a selection would bring to the “true” distribution of RoPs, assumed exponential, in the Appendix. The resulting distribution would be extremely difficult to model exactly in any practical implementation of the LMM. However, we also show that, for the typical level of perimetric variability and a realistic number of tests that would be available in the clinic (4 to 6), the noise in the observed RoPs is such that, in practice, the resulting “true” distribution would only deviate minimally from the original exponential. In other words, such a selection would be largely ineffective in identifying patients truly progressing at a rate within a desired range. Note that this result is valid regardless of the distribution assumed for the “true” RoPs, because it is simply a consequence of the uncertainty affecting the estimation of the true rates from the observed rates.

*exGaussian*, as explained in the Appendix. Modeling the expected distribution under a simulated scenario is relatively straightforward (see Appendix). However, estimating all the parameters for such a distribution from clinical data would be challenging and very likely imprecise. The estimation of a plausible distribution becomes achievable with some simplifications and assumptions, but, in our analyses, it did not provide any improvement in fitting the data over assuming an

*exGaussian*distribution for both intercepts and slopes (based on the WAIC, see Appendix). It should be noted that various choices of the distribution for the intercept, including a simple Gaussian, only had a minimal effect on the estimated distribution of the slopes, the actual parameter of interest.

^{22}and facilitates comparison with standard methodology. Individual variability could be incorporated in the model, for example, by accounting for the level of damage.

^{34}However, global indices such as MD are also largely influenced by global fluctuations in the performance,

^{4}

^{,}

^{34}

^{,}

^{35}which are often subject-dependent. This means that both a systematic and a subject-level element of variability would need to be modeled, considerably increasing the complexity of the model and the number of parameters required. This will be the objective of future work. However, the estimated distribution appeared to accurately describe the observed distribution of RoPs calculated via linear regression (see Fig. 3).

^{22}they pose significant challenges when implementing larger hierarchical structures for more complex distributions, such as in

*exGaussian*-LMMs. Further extension of the model will focus on modeling an additional level in the hierarchy, for example, by modeling the RoP at each location and VF cluster.

^{22}

^{,}

^{27}

^{,}

^{28}This is important for RCTs, because pointwise modeling has been beneficial in some retrospective analyses of trial data to highlight subtle differences

^{27}

^{,}

^{29}in progression. However, a better characterization of fast progressors with the

*exGaussian*-LMMs based on MD might prove sufficient to overcome the limitations of commonly used Gaussian-LMMs and provide similar power to pointwise LMMs. Other methods of analyzing RCT data, such as those based on machine learning,

^{43}

^{–}

^{45}might also provide additional statistical power. However, such empirical techniques have the disadvantage of providing results that are often not directly interpretable. In contrast, our methodology was designed to provide estimates for parameters that had a clear interpretation in terms of disease progression, measurement noise and learning artifact, making it valuable in elucidating the actual impact on disease modification and hopefully improve translation of the results across different studies and datasets.

*exGaussian*-LMM itself could be used for trials to estimate the change in the “true” RoP due to treatment. Although this could be easily implemented for individual trials, sample size calculations would be complicated due to the complex and time-consuming procedure required to estimate the parameters of the

*exGaussian*-LMM. This is mainly because of the relatively complex, but fully interpretable, structure of the model to separate the parameters of the exponential and Gaussian component of the distribution. Defining approximate solution for power calculations with the proposed model will be the objective of future work.

**G. Montesano**, CenterVue-iCare (C), Alcon (C), Relayer, LTD (O), Omikron (C);

**D.F. Garway-Heath**, Carl Zeiss Meditec (C), CenterVue-iCare (C), Heidelberg Engineering (F), Moorfields MDT (P), ANSWERS (P), T4 (P), Visual Field Sensitivity Testing (P), Omikron (C);

**D.M. Wright**, None;

**A. Rabiolo**, None;

**G. Ometto**, Relayer, LTD (O), Alcon (C);

**D.P. Crabb**, ANSWERS (P), T4 (P); Allergan (C), Apellis (C), Roche (C), Thea (C), Santen (F)

*Ophthalmology*. 2023; 130: 469–477. [CrossRef] [PubMed]

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

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

*Am J Ophthalmol*. 2021; 229: 127–136. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2014; 55(7): 4135–4143. [CrossRef] [PubMed]

*Arch Ophthalmol*. 2011; 129(5): 562–568. [CrossRef] [PubMed]

*Ophthalmology*. 2009; 116(12): 2271–2276. [CrossRef] [PubMed]

*Arch Ophthalmol*. 2010; 128(10): 1249–1255. [PubMed]

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

*Transl Vis Sci Technol*. 2015; 4(4): 2. [CrossRef] [PubMed]

*Biom J*. 2015; 57(5): 766–776. [CrossRef] [PubMed]

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

*Optom Vis Sci*. 2008; 85(11): 1043–1048. [CrossRef] [PubMed]

*BMJ Open Ophthalmol*. 2019; 4(1): e000352. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2019; 207: 144–150. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2020; 104(10): 1406–1411. [CrossRef] [PubMed]

*Ophthalmology*. 2017; 124(11): 1612–1620. [CrossRef] [PubMed]

*Acta Ophthalmol Scand*. 2000; 78(5): 519–522. [CrossRef] [PubMed]

*Ophthalmic Physiol Opt*. 2015; 35(2): 225–230. [CrossRef] [PubMed]

*Ophthalmic Physiol Opt*. 2017; 37(1): 82–87. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2012; 96(9): 1185–1189. [PubMed]

*Transl Vis Sci Technol*. 2021; 10(12): 4. [PubMed]

*Quality & Quantity*. 2013; 47(1): 545–552.

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

*Transl Vis Sci Technol*. 2023; 12(10): 20. [PubMed]

*Am J Ophthalmol*. 2023; 251: 143–155. [PubMed]

*Am J Ophthalmol*. 2023; 246: 42–50. [PubMed]

*Ophthalmology*. 2020; 127(10): 1313–1321. [PubMed]

*Stat Comput*. 2017; 27(5): 1413–1432.

*J Am Stat Assoc*. 1951; 46(253): 68–78.

*Front Psychol*. 2019; 10: 2767. [PubMed]

*Lancet*. 2015; 385(9975): 1295–1304. [PubMed]

*Transl Vis Sci Technol*. 2018; 7(3): 22. [PubMed]

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

*J Stat Software*. 2015; 67(1): 1–48.

*J Stat Software*. 2017; 82(13): 1–26.

*Health Technol Assess*. 2018; 22(4): 1–106. [PubMed]

*Invest Ophthalmol Vis Sci*. 2013; 54(2): 1345–1351. [PubMed]

*Trans Am Ophthalmol Soc*. 2017; 115: T4. [PubMed]

*Ophthalmology*. 2014; 121(10): 2023–2027. [PubMed]

*Sci Rep*. 2019; 9(1): 18113. [PubMed]

*PLoS One*. 2019; 14(4): e0214875. [PubMed]

*Invest Ophthalmol Vis Sci*. 2019; 60(1): 365–375. [PubMed]

*Stat Sci*. 1992; 7(4): 457–472.

*Int J Psychol Res*. 2010; 3(1): 68–77.

^{24}The LMMs modelled a linear relationship of MD over time. There were two levels in the hierarchy, the population level and the subject level. The hierarchical LMM modelled individual progression by using random intercepts and slopes. We will describe the implementation of the standard Gaussian LMM and then explain how this was modified to use an

*exGaussian*distribution for the random effects on the individual slopes.

_{0}and σ

_{0}) and the slopes (µ

_{1}and σ

_{1}). The individual intercepts and slopes where random draws from these Gaussian distributions. The residuals error of the predicted MD had a Gaussian distribution with mean 0 and SD σ

_{e}(the residual SE). Note that, in standard LMMs, the random intercepts and slopes are usually modeled as a joint bivariate Gaussian distribution, to model the correlation between intercepts and slopes. We chose to instead model independent Gaussian distributions for the random effects to allow a direct comparison with the

*exGaussian*-LMM, where the joint modeling is not easily implemented.

*exGaussian*-LMM had a similar structure to the Gaussian-LMM. The difference was in the distribution of the random effects, which followed an

*exGaussian*distribution. In this case, the random intercepts and slopes were considered to be the sum of a random draw from a negative exponential distribution and a Gaussian distribution (the noise component). The model estimated two population-level parameters for the

*exGaussian*distribution, the parameter λ for the exponential distribution and the mean (µ

_{n}) of the Gaussian noise component. The SD σ

_{n}of the noise was calculated as the SE of the intercept (\(\sigma _{{{{\rm{\hat \beta }}}_0}}^2\)) and slope (\(\sigma _{{{{\rm{\hat \beta }}}_1}}^2\)) defined as below (

*t*indicates the time vector,

*n*is the number of VF tests):

_{e}rather than estimating the SE for each subject. This also means that the method assumes all subjects to have the same “average” level of variability. Although this is a common assumption in LMMs, it might not be reflective of the data.

^{25}However, our results show that such an approximation does not have a large impact on the accuracy of our estimates (see Fig. 6). The model could be further expanded, in the future, to account for changes in variability due to the level of damage or individual variations, if this is required for specific contexts.

^{46}

*exGaussian*-LMM were set as partially informative for the main calculations. The prior distribution for the parameter λ was a Gamma distribution with shape 1 and rate 1/\(\widehat \lambda \). The value of \(\widehat \lambda \) was calculated using the method of moments for the

*exGaussian*distribution

^{47}from the empirical values of intercepts and slopes obtained through linear regression. This Gamma prior has mean = \(\widehat \lambda \) and variance = 1/\({\widehat \lambda ^2}\). Note that whereas the method of moments can be used to estimate the parameters of an

*exGaussian*distribution and can provide useful prior information, it would fail to establish a link between the observed values and the underlying generating random process, especially for the RoP. The fitting procedure was identical to the Gaussian-LMM.

*exGaussian*-LMM using the

*loo*package for R.

^{30}The package also provides the SE for the difference in WAIC (lower WAIC indicates a preferrable model). As previously explained, because the WAIC can be influenced by the choice of priors, we compared the Gaussian-LMM to an

*exGaussian*-LMM fitted with similarly weakly informative priors. Specifically, the prior distribution for the parameter λ was a Gamma with shape and rate 0.001 (the same as the distribution for the inverse variance of the random effects in the Gaussian-LMM).

*exGaussian*-LMM (192197.5 [725.7]).

*P*value was calculated according to Makowski et al.

^{32}The median of the posterior distribution is determined from the MCMC draws. According to whether the median is above or below 0, the proportion of MCMC draws above or below 0 is the

*P*-direction, respectively. The

*P*-direction can have values between 0.5 (distribution centered on 0) and 1. The Bayesian

*P*value statistic can then be calculated as p = 2 * (1–

*P-direction*)

^{32}and has a similar interpretation to a two-tailed frequentist

*P*value.

*exGaussian*distribution is obtained as the convolution of the Gaussian PDF

_{1}and β

_{2}are the inverse of the rates for the exponential and the Gamma distribution, α is the shape parameter of the Gamma distribution.

*BesselK*indicates a Bessel function of the second kind. These functions can be complex to calculate, but are implemented in most languages, including base R. The effect of perimetric noise and intersubject variability of MD can be modelled with a Gaussian distribution. Similarly to the modeling of the slope, the distribution for the observed intercepts can be obtained by convolving the “true” distribution of the intercepts with this Gaussian distribution. In the case of a simple exponential, the result is an

*exGaussian*. The convolution with the PDF reported above has no close form solution and can be obtained via numerical convolution.

*c*, with an estimated standard error for the observed rates \({\sigma _{{{{\rm{\hat \beta }}}_1}}}\), is:

*c*) can be calculated as: