Using Hierarchical Bayesian Modeling to Enhance Statistical Inference on Contrast Sensitivity

Abstract

Purpose: The purpose of this study is to introduce a nonparametric hierarchical Bayesian model (HBM) that enables advanced statistical inference on contrast sensitivity (CS) both at individual spatial frequencies (SFs) and across multiple SFs in clinical trials, where CS measurements are crucial for assessing safety and efficacy.

Methods: The HBM computes the joint posterior distribution of CS at six Food and Drug Administration–designated SFs across the population, individual, and test levels. It incorporates covariances at both population and individual levels to capture the relationship between CSs across SFs. A Bayesian inference procedure (BIP) is also used to estimate the posterior distribution of CS at each SF independently. Both methods are applied to a quantitative CSF (qCSF) dataset of 112 subjects and compared in terms of precision, test-retest reliability of CS estimates, sensitivity, accuracy, and statistical power in detecting CS changes.

Results: The HBM reveals correlations between CSs in pairs of SFs and provides significantly more precise estimates and higher test-retest reliability compared to the BIP. Additionally, it improves the average sensitivity and accuracy in detecting CS changes for individual subjects, as well as statistical power for detecting group-level CS changes at individual and combinations of multiple SFs between luminance conditions.

Conclusions: The HBM establishes a comprehensive framework to enhance sensitivity, accuracy, and statistical power for detecting CS changes in hierarchical experimental designs.

Translational Relevance: The HBM presents a valuable tool for advancing CS assessments in the clinic and clinical trials, potentially improving the evaluation of treatment efficacy and patient outcomes.

Introduction

The contrast sensitivity function (CSF) is an important measure of spatial vision,^1–3 offering valuable insights into changes in vision associated with the progression of eye diseases and their treatment.^4–9 Various diseases and different stages within a single disease can impact contrast sensitivity (CS) at different spatial frequencies (SF). For instance, glaucoma is particularly known for impairing CS at low SFs,¹⁰ whereas conditions like age-related macular degeneration (AMD),^11–13 cataracts^14,15 (but see reference 16), diabetic retinopathy,^17,18 and idiopathic epiretinal membrane¹⁹ tend to affect CS predominantly at low-to-intermediate SFs. Amblyopia,^20,21 on the other hand, causes CS impairments at intermediate-to-high SFs. Studies have also demonstrated that AMD,²² inherited retinal degeneration,²³ macula-off retinal detachment,²⁴ multiple sclerosis,^25,26 and myopia^27–31 can cause CS impairments across low, intermediate, and high SFs. Consequently, to accurately characterize and monitor disease progression for individual patients in clinical settings and to quantify treatment efficacy for specific treatments in clinical trials, it is imperative to assess CS at multiple SFs. Evaluating CS changes across multiple SFs is crucial for its consideration as a safety and efficacy endpoint in clinical trials.^32,33 

Adaptive procedures have traditionally been used to measure CS independently across a predetermined set of SFs in laboratory settings. In conventional designs, introducing a new SF condition necessitates a minimum number of experimental trials, often ranging from 50 to 100 trials.³⁴ Sampling the CSF at five to 10 SFs typically demands 500 to 1000 trials, translating to 30 to 60 minutes of data collection. This volume of data, although acceptable for experiments focusing on the CSF in a single condition, becomes impractical for measuring multiple CSFs (e.g., for different eyes), even in controlled laboratory environments.³⁵ Conversely, the various CSF charts used in clinical settings (e.g., Arden cards, Vistech, FACT charts; see reference 36 for a review) exhibit limited flexibility and reliability.^37–41 

Recognizing the significance of CSF in both basic and clinical vision and the absence of a precise and efficient assessment instrument prompted the development of the quantitative CSF (qCSF) test.³⁵ The method models the CSF as a four-parameter truncated log parabola function^42,43 and uses a Bayesian active learning algorithm to optimally estimate the posterior distribution of the four parameters by selecting the stimulus for each subsequent trial that maximizes the expected information gain.^35,44 Using a 10-letter identification task, the assessment of the CSF can be achieved in about 20 trials—approximately two minutes—yielding an average standard deviation (SD) of 0.10 log₁₀ units of the CSs estimated from the parametric CSF model.^45,46 A recent analysis showed that qCSF, with the optimal information gain strategy, generated significantly more expected information gain than the Pelli-Robson and the CSV-1000 tests.⁴⁷ 

Since its debut, qCSF has been used to evaluate spatial vision in both normal^48–53 and clinical populations. Applications include aging,^51,54 amblyopia,^20,21,55 cataracts,^15,56 central serous chorioretinopathy, diabetic retinopathy,^17,18,57,58 dysthyroid optic neuropathy,⁵⁹ Fuchs uveitis syndrome,⁵⁶ glaucoma,^10,60–62 idiopathic epiretinal membrane,¹⁹ inherited retinal degeneration,²³ keratoconus,^63,64 multiple sclerosis,^25,26 myopia^27–31, and various maculopathies⁶⁵ such as age-related macular degeneration,^{11–13,22,66–68} geographic atrophy,^12,69 macula-off retinal detachment,²⁴ and retinal vein occlusion.^70,71 However, most studies with qCSF³⁴ have used the area under the log CSF (AULCSF),^{10–12,14,15,17–19,22,24,26–30,46,50,54–57,59,62,63,67,69,71–80} whereas some have used CS at individual SFs^{10–12,14,15,17–19,22,24,26–30,50,51,54,56,59,63,71,73–75,78–80} for statistical inference because of the undesirable mathematical properties associated with the CSs computed from the parametric CSF model.⁴³ Relying on AULCSF or CS at single SFs overlooks the rich information available across SFs in the CSF, thereby limiting the full potential of the CSF test. 

In the four-parameter truncated log parabola CSF model, CSs across all SFs are determined by the mathematical function. Consequently, the CS values generated at different SFs are not independent, given their deterministic relationship dictated by the model. This means that the CS at any SF can be derived if CS values at four specific SFs are known. When computing samples of CSF from the posterior distributions of the model parameters, CSs between pairs of SFs exhibit correlation due to their deterministic connection, reflecting the mathematical properties of the CSF model rather than the true correlation between independently estimated CS values at each SF. To perform meaningful statistical inference on CS at individual SFs or combinations of multiple SFs, it is essential to eliminate the a priori deterministic relationship introduced by the parametric model while considering the true correlations of CS values across SFs. 

On the other hand, estimating CS from qCSF data at individual SFs without a parametric CSF model is highly challenging because of the limited number of trials at each SF in a qCSF assessment (Fig. 1). Figure 1a illustrates the distribution of test stimuli in a typical 25-trial qCSF assessment,⁴⁶ binned at the six SFs designated by the Food and Drug Administration (FDA) (Table). Note that each qCSF trial consists of three optotypes with the same SF characteristics but different contrasts and that SFs near and greater than 18 cpd were labeled as 18 cpd in this study. Attempting to estimate CS at each SF independently with this amount of data is exceptionally challenging and can lead to significant uncertainty. 

Figure 1.

 

(a) Data from a typical subject in a 25-trial qCSF test binned at the six FDA-designated SFs. There were three optotypes with the same SF characteristics but different contrasts in each trial. CS is the reciprocal of contrast threshold: CS = 1/threshold. Log(CS) = −log₁₀(contrast threshold), where contrast threshold is between 0 and 1.⁴⁶ (b) Estimated CSF from the qCSF data in (a) using the BIP, with error bars indicating SDs. (c) Average number of qCSF trials at each of the six SFs based on the data from 112 subjects tested in each of three luminance conditions with 25 qCSF trials.⁴⁶ (d) Average SD of the estimated CSs at the six SFs from the 112 subjects using the BIP, after excluding those with CS < 0.0 log₁₀ units.

Table.

 

Bins for the Six FDA-Designated SFs

 

Figure 1b displays the estimated CSF from the qCSF data in Figure 1a, obtained by fitting a psychometric function to the data at each SF using the Bayesian inference procedure (BIP) described later in this article. With an average SD of 0.23 ± 0.02 log₁₀ units across SFs, which is equal to 70% of the estimated CS value in linear units, the estimated CSF is highly imprecise. 

Figure 1c shows the distribution of the number of qCSF trials at each of the six SFs based on the data from 112 subjects tested in three luminance conditions with 25 qCSF trials.⁴⁶ On average, there are 2.5 ± 0.1, 4.3 ± 0.1, 3.0 ± 0.1, 5.0 ± 0.11, 4.2 ± 0.1, and 6.0 ± 0.1 trials at the six SFs. Figure 1d illustrates the distributions of the SDs of the estimated CS at each of the six SFs based on the data of the 112 subjects using the BIP, after excluding those with CS < 0.0 log₁₀ units. With an average SD of 0.252 ± 0.003 log₁₀ units across all SFs and subjects, which is equal to 79% of the estimated CS values in linear units, the estimated CSs are highly imprecise. 

In this investigation, our aim was to enable advanced statistical inference on changes of CS at individual SFs and across multiple SFs by analyzing the joint posterior distribution of CS across subjects, SFs and experimental conditions. We introduced a nonparametric hierarchical Bayesian model (HBM) structured across population, individual and test levels to compute the joint posterior distribution of CS at the six FDA-designated SFs across all three levels. For comparison, we also implemented the traditional BIP, which computes the posterior distribution of CS at each SF independently. 

Within the Bayesian inference framework, two critical elements influencing the accuracy and precision of the estimated posterior distribution, given a fixed likelihood function, are the prior and the amount of data. The prior represents the probability distribution of the parameters to be estimated before collecting new data. In the ideal scenario, the prior serves as a mathematical description of our knowledge about the parameters before data collection, often expressed as a uniform distribution in cases of little or no prior knowledge (Fig. 2a, column 1), or as a concentrated informative distribution with abundant knowledge (Fig. 2b, column 1). In cases of misinformed priors (Fig. 2c, column 1), bias may be introduced. The more informative the prior, the higher the accuracy and precision of the estimated posterior distribution with the same amount of data (Fig. 2, columns 2–5). To achieve target levels of precision and accuracy, the most informative prior (Fig. 2b) requires the least amount of data, whereas the biased prior demands the most (Fig. 2c). As the amount of data increases, the impact of the prior diminishes (Fig. 2, column 5). This article primarily focuses on constructing informative priors using the HBM because of the limited data available from qCSF assessments. By explicitly modeling the covariance of CSs at the population and individual levels, as well as conditional dependencies across the three levels, the HBM generates an informative prior for CS at each SF in each test by incorporating information across all SFs and tests in a dataset. 

Figure 2.

 

Effects of the prior and the amount of data on the estimated posterior distribution with 2, 5, 10, and 25 trials of data in Bayesian inference. (a) An uninformative uniform prior and estimated posterior distributions. (b) An informative concentrated prior and estimated posterior distributions. (c) A biased prior and estimated posterior distributions. The dotted vertical lines indicate the true CSs.

The HBM functions as a generative model framework that leverages Bayes’ rule to quantify the joint distribution of population-, individual-, and test-level hyperparameters and parameters.^81–85 It can explicitly quantify correlations in the data through covariance parameters.^86–90 By sharing information within and across levels via conditional dependencies, it generates informative priors for each test and can therefore reduce the variance of test-level estimates by shrinking estimated parameters at the lower levels toward the modes of the higher levels when there is insufficient data at the lower levels.^81,84,91 

The HBM has found applications in diverse scientific disciplines, including astronomy,⁹² ecology,^93,94 genetics,⁹⁵ machine learning,⁹⁶ and cognitive science.^{82,84,85,91,97–102} Previous applications demonstrated that the HBM reduces uncertainties in the estimated parameters of the CSF¹⁰³ and visual acuity (VA) behavioral function¹⁰⁴ as well as uncertainties in estimated learning curves^105,106 compared to the BIP. Another development introduced a hierarchical Bayesian joint modeling (HBJM) framework to compute the collective endpoint (CE) from CSF and VA measurements, showing that CE offers more statistical power than the CSF or VA metrics alone.¹⁰⁷ (In the HBJM, CSF was modeled as a three-parameter log parabola function, rather than as CSs at specific SFs in the current HBM. Additionally, the HBJM modeled the relationship between CSF and VA parameters, whereas the current HBM modeled the relationship between CSs at different SFs.) 

In this study, we applied the HBM and BIP to both the first 25 and all 50 trials of a publicly available dataset involving 112 subjects. These subjects underwent 50 qCSF trials in each of three luminance conditions, with two repeated tests in the high luminance condition.⁴⁶ The selected luminance levels (L, M, and H) were designed to induce 0.14, 0.29, and 0.43 log₁₀ units AULCSF changes, mimicking mild, medium, and large CSF deficits in clinical populations.^108–113 After acquisition of posterior distributions of CS from the HBM and BIP, we conducted comparisons across various aspects: (1) mean and uncertainty of the estimated CSs, (2) test-retest reliability of the estimated CSs in the H condition, (3) sensitivity at 95% specificity and accuracy in detecting CS changes between luminance conditions for each subject, and (4) statistical power in detecting CS changes between luminance conditions at individual and across multiple SFs at the group-level. These comparative analyses aimed to provide insights into the performance of the two methods and to highlight the significance of the HBM in enhancing statistical inference on CS at individual SFs and across multiple SFs. 

Methods

Dataset

The dataset comprised 112 subjects who underwent 50 qCSF trials under three luminance conditions (L, M, and H).⁴⁶ Each qCSF trial consisted of three equal-size bandpass-filtered optotypes randomly chosen with replacement from the ten Sloan letters (C, D, H, K, N, O, R, S, V, and Z), with sizes (and therefore center SFs) and contrasts determined by qCSF.^35,45 Subjects were tested once in the L and M conditions and twice in the H condition, resulting in a total of 448 qCSF tests. Written consent was obtained from all subjects, and the study protocol adhered to the tenets of the Declaration of Helsinki, approved by the institutional review board of human subject research of the Ohio State University. 

Apparatus

The qCSF was implemented in MATLAB (MathWorks Corp., Natick, MA, USA) on a 55-inch Samsung UN55FH6030 monitor with a resolution of 1920 × 1080 pixels. A bit-stealing algorithm achieved a nine-bit grayscale resolution.¹¹⁴ Subjects viewed the displays binocularly at a distance of four meters.⁴⁶ Data analysis was conducted using a Dell computer with Intel Xeon W-2145 @ 3.70 GHz CPU (8 cores and 16 threads) and 64GB installed memory (RAM), using MATLAB and JAGS¹¹⁵ in R.¹¹⁶ 

Data Analytic Procedures

Supplementary Materials A contains detailed descriptions of our data analytic procedures. Briefly, the BIP (Fig. 3a) estimates the posterior distribution of CS independently for each individual at each SF in each test. It treats the data in isolation, inferring the posterior distribution of CS at each SF in each test independently, without considering any relationships among CSs across SFs, tests, or individuals. On the other hand, the HBM (Fig. 3b) uses the entire dataset to constrain CS estimates across all tests, individuals, and SFs. It quantifies the relationship between CSs across all SFs at the population and individual levels. Consequently, the estimated CSs from the HBM of all tests, individuals, and SFs mutually constrain each other through the three-level hierarchy and covariances. 

Figure 3.

 

(a) The BIP estimates the posterior distribution of CS independently for each individual at each SF in each test. (b) The HBM uses the entire dataset to constrain CS estimates across all tests, individuals, and SFs through a three-level hierarchy and covariances.

Results

In the main text, we present results from the HBM and BIP based on the first 25 qCSF trials (L = 25). Additionally, the findings from all 50 qCSF trials (L = 50) are generally consistent and reported in Supplementary Material C. 

Posterior Distributions From the HBM and BIP

The HBM quantified the relationships between CSs across SFs. It recovered correlations between CSs in pairs of SFs, with r = 0.66 to 0.90 at the population, r = 0.01 to 0.52 at the individual, and r = −0.01 to 0.45 at the test levels. The large positive correlations across all SFs at the population level are consistent with previous findings that, across subjects, better CS at one SF was associated with better CS at other SFs.^42,43 On the other hand, the positive correlations between CSs at similar SFs at the individual and test levels are consistent with the idea of SF channels underlying the CSF, each selective for a range of SFs.^117–119 The posterior distributions from the HBM and BIP, along with the correlation tables are presented in Supplementary Materials B. 

Goodness of Fit

The Bayesian predictive information criterion (BPIC)^120,121 quantifies the goodness of fit of a model to the trial-by-trial data by evaluating the likelihood of the data based on the joint posterior distribution of the parameters while penalizing model complexity. The smaller the BPIC, the better the fit. The BPIC for the HBM and BIP were 36216 and 48867, respectively. The HBM fit the data significantly better than the BIP (χ²(1928) = 12, 651;  P < 0.001;  Bayes factor = 10¹¹⁴³). 

Mean and Standard Deviation of the Posterior Distributions of the Estimated CSs

Figures 4a–c depict the estimated CSFs of a typical subject from the HBM and BIP in the three luminance conditions. The SD of the posterior distribution of the estimated CS at each SF at the test level (Figs. 4d–f) was used to quantify the precision of the estimates. The smaller the SD, the higher the precision. For this subject, the average SDs of the estimated CSs across all SFs and luminance conditions were 0.106 ± 0.007 log₁₀ units from the HBM and 0.250 ± 0.038 log₁₀ units from the BIP, indicating that the HBM estimates were significantly more precise than those from the BIP (Wilcoxon signed-rank test¹²², P < 0.001). The depressed CSs at mid-low SF (blue curve in Fig. 4c) were produced by the BIP due to limited data and an uninformative prior. These BIP estimates, with large SDs, were unreliable. On the other hand, the HBM (orange curve in Fig. 4c) generated much smoother and more precise estimates. 

Figure 4.

 

(a through c) Estimated CSFs of a typical subject from the HBM (orange) and BIP (blue) in the three luminance conditions. (d through f) Average SDs of the estimated CSs (HBM: orange; BIP: blue) of the included subjects in the three luminance conditions.

Figure 4d–f depict the average SDs of the estimated CSs of the included subjects in the three luminance conditions from the two methods. Across all the subjects, SFs and luminance conditions, the average ± standard error (SE) of SDs were 0.103 ± 0.001 log₁₀ units from the HBM estimates and 0.252 ± 0.004 log₁₀ units from the BIP estimates. Again, the HBM estimates were significantly more precise (Wilcoxon signed-rank test, P < 0.001). An explanation of how the two methods led to different CS priors and posteriors at each SF in each test can be found in Supplementary Material B. 

Test-Retest Reliability in the H Condition

The 95% coefficient of repeatability (CoR)^123,124 and the test-retest reliability coefficient (TRRC; Pearson's correlation coefficient) both quantify the reliability of the test and retest scores across subjects. Smaller CoR and larger TRRC represent higher test-retest reliability. We assessed the CoR of the estimated CSs of the included subjects in the H condition from the HBM and BIP. Across SFs, the average CoR of the estimated CSs were 0.20 ± 0.03 and 0.74 ± 0.07 log₁₀ units from the HBM and BIP, respectively, with corresponding TRRCs of 0.72 ± 0.03 and 0.12 ± 0.07 (Fig. 5). The HBM reduced the average CoR by a factor of 3.7. The improvement was highly significant (Wilcoxon signed-rank test, P = 0.001), as was the improvement on TRRC (Wilcoxon signed-rank test, P = 0.001). 

Figure 5.

 

Bland-Altman plots of the estimated CSs from repeated measures in the H condition at (a) 1, (b) 1.5, (c) 3, (d) 6, (e) 12, and (f) 18 cpd from the HBM (orange) and BIP (blue). The dotted-dashed lines and dashed lines represent mean CS difference and ±  CoR, respectively.

Figure 6 shows Bland-Altman plots of the estimated AULCSFs in the two repeated tests. CoR of the estimated AULCSFs were 0.112 and 0.355 log₁₀ units for the HBM and BIP, respectively, with corresponding TRRCs of 0.85 and 0.27. Again, the HBM improved the CoR of the estimated AULCSFs by a factor of 3.2. 

Figure 6.

 

Bland-Altman plots of the estimated AULCSFs from repeated measures in the H condition from the (a) HBM and (b) BIP. The dotted-dashed lines and dashed lines represent mean AULCSF difference and ±  CoR, respectively.

Subject-Level Sensitivity at 95% Specificity and Discrimination Accuracy

Figure 7 illustrates the baseline and treatment distributions of the estimated CSs at 3 cpd for a typical subject in the three luminance condition pairs. The treatment distribution is the distribution of the difference between estimated CSs in two different luminance conditions, whereas the baseline distribution is constructed by shifting the treatment distribution to the origin (i.e., eliminating the mean difference between CSs in the two luminance conditions). In all three condition pairs, the baseline and treatment distributions from the HBM are narrower and more separated from each other than those from the BIP, leading to higher sensitivity and accuracy in detecting CS changes. Figures 8 and 9 depict the mean and standard error of the subject-level sensitivity at 95% specificity and discrimination accuracy from the HBM and BIP in the three condition pairs, respectively. 

Figure 7.

 

Illustration of baseline (dashed lines) and treatment (solid lines) distributions from the HBM (orange) and BIP (blue) of the estimated CSs at 3 cpd for a subject between (a) L and M, (b) L and H, and (c) M and H conditions.

Figure 8.

 

The mean and standard error (error bars) of the subject-level sensitivity at 95% specificity from the HBM (orange) and BIP (blue) as a function of spatial frequency between the (a) L and M, (b) L and H, and (c) M and H conditions.

Figure 9.

 

The mean and standard error (error bars) of the subject-level discrimination accuracy from the HBM (orange) and BIP (blue) as a function of spatial frequency between the (a) L and M, (b) L and H, and (c) M and H conditions.

Across all included subjects and SFs, the average sensitivity at 95% specificity were 45.9 ± 0.7% from the HBM and 33.2 ± 0.7% from the BIP, respectively, with corresponding average accuracy of 80.3 ± 0.3% and 75.8 ± 0.3%. The CS estimates from the HBM enabled significantly more sensitive and accurate detections of CS changes between luminance conditions than those from the BIP in all except sensitivity and accuracy at 1 and 18 cpd between L and M, and sensitivity and accuracy at 1, 1.5, and 3 cpd between M and H, in which the two methods were not significantly different (Wilcoxon signed-rank test, α = 0.05, with Bonferroni correction). The nonsignificant difference between the HBM and BIP occurred because the very few trials at 1 cpd caused highly variable CS estimates from both the HBM and BIP, estimated CSs at 18 cpd were near the floor at 0 log₁₀ units, and CSF difference between the M and H conditions was small.⁴⁶ 

Group-Level P Values as a Function of Sample Size

At the group level, the HBM captured strong correlations between SFs in the difference distributions, whereas the BIP could not. Using the covariance between CSs at different SFs, the HBM improved statistical inference at individual SFs and combinations of multiple SFs (Figs. 10a–c) relative to the BIP (Figs. 10d–f). 

Figure 10.

 

Average P values with one, two, three, four, five, and six SFs as functions of sample size between the (a, d) L and M, (b, e) L and H, and (c, f) M and H conditions from the (a–c) HBM and (d–f) BIP.

Figure 10 shows the average P values across different combinations of one to six SFs. P values for discriminating CSs between pairs of luminance conditions decreased with increasing sample size, number of SFs, and number of trials (Supplementary Materials B and C). P values of all the SF combinations from the HBM and BIP are presented in Supplementary Materials B. 

Figure 11 shows the average number of subjects needed to detect a significant CS change with one, two, three, four, five, and six SFs from the HBM and BIP. To detect a significant CS change at α = 0.05, on average, the HBM needed 74% fewer subjects than the BIP in all pairs of luminance conditions across all SF combinations. The results suggest that (1) the statistical power in detecting luminance effects on the CSF increases with the number of SFs included in the statistical test, and (2) the HBM significantly improves statistical power for detecting CS changes compared to the BIP. 

Figure 11.

 

Average number of subjects needed to detect a significant CS change (α = 0.05) with one, two, three, four, five, and six SFs from the HBM (orange) and BIP (blue) between the (a) L and M, (b) L and H, and (c) M and H conditions.

Discussion

To enable advanced statistical inference on CS both at individual SFs and across multiple SFs in clinical trials, we addressed the challenge of obtaining precise CS estimates with limited data at individual SFs by developing a nonparametric HBM. With both 25 and 50 qCSF trials, the HBM generated significantly more precise CS estimates (average SDs = 0.103 ± 0.001 and 0.075 ± 0.000 log₁₀ units) than the BIP (average SDs = 0.252 ± 0.004 and 0.137 ± 0.002 log₁₀ units). The HBM generated significantly more precise CS estimates because it provided informative priors (Fig. 2b) based on relationships between CSs across all SFs, tests, and individuals. Most importantly, the joint posterior distribution of CS from the HBM enabled advanced statistical inference both at individual SFs and across multiple SFs, significantly enhancing the statistical power to detect luminance effects on the CSF compared to inference based on the traditional BIP. 

The covariance hyperparameters at the population and individual levels of the HBM quantified the between- and within-individual CS correlations across SFs. The large positive correlations at the population level in the HBM are consistent with previous findings that, across subjects, better CS at one SF was associated with better CS at other SFs.^42,43 The positive correlations between CSs at similar SFs at the individual and test levels are consistent with the idea of SF channels underlying the CSF, each selective for a range of SFs.^117–119 These covariances were determined by the trial-by-trial data. In contrast, the correlations between CSs from the parametric CSF model are influenced by their deterministic relationships. 

In the H condition, the average CoR of CS estimated from the HBM across SFs were 0.20 and 0.17 log₁₀ units with 25 and 50 qCSF trials, respectively. In contrast, the average CoR of estimated CS across SFs were 0.74 and 0.47 log₁₀ units from the BIP, and 0.28 and 0.19 log₁₀ units in the original study⁴⁶ using the parametric CSF model with 25 and 50 qCSF trials, respectively. The nonparametric HBM reduced the average CoR across SFs by a factor of 3.7 and 2.7 compared to the BIP, and a factor of 1.4 and 1.1 compared to the original study with 25 and 50 qCSF trials, respectively. Similar results were obtained for the AULCSF. 

The HBM significantly improved the average sensitivity at 95% specificity (25 qCSF trials: 45.9 ± 0.7% vs. 33.2 ± 0.7% in BIP; 50 qCSF trials: 58.7 ± 0.8% vs. 47.2 ± 0.8% in BIP) and accuracy (25 qCSF trials: 80.3 ± 0.3% vs. 75.8 ± 0.3% in BIP; 50 qCSF trials: 85.4 ± 0.3% and 81.6 ± 0.3% in BIP) in detecting CS changes for individual subjects, as well as statistical power for detecting group-level CS changes at individual SFs and combinations of multiple SFs between pairs of luminance conditions. The reduction in test-retest variability by the HBM is indicative of a significant decrease in the SDs of the measurements. Consequently, the more precise measurements led to higher sensitivity, increased discrimination accuracy, and greater statistical power when we compared CSs under different luminance conditions in this study. These results carry significant implications for diagnosing and staging diseases with frequency-specific CS deficits and the treatment outcome evaluation. 

The HBM improved CS estimates with informative priors by quantifying relationships between CSs across SFs, tests, and individuals in the joint posterior distribution. We demonstrated the utility of the joint posterior distribution for statistical inference at individual SFs and combinations of multiple SFs at the group level after unblinding the data. The baseline and treatment distributions, P(x|G_baseline) and P(x|G_treatment), derived from the joint posterior distribution, accounted for the correlations across subjects between conditions that were not explicitly modeled in the HBM. We defined the P value as the probability of samples from the baseline distribution (difference distribution without change between test and retest) being incorrectly identified as from the treatment distribution. The definition extends the concept of P value in frequentist tests to high-dimensional distributions. The smaller the P value, the more separated the baseline and treatment distributions are and thus the higher the statistical power. We found that the P value for detecting a change of CS between pairs of luminance conditions decreased with increasing sample size, number of SFs, and number of trials. 

We fixed the slope of the psychometric functions across all SFs and subjects in both the HBM and the BIP based on extensive investigations. Initially, we adopted a fixed slope in the qCSF test,³⁵ drawing from findings in the CSF literature^42,43,125 and a simulation study indicating that fixing the slope of the psychometric function primarily affected the variance rather than the mean of threshold estimates.¹²⁶ In several additional studies, we demonstrated that psychometric functions with a fixed slope provided statistically equivalent fits to CSF data compared to those with free-to-vary slopes across SFs.^{51,53,127,128} To further validate the fixed slope assumption, we¹²⁷ conducted a simulation study in which we measured the CSF of a simulated observer with a fixed psychometric function slope of 1.55 across SFs, using the qCSF test that assumed four different fixed slopes (1, 1.55, 2, and 3.5). Our results showed that the estimated CSFs and their precisions overlapped almost completely, indicating that the specific value of the fixed slope did not significantly impact the accuracy and precision of CSF measurements. 

In this study, the HBM was constructed without explicit modeling of experimental conditions. Modeling relationships of CSF between experimental conditions with additional covariance hyperparameters may further improve the statistical power.¹⁰³ The HBM can also be extended to estimate CSs from tests at multiple time points. Covariance hyperparameters can be used to capture CS relationships across multiple time points to improve statistical inference on treatment efficacy at individual and combinations of multiple SFs. Moreover, the HBM can be extended to explicitly model CS relationships not only at different time points but also different treatment arms with additional covariance hyperparameters. The joint posterior distribution from the HBM can be used to make predictions of CS without any data because of the covariance at the population and individual levels (Supplementary Fig. S8 in Supplementary Materials B). In a related study, we showed that modeling the relationship of CSF parameters between experimental conditions with covariance hyperparameters allows the HBM to make accurate predictions of CSFs in yet-to-be-measured conditions.¹²⁹ 

Additionally, the HBM can be applied to model data in other test modalities such as VA,^104,130,131 and model data across different modalities simultaneously via hierarchical Bayesian joint modeling (HBJM). Although we have only demonstrated the use of HBJM to derive a CE with VA and CSF assessments,¹⁰⁷ the HBJM can combine data from multiple tests and modalities, such as CS and retina thickness measured by optical coherence tomography.^132,133 

The hierarchical adaptive design optimization (HADO)^134,135 framework was developed to improve test efficiency with informative priors for new patients based on their group membership and first applied to the qCSF test. Instead of using group membership to generate informative priors, the HBM can generate informative priors for CS tests (1) at untested SFs for patients already tested at other SF(s) based on the within-individual relationship between CS across SFs quantified in the covariance at the individual level, and (2) for new patients based on the between-individual relationship across SFs quantified in the covariance at the population level.¹²⁹ 

The HBM facilitates statistical inference at individual and combinations of multiple SFs, leading to improved statistical power. This holds particular importance in ophthalmic trials, given the diverse clinical conditions that can affect CSs within specific SF ranges. Conditions such as amblyopia,^20,21 cataracts,^14,15 diabetic retinopathy,^17,18 glaucoma,^10,61,62 idiopathic epiretinal membrane,¹⁹ inherited retinal degeneration,²³ keratoconus,⁶⁴ maculopathies,^{11–13,22,24,65} multiple sclerosis,^25,26 and myopia,^27–31 may manifest as CS deficits at different SFs. Moreover, treatments for various eye diseases may exhibit differential efficacies at different SFs.¹³⁶ 

Despite the acknowledged impact of SF-specific CS deficits, most studies have traditionally computed CSs at individual SFs using parametric CSF models. The introduction of the HBM in this paper addresses this limitation, offering a valuable tool for enhancing statistical inference at both individual SFs and combinations of multiple SFs. Notably, the HBM's contribution extends beyond statistical analysis; it has the potential to deepen our understanding of how different diseases and various stages within a single disease influence CSs across SFs. This insight may prove instrumental in refining diagnostic approaches and disease staging based on CSF deficits. 

Studies have consistently demonstrated a strong correlation between CS and quality of life measures in patients with eye diseases as well as decline of the entire CSF in aging.^22,54,68,137 However, the correlations between CSs derived from the parametric CSF model across SFs complicate statistical inference, and prior assessments using CS at a single SF or AULCSF lacked specificity in establishing the relationship between CS deficits and eye diseases. The HBM introduces a promising avenue for improving this specificity, providing a deeper understanding of the intricate interplay between various diseases, different disease stages, and their impact on CSs across SFs. This refined perspective on CS deficits could significantly contribute to the diagnosis and staging of eye diseases. 

As a general framework, the HBM can enhance statistical inference by incorporating relationships across tests, experimental conditions, and individuals. In this study, we demonstrated its utility in estimating CS from limited qCSF data. This framework holds promise for enhancing statistical inference in other structural and functional vision tests in both basic and clinical studies with hierarchical designs or involving multiple conditions. Further research is necessary to explore its potential in various context. 

Moreover, the assumption that hyperparameters can be precisely captured using Gaussian mixtures may not always hold true, particularly for populations with more complicated statistical properties. This assumption could introduce inaccuracies in estimating the posterior distribution and warrants more comprehensive investigation in future studies. 

In conclusion, the HBM improved CS estimates with informative priors by incorporating information across SFs, tests, and individuals in the joint posterior distribution of CS parameters and hyperparameters at the test, individual, and population levels, and enabled evaluation of statistically significant CS differences at both subject and group levels at individual SFs and across multiple SFs with limited data. It provides a general framework to improve statistical inference in multi-condition hierarchical experiment designs. 

Acknowledgments

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

Disclosure: Y. Zhao, None; L.A. Lesmes, Adaptive Sensory Technology, Inc. (E, F, I); M. Dorr, Adaptive Sensory Technology, Inc. (E, F, I); Z.-L. Lu, Adaptive Sensory Technology, Inc. (F, I), Jiangsu Juehua Medical Technology, LTD (F, I)