---
title: "Influence of wheel rim width on rolling resistance and off-road speed in cross-country mountain biking"
author: "Thomas Maier"
date: "February 2018"
output: html_document
---

Welcome to the data analysis code of the corresponding article.

## This document

This is an R Markdown document combining R code with comments.

To execute the code, change the file extension from ".txt" to ".Rmd" and open the file in R Studio <https://www.rstudio.com/>.

Before executing the code, install the necessary packages and place the file *rimwidth_data.csv* in the working directory. This CSV document contains the data for all Cr-trials.

Now press *Run All* or *Cmd/Ctrl + Alt + R*. You can also *knit* this document to an HTML report.

## Loading libraries

```{r setup, message=FALSE, warning = FALSE}
library(tidyverse)
library(forcats)
library(lme4) # mixed modelling
library(modelr)
library(pwr) # power analysis
```

## Power analysis

To find an effect of 1% compared to a typcial error of the measurement of around 3%.

```{r}
pwr.t.test(d = 0.33, sig.level = 0.05 , power = 0.80, type = "paired")
```

## Importing data

```{r message=FALSE}
data <- read_csv("rimwidth_data.csv") %>% 
  mutate(
    wheel_set = factor(wheel_set, levels = c("25baseline", "30samek", "30samep")),
    tyre = factor(tyre),
    participant = factor(participant),
    cr = (cr_inv1+cr_inv2+cr_inv3)/3 # mean cr of the 3 investigators
  )
```

## Custom plot layout

```{r}
theme_custom <- function (base_size = 8, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace% 
    theme(
      axis.ticks.length = unit(-0.1, "cm"),
      axis.ticks = element_line(colour = "black", size = 0.25),
      axis.text.x = element_text(margin = unit(c(0.2, 0.2, 0.0, 0.0), "cm"), 
                                 colour = "black"),
      axis.text.y = element_text(margin = unit(c(0.2, 0.2, 0.0, 0.0), "cm"), 
                                 colour = "black"),
      
      panel.background = element_rect(fill = "white", colour = NA), 
      panel.border = element_rect(fill = NA, colour = "black"), 
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(), 
      strip.background = element_rect(fill = "white", colour = "black"), 
      
      legend.key = element_rect(fill = "white", colour = NA),
      legend.key.width = unit(0.3, "cm"),
      legend.margin = margin(c(0, 0, 0.1, 0), unit = "cm"),
      legend.box.spacing = unit(0.0, "cm"),
      legend.position = "top",
      complete = TRUE)
}
```

## Exclude tests due to wind and suspicious 0-offset

```{r}
# Filter tests for wind
data_excl_wind  <-  data %>% filter(wind_max < 2)
  
# Filter tests for suspicious offset
data_excl_wind_offset <- data_excl_wind %>% 
  add_residuals(lm(offset ~ temperature, data_excl_wind)) %>% 
  filter(resid < 5, resid > -5) %>%
  select(-resid)

# Plot the linear model and excluded tests
ggplot() +
  aes(temperature, offset) +
  geom_point(data = setdiff(data_excl_wind, data_excl_wind_offset), shape = 4, size = 1.5) +
  geom_point(data = data_excl_wind_offset, size = 0.5) +
  geom_smooth(data = data_excl_wind, method = "lm", se = FALSE, color = "grey", size = 0.5) +
  labs(x = "Temperature (°C)", y = "0-Offset (Hz)") +
  theme_custom()

ggsave("figure3.svg", width = 8.4, height = 6, units = "cm")
ggsave("figure3.png", width = 8.4, height = 6, units = "cm", dpi = 1200)
```

## Differences in rolling resistance

### Mixed model with random intercept

```{r}
lmem <-     lmer(cr ~ wheel_set + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
baseline <- lmer(cr ~ 1 +         (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
```

```{r}
summary(lmem)
```
```{r}
anova(baseline, lmem)
```

```{r}
confint_cr <- confint(lmem, level = 0.90, oldNames = FALSE)
```

```{r}
confint_cr
```

### Plot differences by wheel set and participant

```{r}
data_excl_wind_offset %>%
  # compute % difference of each test to trial mean
  group_by(trial) %>% 
  mutate(cr_diff = (cr - mean(cr)) / mean(cr) * 100) %>% 
  ungroup() %>% 
  
ggplot() +
  aes(wheel_set, cr_diff, shape = participant) +
  geom_errorbar(stat = "summary",
                fun.ymin = function(x) mean_cl_normal(x, conf.int = 0.9)$ymin,
                fun.ymax = function(x) mean_cl_normal(x, conf.int = 0.9)$ymax,
                position = position_dodge(width = 0.5), width = 0.25, size = 0.25) +
  
  geom_point(stat = "summary", fun.y = mean, 
             position = position_dodge(width = 0.5), size = 1.5) +
  
  geom_hline(yintercept = 0, linetype = "dotted", size = 0.25) +
  
  scale_x_discrete(labels = c(expression("25"["baseline"]),
                              expression("30"["same"*italic("k")]),
                              expression("30"["same"*italic("p")]))) +
  
  scale_shape_discrete(name = "Rider") +
  labs(x = "Wheel set", y = expression(paste(italic(C)[r]," ", "difference to trial mean (%)"))) +
  coord_flip() +
  theme_custom()

ggsave("figure4.svg", width = 8.4, height = 6, units = "cm")
ggsave("figure4.png", width = 8.4, height = 6, units = "cm", dpi = 1200)
```


### Checking influence of temperature and wind

```{r}
ggplot(data_excl_wind_offset) +
  aes(temperature, cr, color = wind_max) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_custom()
```


```{r}
lmem_full <- lmer(cr ~ temperature + wind_max + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
lmem_temp <- lmer(cr ~ temperature + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
lmem_wind <- lmer(cr ~ wind_max + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
lmem_base <- lmer(cr ~ 1 + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
```

```{r}
anova(lmem_base, lmem_wind)
```
```{r}
anova(lmem_base, lmem_temp)
```

```{r}
anova(lmem_temp, lmem_full)
```

```{r}
summary(lmem_temp)
```

### Checking influence of speed

```{r}
data_excl_wind_offset %>% 
  group_by(participant) %>% 
  summarise(mean_speed = mean(avg_speed_kmh), 
            sd_speed = sd(avg_speed_kmh))
```

```{r}
lmem_speed <- lmer(cr ~ avg_speed_kmh + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
lmem_base <-  lmer(cr ~ 1 + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)

anova(lmem_base, lmem_speed)
```


```{r}
ggplot(data_excl_wind_offset) +
  aes(avg_speed_kmh, cr) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~participant, scales = "free")
```


## Investigator agreement

```{r}
cor(data_excl_wind_offset$cr_inv1, data_excl_wind_offset$cr_inv2)
cor(data_excl_wind_offset$cr_inv1, data_excl_wind_offset$cr_inv3)
cor(data_excl_wind_offset$cr_inv2, data_excl_wind_offset$cr_inv3)
```

```{r}
# limits of agreement
loa <- function(data1, data2){
  diff <- data1-data2
  return(c(mean(diff) + 1.96 * sd(diff), mean(diff) - 1.96 * sd(diff)))
}

loa(data_excl_wind_offset$cr_inv1, data_excl_wind_offset$cr_inv2)
loa(data_excl_wind_offset$cr_inv1, data_excl_wind_offset$cr_inv3)
loa(data_excl_wind_offset$cr_inv2, data_excl_wind_offset$cr_inv3)
```


## Off-road speed

### Mixed model for system weight

```{r}
lmem <- lmer(total_mass ~ wheel_set + (1|participant) + (1|trial) + (1|tyre), data = data_excl_wind_offset)
summary(lmem)
```

```{r warning = FALSE}
confint_m <- confint(lmem, level = 0.9, oldNames = FALSE)
```

```{r}
confint_m
```

### Model flat and uphill speed

```{r}
g = 9.81
CdA = 0.48
rho = 1.05
eta = 0.977

# Function calculates speed from Cr, m, s and P
speed = Vectorize(function(Cr, m, s, P){
  v = Re(polyroot(c(-P*eta, Cr*m*g + s/100*m*g, 0, 0.5*CdA*rho))[1])
  return(v)
})

# Ranges for Cr and m from confints of mixed models
cr_25 <- confint_cr[5,] %>% mean()
cr_ranges <- tibble(cr = c(cr_25, cr_25 + confint_cr[6,], cr_25 + confint_cr[7,]),
                    wheel_set = c("25baseline", "30samek", "30samek", "30samep", "30samep"),
                    limit = c("mean", "low", "high", "low", "high"))

m_25 <- confint_m[5,] %>% mean()
m_ranges  <-  tibble(wheel_set = c("25baseline", "30samek", "30samep"),
                     m = c(m_25, m_25+mean(confint_m[6,]), m_25+mean(confint_m[7,])))

slope <- tibble(slope = c(10,10,10,0,0,0), wheel_set = rep(c("25baseline", "30samek", "30samep"),2))
Power <- tibble(Power = c(400,400,400,300,300,300), wheel_set = rep(c("25baseline", "30samek", "30samep"),2))

all_combinations <- m_ranges %>% full_join(cr_ranges) %>% full_join(slope) %>% full_join(Power)

speed_data <- all_combinations %>% 
  mutate(speed = speed(cr, m, slope, Power)) %>% 
  group_by(slope, Power) %>% 
  mutate(diff_perz_to_25 = (speed / first(speed)-1)*100) %>% 
  arrange(slope, Power)

speed_data %>% 
  group_by(wheel_set) %>% 
  summarise(min = min(diff_perz_to_25),
            max = max(diff_perz_to_25))
```

Most likely trivial as the smallest worthwhile change SWC = 1.2%


```{r}
speed_data %>%
  ungroup() %>% 
  filter(wheel_set != "25baseline") %>% 
  mutate(wheel_set = factor(wheel_set, levels = c("30samek", "30samep"))) %>% 
  select(-cr, -speed) %>% 
  spread(limit, diff_perz_to_25) %>% 
  
ggplot() +
  aes(wheel_set, y = (low+high)/2, ymin = low, ymax = high, shape = factor(Power), color = factor(slope)) +
  geom_errorbar(position = position_dodge(0.5), width = 0.2, size = 0.25) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +

  geom_hline(yintercept = c(-1.2, 0, 1.2), linetype = "dotted", size = 0.25) +
  scale_y_continuous(breaks = seq(-1.8, 1.8, 0.6), 
                     labels = c(-1.8, "-SWC", -0.6, 0, 0.6, "SWC", 1.8)) +
  
  scale_x_discrete(labels = c(expression("30"["same"*italic("k")]),
                              expression("30"["same"*italic("p")])))  +
  
  scale_color_manual(name = "Slope (%)", values = c("darkgrey", "black")) +
  scale_shape_manual(name = "Power output (W)", values = c(16,15)) +
  labs(x = "Wheel set", y = expression(paste("Speed difference to"," ", "25"["baseline"]," ", "(%)"))) +
  coord_flip(ylim = c(-1.8, 1.8)) +
  theme_custom()
  
ggsave("figure5.svg", width = 8.4, height = 6, units = "cm")
ggsave("figure5.png", width = 8.4, height = 6, units = "cm", dpi = 1200)
```