5.1 Types of error

5.2 Variance and covariance

5.3 Ratio plots and product plots

5.4 The variance of CO_{2}production

5.5 Optimal design to minimise the variance

5.6 Variance and the two-point method

5.7 Bias in the estimation of pool size

5.8 Precision and accuracy in the multi-point and two-point methods

5.9 Constancy of errors

5.10 Data and residual plots

5.11 References

Contributor:

Tim Cole

Mike Franklin

Andy Coward

The calculations required to obtain
the pool sizes and flux rates of the two isotopes have been described in Chapter 4. These
quantities together with the relevant fractionation rates, lead to an estimate of the CO_{2}
production rate in individual subjects.

To be of most value, each estimate
of CO_{2} production rate requires an associated standard error, to show how
closely the data satisfy the assumptions made by the method, and to provide reassurance
that the error is small. The calculation of the error assumes that the data provide an
adequate fit to the model of constant flux rate and constant pool size, so that the error
is a measurement of precision. Quite separately there is the possibility of bias in the
estimate, due to imperfections in the model, which defines the accuracy of the estimate.

The fundamental equation for CO_{2}
production rate ignoring fractionation is:

_{}

.........1

where k_{O} and k_{D}
are the flux rates of oxygen and hydrogen respectively, and N_{O} and N_{D }are
the body water pool sizes as estimated by ^{18}O and ^{2}H. (Note that the
use of a prime in r'_{CO2} signifies production rates uncorrected for
fractionation - see Chapter 2). The four quantities are assumed to be constant, and are
estimated from the two regression lines of log isotope enrichment versus time. So long as
there are at least 3 data points the estimates have standard errors, based on the
departures from linearity of the observed log enrichments. It is assumed that the timing
of each sample point is known to full accuracy, although in practice this may not be
strictly true.

The departures from linearity are partly physiological, representing genuine changes in flux from day to day, and partly instrument measurement error. In the multi-point method, the two sources of error are combined in the residual variation about the regression line, whereas in the two-point method there is no residual variation.

__5.2.1 Slope-intercept
covariance__

The regression of log isotope
enrichment on time can take several forms, depending on whether the errors are assumed
constant in enrichment units, or constant in coefficient of variation terms, or somewhere
in between (see Section 5.9 and Chapter 4). In the simplest case a constant coefficient of
variation is assumed, and k_{O} and k_{D} are estimated as the unweighted
linear regression coefficients of log enrichment versus time. Each regression coefficient
has a standard error which is obtained from the regression analysis.

Similarly, estimates of I_{O}
and I_{D} are obtained from the exponentials (or antilog) of the intercepts from
these regressions, and again they have their corresponding standard errors. It is
convenient to express the enrichments of ^{18}O and ^{2}H as the fraction
of dose per mole of body water. Then the intercepts I_{O} and I_{D} are by
definition the reciprocals of the body water pools N_{O} and N_{D}, ie I_{O}
= 1/N_{O} and I_{D} = 1/N_{D}.

Given the standard errors for the
flux rates and pool sizes, it might be thought straightforward to combine the errors,
using Equation 1 to obtain a standard error for the CO_{2} production. However
there are complications, because of the way that the four parameters interact with each
other.

Firstly the estimates of the slope k and intercept log I from each regression line are highly correlated (Fig 5.1). The fitted regression line passes, by definition, through the mean value of the data points, so that the line can pivot about its mean quite independently of the latter's location. Thus the mean log enrichment and the regression slope are uncorrelated with each other. However when the slope changes the intercept changes in the same direction - the slope k and intercept log I are correlated. (It should be emphasised that k is treated here as being positive even though it is actually negative - the negative sign is included in the equations). Thus the standard error of the ratio k/I (or equivalently the product kN) needs to reflect this covariance relationship.

From Fieller's theorem ^{1}
the variance of k/I is given by:

_{}

.........2

where Var(k) and Var(logI) are the variances of k and logI, and r* is the correlation between k and logI. This correlation is given by:

_{}

.........3

where _{} is the arithmetic mean of the sampling times and _{} is the root mean square of the sampling
times (ie square the sampling times, add them together, divide by the number and take the
square root). This correlation lies between zero and +1, and is a function only of the
chosen sampling times, not of the measured enrichments. This takes account of the
covariance between k and I for each isotope, but there is the further problem that after
fitting the two regression lines, the two residuals at a particular time tend to be
correlated.

__5.2.2 Correlated residuals__

If the observed ^{18}O
enrichment is greater than the value predicted by the regression line, there is a high
probability that the ^{2}H enrichment will also be elevated. This correlation
arises in two ways. Firstly, ^{2}H and ^{18}O are both influenced by water
turnover, but only ^{18}O is affected by CO_{2} production. So if the
water turnover from one sampling time to the next is anomalously high say, then the latter
hydrogen enrichment is depressed. The oxygen enrichment is also depressed, but less than
for the hydrogen because only a fraction of the ^{18}O is lost as water. The
converse is true for low water turnovers. So, depending on the fraction of oxygen lost as
water, the deviations of ^{2}H and ^{18}O from their fitted lines will be
correlated to a greater or lesser extent.

The second source of the correlation
is physiological, in that if energy expenditure (ie CO_{2} production) increases
sufficiently, there is likely to be a corresponding increase in water turnover through
sweating. Thus both the CO_{2} and H_{2}O components of ^{18}O may
rise together, depending on the precise nature of the increase in expenditure. In practice
though, this source is likely to be far smaller than that due to the slope-intercept
interaction.

The estimate of error for CO_{2}
production given by Coward *et al* ^{2} fails to take either of these sources
of correlation between the residuals into account. The correlation cannot be predicted on
a theoretical basis, and the best method of incorporating it in the estimate of the
variance of CO_{2} production is by no means clear.

CO_{2} production rate can
be thought of as ^{18}O turnover adjusted for ^{2}H turnover. Thus much of
the information about CO_{2} production is contained in the difference between the
isotope flux rates (k_{O}-k_{D}). This can be estimated directly from the
log regression of the enrichment ratio d_{O}/d_{D} plotted against time. The appeal of this
'ratio plot' regression is that the errors in ^{2}H and ^{18}O caused by
water turnover fluctuations nearly cancel each other out, so that the fit to the line is
much closer. In fact the linearity of the plot is simply a measure of the constancy of CO_{2}
production.

To make use of the ratio plot,
Equation 1 defining CO_{2} production has to be recast in terms of the ratio plot
slope and intercept, k_{r} and I_{r} say. (Note that although k_{r}
= k_{O}-k_{D}, the relationship between I_{r}, I_{O} and I_{D}
is not a simple ratio). Results from the analogous 'product plot', ie d_{O}-d_{D} plotted
on a log scale against time, are also required - here the slope and intercept are k_{p}
and I_{p}.

Although the correlation between
residuals is much reduced by analysing the ratio and product plots, it is not removed
entirely. There may still be the correlation between water turnover and CO_{2}
production discussed earlier, albeit on a much smaller scale.

The assumption of exponential disappearance implies that:

_{}

.........4

where the suffix (*) is either D or O. The corresponding equations for the ratio plot and product are, respectively:

_{}

.........5

and

_{}

.........6

Equation 1 can then be expressed in
terms of I_{r}, he, I_{p} and k_{p} as follows:

_{}

..........7

where

_{}

.........8

and:

_{}

.........9

Using Equation 2 the variance of Equation 7 is given by:

_{}

..........10

where:

_{}.........11

_{}.........12

and r* is the slope-intercept correlation given in Equation 3.

Similarly the Var (*) indicate the
variances (i.e. the squares of the standard errors) for each of the parameters h_{r},
k_{p}, log I_{r} and log I_{p}.

Equation 10 is rather intimidating.
Fortunately it can be reformulated in terms of the residual variances of the two
regression lines _{} and _{}, and the chosen sampling times t_{i},
where i = 1...n:

_{}

.........13

Here _{} is the mean of the sampling times t_{i} and St² is the corrected sum of squares of the t_{1}, where:

_{}

The equation consists of terms which
depend on n, _{} and St² It is clear that if both n and St²
are large, then the variance of CO_{2} production will be small.

Equation 13 gives the required
variance, but in addition it also defines the optimal choice of sampling times to minimise
the variance. For example if the mean of sampling times, _{}, is arranged to exactly equal A_{r}/B_{r},
then the first term in Equation 13 becomes zero; alternatively if t = A_{p}/B_{p}
the second term vanishes.

The ratio A_{r}/B_{r}
is given by (N_{O} + N_{D})/(k_{O}N_{O} + k_{D}N_{D}),
which is the reciprocal of a weighted sum of the flux rates k_{O }and k_{D}.
Similarly Ap/Bp, which is the same as (N_{O} - N_{D})/(k_{O}N_{O}
- k_{D}N_{D}), is a reciprocal of a weighted difference of the two flux
rates.

The importance of Equation 13 cannot
be overstated - it indicates how to design doubly-labelled water experiments optimally,
providing the smallest possible variance for the CO_{2} production. Typical values
for A_{r}, B_{r}, S_{r} and A_{p}, B_{p}, S_{p},
are first required for the subjects under study, and then by suitable choice of the
sampling times, n and _{} and St² can be set to minimise Equation 13.

Based on the 30 sets of data
provided for the workshop, the ratio A_{r}/B_{r} averages about 10 days
while A_{p}/B_{p} is about -1 days. This poses a problem, for in order to
make both the second and third terms in Equation 13 vanish, t needs to be equal to 10 and
-1 simultaneously.

The best choice of t is a compromise
between these two values, the exact choice depending on the relative sizes of (B_{r}S_{r})²
and (B_{p}S_{p})². If they are similar in size, then t should be set to
the average, ie 4 or 5 days. The terms S_{r}² and S_{p}² represent the
noise about the fitted ratio and product lines, and S_{r}² is typically a tenth
of S_{p}² (although it varies widely about this figure). However Br is very much
greater than B_{p}, so that in practice the ratio (B_{r}S_{r})²/(B_{p}S_{p})²
for the workshop datasets is close to 6. This makes the optimal choice for t equal to 7.7
days, as compared with the actual mean of 7.5 days, indicating that the experiments were
close to optimally designed on the whole.

This discussion has assumed that the
flux rate k is known in advance of the experiment, which of course it is not. The optimal
value for t can only be chosen approximately. To ensure that the second and third terms in
Equation 13 are kept small, the value of St² should be as
large as possible, which requires the sampling times to be widely spaced. The possible
range of times is bounded in both directions of course, at time zero at one end and the
time when the signal gets lost in background at the other, so that the optimal design has
one group of samples near the start and another group near the end, arranged symmetrically
around _{}.

The first term of Equation 13 is inversely proportional to the reciprocal of the number n of sampling times in the experiment. In the workshop datasets, where n averages 14, almost exactly 50% of the total variance was contributed by this term. If n were to be reduced to 5 (while keeping the other variables unchanged) then the first term's contribution would be increased 3 times and the total variance doubled. Put the other way, this indicates that using 14 sampling points, as opposed to 5, halves the variance emphasising the improvement in precision that multi-point sampling provides.

Another consideration is the timing
of the first sample. It is argued above that the first sample should be near to zero time,
but in fact this is not a strong requirement from the precision point of view. The term St² is only of secondary importance compared to n and _{}, and there are strong practical reasons
for avoiding a sample immediately after dosing. According to Equation 13, a first sample
as late as 1 or even 2 days after dosing would not have much impact on the variance of CO_{2}
production so long as _{} is chosen
appropriately. Such a delay means that the pool size is being determined through
extrapolation and hence increases the importance of selecting the correct model.

The use of Equations 10 or 13 to
calculate the variance of the CO_{2} production rate requires knowledge of the
residual variances about the ratio and product regression lines. This precludes its
application to two-point data. However it is quite valid to substitute the known
instrument measurement errors into Equation 13, where the errors on ^{2}H and ^{18}O,
d_{D}
and d_{O}
are independent so they combine vectorially. Thus:

_{}

.........14

This variance is only a fraction of
the total variance, and takes no account of day-to-day variation in water and CO_{2}
production. Thus the variance estimated from Equation 13 will be a substantial
underestimate of the true variance, and for this reason it cannot be used for experimental
design purposes.

The calculation of CO_{2}
production rate (Equation 1) assumes that the flux rates of ^{2}H and ^{18}O
are constant, leading to linear plots of log enrichment against time. It is also valid if
the plots are nonlinear so long as the body water pool size N remains constant during the
experiment. Water volume is likely to be constant with adults in non-clinical situations
unless they are dieting or there is serious dehydration, but in small children who are
growing rapidly or in pregnant women there is a chance that the water volume will increase
during the study and introduce a bias. This bias shows itself as a curved plot of log
enrichment versus time.

The pool sizes, obtained from zero
time enrichments, are estimated from the intercepts of the fitted regression lines. In
practice the plots may not be linear, but may oscillate systematically above and below the
fitted lines. This phenomenon is called *serial correlation* or *auto-correlation*.
When this happens, the intercepts are biased estimates of the zero time enrichments,
leading to a corresponding bias in the estimate of CO_{2} production.

Systematic non-linearity in the
data, and the magnitude of the bias, can be seen most easily from residual plots of the
data (see figures in Chapter 11). The plots for ^{18}o and ^{2}H show the
biases for the two water spaces separately, but the biases can also be obtained from the
'product' and 'ratio' plots. The advantage of using the product: and ratio plots is that
the bias is partitioned into an average pool size effect (bias in I_{p}) and an
effect due to the N_{D}/N_{O} ratio (bias in I_{r}). When the
biases in the two water spaces are covariant, as they often are, the resulting bias in CO_{2}
production rate is contained largely in the product plot, while the ratio plot bias is far
smaller.

To see how the biases in I_{p}
and I_{r} relate to the bias in

r'_{CO2.} Equation 7 can be
rearranged to give:

_{}

.........15

Differentiating Equation 15 and rearranging leads to:

_{}

.........16

where dr'_{CO2}, dI_{p} and dI_{r} are the
fractional or percentage biases in r'_{CO2}, I_{p} and I_{r}
respectively. If dI_{r} is zero, Equation 16 shows that dI_{p} needs to be
halved and its sign changed to give the bias r'_{CO2} in CO_{2} production
rate. From the analyses performed prior to the workshop we recommend that the value for dI_{p} is
given by the residual from the product plot about 5 to 8 hours post-dose.

Similarly, the bias in I_{r}
is shown by a series of residuals from the ratio plot between 5 and 8 hours post-dose
which are consistently non-zero. As there is a lot of noise on the ratio plot, a single
large residual during this time should not be treated as bias, it may equally be due to
measurement error. However if bias is present (i.e. dI_{r} is non-zero), its contribution to Equation
16 needs to be taken into account. The value of the multiplier (k_{O}N_{O}
+ k_{D}N_{D})/(k_{O}N_{O} - k_{D}N_{D}) is
typically 9 or a little more, so that dI_{r} should be multiplied by 9/2 to give the bias
in r'_{CO2}.

Thus a bias of say +2% in I_{p}
leads to a -1% bias in r'_{CO2}, whereas a +2% bias in I_{r} biases r'_{CO2}
by -9%. In practice it is difficult to detect bias in the intercept of the ratio plot.
Nevertheless Equation 16 makes clear that if the N_{D}/N_{O} ratio is at
all aberrant, its percentage effect on r'_{CO2} is magnified up some 4 or 5 times.

Unlike the multi-point method, the two-point method explicitly uses an early enrichment to derive the pool size. As a result it avoids the possibility of bias in the estimation of pool size. However this gain in accuracy relative to the multi-point method is offset by a corresponding loss in precision, owing to the smaller number of samples.

Chapter 9 gives further coverage of
the potential errors introduced by systematic and random variations in water and CO_{2}
flow rates, and by errors in estimating N_{D}/N_{O}.