LIMBARE: An Advanced Linear Mixed-Effects Breakpoint Analysis With Robust Estimation Method With Applications to Longitudinal Ophthalmic Studies

Abstract

Purpose: Broken stick analysis is a widely used approach for detecting unknown breakpoints where the association between measurements is nonlinear. We propose LIMBARE, an advanced linear mixed-effects breakpoint analysis with robust estimation, especially designed for longitudinal ophthalmic studies. LIMBARE accommodates repeated measurements from both eyes and over time, and it effectively addresses the presence of outliers.

Methods: The model setup of LIMBARE and the computing algorithm for point and confidence interval estimates of the breakpoint were introduced. The performance of LIMBARE and other competing methods was assessed via comprehensive simulation studies and application to a longitudinal ophthalmic study with 216 eyes (145 subjects) followed for an average of 3.7 ± 1.3 years to examine the longitudinal association between structural and functional measurements.

Results: In simulation studies, LIMBARE showed the smallest bias and mean squared error for estimating the breakpoint, with an empirical coverage probability of corresponding confidence interval estimates closest to the nominal level for scenarios with and without outlier data points. In the application to the longitudinal ophthalmic study, LIMBARE detected two breakpoints between visual field mean deviation (MD) and retinal nerve fiber layer thickness and one breakpoint between MD and cup-to-disc ratio, whereas the cross-sectional analysis approach detected only one and none, respectively.

Conclusions: LIMBARE enhances breakpoint estimation accuracy in longitudinal ophthalmic studies, and the cross-sectional analysis approach is not recommended for future studies.

Translational Relevance: Our proposed method and companion R package provide a valuable computational tool for advancing longitudinal ophthalmology research and exploring the association relationships among ophthalmic variables.

Introduction

Broken stick analysis is a powerful analytical tool for investigating piecewise linear associations between continuous variables, particularly in cases where nonlinear associations are present. The goal of the analysis is to accurately detect breakpoints at which the effect size of the association changes, allowing for piecewise linear association inference to be made before and after these breakpoints. In ophthalmology, broken stick analysis has been used to study the association between structural and functional measurements in glaucoma,^1–3 where varying levels of association between structural measurements, such as average retinal nerve fiber layer (RNFL) thickness, and functional measurements obtained with visual field tests^4–8 have been observed at different stages of disease severity. These differences may be due to the “floor effect” of structural measurements beyond which no further structural damage can be measured, leading to changes in the structural–functional association.⁹ 

Current statistical approaches for broken stick analysis in ophthalmology mainly employ cross-sectional approaches where observations from one eye per subject from a single visit are included in the analysis.⁸ The approach negates the possibility of maximizing available data by using information obtained from both eyes of the same subject and evaluating longitudinal processes such as disease trajectory. Alternatively, aggregates of measurements across all visits from both eyes are used without considering the inter-eye correlations and within-subject correlations between repeated measurements. Without considering the correlations of repeated measurements from two eyes within a subject, the breakpoint estimation may be misleading. 

To address these shortfalls, we propose a new and advanced linear mixed-effects breakpoint analysis with robust estimation (LIMBARE) that can provide accurate broken stick analysis for ophthalmic studies that assess association relationships between measurements longitudinally sampled from both eyes. Specifically, LIMBARE is built upon the segmented mixed model, which accounts for repeated measurements over time and between eyes, and additional confounders commonly encountered in human studies.^10,11 LIMBARE further incorporates the least trimmed squares (LTS) strategy^12,13 to mitigate the impact of outliers and provide a robust and accurate point estimate of breakpoints. We also propose a new standard error estimator for LIMBARE to improve the coverage probability of confidence interval (CI) estimates. With LIMBARE (method workflow is illustrated in Fig. 1), the point and CI estimates of breakpoints, together with the association effect sizes quantified as slopes of the piecewise linear mixed effects model, are reported. 

Figure 1.

 

Schematic framework of the LIMBARE method. The input data are in long format. LIMBARE is comprised of two parts: Algorithms A1 and A2 are used in each part to detect potential outliers and obtain robust point (Output 1) and 95% CI (Output 2) estimates of breakpoints, association slopes, and other unknown parameters. The output is piecewise association inferences between longitudinal variables, which can be applied to identify changes in association patterns for the study of variables of interest.

A simulation study was conducted to evaluate the performance of LIMBARE against competing methods. We also applied these methods to a longitudinal study examining the association between visual field mean deviation (MD) and optical coherence tomography (OCT) RNFL thickness, and between MD and cup-to-disc ratio (CDR), to examine the piecewise association relationship. Our results demonstrate that LIMBARE improves the accuracy of breakpoint estimation, and we recommend the use of this method in future ophthalmic studies. 

Methods

LIMBARE

We used a segmented mixed-effects model (SMM) to model the longitudinal, nonlinear association between two variables of interest. Following the notation of Muggeo,¹⁰ two variables, denoted by x_ijk and y_ijk, are observed from subject i (= 1, …, n) at visit time t_ijk (j = 1, …, n_i) for eye k (= “OD,” “OS”). The model also includes time-dependent covariates, denoted by z_ijk, such as a constant vector 1, visit time t_ijk, and other confounders (see Fig. 1 for an illustration of the data input). Our primary goal is to detect unknown breakpoints in the association between x_ijk and y_ijk. Although the specification of the SMM regression model is straightforward, estimating the unknown parameters is challenging due to a parameter identifiability problem. In the case of the standard linear mixed-effects model, the maximum likelihood estimation (MLE) approach is commonly used to estimate unknown parameters using the derived log-likelihood function.¹⁴ However, the strategy is not applicable for SMM because the corresponding log-likelihood function is not differentiable at the breakpoints. Moreover, MLE is sensitive to outliers, which are common in ophthalmic studies due to measurement error, reading error, and instrumental error, for example. To overcome these challenges, we propose a robust estimation approach that allows for the detection of multiple breakpoints, accounting for potential outliers. Specifically, we incorporate the least trimmed squares (LTS) technique in the estimating algorithms to reduce the impact of outliers in the data on parameter estimation.^15,16 We further propose a companion algorithm for estimating the standard errors (SEs) and CIs of unknown parameters, including breakpoints and their associated slopes. Figure 1 provides a visual overview of the input, algorithm, and output of our approach. For details regarding model specifications and algorithms, please refer to Section A of the Supplementary Material. 

Competing Methods

To evaluate the performance of our proposed method, we compared it with three alternative methods commonly used for estimating breakpoints in longitudinal studies. The first method is the SMM proposed by Muggeo.^10,17 SMM uses an iterative optimization algorithm to estimate the model parameters, but it does not account for outlier data points. Therefore, we also evaluated the effect of outliers on the estimates of breakpoints obtained with this method. The second method is the robust mixed-effects segmented regression (RMESR) model proposed by Zhou et al.¹⁸ RMESR combines the estimates from both segmented linear models and linear mixed-effects models, with the breakpoint estimates based on the segmented linear model even though the data come from a longitudinal design. To accommodate outliers in the dataset, RMESR uses the least trimmed squares technique. The third method is the conventional method, referred to as CROSS, in which longitudinal measurements are aggregated cross-sectionally and analyzed with a segmented linear regression model.¹⁰ The detailed description of each method can be found in Section B of the Supplementary Material. 

Simulation Setup

Our simulation was designed to imitate the structural and functional data from real longitudinal ophthalmic studies. For Scenario 1 (without outliers) and Scenario 2 (with outliers), we simulated two breakpoints (ϕ₁, ϕ₂) for the breakpoint variable x. We also included a time variable to record the time from baseline visit in years and considered random effects between visits of the same subject and between eyes within subjects. The outcome model for variable y is specified in Section C of the Supplementary Material. We simulated data from three or four visits and from two eyes for 60% of the total of 300 subjects. We conducted 1000 replications to evaluate the performance of our proposed method, LIMBARE, and compared it with the SMM, RMESR, and CROSS methods. To simulate outliers under Scenario 2, we randomly selected 5% of observations from the simulated dataset and added [max(y) – min(y)]/±a to the outcome variable y and [max(x) – min(x)]/±b to the breakpoint variable x, where a, b ∼ uniform(4, 20) were randomly sampled from uniform distribution with parameters 4 and 20. We used the finite sample bias and mean squared error (MSE), mean SEs, average length of 95% CIs, and empirical coverage probability (ECP) as comparison metrics for assessing the point and CI estimates of the model parameters, respectively. 

Longitudinal Ophthalmic Study

This longitudinal clinical study examined 216 eyes (163 open-angle glaucoma, 44 glaucoma suspects, and nine healthy eyes) from a cohort of 145 subjects aiming to investigate changes in ocular structure and function over time. The study was approved by the Institutional Review Boards of the University of Pittsburgh and NYU Langone Health, and all participants provided informed consent. Each subject underwent comprehensive ophthalmic examinations, including qualified visual fields and OCT optic nerve head (ONH) and macular scans, in at least five visits over an average period of 3.7 years (Table 1). We considered MD to be the functional parameter of interest, and RNFL thickness and CDR were considered to be the structural parameters of interest. We first conducted breakpoint analyses for each parameter in relation to time from baseline to explore potential nonlinear changes in these parameters over time. We then evaluated the piecewise association between MD and two structural measures: RNFL thickness and CDR. The proportion of data trimming (α) was set to 0.05 for the MD–RNFL association and 0.1 for the MD–CDR association. To control for potential confounding factors, such as time from baseline visit (in years), age at baseline, signal strength of scans, disc area when evaluating CDR, and eye, all breakpoint analysis models were adjusted accordingly. We used the Bayesian information criterion (BIC) to determine the optimal number of breakpoints for each association and each method. The model with the smallest BIC was selected, and the corresponding estimated breakpoints and slopes were reported. 

Table 1.

 

Characteristics of the Longitudinal Clinical Study

 

Results

Simulation Results

We present the estimation performance of all assessed methods under two scenarios, including mean value, empirical bias, and the MSEs for point estimates, mean SEs, average length of 95% CIs, and ECPs for 95% CI estimates of the detected breakpoints in Tables 2 and 3. Additionally, Supplementary Tables S1 and S2 in Section D of the Supplementary Material demonstrate the estimation performance for the estimates of start slope and slope changes for all assessed methods. 

Table 2.

 

Estimation of Breaking Points at −10 and 2 by Different Statistical Methods Using Simulation Scenario 1 (Synthetic Dataset Without Outliers)

 

Table 3.

 

Estimation of Breaking Points at −10 and 2 by Different Statistical Methods Using Simulation Scenario 2 (Synthetic Dataset With Outliers)

 

Under the first scenario, which involves synthetic datasets without outliers (Table 2, Supplementary Table S1), the proposed method, LIMBARE, and SMM are statistically efficient, as indicated by their negligible bias and small MSEs for point estimates. In comparison, the point estimates of RMESR and CROSS have a larger bias and MSE than the other methods. The CI estimates of LIMBARE and SMM have satisfactory ECPs close to the nominal level of 95%. However, the CI estimates of RMESR and CROSS tend to have much lower ECPs than nominal. The under-coverage of RMESR is mainly due to the underestimation of SEs, whereas the poor performance of the CI estimates of CROSS is largely driven by the bias of its point estimate. 

Under the second scenario, which involves synthetic datasets with outliers (Table 3, Supplementary Table S2), the proposed method, LIMBARE, has the markedly smallest bias and MSE for point estimates, indicating its statistical efficiency and superior robustness to outliers. The point estimates of all competing methods have a large bias and MSE, indicating the failure to accurately detect breakpoints and estimate slopes under this scenario. The ECPs of all competing methods are substantially lower than nominal, whereas LIMBARE maintained ECPs closest to nominal. However, we detected a slight over-coverage of the CI estimate of LIMBARE. 

Broken Stick Analysis of the Longitudinal Ophthalmic Study

The histograms in Supplementary Figure S1 show that RNFL thickness follows a relatively normal distribution, whereas the distribution of CDR displays a cluster of data points deviating from the main distribution, as highlighted by the red rectangle, indicating the existence of outliers in the CDR data. No breakpoints were detected in the association between time from baseline and MD, RNFL thickness, and CDR, respectively, using, LIMBARE, SMM, and RMESR. Figure 2 illustrates the estimated rate of change using LIMBARE adjusted for age, eyes, signal strength, and disc area where no breakpoint was detected. Specifically, MDs decreased by 0.29 dB/y, RNFL thickness decreased by 0.69 µm/y, and CDRs increased by 0.005/y. We then conducted a piecewise association analysis between visual field MD and RNFL thickness, as well as MD and CDR, using four different methods: LIMBARE, SMM, RMESR, and CROSS. Figures 3 and 4 present the results of these analyses. 

Figure 2.

 

No breakpoint was detected for MD (left), RNFL thickness (middle), or CDR (right) against the time from baseline. The fitted lines indicate the estimated rate of change per year.

Figure 3.

 

Estimated piecewise associations between visual field MD and RNFL thickness and from all assessed methods.

For the MD–RNFL association (Fig. 3, Table 4), LIMBARE, SMM, and RMESR consistently identified two breakpoints, whereas CROSS detected only one. Specifically, LIMBARE indicated that RNFL thickness decreased significantly and rapidly at a rate of 1.9 µm/dB (95% CI, 1.52–2.29) before reaching the breakpoint of −4.31 dB (Fig. 2, top right to bottom left, which represents the disease spectrum from mild to severe glaucoma). After that, RNFL thickness still significantly decreased but at a slower rate of 1.05 µm/dB (95% CI, 0.81–1.29) until the next breakpoint of −15.70 dB. The association between MD and RNFL approached a plateau with an insignificant rate of change when MD was smaller than −15.70 dB. In contrast, CROSS captured only one breakpoint at MD of −4.08 dB, which is inconsistent with the other three methods. The estimated slope from CROSS indicated that RNFL thickness decreased significantly with MD when the MD was larger than −4.08 dB, which is not typically observed clinically. Noticeably, RMESR reported the narrowest CI estimates for both breakpoints and slopes, which suggests the possibility of under-coverage of the ground truths of these parameters (Table 4). 

Table 4.

 

Parameter Estimates (95% CIs) for Broken Stick Association Analysis Between MD and RNFL

 

For the MD–CRD association (Fig. 4, Table 5), LIMBARE, SMM, and RMESR identified one breakpoint but CROSS detected none. LIMBARE and RMESR, which were designed to accommodate outliers, reported breakpoints close to each other at −8.61 dB and −9.83 dB, respectively. The corresponding slope estimates from LIMBARE and RMESER indicated that CDR increased at the rate of 0.0131/dB (95% CI, 0.0107–0.0155) (Fig. 3, right to left) and 0.0099/dB (95% CI, 0.0089–0.0109) before the MD decreased to the breakpoint. After that, the increasing trend slowed down with a rate of 0.0023/dB (95% CI, 0.0008–0.0038) and 0.0022/dB (95% CI, 0.0014–0.0030). Yet, RMESR reported much narrower CI estimates for both breakpoints and slopes. In comparison, SMM showed that CDR increased at the rate 0.0159/dB (95% CI, 0.0140–0.0178) when the MD decreased before reaching the breakpoint of −12.46 dB. After the breakpoint, CDR decreased at the rate of 0.0019/dB (95% CI, −0.0011 to 0.0048) as the MD declined. The patterns identified by LIMBARE and RMESR agreed with clinical observations, but SMM showed a trend that is less frequently observed. RMESR resulted in the narrowest CI estimates compared with LIMBARE. 

Figure 4.

 

Estimated piecewise associations between visual field MD and CDR from all assessed methods.

Table 5.

 

Parameter Estimates (95% CIs) for Broken Stick Association Analysis Between MD and CDR

 

Overall, the results suggest that LIMBARE is a robust method for identifying breakpoints in the association between MD and RNFL thickness or MD and CDR by providing both point and CI estimates. Conventional methods such as CROSS may miss these breakpoints, leading to inaccurate or incomplete understanding of the relationships between these variables. The identified breakpoints and slopes may provide valuable clinical insights into the study of glaucoma. 

Conclusions

We have proposed a novel advanced LIMBARE for longitudinal data analysis that provides efficient point estimation of breakpoints and other association parameters, largely mitigating the impact of outliers and improving the estimation of CIs for both breakpoints and slopes. Simulation studies show that LIMBARE outperforms three competing methods in terms of negligible bias, smallest mean squared error, and empirical coverage probabilities closest to the nominal level of 95% under the scenarios with and without outliers. 

In the application of broken stick analysis to a longitudinal clinical study, the piecewise associations between MD and RNFL thickness estimated by LIMBARE, SMM, and RMESR were comparable and consistent with clinical observations. Numerous studies have demonstrated the presence of minimal measurable RNFL thickness (floor effect) where no further thinning of RNFL can be detected with OCT.^19,20 As the RNFL approaches this level, the rate of change in its thickness slows down until it levels off, whereas the visual field continues to demonstrate progression. According to our results, when the MD < −15.7 dB, no further change can be expected in the average RNFL thickness. The piecewise associations between the MD and CDR estimated by LIMBARE, SMM, and RMESR showed similar overall patterns, but the results for CROSS were inconsistent with three other methods. Based on our findings, when a subject attains an MD of −8.6 dB or a CDR of 0.76, the rate of change in OCT-measured CDR markedly decelerates, even as the disease appears to be progressing further based on MD assessments. 

Comparing the goodness-of-fit among models continues to be a challenge, primarily stemming from disparities in data processing. Traditional information criteria such as the Akaike information criterion (AIC) and BIC, along with statistical tests such as Vuong's test,²¹ are primarily structured to assess models operating on the same dataset. However, LIMBARE employs the LTS method, a technique that involves trimming a small portion of observations in the dataset to accommodate outliers. Consequently, this approach renders the calculated AIC, BIC, and Vuong's test results unable to be compared with those obtained from other models. However, the simulation study showed that LIMBARE has the best performance in predicting accuracy and statistical inference. 

In conclusion, LIMBARE is recommended for future longitudinal ophthalmic studies as it has exhibited the most favorable results in both simulation scenarios and illustrated patterns that comply with clinical observation. Furthermore, a user-friendly R package, “LIMBARE,” is publicly available at https://github.com/JiyuanHu/. The conventional cross-sectional method is not capable of capturing breakpoints in longitudinal studies, and its use is not recommended. 

Acknowledgments

Supported by grants from the National Institutes of Health (R01-EY013178, R01-EY030770, R01-EY035174, P30-EY013079, U54-MD000538, R33-AG057382) and by an unrestricted grant from Research to Prevent Blindness. 

Disclosure: T. Lee, None; J.S. Schuman, Zeiss (F); M. de los A. Ramos Cadena, None; Y. Zhang, None; G. Wollstein, None; J. Hu, None