Assessing Trends in Functional and Structural Characteristics: A Survey of Statistical Methods With an Example From Ophthalmology

Abstract

Purpose: Clinical decisions on treatment are usually based on short-term records of consecutive measurements of structure and function. Useful models for analyzing average trends and a description of statistical methods for classifying individual subjects on the basis of subject-specific trend progressions are presented.

Methods: Random effects trend models allow intercepts and slopes of the trend regression to vary across subjects around group-specific mean intercepts and mean slopes. Model results assess whether average intercepts and slopes and subject variability in intercepts and slopes are the same across groups. Fisher's discriminant functions are used for classification.

Results: Methods are presented and illustrated on structural visual data from a multiyear perimetry study. Average thickness of the ganglion cell layer from the optical coherence tomography macula scan and of the retinal nerve fiber layer from the optic disc scan for both glaucoma patients on optimal treatment and normal subjects are analyzed. The random effects trend model shows that average intercepts of glaucoma patients and normal subjects are quite different, but that average slopes are the same, and that the subject variability in both intercepts and slopes is larger for the glaucoma group. These findings explain why the subject-specific trend progression is not a good classifier; it is the level of the measurement (intercept or baseline value) that carries useful information in this particular cohort example.

Translational Relevance: Clinicians base decisions on short-term records of consecutive measurements and need simple statistical tools to analyze the information. This paper discusses useful methods for analyzing short time series data. Model results assess whether there exist significant trends and whether average trends are different across groups. The paper discusses whether clinical measures classify patients reliably into disease groups, given their variability. With ever more available data, classification plays a central role of personalized medicine.

Introduction

Clinicians study the progression (i.e., the changes) of functional and structural characteristics over time and need to know how to analyze the resulting data with the appropriate statistical methods. This paper has two objectives: (1) Discuss useful statistical approaches for analyzing groups of short time series and for classifying individual subjects on the basis of subject-specific trend progressions. (2) Illustrate the methodology by analyzing structural optical coherence tomography (OCT) data collected during a multiyear perimetry study conducted at the University of Iowa. 

Methods

The Linear Trend Model for Observations on an Individual Subject

For each patient, n consecutive measurements Display Formula\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\bf{\alpha}}\)\(\def\bupbeta{\bf{\beta}}\)\(\def\bupgamma{\bf{\gamma}}\)\(\def\bupdelta{\bf{\delta}}\)\(\def\bupvarepsilon{\bf{\varepsilon}}\)\(\def\bupzeta{\bf{\zeta}}\)\(\def\bupeta{\bf{\eta}}\)\(\def\buptheta{\bf{\theta}}\)\(\def\bupiota{\bf{\iota}}\)\(\def\bupkappa{\bf{\kappa}}\)\(\def\buplambda{\bf{\lambda}}\)\(\def\bupmu{\bf{\mu}}\)\(\def\bupnu{\bf{\nu}}\)\(\def\bupxi{\bf{\xi}}\)\(\def\bupomicron{\bf{\micron}}\)\(\def\buppi{\bf{\pi}}\)\(\def\buprho{\bf{\rho}}\)\(\def\bupsigma{\bf{\sigma}}\)\(\def\buptau{\bf{\tau}}\)\(\def\bupupsilon{\bf{\upsilon}}\)\(\def\bupphi{\bf{\phi}}\)\(\def\bupchi{\bf{\chi}}\)\(\def\buppsy{\bf{\psy}}\)\(\def\bupomega{\bf{\omega}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\({Y_t}\) are obtained at times Display Formula\(Tim{e_1} \lt Tim{e_2} \lt ... \lt Tim{e_n}\). Times are typically expressed in years, and there is no requirement that times are spaced equally. Times are often expressed relative to the time of the initial visit. For example, Display Formula\(Tim{e_1} = 0\), Display Formula\(Tim{e_2} = 0.5\), Display Formula\(Tim{e_3} = 1.0\), Display Formula\(Tim{e_4} = 1.25\), …, Display Formula\(Tim{e_8} = 4.14\) expresses that eight observations are obtained—at the start of the study, after 6 months, after 12 months, after 15 months, … , and after 4 years and 365(0.14) = 51 days since the start of the study. 

The following trend regression model is considered for the progression of the measurements  

\begin{equation}\tag{1}{Y_t} = \alpha + \beta Tim{e_t} + {\varepsilon _t}{\rm .}\end{equation}

 

The parameter Display Formula\(\alpha \) represents the expected baseline measurement at the time a subject enters the study. The parameter Display Formula\(\beta \) reflects the expected change in the measurement for a unit change in time (which, in our illustration, is 1 year). The model in equation (1) assumes linearity; the mean change over the first year is the same as the mean change in the second year, and so on. The measurement errors Display Formula\({\varepsilon _t}\) are independent, homoscedastic (that is, with constant variance; if the variation increases with the level, a log transformation of the response will stabilize the variance), and normally distributed. Under these assumptions one can obtain the best (least squares) estimates of the parameters, their SEs, confidence intervals, and tests of hypotheses. Least squares estimates of the slope and the baseline are given by Display Formula\(b = {{\sum\nolimits_{t = 1}^n {(Tim{e_t} - \overline {Time} ){Y_t}} } \over {\sum\nolimits_{t = 1}^n {{{(Tim{e_t} - \overline {Time} )}^2}} }}\) and Display Formula\(a = \bar Y - b\overline {Time} \) where Display Formula\(\overline Y = (1/n)\sum\nolimits_{t = 1}^n {{Y_t}} \) and Display Formula\(\overline {Time} = (1/n)\sum\nolimits_{t = 1}^n {Tim{e_t}} \). The estimated SEs are given by Display Formula\(se(b) = {s \over {\sqrt {\sum\nolimits_{t = 1}^n {{{(Time - \overline {Time} )}^2}} } }}\) and Display Formula\(se(a) = s\sqrt {{1 \over n} + {{{{(\overline {Time} )}^2}} \over {\sum\nolimits_{t = 1}^n {{{(Tim{e_t} - \overline {Time} )}^2}} }}} \), where Display Formula\(s = \sqrt {{{\sum\nolimits_{t = 1}^n {{{[{Y_t} - (a + bTim{e_t})]}^2}} } \over {n - 2}}} \) is the unbiased estimate of the measurement variability. 

These results can be found in most regression texts; see, for example, Sections 2.2 to 2.5 of Abraham and Ledolter.¹ 

The Random Effects Trend Model: Modeling Linear Trends for Several Subjects From More Than One Group

Typically, there is not just one single subject, but there are several subjects, which are observed at possibly different times and with different sample sizes. And, typically there is not just a single group of subjects, but several groups. For example, a group of normal/healthy subjects and a group (or groups) of subjects with a certain disease (such as glaucoma in the example in the Results section, or a blast injury when studying traumatic brain injury). 

We denote the observation on the ith subject in group j (j = 1 for normal, j = 2 for disease) at time t with Display Formula\(Y_{it}^{(j)}\)and the time at which the observation is obtained as Display Formula\(Time_{it}^{(j)}\), with Display Formula\(Time_{i1}^{(j)} = 0\) when time is measured from the first visit at the clinic. The model in equation (2) allows for subject variability in the baselines and the slopes. The model  

\begin{equation}\tag{2}Y_{it}^{(j)} = ({\alpha ^{(j)}} + \tilde \alpha _i^{(j)}) + ({\beta ^{(j)}} + \tilde \beta _i^{(j)})Time_{it}^{(j)} + {\rm{\ }}\varepsilon _{it}^{(j)}{\rm{\ }}{\rm }\end{equation}

includes group-specific average baselines and slopes (the fixed effects, Display Formula\({\alpha ^{(j)}}\)and Display Formula\({\beta ^{(j)}}\)) and subject-specific random effects for the baselines and slopes. The random effects Display Formula\(\tilde \alpha _i^{(j)}\) and Display Formula\(\tilde \beta _i^{(j)}\) are independent across subjects, with zero means and group-specific standard deviations Display Formula\(\sigma _\alpha ^{(j)}\) and Display Formula\(\sigma _\beta ^{(j)}\). The model allows for correlation among baseline and slope, Display Formula\(\rho _{\tilde \alpha \tilde \beta }^{(j)} = Corr(\tilde \alpha _i^{(j)},\tilde \beta _i^{(j)})\), in case there is reason to believe that a subject's slope changes with the baseline. The measurement errors Display Formula\(\varepsilon _{it}^{(j)}\) are independent across time and subjects, with group-specific standard deviations Display Formula\(\sigma _\varepsilon ^{(j)}{\rm {.}}\)  

For two groups (in our example, the groups of normal and glaucoma patients), the model in equation (2) can equivalently be written as  

\begin{equation}\tag{3}Y_{it}^{(j)} = [\alpha + {\alpha ^*}IND(Glaucoma)_{it}^{(j)} + \tilde \alpha _i^{(j)}] + [\beta + {\beta ^*}IND(Glaucoma)_{it}^{(j)} + \tilde \beta _i^{(j)}]Time_{it}^{(j)} + \varepsilon _{it}^{(j)}{\rm }\end{equation}

where Display Formula\(IND(Glaucoma)_{it}^{(j)}\) is an indicator variable with value 0 when subject i is from the normal group (Display Formula\(j = 1\)) and value 1 when subject i is from the glaucoma group (Display Formula\(j = 2\)). The parameter Display Formula\({\alpha ^*}\) expresses the difference in the average baselines of glaucoma and normal subjects. The parameter Display Formula\({\beta ^*}\) expresses the difference in the average slopes of glaucoma and normal subjects.  

The models in equations (2) and (3) are known as repeated measurements or mixed linear regression models. They are studied in detail by Searle² and Searle et al.³ Estimation of such models can be carried out with the PROC MIXED procedure of the SAS Statistical Software,⁴ or with the R Statistical Software⁵ through its libraries nlme and lme4. Parameter estimates are usually obtained through restricted maximum likelihood (REML). 

The parameter Display Formula\({\alpha ^*}\) expresses the difference between the average baselines of the two groups while the parameter Display Formula\({\beta ^*}\) assesses whether the average slopes of the two groups are the same. The model is very general, and not all parameters may be needed. Slopes and baselines may be uncorrelated if subject trend changes are unrelated to baselines; that is, Display Formula\(\rho _{\tilde \alpha \tilde \beta }^{(1)} = \rho _{\tilde \alpha \tilde \beta }^{(2)} = 0\). Subject-variability among baseline intercepts, and also among slopes, may be the same for the two groups; that is, Display Formula\(\sigma _\alpha ^{(1)} = \sigma _\alpha ^{(2)}\) and Display Formula\(\sigma _\beta ^{(1)} = \sigma _\beta ^{(2)}\). And, the measurement variability for the two groups may be the same; that is, Display Formula\(\sigma _\varepsilon ^{(1)} = \sigma _\varepsilon ^{(2)}\). 

Time series data for Display Formula\({n_1}\) subjects from the normal group and Display Formula\({n_2}\) subjects from the disease/experimental group are available to estimate the model. The model with the general covariance structure and models with restricted covariance structure and equality constraints among the variances can be compared and assessed using either the traditional log-likelihood ratio test or the Bayesian (BIC) or Akaike (AIC) information criteria, which penalize the likelihoods for the number of parameters; see Chapter 3 of Verbeke and Molenberghs.⁶ 

Random effects models allow us to address important questions. A comparison of the average trend of normal subjects with the average trend of subjects from a comparison group (subjects with a certain disease) is often the main objective of a study, and this objective is addressed by testing the hypothesis whether or not Display Formula\({\beta ^*} = 0\). Another question of interest is whether the variability of slopes in the normal group is the same as the variability of slopes in the disease group, and this question can be answered by testing whether Display Formula\(\sigma _\beta ^{(1)} = \sigma _\beta ^{(2)}\). Hypotheses about baselines (whether average baselines are the same across groups, and whether subject variability of baselines are the same across groups) can be assessed similarly. 

Classifying Subjects Into One of Two (or More) Groups

Another important objective of medical studies is classification. That is, how to use characteristics (such as the slope in a trend regression) to classify a new subject into either one (the healthy) or the other (the disease) group. There may be strong evidence that the difference between the average slopes of the two groups is statistically significant, especially if the sample sizes in the two groups are large. But with large subject variability among the slopes and with considerable overlap of the slopes from the two groups, it will be difficult to classify subjects; the misclassification errors (classifying a true healthy subject as diseased or a diseased subject as healthy) will most likely be unacceptable. Establishing statistical significance of a difference of average slopes is, in general, much easier than the classification of subjects into one or the other group. 

In order to carry out the classification, we adopt a two-step approach that (1) estimates the baseline measurement and slope for each subject and (2) compares and uses for classification the estimated baseline measurements and slopes of all subjects from the two groups. Least squares estimates of the baseline and of the slope of the trend regression are calculated for each subject. This leads to estimates Display Formula\((a_1^{(1)},b_1^{(1)}),(a_2^{(1)},b_2^{(1)}),...,(a_{{n_1}}^{(1)},b_{{n_1}}^{(1)})\) for subjects in group 1 and estimates Display Formula\((a_1^{(2)},b_1^{(2)}),(a_2^{(2)},b_2^{(2)}),...,(a_{{n_2}}^{(2)},b_{{n_2}}^{(2)})\) for subjects in group 2. The reliability of the least squares estimates depends on the number and the time-spacing of the observations. One may want to omit estimates from subjects with only few measurements to make the estimates across subjects comparable. 

In general, each subject is characterized by k features, and information about group membership (healthy group and diseased group) is available on each subject (that is, one deals with supervised learning). If the slope of the progression is the only feature to be considered, then k = 1. If the baseline and the slope of the trend model are used to describe the state of a subject, then k = 2. The baseline may be a useful feature as it measures the subject's starting level of the measured characteristic. The Fisher⁷ linear discriminant analysis classifies a patient with feature vector x into group 2 (disease) if  

\begin{equation}\tag{4}{({\mu _1} - {\mu _2})^T}{\Sigma ^{ - 1}}x \lt (1/2){({\mu _1} - {\mu _2})^T}{\Sigma ^{ - 1}}({\mu _1} + {\mu _2}){\rm .}\end{equation}

 

The superscript Display Formula\(^T\) stands for the vector (matrix) transpose. Here Display Formula\({\mu _1}\) and Display Formula\({\mu _2}\) are the mean vectors (of dimension k) for group 1 (healthy) and group 2 (diseased), and they are estimated with the sample mean vectors of features for subjects in groups 1 and 2. The common covariance matrix of the features Display Formula\(\Sigma \) (a matrix of dimension [k × k]) is estimated with the pooled covariance matrix Display Formula\({{({n_1} - 1){S_1} + ({n_2} - 1){S_2}} \over {({n_1} + {n_2} - 2)}}\), where Display Formula\({S_1}\) and Display Formula\({S_2}\) are the covariance matrices for group 1 and group 2, respectively. The linear discriminant function minimizes the expected cost of misclassification if the features have normal distributions with group mean vectors Display Formula\({\mu _1}\) and Display Formula\({\mu _2}\) and identical group covariance matrices, equal prior probabilities, and equal misclassification costs. See, for example, Chapter 11 of Johnson and Wichern⁸ and Chapter 12 of Ledolter.⁹ If k = 1, Fisher's linear discriminant analysis classifies a subject with feature x into the group with the closest mean. 

The optimal discriminant function becomes quadratic if the covariance matrices for groups 1 and 2 are not the same. The Fisher⁷ quadratic discriminant function classifies a subject with feature vector x into group 2 (disease) if  

\begin{equation}\tag{5} - (1/2){x^T}(\Sigma _1^{ - 1} - \Sigma _2^{ - 1})x + (\mu _1^T\Sigma _1^{ - 1} - \mu _2^T\Sigma _2^{ - 1})x \lt {1 \over 2}\ln {{|{\Sigma _1}|} \over {|{\Sigma _2}|}} + {1 \over 2}(\mu _1^T\Sigma _1^{ - 1}{\mu _1} - \mu _2^T\Sigma _2^{ - 1}{\mu _2}){\rm ,}\end{equation}

with mean and variances estimated by their respective sample estimates.  

Classification uses the known label (that is, the known group association) when constructing the classification rules. Misclassification errors (i.e., the number of times an observation from group 1 is classified into group 2 and the number of times an observation from group 2 is classified into group 1, expressed as a fraction of the number of items classified) can be calculated. Within-sample and out-of-sample misclassification errors can be obtained. For within-sample errors, the data are used “twice”—first for developing the classification rules and then for evaluating the performance. For out-of-sample errors, a subset of the data is used to develop the classification rules while a different subset is used in the evaluation. 

Fisher's discriminant analysis method is one of the simplest approaches to classification. It makes the assumption that the k-variate distributions are normal—an assumption that is often violated. Several nonparametric classification approaches that do not rely on normality are available, among them the Support Vector Machine (SVM) approach. We refer the interested reader to the contributions by Boser et al.¹⁰ and Christianini and Shawe-Taylor,¹¹ and readily available computer software such as Package e1071 in the R Statistical Software⁵ and LIBSVM, an extensive library for SVMs. 

Iowa Ophthalmology Patients: Normal Subjects and Glaucoma Patients

Data on Iowa 105 glaucoma patients and 55 normal subjects are analyzed. Subjects were participants in a multiyear Variability in Perimetry (VIP) study directed by Michael Wall, MD, and funded by the Veterans Administration Rehabilitation Research and Development Division. Adult patients with glaucoma were recruited from the glaucoma clinic at the University of Iowa Department of Ophthalmology and Visual Sciences. The research was approved by the University of Iowa and Veterans Affairs institutional review boards, and informed consent was obtained from all patients after an explanation of the study. Inclusion criteria for cases were as follows: (1) the clinical diagnosis of glaucoma determined by the presence of glaucomatous optic disc changes confirmed through the examination of fundus images and (2) the presence of visual field defects (mean deviation [MD] between −2 and −20 dB on standard automated perimetry, as well as either three adjacent locations in a clinically suspicious area falling outside the deviation limits compared with normative data at P < 0.05 [fifth percentile] or two adjacent locations at P < 0.01 [first percentile]). Cases were not required to have elevated intraocular pressure, but they were excluded if there were cataracts causing visual acuity worse than 20/30, they were younger than 19 years, or they had a pupil size of less than 2.5 mm. If both eyes qualified, one eye was randomly selected for inclusion in the study. 

Healthy adults (controls) were recruited using advertisements inviting participation in a research study. Inclusion criteria for controls were as follows: (1) no history of eye disease, (2) refractive error within a ±5-diopter sphere and ±2-diopter cylinder, (3) no history of diabetes mellitus or systemic arterial hypertension, and (4) a normal ophthalmologic examination result, including 20/30 or better visual acuity. An examination by an ophthalmologist on the day of testing or the results of a complete eye examination within 2 years of the testing date was used to confirm normal ocular health. One eye was randomly selected for the study. If a control subject developed a pattern of visual loss leading to an ophthalmologic diagnosis other than refractive error, the individual was not included in the analysis. This research adhered to the tenets of the Declaration of Helsinki. 

Visual field testing was performed with automated threshold perimetry (Humphrey Visual Field 24-2 SITA Standard test with size III, Zeiss-Meditec, Dublin, CA). Field test results were excluded from the analysis if they were unreliable as defined by excessive false-positive (>10%), false-negative (>33%), or fixation losses (>33%). 

All glaucoma patients were being adequately treated and monitored by the glaucoma service at the University of Iowa Department of Ophthalmology and Visual Sciences during the course of the study. In addition to measurement of visual field sensitivity in one eye with various perimetry testing methods, all subjects underwent spectral domain OCT testing (Cirrus, Zeiss-Meditec) on the same day. A 200 × 200 volume scan of the macula and disc were obtained at each test date and were included if a signal-to-noise ratio was at least 7:10. The thickness of the average ganglion cell layer (GCL) complex derived from the macula scan and the average retinal nerve fiber layer (RNFL) derived from the optic disc scan were measured and analyzed over time. We omitted subjects with three or fewer observations. In addition, we omitted eight glaucoma patients with unreliable segmentation of the GCL. The median number of observations per subject was 10 (mean = 9.01; minimum = 4; maximum = 11) for the glaucoma group and 11 (mean = 10.38; minimum = 5; maximum = 12) for the normal group. The median time span from the first to the last visit (range) per subject was 3.99 years (mean = 3.65; SD = 0.90) for the glaucoma group and 4.12 years (mean = 4.00; SD = 0.80) for the normal group. The median time between consecutive visits was 0.50 years (mean = 0.46, SD = 0.20) for glaucoma patients and 0.48 years (mean = 0.43, SD = 0.23) for normal subjects. 

Results

The Random Effects Trend Model

The model in equation (3) is estimated for RNFL and GCL, and the results are shown in Table 1. The assumption of uncorrelated errors is appropriate for this particular data set as the autocorrelations of the errors in the trend regression are all small (average lag 1 autocorrelation −0.19 for RNFL and −0.15 for GCL). The estimation results (model M1) for RNFL show that the measurement variability is about the same in each group, but that the subject-variability among the baselines and also the slopes is larger in the glaucoma group than in the normal group. The average baseline (RNFL thickness, intercept) in the glaucoma group is significantly thinner (by 25.87 micron units) than the average baseline in the normal group. There is evidence for a negative average slope in the normal group (0.145 micron unit reduction per year, with SE = 0.089; P = 0.107). The average rate of change for the glaucoma group is actually smaller (0.145 – 0.089 = 0.056 micron unit reduction per year), but the difference is not statistically significant. The second model M2 with a common average slope implies a reduction of 0.109 units per year (SE = 0.069; P = 0.11). In summary, we find that the average RNFL baseline of the glaucoma group is considerably thinner than the average RNFL baseline of the normal group (as expected), some evidence for an overall average reduction over time (probability value around 0.10), no statistically significant difference between the average slopes in the normal and the glaucoma groups, and larger subject variability among the slopes in the glaucoma group. 

Table 1

 

Estimation Results

 

The estimation results (model M1) for GCL show that the measurement variability is roughly the same in each group, but that the subject-variability among the baselines and also the slopes is larger in the glaucoma group than in the normal group. The average baseline of GCL thickness in the glaucoma group is significantly thinner (by 15.85 micron units) than the one in the normal group. We find statistically significant evidence for a negative average slope in the normal group (0.291 micron unit reduction per year, with SE = 0.041). The average rate of change for the glaucoma group is slightly larger (0.291 + 0.027 = 0.318 micron unit reduction per year), but the difference is not statistically significant. In summary, we find the average GCL baseline in the glaucoma group very much and significantly smaller than the average GCL baseline of the normal group, evidence for a significant overall average reduction over time (0.297 per year), but no statistically significant difference between the average slopes in the normal and the glaucoma groups. Furthermore, we find larger subject variability among the slopes in the glaucoma group. 

The results—no difference in the average slopes across the two groups and larger variability in the glaucoma group—imply that a strategy that uses just slopes will be unsuccessful in classifying subjects into either the normal or the glaucoma groups. The only distinguishing feature that reliably classifies subjects into one of these two groups is the baseline measurement at the start of the study. 

Classification Using Subject-Specific Least Squares Estimates

We use data on 104 glaucoma patients (we omitted one glaucoma patient with missing GCL measurements) and 55 normal subjects that have been analyzed above. Least squares estimates of baselines and slopes from the regression model in equation (1) are used to obtain the Fisher linear and quadratic discriminant functions (discussed in the Methods section) for classification. Scatter plots of baselines and slopes, for GCL and RFNL, are shown in Figures 1 and 2. Normal subjects are shown as red dots and glaucoma subjects are shown as blue dots. Subjects that are misclassified by the linear (quadratic) discriminant function are displayed with an open circle in a color that reflects the incorrect classification. A red dot with a blue circle indicates a normal subject incorrectly classified as coming from the glaucoma group. A blue dot with a red circle indicates a glaucoma patient incorrectly classified as coming from the normal group. 

Figure 1

 

Scatter plots of intercepts and slopes for GCL. Normal subjects are shown as red dots and glaucoma subjects are shown as blue dots. Subjects that are misclassified by the linear (quadratic) discriminant function are displayed with an open circle in a color that reflects the incorrect classification. A red dot with a blue circle indicates a normal subject incorrectly classified as coming from the glaucoma group. A blue dot with a red circle indicates a glaucoma patient incorrectly classified as coming from the normal group.

Figure 2

 

Scatter plots of intercepts and slopes for RNFL. Normal subjects are shown as red dots and glaucoma subjects are shown as blue dots. Subjects that are misclassified by the linear (quadratic) discriminant function are displayed with an open circle in a color that reflects the incorrect classification. A red dot with a blue circle indicates a normal subject incorrectly classified as coming from the glaucoma group. A blue dot with a red circle indicates a glaucoma patient incorrectly classified as coming from the normal group.

The cross-classification matrix and the misclassification rates of linear and quadratic discriminant functions are summarized in Table 2. Quadratic discriminant functions perform best. The misclassification rates are 10.7% for GCL and 6.9% for RNFL. It is the baseline that helps classify subjects. Quadratic discriminant classification on just the baseline leads to misclassification rates of 9.4% (for both GCL and RNFL), while the misclassification rates when using just slopes are 48.8% (for both GCL and RFNL). For this particular group of treated glaucoma subjects and the group of the age-matched normal subjects, it is the baseline that best classifies the subjects, and not the slope. 

Table 2

 

Cross-Classification Matrix and Misclassification Rates of Linear and Quadratic Discriminant Functions

 

The estimated slopes cannot distinguish the two groups. This is because the average slopes are not significantly different (see results in Table 1), and there is no separation between the two distributions. However, the clinician may not be interested in such classification. Instead, he may be interested in identifying known glaucoma subjects with a reduction in either RNFL or GCL (progressive worsening in retinal structure) that is larger than what is expected for the normal group. Results in Table 1 show that the standard deviation of the slopes in the glaucoma group is about twice as large as the standard deviation of the slopes in the normal group. A rule that uses the 10th percentile of the normal slope distribution as the cutoff will certainly identify more than 10% of glaucoma patients with slopes smaller than that cutoff (negative slope, which then can be targeted for treatment). 

Discussion

Clinicians base decisions on short records of consecutive measurements. This paper discusses useful methods for analyzing such data. Random effects trend models have been applied previously to clinical data sets where progression of disease is the focus, such as glaucoma.^17–22 Model results assess whether there exist significant trends and whether average trends are different across groups. 

The paper also discusses whether clinical measures, such as the baseline and the change of the OCT thickness of the RNFL and the retinal GCL, classify patients reliably into disease groups, given their variability. Interestingly, the glaucoma data set from the Iowa glaucoma study used as an example here did not show a significant change (slope) in structure over time compared to the normal group. This was not unexpected, since the glaucoma patients enrolled were already on maximal treatment and were being carefully monitored in the glaucoma clinic over the study period. An untreated (or not adequately treated sample) would be expected to have very different results, showing more subjects with a significant, negative slope. 

In this study, it was also found that more than 10% of the glaucoma patients enrolled showed slopes larger than the 90th percentile of the normal slope distribution (positive slope, or increasing retinal layer thickness over time). This was an unexpected result, since the thickness of the inner retinal layers is usually stable with glaucoma treatment or becomes less over time in the case of progression and is not traditionally expected to increase in thickness over time in glaucoma. However, this concept may need to be reconsidered. Although measurement variability of retinal layer thickness can occur due to factors that influence segmentation of the retinal layers by software (e.g., signal strength, blood vessels), such factors would not fully explain trend changes resulting in systematic increases in the retinal layer thickness over time (resulting in an increase in slope), which was greater in the glaucoma group compared to age-matched normal subjects. The retinal layers are not static; the retinal GCL complex contains the cell bodies of retinal ganglion cells and dendritic connections. Structural changes associated with neuroplasticity in response to the glaucomatous process may include dendritic sprouting and changes in cytoplasm content and could conceivably cause an increase in thickness.¹² Similarly, the RNFL contains the axons of the retinal ganglion cells, and their axoplasmic contents may be influenced by the characteristics of axoplasmic flow and status of the axon,^12–16 which in turn may cause dynamic changes in the thickness of this layer. The present analysis raises important questions about the increase in variance of the distribution of slopes of the GCL and RNFL in glaucoma patients compared to normal subjects. Specifically, future investigation is warranted to better understand the reason why some subjects with glaucoma show an increase in thickness of the layers of the inner retina over time that the normal group does not. 

The statistical approach presented in this study could also be applied to visual field metrics collected at baseline and monitored over time. The example presented here for classification of individual subjects into the diseased group (glaucoma in this case) or normal group could be applied to visual field data as well as the combination of functional and structural data features for a variety of disorders in addition to glaucoma. A similar approach could also be used to classify patients into a progression versus a nonprogression group. With ever more available data, classification plays a central role of personalized medicine. Classification requires more than significant differences among average trends across groups as the subject variability in the trend progression determines whether the classification will be successful. 

Acknowledgments

Michael Wall, MD, kindly allowed us to make use of the OCT imaging data for this manuscript that was collected as part of his Variability in Perimetry (VIP) study, funded by the Department of Veterans Affairs Rehabilitation, Research, and Development Division. 

Randy Kardon is supported by grants from the Department of Veterans Affairs Rehabilitation Research & Development (RR&D) Division; Center for the Prevention and Treatment of Visual Loss, C9251-C, Chronic Effects of Neurotrauma Consortium (CENC); DOD, CENC0056P; C1786-R; 1 IO1 RX002101 VA-ORD; 2 I01 RX000889-05A2; National Eye Institute (NEI). 

Disclosure: J. Ledolter, None; R.H. Kardon, None