Search

In PK/PD data analysis, statistical models involving random variables are developed. At times an analysis can be carried out by conditioning on certain random events involving these variables. This paper attempts to clarify issues regarding conditioning. In particular, conditioning is examined as it relates to a number of disparate practical matters: missing covariate values, dose titration, BQL data, "no change from baseline" data, and the use of a truncated intraindividual probability distribution for PK observations.

A model can be fit to joint PK/PD data (concentration and effect) either simultaneously or sequentially. The results of a companion paper suggested that when the data-analytic and true models agree, a particular sequential approach is computationally faster than the simultaneous one, yet produces hardly less precise PD parameter estimates, and for suitable designs, about as accurate PD standard error estimates. In this paper, we compare the performance of various methods for the case that the data-analytic model is misspecified. We illustrate these methods by applying them to a set of real data. Using NONMEM, population PK/PD observations were simulated under various study designs according to a one- or two-compartment PK model and direct Emax or sigmoid Emax model. A one-compartment PK model and Emax PD model were fit to the simulated observations by simultaneous and sequential methods. Predictive performance (interpolation and extrapolation) of PD and the type-I error rate of a likelihood ratio test are compared. The real data set consists of PK and (more frequent) PD observations after administration of the muscle relaxant vecuronium. When only the PK data-analytic model is misspecified, the simultaneous method has greater precision than the sequential methods. However a sequential method that uses a non-parametric PK model performs better than both other methods when PK model misspecification is severe. When the PD data-analytic model is misspecified, sequential and simultaneous methods perform similarly. The analysis of the real data shows that the PK fitted with the simultaneous method can be quite sensitive to PD model misspecification, yielding a possible diagnostic for this type of misspecification.

Dose [-concentration]-effect relationships can be obtained by fitting a predictive pharmacokinetic (PK)-pharmacodynamic (PD) model to both concentration and effect observations. Either a model can be fit simultaneously to all the data (“simultaneous” method), or first a model can be fit to the PK data and then a model can be fit to the PD data, conditioning in some way on the PK data or on the estimates of the PK parameters (“sequential” method). Using simulated data, we compare the performance of the simultaneous method with that of three sequential method variants with respect to computation time, estimation precision, and inference. Using NONMEM, under various study designs, observations of one type of PK and one type of PD response from different numbers of individuals were simulated according to a one-compartment PK model and direct Emax PD model, with parameters drawn from an appropriate population distribution. The same PK and PD models were fit to these observations using simultaneous and sequential methods. Performance measures include computation time, fraction of cases for which estimates are successfully obtained, precision of PD parameter estimates, precision of PD parameter standard error estimates, and type-I error rates of a likelihood ratio test. With the sequential method, computation time is less, and estimates are more likely to be obtained. Using the First Order Conditional Estimation (FOCE) method, a sequential approach that conditions on both population PK parameter estimates and PK data, estimates PD parameters and their standard errors about as well as the “gold standard” simultaneous method, and saves about 40% computation time. Type-I error rates of likelihood ratio test for both simultaneous and sequential approaches are close to the nominal rates.

When modeling new data with a complex population pharmacokinetic/pharmacodynamic model, there may not be sufficient information to obtain estimates of all parameters. In this case information from previous studies can also be used to help stabilize estimation. Using simulated data, we explored three different ways to do this. (i) Some parameter values were fixed to estimates obtained from earlier data. (ii) The earlier data were combined with the current data. (iii) The objective function based on the current data was augmented by a penalty function expressing summary information obtained from the earlier data. This last method is similar to the use of a Bayesian prior. It may be particularly useful when either the combined data set of method (ii) is very large and leads to large computation times or when the early data are not readily available. With this method, two different types of penalty functions were used. With our examples, the three methods all resulted in stabilized estimation. Methods (ii) and (iii) gave similar results for parameter and standard error estimation, especially with respect to fixed effects parameters. For hypothesis testing, results obtained with method (i) are very problematic. There are also problems with the results obtained with method (iii), but they are much less severe, and when the design for the earlier data is known, they can be corrected by using a computer-intensive simulation test procedure.

The posterior predictive check (PPC) is a model evaluation tool. It assigns a value (p PPC ) to the probability that the value of a given statistic computed from data arising under an analysis model is as or more extreme than the value computed from the real data themselves. If this probability is too small, the analysis model is regarded as invalid for the given statistic. Properties of the PPC for pharmacokinetic (PK) and pharmacodynamic (PD) model evaluation are examined herein for a particularly simple simulation setting: extensive sampling of a single individual's data arising from simple PK/PD and error models. To test the performance characteristics of the PPC, repeatedly, “real” data are simulated and for a variety of statistics, the PPC is applied to an analysis model, which may (null hypothesis) or may not (alternative hypothesis) be identical to the simulation model. Five models are used here: (PK1) mono-exponential with proportional error, (PK2) biexponential with proportional error, (PK2e) biexponential with additive error, (PD1) E max model with additive error under the logit transform, and (PD2) sigmoid E max model with additive error under the logit transform. Six simulation/analysis settings are studied. The first three, (PK1/PK1), (PK2/PK2), and (PD1/PD1) evaluate whether the PPC has appropriate type-I error level, whereas the second three (PK2/PK1), (PK2e/PK2), and (PD2/PD1) evaluate whether the PPC has adequate power. For a set of 100 data sets simulated/analyzed under each model pair according to a stipulated extensive sampling design, the p PPC is computed for a number of statistics in three different ways (each way uses a different approximation to the posterior distribution on the model parameters). We find that in general; (i) The PPC is conservative under the null in the sense that for many statistics, prob(p PPC ≤α)

Over the past 25 years sophisticated data analytic techniques have been developed which can lead to improved analyses, but at additional computational cost. In particular, this applies to the approach where interindividual random effects are included in a data analytic model for population pharmacokinetic data, which can often lead to substantially improved estimates of fixed-effect parameters. However, there are also commonly occurring situations, notably with some types of pharmacodynamic data, where such improvement is not realized. This study simulates some simple population dichotomous data, and secondarily, some related continuous data. These data are analyzed using both mixed-effect (ME) models that include interindividual random effects and naive (NA) models that do not include interindividual random effects, and it is seen that use of an ME model does not inevitably lead to gains over use of an NA model. In fact, using maximum likelihood estimation with both types of models, the root mean square estimation errors for fixed effect parameters can actually be larger with an ME model than with the corresponding NA model. Using a form of restricted maximum likelihood estimation with the ME model, the two types of models yield root mean square errors which are comparable, but which still do not suggest that there is always marked advantage in using the ME model.

Pharmacokinetic data consist of drug concentration measurements, as well as reports of some measured concentrations being below the quantification limit of the assay (BQL). A pharmacokinetic model may befit to these data, and for this purpose, the BQL observations must be either discarded or handled in a special way. In this paper, seven methods for dealing with BQL observations are evaluated. Both single-subject and population data are simulated from a one-compartment model. A moderate amount of data is simulated for each individual. The actual cv of concentration measurements at the quantification limit is assumed to be no greater than 20%, in accord with the FDA Guidance. The results of this paper should be interpreted in this context. The methods include handling BQL observations as fixed-point censored observations, i.e., by using the likelihoods that these observations are in fact BQL. This method is shown to have some overall statistical advantage. However, the gain in using this method over that of simply discarding the BQL observations is not always much, and this is especially so when the frequency of BQL observations is small. Some simple methods entailing (i) replacing one or more BQL observations with the value 0, or (ii) replacing them with the value QL/2, where QL is the quantification limit, are also included. The first of these two approaches should not be used With population data, use of the second approach can result in some noticeably improved estimation of the typical value of a parameter, but then there is also marked degradation in the estimation of the population variance of the parameter.

In pharmacokinetic data analysis, it is frequently necessary to select the number of exponential terms in a polyexponential expression used to describe the concentration-time relationship. The performance characteristics of several selection criteria, the Akaike Information Criterion (AIC), and the Schwarz Criterion (SC), and the F test (alpha = 0.05), were examined using Monte Carlo simulations. In particular, the ability of these criteria to select the correct model, to select a model allowing estimation of pharmacokinetic parameters with small bias and good precision, and to select a model allowing precise predictions of concentration was evaluated. To some extent interrelationships among these procedures is explainable. Results indicate that the F test tends to choose the simpler model more often than does either the AIC or SC, even when the more complex model is correct. Also, the F test is more sensitive to deficient sampling designs. Clearance estimates are generally very robust to the choice of the wrong model. Other pharmacokinetic parameters are more sensitive to model choice, particularly the apparent elimination rate constant. Prediction of concentrations is generally more precise when the correct model is chosen. The tendency for the F test (alpha = 0.05) to choose the simpler model must be considered relative to the objectives of the study.

The influence of the choice of pharmacokinetic model on subsequent determination of covariate relationships in population pharmacokinetic analysis was studied using both simulated and real data sets. Simulations and data analysis were both performed with the program NONMEM. Data were simulated using a two-compartment model, but at late sample times, so that preferential selection of the two-compartment model should have been impossible. A simple categorical covariate acting on clearance was included. Initially, on the basis of a difference in the objective function values, the two-compartment model was selected over the one-compartment model. Only when the complexity of the one-compartment model was increased in terms of the covariate and statistical models was the difference in objective function values of the two structural models negligible. For two real data sets, with which the two-compartment model was not selected preferentially, more complex covariate relationships were supported with the one-compartment model than with the two-compartment model. Thus, the choice of structural model can be affected as much by the covariate model as can the choice of covariate model be affected by the structural model; the two choices are interestingly intertwined. A suggestion on how to proceed when building population pharmacokinetic models is given.

In a recent paper, Beal (1989, Biometrics 45, 969-977) considers the problem of determining the appropriate sample size when inference about a parameter theta is to be made on the basis of a confidence interval (CI). He suggests that the sample size should be chosen so that the probability that the length of the CI is less than a given value, conditional on the interval including the true theta, is greater than a specified level. In this note, in which we concentrate on two-sided intervals, this suggestion is examined, as is the effect of uncertainty in our knowledge of the population variance sigma 2 on estimates of sample size.

SUMMARY Use of the maximum likelihood method for parameter estimation in a general mixed effects model requires initial estimates for all the model parameters. When the model is complicated, it may not be easy to obtain such estimates, especially estimates of the variances and covariances. A computational strategy is given for obtaining initial estimates. It is based on a general formulation of the method of moments for mixed effects models which allows the variances and covariances to be estimated for fixed values of the regression parameters. A number of numerical algorithms has been proposed for obtaining maximum, or restricted maximum, likelihood estimates of the parameters of a linear mixed effects model (Hemmerle & Hartley, 1973; Jennrich & Sampson, 1976; Corbeil & Searle, 1976; Laird & Ware, 1982). The maximum likelihood method can also be used with a model in which some fixed effects enter nonlinearly, while all random effects enter linearly; see ? 2. With any of the algorithms it is important to have good initial estimates for all the model parameters. One computational strategy sometimes used for finding initial estimates is to first maximize the full likelihood or log likelihood function, referred to here as the objective function, over a mesh, or grid, of points covering the parameter space and then to use the mesh point corresponding to the maximum as the initial estimate. However, with a complicated mixed effects model the dimension of the parameter space can be large enough to make this strategy impractical. Also, this strategy implies that practical bounds can be placed on all variances and covariances of the model, and though this may be possible, it may not be easy to do, especially when the model is complicated. This paper presents a modification of this strategy. It makes use of a general formulation of the method of moments for mixed effects models which is applicable to obtaining estimates of the variances and covariances of such models for fixed values of the regression parameters. Thus initial estimates of the variances-covariances are obtained without recourse to a grid search over at least those dimensions of the parameter space corresponding to the variances and covariances. The method of moments formulation permits estimates to be obtained under various constraints which might be imposed on the variances and covariances. It is also applicable to multivariate mixed effects models, and

Individual pharmacokinetic parameters quantify the pharmacokinetics of an individual, while population pharmacokinetic parameters quantify population mean kinetics, interindividual kinetic variability, and residual variability, including intraindividual variability and measurement error. Individual pharmacokinetics are estimated by fitting a pharmacokinetic model to individual data. Population pharmacokinetic parameters have traditionally been estimated by doing this separately for each individual, and then combining the individual parameter estimates, the Standard Two Stage (STS) approach. Another approach, NONMEM, appropriately pools data across individuals and is therefore less dependent on individual parameter estimates. This study provides further evidence of NONMEM's validity and usefulness by comparing both approaches on simulated routine-type pharmacokinetic data arising from a monoexponential model. The estimates of population parameters (notably those describing interindividual variability) provided by the STS method are poorer than those provided by NONMEM, especially when there is considerable residual error. Further, NONMEM's estimates of population parameters do not require that the data be restricted to special types of routine data such as those obtained only at steady state, or only at peak or trough, nor do the estimates improve with such data. NONMEM's estimates do improve, however, when a data set is enhanced by the addition of single-observation-per-individual type data. Thus, population parameters can be estimated efficiently from data that simulate real clinical pharmacokinetic conditions.

The stochastic formulation of linear kinetic models is elaborated in order to introduce some new concepts and help clarify the meaning and role of residence time moments. Certain conditional moments are introduced. Multicompartment and steady-state dosing within the stochastic context are considered. A general model-independent formula for steady state volume of distribution and a new concept of steady-state moments are presented. A technique for constructing a model of a given topology from its moments is also given.

Consider a finite population, called here the target population, and consider two conditions, each of which we may imagine applying to all the units of this population. Suppose that under each condition, each population unit will exhibit some binary (0-1) response. The problem is to "compare" the set of responses of all the units under one condition with the set of responses of all the units under the second condition. We call these two sets the two population responses. In order to make the comparison, an experiment, in which some responses are actually obtained, must be performed. The population may be large enough to make it impractical to obtain an entire population response, and consequently, part of the experiment may entail determining a sample of units whose response under one or both conditions will be obtained. Also, in a real experiment it may be impossible in general to actually apply both conditions to a unit. (The actual application of a condition to a particular unit may result in altering the unit so that the other condition cannot be applied to it later). Consequently part of the experiment may entail choosing the conditions which shall be applied to particular experimental units. A common experimental design in these situations involves 1) obtaining a random sample from the target population, 2) randomly allocating the sample units to two experimental groups, and 3) applying one condition to the units of one group, while applying the other condition to the units of the other group. Even when both conditions may actually be applied to all units, this design may still be helpful since both groups may undergo experimentation during the same time period, and therefore, total experiment time may be reduced. Fisher's well known hypergeometric test (2,3,8) or a test which is approximate to Fisher's test, such as the chi-squared test with Yates correction (1,2,3,8), is usually used to analyze the responses obtained from an experiment with the above design. When the sample size, N, is sufficiently small compared to the size, M, of the target population, one may regard the responses of the first (second) group as being approximately a sequence of Bernoulli 0-1 outcomes, with probability m1!M (m2!M) for outcome 1, where mi = the number of I's in the ith population response. Thus Fisher's test may be regarded to be approximately testing the null hypothesis m1!M = m2!M (or ml = M2). In this article we explain that the appropriateness of using Fisher's test to test m1 = m2 may be examined without recourse to the Bernoulli approximation. Such examination is especially important when NIM is so large that it is not really meaningful to talk about Bernoulli outcomes. An important requirement of the inferential procedure is that the sample be random. The usual assumption that the sample is indeed random is often rather shaky in application. This is often the case in medical experiments where an accessible group of patients is used for a sample, rather than a group of patients chosen from either a well-defined target population or sampling frame, and where the inference is about some population of interest which is characterized after the sample is determined and which should be identified here with the target population. For a proper interpretation of the experimental result, one may need to redefine the target population in such a way that the randomness assumption becomes more tenable (or, as it is said, in such a way that the sample becomes more representative of the population). In so doing, one may well add restrictions to the original definition of the target population, thereby decreasing the population size. In fact, it may be necessary to define the target population to be, in effect, the set of units comprising the sample; the sample is then undoubtedly random. Of course, when the population is made less encompassing, it may lose some interest value with respect to the conclusion to be drawn about it. However, this loss is often necessary for a valid inference to result. While constricting the population, the ratio of sample size to population size may increase to the point where the inferential argument should not depend on the use of the Bernoulli approximation. These considerations motivated this author's inquiry into the role of Fisher's test in the inference problem to compare the population responses.

It has previously been shown that the extended least squares (ELS) method for fitting pharmacokinetic models behaves better than other methods when there is possible heteroscedasticity (unequal error variance) in the data. Confidence intervals for pharmacokinetic parameters, at the target confidence level of 95%, computed in simulations with several pharmacokinetic and error variance models, using a theoretically reasonable approximation to the asymptotic covariance matrix of the ELS parameter estimator, are found to include the true parameter values considerably less than 95% of the time. Intervals with the ordinary least squares method perform better. Two adjustments to the ELS confidence intervals, taken together, result in better performance. These are: (i) apply a bias correction to the ELS estimate of variance, which results in wider confidence intervals, and (ii) use confidence intervals with a target level of 99% to obtain confidence intervals with actual level closer to 95%. Kineticists wishing to use the ELS method may wish to use these adjustments.

A statistical model is set forth that can be taken as a very general description of most data sets that arise from population pharmacokinetic studies. A (nontraditional) method for estimating the parameters of the model—called the NONMEM method in previously published papers—is described. This method involves a linearization of the model. An investigation of the effect of this linearization is reported. This investigation is an empirical one, based on the simulation of data from special cases of the model and the application of a few different estimation methods to these data. From the limited evidence in this investigation it appears that the linearization per se does not significantly adversely affect the estimates.

It is often difficult to specify weights for weighted least squares nonlinear regression analysis of pharmacokinetic data. Improper choice of weights may lead to inaccurate and/or imprecise estimates of pharmacokinetic parameters. Extended least squares nonlinear regression provides a possible solution to this problem by allowing the incorporation of a general parametric variance model. Weighted least squares and extended least squares analyses of data from a simulated pharmacokinetic experiment were compared. Weighted least squares analysis of the simulated data, using commonly used weighting schemes, yielded estimates of pharmacokinetic parameters that were significantly biased, whereas extended least squares estimates were unbiased. Extended least squares estimates were often significantly more precise than were weighted least squares estimates. It is suggested that extended least squares regression should be further investigated for individual pharmacokinetic data analysis.

An explicit solution to the Michaelis-Menten differential equation with bolus and zero-order input is presented. This solution involves certain simple functions whose values are not readily available. Efficient algorithms for computing values of these functions to any desired degree of accuracy and FORTRAN codes for implementing them are also given.

The precision of pharmacokinetic parameter estimates from several least squares parameter estimation methods are compared. The methods can be thought of as differing with respect to the way they weight data. Three standard methods, Ordinary Least Squares (OLS-equal weighting), Weighted Least Squares with reciprocal squared observation weighting [WLS(y-2)], and log transform OLS (OLS(ln))--the log of the pharmacokinetic model is fit to the log of the observations--are compared along with two newer methods, Iteratively Reweighted Least Squares with reciprocal squared prediction weighting (IRLS,(f-2)), and Extended Least Squares with power function "weighting" (ELS(f-xi)--here xi is regarded as an unknown parameter). The values of the weights are more influenced by the data with the ELS(f-xi) method than they are with the other methods. The methods are compared using simulated data from several pharmacokinetic models (monoexponential, Bateman, Michaelis-Menten) and several models for the observation error magnitude. For all methods, the true structural model form is assumed known. Each of the standard methods performs best when the actual observation error magnitude conforms to the assumption of the method, but OLS is generally least perturbed by wrong error models. In contrast, WLS(y-2) is the worst of all methods for all error models violating its assumption (and even for the one that does not, it is out performed by OLS(ln)). Regarding the newer methods, IRLS(f-2) improves on OLS(ln), but is still often inferior to OLS. ELS(f-xi), however, is nearly as good as OLS (OLS is only 1-2% better) when the OLS assumption obtains, and in all other cases ELS(f-xi) does better than OLS. Thus, ELS(f-xi) provides a flexible and robust method for estimating pharmacokinetic parameters.

Often drug dosage may be chosen rationally by use of plasma concentration (CP) as the "therapeutic" end point. The ability to accurately forecast the CP resulting from a dosage regimen is central to choosing that regimen. Tradionally forecasting has been attempted only by accounting for known influences on pharmacokinetics, such as sex, age, and renal disease. One must also adjust for previously observed CPs. Herein, we discuss and explain an approach to both of these tasks, mainly focusing on the latter. The approach balances observed outcomes against prior expectations taking account of observation CP error. For digoxin, use of 1 measured CP, as opposed to none, improves forecast precision for future CPs by 40% (decrement in variance of forecast error), and 2 CPs improve it by 67%. There is also an increase in forecast accuracy (decrement in mean of forecast error) as the number of CPs used increases. After only 2, forecast accuracy and precision are as good as theoretically possible. Moreover, information from CPs is far more valuable for forecasting than that from observable patient features-sex, age, and the like; use of all the latter information does not improve accuracy and precision as much as only 1 CP.

An “explicit” solution to the Michaelis-Menten differential equation is presented. Instead of involving the exponential function, the solution involves another simple function. A table for this function is presented.