Skip to contents

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

Usage

# 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

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).

data

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

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.

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.

verbose

A logical. TRUE (the default) will return warnings and messages which can be used for data error checking. FALSE will silence these messages. 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).

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

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 half_x half_value
#>  11.9 100.8     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             33.22       5.11
#> 

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