Summary

Time series from 13C breath test recordings are used as a surrogate indicator for gastric emptying. Gastric emptying half time t50 is computed by fitting the series with an exponential beta or a Gamma function. Often, records are too short to obtain stable estimates of the trailing part. One can easily get unrealistic emptying times when single curves are fitted (Ghoos et al. (1993), Maes et al. (1998)).

Population methods can come to the rescue. When these fail, Bayesian methods provide a robust approach to non-linear curve fitting.

Zipped source files and data

Single curve fit

Single curves fits - the standard method used in breath test evaluation - can be computed with function nlsList from package nlme. The exponential beta ExpBetafunction is defined in package D13CBreath; use devtools::install_github("dmenne/d13cbreath") to install. To compute the half-emptying time t50, use function t50BluckCoward or t50Maes from the same packages. The exponential beta function is parameterized as k = 1/tempt, because fits are close to normal distribution when k is estimated instead of tempt. The dose is fixed to 100 (mg); with this choice, factor m can be interpreted as the fraction of 13C metabolized over the complete extrapolated emptying curve.

start = c(m = 30,k = 0.01,beta = 1.8)
d_nls <- nlsList(
  PDR ~ ExpBeta(Time, 100, m, k, beta) | BreathTestRecordID,
  data = d,start = start
)
cf_nls =  coef(d_nls) %>%  
  rownames_to_column("BreathTestRecordID") %>% 
  mutate(
    method = "nls",
    t50BW = t50BluckCoward(.),
    t50Maes = t50Maes(.)
  ) 
cf_nls
   BreathTestRecordID  m      k beta method t50BW t50Maes
1                P435 49 0.0052  1.5    nls  28.7     195
2                P436 53 0.0055  1.6    nls  30.8     189
3                P443 63 0.0115  1.8    nls  24.0     100
4                P444 36 0.0057  1.9    nls  55.4     209
5                P446 NA     NA   NA    nls    NA      NA
6                P455 NA     NA   NA    nls    NA      NA
7                P457 13 0.0107  3.8    nls  84.1     167
8                P475 NA     NA   NA    nls    NA      NA
9                P485 27 0.0107  1.4    nls   6.6      87
10               P487 20 0.0098  2.0    nls  36.8     127
11               P489 NA     NA   NA    nls    NA      NA
12               P496 34 0.0080  2.2    nls  51.8     163

With single curve fits, the underlying function nls does not produce valid fits for 4 records. Packages like MathPad or Matlab might give solutions, but these are most likely associated with huge errors of the coefficients. R uses a much more conservative approach and returns NA as a feature, not a bug (Douglas Bates, personal communication).

Population fit with nlme

The population fit with function nlme from package nlme uses the full information from all curves, giving more stable estimates from "borrowing strength" providing guidance in critical cases such as P443. However, fitting all data with nlme fails; by trial-and-error one can find out that three records must be excluded to force convergence:

d_part = d[!(d$BreathTestRecordID %in% c("P446","P489","P455")), ]
d_nlme = nlme(PDR~ExpBeta(Time, 100, m, k, beta),
           data = d_part,
           fixed = m + k + beta~1,
           random = m + pdDiag(k + beta) ~ 1,
           groups = ~BreathTestRecordID,
           control = nlmeControl(pnlsTol = 0.05),
           start = start) 

In addition to removing the three offending record, the value of pnlsTol has been increased from its default of 0.001. Always use the smallest value that leads to convergence; trial-and-error is required; function nlme_gastempt in my package gastempt automatically adjust pnlsTol until convergence. Values above 0.5 are dubious and indicate that the results might not be trusted.

The computed coefficients from nlme and derived quantity t50 (emptying half time in minutes):

cf_nlme = coef(d_nlme)  %>%  
  rownames_to_column("BreathTestRecordID") %>% 
  mutate(
    method = "nlme",
    t50BW = t50BluckCoward(.),
    t50Maes = t50Maes(.)
  )
cf_nlme
  BreathTestRecordID  m      k beta method t50BW t50Maes
1               P435 48 0.0053  1.6   nlme  29.4     193
2               P436 53 0.0056  1.6   nlme  31.4     187
3               P443 63 0.0115  1.8   nlme  23.9     101
4               P444 35 0.0060  2.0   nlme  56.0     203
5               P457 17 0.0078  2.8   nlme  81.2     195
6               P475 52 0.0029  2.5   nlme 191.6     501
7               P485 27 0.0106  1.4   nlme   6.7      87
8               P487 20 0.0097  2.0   nlme  36.7     127
9               P496 34 0.0079  2.2   nlme  51.6     163

Bayesian fit with Stan

Stan is a comprehensive Open-Source package for Bayesian analysis. For a different evaluation of breath test data using Bayesian methods, see also Bluck et al. (2011).

The following Stan model was used:

# All parameters have hard lower limit at 0
m ~ normal(20,25);
beta ~ normal(2,0.4);
k ~ normal(0.0015,0.001);
sigma ~ cauchy(0,5);

for (n in 1:N){
  rec = Record[n];
  mn  =  m[rec];
  kn = k[rec];
  ktn = kn* Time[n];
  betan = beta[rec];
  # Dose fixed to 100
  pdr = 100*mn*kn*betan*exp(-ktn) * pow(1 - exp(-ktn),(betan -1));
  PDR[n] ~ normal(pdr,sigma);
}

The model uses = for assignment which requires a version of Stan >= 2.10; earlier versions used <-.

Data are passed to Stan as a list.

data = list(N = nrow(d), NRecord = NRecord, Record = d$Record,
  Time  = d$Time, PDR = d$PDR)
List of 5
 $ N      : int 438
 $ NRecord: int 12
 $ Record : int [1:438] 1 1 ...
 $ Time   : num [1:438] 1.5 8.7 ...
 $ PDR    : num [1:438] 4 7.5 ...

The number of data N and the number of records NRecord for indexing of the parameter block must be passed to allocate the data vectors

data{  
  int<lower=0> N; # Number of data
  int<lower=0> NRecord; # Number of records
  int<lower=0> Record[N]; # From BreathTestRecordID
  real<lower=0> Time[N]; 
  real<lower=-30> PDR[N];
}

These following parameters are computed, one for each gastric emptying time curve (=record):

parameters{
  real <lower=0> m[NRecord];
  real <lower=0> k[NRecord];
  real <lower=0> beta[NRecord];
  real <lower=0> sigma;  
}

Stan model compilation and running needs a few minutes, therefore the result is explicitly cached as follows:

stanCache = "BreathTestStan1.rdata"
if (!file.exists(stanCache)) {
  mr.stan = stan("BreathTestStan1.stan", seed = 123, chains = chains, 
                 iter = 200, data = data)
  save(mr.stan,file = stanCache)
} else {
  load(stanCache)
}

4 chains with random initial values are run, each with 100 warmup iterations at a total of 200 iterations.

Stan diagnostics

The left part of the panels with gray background indicate the warmup ("burn-in" in non-Gelman based literatur) range of 100 iterations. Stan has converged after about 40 iteration for all chains

Tracplot of sigma and parameters of the first record

Predictions

All fits and raw data (gray).

To better separate overlapping curves, these are vertically shifted against each other by 0.5 unit.

  • Population fits and Bayes fits are very similar. Bayes and nlme (population) fits do not provide an advantage when the nls fit converges.
  • Since single-curve fits are mostly reliable, the fit is well defined with only little correlation between the parameters.

Things can go wrong

  • For three record (P446, P455, P489), only Stan gives valid fits.
  • For one record (P475), Stan and nlme give valid fits

Records where single fits do not give valid results

  • In extreme cases, only the Bayesian method returns an estimate.
  • The self-corrected estimate of t50 following Bluck-Coward is definitively too low.
  • The estimate of t50 following Maes/Ghoos is too high for some cases (> 10 h).
  • The fact that curves can be fitted does not mean that the fits give medically meaningful results.
method BreathTestRecordID m k beta t50Maes t50BW
Stan P446 117 0.0042 2.1 298 88
Stan P455 61 0.0101 1.3 87 4
Stan P475 58 0.0025 2.3 537 187
Stan P489 40 0.0020 2.0 618 177
nlme P475 52 0.0029 2.5 501 192

Wagner-Nelson to the rescue??

A popular alternative when curve fitting does not work is the Wagner-Nelson method which uses non-parametric approaches for the initial slope. However, it makes a assumption about the terminal slope being the same for all subject (k=0.01/min, 0.65/h). Since this assumption strongly affects the estimate of the half-emptying time, there is little justification to use the Wagner-Nelson method.

Citations

Bluck, Leslie J C., Sarah J. Jackson, Georgios Vlasakakis, and Adrian Mander. 2011. “Bayesian Hierarchical Methods to Interpret the (13)C-Octanoic Acid Breath Test for Gastric Emptying.” Digestion 83 (1-2). MRC Human Nutrition Research, Cambridge, UK. les.bluck@mrc-hnr.cam.ac.uk: 96–107. doi:10.1159/000316823.

Ghoos, Y. F., B. D. Maes, B. J. Geypens, G. Mys, M. I. Hiele, P. J. Rutgeerts, and G. Vantrappen. 1993. “Measurement of Gastric Emptying Rate of Solids by Means of a Carbon-Labeled Octanoic Acid Breath Test.” Gastroenterology 104 (6). Department of Medicine, University Hospital Gasthuisberg, Belgium.: 1640–7.

Maes, B. D., B. J. Geypens, Y. F. Ghoos, M. I. Hiele, and P. J. Rutgeerts. 1998. “13C-Octanoic Acid Breath Test for Gastric Emptying Rate of Solids.” Gastroenterology 114 (4): 856–59.

Previous Post