Fit Kinetics

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,
  x0 = 0,
  window = 30,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  verbose = TRUE,
  ...
)

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

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

# S3 method for class 'peak_slope'
process_kinetics(
  y,
  x = NULL,
  data = NULL,
  x0 = 0,
  window = 30,
  method = c("monoexponential", "sigmoidal", "half_time", "peak_slope"),
  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.

x0

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

window

(Default = 30) A numeric scalar 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 the lesser of either the window or the limits of the data.

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 reoxygenation amplitude.
method = "peak_slope": A non-parametric estimate of the time to reach the peak rolling linear regression slope within a window defined by width.

verbose

A logical. TRUE (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.

width: A numeric scalar defining the window width (in units of the predictor variable x) for rolling slope calculations (only used for method = "peak_slope").
align = c("center", "left", "right"): Specifies the window alignment of width as "center" (the default), "left", or "right". Where "left" is forward looking, and "right" is backward looking by the window width from the current observation (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 L$...:

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 reoxygenation amplitude.

A: text
B: text
half_value: text

"peak_slope": A non-parametric estimate of the time to reach the peak rolling linear regression slope within a window defined by width.

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
#> 1 11.973 98.899 5.436 11.646 17.082     66.921
#> 
#>   Model Diagnostics:
#>       AIC     BIC R2adj  RMSE   MAE  MAPE
#> 1 180.958 188.875 0.994 2.600 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
#> 1 7.929 95.468 5.178 13.919      51.699
#> 
#>   Model Diagnostics:
#>       AIC     BIC R2adj  RMSE   MAE  MAPE
#> 1 199.881 207.798 0.990 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
#> 1 11.663 100.752 16.000     56.207
#> 

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", width = 10)
model
#> 
#> Peak Linear Regression Rolling Slope
#>   data:         data.frame(x = x, y = y)
#>   equation:     peak_slope ~ max(rolling_slope(y ~ x, width = 10))
#> 
#>   Model Coefficients:
#>   peak_slope_x peak_slope_fitted peak_slope
#> 1        8.000            33.223      5.110
#> 

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

Usage

Arguments

Value

Details

See also

Examples