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.2601  38.26802 0.9614559 62.74067  4.211651 38600.89
#> 2  rec_02 494.2307  29.06192 0.7791296 62.74067  4.211651 38600.89
#> 3  rec_03 474.9859  65.59174 0.9788657 62.74067  4.211651 38600.89
#> 4  rec_04 409.2770  70.62793 0.6727204 62.74067  4.211651 38600.89
#> 5  rec_05 364.8374  69.44832 0.4269907 62.74067  4.211651 38600.89
#> 6  rec_06 342.5370 104.75822 0.6046204 62.74067  4.211651 38600.89
#> 7  rec_07 438.8781  48.00173 0.6406705 62.74067  4.211651 38600.89
#> 8  rec_08 370.9166  72.25445 0.4834700 62.74067  4.211651 38600.89
#> 9  rec_09 396.2707  59.60555 0.5053315 62.74067  4.211651 38600.89
#> 10 rec_10 446.1291  82.80581 0.9667218 62.74067  4.211651 38600.89
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.483 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