Fit mNIRS Kinetics Model

Perform parametric curve fitting or non-parametric estimation of mNIRS kinetics on vector data.

Usage

process_multi_kinetics(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  event_channel = NULL,
  event_times = NULL,
  event_labels = NULL,
  event_indices = NULL,
  group_events = list(c("distinct", "ensemble")),
  method = list(c("monoexponential", "sigmoidal", "half_time", "peak_slope")),
  fit_span = list(c(30, 180)),
  end_fit_span = list(c(20)),
  x0 = list(c(0)),
  time_from_zero = TRUE,
  display_mmss = FALSE,
  verbose = TRUE,
  ...
)

# S3 method for class 'monoexponential'
process_kinetics(
  y,
  x = NULL,
  data = NULL,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  end_fit_span = 20,
  x0 = 0,
  verbose = TRUE,
  ...
)

# S3 method for class 'sigmoidal'
process_kinetics(
  y,
  x = NULL,
  data = NULL,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  end_fit_span = 20,
  x0 = 0,
  verbose = TRUE,
  ...
)

# S3 method for class 'half_time'
process_kinetics(
  y,
  x = NULL,
  data = NULL,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  end_fit_span = 20,
  x0 = 0,
  verbose = TRUE,
  ...
)

# S3 method for class 'peak_slope'
process_kinetics(
  y,
  x = NULL,
  data = NULL,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  end_fit_span = 20,
  x0 = 0,
  verbose = TRUE,
  ...
)

Arguments

data

An optional dataframe containing the predictor and response variables named in x and y. Names for x and y must be in quotations.

nirs_channels

A character vector indicating the mNIRS data channels to be processed from your dataframe. Must match data column names exactly. Will be taken from metadata if not defined explicitly.

time_channel

A character string indicating the time or sample channel name. Must match column names in data exactly. Will be taken from metadata if not defined explicitly.

event_channel

A character string indicating the event or lap channel name. Must match column names in data exactly. Will be taken from metadata if not defined explicitly.

event_times

An optional numeric vector corresponding to values of time_channel indicating the start of kinetic events. i.e., by time value or sample number.

event_labels

An optional character vector corresponding to values of event_channel indicating the start of kinetics events. i.e., by an event label such as "end work".

event_indices

An optional numeric vector indicating the starting row indices of kinetics events. i.e., to identify the start of kinetic events by row number.

group_events

Indicates how kinetics events should be analysed. Typically either "distinct" (the default) or "ensemble", but can be manually specified as a list of event numbers (see Details).

method

Indicates which model to evaluate the kinetics event (see Details for method parametrisation).

method = "monoexponential"

A four-parameter monoexponential association function in the form: ifelse(x <= TD, A, A + (B - A) * (1 - exp((TD - x) / tau))).

method = "sigmoidal"

A four-parameter generalised logistic (sigmoidal) function in the form: A + (B - A) / (1 + exp((xmid - x) / scal)).

method = "half_time"

A non-parametric estimate of the time to recover half of the total kinetics amplitude.

method = "peak_slope"

A non-parametric estimate of the peak dynamic rate of change in the mNIRS signal, calculated from rolling local linear regression slope within a window in units of x defined by peak_slope_span.

fit_span

A list of two-element numeric vectors in the form c(before, after) in units of time_channel, defining the window around each kinetics events (in sequence) to include in the model fitting process (default fit_span = list(c(30, 180))). fit_span will be coerced from a vector to a list, and the last vector specified will be carried forward in case more events are observed than fit_span is specified for.

end_fit_span

(Default end_fit_span = 20) A numeric value indicating the local window in units of the predictor variable x after the kinetics extreme (peak or trough) value to look for subsequent greater extremes. The kinetics model will be fit to the data up to the first local extreme with no subsequent greater extremes within a trailing window in units of x defined by end_fit_span, or the limits of the data.

x0

(Default = 0) A numeric value indicating the value of the predictor variable x representing the start of the kinetics event.

time_from_zero

A logical. TRUE (the default) will resample time_channel to start from zero at the kinetic event. FALSE will return the original numeric values of time_channel. Grouping multiple events together will always resample to zero.

display_mmss

A logical to display "mm:ss" on time axis of plots.

verbose

A logical to return (the default) or silence warnings and messages which can be used for data error checking. Abort errors will always be returned.

...

Additional arguments.

peak_slope_span

A numeric value defining the window in units of the predictor variable x for rolling local calculations used for method = "peak_slope".

align = c("center", "left", "right")

Specifies the window alignment of peak_slope_span as "center" (the default), "left", or "right". Where "left" is forward looking, and "right" is backward looking by the window peak_slope_span from the current time. Only used for method = "peak_slope".

fixed parameters

Parameters (coefficients) of the parametric models ("monoexponential" and "sigmoidal") can be defined a priori and fixed, to exclude them from the model fitting optimisation. e.g., A = 10 will define the function SSmonoexp(x, A = 10, B, TD, tau).

y

A numeric vector of the response variable, or the name of the variable within a dataframe.

x

An optional numeric vector of the predictor variable, or the name of the variable within a dataframe. Defaults to using the index of x = seq_along(y).

Value

A list of class mnirs.kinetics with components:

method

The kinetics method used.

model

The model object.

equation

The equation of the kinetics model used.

data

A dataframe of original and fitted model data.

fitted

A vector of fitted values returned by the model.

residuals

A vector of residuals between original and fitted values returned by the model.

x0

The value of the predictor variable indicating the start of kinetics.

coefs

A dataframe of model coefficients, including manually fixed parameters.

diagnostics

A dataframe of model goodness-of-fit metrics (AIC, BIC, R2adj RMSE, MAE, MAPE).

call

The model call.

Details

method %in% c("monoexponential", "sigmoidal") use nls() for nonlinear (weighted) least-squares estimates.

Model Parameterisation

"monoexponential": A four-parameter monoexponential association function in the form: ifelse(x <= TD, A, A + (B - A) * (1 - exp((TD - x) / tau))).

• A: text

• B: text

• TD: text

• tau: text

• MRT: text

"sigmoidal": A four-parameter generalised logistic (sigmoidal) function in the form: A + (B - A) / (1 + exp((xmid - x) / scal)).

• A: text

• B: text

• xmid: text

• scal: text

"half_time": A non-parametric estimate of the time to recover half of the total kinetics amplitude.

• A: text

• B: text

• half_value: text

"peak_slope": A non-parametric estimate of the peak dynamic rate of change in the mNIRS signal, calculated from rolling local linear regression slope within a window in units of x defined by peak_slope_span.

• A: text

• B: text

• x_fitted: text

• peak_slope: text

See also

stats::nls(), stats::SSasymp(), stats::SSfpl(),

Examples

set.seed(13)
x1 <- seq(-10, 60, by = 2)
A <- 10; B <- 100; TD <- 5; tau <- 12
y1 <- monoexponential(x1, A, B, TD, tau) + rnorm(length(x1), 0, 3)

## monoexponential kinetics ===============================
model <- process_kinetics(y1, x1, method = "monoexponential")
model
#> 
#> Monoexponential Association Nonlinear nls() Regression
#>   model:        y ~ SSmonoexp(x, A, B, TD, tau)
#>   data:         data.frame(x = x, y = y)
#>   equation:     y ~ A + (B - A) * (1 - exp((TD - x) / tau))
#> 
#>   Model Coefficients:
#>      A    B    TD   tau   MRT MRT_fitted
#>  11.97 98.9 5.436 11.65 17.08      66.92
#> 
#>   Model Diagnostics:
#>  AIC   BIC RMSE   MAE  MAPE
#>  181 188.9  2.6 2.046 5.465
#> 

if (FALSE) { # \dontrun{
## require(ggplot2)
plot(model, plot_coefs = TRUE, plot_diagnostics = TRUE, plot_residuals = TRUE)
} # }

## sigmoidal kinetics ===============================
model <- process_kinetics(y1, x1, method = "sigmoidal")
model
#> 
#> Sigmoidal Generalised Logistic Nonlinear nls() Regression
#>   model:        y ~ SSfpl(x, A, B, xmid, scal)
#>   data:         data.frame(x = x, y = y)
#>   equation:     y ~ A + (B - A) / (1 + exp((xmid - x) / scal))
#> 
#>   Model Coefficients:
#>      A     B  scal  xmid xmid_fitted
#>  7.929 95.47 5.178 13.92        51.7
#> 
#>   Model Diagnostics:
#>    AIC   BIC  RMSE   MAE  MAPE
#>  199.9 207.8 3.381 2.661 9.035
#> 

if (FALSE) { # \dontrun{
## require(ggplot2)
plot(model, plot_coefs = TRUE, plot_diagnostics = TRUE, plot_residuals = TRUE)
} # }

## half recovery time ===============================
model <- process_kinetics(y1, x1, method = "half_time")
model
#> 
#> Non-parametric Half Response Time
#>   data:         data.frame(x = x, y = y)
#>   equation:     half_x ~ interp_x(y = A + (B - A) / 2)
#> 
#>   Model Coefficients:
#>     A     B B_x half_x half_value
#>  11.9 100.8  52     16      56.32
#> 

if (FALSE) { # \dontrun{
## require(ggplot2)
plot(model, plot_coefs = TRUE, plot_diagnostics = TRUE, plot_residuals = TRUE)
} # }

## peak slope ===============================
model <- process_kinetics(y1, x1, method = "peak_slope", peak_slope_span = 10)
model
#> 
#> Peak Linear Regression Rolling Slope
#>   data:         data.frame(x = x, y = y)
#>   equation:     peak_slope ~ max(rolling_slope(y ~ x, peak_slope_span = 10))
#> 
#>   Model Coefficients:
#>  peak_slope_x peak_slope_fitted peak_slope
#>             8             29.72       5.11
#> 

if (FALSE) { # \dontrun{
## require(ggplot2)
plot(model, plot_coefs = TRUE, plot_diagnostics = TRUE, plot_residuals = TRUE)
} # }