Package {mnirs}

Title:	Muscle Near-Infrared Spectroscopy Processing and Analysis
Version:	0.8.0
Description:	Read, process, and analyse data from muscle near-infrared spectroscopy (mNIRS) devices. Import raw data from file and return time-series data and metadata. Standardised methods for cleaning, filtering, transforming, and analysing mNIRS data. Custom plot theme and colour palette. Intended for mNIRS researchers and practitioners in exercise physiology, sports science, and clinical practice.
License:	MIT + file LICENSE
URL:	https://jemarnold.github.io/mnirs/, https://github.com/jemarnold/mnirs
BugReports:	https://github.com/jemarnold/mnirs/issues
Depends:	R (≥ 4.1)
Imports:	cli, data.table, lifecycle, readxl, rlang, stats, tibble, tidyselect, utils
Suggests:	dplyr, ggplot2, knitr, quarto, scales, signal, testthat (≥ 3.0.0), zoo
VignetteBuilder:	quarto
Config/Needs/website:	quarto, tidyverse
Config/roxygen2/version:	8.1.0
Config/testthat/edition:	3
Encoding:	UTF-8
NeedsCompilation:	no
Packaged:	2026-09-12 23:37:37 UTC; Jem
Author:	Jem Arnold [aut, cre, cph]
Maintainer:	Jem Arnold <jem.arnold@gmail.com>
Repository:	CRAN
Date/Publication:	2026-09-13 00:00:02 UTC

mnirs: Muscle Near-Infrared Spectroscopy Processing and Analysis

Description

Read, process, and analyse data from muscle near-infrared spectroscopy (mNIRS) devices. Import raw data from file and return time-series data and metadata. Standardised methods for cleaning, filtering, transforming, and analysing mNIRS data. Custom plot theme and colour palette. Intended for mNIRS researchers and practitioners in exercise physiology, sports science, and clinical practice.

Author(s)

Maintainer: Jem Arnold jem.arnold@gmail.com (ORCID) [copyright holder]

Authors:

• Jem Arnold jem.arnold@gmail.com (ORCID) [copyright holder]

See Also

Useful links:

• https://jemarnold.github.io/mnirs/

• https://github.com/jemarnold/mnirs

• Report bugs at https://github.com/jemarnold/mnirs/issues

Self-starting biexponential model

Description

Creates initial coefficient estimates for a selfStart wrapper around biexponential(), for use with stats::nls(). Supports both the 5-parameter (A, B, tau, B2, tau2) and 6-parameter forms adding a time delay TD; arity is inferred from the formula passed to stats::nls().

Usage

SSbiexponential(t, A, B, tau, B2, tau2, TD)

Arguments

Details

Model formulas

• 5-parameter: x ~ SSbiexponential(t, A, B, tau, B2, tau2)

• 6-parameter: x ~ SSbiexponential(t, A, B, tau, B2, tau2, TD)

The two phases are weakly identified when tau and tau2 are close, so algorithm = "port" with the time constants bounded non-negative and control = nls.control(warnOnly = TRUE) is recommended. analyse_kinetics() instead fits the phases sequentially, holding the fast phase near a monoexponential estimate.

The 5-parameter form is recommended for small samples or when no obvious time delay is expected, as it converges more reliably. stats::nls() reads the free parameters from the formula right-hand side, so omitting TD incurs no degrees-of-freedom penalty.

Starting estimates are profiled on a coarse grid of tau, tau2 (and TD) with the amplitudes solved by least squares at each grid point, keeping the residual-minimising start. Grid pairs with tau / tau2 > 0.98 are dropped as near-collinear.

The model function returns the analytic gradient for the free parameters as a "gradient" attribute, so stats::nls() does not resort to stats::numericDeriv(). stats::predict() on a fitted model carries the attribute; drop it with as.vector().

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSbiexponential(t, A, B, tau = 5, B2, tau2) holds the fast time constant at 5. Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

biexponential(), analyse_kinetics(), stats::nls(), stats::selfStart(), SSmonoexponential(), SSexponential_drift()

Examples

## create a biexponential excursion-recovery curve with random noise
set.seed(13)
t <- 0:120
x <- biexponential(t, A = 70, B = 40, tau = 5, B2 = 60, tau2 = 40) +
    rnorm(length(t), 0, 0.8)
data <- data.frame(t, x)

## 5-parameter fit
model <- nls(
    x ~ SSbiexponential(t, A, B, tau, B2, tau2),
    data = data,
    algorithm = "port",
    lower = c(-Inf, -Inf, 0, -Inf, 0),
    control = nls.control(warnOnly = TRUE)
)
summary(model)

## fix the fast time constant `tau` at a known value
model_fixed <- nls(
    x ~ SSbiexponential(t, A, B, tau = 5, B2, tau2),
    data = data,
    algorithm = "port",
    lower = c(-Inf, -Inf, -Inf, 0),
    control = nls.control(warnOnly = TRUE)
)
summary(model_fixed)

Self-starting exponential-drift model

Description

Creates initial coefficient estimates for a selfStart wrapper around exponential_drift(), for use with stats::nls(). Supports both the 5-parameter (A, B, tau, slope_B, drift_fraction) and 6-parameter forms adding a time delay TD; arity is inferred from the formula passed to stats::nls().

Usage

SSexponential_drift(t, A, B, tau, slope_B, drift_fraction, TD)

Arguments

Details

Model formulas

• 5-parameter: x ~ SSexponential_drift(t, A, B, tau, slope_B, drift_fraction)

• 6-parameter: x ~ SSexponential_drift(t, A, B, tau, slope_B, drift_fraction, TD)

The hinge at the drift onset TD - tau * log(1 - drift_fraction) is not differentiable, so algorithm = "port" with tau (and TD) bounded non-negative and control = nls.control(warnOnly = TRUE) is recommended.

Starting estimates are profiled on a coarse grid of tau (and TD) with A, B, and slope_B solved by least squares at each grid point, keeping the residual-minimising start.

The model function returns the analytic gradient (one-sided at the hinge) for the free parameters as a "gradient" attribute, so stats::nls() does not resort to stats::numericDeriv(). stats::predict() on a fitted model carries the attribute; drop it with as.vector().

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSexponential_drift(t, A, B, tau, slope_B, drift_fraction = 0.95) holds the drift onset at 95% of the amplitude (TD + 3 * tau). Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

exponential_drift(), analyse_kinetics(), stats::nls(), stats::selfStart(), SSmonoexponential(), SSbiexponential()

Examples

## create an exponential curve with late linear drift and random noise
set.seed(13)
t <- 1:180
x <- exponential_drift(
    t, A = 10, B = 100, tau = 12,
    slope_B = -0.5, drift_fraction = 0.98, TD = 15
) + rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## 6-parameter fit with the drift onset held at 98% of the amplitude
model <- nls(
    x ~ SSexponential_drift(
        t, A, B, tau, slope_B, drift_fraction = 0.98, TD
    ),
    data = data,
    algorithm = "port",
    lower = c(-Inf, -Inf, 0, -Inf, 0),
    control = nls.control(warnOnly = TRUE)
)
summary(model)

Self-starting Gompertz models

Description

Creates initial coefficient estimates for selfStart wrappers around gompertz() and gompertz_left(), for use with stats::nls(). Both wrappers use the same 4-parameter (A, B, xmid, slope) interface.

Usage

SSgompertz(t, A, B, xmid, slope)

SSgompertz_left(t, A, B, xmid, slope)

SSgompertz_left(t, A, B, xmid, slope)

Arguments

Details

Model formulas

• Right-Gompertz: x ~ SSgompertz(t, A, B, xmid, slope)

• Left-Gompertz: x ~ SSgompertz_left(t, A, B, xmid, slope)

Used by analyse_kinetics() with method = "sigmoidal" and shape = "gompertz" or "gompertz_left". Starting estimates locate the inflection from a smoothed first derivative. SSgompertz() masks stats::SSgompertz().

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSgompertz(t, A = 0, B, xmid, slope) fixes the starting asymptote at A = 0. Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

gompertz(), gompertz_left(), analyse_kinetics(), SSlogistic(), stats::nls(), stats::selfStart(), stats::SSgompertz()

Examples

## create a Gompertz curve with random noise
set.seed(15)
t <- 1:60
x <- gompertz(t, A = 10, B = 100, xmid = 30, slope = 4) +
    rnorm(length(t), 0, 2)
data <- data.frame(t, x)

model <- nls(x ~ SSgompertz(t, A, B, xmid, slope), data = data)
summary(model)

## fix the starting asymptote `A` at a known value
model_fixed <- nls(x ~ SSgompertz(t, A = 10, B, xmid, slope), data = data)
summary(model_fixed)

## left-Gompertz
set.seed(16)
x2 <- gompertz_left(t, A = 10, B = 100, xmid = 30, slope = 4) +
    rnorm(length(t), 0, 2)
data2 <- data.frame(t, x = x2)

model_left <- nls(x ~ SSgompertz_left(t, A, B, xmid, slope), data = data2)
summary(model_left)

Self-starting logistic model

Description

Creates initial coefficient estimates for a selfStart wrapper around logistic(), for use with stats::nls(). Supports both the 4-parameter symmetric (A, B, xmid, slope) and 5-parameter asymmetric (A, B, xmid, slope, asym) forms; arity is inferred from the formula passed to stats::nls().

Usage

SSlogistic(t, A, B, xmid, slope, asym)

Arguments

Details

Model formulas

• 4-parameter: x ~ SSlogistic(t, A, B, xmid, slope)

• 5-parameter: x ~ SSlogistic(t, A, B, xmid, slope, asym)

The 4-parameter form is used by analyse_kinetics() with method = "sigmoidal" and shape = "symmetric". The 5-parameter asymmetric form is retained for advanced/experimental use only; analyse_kinetics() instead dispatches to SSgompertz() / SSgompertz_left() for asymmetric shapes, which are more stable. stats::nls() reads the free parameters from the formula right-hand side, so omitting asym incurs no degrees-of-freedom penalty.

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSlogistic(t, A = 0, B, xmid, slope) fixes the starting asymptote at A = 0. Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

logistic(), analyse_kinetics(), stats::nls(), stats::selfStart(), stats::SSfpl(), SSgompertz()

Examples

## create an asymmetric logistic curve with random noise
set.seed(15)
t <- 1:60
x <- logistic(t, A = 10, B = 100, xmid = 30, slope = 4, asym = 0.3) +
    rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## 4-parameter fit
model4 <- nls(x ~ SSlogistic(t, A, B, xmid, slope), data = data)
summary(model4)

## 5-parameter fit on the same data
model5 <- nls(x ~ SSlogistic(t, A, B, xmid, slope, asym), data = data)
summary(model5)

## fix the starting asymptote `A` at a known value
model_fixed <- nls(x ~ SSlogistic(t, A = 10, B, xmid, slope), data = data)
summary(model_fixed)

y4 <- predict(model4, data)
y5 <- predict(model5, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y5, colour = "5-param")) +
            ggplot2::geom_line(ggplot2::aes(y = y4, colour = "4-param"))
    }

Self-starting monoexponential model

Description

Creates initial coefficient estimates for a selfStart wrapper around monoexponential(), for use with stats::nls(). Supports both the 3-parameter (A, B, tau) and 4-parameter (A, B, tau, TD) forms; arity is inferred from the formula passed to stats::nls().

Usage

SSmonoexponential(t, A, B, tau, TD)

Arguments

Details

Model formulas

• 3-parameter: x ~ SSmonoexponential(t, A, B, tau)

• 4-parameter: x ~ SSmonoexponential(t, A, B, tau, TD)

The 3-parameter form is recommended for small samples or when no obvious time delay is expected, as it converges more reliably. stats::nls() reads the free parameters from the formula right-hand side, so omitting TD incurs no degrees-of-freedom penalty.

Starting estimates are profiled on a coarse grid of tau (and TD) with the asymptotes solved by least squares at each grid point, keeping the residual-minimising start.

The model function returns the analytic gradient for the free parameters as a "gradient" attribute, so stats::nls() does not resort to stats::numericDeriv(). stats::predict() on a fitted model carries the attribute; drop it with as.vector().

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSmonoexponential(t, A = 0, B, tau) fixes the baseline at A = 0. Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

monoexponential(), analyse_kinetics(), stats::nls(), stats::selfStart(), stats::SSasymp()

Examples

## create an exponential curve with random noise
set.seed(13)
t <- 1:60
x <- monoexponential(t, A = 10, B = 100, tau = 8, TD = 15) +
    rnorm(length(t), 0, 3)
data <- data.frame(t, x)

## 4-parameter fit
model4 <- nls(x ~ SSmonoexponential(t, A, B, tau, TD), data = data)
summary(model4)

## 3-parameter fit on the same data
model3 <- nls(x ~ SSmonoexponential(t, A, B, tau), data = data)
summary(model3)

## fix the baseline `A` at a known value
model_fixed <- nls(x ~ SSmonoexponential(t, A = 10, B, tau, TD), data = data)
summary(model_fixed)

y4 <- predict(model4, data)
y3 <- predict(model3, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y4, colour = "4-param")) +
            ggplot2::geom_line(ggplot2::aes(y = y3, colour = "3-param"))
    }

Self-starting sigmoidal-drift model

Description

Creates initial coefficient estimates for a selfStart wrapper around sigmoidal_drift(), for use with stats::nls(): a 4-parameter sigmoid (A, B, xmid, slope) with a linear drift slope_B at its ending asymptote from the onset fraction drift_fraction.

Usage

SSsigmoidal_drift(t, A, B, xmid, slope, slope_B, drift_fraction, shape)

Arguments

Details

Model formula

x ~ SSsigmoidal_drift(t, A, B, xmid, slope, slope_B, drift_fraction = 0.95, shape = "gompertz")

drift_fraction should be written as a constant, and shape is a string constant ("symmetric" when omitted); neither is estimated. The hinge at the drift onset is not differentiable, so algorithm = "port" with control = nls.control(warnOnly = TRUE) is recommended.

Starting estimates seed the sigmoid as for SSgompertz(), resolve the drift onset from that seed, and regress the residual past the onset on time to seed slope_B and correct the asymptote B.

Fixing parameters

Any parameter may be held constant by writing a value in place of its name in the formula, e.g. x ~ SSsigmoidal_drift(t, A = 0, B, xmid, slope, slope_B, drift_fraction = 0.95) fixes the starting asymptote at A = 0. Fixed parameters are excluded from estimation and are not returned by stats::coef().

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

sigmoidal_drift(), analyse_kinetics(), stats::nls(), stats::selfStart(), SSlogistic(), SSgompertz(), SSexponential_drift()

Examples

## create a Gompertz curve with late linear drift and random noise
set.seed(13)
t <- 1:120
x <- sigmoidal_drift(
    t, A = 10, B = 100, xmid = 40, slope = 4,
    slope_B = -0.4, drift_fraction = 0.95, shape = "gompertz"
) + rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## fit with the drift onset held at 95% of the amplitude
model <- nls(
    x ~ SSsigmoidal_drift(
        t, A, B, xmid, slope, slope_B,
        drift_fraction = 0.95, shape = "gompertz"
    ),
    data = data,
    algorithm = "port",
    control = nls.control(warnOnly = TRUE)
)
summary(model)

validate_numeric abort message construction

Description

validate_numeric abort message construction

Usage

abort_validation(
  name,
  integer = FALSE,
  msg1 = "",
  msg2 = "",
  env = rlang::caller_env()
)

Accept or reject a non-converged port fit

Description

stats::nls() with algorithm = "port" and warnOnly = TRUE returns a model whose stop certificate failed. It is kept with a warning when ok holds and its coefficients are finite; otherwise it is reported as an error and dropped. The port stop code is reported in prose either way.

Usage

accept_port_fit(model, on_error, ok = TRUE)

Arguments

Value

model or NULL.

Analyse biexponential kinetics across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "biexponential"). Fits a biexponential excursion-recovery curve to each nirs_channel within a single "mnirs" data frame via fit_biexponential(), falling back down the chain in kinetics_fallbacks where the phases are unsupported. See analyse_kinetics() for user-facing documentation.

Usage

analyse_biexponential(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  use_TD = TRUE,
  fix = NULL,
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  tau_flex = 1/3,
  TD_flex = 2,
  A_flex = NULL,
  control = NULL,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, model, A, B, TD, tau, MRT, texc, B2, tau2, MRT_fitted, texc_fitted, plus the columns of the fallback models. Per-channel metadata are attached as attributes:

• "model": an nls model object, or NULL for channels where fitting failed.

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

• "warnings": a data.frame of conditions captured during fitting.

See Also

analyse_kinetics(), biexponential(), SSbiexponential()

Analyse exponential-drift kinetics across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "exponential_drift"). Fits a two-phase monoexponential + linear-drift curve to each nirs_channel within a single "mnirs" data frame via fit_exponential_drift(), falling back to fit_monoexponential() where the drift is unsupported (see kinetics_fallbacks). See analyse_kinetics() for user-facing documentation.

Usage

analyse_exponential_drift(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  use_TD = TRUE,
  drift_fraction = 0.95,
  fix = NULL,
  control = NULL,
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, model, A, B, TD, tau, k, MRT, HRT, texc, slope_B, drift_fraction, MRT_fitted, HRT_fitted, texc_fitted. texc is the excursion point where the drift rate overtakes the decaying primary rate, never before the drift onset (see expdrift_onset()). Per-channel metadata are attached as attributes:

• "model": an nls model object, or NULL for channels where fitting failed.

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

• "warnings": a data.frame of conditions captured during fitting.

See Also

analyse_kinetics(), exponential_drift(), SSexponential_drift()

Analyse kinetics across mNIRS channels and intervals

Description

Fit oxygenation kinetics (response time course) models with various parametric and non-parametric methods.

Usage

analyse_kinetics(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  method = c("response_time", "peak_slope", "monoexponential", "exponential_drift",
    "biexponential", "sigmoidal", "sigmoidal_drift"),
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  group_intervals = "ensemble",
  zero_time = FALSE,
  verbose = TRUE,
  ...,
  response_fraction = 0.5,
  width = NULL,
  span = NULL,
  align = c("centre", "left", "right"),
  partial = FALSE,
  na.rm = FALSE,
  use_TD = TRUE,
  shape = c("symmetric", "gompertz", "gompertz_left"),
  drift_fraction = NULL,
  fix = NULL
)

analyze_kinetics(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  method = c("response_time", "peak_slope", "monoexponential", "exponential_drift",
    "biexponential", "sigmoidal", "sigmoidal_drift"),
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  group_intervals = "ensemble",
  zero_time = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

Data input formats

analyse_kinetics() accepts data in multiple formats:

• A single "mnirs" data frame is processed as a single interval.

• A list of "mnirs" data frames: each interval is processed separately.

• A grouped "mnirs" data frame, e.g. with dplyr::group_by(): the data frame is split by grouping levels and each group is processed as a separate interval.

• A special case for recursive analysis: The results from analyse_kinetics() can be fed into a second call to analyse the results$coefficients table, split into data frames by nirs_channel with one row per interval (see Recursive analysis).

Specified nirs_channels (or channels retrieved from "mnirs" metadata) will be analysed and results returned as a formatted table.

Response start_time and the baseline window

start_time should be specified as the time point separating the pre-response baseline (time_channel <= start_time) from the start of the systematic response fit window (time_channel > start_time). This often corresponds to a stimulus or start/end of an intervention (e.g. start/end of an exercise interval).

For intervals extracted with extract_intervals(), start_time can be retrieved from "mnirs" metadata. Otherwise start_time defaults to 0 or the first positive time_channel value.

All methods are fitted on time elapsed from start_time, so returned time & duration coefficients are relative to response onset (e.g. start_time = 0).

• For "response_time", the baseline window before start_time defines the mean starting amplitude A directly and anchors the start of the response_time parameter.

• For "peak_slope", start_time anchors the start of the peak_slope_time parameter.

• For "exponential"- and "sigmoidal"-family, the baseline window before start_time anchors the starting fitted amplitude A and the start of TD and MRT, or xmid parameters. (see respective method sections below).

The time-delay models ("exponential"-family with use_TD = TRUE) are flat at A before TD, so the pre-onset baseline is included in the fit and anchors A. Their reduced forms (use_TD = FALSE, or a TD fit that failed and fell back) have no such flat region and are fitted only where time_channel >= start_time.

Response direction and the fit end_window

direction is detected automatically by default as either "positive" (upward) or "negative" (downward) response, and can be overwritten manually. end_window is a time span in units of time_channel defining the end of the kinetics fitting window by locating the first extrema (peak/trough, depending on direction) with no greater/lesser values within the subsequent end_window time span. The curve fitting window extends to the end of end_window beyond the detected extrema.

For "exponential"- and "sigmoidal"-family methods, direction also constrains the sign of the fitted amplitude B - A, and the sigmoidal slope. For the "biexponential" method, direction constrains the sign of the fast-phase amplitude B - A. A fit that cannot satisfy the requested direction returns NA coefficients with a warning.

Grouping samples with group_intervals

group_intervals = "ensemble" (the default) analyses every sample of each data frame together as one interval. A list() of sample (row) numbers instead splits each data frame into one interval per group, e.g. for a 20-row data frame:

analyse_kinetics(
    data,
    method = "monoexponential",
    group_intervals = list(trial1 = 1:10, trial2 = 11:20)
)

• List names become interval names; unnamed groups are ⁠interval_<n>⁠.

• Interval names are suffixed ⁠<group>_<df>⁠ (e.g. trial1_A).

• For "mnirs_kinetics" results analysed recursively, the source nirs_channel is prefixed to the analysed coefficient names (e.g. smo2_slope)

• Samples in no group are excluded from analysis (with a message). Samples in more than one group are allowed (with a warning).

• Row-grouped intervals no longer correspond to their extract_intervals() interval_times metadata, which is dropped, so start_time falls back to the first non-negative time_channel value unless supplied explicitly (optionally per-interval, keyed by group name).

• zero_time = TRUE rebases each group's time_channel to its first sample, so start_time then defaults to 0.

• Per-interval arguments key by the group names (see below).

Per-channel and per-interval arguments

Arguments apply globally to all nirs_channels by default. Arguments can instead be uniquely supplied per-channel as a named list() with names matching nirs_channels. For multi-interval input (a list of data frames or a grouped data frame), a named list() can also be keyed by interval name (the list names, group keys, or ⁠interval_<n>⁠) to supply values per-interval, and each per-interval value may itself be a per-channel list(), e.g.

analyse_kinetics(
    data,
    nirs_channels = c(o2hb, hhb),
    method = "peak_slope",
    span = list(10, o2hb = 20),
    direction = list(
        interval_1 = list(hhb = "negative", "auto"),
        interval_2 = "positive"
    )
)

The same rules apply at both levels:

• A non-list value applies to every interval and channel (the default behaviour).

• A list() named by interval or nirs_channels applies to those values per-interval or per-channel.

• A single unnamed value in the list is the fallback applied to any unlisted intervals or channels (e.g. span = list(10, o2hb = 20) gives o2hb 20 and every other channel 10). If no unnamed fallback value in the list, unlisted intervals or channels fall back to the argument's default (i.e. NULL, or may fail with a warning).

• list() names matching neither interval names nor nirs_channels are warned about and ignored.

start_time, direction, and end_window are per-channel and per-interval capable, along with the method-specific arguments except control, which is always global. fix is itself a named list() of model parameters, so a per-channel or per-interval fix is supplied as a list() of list()s keyed by channel or interval name. A plain parameter list applies everywhere:

## fix `A` at 0 for every channel
fix = list(A = 0)

## fix `A` per-channel, leaving unspecified channels free
fix = list(o2hb = list(A = 0), hhb = list(A = 5, B = 20))

## fix `A` per-interval, optionally nested per-channel
fix = list(interval_1 = list(A = 0))
fix = list(interval_1 = list(o2hb = list(A = 0)))

Triple nested list()s is janky, but it works for now!

Limitation: method itself currently only accepts a single value applied globally to all intervals and nirs_channels. So analysing channels or intervals with entirely different kinetics models must be done with independent analyse_kinetics() calls, or other iterative solutions (e.g. lapply() or purrr::map()).

method = "response_time"

Aliases: method = c("response time", "half recovery time", "half time", "HRT").

A non-parametric approach (estimated directly from the observed data without assuming a specific mathematical shape) to estimate the response time at which a signal reaches a specified fraction of its total response amplitude relative to the baseline. e.g. half-response time (response_fraction = 0.5) is the time from response onset to attain 50% of the total amplitude change and approximates the inflection point (xmid of a symmetrical sigmoid function).

response_fraction = 0.632 approximates the time constant (tau; \tau) parameter from a monoexponential function, or the inflection point (xmid) of an asymmetrical left-Gompertz function. response_fraction = 0.368 approximates xmid of a right-Gompertz function. This is a good fallback estimation method if parametric methods are not successfully fit.

The target response value is: fitted = A + (B - A) * response_fraction

Where A is the mean baseline value (time_channel <= start_time) and B is the first local extreme (peak or trough) value with no greater extreme values within end_window. response_value is the first observed sample where the signal is equal to or greater/lesser than the target response_fitted value. response_time is the elapsed time from start_time to response_value. See response_time() for the full algorithm and coefficients.

method = "peak_slope"

Aliases: method = c("peak slope", "slope", "lm").

A semi-parametric approach to estimate the maximum positive or negative local linear slope of a signal using rolling least-squares regression. The steepest local rate of change in NIRS signals can be interpreted as the moment of greatest mismatch between oxygen delivery and extraction. peak_slope_time is the time from response onset start_time to this moment of greatest mismatch.

The local window is defined by either width (number of samples) or span (in units of time_channel). See peak_slope() for window mechanics, partial-window behaviour, and the returned vector-level list.

method = "monoexponential"

Aliases: method = c("monoexp", "exponential", "exp", "tau", "MRT").

A parametric approach fitting a self-starting monoexponential function to the response curve using stats::nls() with SSmonoexponential() for either a 4-parameter (A, B, tau, TD) or 3-parameter (A, B, tau) model.

Model equations:

• 3-parameter: A + (B - A) * (1 - exp(-t / tau))

• 4-parameter: A + (B - A) * (1 - exp(-pmax(t - TD, 0) / tau))

TD is the time delay from start_time to the onset of the exponential response curve. tau is the time constant of the response. The rate constant k is the reciprocal (k = 1 / tau). The mean response time is the time sum MRT = TD + tau. See monoexponential() for the model family and SSmonoexponential() for self-start initialisation.

Any parameter may be held constant with fix, e.g. fix = list(A = 0). This excludes them from the fit optimisation procedure, and effectively reduces the function to a lower-parameter model. TD can only be fixed when use_TD = TRUE and disables the 3-parameter fallback. It is recommended to specify use_TD = FALSE rather than fix TD = 0.

method = "exponential_drift"

Aliases: method = c("exp_drift", "exp_linear", "monoexp_drift").

A parametric approach fitting a self-starting two-phase curve using stats::nls() with SSexponential_drift(). A fast monoexponential() primary response plus a slow linear secondary drift beginning near the primary asymptote.

Model equation:

A + (B - A) * (1 - exp(-pmax(t - TD, 0) / tau)) + slope_B * pmax(t - TD + tau * log(1 - drift_fraction), 0)

A, B, tau, TD, and the derived k, MRT, and HRT are as for "monoexponential". slope_B is the linear drift rate dx/dt. The drift onset is not a free estimate. drift_fraction specifies the fraction (⁠(0.5, 1)⁠) of the primary response amplitude where the drift begins; TD - tau * log(1 - drift_fraction) (default 0.95; TD + 3 * tau).

The excursion point texc is where the drift rate overtakes the decaying primary rate, ⁠TD + tau * log(|B - A| / (|slope_B| * tau))⁠, floored at the drift onset, elapsed from start_time (the same frame as TD and MRT).

The drift component is kept only when the data support it. The model will fall back to "monoexponential" when the fit fails or if the total drift amplitude is below twice the fit RMSE, with a warning recorded in warnings. The model column in coefficients names the final method for each row. A hidden argument model_fallback = FALSE will override the fallback process and retain the more complex model, or return an error.

Parameters may be held constant with fix, e.g. fix = list(A = 0), as above.

method = "biexponential"

Aliases: method = c("biexp", "double exponential").

A parametric approach fitting a self-starting two-phase biexponential excursion-recovery function to the response curve using stats::nls() with SSbiexponential(). A fast primary component driving the initial excursion, and a slow secondary component recovering the response toward a stable plateau.

Model equations:

• 5-parameter: A + (B - A) * (1 - exp(-t / tau)) + (B2 - B) * (1 - exp(-t / tau2))

• 6-parameter, where ts = pmax(t - TD, 0): A + (B - A) * (1 - exp(-ts / tau)) + (B2 - B) * (1 - exp(-ts / tau2))

A is the starting value. B & tau are the asymptote and time constant of the fast response. B2 & tau2 are the asymptote and time constant of the slower response plateau (typically ⁠tau2 >> tau⁠).

Set use_TD = TRUE (default) to specify the time-delay parameter TD. The fast-phase mean response time MRT = TD + tau is reported as for "monoexponential". See biexponential() for the model family and SSbiexponential() for self-start initialisation.

The two phases are fit sequentially.

• Stage 1 fits the fast phase as a "monoexponential" on the supplied end_window window, giving A, tau, and TD (if selected).

• Stage 2 fits the full model to the whole response with A, tau, and TD held within a tight range of their stage-1 values, and B, B2, tau2 free.

Secondary tau2 is floored above the primary tau, so the phases stay separated. tau2 is arbitrarily capped at ten times the fit window timespan, functionally implying the true asymptote is linear not exponential. end_window should be set to isolate the fast phase; by default resolves to 30 sec instead of Inf (recorded in channel_args).

The biexponential fit is kept only when the data support both phases. The model will fall back to "exponential_drift" when the fit fails (e.g. phases not separable), the fitted response is monotonic (no estimable excursion point texc), tau2 exceeds twice the fitted time span (a slow phase the record cannot tell from a linear drift), or the slow-phase amplitude ⁠|B2 - B|⁠ is below twice the fit RMSE.

The exponential-drift fit is in turn subject to its own fallback to "monoexponential" (see above). Each fallback is warned about and recorded in warnings. The model column in coefficients names the final method for each row. A hidden argument model_fallback = FALSE will override the fallback process and retain the more complex model, or return an error.

Parameters may be held constant with fix, e.g. fix = list(A = 0), as above.

method = "sigmoidal"

Aliases: method = c("logistic", "gompertz", "xmid").

A parametric approach fitting a self-starting 4-parameter sigmoidal function to the response curve using stats::nls() in one of three shapes.

Model equations (all 4-parameter):

• shape = "symmetric" (SSlogistic()): A + (B - A) / (1 + exp(-4 * slope * (t - xmid) / (B - A)))

• shape = "gompertz" (SSgompertz()): A + (B - A) * exp(-exp(-k * (t - xmid))) with k = slope * e / (B - A). Early-acceleration; inflection height fixed at A + (B - A) / e; 36.8% of the amplitude.

• shape = "gompertz_left" (SSgompertz_left()): A + (B - A) * (1 - exp(-exp(k * (t - xmid)))) with k = slope * e / (B - A). Late-acceleration; inflection height fixed at A + (B - A) * (1 - 1/e); 63.2% of the amplitude.

xmid is the time from start_time to the inflection point; the steepest point of the response. slope is the response rate dx/dt at the inflection.

A "symmetric" shape is the default when no obvious asymmetry is expected. "gompertz" (right-inflection) growth is appropriate for fast-onset, slow-tail responses. "gompertz_left" for slow-onset, fast-tail responses. See logistic(), gompertz(), and gompertz_left() for the model families and SSlogistic(), SSgompertz(), and SSgompertz_left() for self-start initialisations.

Parameters may be held constant with fix, e.g. fix = list(A = 0), as above.

method = "sigmoidal_drift"

Aliases: method = c("sigmoid_drift", "sig_drift", "sig-lin", "logistic_drift", "gompertz_drift").

A parametric approach fitting a self-starting two-phase curve using stats::nls() with SSsigmoidal_drift(). A fast "sigmoidal" primary response of the given shape plus a slow linear secondary drift beginning near the primary ending asymptote.

Model equation:

S(t) + slope_B * pmax(t - onset, 0)

S(t) and A, B, xmid, and slope are as for "sigmoidal". slope_B is the linear drift rate dx/dt at the asymptote B. The drift is not a free estimate. drift_fraction specifies the fraction (⁠(0.5, 1)⁠) of the primary response amplitude where the drift begins (default 0.95).

The excursion point texc is where the drift rate overtakes the decaying primary rate, ⁠|S'(t)| = |slope_B|⁠, floored at the drift onset, elapsed from start_time (the same frame as xmid).

The drift component is kept only when the data support it. The model will fall back to "sigmoidal" when the fit fails or if the total drift amplitude is below twice the fit RMSE, with a warning recorded in warnings. The model column in coefficients names the final method for each row. A hidden argument model_fallback = FALSE will override the fallback process and retain the more complex model, or return an error.

Parameters may be held constant with fix, e.g. fix = list(A = 0), as above.

Recursive analysis

An "mnirs_kinetics" result may be passed back as data to analyse how coefficients change across intervals, e.g. analyse_kinetics(result, nirs_channels = tau, time_channel = start_time, method = "peak_slope"). nirs_channels and time_channel must name coefficient columns explicitly; no metadata defaults are applied.

Time-point coefficients (response_time, peak_slope_time, TD, MRT, HRT, texc, xmid) are elapsed from each interval's start_time. When one of these is given as time_channel, start_time is added row-wise so the analysis runs on absolute time. start_time itself and duration coefficients (e.g. tau) are unchanged.

Coefficient rows from separate trials can be analysed separately with group_intervals, e.g. 20 occlusion slopes from two trials:

analyse_kinetics(
    result,
    nirs_channels = slope,
    time_channel = peak_slope_time,
    method = "monoexponential",
    group_intervals = list(trial1 = 1:10, trial2 = 11:20)
)

Value

A formatted table of results, with individual elements accessible as a structured list of class "mnirs_kinetics" containing:

See Also

extract_intervals(), response_time(), peak_slope(), monoexponential(), exponential_drift(), biexponential(), logistic(), gompertz(), gompertz_left(), sigmoidal_drift()

Examples

result <- read_mnirs(
    file_path = example_mnirs("train.red"),
    nirs_channels = c(
        smo2_left = "SmO2 unfiltered",
        smo2_right = "SmO2 unfiltered"
    ),
    time_channel = c(time = "Timestamp (seconds passed)"),
    zero_time = TRUE,
    verbose = FALSE
) |>
    resample_mnirs(method = "linear", verbose = FALSE) |>
    extract_intervals(
        group_intervals = "distinct",
        start = by_time(368, 1084),
        span = c(-20, 90),
        zero_time = TRUE,
        verbose = FALSE
    ) |>
    analyse_kinetics(
        nirs_channels = c(smo2_left, smo2_right),
        method = "peak_slope",
        span = 10,          ## 10-second rolling window
        direction = "auto", ## auto-detect slope direction
        verbose = FALSE
    )

## formatted table of results
result

## coefficients are accessible from the result list
result$coefficients

## along with diagnostics and other returned objects
result$diagnostics

## plot results
plot(result)

Process kinetics fits across NIRS channels

Description

Shared per-channel skeleton for all analyse_kinetics() methods. For each channel, resolves the fitting window via find_kinetics_idx(), delegates the method-specific fit to fit_fn, and, for methods listed in kinetics_fallbacks, tests the fit with the method's trigger. A channel with a reason is refit by the reduced method (recursively down the chain) with the arguments it takes, the spec's overrides, and the user-fixed parameters it shares; the fallback is warned about and so recorded in the warnings attribute. A row where every fit in the chain failed reports the last method tried with NA coefficients.

Usage

analyse_kinetics_channels(
  data,
  nirs_channels,
  time_channel,
  per_channel,
  fit_fn,
  verbose = TRUE,
  interval_name = NA_character_,
  extra_args = list(),
  method = NULL,
  fallback = TRUE,
  env = rlang::caller_env()
)

Arguments

Details

Methods with a fallback report the fitting method per row in a model coefficient column and the union of the chain's coefficient columns (NA where a model has no such parameter), so intervals bind regardless of which triggers fire. build_kinetics_results() then drops the columns of fallback models no row resolved to.

Value

A data.frame of coefficients (columns interval, nirs_channels, model for chained methods, and method parameters), one row per channel, with attributes "time_channel" (the resolved time column name), "model" and "fitted_data" (named lists by channel), "diagnostics" and "channel_args" (data frames, one row per channel), and "warnings" (data frame of conditions captured during fitting, regardless of verbose; zero rows when none fire).

Run a kinetics worker over each interval and collate results

Description

Shared skeleton for ⁠analyse_kinetics.*⁠ methods: normalises data to a named list of interval data frames, splits sample groups via split_kinetics_groups(), calls the method worker in kinetics_workers once per interval, and collates results via build_kinetics_results().

Usage

analyse_kinetics_intervals(
  data,
  method,
  worker_args,
  nirs_quo,
  time_quo,
  group_intervals,
  zero_time,
  verbose,
  call,
  env,
  fallback = TRUE
)

Arguments

Value

An "mnirs_kinetics" object from build_kinetics_results().

Analyse logistic kinetics across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "sigmoidal"). Fits a 4-parameter sigmoidal curve to each nirs_channel within a single "mnirs" data frame via fit_sigmoidal() with one of three shapes: "symmetric", "gompertz", or "gompertz_left". See analyse_kinetics() for user-facing documentation.

Usage

analyse_logistic(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  shape = c("symmetric", "gompertz", "gompertz_left"),
  fix = NULL,
  control = NULL,
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, A, B, xmid, slope, xmid_fitted. Per-channel metadata are attached as attributes:

• "model": an nls model object, or NULL for channels where fitting failed.

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

See Also

analyse_kinetics(), logistic(), SSlogistic(), gompertz(), gompertz_left(), SSgompertz(), SSgompertz_left()

Analyse monoexponential kinetics across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "monoexponential"). Fits a monoexponential curve to each nirs_channel within a single "mnirs" data frame via fit_monoexponential(). See analyse_kinetics() for user-facing documentation.

Usage

analyse_monoexponential(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  use_TD = TRUE,
  fix = NULL,
  control = NULL,
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, A, B, TD, tau, k, MRT, HRT, MRT_fitted, HRT_fitted. Per-channel metadata are attached as attributes:

• "model": an nls model object, or NULL for channels where fitting failed.

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

See Also

analyse_kinetics(), monoexponential(), SSmonoexponential()

Analyse peak linear slope across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "peak_slope"). Computes the maximum local linear slope for each nirs_channel within a single "mnirs" data frame. See analyse_kinetics() for user-facing documentation.

Usage

analyse_peak_slope(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  start_time = NULL,
  width = NULL,
  span = NULL,
  align = c("centre", "left", "right"),
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  partial = FALSE,
  na.rm = FALSE,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, slope, intercept, y, peak_slope_time, idx. Per-channel metadata are attached as attributes:

• "model": a linear regression model object via stats::lm().

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

See Also

analyse_kinetics(), peak_slope()

Analyse fractional kinetics response time across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "response_time"). Computes the fractional response time for each nirs_channel within a single "mnirs" data frame. See analyse_kinetics() for user-facing documentation.

Usage

analyse_response_time(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  start_time = NULL,
  response_fraction = 0.5,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel per response_fraction and columns nirs_channels, response_fraction, A, B, response_time, response_value, fitted, idx. Per-channel metadata are attached as attributes:

• "model": NULL (no parametric model is fitted).

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted, containing the baseline, response, and extreme key points.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

See Also

analyse_kinetics(), response_time()

Analyse sigmoidal-drift kinetics across NIRS channels

Description

Internal channel-level dispatch for analyse_kinetics(method = "sigmoidal_drift"). Fits a two-phase sigmoidal + linear-drift curve to each nirs_channel within a single "mnirs" data frame via fit_sigmoidal_drift(), falling back to fit_sigmoidal() where the drift is unsupported (see kinetics_fallbacks). See analyse_kinetics() for user-facing documentation.

Usage

analyse_sigmoidal_drift(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  shape = c("symmetric", "gompertz", "gompertz_left"),
  drift_fraction = 0.95,
  fix = NULL,
  control = NULL,
  start_time = NULL,
  direction = c("auto", "positive", "negative"),
  end_window = Inf,
  verbose = TRUE,
  ...,
  env = rlang::caller_env()
)

Arguments

Value

A data.frame with one row per nirs_channel and columns nirs_channels, model, A, B, xmid, slope, texc, slope_B, drift_fraction, xmid_fitted, texc_fitted. texc is the excursion point where the drift rate overtakes the decaying sigmoid rate, never before the drift onset (see sigdrift_texc()). Per-channel metadata are attached as attributes:

• "model": an nls model object, or NULL for channels where fitting failed.

• "fitted_data": a named list of per-channel data frames with columns window_idx and fitted.

• "diagnostics": a data.frame with one row per nirs_channel containing model fit diagnostics.

• "channel_args": a data.frame with one row per nirs_channel recording the resolved arguments used.

• "warnings": a data.frame of conditions captured during fitting.

See Also

analyse_kinetics(), sigmoidal_drift(), SSsigmoidal_drift()

Apply grouping to intervals

Description

Apply grouping to intervals

Usage

apply_interval_groups(
  df_list,
  group_channels,
  metadata,
  group_intervals,
  zero_time = FALSE,
  verbose = TRUE,
  env = rlang::caller_env()
)

Arguments

apply span to resolved times and build interval_spec data frame

Description

apply span to resolved times and build interval_spec data frame

Usage

apply_span(
  interval_list,
  t_vec,
  span,
  verbose = TRUE,
  env = rlang::caller_env()
)

Arguments

10 Hz Artinis Oxysoft export recorded with Oxymon MKIII

Description

Exported from Artinis Oxysoft, recorded on Oxymon MKIII at 50 Hz and exported at 10 Hz. Containing two 5-minute cycling work intervals and an ischaemic occlusion, placed on the vastus lateralis muscle site.

Format

.xlsx file with header metadata and five columns and 20919 rows:

Column 1

Sample index (divide by sample rate for seconds).

Column 2

O2Hb: oxyhaemoglobin concentration change (\muM).

Column 3

HHb: deoxyhaemoglobin concentration change (\muM).

Column 4

Event marker (character).

Column 5

Unmarked event label (character).

Channels are detected automatically from the file legend, or can be specified explicitly for read_mnirs():

• nirs_channels = c(O2Hb = 2, HHb = 3)

• time_channel = c(sample = 1)

• event_channel = c(event = 4)

• ⁠interval_times = list( ## two intervals, post-exercise occlusion start = c(158, 999, 1750), end = c(493, 1333, 1961) )⁠

Source

Artinis Medical Systems. Oxymon MKIII, exported via Oxysoft desktop software (https://artinis.com/)

See Also

read_mnirs(), example_mnirs()

Examples

example_mnirs("artinis_intervals")

Coerce data input to a named list of data frames

Description

Accepts a single or grouped data frame, a list of data frames, or an "mnirs_kinetics" object, whose coefficients are split by nirs_channels into one data frame per channel (a row per interval) for recursive analysis of coefficients.

Usage

as_data_list(data, env = rlang::caller_env())

Arguments

coerce raw values to mnirs_interval objects

Description

coerce raw values to mnirs_interval objects

Usage

as_mnirs_interval(x, arg = "start", env = rlang::caller_env())

Arguments

Validate and bind a list of mnirs data frames for plotting

Description

Validate and bind a list of mnirs data frames for plotting

Usage

as_plot_data(x, env = rlang::caller_env())

Arguments

Value

For a single-element list, that element unchanged. Otherwise a row-bound data.frame with an interval factor column, carrying attributes nirs_channels (the union across elements), time_channel, and channel_map – a named list mapping each channel to the interval names whose source element declares it, so plot.mnirs() draws each channel only in its own panels.

Biexponential model with gradient

Description

biexp_core() evaluates the curve and its partial derivatives on the canonical parameters. biexp_model() is the model function of SSbiexponential(): biexponential() plus the gradient for the parameters written as bare symbols in the call (see free_params()), so stats::nls() skips stats::numericDeriv().

Usage

biexp_core(t, A, B, tau, B2, tau2, TD = NULL)

biexp_model(t, A, B, tau, B2, tau2, TD = NULL)

Arguments

Value

biexp_core(): a list of the curve val and the partial derivatives by parameter name. biexp_model(): a numeric vector of predicted values with a "gradient" attribute when any parameter is free.

Initiate self-starting biexponential model

Description

biexp_init(): Returns initial values for the parameters in a selfStart model.

Usage

biexp_init(mCall, data, LHS, ...)

Arguments

Value

biexp_init(): Initial starting estimates for parameters in the model called by SSbiexponential().

Grid-profiled starting estimates for the biexponential model

Description

Vector-level initialiser behind biexp_init(), called directly by the kinetics worker with tau and TD held at their stage-1 values to seed the slow phase. Profiles the time constants (and TD) on a coarse grid and keeps the RSS-minimising start (cf. expdrift_start()). The model is linear in A, B, and B2 once tau, tau2, and TD are held, so those are solved by least squares at every grid point at once: the Gram entries of the bases e1, e2 - e1, 1 - e2 for every ⁠(tau, tau2)⁠ pair follow from the column products of the two exponential matrices, and solve_grid3() solves the pairs in one pass. User-fixed values narrow the grids; the amplitudes are always solved free, as this is only a seed. Pairs with tau / tau2 > 0.98 are dropped as their bases are near-collinear, unless both time constants are fixed.

Usage

biexp_start(x, t, fixed = list(), has_TD = FALSE)

Arguments

Value

A named numeric vector of starting estimates in model order.

Biexponential function

Description

Calculate a two-phase curve: a fast monoexponential primary response toward B and a slow monoexponential secondary response from B toward a stable plateau at B2, both clocked from the response onset and summed. Model family fit by analyse_kinetics() with method = "biexponential", and by stats::nls() via the self-starting wrapper SSbiexponential().

Usage

biexponential(t, A, B, tau, B2, tau2, TD = NULL)

Arguments

Details

Model equations

• 5-parameter: A + (B - A) * (1 - exp(-t / tau)) + (B2 - B) * (1 - exp(-t / tau2))

• 6-parameter, where ts = pmax(t - TD, 0): A + (B - A) * (1 - exp(-ts / tau)) + (B2 - B) * (1 - exp(-ts / tau2))

A, B, and B2 are all values on the response scale. The fast component is a monoexponential() response from A toward B with amplitude B - A; the slow component runs concurrently from the same onset with amplitude B2 - B. The curve starts at A, approaches B2 as t grows, and is smooth throughout. If B = B2, the curve reduces to a monoexponential() with time constant tau and asymptote B2.

Excursion point

The expected response is a fast excursion toward a minimum or maximum short of B, followed by a slow recovery back to a stable plateau at B2. The excursion point texc occurs where the two phase rates cancel: texc = TD + log(ratio) / (1 / tau - 1 / tau2) with ratio = -(B - A) * tau2 / ((B2 - B) * tau), which exists only when the amplitudes oppose in sign and the fast phase dominates at the onset (ratio > 1). If B is between A and B2, the response is monotonic but still two-phase.

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

analyse_kinetics(), SSbiexponential(), monoexponential(), exponential_drift()

Examples

## create a biexponential excursion-recovery curve with random noise
set.seed(1)
t <- 0:120
x <- biexponential(t, A = 70, B = 40, tau = 5, B2 = 60, tau2 = 40) +
    rnorm(length(t), 0, 0.8)
data <- data.frame(t, x)

## 5-parameter fit with the self-starting wrapper
model <- nls(
    x ~ SSbiexponential(t, A, B, tau, B2, tau2),
    data = data,
    algorithm = "port",
    lower = c(-Inf, -Inf, 0, -Inf, 0),
    control = nls.control(warnOnly = TRUE)
)
summary(model)

y <- predict(model, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y))
    }

Breaks for time span data

Description

Pretty time span breaks for plotting in units of 5, 15, 30, 60 sec, etc. Modified from scales::breaks_timespan().

Usage

breaks_timespan(unit = c("secs", "mins", "hours", "days", "weeks"), n = 5)

Arguments

Value

Returns a function for generating breaks.

Examples

x <- 0:120
y <- sin(2 * pi * x / 15) + rnorm(length(x), 0, 0.2)

ggplot2::ggplot(data.frame(x, y), ggplot2::aes(x, y)) +
    theme_mnirs() +
    ggplot2::scale_x_continuous(breaks = breaks_timespan()) +
    ggplot2::geom_line()

Assemble a fitted channel result

Description

Counterpart of build_na_results() for a successful fit: the coefs/model/fitted_data/diag list expected by analyse_kinetics_channels(), with fitted values and diagnostics derived from model on the rows in keep.

Usage

build_fit_results(
  coefs,
  model,
  x_fit,
  t_fit,
  valid,
  keep = TRUE,
  env = rlang::caller_env()
)

Arguments

Value

A named list with elements coefs, model, fitted_data, and diag.

Gather per-interval mnirs_kinetics into results structure

Description

Shared helper for ⁠analyse_kinetics.*⁠ methods. Takes a list of per-interval (per-data frame) kinetics results data frames (each carrying "fitted_data", "channel_args", and "diagnostics" attributes) and the original data_list. Interval names are taken from names(data_list). Where rows carry a model column, coefficient columns owned only by fallback models no row resolved to are dropped.

Usage

build_kinetics_results(data_list, result_list, method, call)

Arguments

Value

A named list with: method, model, coefficients, data, interval_times, diagnostics, channel_args, warnings, call.

Build a standardised NA result for a failed channel

Description

Returns the method-specific coefs/model/fitted_data/diag list expected by analyse_kinetics_channels() when a model fit fails, populated with NA/NULL values.

Usage

build_na_results(na_coefs)

Arguments

Value

A named list with elements coefs, model, fitted_data, and diag.

Build a self-start model formula with optional fixed parameters

Description

Constructs x ~ fn(t, ...) with each free parameter as a bare symbol and each fixed parameter substituted as its constant value.

Usage

build_ss_formula(fn, params, fix = list(), x, t)

Arguments

Value

A two-sided formula on x and t.

Specify interval boundaries by time, label, lap, or sample

Description

Helper functions to define interval start or end boundaries for extract_intervals().

Usage

by_time(...)

by_label(..., ignore_case = FALSE, fixed = FALSE)

by_lap(...)

by_sample(...)

Arguments

Details

These helpers can be used explicitly for arguments start/end, or raw values can be passed directly:

• Numeric -> by_time()

• Character -> by_label(),

• Explicit integer (e.g. 2L) -> by_lap().

• Use by_sample() explicitly for sample indices.

Multiple specification types can be combined for a single boundary with list() (e.g. list(by_time(30), by_label("go"))). Resolved boundary times are concatenated in the order supplied. Combined specifications must use the by_ helpers directly: raw values are ignored with a warning.

Value

An object of class "mnirs_interval" for use with the start and end arguments of extract_intervals().

Examples

## read example data
data <- read_mnirs(
    example_mnirs("train.red"),
    nirs_channels = c(
        smo2_left = "SmO2 unfiltered",
        smo2_right = "SmO2 unfiltered"
    ),
    time_channel = c(time = "Timestamp (seconds passed)"),
    event_channel = c(lap = "Lap/Event"),
    zero_time = TRUE,
    verbose = FALSE
)

## start and end by time
extract_intervals(data, start = by_time(66), end = by_time(357))

## start by lap
extract_intervals(data, start = by_lap(2, 4), span = 0)

## combine multiple specification types
extract_intervals(
    data,
    start = list(by_lap(2), by_time(400)),
    end = by_sample(1500)
)

## simulate event_channel with character label match
data$event <- NA_character_
data$event[c(1000, 1001)] <- c("start", "lap.1")
data <- create_mnirs_data(data, event_channel = "event")

## case-insensitive label match
extract_intervals(data, start = by_label("START", ignore_case = TRUE))

## literal-string label match (regex metacharacters treated as text)
extract_intervals(data, start = by_label("lap.1", fixed = TRUE))

Clean legend trace names to syntactic column names

Description

Clean legend trace names to syntactic column names

Usage

clean_channel_names(x)

Flatten a captured condition message

Description

Drops cli bullet glyphs and console-width wrapping from a condition message so the stored warnings text reads as plain sentences.

Usage

clean_cnd_message(cnd)

Arguments

Value

A single character string.

Compute model diagnostics

Description

Compute model diagnostics

Usage

compute_diagnostics(x, t, fitted, n_params = 1L, env = rlang::caller_env())

Arguments

Details

r2

Squared Pearson correlation between observed and fitted values. Equals the classic 1 - SSres / SStot for OLS linear fits (matches summary(lm)$r.squared); a bounded ⁠[0, 1]⁠ pseudo-R^2 for non-linear fits such as "monoexponential" and "sigmoidal".

adj_r2

Adjusted R^2 penalised by n_params. Appropriate for OLS linear models; interpret with caution for non-linear fits.

aic, aicc, bic

Information criteria derived from a Gaussian log-likelihood with the maximum-likelihood residual variance sigma_hat^2 = SSres / n_obs. The effective parameter count is k = n_params + 1 (the +1 accounts for the estimated residual variance). Values match stats::AIC() and stats::BIC() for lm and nls fits. aicc is the small-sample correction and is NA when n_obs - k - 1 <= 0.

Value

A 1-row data.frame with columns n_obs, n_params, r2, adj_r2, rmse, cv_rmse, snr, aic, aicc, and bic.

Computes rolling local values

Description

compute_window_bounds(): Compute the start and end indices of rolling windows along a time variable t.

window_sums(): Windowed sums by cumulative-sum differencing.

window_min_obs(): Minimum number of samples spanned by a complete window.

compute_local_mean(): Compute rolling means from window bounds.

compute_local_fun(): Compute a rolling function along x from a list of rolling sample windows.

median_no_na(): Fast median for numeric vectors. Strips NAs and replicates median.default arithmetic without S3 dispatch.

compute_col_medians(): Column medians of an NA-padded numeric matrix via a single radix sort. NAs sort last per column; medians indexed from per-column valid counts. Matches median(w, na.rm = TRUE).

compute_outliers(): Computes a vector of local medians and logicals indicating outliers of x within rolling windows defined by width or span.

compute_valid_neighbours(): Compute a list of rolling window indices along x to either side of NAs.

Usage

compute_window_bounds(
  t,
  idx = seq_along(t),
  width = NULL,
  span = NULL,
  align = c("centre", "left", "right"),
  env = rlang::caller_env()
)

window_sums(v, bounds)

window_min_obs(width, span, t, min_n = 1L, env = rlang::caller_env())

compute_local_mean(x, bounds, na.rm = FALSE, min_obs = 1L)

compute_local_fun(x, window_idx, fn, ...)

median_no_na(w)

compute_col_medians(m)

compute_outliers(
  x,
  t,
  outlier_cutoff,
  width = NULL,
  span = NULL,
  env = rlang::caller_env()
)

compute_valid_neighbours(
  x,
  t = seq_along(x),
  width = NULL,
  span = NULL,
  verbose = TRUE,
  env = rlang::caller_env()
)

Arguments

Details

The local rolling window can be specified by either width as the number of samples, or span as the time span in units of t. Specifying width is often faster than span.

align defaults to "centre" the local window around idx between ⁠[idx - floor((width-1)/2),⁠ ⁠idx + floor(width/2)]⁠ when width is specified. Even width values will bias align to "left", with the unequal sample forward of idx, effectively returning NA at the last sample index. When span is specified, the local window is between ⁠[t - span/2, t + span/2]⁠.

window_min_obs() converts span to a sample count via the estimated sample rate, less two samples to buffer irregular t at the start and end of each window.

compute_local_mean() computes all window means in O(n) via window_sums(). Values are centred first so the cumulative-sum differencing error stays around eps * sqrt(n) * sd, far below measurement resolution.

Value

compute_window_bounds(): A list() with start and end integer vectors the same length as idx, giving the inclusive window bounds at each index.

window_sums(): A numeric vector the same length as bounds$start.

window_min_obs(): An integer value.

compute_local_mean(): A numeric vector the same length as bounds$start.

compute_local_fun(): A numeric vector the same length as x.

median_no_na(): A numeric value.

compute_col_medians(): A numeric vector of length ncol(m).

compute_outliers(): A list() with vectors the same length as x for with numeric local medians and logical identifying where is_outlier.

compute_valid_neighbours(): A list the same length as the NA values in x with numeric vectors of sample indices of length width samples or span units of time t for valid values neighbouring split to either side of the invalid NAs.

Coerce column types by role: nirs numeric, event integer, others detected

Description

Coerce column types by role: nirs numeric, event integer, others detected

Usage

convert_type(data, channels, verbose = TRUE, env = rlang::caller_env())

Arguments

Correct for blood volume changes

Description

Normalises mNIRS channels for the effects of blood volume changes, following the sample-wise iterative method of Beever & Tripp et al, 2020.

Usage

correct_blood_volume(
  data,
  oxy_channel = NULL,
  deoxy_channel = NULL,
  total_channel = NULL,
  verbose = TRUE
)

Arguments

Details

Specify NIRS component channels

At least two of oxy_channel, deoxy_channel, and total_channel must be specified to calculate the blood volume correction factor. Best practice is to specify all existing channels in data. Missing channels are derived from the specified pair before the correction is applied.

• total = oxy + deoxy

• oxy = total - deoxy

• deoxy = total - oxy

Multiple channel pairs can be corrected in one call by passing equal-length vectors, with each element number forming a pair (e.g. ⁠oxy_channel = c(o2hb_1, o2hb_2), deoxy_channel = c(hhb_1, hhb_2)⁠).

NOTE: the returned data frame will ONLY include corrected values for the specified channels. Non-specified channels will remain uncorrected and will therefore no longer be comparable to corrected channels. Best practice is to specify all existing channels in data.

Compute blood volume correction

If any NIRS channels have negative values, all specified channels will be ensemble-shifted by a common offset so that all channels contain only positive values. Relative scaling across channels is preserved. This is modified from the method in Beever & Tripp et al, 2020 to properly calculate total[haem] and the blood volume correction factor beta when there are negative NIRS values.

The correction factor beta is effectively the single-channel fractional (%) oxygen saturation used to normalise oxy[haem] and deoxy[haem] relative to an adjusted invariant total[haem]. This is computed as the cumulative sum of adjusted incremental differences:

\Delta\text{O2Hb}_c = \Delta\text{O2Hb} - \beta \cdot \Delta\text{THb}

\Delta\text{HHb}_c = \Delta\text{HHb} - (1 - \beta) \cdot \Delta\text{THb}

After correction, total[haem] is zero (blood volume changes are normalised).

Value

A tibble of class "mnirs" with blood volume-corrected channels written back to the specified columns, and with metadata available with attributes(). For list or grouped data frame input, returns a named list of "mnirs" tibbles, one per interval.

Data input formats

mnirs processing functions accept data in multiple formats:

• A single "mnirs" data frame is processed and returned directly.

• A list of "mnirs" data frames: each interval is processed separately and returned as a named list.

• A grouped "mnirs" data frame, e.g. with dplyr::group_by(): the data frame is split by grouping levels and each group is processed as a separate interval, returned as a named list.

References

Beever AT, Tripp TR, Zhang J, MacInnis MJ (2020) Nirs-Derived Skeletal Muscle Oxidative Capacity Is Correlated with Aerobic Fitness and Independent of Sex. J Appl Physiol (1985). doi:10.1152/japplphysiol.00017.2020

Ryan TE, Erickson ML, Brizendine JT, et al. (2012) Noninvasive Evaluation of Skeletal Muscle Mitochondrial Capacity with near-Infrared Spectroscopy: Correcting for Blood Volume Changes. J Appl Physiol (1985). doi:10.1152/japplphysiol.00319.2012

Examples

data <- read_mnirs(
    file_path = example_mnirs("artinis"),
    nirs_channels = c(o2hb = 2, hhb = 3),
    time_channel = c(sample = 1),
    verbose = FALSE,
)

plot(data)

result <- correct_blood_volume(
    data,
    oxy_channel = "o2hb",
    deoxy_channel = "hhb", ## thb will be derived from o2hb + hhb
)

plot(result)

Count maximum decimal places across a numeric vector

Description

Returns the largest number of decimal places present in any finite, non-NA element of x. Used internally by signif_trailing() for format = "digits".

Usage

count_decimals(x)

Arguments

Value

A single non-negative integer.

Count maximum significant figures across a numeric vector

Description

Returns the largest number of significant figures present in any finite, non-NA element of x. Used internally by signif_trailing() for format = "signif".

Usage

count_sigfigs(x)

Arguments

Value

A single positive integer (minimum 1).

Create an mnirs data frame with metadata

Description

Manually add class "mnirs" and metadata to an existing data frame.

Usage

create_mnirs_data(data, ...)

Arguments

Details

Intended primarily for internal use, but can be used to inject mnirs metadata into any data frame.

Value

A tibble of class "mnirs". Metadata are stored as attributes and can be accessed with attributes(data).

Examples

data <- data.frame(
    A = 1:3,
    B = seq(10, 30, 10),
    C = seq(11, 33, 11)
)

attributes(data)

## inject metadata
nirs_data <- create_mnirs_data(
    data,
    nirs_channels = c("B", "C"),
    time_channel = "A",
    sample_rate = 1
)

attributes(nirs_data)

## rename channels and update metadata
create_mnirs_data(
    nirs_data,
    nirs_channels = c(smo2 = "B", thb = "C"),
    time_channel = c(time = "A")
)

Detect the direction of a response signal

Description

Resolves whether a signal responds upward ("positive") or downward ("negative") by comparing the excursions of x above and below its initial baseline, taken as the median of the earliest samples ordered by t. The dominant excursion captures the primary response direction even when a fast initial component partially recovers over most of the record (e.g. biexponential drop-recovery), where a net slope would misreport the trend. Used internally to disambiguate peak (maximum) from trough (minimum) detection when direction = "auto".

Usage

detect_direction(
  x,
  t = seq_along(x),
  fallback = x,
  direction = c("auto", "positive", "negative")
)

Arguments

Value

A character string: "positive" or "negative".

Detect the first dttm_opts format matching a character vector

Description

Tested on the first non-empty value only, in UTC (local time zone parsing is slow on Windows).

Usage

detect_dttm_format(x)

Value

A format string, or NULL when none match.

Report warnings for unbalanced time_channel samples

Description

Report warnings for unbalanced time_channel samples

Usage

detect_irregular_samples(
  x,
  time_channel,
  verbose = TRUE,
  env = rlang::caller_env()
)

Arguments

Detect mnirs device from file metadata

Description

Detect mnirs device from file metadata

Usage

detect_mnirs_device(data, chunk = 200L)

Detect time_channel from column names or time-formatted values

Description

Detect time_channel from column names or time-formatted values

Usage

detect_time_channel(data, verbose = TRUE, env = rlang::caller_env())

Arguments

Known channel names and detection patterns for supported mNIRS devices

Description

Per device: pattern strings which must all match one header row (fixed regex flag); default time_channel and event_channel names; extra_channels companion columns (e.g. sample index, numeric tag) returned with keep_all = TRUE.

Usage

device_patterns

Datetime format strings for POSIXct parsing

Description

Time-only format must stay first: parse_dttm() treats dttm_opts[1L] as a relative time of day and the rest as absolute date-times.

Usage

dttm_opts

Make an nls model call self-contained

Description

stats::nls() stores its call arguments as the expressions it was called with (formula, .data, start[free]), resolvable only in the fitting frame. Evaluating them in that frame and storing the values lets stats::update(), stats::profile(), and insight::get_data().

Usage

embed_fit_call(model, env = parent.frame())

Arguments

Value

model with every call argument replaced by its value.

Enforce the requested direction on a converged parametric fit

Description

Direction is the sign of the primary amplitude D = B - A (B - A for the biexponential), where the primary asymptote follows A in model parameter order. A fit is satisfied when D and every free parameter lie inside the refit box (D sign-constrained; lower/ upper for the rest, e.g. a sigmoid slope sign floor) and is returned unchanged. Otherwise the model is refit on D via nls(algorithm = "port") with D sign-bounded and its magnitude floored strictly above zero (sigmoid models divide by D), then re-expressed in the original parameterisation from that optimum so the returned model reports consistent coefficient names. A refit that fails, pins a sign-floored coefficient (degenerate flat fit), or loses the requested sign on re-expression warns and returns NULL. Parameters in fix are held constant: a fixed A or B is substituted into the amplitude reparameterisation; with both fixed the amplitude sign is predetermined and no refit is possible.

Usage

enforce_direction(
  model,
  coefs,
  fit_data,
  direction,
  amp_fn,
  fn = sub("^(SS)?", "SS", as.character(amp_fn)),
  lower = NULL,
  upper = NULL,
  floor_params = NULL,
  fix = list(),
  control = NULL,
  .nirs,
  interval_name,
  env = rlang::caller_env()
)

Arguments

Value

A named list list(model, coefs) with coefs a named numeric vector in ⁠(A, B, ...)⁠ space including fixed values, or NULL when the direction cannot be satisfied (caller returns build_na_results()).

Ensemble average multiple intervals

Description

group_channels is a character vector applied to every interval, or a list of per-interval channel vectors; channels excluded from an interval do not contribute to that channel's ensemble-mean.

Usage

ensemble_intervals(
  df_list,
  group_channels,
  metadata,
  verbose = TRUE,
  env = rlang::caller_env()
)

Arguments

Get path to mnirs example files

Description

Get path to mnirs example files

Usage

example_mnirs(file = NULL)

Arguments

Value

A file path character string for selected example files stored in this package.

Examples

## lists all files
example_mnirs()

## partial matching will error if matches multiple
try(example_mnirs("moxy"))

example_mnirs("moxy_ramp")

Initiate self-starting exponential-drift model

Description

expdrift_init(): Returns initial values for the parameters in a selfStart model.

Usage

expdrift_init(mCall, data, LHS, ...)

Arguments

Value

expdrift_init(): Initial starting estimates for parameters in the model called by SSexponential_drift().

Exponential-drift model with gradient

Description

Model function of SSexponential_drift(): exponential_drift() plus the partial derivatives for the parameters written as bare symbols in the call (see free_params()), so stats::nls() skips stats::numericDeriv(). The hinge derivatives are one-sided at the drift onset.

Usage

expdrift_model(t, A, B, tau, slope_B, drift_fraction, TD = NULL)

Arguments

Value

A numeric vector of predicted values with a "gradient" attribute when any parameter is free.

Drift onset time of the exponential-drift model

Description

The time at which a monoexponential response reaches the drift_fraction fraction of its amplitude, by the analytic inverse TD - tau * log(1 - drift_fraction) (see exponential_drift()).

Usage

expdrift_onset(tau, drift_fraction, TD = NULL)

Arguments

Value

A numeric vector of onset times, TD = 0 when NULL.

Grid-profiled starting estimates for the exponential-drift model

Description

Vector-level initialiser behind expdrift_init(), called directly by the kinetics worker on the fit window. Profiles tau (and TD) on a coarse grid and keeps the RSS-minimising start (cf. monoexp_start()). The model is linear in A, B, and slope_B once tau and TD are held, so those are solved by least squares at every grid point at once via solve_grid3(). User-fixed tau, TD, and drift_fraction narrow the grids; the linear parameters are always solved free, as this is only a seed. tau is capped so the drift onset stays inside the record; a grid point whose hinge has no support is singular and skipped.

Usage

expdrift_start(x, t, fixed = list(), has_TD = FALSE)

Arguments

Value

A named numeric vector of starting estimates in model order.

Exponential-drift function

Description

Calculate a two-phase curve: a fast monoexponential() primary response plus a slow linear secondary drift beginning near the primary asymptote. Model family fit by analyse_kinetics() with method = "exponential_drift", and by stats::nls() via the self-starting wrapper SSexponential_drift().

Usage

exponential_drift(t, A, B, tau, slope_B, drift_fraction, TD = NULL)

Arguments

Details

Model equations

• 5-parameter: A + (B - A) * (1 - exp(-t / tau)) + slope_B * pmax(t + tau * log(1 - drift_fraction), 0)

• 6-parameter: A + (B - A) * (1 - exp(-pmax(t - TD, 0) / tau)) + slope_B * pmax(t - TD + tau * log(1 - drift_fraction), 0)

A, B, tau, and TD are as for monoexponential(). The drift onset is not a free estimate: the secondary drift is exactly zero before TD - tau * log(1 - drift_fraction) (TD = 0 when absent), and drift_fraction = 0.95 places the onset at TD + 3 * tau.

The excursion point texc is where the drift rate overtakes the decaying primary rate, ⁠TD + tau * log(|B - A| / (|slope_B| * tau))⁠, floored at the drift onset.

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

analyse_kinetics(), SSexponential_drift(), monoexponential(), biexponential(), sigmoidal_drift()

Examples

## create an exponential curve with late linear drift and random noise
set.seed(13)
t <- 1:180
x <- exponential_drift(
    t, A = 10, B = 100, tau = 12,
    slope_B = -0.5, drift_fraction = 0.95, TD = 15
) + rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## the drift onset fraction is held constant in the formula
model <- nls(
    x ~ SSexponential_drift(
        t, A, B, tau, slope_B, drift_fraction = 0.95, TD
    ),
    data = data,
    algorithm = "port",
    lower = c(-Inf, -Inf, 0, -Inf, 0),
    control = nls.control(warnOnly = TRUE)
)
summary(model)

y <- predict(model, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y))
    }

Extract interval data by time range

Description

Extract interval data by time range

Usage

extract_df_list(data, t_vec, interval_spec, group_channels)

Extract intervals from mnirs data

Description

Extract intervals from "mnirs" time series data, specifying interval start and end boundaries by time value, event label, lap number, or sample index.

Usage

extract_intervals(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  event_channel = NULL,
  sample_rate = NULL,
  group_intervals = c("distinct", "ensemble"),
  group_channels = NULL,
  start = NULL,
  end = NULL,
  span = list(c(-60, 60)),
  zero_time = FALSE,
  verbose = TRUE,
  event_groups = deprecated()
)

Arguments

Details

Interval specification

Interval start and end boundaries are specified using helper functions, or by passing raw values directly:

by_time()

Time values in units of time_channel.

by_label()

Strings to match in event_channel. All matching occurrences are returned.

by_lap()

Lap numbers to match in event_channel. Resolves to the first sample of each lap for start, and the last lap sample for end

by_sample()

Integer sample indices (row numbers).

Raw values supplied to start/end are auto-coerced:

• Numeric -> by_time()

• Character -> by_label(),

• Explicit integer (e.g. 2L) -> by_lap().

• Use by_sample() explicitly for sample indices.

start and end can use different specification types (e.g., start by label, end by time). When lengths differ, the shorter is recycled.

Multiple specification types can be combined for a single boundary with list() (e.g. start = list(by_time(30), by_label("go"))). Resolved boundary times are concatenated in the order supplied. Combined specifications must use the by_ helpers directly: raw values are ignored with a warning.

Time span window

span additively expands the time span window around interval boundaries.

• A two-value vector expands the start and end, respectively: span = c(-60, 60) expands the start earlier by 60, and the end later by 60. For example, ⁠start = by_time(30), end = by_time(60), span = c(-5, 10)⁠ returns an interval of ⁠[25, 70]⁠.

• A single numeric value is recycled according to the sign: span = -60 becomes c(-60, 0) to expand the start earlier. span = 60 becomes c(0, 60) to expand the end later.

• If only start is specified alone, both span values expand the single boundary window: ⁠start = by_time(30), span = c(-5, 60)⁠ returns ⁠[25, 90]⁠.

Per-group channel selection with group_channels

When group_intervals = "ensemble" or a list of numeric grouped intervals, group_channels can be specified as a list of column names to override which channels are ensemble-averaged within each group. For example, to exclude a channel from one interval:

group_channels = list(
  c(A, B, C),
  c(A, C) ## channel "B" data are excluded from the second interval
)

If all grouped intervals include all nirs_channels, group_channels can be left as NULL (the default) and all channels are ensemble-averaged within every group.

Grouping intervals

group_intervals controls whether extracted intervals are returned as distinct data frames or ensemble-averaged.

"distinct"

The default. Extract each interval and return a list of independent data frames.

"ensemble"

Ensemble-average each specified nirs_channel across all detected intervals and return a one-item list with a single data frame.

list(c(1, 2), c(3, 4))

Ensemble-average each specified nirs_channel within each group and return a list with one data frame for each group. Any intervals detected but not specified in group_intervals are returned as distinct.

group_intervals lists can be named (e.g. list(low = c(1, 2), high = c(3, 4))) and will pass those names to the returned list of data frames.

When group_intervals is a list of numeric interval numbers, list items in group_channels and span are recycled to the number of groups. If lists are only partially specified, the final item is recycled forward as needed. Extra items are ignored.

Value

A named list() of tibbles of class "mnirs", each with metadata available via attributes(). When data is a list of data frames, results are flattened into a single-layer list with interval names indicating nested number of data frame and interval as ⁠interval_<df>.<interval>⁠. Other interval names (e.g. "ensemble" or custom group names) are suffixed as ⁠<name>_<df>⁠.

Examples

## read example data
data <- read_mnirs(
    example_mnirs("train.red"),
    nirs_channels = c(
        smo2_left = "SmO2 unfiltered",
        smo2_right = "SmO2 unfiltered"
    ),
    time_channel = c(time = "Timestamp (seconds passed)"),
    zero_time = TRUE,
    verbose = FALSE
) |>
    ## avoid issues ensemble-averaging irregular samples
    resample_mnirs(method = "linear", verbose = FALSE)

## ensemble-average across multiple intervals
interval_list <- extract_intervals(
    data,                         ## channels recycled to all intervals by default
    nirs_channels = c(smo2_left, smo2_right),
    group_intervals = "ensemble", ## ensemble-average across two intervals
    start = by_time(368, 1084),   ## manually identified interval start times
    span = c(-20, 90),            ## include the last 180-sec of each interval (recycled)
    zero_time = TRUE              ## re-calculate common time to start from `0`
)

interval_list[[1L]]


  if (requireNamespace("ggplot2", quietly = TRUE)) {
    plot(interval_list, time_labels = TRUE) +
      ggplot2::geom_vline(xintercept = 0, linetype = "dotted")
  }

Extract earliest POSIXct value from file header metadata

Description

Extract earliest POSIXct value from file header metadata

Usage

extract_start_timestamp(file_header)

Apply a Butterworth digital filter

Description

Apply a Butterworth digital filter to vector data with signal::butter() and signal::filtfilt() which handles 'edges' better at the start and end of the data.

Usage

filter_butterworth(
  x,
  order = 2L,
  W,
  type = c("low", "high", "stop", "pass"),
  edges = c("rev", "rep1", "none"),
  na.rm = FALSE,
  ...
)

filter_butter(
  x,
  order = 2L,
  W,
  type = c("low", "high", "stop", "pass"),
  edges = c("rev", "rep1", "none"),
  na.rm = FALSE,
  ...
)

Arguments

Details

Applies a centred (two-pass symmetrical) Butterworth digital filter from signal::butter() and signal::filtfilt().

Filter type defines how the desired signal frequencies are either passed or rejected from the output signal. Low-pass and high-pass filters allow only frequencies lower or higher than the cutoff frequency W to be passed through as the output signal, respectively. Stop-band defines a critical range of frequencies which are rejected from the output signal. Pass-band defines a critical range of frequencies which are passed through as the output signal.

The filter order (number of passes) is defined by order, typically in the range ⁠order = [1, 10]⁠. Higher filter order tends to capture more rapid changes in amplitude, but also causes more distortion around those change points in the signal. General advice is to use the lowest filter order which sufficiently captures the desired rapid responses in the data.

The critical (cutoff) frequency is defined by W, a numeric value for low-pass and high-pass filters, or a two-element vector c(low, high) defining the lower and upper bands for stop-band and pass-band filters. W represents the desired fractional cutoff frequency in the range ⁠W = [0, 1]⁠, where 1 is the Nyquist frequency, i.e., half the sample rate of the data in Hz.

Missing values (NA) in x will cause an error unless na.rm = TRUE. Then NAs will be ignored and passed through to the returned vector.

Value

A numeric vector the same length as x.

See Also

signal::filtfilt(), signal::butter()

Examples

set.seed(13)
sin <- sin(2 * pi * 1:150 / 50) * 20 + 40
noise <- rnorm(150, mean = 0, sd = 6)
noisy_sin <- sin + noise
without_edge_detection <- filter_butterworth(
    x = noisy_sin,
    order = 2,
    W = 0.1,
    edges = "none"
)
with_edge_detection <- filter_butterworth(
    x = noisy_sin,
    order = 2,
    W = 0.1,
    edges = "rep1"
)


ggplot2::ggplot(data.frame(), ggplot2::aes(x = seq_along(noise))) +
    theme_mnirs() +
    scale_colour_mnirs(name = NULL) +
    ggplot2::geom_line(ggplot2::aes(y = noisy_sin)) +
    ggplot2::geom_line(
        ggplot2::aes(y = without_edge_detection, colour = "without")
    ) +
    ggplot2::geom_line(
        ggplot2::aes(y = with_edge_detection, colour = "with")
    )

Filter a data frame

Description

Apply digital filtering/smoothing to numeric vector data within a data frame using either:

1. A cubic smoothing spline.

2. A Butterworth digital filter.

3. A simple moving average.

Note the method-specific arguments below.

Usage

filter_mnirs(
  data,
  nirs_channels = NULL,
  time_channel = NULL,
  method = c("smooth_spline", "butterworth", "moving_average"),
  na.rm = FALSE,
  verbose = TRUE,
  ...,
  spar = NULL,
  order = 2L,
  W = NULL,
  fc = NULL,
  sample_rate = NULL,
  type = c("low", "high", "stop", "pass"),
  edges = c("rev", "rep1", "none"),
  width = NULL,
  span = NULL,
  partial = FALSE
)

Arguments

Details

method = "smooth_spline"

Aliases: method = c("smooth spline", "spline")

Applies a non-parametric cubic smoothing spline from stats::smooth.spline(). Smoothing is defined by the parameter spar, which can be left as NULL and automatically determined via penalised log likelihood. This usually works well for responses occurring on the order of minutes or longer. spar can be specified typically, but not necessarily, in the range ⁠spar = [0, 1]⁠.

method = "butterworth"

Aliases: method = c("butter")

Applies a centred (two-pass symmetrical) Butterworth digital filter from signal::butter() and signal::filtfilt().

Filter type defines how the desired signal frequencies are either passed or rejected from the output signal. Low-pass and high-pass filters allow only frequencies lower or higher than the cutoff frequency, respectively to be passed through to the output signal. Stop-band defines a critical range of frequencies which are rejected from the output signal. Pass-band defines a critical range of frequencies which are passed through as the output signal.

The filter order (number of passes) is defined by order, typically in the range ⁠order = [1, 10]⁠. Higher filter order tends to capture more rapid changes in amplitude, but also causes more distortion around those change points in the signal. General advice is to use the lowest filter order which sufficiently captures the desired rapid responses in the data.

The critical (cutoff) frequency can be defined by W, a numeric value for low-pass and high-pass filters, or a two-element vector c(low, high) defining the lower and upper bands for stop-band and pass-band filters. W represents the desired fractional cutoff frequency in the range ⁠W = [0, 1]⁠, where 1 is the Nyquist frequency, i.e., half the sample_rate of the data in Hz.

Alternatively, the cutoff frequency can be defined by fc and sample_rate together. fc represents the desired cutoff frequency directly in Hz, and sample_rate is the sample rate of the recorded data in Hz. Where W = fc / (sample_rate / 2).

Only one of either W or fc should be defined. If both are defined, W will be preferred over fc.

method = "moving_average"

Aliases: method = c("moving average", "ma")

Applies a centred (symmetrical) moving average filter in a local window, defined by either width as the number of samples around idx between ⁠[idx - floor(width/2), idx + floor(width/2)]⁠. Or by span as the timespan in units of time_channel between ⁠[t - span/2, t + span/2]⁠.

Missing values

Missing values (NA) in nirs_channels will cause an error for method = "smooth_spline" or "butterworth", unless na.rm = TRUE. Then NAs will be ignored and passed through to the returned data.

For method = "moving_average", na.rm controls whether NAs within each local window are either propagated to the returned vector when na.rm = FALSE (the default), or ignored before processing if na.rm = TRUE.

Value

A tibble of class "mnirs" with metadata available with attributes(). For list or grouped data frame input, returns a named list of "mnirs" tibbles, one per interval.

Data input formats

mnirs processing functions accept data in multiple formats:

• A single "mnirs" data frame is processed and returned directly.

• A list of "mnirs" data frames: each interval is processed separately and returned as a named list.

• A grouped "mnirs" data frame, e.g. with dplyr::group_by(): the data frame is split by grouping levels and each group is processed as a separate interval, returned as a named list.

Per-channel arguments

Arguments apply globally to all nirs_channels by default. Relevant arguments can instead be supplied uniquely per-channel as a named list(), with names matching nirs_channels, e.g.

replace_mnirs(
    data,
    nirs_channels = c(hhb, smo2),
    invalid_values = list(hhb = -1, smo2 = c(0, 100)),
    invalid_above = list(hhb = 10),
    span = list(3, hhb = 5)
)

• A non-list value applies to every channel (the default behaviour).

• A list() named by nirs_channels applies per-channel values.

• A single unnamed value in the list will be applied to unlisted channels (e.g. span = list(3, hhb = 5) gives hhb 5 and every other channel 3). If no unnamed fallback value in the list, channels not named in the list will be returned un-processed (e.g. span = list(hhb = 5) will only process hhb).

• list() names not matching nirs_channels are warned about and ignored.

Examples

## read example data and clean for outliers
data <- read_mnirs(
    file_path = example_mnirs("moxy_ramp"),
    nirs_channels = c(smo2 = "SmO2 Live"),
    time_channel = c(time = "hh:mm:ss"),
    verbose = FALSE
) |>
    replace_mnirs(
        invalid_values = c(0, 100),
        outlier_cutoff = 3,
        width = 7,
        verbose = FALSE
    )

data

data_filtered <- filter_mnirs(
    data,                   ## blank channels will be retrieved from metadata
    method = "butterworth", ## Butterworth digital filter is a common choice
    order = 2,              ## filter order number
    W = 0.02,               ## filter fractional critical frequency `[0, 1]`
    type = "low",           ## specify a "low-pass" filter
    na.rm = TRUE            ## explicitly ignore NAs
)

## note the smoothed `smo2` values
data_filtered


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ## plot filtered data on top of raw to compare
        plot(data_filtered, time_labels = TRUE) +
            ggplot2::geom_line(
                data = data,
                ggplot2::aes(y = smo2, colour = "smo2"), alpha = 0.4
            )
    }

Apply a moving average filter

Description

Apply a simple moving average smoothing filter to vector data. filter_ma() is an alias of filter_moving_average().

Usage

filter_moving_average(
  x,
  t = seq_along(x),
  width = NULL,
  span = NULL,
  partial = FALSE,
  na.rm = FALSE,
  verbose = TRUE,
  ...
)

filter_ma(
  x,
  t = seq_along(x),
  width = NULL,
  span = NULL,
  partial = FALSE,
  na.rm = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

Rolling window

Applies a centred (symmetrical) moving average filter in a local window, defined by either width as the number of samples around idx between ⁠[idx - floor(width/2), idx + floor(width/2)]⁠. Or by span as the timespan in units of time_channel between ⁠[t - span/2, t + span/2]⁠.

Partial windows

The default partial = FALSE requires a complete number of samples specified by width or span (estimated from the sample rate of t when span is used). NA is returned if fewer samples are present in the local window.

Setting partial = TRUE allows computation with only a single valid sample, such as at edge conditions. But these values will be more sensitive to noise and should be used with caution.

Missing values

na.rm controls whether missing values (NAs) within each local window are either propagated to the returned vector when na.rm = FALSE (the default), or ignored before processing if na.rm = TRUE.

Value

A numeric vector the same length as x.

Examples

x <- c(1, 3, 2, 5, 4, 6, 5, 7)
t <- c(0, 1, 2, 4, 5, 6, 7, 10)  ## irregular time with gaps

## width: centred window of 3 samples
filter_moving_average(x, width = 3)

## partial = TRUE fills edge values with a narrower window
filter_moving_average(x, width = 3, partial = TRUE)

## span: centred window of 2 time-units (accounts for irregular sampling)
filter_moving_average(x, t, span = 2)

## na.rm = FALSE (default): any NA in the window propagates to the result
x_na <- c(1, NA, 3, 4, 5, NA, 7, 8)
filter_moving_average(x_na, width = 3, na.rm = FALSE)

## na.rm = TRUE: skip NAs and return the local mean of local valid values
filter_moving_average(x_na, width = 3, partial = TRUE, na.rm = TRUE)

Find the header row containing all nirs_channels

Description

Find the header row containing all nirs_channels

Usage

find_header_row(raw, nirs_channels, start = 1L, env = rlang::caller_env())

Arguments

resolve a single mnirs_interval object to time values

Description

resolve a single mnirs_interval object to time values

Usage

find_interval_time(
  interval,
  t_vec,
  event_vec = NULL,
  position = c("first", "last"),
  env = rlang::caller_env()
)

Arguments

Find valid model-fitting indices up to the first extreme

Description

Filters x and t to valid finite values, locates the first valid peak (maximum) or trough (minimum) where t >= 0, and returns the integer indices of all finite observations up to end_window past that extreme.

Usage

find_kinetics_idx(
  x,
  t = seq_along(x),
  end_window = Inf,
  direction = c("auto", "positive", "negative"),
  ...,
  env = rlang::caller_env()
)

Arguments

Details

Direction detection

When direction = "auto", the excursions of x above and below its initial baseline (the median of the earliest samples) are compared via detect_direction(). If the upward excursion dominates, the function searches for a peak (maximum); if the downward excursion dominates, a trough (minimum). When the excursions tie, the direction is determined by comparing abs(max(x)) to abs(min(x)), with ties defaulting to "positive".

Negative time handling

Only samples where t >= 0 are used for detecting the extreme, allowing pre-baseline (negative time) data to be excluded from the search. However, indices where t < 0 are included in the returned vector provided they are finite.

Value

A named list with three elements:

direction

Character; the resolved direction used – "positive" (peak) or "negative" (trough).

extreme

Integer or NULL; the index of the first qualifying peak or trough in original x space, or NULL if no qualifying extreme was found (monotonic, horizontal, or degenerate input).

idx

Integer vector of all valid finite indices, truncated at t[extreme] + end_window.

Fit a biexponential model to one channel

Description

Channel fitter of analyse_biexponential() (see analyse_kinetics_channels()), in two stages. Stage 1 fits the fast phase as a monoexponential on the end_window window (fit_monoexponential()). Stage 2 fits the full SSbiexponential() model on the whole response via stats::nls() with algorithm = "port", A, tau, and TD box-bounded about their stage-1 values by the ⁠*_flex⁠ half-widths and B, B2, tau2 free, seeded by biexp_start() with the fast phase held. A failed stage returns NA, and the fallback chain resolves the row upstream.

Usage

fit_biexponential(x, t, valid, .a, ctx)

Arguments

Value

The coefs/model/fitted_data/diag list of build_fit_results(), or build_na_results() when a stage fails.

Merge user nls control over a fit's internal defaults

Description

Merge user nls control over a fit's internal defaults

Usage

fit_control(control, ...)

Arguments

Value

A control list for stats::nls().

Fit an exponential-drift model to one channel

Description

Channel fitter of analyse_exponential_drift() (see analyse_kinetics_channels()). Self-starting SSexponential_drift() via stats::nls() with algorithm = "port", seeded by expdrift_start(); a failed 6-parameter fit falls back to the 5-parameter model (fit_td_fallback()), and the requested direction is enforced on B - A (enforce_direction()).

Usage

fit_exponential_drift(x, t, valid, .a, ctx)

Arguments

Value

The coefs/model/fitted_data/diag list of build_fit_results(), or build_na_results() when the fit fails.

Fit a monoexponential model to one channel

Description

Channel fitter of analyse_monoexponential() (see analyse_kinetics_channels()), also the fast-phase (stage 1) fit of fit_biexponential() and the fallback of fit_exponential_drift(). Self-starting SSmonoexponential() via stats::nls(); a failed 4-parameter fit falls back to the 3-parameter model (fit_td_fallback()), and the requested direction is enforced on B - A (enforce_direction()).

Usage

fit_monoexponential(x, t, valid, .a, ctx)

Arguments

Value

The coefs/model/fitted_data/diag list of build_fit_results(), or build_na_results() when the fit fails.

Alias fit column names that collide with model parameters

Description

stats::nls() formula symbols must be disjoint: a name cannot be both a data column and a parameter. The fit data frame carries the channel names so the model predicts on them, so a channel named after a model parameter (or D, the amplitude used by enforce_direction()) is prefixed with . in the fit and in the stored model formula. Coefficients and results keep the original names.

Usage

fit_names(x, t, params)

Arguments

Value

A length-2 character vector of the response and time column names to fit on.

Fit a sigmoidal model to one channel

Description

Channel fitter of analyse_logistic() (see analyse_kinetics_channels()), also the fallback of fit_sigmoidal_drift(). Self-starting SSlogistic(), SSgompertz(), or SSgompertz_left() per the channel shape via stats::nls(), with the requested direction enforced on B - A and the sign of slope (enforce_direction()).

Usage

fit_sigmoidal(x, t, valid, .a, ctx)

Arguments

Value

The coefs/model/fitted_data/diag list of build_fit_results(), or build_na_results() when the fit fails.

Fit a sigmoidal-drift model to one channel

Description

Channel fitter of analyse_sigmoidal_drift() (see analyse_kinetics_channels()). Self-starting SSsigmoidal_drift() of the channel shape via stats::nls() with algorithm = "port", seeded by sigdrift_start(), with the requested direction enforced on B - A and the sign of slope (enforce_direction()).

Usage

fit_sigmoidal_drift(x, t, valid, .a, ctx)

Arguments

Value

The coefs/model/fitted_data/diag list of build_fit_results(), or build_na_results() when the fit fails.

Fit a self-start model with time-delay fallback

Description

Shared attempt skeleton for the nls-based kinetics workers. The TD model is flat at A before TD, so the pre-onset baseline anchors A; the reduced model has no such region and diverges at t < 0, so it is fit from start_time onward. An under-determined attempt is rejected before it reaches fitter, and a failed TD fit falls back to the reduced model without TD unless TD is user-fixed. Every failure is reported through warn_fit_failed().

Usage

fit_td_fallback(x_fit, t_fit, params, .a, fitter, fn, ctx)

Arguments

Value

A list with model (or NULL), the params actually fit, the logical row filter keep, and the fit data frame.

Format time span data as h:mm:ss

Description

Convert numeric time span data to h:mm:ss format for pretty plotting. Inspired by ggplot2::scale_x_time().

Usage

format_hmmss(x)

Arguments

Details

If all values are less than 3600 (1 hour), then format is returned as mm:ss. If any value is greater than 3600, format is returned as h:mm:ss with leading zeroes.

Value

A character vector the same length as x.

Examples

x <- 0:120
y <- sin(2 * pi * x / 15) + rnorm(length(x), 0, 0.2)

ggplot2::ggplot(data.frame(x, y), ggplot2::aes(x, y)) +
    theme_mnirs() +
    ggplot2::scale_x_continuous(
        breaks = breaks_timespan(),
        labels = format_hmmss
    ) +
    ggplot2::geom_line()

Free parameters of a self-start model call

Description

A parameter written as a bare symbol in the model call is free (fitted by stats::nls()); one written as a constant or expression is fixed. Used by the model functions to return gradient columns for the free parameters only, in call order, as stats::nls() indexes the "gradient" attribute by position.

Usage

free_params(mCall, params)

Arguments

Value

A character vector; the subset of params that are free.

Combine fitted and fixed coefficients into the full parameter vector

Description

Combine fitted and fixed coefficients into the full parameter vector

Usage

full_coefs(model, params, fix = list())

Arguments

Value

A named numeric vector ordered by params containing the fitted coefficients with fixed values merged in.

Gompertz growth functions

Description

Calculate 4-parameter Gompertz (asymmetric sigmoidal) curves. Model families fit by analyse_kinetics() with method = "sigmoidal" and shape = "gompertz" or "gompertz_left", and by stats::nls() via the self-starting wrappers SSgompertz() and SSgompertz_left().

Usage

gompertz(t, A, B, xmid, slope)

gompertz_left(t, A, B, xmid, slope)

Arguments

Details

gompertz() (right-Gompertz) is asymmetric with the inflection point xmid closer to the starting asymptote A: early acceleration away from A, and a slow approach to the ending asymptote B. Appropriate for fast-onset, slow-tail responses.

gompertz_left() (left-Gompertz) has the inflection point closer to the ending asymptote B: slow departure from A, and late acceleration toward B. Appropriate for slow-onset, fast-tail responses.

Model equations

Both forms are re-parameterised so xmid is the time at inflection and slope is the response rate dx/dt at the inflection, with k = slope * e / (B - A).

• gompertz(): A + (B - A) * exp(-exp(-k * (t - xmid))). Inflection height fixed at A + (B - A) / e; 36.8% of the amplitude.

• gompertz_left(): A + (B - A) * (1 - exp(-exp(k * (t - xmid)))). Inflection height fixed at A + (B - A) * (1 - 1/e); 63.2% of the amplitude.

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

analyse_kinetics(), SSgompertz(), SSgompertz_left(), logistic(), sigmoidal_drift()

Examples

## create a Gompertz curve with random noise
set.seed(15)
t <- 1:60
x <- gompertz(t, A = 10, B = 100, xmid = 30, slope = 4) +
    rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## fit with the self-starting wrapper
model <- nls(x ~ SSgompertz(t, A, B, xmid, slope), data = data)
summary(model)

y <- predict(model, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y))
    }

Initiate self-starting Gompertz model

Description

gompertz_init(): Returns initial values for the parameters in a selfStart model. Used by both SSgompertz() and SSgompertz_left(); the symmetric logistic linearisation does not apply to Gompertz forms, so initialisation is derivative-based via init_inflection().

Usage

gompertz_init(mCall, data, LHS, ...)

Arguments

Value

gompertz_init(): Initial starting estimates for parameters in the model called by SSgompertz() or SSgompertz_left().

Convert ⁠H:MM(:SS.fff)⁠ strings to seconds of day

Description

Convert ⁠H:MM(:SS.fff)⁠ strings to seconds of day

Usage

hms_to_seconds(x)

Estimate baseline and asymptote from the first/last quintile of x

Description

Shared helper used by self-start initialisers for logistic / Gompertz model families.

Usage

init_asymptotes(x, n = length(x))

Arguments

Value

A list with elements A (starting asymptote estimate) and B (ending asymptote estimate).

Wrap a self-start initialiser to support fixed parameters

Description

Decorates a selfStart initial function so parameters supplied as values in the model formula (e.g. SSmonoexponential(t, A = 0, B, tau)) are excluded from the returned start vector. stats::nls() reads the free parameters from the names of that vector, so excluded parameters are treated as constants in the formula. Fixed values are forwarded to the wrapped initialiser as a fixed list argument to seed the remaining free estimates.

Usage

init_fixed(init, params)

Arguments

Value

A function suitable for the initial argument of stats::selfStart().

Estimate inflection point from a smoothed first derivative

Description

Shared helper that locates the empirical inflection (peak of ⁠|dx/dt|⁠ after smoothing) and returns the corresponding xmid and slope initial values. Falls back to the half-response point and a mean-rate slope when the derivative is degenerate.

Usage

init_inflection(x, t, A_init, B_init)

Arguments

Value

A list with elements idx (integer index into x), xmid (numeric t value at the inflection), and slope (numeric dx/dt at the inflection).

Classify a per-channel/per-interval argument map

Description

An argument is a map when it is a list() with at least one named element and at most one unnamed element (the fallback for unlisted keys). Shared by resolve_channel_args() and resolve_interval_args().

Usage

is_arg_map(x)

Arguments

Value

A logical scalar.

Detect empty or NA strings

Description

Detect empty or NA strings

Usage

is_empty(x)

Build per-panel kinetics marker and label annotations

Description

Maps a fitted mnirs_kinetics method to its key coefficient markers (xval, yval) and formatted label lines for plot.mnirs_kinetics(). Marker rows are one per nirs_channel per key point per interval, with x-coordinates the resolved onset plus the method's time coefficient. Label rows are one per label line, anchored (xval = Inf, yval = -Inf or Inf) at the corner of the panel's right edge vacated by the observed signal: the bottom corner when the median of all channels over the right half of the interval sits above the y-axis midpoint, otherwise the top. vjust stacks the lines inward from the corner in channel order within each interval.

Usage

kinetics_annotations(x, free_y = FALSE)

Arguments

Value

A data.frame with columns interval, nirs_channels, xval, yval, label, and vjust. Marker rows have an empty label and NA vjust; label rows have infinite xval/yval. Rows are annotated by the model that fit them (the model coefficient column where the method has a fallback chain, else the method). NULL for a method with no annotation spec, in which case plot.mnirs_kinetics() draws the fitted curve alone.

Zero-row kinetics warnings scaffold

Description

Stable column template for captured fit conditions, so binding and the returned warnings element keep consistent columns when none fire.

Usage

kinetics_warnings_df()

Value

A zero-row data.frame with columns interval, nirs_channels, type, and message.

Generalised logistic function

Description

Calculate a 4- or 5-parameter logistic (sigmoidal) curve. The 4-parameter symmetric form is fit by analyse_kinetics() with method = "sigmoidal" and shape = "symmetric" (default), and by stats::nls() via the self-starting wrapper SSlogistic().

Usage

logistic(t, A, B, xmid, slope, asym = NULL)

Arguments

Details

The 5-parameter Richards form is exported for advanced use directly with stats::nls() but is not used by analyse_kinetics() due to convergence instability. For asymmetric responses, prefer gompertz() / gompertz_left(), which are more stable.

Model equations

Both forms are re-parameterised from the Richards generalised logistic model so xmid is the time at inflection and slope is the response rate dx/dt at the inflection.

• 4-parameter (symmetric): A + (B - A) / (1 + exp(-4 * slope * (t - xmid) / (B - A)))

• 5-parameter (asymmetric): A + (B - A) / (1 + exp(-k * (t - xmid)))^(1 / v) with v = -log(2) / log(asym) and k = 2 * slope * v / ((B - A) * asym).

The inflection is at t = xmid with dx/dt = slope and y(xmid) = A + (B - A) * asym for any asym in ⁠(0, 1)⁠:

• asym = 0.5 (v = 1) collapses to the 4-parameter form.

• asym -> 0 gives an early-acceleration curve (inflection near A).

• asym -> 1 gives a late-acceleration curve (inflection near B).

• asym = 0.368 (1/e) approximates a right-inflection gompertz() curve.

• asym = 0.632 (1 - 1/e) approximates a left-inflection gompertz_left() curve.

Value

A numeric vector of predicted values the same length as the predictor variable t.

See Also

analyse_kinetics(), SSlogistic(), gompertz(), gompertz_left(), sigmoidal_drift(), monoexponential()

Examples

## create an asymmetric logistic curve with random noise
set.seed(15)
t <- 1:60
x <- logistic(t, A = 10, B = 100, xmid = 30, slope = 4, asym = 0.3) +
    rnorm(length(t), 0, 2)
data <- data.frame(t, x)

## 5-parameter fit with the self-starting wrapper
model <- nls(x ~ SSlogistic(t, A, B, xmid, slope, asym), data = data)
summary(model)

y <- predict(model, data)


    if (requireNamespace("ggplot2", quietly = TRUE)) {
        ggplot2::ggplot(data, ggplot2::aes(t, x)) +
            theme_mnirs() +
            ggplot2::geom_point() +
            ggplot2::geom_line(ggplot2::aes(y = y))
    }

Initiate self-starting logistic model

Description

logistic_init(): Returns initial values for the parameters in a selfStart model.

Usage

logistic_init(mCall, data, LHS, ...)

Arguments

Value

logistic_init(): Initial starting estimates for parameters in the model called by SSlogistic().

Apply an mnirs function over each interval of a multi-interval input

Description

Shared entry point for ⁠*_mnirs()⁠ transformer functions accepting a list of data frames or a grouped data frame. Normalises data to a named list via as_data_list(), then re-evaluates the captured user-facing call once per interval with data swapped, so all arguments (including NSE channel expressions) are forwarded verbatim.

Usage

map_mnirs_intervals(data, call, eval_env, env = rlang::caller_env())

Arguments

Value

A named list of class "mnirs" containing processed "mnirs" data frames, one per interval.

Data input formats

mnirs processing functions accept data in multiple formats:

• A single "mnirs" data frame is processed and returned directly.

• A list of "mnirs" data frames: each interval is processed separately and returned as a named list.

• A grouped "mnirs" data frame, e.g. with dplyr::group_by(): the data frame is split by grouping levels and each group is processed as a separate interval, returned as a named list.

Metadata names of class "mnirs", retrieved with attr()

Description

Metadata names of class "mnirs", retrieved with attr()

Usage

mnirs_metadata

Initiate self-starting monoexponential model

Description

monoexp_init(): Returns initial values for the parameters in a selfStart model.

Usage

monoexp_init(mCall, data, LHS, ...)

Arguments

Value

monoexp_init(): Initial starting estimates for parameters in the model called by SSmonoexponential().

Monoexponential model with gradient

Description

Model function of SSmonoexponential(): monoexponential() plus the partial derivatives for the parameters written as bare symbols in the call (see free_params()), so stats::nls() skips stats::numericDeriv() and a parameter fixed as a constant in the formula contributes no gradient column.

Usage

monoexp_model(t, A, B, tau, TD = NULL)

Arguments

Value

A numeric vector of predicted values with a "gradient" attribute when any parameter is free.

Grid-profiled starting estimates for the monoexponential model

Description

Vector-level initialiser behind monoexp_init(), called directly by the kinetics worker on the fit window. Profiles tau (and TD for the 4-parameter model) on a coarse grid and keeps the RSS-minimising start (cf. biexp_start()). The model is linear in A and B once tau and TD are held, so the asymptotes are solved by least squares at every grid point at once. Point estimates from derivative changepoints or log-linearisation are too sensitive to noise, overshoot, and plateau data on real NIRS signals, and can strand nls with a singular gradient.

Usage

monoexp_start(x, t, fixed = list(), has_TD = FALSE)

Arguments