Minimal model SI 0 problem in NIDDM subjects ... - Semantic Scholar
Am J Physiol Endocrinol Metab 282: E564–E573, 2002; 10.1152/ajpendo.00576.2000.
Minimal model SI⫽0 problem in NIDDM subjects: nonzero Bayesian estimates with credible confidence intervals GIANLUIGI PILLONETTO,1 GIOVANNI SPARACINO,1 PAOLO MAGNI,2 RICCARDO BELLAZZI,2 AND CLAUDIO COBELLI1 1 Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, 35131 Padova; and 2 Dipartimento di Informatica e Sistemistica, Universita` degli Studi di Pavia, 27100 Pavia, Italy Received 14 December 2000; accepted in final form 10 October 2001
Pillonetto, Gianluigi, Giovanni Sparacino, Paolo Magni, Riccardo Bellazzi, and Claudio Cobelli. Minimal model SI⫽0 problem in NIDDM subjects: nonzero Bayesian estimates with credible confidence intervals. Am J Physiol Endocrinol Metab 282: E564–E573, 2002; 10.1152/ajpendo. 00576.2000.—The minimal model of glucose kinetics, in conjunction with an insulin-modified intravenous glucose tolerance test, is widely used to estimate insulin sensitivity (SI). Parameter estimation usually resorts to nonlinear least squares (NLS), which provides a point estimate, and its precision is expressed as a standard deviation. Applied to type 2 diabetic subjects, NLS implemented in MINMOD software often predicts SI⫽0 (the so-called “zero” SI problem), whereas general purpose modeling software systems, e.g., SAAM II, provide a very small SI but with a very large uncertainty, which produces unrealistic negative values in the confidence interval. To overcome these difficulties, in this article we resort to Bayesian parameter estimation implemented by a Markov chain Monte Carlo (MCMC) method. This approach provides in each individual the SI a posteriori probability density function, from which a point estimate and its confidence interval can be determined. Although NLS results are not acceptable in four out of the ten studied subjects, Bayes estimation implemented by MCMC is always able to determine a nonzero point estimate of SI together with a credible confidence interval. This Bayesian approach should prove useful in reanalyzing large databases of epidemiological studies.
in conjunction with the insulin-modified intravenous glucose tolerance test (IM-IVGTT), is widely used to measure insulin sensitivity in subjects with impaired glucose tolerance and type 2 diabetes in both clinical and epidemiological studies (1, 14, 17, 23). Although some of its assumptions, in particular the single compartment approximation to describe glucose kinetics, have been recently reexamined (7, 8), it is important to understand the performance of the classical minimal model when it is identified with the most advanced estimation techniques. Virtually all minimal model
identification strategies employed in the literature resort to nonlinear least squares (NLS) estimation. NLS provides a point estimate of each model parameter and, by means of the Fisher information matrix, a measure of its uncertainty expressed as standard deviation (SD) can be obtained, from which a confidence interval can be determined. However, this approach has difficulties in handling possible asymmetries in the probability distribution of the estimates. In providing the estimates, NLS exploits only the experimental data (e.g., plasma glucose concentration time series in our case) and the knowledge of the measurement error statistics. NLS is thus a purely data-driven approach. One reported problem with this approach in minimal model studies is that, in a number of subjects, especially those with type 2 diabetes, insulin sensitivity (SI) is calculated as SI⫽0 with the minimal model program MINMOD (19). For instance, Saad et al. (23) found that SI⫽0 occurred in 12 of 24 type 2 diabetic subjects, i.e., a prevalence of 50%. As a consequence, in population studies, histograms with a likely artificial bimodal pattern (i.e., a peak at SI⫽0 and another peak at a positive SI value) are obtained. At the same time, with more general purpose software packages, e.g., SAAM II or ADAPT (2, 10), SI is estimated to be nonzero, and thus physically sound, but very small and with a very large uncertainty. As a result, negative values are included in the confidence interval. The above interpretative difficulties, often referred to as the “zero” SI problem, call for the adoption of parameter estimation techniques more sophisticated than NLS. The so-called Bayes approach is a methodology that can be employed to attack the parameter estimation problem as an alternative to NLS, but it is less adopted in practice. Bayesian estimation techniques exploit not only the experimental data but also the a priori information (also denoted in the following by the term “prior,” as is usual in the statistics literature) on the unknown parameters of the model (27). In physiological modeling, this a priori information, e.g., nonnegativity and range of variability, can typically be made available from
Address for reprint requests and other correspondence: C. Cobelli, Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, Via Gradenigo, 6a, 35131 Padova, Italy (E-mail: [email protected]).
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked ‘‘advertisement’’ in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
insulin sensitivity; insulin resistance; mathematical model; parameter estimation; type 2 diabetes
THE MINIMAL MODEL OF GLUCOSE KINETICS,
E564
0193-1849/02 $5.00 Copyright © 2002 the American Physiological Society
http://www.ajpendo.org
E565
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
physical considerations and population studies (see Refs. 8 and 24 for two recent applications of Bayesian estimation in metabolism). The key advantage of Bayes estimation over NLS is that the former returns the entire a posteriori probability distribution of the model parameters, from which the marginal probability distribution of each parameter can be obtained. From this probabilistic quantity, a point estimate and its confidence interval can be computed. Notably, with use of Bayes estimation, confidence intervals are in general allowed to be asymmetric with respect to the point parameter estimates, in contrast to NLS, where the implicit assumption of parameter estimates to be Gaussian leads to confidence intervals symmetric by construction. The limited popularity of Bayes estimation in physiological modeling is likely due to its heavy computational burden, because often the a posteriori probability distribution of the model parameters must be numerically computed by the Markov chain Monte Carlo (MCMC) simulation approach (13). In this article we develop a Bayesian parameter estimation strategy to identify the minimal model of glucose kinetics in diabetic subjects and implement it by MCMC. First, we only supply to the estimator a rough prior, i.e., the unknown model parameters are positive. Then, because more refined a priori information on SI can be made available from the literature, we reidentify the model and show improved results. The database consists of IM-IVGTT performed in ten subjects. Whereas NLS provides SI estimates and confidence intervals that are not acceptable in four cases, Bayes estimation always determines a nonzero estimate of SI and, when resorting to the refined prior model, a credible confidence interval. The Bayesian approach we propose should prove useful to reanalyze large databases of epidemiological studies, such as GENNID (20), the Finland-US investigation of NIDDM genetics (FUSION) (26), and the Insulin Resistance Atherosclerosis Study (IRAS) (17), where the numbers of subjects are in the range of thousands.
MATERIALS AND METHODS
Subjects, Protocol, and Data An IM-IVGTT (glucose dose of 300 mg/kg at time 0 plus 0.05 U/kg of insulin given as a square wave between 20 and 25 min) was performed in 10 type 2 diabetic subjects. The data of 7 subjects have already been published in Ref. 1, to which we refer the reader for details on protocol and measurement. Plasma glucose and insulin concentrations were frequently measured for 4 h. Minimal Model of Glucose Kinetics The minimal model of glucose kinetics (4) during an IVGTT (Fig. 1) is given by ˙ 共t兲 ⫽ ⫺共SG ⫹ X共t兲兲G共t兲 ⫹ SGGb G共0兲 ⫽ G0 G ˙ 共t兲 ⫽ ⫺p2共X共t兲 ⫺ SI(I共t兲 ⫺ Ib兲) X
X共0兲 ⫽ 0
(1) (2)
where G (mg/dl) and I (U/ml) are plasma glucose and insulin concentrations, respectively, with Gb and Ib representing their baseline values; X is remote insulin equal to (k4 ⫹ k6)I(t); SG (min⫺1) is glucose effectiveness equal to k1⫹ k5; SI (min⫺1/U ml⫺1) is insulin sensitivity equal to k2,(k4 ⫹ k5)/k3; p2 (min⫺1) is a rate parameter equal to k3; and G0 is the value of G extrapolated at time 0. The four model parameters SG, SI, p2, and G0 are a priori uniquely identifiable (6, 9) from the glucose data, once I(t) in Eq. 2 is assumed a known input by linearly interpolating the measured plasma concentrations; [neglecting the insulin measurement errors is without significant consequences on parameter estimates, because the key model variable is insulin action X, which is obtained by integrating I (3)]. For what follows, it is useful to consider the unknown parameter vector ⫽ [SG, p2, G0, SI]T. NLS Identification From the model of Eqs. 1 and 2, glucose concentration at time ti, i.e., G(ti), is predicted by a function of the unknown parameter vector , i.e., h(ti,). If we denote by yi the ith glucose plasma concentration noisy measurement, one has y i ⫽ h共ti, 兲 ⫹ vi i ⫽ 1, . . . , N
(3)
where vi is the measurement error, assumed to be additive, independent, and Gaussian with zero mean and a 2% coefficient of variation.
Fig. 1. Cold minimal model of glucose disappearance. NHGB, net hepatic glucose balance; D␦(t), cold glucose dose; Ib, basal insulin dose; G(t), glucose concentration at time t; Q(t), glucose mass at time t; Up(t), peripheral glucose disposal at time t; ki, rate constants characterizing material fluxes (solid lines) or control actions (dashed lines).
AJP-Endocrinol Metab • VOL
282 • MARCH 2002 •
www.ajpendo.org
E566
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 2. Glucose (units at left, interpolated with solid lines) and insulin (units at right, interpolated with dashed lines) plasma concentrations of the 10 subjects.
As in Ref. 1, we estimate, in each of the 10 subjects of our database, the unknown parameter vector from the data vector y ⫽ [y1, y2,. . ., yN]T by weighted NLS, i.e.
兺冉 N
ˆ ⫽ arg min
i⫽1
冊
yi ⫺ h共ti, 兲 wi
2
(4)
where the weights wi are optimally chosen, i.e., equal to the known measurement error SD (6, 9). The precision of each AJP-Endocrinol Metab • VOL
estimate, expressed by its SD, is obtained from the square root of the diagonal elements of the inverse of the Fisher information matrix. Bayesian Identification Bayes estimation is based on the concept of a priori information on the unknown parameter vector , mathematically specified by its a priori probability density function f(). In this context, a priori means “before having seen the data y.”
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Table 1. NLS parameter estimates Subject Nos.
SI, 10⫺4 min⫺1/U ml⫺1
SG, 10⫺2 min⫺1
1 2 3 4 5 6 7 8 9 10
0.92 (20) 1.51 (3) 0.22 (148) 0.86 (23) 2.1 (5) 0.14 (186) 0.98 (353) 0.99 (7.6) 0.52 (128) 0.58 (10)
1.54 (14) 2.1 (10) 1.33 (19) 0.97 (17) 0.87 (29) 1.33 (12) 1.5 (10) 0.9 (16) 1.7 (4.9) 0.84 (19)
p2 , min⫺1
Go , mg/dl
0.015 (21) 297.9 (2.3) 0.072 (19) 381.6 (23) 0.00495 (402) 444.13 (1.3) 0.007 (37) 315.6 (2.2) 0.098 (14) 409 (2.3) 0.0074 (176) 563 (1.6) 0.0016 (460) 315 (1.5) 0.3 (30) 323.8 (1.7) 0.0011 (147) 353.9 (1.4) 0.027 (8) 263 (2.2)
NLS, nonlinear least squares; SI, insulin sensitivity; SG, glucose effectiveness; p2, rate parameter; Go, glucose at time 0. Nos. in parentheses show precision of the estimate expressed as a percent coefficient of variation.
A posteriori, i.e., “after having seen the data y,” the probability density function from which we expect to be taken obviously changes. This function goes under the name of a posteriori probability density function and is denoted by f y(兩y) (the symbol 兩y reads as “ given y”). The a posteriori probability function of contains extremely insightful information. For instance, a point estimate, called the posterior mean of the unknown parameter vector, can be determined as the expected value of the vector 兩y ˆ ⫽ E关円y兴 ⫽
兰
f円y共円y兲d
E567
clared. In formal terms, we state that each value of SG, p2, and Go is a priori equally probable in [0, a] with a3⬁. The same poorly informative prior could obviously also be applied to SI. One can note from Eq. 6 that, in this way, the posterior density function coincides, except for a constant factor, with the likelihood function in the domain of interest. As will be illustrated in RESULTS, this prior, albeit rough, allows a solution of the “zero” SI problem, because a confidence interval, including only nonnegative values, is obtained. However, because one additionally experienced problem with NLS is that the confidence interval of SI tends to be large, it is worthwhile trying to incorporate more knowledge into the SI prior, as we will explain. A refined SI prior. Ideally, to define the a priori probability density function of SI, one would need a large database, e.g., containing thousands of subjects, to extract it directly. Unfortunately, such a database is not available to us. However, several insulin sensitivity NLS studies in diabetic subjects are available in the literature that allow one to obtain an approximation of the actual distribution of the true SI values. From these investigations, one knows that it is unlikely that SI will exceed certain values. For instance, in Ref. 18, 30 type 2 diabetic subjects were studied, and the average SI was 0.7 (10⫺4 min⫺1/U ml⫺1) with a standard error of 0.1. Of note is
(5)
The a posteriori probability function of also allows the derivation of the 95% confidence intervals [L,H], e.g., the interval between the 0.025 and the 0.975 quantile of the a posteriori probability distribution of . The a posteriori probability density function of the parameter vector can be obtained from the Bayes theorem as f 円y共円y兲 ⫽
fy円共y円兲f共兲 fy共y兲
(6)
In the right-hand side of Eq. 6, fy(y) does not depend on and, for us, is just a scale factor we are not interested in (13), f() is given (see the next paragraph), and fy円(y兩) (often referred to as the likelihood function of the data) can be obtained from Eq. 3 by exploiting the knowledge of measurement error statistics. For instance, when, as in our case study, v ⫽ [v1, v2, . . ., vN]T is Gaussian with zero mean and covariance matrix ¥v, the random vector y兩 is Gaussian with covariance matrix ¥v and mean given by [h(t1,), h(t2,), . . ., h(tN,)]T. A nonnegativity prior on the four parameters. Prior information is a key additional ingredient of Bayes estimation, i.e., one needs to specify the a priori probability density function f() of Eq. 6, which summarizes our beliefs on model parameters “before having seen the data.” From Eq. 6 one can note that the Bayesian approach coincides with maximum likelihood, if a uniform distribution from ⫺⬁ and ⫹⬁ is assumed as a prior for all of the unknown parameters. In our minimal model studies, we know in advance that all the parameters (SG, p2, Go, and SI) are intrinsically nonnegative. Therefore, it is reasonable to choose a prior f() that excludes negative parameter values. For SG, p2, and Go, nonnegativity is the only information embedded into the prior distribution, i.e., an expected range of variability is not deAJP-Endocrinol Metab • VOL
Fig. 3. Nonlinear least squares (NLS) and Bayes insulin sensitivity (SI) estimates in 10 type 2 diabetic subjects. Vertical line, point estimates; horizontal bars, 95% confidence intervals; NLS estimates were assumed Gaussian to determine their confidence interval as ⫾1.96 SD; for Bayes estimation, point estimates and confidence intervals were obtained from information depicted in Fig. 5 (see RESULTS, A refined SI prior).
282 • MARCH 2002 •
www.ajpendo.org
E568
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
that, in this data set, no “zero” SI problems were present. Thus one expects that SI values greater than a certain threshold are less and less probable. At the same time, because SI can be very small, although we do not know how small it can be in a certain subject before seeing his or her data, it is wise not to incorporate in the prior any beliefs about its variability at low levels. Thus we have assumed for SI an a priori probability density function in which SI values ⬍2 ⫻ 10⫺4 (min⫺1/Uml⫺1) are equally probable, while SI values ⬎2 ⫻ 10⫺4 are less and less probable, according to a decreasing exponential law (with exponent equal to 10⫺4 U ml⫺1/min⫺1). The chosen threshold value of 2 ⫻ 10⫺4 min⫺1/U ml⫺1, even if not directly obtained from a distribution of SI estimates in NIDDM subjects, is realistic and is supported by the literature, where it is unlikely to find insulin sensitivity values in NIDDM subjects ⬎2 ⫻ 10⫺4 min⫺1/U ml⫺1. The MCMC method. The goal is to obtain f兩y(兩y), i.e., the joint probability density function of the unknown minimal model parameters, given the glucose data. This task is analytically intractable, both because of the complex relationships between parameters and data and because of the shape of the probability distributions involved. A solution suggested in the literature is to derive f兩y(兩y) in sampled form by adopting a simulation strategy known as Markov chain Monte Carlo (13). MCMC methods are based on two steps. In the first step, a suitable Markov chain that converges (in distribution) to a target distribution, f兩y(兩y) in our case, is generated. Then, a Monte Carlo integration step is performed to numerically compute the integrals of interest, e.g., Eq. 5. There are many MCMC methods in the literature; see Ref. 13 for a review. They differ from each other in the way the Markov chain is generated. However, all of the different strategies proposed in the literature can be traced back to the Metropolis-Hastings algorithm (15). In this work, we generated the Markov chain by a symmetric transition kernel that extracts a sample of the chain around the previous one (this scheme is usually referred to as “random-walk Metropolis”) by use of independent uniform probability densities. The variance of such densities was chosen case by case. As far as implementation issues of the MCMC method are concerned, we note that convergence of the chain was assessed by the Raftery criterion (21), which is also known in the literature as binary-control. To compute the number of iterations necessary to estimate a posterior quantile from a single run of a Markov chain, a pilot analysis of output samples was used to fit a two-state Markov chain model, and from it one can calculate the length of the burn-in period and, then, the number of further iterations required to estimate quantiles of interest with the required accuracy. In the sequel, we have
been required to estimate quantiles 0.025, 0.25, 0.5, 0.75, 0.975 with precision factors of 0.02, 0.05, 0.05, 0.05, and 0.02, respectively, with a probability of 0.95. From these quantiles, numerically robust confidence intervals can be calculated. It is worth stressing that, from the practical point of view, switching from NLS to Bayes estimation is not without a price. To give an idea, minutes instead of seconds are required to identify the model in each individual. The reason for this increased time is that Markov chain generation and Monte Carlo integration, needed to handle Eqs. 5 and 6, are computationally demanding. In particular, to obtain the results we present in the next section, a number of iterations varying from 35,000 to 80,000, depending on the subject under study, were performed. RESULTS
Data Plasma glucose and insulin concentrations of the 10 subjects are shown in Fig. 2. NLS Identification Results of NLS for SG, p2,, G0, and SI, obtained with SAAM II (2), are reported in Table 1 together with the precision of the estimates. The SI estimates in all of the subjects are shown in Fig. 3, together with their 95% confidence intervals, obtained by adding and subtracting to the point estimate the quantity 1.96 SD (results of subjects 1–7 are those already reported in Table 2 of Ref. 1). As a consequence of the implicit Gaussian assumption, the confidence interval is symmetric and centered around the point estimate. As apparent from the picture, SI estimates are not equally satisfactory in all of the subjects; whereas in subjects 1, 2, 4, 5, 8, and 10 a narrow confidence interval lying on the nonnegative axis is obtained, in subjects 3, 6, 7, and 9, the low value of the SI estimate, combined with its poor precision, produce an unrealistic confidence interval. In fact, in each of these four subjects, the SI confidence interval tells us that the true SI can fall with nonzero probability in the negative portion of the x-axis. Note that, because the covariance matrix of the measurement error is assumed to be known from the data, NLS coincides with maximum likelihood (ML). If it were assumed to be dependent on the model, then NLS and ML are no longer equivalent on a theoretical basis.
Table 2. Bayesian parameter estimates with the refined SI prior Subject Nos.
SI, 10⫺4 min⫺1/U ml⫺1
SG, 10⫺2 min⫺1
p2, min⫺1
Go , mg/dl
1 2 3 4 5 6 7 8 9 10
1.2 (0.61–2.7) 1.5 (1.4–1.6) 1.7 (0.4–4.1) 0.77 (0.58–1.1) 2.1 (1.8–2.2) 0.39 (0.032–2.3) 1.51 (0.36–4.1) 0.98 (0.8–1.11) 0.8 (0.18–2.7) 0.6 (0.46–0.72)
1.59 (1.07–1.98) 2.1 (1.8–2.38) 1.06 (0.78–1.35) 1.02 (0.76–1.25) 0.89 (0.48–1.32) 1.35 (1–1.62) 1.48 (1.32–1.62) 0.92 (0.66–1.19) 1.61 (1–1.85) 0.84 (0.52–1.18)
0.013 (0.0013–0.033) 0.071 (0.057–0.84) 0.0018 (0.0004–0.0068) 0.008 (0.0033–0.013) 0.098 (0.078–0.12) 0.0070 (0.00014–0.07) 0.0019 (0.00033–0.008) 0.35 (0.2–0.68) 0.006 (0.00034–0.024) 0.025 (0.021–0.03)
299.7 (284.1–313.3) 382 (369.9–394.7) 434.6 (421.9–448.4) 318 (307.4–329.2) 411 (395.6–428.4) 563 (543.3–581.3) 314 (304.8–323.5) 325.4 (315.3–336.2) 350.9 (335.4–363.2) 263 (252.5–275.2)
Nos. in parentheses show 95% confidence intervals. AJP-Endocrinol Metab • VOL
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
E569
Fig. 4. The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a nonnegativity prior on the 4 model parameters. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The y-axis quantity is unitless.
However, in this case, the two results are virtually superimposable (unpublished data). Bayesian Identification A nonnegativity prior on the four parameters. The a posteriori probability density function of SI estimated in each subject is reported in Fig. 4. As somewhat expected, in subjects 1, 2, 4, 5, 8, and 10, where NLS was successful, the posterior provided (in sampled form) by MCMC is concentrated around the NLS estimate. However, the new method reveals asymmetries in the distribution of the estimates, particAJP-Endocrinol Metab • VOL
ularly in subjects 1 and 4. This is not surprising, because the symmetry of NLS confidence intervals reflects only the Gaussian assumption implicitly made on the estimates. Moreover, in subjects 3, 6, 7, and 9, MCMC analysis reveals that the marginal posterior densities have a high peak located at low and realistic SI values but that they have a long tail. The confidence intervals, however, would be too pessimistic, because it is known from the literature that such high SI values are unrealistic, e.g., it is impossible that SI is, say, equal to 100 ⫻ 10⫺4 min⫺1/U ml⫺1 (see subject 3 in Fig. 4). So, to obtain a more
282 • MARCH 2002 •
www.ajpendo.org
E570
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 5. The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a refined SI prior. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The yaxis quantity is unitless.
credible confidence interval, it is necessary to adopt a more realistic SI prior. A refined SI prior. Results for SG, p2, G0, and SI are reported in Table 2, together with the 95% confidence interval. Results for SI are also displayed in Fig. 5, where the estimated a posteriori probability density function of SI in each subject is reported (the axes are now different from those used in Fig. 4), and in Fig. 3, where a summary of Fig. 5 results is depicted to make the comparison with NLS easier. One can note, by AJP-Endocrinol Metab • VOL
comparing Figs. 5 and 4, how the threshold influences the tail of the posterior distribution of some subjects. One can also note that the posterior probability density function is zero in an interval at the right of SI⫽0. This allows us to exclude, on a probabilistic basis, that SI is greater than a chosen threshold. From the a posteriori density function of SI, one can obtain the interval where the true value of SI falls with 95% of confidence as the interval between the quantiles 0.025 and 0.975. This interval is displayed in Fig. 3, together with the
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
posterior mean of SI. The novel features of the approach we propose emerges in subjects 3, 6, 7, and 9, where NLS results were largely unsatisfactory with noncredible confidence intervals. One can see that the marginal a posteriori probability density function of SI does not collapse to zero (this also applies to Fig. 4, but it cannot be visually appreciated due to the larger range of SI values shown in the x-axis). In all of these subjects, the posterior mean of SI is higher than the NLS estimate (markedly in subjects 3 and 7), because the minimum variance estimate takes into account the entire shape of the posterior and, so, also the presence in some subjects of the tails. More important, in each subject the confidence interval now lies entirely in the positive axis and allows the investigator to infer that the probability that the true SI lies at the left of the displayed bar is 2.5%. In summary, a nonzero SI value is detected with high probability by Bayesian identification in these very insulin-resistant individuals. This result is not always possible with NLS, where, given the uncertainty of the SI estimate, one can only conclude that SI is indistinguishable from zero (Fig. 3, e.g., subject 7). A credible confidence interval is obtained, which allows the minimal model user to determine with 97.5% of confidence a (nonzero) lower-bound number for SI. DISCUSSION
In this article, we have attacked the identification of the minimal model of glucose kinetics in type 2 diabetic subjects in a Bayesian context. The motivation for this stems from the interpretative difficulties generated from what is usually referred to as the “zero SI” problem, i.e., in a significant number of cases NLS estimation returns an SI⫽0 when the minimal model program MINMOD (19) is used, or very small values with an unrealistic confidence interval including negative values with software packages like SAAM II (2) or ADAPT (10). In this case study, the Bayesian approach allowed the derivation (albeit only numerically through MCMC) of the joint a posteriori probability distribution of the model parameters, and thus the determination of a nonnegative point estimate of SI, together with a plausible confidence interval. The differences of the point estimates between the Bayesian approach and NLS, at least in most subjects, are not large, but it is when the confidence intervals of the parameter estimates are considered that the Bayesian approach emerges as superior. In fact, because we think that the chosen prior as well as the confidence intervals employed are both reasonable, one can exclude on a probabilistic basis that SI in the critical subjects is higher than a given threshold. Moreover, from the confidence interval obtained by Bayes estimation, one can exclude with 97.5% of probability that the true SI lies at the left of the displayed variability range. In summary, the Bayesian approach is able to cope with those situations in which the confidence interval, obtained as two times the SD of the error given by the Fisher information matrix, is not realistic, because such a method is unAJP-Endocrinol Metab • VOL
E571
able to cope with asymmetries in the distribution of the estimates, like those present in this study. A critical point in Bayes estimation is the definition with a probability density function of the information available on the unknown model parameters before seeing the data. As far as SG, p2, and G0 are concerned, we have limited ourselves to express the poorly informative knowledge that they are nonnegative, i.e., a uniform distribution along the positive axis was chosen. A more sophisticated prior was chosen for SI, the very critical parameter of the model, based on the observation that, for values of SI greater than a certain threshold, the higher the candidate SI, the lower, in all likelihood, its probability. We have translated these beliefs about SI in statistical terms by describing its a priori probability density function as a uniform distribution between zero and a certain threshold (i.e., 2 ⫻ 10⫺4 min⫺1/U ml⫺1, a value that appears as prudential on the basis of published studies on insulin resistance in type 2 diabetics) and an exponential decay for SI greater than this threshold. It is important to note that, in the prior chosen for SI, we do not include beliefs on how small it can be, i.e., where the whole “zero SI” problem began. Choosing also for SI a uniform distribution prior, as in Bayesian identification with only the nonnegativity prior, one would have obtained an SI a posteriori probability density function lower than the one presented. This means that low SI values are less probable and supports the robustness in each subject of the 0.025 quantile of Fig. 1, which allows us to exclude with 97.5% of probability that the true SI lies at the left of the displayed range. Another consequence of using a uniform prior for SI along the positive axis is represented by a gratuitous shift to the right of the 0.975 quantile. In fact, the new technique permits us to show that, in some cases, the model defines likelihood functions which assign a nonnegligible probability to high SI values that are largely unrealistic. This kind of information was impossible to obtain so clearly by the standard Fisherian approaches. The rationale of choosing for SI a prior more sophisticated than that employed for SG, p2, and G0 was to determine with 95% confidence an interval as narrow as possible where SI lies. The role played by this a priori knowledge is not hidden but transparent in terms of its effect on the results. For example, a strategy to obtain nonnegative confidence intervals by NLS could consist in reparameterizing the model with the aim of estimating the logarithm of the parameters. However, this would be obtained by approximating the likelihood as a function of the logarithm of the parameters, with a Gaussian distribution. This is questionable, because it would not be possible to evaluate the quality of such approximation. As already mentioned, switching from NLS to Bayes estimation is not without a price, because the MCMC results are computationally more demanding to obtain. However, this strategy offers a new tool that could also become important in classification studies. In fact, it would be interesting to reanalyze large databases of epidemiological studies in which the standard NLS
282 • MARCH 2002 •
www.ajpendo.org
E572
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 6. A: a qualitative distribution of SI estimates that would be obtained in a study in which the SI⫽0 problem is present. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The y-axis shows the no. of SI estimates falling in equally spaced containers. B: qualitative distribution of SI estimates that would be obtained by the Markov chain Monte Carlo (MCMC) method in the same study.
approach has produced SI histograms with a likely artificial bimodal pattern (with one peak at SI⫽0 and another peak at a positive SI value), like that qualitatively displayed in Fig. 6, which are difficult to interpret. The value SI⫽0, which Fisherian techniques return as the joint mode of a suitable objective function, can be a mathematical artifact, because it can be far from the minimum variance estimate provided by MCMC. Therefore, we believe that if the new strategy proposed in this paper were applied in epidemiological studies like those mentioned above, one would obtain a more reliable distribution of SI estimates, e.g., changing the qualitative profile of Fig. 6A into that shown in Fig. 6B. In this way, more reliable information on NIDDM would be made available. Moreover, a reanalysis of databases such as GENNID (20), FUSION (26), and IRAS (16, 17), where the numbers of subjects are in the range of thousands, could also allow a refinement of our MCMC method, because a more accurate description of the a priori probability density function of SI will be possible. It is clear, in fact, that the refined prior used in Fig. 5 provides better results than those obtained with the positivity prior and shown in Fig. 4. However, the refined prior still has a certain degree of roughness. A future evolution of the method could also envisage a population approach with the aim of estimating, jointly with model parameters, the common prior that has generated them. Finally, it is worth stressing that important questions concerning possible problems introduced by the structure of the considered model still remain open also in the light of our MCMC results. In fact, one can clearly note that the a posteriori probability density functions of SI of subjects 3 and 7 are quite close to the AJP-Endocrinol Metab • VOL
prior. As a consequence, the point estimates are close to the mean of the prior, equal to 1.67 10⫺4 min⫺1/U ml⫺1, suggesting that sometimes the model extracts little information on SI from IVGTT data (in fact the posterior distribution approaches the prior as the data contain decreasing information about the parameters). Moreover, the presence of long tails in the a posteriori probability density of SI of some subjects is a critical issue that will have to be clarified in future studies. Our Bayesian analysis, including reasonable a priori information on SI, improves the estimation process, reducing this kind of problem. However, it is clear that no single parameter estimation technique, however sophisticated, can be a panacea, able to compensate possible structural defects of a model. So, further investigations are necessary to clarify whether the cause of these problems is the fact that the model is structurally inadequate for the estimation of a low SI. We are inclined to believe that, in many of such critical situations, MCMC is really a more robust estimator than NLS or ML and can significantly improve the estimation process. However, we also believe that in the future it will also be worth exploring a different description of the same process, like the two-compartment model described in Refs. 7 and 8, which may overcome these problems. Maybe, after including some suitable a priori information, which is necessary because of a priori identifiability issues, the problems mentioned here may be solved by resorting to a maximum a posteriori estimator. In any case, it will be interesting to identify the two-compartment model by resorting to the same MCMC approach described here. In fact, even if more computationally demanding, this technique is likely to provide a clearer picture of the
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
advantages of the two- vs. the one-compartment description. This work was in part supported by National Institutes of Health Grant RR-12609, by the MURST COFIN 2000 project, “Estimation of nonaccessibile parameters in physiological systems,” by a University of Padova grant (STIM-PET) to G. Sparacino, and by a University of Pavia grant to P. Magni.
15. 16.
17.
REFERENCES 1. Avogaro A, Vicini P, Valerio A, Caumo A, and Cobelli C. The hot but not the cold minimal model allows precise assessment of insulin sensitivity in NIDDM subjects. Am J Physiol Endocrinol Metab 270: E532–E540, 1996. 2. Barrett PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P, and Foster D. SAAM II: simulation, analysis and modeling software for tracer and pharmacokinetic studies. Metabolism 47: 484–492, 1998. 3. Bergman RN, Bortolan G, Cobelli C, and Toffolo G. Identification of a minimal model of glucose disappearance for estimating insulin sensitivity. In: Proc IFAC Symposium Identification and System Parameter Estimation. Darmstadt, Germany: IFAC, 1979, vol. 2, p. 883–890. 4. Bergman RN, Finegood DT, and Ader M. Assessment of insulin sensitivity in vivo. Endocrinol Rev 6: 45–46, 1985. 5. Bergman RN, Ider YZ, Bowden CR, and Cobelli C. Quantitative estimation of insulin sensitivity. Am J Physiol Endocrinol Metab Gastrointest Physiol 236: E667–E677, 1979. 6. Carson ER, Cobelli C, and Finkelstein L. The Mathematical Modelling of Metabolic and Endocrine Systems. New York: Wiley, 1983. 7. Caumo A, Vicini P, Zachwieja JJ, Avogaro A, Yarasheski K, Bier DM, and Cobelli C. Undermodeling affects minimal model indexes: insights from a two-compartment model. Am J Physiol Endocrinol Metab 276: E1171–E1193, 1999. 8. Cobelli C, Caumo A, and Omenetto M. Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian two-compartment model. Am J Physiol Endocrinol Metab 277: E481–E488, 1999. 9. Cobelli C, Foster D, and Toffolo G. Tracer Kinetic in Biomedical Research: from Data to Model. New York: Kluwer Academic/ Plenum, 2000. 10. D’Argenio D and Schumitzky A. ADAPT II Users Guide: Pharmacokinetic/Pharmacodynamic Systems Analysis Software. Los Angeles, CA: Biomedical Simulations Resource, University of Southern California, 1997. 11. DeFronzo RA, Bonadonna RC, and Ferranini E. Pathogenesis of NIDDM: a balanced overview. Diabetes Care 15: 318–368, 1992. 12. DeFronzo RA, Tobin JD, and Andres R. Glucose clamp technique: a method for quantifying insulin secretion and resistance. Am J Physiol Endocrinol Metab Gastrointest Physiol 237: E214–E223, 1979. 13. Gilks WR, Richardson S, and Spiegelhalter DJ. Markov Chain Monte Carlo in Practice. London: Chapman & Hall, 1996. 14. Haffner SM, D’Agostino R Jr, Mykkanen L, Tracy R, Howard B, Rewers M, Selby J, Savage PJ, and Saad MF. Insulin sensitivity in subjects with type 2 diabetes. Relationship
AJP-Endocrinol Metab • VOL
18. 19.
20.
21.
22. 23.
24.
25. 26.
27. 28.
E573
to cardiovascular risk factors: the Insulin Resistance Atherosclerosis Study. Diabetes Care 22: 562–568, 1999. Hastings WK. Monte Carlo sampling methods using Markov chain and their applications. Biometrika 57: 97–109, 1970. Howard G, Bergman R, Wagenknecht LE, Haffner SM, Savage PJ, Saad MF, Laws A, and D’Agostino RB Jr. Ability of alternative indices of insulin sensitivity to predict cardiovascular risk: comparison with the ’minimal model.’ Ann Epidemiol 8: 358–369, 1998. Mykkanen L, Zaccaro DJ, Hales CN, Festa A, and Haffner SM. The relation of proinsulin and insulin to insulin sensitivity and acute insulin response in subjects with newly diagnosed type II diabetes: the Insulin Resistance Atherosclerosis Study. Diabetologia 42: 1060–1066, 1999. Owens DR, Luzio SD, and Coates PA. Insulin secretion and sensitivity in newly diagnosed NIDDM Caucasians in the UK. Diabetes Med 13: S19–S24, 1996. Pacini G and Bergman RN. MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsively from the frequently sampled intravenous glucose tolerance test. Comput Methods Prog Biomed 23: 113–122, 1986. Raffel LJ, Robbins DC, Norris JM, Boerwinkle E, DeFronzo RA, Elbein RC, Fujimoto W, Hanis CL, Kahn SE, Permutt MA, Chiu KC, Cruz J, Ehrmann DA, Robertson RP, Rotter JI, and Buse J. The GENNID study. A resource for mapping the genes that cause NIDDM. Diabetes Care 19: 864– 872, 1996. Raftery AE and Lewis SM. Implementing MCMC. In: Markov Chain Monte Carlo in Practice, edited by WR Gilks, S Richardson, and DJ Spiegelhalter. London: Chapman & Hall, 1996, p. 115–130. Reaven GM. Banting lecture 1988. Role of insulin resistance in human disease. Diabetes 37: 1595–1607, 1988. Saad MF, Anderson RL, Laws A, Watanabe RM, Kades WW, Chen YDI, Sands RE, Pei D, Savage PJ, and Bergman RN. A comparison between the minimal model and the glucose clamp in the assessment of insulin sensitivity across the spectrum of glucose tolerance. Diabetes 43: 1114–1121, 1994. Sparacino G, Tombolato C, and Cobelli C. Maximum-likelihood vs. maximum a posteriori parameter estimation of physiological system models: the C-peptide impulse response case study. IEEE Trans Biomed Eng 47: 801–811, 2000. Tierney L. Markov chains for posterior distributions (with discussion). Ann Statistics 22: 1701–1762, 1994. Valle T, Tuomilehto J, Bergman RN, Gosh S, Hauser ER, Eriksson J, Nylund SJ, Kohtamaki K, Toivanen L, Vidfren G, Tuomilehto-Wolf E, Ehnolm C, Blaschal J, Langefeld CD, Watanabe RM, Magnuson V, Ally DS, Hagopian WA, Ross E, Buchanan TA, Collins F, and Boehnke M. Mapping genes for NIDDM design of the Finland-United States investigation of NIDDM genetics (FUSION) study. Diabetes Care 21: 949–956, 1998. Walter E and Pronzato L. Identification of Parametric Models from Experimental Data. Paris: Springer, 1997. Welch S, Gebhart SS, Bergman RN, and Phillips LS. Minimal model analysis of intravenous glucose tolerance test-derived insulin sensitivity in diabetic subjects. J Clin Endocrinol Metab 71: 1508–1518, 1990.
282 • MARCH 2002 •
www.ajpendo.org
Minimal model SI⫽0 problem in NIDDM subjects: nonzero Bayesian estimates with credible confidence intervals GIANLUIGI PILLONETTO,1 GIOVANNI SPARACINO,1 PAOLO MAGNI,2 RICCARDO BELLAZZI,2 AND CLAUDIO COBELLI1 1 Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, 35131 Padova; and 2 Dipartimento di Informatica e Sistemistica, Universita` degli Studi di Pavia, 27100 Pavia, Italy Received 14 December 2000; accepted in final form 10 October 2001
Pillonetto, Gianluigi, Giovanni Sparacino, Paolo Magni, Riccardo Bellazzi, and Claudio Cobelli. Minimal model SI⫽0 problem in NIDDM subjects: nonzero Bayesian estimates with credible confidence intervals. Am J Physiol Endocrinol Metab 282: E564–E573, 2002; 10.1152/ajpendo. 00576.2000.—The minimal model of glucose kinetics, in conjunction with an insulin-modified intravenous glucose tolerance test, is widely used to estimate insulin sensitivity (SI). Parameter estimation usually resorts to nonlinear least squares (NLS), which provides a point estimate, and its precision is expressed as a standard deviation. Applied to type 2 diabetic subjects, NLS implemented in MINMOD software often predicts SI⫽0 (the so-called “zero” SI problem), whereas general purpose modeling software systems, e.g., SAAM II, provide a very small SI but with a very large uncertainty, which produces unrealistic negative values in the confidence interval. To overcome these difficulties, in this article we resort to Bayesian parameter estimation implemented by a Markov chain Monte Carlo (MCMC) method. This approach provides in each individual the SI a posteriori probability density function, from which a point estimate and its confidence interval can be determined. Although NLS results are not acceptable in four out of the ten studied subjects, Bayes estimation implemented by MCMC is always able to determine a nonzero point estimate of SI together with a credible confidence interval. This Bayesian approach should prove useful in reanalyzing large databases of epidemiological studies.
in conjunction with the insulin-modified intravenous glucose tolerance test (IM-IVGTT), is widely used to measure insulin sensitivity in subjects with impaired glucose tolerance and type 2 diabetes in both clinical and epidemiological studies (1, 14, 17, 23). Although some of its assumptions, in particular the single compartment approximation to describe glucose kinetics, have been recently reexamined (7, 8), it is important to understand the performance of the classical minimal model when it is identified with the most advanced estimation techniques. Virtually all minimal model
identification strategies employed in the literature resort to nonlinear least squares (NLS) estimation. NLS provides a point estimate of each model parameter and, by means of the Fisher information matrix, a measure of its uncertainty expressed as standard deviation (SD) can be obtained, from which a confidence interval can be determined. However, this approach has difficulties in handling possible asymmetries in the probability distribution of the estimates. In providing the estimates, NLS exploits only the experimental data (e.g., plasma glucose concentration time series in our case) and the knowledge of the measurement error statistics. NLS is thus a purely data-driven approach. One reported problem with this approach in minimal model studies is that, in a number of subjects, especially those with type 2 diabetes, insulin sensitivity (SI) is calculated as SI⫽0 with the minimal model program MINMOD (19). For instance, Saad et al. (23) found that SI⫽0 occurred in 12 of 24 type 2 diabetic subjects, i.e., a prevalence of 50%. As a consequence, in population studies, histograms with a likely artificial bimodal pattern (i.e., a peak at SI⫽0 and another peak at a positive SI value) are obtained. At the same time, with more general purpose software packages, e.g., SAAM II or ADAPT (2, 10), SI is estimated to be nonzero, and thus physically sound, but very small and with a very large uncertainty. As a result, negative values are included in the confidence interval. The above interpretative difficulties, often referred to as the “zero” SI problem, call for the adoption of parameter estimation techniques more sophisticated than NLS. The so-called Bayes approach is a methodology that can be employed to attack the parameter estimation problem as an alternative to NLS, but it is less adopted in practice. Bayesian estimation techniques exploit not only the experimental data but also the a priori information (also denoted in the following by the term “prior,” as is usual in the statistics literature) on the unknown parameters of the model (27). In physiological modeling, this a priori information, e.g., nonnegativity and range of variability, can typically be made available from
Address for reprint requests and other correspondence: C. Cobelli, Dipartimento di Elettronica e Informatica, Universita` degli Studi di Padova, Via Gradenigo, 6a, 35131 Padova, Italy (E-mail: [email protected]).
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked ‘‘advertisement’’ in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
insulin sensitivity; insulin resistance; mathematical model; parameter estimation; type 2 diabetes
THE MINIMAL MODEL OF GLUCOSE KINETICS,
E564
0193-1849/02 $5.00 Copyright © 2002 the American Physiological Society
http://www.ajpendo.org
E565
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
physical considerations and population studies (see Refs. 8 and 24 for two recent applications of Bayesian estimation in metabolism). The key advantage of Bayes estimation over NLS is that the former returns the entire a posteriori probability distribution of the model parameters, from which the marginal probability distribution of each parameter can be obtained. From this probabilistic quantity, a point estimate and its confidence interval can be computed. Notably, with use of Bayes estimation, confidence intervals are in general allowed to be asymmetric with respect to the point parameter estimates, in contrast to NLS, where the implicit assumption of parameter estimates to be Gaussian leads to confidence intervals symmetric by construction. The limited popularity of Bayes estimation in physiological modeling is likely due to its heavy computational burden, because often the a posteriori probability distribution of the model parameters must be numerically computed by the Markov chain Monte Carlo (MCMC) simulation approach (13). In this article we develop a Bayesian parameter estimation strategy to identify the minimal model of glucose kinetics in diabetic subjects and implement it by MCMC. First, we only supply to the estimator a rough prior, i.e., the unknown model parameters are positive. Then, because more refined a priori information on SI can be made available from the literature, we reidentify the model and show improved results. The database consists of IM-IVGTT performed in ten subjects. Whereas NLS provides SI estimates and confidence intervals that are not acceptable in four cases, Bayes estimation always determines a nonzero estimate of SI and, when resorting to the refined prior model, a credible confidence interval. The Bayesian approach we propose should prove useful to reanalyze large databases of epidemiological studies, such as GENNID (20), the Finland-US investigation of NIDDM genetics (FUSION) (26), and the Insulin Resistance Atherosclerosis Study (IRAS) (17), where the numbers of subjects are in the range of thousands.
MATERIALS AND METHODS
Subjects, Protocol, and Data An IM-IVGTT (glucose dose of 300 mg/kg at time 0 plus 0.05 U/kg of insulin given as a square wave between 20 and 25 min) was performed in 10 type 2 diabetic subjects. The data of 7 subjects have already been published in Ref. 1, to which we refer the reader for details on protocol and measurement. Plasma glucose and insulin concentrations were frequently measured for 4 h. Minimal Model of Glucose Kinetics The minimal model of glucose kinetics (4) during an IVGTT (Fig. 1) is given by ˙ 共t兲 ⫽ ⫺共SG ⫹ X共t兲兲G共t兲 ⫹ SGGb G共0兲 ⫽ G0 G ˙ 共t兲 ⫽ ⫺p2共X共t兲 ⫺ SI(I共t兲 ⫺ Ib兲) X
X共0兲 ⫽ 0
(1) (2)
where G (mg/dl) and I (U/ml) are plasma glucose and insulin concentrations, respectively, with Gb and Ib representing their baseline values; X is remote insulin equal to (k4 ⫹ k6)I(t); SG (min⫺1) is glucose effectiveness equal to k1⫹ k5; SI (min⫺1/U ml⫺1) is insulin sensitivity equal to k2,(k4 ⫹ k5)/k3; p2 (min⫺1) is a rate parameter equal to k3; and G0 is the value of G extrapolated at time 0. The four model parameters SG, SI, p2, and G0 are a priori uniquely identifiable (6, 9) from the glucose data, once I(t) in Eq. 2 is assumed a known input by linearly interpolating the measured plasma concentrations; [neglecting the insulin measurement errors is without significant consequences on parameter estimates, because the key model variable is insulin action X, which is obtained by integrating I (3)]. For what follows, it is useful to consider the unknown parameter vector ⫽ [SG, p2, G0, SI]T. NLS Identification From the model of Eqs. 1 and 2, glucose concentration at time ti, i.e., G(ti), is predicted by a function of the unknown parameter vector , i.e., h(ti,). If we denote by yi the ith glucose plasma concentration noisy measurement, one has y i ⫽ h共ti, 兲 ⫹ vi i ⫽ 1, . . . , N
(3)
where vi is the measurement error, assumed to be additive, independent, and Gaussian with zero mean and a 2% coefficient of variation.
Fig. 1. Cold minimal model of glucose disappearance. NHGB, net hepatic glucose balance; D␦(t), cold glucose dose; Ib, basal insulin dose; G(t), glucose concentration at time t; Q(t), glucose mass at time t; Up(t), peripheral glucose disposal at time t; ki, rate constants characterizing material fluxes (solid lines) or control actions (dashed lines).
AJP-Endocrinol Metab • VOL
282 • MARCH 2002 •
www.ajpendo.org
E566
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 2. Glucose (units at left, interpolated with solid lines) and insulin (units at right, interpolated with dashed lines) plasma concentrations of the 10 subjects.
As in Ref. 1, we estimate, in each of the 10 subjects of our database, the unknown parameter vector from the data vector y ⫽ [y1, y2,. . ., yN]T by weighted NLS, i.e.
兺冉 N
ˆ ⫽ arg min
i⫽1
冊
yi ⫺ h共ti, 兲 wi
2
(4)
where the weights wi are optimally chosen, i.e., equal to the known measurement error SD (6, 9). The precision of each AJP-Endocrinol Metab • VOL
estimate, expressed by its SD, is obtained from the square root of the diagonal elements of the inverse of the Fisher information matrix. Bayesian Identification Bayes estimation is based on the concept of a priori information on the unknown parameter vector , mathematically specified by its a priori probability density function f(). In this context, a priori means “before having seen the data y.”
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Table 1. NLS parameter estimates Subject Nos.
SI, 10⫺4 min⫺1/U ml⫺1
SG, 10⫺2 min⫺1
1 2 3 4 5 6 7 8 9 10
0.92 (20) 1.51 (3) 0.22 (148) 0.86 (23) 2.1 (5) 0.14 (186) 0.98 (353) 0.99 (7.6) 0.52 (128) 0.58 (10)
1.54 (14) 2.1 (10) 1.33 (19) 0.97 (17) 0.87 (29) 1.33 (12) 1.5 (10) 0.9 (16) 1.7 (4.9) 0.84 (19)
p2 , min⫺1
Go , mg/dl
0.015 (21) 297.9 (2.3) 0.072 (19) 381.6 (23) 0.00495 (402) 444.13 (1.3) 0.007 (37) 315.6 (2.2) 0.098 (14) 409 (2.3) 0.0074 (176) 563 (1.6) 0.0016 (460) 315 (1.5) 0.3 (30) 323.8 (1.7) 0.0011 (147) 353.9 (1.4) 0.027 (8) 263 (2.2)
NLS, nonlinear least squares; SI, insulin sensitivity; SG, glucose effectiveness; p2, rate parameter; Go, glucose at time 0. Nos. in parentheses show precision of the estimate expressed as a percent coefficient of variation.
A posteriori, i.e., “after having seen the data y,” the probability density function from which we expect to be taken obviously changes. This function goes under the name of a posteriori probability density function and is denoted by f y(兩y) (the symbol 兩y reads as “ given y”). The a posteriori probability function of contains extremely insightful information. For instance, a point estimate, called the posterior mean of the unknown parameter vector, can be determined as the expected value of the vector 兩y ˆ ⫽ E关円y兴 ⫽
兰
f円y共円y兲d
E567
clared. In formal terms, we state that each value of SG, p2, and Go is a priori equally probable in [0, a] with a3⬁. The same poorly informative prior could obviously also be applied to SI. One can note from Eq. 6 that, in this way, the posterior density function coincides, except for a constant factor, with the likelihood function in the domain of interest. As will be illustrated in RESULTS, this prior, albeit rough, allows a solution of the “zero” SI problem, because a confidence interval, including only nonnegative values, is obtained. However, because one additionally experienced problem with NLS is that the confidence interval of SI tends to be large, it is worthwhile trying to incorporate more knowledge into the SI prior, as we will explain. A refined SI prior. Ideally, to define the a priori probability density function of SI, one would need a large database, e.g., containing thousands of subjects, to extract it directly. Unfortunately, such a database is not available to us. However, several insulin sensitivity NLS studies in diabetic subjects are available in the literature that allow one to obtain an approximation of the actual distribution of the true SI values. From these investigations, one knows that it is unlikely that SI will exceed certain values. For instance, in Ref. 18, 30 type 2 diabetic subjects were studied, and the average SI was 0.7 (10⫺4 min⫺1/U ml⫺1) with a standard error of 0.1. Of note is
(5)
The a posteriori probability function of also allows the derivation of the 95% confidence intervals [L,H], e.g., the interval between the 0.025 and the 0.975 quantile of the a posteriori probability distribution of . The a posteriori probability density function of the parameter vector can be obtained from the Bayes theorem as f 円y共円y兲 ⫽
fy円共y円兲f共兲 fy共y兲
(6)
In the right-hand side of Eq. 6, fy(y) does not depend on and, for us, is just a scale factor we are not interested in (13), f() is given (see the next paragraph), and fy円(y兩) (often referred to as the likelihood function of the data) can be obtained from Eq. 3 by exploiting the knowledge of measurement error statistics. For instance, when, as in our case study, v ⫽ [v1, v2, . . ., vN]T is Gaussian with zero mean and covariance matrix ¥v, the random vector y兩 is Gaussian with covariance matrix ¥v and mean given by [h(t1,), h(t2,), . . ., h(tN,)]T. A nonnegativity prior on the four parameters. Prior information is a key additional ingredient of Bayes estimation, i.e., one needs to specify the a priori probability density function f() of Eq. 6, which summarizes our beliefs on model parameters “before having seen the data.” From Eq. 6 one can note that the Bayesian approach coincides with maximum likelihood, if a uniform distribution from ⫺⬁ and ⫹⬁ is assumed as a prior for all of the unknown parameters. In our minimal model studies, we know in advance that all the parameters (SG, p2, Go, and SI) are intrinsically nonnegative. Therefore, it is reasonable to choose a prior f() that excludes negative parameter values. For SG, p2, and Go, nonnegativity is the only information embedded into the prior distribution, i.e., an expected range of variability is not deAJP-Endocrinol Metab • VOL
Fig. 3. Nonlinear least squares (NLS) and Bayes insulin sensitivity (SI) estimates in 10 type 2 diabetic subjects. Vertical line, point estimates; horizontal bars, 95% confidence intervals; NLS estimates were assumed Gaussian to determine their confidence interval as ⫾1.96 SD; for Bayes estimation, point estimates and confidence intervals were obtained from information depicted in Fig. 5 (see RESULTS, A refined SI prior).
282 • MARCH 2002 •
www.ajpendo.org
E568
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
that, in this data set, no “zero” SI problems were present. Thus one expects that SI values greater than a certain threshold are less and less probable. At the same time, because SI can be very small, although we do not know how small it can be in a certain subject before seeing his or her data, it is wise not to incorporate in the prior any beliefs about its variability at low levels. Thus we have assumed for SI an a priori probability density function in which SI values ⬍2 ⫻ 10⫺4 (min⫺1/Uml⫺1) are equally probable, while SI values ⬎2 ⫻ 10⫺4 are less and less probable, according to a decreasing exponential law (with exponent equal to 10⫺4 U ml⫺1/min⫺1). The chosen threshold value of 2 ⫻ 10⫺4 min⫺1/U ml⫺1, even if not directly obtained from a distribution of SI estimates in NIDDM subjects, is realistic and is supported by the literature, where it is unlikely to find insulin sensitivity values in NIDDM subjects ⬎2 ⫻ 10⫺4 min⫺1/U ml⫺1. The MCMC method. The goal is to obtain f兩y(兩y), i.e., the joint probability density function of the unknown minimal model parameters, given the glucose data. This task is analytically intractable, both because of the complex relationships between parameters and data and because of the shape of the probability distributions involved. A solution suggested in the literature is to derive f兩y(兩y) in sampled form by adopting a simulation strategy known as Markov chain Monte Carlo (13). MCMC methods are based on two steps. In the first step, a suitable Markov chain that converges (in distribution) to a target distribution, f兩y(兩y) in our case, is generated. Then, a Monte Carlo integration step is performed to numerically compute the integrals of interest, e.g., Eq. 5. There are many MCMC methods in the literature; see Ref. 13 for a review. They differ from each other in the way the Markov chain is generated. However, all of the different strategies proposed in the literature can be traced back to the Metropolis-Hastings algorithm (15). In this work, we generated the Markov chain by a symmetric transition kernel that extracts a sample of the chain around the previous one (this scheme is usually referred to as “random-walk Metropolis”) by use of independent uniform probability densities. The variance of such densities was chosen case by case. As far as implementation issues of the MCMC method are concerned, we note that convergence of the chain was assessed by the Raftery criterion (21), which is also known in the literature as binary-control. To compute the number of iterations necessary to estimate a posterior quantile from a single run of a Markov chain, a pilot analysis of output samples was used to fit a two-state Markov chain model, and from it one can calculate the length of the burn-in period and, then, the number of further iterations required to estimate quantiles of interest with the required accuracy. In the sequel, we have
been required to estimate quantiles 0.025, 0.25, 0.5, 0.75, 0.975 with precision factors of 0.02, 0.05, 0.05, 0.05, and 0.02, respectively, with a probability of 0.95. From these quantiles, numerically robust confidence intervals can be calculated. It is worth stressing that, from the practical point of view, switching from NLS to Bayes estimation is not without a price. To give an idea, minutes instead of seconds are required to identify the model in each individual. The reason for this increased time is that Markov chain generation and Monte Carlo integration, needed to handle Eqs. 5 and 6, are computationally demanding. In particular, to obtain the results we present in the next section, a number of iterations varying from 35,000 to 80,000, depending on the subject under study, were performed. RESULTS
Data Plasma glucose and insulin concentrations of the 10 subjects are shown in Fig. 2. NLS Identification Results of NLS for SG, p2,, G0, and SI, obtained with SAAM II (2), are reported in Table 1 together with the precision of the estimates. The SI estimates in all of the subjects are shown in Fig. 3, together with their 95% confidence intervals, obtained by adding and subtracting to the point estimate the quantity 1.96 SD (results of subjects 1–7 are those already reported in Table 2 of Ref. 1). As a consequence of the implicit Gaussian assumption, the confidence interval is symmetric and centered around the point estimate. As apparent from the picture, SI estimates are not equally satisfactory in all of the subjects; whereas in subjects 1, 2, 4, 5, 8, and 10 a narrow confidence interval lying on the nonnegative axis is obtained, in subjects 3, 6, 7, and 9, the low value of the SI estimate, combined with its poor precision, produce an unrealistic confidence interval. In fact, in each of these four subjects, the SI confidence interval tells us that the true SI can fall with nonzero probability in the negative portion of the x-axis. Note that, because the covariance matrix of the measurement error is assumed to be known from the data, NLS coincides with maximum likelihood (ML). If it were assumed to be dependent on the model, then NLS and ML are no longer equivalent on a theoretical basis.
Table 2. Bayesian parameter estimates with the refined SI prior Subject Nos.
SI, 10⫺4 min⫺1/U ml⫺1
SG, 10⫺2 min⫺1
p2, min⫺1
Go , mg/dl
1 2 3 4 5 6 7 8 9 10
1.2 (0.61–2.7) 1.5 (1.4–1.6) 1.7 (0.4–4.1) 0.77 (0.58–1.1) 2.1 (1.8–2.2) 0.39 (0.032–2.3) 1.51 (0.36–4.1) 0.98 (0.8–1.11) 0.8 (0.18–2.7) 0.6 (0.46–0.72)
1.59 (1.07–1.98) 2.1 (1.8–2.38) 1.06 (0.78–1.35) 1.02 (0.76–1.25) 0.89 (0.48–1.32) 1.35 (1–1.62) 1.48 (1.32–1.62) 0.92 (0.66–1.19) 1.61 (1–1.85) 0.84 (0.52–1.18)
0.013 (0.0013–0.033) 0.071 (0.057–0.84) 0.0018 (0.0004–0.0068) 0.008 (0.0033–0.013) 0.098 (0.078–0.12) 0.0070 (0.00014–0.07) 0.0019 (0.00033–0.008) 0.35 (0.2–0.68) 0.006 (0.00034–0.024) 0.025 (0.021–0.03)
299.7 (284.1–313.3) 382 (369.9–394.7) 434.6 (421.9–448.4) 318 (307.4–329.2) 411 (395.6–428.4) 563 (543.3–581.3) 314 (304.8–323.5) 325.4 (315.3–336.2) 350.9 (335.4–363.2) 263 (252.5–275.2)
Nos. in parentheses show 95% confidence intervals. AJP-Endocrinol Metab • VOL
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
E569
Fig. 4. The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a nonnegativity prior on the 4 model parameters. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The y-axis quantity is unitless.
However, in this case, the two results are virtually superimposable (unpublished data). Bayesian Identification A nonnegativity prior on the four parameters. The a posteriori probability density function of SI estimated in each subject is reported in Fig. 4. As somewhat expected, in subjects 1, 2, 4, 5, 8, and 10, where NLS was successful, the posterior provided (in sampled form) by MCMC is concentrated around the NLS estimate. However, the new method reveals asymmetries in the distribution of the estimates, particAJP-Endocrinol Metab • VOL
ularly in subjects 1 and 4. This is not surprising, because the symmetry of NLS confidence intervals reflects only the Gaussian assumption implicitly made on the estimates. Moreover, in subjects 3, 6, 7, and 9, MCMC analysis reveals that the marginal posterior densities have a high peak located at low and realistic SI values but that they have a long tail. The confidence intervals, however, would be too pessimistic, because it is known from the literature that such high SI values are unrealistic, e.g., it is impossible that SI is, say, equal to 100 ⫻ 10⫺4 min⫺1/U ml⫺1 (see subject 3 in Fig. 4). So, to obtain a more
282 • MARCH 2002 •
www.ajpendo.org
E570
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 5. The a posteriori probability density function of SI for the 10 subjects (area under each curve is unitary) from Bayesian identification, with a refined SI prior. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The yaxis quantity is unitless.
credible confidence interval, it is necessary to adopt a more realistic SI prior. A refined SI prior. Results for SG, p2, G0, and SI are reported in Table 2, together with the 95% confidence interval. Results for SI are also displayed in Fig. 5, where the estimated a posteriori probability density function of SI in each subject is reported (the axes are now different from those used in Fig. 4), and in Fig. 3, where a summary of Fig. 5 results is depicted to make the comparison with NLS easier. One can note, by AJP-Endocrinol Metab • VOL
comparing Figs. 5 and 4, how the threshold influences the tail of the posterior distribution of some subjects. One can also note that the posterior probability density function is zero in an interval at the right of SI⫽0. This allows us to exclude, on a probabilistic basis, that SI is greater than a chosen threshold. From the a posteriori density function of SI, one can obtain the interval where the true value of SI falls with 95% of confidence as the interval between the quantiles 0.025 and 0.975. This interval is displayed in Fig. 3, together with the
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
posterior mean of SI. The novel features of the approach we propose emerges in subjects 3, 6, 7, and 9, where NLS results were largely unsatisfactory with noncredible confidence intervals. One can see that the marginal a posteriori probability density function of SI does not collapse to zero (this also applies to Fig. 4, but it cannot be visually appreciated due to the larger range of SI values shown in the x-axis). In all of these subjects, the posterior mean of SI is higher than the NLS estimate (markedly in subjects 3 and 7), because the minimum variance estimate takes into account the entire shape of the posterior and, so, also the presence in some subjects of the tails. More important, in each subject the confidence interval now lies entirely in the positive axis and allows the investigator to infer that the probability that the true SI lies at the left of the displayed bar is 2.5%. In summary, a nonzero SI value is detected with high probability by Bayesian identification in these very insulin-resistant individuals. This result is not always possible with NLS, where, given the uncertainty of the SI estimate, one can only conclude that SI is indistinguishable from zero (Fig. 3, e.g., subject 7). A credible confidence interval is obtained, which allows the minimal model user to determine with 97.5% of confidence a (nonzero) lower-bound number for SI. DISCUSSION
In this article, we have attacked the identification of the minimal model of glucose kinetics in type 2 diabetic subjects in a Bayesian context. The motivation for this stems from the interpretative difficulties generated from what is usually referred to as the “zero SI” problem, i.e., in a significant number of cases NLS estimation returns an SI⫽0 when the minimal model program MINMOD (19) is used, or very small values with an unrealistic confidence interval including negative values with software packages like SAAM II (2) or ADAPT (10). In this case study, the Bayesian approach allowed the derivation (albeit only numerically through MCMC) of the joint a posteriori probability distribution of the model parameters, and thus the determination of a nonnegative point estimate of SI, together with a plausible confidence interval. The differences of the point estimates between the Bayesian approach and NLS, at least in most subjects, are not large, but it is when the confidence intervals of the parameter estimates are considered that the Bayesian approach emerges as superior. In fact, because we think that the chosen prior as well as the confidence intervals employed are both reasonable, one can exclude on a probabilistic basis that SI in the critical subjects is higher than a given threshold. Moreover, from the confidence interval obtained by Bayes estimation, one can exclude with 97.5% of probability that the true SI lies at the left of the displayed variability range. In summary, the Bayesian approach is able to cope with those situations in which the confidence interval, obtained as two times the SD of the error given by the Fisher information matrix, is not realistic, because such a method is unAJP-Endocrinol Metab • VOL
E571
able to cope with asymmetries in the distribution of the estimates, like those present in this study. A critical point in Bayes estimation is the definition with a probability density function of the information available on the unknown model parameters before seeing the data. As far as SG, p2, and G0 are concerned, we have limited ourselves to express the poorly informative knowledge that they are nonnegative, i.e., a uniform distribution along the positive axis was chosen. A more sophisticated prior was chosen for SI, the very critical parameter of the model, based on the observation that, for values of SI greater than a certain threshold, the higher the candidate SI, the lower, in all likelihood, its probability. We have translated these beliefs about SI in statistical terms by describing its a priori probability density function as a uniform distribution between zero and a certain threshold (i.e., 2 ⫻ 10⫺4 min⫺1/U ml⫺1, a value that appears as prudential on the basis of published studies on insulin resistance in type 2 diabetics) and an exponential decay for SI greater than this threshold. It is important to note that, in the prior chosen for SI, we do not include beliefs on how small it can be, i.e., where the whole “zero SI” problem began. Choosing also for SI a uniform distribution prior, as in Bayesian identification with only the nonnegativity prior, one would have obtained an SI a posteriori probability density function lower than the one presented. This means that low SI values are less probable and supports the robustness in each subject of the 0.025 quantile of Fig. 1, which allows us to exclude with 97.5% of probability that the true SI lies at the left of the displayed range. Another consequence of using a uniform prior for SI along the positive axis is represented by a gratuitous shift to the right of the 0.975 quantile. In fact, the new technique permits us to show that, in some cases, the model defines likelihood functions which assign a nonnegligible probability to high SI values that are largely unrealistic. This kind of information was impossible to obtain so clearly by the standard Fisherian approaches. The rationale of choosing for SI a prior more sophisticated than that employed for SG, p2, and G0 was to determine with 95% confidence an interval as narrow as possible where SI lies. The role played by this a priori knowledge is not hidden but transparent in terms of its effect on the results. For example, a strategy to obtain nonnegative confidence intervals by NLS could consist in reparameterizing the model with the aim of estimating the logarithm of the parameters. However, this would be obtained by approximating the likelihood as a function of the logarithm of the parameters, with a Gaussian distribution. This is questionable, because it would not be possible to evaluate the quality of such approximation. As already mentioned, switching from NLS to Bayes estimation is not without a price, because the MCMC results are computationally more demanding to obtain. However, this strategy offers a new tool that could also become important in classification studies. In fact, it would be interesting to reanalyze large databases of epidemiological studies in which the standard NLS
282 • MARCH 2002 •
www.ajpendo.org
E572
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
Fig. 6. A: a qualitative distribution of SI estimates that would be obtained in a study in which the SI⫽0 problem is present. On the x-axis, SI is measured in 10⫺4 min⫺1/U ml⫺1. The y-axis shows the no. of SI estimates falling in equally spaced containers. B: qualitative distribution of SI estimates that would be obtained by the Markov chain Monte Carlo (MCMC) method in the same study.
approach has produced SI histograms with a likely artificial bimodal pattern (with one peak at SI⫽0 and another peak at a positive SI value), like that qualitatively displayed in Fig. 6, which are difficult to interpret. The value SI⫽0, which Fisherian techniques return as the joint mode of a suitable objective function, can be a mathematical artifact, because it can be far from the minimum variance estimate provided by MCMC. Therefore, we believe that if the new strategy proposed in this paper were applied in epidemiological studies like those mentioned above, one would obtain a more reliable distribution of SI estimates, e.g., changing the qualitative profile of Fig. 6A into that shown in Fig. 6B. In this way, more reliable information on NIDDM would be made available. Moreover, a reanalysis of databases such as GENNID (20), FUSION (26), and IRAS (16, 17), where the numbers of subjects are in the range of thousands, could also allow a refinement of our MCMC method, because a more accurate description of the a priori probability density function of SI will be possible. It is clear, in fact, that the refined prior used in Fig. 5 provides better results than those obtained with the positivity prior and shown in Fig. 4. However, the refined prior still has a certain degree of roughness. A future evolution of the method could also envisage a population approach with the aim of estimating, jointly with model parameters, the common prior that has generated them. Finally, it is worth stressing that important questions concerning possible problems introduced by the structure of the considered model still remain open also in the light of our MCMC results. In fact, one can clearly note that the a posteriori probability density functions of SI of subjects 3 and 7 are quite close to the AJP-Endocrinol Metab • VOL
prior. As a consequence, the point estimates are close to the mean of the prior, equal to 1.67 10⫺4 min⫺1/U ml⫺1, suggesting that sometimes the model extracts little information on SI from IVGTT data (in fact the posterior distribution approaches the prior as the data contain decreasing information about the parameters). Moreover, the presence of long tails in the a posteriori probability density of SI of some subjects is a critical issue that will have to be clarified in future studies. Our Bayesian analysis, including reasonable a priori information on SI, improves the estimation process, reducing this kind of problem. However, it is clear that no single parameter estimation technique, however sophisticated, can be a panacea, able to compensate possible structural defects of a model. So, further investigations are necessary to clarify whether the cause of these problems is the fact that the model is structurally inadequate for the estimation of a low SI. We are inclined to believe that, in many of such critical situations, MCMC is really a more robust estimator than NLS or ML and can significantly improve the estimation process. However, we also believe that in the future it will also be worth exploring a different description of the same process, like the two-compartment model described in Refs. 7 and 8, which may overcome these problems. Maybe, after including some suitable a priori information, which is necessary because of a priori identifiability issues, the problems mentioned here may be solved by resorting to a maximum a posteriori estimator. In any case, it will be interesting to identify the two-compartment model by resorting to the same MCMC approach described here. In fact, even if more computationally demanding, this technique is likely to provide a clearer picture of the
282 • MARCH 2002 •
www.ajpendo.org
MINIMAL MODEL “ZERO” SI SOLVED BY BAYESIAN APPROACH
advantages of the two- vs. the one-compartment description. This work was in part supported by National Institutes of Health Grant RR-12609, by the MURST COFIN 2000 project, “Estimation of nonaccessibile parameters in physiological systems,” by a University of Padova grant (STIM-PET) to G. Sparacino, and by a University of Pavia grant to P. Magni.
15. 16.
17.
REFERENCES 1. Avogaro A, Vicini P, Valerio A, Caumo A, and Cobelli C. The hot but not the cold minimal model allows precise assessment of insulin sensitivity in NIDDM subjects. Am J Physiol Endocrinol Metab 270: E532–E540, 1996. 2. Barrett PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P, and Foster D. SAAM II: simulation, analysis and modeling software for tracer and pharmacokinetic studies. Metabolism 47: 484–492, 1998. 3. Bergman RN, Bortolan G, Cobelli C, and Toffolo G. Identification of a minimal model of glucose disappearance for estimating insulin sensitivity. In: Proc IFAC Symposium Identification and System Parameter Estimation. Darmstadt, Germany: IFAC, 1979, vol. 2, p. 883–890. 4. Bergman RN, Finegood DT, and Ader M. Assessment of insulin sensitivity in vivo. Endocrinol Rev 6: 45–46, 1985. 5. Bergman RN, Ider YZ, Bowden CR, and Cobelli C. Quantitative estimation of insulin sensitivity. Am J Physiol Endocrinol Metab Gastrointest Physiol 236: E667–E677, 1979. 6. Carson ER, Cobelli C, and Finkelstein L. The Mathematical Modelling of Metabolic and Endocrine Systems. New York: Wiley, 1983. 7. Caumo A, Vicini P, Zachwieja JJ, Avogaro A, Yarasheski K, Bier DM, and Cobelli C. Undermodeling affects minimal model indexes: insights from a two-compartment model. Am J Physiol Endocrinol Metab 276: E1171–E1193, 1999. 8. Cobelli C, Caumo A, and Omenetto M. Minimal model SG overestimation and SI underestimation: improved accuracy by a Bayesian two-compartment model. Am J Physiol Endocrinol Metab 277: E481–E488, 1999. 9. Cobelli C, Foster D, and Toffolo G. Tracer Kinetic in Biomedical Research: from Data to Model. New York: Kluwer Academic/ Plenum, 2000. 10. D’Argenio D and Schumitzky A. ADAPT II Users Guide: Pharmacokinetic/Pharmacodynamic Systems Analysis Software. Los Angeles, CA: Biomedical Simulations Resource, University of Southern California, 1997. 11. DeFronzo RA, Bonadonna RC, and Ferranini E. Pathogenesis of NIDDM: a balanced overview. Diabetes Care 15: 318–368, 1992. 12. DeFronzo RA, Tobin JD, and Andres R. Glucose clamp technique: a method for quantifying insulin secretion and resistance. Am J Physiol Endocrinol Metab Gastrointest Physiol 237: E214–E223, 1979. 13. Gilks WR, Richardson S, and Spiegelhalter DJ. Markov Chain Monte Carlo in Practice. London: Chapman & Hall, 1996. 14. Haffner SM, D’Agostino R Jr, Mykkanen L, Tracy R, Howard B, Rewers M, Selby J, Savage PJ, and Saad MF. Insulin sensitivity in subjects with type 2 diabetes. Relationship
AJP-Endocrinol Metab • VOL
18. 19.
20.
21.
22. 23.
24.
25. 26.
27. 28.
E573
to cardiovascular risk factors: the Insulin Resistance Atherosclerosis Study. Diabetes Care 22: 562–568, 1999. Hastings WK. Monte Carlo sampling methods using Markov chain and their applications. Biometrika 57: 97–109, 1970. Howard G, Bergman R, Wagenknecht LE, Haffner SM, Savage PJ, Saad MF, Laws A, and D’Agostino RB Jr. Ability of alternative indices of insulin sensitivity to predict cardiovascular risk: comparison with the ’minimal model.’ Ann Epidemiol 8: 358–369, 1998. Mykkanen L, Zaccaro DJ, Hales CN, Festa A, and Haffner SM. The relation of proinsulin and insulin to insulin sensitivity and acute insulin response in subjects with newly diagnosed type II diabetes: the Insulin Resistance Atherosclerosis Study. Diabetologia 42: 1060–1066, 1999. Owens DR, Luzio SD, and Coates PA. Insulin secretion and sensitivity in newly diagnosed NIDDM Caucasians in the UK. Diabetes Med 13: S19–S24, 1996. Pacini G and Bergman RN. MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsively from the frequently sampled intravenous glucose tolerance test. Comput Methods Prog Biomed 23: 113–122, 1986. Raffel LJ, Robbins DC, Norris JM, Boerwinkle E, DeFronzo RA, Elbein RC, Fujimoto W, Hanis CL, Kahn SE, Permutt MA, Chiu KC, Cruz J, Ehrmann DA, Robertson RP, Rotter JI, and Buse J. The GENNID study. A resource for mapping the genes that cause NIDDM. Diabetes Care 19: 864– 872, 1996. Raftery AE and Lewis SM. Implementing MCMC. In: Markov Chain Monte Carlo in Practice, edited by WR Gilks, S Richardson, and DJ Spiegelhalter. London: Chapman & Hall, 1996, p. 115–130. Reaven GM. Banting lecture 1988. Role of insulin resistance in human disease. Diabetes 37: 1595–1607, 1988. Saad MF, Anderson RL, Laws A, Watanabe RM, Kades WW, Chen YDI, Sands RE, Pei D, Savage PJ, and Bergman RN. A comparison between the minimal model and the glucose clamp in the assessment of insulin sensitivity across the spectrum of glucose tolerance. Diabetes 43: 1114–1121, 1994. Sparacino G, Tombolato C, and Cobelli C. Maximum-likelihood vs. maximum a posteriori parameter estimation of physiological system models: the C-peptide impulse response case study. IEEE Trans Biomed Eng 47: 801–811, 2000. Tierney L. Markov chains for posterior distributions (with discussion). Ann Statistics 22: 1701–1762, 1994. Valle T, Tuomilehto J, Bergman RN, Gosh S, Hauser ER, Eriksson J, Nylund SJ, Kohtamaki K, Toivanen L, Vidfren G, Tuomilehto-Wolf E, Ehnolm C, Blaschal J, Langefeld CD, Watanabe RM, Magnuson V, Ally DS, Hagopian WA, Ross E, Buchanan TA, Collins F, and Boehnke M. Mapping genes for NIDDM design of the Finland-United States investigation of NIDDM genetics (FUSION) study. Diabetes Care 21: 949–956, 1998. Walter E and Pronzato L. Identification of Parametric Models from Experimental Data. Paris: Springer, 1997. Welch S, Gebhart SS, Bergman RN, and Phillips LS. Minimal model analysis of intravenous glucose tolerance test-derived insulin sensitivity in diabetic subjects. J Clin Endocrinol Metab 71: 1508–1518, 1990.
282 • MARCH 2002 •
www.ajpendo.org
