Abstract
Nonlinear Bayesian regression permits curve fitting to a group of subjects simultaneously rather than individually. We evaluated this approach for interpreting plasma clearance curves with the goal of reducing curve-fitting failures and dealing objectively with problem datasets that may arise in clinical settings. Methods: 99mTc-Diethylenetriaminepentaacetic acid plasma clearance curves from 79 subjects were analyzed. The data typically comprised 7–9 samples obtained from 5–10 to 180–240 min after injection. A 2-compartment model was fitted by Bayesian regression to yield compartmental hyperparameters V1, L21, and L12 corresponding to the volume of the compartment into which tracer was injected and the transfer rates from compartment 1 to compartment 2 and from compartment 2 to compartment 1, respectively. This also yielded a clearance estimate for each subject. Results: Estimated hyperparameters were V1 = 8.9 L, L21 = 0.026 min−1, and L12 = 0.040 min−1. Conventional methods led to fitting failures in 2 of the 79 subjects but there were no failures with the Bayesian method. The hyperparameters were used to calculate the glomerular filtration rate for each subject from a single plasma sample with a root-mean-square error of 7.3 mL/min, which was not significantly different from the widely used Christensen—Groth formula. Conclusion: Fewer fitting failures were encountered than with conventional methods, offering an objective means of dealing with problem data. This conceptually simple model can be used directly to calculate clearance from a single plasma sample. It requires only the 3 parameters described above, whereas the Christensen—Groth method requires 6 parameters.
To calculate tracer clearance from a plasma time—activity curve, the curve is commonly fitted numerically to a sum of exponentials. However, this mathematic problem is ill conditioned or ill posed, meaning that small errors in the data can lead to large errors in calculated quantities such as renal clearance (1,2). This can lead to negative parameter estimates and can cause convergence failure with iterative curve-fitting algorithms. This problem has been handled traditionally by manual intervention, but with no assurance of a correct result. Modern nonlinear regression programs better reveal the problem. Anyone who handles many plasma clearance curves by modern methods will encounter occasional curve-fitting problems.
When a problem is ill conditioned, one seeks more information than the data alone can provide. The observed data must be supplemented by additional information that restricts the set of possible solutions. As an example, when filtered backprojection is used in computerized tomography, conditioning is improved by means of a smoothing filter that imposes a smoothness restriction on the solution. In the case of plasma clearance curves, conditioning can be improved by restricting compartmental parameters to physiologically reasonable ranges, which can be regarded as a kind of smoothing. Bayesian regression automatically applies such a restriction and, thus, can improve conditioning.
Bayesian regression treats all subjects simultaneously, in contrast to the usual practice of fitting each subject separately (3). Parameters for individual subjects are regarded as drawn from statistical frequency distributions, each described by its own hyperparameters. For example, the volume of the first compartment for a given subject, a conventional parameter, is regarded as a random variable, and the median of this random variable is a hyperparameter that describes the subject population as a whole. The solution space is restricted, which improves the conditioning, by the implicit requirement that parameters for individual subjects cluster near the group averages. The clustering need not be tight. For example, the relative SD for the volume of the first compartment was found to be 22% (Table 1). This imposes only a weak restriction.
Compartmental Hyperparameter Estimates
In 1985, we showed that the parameters of a 2-compartment model could be used directly to estimate clearance from a single plasma sample (4). This direct method attracted less attention than other methods that were simpler to program and that were presented in the same publication. However, the elegance and conceptual simplicity of dealing with just the 3 parameters of the conventional model were attractive. With our data, it worked as well as other methods.
In 1986, Christensen and Groth (5) reported a 1-sample method based on the 2-compartment model that incorporated an empiric formula for extracellular fluid volume, providing for variation in body size. That has become the standard 1-sample method. It fits our own data better than our 1985 formula and slightly (but not significantly) better than even the method reported in this article. But it lacks theoretic elegance: It is a 6-parameter method, none of which correspond to the parameters of the standard model, and pediatric use requires a completely different set of parameters.
Was the Christensen—Groth method better than our 1985 method because of the scaling for patient size implicit in use of the extracellular volume? That is possible, but we have never found scaling for size to improve significantly the results of a 1-sample method unless children were included (6,7). More likely, our estimates of the compartmental parameters were not good enough. Until now, our attempts to obtain better estimates of these parameters have been defeated by ill conditioning of the curve-fitting problem. (With conventional formulation of the problem, different combinations of the 3 parameters could lead to nearly the same residual variance.) With the Bayesian method, we have now found parameters for the simple 2-compartment model that yield clearance estimates rivaling the best empiric methods.
MATERIALS AND METHODS
Patient Population
Volunteers were recruited from adult patients referred for evaluation of renal function to the nuclear medicine clinics at the University of Alabama in Birmingham (63 subjects) or Emory University (16 subjects). Informed consent was obtained in accord with institutionally approved procedures. Patients were recruited to include a broad range of renal function from low to high to ensure that the model would fit the full range encountered clinically.
Plasma Clearance Curves
After administration of an imaging dose of 185 MBq 99mTc-diethylenetriaminepentaacetic acid (DTPA), 7–9 blood samples were obtained from an indwelling catheter over a time interval beginning 5–10 min after injection and continuing for 180–240 min. The detailed methods have been described (4). All data were corrected for protein-bound impurities either directly (4,8) or by counting only the plasma ultrafiltrate (9). If significant dose infiltration was seen on images of the injection site, data from that subject were excluded.
Scaling for Patient Size
The scaling used here has been used by physiologists for cross-species scaling of physiologic parameters (e.g., plasma volume and glomerular filtration rate [GFR] from mouse to elephant): volume parameters are scaled by body weight and flux parameters (clearance, intercompartmental flows) are scaled by body surface area. We have discussed this scaling method and shown its validity for human subjects elsewhere (6). (Other investigators, reasoning from the fact that measurement of extracellular fluid volume in children tends to correlate better with surface area than with weight, have obtained useful results scaling volumes by surface area rather than weight (10). We have been unable, thus far, to show a statistically significant difference between these 2 scaling methods in the context of clearance estimation.) Good results with the scaling method we used earlier (6) are reported in this article.
Bayesian Statistical Model
Overview.
A conventional linear 2-compartment model was assumed, with injection into and excretion from compartment 1 (having virtual volume V1) and with transfer of activity from compartment 1 to compartment 2 (having transfer coefficient L21) and in the reverse direction (having transfer coefficient L12). This is only a model, in which the volumes are virtual volumes that may not have meaningful physiologic correlates. (The physiologically meaningful volume is the sum of V1 and V2, which approximates the extracellular fluid volume.)
The 2-compartment model was fitted by Bayesian regression. The values of V1, L21, and L12 were allowed to vary from subject to subject. (Actually, the logarithms of these parameters were used to ensure that negative values would not arise in the curve-fitting procedure. For nonnegative physiologic quantities having large variance, such logarithmic transformation frequently leads to approximately normal distributions.) The clearance for each subject was not treated as a random variable but as a unique parameter describing that subject. The likelihood function for this nonlinear model was derived in the manner described by Gelman et al. (3) for the closely related linear model. A noninformative uniform prior distribution was used except for the variance parameters, which were restricted to a physiologically reasonable range. To obtain hyperparameter estimates from the posterior distribution (which is simply the product of the likelihood function and the prior distribution), the mode of the posterior distribution was calculated using the program of Bunch et al. (11).
Equations.
Consider a time—activity curve from a 2-compartment model with N plasma samples designated by index i. At sample time ti, let the observed activity be yi and the predicted activity fi. Each fi is a function of the sample time ti and the vector β of compartmental parameters having M components βj. (In this case M = 3, and the components of vector β are the 3 parameters V1, L12, and L12.) Assume that yi = fi + εi, where the measurement errors εi are normally distributed with mean 0 and variance σ2. This is the standard nonlinear regression model for a compartmental model, so that fi can be calculated from ti and β by standard methods.
Now extend this model from a single subject to a group of K subjects, each designated by index k, and extend the notation by replacing each subscript i by ik. Assume that βjk = β0j + γjk, where b0j is constant and γjk is normally distributed with mean 0 and variance σj2. (Here j indexes a component of the parameter vector, and k indexes a subject.) Following the reasoning given by Gelman et al. (3) for the closely related linear regression model leads one to the following expression for the likelihood function:
where ln(lik) is the log likelihood, δyik is a residual (the difference yik − fik between observed and predicted values), and δβjk is the difference βjk − β0j between βjk and its expected value β0j. Calculation of the residuals δyik, which involves calculating predicted activities from the parameters, can be done by various methods for this standard linear compartmental model. We used a matrix diagonalization method in order not to limit the computer program to 2 compartments, although only 2 compartments are used in this article (12).
The Bayesian posterior distribution is the product of the above likelihood function and the Bayesian prior distribution. We assumed a noninformative uniform prior distribution (simply a factor of 1) for each parameter other than the variance parameters, and so these factors can be ignored. The variance parameters (squares of the SDs) were treated differently because simple noninformative prior distributions for variance parameters place too much mass at the origin when using hierarchic models. (Some parameters may shrink to 0 when finding the mode numerically.) As suggested by Gelman et al. (3), the prior distribution for each variance parameter was taken as a scaled inverse χ2 distribution that embraced a physically reasonable range for the variance parameters. Thus,
where f(θ) is the probability density function for parameter θ, ν and s are parameters of the distribution, and Γ is the mathematic γ-function. The variance hyperparameter for measuring plasma activity was forced into a range of approximately 1%–10% (expressed as relative SD), a range chosen from experience with similar measurements, by selecting values of 3 and 0.0004 for ν and s, respectively. The variance hyperparameters for compartmental parameters were forced into a range of approximately 10%–40% by selecting values of 7 and 0.04 for ν and s, respectively. The posterior distribution was obtained by multiplying the above likelihood function by the product of these inverse χ2 factors, one for each variance hyperparameter. The compartmental hyperparameters were then estimated by finding the mode of the posterior distribution numerically.
Conventional Statistical Model
As in our previous work (7,13), weighted regression was used on the basis of the premise that measurement error arises chiefly from volumetric laboratory manipulations. Measurement error was thus assumed proportional to the measured activity. (For the data presented in this article, unweighted regression was found to increase the number of fitting failures and would have led to an even more favorable comparison for the Bayesian model.)
RESULTS
The estimates of the Bayesian hyperparameters are shown in Table 1. Because logarithms of the parameters were used, Table 1 presents the median and relative SD (relative to the measured value), which take simple forms under the logarithmic transformation.
Although conventional methods led to fitting failures in 2 subjects, no failures occurred with the Bayesian model so that all 79 subjects were included in the results shown in Table 1. The residual SD attributable to measurement error was 3.4%, close to the value of 3% determined by direct measurement of pipeting error (13). For the 77 subjects for whom the data could be fitted by conventional methods, the root-mean-square (rms) difference between Bayesian and conventional clearance estimates was small, only 3.4 mL/min.
The parameter values and variances indicated above correspond to the joint mode of the posterior distribution. (Note that, from the Bayesian viewpoint, the statistical parameters themselves have a frequency distribution, so that we can speak of the mode of the median, and there is no discrepancy with the headings in Table 1.) A full Bayesian analysis can describe the entire distribution. Thus, we were able to apply the Metropolis algorithm (3) to obtain a pseudorandom sample from the posterior distribution, but the size of the problem (441 parameters, including latent parameters) meant difficulty in scaling and very slow convergence, measured in days. The sample medians calculated in this way for V1 (8.3 L), L21 (0.033 min−1), and L12 (0.048 min−1) differed somewhat from the modal estimates shown in Table 1, not surprising when the number of parameters is large in proportion to the number of data points. (Both methods estimated the median of the population, but one was based on the mode of the posterior distribution and the other was based on the median of a sample drawn from the posterior distribution by means of the Metropolis algorithm.) However, clearance estimates for individual subjects were quite precise: For these, the mean interquartile range, or mean interval from 25 percentile to 75 percentile (a measure of error in the clearance measurement), was 3.3 mL/min/1.73 m2. Furthermore, for single-subject clearance, the Metropolis and modal estimates were in close agreement (rms difference, 1.7 mL/min/1.73 m2). We conclude that the modal parameter estimates of Table 1 are satisfactory for calculating clearance. In a larger study, the Metropolis method would not be practical, and so it is fortunate that the simpler modal estimates suffice.
Individualized error estimates can be calculated for the clearance of each subject (from an individual multisample time—activity curve) by means of the Metropolis algorithm. The largest interquartile range for a single subject in this study was 11 mL/min/1.73 m2, nearly 4 times the average value of 3.3 mL/min/1.73 m2. The clearance itself can be calculated quickly, but calculating the error estimate can take as much as 1 h with our present software and is not practical for routine use.
The hyperparameters in Table 1 can also be used as a prior distribution for direct estimation of clearance from 1 or 2 (or more) plasma samples. In the 1-sample case, this is mathematically equivalent to the Newton’s method approach that we have described (4), so that programming is simple and execution is fast. (Only the clearance is unknown because the other single-patient parameters are the same as the group hyperparameters.) The results are shown in Figure 1. We found an rms error of 7.3 mL/min/m2 using a single sample at 180 min and 4.5 mL/min/1.73 m2 using 2 samples at 60 and 180 min. For comparison, using the formula of Christensen and Groth (5) with our data, the rms error was 6.8 mL/min/1.73 m2. The difference between the 2 single-sample methods was not statistically significant (F test; P < 0.05). Another comparison is shown in Figure 2, which shows the GFR calculated by each method versus the effective volume of distribution at 3 h for an adult of normal height and weight. Note the close agreement between the 2 methods.
Correlation of Bayesian 1- and 2-sample methods and Christensen—Groth 1-sample method (5) with multisample reference method. Sampling was at approximately 3 h for 1-sample methods and at 1 and 3 h for 2-sample method.
GFR vs. effective volume (Veff) of distribution at 3 h after injection, for adult of average size, calculated from Bayesian parameter estimates and from Christensen—Groth formula (5).
DISCUSSION
In contrast to the Bayesian method, the conventional method was more susceptible to fitting failures, which occurred in 2 of the 79 subjects. In a research study, such datasets can be deleted, but in the clinic, fitting failure can be a serious problem. The Bayesian method reduced the number of such failures, providing an objective means of dealing with problem datasets.
The method used to scale for height and weight has been used successfully in the past to extrapolate from adult to pediatric data (6), so it is of interest to try such extrapolation with the current model. Figure 3 compares this Bayesian model with pediatric formulas proposed by Groth and Aasted (10) and by Ham and Piepsz (14). GFR as calculated by all 3 methods is plotted against effective volume of distribution at 120 min for a 1-y-old child of average height and weight. (The 2 pediatric methods were derived for 51Cr-ethylenediaminetetraacetic acid [EDTA], not 99mTc-DTPA, but the properties of these agents are known to be similar.) Note the close agreement between the Bayesian 2-compartment model and the Groth—Aasted pediatric formula. The same Bayesian model (with identical parameters) was used for adults and children, unlike the Groth formulas, which differ for adults (Christensen—Groth) (5) and children (Groth—Aasted). Although Figure 3 is highly suggestive, it must be emphasized that this method has not been tested directly in children. The Bayesian model clearly fits published pediatric data for 51Cr-EDTA, but it should be tested with pediatric data for 99mTc-DTPA as well.
For adult subjects this method led to results that did not differ significantly from the standard Christensen—Groth method (5). Note that this model has only 3 adjustable parameters compared with 6 in the Christensen—Groth formula (including the 2 parameters required to calculate extracellular volume). The Christensen—Groth formula cannot be used for children; instead, there is a separate Groth—Aasted formula (10). We suggest direct use of the 2-compartment model with the same 3 parameters for adults and children.
CONCLUSION
A Bayesian 2-compartment model was applied to the plasma clearance of 99mTc-DTPA. Compartmental parameters were estimated by fitting a single global model simultaneously to data from multiple subjects. The result was a conceptually simple 2-compartment model having only 3 adjustable parameters that can be used in lieu of empiric formulas for 1-sample and 2-sample clearance estimation. This approach can also be used for more robust fitting of multisample clearance curves, with fewer fitting failures than conventional methods, to provide a means of dealing with problem datasets.
Footnotes
Received Jul. 17, 2001; revision accepted Jan. 31, 2002.
For correspondence or reprints contact: Charles D. Russell, MD, PhD, Division of Nuclear Medicine, University of Alabama Hospital, 619 S. 19th St., Birmingham, AL 35233.
E-mail: crussell{at}uab.edu
