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.
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
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 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).
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.
nls
fit
converges.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.