Improved Statistical Methods for Quantal Assay

This letter will argue that a recent report in Anesthesiology would have been improved by different statistical methods; this letter is not a criticism of other research methods in the report. Ginosar et al.  1reported a dose–response study of intrathecal hyperbaric bupivacaine (with adjuvant intrathecal opioids) administered for cesarean delivery in patients having a combined spinal–epidural technique. Seven groups of parturients—six patients per group—were randomly assigned to linearly spaced doses of bupivacaine (6, 7, 8, 9, 10, 11, or 12 mg). Success of the block (binary yes/no) was declared at 10 min after administration if a bilateral T6 pinprick sensory level existed (success^induction); success^operationwas recorded if surgery proceeded without administration of any supplemental epidural local anesthetics. The observed response rates with increasing dose were success^induction(3/6, 2/6, 6/6, 6/6, 6/6, 6.6) and success^operation(1/6, 2/6, 4/6, 3/6, 6/6, 6/6, 6/6). Ginosar et al.  used a version of the Hill equation (also known as the quantal sigmoid Emax model) to relate probability of success to dose of bupivacaine with two parameters, γ and dose⁵⁰:

This was further described as logistic regression analysis of naïve-pooled data (one observation per patient) using Laplacian estimation routines of NONMEM (version V) statistical software (NONMEM Project Group, University of California, San Francisco, CA). Estimates of ED⁵⁰with SEs were reported for success^induction(6.7 ± 0.6 mg) and success^operation(7.6 ± 0.4 mg); point estimates of ED⁹⁵without standard errors were also reported (11.0 and 11.2 mg, respectively).

The community of statisticians has produced an extensive repertory of methods for the analysis of quantal response data.2It is usually assumed that each individual of the relevant population has a dose tolerance or threshold for the particular substance being tested; a descriptive model characterizes the distribution of tolerances. It is possible, but not necessary, to use a logarithmic dose transformation. Assuming large sample properties, sigmoidicity, symmetry, and homoscedasticity, the most common simple model is called the logit or logistic regression:

where α and β are location and scale parameters. Further, assuming a logistic probability density function, maximum likelihood estimation routines for α and β by iteratively reweighted linear regression are available in most statistical software packages; using the δ method, ED^xwith standard errors may be estimated for any x (0 < x < 100). The ED⁵⁰represents the median value of the distribution of tolerances in the population. Counts for success and failure being available in their figure 1, the data of Ginosar et al.  1can be reanalyzed assuming linear spacing of doses. Using the open software R statistical computing and graphics package* (version 1.8.1) with the base and MASS libraries, the estimates for success^inductionare ED⁵⁰(6.5 ± 0.4 mg) and ED⁹⁵(8.6 ± 0.7 mg); the estimates for success^operationare ED⁵⁰(7.7 ± 0.4 mg) and ED⁹⁵(10.6 ± 0.9 mg). The poor precision of the ED⁹⁵estimates is apparent; the 95% confidence interval for success^inductionis 7.2–10.0 mg and for success^operationis 8.8–12.4 mg.

Objections to the statistical methods can be briefly summarized. First, the Hill equation originated in studies of multiple ligand binding to allosteric proteins, in particular hemoglobin, the exponent γ (a slope parameter) being interpreted in a mechanistic way to reflect the cooperativity (interaction of ligands) in binding. In quantal assay, the γ parameter is sometimes described as the steepness of the probability of effect curve for an individual patient.3The complexity of general and spinal anesthesia allow considerable skepticism that the single parameter γ can have any deterministic/mechanistic interpretation. The more conservative approach is to consider any slope estimate of the anesthetic dose–response curve as purely descriptive of the distribution of thresholds in the patient population.

Second, as an extension of concepts developed in the mixed effects modeling of population pharmacokinetic data (multiple observations per subject), it has been argued that mixed effects modeling can be used on single response data to estimate an intraindividual and an interindividual variance.4However, a fundamental flaw was demonstrated in assumptions for such methods5: Data with only one observation per subject cannot be used to estimate an intraindividual variance. Thus, the nomenclature “naive-pooled data,” commonly used in population pharmacokinetic modeling, is an incorrect description of the data structure in this experiment.

Third, the NONMEM statistical package is one of the prominent software tools developed for mixed effects modeling. Mathematical calculations within NONMEM such as Laplacian estimation are extremely complex, involving many assumptions and approximations; there is no consensus about the optimal estimation routines for mixed effects software.6Standard logistic regression software uses commonly accepted routines that allow the estimation of confidence intervals for arbitrary ED^x, missing in the NONMEM output. It should be emphasized that even standard logistic regression analysis gives less precision for ED^xvalues at the upper and lower edges of the sigmoid curve.

Fourth, simple algebraic manipulation shows that the quantal sigmoid Emax model can be rewritten in a logistic format:

This restatement reveals that the quantal sigmoid Emax model enforces a logarithmic transformation of dose. Such a transformation may or may not be desirable. The standard logistic model leaves this choice to the modeler. In this experiment, bupivacaine doses were linearly spaced, and a logarithmic transform seems unnecessary. However, using a logarithmic dose transformation, the logistic regression estimates (with 95% confidence intervals) are ED⁵⁰(6.5 [5.8–7.3] mg) and ED⁹⁵(8.7 [7.1–10.5] mg) for success^inductionand ED⁵⁰(7.6 [6.8–8.5] mg) and ED⁹⁵(11.0 [8.8–13.6] mg) for success^operation. The NONMEM estimate for the ED⁹⁵of success^inductionis 11.0 mg, a value extremely discordant with the observed response rates (fig. 1, Ginosar et al.  1); standard logistic regression gives estimates (linear dose, 8.6 mg; logarithmic dose, 8.7 mg) consistent with the observed response rate.

The observed response rates for this experiment were nonmonotonic, decreasing at some intermediate doses. This creates difficulty in reliable statistical estimation. Nonparametric methods of obtaining doses for ED^xare available and could have been considered.2If a simple parametric model were desired, standard logistic regression analysis would have been preferred.

University of Utah, Salt Lake City, Utah. n.l.pace@utah.edu