A Technique for Population Pharmacodynamic Analysis of Concentration-Binary Response Data 

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⁵⁰] is the area under a Gaussian curve (with standard deviation of 1) from minus infinity to gamma lnC-gamma lnC⁵⁰(see appendix 1 for details). Thus Equation 2resembles Equation 1(logistic regression) in that it is characterized by two parameters, gamma and C⁵⁰, 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⁵⁰, 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⁵⁰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⁵⁰>, the mean value of C⁵⁰in the population; gamma; and omega², the variance of In C⁵⁰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⁵⁰, 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²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⁵⁰, value and a gamma value were assigned to each hypothetical “patient” using the Excel random-number generator and assuming that both C⁵⁰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⁵⁰and gamma values. For all simulations we assumed that the mean C⁵⁰(denoted <C⁵⁰>) 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⁵⁰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⁵⁰, (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⁵⁰, gamma, and epsilon and the sum gamma lnC-gamma lnC⁵⁰+ epsilon, which determined R. In one of the simulations, patient 6 was assigned a drug concentration of 80 units/ml. The values of C⁵⁰, 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⁵⁰, gamma, and epsilon returned by random-number generation were 115, 5.4, and 1.4, respectively. For this patient, gamma lnC-gamma lnC⁵⁰+ 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²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⁵⁰), 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⁵⁰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⁵⁰. For sufentanil, omega was directly calculated from reported values of C⁵⁰. [16]

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⁵⁰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⁵⁰>, gamma, and omega for two of these simulations. For Figure 2, the simulated data were generated assuming that <C⁵⁰> = 100, gamma = 10, and omega = 1. Qualitatively, this represents steep individual concentration-response curves with substantial interpatient variability in the C⁵⁰value. For Figure 3, the target values are <C⁵⁰> = 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⁵⁰>, 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⁵⁰⁹> = 100, so that the “true” curve can also be illustrated (this is necessary because individual C⁵⁰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).

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⁵⁰variability. Although many pharmacodynamic studies concentrate on estimating C⁵⁰, 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⁵⁰is important, then we should devote the same effort to estimating gamma, because it will determine C⁹⁰, or C⁹⁵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⁵⁰>, gamma, and omega. <C⁵⁰> is the mean C⁵⁰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²is a measure of interpatient variability, the variance of ln C⁵⁰, 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⁵⁰. 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⁵⁰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⁵⁰values (in this case, a standard deviation for ln C⁵⁰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⁵⁰>)/radical(1 + omega²gamma²). 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⁵⁰> because x =(gamma lnC - gamma ln<C⁵⁰>)/radical(1 + omega²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⁵⁰>)/radical(1 + omega²gamma²) 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⁵⁰, 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⁵⁰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⁵⁰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⁵⁰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⁵⁰. 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⁵⁰. Accounting for gamma variability would add a parameter to the model (two additional parameters would be required if covariance between C⁵⁰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²naive” approach, logistic regression, is a good estimator of <C⁵⁰> (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⁵⁰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⁵⁰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⁵⁰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⁵⁰and gamma. If these assumptions are plausible, our simulations indicate that this technique can accurately estimate <C⁵⁰>, 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⁵⁰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⁵⁰; that is, ln C⁵⁰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⁵⁰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⁵⁰, a parameter gamma that influences the steepness of the concentration-response curve, and omega², the variance of lnC⁵⁰in the population. The model assumes a log-normal distribution of C⁵⁰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.

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⁵⁰and gamma), we can arbitrarily assume that the logistic distribution of epsilon has a mean of zero and a standard deviation of pi²/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⁵⁰). 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⁵⁰(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⁵⁰has a log-normal distribution; that is, C⁵⁰= <C⁵⁰> exp (deltaⁱ) where <C⁵⁰> is the mean C⁵⁰value and deltaⁱis normally distributed with a mean of zero and in unknown variance omega sup 2 and is independent of epsilonⁱ. 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ⁱis a normally distributed variable with variance equal to 1 +(gamma omega)². Because we observe a drug effect if yⁱ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⁵⁰>, the mean value of C⁵⁰for individual patients, gamma a measure of the steepness of the concentration-response curve, and omega², the variance of C⁵⁰values in the population.

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⁵⁰> 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⁵⁰>, 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⁵⁰>, gamma, and omega to have positive values. This option is available in the Solver menu.