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
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.
morph$date <-mdy(morph$date) #convert the date column from characters to true datemorph <- 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 treatmentmean_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.
# To visualize the mean health scores using ggplot2ggplot(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 graphggplot(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 factorsmorph_fact$treatment_mg_per_L <-as.factor(morph_fact$treatment_mg_per_L)morph_fact$hpe <-as.factor(morph_fact$hpe)## Perform one-way ANOVAanova_result <-aov(health ~ treatment_mg_per_L + hpe + treatment_mg_per_L:hpe , data = morph_fact)# Summarize the ANOVA resultssummary(anova_result)## Meeting the ANOVA assumptionsleveneTest(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 testshapiro_test <-shapiro.test(residuals(anova_result))print(shapiro_test)#meets assumption of normality, p value > 0.05## Perform Tukey's HSD test for significant interactiontukey_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 comparisonsggplot(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 gridlinespanel.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/Lgeom_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/Lgeom_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 )