Giant African Snails (Achatina fulica) as Feeds for Mud Crabs (Scylla serrata) Statistical Test
Kyle Kenneth Ruaya
2024-05-02
library(readxl)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.4.4 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ car::recode() masks dplyr::recode()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggpubr)
library(rstatix)##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filterData
library(readxl)
Data <- read_excel("D:/SNSHS STATS/Hannah/HANNAH.xlsx")
View(Data)library(rmarkdown)
paged_table(Data)Normality Test
res_aov <- aov(Weight ~ Treatment,
data = Data
)shapiro.test(res_aov$residuals)##
## Shapiro-Wilk normality test
##
## data: res_aov$residuals
## W = 0.97177, p-value = 0.6491res_aov <- aov(Length ~ Treatment,
data = Data
)shapiro.test(res_aov$residuals)##
## Shapiro-Wilk normality test
##
## data: res_aov$residuals
## W = 0.86298, p-value = 0.002104res_aov <- aov(Width ~ Treatment,
data = Data
)shapiro.test(res_aov$residuals)##
## Shapiro-Wilk normality test
##
## data: res_aov$residuals
## W = 0.94868, p-value = 0.1991Homogeneity of Variance
leveneTest(`Weight` ~ Treatment,
data = Data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2382 0.3078
## 24The p-value is greater than the 0.05 level of significance. Thus, the homogeneity assumption of the variance is met.
leveneTest(`Length` ~ Treatment,
data = Data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.8006 0.1868
## 24The p-value is greater than the 0.05 level of significance. Thus, the homogeneity assumption of the variance is met.
leveneTest(`Width` ~ Treatment,
data = Data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 6.3509 0.006113 **
## 24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Weight
group_by(Data, Treatment) %>%
summarise(
count = n(),
mean = mean(Weight, na.rm = TRUE),
sd = sd(Weight, na.rm = TRUE)
)## # A tibble: 3 × 4
## Treatment count mean sd
## <chr> <int> <dbl> <dbl>
## 1 Treatment 1 9 90.7 14.0
## 2 Treatment 2 9 72.7 11.1
## 3 Treatment 3 9 68.3 6.85library("ggpubr")
ggboxplot(Data, x = "Treatment", y = "Weight",
color = "Treatment", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Weight", xlab = "Treatment")
ggline(Data, x = "Treatment", y = "Weight",
add = c("mean_se", "jitter"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Weight", xlab = "Treatment")
library("gplots")##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplotmeans(Weight ~ Treatment, data = Data, frame = TRUE,
xlab = "Treatment", ylab = "Weight",
main="Mean Plot with 95% CI")## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter## Warning in axis(1, at = 1:length(means), labels = legends, ...): "frame" is not
## a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter
One-way ANOVA
res.aov <- aov(Weight ~ Treatment, data = Data)
summary(res.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 2 2532 1265.9 10.39 0.000562 ***
## Residuals 24 2924 121.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since p-value = 0.000562 < 0.05, it is conclusive that we reject the null hypothesis, that is, there is significant difference between the treatments in terms of weight.
TukeyHSD(res.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Weight ~ Treatment, data = Data)
##
## $Treatment
## diff lwr upr p adj
## Treatment 2-Treatment 1 -17.955556 -30.95036 -4.960752 0.0056880
## Treatment 3-Treatment 1 -22.400000 -35.39480 -9.405197 0.0006878
## Treatment 3-Treatment 2 -4.444444 -17.43925 8.550359 0.6736899Length
ggboxplot(Data, x = "Treatment", y = "Length",
color = "Treatment", palette = c("#FC4E07", "#00FF00", "#4E07FC"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Length", xlab = "Treatment")
ggline(Data, x = "Treatment", y = "Length",
add = c("mean_se", "jitter"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Length", xlab = "Treatment")
library("gplots")
plotmeans(`Length` ~ Treatment, data = Data, frame = TRUE,
xlab = "Treatment", ylab = "Length",
main="Mean Plot with 95% CI")## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter## Warning in axis(1, at = 1:length(means), labels = legends, ...): "frame" is not
## a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter
PlantGrowth <- Data %>%
reorder_levels(Treatment, order = c("Treatment 1", "Treatment 2", "Treatment 3"))PlantGrowth %>%
group_by(Treatment) %>%
get_summary_stats(`Length`, type = "common")## # A tibble: 3 × 11
## Treatment variable n min max median iqr mean sd se ci
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Treatment 1 Length 9 5.2 5.74 5.46 0.1 5.43 0.164 0.055 0.126
## 2 Treatment 2 Length 9 2.82 6.02 4.74 0.72 4.84 0.894 0.298 0.687
## 3 Treatment 3 Length 9 4.74 6.12 4.98 0.88 5.24 0.532 0.177 0.409ggboxplot(PlantGrowth, x = "Treatment", y = "Length", fill="Treatment")
Kruskal-wallis Test
res.kruskal <- PlantGrowth %>% kruskal_test(`Length` ~ Treatment)
res.kruskal## # A tibble: 1 × 6
## .y. n statistic df p method
## * <chr> <int> <dbl> <int> <dbl> <chr>
## 1 Length 27 3.78 2 0.151 Kruskal-WallisBased on the p-value, there is no significant difference was observed between the group pairs in terms of length.
Pairwise Comparisons
res1<- PlantGrowth %>%
dunn_test(`Length` ~ Treatment, p.adjust.method = "bonferroni")
res1## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Length Treatment 1 Treatment 2 9 9 -1.94 0.0520 0.156 ns
## 2 Length Treatment 1 Treatment 3 9 9 -0.927 0.354 1 ns
## 3 Length Treatment 2 Treatment 3 9 9 1.02 0.309 0.928 nsWidth
ggboxplot(Data, x = "Treatment", y = "Width",
color = "Treatment", palette = c("#FC4E07", "#00FF00", "#4E07FC"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Width", xlab = "Treatment")
ggline(Data, x = "Treatment", y = "Width",
add = c("mean_se", "jitter"),
order = c("Treatment 1", "Treatment 2", "Treatment 3"),
ylab = "Width", xlab = "Treatment")
library("gplots")
plotmeans(`Width` ~ Treatment, data = Data, frame = TRUE,
xlab = "Treatment", ylab = "Width",
main="Mean Plot with 95% CI")## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter## Warning in axis(1, at = 1:length(means), labels = legends, ...): "frame" is not
## a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "frame" is not a
## graphical parameter
PlantGrowth <- Data %>%
reorder_levels(Treatment, order = c("Treatment 1", "Treatment 2", "Treatment 3"))PlantGrowth %>%
group_by(Treatment) %>%
get_summary_stats(`Width`, type = "common")## # A tibble: 3 × 11
## Treatment variable n min max median iqr mean sd se ci
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Treatment 1 Width 9 7.58 8.52 8.06 0.22 8.05 0.296 0.099 0.228
## 2 Treatment 2 Width 9 7.02 8.6 7.02 0.96 7.42 0.617 0.206 0.474
## 3 Treatment 3 Width 9 2.76 9.1 7.36 3.42 6.17 2.19 0.73 1.68ggboxplot(PlantGrowth, x = "Treatment", y = "Width", fill="Treatment")
Welch Test
oneway.test(Width ~ Treatment, var.equal = FALSE, data = Data)##
## One-way analysis of means (not assuming equal variances)
##
## data: Width and Treatment
## F = 6.3665, num df = 2.000, denom df = 12.735, p-value = 0.01212Based on the p-value, there is significant difference was observed between the group pairs in terms.
Pairwise Comparisons
res1<- PlantGrowth %>%
dunn_test(`Width` ~ Treatment, p.adjust.method = "bonferroni")
res1## # A tibble: 3 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 Width Treatment 1 Treatment 2 9 9 -2.06 0.0397 0.119 ns
## 2 Width Treatment 1 Treatment 3 9 9 -2.46 0.0139 0.0417 *
## 3 Width Treatment 2 Treatment 3 9 9 -0.402 0.687 1 ns