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

Load data and add speed and time variables

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

library(SiZer)
(fit <- piecewise.linear(x=x1$speed, y=x1$heart_rate, CI = TRUE))
## [1] "Threshold alpha: 12.2688986807056"
## [1] ""
## [1] "Model coefficients: Beta[0], Beta[1], Beta[2]"
## (Intercept)           x           w 
##   68.366210    8.766552   -5.779468 
##       Change.Point Initial.Slope Slope.Change Second.Slope
## 2.5%      11.70218      8.475224    -6.387856     2.399742
## 97.5%     12.51569      9.218629    -4.888003     3.953101

Plot HR vs time, with ALL HR measurements

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

x2 <- group_by(x1, speed) %>%
  mutate(heart_rate = mean(tail(heart_rate, 5))) %>%
  select(speed, heart_rate) %>%
  unique()
x2
## # A tibble: 10 × 2
## # Groups:   speed [10]
##    speed heart_rate
##    <dbl>      <dbl>
##  1     6       126.
##  2     7       133.
##  3     8       143 
##  4     9       148.
##  5    10       159 
##  6    11       167.
##  7    12       174.
##  8    13       179.
##  9    14       182.
## 10    15       185.
(fit2 <- piecewise.linear(x=x2$speed, y=x2$heart_rate, CI = TRUE))
## [1] "Threshold alpha: 12.3745900390539"
## [1] ""
## [1] "Model coefficients: Beta[0], Beta[1], Beta[2]"
## (Intercept)           x           w 
##   76.942904    8.142851   -5.542891 
##       Change.Point Initial.Slope Slope.Change Second.Slope
## 2.5%      9.089871      7.563235   -7.8333868    0.3369538
## 97.5%    13.005202      8.717926   -0.1769858    7.8646193
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

x3 <- data.frame(speed = seq(6, 14),
                 heart_rate = c(131, 142, 153, 162, 170, 177, 183, 187, 192))
x3
##   speed heart_rate
## 1     6        131
## 2     7        142
## 3     8        153
## 4     9        162
## 5    10        170
## 6    11        177
## 7    12        183
## 8    13        187
## 9    14        192
(fit3 <- piecewise.linear(x=x3$speed, y=x3$heart_rate, CI = TRUE))
## [1] "Threshold alpha: 9.59998916911483"
## [1] ""
## [1] "Model coefficients: Beta[0], Beta[1], Beta[2]"
## (Intercept)           x           w 
##   68.999918   10.400012   -5.000005 
##       Change.Point Initial.Slope Slope.Change Second.Slope
## 2.5%      8.499997      8.000007    -6.500003     3.999999
## 97.5%    11.353481     11.000024    -3.326153     7.000001
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

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))
## Joining, by = c("speed", "heart_rate", "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))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.