Modeling normative kinetic perimetry isopters using mixed-effects quantile regression

Abstract

Kinetic perimetry is used to quantify visual field size/sensitivity. Clinically, perimetry can be used to diagnose and monitor ophthalmic and neuro-ophthalmic disease. Normative data are integral to the interpretation of these findings. However, there are few computational developments that allow clinicians to collect and analyze normative data from kinetic perimeters. In this article we describe an approach to fitting kinetic responses using linear quantile mixed models. Analogously to traditional linear mixed-effects models for the mean, linear quantile mixed models account for repeated measurements taken from the same individual, but differently from linear mixed-effects models, they are more flexible as they require weaker distributional assumptions and allow for quantile-specific inference. Our approach improves on parametric alternatives based on normal assumptions. We introduce the R package kineticF, a freely available and open-access resource for the analysis of perimetry data. Our proposed approach can be used to analyze normative data from further studies.

Introduction

Sensitivity to light reduces with increasing eccentricity from the center of the visual field. Kinetic perimetry assesses the location at which moving light of a fixed size/intensity can be seen and is used to quantify visual field sensitivity. The location of each test point is recorded as a polar coordinate, and points are joined to form an isopter, denoting a line of differential light sensitivity, within which light of a particular size/intensity can be perceived. 

Constructing robust normative perimetry standards provides an evidence base for interpretation of clinical/research findings—aiding detection of visual field defects. Analysis of normative kinetic perimetry data requires an understanding of anatomical and physiological characteristics, and the correct application of appropriate statistical methods. However, the analytical methods currently employed in the literature do not take into account the repeated-measures structure of a kinetic isopter, nor do they adequately address the lack of normality in the empirical distribution of data points. 

Niederhauser and Mojon (2002) summarized kinetic isopters at each meridian with the mean and 95% confidence intervals. Similar methods have been used by other researchers (Egge, 1984; Wilscher, Wabbels, and Lorenz, 2010) though they assume symmetric distributions along test meridians. There are two issues with this approach: It can define poorly fitted regions of uncertainty that may fall outside the edge of the perimeter test surface or cause crossing of confidence bands; also, as it is not model-based, it cannot be used for inferential purposes (e.g., to compare goodness-of-fit or to determine whether the resulting normative curves and confidence regions need further adjustment by relevant covariates, such as age). Other techniques for summarizing normative fields, such as total isopter area (Bjerre, Codina, & Griffiths, 2014; Quinn, Fea, & Minguini, 1991) and visual field extent along meridians (Wilson, Quinn, Dobson, & Breton, 1991) do not provide reference values that are easily interpretable in a clinical setting. 

Isopter data may show conditional distributions that do not conform to normality assumptions, such as symmetry and constant variance (homoscedasticity), either on the observed or the transformed (e.g., logarithmic) scale of the outcome. For these reasons, modeling based on mean (normal) regression can lead to incorrect inference. Even when approximate normality is achieved after transformation, back-transformation of conditional expectations is troublesome as it may lead to estimates and/or confidence regions outside the admissible range of the outcome. Moreover, isopter data that are collected repeatedly on the same subject are correlated by design. While mixed-effects models for the mean account for the clustered design, they are still subject to strong distributional assumptions and back-transformation issues. 

As an alternative to mean regression, we consider quantile regression (QR; Koenker, 2005), which introduces weak assumptions on the distribution of the error and therefore is robust to deviations from normality. For this reason, no transformation is required—note, however, that even when transformations are introduced to achieve linearity, back-transforming quantiles is a simple task (see for example, Geraci & Jones, 2015). More specifically, in this article we address the analytical challenges of modeling isopter data using mixed-effects models for conditional quantiles called linear quantile mixed models (LQMMs; Geraci & Bottai, 2014). The inclusion of subject-specific effects in mixed-effects models allows for within-subjects correlation resulting from repeated measurements. This approach yields clinically appropriate estimates within expected clinical ranges and correct inference from kinetic perimetry data. 

Here we report the application of LQMMs to normative kinetic isopter data obtained from healthy children (Patel, Cumberland, Walters, Russell-Eggitt, Cortina-Borja et al., 2015) using the R (The R project for Statistical Computing, version 3.1.2; R Core Team, 2016) package kineticF (Patel & Cortina-Borja, 2015). The kineticF package is a collection of functions for cleaning, processing, visualization, and analysis of manual (Goldmann) and automated (Octopus 900) kinetic visual field data, which depends on the package lqmm (Geraci, 2014) used to fit QR models with random effects. 

Methods

Data

The data used in the examples here were collected at Moorfields Eye Hospital, London, as part of a wider program of research, which has been described elsewhere (Patel, Cumberland, Walters, Russell-Eggitt, & Rahi, 2015). Briefly, we recruited children without an ophthalmic condition that affects the visual field to generate normative data and explore visual field development in childhood. Informed written consent was obtained from the children's parents/guardians. This study was approved by the National Health Service Research Ethics Committee for London–Bloomsbury and conforms to the tenets of the Declaration of Helsinki. 

Children aged 5 to 15 years were examined by one clinician (DEP) under clinical conditions using Goldmann and Octopus kinetic perimeters (Haag-Streit, Bern, Switzerland). Manual Goldmann data were scanned and digitized using Engauge digitizer software (open-source, www.digitizer.sourceforge.net). Octopus data were extracted from the perimeter and processed using the R package kineticF (Patel & Cortina-Borja, 2015) available from the Comprehensive R Archive Network repository. 

Statistical methods

Data were collected in polar coordinates (r, θ) defining points along meridians. Due to the cyclical nature of the isopters (Figure 1A), we used harmonic linear predictors in our regression models. For these data, a simple model with sine and cosine terms of periods π and π/2 could be specified as follows:  

 

where r is the isopter value corresponding to meridian θ and the βs are the model's parameters. Of course, Model 1 can include interactions between the sine and cosine terms, or higher frequency harmonics, and be adjusted for confounders such as age and sex. Typically, the error term ε is assumed to be normal, with zero mean and constant variance.  

Figure 1

 

Subject lines (A), points (B), and methods of summarizing normative data; mean and 95% confidence interval (C), mean with 2.5% and 97.5% quantiles (D), and a linear quantile mixed-effects model without (E) and with (F) subject data points.

Since the goal is to produce normative standards (i.e., reference charts), we need to predict specific isopter quantiles. Under normal assumptions, the 100 · pth centile of the conditional distribution of r from Equation 1 is defined as  

 

where β_0,p is an intercept that depends on the quantile p. However, note that the slopes in Equation 2 are constrained by the normal model to be the same for all quantiles. This means that all the individuals in the sample are assumed to have the same distribution, except for a shift β_0,p, which is determined by the normal distribution.  

In QR, individual conditional distributions are allowed to differ from subject to subject by letting all the parameters, including the slopes, vary with p. Thus, the location, scale, and shape of the error term ε are not determined by a theoretical distribution (e.g., normal) but are modeled empirically. 

A QR approach is particularly appropriate to model isopter data. The distribution of points along meridians (as observed from multiple subjects) is skewed and leptokurtotic (Figure 1B). Rather than use a transformation (e.g., logarithmic) of these values to achieve normality, we used QR models, which do not require specific distributional assumptions on the error (distribution-free), adequately describe skewed and kurtotic distributions, and operate on the data's original scale. In these models, we define the regression parameters in relation to a set of quantile levels, namely 2.5%, 50%, and 97.5%, which describe the central tendency and the central 95% of the conditional distribution. Moreover, since by design the isopter data are clustered within individuals, we used LQMMs as proposed by Geraci and Bottai (2014) to account for the within-subjects correlation. For our data, we considered several specifications of LQMMs. 

In our model selection strategy, we first used a forward stepwise approach in which terms were included if their parameter estimates reached a 10% significance level. We only considered random effects terms on the intercept of each model given that, for all models fitted, adding other random effects in the models did not improve the trade-off between goodness-of-fit and number of parameters as measured with the Bayesian Information Criterion (Kuha, 2004). 

The final LQMM  

 

included a random intercept (b_0,p) and a fixed effect for the interaction between cosine and sine terms (β_5,p), along with fixed effects for individual cosine and sine terms. Note that all regression coefficients depend on the particular centile p. 

The R package kineticF makes use of the routines from the package lqmm (Geraci, 2014) to fit QR models with random effects. The output from lqmm is subsequently taken by kineticF to provide graphical and statistical summaries matching those produced by common perimeters, thus immediately facilitating clinical interpretation. Cubic smoothing spline models were used to smooth the predicted isopters between consecutive meridians by interpolating the fitted quantile values in each meridian. The smoothing parameter of these spline models was chosen using a cross-validation procedure (Green & Silverman, 1994). This interpolating process is included in kineticF and produces smooth bands for specific quantiles of the isopters. These predicted normative bands avoid assuming a constant standard deviation and symmetry across meridians. 

These outputs were then formatted as templates for use in both clinical and research settings (available from http://e-lucid.com/i/video_and_images/optic_templates.html). All data management and statistical modeling was performed using Microsoft Excel 2010® and R. 

Results

Data were collected from 154 children. The examples shown here are from 30 subjects aged 8 to 11 years, plotting isopter I4e. 

Figure 1A displays the trajectories of individuals, demonstrating the need for modeling the within-subject structure present in the dataset. (Note that the lack of data points along the 0° meridian corresponds to the void [hashed] area within the Goldmann perimeter bowl—that is, an area where no points can be plotted.) Figure 1B shows the sample data points along 24 meridians located at 15° intervals; the asymmetry of these observations along each meridian is apparent. For example, along the 150° meridian, there is a concentration of points at 50° eccentricity (near the limit of the nasal field), with a long lower tail extending to approximately 35°. Figure 1C and D show these data's mean and 95% confidence interval, and the 2.5th, 50th, and 97.5th empirical centiles along each meridian, which are methods that have been used previously in the literature. Finally, Figure 1E and F show the fitted mixed-effects QR model as defined in Equation 3 with predicted and smoothed curves for the 2.5%, 50%, and 97.5% quantiles, without and with observed data points, respectively. 

Discussion

We have applied robust statistical methods to the modeling of normative kinetic perimetry data sampled from multiple subjects and developed a collection of visual and analytic software routines to carry out similar analyses. Our approach is robust against deviations from normality and homoscedasticity, and represents a substantial advance to current practice in that it lays the basis for statistical model comparisons based on goodness-of-fit criteria, and improves clinical interpretation by providing smooth quantile curves representing the variability of normative visual fields. Since these two aspects have not been fully addressed in the literature, we cannot compare our approach with other modeling strategies (Vonthein et al., 2007). However, Figure 1F shows a substantial improvement with respect to other approaches that do not use regression models, namely Figure 1C and D. Both are flawed: The former assumes that the mean and a constant dispersion measure are enough to model the conditional distribution; the latter estimates the quantiles of interest for each isopter without reference to the whole data set or to the within-subjects correlation, and does not produce smooth normative bands. 

A potential limitation of methods based on QR is that small samples (data sparsity) may lead to erratic estimates of tail quantiles and, in the worst case, to quantile crossing. This can be avoided by ensuring a balanced design and an adequately sized sample. 

As noted, we used a cross-validation criterion to determine the smoothness of the cubic splines that interpolated the quantiles of the isopters predicted by our model. This could have been achieved in a Bayesian setting assuming a prior distribution for the smoothing parameters and estimating them with, for instance, the median of its posterior distribution. 

Our approach can be applied to clinical data from control subjects thus providing normative data against which observations from case subjects can be compared. Without defining these normative standards, it is impossible to accurately assess disease status. The mixed-effects quantile models described in this paper can be readily applied to perimetric data using packages lqmm and kineticF, which are freely available from public repositories. Our findings support the use of the models reported here when generating normative kinetic perimetry standards. 

Acknowledgments

We would like to acknowledge the Study of Optimal Perimetric Testing In Children (OPTIC) group members: Peng Tee Khaw, Bronwen Walters, Phillippa Cumberland, Isabelle Russell-Eggitt, Chris Timms, John Brookes, Anthony Moore, Maria Papadopoulos, David Garway-Heath, Ananth Viswanathan, Alki Liasis, David Crabb, Mario Cortina-Borja, Dipesh Patel, and Jugnoo Rahi. The OPTIC study was supported by The Guide Dogs for the Blind (Grant: OR2009-04e, http://www.guidedogs.org.uk/). 

Commercial relationships: None. 

Corresponding author: Mario Cortina-Borja. 

Email: [e-mail]. 

Address: Clinical Epidemiology, Nutrition and Biostatistics Section, UCL Institute of Child Health, London, UK.