Many pharmacologic studies record data as binary, yes-or-no, variables with analysis using logistic regression. In a previous study, it was shown that estimates of C50, the drug concentration associated with a 50% probability of drug effect, were unbiased, whereas estimates of gamma, the term describing the steepness of the concentration-effect relationship, were biased when sparse data were naively pooled for analysis. In this study, it was determined whether mixed-effects analysis improved the accuracy of parameter estimation.
Pharmacodynamic studies with binary, yes-or-no, responses were simulated and analyzed with NONMEM. The bias and coefficient of variation of C50 and gamma estimates were determined as a function of numbers of patients in the simulated study, the number of simulated data points per patient, and the "true" value of gamma. In addition, 100 sparse binary human data sets were generated from an evaluation of midazolam for postoperative sedation of adult patients undergoing cardiac surgery by random selection of a single data point (sedation score vs. midazolam plasma concentration) from each of the 30 patients in the study. C50 and gamma were estimated for each of these data sets by using NONMEM and were compared with the estimates from the complete data set of 656 observations.
Estimates of C50 were unbiased, even for sparse data (one data point per patient) with coefficients of variation of 30-50%. Estimates of gamma were highly biased for sparse data for all values of gamma greater than 1, and the value of gamma was overestimated. Unbiased estimation of gamma required 10 data points per patient. The coefficient of variation of gamma estimates was greater than that of the C50 estimates. Clinical data for sedation with midazolam confirmed the simulation results, showing an overestimate of gamma with sparse data.
Although accurate estimations of C50 from sparse binary data are possible, estimates of gamma are biased. Data with 10 or more observations per patient is necessary for accurate estimations of gamma.
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. 1This 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. 13We 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. 13Excel (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. 14This was performed by using the software package NONMEM (University of California, San Francisco, San Francisco, CA) version 5. 15The 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. 17The 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 . 18This 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 γ.
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.
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 γ.
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.
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. 13In 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. 18When 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, 13whereas 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).
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. 19administered 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. 20Although the number of observations is not cited, inspection of figure 1of reference 20clearly 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. 13The agreement of the estimate of γ from our mixed-effects analysis and the naively pooled data analysis of Barr et al . 19suggests 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 . 20The comparable values for propofol are 0.74 and 0.5 μg/ml. 19The 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.
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.