## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  out.width = "100%"
)

## ----setup--------------------------------------------------------------------
library(MFO)

## ----data---------------------------------------------------------------------
head(basal_df)   # resting measurement, 45 records = 15 min
head(MFO_df)     # graded MFO test, 72 records = 8 steps of 3 min
head(VO2max_df)  # maximal graded test

## ----vo2max-------------------------------------------------------------------
step_time <- 20
n_vo2max <- ceiling(30 / step_time)                         # 2 records
VO2max <- max(stats::filter(VO2max_df$VO2, rep(1 / n_vo2max, n_vo2max), sides = 1),
              na.rm = TRUE)
c(VO2max = VO2max, instantaneous_max = max(VO2max_df$VO2))

## ----basal--------------------------------------------------------------------
window <- 300 / step_time                                   # 15 records = 5 min
rest <- basal_df[-seq_len(window), ]                        # drop the first 5 min
cv <- sapply(seq_len(nrow(rest) - window + 1), function(i) {
  x <- rest$RER[i:(i + window - 1)]
  sd(x) / mean(x)
})
best <- which.min(cv)
basal_window <- rest[best:(best + window - 1), c("VO2", "VCO2", "RER", "HR")]
basal <- colMeans(basal_window)
round(basal, 3)

## ----steps--------------------------------------------------------------------
records_per_step <- 180 / step_time                         # 9
records_steady <- 60 / step_time                            # 3
mfo <- data.frame(MFO_df)
mfo$step <- rep(seq_len(nrow(mfo) %/% records_per_step), each = records_per_step)
position <- ave(mfo$step, mfo$step, FUN = seq_along)
steady <- mfo[position > records_per_step - records_steady, ]
steps <- aggregate(cbind(VO2, VCO2, RER, HR) ~ step, data = steady, FUN = mean)
steps$porc_VO2 <- steps$VO2 / VO2max * 100
round(steps, 2)

## ----equation-----------------------------------------------------------------
frayn <- function(VO2, VCO2) {
  data.frame(CHO = 4.55 * VCO2 / 1000 - 3.21 * VO2 / 1000,
             FAT = 1.67 * VO2 / 1000 - 1.67 * VCO2 / 1000)
}
steps <- cbind(steps, frayn(steps$VO2, steps$VCO2))
basal_ox <- frayn(basal[["VO2"]], basal[["VCO2"]])
round(steps[, c("step", "porc_VO2", "RER", "CHO", "FAT")], 3)
round(basal_ox, 3)

## ----negative-----------------------------------------------------------------
steps[steps$FAT < 0, c("step", "RER", "FAT")]
steps$FAT[steps$FAT < 0] <- 0                                # negative_FAT = "truncate"

## ----fit----------------------------------------------------------------------
p3 <- lm(FAT ~ porc_VO2 + I(porc_VO2^2) + I(porc_VO2^3), data = steps)
b <- coef(p3)
roots <- Re(polyroot(c(b[2], 2 * b[3], 3 * b[4])))           # FAT'(x) = 0
roots <- roots[roots >= min(steps$porc_VO2) & roots <= max(steps$porc_VO2)]
fat_at <- function(x) predict(p3, newdata = data.frame(porc_VO2 = x))
Fatmax <- roots[which.max(fat_at(roots))]
MFO_hand <- unname(fat_at(Fatmax))
c(MFO = MFO_hand, Fatmax = Fatmax, R2 = summary(p3)$r.squared)

## ----outputs------------------------------------------------------------------
grid <- seq(min(steps$porc_VO2), max(steps$porc_VO2), length.out = 1000)
zone <- range(grid[fat_at(grid) >= 0.9 * MFO_hand])
diff_ox <- steps$FAT - steps$CHO
i <- which(diff(sign(diff_ox)) != 0)[1]
crossover <- approx(diff_ox[i:(i + 1)], steps$porc_VO2[i:(i + 1)], xout = 0)$y
HR_Fatmax <- approx(steps$porc_VO2, steps$HR, xout = Fatmax)$y
round(c(zone_lower = zone[1], zone_upper = zone[2],
        crossover = crossover, HR_Fatmax = HR_Fatmax), 2)

## ----mfo-call-----------------------------------------------------------------
res <- MFO(step_time = 20, db_MFO = MFO_df, db_basal = basal_df,
           db_graded = VO2max_df, cv_var = "RER", author = "Frayn",
           energy_equivalents = "classic")
hand <- c(VO2max, MFO_hand, Fatmax, zone, crossover, HR_Fatmax,
          basal_ox$CHO, basal_ox$FAT)
package <- c(res$VO2max, res$MFO, res$FAT_MAX, res$FatMax_zone_lower,
             res$FatMax_zone_upper, res$crossover_point, res$FAT_MAX_HR,
             res$x_CHO, res$x_FAT)
all.equal(unname(hand), unname(package))

## ----mfo-db, fig.alt = "Fat and carbohydrate oxidation (g/min) against %VO2max for the example participant: observed steps, fitted P3 curves, shaded Fatmax zone and vertical lines at Fatmax and at the crossover point."----
res$MFO_db
res$MFO_plot

## ----helper-------------------------------------------------------------------
run_MFO <- function(...) {
  warnings <- character(0)
  out <- withCallingHandlers(
    suppressMessages(MFO(step_time = 20, db_MFO = MFO_df, db_basal = basal_df,
                         db_graded = VO2max_df, cv_var = "RER", author = "Frayn",
                         energy_equivalents = "classic", ...)),
    warning = function(w) {
      warnings <<- c(warnings, conditionMessage(w))
      invokeRestart("muffleWarning")
    })
  out$warnings <- warnings
  out
}

## ----methods------------------------------------------------------------------
methods <- c("P3", "P2", "MV", "MRER")
by_method <- lapply(methods, function(m) run_MFO(method = m))
knitr::kable(data.frame(
  method = methods,
  MFO = sapply(by_method, `[[`, "MFO"),
  Fatmax = sapply(by_method, `[[`, "FAT_MAX"),
  zone_lower = sapply(by_method, `[[`, "FatMax_zone_lower"),
  zone_upper = sapply(by_method, `[[`, "FatMax_zone_upper"),
  zone_method = sapply(by_method, `[[`, "FatMax_zone_method"),
  R2 = sapply(by_method, `[[`, "R_squared"),
  RMSE = sapply(by_method, `[[`, "RMSE")
), digits = 3, caption = "MFO (g/min) and Fatmax and zone (%VO2max) by method.")

## ----sin, fig.alt = "Fat oxidation as %MFO against %VO2max: observed steps, the basic SIN curve in grey and the fitted SIN curve in orange between its two zeros; d, s, t, Fatmin and R² in the subtitle."----
kin <- MFO_kinetics(res$MFO_db)
knitr::kable(data.frame(
  parameter = c("d", "s", "t"),
  estimate = c(kin$d, kin$s, kin$t),
  std_error = c(kin$d_se, kin$s_se, kin$t_se)
), digits = 3)
round(c(FatMax_SIN = kin$FatMax_SIN, Fatmin = kin$Fatmin, R2 = kin$R_squared), 2)
kin$MFO_kinetics_plot

## ----sin-sensitivity----------------------------------------------------------
sin_row <- function(label, fit) {
  k <- suppressWarnings(MFO_kinetics(fit$MFO_db))
  data.frame(setting = label, steps = nrow(fit$MFO_db), d = k$d, s = k$s,
             t = k$t, FatMax_SIN = k$FatMax_SIN, Fatmin = k$Fatmin)
}
knitr::kable(rbind(
  sin_row("negative_FAT = \"truncate\"", res),
  sin_row("negative_FAT = \"exclude\"", run_MFO(negative_FAT = "exclude")),
  sin_row("RER_cutoff = 1", run_MFO(RER_cutoff = 1))
), digits = 3, caption = "SIN parameters; Fatmax_SIN and Fatmin in %VO2max.")

## ----rer-cutoff---------------------------------------------------------------
cutoffs <- list(none = NULL, classic = "classic", safe = "safe",
                sedentary = "sedentary", trained = "trained")
by_cutoff <- lapply(cutoffs, function(cut) run_MFO(RER_cutoff = cut))
knitr::kable(data.frame(
  RER_cutoff = names(cutoffs),
  steps = sapply(by_cutoff, function(r) nrow(r$MFO_db)),
  MFO = sapply(by_cutoff, `[[`, "MFO"),
  Fatmax = sapply(by_cutoff, `[[`, "FAT_MAX"),
  R2 = sapply(by_cutoff, `[[`, "R_squared"),
  warnings = sapply(by_cutoff, function(r) length(r$warnings)),
  row.names = NULL
), digits = 3)
by_cutoff$trained$warnings

## ----include-basal------------------------------------------------------------
by_basal <- list(`FALSE` = res, `TRUE` = run_MFO(include_basal = TRUE))
knitr::kable(data.frame(
  include_basal = names(by_basal),
  rows = sapply(by_basal, function(r) nrow(r$MFO_db)),
  MFO = sapply(by_basal, `[[`, "MFO"),
  Fatmax = sapply(by_basal, `[[`, "FAT_MAX"),
  R2 = sapply(by_basal, `[[`, "R_squared"),
  FatMax_SIN = sapply(by_basal, function(r)
    suppressWarnings(MFO_kinetics(r$MFO_db))$FatMax_SIN),
  row.names = NULL
), digits = 3)

## ----report-values------------------------------------------------------------
HRmax <- max(VO2max_df$HR)
final <- run_MFO(HRmax = HRmax)
n_steps <- sum(final$MFO_db$Load != "basal")

## ----report, results = "asis", echo = FALSE-----------------------------------
paragraph <- sprintf(paste(
  "Substrate oxidation was computed from the mean VO2 and VCO2 of the last %d s",
  "of each %d-min step with the equations of Frayn (1983), with urinary nitrogen",
  "set to zero. VO2max (%.0f mL/min) was the highest %d-s rolling mean of VO2 in",
  "the maximal test. The resting value was the mean of the 5-min window with the",
  "lowest coefficient of variation of RER after discarding the first 5 min.",
  "Negative fat oxidation rates (RER > 1) were set to zero and the resting point",
  "was not included in the fit. MFO and Fatmax were determined with a",
  "third-degree polynomial of fat oxidation on %%VO2max fitted to %d steps",
  "(R² = %.2f). The Fatmax zone was defined as the intensities eliciting at least",
  "90%% of MFO (Achten et al., 2002). Fat oxidation kinetics were described with",
  "the SIN model (Chenevière et al., 2009). MFO was %.2f g/min at %.1f%% VO2max",
  "(%.0f bpm, %.1f%% HRmax; Fatmax zone %.1f–%.1f%% VO2max); the SIN parameters",
  "were d = %.2f, s = %.2f and t = %.2f, with Fatmin at %.1f%% VO2max. Analyses",
  "were run with the R package MFO, version %s, [protocol: ergometer, initial",
  "load and increment, fasting and pre-test conditions]."),
  60, 3, final$VO2max, final$VO2max_window, n_steps, final$R_squared,
  final$MFO, final$FAT_MAX, final$FAT_MAX_HR, final$FAT_MAX_porcHRmax,
  final$FatMax_zone_lower, final$FatMax_zone_upper, kin$d, kin$s, kin$t,
  kin$Fatmin, as.character(utils::packageVersion("MFO")))
cat("> ", paragraph, "\n")

## ----batch-data---------------------------------------------------------------
root <- file.path(tempdir(), "mfo-batch")
dir.create(root, showWarnings = FALSE)
scale_gas <- function(db, k) transform(db, VO2 = VO2 * k, VCO2 = VCO2 * k)
for (id in c("P01", "P02")) {
  k <- if (id == "P01") 1 else 1.1
  dir.create(file.path(root, id), showWarnings = FALSE)
  openxlsx::write.xlsx(scale_gas(basal_df, k), file.path(root, id, "basal_df.xlsx"))
  openxlsx::write.xlsx(scale_gas(MFO_df, k), file.path(root, id, "MFO_df.xlsx"))
  openxlsx::write.xlsx(scale_gas(VO2max_df, k), file.path(root, id, "VO2max_df.xlsx"))
}
list.files(root, recursive = TRUE)

## ----batch--------------------------------------------------------------------
batch <- MFOs(from = "folder", path = root,
              db_basal_name = "basal_df", db_MFO_name = "MFO_df",
              db_graded_name = "VO2max_df", step_time = 20, cv_var = "RER",
              author = "Frayn", energy_equivalents = "classic",
              save_plot = FALSE, save_result = FALSE)
knitr::kable(batch[, c("ID", "MFO_g_min", "FAT_MAX_perVO2max", "VO2max",
                       "FatMax_zone_lower", "FatMax_zone_upper", "d", "s", "t",
                       "Fatmin")], digits = 3)

## ----cleanup, include = FALSE-------------------------------------------------
unlink(root, recursive = TRUE)