Free
Education  |   December 2003
Reliability of Pharmacodynamic Analysis by Logistic Regression: Mixed-effects Modeling
Author Affiliations & Notes
  • Wei Lu, Ph.D.
    *
  • James G. Ramsay, M.D.
  • James M. Bailey, M.D., Ph.D.
  • *Research Fellow, †Professor, ‡Associate Professor.
  • Received from the Department of Anesthesiology, Emory University School of Medicine, Atlanta, Georgia.
Article Information
Education
Education   |   December 2003
Reliability of Pharmacodynamic Analysis by Logistic Regression: Mixed-effects Modeling
Anesthesiology 12 2003, Vol.99, 1255-1262. doi:
Anesthesiology 12 2003, Vol.99, 1255-1262. doi:
A GOAL of pharmacodynamic modeling is to quantitatively describe the response to pharmacologic agents in terms of a small number of parameters. For example, one of the more common pharmacodynamic models assumes that the effect of a drug (E) is described by the sigmoid Emax, or Hill, equation, E = Emax[Cγ/(Cγ+ C50γ)], where Emaxis the maximal effect; C is the drug concentration; C50is the drug concentration associated with an effect equal to 50% of Emax; and γ is a measure of how steep the drug concentration–effect relationship is. 1 This model is useful when the effect is a continuous variable. In anesthesia research, however, the response to a drug is often not a continuous, but rather a binary variable. For example, the patient may or may not be responsive after administration of a hypnotic agent. In this case, Emaxis equal to 1 and the effect becomes the probability of a particular response, given by P = Cγ/(Cγ+ C50γ). The concentration–effect relationship is thus determined by two parameters, C50and γ. 2–12 
This article is accompanied by an Editorial View. Please see: Sani O, Shafer SL: MAC attack? Anesthesiology 2003; 99:1249–50.
Pharmacodynamic data are often sparse. That is, there may be many data points, but there are few data points from the same patient. This is a reflection of the difficulty of repetitively observing the response to drug administration in the same patient. In a previous study, we investigated the accuracy of parameter estimates when data consisting of one data point per patient was analyzed in a naive manner, i.e.  , an analysis in which the interpatient variability of C50was ignored. 13 We found that naive analysis of sparse data resulted in relatively accurate estimates of C50. However, there was a substantial bias in the estimate of γ. In the current study, we investigated whether a population analysis in which interpatient C50variability is taken into consideration, results in less biased estimates of γ.
Methods
Simulations
Simulations were performed as previously described. 13 Excel (Microsoft, Redmonds, WA) was used to generate random values of C50and γ, assuming that both parameters had a log-normal distribution, i.e.  , that the parameters are distributed as P = PTVexp(η), in which P denotes the parameter (C50or γ); PTVdenotes the typical value of the parameter; and η has a normal distribution with mean value of zero and SD of 0.3. We assumed that PTVfor C50was 100 and that PTVfor γ was 1, 1.5, 3, 4.5, or 6. Simulations were performed for varying numbers of patients (n = 10, 20, 30, 40, 50, 75, and 100) and varying numbers of simulated data points for each patient (m = 1, 3, 5, and 10). To begin each simulation, data points were generated by randomly selecting n (number of patients) × m (number of data points per patient) drug concentrations distributed uniformly on a logarithmic scale from 25 to 400 units by using the random number generator of an Excel spreadsheet. This corresponds, in a human study, to the investigator assigning a drug dose to each patient enrolled in the study, with variability in the plasma concentrations resulting from the dose. At this point, if the concentration–effect relationship were totally deterministic (i.e.  , γ were infinite), a positive drug effect would be observed if the drug concentration, C, assigned to the data point exceeded C50. However, because we assess the concentration–effect relationship in terms of the probability of the effect, there is obviously an element of intrapatient randomness in the concentration–effect relationship, and this is reflected in the fact that γ is finite. To take into account this randomness, a uniformly distributed random variable from 0 to 1 was generated for each data point, again by using the random number generator of an Excel spreadsheet. If this number was less than Cγ/(C50γ+ Cγ), the simulated patient was assumed to have a positive drug effect; the response variable, denoted R, was given a value of 1. Otherwise, it was assumed that R is equal to 0. In this manner, responses were obtained for a range of concentrations, and spreadsheets consisting of columns of C, the drug concentration, and R, the response variable, were generated. For each simulation, the parameters (C50and γ) were then estimated by maximal likelihood estimation by using mixed-effects modeling. 14 This was performed by using the software package NONMEM (University of California, San Francisco, San Francisco, CA) version 5. 15 The Laplacian estimation technique was used. To minimize the number of parameters estimated, the initial analysis was performed by assuming that C50had a log-normal distribution, whereas γ did not vary from patient to patient. A subsequent analysis was performed by assuming that C50and γ had log-normal distributions. Each simulation was repeated 100 times. Simulations were performed for starting values of C50, γ, and <η2> equal to the true values, to 50% of the true values, or to 200% of the true va(lues. Bias was defined as the estimate minus the true value with the difference normalized to the true value of the parameter estimate.
C50, γ, and √(<η2>) bias were calculated, as were the coefficients of variation (SDs of parameter estimates normalized to the mean estimates and denoted CV).
Human Studies
After obtaining approval from the Emory University School of Medicine Human Investigations Committee and written informed consent, 30 patients were enrolled in a study of postoperative sedation of adult patients undergoing cardiac surgery. The study was conducted in the intensive care unit after elective coronary artery bypass graft surgery. Patients with decreased left ventricular function (ejection fraction of <40%), chronic pulmonary disease, neuromuscular disease, or a history or laboratory evidence of renal, hepatic, or hematologic disease were excluded. Patients were also excluded from the study if they required an intraaortic balloon pump, showed evidence of active bleeding, or were expected to be intubated for more than 72 h.
All patients were premedicated with diazepam, 0.1 mg/kg; morphine, 0.1 mg/mg; and scopolamine, 0.05 mg/kg. Anesthesia was induced with sodium thiopental and was maintained with fentanyl, 50 mcg/kg, or sufentanil, 10 mcg/kg, supplemented with midazolam (up to 10 mg) or sodium thiopental (up to 4 mg/kg) during cardiopulmonary bypass. Muscle relaxation was achieved with pancuronium or vecuronium.
After arrival in the intensive care unit, the level of sedation was evaluated by a research nurse (not a caregiver) using the following scale 16 : 1, no response to pain; 2, responds to pain only; 3, eyes closed, calm, responds only to loud verbal or physical stimulus; 4, eyes closed, calm, responds to verbal stimulus; 5, awake; 6, agitated.
Once patients had a sedation score of 3 or higher, midazolam, 0.015 mg/kg, was administered every 2 min until a sedation score of 2 was achieved or until a maximum of 0.2 mg/kg was given. After this first dose or series of doses, when a sedation score of 4 was observed, midazolam, 0.015 mg/kg, was given to achieve a sedation score of 3. If patients were coughing, despite achieving the target sedation score, morphine was given in 2-mg increments to suppress coughing. If patients indicated pain in response to a direct question, 2–5 mg morphine was given as needed. The last dose of midazolam was given after midnight on the day of surgery, but before 5:00 am of the first postoperative day at the discretion of the research nurse.
Arterial blood samples were drawn 5 min after each dose or series of doses. Plasma midazolam concentrations were assayed by high-performance gas chromatography and mass spectroscopy. 17 The interassay precision was 4.1%, and the intraassay precision was 4.4%. The limit of detection was 10 ng/ml.
To compare our ordinal clinical data with the results of simulations of a binary effect, we converted the ordinal clinical data to binary data by assuming that a positive drug effect was a sedation score of 3 or less, whereas a score of 4 or greater indicated a negative drug effect. Data were analyzed with NONMEM by using the Laplacian estimation technique. It was assumed that C50had a log-normal distribution in the population. Interpatient variability in γ was ignored. Furthermore, because these data were collected at a time when high-dose opioid and heavy premedication was used for cardiac anesthesia, we considered a model that accounted for the residual and decaying effects of opioids, premedication, and effects of cardiopulmonary bypass on the level of consciousness by including a “virtual” drug, as first introduced by Somma et al  . 18 This virtual drug was modeled by assuming that the effective midazolam concentration is equal to the measured concentration plus a monoexponential term, A exp(−kt), in which t is the time in the intensive care unit and A and k are estimated by nonlinear regression.
A total of 656 data points were available. To simulate a population analysis of sparse data, a random number generator was used to select one data point for each patient. The set of data points for all 30 patients was then analyzed with NONMEM as described above and incorporating the virtual drug effect. This random selection of a sparse data set was repeated 100 times. The results were compared with the results of a NONMEM analysis of all 656 data points.
Results
Simulations
The bias of C50estimates was small (within 5% of the true value), even for sparse data (one data point per simulated patient), as long as the total number of patients in the simulation was 20 or more (fig. 1). This was true for all values of γ, although we present only the results for γ= 1, 3, and 6 and omit the results for γ= 1.5 and 4.5. This was also true whether interpatient variability of γ was ignored in the estimation procedure and regardless of the parameter value used to start the process. The results shown in fig. 1and subsequent figures present results in which γ interpatient variability was ignored and the starting parameter estimates were equal to the true value. Despite the insignificant bias, we found that the coefficient of variation of C50estimates was greater for sparse data (one data point per patient) than for dense data (10 data points per patient). Figure 2shows graphs for sparse and dense data of the coefficient of variation of C50estimates as a function of the number of simulated patients and γ.
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel 
	shows the results for one data point per patient, and the lower panel 
	shows the results for 10 data points per simulated patient.
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
×
In contrast to the insignificant bias of C50estimates, there was a large bias in the estimates of γ for sparse data when the true value of γ was 3 or greater. Figure 3shows γ estimate bias as a function of numbers of simulated patients and number of data points per patient for γ= 1, 3, and 6. Bias is less pronounced for a low value of γ. For larger values of γ (greater than 1), the bias for γ estimates was much higher than the bias for C50estimates, unless the data were dense (10 data points per patient) and there were at least 20 patients in the simulated study. The results were virtually identical regardless of whether interpatient γ variability was ignored in the estimation process. The results in figure 3were obtained by using starting values in the estimation procedure equal to the true value. When lower or greater starting values are used, the bias is even greater.
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
The potential for bias in estimation of γ when there was only one data point per simulated patient was not always indicated by the SE of the estimate. A large SE would indicate the possibility of significant bias. The SE of the estimate of γ was large enough to encompass the true value in only 30–60% (depending on the number of simulated patients) of the simulations for a true value of γ= 6 and in 50–80% of the simulations for a true value of γ= 3.
In general, variability of γ estimates was also greater than that of C50estimates. Figure 4illustrates this for sparse data (one data point per patient) and dense data (10 data points per patient) with graphs of the coefficient of variation as a function of the number of simulated patients and γ.
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel 
	shows the results for one data point per patient, and the lower panel 
	shows the results for 10 data points per simulated patient.
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
×
Although estimation of C50was performed with minimal bias, this was not true for the estimation of the population variance of C50, denoted ETA. Figure 5shows the bias of ETA estimates as a function of the number of simulated patients and number of data points per patient for γ= 1, 3, and 6. In contrast to the bias of γ estimates, the bias of ETA estimates was greatest for a true value of γ equal to 1. In general, ETA bias decreased with the number of simulated patients and the number of data points per simulated patient.
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
Clinical Study
There were at least 10 data points for each patient in the clinical study. The basic model that ignored a virtual drug effect had an objective function (−2 times the logarithm of the likelihood of the results) of 710, and the estimate of C50was 41.6 ± 3.8 and the estimate of γ was 2.69 ± 0.43. Inclusion of a virtual drug in the model decreased the objective function for the complete data set by 99 units, which is a highly significant improvement in the quality of the fit. The estimate of C50was 47.5 ± 4.7 ng/ml; the estimate of γ was 4.17 ± 0.59; the coefficient of the virtual drug effect was 63.5 ± 25.8 ng/ml; and the rate constant of the virtual drug effect was 0.00903 ± 0.00416 min−1. When sparse data sets were generated from the complete data set, the mean estimate of C50was 67.2 ± 53.5 ng/ml. The mean estimate of γ from the 100 sparse data sets was 14.4 ± 7.88.
Discussion
The results of the simulations of this study confirm our previous finding that C50may be estimated with insignificant bias with sparse data. 13 In the previous study, we examined the use of naively pooled data analysis, whereas the current study looked at population analysis, implemented with NONMEM. It is, of course, not surprising that the mathematically sophisticated analytic technique, NONMEM, produced unbiased estimates, given our previous observation that estimates derived from naively pooled data analysis were unbiased. We observed that the coefficient of variation of C50estimates with NONMEM was generally less than we reported previously for naively pooled data, but the difference was small. It appears that the use of population methods does not significantly enhance our ability to estimate C50with sparse data. The coefficient of variation of C50estimates was greater for sparse data (one data point per patient) than for dense data (10 data points per patient), as expected, because the variance of an estimator of a typical value will be of the order of 1/(nm).
In contrast to C50, estimates of γ by NONMEM were usually highly biased with sparse data. Bias was less pronounced for the lowest true value of γ= 1 and was particularly notable for γ= 6. In general, NONMEM overestimated the value of γ. This is an interesting contrast to our previous observation that naively pooled data analysis underestimates γ. The coefficient of variation for estimates of γ reflected this bias and was greater than we previously reported for naively pooled data analysis. Bias of γ estimates decreased as the number of data points per patient increased, but bias became insignificant only if there were 10 data points per patient.
The estimation of the parameters used to generate figures 1–5were performed by assuming that each simulated patient had the same value of γ; that is, we ignored interpatient variability. We used this approach for our initial estimations to minimize the number of parameters that needed to be estimated. However, when we assumed that γ had a log-normal distribution in the estimation step, the results were essentially identical (data not shown). Also, the results presented in figures 1–5were done by using starting parameter estimates equal to the true values. We did this to present the best case scenario for mixed-effects modeling. When we used lower or higher starting values, the bias in the estimates of γ were often dramatically worse.
We further investigated the reliability of C50and γ estimation by using population analysis using human data on the sedation of patients with midazolam after cardiac surgery. The data set consisted of 656 matched observations of sedation score and midazolam plasma concentration from 30 patients. The original six-point ordinal sedation score was converted to binary data by defining a positive drug effect as a patient who was unresponsive or responsive only to pain or loud verbal or physical stimulus. These data were collected in an era when patients undergoing cardiac surgery received high doses of opioids and heavy premedication, and our modeling confirmed the conclusion of Somma et al.  that postulating a virtual drug with monoexponential kinetics to account for the decaying effects of opioids and premedication improved the fit of the model to the data. 18 When the full data set, which contained more than 10 data points per patient, was analyzed, the estimate of C50was 47.5 ng/ml and the estimate of γ was 4.17. We then randomly generated 100 sparse data sets (one data point per patient) and analyzed each of these. The mean estimate of C50from sparse data was 67.2 ng/ml, close to the result from the full data set. However, the mean estimate of γ was 14.4, which was noticeably higher than the result from the full data set. This confirms our observation that population analysis, at least as implemented by NONMEM, overestimates the true value of γ with sparse data.
In our previous study of naively pooled data, we found that γ was underestimated, 13 whereas in this study, we report overestimation of γ. It is instructive to consider why NONMEM overestimates γ, in contrast to naively pooled data. With naively pooled data, interpatient variability in C50is ignored. The response of the population does not resemble any individual. However, with mixed-effects modeling, each individual is allowed his or her own C50. Now consider the situation in which we make a single observation, a positive drug response at a drug concentration C. The probability of that response is P = Cγ/(Cγ+ C50γ), and the likelihood of the observed response is maximized by any estimate of C50less than C, as long as the estimate of γ=∞. Similarly, if the single observation were a negative response, the probability of the observation is P = Cγ/(Cγ+ C50γ), and again, the likelihood of the observed result is maximized by any estimate of C50greater than C and an estimate of γ=∞. When data from multiple patients is pooled, NONMEM estimates <C50>, the typical value of C50in the population. If the drug concentration, C, exceeds <C50> and a positive drug response is observed or if the drug concentration, C, is less than <C50> and a positive drug response is not observed, NONMEM could maximize likelihood with an infinite estimate of γ (fig. 6) without assuming any interpatient C50variability. Of course, infinite estimates are not returned, and the reason is that, because of intrapatient variability, in any real study, we would expect to observe cases in which C is less than <C50>, but there is a positive drug response and vice versa  . This intrapatient variability is reflected in a finite value of γ. Thus, there are two elements of the variance of the model, interpatient C50variability and intrapatient variability. NONMEM estimates parameters by finding those parameters that minimize an objective function (−2 times the logarithm of the likelihood of the observed results), which applies a penalty to the model variance. If γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50(fig. 7).
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve  illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve 
	illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve  illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
×
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves  illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve  ) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves 
	illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve 
	) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves  illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve  ) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
×
Because of our observation that γ may not be estimated accurately from sparse data, it is interesting to ask how many drugs commonly used in anesthesia have been studied with sufficiently dense data to estimate γ accurately. Our simulations suggest that 10 data points per patient are necessary. As noted earlier, it is difficult to generate data sets this dense in anesthesia research. For example, in the study of a hypnotic used for intravenous induction of anesthesia, it is almost impossible to repetitively observe the response (i.e.  , loss of consciousness) repetitively in the same patient. However, there is one situation in which this type of research is more feasible and that is the study of sedation in an intensive care unit. Barr et al.  19 administered propofol for intensive care unit sedation and determined C50and γ values for modified Ramsay sedation scores similar to those used in our human study. Their data set comprised 643 observations from 20 patients, which is dense by the criteria we have identified. Thus, we think their estimate of γ by mixed-effects modeling, 1.7, is accurate and unbiased. The same group also studied intensive care unit sedation with midazolam and lorazepam, again determining C50and γ for modified Ramsay sedation scores. 20 Although the number of observations is not cited, inspection of figure 1of reference 20 clearly indicates that the data set was dense (greater than 10 observations per subject). The reported values of γ estimated by the naively pooled data approach was 4.5 for midazolam and 3.5 for lorazepam. This estimate for midazolam is within the SE for the estimate from our dense data set using mixed-effects modeling. We have previously reported that estimation with naively poled data tends to underestimate γ for sparse data, but we have performed no simulations of this question for dense data. 13 The agreement of the estimate of γ from our mixed-effects analysis and the naively pooled data analysis of Barr et al  . 19 suggests that naively pooled data may be accurate for dense balanced data sets.
If the generation of a dense data set is not feasible, as in a typical minimal alveolar concentration study, there seem to be only two alternatives for analysis. One could use the naively pooled data approach, realizing that it will underestimate γ but will avoid the infinite bias of a single observation. Alternatively, one could fix the variance of C50based on previous data, such as the variance in the minimal analgesic concentrations.
One could ask how important the estimation of γ is. In pharmacodynamic modeling, we attempt to describe the concentration–effect relationship in terms of a small number of parameters. The midpoint of this relationship is C50, which can be viewed as a location parameter. The determination of C50values is a major goal of pharmacodynamic modeling. However, the value of C50tells us nothing about the shape or, alternatively, the scale of the concentration–effect relationship. Estimation of C50only is incomplete knowledge. As an example of the clinical importance of the shape of the concentration–effect relationship (γ), consider two hypnotic agents, both of which have C50values of 100 units/ml, but for drug A, γ equals 1, and for drug B, γ equals 6. We make the assumption that, during maintenance of hypnosis, the practitioner would not be satisfied with the patient having a 50% chance of awareness and will assume arbitrarily that it would be more appropriate to maintain drug levels that guarantee at least a 95% chance of hypnosis. We also assume that recovery will not be complete until there is less than a 5% chance of hypnosis. Thus, the concentration decrement defining recovery is the difference between C95and C5. This concentration decrement may be remarkably different for drug A and drug B (figure 8). The concentration decrement for drug B is so much smaller than that for drug A that recovery from drug B could be faster even if the kinetics were slower. We think this hypothetical example illustrates the importance of understanding the shape or scale of the concentration–effect relationship and its midpoint. This point can be illustrated using known data for midazolam and propofol. The C50values for deep sedation (asleep, unresponsive to commands, able to be aroused with shoulder tap or loud verbal stimulation) and moderate sedation (asleep but able to be aroused with simple verbal commands) are 208 ng/ml and 101 ng/ml, respectively, for midazolam using the optimal model of Barr et al  . 20 The comparable values for propofol are 0.74 and 0.5 μg/ml. 19 The respective values of γ are 1.7 and 4.5 for propofol and midazolam. One may reasonably postulate that recovery from deep sedation could be defined as the time needed for a decrement in the drug concentration associated with a 95% or better probability of deep sedation to the drug concentration defined by a 5% or less chance of moderate sedation. By using the reported parameter estimates for midazolam, this would require a decrease in plasma concentration from 400 to 50 ng/ml, an eightfold decrease. By using the reported parameter estimates for propofol, this would require a decrease in plasma concentration from 4 to 0.09 μg/ml, a nearly 40-fold decrease. Thus, the relative decrement in concentration defining recovery is much less for midazolam than propofol, although the difference in C50values for the two sedation scores is less for propofol. This concretely shows how the value of γ affects the concentration decrement associated with recovery. However, it must be emphasized that this argument is made to illustrate the importance of the parameter, λ, and not to actually argue that recovery from midazolam is faster than recovery from propofol. The argument could be more concrete only if the estimates of λ are accurate and if the definitions of recovery are well delineated.
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines  (♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line  ) and drug B (lower line  ).
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines 
	(♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line 
	) and drug B (lower line 
	).
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines  (♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line  ) and drug B (lower line  ).
×
In conclusion, the results of simulations, the study of human subjects, and theoretical considerations each indicate that, although estimation of C50is accurate with sparse population data, the estimation of γ is not. Accurate estimation of γ requires dense data sets with 10 data points per patient.
References
Hill AV: The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol 1910; 40: iv–viiHill, AV
Waud DR: On biological assays involving quantal responses. J Pharmacol Exp Ther 1972; 183: 577–607Waud, DR
De Jong RE, Eger EI II: MAC expanded: AD50and AD90values of common inhalation anesthetics in man. A nesthesiology 1975; 42: 384–9De Jong, RE Eger, EI
Fisher DM, Zwass MS: MAC desflurane in 60% nitrous oxide in infants and children. A nesthesiology 1992; 76: 354–6Fisher, DM Zwass, MS
Chortkoff BS, Bennett HL, Eger EI II: Subanesthetic concentrations of isoflurane suppress learning as defined by the category-example task. A nesthesiology 1993; 79: 16–22Chortkoff, BS Bennett, HL Eger, EI
Ausems ME, Hug CC Jr, Stanski DR, Burm AGL: Plasma concentrations of alfentanil required to supplement nitrous oxide anesthesia for general surgery. A nesthesiology 1986; 65: 362–73Ausems, ME Hug, CC Stanski, DR Burm, AGL
Wessen A, Persson PM, Nilsson A, Hartvig P: Concentration-effect relationship of propofol after total intravenous anesthesia. Anesth Analg 1993; 77: 1000–7Wessen, A Persson, PM Nilsson, A Hartvig, P
Vuyk J, Lim T, Engbers FHM, Burm AGL, Vleter AA, Bovill JG: Pharmacodynamics of alfentanil as a supplement to propofol or nitrous oxide for lower abdominal surgery in female patients. A nesthesiology 1993; 78: 1036–45Vuyk, J Lim, T Engbers, FHM Burm, AGL Vleter, AA Bovill, JG
Glass PSA, Doherty M, Jacobs JR, Goodman D, Smith LR: Plasma concentration of fentanyl, with 70% nitrous oxide, to prevent movement at skin incision. A nesthesiology 1993; 78: 842–7Glass, PSA Doherty, M Jacobs, JR Goodman, D Smith, LR
Vuyk J, Lim T, Engbers FHM, Burm AGL, Vletter AA, Bovill JG: The pharmacodynamic interaction of propofol and alfentanil during lower abdominal surgery in women. A nesthesiology 1995; 83: 8–22Vuyk, J Lim, T Engbers, FHM Burm, AGL Vletter, AA Bovill, JG
Jacobs JR, Reves JG, Marty J, White WD, Bai SA, Smith LR: Aging increases pharmacodynamic sensitivity to the hypnotic effects of midazolam. Anesth Analg 1995; 80: 143–8Jacobs, JR Reves, JG Marty, J White, WD Bai, SA Smith, LR
Bailey JM, Schwieger IM, Hug CC Jr: Evaluation of sufentanil anesthesia obtained by a computer-controlled infusion for cardiac surgery. Anesth Analg 1993; 76: 247–52Bailey, JM Schwieger, IM Hug, CC
Lu W, Bailey JM: Reliability of pharmacodynamic analysis by logistic regression: A computer simulation study. A nesthesiology 2000; 92: 985–92Lu, W Bailey, JM
Davidian M, Giltinan DM: Nonlinear Models for Repeated Measurement Data. Boca Raton, Florida, Chapman and Hall, 1995
Beal S, Sheiner L: NONMEM User's Guide. San Francisco, University of California, San Francisco, 1992.
Ramsay MA, Savage TM, Simpson BR, Goodwin R: Controlled sedation with alphaxalone-alphadone. BMJ 1974; 2: 656–9Ramsay, MA Savage, TM Simpson, BR Goodwin, R
Arendt RM, Greenblatt DJ, Garland WA: Quantification by gas chromatography of the 1- and 3-hydroxy metabolites of midazolam in human plasma. Pharmacology 1984; 29: 158–64Arendt, RM Greenblatt, DJ Garland, WA
Somma J, Donner A, Zomorodi K, Sladen R, Ramsay J, Geller E, Shafer SL: Population pharmacodynamics of midazolam administered by target-controlled infusion in SICU patients after CABG surgery. A nesthesiology 1998; 89: 1430–43Somma, J Donner, A Zomorodi, K Sladen, R Ramsay, J Geller, E Shafer, SL
Barr J, Egan TD, Sandovai N, Zomorodi K, Cohane C, Gambus PL, Shafer SL: Propofol dosing regimens for ICU sedation based upon an integrated pharmacokinetic-pharmacodynamic model. A nesthesiology 2001; 95: 324–33Barr, J Egan, TD Sandovai, N Zomorodi, K Cohane, C Gambus, PL Shafer, SL
Barr J, Zomorodi K, Bertaccini EJ, Shafer SL, Geller E: A double-blind randomized comparison of IV lorazepam versus midazolam for sedation of ICU patients via a pharmacologic model. A nesthesiology 2001; 95: 286–98Barr, J Zomorodi, K Bertaccini, EJ Shafer, SL Geller, E
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 1. The bias of simulated C50estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel 
	shows the results for one data point per patient, and the lower panel 
	shows the results for 10 data points per simulated patient.
Fig. 2. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of C50estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
×
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 3. The bias of simulated γ estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel 
	shows the results for one data point per patient, and the lower panel 
	shows the results for 10 data points per simulated patient.
Fig. 4. The coefficient of variation (the square root of the mean squared difference between the true value and the estimate, normalized to the true value) of γ estimates as a function of the number of patients for true γ= 1 (♦), true γ= 1.5 (▵), true γ= 3 (▪), true γ= 4.5 (□), and true γ= 6 (▴). The upper panel  shows the results for one data point per patient, and the lower panel  shows the results for 10 data points per simulated patient.
×
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel 
	shows the results for true γ= 1; the middle panel 
	shows the results for true γ= 3; and the lower panel 
	shows the results for true γ= 6.
Fig. 5. The bias of simulated ETA estimates (the true value minus the mean of the estimate, normalized by the true value) for different numbers of subjects and 1 (□), 3 (▪), 5 (▴), or 10 (▵) data points per patient. The upper panel  shows the results for true γ= 1; the middle panel  shows the results for true γ= 3; and the lower panel  shows the results for true γ= 6.
×
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve  illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve 
	illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
Fig. 6. An illustration of how data can be perfectly fit with an infinite estimate of γ. Two data points are shown, one at a concentration of 75 units/ml, for which a positive drug effect is not observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is observed. The solid curve  illustrates that a model with a single value of C50(100 units/ml) and an infinite value of γ provides a perfect fit to both observations.
×
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves  illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve  ) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves 
	illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve 
	) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
Fig. 7. An illustration of a hypothetical study with two data points, one at a concentration of 75 units/ml, for which a positive drug effect is observed, and the other at a concentration of 125 units/ml, for which a positive drug effect is not observed. The dashed curves  illustrate that a model with an infinite value of γ requires that we postulate two values of C50, one just slightly less than 75 units/ml and the other just slightly greater than 125 units/ml, to explain the observed data. In contrast, a model with a smaller value of γ= 3 (solid curve  ) can provide a fit to the data by postulating a single value of C50equal to100 units/ml. This example shows that if γ is fixed at infinity, NONMEM could only account for the variance of observations in which C is on the wrong side of C50by assuming a higher variance on C50.
×
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines  (♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line  ) and drug B (lower line  ).
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines 
	(♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line 
	) and drug B (lower line 
	).
Fig. 8. The probability of drug effect for two hypothetical drugs each with C50= 100. The value of γ for drug A (the steeper curve) is 6 and the value of γ for drug B is 2. The horizontal lines  (♦) show the concentration decrement needed for a decrease in the probability of drug effect from 95% to 5% for drug A (upper line  ) and drug B (lower line  ).
×