Analysis of “Ground reaction force metrics are not strongly correlated with tibial bone load when running across speeds and slopes: Implications for science, sport and wearable tech”

data_from_matlab_gui.csv

GRF_Tibia_Load.R

library(ggplot2)
library(dplyr)
library(Metrics)

## Import discrete variables ----
data <- read.csv('data_from_matlab_gui.csv', header = TRUE) # csv compiled from GUI data
data$sub <- as.factor(data$sub)

## Plot Correlation across Slopes ----
ggplot(data = data, 
       aes(x = GRF_active, 
           y = tibia_max,
           fill= sub, 
           color = sub, 
           group = sub))+
  geom_smooth(method = 'lm', 
              se = F) +
  geom_smooth(aes(group = as.factor(incline)), 
              method = 'lm', 
              se = F, 
              color = 'black') +
  geom_point() +
  theme_classic() +
  ggtitle('Peak TBL ~ Peak GRF across Speeds & Slopes')+
  facet_grid( ~incline)

## Calculate Correlations ----
library(plyr)  # Throws warning about conflict with dplyr. Detach plyr after use.
corfun <- function(x) {
  return(data.frame(COR = cor(x$GRF_active, x$tibia_max)))
}

# Not grouped by subject:
ddply(data, .(incline), corfun)

# Grouped by subject:
correlations <- matrix(nrow=7, ncol = 3)
all_correlations <- matrix(nrow=0, ncol = 3)

for (i in 1:10){
  
  r_values <- ddply(data[data$sub == i,], .(incline), corfun)$COR  # by incline
  incline <- ddply(data[data$sub == i,], .(incline), corfun)$incline  # by incline
  
  correlations[,1] = i
  correlations[,2] = incline
  correlations[,3] = r_values
  
  all_correlations <- rbind(all_correlations, correlations)
}

all_correlations <- data.frame(all_correlations)
colnames(all_correlations) <- c('sub', 'incline', 'r')
detach(package:plyr)

all_correlations %>% group_by(incline) %>% summarize(avg_r = mean(r), sd_r = sd(r))

## Fit a Linear Regression ----
# fit model on subjects 1 - 9, test on subject 10
formula_lr <- tibia_max ~ GRF_active + speed + incline
model <- lm(formula_lr, data = data[data$sub != 10,], na.action = na.omit)
summary(model)

#predict data for subject 10
data$pred_peak[data$sub == 10] <- predict(model, newdata = data[data$sub == 10,], level = 0)

# RMSE 
rmse(data$pred_peak[data$sub == 10], data$tibia_max[data$sub == 10])
# 0.370 BW

# plot prediction vs measured peak tibia bone load
ggplot(data[data$sub == 10,], aes(x = pred_peak, y = tibia_max)) +
  geom_point(aes(color = as.factor(incline)))+
  geom_abline() +
  theme_classic()+
  coord_fixed(xlim = c(5, 9.5), ylim = c(5, 9.5))