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.

The methods described here are implemented in packages breathtestcore and breathtestshiny, and can be tested online here

\===

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
(**???**).

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 |

Estimated parameters from Stan for the 3 records with invalid fits. t50BW: t50 self-corrected after Bluck-Coward; t50Maes: t50 following Ghoos/Maes.

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.

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): 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.