Estimation of doubly labeled water energy expenditure with confidence intervals

Richard H. Jones, Bakary J. Sonko, Leland V. Miller, Patti J. Thureen, and Paul V. Fennessey

Departments of Preventive Medicine and Biometrics and of Pediatrics, School of Medicine, University of Colorado Health Sciences Center, Denver, Colorado 80262


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

Bivariate regression is used to estimate energy expenditure from doubly labeled water data. Two straight lines are fitted to the logarithms of the enrichments of oxygen-18 and deuterium simultaneously as a bivariate regression, so that the correlations between the oxygen and deuterium regression coefficients can be estimated. Maximum likelihood methods are used to extend bivariate regression to unbalanced situations caused by missing observations and to include replicate laboratory determination from the same urine samples, even if one of the replicates is missing. Use of maximum likelihood allows the determination of a confidence interval for the energy expenditure based on the log likelihood surface rather than use of the propagation of variance methods for nonlinear transformations. The model is extended to include the subject's deviations from the two lines as a bivariate continuous-time first-order autoregression to allow for serial correlation in the observations. The analysis of data from two subjects, one without apparent serial correlation and one with serial correlation, is presented.

deuterium; oxygen-18; maximum likelihood


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

THE DOUBLY LABELED WATER TECHNIQUE, pioneered by Lifson (7, 8) for energy expenditure measurement in animals, is now frequently used for precise estimation of free living energy expenditure in humans (2, 13). The technique uses two stable isotopes, oxygen-18 (18O) and deuterium (2H), which can exist in water form as H218O and 2H2O, respectively. The mixed isotopes are given to the subject to drink, and after the isotopes are allowed to equilibrate with total body water, samples of body fluid, e.g., urine or saliva, are taken to determine isotopic enrichment over the measurement period, 14-21 days for adults and 5-8 days for infants (10). The elimination of the isotopes from the body describes a monoexponential function in which the 2H is eliminated only as water, and the 18O is eliminated both as water and carbon dioxide (CO2) (10). The difference in the elimination rates of the two isotopes gives an estimate of the rate of CO2 production in the body. The CO2 production can be converted to energy expenditure by using either Weir's equation (14) or Elia's equation (5), if the respiratory quotient (RQ) of the subject over the study period is known. In an "ideal" experiment, the isotope elimination curve would consist of as many points (samples) as could be conveniently collected during the course of the experiment. Under more practical conditions, the number of points in each curve is limited by both the sample collection and sample processing constraints, as well as the cost of 18O. Current practice, as determined from the various methods sections of publications, records as few as two (two-point method) data points to as many as 14-22 (multipoint method) data points (10) used to determine the rate of the isotope elimination.

Cole and Coward (1) have reported a statistics technique to determine confidence intervals for these approaches. Their model, in their notation, follows the recommendation of Prentice (10)
log E<SUB>D<SUB><IT>i</IT></SUB></SUB> = log <IT>I</IT><SUB>D</SUB> − <IT>k</IT><SUB>D</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + error term (1)

log E<SUB>O<SUB><IT>i</IT></SUB></SUB> = log <IT>I</IT><SUB>O</SUB> − <IT>k</IT><SUB>O</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + error term (2)
for i = 1...n, where the subscript D denotes deuterium, and O denotes oxygen. Ei is the isotope enrichment after subtraction of background expressed as a fraction of dose per mole of body water, ti is the time (in units of days), I is the enrichment at time 0, and k is the flux rate. The error terms are assumed to have zero means and constant variances. The multipoint method estimates log I and k in each equation by linear regression, which also produces estimated standard errors of the estimates. This model is appropriate when the original (untransformed) data have a constant coefficient of variation, i.e., the standard deviation of the observations, EDi and EOi are proportional to the means, ID exp(-kDti) and IO exp(-kOti). In this case, the log transformation is a variance-stabilizing transformation. This assumption can be checked by plotting the residuals after fitting the lines to the transformed data (the differences between the observed log enrichment and the predicted value based on the estimated regression coefficients). These residuals should show a fairly random pattern. It is known that the error terms between the two models are correlated; therefore, a bivariate regression is appropriate where the observations are treated as pairs.

The CO2 production rate, with fractionation ignored, is
r′<SUB>CO<SUB>2</SUB></SUB> = (<IT>k</IT><SUB>O</SUB><IT>N</IT><SUB>O</SUB> − <IT>k</IT><SUB>D</SUB><IT>N</IT><SUB>D</SUB>)/2 (3)
where NO and ND are the body water pool sizes obtained as
<IT>N</IT><SUB>O</SUB> = 1/<IT>I</IT><SUB>O</SUB> <IT>N</IT><SUB>D</SUB> = 1/<IT>I</IT><SUB>D</SUB>
The main statistical problem in obtaining the estimate of the CO2 production rate in Eq. 3 is that this is a highly nonlinear function of the estimated parameters in the two linear regressions, and the covariance between the estimates of kONO and kDND is unknown. Cole and Coward (1) argue that this covariance problem is much reduced by working with two different linear regression equations, the difference between Eqs. 1 and 2, and the sum of Eqs. 1 and 2
log(E<SUB>O<SUB><IT>i</IT></SUB></SUB>/E<SUB>D<SUB><IT>i</IT></SUB></SUB>) = log <IT>I</IT><SUB>r</SUB> − <IT>k</IT><SUB>r</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + error term

log(E<SUB>O<SUB><IT>i</IT></SUB></SUB>E<SUB>D<SUB><IT>i</IT></SUB></SUB>) = log <IT>I</IT><SUB>p</SUB> − <IT>k</IT><SUB>p</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + error term
where Ir = IO/ID (r denotes ratio), Ip = IOID (p denotes product), kr = kO - kD, and kp = kO + kD. The error terms in these transformed equations will be uncorrelated if the variances in Eqs. 1 and 2 are equal, even if the errors in the untransformed equations are correlated. The CO2 production rate is now expressed as a function of the four parameters in these last two linear regressions. The approximate variance of the CO2 production rate is calculated using the delta method, where the expression for the CO2 production rate is expanded in a Taylor series (keeping only the first-order terms) about the estimated values of the four parameters. The expression for the approximate variance of the CO2 production rate is given as Eq. 13 in Cole and Coward (1).

In this article, we use maximum likelihood estimation based on the multipoint method. Confidence limits for energy expenditure are obtained from the log likelihood surface so that methods of propagation of error variances through nonlinear functions are not necessary. Also, instead of treating the logs of the isotope enrichments as univariate regressions, the pair of observations 2H and 18O are treated as a bivariate regression, with correlation between the two observations in the pair. In the analysis of doubly labeled water data, the problem can be somewhat more complicated if one of the measurements at a given time is missing. The other observation can be discarded, but this is a waste of information. If each urine sample is split and analyzed twice, it is possible to estimate the measurement error (analytical variance of the laboratory) as well as the physiological variance of the subject. Prentice (Appendix 4.2, p. 290-293 of Ref. 10) discusses the effect of autocorrelated errors on linear regression and states that "the effect of auto-correlation is to reduce the advantage of the multi-point method." It is also stated that

Because autocorrelation between errors may arise naturally or may arise from attempting to fit an incorrect model, we recommend that if there are ten or more equally-spaced data points, the autocorrelation coefficient for the residuals is obtained. A high absolute value may help to identify wrongly specified models.

Autocorrelation is also referred to as serial correlation. Serial correlation can be observed in residual plots when the residuals seem to have systematic excursions rather than a random appearance. Serial correlation does not affect the position of the fitted line very much, but it does affect confidence intervals of estimates in a nonconservative way. If serial correlation exists and is ignored in the model, the confidence intervals will be too narrow. In this paper we take the approach that serial correlation in the errors can be modeled, even for unequally spaced data, and tested using likelihood ratio tests to determine whether the serial correlation is statistically significant (6). If serial correlation is present in the data, the error model gives more appropriate statistical tests and confidence intervals.


    METHODS
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

The methods presented here are based on bivariate regression (a special case of multivariate regression) in which there are two correlated response variables. Classically, these methods require balanced data, so there can be no missing observations, and, if there are replicate analytical determinations at each time, there must be the same number at every time, and the means of these replications are used as data. Models of this type can be fitted using SAS (11) PROC GLM by using the MANOVA statement. However, replicates must be averaged, and, if an observation is missing at some time point, that time point is discarded. This is due to the assumption of equal covariance matrices at different times in the bivariate regression. Because there are two components of variance, the analytical variance of the laboratory and the physiological variation of the subject, it is not possible to weight the regression by the number of replications. A proper weighting would be a function of these two components of variance.

Physiological variation causes both log enrichments to vary in a similar way about the fitted lines. This is where serial correlation may come into the model. This variation about the lines may not be independent from observation time to observation time, as is required by standard regression analysis. If an enrichment is above the line at one time point, it is likely to be above the line at the next time point. This variation about the line in a dependent fashion is serial correlation, and the physiological variation is a random process. A physiological model that fits the doubly labeled water data well is
<IT>y</IT><SUB>D</SUB>(<IT>t</IT>) = &agr;<SUB>D</SUB> + &bgr;<SUB>D</SUB><IT>t</IT> + g<SUB>D</SUB>&ggr;(<IT>t</IT>)

<IT>y</IT><SUB>O</SUB>(<IT>t</IT>) = &agr;<SUB>O</SUB> + &bgr;<SUB>O</SUB><IT>t</IT> + g<SUB>O</SUB>&ggr;(<IT>t</IT>)
where yD(t) = logED(t), yO(t) = logEO(t), alpha D = logID, alpha O = log IO, beta D = -kD, beta O = -kO, and gamma (t) represents the physiological variation about the lines as a function of time and is assumed to have zero mean and unit variance. This random physiological variation enters both equations with different weights, gD and gO, that are estimated from the subject's data; gD and gO are the physiological standard deviations of the two equations. This single random input entering both equations produces correlation between the two equations. The model for a single observational pair is
<IT>y</IT><SUB>D<SUB><IT>j</IT></SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>) = &agr;<SUB>D</SUB> + &bgr;<SUB>D</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + g<SUB>D</SUB>&ggr;(<IT>t</IT><SUB><IT>i</IT></SUB>) + &egr;<SUB>D<SUB><IT>j</IT></SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)

<IT>y</IT><SUB>O<SUB><IT>j</IT></SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>) = &agr;<SUB>O</SUB> + &bgr;<SUB>O</SUB><IT>t</IT><SUB><IT>i</IT></SUB> + g<SUB>O</SUB>&ggr;(<IT>t</IT><SUB><IT>i</IT></SUB>) + &egr;<SUB>O<SUB><IT>j</IT></SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>) (4)
In this model, ti denotes the time of the observations, the subscript j denotes replications, and epsilon Dj(ti) and epsilon Oj(ti) are the laboratory errors for replications and are assumed to have means of zero, two different variances, sigma 2D and sigma 2O, and to be uncorrelated. The errors are also assumed to be independent at different observation times and independent of the physiological variation. In the Cole and Coward (1) model, the error term in Eq. 1 would be gDgamma (ti) + epsilon D(ti) if there are no replicates, and has variance g2D + sigma 2D. If two replicates are averaged, the error term would have variance g2D + sigma 2D/2. This produces different error variances if some of the replicates are missing.

Estimating the variance of the laboratory precision. The variance of the laboratory precision can be estimated directly from the replication data. If there are mi replications at observation time i, these replications can be averaged and the sum of squares of the deviations from the estimated mean can be calculated. This estimate has mi - 1 degrees of freedom. The sums of squares and degrees of freedom are summed over observation times, and the variance is estimated as the final sum of squares divided by the total number of degrees of freedom. This must be carried out separately for D and O, because the laboratory precision may be different for deuterium and oxygen-18 determination. Because these two variances can be estimated first, it eliminates two parameters from the maximum likelihood estimation. The standard deviations of the laboratory precision are estimated as the square roots of these variances.

Modeling the physiological variation. The simplest model for physiological variation is that gamma (ti) is statistically independent at different observation times and has a Gaussian distribution with zero mean and unit variance. The variances of the two physiological inputs are then g2D and g2O, and the correlation is 1. This assumption means that the random input due to physiological variation is uncorrelated from time point to time point, but at each time, the random deviation from the deuterium line is proportional to the random deviation from the oxygen line. Because of the laboratory error, the actual observations are not perfectly correlated.

This model has four parameters in the linear regression, two standard deviation parameters for the laboratory precision, sigma D and sigma O, and two nonlinear parameters, gD and gO, that are the standard deviations of the physiological variation. Details of the maximum likelihood estimation, the method of determining confidence intervals on the energy expenditure, and the extension to models with physiological serial correlation are shown in the APPENDIX.


    EXAMPLES
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

Table 1 shows the enrichment data for a subject without significant serial correlation (subject 1). Missing values are indicated by periods. The details for the calculation of the scale factors are reported by Prentice (p. 217, Ref. 10). The enrichments are represented as a fraction of the initial dose giving a normalized enrichment by use of the formula
<FR><NU>&dgr;<SUB>s</SUB> − &dgr;<SUB>p</SUB></NU><DE>&dgr;<SUB>a</SUB> − &dgr;<SUB>t</SUB></DE></FR> × <FR><NU>18.02a</NU><DE>WA</DE></FR>
where delta s is the enrichment of the sample, delta p is the predose enrichment of the subject, delta a is the enrichment of the dose, delta t is the enrichment in the tap water, a is the amount of the dose diluted for analysis in grams, W is the amount of water used to dilute the dose (g), A is the amount of dose administered (g), and 18.02 converts grams of water into moles. To begin, the average of the two baseline O values is subtracted from each O in the raw data, and the average of the two baseline D values is subtracted from each D. This is the delta s - delta p part of the conversion. The other parameters for subject 1 are delta a(O) = 219.51, delta a(D) = 847.97, delta t(O) -12.20, delta t(D) = -77.416, a = 0.0121, W = 2.3845, and A = 173.5. The RQ for this subject, needed to convert CO2 production to energy expenditure, is estimated from food intake data to be RQ = 0.86. 

                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Data for subject 1 

The natural logarithms (base e) of the scaled data are calculated, and the standard deviations of the laboratory precision are estimated from the replications of the log data. For this data set, the standard deviations of the laboratory precisions are estimated as delta O = 0.00073 and delta D = 0.01166. It should be noted that the standard deviation of the laboratory precision based on the logs of the data is the coefficient of variation (standard deviation divided by the mean) of the original data.

After the lines are fitted to the data, the predicted values of the lines at each time point are calculated as well as the residuals. The normalized residuals are calculated by dividing each residual by its standard deviation. Plots of the fitted lines and the normalized residuals are shown in Fig. 1 (left). A normalized residual that is >2.5 in absolute value should be investigated as a possible error. Various investigators use values of normalized residuals in the range of 2.0-3.0 for detecting errors. The second D normalized residual at 7.06 days is -2.93, whose absolute value is >2.5. Going back to the raw data, we see that the value 521.742 appears to be too low relative to the other data. Setting this value to missing produces normalized residuals, where the largest in absolute value is the second D normalized residual of -2.06 at 12.06 days. This is considered to be acceptable. The estimated coefficient of variation of the measurement precision is decreased slightly to <A><AC>&sfgr;</AC><AC>ˆ</AC></A>D = 0.00969. 


View larger version (24K):
[in this window]
[in a new window]
 
Fig. 1.   Results for 2 subjects. Subject 1 (top left) has no apparent serial correlation. Plot shows log data points and fitted straight lines for subject 1. Plot (bottom left) shows normalized residuals for subject 1. It was determined that the one normalized residual at 7.06 days with a value of -3.03 was probably an error, so this data point was declared missing in the final analysis. Results for subject 2, with significant serial correlation, are on right. Normalized residuals show systematic trends. Recursive residuals use serial correlation to predict a residual from the previous residuals and should show a more random pattern. The low set of recursive residuals at 5 days is possibly a timing error, but these data were included in the analysis.

From the regression analysis, 95% confidence limits on the slopes and intercepts can be directly calculated and inverted to get confidence limits on NO and ND
<IT>k</IT><SUB>O</SUB> = 0.0820 (0.0807, 0.0832)/day

<IT>k</IT><SUB>D</SUB> = 0.0574 (0.0562, 0.0586)/day

<IT>N</IT><SUB>O</SUB> = 1879 (1861, 1898) moles

<IT>N</IT><SUB>D</SUB> = 1941 (1922, 1960) moles
The ND-to-NO ratio is estimated as ND/NO = 1.033 (1.023, 1.042). This is in the range of physiological values. The uncorrected CO2 production is estimated as (NOkO - NDkD)/2 = 21.33 (20.81, 21.90) mol/day. The confidence intervals on CO2 and energy expenditure are obtained by using the log likelihood surface, as explained in APPENDIX. A fractionation correction is then applied to correct for fractionation loss of both isotopes in breath and transcutaneous water, and 18O in CO2. Let the fractionation factors for 2H and 18O in breath and transcutaneous water be f1 and f2, respectively, and the fractionation factor for 18O in CO2 be f3. These have values: f1 = 0.941, f2 = 0.992, and f3 = 1.040. The corrected CO2 production is calculated (see p. 101 in Ref. 10) as
<FR><NU><IT>N</IT><SUB>O</SUB><IT>k</IT><SUB>O</SUB> − <IT>N</IT><SUB>D</SUB><IT>k</IT><SUB>D</SUB> − 27(<IT>f</IT><SUB>2</SUB> − <IT>f</IT><SUB>1</SUB>)</NU><DE>2<IT>f</IT><SUB>3</SUB> + 1.1(<IT>f</IT><SUB>2</SUB> − <IT>f</IT><SUB>1</SUB>)</DE></FR> = <FR><NU><IT>N</IT><SUB>O</SUB><IT>k</IT><SUB>O</SUB> − <IT>N</IT><SUB>D</SUB><IT>k</IT><SUB>D</SUB> − 1.377</NU><DE>2.1361</DE></FR>
giving a corrected CO2 estimate of 19.32 (18.84, 19.86) mol/day. The estimated CO2 production is multiplied by 22.4 l/mol to convert to liters/day of CO2. To convert from liters of CO2 per day to total energy expenditure (kilojoules per day), liters per day is multiplied by (5.55 + 15.48/RQ) (see p. 198 of Ref. 10), where RQ is the subject's respiratory quotient. The estimated total energy expenditure is 10,194 (9,938, 10,478) kJ/day. This subject has a nonsignificant test for serial correlation by use of a likelihood ratio test (see APPENDIX), and the normalized residuals shown at the bottom left of Fig. 1 do not appear to have any systematic trends.

The data for a second subject, who had significant serial correlation, are shown in Table 2. The parameters for subject 2 are delta a(O) = 173.88, delta a(D) = 649.969, delta t(O) = -16.75, delta t(D) = -117.745, a = 0.0121, W = 4.9257, and A = 211.52. The RQ for this subject is estimated to be 0.857. 

                              
View this table:
[in this window]
[in a new window]
 
Table 2.   Data for subject 2 

The results of the regression for subject 2 are
<IT>k</IT><SUB>O</SUB> = 0.0930 (0.0875, 0.0986)/day

<IT>k</IT><SUB>D</SUB> = 0.0707 (0.0658, 0.0756)/day

<IT>N</IT><SUB>O</SUB> = 2963 (2835, 3097) moles

<IT>N</IT><SUB>D</SUB> = 3027 (2909, 3150) moles
The ND-to-NO ratio is estimated as ND/NO = 1.022 (1.010, 1.033). This is also in the range of physiological values. The uncorrected CO2 production is estimated as (NOkO - NDkD)/2 = 30.78 (29.79, 31.85) mol/day. The corrected CO2 production is calculated as 28.17 (27.25, 29.17) mol/day. The estimate of energy expenditure is 14,902 (14,413, 15,431) kJ/day.

The normalized residuals for this subject are shown in Fig. 1 (upper right). These residuals show a more systematic behavior than those from the first subject. This second subject has significant serial correlation based on a likelihood ratio test. The change in -2 ln likelihood when the serial correlation parameter is included is 4.93 and is tested as chi 2 with one degree of freedom (P = 0.014). When serial correlation is modeled, the physiological deviations from the lines are related at different times. The estimated deviation at one time is used to predict the deviation at the next time. The difference between the predicted deviation and observed deviation can be called a recursive residual. These recursive residuals should be nearly uncorrelated if the error structure is properly modeled. Figure 1 (bottom right) shows the normalized recursive residuals for the second subject. Modeling the serial correlation appears to have removed the systematic behavior. The residual at 5 days looks suspicious, because the data points are all slightly below -2.0. This could be a timing error, possibly caused by an improper urine collection. These data points were not removed from the analysis.

Table 3 shows a comparison of the estimated standard errors of the estimated energy expenditure obtained from the -2 ln likelihood surface and by the method of Cole and Coward (1) for 10 subjects. The standard error of the likelihood method is obtained by dividing the 95% confidence interval by 4 (±2 SE). A computer program in FORTRAN is available from the first author. The data input part of the program is specific for our data format from the spreadsheet that we use. This part of the program would need to be modified by a user with a different data format. The program is fairly large because it includes nonlinear optimization subroutines used to search for maximum likelihood estimates of the nonlinear parameters.

                              
View this table:
[in this window]
[in a new window]
 
Table 3.   Estimates of energy expenditure with estimated standard errors for 10 subjects by the likelihood method and the method of Cole and Coward


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

We have presented a new multipoint method for estimating energy expenditure using doubly labeled water. The method is based on bivariate regression with two components of variation, physiological variation and the variances of the laboratory precision. Unbalanced data caused by missing observations are handled using a maximum likelihood approach. Confidence intervals for the estimated energy expenditure are obtained directly from the log likelihood surface.

Table 3 shows that the estimated energy expenditures are very close to the estimates obtained by the Cole and Coward (1) method. The standard errors for the Cole and Coward method tend to be somewhat larger than those of the likelihood method. This can be partly explained by the fact that the likelihood method uses all available data when a replicate or a replicate pair are missing. It is possible that the approximations used when propagating variances through nonlinear functions may tend to increase the estimated standard errors.

The model is extended to handle serial correlation in the physiological variation. This model is an improvement over the model in the book by Jones (see p. 181 in Ref. 6). The new error model based on a continuous-time first-order autoregression fits the data better.

The variances of the logs of the laboratory precision are estimated from the replicate data, where the urine samples are divided into two samples and analyzed separately. Each replicate pair produces one degree of freedom for the variance estimate. If one of the replicates is missing, there is no information about the variance from that pair. Two variances are estimated, one for deuterium and one for oxygen-18. These unbiased estimates of the two variances are then held fixed in the maximum likelihood estimates of the coefficients of the two lines and the parameters of the physiological variation.

Normalized residuals from both the estimation of the variances of the laboratory precisions and the fitting of the lines are used for error detection. Based on the standard normal distribution, there is some question about how to choose a cutoff value for detecting an outlier. In our data we usually have ~14 observations on both deuterium and oxygen-18. If there is an outlier in the data, it tends to pull the fitted line toward itself, making the residual smaller. We recommend that an outlier with an absolute value >2.5 be considered as a possible error. A case can be made for a cutoff of anything in the range of 2.5-3, because the number of observations increases the probability of an observation in the range of 2-3 in absolute value. In engineering, the value of 3 is often used, but there may be thousands of data points.

If there is serial correlation in the physiological variation, it is important that it be modeled in the error structure. Modeling the serial correlation error structure affects the estimation of the lines very little, but it does affect the width of the confidence intervals. In general, confidence intervals are shorter if the data have serial correlation and it is not modeled, and these shorter intervals are not correct. It is possible that our model can underestimate the width of the confidence intervals. Some of the parameters, such as RQ, used in the model are assumed to be known, when they are actually estimates.


    APPENDIX
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

Maximum Likelihood Estimation

To handle laboratory replications, physiological variation, and missing observations, it is necessary to use maximum likelihood estimation. These estimates are usually obtained by finding the maximum of the log likelihood or the minimum of -2 log likelihood with respect to the unknown parameters. The reason -2 log likelihood is often used is that chi 2 likelihood ratio tests for significance of parameters can be calculated from the change in -2 log likelihood when parameters are added or eliminated from the model. In textbook problems, maximum likelihood estimates are usually obtained by differentiating the log likelihood with respect to the unknown parameters, setting the derivatives to zero, and solving a system of equations. In the real world there is often no closed-form solution, and numerical searches need to be used.

In matrix notation, the regression model at a single time point, ti, with a replicate for both the 2H and 18O, is
<FENCE><AR><R><C><IT>y</IT><SUB>D<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>D<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>O<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>O<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R></AR></FENCE> = <FENCE><AR><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>0</C><C>0</C><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C></R><R><C>0</C><C>0</C><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C></R></AR></FENCE> <FENCE><AR><R><C>&agr;<SUB>D</SUB></C></R><R><C>&bgr;<SUB>D</SUB></C></R><R><C>&agr;<SUB>O</SUB></C></R><R><C>&bgr;<SUB>O</SUB></C></R></AR></FENCE> + <FENCE><AR><R><C>g<SUB>D</SUB></C></R><R><C>g<SUB>D</SUB></C></R><R><C>g<SUB>O</SUB></C></R><R><C>g<SUB>O</SUB></C></R></AR></FENCE> &ggr;(<IT>t</IT><SUB><IT>i</IT></SUB>) + <FENCE><AR><R><C>&egr;<SUB>D<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C>&egr;<SUB>D<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C>&egr;<SUB>O<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C>&egr;<SUB>O<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R></AR></FENCE> (5)
If the replication error variances are sigma 2D and sigma 2O, the total covariance matrix of the errors in this model is
<B>V</B><SUB>i</SUB> = <FENCE><AR><R><C>g<SUP>2</SUP><SUB>D</SUB> + &sfgr;<SUP>2</SUP><SUB>D</SUB></C><C>g<SUP>2</SUP><SUB>D</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C></R><R><C>g<SUP>2</SUP><SUB>D</SUB></C><C>g<SUP>2</SUP><SUB>D</SUB> + &sfgr;<SUP>2</SUP><SUB>D</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C></R><R><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUP>2</SUP><SUB>O</SUB> + &sfgr;<SUP>2</SUP><SUB>O</SUB></C><C>g<SUP>2</SUP><SUB>O</SUB></C></R><R><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUP>2</SUP><SUB>O</SUB></C><C>g<SUP>2</SUP><SUB>O</SUB> + &sfgr;<SUP>2</SUP><SUB>O</SUB></C></R></AR></FENCE> (6)
If we let
<B>y</B><SUB><B>i</B></SUB> = <FENCE><AR><R><C><IT>y</IT><SUB>D<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>D<SUB>2</SUB></SUB>(<IT>t</IT><SUB>i</SUB>)</C></R><R><C><IT>y</IT><SUB>O<SUB>1</SUB></SUB>(<IT>t</IT><SUB>i</SUB>)</C></R><R><C><IT>y</IT><SUB>O<SUB>2</SUB></SUB>(<IT>t</IT><SUB>i</SUB>)</C></R></AR></FENCE> <B>X</B><SUB><B>i</B></SUB> = <FENCE><AR><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>0</C><C>0</C><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C></R><R><C>0</C><C>0</C><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C></R></AR></FENCE> <B>&bgr; = </B><FENCE><AR><R><C><B>&agr;<SUB>D</SUB></B></C></R><R><C><B>&bgr;<SUB>D</SUB></B></C></R><R><C><B>&agr;<SUB>O</SUB></B></C></R><R><C><B>&bgr;<SUB>O</SUB></B></C></R></AR></FENCE>
and if we assume that the errors are Gaussian, the contribution to -2 log likelihood from this time point (dropping the constant terms involving pi ), is (6)
ℓ<SUB><IT>i</IT></SUB> = log ‖<B>V<SUB>i</SUB></B>‖ + (<B>y</B><SUB><B>i</B></SUB> − <B>X</B><SUB><B>i</B></SUB>&bgr;)′<B>V</B><SUP>−1</SUP><SUB><B>i</B></SUB>(<B>y</B><SUB><B>i</B></SUB> − <B>x</B><SUB><B>i</B></SUB>&bgr;)
where |Vi| denotes the determinant, and ' denotes the transposed vector. Assuming that the physiological variation does not have serial correlation, i.e., it is statistically independent at different observation times, the -2 log likelihood from all observation times is
ℓ = <LIM><OP>∑</OP><LL>i</LL></LIM> [log ‖<B>V</B><SUB><B>i</B></SUB>‖ + (<B>y</B><SUB><B>i</B></SUB> − <B>X</B><SUB><B>i</B></SUB>&bgr;)′<B>V</B><SUP>−1</SUP><SUB><B>i</B></SUB> (<B>y</B><SUB><B>i</B></SUB> − <B>X</B><SUB><B>i</B></SUB>&bgr;)] (7)
To obtain maximum likelihood estimates, it is necessary to find the minimum of ell  with respect to the unknown parameters. These unknown parameters are the four intercepts and slopes in the vector beta , gD and gO.

The search for a minimum of -2 log likelihood can be simplified, because for given values of gD, gO, sigma 2D and sigma 2O, beta  has a closed form solution that minimizes -2 log likelihood
<A><AC>&bgr;</AC><AC>ˆ</AC></A> = <FENCE><LIM><OP>∑</OP><LL><IT>i</IT></LL></LIM> <B>X</B>′<SUB><B>i</B></SUB> <B>V</B><SUB><B>i</B></SUB><SUP>−1</SUP> <B>X</B><SUB><B>i</B></SUB></FENCE><SUP>−1</SUP> <FENCE><LIM><OP>∑</OP><LL><IT>i</IT></LL></LIM> <B>X</B>′<SUB><B>i</B></SUB> <B>V</B><SUP>−1</SUP><SUB><B>i</B></SUB> <B>y</B><SUB><B>i</B></SUB></FENCE>
Substituting this into Eq. 7
ℓ = <LIM><OP>∑</OP><LL><IT>i</IT></LL></LIM> [log ‖<B>V</B><SUB><B>i</B></SUB>‖ + (<B>y</B><SUB><B>i</B></SUB> − <B>X</B><SUB><B>i</B></SUB><A><AC>&bgr;</AC><AC>ˆ</AC></A>)′<B>V</B><SUP>−1</SUP><SUB><B>i</B></SUB> (<B>y</B><SUB><B>i</B></SUB> − <B>X</B><SUB><B>i</B></SUB><A><AC>&bgr;</AC><AC>ˆ</AC></A>)]
which is now a function of only gD and gO. The nonlinear search need only be carried out with respect to these two parameters.

In the unbalanced case, consider the modification when one of the replications is missing. In Eq. 5, suppose the second oxygen replication is missing. The resulting equation is
<FENCE><AR><R><C><IT>y</IT><SUB>D<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>D<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C><IT>y</IT><SUB>O<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R></AR></FENCE> = <FENCE><AR><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C><C>0</C><C>0</C></R><R><C>0</C><C>0</C><C>1</C><C><IT>t</IT><SUB><IT>i</IT></SUB></C></R></AR></FENCE> <FENCE><AR><R><C>&agr;<SUB>D</SUB></C></R><R><C>&bgr;<SUB>D</SUB></C></R><R><C>&agr;<SUB>O</SUB></C></R><R><C>&bgr;<SUB>O</SUB></C></R></AR></FENCE> + <FENCE><AR><R><C>g<SUB>D</SUB></C></R><R><C>g<SUB>D</SUB></C></R><R><C>g<SUB>O</SUB></C></R></AR></FENCE> &ggr;<SUB>D</SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>) + <FENCE><AR><R><C>&egr;<SUB>D<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C>&egr;<SUB>D<SUB>2</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R><R><C>&egr;<SUB>O<SUB>1</SUB></SUB>(<IT>t</IT><SUB><IT>i</IT></SUB>)</C></R></AR></FENCE>
The covariance matrix in Eq. 6 becomes
<B>V</B><SUB><B>i</B></SUB> = <FENCE><AR><R><C>g<SUP>2</SUP><SUB>D</SUB> + &sfgr;<SUP>2</SUP><SUB>D</SUB></C><C>g<SUP>2</SUP><SUB>D</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C></R><R><C>g<SUP>2</SUP><SUB>D</SUB></C><C>g<SUP>2</SUP><SUB>D</SUB> + &sfgr;<SUP>2</SUP><SUB>D</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C></R><R><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUB>D</SUB>g<SUB>O</SUB></C><C>g<SUP>2</SUP><SUB>O</SUB> + &sfgr;<SUP>2</SUP><SUB>O</SUB></C></R></AR></FENCE>
The estimation then continues with these smaller arrays. This method can be used for any pattern of missing observations, including that pattern when all replications of deuterium or oxygen are missing. The pattern of missing observations can be different at different observation times. This methodology assumes that the missing observations are missing at random (9). This means that the reason an observation is missed is not related to the value of the missed observation.

Confidence Intervals From the Log Likelihood Surface

Maximum likelihood estimation also allows confidence intervals for the estimated parameters to be estimated directly from the -2 log likelihood surface. This method is sometimes referred to as "likelihood ratio-based confidence intervals" (see p. 289 of Ref. 12) or "profile likelihood confidence intervals" (see p. 43 of Ref. 4).

Model Eq. 4 has four linear parameters, alpha D, kD, alpha O, and kO. The main quantity of interest is
2r′<SUB>CO<SUB>2</SUB></SUB> = <IT>k</IT><SUB>O</SUB><IT>N</IT><SUB>O</SUB> − <IT>k</IT><SUB>D</SUB><IT>N</IT><SUB>D</SUB> = <IT>k</IT><SUB>O</SUB><IT>e</IT><SUP>−&agr;<SUB>O</SUB></SUP> − <IT>k</IT><SUB>D</SUB><IT>e</IT><SUP>−&agr;<SUB>D</SUB></SUP>
which is a nonlinear function of the alpha -values, and the two terms are correlated. To obtain a confidence interval on 2r'CO2, the maximum likelihood estimates of the parameters are obtained first by minimizing -2 log likelihood. The value of 2r'CO2 is varied to find the point at which -2 log likelihood increases by the significance level for chi 2 with one degree of freedom. For a 5% level of significance, this is 3.84. At each value of 2r'CO2 that is tried in the search, -2 log likelihood is minimized with respect to kO, alpha O, kD, and alpha D, under the constraint that
<IT>k</IT><SUB>O</SUB><IT>e</IT><SUP>−&agr;<SUB>O</SUB></SUP> − <IT>k</IT><SUB>D</SUB><IT>e</IT><SUP>−&agr;<SUB>D</SUB></SUP> = 2r′<SUB>CO<SUB>2</SUB></SUB>

Serial Correlation

The physiological variation can be modeled as a continuous-time first-order autoregression (see p. 92 of Ref. 6). Let delta i ti - ti - 1 be the time interval between observation time ti - 1 and observation time ti. The model is
&ggr;(<IT>t</IT><SUB><IT>i</IT></SUB>) = &ggr;(<IT>t</IT><SUB><IT>i</IT> − 1</SUB>)<IT>e</IT><SUP>−a&dgr;<SUB><IT>i</IT></SUB></SUP> + &eegr;(<IT>t</IT><SUB><IT>i</IT></SUB>)
In this model the parameter must be positive and represents the rate of decay of the physiological components toward the lines. There is one random input to the model, eta (ti) that represents the unpredictability of the physiological component between the two times. The random input has mean zero and variance, - exp(-2adelta i), that depend on the time interval. The variance of the autoregressive process is 1. Notice that the variance of the random input approaches the variance of the process sigma 2 as delta i becomes large, and approaches zero as delta i approaches zero.

Given this model for the physiological variation, the parameters of the model can be estimated using the state space methodology, as in Jones (6). To obtain the likelihood ratio test for whether the inclusion of the parameter for serial correlation is significant, the models with and without serial correlation are fitted to the data, and the change in -2 ln likelihood is tested as chi 2 with one degree of freedom. The degrees of freedom are the number of additional parameters included in the model, in this case the single parameter "a."


    ACKNOWLEDGEMENTS

We thank the reviewers for their comments.


    FOOTNOTES

The first author is partially supported by National Institute of General Medical Studies Grant GM-38519. The mass spectrometry core is supported in part by the following National Institutes of Health grants: HD-04024, HD-20716, DK-48520, and DK-49181. The Center for Human Nutrition has supported one author (B. J. Sonko) for a pilot study.

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. §1734 solely to indicate this fact.

Address for reprint requests and other correspondence: R. H. Jones, Dept. of Preventive Medicine and Biometrics, School of Medicine, Box B-119, Univ. of Colorado Health Sciences Center, 4200 East 9th Ave., Denver, CO 80262 (E-mail: rhj{at}times.uchsc.edu).

Received 28 May 1999; accepted in final form 1 October 1999.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
METHODS
EXAMPLES
DISCUSSION
REFERENCES
APPENDIX

1.   Cole, T. J., and W. A. Coward. Precision and accuracy of doubly labeled water energy expenditure by multipoint and two-point methods. Am. J. Physiol. Endocrinol. Metab. 263: E965-E973, 1992[ISI][Medline].

2.   Coward, W. A. The doubly-labeled-water (2H218O) method: principles and practice. Proc. Nutr. Soc. 47: 209-218, 1988[ISI][Medline].

3.   Coward, W. A., S. B. Roberts, and T. J. Cole. Theoretical and practical considerations in the doubly-labelled water (2H218O) method for the measurement of carbon dioxide production in man. Eur. J. Clin. Nutr. 42: 207-212, 1988[ISI][Medline].

4.   Cox, D. R., and D. Oakes. Analysis of Survival Data. London: Chapman & Hall, 1984.

5.   Elia, M. Converting carbon dioxide production to energy expenditure. In: The Doubly-Labeled Water Method for Measuring Energy Expenditure. Technical Recommendations for Use in Humans, edited by A. M. Prentice. Vienna: IDECG/IAEA, 1990.

6.   Jones, R. H. Longitudinal Data with Serial Correlation: A State-Space Approach. London: Chapman & Hall, 1993.

7.   Lifson, N. Theory of the turnover rates of body water for measuring energy and material balance. J. Theoret. Biol. 12: 46-74, 1966[ISI][Medline].

8.   Lifson, N., G. B. Gordon, and R. McClintock. Measurement of total carbon dioxide production by means of D2O18. J. Appl. Physiol. 7: 704-710, 1955[Free Full Text].

9.   Little, R. J. A. Modeling the drop-out mechanism in repeated-measures studies. J. Am. Statis. Soc. 90: 1112-1121, 1995.

10.   Prentice, A. M. (Editor). The Doubly-Labeled Water Method for Measuring Energy Expenditure. Technical Recommendations for Use in Humans. Vienna: IDECG/IAEA, 1990.

11.   SAS Institute, Inc. SAS/STAT User's Guide, Release 6.03 Edition. Cary, NC: SAS, 1988, p. 1028.

12.   SAS Institute, Inc. Software: Changes and Enhancements through Release 6.12. Cary, NC: SAS, 1997, p. 1167.

13.   Schoeller, D. A. Measurement of energy expenditure in free living humans by using doubly labeled water. J. Nutr. 118: 1278-1289, 1988[ISI][Medline].

14.   Wier, J. B. de V. New methods for calculating metabolic rate with special reference to protein metabolism. J. Physiol. (Lond.) 109: 1-9, 1949[ISI].


Am J Physiol Endocrinol Metab 278(3):E383-E389
0193-1849/00 $5.00 Copyright © 2000 the American Physiological Society