The effect of changes in intraocular pressure on the risk of primary openangle glaucoma in patients with ocular hypertension: an application of latent class analysis
 Feng Gao^{1}Email author,
 J Philip Miller^{1},
 Stefano Miglior^{3},
 Julia A Beiser^{2},
 Valter Torri^{4},
 Michael A Kass^{2} and
 Mae O Gordon^{1, 2}
DOI: 10.1186/1471228812151
© Gao et al.; licensee BioMed Central Ltd. 2012
Received: 21 October 2011
Accepted: 25 September 2012
Published: 4 October 2012
Abstract
Background
Primary openangle glaucoma (POAG) is one of the leading causes of blindness in the United States and worldwide. While lowering intraocular pressure (IOP) has been proven to be effective in delaying or preventing the onset of POAG in many largescale prospective studies, one of the recent hot topics in glaucoma research is the effect of IOP fluctuation (IOP lability) on the risk of developing POAG in treated and untreated subjects.
Method
In this paper, we analyzed data from the Ocular Hypertension Treatment Study (OHTS) and the European Glaucoma Prevention Study (EGPS) for subjects who had at least 2 IOP measurements after randomization prior to POAG diagnosis. We assessed the interrelationships among the baseline covariates, the changes of postrandomization IOP over time, and the risk of developing POAG, using a latent class analysis (LCA) which allows us to identify distinct patterns (latent classes) of IOP trajectories.
Result
The IOP change in OHTS was best described by 6 latent classes differentiated primarily by the mean IOP levels during followup. Subjects with high postrandomization mean IOP level and/or large variability were more likely to develop POAG. Five baseline factors were found to be significantly predictive of the IOP classification in OHTS: treatment assignment, baseline IOP, gender, race, and history of hypertension. In separate analyses of EGPS, LCA identified different patterns of IOP change from those in OHTS, but confirmed that subjects with high mean level and large variability were at high risk to develop POAG.
Conclusion
LCA provides a useful tool to assess the impact of postrandomization IOP level and fluctuation on the risk of developing POAG in patients with ocular hypertension. The incorporation of postrandomization IOP can improve the overall predictive ability of the original model that included only baseline risk factors.
Keywords
Latent class analysis Longitudinal data Timedependent covariate Prediction model Survival data Primary openangle glaucoma Intraocular pressure fluctuationBackground
Ocular hypertension is a leading risk factor for the development of primary openangle glaucoma (POAG) which remains one of the major causes of blindness in the United States and worldwide [1–5]. It is estimated that approximately 4%  7% of the population over the age of 40 years have ocular hypertension without detectable glaucomatous damage using standard clinical tests, and thus as many as 3 to 6 million Americans are at risk for developing glaucoma because of ocular hypertension [6–8]. Intraocular pressure (IOP) is the only known modifiable risk factor for POAG. Lowering the level of IOP has been shown to effectively delay or prevent glaucomatous visual damage in different phases of disease progression by many largescale multicenter clinical trials, including the Ocular Hypertension Treatment Study (OHTS) [9], the Early Manifest Glaucoma Trial [10], and the Advanced Glaucoma Intervention Study [11].
In recent years, one of the hot topics in glaucoma research has been the effect of IOP fluctuation (IOP lability), both within a single day (shortterm fluctuation) and from visit to visit (longterm fluctuation) on POAG [12, 13]. Measures of IOP fluctuation have included a wide range of quantities  peak, trough, variance, and range, etc. [13] However, since subjects with high mean IOP often show large IOP variability over time, it is challenging to disentangle the effect of fluctuation from mean IOP. A recently emerged technique for longitudinal data analysis, latent class analysis (LCA) [14], provides an appealing approach to this question. Rather than dealing with individual measures of fluctuation, LCA identifies distinct patterns of longitudinal profiles based on the combination of summary statistics (i.e., mean level and variability) and hence provides information complementary to the conventional methods. LCA uses the patterns of serial biomarker readings available for subjects, together with baseline covariates and disease outcomes, to divide subjects into a number of mutually exclusive subpopulations (classes). The class membership is unobserved (latent) and determined by the classspecific parameters in a datadriven basis.
In this paper, we used LCA to model the postrandomization IOP in the OHTS. For each class, the change of IOP was characterized by 4 parameters: the initial IOP level (I), the linear (L) and quadratic (Q) trend over time, and the variance of IOP (V). We used data from the European Glaucoma Prevention Study (EGPS) [15], another largescale multicenter randomized clinical trial of patients with ocular hypertension, for external independent validation. We first fit an unconditional (without any covariates) LCA to determine the optimal number of distinct patterns that best described the IOP change for each study. Then a conditional model was constructed by adding baseline covariates as the antecedents (predictors) of IOP change and time to POAG as a consequence (outcome) of IOP change [16]. This analysis enhanced our understanding of the interrelationships among the IOP change, the baseline covariates, and the risk of developing POAG. This also provided evidence towards our ultimate goal to improve the prediction of POAG in patients with ocular hypertension.
Methods
Study cohort
Our study used data from OHTS and EGPS, the two largest randomized trials to test safety and efficacy of topical hypotensive medication in preventing the development of POAG. In OHTS, 1636 subjects were randomized to either observation or treatment with ocular hypotensive medication and followed for a median of 78 months [9]. In EGPS, 1077 subjects were randomized to either placebo or an active treatment (dorzolamide) and followed for a median of 55 months [15]. The two studies shared many key similarities in the study protocol and generated data of high quality. In both studies, for example, the outcome ascertainment was performed by specialized resource centers where readers were masked as to randomization assignment and information about the participant’s clinical status, and the attribution of abnormality due to POAG was performed by a masked Endpoint Committee. Detailed information on the similarity and discrepancy between OHTS and EGPS as described by Gordon et al. [17]. This study was approved by the Institutional Review Boards of Washington University in St. Louis and the University Bicocca of Milan.
Summary statistics of baseline predictors and followup data for OHTS and EGPS, where categorical variables are summarized as counts and proportions, while the continuous variables are summarized in means and standard deviations (SD).
Variables  OHTS (N = 1600)  EGPS (N = 971) 

Baseline predictors  
TRT  795 (49.7%)  487 (50.2%) 
Male  687 (42.9%)  445 (45.8%) 
Black  396 (24.8%)  1 (0.1%) 
AGE (decades)  5.56 (0.96)  5.70 (1.02) 
IOP0 (mmHg)  24.9 (2.69)  23.4 (1.62) 
CCT (μm)  572.6 (38.5)  573.3 (37.5) 
PSD (dB)  1.91 (0.21)  2.00 (0.52) 
VCD  0.39 (0.19)  0.32 (0.14) 
BB  71 (4.4%)  64 (6.6%) 
CHB  190 (11.9%)  66 (6.8%) 
DM  188 (11.8%)  55 (5.7%) 
Heart  99 (6.2%)  109 (11.2%) 
HBP  606 (37.9%)  279 (28.7%) 
Postrandomization IOP  
Mean (mmHg)  21.44 (3.45)  19.73 (2.57) 
SD (mmHg)  2.27 (1.04)  2.22 (1.03) 
Median #visits (minmax)  13 (2–16)  9 (2–10) 
POAG  146 (9.1%)  107 (11.0%) 
Statistical analysis
Unconditional LCA
where $\mathrm{Pr}\left({C}_{i}=g\right)$ represented the size (mixing proportion) of g^{th} latent class in the mixture and $f\left({Y}_{i}{C}_{i}=g;{\theta}_{g}\right)$ was the classspecific distribution of Y_{ i } as detailed below.

The mixing probability $\mathrm{Pr}\left({C}_{i}=g\right)$ was modeled as a multinomial logistic regression, $\mathrm{Pr}\left({C}_{i}=g\right)=\frac{exp\left({\alpha}_{0g}\right)}{{\displaystyle \sum _{h=1}^{G}exp\left({\alpha}_{0g}\right)}}$, where α_{0g} represented the log odds of membership in the g^{th} class relative to a reference class (class 1, say), with the parameter in the reference being 0 for identification.

The specification of $f\left({Y}_{i}{C}_{i}=g;{\theta}_{g}\right)$ was aided by our previous experience on the joint modeling of longitudinal IOP and time to POAG in OHTS [19]. The joint model identified IOP variability as an independent predictor for POAG and also revealed that the IOP change can be better fit by a quadratic functional form. Therefore, we set $f\left({Y}_{i}{C}_{i}=g;{\theta}_{g}\right)={I}_{g}+{L}_{g}{t}_{i}+{Q}_{g}{t}_{i}^{2}+{\epsilon}_{i}$_{,} with ${\epsilon}_{i}{\text{~N(0, V}}_{\text{g}})\phantom{\rule{0.5em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}{\theta}_{g}=\left\{{I}_{g},{L}_{g},{Q}_{g},{V}_{g}\right\}.$ Because high IOP was an eligibility criterion in both OHTS and EGPS, the estimated initial level (intercept I_{g}) may be influenced by “regression to the mean”. To address this concern, we reset the time 0 and the intercept was actually estimated at 1year after randomization. We also assumed that followup IOPs were measured regularly every 6 months according to the protocol, i.e., with timing t_{i} = {−0.5, 0, 0.5, 1, …}. Figure 1A showed the diagram of an unconditional LCA for the OHTS data.

Given the estimated parameters ${\stackrel{\u2322}{\theta}}_{g}$ and the observed IOP, each individual can be assigned to the most likely class based on the probability of class membership (often termed as posterior class probability) [14],${p}_{\mathit{ig}}=\frac{\stackrel{\u2322}{P}r\left({C}_{i}=g\right)\xb7f\left({Y}_{i}{C}_{i}=g;{\stackrel{\u2322}{\theta}}_{g}\right)}{{\displaystyle \sum _{h=1}^{G}\{}\stackrel{\u2322}{P}r\left({C}_{i}=h\right)\xb7f\left({Y}_{i}{C}_{i}=h;{\stackrel{\u2322}{\theta}}_{h}\right)\}}\text{.}$(2)
The best unconditional LCA was selected by enumerating and comparing a set of competing models differing only in the number of classes. In this paper, the model comparison was based primarily on the log likelihood values, including the Bayesian Information Criteria (BIC, with a smaller BIC indicating a better fit) and the LoMendellRubin adjusted likelihood ratio test (LMRLRT) [20]. A significant test of LMRLRT indicated that the model with G1 classes should be rejected in favor of the Gclass LCA. In addition to the above statistical criteria, we also specified a minimum size for each class (with at least 5% participants in OHTS or 10% participants in EGPS) to ensure reliable withinclass estimation. Once an optimal LCA was developed, a bootstrap method was used to assess whether patients with different patterns of IOP change have different susceptibility to POAG. Specifically, a class membership was generated for each individual from a multinomial distribution using the posterior class probability, and then a Cox model was fit to assess the effects of latent classes on POAG. Summary statistics such as hazard ratios and their 95% confidence intervals were estimated by repeating the above procedure 1,000 times.
Conditional LCA

Similar to Model (1), the term $\mathrm{Pr}\left({C}_{i}=g;{\alpha}_{g}\right)=\frac{exp\left({\alpha}_{0g}+{\alpha}_{g}{X}_{i}\right)}{{\displaystyle \sum _{h=1}^{G}exp(}{\alpha}_{0h}+{\alpha}_{h}{X}_{i})}$ represented the size of g^{th} class in the mixture distribution and $f\left({Y}_{i}{C}_{i}=g;{\theta}_{g}\right)$ described the withinclass IOP change.

The term $\lambda \left({T}_{i}{C}_{i}=g;\beta \right)={\lambda}_{0g}(\text{t)}\xb7\text{exp}\left(\beta {X}_{i}\right)$ described the risk of developing POAG in g^{th} class and $S\left({T}_{i}{C}_{i}=g;\beta \right)=exp(\int {\lambda}_{0g}(\text{t)}\xb7\text{exp}(\beta {X}_{i}\left(\text{dt}\right)$ was the corresponding cumulative POAGfree probability, where λ_{0g}(t) was the classspecific baseline hazard with all covariates being 0. In this paper, λ_{0g}(t) was approximated by a piecewise stepfunction with a 6month interval. Following the conventional practice in joint latent class modeling [21, 22], we assumed that the association between IOP change and time to POAG was introduced exclusively via λ_{0g}(t), so that the longitudinal process and survival process were completely independent given the class membership. Therefore, neither timedependent IOP values nor random effects of IOP were included in the survival function. We also assumed that the effects of covariates on POAG were common across latent classes.
where $({\widehat{\theta}}_{g},{\widehat{\alpha}}_{g},\widehat{\beta},\text{and}\phantom{\rule{0.25em}{0ex}}{\widehat{\lambda}}_{0g}(\text{t}))$ were the estimated parameters from the conditional LCA. In this paper, the parameter estimation for LCA was implemented using statistical package Mplus [23], while all the other analyses were performed using statistical package R [24].
Results
Unconditional LCA
Fitting statistics of 7 competing models that are only different in the number of latent classes
# latent classes (G)  OHTS  EGPS  

BIC  LMRLRT*  Minimal class size  BIC  LMRLRT  Minimal class size  
2  97097  <0.001  47%  39235  <0.001  44% 
3  94219  0.002  24%  38395  0.001  14% 
4  92922  0.609  14%  38109  0.005  11% 
5  92107  0.003  13%  37870  0.009  12% 
6  91644  0.042  10%  37760  0.452  9% 
7  91289  0.147  7%  37682  0.060  5% 
8  91045  0.406  6%  37608  0.011  4% 
Observed proportions of POAG, as well as estimated hazard ratios (HR) and 95% confidence intervals (CI) for POAG development in the unconditional LCAs for the OHTS and EGPS data, where the HR and 95% CI were based on 1000 bootstrapping samples to account for the uncertainty in the latent class membership
Latent class  OHTS  EGPS  

POAG%  HR  95% CI  POAG%  HR  95% CI  
1  5.9%  1.00    8.3%  1.00   
2  3.9%  0.59  0.37  0.88  10.2%  1.28  0.76  2.06 
3  4.3%  0.83  0.57  1.14  8.7%  1.13  0.73  1.65 
4  10.1%  1.87  1.32  2.57  10.5%  1.40  0.85  2.18 
5  11.4%  1.93  1.50  2.46  19.4%  2.66  1.92  3.69 
6  31.2%  5.61  4.46  7.08 
Distribution of the randomization groups across latent classes, where the latent classes were based on the most likely posterior class probability from the optimal unconditional LCAs for the OHTS and EGPS data
Latent class  OHTS  EGPS  

Observation  Treatment  Placebo  Treatment  
1  15 (1.9%)  191 (24.0%)  113 (23.3%)  143 (29.4%)  
2  67 (8.3%)  385 (48.4%)  69 (14.3%)  112 (23.0%)  
3  226 (28.1%)  106 (13.3%)  162 (33.5%)  136 (27.9%)  
4  55 (6.8%)  84 (10.6%)  64 (13.2%)  63 (12.9%)  
5  279 (34.7%)  19 (2.4%)  76 (15.7%)  33 (6.8%)  
6  163 (20.2%)  10 (1.3%)  
Total  805 (100%)  795 (100%)  484 (100%)  487 (100%) 
Conditional LCA
Estimated parameters of the conditional LCAs in the OHTS and EGPS data
OHTS  

Variables  Parameters for IOP change and the effects of covariates on class membership  Effects on POAG  
Class 1  Class2 (Ref.)  Class 3  Class 4  Class 5  Class 6  
Class Size IOP Change  14.2%  27.1%  21.1%  9.1%  17.6%  10.9%  
I  17.58(0.20)^{#}  19.83(0.21)^{#}  22.30(0.21)^{#}  22.82(0.79)^{#}  24.74(0.22)#  27.70(0.28)^{#}   
L  −0.57(0.08)^{#}  −0.53(0.06)^{#}  −0.47(0.09)^{#}  −0.95(0.30)^{#}  −0.20(0.09)*  0.16(0.17)   
Q  0.06(0.01)^{#}  0.05(0.01)^{#}  0.05(0.01)^{#}  0.07(0.04)  0.02(0.02)  −0.05(0.03)   
V  4.36(0.27)^{#}  4.30(0.32)^{#}  4.80(0.25)^{#}  16.07(1.46)^{#}  4.66(0.31)#  12.15(1.15)^{#}   
Covariates  
Intercept  −2.64(0.63)^{#}  2.08(0.47)^{#}  −0.06(0.47)  2.35(0.54)^{#}  1.04(0.61)    
TRT  1.60(0.66)^{*}    −3.25(0.35)^{#}  −2.19(0.63)^{#}  −5.98(0.56)^{#}  −6.46(0.70)^{#}  0.16(0.29) 
MALE  0.25(0.23)    −0.99(0.24)^{#}  0.24(0.27)  −0.68(0.28)^{*}  −0.22(0.30)  0.23(0.19) 
RACEB  −0.10(0.27)    −0.37(0.30)  0.75(0.31)^{*}  −0.91(0.34)^{*}  0.05(0.37)  −0.05(0.24) 
AGE  0.06(0.12)    −0.01(0.12)  0.08(0.16)  −0.05(0.14)  0.13(0.15)  0.18(0.09) ^{*} 
IOP0  −0.79(0.18)^{#}    0.21(0.22)  1.03(0.36)^{*}  0.89(0.18)^{#}  1.73(0.22)^{#}  −0.10(0.11) 
CCT  −0.35(0.12)^{*}    0.18(0.11)  −0.08(0.17)  0.14(0.13)  0.09(0.16)  −0.64(0.13)^{#} 
PSD  0.13(0.11)    0.04(0.11)  0.08(0.13)  0.17(0.13)  0.12(0.15)  0.23(0.10)^{*} 
VCD  0.06(0.11)    −0.08(0.12)  −0.15(0.17)  −0.10(0.13)  −0.09(0.15)  0.60(0.10)^{#} 
BB  −0.60(0.48)    −0.37(0.61)   ^{**}  −0.30(0.62)  −1.21(0.77)  0.19(0.57) 
CHB  −0.19(0.37)    −0.32(0.42)  0.47(0.45)  −0.42(0.44)  −0.50(0.49)  0.09(0.31) 
DM  −0.23(0.35)    0.22(0.32)  −0.71(0.45)  0.23(0.36)  0.64(0.40)  −1.67(0.53)^{*} 
HEART  0.70(0.44)    0.10(0.49)  0.43(0.56)  −0.09(0.56)  −0.48(0.73)  0.71(0.29)^{*} 
HBP  0.47(0.24)^{*}    0.40(0.27)  0.06(0.34)  0.41(0.30)  0.66(0.33)^{*}  0.08(0.22) 
EGPS  
Variables  Parameters for IOP change and the effects of covariates on class membership  Effects on POAG  
Class1 (Ref.)  Class 2  Class 3  Class 4  Class 5  
Class Size IOP Change  26.3%  20.1%  29.3%  12.9%  11.4%  
I  18.66(0.25) ^{#}  18.24(0.18)^{#}  21.24(0.20)^{#}  21.85(0.47)^{#}  24.33(0.34)^{#}    
L  −0.34(0.14)^{*}  −1.29(0.18)^{#}  −0.25(0.11)^{*}  −1.76(0.28)^{#}  0.06(0.26)    
Q  0.05(0.03)  0.11(0.04)*  0.02(0.03)  0.02(0.07)  −0.08(0.07)    
V  3.79(0.23)^{#}  3.75(0.25)^{#}  4.59(0.24)^{#}  7.12(0.85)^{#}  12.17(1.32)^{#}    
Covariates  
Intercept    −0.84(0.33)^{*}  0.30(0.28)  −0.85(0.58)  −0.71(0.34)^{*}    
TRT    0.36(0.27)  −0.65(0.23)^{*}  −0.58(0.37)  −1.79(0.37)^{#}  −0.01(0.21)  
MALE    −0.35(0.29)  0.14(0.23)  0.18(0.40)  0.38(0.33)  −0.24(0.22)  
Black              
AGE    −0.09(0.16)  0.01(0.13)  0.40(0.23)  0.45(0.20)^{*}  0.16(0.10)  
IOP0    −0.61(0.24)^{*}  0.82(0.17)^{#}  1.24(0.24)^{#}  1.79(0.23)^{#}  0.11(0.13)  
CCT    −0.33(0.13)^{*}  −0.14(0.12)  −0.43(0.15)^{*}  0.09(0.15)  −0.36(0.12)^{*}  
PSD    0.27(0.16)  −0.23(0.14)  −0.18(0.24)  −0.41(0.23)  0.18(0.09)^{*}  
VCD    −0.13(0.17)  −0.03(0.13)  0.72(0.26)^{*}  0.17(0.16)  0.46(0.12)^{#}  
BB    −0.17(0.50)  −0.82(0.52)   ^{**}  −0.58(0.69)  −0.07(0.41)  
CHB    −0.10(0.56)  −0.97(0.52)  0.89(1.03)  −1.22(0.79)  −0.28(0.47)  
DM    −0.46(0.52)  0.12(0.45)  −1.32(1.27)  0.82(0.81)  −0.18(0.54)  
HEART    0.78(0.41)  0.11(0.44)  −0.83(0.77)  −0.79(0.61)  0.74(0.32)^{*}  
HBP    0.02(0.36)  0.53(0.30)  −1.27(0.90)  0.11(0.54)  0.24(0.26) 
Effects of the baseline covariates on IOP classification

To identify baseline predictors for IOP classification, we only focused on factors that were significantly associated with the high risk groups (Classes 4, 5, 6 in OHTS, and Class 5 in EGPS), while treating the lowest risk group (Class 2 in OHTS and Class 1 in EGPS) as the reference. In OHTS, treatment assignment and baseline IOP were two most important predictors for IOP classification. Subjects randomized to treatment group had a much lower chance of inclusion in the high risk groups (with OR = 0.11, 0.003, and 0.002 for Classes 4, 5, and 6, respectively), while these with a higher baseline IOP were more likely to be in the Classes 4, 5, or 6 (with OR = 2.80, 2.44, and 5.64 respectively). The results also showed that male subjects were less likely to be in Class 5 (OR = 0.51), the black subjects were more likely to be in Class 4 (OR = 2.12) but with a lower chance in Class 5 (OR = 0.40), and subjects with a history of hypertension were more likely in Class 6 (OR = 1.93). In EGPS, the results confirmed that treatment assignment (OR = 0.17) and baseline IOP level (OR = 5.99) were important predictors for Class 5. The result also showed that older age (OR = 1.57) was significantly associated with Class 5.
Effects of the baseline covariates on the risk of POAG development

As expected, the effects of baseline covariates on the risk of developing POAG from the conditional LCA reached consistent conclusions as previous analyses using Cox models [17, 25]. In OHTS, subjects with older age (HR = 1.20), higher PSD (HR = 1.26), large VCD (HR = 1.82), and history of heart diseases (HR = 2.03) had a higher risk of developing POAG, while thicker CCT (HR = 0.53) and history of diabetes (HR = 0.19) reduced the risk of developing POAG. Interestingly, despite marked differences between OHTS and EGPS in the patterns of IOP change, the EGPS confirmed 4 of the 6 predictors (except age and history of diabetes) identified in OHTS. In both studies, baseline IOP and treatment assignment were not significantly associated with POAG directly, but appeared to affect the risk indirectly through their strong influence on the classification of IOP change.
Sensitivity analysis comparing the overall predictive accuracy (measured as Cindex and Calibration Chisquare statistics) for LCAs with different model specifications
Model  Model Features  C index  Calibration Chisquare  

OHTS  EGPS  OHTS  EGPS  
0  Cox mode with baseline factors only  0.778  0.706  5.0  2.1 
1  LCA with a quadratic withinclass functional form  0.821  0.719  11.3  7.0 
2  LCA with a linear withinclass functional form  0.825  0.720  10.2  4.9 
3  LCA with a constant withinclass functional form  0.823  0.727  10.5  13.5 
Sensitivity analyses
As in all the statistical models, LCAs were inevitably based on certain assumptions. One assumption of our LCA was that the trajectories of IOP followed a quadratic functional form. It is known that the parameter estimates, class sizes, and interpretation of latent classes could be heavily influenced by the withinclass distribution of longitudinal data [16]. In this section, first we assessed the sensitivity of risk prediction to different LCA specifications. Table 6 presented the Cindex and calibration statistics for LCAs after removing the quadratic term (Model 2) or removing both quadratic and linear terms (Model 3). The results showed that LCAs had a robust performance in terms of predictive accuracy for POAG development.
Finally, our LCA also made an implicit assumption that the baseline covariates influenced the IOP change exclusively through their effects on the class membership (i.e., no direct effects on the withinclass growth parameters). The validity of this assumption can be checked by comparing the conditional LCAs with the unconditional models. The assumption violation is often signified by a dramatic shifting in the meaning and size of latent classes when the baseline predictors are added to the unconditional LCA [16]. Based on the estimated classspecific parameters (Table 5, Figures 2 and 3), this assumption was well satisfied in both studies.
Discussion
In recent years, one of the hot topics in glaucoma research has been the effect of IOP fluctuation on POAG. Although more and more studies have confirmed that a decrease in the mean IOP level can reduce the risk of developing POAG, the findings from major prospective clinical trials about the impact of IOP fluctuation on POAG remain controversial [25, 27–30]. In this paper, we analyzed the postrandomization IOPs from OHTS and EGPS taking a latent class analysis (LCA) approach. The LCA allows us to identify distinct patterns of IOP change over time and then associates the changes in IOP with the risk of POAG. The results from both studies showed that different patterns of IOP change could markedly affect the risk of POAG (irrespective of their baseline, prerandomization IOP levels). In OHTS, the change in IOP was best described by 6 distinct patterns. The model identified a subset of participants in whom IOP variability also played an important role in predicting POAG. This subgroup showed the highest IOP variability and had a higher risk than those with a comparable IOP mean. Comparing to the reference class, these participants were less likely from treatment group (OR = 0.11), more likely selfclassified as being black (OR = 2.12), and had relatively higher baseline IOP (OR = 2.80). However, the subgroup only accounted for about 10% of the OHTS sample, and this may partially explain our finding that IOP variability was an independent risk factor in the OHTS but had little impact on the overall predictive accuracy for POAG (manuscript in progression). In a sensitivity analysis using the nonBlack participants, the LCA identified similar patterns of IOP change as in the whole OHTS dataset. This result was consistent with a treebased model in the OHTSEGPS metaanalysis which showed that race was no longer an important predictor for POAG development after considering other risk factors [17]. In EGPS, LCA identified 5 distinct latent classes and confirmed that those subjects with the highest mean IOP were most likely to develop POAG. However, it failed to disentangle the effect of fluctuation from mean because these participants with the highest mean level also had the largest IOP variability. Interestingly, despite the marked differences between EGPS and OHTS in the treatment intervention and magnitude of IOP lowering achieved, both studies showed that adding IOP change into the baseline model improved the overall predictive accuracy for POAG development.
Conventionally the change of longitudinal data is described using linear mixed models with random coefficients [31]. Though the mixed model recognizes the heterogeneous nature of the data by allowing each individual to have his/her own intercept and slope, it assumes that all individuals come from a single population and uses an average trajectory for the entire population. A LCA analyzes data from a rather different perspective. The model approximates the unknown heterogeneity in the distribution of longitudinal outcome using a finite number of polynomial functions each describing a unique subpopulation [14, 32]. It classifies individuals into distinct groups based on the patterns of longitudinal outcome, so that individual within a group are more similar than those between different groups. This LCA possesses some unique advantages as comparing to conventional methods. First, the model lends itself directly to a set of well characterized subpopulations and also provides a formal statistical procedure to determine the appropriate number of subpopulations. It thus enables the discovery of unexpected yet potentially meaningful subpopulations that may be otherwise missed with conventional methods. Second, the method permits one to relate the developmental patterns of longitudinal data to its antecedents (predictors or covariates) and consequences (clinical outcomes), and thus allows estimation of both direct and indirect (via longitudinal data) effects of a covariate on the distant outcome [16, 23]. Finally, the recent advances of the dual trajectory modeling also allow investigators to assess the joint evolution of multiple longitudinal processes, which may evolve contemporaneously or over different time periods [32].
LCA also provides an attractive alternative for making prediction with timedependent covariates [21, 22]. A LCA takes a joint modeling approach to assess the association between longitudinal and survival data and thus uses information more efficiently, resulting less biased estimates. Unlike the conventional joint models that assess the association via shared random effects [19, 33, 34], a LCA relates the longitudinal process to survival process by latent classes and assumes the two stochastic processes independent given the class membership [22]. Therefore, neither timedependent covariates nor random effects of the longitudinal data are needed in the survival submodel. Such a model specification will avoid the intensive computation to obtain the random effects for new subjects and hence facilitates a realtime individualized prediction [21]. The key to build an accurate prediction in a LCA setting is to have a reliable classification given the observed data. Generally speaking, the more the available serial biomarker readings, the more reliable a classification is. To this consideration, the impact of followup IOP on POAG may be overestimated in OHTS because an average length of 6.5year IOP readings was used to calculate the 5year POAGfree rate. To solve this dilemma, which is rather common in all predictions involving timedependent covariates, one of the most frequently used approaches in medical literature is a landmark analysis that consists of fitting a serial of survival models only to the subjects still at risk, that is, computation of the predictive distribution at a certain time given the history of event and covariates until that moment [35]. In a LCA setting, such a dynamic prediction can be conveniently implemented because the conditional survival probability at any time can be calculated analytically from a single LCA once the parameters are estimated [21].
Despite its advantages, the LCA has several limitations. First, the computational load of LCA can be high, especially for models with complexity structures. In OHTS data (N = 1600), for example, it ran less than 10 minutes for an unconditional 6class LCA, but it took more than 30 minutes to develop the full conditional model. Because of the exploratory nature of data analysis with LCA, the cumulative time can be substantial. For this consideration, in practice the best LCA model is often constructed taking a twostep approach as in this paper. Another issue in LCA is that the loglikelihood function may end up at local rather than global maxima. Fortunately this issue has been taken into consideration by the statistical package Mplus which automatically uses 10 sets of randomly generated starting values for estimation. The program also allows investigators to rerun and compare the estimates from user specified starting values if necessary [23].
Conclusion
LCA provides a useful alternative to understand the interrelationship among the baseline covariates, the change in followup IOP, and the risk of developing POAG. The inclusion of postrandomization IOP can improve the predictive ability of the original prediction model that only included baseline risk factors.
Abbreviations
 IOP:

Intraocular pressure
 POAG:

Primary openangle glaucoma
 LCA:

Latent class analysis
 OHTS:

The ocular hypertension treatment study
 EGPS:

The european glaucoma prevention study
 HR:

Hazards ratio
 OR:

Odds ratio
 CCT:

Central corneal thickness
 PSD:

Pattern standard deviation
 VCD:

Vertical cup/disc ratio.
Declarations
Acknowledgments
This study is supported by grants from the National Eye Institute of Health (EY091369 and EY09341) and Research to Prevent Blindness (RPB).
Authors’ Affiliations
References
 Quigley HA: Number of people with glaucoma worldwide. Br J Ophthalmol. 1996, 80 (5): 389393. 10.1136/bjo.80.5.389.PubMed CentralView ArticlePubMed
 Sommer A, Tielsch JM, Katz J, et al: Racial differences in the causespecific prevalence of blindness in east Baltimore. N Engl J Med. 1991, 325 (20): 14121417. 10.1056/NEJM199111143252004.View ArticlePubMed
 US Department of Health, Education, and Welfare: Statistics on blindness in the model reporting area 1969–1970. 1973, Washington, DC: US Government Printing Office
 Quigley HA, Vitale S: Models of openangle glaucoma prevalence and incidence in the United States. Invest Ophthalmol Vis Sci. 1997, 38 (1): 8391.PubMed
 Hyman L, Wu SY, Connell AM, et al: Prevalence and causes of visual impairment in the Barbados eye study. Ophthalmology. 2001, 108 (10): 17511756. 10.1016/S01616420(01)005905.View ArticlePubMed
 Leibowitz HM, Krueger DE, Maunder LR, et al: The Framingham Eye study monograph: an ophthalmological and epidemiological study of cataract, glaucoma, diabetic retinopathy, macular degeneration, and visual acuity in a general population of 2631 adults, 1973–1975. Surv Ophthalmol. 1980, 24 (suppl): 335610.PubMed
 Armaly MF, Kreuger DE, Maunder L, et al: Biostatistical analysis of the collaborative glaucoma study. I. Summary report of the risk factors for glaucomatous visualfield defects. Arch Ophthalmol. 1980, 98: 21632171. 10.1001/archopht.1980.01020041015002.View ArticlePubMed
 Quigley HA, Enger C, Katz J, et al: Risk factors for the development of glaucomatous visual field loss in ocular hypertension. Arch Ophthalmol. 1994, 112: 644649. 10.1001/archopht.1994.01090170088028.View ArticlePubMed
 Kass MA, Heuer DK, Higginbotham EJ, et al: The ocular hypertension treatment study: a randomized trial determines that topical ocular hypotensive medication delays or prevents the onset of primary openangle glaucoma. Arch Ophthalmol. 2002, 120 (6): 701713. discussion 829–830View ArticlePubMed
 Heijl A, Leske MC, Bengtsson B, et al: Reduction of intraocular pressure and glaucoma progression: results from the early manifest glaucoma trial. Arch Ophthalmol. 2002, 120 (10): 12681279.View ArticlePubMed
 The Advanced Glaucoma Intervention Study Group: The advanced glaucoma intervention study (AGIS): 7. The relationship between control of intraocular pressure and visual field deterioration. Am J Ophthalmol. 2000, 130: 429440.View Article
 Singh K, Shrivastava A: Intraocular pressure fluctuations: how much do they matter?. Curr Opin Ophthalmol. 2009, 20: 8487. 10.1097/ICU.0b013e328324e6c4.View ArticlePubMed
 Sultan MB, Mansberger SL, Lee PP: Understanding the importance of IOP variables in glaucoma: a systematic review. Surv Ophthalmol. 2009, 54: 643662. 10.1016/j.survophthal.2009.05.001.View ArticlePubMed
 Nagin DS: Groupbased modeling of development. 2005, Cambridge: MA. Harvard University PressView Article
 Miglior S, Pfeiffer N, Torri V, Zeyen T, CunhaVaz J, Adamsons I: Predictive factors for openangle glaucoma among subjects with ocular hypertension in the European glaucoma prevention study. Ophthalmology. 2007, 114 (1): 39.View ArticlePubMed
 Petras H, Masyn K: General growth mixture analysis with antecedents and consequences of change. 2009, Springer: Handbook of Quantitative Criminology
 Gordon MO, Torri V, Miglior S, Beiser JA, Floriani I, Miller JP, Gao F, Adamsons I, Poli D, D’Agostino RB, Kass MA: The ocular hypertension treatment study group and the European glaucoma prevention study group. A validated prediction model for the development of primary open angle glaucoma in individuals with ocular hypertension. Ophthalmology. 2007, 114 (1): 1019.View ArticlePubMed
 Gordon MO, Beiser JA, Brandt JD, Heuer D, Higginbotham E, Johnson C, et al: The ocular hypertension treatment study: baseline factors that predict the onset of primary openangle glaucoma. Arch Ophthalmol. 2002, 120: 714720.View ArticlePubMed
 Gao F, Miller JP, Miglior S, Beiser JA, Torri V, Kass MA, Gordon MO: A joint model for prognostic effect of biomarker variability on outcomes: longterm intraocular pressure (IOP) fluctuation on the risk of developing primary openangle glaucoma (POAG). JP Journal of Biostatistics. 2011, 5 (2): 7396.PubMed CentralPubMed
 Lo Y, Mendell NR, Rubin DB: Testing the number of components in a normal mixture. Biometrika. 2001, 88: 767778. 10.1093/biomet/88.3.767.View Article
 ProustLima C, Taylor JMG: Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of posttreatment PSA: a joint modeling approach. Biostatistics. 2009, 10: 535549. 10.1093/biostatistics/kxp009.PubMed CentralView ArticlePubMed
 Lin H, Turnbull BW, Mcculloch CE, Slate EH: Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostatespecific antigen readings and prostate cancer. J Am Stat Assoc. 2002, 97: 5365. 10.1198/016214502753479220.View Article
 Muthén LK, Muthén BO: Mplus User’s guide. 1998–2010, Los Angeles, CA: Muthén & Muthén, Sixth
 R Core Team: R: a language and environment for statistical computing. 2012, Vienna, Austria: R Foundation for Statistical Computing, http://www.Rproject.org/ ISBN 3900051070
 Miglior S, Torri V, Zeyen T, Pfeiffer N, CunhaVaz J, Adamsons I: Intercurrent factors associated with the development of openangle glaucoma in the European glaucoma prevention study. Am J Ophthalmol. 2007, 144: 266275. 10.1016/j.ajo.2007.04.040.View ArticlePubMed
 Harrell FE, Lee KL, Mark DB: Tutorial in biostatistics multivariate prognostic models: Issues in developing models, evaluation assumptions and adequacy, and measuring and reducing errors. Stat Med. 1996, 15: 361387. 10.1002/(SICI)10970258(19960229)15:4<361::AIDSIM168>3.0.CO;24.View ArticlePubMed
 NouriMahdavi K, Hoffman D, Coleman AL, et al: Advanced glaucoma intervention study. Predictive factors for glaucomatous visual field progression in the advanced glaucoma intervention study. Ophthalmology. 2004, 111 (9): 16271635. 10.1016/j.ophtha.2004.02.017.View ArticlePubMed
 Bengtsson BL, Leske MC, Hyman L, Heijl A, Early Manifest Glaucoma Trial Group: Fluctuation of intraocular pressure and glaucoma progression in the early manifest glaucoma trial. Ophthalmology. 2007, 114: 205209. 10.1016/j.ophtha.2006.07.060.View ArticlePubMed
 Caprioli J, Coleman AL: Intraocular pressure fluctuation: a risk factor for visual field progression at low intraocular pressures in the advanced glaucoma intervention study. Ophthalmology. 2008, 115: 11231129. 10.1016/j.ophtha.2007.10.031.View ArticlePubMed
 Medeiros FA, Weinreb RN, Zangwill LM, et al: Longterm intraocular pressure fluctuations and risk of conversion from ocular hypertension to glaucoma. Ophthalmology. 2008, 115: 934940. 10.1016/j.ophtha.2007.08.012.PubMed CentralView ArticlePubMed
 Laird NM, Ware JH: Randomeffects models for longitudinal data. Biometrics. 1982, 38: 963974. 10.2307/2529876.View ArticlePubMed
 Jones BL, Nagin DS: Advances in groupbased trajectory modeling and a SAS procedure for estimating them. Sociological Methods & Research. 2007, 35 (3): 542571.View Article
 Henderson R, Diggle P, Dobson A: Joint modeling of longitudinal measurements and event time data. Biostatistics. 2000, 4: 465480.View Article
 Vonesh EF, Greene T, Schluchter MD: Shared parameter models for the joint analysis of longitudinal data and event times. Statistics in Medicine. 2006, 25: 143163. 10.1002/sim.2249.View ArticlePubMed
 Van Houwelingen HC: Dynamic prediction by landmarking in event history analysis. Scand J Statist. 2007, 34: 7085. 10.1111/j.14679469.2006.00529.x.View Article
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14712288/12/151/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.