Simplified population fit of gastric emptying data

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
#> # A tibble: 10 × 7
#> # Rowwise: 
#>    record    v0 tempt kappa   t50 slope_t50    auc
#>    <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl>  <dbl>
#>  1 rec_01  514.  38.3 0.961  62.7      4.21 38601.
#>  2 rec_02  494.  29.1 0.779  42.1      5.39 25554.
#>  3 rec_03  475.  65.6 0.979 109.       2.27 61652.
#>  4 rec_04  409.  70.6 0.673  94.2      1.87 48352.
#>  5 rec_05  365.  69.4 0.427  74.3      1.86 36156.
#>  6 rec_06  343. 105.  0.605 132.       1.07 57579.
#>  7 rec_07  439.  48.0 0.641  62.3      2.97 34564.
#>  8 rec_08  371.  72.3 0.484  81.5      1.76 39757.
#>  9 rec_09  396.  59.6 0.505  68.7      2.26 35557.
#> 10 rec_10  446.  82.8 0.967 136.       1.69 72655.
coef(fit_d, signif = 3) # Can also set number of digits
#> # A tibble: 10 × 7
#> # Rowwise: 
#>    record    v0 tempt kappa   t50 slope_t50   auc
#>    <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl> <dbl>
#>  1 rec_01   514  38.3 0.961  62.7      4.21 38600
#>  2 rec_02   494  29.1 0.779  42.1      5.39 25600
#>  3 rec_03   475  65.6 0.979 109        2.27 61700
#>  4 rec_04   409  70.6 0.673  94.2      1.87 48400
#>  5 rec_05   365  69.4 0.427  74.3      1.86 36200
#>  6 rec_06   343 105   0.605 132        1.07 57600
#>  7 rec_07   439  48   0.641  62.3      2.97 34600
#>  8 rec_08   371  72.3 0.484  81.5      1.76 39800
#>  9 rec_09   396  59.6 0.505  68.7      2.26 35600
#> 10 rec_10   446  82.8 0.967 136        1.69 72700
# 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