Pharmacodynamic data frequently consist of the binary assessment (a "yes" or "no" answer) of the response to a defined stimulus (verbal stimulus, intubation, skin incision, and so on) for multiple patients. The concentration-effect relation is usually reported in terms of C50, the drug concentration associated with a 50% probability of drug effect, and a parameter the authors denote gamma, which determines the shape of the concentration-probability of effect curve. Accurate estimation of gamma, a parameter describing the entire curve, is as important as the estimation of C50, a single point on this curve. Pharmacodynamic data usually are analyzed without accounting for interpatient variability. The authors postulated that accounting for interpatient variability would improve the accuracy of estimation of gamma and allow the estimation of C50 variability.

A probit-based model for the individual concentration-response relation was assumed, characterized by two parameters, C50 and gamma. This assumption was validated by comparing probit regression with the more commonly used logistic regression of data from individual patients found in the anesthesiology literature. The model was then extended to analysis of population data by assuming that C50 has a log-normal distribution. Population data were analyzed in terms of three parameters, (C50), the mean value of C50 in the population; omega, the standard deviation of the distribution of the logarithm of C50; and gamma. The statistical characteristics of the technique were assessed using simulated data. The data were generated for a range of gamma and omega values, assuming that C50 and gamma had a log-normal distribution.

The probit-based model describes data from individual patients and logistic regression does. Population analysis using the extended probit model accurately estimated (C50), gamma, and omega for a range of values, despite the fact that the technique accounts for C50 variability but not gamma variability.

A probit-based method of pharmacodynamic analysis of pooled population data facilitates accurate estimation of the concentration-response curve.

Pharmacodynamic data are often recorded as binary variables; that is, the patient does or does not respond to a command, can or cannot maintain adequate spontaneous ventilation, does or does not have a hemodynamic or somatic response to surgical stimulus, and so forth. In this situation, the most common technique of data analysis is logistic regression, in which the probability of drug effect is evaluated as a function of C, the drug concentration in plasma (or at the effect site), with the equation Equation 1

C^{50}is the concentration at which the probability of drug effect is 50% and gamma is a measure of the steepness of the concentration-effect curve (throughout this article P refers to the probability of drug effect, such as the probability of ablating the response to some stimulus, such as skin incision). Logistic regression has been used for the analysis of the pharmacodynamics of inhaled and intravenous anesthetics, [1–10] usually with a primary focus of determining C^{50}values. However, logistic regression not only estimates the median point of the concentration-response curve (C^{50}), but it also provides information on the shape of this curve by estimating the steepness parameter gamma. Note that Equation 1, when applied to a single patient, implies that there is a fundamental element of intrapatient variability in the concentration-effect relation. There is a finite, albeit small, probability of drug effect at very low drug concentrations, and at high drug concentrations there is still some probability that a drug effect will not be observed. This reflects what many studies have shown: specifically that the drug concentration needed to block the response of an individual to a repetitive and unchanging stimulus varies randomly around some mean value during the study. This also is consistent with the clinical observation that anesthetic requirements vary during a procedure, even when the stimulus is seemingly unchanging. The parameter gamma is a measure of this intrapatient variability and has substantial clinical significance. When gamma is large (> 6), the probability of effect is very close to zero when C is less than C^{50}(even if only slightly less) and very close to 1 when C is greater than C^{50}(even if only slightly greater). In other words, if gamma is large there is a well-defined threshold for drug effect. In contrast, when gamma is small, the drug effect-concentration curve is not as steep, and it is difficult to clearly define a threshold for drug effect. For example, if gamma = 1, Equation 1implies that the concentration must increase by a factor of 9 to increase from a 10% chance to a 90% chance of drug effect. In contrast, if gamma = 10, the concentration need only increase by a factor of 1.5 to change from a 10% to 90% chance of drug effect. Understanding the shape of the concentration-effect curve is important for clinical practice. [11,12] Although C^{50}is a single point of the concentration-effect curve, the parameter gamma characterizes the entire relation. For example, given C sub 50, the drug concentration associated with a 90% or 95% probability of drug effect (C^{90}or C^{95}) depends on gamma. Given that most anesthesiologists would prefer to maintain the anesthetic concentrations at C^{90}, or C^{95}rather than C^{50}, it seems that accurate estimation of gamma could be viewed as important as determining C^{50}.

Ideally, logistic regression would be applied to data taken from individual patients for evaluation of C^{50}and gamma and then these individual values could be used to generate a measure of central tendency (mean, median, and so on) for the population. But in reality, many important stimuli, such as tracheal intubation, loss of consciousness, skin incision, or sternotomy, can only be evaluated once during the study of any one patient and logistic regression is most often applied to pooled data of multiple patients with one data point per patient. Analysis of pooled data should account for interpatient variability around the typical concentration-response curve. Evaluating interpatient variability has proved useful in pharmacokinetic investigations, [13–15] and we believe it is important to develop methods to evaluate pharmacodynamic interpatient variability.

Although quantifying interpatient variability is of interest in its own right, it is even more important to account for interpatient variability to evaluate more accurately the shape of the individual concentration-effect curves. Evaluation of the steepness of the individual concentration-effect relation may be confounded by the lack of distinction between intrapatient and interpatient variability. This is illustrated conceptually in Figure 1, where the concentration-effect relations of nine hypothetical patients are presented. Although each curve is steep, the curve generated by one point from each patient is much flatter (and, in fact, the curve would have appeared even flatter if we had not selected low probability values from the leftward curves and high probability values from the rightward curves). Although this is a contrived example, it nevertheless illustrates how intrapatient variability may be overestimated if interpatient variability is significant. The resultant concentration-effect curve may appear artificially flat while the curve for any one patient may be steep. No readily accessible method is available to evaluate interpatient C^{50}variability or to distinguish interpatient C^{50}variability in pooled data analysis from flat individual concentration-effect curves. This article presents a technique to achieve these goals.

**Methods**

Appendix 1 presents the mathematical basis for this probit-based approach to pharmacodynamic analysis. The technique is based on the assumption that for individual patients the probability of a drug effect is given by Equation 2where Phi denotes the standardized cumulative normal distribution; that is, Phi[gamma lnC-gamma lnC^{50}] is the area under a Gaussian curve (with standard deviation of 1) from minus infinity to gamma lnC-gamma lnC^{50}(see appendix 1 for details). Thus Equation 2resembles Equation 1(logistic regression) in that it is characterized by two parameters, gamma and C^{50}, and the probability increases with increasing C.

In two separate studies, Vuyk et al. [7,9] collected data during gynecologic surgery on the response to repetitive intra-abdominal stimuli (electric cautery) of individual patients as a function of alfentanil plasma concentration when supplemented by nitrous oxide or various concentrations of propofol. Because this is, to our knowledge, the largest collection of individual concentration-response data, we used it to assess the validity of probit analysis for individual patients (Equation 2) by comparing the results to those of logistic regression (Equation 1). Observed responses and drug concentrations were taken from the figures of these reports. The parameters (C^{50}, gamma) in Equation 1and Equation 2were evaluated by maximum likelihood estimation. The logarithm of the likelihood of observed results Equation 3

(P is given by Equation 1or Equation 2and R = 1 if a drug effect is observed and R = 0 otherwise and I indexes patients) was maximized as a function of C^{50}and gamma. The estimation procedure was easily implemented using an Excel spreadsheet (Microsoft, Redmond, WA), taking advantage of its built-in function that evaluates the error function and the Solver function for the maximization step.

Appendix 1 describes generalization of Equation 2to account for interpatient variability. The probability of drug effect when concentration-response data are collected from multiple patients with one observation per patient is given by Equation 4where Phi again denotes the standardized cumulative normal distribution, but now the probability of drug effect is determined by three parameters, <C^{50}>, the mean value of C^{50}in the population; gamma; and omega^{2}, the variance of In C^{50}values in the population. Derivation of Equation 4is contingent on the assumption that the probability of drug effect for individual patients is given by Equation 2and that the distribution of C^{50}, values in the population has a log-normal distribution.

Direct validation of Equation 4was not possible because we are not aware of a sufficiently large data set comprised of single data points from multiple patients for which the interpatient variability of C sub 50, (omega^{2}in Equation 4) is also known. Instead, computer simulation was used for indirect validation. Each simulated study consisted of 30 different participants (a number fairly typical of many anesthesiology pharmacodynamic studies). A C^{50}, value and a gamma value were assigned to each hypothetical “patient” using the Excel random-number generator and assuming that both C^{50}and gamma had log-normal distributions. This is analogous to the enrollment of patients, during a true human participant study, from a population with randomly varying C^{50}and gamma values. For all simulations we assumed that the mean C^{50}(denoted <C^{50}>) was 100 units per milliliter, that the standard deviation of the distribution of log gamma was 0.5 (see Discussion to follow), and the mean value of gamma and the variance of C sub 50, were varied, as will be described. Each simulated patient was also assigned a single drug concentration from one of 30 values distributed uniformly from 10 to 300 units/ml. This corresponds in a real study to the investigator assigning a drug dose to each patient enrolled in the study. At this point, if the concentration-effect relation were completely deterministic (no intrapatient variability) then a positive drug effect would be observed if the drug concentration, C, assigned to the patient exceeded the value of C^{50}for that patient (see appendix 1). However, as noted earlier, there is an element of randomness even in the individual concentration-effect relation. To account for this intrapatient variability, a normally distributed random variable (denoted epsilon) with mean value of 0 and standard deviation of 1, was generated for each “patient” and added to gamma lnC-gamma lnC^{50}, (see equation A-1 in the appendix 1). If this sum exceeded zero, the simulated patient was assumed to have a positive drug effect (the response variable R was given a value of 1).

This simulation technique is best illustrated with an example. Each simulation was done on a spreadsheet and each row corresponded to a single “patient” and had entries for the assigned value of C and the randomly generated values of C^{50}, gamma, and epsilon and the sum gamma lnC-gamma lnC^{50}+ epsilon, which determined R. In one of the simulations, patient 6 was assigned a drug concentration of 80 units/ml. The values of C^{50}, gamma, and epsilon returned by random-number generation were 95, 9.4, and 0.9, respectively. Thus gamma lnC-gamma lnC sub 50 + epsilon was equal to -0.72, and this indicates that a drug effect was not observed (i.e., R = 0). In contrast, patient 13 was assigned a drug concentration of 90 units/ml. The values of C^{50}, gamma, and epsilon returned by random-number generation were 115, 5.4, and 1.4, respectively. For this patient, gamma lnC-gamma lnC^{50}+ epsilon was equal to 0.08 and R = 1.

Six different simulations were conducted by assigning gamma a mean value of either 10 (steep concentration-response curves), 5 (intermediate concentration-response curves), or 2 (flat concentration-response curves) and assigning omega^{2}values of either 1 or 0.04. Twenty-five repetitions of each simulation were conducted to assess the statistical properties of this technique. Maximum likelihood estimates of (C^{50}), gamma, and omega were calculated using Equation 3and Equation 4and implemented with the Excel spreadsheet. We also applied “naive” probit analysis (ignoring interpatient variability) to the simulated data sets by assigning omega the value of 0.

Appendix 2 provides a detailed description of how to use an Excel spreadsheet to implement this technique.

This technique was applied to real data presented by Ausems et al. [5] on the relation among alfentanil (when supplemented by nitrous oxide), plasma concentration, and the responses to intubation and skin incision, with data retrieved from the published graphs of Jacobs et al. [10] on the relation between midazolam plasma concentrations and loss of responsiveness, and of Bailey et al. [16] on the relation between sufentanil plasma concentrations and the responses to intubation, skin incision, and sternotomy in cardiac surgical patients, using the original data. For alfentanil, an indirect and approximate estimate of omega was made by equating it to the standard deviation of In C^{50}values for the electroencephalographic effect reported by Egan et al. [17] For midazolam, the standard deviation of patient ages in the study of Jacobs et al. [10] was used as a measure of omega, because these investigators found that age was an important covariate of ln C^{50}. For sufentanil, omega was directly calculated from reported values of C^{50}. [16]

**Results**

Population pharmacodynamic analysis using Equation 4is based on the assumption that a probit model (Equation 2) accurately describes the concentration-response relation for individual patients. We assessed this assumption using the data of Vuyk and colleagues, who have studied extensively the pharmacodynamics of alfentanil in combination with varying doses of propofol or nitrous oxide in healthy women undergoing lower abdominal gynecologic surgery. [7,9] They have presented concentration-response data from 38 patients. In seven cases the data were not amenable to either logistic or probit analysis (because a concentration was found above which there was always a positive drug effect and below which there was never a positive drug effect). For the remaining 31 patients, the log-likelihoods of the “best”(maximum likelihood) description of the data by either logistic regression or probit analysis are shown in Table 1, along with C^{50}values. There was no difference between the two approaches in their ability to describe these data.

The precision and bias of population analysis using Equation 4was evaluated by computer simulations. Figure 2and Figure 3show the distribution of population estimates of <C^{50}>, gamma, and omega for two of these simulations. For Figure 2, the simulated data were generated assuming that <C^{50}> = 100, gamma = 10, and omega = 1. Qualitatively, this represents steep individual concentration-response curves with substantial interpatient variability in the C^{50}value. For Figure 3, the target values are <C^{50}> = 100, gamma = 2, and omega = 0.2, which qualitatively represents flat individual concentration-response curves (large intrapatient variability) with relatively insignificant interpatient variability. The statistical fidelity of the method is indicated by how close the individual estimates are to these target values (shown as horizontal lines in the figures). For quantitative appraisal, the mean estimates of <C^{50}>, gamma, and omega for 25 repetitions of six different simulations (varying in the values of gamma and omega) are shown in Table 2, along with standard deviations of the estimates.

(Figure 4) compares the predicted concentration-response curves using parameters estimated by “naive” probit regression (Equation 2) and population probit regression (Equation 4) for simulated data from a population with <C^{509}> = 100, so that the “true” curve can also be illustrated (this is necessary because individual C^{50}values vary in the population). By normalizing the curve, we can use a single curve to represent the underlying population).

The results of applying Equation 4to real data are shown in Table 3. Also shown are independent but indirect and approximate estimates of omega. Also shown are estimates of gamma made accounting for interpatient variability (Equation 4) and by naive estimation (assuming that omega = 0 in Equation 4).

**Discussion**

We have presented a technique for population pharmacodynamic analysis of binary data, our goal being to improve the accuracy of estimating the parameter characterizing the shape of the individual concentration-effect relation and to provide a simple method to assess interpatient C^{50}variability. Although many pharmacodynamic studies concentrate on estimating C^{50}, this is only one point of the concentration-effect curve. In this study (as with logistic regression), the shape of the concentration-effect curve is largely determined by a single parameter, which we denoted gamma. Drugs with a relatively large value of 3 have a steep concentration-effect relation with well-defined thresholds for effect. Furthermore, if we believe that accurate estimation of C^{50}is important, then we should devote the same effort to estimating gamma, because it will determine C^{90}, or C^{95}and most clinicians would prefer to maintain patients at this level rather than C sub 50 (viewing a 50% chance of awareness, for example, as undesirable). In short, accurate assessment of this aspect of pharmacodynamics is essential for the clinical use of the drug.

The technique we describe characterizes the concentration-response curve with three parameters, <C^{50}>, gamma, and omega. <C^{50}> is the mean C^{50}value for the population and 3 is a measure of how steep the individual concentration-response curve is. Because this curve is expressed as the probability of drug effect as a function of drug concentration, gamma is a measure of intrapatient variability. In contrast, omega^{2}is a measure of interpatient variability, the variance of ln C^{50}, values in the population. By accounting for interpatient variability, we hope to avoid confounding intrapatient and interpatient variability, leading to artificially low estimates of gamma.

In a previous report, we defined a measure of recovery, denoted mean effect time, which depended closely on the steepness of the concentration-effect relation as measured by gamma. [11] In an editorial comment, Schnider and Shafer [12] noted that if estimates of gamma were artificially low, then the calculated mean effect time would grossly exaggerate recovery time. A source of error in the estimation of gamma is the failure to consider interpatient variability in C^{50}. Methods of data analysis that use pooled data from a population and do not account for the fact that different patients will have different C^{50}values may lead to inaccurate estimates of gamma. This is analogous to the situation in pharmacokinetic analysis [13–15] and is illustrated by Figure 4, which shows that if there is substantial interpatient variability of C^{50}values (in this case, a standard deviation for ln C^{50}of 1), the concentration-response curve estimated by naive application of logistic regression to pooled data is very flat, even when the individual curves are steep (in this case, a mean gamma value of 10). In our applications of this technique to real data (Table 3), we see that “naive” estimates of gamma (ignoring interpatient variability and assuming that omega = 0) are considerably lower than those that account for interpatient variability, as predicted by the simulation illustrated in Figure 4.

We can illustrate these points quantitatively by calculating mean effect time for midazolam, as an example, using “naive” versus population estimates of gamma (Table 3). The technique for calculating the mean effect time is described by Bailey, [11] with the substitution of Equation 2, rather than Equation 1, for the probability of drug effect, and we consider the specific case of a duration of administration of 100 min. After discontinuing administration, the mean effect time based on the naive estimate of gamma (1.31) is 42 min. Using the population-based estimate of gamma (5.05), the estimate is 14 min. The naive estimate exaggerates the estimate of recovery time. However, we must also emphasize that accurate estimation of gamma is important for understanding the pharmacology of the drug only if the clinician titrates the drug to effect. If drugs are administered without end points for titration, the naive estimate of gamma may actually be more informative, although this would be grossly misleading about the concentration-effect curve of the typical patient. For example, if we seek the drug concentration that ensures that an effect will be observed in 90% of all patients, then the naive estimate of gamma should be used for the calculation. However, this drug concentration may be a gross overdose for any one patient if we can titrate to effect. The population-based estimate of gamma is a better reflection of the concentration-effect curve in this situation, although both estimates are useful and we recommend using both types of analysis in pharmacodynamic research.

It is also of some interest that the mean effect time calculated using a naive estimate of gamma based on logistic regression (Equation 1) was considerably larger than the estimate using the naive probit (Equation 2)(89 versus 42 min). This presumably reflects the fact that the logistic distribution has longer tails than the normal distribution and may confirm the speculation by Schnider and Shafer [12] that calculations of mean effect time may be inflated due to the tails of the logistic distribution.

The equation we derive to describe the probability of drug effect (Equation 4) is not as straightforward and easily understood as the one used for “naive” logistic regression of pooled data (Equation 1). Thus a discussion is warranted. Equation 4relates the probability of drug effect for the population to the cumulative normal distribution, by convention denoted phi(x). If we plot a standardized normal curve (a Gaussian or bell-shaped curve with a standard deviation of 1) on the x-y axis, centered around x = 0, phi(x) is the area under this curve to the left of x. Note that x can be either positive or negative; if x = 0, then phi(x)= 0.5, because x = 0 is the midpoint of the curve. In Equation 4, x =(gamma lnC - gamma ln<C^{50}>)/radical(1 + omega^{2}gamma^{2}). So as C increases, phi(x) also increases, which is logical because the probability of drug effect should increase as the drug concentration increases. The probability of drug effect is 50% when C = <C^{50}> because x =(gamma lnC - gamma ln<C^{50}>)/radical(1 + omega^{2}gamma sup 2)= 0 at this point and the normal curve is centered around x = 0. Also, note that as gamma becomes larger, x =(gamma lnC - gamma ln<C^{50}>)/radical(1 + omega^{2}gamma^{2}) increases more per unit increase in concentration; that is, the concentration-effect curve is steeper. Finally, as omega becomes larger (greater interpatient variability), the increase in x and phi(x) per unit increase in drug concentration decreases; that is, the concentration-curve for the population is flatter. This reflects the fact that with greater interpatient variability there will be more patients with extremes of C^{50}, both low and high, so that for the population the probability of drug effect at low drug concentrations is increased while it is decreased at high concentrations, flattening the curve.

Within the framework of the model and its assumptions, the mathematical derivation of Equation 4is exact. Thus the utility of this technique is determined solely by the validity of the assumptions of the model. The critical assumptions are (1) the concentration-response relation for individual patients is described by Equation 2, a probit relation, rather than by the more commonly used logistic equation (Equation 1);(2) the distribution of C^{50}values in the population conforms to a log-normal relation; and (3) we can ignore the interpatient variability of gamma without substantial loss of accuracy.

The first assumption, that individual concentration-response curves are appropriately described by a probit model, was directly evaluated using data for individual patients and was compared with logistic regression. Although logistic regression has been commonly used for binary response pharmacodynamic analysis, we know of no previous effort to determine whether intrapatient variability is well described by a logistic distribution. In our analysis of the data of Vuyk et al., [7,9] we found that the probit model (which uses a normal distribution to describe intrapatient variability) fit the data as well as the logistic model. This reflects the fact that the logistic and normal distributions are not that dissimilar, both being symmetric and sigmoid. We cannot conclude that the probit model is superior to the logistic model for individual patients, but it seems apparent that it is just as appropriate.

The second assumption of the model is that the distribution of C^{50}values conforms to a log-normal distribution. This type of assumption is commonly used in both pharmacokinetic and pharmacodynamic modeling and conforms to observed distributions of pharmacokinetic and pharmacodynamic parameters. [13] However, if the actual distribution of C^{50}values are multimodal, this assumption is clearly incorrect.

The final critical assumption of this model is that it is “gamma naive”; that is, we have not accounted for variability in gamma values. Although we have referred to gamma as the parameter that determines the steepness of the concentration-response relation, the slope of this curve is also influenced by C^{50}. Thus ignoring gamma variability is not equivalent to assuming that the shape of the concentration-response curve is the same for all patients. In addition, interpatient gamma variability could be incorporated into the model if the distribution of gamma were normal. As we noted previously, a log-normal distribution seems more plausible. However, if the variance is not too large, a normal distribution could approximate a log-normal, as a “first-order” approximation similar to that used in NONMEM. [13,14] and the model could be expanded to consider interpatient variability in gamma. However, we have not extended the model in this manner because our primary focus is the variability of C^{50}. Accounting for gamma variability would add a parameter to the model (two additional parameters would be required if covariance between C^{50}and gamma is considered). Few anesthesiology data sets of this type are large enough to support accurate estimation of four (or five) parameters. In our simulations, we have assumed a log-normal distribution of gamma values in the population, with a standard deviation equal to 50% of the mean value, a substantial variability. We chose this value because the mean standard deviation of ln gamma values among the separate patient subsets of the data from Vuyk et al. [7,9](after “trimming” the highest and lowest gamma value in each group) was 0.477. Thus a ln gamma standard deviation of 0.5 seems realistic. We found that our estimates of the mean value of gamma were accurate despite our “gamma naive” estimating equation (Equation 4). This seems to parallel the observation that the “omega^{2}naive” approach, logistic regression, is a good estimator of <C^{50}> (see Figure 4) but may fail in the estimation of gamma, a measure of variability (intrapatient). Similarly, our “gamma naive” technique estimates <gamma> reasonably accurately but provides no information about its variability. Although our model does not fully account for interpatient variability, we believe it is an improvement over the fully “naive” logistic regression approach.

A possible concern with this project is that we cannot validate our technique directly. There are simply no data sets available consisting of single concentration-response observations from multiple patients for whom the variance of C^{50}in the population is known. Consequently, we evaluated the statistical properties of the technique using computer simulation. We “created” 30 patients for each simulation by randomly generating values of C^{50}and gamma, analogous to the experimental investigator enrolling patients from a population with log-normal distributions of these variables, and we assigned a single drug concentration to each patient, analogous to the experimental investigator selecting a drug dose for each patient in the study. We then accounted for intrapatient variability by letting the response variable, R, be determined by the sum of gamma lnC - gamma lnC^{50}and a term representing random intrapatient variability, generated by a random-number generator. Thus the simulations emulated a typical pharmacodynamic study. We wish to emphasize that computer simulation is a commonly used validation procedure for statistical methods, and our approach parallels the validation of widely used population pharmacokinetic techniques. [14] We believe our simulations are realistic. They are based on a probit model (Equation 2) for the individual concentration-response relations, which we have shown is a reasonable assumption, and on the assumption of log-normal distributions for C^{50}and gamma. If these assumptions are plausible, our simulations indicate that this technique can accurately estimate <C^{50}>, gamma, and omega.

We illustrated this technique by applying it to three sets of actual patient data, that reported by Ausems et al. [5] for the relation between alfentanil (supplemented with nitrous oxide) concentration and the response to intubation or skin incision, that reported by Jacobs et al. [10] for the loss of responsiveness and midazolam plasma concentration, and data for sufentanil reported by Bailey et al. [16] In each of these cases we have approximate, indirect measures of omega, and the estimates with our technique and these indirect measures are in reasonable agreement. But this observation must be tempered by the approximate nature of the indirect measures. In the case of alfentanil, the indirect measure is the standard deviation of electroencephalographic measures of ln C^{50}reported by Egan et al. [17] Whether this has an important relation to the somatic or sympathetic response to intubation or incision is uncertain. In the case of midazolam, the investigators found that patient age was a significant covariate in the estimation of ln C^{50}; that is, ln C^{50}was proportional to age. [10] The independent estimate of omega in the table is simply the standard deviation of the patient ages in the study (multiplied by the coefficient of proportionality reported by the authors) and so should probably be viewed as a lower bound to the true standard deviation of C^{50}in the population. For the sufentanil data, the estimate of omega is actually direct but must be viewed as crude because the number of patients used for the estimate was small (n = 5) and the responses to different stimuli (intubation, incision, and sternotomy) were pooled for analysis.

There is a considerable statistical literature on the analysis of binary or ordered categorical data (see McCullagh and Nelder [18] and Zeger and Liang [19] for overviews) and, in particular, Sheiner [20] has described the analysis of categorical data within the framework of the NONMEM program. We have not examined the analysis of binary data with NONMEM, primarily because the probit technique is so easily implemented with an Excel spreadsheet (as detailed in appendix 2). However, in the analysis of binary data it is necessary to evaluate complicated integrals. In our probit model, these integrals (the cumulative normal distribution) are provided by a built-in function with a high degree of accuracy. With NONMEM, the integrals are evaluated by a technique known as the Laplace approximation. [20,21] To our knowledge, application of the Laplace approximation to the analysis of sparse binary data (one data point per patient) has not been validated in simulations such as those we used here.

We have presented a method for population analysis of binary, “yes/no,” pharmacodynamic data. The technique estimates the mean value of C^{50}, a parameter gamma that influences the steepness of the concentration-response curve, and omega^{2}, the variance of lnC^{50}in the population. The model assumes a log-normal distribution of C^{50}values. The underlying model for the individual concentration-response relations was validated by analysis of individual patient data. The statistical properties of the technique were evaluated by computer simulation.

**Appendix 1**

Although logistic regression is commonly used in anesthesiology research, the underlying statistical model is seldom described. A useful general reference to these types of statistical models is the textbook by McCullagh and Nelder. [18] This model assumes that anesthetic effect is measured by an underlying, continuous variable, denoted y, that is related to the drug concentration by the equation Equation 5where epsilon is a random variable that is assumed to have a logistic distribution. Because we describe the concentration-effect curve by two parameters (C^{50}and gamma), we can arbitrarily assume that the logistic distribution of epsilon has a mean of zero and a standard deviation of pi^{2}/3. The model is applicable to binary, “yes/no” data because we assume that a positive drug effect is observed only if y greater or equal to 0 or equivalently if epsilon greater or equal to -(gamma lnC - gamma lnC^{50}). Thus the probability of drug effect is Equation 6.

(The integrand is the logistic density function.)

As an alternative, we propose using probit analysis. The basic model is again Equation 7but epsilon now has a normal distribution with a mean of zero and a variance of 1. Again we assume that a drug effect is observed only if y greater or equal to 0. Also note that again the concentration-effect curve is characterized by two parameters, gamma and C^{50}(which is why the normal distribution may be arbitrarily assumed to have a variance of 1). With this model the probability of observing a drug effect is Equation 8where phi denotes the cumulative standardized normal distribution.

(Equation 8) is applicable to data from a single patient. For the analysis of pooled population data with one observation per individual, this expression must be generalized. We rewrite Equation 7as Equation 9where the subscript I indexes patients. Note that we assume that gamma is the same for all patients. We also assume that C^{50}has a log-normal distribution; that is, C^{50}= <C^{50}> exp (delta^{i}) where <C^{50}> is the mean C^{50}value and delta^{i}is normally distributed with a mean of zero and in unknown variance omega sup 2 and is independent of epsilon^{i}. With these assumptions we can write Equation 10.

Because the sum of two independent normally distributed variables is itself a normally distributed variable whose variance is the sum of the variances of its individual components, in this model y^{i}is a normally distributed variable with variance equal to 1 +(gamma omega)^{2}. Because we observe a drug effect if y^{i}greater or equal to 0, the probability of a drug effect in the ith individual for population analysis is Equation 11.

(For convenience we have dropped the index I). Note that the probability of drug effect for this population model is now characterized by three parameters: <C^{50}>, the mean value of C^{50}for individual patients, gamma a measure of the steepness of the concentration-response curve, and omega^{2}, the variance of C^{50}values in the population.

**Appendix 2**

In this appendix we describe the implementation of population probit analysis using an Excel spreadsheet. This implementation is illustrated in Figure 5. These instructions assume a computer with a mouse. For personal computers, cells are highlighted by using the mouse to align the “[center dot]” over a cell and striking the mouse button. Options from the menu above the spreadsheet itself are selected by using the mouse to align the “[center dot]” over the option and clicking the mouse button.

1. Enter data in columns A and B. The concentrations at which observations of drug effect were made are entered in column A (cells A1:A33 in this example). The response variables are entered in column B (enter a “1” if a drug effect was observed and “0” otherwise).

2. Enter initial estimates, or guesses, of <C^{50}> in cell C1, of gamma in cell C2, and of gamma in cell C3.

3. Highlight cell D1. Type “=($c$2*ln(a1)-$c$2*ln($c$1))/((1 +($c$2*$c$3)[center dot] 2)[center dot] .5)”.(The quotation marks are not included in the expression to be typed). Strike the Enter key. Now position the mouse over the lower right corner of cell D1 until the “+” symbol appears (distinguish this from the “[center dot]” symbol). While depressing the mouse button, drag the mouse down column D, in this example to row 33. This will fill in cells D2:D33 with the above expression except that the drug concentrations in cells A2:A33 will be used (in place of the value in cell A1) for cells D2:D33.

4. Highlight cell E1. Select the “function wizard,” which is the box labeled “f^{x}” in the menu above the spreadsheet. A menu of categories of functions will appear on the screen. Select “Statistical.” A menu of statistical functions will appear. Scroll down to “NORMS-DIST” and select this option. Then select the “Next” option. A line for inserting the argument of the function will appear. Type “d1” and strike the Enter key. This will cause the value of the cumulative standardized normal distribution of the value in cell D1 to appear in cell E1. Place the mouse over the lower right-hand corner of cell E1 until the “+” symbol appears. Keeping the mouse button depressed, drag the mouse down column E to fill in cells E2:E33 with the values of the normal distribution corresponding to cells C2:C33 and D2:D33.

5. Highlight cell F1. Type “= b1*log(e1)+(1 - b1)*log(1 - e1)” and strike the Enter key. Place the mouse over the lower right-hand corner of cell F1 until the “+” symbol appears, and while keeping the mouse button depressed drag the mouse down column F to fill in cells F2:F33.

6. Highlight the cell below the last cell entered in column F. In this example, this is cell F34. Select the icon labeled “Sigma” in the menu above the spreadsheet and strike the Enter key. This will enter the sum of the values in column F (the log-likelihood values).

7. Select the “Tools” option in the menu at the top of the spreadsheet and then select the “Solver” option. Select “Set Target Cell” and designate the cell containing the sum of the log-likelihood values, in this case typing “f34”. Select “Max” and then select “By Changing Cells” and type “c1,c2,c3”. Strike the Enter key. This will result in maximum likelihood estimation by varying the values of <C^{50}>, gamma, and omega in cells C1:C3. This process should be repeated using multiple starting values in cells C1:C3. Convergence of estimation can be improved by constraining <C^{50}>, gamma, and omega to have positive values. This option is available in the Solver menu.