Time series from ^{13}C 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.

Single curves fits - the standard method used in breath test evaluation - can be computed with function `nlsList`

from package `nlme`

. The exponential beta `ExpBeta`

function 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 ^{13}C 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).

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
```

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.

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

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.

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

- 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 |

The Wagner-Nelson method does not directly use the fit, but could come to the rescue here: it either makes an strong assumption of the terminal slope (`k=0.01/min, 0.65/h`

) or needs an estimate from a range of the curve far outside the recording.

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.