#Load these libraries (citations at the end)
library(lme4)
library(ggplot2)
library(dplyr)
library(viridis)
data <- read.csv("Puglisi_2021_Data.csv")
# Histogram to analyze normality
hist(data$snailAbund)
#regression model will not be not normal distribution
#most likely poisson
m <- glmer(snailAbund ~ pH + Cd + d.o. + temp + (1|site ), data = data, family = "poisson")
summary(m)
#found strong correlations when plotted with "glmer" but could be due to random site effects
#Basic plots to show relationships
plot(data$snailAbund ~ pH, data=data)
#pH to abundance does not seem to have much of a correlation, indicates a potential wide range
plot(data$snailAbund ~ temp, data=data)
#temp seems to have some correlation, but most likely due to data being gathered in different air temperature and potentially swayed from how abundant sites 1 and 2 were
plot(data$snailAbund ~ Cd, xlab="Conductivity (Âµs/cm)", ylab="Snail Abundance", data=data)
#cD seems to influence snail abundance as well through a negative relationship
plot(data$snailAbund ~ d.o., xlab="Dissolved Oxygen (mg/L)", ylab="Snail Abundance", data=data)
#best correlation so far and seems to show a positive relationship on abundance
m2 <- glm(snailAbund ~ pH + Cd + d.o. + temp, data = data, family = "poisson")
summary(m2)
#shows that the data was not randomly significant without the random effects included
#negative z value indicates a negative relationship and a positive indicates a positive relationship
sitelabs <- c("1","2","3","4","5","6","7","8","9","10")
#changed the x axis tickmarks for plots
p1 <- ggplot(data, aes(site, snailAbund, color=site))+
geom_boxplot()+labs(x="Site", y="Snail Abundance")+
scale_fill_viridis(option="turbo")+theme_classic(base_size=20)+ scale_x_discrete(labels=sitelabs)
p1
#Overall graph showing abundance at each site
p2 <- ggplot(data, aes(d.o., snailAbund, color=site))+
geom_point(shape=17, size=2)+
labs(x="Dissolved Oxygen (mg/L)", y="Snail Abundance")+
scale_fill_viridis(option="turbo")+theme_classic(base_size=20)
p2
#dissolved oxygen relationship plot with color
p3 <- ggplot(data, aes(Cd, snailAbund, color=site))+
geom_point(shape=17, size=2)+
labs(x="Conductivity (µs/cm)", y="Snail Abundance")+
scale_fill_viridis(option="turbo")+theme_classic(base_size=20)
p3
#conductivity relationship shown to be negative with a colored plot
p4 <- ggplot(data, aes(Cd, d.o., color=site))+
geom_point(shape=17, size=2)+
labs(x="Conductivity (µs/cm)", y="Dissolved Oxygen (mg/L)")+
scale_fill_viridis(option="turbo")+theme_classic(base_size=18.5)
p4
#p4 is a predictor vs predictor graph showing the inverse relationship they have and the effect on snail abundance
p5 <- ggplot(data, aes(pH, snailAbund, color=site))+
geom_point(shape=17, size=2)+
labs(x="pH", y="Snail Abundance")+
scale_fill_viridis(option="turbo")+theme_classic(base_size=20)
p5
#p5 shows the relationship between pH and abundance. This does not effect abundance, but shows snails can function in a range of pH levels.
#Citations
citation("dplyr")
citation("lme4")
citation("ggplot2")
citation("viridis")