---
title: "2_analyze_thread_strength"
author: "Matt George"
date: "2023-04-25"
output: html_document
---

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

### LDOd packages; set ggplot theme
```{r message=FALSE, warning=FALSE}
## clear workspace
rm(list=ls())

## LDOd Packages
lDOd.lib<-c("readxl","tidyverse","bestNormalize","agricolae","nlme","multcomp","rstatix","lme4","lmerTest") # List of required packages
install.lib <- lDOd.lib[!lDOd.lib %in% installed.packages()] # Select missing packages
for(lib in install.lib) install.packages(lib,dependencies=TRUE) # Install missing packages + dependencies
sapply(lDOd.lib,require,character=TRUE) # LDOd all packages.

## Set ggplot theme
my_theme <- theme(line              = element_line(size=1.2),
                  rect              = element_rect(size=1.2),
                  text              = element_text(size=14,color="black"),
                  panel.background  = element_blank(),
                  panel.grid.major  = element_blank(), 
                  panel.grid.minor  = element_blank(),
                  axis.text.x       = element_text(size=16,color="black"), #,angle=90),
                  axis.text.y       = element_text(size=16,color="black"),
                  axis.title.x      = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
                  axis.title.y      = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
                  axis.ticks.x      = element_line(color="black"),
                  axis.ticks.y      = element_line(color="black"),
                  # axis.line         = element_line(color = "black", size = 0.1),
                  panel.border      = element_rect(color = "black", fill=NA, size=1.2),
                  legend.key        = element_blank(), # removes background of legend bullets
                  legend.position   = "none"
                  ) 
```

```{bash}
mkdir analyses
mkdir analyses/thread_strength

```


### LOAD data set
```{r}
# Load datasets
thread_plot       <- read_excel("thread_strength/summarized_data/summary.xlsx", 
                              sheet = "data", col_names = TRUE)
# trt_list        <- read_excel("data/respirometry/respirometry.xlsx", 
#                               sheet = "trt_list", col_names = TRUE)
thread_plot$treatment  <- factor(thread_plot$treatment,  levels = c("control","OW","OA","DO"), ordered=TRUE)
thread_plot$mussel     <- factor(thread_plot$mussel)
thread_plot$thread     <- factor(thread_plot$thread)

# # set alpha for boxplots
# alpha_set = 0.9
# 
# # summarize by ploidy * timepoint * treatment
# resp_plot %>% group_by(trt_list) %>% summarise(mean=mean(SMR_weight), sd=sd(SMR_weight), count=n())
# 
# # define SE
# summary_se <- function(.data, measure_var, ..., .ci = 0.95, na.rm = FALSE) {
#   
#   measure_var <- dplyr::enquo(measure_var)
#   group_var <- dplyr::enquos(...)
#   
#   .data %>%
#     group_by(!!! group_var) %>%
#     summarise(mean = mean(!! measure_var, na.rm = na.rm),
#               sd = sd(!! measure_var, na.rm = na.rm),
#               n = n(),
#               se = sd/sqrt(n),
#               ci = se * qt(.ci/2 + 0.5, n-1)) %>%
#     ungroup()
#   
# }

```



### boxplot, all
```{r}

# filter the data to only include the mussels w/ repeated observations
thread_plot_repeated <- thread_plot %>%
  group_by(mussel) %>%
  filter(n_distinct(treatment) == 2)

# create a new dataframe that groups the data by mussel and treatment and then summarizes the data by calculating the mean and standard deviation of the max_force, pad_area, and adhesion_kpa columns
thread_plot_all_mean <- thread_plot %>%
  group_by(mussel, treatment) %>% 
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))

# average across thread for each mussel * treatment combination, only mussels with repeated measurements
# create a new dataframe called thread_plot_repeated_mean
thread_plot_repeated_mean <- thread_plot %>%
  # group the data by mussel
  group_by(mussel) %>%
  # filter the data to only include mussels that have two treatments
  filter(n_distinct(treatment) == 2) %>%
  # group the data by mussel and treatment
  group_by(mussel, treatment) %>% 
  # summarize the data by calculating the mean and sd of max_force, pad_area, and adhesion_kpa
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))



# Plot all data
p1 <- ggplot(thread_plot, aes(x=treatment, y=adhesion_kpa)) +
       geom_boxplot(colour = "black", size = 0.8,outlier.colour="black", outlier.shape = 16,
                           outlier.size=1, notch=FALSE, varwidth = TRUE, alpha = 0.8) +
      scale_y_continuous(breaks = seq(0, 200, 25), limits = c(0, 210)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p1


# plot only mussels with repeated measurements
p2 <- ggplot(thread_plot_repeated, aes(x=treatment, y=adhesion_kpa)) +
       geom_boxplot(colour = "black", size = 0.8,outlier.colour="black", outlier.shape = 16,
                           outlier.size=1, notch=FALSE, varwidth = TRUE, alpha = 0.8) +
      scale_y_continuous(breaks = seq(0, 200, 25), limits = c(0, 210)) +
      xlab(NULL) + ylab(NULL) + my_theme

p2


# plot data from all mussels, averaging across threads produced by the mussel at each timepoint
p3 <- ggplot(thread_plot_all_mean, aes(x=treatment, y=adhesion_kpa_mean)) +
       geom_boxplot(colour = "black", size = 0.8,outlier.colour="black", outlier.shape = 16,
                           outlier.size=1, notch=FALSE, varwidth = TRUE, alpha = 0.8) +
      scale_y_continuous(breaks = seq(0, 200, 25), limits = c(0, 210)) +
      xlab(NULL) + ylab(NULL) + my_theme

p3

# plot only mussels with repeated measures, averaging across threads produced by the mussel at each timepoint
p4 <- ggplot(thread_plot_repeated_mean, aes(x=treatment, y=adhesion_kpa_mean)) +
       geom_boxplot(colour = "black", size = 0.8,outlier.colour="black", outlier.shape = 16,
                           outlier.size=1, notch=FALSE, varwidth = TRUE, alpha = 0.8) +
      scale_y_continuous(breaks = seq(0, 200, 25), limits = c(0, 210)) +
      xlab(NULL) + ylab(NULL) + my_theme

p4


ggsave("analyses/thread_strength/BP_AllMussels_AllThreads.png",
       plot   = p1,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/BP_OnlyRepeatedMussels_AllThreads.png",
       plot   = p2,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/BP_AllMussels_MeanThreads.png",
       plot   = p3,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/BP_OnlyRepeatedMussels_MeanThreads.png",
       plot   = p4,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

```

### linegraph, DO only
```{r}

# ChatGPT prompts:
# create a new dataframe called "thread_plot_DO" that is only data from thread_plot when the column "treatment" is control or DO and the same "mussel" ID is represented in both treatments (as noted in the "mussel" column)
# create a new dataframe called "thread_plot_DO_mean" from "thread_plot" that produces the mean and standard deviation of the columns "max_force", "pad_area", and "adhesion_kpa"  across unique combinations of "mussel" and "treatment".


# filter the data to only include the control and DO treatment
thread_plot_DO_all <- thread_plot %>% filter(treatment == "control" | treatment == "DO")

# filter the data to only include the control and DO treatments
thread_plot_DO_repeated <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "DO")) %>%
  filter(n_distinct(treatment) == 2)

# create a new data frame called thread_plot_DO_mean
# group the data by mussel
# filter the data to only include the control and DO treatments
# group the data by mussel and treatment
# summarize the data by calculating the mean and sd for max_force, pad_area, and adhesion_kpa
# rename the columns to include the statistic calculated
thread_plot_DO_mean <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "DO")) %>%
  group_by(mussel, treatment) %>%
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))

thread_plot_DO_repeated_mean <- thread_plot %>%
  # group the data by mussel
  group_by(mussel) %>%
  # filter the data to only include the control and DO treatments
  filter(treatment %in% c("control", "DO")) %>%
  # filter the data to only include mussels that have two treatments
  filter(n_distinct(treatment) == 2) %>%
  # group the data by mussel and treatment
  group_by(mussel, treatment) %>% 
  # summarize the data by calculating the mean and sd of max_force, pad_area, and adhesion_kpa
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))



p1 <- ggplot(thread_plot_DO_all, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p1

p2 <- ggplot(thread_plot_DO_repeated, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p2

p3 <- ggplot(thread_plot_DO_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p3

p4 <- ggplot(thread_plot_DO_repeated_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p4



ggsave("analyses/thread_strength/LineGraph_DO_AllMussels_MeanThreads.png",
       plot   = p3,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/LineGraph_DO_OnlyRepeatedMussels_MeanThreads.png",
       plot   = p4,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")


```

### linegraph, OW only
```{r}

# ChatGPT prompts:
# create a new dataframe called "thread_plot_OW" that is only data from thread_plot when the column "treatment" is control or OW and the same "mussel" ID is represented in both treatments (as noted in the "mussel" column)
# create a new dataframe called "thread_plot_OW_mean" from "thread_plot" that produces the mean and standard deviation of the columns "max_force", "pad_area", and "adhesion_kpa"  across unique combinations of "mussel" and "treatment".


# filter the data to only include the control and OW treatment
thread_plot_OW_all <- thread_plot %>% filter(treatment == "control" | treatment == "OW")

# filter the data to only include the control and OW treatments
thread_plot_OW_repeated <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "OW")) %>%
  filter(n_distinct(treatment) == 2)

# create a new data frame called thread_plot_OW_mean
# group the data by mussel
# filter the data to only include the control and OW treatments
# group the data by mussel and treatment
# summarize the data by calculating the mean and sd for max_force, pad_area, and adhesion_kpa
# rename the columns to include the statistic calculated
thread_plot_OW_mean <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "OW")) %>%
  group_by(mussel, treatment) %>%
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))

thread_plot_OW_repeated_mean <- thread_plot %>%
  # group the data by mussel
  group_by(mussel) %>%
  # filter the data to only include the control and OW treatments
  filter(treatment %in% c("control", "OW")) %>%
  # filter the data to only include mussels that have two treatments
  filter(n_distinct(treatment) == 2) %>%
  # group the data by mussel and treatment
  group_by(mussel, treatment) %>% 
  # summarize the data by calculating the mean and sd of max_force, pad_area, and adhesion_kpa
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))



p1 <- ggplot(thread_plot_OW_all, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p1

p2 <- ggplot(thread_plot_OW_repeated, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p2

p3 <- ggplot(thread_plot_OW_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p3

p4 <- ggplot(thread_plot_OW_repeated_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p4



ggsave("analyses/thread_strength/LineGraph_OW_AllMussels_MeanThreads.png",
       plot   = p3,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/LineGraph_OW_OnlyRepeatedMussels_MeanThreads.png",
       plot   = p4,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")


```
### linegraph, DO only
```{r}

# ChatGPT prompts:
# create a new dataframe called "thread_plot_DO" that is only data from thread_plot when the column "treatment" is control or DO and the same "mussel" ID is represented in both treatments (as noted in the "mussel" column)
# create a new dataframe called "thread_plot_DO_mean" from "thread_plot" that produces the mean and standard deviation of the columns "max_force", "pad_area", and "adhesion_kpa"  across unique combinations of "mussel" and "treatment".


# filter the data to only include the control and DO treatment
thread_plot_DO_all <- thread_plot %>% filter(treatment == "control" | treatment == "DO")

# filter the data to only include the control and DO treatments
thread_plot_DO_repeated <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "DO")) %>%
  filter(n_distinct(treatment) == 2)

# create a new data frame called thread_plot_DO_mean
# group the data by mussel
# filter the data to only include the control and DO treatments
# group the data by mussel and treatment
# summarize the data by calculating the mean and sd for max_force, pad_area, and adhesion_kpa
# rename the columns to include the statistic calculated
thread_plot_DO_mean <- thread_plot %>%
  group_by(mussel) %>%
  filter(treatment %in% c("control", "DO")) %>%
  group_by(mussel, treatment) %>%
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))

thread_plot_DO_repeated_mean <- thread_plot %>%
  # group the data by mussel
  group_by(mussel) %>%
  # filter the data to only include the control and DO treatments
  filter(treatment %in% c("control", "DO")) %>%
  # filter the data to only include mussels that have two treatments
  filter(n_distinct(treatment) == 2) %>%
  # group the data by mussel and treatment
  group_by(mussel, treatment) %>% 
  # summarize the data by calculating the mean and sd of max_force, pad_area, and adhesion_kpa
  summarize(across(c(max_force, pad_area, adhesion_kpa), list(mean = mean, sd = sd), .names = "{.col}_{.fn}"))



p1 <- ggplot(thread_plot_DO_all, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p1

p2 <- ggplot(thread_plot_DO_repeated, aes(x=treatment, y=adhesion_kpa, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p2

p3 <- ggplot(thread_plot_DO_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p3

p4 <- ggplot(thread_plot_DO_repeated_mean, aes(x=treatment, y=adhesion_kpa_mean, group=as.factor(mussel))) +
      geom_point() +
      geom_line(linewidth=1, alpha = 1) +
      # facet_wrap(~ID) +
      # scale_color_manual(values=c("royalblue4","royalblue1")) +
      scale_y_continuous(breaks = seq(0, 150, 50), limits = c(0, 155)) +
      # scale_x_continuous(breaks = seq(0, 30, 5), limits = c(0, 32)) +
      xlab(NULL) + ylab(NULL) + my_theme

p4



ggsave("analyses/thread_strength/LineGraph_DO_AllMussels_MeanThreads.png",
       plot   = p3,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")

ggsave("analyses/thread_strength/LineGraph_DO_OnlyRepeatedMussels_MeanThreads.png",
       plot   = p4,
       dpi    = 1200,
       device = "png",
       width  = 4,
       height = 4,
       units  = "in")


```


### Statistical Testing - all
```{R}
# Define dataset for analysis, factors
dat_stat <- thread_plot_repeated_mean
test_me  <- dat_stat$adhesion_kpa_mean

dat_stat$ID         <- factor(dat_stat$mussel)
dat_stat$trt        <- factor(dat_stat$treatment)
# dat_stat$thread     <- factor(dat_stat$thread)
# dat_stat$pad        <- factor(dat_stat$pad_area)

# Test for normality
qqnorm(test_me) # check linearity of transformed data
qqline(test_me)
norm_test <- shapiro.test(test_me) # p-value fail = good, don't need transformation
print(norm_test$p.value)

# Normalize data if normality test failed
if(norm_test$p.value<0.05)     {
        normalized <- bestNormalize(test_me)
        test_me <- normalized$x.t # overwrite
        qqnorm(test_me) # check linearity of transformed data
        qqline(test_me)
        print(shapiro.test(test_me))
        print("transformed!",quote=FALSE)}

dat_stat$response <- test_me

# Run repeated measures anova using oyster ID has random effect
# https://www.r-bloggers.com/2015/08/two-way-anova-with-repeated-measures/
my_test <- aov(response ~ trt + Error(ID), # Error(1|(factor(ID)/factor(trt))),
                       data = dat_stat)

my_test_summary <- summary(my_test)
summary(my_test)

# Significant interaction observed between timepoint:trt:ploidy. Proceed with HSD test.
tx <- with(dat_stat, interaction(trt))
amod <- aov(response ~ tx, data = dat_stat)
mult_comp <- HSD.test(amod, "tx", group=TRUE, console=TRUE)
# 
# # save output
# write.table(my_test_summary$`Error: Within`[[1]], file = "stats/respirometry/AOV_respirometry.csv",      row.names = TRUE)
# write.table(mult_comp$groups,                     file = "stats/respirometry/HSD_test_respirometry.csv", row.names = TRUE)


```



### Statistical Testing - OA
```{R}
# Define dataset for analysis, factors
dat_stat <- thread_plot_OA_all
test_me  <- dat_stat$adhesion_kpa

dat_stat$ID         <- factor(dat_stat$mussel)
dat_stat$trt        <- factor(dat_stat$treatment)
dat_stat$thread     <- factor(dat_stat$thread)
dat_stat$pad        <- factor(dat_stat$pad_area)

# Test for normality
qqnorm(test_me) # check linearity of transformed data
qqline(test_me)
norm_test <- shapiro.test(test_me) # p-value fail = good, OAn't need transformation
print(norm_test$p.value)

# Normalize data if normality test failed
if(norm_test$p.value<0.05)     {
        normalized <- bestNormalize(test_me)
        test_me <- normalized$x.t # overwrite
        qqnorm(test_me) # check linearity of transformed data
        qqline(test_me)
        print(shapiro.test(test_me))
        print("transformed!",quote=FALSE)}

dat_stat$response <- test_me


# ChatGPT prompt:
# All data is in the "thread_plot_OA_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov funtion.# All data is in the "thread_plot_OA_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov function.


# This will generate an ANOVA model with "treatment" as the fixed effect and "mussel" and "thread" as ranOAm grouping factors.
# The "Error" term specifies that you want a within-subjects design, with the "mussel" and "thread" factors nested within each other.
my_test <- aov(response ~ trt + Error(ID/thread), data = dat_stat)

my_test_summary <- summary(my_test)
summary(my_test)

# Significant interaction observed between timepoint:trt:ploidy. Proceed with HSD test.
tx <- with(dat_stat, interaction(trt))
amod <- aov(response ~ tx, data = dat_stat)
mult_comp <- HSD.test(amod, "tx", group=TRUE, console=TRUE)
# 
# # save output
# write.table(my_test_summary$`Error: Within`[[1]], file = "stats/respirometry/AOV_respirometry.csv",      row.names = TRUE)
# write.table(mult_comp$groups,                     file = "stats/respirometry/HSD_test_respirometry.csv", row.names = TRUE)

```


### Statistical Testing - OW
```{R}
# Define dataset for analysis, factors
dat_stat <- thread_plot_OW_all
test_me  <- dat_stat$adhesion_kpa

dat_stat$ID         <- factor(dat_stat$mussel)
dat_stat$trt        <- factor(dat_stat$treatment)
dat_stat$thread     <- factor(dat_stat$thread)
dat_stat$pad        <- factor(dat_stat$pad_area)

# Test for normality
qqnorm(test_me) # check linearity of transformed data
qqline(test_me)
norm_test <- shapiro.test(test_me) # p-value fail = good, don't need transformation
print(norm_test$p.value)

# Normalize data if normality test failed
if(norm_test$p.value<0.05)     {
        normalized <- bestNormalize(test_me)
        test_me <- normalized$x.t # overwrite
        qqnorm(test_me) # check linearity of transformed data
        qqline(test_me)
        print(shapiro.test(test_me))
        print("transformed!",quote=FALSE)}

dat_stat$response <- test_me


# ChatGPT prompt:
# All data is in the "thread_plot_OW_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov funtion.# All data is in the "thread_plot_OW_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov function.


# This will generate an ANOVA model with "treatment" as the fixed effect and "mussel" and "thread" as random grouping factors.
# The "Error" term specifies that you want a within-subjects design, with the "mussel" and "thread" factors nested within each other.
my_test <- aov(response ~ trt + Error(ID/thread), data = dat_stat)

my_test_summary <- summary(my_test)
summary(my_test)

# Significant interaction observed between timepoint:trt:ploidy. Proceed with HSD test.
tx <- with(dat_stat, interaction(trt))
amod <- aov(response ~ tx, data = dat_stat)
mult_comp <- HSD.test(amod, "tx", group=TRUE, console=TRUE)
# 
# # save output
# write.table(my_test_summary$`Error: Within`[[1]], file = "stats/respirometry/AOV_respirometry.csv",      row.names = TRUE)
# write.table(mult_comp$groups,                     file = "stats/respirometry/HSD_test_respirometry.csv", row.names = TRUE)

```

### Statistical Testing - DO
```{R}
# Define dataset for analysis, factors
dat_stat <- thread_plot_DO_all
test_me  <- dat_stat$adhesion_kpa

dat_stat$ID         <- factor(dat_stat$mussel)
dat_stat$trt        <- factor(dat_stat$treatment)
dat_stat$thread     <- factor(dat_stat$thread)
dat_stat$pad        <- factor(dat_stat$pad_area)

# Test for normality
qqnorm(test_me) # check linearity of transformed data
qqline(test_me)
norm_test <- shapiro.test(test_me) # p-value fail = good, don't need transformation
print(norm_test$p.value)

# Normalize data if normality test failed
if(norm_test$p.value<0.05)     {
        normalized <- bestNormalize(test_me)
        test_me <- normalized$x.t # overwrite
        qqnorm(test_me) # check linearity of transformed data
        qqline(test_me)
        print(shapiro.test(test_me))
        print("transformed!",quote=FALSE)}

dat_stat$response <- test_me


# ChatGPT prompt:
# All data is in the "thread_plot_DO_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov funtion.# All data is in the "thread_plot_DO_all" dataframe. Multiple mussels were repeatedly sampled. Mussels were given a unique ID, reported in the "mussel" column. Each Mussel produced multiple threads before and after exposure to a "treatment". Each thread is unique and given a label in the "thread" column. The strength of threads were tested and reported as "adhesion_kpa". Generate a repeated measures anova that tests the impact of "treatment" on "adhesion_kpa". Use the aov function.


# This will generate an ANOVA model with "treatment" as the fixed effect and "mussel" and "thread" as random grouping factors.
# The "Error" term specifies that you want a within-subjects design, with the "mussel" and "thread" factors nested within each other.
my_test <- aov(response ~ trt + Error(ID/thread), data = dat_stat)

my_test_summary <- summary(my_test)
summary(my_test)

# Significant interaction observed between timepoint:trt:ploidy. Proceed with HSD test.
tx <- with(dat_stat, interaction(trt))
amod <- aov(response ~ tx, data = dat_stat)
mult_comp <- HSD.test(amod, "tx", group=TRUE, console=TRUE)
# 
# # save output
# write.table(my_test_summary$`Error: Within`[[1]], file = "stats/respirometry/AOV_respirometry.csv",      row.names = TRUE)
# write.table(mult_comp$groups,                     file = "stats/respirometry/HSD_test_respirometry.csv", row.names = TRUE)

```
