Combining optical coherence tomography with visual field data to rapidly detect disease progression in glaucoma: a diagnostic accuracy study

Overview

Evaluation of ONH and RNFL imaging measurements as potential clinical trial outcomes was a secondary objective of the UKGTS. 60 The trial protocol for data collection was established before the index and reference standard tests were performed.

This work draws together a number of different statistical strands to establish statistical tools to analyse glaucomatous disease deterioration. The main emphasis was to investigate ‘trend’ analyses, which use several observations taken over a period of time to estimate a rate of change in the outcome over that time period. A related secondary set of methods that were investigated were ‘event’-based analyses, in which patients are declared to have experienced an event if they have deteriorated by a certain amount compared with their baseline measurements. With the reference method (see below), an event is identified when follow-up measurements differ by a certain amount from the baseline measurements. With the index methods, the event is based on the ‘trend’ analysis, so that the event is the point in time at which a trend becomes statistically significant.

Before evaluating methods that combine the VF and OCT measurements, we evaluated some data attributes required to assess the suitability of the imaging outcome as an alternative or additional outcome to VF testing. To do this we compared the distribution of individual rates of change in VF mean sensitivity and OCT mean RNFL thickness between treatment arms of the UKGTS. Potential surrogate outcomes also need to be shown to capture the effect of a treatment intervention on the clinically relevant outcome. 64,65 We therefore assessed the association of the rate of RNFL loss in individual participants with the time to VF event in a Cox proportional hazards model.

The reference method to identify VF deterioration, for comparison with the index methods, was the event-based method available in the VF testing instrument [Humphrey Field Analyzer™ (HFA) II-i Guided Progression Analysis™ (GPA) software (Carl Zeiss Meditec Inc., Dublin, CA, USA)] with the criterion for deterioration that was applied in the UKGTS. 60,67 The GPA method compares each follow-up VF with the baseline pair of VFs and ‘flags’ deterioration when the follow-up measurement at a VF location has a sensitivity that is lower than the 5% limit predicted from population test–retest data. Criteria for deterioration are based on various combinations of the number of flagged location and the number of consecutive times that the location is flagged. 68

Two types of index analysis methods were considered: those that analyse VF measurements alone (‘VF-only index methods’) and those that analyse both VF and OCT outcomes (‘VF + OCT index methods’). These methods were compared against the reference method.

Analysis methods were compared as detailed in the following sections.

Review of models to identify visual field deterioration and incorporate imaging outcomes

Testing of the VF, using instruments such as the HFA, has been used as both a diagnostic tool and an outcome measure in randomised clinical trials (RCTs) in glaucoma for many years. Bosworth et al. 71 reviewed current VF testing practice. There has been much research into the characteristics of VF measurements from the 24-2 pattern acquired with the Swedish interactive thresholding algorithm (SITA) of the HFA. 72 Many papers have described the large variability in VF measurements and the increasing variability as sensitivity declines. 7–9 Gardiner et al. 73 found that sensitivity values below a cut-off value of between 15 and 19 dB are particularly unreliable and that values at these levels of sensitivity agree poorly with estimates obtained from ‘frequency of seeing’ curves. Gardiner et al. 73 also discussed the concept that sensitivity values may not be truly defined at such low levels, as the retinal ganglion cells remaining cannot be stimulated sufficiently to produce a 50% response rate to stimuli, regardless of the intensity of the presentation. Gardiner et al. 74 found that the proportion of eyes identified as deteriorating was not decreased by truncating sensitivities at 10 dB (i.e. setting all sensitivities below 10 dB equal to 10 dB) in the two cohorts studied. Furthermore, in one of the cohorts, no reduction in the proportion deteriorating was observed when truncating at 15 dB; in the other, only a small reduction was observed.

In addressing the increased variability in sensitivity estimates with declining sensitivity, Ibáñez and Simó75 used geostatistical methods to model the spatiotemporal characteristics of VF data taken from healthy eyes. They noted the relationship between variability, mean sensitivity and distance from the centre of the VF and considered the use of a Box–Cox transformation to adjust for the heteroskedasticity of the data.

Bryan et al. 76 compared pointwise simple linear regression (SLR) with exponential regression, censored regression and median-based regression to describe VF deterioration. They found that all models performed similarly, with SLR performing slightly better in terms of prediction errors and model fit, but that all models gave large prediction errors. However, the authors noted a need for more complex models that explore the spatiotemporal relationships of VF data to improve predictions further. For the main part of their paper, Bryan et al. 76 considered censored regression models (also known as Tobit models) with censoring at 0 dB, which is the limit of the VF machine’s dynamic range. In the discussion section, the authors briefly mention also considering censoring at 10 dB, but that SLR still performed better in terms of model fit and prediction errors. The authors acknowledged the difficulties of using prediction errors as a metric to compare models given the high level of variability in VF measurements, particularly at low dB values.

Bryan et al. 77 explored the use of mixed-effects models in a Bayesian framework in the context of analysis of glaucoma data. In their model they included levels for person, eye and hemi-field within eye. In such a model, the VF locations are assumed to have an exchangeable correlation structure within the hemi-field. The authors also included a ‘global visit effect’ in their model (to account for visit-by-visit differences in performance) and allowed the log of the SD in their model to vary linearly with the sensitivity. Bryan et al. 77 used a two-stage Bayesian Monte Carlo procedure to obtain estimates from such a model. They found that by taking into account the global visit effect and including the relationship between variability and sensitivity, the model fit improved, as well as future VF prediction accuracy.

Acknowledging the heteroskedastic nature of VF variability, O’Leary et al. 78 considered the use of permutation methods. They used SLR to obtain a p-value at each VF location. They then used a generalisation of Fisher’s method to combine the p-values into a summary statistic, which they permuted. With this method, which they named permutation analyses of pointwise linear regression (PoPLR), they obtained a superior ‘hit rate’, with specificity held at 5%, to the previously described pointwise linear regression. 79 A particular advantage of this method over the GPA facility in the HFA is that the statistical significance of the slope of change is estimated from the patient’s own data rather than from population statistics.

Zhu et al. 80 proposed an approach to identify deterioration that formally modelled the variability characteristics at each level of VF sensitivity; this was termed ‘Analysis with Non-Stationary Weibull Error Regression and Spatial enhancement’ (ANSWERS). The methodology was developed from a data set that included 12 VF tests on each of 30 subjects. 81 To develop their methodology, Zhu et al. 80 formed the 1980 (30 × 12 × 11/2) pairs of VFs and treated these as independent representations of within-subject variability. However, these pairs were not truly independent as they came from only 30 subjects; it is not clear to what extent this non-independence may impact on the validity of the method, but the sensitivity-determined variability characteristics of the model were consistent with previous publications. 7,8 The ‘hit rate’ and prediction accuracy of ANSWERS was compared with the SLR of a VF summary measure [mean deviation (MD)] and with PoPLR. The hit rate was significantly better than the SLR of the MD and PoPLR, especially in short time series. The prediction accuracy of ANSWERS was also better than that of PoPLR, particularly in shorter series; 75% of VF series were better predicted by ANSWERS than by PoPLR and the average prediction error of ANSWERS was 15% lower than that of PoPLR. 82

Medeiros et al. 83 proposed a Bayesian framework to integrate event-based (GPA) and trend-based (SLR of VFI) analyses. The results from the event-based analysis were incorporated into the prior for the trend-based analysis. Although this approach identified more subjects as deteriorating than the GPA alone, at the same level of estimated specificity, the approach appears to use the same data twice, which may lead to incorrect estimates of variability when making inferences.

Recognising the potential of structural measures to support VF testing, several groups have proposed methods to combine VF and imaging data for glaucoma diagnosis and disease staging84–86 that are reported to perform better than the separate structure or function tests. Medeiros et al. 87 applied their approach to identify glaucoma deterioration. The analysis is effectively SLR on combined summary measures of glaucoma damage and does not address the heteroskedasticity in the VF data. Analysis of the combined data provided a better hit rate than the VF and structure data alone, for the same level of specificity. This relatively simple methodology could be more statistically efficient than more complex approaches; however, it is probably more appropriate to exploit the multivariate nature of the repeated measures of the imaging outcome together with the data from all locations in the VF.

Two further studies by Medeiros et al. 62,88 used a Bayesian approach to combine structural and functional information. In the first study,62 a Bayesian hierarchical model was used to integrate information from the longitudinal VF (VFI) and structure (RNFL thickness estimated with the scanning laser polarimeter) measurements to classify individual eyes as progressing or not. The output was compared with SLR of the component VF and structural data; more eyes were identified as deteriorating with the Bayesian approach. ‘Specificity’ was estimated from healthy eyes. In the second study,88 a Bayesian joint regression model was used to integrate structural (rim area measurements from a scanning laser tomograph) and functional (VF MD) information. Compared with SLR on either the structural or the functional components alone, the deterioration slopes were less variable and prediction accuracy was better. Both approaches use a summary VF outcome measure (VFI or MD) and neither formally addresses heteroskedasticity in the data. Furthermore, the use of healthy eyes to evaluate criterion specificity is not appropriate because the range of sensitivity values and associated variability differs.

Russell et al. 89 proposed a Bayesian framework in which the rate of neural rim change measured by scanning laser tomography formed the prior for SLR of a VF summary measure (mean sensitivity). The Bayesian approach incorporating the structural data resulted in significantly better predictions of future VF mean sensitivity than SLR of mean sensitivity alone. As with other published methods, this approach used a summary VF outcome measure and did not addresses heteroskedasticity in the VF data. However, this, and the other approaches, demonstrate the potential benefit of including structure outcomes in the analysis of serial VF data to identify deterioration.

Developing analysis approaches for the visual field and imaging outcomes

In this project we made use of a variety of techniques that are used in mainstream medical statistics in many disease areas, as well as building on some methods that have also been used in a glaucoma setting; these are detailed in Chapter 6 (see Index methods: newly developed). For example, we considered a variety of transformations, including those from the two-parameter Box–Cox family of transformations,90 to investigate whether or not one could be found under which the transformed VF values would follow an approximate normal distribution. Ibáñez and Simó75 considered the use of Box–Cox transformations and transformations more generally have been used in a wide variety of settings to allow normality assumptions to be made.

We considered the use of censored regression models (also known as Tobit models) in this project, using a cut-off value of 15 dB. This relates to work by Bryan et al. ,76 among others, who considered regression models that assume censoring below 0 dB. We considered a higher cut-off value, as described in Chapter 6 (see Censored regression), which has the dual advantage of allowing normality and homoskedasticity assumptions to be made and avoids the need to account for zero weighting caused by the truncation at 0 dB. The recent study by Gardiner et al. 74 considered ‘censoring’ at a cut-off value above 0 dB, although in their work they simply truncated at the cut-off value and replaced all values below the cut-off value with the cut-off value, instead of using a censored regression model as proposed.

Another method we have made use of is the permutation test. These tests are commonly used in many areas of statistics to compute p-values when comparing levels of a particular outcome variable between two groups when there are concerns over the assumptions made by a test such as a t-test. 91 Such permutation tests have been utilised for the analysis of serial VF data to identify deterioration. 78

Multiple imputation, and specifically the use of chained equations,92 is frequently used to handle missing data in a principled way in many disease areas. A similar approach to the model that we propose in Chapter 6 (see Multiple imputation) has been previously used to analyse longitudinal dietary data. 93

Linear mixed-effects models (also known as hierarchical models or mixed models) are commonly used in other disease areas to analyse longitudinal and/or hierarchical data. 94 They have also been applied to glaucoma data by Bryan et al. 77 in a Bayesian framework. The specific use of a Kronecker product to specify the residual error structure is described by O’Brien et al. 95 Generalised estimating equations (GEEs) are also a standard methodology that is frequently used in other disease areas. 94

We also make use of the known spatial relationship between regions of the VF and RNFL sectors around the ONH. 22

United Kingdom Glaucoma Treatment Study

The UKGTS was a RCT that compared the effects of latanoprost, a topical treatment to lower IOP, with those of placebo on survival from VF deterioration. In total, 516 patients with newly diagnosed OAG were enrolled, with 777 eyes eligible for entry into the study; details of the eligibility criteria and baseline characteristics have been published elsewhere. 60,96 The study was undertaken in accordance with good clinical practice guidelines97 and adhered to the Declaration of Helsinki. 98 The trial was approved by the Moorfields and Whittington Research Ethics Committee on 1 June 2006 (reference 09/H0721/56). All patients provided written informed consent before the screening investigations were carried out. An independent Data and Safety Monitoring Committee (DSMC) was appointed by the Trial Steering Committee. The trial manager monitored adverse events, which were reported immediately to the operational DSMC at Moorfields Eye Hospital. Serious adverse events were reported to the Medicines and Healthcare products Regulatory Agency. The trial is registered as ISRCTN96423140.

The principal baseline characteristics of the participants are presented in Table 1.

Participants were followed up every 2–3 months after eye drop therapy was initiated, for up to 11 scheduled visits (see Appendix 1). Participants attended for additional visits if VF deterioration was identified according to certain preset criteria, at which VF testing and imaging were repeated. The baseline visit was the first visit, 6 weeks after randomised therapy (drops) was started. Visual function was monitored by VF testing (see Data types, Visual field measurements) and ONH structure was monitored with the Heidelberg Retina Tomograph at all study sites and with TD Stratus OCT™ (Carl Zeiss Meditec Inc., Dublin, CA, USA) (software version 5.0; see Data types, Optical coherence tomography measurements) and GDxECC Nerve Fiber Analyzer (Carl Zeiss Meditec Inc., Dublin, CA, USA) at study sites with those devices. The sample size in the UKGTS was determined for a two-sided error at a significance level (α) of 0.05 to detect the difference between 24% and 11% incident deterioration over a 24-month follow-up at 90% power and assuming a 25% attrition rate. The subset of UKGTS participants with both VF testing and OCT imaging was used in this work.

The primary outcome for the trial was glaucomatous VF progression (deterioration) within 24 months. Details of the method for determining progression in the VFs (see Chapter 3, Visual field measurements) have been published previously. 60,67 If tentative deterioration was identified, participants returned for confirmation tests within 1 month. At this confirmation visit, two VF tests were performed; if the deterioration was confirmed, then participants were considered to have progressed. Participants who were deemed to have progressed left the trial and treatment was adjusted as appropriate. Participants leaving the trial were invited to an ‘exit visit’ before treatment adjustment. Participants found to not be progressing at the confirmation visit were returned to the standard visit schedule.

The primary outcome was analysed as survival in a Cox model, which found an adjusted treatment hazard ratio (HR) of 0.44 [95% confidence interval (CI) 0.28 to 0.69, p < 0.001]. Details of the covariates used in the outcome model have been published previously. 67

RAPID study

Eighty-two glaucoma patients under standard treatment were recruited to a test–retest study. Seventy-seven (148 eyes) of the patients recruited attended for up to 10 visits within a 3-month period, totalling 1256 patient-eye visits. This data set was taken to represent a ‘stable glaucoma’ cohort; assumptions made include that, over such a short length of time, no clinically meaningful changes in the VF or RNFL structure would occur and that the variability in characteristics of the VF and RNFL measurements are similar to those seen in clinical practice over longer periods of time.

The study was undertaken in accordance with good clinical practice guidelines97 and adhered to the Declaration of Helsinki. 98 The study was approved by the North of Scotland National Research Ethics Service committee on 27 September 2013 (reference 13/NS/0132) and NHS Permissions for Research was granted by the Joint Research Office at University College London Hospitals NHS Foundation Trust on 3 December 2013. All patients provided written informed consent before the screening investigations were carried out.

Recruitment criteria were based on those for the UKGTS. 60 Patients were required to have reproducible VF loss with corresponding damage to the ONH and no other condition that could lead to VF loss, be aged > 18 years and have a visual acuity of ≥ 6/12, a refractive error within ±8 dioptres and an IOP of ≤ 30 mmHg. The VF MD had to be better than –16 dB in the worse eye and better than –12 dB in the better eye. VF loss was defined as a reduction in sensitivity at two or more contiguous locations with p < 0.01 loss or more, three or more contiguous locations with p < 0.05 loss or more, or a 10-dB difference across the nasal horizontal midline at two or more adjacent locations in the total deviation plot.

Participants attended approximately once a week for 10 visits, with VF testing and OCT imaging carried out twice at the first visit and once at each subsequent visit. The mean time between visits was 8.2 days (range 3–63 days). VF testing was undertaken with the HFA, as detailed below, and OCT imaging was carried out using TD Stratus OCT and Spectralis^® SD OCT (Heidelberg Engineering, Heidelberg, Germany); we present the TD OCT results in this report because this was the form of OCT imaging used in the UKGTS.

Halifax study

Initial statistical modelling was undertaken while the RAPID study was still in the data collection phase. Therefore, we made use of a similar test–retest data set. 81 This consisted of 30 glaucoma patients who took 12 VF tests over a 3-month period, approximately once a week, in the Department of Ophthalmology, Dalhousie University, Halifax, NS, Canada. In accordance with the Declaration of Helsinki,98 the institutional research ethics board approved the protocol and all patients gave written informed consent.

Visual field measurements

In a VF test, a patient fixates with the eye to be tested (one at a time) on a central (‘fixation’) point and is then presented with flashes of light of varying intensity at locations at various distances (eccentricity) from the fixation point. The patient is provided with a button and instructed to click this when he or she can see a flash of light. In the standard 24-2 VF pattern, 54 locations in a regular grid 6° apart are tested and a sensitivity level is recorded for each. The locations above and below the physiological blind spot (see Figure 1) are discarded so that 52 locations are analysed. At a sensitive retinal location a dim flash can be seen and at a location with poor sensitivity only bright flashes can be seen. The unit of sensitivity measurement (decibel) is 10 times the log of the reciprocal of the dimmest intensity seen, so that a sensitive location has a high decibel value. A sensitivity of < 0 dB implies that the patient has been unable to see the brightest light that the machine can produce. Thus, the measurements are bounded at 0 dB. They are also heteroskedastic (the variability associated with the measurement depends on the mean level of the measurement), so that the variability of the VF sensitivities increases as sight deteriorates (at low decibel values).

All VF tests were performed with the HFA II (or II-i) and the SITA standard 24-2 program. A reliable VF was one with a false-positive rate of < 15% and < 20% fixation losses (for fixation losses of > 20%, reliability was based on the subjective judgement of the technician supervising the test and the clinician reading the test, including an assessment of the eye tracker trace). Unreliable tests were repeated, either on the same day (with a break of at least 30 minutes) or on a subsequent occasion.

The reference standard analysis for VF deterioration was that used for the outcome of the UKGTS and was undertaken with the HFA II-i GPA software (version 5.1.1). The criterion for tentative deterioration (progression) was three locations worse than baseline in two consecutive VFs (three half-shaded or full-shaded locations). Definite deterioration was identified if the same criterion of three half-shaded or full-shaded locations was satisfied in the next two VF tests. 60,67

Optical coherence tomography measurements

During OCT image acquisition, a patient sits in front of the image device, which is aligned to the patient’s pupil. The patient views a fixation light and the instrument scans the retina, acquiring reflectivity measurements around a circle centred on the ONH. OCT imaging was performed with the TD Stratus OCT. Images were acquired through dilated pupils with the fast RNFL thickness (3.4) scanning protocol and using the landmark function. With this scanning protocol, three images are acquired in quick succession. The OCT instrument software averages the measurements from these three images. A signal strength of ≥ 7 was required; images were retaken if necessary to acquire adequate-quality images. Images of lower quality, or those with a software alert, were not included in the analyses.

Retinal nerve fibre layer measurements are provided as means (average RNFL around the ONH) and in clock-hour sectors.

Structure/function mapping

There is an established anatomical correspondence between regions of the VF and sectors of the ONH (Figure 4). 22 In some analyses we make use of this mapping, so that either regions of the VF or individual VF locations are associated with corresponding RNFL measurements.

FIGURE 4.

Relationship between VF regions and RNFL sectors around the ONH. The optics of the eye result in inversion of the VF compared with the anatomy. The VF test locations are shown in ‘retinal view’ (the VF has been inverted, so that inferior locations are shown at the top of the image, to allow direct correspondence between the anatomy and the VF). The RNFL sectors are numbered as follows: (1) superior 30° sector corresponding to the inferior VF, (2) superotemporal 30° sector corresponding to the inferior arcuate VF region, (3) nasal 120° sector corresponding to the temporal VF, (4) temporal 120° sector corresponding to the central VF, (5) inferior 30° sector corresponding to the superior VF and (6) inferotemporal 30° sector corresponding to the superior arcuate VF region.

Data

A subset of 528 eyes of 361 UKGTS participants had OCT imaging of adequate quality acquired during the follow-up, 178 participants in the placebo group and 183 participants in the latanoprost group. The principal characteristics of these participants (Table 2) were very similar to those of the complete UKGTS cohort (see Table 1).

Figure 5 shows the process for identifying subjects with VF and OCT image data of adequate quality. Confirmation and exit visits were included in the analyses provided that they occurred within 2 years of the baseline visit.

FIGURE 5.

Flow chart illustrating the process for identifying patients with both VF and OCT data of adequate quality.

Data for the analyses in the following two sections came from the 284 UKGTS participants (141 in the latanoprost arm and 143 in the placebo arm; see Figure 5) with adequate-quality VF and OCT data, with > 6 months of follow-up, who underwent five or more visits and with data for both VFs and OCT at the baseline visit.

Rate of change of visual field mean sensitivity and mean retinal nerve fibre layer thickness

The rate of change of VF mean sensitivity (mean of the 52 locations) and mean RNFL thickness were calculated by SLR and the individual rates of change were compared between treatment arms using a Mann–Whitney two-tailed test.

Association of retinal nerve fibre layer thickness change with time to visual field progression

To identify whether the rate of change in OCT RNFL thickness was associated with VF progression (reference analysis), a Cox proportional hazards model was fitted to the data for factors potentially associated with survival failure (treatment allocation, age, baseline IOP, baseline VF MD and occurrence of a disc haemorrhage in either eye during follow-up)99 and the slope of change in OCT RNFL thickness. Calculations were performed using MedCalc Statistical Software version 17.1 [MedCalc Software bvba, Ostend, Belgium; see www.medcalc.org (accessed 7 November 2017)].

Data

Reference analysis

Progression was defined according to the reference GPA (detailed in Chapter 3, Visual field measurements). Progression-free VF survival for the treatment and placebo arms was assessed using Kaplan–Meier survival analysis in two subsets: the 284 UKGTS participants in Visual field and imaging outcomes in the United Kingdom Glaucoma Treatment Study and the 320 participants in Index analyses. Calculations were performed with MedCalc Statistical Software version 17.1.

The specificity of the reference analysis was evaluated in the RAPID study participants; the first two VFs formed the baseline and each subsequent VF test was compared with the baseline pair. The RAPID data were not permuted for this analysis because this is not possible using the HFA II-i GPA software. Estimates of criterion specificity took into account the manner in which tests for progression are applied in clinical practice and in clinical trials. Tests for progression are applied each time a patient has a new VF or OCT test and so there is an opportunity for false-positive identification of progression. Therefore, in every generated series, the ‘progression test’ is applied at each time point and a series is flagged as progressing if the progression criterion is met at any time point.

Index analyses

Structure-guided ANSWERS

Structure-guided ANSWERS is a modification of ANSWERS in which there is a two-layered hierarchical Bayesian model (Bayesian belief network); the prior distribution of the VF progression rate at each VF location is set by the slopes and variance of the rate of change of RNFL thickness measurements (Figure 6). This is similar to the approach described previously to incorporate scanning laser ophthalmoscope rim area measurement slopes into VF progression analysis. 89

FIGURE 6.

Structure of sANSWERS: the structure measure S and function measure F are dependent when the function progression parameter w_f is not conditioned on, but become independent only when w_f is observed (conditioned). F, function measure; w_f, the prior distribution of the slopes and intercepts of VF progression rate; ∑_{f_s}, The prior distribution of the slopes of the VF progression rate; ∑_{f_i}, The prior distribution of the intercepts of the VF progression rate; w_s, The distribution of the rate of change of RNFL thickness measurements; ∑_space, The spatial correlation between each location and all other locations in the visual field; S, structure measure.

As the spatial correspondence of peripapillary circle sectors and VF locations is known,22 each VF location was mapped to one of 12 peripapillary RNFL sector measurements (see Figure 4); the slope and variance of RNFL thickness over time formed the Bayesian prior for the VF slope.

Assessment

Data

Methods and assessment

Visual field test and transformations

As can be seen in Figures 7 and 8, which display data from the Halifax study,81 VF measures are highly non-normal. They are truncated at zero and are heteroskedastic, with greater variation as vision deteriorates (i.e. as VF values decrease). Figure 9 illustrates the complex spatial correlation structure of VF measures. Mixed-effects models, also known as hierarchical models, are commonly used in other disease areas to analyse multivariate, hierarchical data and can handle complex correlation structures such as those observed in VF data. However, they are easiest to implement when outcomes approximately follow a normal distribution. We therefore explored a variety of transformations, including those from the two-parameter Box–Cox family of transformations, to investigate whether one could be found under which the transformed values would follow an approximate normal distribution. 90 Examples of the transformations that we considered can be found in Figures 10–13. None of the transformations that we considered appeared to stabilise the variance sufficiently for the transformed values to be reliably modelled under a normality assumption.

FIGURE 7.

Visual field sensitivity values at each test location from 30 glaucoma patients from the Halifax study. There is one plot for each of the 54 VF test locations tested in a 24-2 pattern (two of the locations are adjacent to the physiological blind spot). Each line on the location-specific plots represents the trajectory of a single patient across 12 repeated visits. Differential light sensitivity (DLS) measured in decibels is plotted against visit number.

FIGURE 8.

Plots of within-person variance vs. within-person mean at each VF location for the 30 patients from the Halifax study.

FIGURE 9.

Correlations between mean light sensitivity values across all VF test locations. This figure shows the correlations between test locations (the mean of the 12 repeat VF measurements at each location) in the 30 patients from the Halifax study. Each plot represents one VF location, which is highlighted in dark blue. The other pixels represent the strength of the correlations with other locations, with darker colours indicating larger correlations. The lightest green shading represents correlations of between –0.6 and 0.4, the mid-green shading indicates correlations between 0.4 and 0.6, the dark green shading indicated correlation between 0.6 and 0.8 and the darkest green represents correlations between 0.8 and 1.

FIGURE 10.

Plots of trajectories over time for log(45 – DLS)-transformed Halifax study data. Log(45 – DLS) are plotted against visit number.

FIGURE 11.

Plots of within-person variance vs. within-person mean for log(45 – DLS) using transformed Halifax study data.

FIGURE 12.

Plots of trajectories over time for Halifax study data with a two-parameter Box–Cox transformation: [(35 – DLS)^0.56 – 1]/0.56. The transformed data are plotted against visit number.

FIGURE 13.

Plots of within-person variance vs. within-person mean for Halifax study data with a two-parameter Box–Cox transformation: [(35 – DLS)^0.56 – 1]/0.56.

Censored regression

An alternative to the transformations is to evaluate the normality of the test–retest distributions over part of the sensitivity range. To do this we utilised a variation of a quantile–quantile plot in which values are plotted against the values that would be expected under a normal distribution given their rank. Such a standard quantile–quantile plot would generally be used for a variable with independent values, for example with each value obtained from a different person. However, in this setting we have multiple values obtained over time, and from several locations across the eye, for each person. Each person and location has values distributed around their own mean level. To avoid the complexity of multiple locations we plotted each separately. To account for every person having their own mean we used the following variation of a quantile–quantile plot:

• For each location separately, a censored regression model was fitted to the 12 repeat measurements for the 30 subjects. This model allowed a different mean for each subject, with values of < 15 considered as censored and a common within-subject variance. As a separate model was fitted for each location, the variance was free to vary between different spatial locations across the eye.

• 12 values were drawn for each patient from a normal distribution centred on the patient’s mean and using the common variance from the censored regression model.

• These were ordered according to their size.

• The previous two steps were repeated 500 times.

• An average was taken across each of these repeats to obtain an accurate estimate of the expected value for each rank.

• The average, or expected, value for each rank was plotted against the actual value.

If the values are in fact normally distributed, or close to normally distributed, then they would be expected to fall close to the line of equality. As can be seen in Figure 14, the VF data appear to be reasonably well described by a mixture of normal distributions above 15 dB, although at a few locations there is some minor departure from normality in the tail of the distribution. Such departures were appreciably more marked using a cut-off value of 10 dB (Figure 15) and not completely eliminated using a cut-off value of 20 dB (Figure 16), providing some justification for our choice of a 15-dB cut-off value. In addition, using a higher cut-off value such as 20 dB results in a higher proportion of censored observations and might lead to a loss of important information.

FIGURE 14.

Joint quantile–quantile plots with a cut-off value of 15 dB for the Halifax study data set.

FIGURE 15.

Joint quantile–quantile plots using a cut-off value of 10 dB for the Halifax study data set.

FIGURE 16.

Joint quantile–quantile plots using a cut-off value of 20 dB for the Halifax test-retest data set.

Visual field

The distribution of the rate of VF mean sensitivity change is shown in Figure 17. The rate of loss was taken from the eye with the worse baseline MD or the eye demonstrating incident VF loss. It can be seen clearly in the histogram that the placebo group has faster rates of deterioration than the latanoprost group (data shifted to the left). The mean (SD) rates of change were –0.36 (2.20) dB per year in the placebo group and –0.02 (0.96) dB per year in the latanoprost group. The d’Agostino–Pearson test for normal distribution of individual rates of change rejected normality (p < 0.0001). The median rate of change was –0.21 (5th to 95th percentile –2.84 to +1.40) dB per year in the placebo group and +0.12 (5th to 95th percentile –1.95 to +1.39) dB per year in the latanoprost group. A Mann–Whitney two-tailed test (independent samples) identified that the distribution of slopes was significantly different (p = 0.0034).

FIGURE 17.

Distribution of the rate of VF mean sensitivity (MS) change in decibels per year for the subset of UKGTS participants with OCT images (placebo, n = 143 participants; latanoprost, n = 141 participants).

Retinal nerve fibre layer thickness

The distribution of rate of RNFL thickness change is shown in Figure 18. The rate of loss was taken from the eye with the worse baseline MD or the eye demonstrating incident VF loss. Similarly to the VF data, the placebo group had faster rates of deterioration than the latanoprost group (data shifted to the left). The d’Agostino–Pearson test for normal distribution rejected normality (p = 0.0026). The median rate of change was –1.56 (5th to 95th percentile –9.17 to +6.16) µm per year in the placebo group and –1.05 (5th to 95th percentile –8.05 to +3.89) µm per year in the latanoprost group. The difference in distribution of slopes was not statistically significant (Mann–Whitney two-tailed test, p = 0.18).

FIGURE 18.

Distribution of the rate of OCT RNFL thickness change for the subset of UKGTS participants with OCT images (placebo, n = 143 participants; latanoprost, n = 141 participants).

‘Hit rate’ compared with specificity

Figure 20 illustrates the proportion of eyes of the 320 UKGTS participants with five or more time points with both VF tests and OCT images available that were identified as progressing (the ‘hit rate’, equivalent to the true positives plus the false positives) plotted against the false-positive frequency (proportion of eyes of RAPID subjects identified as deteriorating) as the criterion for flagging an eye as deteriorating is varied. This is a larger subset than other analyses because baseline OCT images were not required (see Figure 5).

FIGURE 20.

The proportion of participants in the UKGTS identified as progressing (hit rate) plotted against the false-positive frequency as the criterion for progression is varied. The ‘hit rate’ is the proportion of eyes of UKGTS participants identified as deteriorating at criterion false-positive rates between 0% and 15%. Analyses are shown for the ANSWERS, PoPLR and sANSWERS models. Data are shown for series intervals (baseline to final observation) of up to (a) 7 months (404 eyes), (b) 13 months (378 eyes), (c) 18 months (286 eyes) and (d) 22 months (230 eyes). The shorter series are subsets of the longer series so that an eye identified as ‘progressed’ earlier in the series is carried forward as ‘progressed’ in the later series. Data are shown for 404 eyes of 320 participants.

At a 5% false-positive frequency and after 22 months’ observation, the hit rate for the ANSWERS and PoPLR methods was similar, at around 20%. For comparison, the hit rate with the GPA criterion applied in the UKGTS in this subset of eyes with OCT data was 87/394 eligible eyes (22%). The hit rate for the sANSWERS method was considerably greater, at > 40%, suggesting that, for the same false-positive frequency, sANSWERS is much more sensitive at identifying a progressing eye. A similar pattern was seen for shorter follow-up durations, but with the ANSWERS method showing greater sensitivity than the PoPLR method for short follow-up durations.

Prediction of future visual field state

The sensitivity at each location in the last VF in a UKGTS participant’s series was predicted by projecting the trend for sensitivity change observed over the first five visits to the time of the last VF. The period over which the initial trend line was fitted was a mean (SD) of 42.4 (6.2) weeks and the mean (SD) interval from the initial period to the predicted VF was 49.2 (19.8) weeks.

The mean prediction errors (average error across the 52 locations in the VF) across subjects were not normally distributed (D’Agostino–Pearson test, p < 0.0001). The median mean prediction error across subjects was 3.8 (5th to 95th centile 1.7 to 7.6) dB for SLR, 3.0 (5th to 95th centile 1.5 to 5.7) dB for ANSWERS and 2.3 (5th to 95th centile 1.3 to 4.5) dB for sANSWERS. The difference between methods was evaluated with the Wilcoxon signed-rank test; all pairs of comparisons were significantly different at the p < 0.0001 level.

Survival analyses

The following survival analyses were conducted on data from the 320 subjects who had five or more time points with both VF tests and OCT images available (see Figure 5) (in contrast to Reference method: Guided Progression Analysis, which describes the analysis of patients with only OCT data at the baseline visit). Survival was determined at the patient level, with the first UKGTS-eligible eye showing progression classifying the patient as ‘progressed’.

Guided Progression Analysis (reference test)

The survival analysis employing the GPA progression criterion applied in the UKGTS is shown in Figure 21. The number of participants with incident progression was 42 (27.1%) in the placebo group, 29 (17.6%) in the latanoprost group and 71 (22.2%, 95% CI 17.8% to 27.2%) overall. The overall mean survival time was 93.6 weeks (95% CI 91.3 to 96.0 weeks). The HR was 0.566 (95% CI 0.354 to 0.903) and the log-rank test to compare the survival curves was significant at p = 0.016.

FIGURE 21.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the GPA criterion for progression.

ANSWERS

The survival analysis employing the ANSWERS survival criterion giving a 5% false-positive rate when applied sequentially to each VF in a participant’s VF series in the RAPID data set is shown in Figure 22. The number of participants with incident progression was 90 (58.0%) in the placebo group, 82 (49.7%) in the latanoprost group and 172 (53.8%, 95% CI 48.1% to 59.3%) overall. The proportion identified as progressing was statistically significantly greater than that identified with the GPA criterion (McNemar test, p < 0.0001). The overall mean survival time was 72.0 weeks (95% CI 68.6 to 75.5 weeks). The HR was 0.758 (95% CI 0.561 to 1.023) and the log-rank test to compare the survival curves was borderline statistically significant at p = 0.065.

FIGURE 22.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the ANSWERS criterion for progression.

Permutation analyses of pointwise linear regression

The survival analysis employing the PoPLR survival criterion, giving a 5% false-positive rate when applied sequentially to each VF in a participant’s VF series in the RAPID data set, is shown in Figure 23. The number of participants with incident progression was 76 (49.0%) in the placebo group, 57 (34.6%) in the latanoprost group and 133 (41.6%, 95% CI 36.1% to 47.2%) overall. The proportion identified as progressing was statistically significantly greater than that identified with the GPA criterion (McNemar test, p < 0.0001), but statistically significantly less than that identified with the ANSWERS criterion (McNemar test, p < 0.0001). The overall mean survival time was 82.5 weeks (95% CI 79.5 to 85.5 weeks). The HR was 0.590 (95% CI 0.419 to 0.831) and the log-rank test to compare the survival curves was significant at p = 0.002.

FIGURE 23.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the PoPLR criterion for progression.

Structure-guided ANSWERS

The survival analysis employing the sANSWERS survival criterion, giving a 5% false-positive rate when applied sequentially to each VF in a participant’s VF series in the RAPID data set, is shown in Figure 24. The number of participants with incident progression was 103 (66.5%) in the placebo group, 93 (56.4%) in the latanoprost group and 196 (61.3%, 95% CI 55.7% to 66.6%) overall. The proportion identified as progressing was statistically significantly greater than that identified with the GPA and PoPLR criteria (McNemar test, p < 0.0001) and statistically significantly greater than that identified with the ANSWERS criterion (McNemar test, p = 0.042). The overall mean survival time was 69.1 weeks (95% CI 65.7 to 72.4 weeks). The HR was 0.704 (95% CI 0.531 to 0.933) and the log-rank test to compare the survival curves was significant at p = 0.012.

FIGURE 24.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the sANSWERS criterion for progression.

Post hoc survival analyses

Application of the progression criterion for ANSWERS and sANSWERS resulted in the identification of 53.8% and 61.3%, respectively, of UKGTS participants as progressing, substantially more than the 22.2% identified as progressing with the reference standard GPA criterion. The HRs resulting from application of the ANSWERS and sANSWERS criteria (0.758 and 0.704, respectively) suggested a smaller treatment effect than the HR resulting from application of the GPA criterion (0.566). A possible explanation for this finding is that the treatment is more effective in more rapidly progressing eyes. To test this hypothesis, in a post hoc analysis, to label an eye as progressing, an additional ‘rate of change’ in MD criterion was added to the requirement for a negative slope, with statistical significance sufficient for 95% specificity. For sANSWERS, a rate of change cut-off value at approximately the 50th centile of significantly progressing participants was selected; the rate of change in MD at this cut-off value was –0.35 dB per year. Residual life expectancy at diagnosis for glaucoma patients is a median of about 16 years. 106 For a patient presenting with relatively early VF loss in the better eye of MD –4 dB, –0.35 dB per year would result in a MD of nearly –10 dB; this change in MD would alter the probability of failing the vision requirements for driving from about 24% to about 70%. 107 An MD change of –0.35 dB per year is therefore clinically meaningful.

The rates of change estimated by sANSWERS were more accurate than those estimated by SLR or ANSWERS (see Prediction of future visual field state). Therefore, for comparison with sANSWERS, we selected a rate of change criterion for ANSWERS that identified the same number of subjects as progressing as the sANSWERS criterion; this was a MD change of –1.05 dB per year. Survival analyses for ANSWERS and sANSWERS with the additional rate criterion were calculated.

ANSWERS

The survival analysis employing the ANSWERS survival criterion including the rate of change is shown in Figure 25. The number of participants with incident progression was 60 (38.7%) in the placebo group, 44 (26.7%) in the latanoprost group and 104 (32.5%, 95% CI 27.5% to 38.0%) overall. The proportion identified as progressing was statistically significantly greater than that identified with the GPA criterion (McNemar test, p < 0.005). The overall mean survival time was 85.0 weeks (95% CI 81.2 to 88.1 weeks). The HR was 0.593 (95% CI 0.403 to 0.873) and the log-rank test to compare the survival curves was borderline statistically significant at p = 0.0075.

FIGURE 25.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the criterion for progression for ANSWERS including a MD rate of progression of > –1.05 dB per year.

Structure-guided ANSWERS

The survival analysis employing the sANSWERS survival criterion including the rate of change is shown in Figure 26. The number of participants with incident progression was 63 (40.7%) in the placebo group, 41 (24.9%) in the latanoprost group and 104 (32.5%, 95% CI 27.5% to 38.0%) overall. The proportion identified as progressing was statistically significantly greater than that identified with the GPA criterion (McNemar test, p < 0.005). The overall mean survival time was 83.7 weeks (95% CI 80.3 to 87.0 weeks). The HR was 0.573 (95% CI 0.390 to 0.843) and the log-rank test to compare the survival curves was borderline statistically significant at p = 0.0047.

FIGURE 26.

Kaplan–Meier survival curves for the subset of UKGTS participants with OCT images applying the criterion for progression for sANSWERS including a MD rate of progression of > –0.35 dB per year.

Sample size calculations

The purpose of the sample size calculations was to estimate, for each analysis method and the HR estimated for the analysis method (see Survival analyses), the sample size required for a clinical trial of a treatment with the same effect size as latanoprost in the UKGTS and for an observation period of 18 months per participant. A second set of sample size calculations was carried out to assess a HR of 0.6, based on the GPA [see Guided Progression Analysis (reference test)] and PoPLR (section 7.4.3.3) criteria and the ANSWERS and sANSWERS post hoc criteria (see Post hoc survival analyses).

The sample size calculations were based on the survival curves in the preceding section. The number of events required was determined from the observed HR for each method and a HR of 0.60, a type 1 error of 5%, a type 2 error of 10% and equal allocation between treatment groups. In the UKGTS, progression (deterioration) events were observed from 4.5 months onwards (once sufficient data had been collected to enable application of the GPA criterion).The sample size was estimated for an 18-month (4.5 plus 13.5 month) trial from the number of events required, the observed event rate in the placebo group for each method and the censoring rate observed in the UKGTS. Event rates were calculated over the 13.5-month interval, after the initial 4.5 months, for an 18-month trial. The loss to follow-up (censoring) rate was about 10% per year in the UKGTS.

The sample size calculations were made with an online calculator. 108,109

Analyses

The analysis methods applied were as follows:

• PERM:

  • with the timescale as either visit number or time since baseline

  • using individual VF location sensitivity or VF region mean sensitivity

  • when more than one VF test was taken at a visit, taking either the mean values or the values in the initial test at the visit

• MaHMIC

• MaGIC.

All analyses were performed using Stata^® software (version 14; StataCorp LP, College Station, TX, USA) or R.

Permutation test

Participants with three or fewer visits were not included in the permutation analyses as there are only six permutations of three objects and this is not sufficient to produce significance at the 5% level (there is a lack of power – even if the original ordering of the visits is such that the slope is the most extreme that it could be, the one-sided p-value would still be only 1/6 = 0.167). With four visits, if the observed data give a test statistic that is more extreme than those from all other permutations, the one-sided p-value is 1/24 = 0.042, which is just statistically significant using a conventional cut-off value of 5% (although not if corrections for multiple testing are carried out). For four visits, all 23 permutations plus the original ordering of visits were used; for five visits, all 119 permutations plus the original were used; for six visits, all 719 permutations plus the original were used; and, for seven or more visits, 719 permutations were selected at random for use in the PERM in addition to the original ordering of visits.

Visual field sensitivity at each location (or region mean) and the average RNFL thickness from OCT were regressed against time (visit number or time from baseline) using interval-censored regression with VF sensitivities censored at 15 dB (note that imaging outcomes were not censored). Locations with one or fewer uncensored observations were not included in the regression model; the rationale for this is that these locations contain very little information on whether the patient is progressing or not, but will add uncertainty to the model estimates. When more than one VF test was taken at the same visit, either the mean value at each location was used (see RAPID data set, Permutation with visit number as the timescale and Permutation with time since baseline as the timescale, and United Kingdom Glaucoma Treatment Study, Permutation with visit number as the timescale and Permutation with time since baseline as the timescale) or only the first VF of the visit was used (see RAPID data set, Permutations using one visual field per visit, and United Kingdom Glaucoma Treatment Study, Permutations using one visual field per visit).

The regression model included both individual means and slopes over time for each VF location and for the imaging outcome. The model also allowed heterogeneity of the residual variance according to measurement type (separate variances were estimated for the VF and OCT measurements). The regression model was applied to the original ordering of the visits and to each of the (up to 719) permutations described above.

After each model was run, the following parameters were stored:

• Simple average of the slopes from the individual slopes from the VF locations. Note that this is the average of 52 slopes for people with no locations that have one or fewer uncensored observations (i.e. two or more values of ≥ 15 dB at that location). For people with locations that have one or fewer uncensored values, the average will be taken across fewer than 52 slopes.

• Chi-squared test statistic for the simple average of the VF slopes.

• Slope of the imaging outcome.

• Chi-squared test statistic for the imaging outcome slope.

• Chi-squared test statistic for a joint test of all of the individual VF slopes.

• Chi-squared test statistic for a joint test of all of the individual VF slopes and the imaging outcome slope. This test statistic is used in two different ways:

  • requiring the mean VF slope to be negative

  • requiring both the mean VF slope and the imaging slope to be negative.

• Chi-squared test statistic for a joint test of the simple average of the VF slopes and the imaging slope. This test statistic is used in two different ways:

  • requiring the mean VF slope to be negative

  • requiring both the mean VF slope and the imaging slope to be negative.

• Mean slopes within the six sectors defined in the structure/function map (see Figure 4). 22

PERM was then applied to each of these parameters in turn to give a p-value. Progressing eyes were then defined in the following way:

• For the simple average of VF slopes, the imaging slope and the test statistics for both, a one-sided p-value at the 5% significance level was used so that all progressors have a negative slope over time. A one-sided test was used, as clinically a decision to intensify treatment will be made on the basis of deteriorating VF and imaging values, not improvement.

• For the other test statistics, a two-sided p-value was used as it is not clear from a chi-squared test statistic which tail the statistic belongs to. However, to declare progression a negative mean VF slope or both a negative mean VF slope and a negative imaging slope was required. The significance level was set such that a type 1 error of close to 5% was achieved in the RAPID data.

• The mean from each sector was permuted separately and then progressing eyes were defined as those with at least one sector-specific p-value below a critical value. A Bonferroni correction was used to account for multiple testing (described in Mean slopes within visual field sector). Given that such a correction is likely to be overly conservative, other corrections were also considered.

In computing hit rates and specificities using data that accumulate over time, there are two possible approaches. One is to compute the hit rate and specificity at a single visit, using data from that visit and earlier ones to assess trends. The second is to compute the hit rate and specificity cumulatively so that they represent the probability of a positive result at the current and previous visits combined.

There are advantages and disadvantages of both approaches. Using the former approach it is possible to constrain the specificity to take a particular value (such as 95%) at each visit, whereas using the latter approach the specificity will inevitably increase with an increasing number of visits. The latter approach might be preferable if one was considering a particular trial design in which every patient has a fixed maximum number of (say) six visits, whereas the former is arguably preferable for continuous monitoring when there is no upper limit to the number of visits. In this section we have opted for the former approach because we are focusing on the context of a single patient being seen repeatedly in clinical practice.

RAPID data set

In the following sections we consider whether the number of eyes identified as progressing is consistent with a false-positive error (type 1 error) of 5%. This is equivalent to considering whether the specificity, or true negative rate, is consistent with 95%.

United Kingdom Glaucoma Treatment Study data set

Multiple imputation

Multiple imputation models applied to United Kingdom Glaucoma Treatment Study data

Based on the results of the simulation study, we started by applying a multiple imputation model with three spatial neighbours at the same time point, within the same VF region,22 and the nearest time neighbours at the same location as predictors in the model. Instead of adding an extra level of hierarchy into an already complex mixed-effects model, which may have led to problems with convergence in the mixed-effects model, we chose to average across repeated measures recorded at the same visit. However, we found that using such a method resulted in lack of convergence for the multiple imputation process. This behaviour can be seen in the chain plots in Figure 27. This figure shows the burn-in chains for the imputed variables at the first time point for the first 100 iterations, with one plot for each VF location. To produce such a plot, the mean and SD of the imputed values for each variable are calculated after each iteration during the burn-in period. If the multiple imputation process were converging, then one would expect an initial move away from any effect of the starting values and for the chains to then not exhibit a trend as the number of iterations increases. However, there are clear trends observable at some locations in Figure 27, with means tending to decrease over time and SDs correspondingly increasing.

FIGURE 27.

Chains for the means (a) and SDs (b) of the imputed VF variables at time point 1, using the average across VF repeats.

To improve the convergence of these chains, we therefore moved to using models in which a single VF repeat is used, instead of the average across VF repeats. An example of the improved convergence for such a model is shown in Figure 28. This model also uses three spatial neighbours at the same time point and the nearest time neighbours at the same location as predictors. One possible explanation of the improved convergence behaviour in the model using a single VF scan is that it makes full use of the spatial correlations to predict the imputed values. In contrast, by averaging across VF repeats, some of that local correlation structure is being lost or ‘smoothed out’.

FIGURE 28.

Chains for the means (a) and SDs (b) of the imputed VF variables at time point 1, using a single VF repeat.

In addition, we also considered a model that uses a three-stage approach to imputing the data. The data that are imputed in this data set are a mixture of values that are censored (i.e. that lie below 15 dB) and values that are missing (i.e. visits that do not appear in the data set, e.g. because a patient had left the study or did not attend the visit or because the visit did not have both a VF and an OCT measurement). In the first model discussed, we imputed both the censored and the missing data in the same process. In the three-stage approach, we first imputed the censored values only using a constant model (i.e. no predictors were used). This step is performed to obtain some values that are approximately normally distributed before moving on to the second stage, which is to impute the missing visit values using the imputed values from the first stage. SLR models are used for this second step, with three spatial neighbours in the same sector and the nearest time neighbours as predictors. In the third stage, all values that are < 15 dB were imputed using censored regressions, again with three spatial neighbours in the same sector and the nearest time neighbours as predictors.

We used a burn-in period of 50 iterations for the first multiple imputation model and a burn-in of 20 iterations for the third stage of the second model. Ten imputations for each model were produced. Burn-in chains for the imputed data sets at time point 1 are provided in Figures 29 and 30.

FIGURE 29.

Chains for the means (a) and SDs (b) over a burn-in of 50 iterations for the 10 imputations produced by the first model that uses a single VF repeat, three spatial neighbours within a sector and the nearest time neighbours as predictors.

FIGURE 30.

Chains for the means (a) and SDs (b) over a burn-in of 20 iterations for the 10 imputations produced by the second model that uses a single VF repeat and a three-stage approach. The burn-in chains are for the third stage.

Kronecker hierarchical models

The Kronecker hierarchical models, which form part of the MaHMIC method, were performed in R software.

MaHMIC model applied to United Kingdom Glaucoma Treatment Study visual field and imaging data

The results for the fixed effects from the model outlined in Chapter 4 (see Index methods: newly developed, Kronecker model) are provided in Table 22 for the first set of 10 imputations described in Chapter 4 (see Index methods: newly developed, Multiple imputation). These imputations are from the one-stage imputation model that uses the nearest time neighbour and three spatial neighbours in the same sector as predictors in the model. The results from the individual imputations were combined using Rubin’s rules. 102

The slope over time for the imaging outcome was estimated as negative and was statistically significant at the 5% level. However, the slope over time for the VF outcomes was not statistically significant and the point estimate was actually positive. This could be because of learning effects (patients can become more adept at taking the VF test over time), but it could also be because of an inability of the multiple imputation process to account for the selectively missing data. In the original UKGTS, patients were monitored for progression on the basis of their VF values and patients with confirmed progression were removed from the study. Subsequent data are therefore not observed. In a cohort with positive and negative slopes, removing subjects with negative slopes over time will cause a shift to a more positive average for the slopes over time. In a missing at random (MAR) setting, in which all information needed to predict the missing values is present in the observed data set, and with a correctly specified model, it would be possible to adjust for this artefactual positive slope. However, in this setting, with the restricted data set, it is possible that there is a missing not at random (MNAR) scenario. Alternatively, a MAR mechanism may be involved but the proposed multiple imputation model is not using all of the relevant information.

Neither of the time-by-treatment interactions is statistically significant, although both are in the direction of benefit for the trial drug, latanoprost. The VF time-by-treatment interaction also has a lower confidence limit that is not far from zero. The joint test of the time-by-treatment interactions for VF and imaging is not statistically significant (p = 0.24).

The results for the second set of 10 imputations are provided in Table 23. These imputations are from the multiple imputation described in Chapter 4 (see Index methods: newly developed, Multiple imputation), a three-stage approach that first imputes the censored observations from a constant model, then imputes the missing visits using regression models with the nearest time neighbour and three spatial locations within the same sector and finally re-imputes all censored observations using the same predictors.

Using these imputations, the slope over time for VF is again positive. In this case the effect appears more extreme, with a statistically significant positive slope. Again, both time-by-treatment interactions directionally favour the latanoprost group, but neither is statistically significant. In these imputations neither treatment effect is close to being statistically significant. The joint test of the time-by-treatment interactions for VF and imaging is again not statistically significant (p = 0.54).

MaHMIC model applied to United Kingdom Glaucoma Treatment Study visual field data only

In this section we provide the results for the MaHMIC model applied to data sets containing only the VF data, without the addition of imaging data. The model for VF data only is outlined in Chapter 4 (see Index methods: newly developed, Kronecker model).

The results for the first set of imputations are provided in Table 24 and for the second set of imputations are provided in Table 25. Both sets of results have been combined using Rubin’s rules.

The results obtained from the VF-only models are very similar to those obtained from the joint VF and imaging models.

MaGIC model applied to United Kingdom Glaucoma Treatment Study visual field and imaging data

In this section we provide the results for the MaGIC model, described in Chapter 4 (see Index methods: newly developed, Generalised estimating equations), applied to the imputed data sets using Stata software. As for the MaHMIC models, the results from each imputation were combined using Rubin’s rules. The results from the first and second sets of imputations are provided in Tables 26 and 27 respectively.

The parameter estimates obtained from the MaGIC models were quite similar to those obtained from the MaHMIC models, although the CIs were generally slightly wider. There was a slight discrepancy in the imaging time-by-treatment interaction estimates, with the MaGIC model giving small negative estimates and the MaHMIC model producing positive estimates. However, given the wide CIs for the MaGIC estimates, these estimates are still broadly in agreement. In both models, the imaging time-by-treatment interaction was not statistically significant. Similarly, the VF time-by-treatment interaction estimated in the second set of imputations differed in sign to that estimated by the MaHMIC model, although in both models this term was not statistically significant.

As for the MaHMIC models, the joint test of the time-by-treatment interaction for VF and imaging was not statistically significant in either set of imputations (p = 0.76 for the first set of imputations and p = 0.88 for second set of imputations).

We also applied the first imputation model to the two treatment groups separately and combined the imputed data sets for analysis (as opposed to applying the imputation model to both treatment groups combined). When analysing the imputations from the two treatment groups with the MaGIC model, substantively similar results to those presented in Table 26 were obtained.

MaGIC model applied to United Kingdom Glaucoma Treatment Study visual field data only

In this section we present the results of the MaGIC model applied to data sets containing only the VF data, without the addition of imaging data. The model for VF data only is outlined in Chapter 4 (see Index methods: newly developed, Generalised estimating equations). The results are presented in Tables 28 and 29 for the two sets of imputations respectively.

As for the joint model, these results broadly agreed with the estimates produced by the MaHMIC models. The VF-only model also provided similar results to the joint VF and imaging models.