Variation of mussel
population
Mean population
density
# Calculate mean
a=4
mean= (population$i1+ population$i2+ population$i3)/(3*a)
# tiff("docs/Fig_3.tiff", height=3, width=5, units='in', res=600)
par(mar=c(5, 4.5, 1, 1))
plot(mean,
type='b', col='blue',
xlab='', ylab=expression(Population~density~(no./m^2)), xlim=c(1,12), xaxt='n', lwd=2)
# legend('top', legend=c('Mean'), col=c('blue'), lwd=c(3), lty=1)
axis(1, at=seq(1,12), population$Month, las=2, cex.axis = 1)
# title('Monthly variation in population density of L. marginalis')
# dev.off()
Individual quadrat
population
plot(population$i1,
type='b', col='blue',
xlab='', xlim=c(1,12), xaxt='n')
axis(1, at=seq(1,12), population$Month, las=2)
plot(population$i2,
type='b', col='blue',
xlab='', xlim=c(1,12), xaxt='n')
axis(1, at=seq(1,12), population$Month, las=2)
plot(population$i3,
type='b', col='blue',
xlab='', xlim=c(1,12), xaxt='n')
axis(1, at=seq(1,12), population$Month, las=2)
Analyze quadrat
samples using Mann-Whitney U test
wilcox.test(population$i1, population$i2)
## Warning in wilcox.test.default(population$i1, population$i2): cannot compute
## exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: population$i1 and population$i2
## W = 100.5, p-value = 0.1041
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(population$i2, population$i3)
## Warning in wilcox.test.default(population$i2, population$i3): cannot compute
## exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: population$i2 and population$i3
## W = 48.5, p-value = 0.1783
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(population$i3, population$i1)
## Warning in wilcox.test.default(population$i3, population$i1): cannot compute
## exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: population$i3 and population$i1
## W = 59.5, p-value = 0.4858
## alternative hypothesis: true location shift is not equal to 0
Variation of water
quality parameters
water=read.csv("data/water-quality.csv")
for (w in c('Depth','Temperature','pH')) {
arr=water[[w]]
cat(sprintf('%15s - mean: %f, sd: %f min: %f, max: %f\n', w, mean(arr), sd(arr),
min(arr), max(arr)))
}
## Depth - mean: 1.506667, sd: 0.295030 min: 1.150000, max: 2.100000
## Temperature - mean: 23.589167, sd: 3.712300 min: 15.670000, max: 27.350000
## pH - mean: 7.367500, sd: 0.452611 min: 6.410000, max: 7.720000
Temperature
library(ggplot2)
# Month needs to be factorized to facilitate plotting and preserve order
Month1=factor(water$Month, levels=water$Month)
df= data.frame(Month=Month1, Depth=water$Depth,
Temperature=water$Temperature, pH=water$pH, mean=mean)
scale=6
ggplot(data = df, aes(x = Month, group = 1)) +
geom_line(aes(y = mean, color = "black"), linewidth = 1.3) +
geom_point(aes(y = mean, color = "black"), shape=1, size=3) +
geom_line(aes(y = Temperature / scale, color = "blue"), linewidth = 1.3) +
geom_point(aes(y = Temperature / scale, color = "blue"), shape=1, size=3) +
scale_y_continuous(name = "Population / m^2",
sec.axis = sec_axis(transform = ~ . * scale, name = "Temperature °C")) +
scale_color_manual('', labels= c('Population / m^2','Temperature °C'),
values = c("black", "blue")) +
theme_bw() +
theme(legend.position = "top", axis.text.x= element_text(size=11, angle = 90, vjust = 0.5),
axis.text.y.right= element_text(size=11, color='blue'),
axis.title.y.right = element_text(color='blue'))
Depth
scale=1
ggplot(data = df, aes(x = Month, group = 1)) +
geom_line(aes(y = mean, color = "black"), linewidth = 1.3) +
geom_point(aes(y = mean, color = "black"), shape=1, size=3) +
geom_line(aes(y = Depth / scale, color = "blue"), linewidth = 1.3) +
geom_point(aes(y = Depth / scale, color = "blue"), shape=1, size=3) +
scale_y_continuous(name = "Population / m^2",
sec.axis = sec_axis(transform = ~ . * scale, name = "Depth")) +
scale_color_manual('', labels= c('Population / m^2','Depth'),
values = c("black", "blue")) +
theme_bw() +
theme(legend.position = "top", axis.text.x= element_text(size=11, angle = 90, vjust = 0.5),
axis.text.y.right= element_text(size=11, color='blue'),
axis.title.y.right = element_text(color='blue'))
pH
scale=2
ggplot(data = df, aes(x = Month, group = 1)) +
geom_line(aes(y = mean, color = "black"), linewidth = 1.3) +
geom_point(aes(y = mean, color = "black"), shape=1, size=3) +
geom_line(aes(y = pH / scale, color = "blue"), linewidth = 1.3) +
geom_point(aes(y = pH / scale, color = "blue"), shape=1, size=3) +
scale_y_continuous(name = "Population / m^2",
sec.axis = sec_axis(transform = ~ . * scale, name = "pH")) +
scale_color_manual('', labels= c('Population / m^2','pH'),
values = c("black", "blue")) +
theme_bw() +
theme(legend.position = "top", axis.text.x= element_text(size=11, angle = 90, vjust = 0.5),
axis.text.y.right= element_text(size=11, color='blue'),
axis.title.y.right = element_text(color='blue'))
Correlation between
mussel population density and water quality parameters
Plot data points
water=read.csv("data/water-quality.csv")
str(water)
## 'data.frame': 12 obs. of 4 variables:
## $ Month : chr "Apr 2023" "May 2023" "Jun 2023" "Jul 2023" ...
## $ Depth : num 1.15 1.35 1.56 1.77 1.93 2.1 1.55 1.46 1.42 1.35 ...
## $ Temperature: num 25.8 26.1 26.3 26.2 27.1 ...
## $ pH : num 7.63 7.65 7.67 7.72 7.64 7.7 7.44 6.81 6.76 6.41 ...
plot(water$Depth, mean,
type='p', col='blue',
xlab='Depth (m)', ylab='Population / m^2')
plot(water$Temperature, mean,
type='p', col='blue',
xlab='Temperature °C', ylab='Population / m^2')
plot(water$pH, mean,
type='p', col='blue',
xlab='pH', ylab='Population / m^2')
Reference: https://www.sthda.com/english/wiki/correlation-test-between-two-variables-in-r#compute-correlation-in-r
Compute
correlation
for (i in colnames(df)[2:4]) {
print(i)
print(cor.test(df[[i]], mean, method='spearman'))
cat(sprintf('\n'))
}
## [1] "Depth"
## Warning in cor.test.default(df[[i]], mean, method = "spearman"): Cannot compute
## exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: df[[i]] and mean
## S = 559.94, p-value = 9.792e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.9578222
##
##
## [1] "Temperature"
## Warning in cor.test.default(df[[i]], mean, method = "spearman"): Cannot compute
## exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: df[[i]] and mean
## S = 476.67, p-value = 0.0179
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.6666708
##
##
## [1] "pH"
## Warning in cor.test.default(df[[i]], mean, method = "spearman"): Cannot compute
## exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: df[[i]] and mean
## S = 414.45, p-value = 0.143
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4491256
Visualize
correlation
library("ggpubr")
## Warning: package 'ggpubr' was built under R version 4.4.3
library(cowplot)
## Warning: package 'cowplot' was built under R version 4.4.3
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
##
## get_legend
# tiff("docs/Fig_4.tiff", height=9, width=5, units='in', res=600)
figA= ggscatter(df, x = "Depth", y = "mean",
add = "reg.line", conf.int = TRUE,
xlab = "Depth (m)",
title = 'A') +
labs(y=expression(Population~density~(no./m^2))) +
stat_cor(cor.coef.name = 'rho', label.x=1.5, label.y=3.5, method='spearman',
size=5, p.accuracy=0.001)
figB= ggscatter(df, x = "Temperature", y = "mean",
add = "reg.line", conf.int = TRUE,
xlab = "Temperature (°C)",
title = 'B') +
labs(y=expression(Population~density~(no./m^2))) +
stat_cor(cor.coef.name = 'rho', label.x=20, label.y=4.3, method='spearman',
size=5, p.accuracy=0.001)
figC= ggscatter(df, x = "pH", y = "mean",
add = "reg.line", conf.int = TRUE,
xlab = "pH",
title = 'C') +
labs(y=expression(Population~density~(no./m^2))) +
stat_cor(cor.coef.name = 'rho', label.x=6.9, label.y=4.3, method='spearman',
size=5, p.accuracy=0.001)
plot_grid(figA, figB, figC, ncol=1, align = "v")
# dev.off()