**Purpose**:
We develop a longitudinal statistical model describing best-corrected visual acuity (BCVA) changes in anti-VEGF therapy in relation to imaging data, and predict the future BCVA outcome for individual patients by combining population-wide trends and initial subject-specific time points.

**Methods**:
Automatic segmentation algorithms were used to measure intraretinal (IRF) and subretinal (SRF) fluid volume on monthly spectral-domain optical coherence tomography scans of eyes with central retinal vein occlusion (CRVO) receiving standardized anti-VEGF treatment. The trajectory of BCVA over time was modeled as a multivariable repeated-measure mixed-effects regression model including fluid volumes as covariates. Subject-specific BCVA trajectories and final treatment outcomes were predicted using a population-wide model and individual observations from early follow-up.

**Results**:
A total of 193 eyes (one per patient, 12-month follow-up, 2420 visits) were analyzed. The population-wide mixed model revealed that the impact of fluid on BCVA is highest for IRF in the central millimeter around the fovea, with −31.17 letters/mm^{3} (95% confidence interval [CI], −39.70 to −23.32), followed by SRF in the central millimeter, with −17.50 letters/mm^{3} (−31.17 to −4.60) and by IRF in the parafovea, with −2.87 letters/mm^{3} (−4.71 to −0.44). The influence of SRF in the parafoveal area was −1.24 letters/mm^{3} (−3.37–1.05). The conditional *R*^{2} of the model, including subject-specific deviations, was 0.887. The marginal *R*^{2} considering the population-wide trend and fluid changes was 0.109. BCVA at 1 year could be predicted for an individual patient after three visits with a mean absolute error of six letters and a predicted *R*^{2} of 0.658 using imaging information.

**Conclusions**:
The mixed-effects model revealed that retinal fluid volumes and population-wide trend only explains a small proportion of the variation in BCVA. Individual BCVA outcomes after 1 year could be predicted from initial BCVA and fluid measurements combined with the population-wide model. Accounting for fluid in the predictive model increased prediction accuracy.

^{1}In clinical practice, the functional response to anti-VEGF treatment is markedly heterogeneous and ranges from substantial vision gains with subsequent stability, to loss of vision despite aggressive therapy.

^{2}The individual functional potential to respond to VEGF inhibition is determined by several factors, including baseline visual acuity (VA) and variables obtained by high-resolution imaging such as optical coherence tomography (OCT).

^{3}While some changes observed within the retina are amenable to treatment (such as exudative fluid), others may become permanent despite therapeutic intervention (e.g., photoreceptor damage).

^{3}However, the extent of each factor's contribution is unclear and the interplay between imaging biomarkers, visual function, and treatment response remains poorly understood. Conventional statistical methods reporting outcomes (e.g., mean change in best-corrected visual acuity [BCVA]) cannot describe the longitudinal trajectories of visual function or the relation between vision and retinal morphology.

^{4}In this study, we integrated longitudinal clinical and imaging data to model the influence of pathologic changes in the retina on the clinical outcome variable and to predict individual future therapeutic outcomes based on the current state. We demonstrated our method on automatically segmented OCT characteristics and BCVA letter score as a clinical variable in a large population of patients with macular edema secondary to central retinal vein occlusion (CRVO) receiving anti-VEGF therapy (Fig. 1). We assumed that the individual VA and its development during therapy are affected by two major factors: reversible damage from which the retina recovers over time (where fluid in the retina is the major contributor) and irreversible damage. Hence, we modeled the population-wide temporal BCVA trajectory as a function of time and retinal fluid. By introducing random effects into the model, we accounted for the subject-specific deviation from the population-wide BCVA trajectory due to variation in reversible and irreversible damage as well as speed of recovery. We demonstrated that the combination of imaging and clinical data allows accurate modeling as well as prediction of treatment outcomes.

^{5}Therefore, only procedures immediately relevant to the current analysis are described herein. The study was conducted in compliance with the Declaration of Helsinki and all participants provided written informed consent before inclusion. Approval was obtained from the ethics committee at the Medical University of Vienna and each study site participating in the CRYSTAL trial.

*n*= 193, 13 visits) were included in our retrospective analysis of the prospective study. All patients underwent standardized SD-OCT imaging by Cirrus HD-OCT (

*n*= 58; Carl Zeiss Meditec, Dublin, CA, USA) and Spectralis OCT (

*n*= 135; Heidelberg Engineering, Dossenheim, Germany) at monthly visits by certified operators. A 6 × 6 mm macular cube scan pattern with a resolution of 200 × 200 for Cirrus and 512 × 49 for Spectralis was used. A total of 89 scans that failed preprocessing steps described below due to poor image quality were removed from the analysis (see Fig. 2 for examples). Overall, 2420 scans were analyzed.

^{6}First, follow-up patient scans were aligned rigidly to match the segmented vessel structures of the two images. Interpatient affine registration was obtained by aligning the center of the retina (fovea) and optic nerve head (ONH) center within patients (Fig. 3A).

^{7}(Fig. 3B). The segmentations were transformed in the joint coordinate system and the total volume of IRF and SRF was computed within the central 1 mm region around the fovea and parafoveal region in the 1- to 3-mm radius ring (Fig. 3C). In the following, we denote the amount of foveal fluid volume measured as

^{8}with a quadratic growth term. BCVA was the outcome variable, time and fluid volume covariates were fixed effects, and individually varying intercept and slope were modeled as random effects. The BCVA value of a single visit

^{8}

^{2}, fluid volumes) as columns, and visits as rows. The coefficients of the fixed effects are subsumized in the

^{9}where ML is known to produce biased estimates. On the other hand, doing likelihood-based model comparisons does not provide meaningful results for REML models with different fixed effects structures.

^{10}Thus, in our work we used the ML estimates for likelihood-based model selection and comparison only, and the REML estimates otherwise.

^{11}

*t*using Equation 4, and finally the BCVA using the predicted volumes in Equation 1.

^{12,13}

*t*is then the interval between the

_{ij}^{9}from the statistics software R. We used the EM estimate to compare the models.

^{14}and the Akaike information criterion (AIC).

^{15}As a goodness-of-fit measure, we used the coefficient of determination (

*R*

^{2}) for MRMs,

^{16,17}which is divided into the marginal

*R*

^{2}that describes the proportion of variance explained by the fixed factors only and the conditional

*R*

^{2}that considers fixed and random effects, thus measuring the variation explained by all effects. The coefficients were determined using REML for the final estimates.

*T*-statistics and

*P*values were calculated based on Satterthwaite's approximations and 95% confidence intervals (CIs) were obtained by bootstrap sampling with 500 simulations.

*R*

^{2}, which measures how well unseen samples are likely to be predicted by the model, and the mean absolute error (MAE) from the observed and predicted BCVA values. Because patients may suffer from recurring fluid at month 12 causing a drop in the VA, we used the median BCVA from months 10 to 12 as the prediction target to obtain a more reliable treatment outcome variable that is less affected by temporary events. Trajectories for each patient in the test-sets were predicted as described in Section 2.3 by first estimating the model parameters

^{18,19}with the alternative hypothesis that the prediction results of the nested model is less accurate than the full model in terms of MAE.

*R*

^{2}indicate that adding fluid volume improves the general model fit. The conditional

*R*

^{2}of 0.887 indicates a general good model fit. Most of the variance is covered by the subject-specific random intercept and slope, while a marginal

*R*

^{2}of 0.109 shows that only a rather small proportion of the variance in the data is explained by the fixed effects (intercept, slopes, and fluid volumes).

^{3}, which causes a drop of 5.30 letters, whereas the median

^{3}causes a drop of 1.38 letters. Similarly, the baseline median

^{3}results in a drop of 0.98 letters and

^{3}reduces BCVA by 0.17 letters.

*R*

^{2}is given in Figure 5. Furthermore, we exemplarily visualize the individual BCVA trajectories and prediction intervals for four patients in Figure 6, where the four models are compared in Figures 6A and 6B. How the predicted trajectory adapts to an increasing number of time points available for prediction is shown in Figures 6C and 6D.

^{20}provide a rich statistical framework for approaching these issues. The linear relations between a univariate longitudinal measurement and time points of observations as well as population-wide and individual effects can be modeled by using the particular random slope and intercept model proposed by Laird and Ware.

^{8}Intersubject variability is accounted for by allowing random deviations from the population-wide mean trajectory on an individual patient level by incorporating random intercepts and slopes. The superior statistical power of MRM compared to alternative methods using longitudinal imaging data of neurodegenerative diseases has been demonstrated.

^{21}Furthermore, MRM models have been shown to be capable of predicting the future development of early infant brain maturation by estimating the growth trajectories from magnetic resonance diffusion tensor images.

^{13}In previous work,

^{6}we used SD-OCT data combined with sparse logistic regression to predict disease recurrence from early time-point observations. However, this method is limited to balanced data without missing time points.

*R*

^{2}= 0.1), and foveal fluid and IRF have more influence on VA than parafoveal fluid and SRF. Thus, effects other than fluid cause a large fraction of vision. We subsumed these effects in the random effects that models baseline irreversible damage as subject-specific intercept deviation and reversible damage as subject-specific slope deviation. Note that we model only the immediate effects of fluid on BCVA explicitly, and not the secondary effects caused by fluid, such as damage to photoreceptor cells

^{22,23}or conductive elements of the retina.

^{24}

^{25}and other biomarkers that may be quantified in the future. The limited explanatory value of the fluid parameters on VA emphasizes the necessity of incorporating additional possible biomarkers. The proposed model has the capability to assess and quantify the influence of these biomarkers on VA similarly to the fluid model, by marginal/conditional

*R*

^{2}, and coefficient effect size and significance, and ultimately it may help to obtain additional insights in the disease, its progression, and outcome. Moreover, the specific effects of different pharmacologic agents or other treatments may be compared efficiently. The prediction model also could be applied to other diseases, such as diabetic maculopathy or age-related macular degeneration.

**W.-D. Vogl**, None;

**S.M. Waldstein**, None;

**B.S. Gerendas**, None;

**T. Schlegl**, None;

**G. Langs**, None;

**U. Schmidt-Erfurth**, None

*. 2016; 123: S78–S88.*

*Ophthalmology**. 2014; 121: 1092–101.*

*Ophthalmology**. 2016; 50: 1–24.*

*Prog Retin Eye Res**. Hoboken, NJ: John Wiley & Sons; 2011.*

*Applied Longitudinal Analysis*. 2nd ed*. 2016; 123: 1101–1111.*

*Ophthalmology**. 2015; 24: 152–163.*

*Inf Proc Med Imaging**. 2015: 24: 437–448.*

*Inf Proc Med Imaging**. 1982; 38: 963–974.*

*Biometrics**. 2015; 67: 1–48.*

*J Stat Statistical Software**. New York: Springer Science & Business Media; 2009.*

*Linear Mixed Models for Longitudinal Data**. 1991; 6: 15–32.*

*Stat Sci**. New York: Marcel Dekker, Inc.; 1997.*

*Linear and Nonlinear Models for the Analysis of Repeated Measurements**Medical Image Computing and Computer-Assisted Intervention – MICCAI 2014: 17th International Conference*.

*Springer International Publishing;*2014: 33–40.

*. 1946; 2: 110–114.*

*Biometrics**. 1974; 19: 716–723.*

*IEEE Trans Auto Control*^{2}from generalized linear mixed-effects models.

*. 2013; 4: 133–142.*

*Methods Ecol Evol**. 2014; 5: 944–946.*

*Methods Ecol Evol**. 1995; 13: 253–263.*

*J Bus Econ Stat**. 1997; 13: 281–291.*

*Int J Forcastt**. New York: Springer Science & Business Media; 2006.*

*Mixed-Effects Models in S and S-PLUS**. 2013; 66: 249–260.*

*Neuroimage**. 2011; 249: 183–192.*

*Graefe's Arch Clin Exp Ophthalmol**. 2009; 29: 1242–1248.*

*Retina**. 2011; 52: 2741–2748.*

*Invest Ophthalmol Vis Sci**. 2017; 177: 116–125.*

*Am J Ophthalmol*