library(knitr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidyr)
library(dplyr)
library(hrbrthemes)
library(ggplot2)
library(car)

Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

library(RColorBrewer)
library(ggpubr)
knitr::opts_chunk$set(echo = TRUE,
                      eval = FALSE)

Objective and Background

Data below was collected from a pilot study conducted in Autumn 2023 where field-collected individual Botryllus schlosseri were exposed to a low and high treatment of nickel chloride and then were assessed for changes in morphology. Refer to Pilot 1.1 Nickel Exposure notebook post for experimental design details.

Retrieving Data from Google Sheets

Make sure you have made your Google sheet publicly available to anyone that has the link. If you make any updates to the sheet just re-curl the data, meaning just re-run the code below.

curl -L https://docs.google.com/spreadsheets/d/10uM3N3PD9xIP4yUnadfhkcXa8TPMRMD-adOIKbmkYzY/export?exportFormat=csv | tee ~/Documents/GitHub/Git/bestblogever/posts/morph-nicl2-botryllus-AU23/data/morph.csv

Read in the data to your local R environment.

morph <- read.csv(file = "data/morph.csv")

Cleaning up Data

morph$date <- mdy(morph$date) #convert the date column from characters to true date

morph <- morph %>%
  separate(jar_id, c("treatment_mg_per_L", "replicate"), sep = "-") #create two new columns, treatment and replicate from jar id 

morph_fact <- morph %>%
  mutate(treatment_mg_per_L = as.factor(treatment_mg_per_L)) %>%
  mutate(stage = as.factor(stage)) %>%
  mutate(animal_id = as.factor(animal_id)) %>%
  mutate(date = as.factor(date)) %>%
  mutate(treatment_order = factor(paste(treatment_mg_per_L, animal_id))) # Create a new variable for ordering by treatment

summary(morph_fact)

# Calculate the mean health score binned by hours post-exposure and by treatment
mean_health_score <- morph_fact %>%
  group_by(hpe, treatment_mg_per_L) %>%
  summarise(mean_health = round(mean(health, na.rm = TRUE)))

`summarise()` has grouped output by 'hpe'. You can override using the `.groups`
argument.

# Calculate SD
SD_health_score <- morph_fact %>%
  group_by(hpe, treatment_mg_per_L) %>%
  summarise(SD_health = round(sd(health, na.rm = TRUE)))

`summarise()` has grouped output by 'hpe'. You can override using the `.groups`
argument.

n <- 14

SE_health_score <- SD_health_score$SD_health/ sqrt(n)
# Print the calculated mean health scores
print(mean_health_score)

# A tibble: 6 × 3
# Groups:   hpe [2]
    hpe treatment_mg_per_L mean_health
  <int> <fct>                    <dbl>
1     0 0                            8
2     0 05                           8
3     0 45                           9
4    24 0                            8
5    24 05                           6
6    24 45                           6

print(SE_health_score)

[1] 0.2672612 0.2672612 0.2672612 0.2672612 0.2672612 0.5345225

Exploratory Plots

# To visualize the mean health scores using ggplot2
ggplot(mean_health_score, aes(x = hpe, y = mean_health, color = treatment_mg_per_L)) +
  geom_line() +
  geom_point() +
  labs(title = "Mean Health Score by Days Post-Exposure and Treatment",
       x = "Hours Post-Exposure",
       y = "Mean Health Score",
       color = "Concentration of Nickel Chloride (mg/L)") +
  theme_minimal()

Just the plot alone with no added significance bars and not adjusted for printing.

# Create the clustered bar graph
ggplot(mean_health_score, aes(x = factor(hpe), y = mean_health, fill = as.factor(treatment_mg_per_L))) +
  geom_bar(stat = "identity", position = "dodge", width = 0.8) +
scale_fill_brewer(palette = "Green", name = "Nickel Chloride (mg/L)", 
                    labels = c("0" = "0", "05" = "5", "45" = "45")) +
  labs(title = "Health Score by Hours Post-Exposure and Treatment",
       x = "Hours Post-Exposure",
       y = "Health Score") +
  theme_minimal() +
  scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10), limits = c(0, 10))

Statistically evaluate differences in health score

## Convert Treatment and HoursPostExposure to factors if they are not already factors
morph_fact$treatment_mg_per_L <- as.factor(morph_fact$treatment_mg_per_L)
morph_fact$hpe <- as.factor(morph_fact$hpe)


## Perform one-way ANOVA
anova_result <- aov(health ~ treatment_mg_per_L + hpe + treatment_mg_per_L:hpe , data = morph_fact)

# Summarize the ANOVA results
summary(anova_result)


## Meeting the ANOVA assumptions
leveneTest(health ~ treatment_mg_per_L*hpe, data = morph_fact)
#p-value > 0.05, therefore equal-variances are met

# Check normality of residuals using Shapiro-Wilk test
shapiro_test <- shapiro.test(residuals(anova_result))
print(shapiro_test)
#meets assumption of normality, p value > 0.05

## Perform Tukey's HSD test for significant interaction

tukey_result <- TukeyHSD(x= anova_result, ordered = TRUE, conf.level = 0.95)
print(tukey_result)

Printing out a png image of graph. Make sure to adjust pixel dimensions for the ideal DPI and desired printing size.

Expository Graphs

#png(filename = "bar_plot.png", width = 800, height = 800) # good for making a 3x3 in image for printing ~150 dpi

# Add significance bars based on specific comparisons
ggplot(mean_health_score, aes(x = factor(hpe), y = mean_health, fill = as.factor(treatment_mg_per_L))) +
  geom_bar(stat = "identity", position = "dodge", width = 0.8) +
scale_fill_brewer(palette = "Dark2", name = "Treatment", 
                    labels = c("0" = "Control", "05" = "5 mg/L", "45" = "45 mg/L")) +
  labs(
       x = "Hours Post-Exposure",
       y = "Mean Health Score", ) +
  theme_minimal() +
  scale_y_continuous(breaks = c(0, 5, 10), limits = c(0, 10)) +
  geom_errorbar(aes(ymin=mean_health_score$mean_health-SE_health_score, ymax=mean_health_score$mean_health+SE_health_score), width=.2,
                 position=position_dodge(0.81)) +
  theme(panel.grid.major.x = element_blank(),  # Remove major vertical gridlines
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(color = "darkgrey", size = 1),
        panel.grid.minor.y = element_line(color = "darkgrey", size = 0.5), 
        axis.text=element_text(size=30),
        axis.title=element_text(size=40),
        legend.text = element_text(size = 40),
        legend.title = element_text(size = 40)) + # Adjust the font size of the legend text 
  geom_text(aes(x = 1.98, y = 9.25, label = "**"), size = 14, vjust = -0.5) +  # Asterisk for significant difference at 24 hours for Control vs. 45 mg/L
    geom_segment(aes(x = 1.65, xend = 2.27, y = 9.5, yend = 9.5), linetype = "dashed", size = 1 ) +  # Line below the asterisk for Control vs. 45 mg/L
  geom_segment(aes(x = 2.27, xend = 2.27, y = 6.7, yend = 9.5), linetype = "dashed", size = 1 ) +
  geom_segment(aes(x = 1.65, xend = 1.65, y = 8.4, yend = 9.5), linetype = "dashed", size = 1 ) +
  geom_text(aes(x = 1.9, y = 8.25, label = "**"), size = 14, vjust = -0.5) +   #asterisk for significant difference at 24 hours for control vs 5 mg/L 
  geom_segment(aes(x = 1.8, xend = 2, y = 8.5, yend = 8.5), linetype = "dashed", size = 1 ) +
  geom_segment(aes(x = 1.8, xend = 1.8, y = 8.4, yend = 8.5), linetype = "dashed", size = 1 ) +
  geom_segment(aes(x = 2, xend = 2, y = 6.5, yend = 8.5), linetype = "dashed", size = 1 )