Analysis of Garmin tcx heart rate data from a stepwise treadmill test to identify lactate threshold

stepwise_treadmill_test.Rmd

---
title: "Stepwise treadmill test"
author: "Levi Waldron"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

This code:
* removes observations from before and after the stepwise test
* assumes starting at 6km/h and increasing 1km/h each 90 seconds
* fits a change-point linear model where the change-point is the estimate of the lactate threshold
* plots results

```{r}
suppressPackageStartupMessages({
  library(trackeR)
  library(dplyr)
  library(ggplot2)
})
```

# Load data and add speed and time variables

```{r}
x <- readTCX("activity_10318883291.tcx")
x1 <- filter(x, cadence_running > 65) %>%
  slice_tail(n = nrow(.) - 5) %>%
  slice_head(n = nrow(.) - 5) %>%
  arrange(time) %>%
  mutate(minutes = as.numeric((time - time[1]) / 60)) %>%
  mutate(n = 1:nrow(.)) %>%
  mutate(speed = as.numeric(cut(.$minutes, breaks = seq(-0.001, 14.999, by = 1.5), labels = 5:14)) + 5)
```

Fit change-point model

```{r}
library(SiZer)
(fit <- piecewise.linear(x=x1$speed, y=x1$heart_rate, CI = TRUE))
```
Plot HR vs time, with ALL HR measurements

```{r, message=FALSE}
fitdat <- data.frame(speed = fit$x, heart_rate = fit$y)
ggplot(x1, aes(x=speed, y=heart_rate)) +
  geom_point() +
#  geom_vline(xintercept=seq(1.5, 15, by = 1.5)) + 
  xlim(5, 15) + ylim(90, 200) +
  scale_x_continuous(breaks=seq(5, 15)) +
  scale_y_continuous(breaks=seq(90, 200, by=10)) +
  geom_line(data = fitdat, aes(x=speed, y = heart_rate), col = "red", linewidth = 1)
```

# Using only the average of the last 5 HR measurements in each step

```{r}
x2 <- group_by(x1, speed) %>%
  mutate(heart_rate = mean(tail(heart_rate, 5))) %>%
  select(speed, heart_rate) %>%
  unique()
x2
(fit2 <- piecewise.linear(x=x2$speed, y=x2$heart_rate, CI = TRUE))
```

```{r, message=FALSE}
fitdat2 <- data.frame(speed = fit2$x, heart_rate = fit2$y)
ggplot(x2, aes(x=speed, y=heart_rate)) +
  geom_point() +
  xlim(5, 15) + ylim(90, 200) +
  scale_x_continuous(breaks=seq(5, 15)) +
  scale_y_continuous(breaks=seq(90, 200, by=10)) +
  geom_line(data = fitdat2, aes(x=speed, y = heart_rate), col = "red", linewidth = 1)
```

# Using January 2022 data 

```{r}
x3 <- data.frame(speed = seq(6, 14),
                 heart_rate = c(131, 142, 153, 162, 170, 177, 183, 187, 192))
x3
```
```{r, warning=FALSE}
(fit3 <- piecewise.linear(x=x3$speed, y=x3$heart_rate, CI = TRUE))
```

```{r, message=FALSE}
fitdat3 <- data.frame(speed = fit3$x, heart_rate = fit3$y)
ggplot(x3, aes(x=speed, y=heart_rate)) +
  geom_point() +
  xlim(5, 15) + ylim(90, 200) +
  scale_x_continuous(breaks=seq(5, 15)) +
  scale_y_continuous(breaks=seq(90, 200, by=10)) +
  geom_line(data = fitdat3, aes(x=speed, y = heart_rate), col = "red", linewidth = 1)
```

# 2022-23 together

```{r}
x2$year = 2023L
x3$year = 2022L
yr <- factor(c(2022, 2023))
fitdat2$year = yr[2]
fitdat3$year = yr[1]
xall <- full_join(x=x2, y=x3) %>%
  mutate(year = factor(year))
ggplot(data=xall, mapping = aes(x=speed, y=heart_rate, color = year, pch = year)) +
  geom_point(size = 3) +
  xlim(5, 15) + ylim(90, 200) +
  scale_x_continuous(breaks=seq(5, 15)) +
  scale_y_continuous(breaks=seq(90, 200, by=10)) +
  geom_line(data = fitdat3, aes(x=speed, y = heart_rate, color =  year)) +
  geom_line(data = fitdat2, aes(x=speed, y = heart_rate, color =  year))
```