**Purpose**:
The goal of this study is to develop a hierarchical Bayesian model (HBM) to better quantify uncertainty in visual acuity (VA) tests by incorporating the relationship between VA threshold and range across multiple individuals and tests.

**Methods**:
The three-level HBM consisted of multiple two-dimensional Gaussian distributions of hyperparameters and parameters of the VA behavioral function (VABF) at the population, individual, and test levels. The model was applied to a dataset of quantitative VA (qVA) assessments of 14 eyes in 4 Bangerter foil conditions. We quantified uncertainties of the estimated VABF parameters (VA threshold and range) from the HBM and compared them with those from the qVA.

**Results**:
The HBM recovered covariances between VABF parameters and provided better fits to the data than the qVA. It reduced the uncertainty of their estimates by 4.2% to 45.8%. The reduction of uncertainty, on average, resulted in 3 fewer rows needed to reach a 95% accuracy in detecting a 0.15 logMAR change of VA threshold or both parameters than the qVA.

**Conclusions**:
The HBM utilized knowledge across individuals and tests in a single model and provided better quantification of the uncertainty of the estimated VABF, especially when the number of tested rows was relatively small.

**Translational Relevance**:
The HBM can increase the accuracy in detecting VA changes. Further research is necessary to evaluate its potential in clinical populations.

^{1}Precise assessment of VA is of paramount importance for accurate diagnosis of eye diseases, monitoring disease progression, and evaluating treatment efficacy, as well as setting classification and qualification standards in many sports and professions.

^{2}

^{–}

^{8}However, VA measurements are intrinsically noisy and the estimated VA scores carry a lot of uncertainty. This is because human VA behavior is probabilistic in nature.

^{9}

^{–}

^{12}Given the exact same optotypes, the responses from a human observer may vary from trial to trial because of the various noises in the perceptual system.

^{13}

^{–}

^{16}In this study, we apply the Hierarchical Bayesian Modeling (HBM) approach to improve quantification of uncertainty in VA assessments.

^{17}

^{–}

^{20}Although traditional VA assessments have only focused on estimating VA threshold, both VABF parameters may vary across individuals and disease stages.

^{18}

^{,}

^{21}

^{–}

^{23}A complete characterization of the VABF requires a joint estimate of both parameters.

^{24}

^{,}

^{25}

^{26}Given the results from a single VA assessment (Fig. 1D), the probability distribution of VABF can be estimated by their likelihood or Bayesian inference (Fig. 1E). The corresponding parameters and the associated probabilities form a joint distribution of the estimated VABF parameters (Fig. 1F). The mean and spread of the joint distribution can be used to quantify the estimated parameters and their uncertainties,

^{26}

^{–}

^{30}which is often quantified by the half width of the 68.2% credible interval (68.2% HWCI) of the distribution,

^{26}

^{,}

^{27}defined as the smallest interval that contains the true value of the estimated quantity with 68.2% probability. The estimated uncertainties from repeated tests (see Fig. 1C) and a single measurement (see Fig. 1F) are equivalent under ideal experimental conditions.

^{26}

^{,}

^{30}

^{,}

^{31}

^{32}

^{–}

^{34}The method of limits and heuristic termination rules are used to determine a single VA score (but see Ref. 35). Although uncertainty is usually not estimated, studies with repeated tests found large uncertainties (approximately 0.07–0.34 logMAR)

^{33}

^{,}

^{36}

^{–}

^{43}because of the coarse sampling of the VABF as well as the use of heuristic termination rules that depend on the probabilistic behavior of the human observer. In addition, the exact performance level (probability of correct optotype identification) associated with each heuristic termination rule depends on the range parameter of the VABF.

^{24}

^{,}

^{25}We cannot meaningfully compare VA scores from patients with different VA range parameters or interpret changes of VA scores of the same patient with range parameter changing between assessments.

^{18}

^{,}

^{19}

^{24}

^{,}

^{25}is a Bayesian inference procedure developed to characterize the full VABF and quantify the uncertainty in each assessment of individual

*i*in test

*j*. In the qVA, the VABF (see Fig. 1) is characterized with parameter θ

_{ij}= (α

_{ij}, β

_{ij}), where α

_{ij}is VA threshold, corresponding to the \(d^{\prime}\) = 2 performance level, β

_{ij}is the range parameter of the function that covers \(d^{\prime}\) = 1 to \(d^{\prime}\) = 4 performance levels.

^{24}

^{,}

^{25}The qVA incorporates high density optotype size sampling and starts with a broad joint prior distribution of θ

_{ij}, representing existing knowledge of the general pattern of the VABF (Fig. 2A). An active learning approach is used to optimize stimulus selection to reduce the expected uncertainty of the posterior distribution of θ

_{ij}in each trial (see Fig. 2). Bayes’ rule is used to compute the joint posterior distribution of θ

_{ij }(Figs. 2B, 2C), which allows us to quantify not only the VABF parameters but also their uncertainties from a single measurement. In computer simulations, we showed that the qVA could assess VABF parameters with virtually no bias, and very small uncertainty in α

_{ij}(HWCI = 0.028 logMAR) and small uncertainty in β

_{ij}(HWCI = 0.092 logMAR), reflecting 52.5% and 49.5% uncertainty reductions of the estimated α

_{ij}relative to the electronic-early treatment diabetic retinopathy study (E-ETDRS) and Freiburg Visual Acuity and Contrast Test (FrACT) methods, respectively.

^{25}The results were confirmed in a psychophysical study: estimated θ

_{ij}from the qVA exhibited very small (α

_{ij}= 0.019 logMAR HWCI) and relatively small (β

_{ij}= 0.062 logMAR HWCI) uncertainty. In addition, we found a significant correlation (r = 0.412,

*P*< 0.001) between estimated α

_{ij}and β

_{ij}across individuals and tests in the psychophysical experiment.

^{25}

_{ij}and β

_{ij},

^{24}

^{,}

^{25}which could be useful in further constraining the estimated θ

_{ij}s and thereby reducing their uncertainties. In other words, the broad prior and independent treatment of each qVA test may overestimate the uncertainty of the estimated VABF. Although the qVA has greatly reduced the uncertainty of α

_{ij}estimates, the uncertainty of β

_{ij}estimates is still relatively large.

^{25}

_{ij}and β

_{ij}across individuals and tests in our previous study

^{25}motivated us to develop a three-level HBM (Fig. 3B) to quantify uncertainties at the population, individual, and test levels in a single model. At the population level,

*p*(η), the joint distribution of VABF hyperparameter η across all the individuals, is modeled as a mixture of 2-dimensional Gaussian distributions with mean µ and covariance Σ, which have distributions

*p*(µ) and

*p*(Σ). At the individual level,

*p*(τ

_{i}|η), the joint distribution of VABF hyperparameter τ

_{i}of individual

*i*across all the tests performed by the individual, is modeled as a mixture of 2-dimensional Gaussian distributions with mean ρ

_{i}and covariance φ, which have distributions

*p*(ρ

_{i}|η) and

*p*(φ).

*p*(ρ

_{i}|η) denotes that ρ

_{i}is conditioned on η. At the test level,

*p*(θ

_{ij}|τ

_{i}), the joint distribution of the VABF parameters of individual

*i*in test

*j*, θ

_{ij}, is conditioned on τ

_{i}. The cross- and within-individual regularities are modeled as covariances Σ and φ in the HBM.

^{44}

^{–}

^{48}Bayes’ rule is then used to update the joint prior distribution of all the parameters and hyperparameters to their joint posterior distribution,

^{49}

^{–}

^{53}which allows us to quantify the estimated hyperparameters and parameters as well as their uncertainties at the population, individual and test levels.

^{49}

^{,}

^{54}

^{,}

^{55}

^{56}ecology,

^{57}

^{,}

^{58}genetics,

^{59}and cognitive science.

^{50}

^{,}

^{51}

^{,}

^{53}

^{,}

^{60}

^{–}

^{65}By decomposing the variability of an entire dataset into distributions at multiple levels of the hierarchy, it can better quantify uncertainty at each level.

^{59}

^{,}

^{52}

^{,}

^{65}Here, we developed an HBM and evaluated its performance relative to that of the qVA using the qVA data in 14 eyes tested in 4 Bangerter foil conditions.

^{25}Our goal was to improve quantification of uncertainty at the test level. We hypothesized that the HBM would account for the data better than the qVA and reduce the uncertainties of estimated VABF parameters in single VA tests relative to the qVA.

^{66}

^{,}

^{67}and without foil.

^{25}Each qVA test consisted of 45 rows (

*K*= 45) of 3 optotypes (Fig. 4). In each qVA trial, 3 letters on a row, randomly sampled without replacement from the 10 Sloan letters (C, D, H, K, N, O, R, S, V, and Z) and with their size determined by the qVA, were presented to the subject. The subject was instructed to report the identity of the letters verbally. The joint posterior distributions of VABF parameters in each qVA test after

*k*= 5, 10, 15, 20, 25, 30, 35, 40, and 45 rows were used to conduct statistical analysis. The details of the qVA can be found in prior publications.

^{24}

^{,}

^{25}

*i*’s response

*r*(the number of correct identifications: 0, 1, 2, or 3) of the 3 optotypes with optotype size

_{ijk}*os*in row

_{ijk}*k*of test

*j*:

*os*,θ

_{ijk}_{ij}) is the VABF of correct identification of a single optotype with size

*os*and model parameters θ

_{ijk}_{ij}= (α

_{ij}, β

_{ij});

*f*(g(

*os*, θ

_{ijk}_{ij}),

*n*) is the probability of observing the number of correct identifications given g(

*os*, θ

_{ijk}_{ij}) and a specific chart design (e.g.

*n*= 3 optotypes on a row). The details of functions

*f*and

*g*can be found in Supplementary Materials A and previous publications.

^{24}

^{,}

^{43}The probability of obtaining the observed responses in all rows in test

*j*from individual

*i*is the product of the probability of responses in all rows (Equation 1a):

*p*(η), the joint distribution of hyperparameter η, is modeled as a mixture of two-dimensional Gaussian distributions \({\cal N}\) with mean µ and covariance Σ, which have distributions

*p*(µ) and

*p*(Σ):

*p*(τ

_{i}|η), the joint distribution of hyperparameter τ

_{i}of individual

*i*, is modeled as a mixture of two-dimensional Gaussian distributions \({\cal N}\) with mean ρ

_{i}and covariance φ, which have distributions

*p*(ρ

_{i}|η) and

*p*(φ):

*p*(ρ

_{i}|η) denotes that ρ

_{i}is conditioned on η.

*p*(θ

_{ij}|τ

_{i}), the joint distribution of parameter θ

_{ij}is conditioned on τ

_{i}.

*X*= (θ

_{1:I, 1:J}, ρ

_{1:I}, φ, µ, Σ) are all the parameters and hyperparameters in the HBM,

*I*is the total number of individuals, and

*J*is the total number of tests on each individual.

_{0,min}= (−0.5, 0.1) logMAR and µ

_{0,max}= (1.3, 1.5) logMAR; precision matrices φ

^{−1}and Σ

^{−1}are the inverse of covariance φ and Σ; \({\cal W}( {\frac{{\rm{Y}}}{{\rm{\nu }}},{\rm{\nu }}} )\) denotes a Wishart distribution with expected mean precision matrix Y and degrees of freedom ν = 2; \(\phi _{qVA}^{ - 1}\) is the inverse of the average covariance matrix φ

_{qVA}computed from the joint posterior distributions across all qVA tests in the dataset; Σ

_{qVA}

^{−1}, is the inverse of the covariance matrix Σ

_{qVA}computed from the estimated θ

_{ij}across all qVA tests.

*X*and is a constant for a given dataset and HBM.

*I = 56*. In addition, each individual was only tested once (

*J = 1*).

^{68}package in R

^{69}to evaluate the joint posterior distribution (see Equation 6). JAGS uses an efficient sampling algorithm to generate representative samples of the joint posterior distribution of all the parameters and hyperparameters in the HBM via Markov Chain Monte Carlo (MCMC). Each MCMC generated 15,000 samples via a random walk process (see Supplementary Materials B for details). Because the initial part of the random walk was largely influenced by the arbitrary starting position, steps in the burn-in and adaptation (further optimization of the sampling algorithm by JAGS) phases were discarded and not included in the analysis. The exact number of steps discarded differ for each model. In this study, 20,000 and 100,000 steps were used for burn-in and adaptation in each HBM model based on results from pilot studies. Convergence of each parameter was evaluated with Gelman and Rubin's diagnostic rule

^{70}based on the ratio of between- and within-MCMC variances along each dimension, that is, the variance of the samples across MCMC processes divided by the variance of the samples in each MCMC process. The HBM was deemed “converged” when the ratios for all the parameters were smaller than 1.05. The data with

*k*= 5, 10, 15, 20, 25, 30, 35, 40, and 45 rows were fit by the HBM separately.

_{i1}obtained from the HBM fit to the real data. Each qVA test consisted of 45 three-optotype rows, and the row-by-row responses were determined by the VABF specified with the parameters of the simulated observer. The HBM was fit to the simulated dataset of 56 tests. The procedure was repeated 10 times.

^{71}

^{,}

^{72}was used to quantify and compare the goodness-of-fit of the HBM against that of the qVA. The BPIC quantifies the likelihood of the data based on the joint posterior distribution of the parameters of a model and penalizes model complexity.

_{i1}from the qVA and HBM. The uncertainty of each estimated parameter was quantified by the 68.2% HWCI of its marginal posterior distribution.

^{26}

^{,}

^{27}Paired

*t*-test (R function

*t.test*

^{69}) was used to compare the 68.2% HWCIs between the two methods.

_{i1}) change, a 0.15 logMAR range (β

_{i1}) change, and a 0.15 logMAR change of both parameters. We chose 0.15 logMAR as the magnitude of α

_{i1}change because a greater-than-15-letter VA improvement (0.3 logMAR) is considered by the US Food & Drug Administration (FDA) as an acceptable end point of a clinical trial, although a greater-than-10-letter VA improvement (0.2 logMAR) has been used when the benefits can outweigh the safety risks of the proposed method or product.

^{73}We chose 0.15 logMAR as the magnitude of β

_{ij}change, because a similar change has been reported in eyes with degraded vision.

^{18}

^{,}

^{25}

^{74}Whereas sensitivity and specificity depend on the change criterion, the AUROC in signal detection theory provides a criterion-free estimate of the accuracy of a test to detect a change.

^{74}It only depends on the discriminability \(d^{\prime}\) that weighs the evidence against uncertainty and quantifies the signal (change) to noise (HWCI) ratio of two probability distributions of the estimates in two conditions.

^{74}:

^{75}:

*cov*(Δ) is the covariance matrix of the VABF parameter difference distribution,

*cov*(Δ)

^{−1}is the inverse of

*cov*(Δ), Δ

^{T}is the transpose of Δ, * represents matrix multiplication, and

*corr*is the correlation coefficient at the test level in the HBM.

_{1: I, 1}, it did not significantly change the expected values of the estimated θ

_{1:I, 1}: the average absolute differences of the estimates from the qVA and HBM across all individuals were 0.002 and 0.022 logMAR, respectively, indicating that HBM and qVA estimates exhibited virtually the same central tendency.

^{1}and η

^{2}was 0.662. (

*P*< 0.001). Figure 7 illustrates the joint posterior distributions of hyperparameters φ, and ρ

_{i}, and τ

_{i}for 8 out of the 56 individuals after 45 rows. Whereas ρ

_{i}and φ were from the HBM fits to the data (see Equation 6), τ

_{i}was constructed from the distributions of ρ

_{i}and φ (see Equation 3). Figure 8A illustrates the joint posterior distributions of parameters θ

_{i1}from the 8 corresponding individuals in Figure 7 after 45 rows. Figure 8B shows the mean θ

_{i1}of all the 56 tests.

_{i1}Estimates

_{i1}after 5, 15, and 45 rows for one individual, obtained from both the qVA and HBM. Across all the individuals and tests, the average 68.2% HWCI of the estimated θ

_{i1}decreased with the number of rows in both methods (Fig. 10, Table 1). With 5 rows, the 68.2% HWCIs of the estimated α

_{i1}(t(55) = 14.8,

*P*< 0.001) and β

_{i1}(t(55) = 13.7,

*P*< 0.001) were significantly different between the methods. With 45 rows, the 68.2% HWCIs of the estimated α

_{i1}(t(55) = 4.63,

*P*< 0.001) and β

_{i1}(t(55) = 6.74,

*P*< 0.001) were significantly different between the methods. Relative to the qVA, the HBM reduced the HWCI of the estimated α

_{i1}by 4.2% to 26.9%, and β

_{i1}by 20.8% to 45.8%.

_{i1}. Because the trial-by-trial data and distribution of τ

_{i}were different for each individual and test, the uncertainties were different across individuals and tests in the HBM (Fig. 11).

_{i1}(−0.0024 logMAR) and β

_{i1}(0.0057 logMAR) across all the simulations. The average root–mean-square errors (RMSEs) between the HBM estimates and the true values (0.019 and 0.047 logMAR for α

_{i1}and β

_{i1}) were comparable to the HWCIs of the parameter estimates (0.019 and 0.050 logMAR for α

_{i1}and β

_{i1}), indicating good convergence of the HBM fits. In addition, the HWCIs from the simulations were also comparable to those from the experimental data (0.018 and 0.048 logMAR for α

_{i1}and β

_{i1}), suggesting that the simulations captured the uncertainties in the real experiment.

_{i1}change with the qVA and HBM, and 10 and 6 rows were needed to detect a 0.15 logMAR change of both α

_{i1}and β

_{i1}, respectively.

_{i1}estimates by 26.9%, 12.2%, and 4.2%, and β

_{i1}estimates by 45.8%, 35.3%, and 20.8% for tests with 5, 15, and 45 rows, respectively. Simulations showed that the HBM recovered α

_{i1}spanning about 0.8 logMAR (−0.15 to approximately 0.68 logMAR) and β

_{i1}spanning about 0.4 logMAR (approximately 0.11 to 0.49 logMAR). The simulation results also provided validity for the discrimination accuracy analysis.

_{i1}and improved the discrimination accuracy of the tests. The reduction of uncertainty, on average, resulted in 3 fewer rows needed to reach 95% accuracy in detecting changes of 0.15 logMAR in α

_{i1}or both α

_{i1}and β

_{i1}than the qVA (see Table 2, Fig. 12). The increased accuracy in detecting β

_{i1}changes in the HBM may enable a new range endpoint in functional vision,

^{5}because it has been documented that poor VA thresholds were accompanied by wider VA ranges.

^{25}In addition, using the joint distribution of θ

_{i1}can further increase the accuracy in detecting changes of VA behavior because the two-dimensional distributions may contain more information than the one-dimensional marginal distributions.

^{76}

^{–}

^{78}

_{ij}of all the tests conditioned on τ

_{i}.

_{i1}as the original HBM fit to the 14 eyes, with slightly higher HWCIs of the estimated β

_{i,1}(0.018 ± 0.00082 and 0.050 ± 0.00184 logMAR for α

_{i1}and β

_{i1}in the 7 eye fits, versus 0.018 and 0.048 logMAR in the 14 eye fit). The RMSEs between the HBM fits to the 7 and 14 eyes were 0.00078 and 0.00793 logMAR for α

_{i1}and β

_{i1}, and 0.00036 and 0.00350 logMAR for the corresponding HWCIs. The results also suggest that the original experimental dataset with 14 eyes tested in 4 conditions was sufficient to constrain the HBM because cutting the amount of data to half generated virtually the same estimates.

_{i}at the population and individual levels from the HBM can also be used to construct informative priors (see Figs. 6C, 7C) within the hierarchical adaptive design optimization (HADO) framework

^{79}

^{,}

^{80}for new individuals (see Fig. 6C) and repeated tests of the same individual (see Fig. 7C), respectively. In the HADO framework, the joint distribution in the HBM is updated after each test and used as an informative prior in the next test to improve test efficiency for new individuals or repeated tests of the same individual. The 6-hour computation time of the HBM on a desktop computer makes it realistic to incorporate it within the HADO framework with daily updates of the joint posterior distribution.

**Y. Zhao**, None;

**L.A. Lesmes**, (I, E, P);

**M. Dorr**, (I, E, P);

**Z.-L. Lu**, (I, P)

*Ophthalmol Clin North Am*. 2003; 16: 155–170. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2010; 94: 400–405. [CrossRef] [PubMed]

*Eye.*2003; 17(5): 579–582. [CrossRef]

*Clinical Methods: The History, Physical, and Laboratory Examinations*(3rd ed.). Boston, US: Butterworths; 1990.

*Br Med J.*2006; 332(7545): 820–824. [CrossRef]

*Invest Ophthalmol Vis Sci.*1991; 32: 422–432. [PubMed]

*Optom Vis Sci.*2001; 78: 529–538. [CrossRef] [PubMed]

*Ophthalmic Physiol Opt.*1988; 8: 363–370. [CrossRef] [PubMed]

*Optom Vis Sci.*1998; 75: 342–348. [CrossRef] [PubMed]

*J Optical Society of Am*. 1956; 46: 634–639. [CrossRef]

*Nat Rev Neurosci*. 2008; 9: 292–303. [CrossRef] [PubMed]

*Psychol. Rev.*2008; 115: 44–82. [CrossRef] [PubMed]

*Psychon Bull Rev*. 2010; 17: 802–808. [CrossRef] [PubMed]

*Vision Res*. 1997; 37: 813–819. [CrossRef] [PubMed]

*Optom Vis Sci*. 2001; 78(2): 113–121. [CrossRef] [PubMed]

*J Optom*. 1985; 62: 895–900. [CrossRef]

*Percept Psychophys.*2001; 63: 1348–1355. [CrossRef] [PubMed]

*Graefes Arch Clin Exp Ophthalmol*. 2007; 245: 965–971. [CrossRef] [PubMed]

*Ophthalmic Physiol Opt*. 2017; 37(2): 118–127. [CrossRef] [PubMed]

*Optom Vis Sci*. 2002; 79(12): 788–792. [CrossRef] [PubMed]

*Proceedings of the 2nd International Conference on Applications of Intelligent Systems*. APPIS ’19. New York, NY: Association for Computing Machinery; 2019: 1–6.

*Trans Vis Sci Tech.*2021; 10: 1. [CrossRef]

*Visual Psychophysics: From Laboratory to Theory*. Cambridge, MA: MIT Press; 2013.

*Statistical models in epidemiology*. Oxford, UK: Oxford University Press; 1993.

*Psychol Rev*. 1963; 70(3): 193–242. [CrossRef]

*Percept Psychophys*. 2001; 63: 1293–1313. [CrossRef] [PubMed]

*Percept Psychophys*. 2001; 63: 1314–1329. [CrossRef] [PubMed]

*J Vis*. 2015; 15: 2. [CrossRef] [PubMed]

*Am J Optom Physiol Opt*. 1976; 53: 740–745. [CrossRef] [PubMed]

*Am J Ophthalmol*. 2003; 135: 194–205. [CrossRef] [PubMed]

*Am J Ophthalmol*. 1982; 94: 91–96. [CrossRef] [PubMed]

*Optom Vis Sci*. 1996; 73: 49–53. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 1993; 34: 120–129. [PubMed]

*Eye*. 2015; 29: 1085–109. [CrossRef] [PubMed]

*Ophthalmic Physiol Opt*. 1988; 8: 397–401. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2003; 87(10): 1232–1234. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2008; 92: 241–244. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci.*2008; 49: 4347–4352. [CrossRef] [PubMed]

*Br J Ophthalmol*. 2001; 85: 432–436. [CrossRef] [PubMed]

*Eye*. 1997; 11: 411–417. [CrossRef] [PubMed]

*J Am Stat Assoc*. 1999; 94(448): 1254 [CrossRef]

*Front Psychol*. 2019; 10: 1675. [CrossRef] [PubMed]

*Stat Med.*2003; 22(5): 763–780. [CrossRef] [PubMed]

*Stat Methods Med Res*. 2020; 29(4): 1112–1128. [CrossRef] [PubMed]

*Bayesian Anal*. 2011; 11(3): 649–670.

*Doing Bayesian data analysis: a tutorial with R, JAGS, and Stan*. San Diego, CA: Academic Press; 2014.

*Cogn Sci*. 2006; 30(3): 1–26. [CrossRef] [PubMed]

*J Math Psychol*. 2011; 55(1): 1–7. [CrossRef]

*Psychon Bull Rev*. 2005; 12(4): 573–604 [CrossRef] [PubMed]

*Comput Brain Behav*. 2020; 3: 384–399. [CrossRef]

*Front Psychol*. 2019; 10: 1675. [CrossRef] [PubMed]

*Psychon Bull Rev*. 2008; 15(1): 1–15. [CrossRef] [PubMed]

*Publ Astron Soc Aust*. 2019; 36: e010. [CrossRef]

*Mar Biol*. 2015; 162(2): 469–478. [CrossRef]

*Ecology*. 2003; 84(6): 1382–1394. [CrossRef]

*Evolution*. 2002; 56(1): 154–166. [CrossRef] [PubMed]

*J Neurosci Psychol Econ*. 2011; 4(2): 95–110. [CrossRef] [PubMed]

*J Opt Soc Am A Opt Image Sci Vis*. 2003; 20(7): 1434–1448. [CrossRef] [PubMed]

*J Math Psychol*. 2011; 55(1): 57–67. [CrossRef]

*Comput Brain Behav*. 2018; 1(2): 184–213. [CrossRef]

*J Cogn Neurosci*. 2019; 31(12): 1976–1996. [CrossRef] [PubMed]

*Psychometrika*. 2003; 68(4): 589–606. [CrossRef]

*Invest Ophthalmol Vis Sci*. 2010; 51: 609–613. [CrossRef] [PubMed]

*J AAPOS*. 2008; 12(6): 555–559. [CrossRef] [PubMed]

*Proceedings of the 3rd international workshop on distributed statistical computing*. 2003. Available at: https://www.r-project.org/nosvn/conferences/DSC-2003/Drafts/Plummer.pdf.

*R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing; 2003. Available at: https://www.R-project.org/.

*Stat Sci*. 1992; 7: 457–511.

*Biometrika*. 2007; 94: 443–458. [CrossRef]

*Am J Math Manag Sci*. 2011; 31: 13–38.

*Invest Ophthalmol Vis Sci*. 2017; 58(9): 3456–3463. [CrossRef] [PubMed]

*Signal Detection Theory and Psychophysics*. New York, NY: John Wiley & Sons; 1966.

*Psychol Rev*. 1986; 93(2): 154–179. [CrossRef] [PubMed]

*J Clin Epidemiol*. 2007; 60(12): 1234–1238. [CrossRef] [PubMed]

*BMC Med Res Methodol*. 2014; 14: 49. [CrossRef] [PubMed]

*Modern Statistics for the Social and Behavioral Sciences: A Practical Introduction*(pp. 101–102). Boca Raton, FL: CRC Press; 2012.

*J Vis*. 2016; 16(6): 15. [CrossRef] [PubMed]

*Neural Comput*. 2014; 26(11): 2465–2492. [CrossRef] [PubMed]