Skip to contents

Simplified population fit of gastric emptying data

Source: R/fitnlme.R
nlme_gastempt.Rd

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.2572  38.26755 0.9614860 62.74108  4.211677 38600.8
#> 2  rec_02 494.2277  29.06142 0.7791640 62.74108  4.211677 38600.8
#> 3  rec_03 474.9852  65.59146 0.9788760 62.74108  4.211677 38600.8
#> 4  rec_04 409.2757  70.62709 0.6727428 62.74108  4.211677 38600.8
#> 5  rec_05 364.8369  69.44775 0.4270034 62.74108  4.211677 38600.8
#> 6  rec_06 342.5338 104.75216 0.6047117 62.74108  4.211677 38600.8
#> 7  rec_07 438.8791  48.00223 0.6406530 62.74108  4.211677 38600.8
#> 8  rec_08 370.9151  72.25262 0.4835091 62.74108  4.211677 38600.8
#> 9  rec_09 396.2782  59.61221 0.5051639 62.74108  4.211677 38600.8
#> 10 rec_10 446.1316  82.80695 0.9666869 62.74108  4.211677 38600.8
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.961 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.3 0.484 62.7      4.21 38600
#> 9  rec_09 396  59.6 0.505 62.7      4.21 38600
#> 10 rec_10 446  82.8 0.967 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