catwalk_statistics.rmd

---
title: "R Notebook"
output: html_notebook
---

```{r}
library(tidyverse)
library(ez)
library(car)
library(broom)
library(readxl)
library(writexl)
library(lsmeans)
```


```{r, eval = F}
catstat_df <- raw_catwalk %>% 
  filter(gt != "1" & uid != 928211 & uid != 841350) %>% 
  mutate(age = factor(age, labels = c("YG", "AD", "OLD")),
         gt = factor(gt, labels = c("WT", "KO")),
         sex = factor(sex, labels = c("M", "F")))

catstat_df %>% 
  group_by(sex) %>% 
  count()
```

```{r}
catstat_df <- read.csv("catstat.csv")
catstat_df <- catstat_df %>% 
  mutate(age = factor(age, labels = c("YG", "AD", "OLD")),
         gt = factor(gt, labels = c("WT", "KO")),
         sex = factor(sex, labels = c("M", "F")))
  
```



```{r, eval = F}
catstat_df <- raw_catwalk %>% 
  filter(gt != "1" & uid != 928211 & uid != 841350) %>% 
  mutate(age = factor(age, labels = c("YG", "AD", "OLD")),
         gt = factor(gt, labels = c("WT", "KO")),
         sex = factor(sex, labels = c("M", "F")))

catstat_df %>% 
  group_by(sex) %>% 
  count()
```

```{r}
# have to set to helmert contrasts
options(contrasts = c("contr.helmert", "contr.poly"))
```

```{r}
# custom function to generate a tidy type 3 ANOVA
type_3_aov <- function(x, dataframe) {
  broom::tidy(car::Anova(lm(x ~ sex * age * gt, 
                            data=dataframe), 
                         type=3)) %>% 
    mutate(p.value = round(p.value, 3))
}

# custom function to generate a tukey post-hoc
post_hoc_lm <- function(x, dataframe){
  summary(lsmeans::lsmeans(lm(x ~ sex * age * gt, data = dataframe), 
                           pairwise ~ sex:age:gt, 
                           adjust = "tukey"))$contrasts %>% 
    mutate(p.value = round(p.value, 3)) %>% 
    filter(p.value < 0.051)
}

```

```{r, eval = F, echo = F}
# purr functions for calculating all things of interest at once

# models well
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ car::Anova(lm(.x ~ catstat_df$sex * catstat_df$age * catstat_df$gt))) %>% 
  map(~ broom::tidy(.)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5)))  %>% 
  map( ~ mutate(.x, sig = case_when(p.value <= 0.05 ~ "*",
                                p.value > 0.05 ~ " "))) %>% 
  write_xlsx(., "test.xlsx")

  
# tukey posthocs
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ lm(.x ~ sex * age * gt, data = catstat_df)) %>% 
  map(~ lsmeans::lsmeans(.x, pairwise ~ sex:age:gt, data = catstat_df, adjust = "tukey")$contrasts) %>% 
  map(~ tidy(.x)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5))) %>% 
  map(~ filter(.x, p.value < 0.051)) %>% 
  write_xlsx(., "test.xlsx")

# uncorrected lsd posthocs
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map( ~ pairwise.t.test(.x,
                  catstat_df$sex:catstat_df$age:catstat_df$gt, 
                  p.adj = "none")) %>% 
  map( ~ tidy(.x)) %>% 
  map( ~ mutate(.x, p.value = round(p.value, 5))) %>%
  map( ~ filter(.x, p.value < 0.051)) %>% 
  map( ~ select(.x, group2, group1, p.value)) %>% 
  map( ~ mutate(.x, sig = case_when(p.value <= 0.05 ~ "*",
                                p.value > 0.05 ~ " ")))

```


Speed ANOVA
```{r}
speed_aov <- aov(run_speed~sex*age*gt, data = catstat_df)
summary(speed_aov)
TukeyHSD(speed_aov, which = 'sex:age:gt')
# call and save the pair.t.test
speed_lsd <- tidy(pairwise.t.test(catstat_df$run_speed,
                             catstat_df$sex:catstat_df$age:catstat_df$gt, 
                             p.adj = "none")) %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
                
speed_lsd

# type 3 ANOVA and corresponding post-hoc
speed_aov_3 <- type_3_aov(catstat_df$run_speed, catstat_df)
speed_aov_3
speed_posthoc <- post_hoc_lm(catstat_df$run_speed, catstat_df)
```

front_bos ANOVA
```{r}
front_bos_aov <- aov(front_bos~age*gt, data = catstat_df)
summary(front_bos_aov)
TukeyHSD(front_bos_aov)

# call and save the pair.t.test
front_bos_lsd <- tidy(pairwise.t.test(catstat_df$front_bos,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
front_bos_lsd
```

hind_bos ANOVA
```{r}
hind_bos_aov <- aov(hind_bos~age*gt, data = catstat_df)
summary(hind_bos_aov)
TukeyHSD(hind_bos_aov)

# call and save the pair.t.test
hind_bos_lsd <- tidy(pairwise.t.test(catstat_df$hind_bos,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
hind_bos_lsd
```

front_bos ANOVA
```{r}
front_bos_aov <- aov(front_bos~age*gt, data = catstat_df)
summary(front_bos_aov)
TukeyHSD(front_bos_aov)

# call and save the pair.t.test
front_bos_lsd <- tidy(pairwise.t.test(catstat_df$front_bos,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
front_bos_lsd
```

front_strlen ANOVA
```{r}
front_strlen_aov <- aov(front_strlen~age*gt, data = catstat_df)
summary(front_strlen_aov)
TukeyHSD(front_strlen_aov)

# call and save the pair.t.test
front_strlen_lsd <- tidy(pairwise.t.test(catstat_df$front_strlen,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
front_strlen_lsd
```

hind_strlen ANOVA
```{r}
hind_strlen_aov <- aov(hind_strlen~age*gt, data = catstat_df)
summary(hind_strlen_aov)
TukeyHSD(hind_strlen_aov)

# call and save the pair.t.test
hind_strlen_lsd <- tidy(pairwise.t.test(catstat_df$hind_strlen,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
hind_strlen_lsd

```
front_swingspd ANOVA
```{r}
front_swingspd_aov <- aov(front_swingspd~age*gt, data = catstat_df)
summary(front_swingspd_aov)
TukeyHSD(front_swingspd_aov)

# call and save the pair.t.test
front_swingspd_lsd <- tidy(pairwise.t.test(catstat_df$front_swingspd,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
front_swingspd_lsd
```

hind_swingspd ANOVA
```{r}
hind_swingspd_aov <- aov(hind_swingspd~age*gt, data = catstat_df)
summary(hind_swingspd_aov)
TukeyHSD(hind_swingspd_aov)

# call and save the pair.t.test
hind_swingspd_lsd <- tidy(pairwise.t.test(catstat_df$hind_swingspd,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
hind_swingspd_lsd
```

front_stepcyc ANOVA
```{r}
front_stepcyc_aov <- aov(front_stepcyc~age*gt, data = catstat_df)
summary(front_stepcyc_aov)
TukeyHSD(front_stepcyc_aov)

# call and save the pair.t.test
front_stepcyc_lsd <- tidy(pairwise.t.test(catstat_df$front_stepcyc,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
front_stepcyc_lsd
```

hind_stepcyc ANOVA
```{r, eval = F, echo = F}
hind_stepcyc_aov <- aov(hind_stepcyc~age*gt, data = catstat_df)
summary(hind_stepcyc_aov)
TukeyHSD(hind_stepcyc_aov)

# call and save the pair.t.test
hind_stepcyc_lsd <- tidy(pairwise.t.test(catstat_df$hind_stepcyc,
                             catstat_df$age:catstat_df$gt, 
                             p.adj = "none"))  %>% 
              mutate(p.value = round(p.value, digits = 3)) %>% 
              filter(p.value < 0.051)
hind_stepcyc_lsd

```

```{r}
# calculate catwalk ANOVAs and write to Excel
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ car::Anova(lm(.x ~ sex * age * gt, data = catstat_df))) %>% 
  map(~ broom::tidy(.)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5)), .id = 'source') %>% 
  write_xlsx(., "catwalk_anovas.xlsx")

# calculate LSD post-hocs and write to excel
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map( ~ pairwise.t.test(.x,
                  catstat_df$sex:catstat_df$age:catstat_df$gt, 
                  p.adj = "none")) %>% 
  map( ~ tidy(.x)) %>% 
  map( ~ mutate(.x, p.value = round(p.value, 5))) %>%
  map( ~ filter(.x, p.value < 0.051)) %>% 
  write_xlsx(., "catwalk_lsds.xlsx")
```


```{r, echo = F, eval = F}
# purr functions for calculating all things of interest at once

# this works
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ car::Anova(lm(.x ~ catstat_df$sex * catstat_df$age * catstat_df$gt))) %>% 
  map(~ broom::tidy(.)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5)))

catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ car::Anova(lm(.x ~ sex * age * gt, data = catstat_df))) %>% 
  map(~ broom::tidy(.)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5)), .id = 'source') %>% 
  write_xlsx(., "catwalk_anovas.xlsx")

# this works!  
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ lm(.x ~ catstat_df$sex * catstat_df$age * catstat_df$gt)) %>% 
  map(~ lsmeans::lsmeans(.x, pairwise ~ catstat_df$sex:catstat_df$age:catstat_df$gt, adjust = "tukey")$contrasts) %>% 
  map(~ tidy(.x)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5))) %>% 
  map(~ filter(.x, p.value < 0.051))

catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ lm(.x ~ sex * age * gt, data = catstat_df)) %>% 
  map(~ lsmeans::lsmeans(.x, pairwise ~ sex:age:gt, adjust = "tukey")$contrasts) %>% 
  map(~ tidy(.x)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5))) %>% 
  map(~ filter(.x, p.value < 0.051))

catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map(~ lm(.x ~ sex * age * gt, data = catstat_df)) %>% 
  map(~ lsmeans::lsmeans(.x, pairwise ~ sex:age:gt, adjust = "tukey", data = catstat_df)$contrasts) %>% 
  map(~ tidy(.x)) %>% 
  map(~ mutate(.x, p.value = round(p.value, 5))) %>% 
  map(~ filter(.x, p.value < 0.051))

#this works!
catstat_df %>% 
  select(run_speed:hind_stepcyc) %>% 
  map( ~ pairwise.t.test(.x,
                  catstat_df$sex:catstat_df$age:catstat_df$gt, 
                  p.adj = "none")) %>% 
  map( ~ tidy(.x)) %>% 
  map( ~ mutate(.x, p.value = round(p.value, 5))) %>%
  map( ~ filter(.x, p.value < 0.051)) %>% 
  write_xlsx(., "catwalk_lsds.xlsx")

lsd_post$run_speed

```