Collective endpoint of visual acuity and contrast sensitivity function from hierarchical Bayesian joint modeling

Abstract

Clinical trials typically analyze multiple endpoints for signals of efficacy. To improve signal detection for treatment effects using the high-dimensional data collected in trials, we developed a hierarchical Bayesian joint model (HBJM) to compute a five-dimensional collective endpoint (CE_5D) of contrast sensitivity function (CSF) and visual acuity (VA). The HBJM analyzes row-by-row CSF and VA data across multiple conditions, and describes visual functions across a hierarchy of population, individuals, and tests. It generates joint posterior distributions of CE_5D that combines CSF (peak gain, peak frequency, and bandwidth) and VA (threshold and range) parameters. The HBJM was applied to an existing dataset of 14 eyes, each tested with the quantitative VA and quantitative CSF procedures in four Bangerter foil conditions. The HBJM recovered strong correlations among CE_5D components at all levels. With 15 qVA and 25 qCSF rows, it reduced the variance of the estimated components by 72% on average. Combining signals from VA and CSF and reducing noises, CE_5D exhibited significantly higher sensitivity and accuracy in discriminating performance differences between foil conditions at both the group and test levels than the original tests. The HBJM extracts valuable information about covariance of CSF and VA parameters, improves precision of the estimated parameters, and increases the statistical power in detecting vision changes. By combining signals and reducing noise from multiple tests for detecting vision changes, the HBJM framework exhibits potential to increase statistical power for combining multi-modality data in ophthalmic trials.

Introduction

As research studies that test new methods of prevention, detection, treatment, and care in healthcare, clinical trials often consist of five phases, early phase I and phases I–IV, with each phase designed to test a different aspect of the new drug, device, or method (Rosenberger & Haines, 2002; Thall, 2008; Torres-Saavedra & Winter, 2022; Umscheid, Margolis, & Grossman, 2011). Although early phase I trials test the initial safety and potency of the new treatment, phase I trials aim to determine how effective the new drug, treatment, or device is and how the human body processes it, phase II trials investigate if the new treatment leads to prolonged life or better health outcomes, phase III trials compare the new treatment to existing standard treatment, and phase IV trials evaluate the new treatment after approval from the Food and Drug Administration. Although clinical trials often involve hundreds or even thousands of subjects, the sample size in most early phase I and phase I trials is about 10 to 20 (Cook, Hansen, Siu, & Abdul, 2015; Rosenberger & Haines, 2002; Storer, 1989; Umscheid et al., 2011). 

Endpoints are measured in clinical trials to evaluate the efficacy of a new treatment. An endpoint is an event or outcome that can be measured objectively to determine whether the treatment being studied is beneficial. In ophthalmic trials, the endpoint could be a clinical event (e.g., retinal hemorrhage, stroke), patient symptom (e.g., pain), visual function (e.g., visual acuity, contrast sensitivity), retinal structure (e.g., optical coherence tomography, fundus photography), or an equivalent surrogate (McLeod et al., 2019; US Department of Health and Human Services et al., 1998; US Department of Health and Human Services et al., 2022). To obtain multi-faceted assessment of treatment outcomes, most clinical trials collect multiple endpoints and analyze them independently (but see composite endpoint; Kramer et al., 2006; McLeod et al., 2019; US Department of Health and Human Services et al., 1998; US Department of Health and Human Services et al., 2022). 

In this study, we attempted to address an important issue in early-stage retina trials with multiple functional vision endpoints, in which statistical power for detecting vision gains based on each endpoint is limited by the variance caused by small sample size, noisy measurements, and variable treatment response (Abou-Hanna, Andrews, Khan, Musch, & Jayasundera, 2021; Allingham, Mettu, & Cousins, 2022; Cukras et al., 2018; da Cruz et al., 2018; de Guimaraes, Georgiou, Bainbridge, & Michaelides, 2021; Reichel et al., 2022). We developed a hierarchical Bayesian joint model (HBJM) to exploit the relationship between two endpoints, visual acuity (VA) and contrast sensitivity function (CSF), and computed a five-dimensional collective endpoint (CE_5D) that consists of two parameters of the visual acuity behavioral function (VABF) (Lesmes & Dorr, 2019) and three parameters of contrast sensitivity function (Watson & Ahumada, 2005). CE_5D is considered as a “collective endpoint” because it combines information from two endpoint assessments. We demonstrated that the high dimensional collective endpoint significantly improved the statistical power in a study with a small sample size (N = 14 eyes). 

Visual acuity and contrast sensitivity function tests are designed to assess different aspects of vision. While VA focuses on the spatial resolution of vision that is important for resolving details in tasks such as reading (Whittaker & Lovie-Kitchin, 1993), CSF measures the ability to detect subtle light changes against the background that is important for tasks such as mobility (Lord, 2006) and driving (Owsley & McGwin, 2010). It's also well known that many visual diseases affect CSF and VA differentially, with reports of CSF deficits but normal or near-normal VA in early stages of many clinical conditions (Elliott & Situ, 1998; Hawkins, Szlyk, Ardickas, Alexander, & Wilensky, 2003; Hou et al., 2010; Jia, Zhou, Lu, Lesmes, & Huang, 2015; Lesmes et al., 2016; Midena, Angeli, Blarzino, Valenti, & Segato, 1997; Onal, Yenice, Cakir, & Temel, 2008; Ou, Lesmes, Christie, Denlar, & Csaky, 2021; Ramulu, Dave, & Friedman, 2015; Shandiz et al., 2011; Stellmann, Young, Pottgen, Dorr, & Heesen, 2015; Vingopoulos et al., 2021; Wai et al., 2021; Wilensky & Hawkins, 2001; Zimmern, Campbell, & Wilkinson, 1979). On the other hand, VA and CS both reflect properties of the same visual system for an individual, and there are many reports of significant correlations between VA and CSF in normal older adults and clinical populations (Alexander, Derlacki, & Fishman, 1995; Bellmann, Unnebrink, Rubin, Miller, & Holz, 2003; Brown & Lovie-Kitchin, 1989; Cormack, Tovee, & Ballard, 2000; Elliott & Situ, 1998; Haegerstrom-Portnoy Schneck, & Brabyn, 1999; Kiser, Mladenovich, Eshraghi, Bourdeau, & Dagnelie, 2005; Roark & Stringham, 2019; Rubin, Roche, Prasada-Rao, & Fried, 1994; West et al., 2002). Although one could ultimately construct a unified model in which VA and CSF are determined by both shared and independent parameters that are affected differentially by different diseases and in different disease stages, much more basic and clinical research is necessary. Here, we take an alternative approach based on the covariance joint modeling framework (Palestro et al., 2018; Turner, 2015; Turner et al., 2013; Turner, Van Maanen, & Forstmann, 2015; Turner, Rodriguez, Norcia, McClure, & Steyvers, 2016; Turner, Wang, & Merkle, 2017) in which VA and CSF are modeled separately but with covariances between parameters of the models. 

Most experimental designs follow a hierarchical structure (Kim, Pitt, Lu, Steyvers, & Myung, 2014; Yin, Qin, Sargent, Erlichman, & Shi, 2018). At the top, there is the study population, which can be divided into different treatment arms. Each arm consists of multiple subjects that are tested with multiple assessments, and each assessment could be repeated multiple times. While many common statistical methods assume either homogeneity (e.g., fixed effects) or complete independence across subjects or tests, hierarchical modeling is a statistical method that has been developed to quantitatively and coherently combine heterogeneous information (Kruschke, 2014; Rouder & Lu, 2005). Constructed by combining sub-models and probability distributions at multiple levels that represent relationships between information arising within and between levels, hierarchical models are often computed with Bayes’ theorem from all the observed data in a study (rather than data from a subject or a test) (Kruschke, 2014; Kruschke & Liddell, 2018). The result is an updated joint posterior distribution of all the hyperparameters and parameters (Gu et al., 2016; Kruschke, 2014; Zhao, Lesmes, Dorr, & Lu, 2021). By sharing information within and across levels via conditional dependencies, hierarchical Bayesian modeling (HBM) can reduce the variance of the estimated posterior distributions through (1) decomposition of variabilities from different sources with parameters and hyperparameters (Song, Behmanesh, Moaveni, & Papadimitriou, 2020), and (2) shrinkage of the estimated parameters at the lower levels toward the modes of the higher levels when there is not sufficient data at the lower level (Kruschke, 2014; Rouder, Sun, Speckman, Lu, & Zhou, 2003; Rouder & Lu, 2005). 

Hierarchical Bayesian models (HBMs) have been developed and widely applied in cognitive science, including signal detection (Rouder et al., 2003; Rouder & Lu, 2005), decision making (Lee, 2006; Merkle, Smithson, & Verkuilen, 2011), and functional magnetic resonance imaging (Ahn, Krawitz, Kim, Busmeyer, & Brown, 2011; Wilson, Cranmer, & Lu, 2020). HBJM has been developed in cognitive neuroscience to investigate brain-behavior relationships (Molloy et al., 2018; Molloy, Bahg, Lu, & Turner, 2019; Palestro et al., 2018). As an alternative to the two-stage correlation approach in cognitive neuroscience, where estimated parameters of a cognitive model on human behavior are simply correlated with a neural measure of interest, joint modeling enforces constraints on the parameters of the cognitive model using the neural data. It has been used to perform inference on brain-behavior relationships, improve the precision of the estimated model parameters, and make predictions about either neural or behavioral data (Cassey, Gaut, Steyvers, & Brown, 2016; Turner, 2015; Turner et al., 2013; Turner et al., 2015; Turner et al., 2016; Turner et al., 2017). Although we have previously developed HBMs to estimate parameters of the VABF and CSF models within each type of measurement and significantly improved the accuracy and precision of the estimates (Zhao, Lesmes, Dorr, & Lu, 2021; Zhao, Lesmes, Hou, & Lu, 2021), we used HBJM to model the relationships between VA and CSF in this study. Our goals were to (1) compute a multi-dimensional collective endpoint (CE_5D) of CSF and VA by modeling row-by-row VA and CSF data simultaneously, and (2) exploit the statistically reciprocal relationships between VA and CSF parameters to combine signals from VA and CSF tests and reduce the variability of the estimated endpoints and therefore improve the power of statistical inference (Figure 1). 

Figure 1.

 

Illustration of a two-dimensional CE. (A) The gray iso-density contours represent two-dimensional distributions of a VA-CSF CE in two conditions. The blue and purple curves represent VA and CSF distributions computed from the BIP. The black curves represent the marginal distributions of CE components on the VA and CSF axes, and the projection of CE on the decision axis. The p values for detecting the change between the two conditions are 0.15 based on VA and CSF from BIP, 0.10 based on VA and CSF from CE, and 0.006 based on CE. (B) Signals from BIP analyses of VA (blue) and CSF (purple), and from CE components and CE (black). (C) Noise from BIP analyses of VA (blue) and CSF (purple), and from CE components and CE (black). (D) Signal-to-noise ratios from BIP analyses of VA (blue) and CSF (purple), and from CE components and CE (black).

We began with generative models of VA and CSF that, with appropriate parameters, can reproduce the row-by-row data from the measurements. We used the VABF with two parameters, VA threshold and VA range, as the generative model of observer performance in VA assessment (Figure 2A) (Lesmes & Dorr, 2019), and a log parabola with three parameters, peak gain (PG), peak spatial frequency (PF) and bandwidth at half-height (BH) (Figure 2C) (Watson & Ahumada, 2005), as the generative model of observer performance in CSF assessment. Both models have previously been used to model and reproduce row-by-row data in their respective measurements (Figures 2B, 2E) (Hou, Lesmes, Bex, Dorr, & Lu, 2015; Hou et al., 2016; Lesmes, Lu, Baek, & Albright, 2010; Zhao, Lesmes, Dorr, Bex, & Lu, 2021; Zhao, Lesmes, Dorr, & Lu, 2021; Zhao, Lesmes, Hou, et al., 2021). We defined CE_5D as a five-dimensional vector consisting of VA threshold, VA range, PG, PF, and BH. 

Figure 2.

 

(A) The VABF with two parameters: threshold \(\theta _{ij,Threshold}^{VA}\;\)and range \(\theta _{ij,Range}^{VA}\). (B) Simulated row-by-row data generated by the VABF model in a VA assessment. (C) The log parabola model of the CSF with three parameters: peak gain \(\theta _{ij,PG}^{CSF}\), peak spatial frequency \(\theta _{ij,PF}^{CSF}\), and bandwidth \(\theta _{ij,BH}^{CSF}\). (D) Psychometric functions at different spatial frequencies in a CSF test. (E) Simulated row-by-row data generated by the CSF model in a CSF assessment.

Next, we developed an HBJM with population, individual, and test levels to model the row-by-row data from both VA and CSF assessments on all individuals in an experiment within a single structure and compute the CE parameter and hyperparameter distributions at all three levels, taking into account of the hierarchical experimental design, covariance of CE components, and conditional dependency across different levels (Kruschke, 2014; Zhao, Lesmes, Dorr, & Lu, 2021; Zhao, Lesmes, Hou, et al., 2021). The new development enables us to compute the joint distributions of the five-dimensional CE_5D from the row-by-row data of human observers in all the tests and make statistical inferences at the population, individual, and test levels. 

We applied the HBJM to an existing experimental dataset (Zhao, Lesmes, Dorr, Bex, et al., 2021; Zhao, Lesmes, Dorr, & Lu, 2021) with 14 eyes, each tested with the quantitative VA (qVA) (Lesmes & Dorr, 2019; Zhao, Lesmes, Dorr, & Lu, 2021) and quantitative CSF (qCSF) (Hou et al., 2015; Lesmes et al., 2010) procedures in four Bangerter foil conditions, and computed CE_5D distribution for each eye in each condition. We quantified the uncertainty, sensitivity, and discriminability of the estimated CE_5D, and compared them with metrics from Bayesian inference procedure (BIP) analyses of the tests, VA_1D and AULCSF_1D. 

Methods

Data

The dataset consisted of 14 eyes of seven subjects with normal or corrected-to-normal vision, each tested with qVA and qCSF in four conditions: no foil (F0) and three levels of Bangerter foils (Ryser Ophtalmologie, St. Gallen, Switzerland) with nominal acuities of 20/25 (F1), 20/30 (F2), and 20/100 (F3). The qVA data were included in previous publications (Zhao, Lesmes, Dorr, Bex, et al., 2021; Zhao, Lesmes, Dorr, & Lu, 2021). The qCSF data have not been published. The test order of the four conditions was randomized across subjects. Each qVA test consisted of 45 rows (M = 45), each with three equal-size high-contrast optotypes (Figure 3A) randomly sampled without replacement from the 10 Sloan letters (C, D, H, K, N, O, R, S, V, and Z). The size of the optotypes was determined by qVA (Lesmes & Dorr, 2019; Zhao, Lesmes, Dorr, Bex, et al., 2021). Each qCSF test consisted of 50 rows (K = 50), each with three equal-size bandpass-filtered optotypes (Figure 3B) randomly sampled with replacement from the ten Sloan letters (C, D, H, K, N, O, R, S, V, and Z). The sizes (and therefore center spatial frequencies) and contrasts of the optotypes were determined by qCSF (Hou et al., 2015; Lesmes et al., 2010). In both tests, the stimuli stayed on the screen until the subject verbally reported the identities of the letters presented on the screen, which were entered into the computer by the experimenter. 

Figure 3.

 

Illustrations of the test stimuli in one row used in the (A) qVA and (B) qCSF.

Apparatus

Both qVA and qCSF were implemented in MATLAB (MathWorks Corp., Natick, MA, USA). A 24-inch Dell P2415Q liquid-crystal display monitor at a 3840 × 2160 pixel resolution and a 55-inch Samsung UN55FH6030 monitor at a 1920 × 1080 pixel resolution were used for the qVA and qCSF tests, respectively. Subjects viewed the displays monocularly at four meters. The untested eye was covered with an opaque patch. A Dell computer with Intel Xeon W-2145 @ 3.70GHz CPU (8 cores and 16 threads) and 64GB installed memory (RAM) was used to conduct all analyses with JAGS (Plummer, 2003) in R (R Core Team, 2003). 

Hierarchical Bayesian joint modeling of VA and CSF assessments

In this section, we first describe the generative models of VA and CSF, and the corresponding BIP to infer the parameters of the generative models for an individual in a single test. We then introduce the three-level HBJM to account for the row-by-row data from both VA and CSF assessments on all individuals within a single structure and compute the joint distributions of the collective endpoint (CE_5D), comprised of parameters of the two generative models, across all tests and individuals and at the population level. Finally, we describe the computational and statistical inference procedures. In the current development, we treat each unique combination of eye and Bangerter foil condition as an individual. 

Generative model of visual acuity and Bayesian inference procedure

With parameters \(\theta _{ij}^{VA} = ( {\theta _{ij,Threshold}^{VA},\theta _{ij,Range}^{VA}} )\), where \(\theta _{ij,Threshold}^{VA}\) is VA threshold corresponding to the dʹ = 2 performance level, and \(\theta _{ij,Range}^{VA}\) is VA range defined as the difference between VA thresholds at dʹ = 1 and dʹ = 4 performance levels, the generative model of visual acuity (Figure 4A) specifies the probability of obtaining r_ijm (=0, 1, 2, or 3) correct responses from individual i tested with optotype size os_ijm in row m of the j^th test (Figure 2A) (Lesmes & Dorr, 2019; Zhao, Lesmes, Dorr, Bex, et al., 2021):  

\begin{eqnarray}p\left( {r_{ijm}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ijm}}} \right) = {\rm{G}}\left( {o{s_{ijm}},\theta _{ij}^{VA}} \right),\quad\end{eqnarray}

(1)

where \({\rm{G}}( {o{s_{ijm}},\theta _{ij}^{VA}} )\) is the probability of the number of correct identifications in a row (Supplementary Material A). The probability of observing the responses of individual i in all M rows in the j^th test is the product of \(p( {r_{ijm}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ijm}}} )\) in all rows:  

\begin{eqnarray} && p\left( {r_{ij,1:M}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ij,1:M}}} \right)\nonumber\\ && \quad = \mathop \prod \limits_{m = 1}^M p\left( {r_{ijm}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ijm}}} \right).\quad\end{eqnarray}

(2)

 

Figure 4.

 

(A) The generative model of visual acuity specifies the VABF of individual i in test j with two parameters: \(\theta _{ij}^{VA} = ( {\theta _{ij,Threshold}^{VA},\theta _{ij,Range}^{VA}} )\), where \(\theta _{ij,Threshold}^{VA}\) and \(\theta _{ij,Range}^{VA}\) are the VA threshold and VA range. (B) The generative model of CSF specifies the CSF of individual i in test j with three parameters\(:\;\theta _{ij}^{CSF} = ( {\theta _{ij,PG}^{CSF},\theta _{ij,PF}^{CSF},\theta _{ij,BH}^{CSF}} )\), where \(\theta _{ij,PG}^{CSF}\;\)is the peak sensitivity, \(\theta _{ij,PF}^{CSF}\) is the peak spatial frequency (cycles/degree), and \(\theta _{ij,BH}^{CSF}\) is the bandwidth (octaves) at half of the peak sensitivity. (C) A three-level HBJM of the VABF and CSF across multiple individuals and tests. At the population level, μ^CE and Σ are the mean and covariance hyperparameters of the population. At the individual level \(\rho _i^{{\rm{CE}}}\), ϕ_CSF, and ϕ_VA are the mean and covariance hyperparameters of individual i. At the test level, \(\theta _{ij}^{CE}\) are the parameters of individual i in test j.

In other words, with known \(\theta _{ij}^{VA}\), the VA generative model allows us to compute the probability of the number of correct optotype identifications in each row as well as the probability of obtaining the entire dataset in a single test. It can be used to generate the row-by-row responses in VA tests (Figure 2B). 

On the flipside, we can infer the posterior distribution of \(\theta _{ij}^{VA}\) from the row-by-row data \(Y_{ij}^{VA} = \{ {( {r_{ij,1:M}^{VA},o{s_{ij,1:M}}} )} \}\;\)of individual i in the j^th test with the BIP (Lesmes et al., 2010; Zhao, Lesmes, Dorr, Bex, et al., 2021; Zhao, Lesmes, Dorr, & Lu, 2021):  

\begin{eqnarray} && p\left( {\theta _{ij}^{VA}{\rm{|}}Y_{ij}^{VA}} \right)\nonumber\\ && \quad = \frac{{\mathop \prod \nolimits_{m = 1}^M p\left( {r_{ijm}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ijm}}} \right){p_0}\left( {\theta _{ij}^{VA}} \right)}}{{\smallint \mathop \prod \nolimits_{m = 1}^M p\left( {r_{ijm}^{VA}{\rm{|}}\theta _{ij}^{VA},o{s_{ijm}}} \right){p_0}\left( {\theta _{ij}^{VA}} \right)d\theta _{ij}^{VA}}},\quad\end{eqnarray}

(3)

where \({p_0}( {\theta _{ij}^{VA}} )\) is the prior probability distribution of \(\theta _{ij}^{VA}\), and the denominator is the integral across all possible values of \(\theta _{ij}^{VA}\). 

Generative model of contrast sensitivity function and Bayesian inference procedure

The generative model of the CSF (Figure 4B) specifies contrast sensitivity \(S( {{f_{ijk}},\theta _{ij}^{CSF}} )\) at spatial frequency f_ijk for individual i in row k of the j^th test using a log parabola function with three parameters \(\theta _{ij}^{CSF} = ( {\theta _{ij,PG}^{CSF},\theta _{ij,PF}^{CSF},\theta _{ij,BH}^{CSF}} )\) (Figure 2C) (Lesmes et al., 2010; Rohaly & Owsley, 1993; Watson & Ahumada, 2005; Zhao, Lesmes, Hou, et al., 2021):  

\begin{eqnarray}&&{\rm{lo}}{{\rm{g}}_{10}}\left( {S\left( {{f_{ijk}},\theta _{ij}^{CSF}} \right)} \right) = \;{\rm{lo}}{{\rm{g}}_{10}}\left( {\theta _{ij,PG}^{CSF}} \right)\nonumber\\ && \quad - \frac{4}{{{\rm{lo}}{{\rm{g}}_{10}}\left( 2 \right)}}{\left( {\frac{{{\rm{lo}}{{\rm{g}}_{10}}\left( {{f_{ijk}}} \right) - {\rm{lo}}{{\rm{g}}_{10}}\left( {\theta _{ij,PF}^{CSF}} \right)}}{{\theta _{ij,BH}^{CSF}}}} \right)^2},\quad\end{eqnarray}

(4)

where \(\theta _{ij,PG}^{CSF}\;\)is the peak gain, \(\theta _{ij,PF}^{CSF}\) is the peak spatial frequency (cycles/degree), and \(\theta _{ij,BH}^{CSF}\) is the bandwidth (octaves) at half of the peak sensitivity. The probability of obtaining a correct response is described with a psychometric function (Figure 2D) (Hou et al., 2015):  

\begin{eqnarray} && p\left( {r_{ijkl}^{CSF} = 1{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}\;} \right) = g + \left( {1 - g - \frac{\lambda }{2}} \right)\nonumber\\ && \quad {\rm{\Phi }}\left( {\frac{{{\rm{lo}}{{\rm{g}}_{10}}\left( {{c_{ijkl}}} \right) + {\rm{lo}}{{\rm{g}}_{10}}\left( {S\left( {{f_{ijk}},\theta _{ij}^{CSF}} \right)} \right)}}{\sigma }} \right),\end{eqnarray}

(5)

where \(r_{ijkl}^{CSF}\)and c_ijkl are the response to and contrast of the l^th (=1, 2, 3) optotype in row k with a total of L optotypes; f_ijk is the spatial frequency for all three optotypes in the row; g is the guessing rate, λ, set to 0.04, is the lapse rate (Hou et al., 2015; Lesmes et al., 2010; Wichmann & Hill, 2001), Φ is the standard cumulative Gaussian function, and σ determines the steepness of the psychometric function. The probability of an incorrect response is:  

\begin{eqnarray}&& p\left( {r_{ijkl}^{CSF} = 0{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}} \right)\nonumber\\ && \quad = 1 - p\left( {r_{ijkl}^{CSF} = 1{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}\;} \right).\quad\end{eqnarray}

(6)

 

The likelihood function, defined by Equations 5 and 6, quantifies the probability of a correct or incorrect response given the stimulus and CSF parameters in a row. The probability of observing the particular responses of individual i in all K rows in the j^th test is the product of \(p( {r_{ijkl}^{CSF}{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}} )\) in all rows:  

\begin{eqnarray} && p\left( {r_{ij,1:K,1:L}^{CSF}{\rm{|}}\theta _{ij}^{CSF},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} \right)\nonumber \\ && \quad = \mathop \prod \limits_{k = 1}^K \mathop \prod \limits_{l = 1}^L p\left( {r_{ijkl}^{CSF}{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}} \right).\quad\end{eqnarray}

(7)

 

In other words, with known \(\theta _{ij}^{CSF}\), the CSF generative model allows us to compute the probability of a correct or incorrect response to each optotype in every row as well as the probability of obtaining the entire dataset in a single test. It can be used to generate the row-by-row responses in CSF tests (Figure 2E). 

On the flipside, we can infer the posterior distributions of \(\theta _{ij}^{CSF}\) from the row-by-row data \(Y_{ij}^{CSF} = \{ {( {r_{ij,1:K,1:L}^{CSF},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} )} \}\;\)of individual i in the j^th test with the BIP (Hou et al., 2015; Hou et al., 2016; Lesmes et al., 2010):  

\begin{eqnarray} && p\left( {\theta _{ij}^{CSF}{\rm{|}}Y_{ij}^{CSF}} \right)\nonumber\\ && \quad = \frac{{\mathop \prod \nolimits_{k = 1}^K \mathop \prod \nolimits_{l = 1}^L p\left( {r_{ijkl}^{CSF}{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}} \right){p_0}\left( {\theta _{ij}^{CSF}} \right)}}{{\smallint \mathop \prod \nolimits_{k = 1}^K \mathop \prod \nolimits_{l = 1}^L p\left( {r_{ijkl}^{CSF}{\rm{|}}\theta _{ij}^{CSF},{f_{ijk}},{c_{ijkl}}} \right){p_0}\left( {\theta _{ij}^{CSF}} \right)d\theta _{ij}^{CSF}}},\quad\end{eqnarray}

(8)

where \({p_0}( {\theta _{ij}^{CSF}} )\) is the prior probability distribution of \(\theta _{ij}^{CSF}\), and the denominator is the integral across all possible values of \(\theta _{ij}^{CSF}\). 

HBJM

Although the BIP can be used to infer the parameters of the generative model of an individual in a single test, the prior distributions in the procedure are usually uninformative, and the posterior parameters distributions are computed independently for each assessment in each test of an individual, without considering potential conditional dependencies and covariance of the parameters across assessments, tests and individuals. Here, we developed a three-level HBJM to model the row-by-row data from multiple endpoint assessments of all individuals with a joint parameter and hyperparameter distribution of the CE_5D, comprised of the parameters of the generative models of all the endpoint assessments. In the current development, the five-dimensional collective endpoint \(\theta _{ij}^{{\rm{CE}}} = ( {\theta _{ij,{\rm{VA}}}^{{\rm{CE}}},\;\theta _{ij,{\rm{CSF}}}^{{\rm{CE}}}} )\;\)for individual i in the j^th test consists of two-dimensional \(\theta _{ij,{\rm{VA}}}^{{\rm{CE}}} = ( {\theta _{ij,Threshold}^{{\rm{CE}}},\theta _{ij,Range}^{{\rm{CE}}}} )\) and three-dimensional \(\theta _{ij,CSF}^{CE} = ( {\theta _{ij,PG}^{CE},\theta _{ij,PF}^{CE},\theta _{ij,BH}^{CE}} )\) components. The CE_5D hyperparameters at the population and individual levels are five-dimensional random variables η^CE and \(\tau _i^{{\rm{CE}}}.\) In the HBJM, the distribution at a higher level sets the prior for the means of the distributions at a lower level, i.e., the distribution of the population hyperparameter η^CE sets the prior for \(\rho _i^{{\rm{CE}}}\), the mean of the distribution of the hyperparameter \(\tau _i^{{\rm{CE}}}\) of individual i, which in turn sets the prior for the mean of the distribution of parameter \(\theta _{ij}^{{\rm{CE}}}\) of individual i in the j^th test. 

The probability distribution p(η^CE) of the CE_5D hyperparameter \({\eta ^{{\rm{CE}}}} = (\eta _{VA}^{{\rm{CE}}},\) \(\eta _{CSF}^{{\rm{CE}}})\) with \(\eta _{VA}^{{\rm{CE}}} = ( {\eta _{Threshold}^{{\rm{CE}}},\eta _{Rang{\rm{e}}}^{{\rm{CE}}}} )\) and \(\eta _{CSF}^{{\rm{CE}}} = ( {\eta _{PG}^{{\rm{CE}}},\eta _{PF}^{{\rm{CE}}},\eta _{BH}^{{\rm{CE}}}} )\) at the population level is modeled as a mixture of five-dimensional Gaussian distributions \(\mathcal{N}( {{\eta ^{{\rm{CE}}}},{\mu ^{{\rm{CE}}}},\Sigma } )\):  

\begin{eqnarray}p\left( {{\eta ^{{\rm{CE}}}}} \right) = \mathcal{N}\left( {{\eta ^{{\rm{CE}}}},{\mu ^{{\rm{CE}}}},\Sigma } \right)p\left( {{\mu ^{{\rm{CE}}}}} \right)p\left( \Sigma \right),\quad\end{eqnarray}

(9)

where μ^CE is the mean of η^CE, with distribution p(μ^CE), and Σ is the covariance of η^CE, with distribution p(Σ). 

The probability distribution \(p( {\tau _i^{{\rm{CE}}}{\rm{|}}{\eta ^{{\rm{CE}}}}} )\;\)of the CE_5D hyperparameter \(\tau _i^{\rm CE} = (\tau _{i,VA}^{\rm CE},\) \(\tau _{i,CSF}^{\rm CE})\) of individual i is modeled as the product of two-dimensional \(\tau _{i,VA}^{\rm CE} = ( {\tau _{i,Threshold}^{\rm CE},\tau _{j,Range}^{\rm CE}} )\) and three-dimensional \(\tau _{i,CSF}^{\rm CE} = ( {\tau _{i,PG}^{\rm CE},\tau _{i,PF}^{\rm CE},\tau _{i,BH}^{\rm CE}} )\) Gaussian distributions:  

\begin{eqnarray}&& p\left( {\tau _i^{{\rm{CE}}}{\rm{|}}{\eta ^{{\rm{CE}}}}} \right) = \mathcal{N}\left( {\tau _{i,VA}^{{\rm{CE}}},\rho _{i,VA}^{{\rm{CE}}},{\phi _{VA}}} \right)p\left( {{\phi _{VA}}} \right)\nonumber\\ && \mathcal{N}\left( {\tau _{i,CSF}^{{\rm{CE}}},\rho _{i,CSF}^{{\rm{CE}}},{\phi _{CSF}}} \right)p\left( {{\phi _{CSF}}} \right)p\left( {\rho _i^{{\rm{CE}}}|{\eta ^{{\rm{CE}}}}} \right),\quad\end{eqnarray}

(10)

where \(\rho _i^{{\rm{CE}}} = ( {\rho _{i,{\rm{VA}}}^{{\rm{CE}}},\;\rho _{i,{\rm{CSF}}}^{{\rm{CE}}}} )\) is the mean of the \(\tau _i^{{\rm{CE}}}\;\)distribution and conditioned on η^CE, ϕ_VA and  ϕ_CSF are 2  ×  2 and 3  ×  3 covariances of \(\tau _{i,{\rm{VA}}}^{{\rm{CE}}}\;\)and \(\tau _{i,{\rm{CSF}}}^{{\rm{CE}}}\), with distributions p(ϕ_VA) and p(ϕ_CSF). 

The probability distribution \(p(\theta _{ij}^{{\rm{CE}}}|\tau _i^{{\rm{CE}}})\) of the CE_5D parameter \(\theta _{ij}^{{\rm{CE}}}\) for individual i in the j^th test is the product of two distributions conditioned on \(\tau _i^{{\rm{CE}}}\):  

\begin{eqnarray} p(\theta _{ij}^{{\rm{CE}}}|\tau _i^{{\rm{CE}}}) = p(\theta _{ij,VA}^{{\rm{CE}}}|\tau _{i,VA}^{{\rm{CE}}})p(\theta _{ij,CSF}^{{\rm{CE}}}|\tau _{i,CSF}^{{\rm{CE}}}).\quad\end{eqnarray}

(11)

 

The probability of obtaining the entire dataset, including all the row-by-row data from both VA and CSF assessments of all the individuals, is computed by probability multiplication:  

\begin{eqnarray}\begin{array}{@{}l@{}} p\left( {r_{1:I,1:J,1:M}^{VA},r_{1:I,1:J,1:K,1:L}^{CSF},{\rm{|}}X,o{s_{1:I,1:J,1:M}},{f_{1:I,1:J,1:K}},{c_{1:I,1:J,1:K,1:L}}} \right)\\ = \mathop \prod \limits_{i = 1}^I \mathop \prod \limits_{j = 1}^J p\left( {r_{ij,1:M}^{VA}{\rm{|}}\theta _{ij,VA}^{{\rm{CE}}},o{s_{ij,1:M}}} \right)p\left( {r_{ij,1:K,1:L}^{CSF}{\rm{|}}\theta _{ij,CSF}^{{\rm{CE}}},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} \right)\\ p\left( {\theta _{ij}^{{\rm{CE}}}{\rm{|}}\tau _i^{{\rm{CE}}}} \right)p\left( {\tau _i^{{\rm{CE}}}{\rm{|}}{\eta ^{{\rm{CE}}}}} \right)p\left( {{\eta ^{{\rm{CE}}}}} \right) = \mathop \prod \limits_{i = 1}^I \mathop \prod \limits_{j = 1}^J p\left( {r_{ij,1:M}^{VA}{\rm{|}}\theta _{ij,VA}^{{\rm{CE}}},o{s_{ij,1:M}}} \right)\\ p\left( {\theta _{ij,VA}^{{\rm{CE}}}{\rm{|}}\tau _{i,VA}^{{\rm{CE}}}} \right)p\left( {r_{ij,1:K,1:L}^{CSF}{\rm{|}}\theta _{ij,CSF}^{{\rm{CE}},},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} \right)p\left( {\theta _{ij,CSF}^{{\rm{CE}}}{\rm{|}}\tau _{i,CSF}^{{\rm{CE}}}} \right)\\ \mathcal{N}\left( {\tau _{i,VA}^{{\rm{CE}}},\rho _{i,VA}^{{\rm{CE}}},{\phi _{VA}}} \right)\mathcal{N}\left( {\tau _{i,CSF}^{{\rm{CE}}},\;\rho _{i,CSF}^{{\rm{CE}}},{\phi _{CSF}}} \right)p\left( {{\phi _{VA}}} \right) p\left( {{\phi _{CSF}}} \right)\\ p\left( {\rho _i^{{\rm{CE}}}|{\eta ^{{\rm{CE}}}}} \right)\mathcal{N}\left( {{\eta ^{{\rm{CE}}}},{\mu ^{{\rm{CE}}}},\Sigma } \right)p\left( {{\mu ^{{\rm{CE}}}}} \right)p\left( \sum \right) \end{array}\quad\end{eqnarray}

(12)

where \(X = ( {\theta _{1:I,1:J}^{{\rm{CE}}},\rho _{1:I}^{{\rm{CE}}},{\phi _{VA}},{\phi _{CSF}},{\mu ^{{\rm{CE}}}},\Sigma } )\) are all the parameters and hyperparameters in the HBJM, I is the total number of individuals, and J is the total number of tests on each individual. 

We start with prior distributions of μ^CE, Σ, ϕ_CSF, and ϕ_VA. 

\begin{eqnarray}{p_0}\left( {{\mu ^{{\rm{CE}}}}} \right) = \mathcal{U}\left( {\mu _{0,min}^{{\rm{CE}}},\mu _{0,max}^{{\rm{CE}}}} \right),\qquad\end{eqnarray}

(13a)

 

\begin{eqnarray}{p_0}\left( {{\Sigma ^{ - 1}}} \right) = \mathcal{W}\left( {\Sigma _{BIP}^{ - 1}/{\rm{\nu }},{\rm{\nu }}} \right),\qquad\end{eqnarray}

(13b)

 

\begin{eqnarray}{p_0}\left( {\phi _{VA}^{ - 1}} \right) = \mathcal{W}\left( {\phi _{qVA}^{ - 1}/{{\rm{\nu }}_{VA}},{{\rm{\nu }}_{VA}}} \right),\qquad\end{eqnarray}

(13c)

 

\begin{eqnarray}{p_0}\left( {\phi _{CSF}^{ - 1}} \right) = \mathcal{W}\left( {\phi _{qCSF}^{ - 1}/{{\rm{\nu }}_{CSF}},{{\rm{\nu }}_{CSF}}} \right),\qquad\end{eqnarray}

(13d)

where \(\mathcal{U}\) is a five-dimensional uniform distribution with lower bound \(\mu _{0,min}^{{\rm{CE}}} = ( {1.05,{\rm{\;}}0.1,{\rm{\;}}1,{\rm{\;}} - 0.5,{\rm{\;}}0.1} )\) and upper bound \(\mu _{0,max}^{{\rm{CE}}} = ( {1050,\;20,\;9,\;1.3,\;1.5} )\); precision matrices Σ⁻¹,\(\;\phi _{VA}^{ - 1}\), and \(\phi _{CSF}^{ - 1}\) are the inverse of covariances Σ, ϕ_VA, and ϕ_CSF; \(\mathcal{W}( {\frac{{\rm{Y}}}{{{{\rm{\nu }}_{\rm{Y}}}}},{{\rm{\nu }}_{\rm{Y}}}} )\) denotes a Wishart distribution with mean precision matrix Y and degrees of freedom ν_Y, with ν = 5, ν_VA = 2 , and ν_CSF = 3;\(\;\phi _{qVA}^{ - 1}\) and \(\phi _{qCSF}^{ - 1}\) are the inverse of the average covariance matrices ϕ_qVA and ϕ_qCSF computed from the BIP posterior distributions across all qVA and qCSF tests; and Σ⁻¹ is the inverse of the covariance matrix Σ computed from the estimated \(\theta _{ij}^{VA}\) and \(\theta _{ij}^{CSF}\) across all qVA and qCSF tests. 

The joint posterior distribution of all the parameters and hyperparameters in the HBJM is computed by Bayes’ rule:  

\begin{eqnarray}\begin{array}{@{}l@{}} p\left( {X{\rm{|}}Y_{1:I,1:J,1:M}^{VA},Y_{1:I,1:J,1:K,1:L}^{CSF}} \right)\\ = \frac{{\mathop \prod \nolimits_{i = 1}^I \mathop \prod \nolimits_{j = 1}^J p\left( {r_{ij,1:M}^{VA}{\rm{|}}\theta _{ij,VA}^{{\rm{CE}}}o{s_{ij,1:M}}} \right)p\left( {\theta _{ij}^{{\rm{CE}}}{\rm{|}}\tau _i^{{\rm{CE}}}} \right)p\left( {r_{ij,1:K,1:L}^{CSF}{\rm{|}}\theta _{ij,CSF}^{{\rm{CE}}},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} \right)}}{{\smallint \mathop \prod \nolimits_{i = 1}^I \mathop \prod \nolimits_{j = 1}^J p\left( {r_{ij,1:M}^{VA}{\rm{|}}\theta _{ij,VA}^{{\rm{CE}}},o{s_{ij,1:M}}} \right)p\left( {\theta _{ij}^{{\rm{CE}}}{\rm{|}}\tau _i^{{\rm{CE}}}} \right)p\left( {r_{ij,1:K,1:L}^{CSF}{\rm{|}}\theta _{ij,CSF}^{{\rm{CE}}},{f_{ij,1:K}},{c_{ij,1:K,1:L}}} \right)}}\\ \frac{{\mathcal{N}\left( {\tau _{i,VA}^{{\rm{CE}}},\rho _{i,VA}^{{\rm{CE}},},{\phi _{VA}}} \right)\mathcal{N}\left( {\tau _{i,CSF}^{{\rm{CE}}},\rho _{i,CSF}^{{\rm{CE}}},{\phi _{CSF}}} \right){p_0}\left( {{\phi _{VA}}} \right){p_0}\left( {{\phi _{CSF}}} \right)p\left( {\rho _i^{{\rm{CE}}}|{\eta ^{{\rm{CE}}}}} \right)\mathcal{N}\left( {{\eta ^{{\rm{CE}}}},{\mu ^{{\rm{CE}}}},\Sigma } \right){p_0}\left( {{\mu ^{{\rm{CE}}}}} \right){p_0}\left( \Sigma \right)}}{{\mathcal{N}\left( {\tau _{i,VA}^{{\rm{CE}}},\rho _{i,VA}^{{\rm{CE}}},{\phi _{VA}}} \right)\mathcal{N}\left( {\tau _{i,CSF}^{{\rm{CE}}},\rho _{i,CSF}^{{\rm{CE}}},{\phi _{CSF}}} \right){p_0}\left( {{\phi _{VA}}} \right){p_0}\left( {{\phi _{CSF}}} \right)p\left( {\rho _i^{{\rm{CE}}}|{\eta ^{{\rm{CE}}}}} \right)\mathcal{N}\left( {{\eta ^{{\rm{CE}}}},{\mu ^{{\rm{CE}}}},\Sigma } \right){p_0}\left( {{\mu ^{{\rm{CE}}}}} \right){p_0}\left( \Sigma \right)dX}}, \end{array}\quad\end{eqnarray}

(14)

where the denominator is the integral of the probability of obtaining the entire dataset across all possible values of X and is a constant for a given dataset and HBJM. 

Computing the joint posterior distribution

We used R (R Core Team, 2003) function run.jags in JAGS (Plummer, 2003) to compute representative samples of the joint posterior distribution of \(\theta _{i1}^{{\rm{CE}}}\;\)(5 parameters × 56 individuals=280 parameters), \(\rho _i^{{\rm{CE}}}\) (5 parameters × 56 subjects = 280 parameters), ϕ_CSF and ϕ_VA (9 parameters), μ^CE(5 parameters), and Σ (15 parameters) in three Markov Chain Monte Carlo (MCMC) chains. The MCMC efficiently samples the joint posterior distribution in a random walk (Kruschke, 2014). It started at a position randomly selected in the 589-dimensional parameter space. In each step, one of the 589 parameters was randomly selected. The values of all the other 588 parameters at the current position were fixed and the one-dimensional conditional probability distribution of the selected parameter was evaluated, from which a new value of the selected parameter was generated. The process reiterated until a pre-defined convergence criterion was reached so that the probability of visiting a position in the random walk approximated the joint posterior distribution of all the 589 parameters in Equation 14. 

Because the arbitrary starting position influenced the initial part of the random walk, steps in the adaptation and burn-in (optimization of the sampling algorithm by JAGS) phases were discarded and not included in the analysis. The exact number of steps discarded depends on each model and data. Convergence of each parameter was evaluated with Gelman and Rubin's diagnostic rule (Gelman & Rubin, 1992) based on the ratio of between- and within-MCMC variances along each dimension, that is, the variance of the samples across MCMC chains divided by the variance of the samples within each MCMC chain. The HBJM was considered “converged” when the Gelman and Rubin's statistics for all the parameters were smaller than 1.05. In this study, 500,000 and 20,000 steps were used for adaptation and burn-in in each MCMC chain based on pilot studies. The random walk continued until 2,000,000 additional samples were generated in each MCMC chain. 20,000 of the 2,000,000 samples (1 of every 100 sample) were stored (thinning ratio = 100). We fit the HBJM to three test lengths: short (the first 5 qVA and first 15 qCSF rows; M = 5, K = 15), medium (the first 15 qVA and first 25 qCSF rows; M = 15, K = 25), and long (all qVA and qCSF rows; M = 45, K = 50), separately. 

Similarly, we used R (R Core Team, 2003) function run.jags in JAGS (Plummer, 2003) to compute representative samples of the posterior distribution of \(\theta _{i1}^{VA}\) and \(\theta _{i1}^{CSF}\) for each qVA and qCSF test independently using their respective BIP. Five thousand and 5000 steps were used for adaptation and burn-in in each MCMC chain based on pilot studies. The random walk continued until 10,000 were stored in each MCMC chain with thinning ratio = 150, 50, and 50 for short, medium, and long test lengths, respectively. 

Statistical analysis

We evaluated the goodness of fit of the BIP and HBJM to the entire dataset using Bayesian predictive information criterion (BPIC) (Ando, 2007; Ando, 2011), which quantifies the likelihood of obtaining the observed data given the joint probability distribution of all the parameters of a model and penalizes model complexity. Smaller BPIC value indicates better fit. 

We then compared the uncertainties of the estimated \(\theta _{i1}^{VA},\theta _{i1}^{CSF}\), and \(\theta _{i1}^{{\rm{CE}}}\). The uncertainty was quantified as the variance of the marginal probability distribution of each parameter (Clayton & Hills, 1993; Edwards et al., 1963). Hotelling's T-squared test (Anderson, 2003) was used to compare the variances of \(\theta _{i1}^{{\rm{BIP}}} = (\theta _{i1}^{VA},\theta _{i1}^{CSF})\) and \(\theta _{i1}^{{\rm{CE}}}\) with R function HotellingsT2 (Nordhausen, Sirkia, Oja, & Tyler, 2018; R Core Team, 2003). 

Finally, we unblinded the data and computed sensitivity at 95% specificity, accuracy, and p value between the six Bangerter foil condition pairs (F0-F1, F0-F2, F0-F3, F1-F2, F1-F3, F2-F3) for (1) each eye based on distributions of \(\theta _{i1}^{VA}\), \(\theta _{i1}^{CSF}\), and \(\theta _{i1}^{{\rm{CE}}}\), and (2) group-distributions constructed by random sampling from distributions of \(\theta _{i1}^{VA}\), \(\theta _{i1}^{CSF}\), and \(\theta _{i1}^{{\rm{CE}}}\). 

At the test level for individual eyes, the distributions from the same eye in each condition pair were used for analysis based on VA threshold from the BIP (VA_1D), area under the log CSF from the BIP (AULCSF_1D), and the collective endpoint (CE_5D) from the HBJM. AULCSF_1D was constructed by computing the AULCSF between 1.5 and 18 cycles per degree. 

The group level null and difference distributions of VA_1D, AULCSF_1D, and CE_5D between a condition pair, p(x|G_null) and p(x|G_diff), where x = VA_1D, AULCSF_1D, or CE_5D, were constructed from the corresponding individual test-level distributions in the two conditions: (1) For each eye, a pair of random samples were drawn from the test-level distributions in the two conditions. The difference between the two samples was computed. (2) Step (1) was repeated for all 14 eyes. (3) The mean difference across the 14 eyes was added as a new sample in p(x|G_diff). (4) By repeating (1) ∼(3), p(x|G_diff) was constructed. And (5) the baseline distribution p(x|G_null) was constructed by shifting the mean of p(x|G_diff) to the origin. 

Similarly, for the power analysis, the group level null and difference distributions of VA_1D, AULCSF_1D, and CE_5D for sample sizes N = 1∼14 between each condition pair, p(x|G_null) and p(x|G_diff), where x = VA_1D, AULCSF_1D, or CE_5D, were constructed from the corresponding individual test-level distributions in the two conditions: (1) N eyes were randomly selected with replacement from the 14 eyes. (2) For each selected eye, a pair of random samples were drawn from the test-level distributions in the two conditions. The difference between the two samples was computed. (3) Step (2) was repeated for the N randomly selected eyes. (4) The mean difference across the N eyes was added as a new sample in p(x|G_diff). (5) By repeating (1) ∼(4) 10,000 times, p(x|G_diff) was constructed. (6) the baseline distribution p(x|G_null) was constructed by shifting the mean of p(x|G_diff) to zero (the origin). 

Probability density function (pdf) was obtained by fitting two-, three- and five- dimensional Gaussian mixture models (GMMs) to the MCMC samples of \(\theta _{i1}^{VA}\), \(\theta _{i1}^{CSF}\), and \(\theta _{i1}^{{\rm{CE}}}\), respectively, and the samples in the corresponding difference distributions. MATLAB function fitgmdist was used to fit GMMs with one to six Gaussian components. The Bayes information criterion (BIC) was computed for each GMM fit. The best fitting GMM was determined from the elbow of the BIC versus number of Gaussian component function. 

Given a baseline distribution p(x|B) and a treatment distribution p(x|T) (Figure 5A), specificity is the probability of correctly identifying a random sample x_s,B from p(x|B), and sensitivity is the probability of correctly identifying a random sample x_s,T from p(x|T) (Green & Swets, 1966; Swets & Pickett, 1982; US Department of Health and Human Services et al., 2007). We used the likelihood ratio \(\frac{{pdf( {x{\rm{|}}B} )}}{{pdf( {x{\rm{|}}T} )}}\) to define the criterion and computed sensitivity and specificity over a wide range of criteria to construct the receiver operating characteristic curves (ROC): (1-specificity) versus sensitivity. Accuracy is quantified by the area under ROC (AUROC), which is a criterion-free measure in signal detection theory (Figure 5B). (Green & Swets, 1966; Swets & Pickett, 1982; US Department of Health and Human Services et al., 2007) For a given pair of p(x|B) and p(x|T), the p value is defined as the probability of a random sample x_s,B from p(x|B) being incorrectly identified as from p(x|T) with \(\frac{{pdf( {x{\rm{|}}B} )}}{{pdf( {x{\rm{|}}T} )}} \le 1\) (Figure 5A):  

\begin{eqnarray}P = \mathop \smallint \nolimits_{pdf({x_{s,B}}|B) \le pdf({x_{s,B}}|T)}^{} p\left( {{x_{s,B}}{\rm{|}}B} \right)d{x_{s,B}}.\quad\end{eqnarray}

(15)

 

Figure 5.

 

(A) Scatter plots of x_s,B (red dots) from the base-line distribution p(x|B) and x_s,T (blue dots) from the treatment distribution p(x|T). Black dashed line represents a criterion of \(\frac{{pdf( {x{\rm{|}}B} )}}{{pdf( {x{\rm{|}}T} )}} = 1\) and the iso-density contours illustrate the best fitting Gaussian mixture models to the two distributions. (B) The bold black curve is the ROC curve. The grey area under the ROC is the AUROC, a criterion-free estimate of the test accuracy. Black dashed line represents chance performance.

In this article, p(x|B) and p(x|T) are test-level distributions in two conditions for individual eyes and p(x|G_null) and p(x|G_diff) at the group level. We also computed the p values (Equation 15) between \(\theta _{i1}^{VA}\) and the VA component of \(\theta _{i1}^{{\rm{CE}}}\), and between \(\theta _{i1}^{CSF}\) and the CSF parameter components of \(\theta _{i1}^{{\rm{CE}}}\) to quantify the agreement between BIP and HBJM estimates in each test, with BIP as the baseline. 

Results

We report the detailed results based on the medium test length in the main text. The results from the short and long test lengths were largely consistent and can be found in Supplementary Materials B and C. 

Goodness-of-fit

The HBJM (BPIC_HBJM = 4591) fit the data better than the BIP (BPIC_BIP = 4617). 

Posterior distributions from the HBJM and the BIP

Figure 6 illustrates the two-dimensional marginal posterior distributions of η^CE at the population level. We found large correlations among η^CE components (Table 1). Specifically, \(\eta _{Threshold}^{{\rm{CE}}}\) was positively correlated with \(\eta _{Range}^{{\rm{CE}}}\) (r = 0.516), but negatively correlated with \(\eta _{PG}^{{\rm{CE}}}\) (r = −0.486), \(\eta _{PF}^{{\rm{CE}}}\) (r = −0.583), and \(\eta _{BH}^{{\rm{CE}}}\) (r = −0.288); \(\eta _{range}^{{\rm{CE}}}\) was negatively correlated with \(\eta _{PF}^{{\rm{CE}}}\) (r = −0.494); \(\eta _{BH}^{{\rm{CE}}}\) was positively correlated with \(\eta _{PG}^{{\rm{CE}}}\) (r = 0.577), but was negatively correlated with \(\eta _{PF}^{{\rm{CE\;}}}\)(r = −0.536). The positive correlation between \(\eta _{Threshold}^{{\rm{CE}}}\) and \(\eta _{Range}^{{\rm{CE}}}\) is consistent with our previous observation that worse VA thresholds were associated with wider ranges (shallower slope of VABF) across individuals (Zhao, Lesmes, Dorr, Bex, et al., 2021; Zhao, Lesmes, Dorr, & Lu, 2021). The negative correlations between \(\eta _{VA}^{{\rm{CE}}}\) and \(\eta _{CSF}^{{\rm{CE}}}\) components were expected because increased VA threshold is associated with worse vision, and increased CSF is associated with better vision. The correlations among \(\eta _{CSF}^{{\rm{CE}}}\) components indicated that individuals with wider CSF bandwidths tended to have higher peak gains but lower peak frequencies. 

Figure 6.

 

Two-dimensional marginal posterior distributions of η^CE.

Table 1.

 

Correlations of η^CEcomponents at the population level.

 

Figure 7 illustrates two-dimensional marginal posterior distributions of \(\tau _i^{{\rm{CE}}}\) for one individual. Table 2 shows the correlations of \(\tau _i^{{\rm{CE}}}\) components, averaged across all individuals. \(\tau _{i,Threshold}^{{\rm{CE}}}\) and \(\tau _{i,Range}^{{\rm{CE}}}\) did not exhibit large correlation, consistent with our previous publication (Zhao, Lesmes, Dorr, & Lu, 2021), because they quantified two independent properties of the VABF, namely the position and shape of the psychometric function at the individual level (Figure 2A). We found large correlations among \(\tau _{i,CSF}^{{\rm{CE}}}\;\)components: \(\tau _{i,PF}^{{\rm{CE}}}\) was negatively correlated with \(\tau _{i,PG}^{{\rm{CE}}}\;\)(r = −0.757 ± 0.028) and \(\tau _{i,BH}^{{\rm{CE}}}\) (r = −0.900 ± 0.010), whereas \(\tau _{i,PG}^{{\rm{CE}}}\) was positively correlated with \(\tau _{i,BH}^{{\rm{CE}}}\;\)(r = 0.551 ± 0.047). The correlations among \(\tau _{i,CSF}^{{\rm{CE}}}\) components resulted from the properties of the generative model of CSF: The three CSF model parameters jointly determine contrast sensitivities at all the spatial frequencies and different combinations of parameter values may result in similar probabilities of correct responses in some experimental conditions. On the other hand, correlations between \(\tau _{i,VA}^{{\rm{CE}}}\) and \(\tau _{i,CSF\;}^{{\rm{CE}}}\)were small because VA and CSF assessments were conducted separately. 

Figure 7.

 

Two-dimensional marginal posterior distributions of \(\tau _i^{{\rm{CE}}}\) for one individual from the HBJM.

Table 2.

 

Correlations of \(\tau _i^{{\rm{CE}}}\) components, averaged across all individuals and with SD in parentheses.

 

The two-dimensional marginal posterior distributions of \(\theta _{i1}^{{\rm{CE}}}\) at the test level for the same individual are illustrated in Figure 8. The posterior distributions of \(\theta _{i1}^{VA}\) and \(\theta _{i1}^{CSF}\) from the BIP are illustrated in Figure 9. The corresponding correlations, averaged across all tests from the two analyses, are listed in Tables 3 and 4. We found large correlations among \(\theta _{i1,{\rm{CSF}}}^{{\rm{CE}}}\) components and among \(\theta _{i1}^{CSF}\) components. \(\theta _{i1,PF}^{{\rm{CE}}}\) was negatively correlated with \(\theta _{i1,PG\;}^{{\rm{CE}}}\)(r = −0.730 ± 0.210) and \(\theta _{i1,BH}^{{\rm{CE}}}\) (r = −0.930 ± 0.024), whereas \(\theta _{i1,PG}^{{\rm{CE}}}\) was positively correlated with \(\theta _{i1,BH\;}^{{\rm{CE}}}\)(r = 0.520 ± 0.271). \(\theta _{i1,PF}^{CSF}\) was negatively correlated with \(\theta _{i1,PG}^{CSF}\) (r = −0.762 ± 0.271) and \(\theta _{i1,BH}^{CSF}\) (r = −0.931 ± 0.026), whereas \(\theta _{i1,PG}^{CSF}\) was positively correlated with \(\theta _{i1,BH}^{CSF}\) (r = 0.562 ± 0.355). These patterns are consistent with those at the individual level and our previous publications (Zhao, Lesmes, Dorr, Bex, et al., 2021; Zhao, Lesmes, Dorr, & Lu, 2021; Zhao, Lesmes, Hou, et al., 2021). 

Figure 8.

 

Two-dimensional marginal posterior distributions of \(\theta _{i1}^{{\rm{CE}}}\) at the test level for the same individual in Figure 7 from the HBJM.

Figure 9.

 

Two-dimensional marginal posterior distributions of \(\theta _{i1}^{VA}\) and \(\theta _{i1}^{CSF}\) for the same individual in Figure 7 from the BIP.

Table 3.

 

Correlations of \(\theta _{i1}^{{\rm{CE}}}\) components, averaged across all tests and with SD in parentheses.

 

Table 4.

 

Correlations of \(\theta _{i1}^{CSF}\) and\(\;\theta _{i1}^{VA}\) components, averaged across all tests and with SD in parentheses.

 

Mean and uncertainty of test parameters

The BIP and HBJM estimates of the parameters at the test level showed excellent agreement. The p values for the difference between the posterior distributions of \(\theta _{i1}^{VA}\) and \(\theta _{i1}^{CSF}\) and those of the corresponding \(\theta _{i1}^{{\rm{CE}}}\) components were 0.20 to 0.41 and 0.20 to 0.45. The average variance of the estimated \(\theta _{i1}^{{\rm{CE}}}\) components (mean = 0.0091 log10 units; range 0.00026∼0.059 log10 units) was 71.5% less than that of the corresponding \(\theta _{i1}^{{\rm{BIP}}}\) components (mean = 0.032 log10 units; range 0.00026∼0.30 log10 units) (t² (5, 51) = 46, p < 0.001), consistent with the well-known variance shrinkage effect of HBMs (Kruschke, 2014). Figure 10 shows histograms of the variance of \(\theta _{i1}^{{\rm{BIP}}}\) and \(\theta _{i1}^{{\rm{CE}}}\) components. 

Figure 10.

 

Histograms of the variance of \(\theta _{i1}^{{\rm{BIP}}}\) and \(\theta _{i1}^{{\rm{CE}}}\) components.

Signal-to-noise analysis

The mean VA (Figures 11A–C) and AULCSF (Figures 11D–F) changes between pairs of foil conditions were virtually equal across BIP and CE analyses in all test durations. The recovery of strong covariances by CE analysis reduced the variance of CE estimates relative to BIP estimates (Figures 11G–L), especially for the short testing duration, consistent with our previous HBM results in each single modality of CSF and VA (Zhao, Lesmes, Dorr, & Lu, 2021; Zhao, Lesmes, Hou, et al., 2021). In addition, CE_5D combines signals from VA and CSF. The increased signal and reduced noise of CE_5D improved signal detection. At the test level for individual eyes, median p values only for CE_5D from the five-dimensional HBJM analysis were all significant for short (p < 0.02), medium (p < 0.002), and long testing (p < 0.0004) (Figure 12). The advantage provided by the HBJM was stronger when less information was available with the short and medium test lengths relative to those with the long test length. 

Figure 11.

 

(A–C) Average signal between F1 and F2, F1 and F3, and F2 and F3 for VA_1D and the VA component of CE_5D across 14 eyes. (D–F) Average signal between F1 and F2, F1 and F3, and F2 and F3 for AULCSF_1D and AULCSF based on CE_5D. (G–I) Average variance of the difference distributions between F1 and F2, F1 and F3, and F2 and F3 for VA_1D and the VA component of CE_5D. (J–L) Average variance of the difference distributions between F1 and F2, F1 and F3, and F2 and F3 for AULCSF_1D and AULCSF based on CE_5D. (1dB = 0.1 log10 units).

Figure 12.

 

Medium p values between F1 and F2 (A), F1 and F3 (B,), and F2 and F3 (C) for VA_{1D ,} AULCSF_1D, and CE_5D in the 14 eyes.

Sensitivity at 95% specificity, accuracy, and p value

Since the differences between condition pairs F0-F1, F0-F2, and F0-F3 were large, sensitivity at 95% specificity and accuracy were all close to 100% and the p values were close to 0 for VA_1D, AULCSF_1D, and CE_5D at both the individual and group levels. CE_5D exhibited advantage over VA_1D and AULCSF_1D between condition pairs with smaller differences (F1-F2, F2-F3, F1-F3), and the advantage was more profound when less information was collected with shorter test lengths or smaller numbers of eyes (see results for the short and long test lengths in Supplementary Materials B and C). 

Test level

Tables 5 and 6 show the average sensitivity at 95% specificity and accuracy across all eyes between the six condition pairs at the test level, with standard error (SE) in parentheses. CE_5D exhibited significantly higher sensitivity at 95% specificity than VA_1D between F1 and F2 (t(13) = 2.9, p = 0.006), between F1 and F3 (t(13) = 1.9, p = 0.04), and between F2 and F3 (t(13) = 3.5, p = 0.002), and AULCSF_1D between F1 and F2 (t(13) = 2.7, p = 0.009), and between F2 and F3 (t(13) = 2.3, p = 0.02). CE_5D exhibited significantly higher accuracy than VA_1D between F1 and F2 (t(13) = 2.6, p = 0.01), and between F2 and F3 (t(13) = 3.2, p = 0.004), and AULCSF_1D between F1 and F2 (t(13) = 2.2, p = 0.02), and between F2 and F3 (t(13) = 2.0, p = 0.03). Table 7 shows that larger number of eyes reached p < 0.01 with CE_5D relative to VA_1D and AULCSF_1D between F1 and F2, between F1 and F3, and between F2 and F3, indicating highest discriminative power for CE_5D. 

Table 5.

 

Average (SE) sensitivity (%) at 95% specificity of VA_1D, AULCSF_1D, and CE_5D across all eyes between six condition pairs.

 

Table 6.

 

Average (SE) accuracy (%) of VA_1D, AULCSF_1D, and CE_5D across all eyes between six condition pairs.

 

Table 7.

 

Number of eyes required to reach significant differences (p < 0.01) based on VA_1D, AULCSF_1D, and CE_5D.

 

Group level

At the group level, sensitivity at 95% specificity and accuracy between all condition pairs based on VA_1D, AULCSF_1D, and CE_5D were all 100% and the p values were all <0.001, except for those of VA_1D between F2 and F3 with sensitivity and p values of 99% and 0.03. 

Power analysis

CE_5D exhibited much more statistical power than VA_1D and AULCSF_1D (Figure 13). Only five, two, and five eyes were needed for CE_5D to reached 95% sensitivity at 95% specificity between F1 and F2, F1 and F3, and F2 and F3. Many more eyes were needed for VA_1D and AULCSF_1D to reach the same level of specificity. Similar results were obtained in terms of accuracy and p values. 

Figure 13.

 

(A–C) Sensitivity at 95% specificity between F1-F2, F1-F3, and F2-F3 condition pairs as functions of sample size for VA_1D, AULCSF_1D, and CE_5D. (D–F) Accuracy between F1-F2, F1-F3, and F2-F3 condition pairs. (G–I) p value between F1-F2, F1-F3, and F2-F3 condition pairs.

Discussion

In this article, we developed a three-level HBJM framework to compute a collective endpoint (CE_5D) from multiple assessment modalities by modeling data from two test modalities. Using information across all individuals and test modalities, the framework explicitly quantifies the statistical relationships between multiple test modalities and generates a multidimensional collective endpoint that can improve statistical power. We applied the framework to compute the five-dimensional CE_5D, consisting of VA threshold, VA range, peak gain, peak spatial frequency, and bandwidth, from trial-by-trial qVA and qCSF data. The HBJM revealed strong correlations among CE_5D components at the population, individual and test levels. Compared with the one-dimensional endpoints from BIP analysis (VA_1D and AULCSF_1D), the five-dimensional CE_5D from the HBJM generated the highest sensitivity at 95% specificity, the highest accuracy, and the most significant p values at both the test and group levels. The higher sensitivity and accuracy and more significant p values of CE_5D resulted mainly from the HBJM, which optimally combines information from test modalities, individuals, and tests to mutually constrain the estimates to reduce the uncertainties of the parameters and hyperparameters of the VA and CSF models (Figures 11 and 12). 

Because of the strong correlations between the VA and CSF components of CE_5D at the population level, the HBJM can be used to generate informative priors for qCSF based on qVA assessment or priors for qVA based on qCSF assessment by treating the to-be-observed data as missing data in the model. We previously developed the hierarchical adaptive design optimization (Gu et al., 2016; Kim et al., 2014) algorithm to simultaneously categorize each new patient into existing groups and provide the most informative prior for them based on their group membership for qCSF. A combination of hierarchical adaptive design with the collective endpoint analysis developed in this article may further improve the efficiency of adaptive tests in multiple assessment modalities. 

Although it was applied to a dataset with only two assessment modalities in this study, the HBJM can be extended to compute CE from more than two assessment modalities. A generative model of each endpoint is required in the HBJM framework. In the current study, generative models of two functional vision endpoints (VA and CSF) enabled our development of CE_5D of VA and CSF. Generative models of other endpoints including clinical events (Yin et al., 2018) and structural endpoints can be developed and incorporated into the HBJM framework. They can be either statistical models that capture the probability of certain aspects of a human condition or mechanistic models based on human biology or physiology. For example, optical coherence tomography (OCT) (Huang et al., 1991) uses light to generate cross-section retina images. Anatomical changes (e.g., thickness of a retina layer) identified with OCT have been related to the severities of eye diseases (Vujosevic et al., 2022). A generative model can be constructed by specifying the probability of an eye disease and its severity as a function of the anatomical features of the retina, such as the thickness of retinal layers in different locations. The HBJM can be extended to compute CE from OCT, VA, and CSF and further facilitate disease staging and treatment monitoring and improve statistical evaluation in clinical trials. 

In this article, we constructed a three-level HBJM without considering any hierarchy related to treatment conditions because our goal was to compute the CE for each individual from blinded data at one point of assessment. The HBJM can be extended to compute CE from blinded data at different time points. In this case, the CE includes the parameters of all the assessment modalities at multiple time points to capture treatment effects over time and can be used to improve statistical evaluation of treatment efficacy. The HBJM can also be extended to compute CE from unblinded data in clinical trials with multiple endpoints at multiple time points. In this case, the HBJM takes into account of the treatment arms. Further research is necessary to investigate the utilities of these two forms of HBJM in clinical trials. 

In conclusion, the HBJM extracts and optimally combines valuable information about covariance of CS and VA parameters and improves precision of the estimated parameters. By combining signals and reducing noise from multiple tests for detecting vision changes, this framework exhibits potential to improve statistical power for combining multimodality data in ophthalmic trials. 

Acknowledgments

Supported by the National Eye Institute (EY017491 and EY032125). 

Commercial relationships: L.A.L and M.D. have intellectual and equity interests, and hold employment in Adaptive Sensory Technology, Inc. (San Diego, CA). L.A.L and M.D. hold employment in AST. Z.-L.L. has intellectual and equity interests in Adaptive Sensory Technology, Inc. (San Diego, CA) and Jiangsu Juehua Medical Technology, LTD (Jiangsu, China). Y.Z. has no competing interest. 

Corresponding author: Zhong-Lin Lu. 

Email: [e-mail]. 

Address: Center for Neural Science, New York University, 4 Washington Place, New York, NY 10003, USA.