A joint multivariate mixed-effects model was implemented within a Bayesian hierarchical modeling framework to evaluate the relationship between the two longitudinal measures obtained over time (i.e., the GDx RNFL measurements and SAP VFI).
32 –35 Linear mixed models were used to evaluate the evolution of each response over time. In these models, the average evolution of a specific response is described using some function of time, and subject-specific deviations from this average evolution are introduced by random intercepts and random slopes, allowing for different baseline values and different rates of change for each patient. Linear mixed models are a natural extension of Bayesian hierarchical models, where the first level of hierarchy corresponds to within-patient variation and the second level to between-patient variation. The Bayesian framework, however, allows more flexible specification of the model assumptions and easier implementation of the joint modeling approach, as described below.
36 –42
In a joint modeling approach using mixed models, random-effects are assumed for each response process and the different processes are associated by imposing a joint multivariate distribution on the random effects. This approach has great flexibility and allows the joining of models for responses of different types and also with a different number of observations, as in our application.
32
Each response process is described using the following linear mixed-effects model
where
Y1 and
Y2 represent the longitudinal measurements of SLP and SAP, respectively, for a subject
i taken at time
t, and
m1(
t) and
m2(
t) represent the average linear evolution of each response over time, that is, the average rate of change of SLP measurements and the average rate of change in SAP sensitivity. Both response trajectories are tied together through a multivariate distribution of the random effects, where
a1 and
b1 correspond to random intercepts and random slopes, respectively, for process 1 (SLP), and
a2 and
b2 correspond to random intercepts and random slopes for process 2 (SAP). By adding an additional level for eye nested within patient, it is also possible to use information from both eyes of the same patient taking into account the correlation between them.
Conventional linear mixed models assume a normal distribution of random effects. However, previous studies have shown that the general assumption of normally distributed random effects may lead to biased estimates of individual change parameters when there is heterogeneity in the population.
43,44 In the present application, heterogeneity is to be expected, as only a proportion of eyes will show progression over time. Further, in the progressing group, only a small proportion is expected to have relatively fast progression. This situation can induce considerable nonnormality, or skewness, in the random effects distribution. To address the problem, we used a multivariate skew
t distribution to model random effects. Both the
t and skew
t distributions allow greater flexibility in the distribution of random effects compared with the normal distribution and have been successfully applied to model nonnormal random effects.
41,45 The probability density function for a random effect
bj that is
t distributed
t(μ,
v,
k) is
where μ is the mean,
v the scale parameter, and
k the degrees of freedom determining the weight of the tails, giving variance σ
2 = [
k/(
k − 2)]
v. The
t distribution can be extended to accommodate the multivariate case.
46 To allow for skewing, Fernandez and Steel
47 proposed the introduction of skewness to any unimodal distribution symmetric ∼0, using a scale factor γ on each side of 0. For a random variable −∞ <
x < ∞, this method gives:
where
I[c,d] is the indicator function for
x, being between
c and
d, and γ controls the mass on each side of 0, representing the skewness.
Therefore, in our Bayesian model, the intercepts and slopes were assumed to follow a multivariate skew
t distribution. This multivariate distribution acted as a “prior” for the estimation of the intercepts and slopes for each eye and the parameters of the multivariate distribution were themselves estimated from the data. However, to perform Bayesian inference, we also need to specify priors for the parameters of the multivariate distribution (i.e., the hyperparameters). Unless there are strong prior beliefs in the hyperparameter values, it is usually desirable to use noninformative priors. Therefore, we used normal (0,1000) for μ and uniform (0,100) for σ.
48 An exponential (0.1) prior was used for
k in the
t distribution,
47 but restricted to
k > 2.5 to prevent problems with undefined variance when
k ≤ 2. For the skew parameter γ, we used a γ(0.5, 0.318) prior for φ = γ
2. These are general priors used in previous applications of Bayesian hierarchical models using
t and skew
t distributions.
45
Estimates of the posterior distributions of the parameters of interest (i.e., the random effects, were obtained by Markov chain Monte Carlo (MCMC) procedures. The MCMC sampler was implemented in WinBUGS software.
49 We used 10,000 iterations after discarding the initial 5,000 iterations for burn-in. Convergence of the generated samples was assessed by standard tools in WinBUGS (trace plots and autocorrelation function [ACF] plots) as well as Gelman-Rubin convergence diagnostics. After the posterior distributions were estimated, summary measures were calculated, such as mean and credible intervals. For a specific test, we considered that progression had occurred if the upper limit of the 95% credible interval for the slope was less than 0.
Importantly, no information on the classification of subjects (i.e., whether there was optic disc progression on stereophotographs or not) was inserted into the model. Also, the likelihood information from healthy eyes was not allowed to influence calculations of model parameters. This was done to have a model that would reflect a population that is commonly followed under a clinical scenario, that is, composed of eyes with glaucoma or suspected of having disease. We calculated sensitivity, specificity, and receiver operating characteristic (ROC) curves as measures of discrimination. To estimate the effect of the sample composition on the predictive distribution, we ran 200 bootstrap resamples and re-estimated sensitivity and specificity. To avoid potential biases in the estimation of sensitivity and specificity, we also used a leave-one-out approach to calculate the estimates for individual eyes. That is, the model was built with information from
n − 1 eyes and applied to obtain the estimates of progression in the
nth eye. We compared our algorithm for detection of progression to standard approaches reported on the printouts of the instruments and used in previous studies.
8,10,25,50 –53 The standard approach is based on OLS linear regression of measurements over time, so that an eye is considered to have progressed if a negative OLS regression slope is significantly different from 0, with
P < 0.05. OLS regressions were performed separately for each eye and for each test (i.e., SAP VFI and SLP TSNIT average).