Skip to contents

Compute coefficients v0, tempt and kappa of a mixed model fit to a linexp function with one grouping variable

Usage

nlme_gastempt(d, pnlsTol = 0.001, model = linexp, variant = 1)

Arguments

d

A data frame with columns

  • record Record descriptor as grouping variable, e.g. patient ID

  • minute Time after meal or start of recording.

  • vol Volume of meal or stomach

pnlsTol

The value of pnlsTol at the initial iteration. See nlmeControl When the model does not converge, pnlsTol is multiplied by 5 and the iteration repeated until convergence or pnlsTol >= 0.5. The effective value of pnlsTol is returned in a separate list item. When it is known that a data set converges badly, it is recommended to set the initial pnlsTol to a higher value, but below 0.5, for faster convergence.

model

linexp (default) or powexp

variant

For both models, there are 3 variants

  • variant = 1 The most generic version with independent estimates of all three parameters per record (random = v0 + tempt + kappa ~ 1 | record). The most likely to fail for degenerate cases. If this variant converges, use it.

  • variant = 2 Diagonal random effects (random = pdDiag(v0 + tempt + kappa) ~ 1; groups = ~record ). Better convergence in critical cases. Note: I never found out why I have to use the groups parameter instead of the |; see also p. 380 of Pinheiro/Bates.

  • variant = 3 Since parameters kappa and beta respectively are the most difficult to estimate, these are fixed in this variant (random = v0 + tempt ~ 1). This variant converges in all reasonable cases, but the estimates of kappa and beta cannot be use for secondary between-group analysis. If you are only interested in t50, you can use this safe version.

Value

A list of class nlme_gastempt with elements coef, summary, plot, pnlsTol, message

  • coef is a data frame with columns:

    • record Record descriptor, e.g. patient ID

    • v0 Initial volume at t=0

    • tempt Emptying time constant

    • kappa Parameter kappa for model = linexp

    • beta Parameter beta for model = powexp

    • t50 Half-time of emptying

    • slope_t50 Slope in t50; typically in units of ml/minute

    On error, coef is NULL

  • nlme_result Result of the nlme fit; can be used for addition processing, e.g. to plot residuals or via summary to extract AIC. On error, nlme_result is NULL.

  • plot A ggplot graph of data and prediction. Plot of raw data is returned even when convergence was not achieved.

  • pnlsTol Effective value of pnlsTo after convergence or failure.

  • message String "Ok" on success, and the error message of nlme on failure.

Examples

suppressWarnings(RNGversion("3.5.0"))
set.seed(4711)
d = simulate_gastempt(n_record = 10, kappa_mean = 0.9, kappa_std = 0.3,
                      model = linexp)$data
fit_d = nlme_gastempt(d)
# fit_d$coef # direct access
coef(fit_d) # better use accessor function
#>    record       v0     tempt     kappa      t50 slope_t50      auc
#> 1  rec_01 514.2462  38.26585 0.9615993 62.74269  4.211766 38600.49
#> 2  rec_02 494.2109  29.05876 0.7793491 62.74269  4.211766 38600.49
#> 3  rec_03 474.9858  65.59160 0.9788704 62.74269  4.211766 38600.49
#> 4  rec_04 409.2693  70.62264 0.6728595 62.74269  4.211766 38600.49
#> 5  rec_05 364.8323  69.44100 0.4271404 62.74269  4.211766 38600.49
#> 6  rec_06 342.5194 104.72522 0.6051172 62.74269  4.211766 38600.49
#> 7  rec_07 438.8806  48.00304 0.6406247 62.74269  4.211766 38600.49
#> 8  rec_08 370.9063  72.24122 0.4837465 62.74269  4.211766 38600.49
#> 9  rec_09 396.3088  59.63892 0.5044898 62.74269  4.211766 38600.49
#> 10 rec_10 446.1476  82.81441 0.9664679 62.74269  4.211766 38600.49
coef(fit_d, signif = 3) # Can also set number of digits
#>    record  v0 tempt kappa  t50 slope_t50   auc
#> 1  rec_01 514  38.3 0.962 62.7      4.21 38600
#> 2  rec_02 494  29.1 0.779 62.7      4.21 38600
#> 3  rec_03 475  65.6 0.979 62.7      4.21 38600
#> 4  rec_04 409  70.6 0.673 62.7      4.21 38600
#> 5  rec_05 365  69.4 0.427 62.7      4.21 38600
#> 6  rec_06 343 105.0 0.605 62.7      4.21 38600
#> 7  rec_07 439  48.0 0.641 62.7      4.21 38600
#> 8  rec_08 371  72.2 0.484 62.7      4.21 38600
#> 9  rec_09 396  59.6 0.504 62.7      4.21 38600
#> 10 rec_10 446  82.8 0.966 62.7      4.21 38600
# Avoid ugly ggplot shading (not really needed...)
library(ggplot2)
theme_set(theme_bw() + theme(panel.spacing = grid::unit(0,"lines")))
# fit_d$plot  # direct access is possible
plot(fit_d) # better use accessor function