#Rcode belonging to Smulders et al. 2022. Data underlying publication: Green turtles shape the seascape #through grazing patch formation around habitat features: experimental evidence # Clean workspace rm(list=ls()) # Load libraries library(readr) library(dplyr) library(ggplot2) library(plyr) library(lattice) library(reshape) library(tidyverse) library(dbplyr) library(ggpubr) library(car) library(lawstat) library(tidyverse) library(rstatix) library(ggpubr) # Set appropriate working directory and load datasets: setwd("~/Dropbox/Dropbox Fee/PhD/Projects/Bioturbation experiment/Data management") Experiment <- read.csv("Pilot_experiment.csv") Transects <- read.csv("Transects.csv") ### 1. Analysis pilot experiment # Select endpoint End<-subset(Experiment, Days =="45") # Explore normality assumption with QQ plots ggqqplot(End$Hs_growth) ggqqplot(End$Jellyfish) ggqqplot(End$Thalassia_shoots) # Test normality assumption with Shapiro-Wilk tests shapiro.test(End$Hs_growth) #normally distributed shapiro.test(End$Jellyfish) # not normally distributed shapiro.test(End$Thalassia_shoots) #normally distributed # Test homogeneity of variances using Levene's test levene_test(Hs_growth ~ Treatment, data=End) #>0.05 equal variances levene_test(Jellyfish~ Treatment, data=End) #>0.05 equal variances levene_test(Thalassia_shoots ~ Treatment, data=End) #>0.05 equal variances # Summarize and compute averages, first for Halophila stipulacea new shoots grown per day ggboxplot(End, x = "Treatment", y = "Hs_growth", color="Treatment", palette = c("#00AFBB", "#E7B800","#E7B800"),ylab="Hs shoots)",xlab="") sumHs<-ddply(End, .(Treatment), summarize, Avg = mean(Hs_growth), sd = sd(Hs_growth), N = length(Hs_growth), se = sd/sqrt(N)) # Test for differences between experimental treatments res.aov <- aov(Hs_growth ~ Treatment, data = End) summary(res.aov) # Summarize and compute averages for number of Cassiopea jellyfish at the end of the experiment ggboxplot(End, x = "Treatment", y = "Jellyfish", color="Treatment", palette = c("#00AFBB", "#E7B800","#E7B800"),ylab="Turtle density (ha m-1)",xlab="") sumJelly<-ddply(End, .(Treatment), summarize, Avg = mean(Jellyfish), sd = sd(Jellyfish), N = length(Jellyfish), se = sd/sqrt(N)) # Test for differences between experimental treatments, kruskal.test(Jellyfish ~ Treatment, data = End) ### 2. Analysis Bonaire transect data # Check normality ggqqplot(Transects$Halophila_shoots_m2) ggqqplot(Transects$Mounds) # Test for normality using Shapiro-Wilk tests shapiro.test(Transects$Halophila_shoots_m2) #not normally distributed shapiro.test(Transects$Mounds) # normally distributed # Test for differences between inward and outward habitat using a paired t-test for the number of mounds # and a paired Wilcoxon test for the number of Halophila stipulacea shoots res <- t.test(Mounds ~ Type, data = Transects, paired = TRUE) res res <- wilcox.test(Halophila_shoots_m2 ~ Type, data=Transects,alternative = "two.sided",exact = FALSE, paired=TRUE) res