In fmmattioni/whippr: Tools for Manipulating Gas Exchange Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  fig.align = "center",
  dpi = 300
)

# from https://github.com/r-lib/crayon/issues/24#issuecomment-581068792
# crayon needs to be explicitly activated in Rmd
options(crayon.enabled = TRUE)
# Hooks needs to be set to deal with outputs
# thanks to fansi logic
old_hooks <- fansi::set_knit_hooks(knitr::knit_hooks,
                                   which = c("output", "message", "error"))

```{block, type = 'rmdinfo'} Here you may find a walk-through on how to perform VO2 kinetics analysis in the moderate-intensity domain.

Functions for analysis in the heavy- and severe-intensity domains will be added in the near future.

For making everyone's life easier, the general function `vo2_kinetics()` was created. This function calls smaller separate functions that **fully** automate the VO2 kinetics data analysis. The following interactive tree diagram shows how each function is called:

```r
library(dplyr, warn.conflicts = FALSE)
library(collapsibleTree)

normalize_transitions <- tibble(
  a = "vo2_kinetics()",
  b = "detect_outliers()",
  c = "normalize_transitions()",
  d = "normalize_first_breath()"
)

predict_bands <- tibble(
  a = "vo2_kinetics()",
  b = "detect_outliers()",
  c = "predict_bands()",
  d = c("predict_bands_baseline()", "predict_bands_transition()")
)

plot_outliers <- tibble(
  a = "vo2_kinetics()",
  b = "plot_outliers()"
)

process_data <- tibble(
  a = "vo2_kinetics()",
  b = "process_data()",
  c = c("interpolate()", "perform_average()", "normalize_time()")
)

perform_kinetics <- tibble(
  a = "vo2_kinetics()",
  b = "perform_kinetics()",
  c = c("get_residuals()", "model_diagnostics()")
)

bind_rows(normalize_transitions, predict_bands, plot_outliers, process_data, perform_kinetics) %>% 
  collapsibleTree(
    hierarchy = c("a", "b", "c", "d"),
    root = "VO2 kinetics analysis",
    width = 800,
    height = 300,
    zoomable = FALSE
  )

Read the data

The first step is to read the raw data with the read_data() function. Here we are going to use the example file that comes with {whippr}, which is a file exported from the COSMED metabolic cart.

library(whippr)

raw_data <- read_data(path = system.file("example_cosmed.xlsx", package = "whippr"), metabolic_cart = "cosmed", time_column = "t")

raw_data

As you can see in the following graph, this is a protocol where 3 transitions from a baseline exercise intensity to to an exercise intensity below the gas exchange threshold. Therefore, this is a VO2 kinetics test in the moderate-intensity domain. In this specific case, the following was done:

• 3 x 6-min baseline periods at 20 W.
• 3 x 6-min transition periods at the power output associated with 90% of the gas exchange threshold.

library(ggplot2)

raw_data %>% 
  ggplot(aes(t, VO2)) +
  geom_point(shape = 21, size = 3, fill = "white") +
  theme_whippr()

Perform the analysis

After reading the raw data, we can move directly to performing the VO2 kinetics analysis with vo2_kinetics(). This function will:

• Recognize each baseline and transition phase
• Normalize the first breath in each transition in a safe way to prevent time misalignment
• Recognize outliers
• Remove outliers
• Interpolate each transition
• Time-align the data
• Ensemble-average the transitions
• Perform the chosen bin-average
• Fit the final mono-exponential model from VO2 kinetics from the options chosen
• Calculate residuals

```{block, type = 'rmdinfo'} For modeling VO2 kinetics analysis in the moderate-intensity domain, a mono-exponential model is used:

$$VO_2\left(t\right)=baseline+amplitude\cdot\left(1-e^{^{-\frac{\left(t-TD\right)}{τ}}}\right)$$

where:

• VO2(t) = the oxygen uptake at any given time.
• baseline = the oxygen uptake associated with the baseline phase.
• amplitude = the steady-state increase increase in oxygen uptake above baseline.
• TD = the time delay.
• τ = the time constant defined as the duration of time for the oxygen uptake to increase to 63% of the steady-state increase.

### Important options

In `vo2_kinetics()` you must set important options before continuing.

**Protocol-related options:**

* `protocol_n_transitions` = Number of transitions performed.
* `protocol_baseline_length` = The length of the baseline (in seconds).
* `protocol_transition_length` = The length of the transition (in seconds).

**Data cleaning-related options:**

* `cleaning_level` = A numeric scalar between 0 and 1 giving the confidence level for the intervals to be calculated during the data cleaning process. Breaths lying outside the prediction bands will be excluded.
* `cleaning_baseline_fit` = A vector of the same length as the number in `protocol_n_transitions`, indicating what kind of fit to perform for each baseline. Either *linear* or *exponential*.

**Fitting-related options:**

* `fit_level` = A numeric scalar between 0 and 1 giving the confidence level for the parameter estimates in the final VO2 kinetics fit.
* `fit_bin_average` = The bin average to be performed for the final fit.
* `fit_phase_1_length` = The length of the phase I that you wish to exclude from the final exponential fit, in seconds.
* `fit_baseline_length` = The length the baseline to perform the final linear fit, in seconds. 
* `fit_transition_length` = The length of the transition to perform the final exponential fit, in seconds.

The analysis is performed like the following:

```r
results <- vo2_kinetics(
  .data = raw_data,
  intensity_domain = "moderate",
  vo2_column = "VO2",
  protocol_n_transitions = 3,
  protocol_baseline_length = 360,
  protocol_transition_length = 360,
  cleaning_level = 0.95,
  cleaning_baseline_fit = c("linear", "exponential", "exponential"),
  fit_level = 0.95,
  fit_bin_average = 5,
  fit_phase_1_length = 20,
  fit_baseline_length = 120,
  fit_transition_length = 240,
  verbose = TRUE
)

results

Fit parameters

Fit parameters and confidence intervals may be accessed through the model_summary column.

results$model_summary[[1]]

Fit plot

The fit plot may be accessed through the plot_model column.

results$plot_model[[1]]

Checking what was done during data cleaning

The data cleaning process may be accessed through the plot_outliers column.

results$plot_outliers[[1]]

Model diagnostics

Model residuals plot may be accessed through the plot_residuals column.

results$plot_residuals[[1]]

Additional columns

Raw data with detected outliers

The raw data with additional columns from the data cleaning process may be accessed through the data_outliers column.

results$data_outliers[[1]]

Processed data

The processed data (cleaned, interpolated, time-aligned, ensemble-averaged, and bin-averaged) may be accessed through the data_processed column.

results$data_processed[[1]]

Fitted data

The data from the baseline and transition fits may be accessed through the data_fitted column.

results$data_fitted[[1]]

Model

The model used for fitting the mono-exponential model may be accessed through the model column.

results$model[[1]]

summary(results$model[[1]])

Residuals data

The model residuals data may be accessed through the model_residuals column.

results$model_residuals[[1]]

fmmattioni/whippr documentation built on Feb. 23, 2024, 11:20 a.m.