library (knitr)
## Warning: package 'knitr' was built under R version 4.0.5
###Getting the directory###
getwd()
## [1] "C:/Users/Acer/Desktop"
####setting up working directory####
setwd("C:/Users/Acer/Desktop")
####Creating data frame####
###We name the dataframe RPractice###
Practice<-read.csv("C:/Users/Acer/Desktop/RPractrice.csv")
###Verifying data frame###
names(Practice)
## [1] "Treatments" "Block" "Samples" "Response1" "Response2"
## [6] "Response3" "Response4" "Response5" "Response6"
head (Practice)
## Treatments Block Samples Response1 Response2 Response3 Response4 Response5
## 1 T1 R1 S1 17.62 23.25 22.47 25.51 17.68
## 2 T1 R1 S2 13.05 16.43 12.47 16.70 11.28
## 3 T1 R1 S3 13.91 13.68 16.56 24.76 15.59
## 4 T1 R1 S4 14.31 18.05 29.68 32.43 32.53
## 5 T2 R1 S1 18.11 26.46 23.26 29.65 42.90
## 6 T2 R1 S2 14.45 22.31 23.60 31.44 43.81
## Response6
## 1 17.40
## 2 11.95
## 3 17.98
## 4 21.37
## 5 46.78
## 6 49.16
str(Practice)
## 'data.frame': 60 obs. of 9 variables:
## $ Treatments: chr "T1" "T1" "T1" "T1" ...
## $ Block : chr "R1" "R1" "R1" "R1" ...
## $ Samples : chr "S1" "S2" "S3" "S4" ...
## $ Response1 : num 17.6 13.1 13.9 14.3 18.1 ...
## $ Response2 : num 23.2 16.4 13.7 18.1 26.5 ...
## $ Response3 : num 22.5 12.5 16.6 29.7 23.3 ...
## $ Response4 : num 25.5 16.7 24.8 32.4 29.6 ...
## $ Response5 : num 17.7 11.3 15.6 32.5 42.9 ...
## $ Response6 : num 17.4 11.9 18 21.4 46.8 ...
summary(Practice)###quick stats of dataframe##
## Treatments Block Samples Response1
## Length:60 Length:60 Length:60 Min. : 9.14
## Class :character Class :character Class :character 1st Qu.:14.94
## Mode :character Mode :character Mode :character Median :16.40
## Mean :17.40
## 3rd Qu.:19.03
## Max. :31.88
## Response2 Response3 Response4 Response5
## Min. : 7.05 Min. : 1.12 Min. : 0.00 Min. : 0.00
## 1st Qu.:15.22 1st Qu.:15.94 1st Qu.:18.75 1st Qu.:15.49
## Median :18.18 Median :19.64 Median :24.28 Median :20.50
## Mean :18.57 Mean :19.78 Mean :24.74 Mean :28.47
## 3rd Qu.:22.07 3rd Qu.:23.60 3rd Qu.:31.52 3rd Qu.:43.92
## Max. :31.60 Max. :37.65 Max. :42.10 Max. :59.04
## Response6
## Min. : 0.00
## 1st Qu.:17.61
## Median :21.54
## Mean :32.36
## 3rd Qu.:54.46
## Max. :67.08
#####Installing prerequesites####
library(agricolae)
## Warning: package 'agricolae' was built under R version 4.0.5
###Perform ANOVA###
PracAnov<-aov(Response1~Treatments+Block,data = Practice)
summary(PracAnov) ###ouput of ANOVA##
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatments 4 284.8 71.21 4.614 0.00286 **
## Block 2 97.8 48.90 3.168 0.05014 .
## Residuals 53 817.9 15.43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Check histogram of data###
hist(Practice$Response1)
###Assumption Checking of Data###
plot(PracAnov)
## hat values (leverages) are all = 0.1166667
## and there are no factor predictors; no plot no. 5
##Check normal distribution of data##
shapiro.test(PracAnov$residuals)
##
## Shapiro-Wilk normality test
##
## data: PracAnov$residuals
## W = 0.97975, p-value = 0.4181
###In case needed, Bartlett may also be perform##
bartlett.test(Response1~Treatments, data=Practice)
##
## Bartlett test of homogeneity of variances
##
## data: Response1 by Treatments
## Bartlett's K-squared = 13.165, df = 4, p-value = 0.0105
##Post hoc Analysis HSD###
TukeyHSD(PracAnov, which ="Treatments") ###Lwr and upr limit##
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Response1 ~ Treatments + Block, data = Practice)
##
## $Treatments
## diff lwr upr p adj
## T2-T1 0.76833333 -3.7605573 5.297224 0.9889480
## T3-T1 5.09583333 0.5669427 9.624724 0.0200700
## T4-T1 -1.41666667 -5.9455573 3.112224 0.9017957
## T5-T1 0.73166667 -3.7972240 5.260557 0.9908155
## T3-T2 4.32750000 -0.2013906 8.856391 0.0676039
## T4-T2 -2.18500000 -6.7138906 2.343891 0.6539110
## T5-T2 -0.03666667 -4.5655573 4.492224 0.9999999
## T4-T3 -6.51250000 -11.0413906 -1.983609 0.0014690
## T5-T3 -4.36416667 -8.8930573 0.164724 0.0640516
## T5-T4 2.14833333 -2.3805573 6.677224 0.6681014
Tukey1<-HSD.test(PracAnov, trt = "Treatments")##mean separation##
Tukey1 ##output of Tukey##
## $statistics
## MSerror Df Mean CV MSD
## 15.4325 53 17.40333 22.57281 4.528891
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Treatments 5 3.993598 0.05
##
## $means
## Response1 std r Min Max Q25 Q50 Q75
## T1 16.36750 2.430979 12 13.05 21.11 14.9175 15.735 17.9625
## T2 17.13583 3.838604 12 9.14 23.03 14.3850 17.985 19.6225
## T3 21.46333 6.104727 12 13.82 31.88 17.0100 19.715 27.4825
## T4 14.95083 2.419634 12 10.67 18.80 13.6300 15.175 16.3000
## T5 17.09917 4.413612 12 9.31 28.03 15.4450 16.615 18.7550
##
## $comparison
## NULL
##
## $groups
## Response1 groups
## T3 21.46333 a
## T2 17.13583 ab
## T5 17.09917 ab
## T1 16.36750 b
## T4 14.95083 b
##
## attr(,"class")
## [1] "group"
plot(Tukey1) ##Plot of Tukey##
###Post Hoc LSD###
###p.adjustment.methods##
##c("none", "bonferroni", "sidak", "holm", "hs", "hochberg", "bh", "by")
###Choose on the p adjustment methods###
LSD<- LSD.test(PracAnov,"Treatments",alpha = 0.05, p.adj="fdr")
LSD ###output of LSD###
## $statistics
## MSerror Df Mean CV t.value MSD
## 15.4325 53 17.40333 22.57281 2.005746 3.216759
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD fdr Treatments 5 0.05
##
## $means
## Response1 std r LCL UCL Min Max Q25 Q50 Q75
## T1 16.36750 2.430979 12 14.09291 18.64209 13.05 21.11 14.9175 15.735 17.9625
## T2 17.13583 3.838604 12 14.86124 19.41043 9.14 23.03 14.3850 17.985 19.6225
## T3 21.46333 6.104727 12 19.18874 23.73793 13.82 31.88 17.0100 19.715 27.4825
## T4 14.95083 2.419634 12 12.67624 17.22543 10.67 18.80 13.6300 15.175 16.3000
## T5 17.09917 4.413612 12 14.82457 19.37376 9.31 28.03 15.4450 16.615 18.7550
##
## $comparison
## NULL
##
## $groups
## Response1 groups
## T3 21.46333 a
## T2 17.13583 b
## T5 17.09917 b
## T1 16.36750 b
## T4 14.95083 b
##
## attr(,"class")
## [1] "group"
plot(LSD) ##Plotting the LSD##
#### Data Visualization Simple Boxplot###
boxplot(Response1~Treatments, data = Practice)
#### Data Visualization Simple Colored Boxplot###
library(ggplot2)
ggplot(Practice, aes(x=Treatments, y=Response1)) +
geom_boxplot(color="red", fill="blue", lwd=1, alpha=0.5)
#### Data Visualization Boxplot with jitters###
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.0.5
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1
## v purrr 0.3.4
## Warning: package 'tibble' was built under R version 4.0.5
## Warning: package 'tidyr' was built under R version 4.0.5
## Warning: package 'readr' was built under R version 4.0.5
## Warning: package 'dplyr' was built under R version 4.0.5
## Warning: package 'forcats' was built under R version 4.0.5
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.0.5
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.0.5
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(broom)
## Warning: package 'broom' was built under R version 4.0.5
ggboxplot(Practice, x = "Treatments", y = "Response1", add = "jitter")
####Data visualization with jitters data effects###
ggplot(Practice) +
aes(
x = Treatments,
y = Response1,
fill = Treatments,
colour = Response1,
size = Response1
) +
geom_boxplot() +
geom_jitter() +
scale_fill_hue(direction = 1) +
scale_color_gradient() +
theme_minimal()
####Data Visualization Rain Cloud Plot###
library(ggforce)
## Warning: package 'ggforce' was built under R version 4.0.5
library(ggdist)
## Warning: package 'ggdist' was built under R version 4.0.5
library(gghalves)
## Warning: package 'gghalves' was built under R version 4.0.5
ggplot(Practice, aes(Treatments, Response1)) +
ggdist::stat_halfeye(adjust = .3, width = .7, .width = 0, justification = -.2, fill="Dark Green", point_colour = "blue") +
geom_boxplot(width = .2, outlier.shape = 18) +
geom_jitter(width = .05, alpha = .3)
### In case needed, Non Parametric Analysis##
kruskal.test(Response1 ~ Treatments, data = Practice)
##
## Kruskal-Wallis rank sum test
##
## data: Response1 by Treatments
## Kruskal-Wallis chi-squared = 10.351, df = 4, p-value = 0.03491
###Dunns Test###
library(dunn.test)
dunn.test (Practice$Response1, g=Practice$Treatments, method="sidak", kw=TRUE, label=TRUE,
wrap=FALSE, table=TRUE, list=FALSE, rmc=FALSE, alpha=0.05, altp=FALSE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 10.3513, df = 4, p-value = 0.03
##
##
## Comparison of x by group
## (Šidák)
## Col Mean-|
## Row Mean | T1 T2 T3 T4
## ---------+--------------------------------------------
## T2 | -0.672085
## | 0.9443
## |
## T3 | -2.144828 -1.472743
## | 0.1488 0.5181
## |
## T4 | 0.981828 1.653914 3.126657
## | 0.8314 0.3954 0.0088*
## |
## T5 | -0.502602 0.169482 1.642225 -1.484431
## | 0.9747 0.9965 0.4030 0.5100
##
## alpha = 0.05
## Reject Ho if p <= alpha/2