Loading Libraries
library(readxl)
library(ExpDes)
library(agricolae)
##
## Attaching package: 'agricolae'
## The following objects are masked from 'package:ExpDes':
##
## lastC, order.group, tapply.stat
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggstatsplot)
## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(patchwork)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:patchwork':
##
## area
##
## The following object is masked from 'package:dplyr':
##
## select
##
## The following object is masked from 'package:ExpDes':
##
## ginv
##
##
## Attaching package: 'TH.data'
##
## The following object is masked from 'package:MASS':
##
## geyser
library(dplyr)
library(readxl)
#Import the excel file into R
fact_crd <- read_excel("Factorial.xlsx",
sheet = "CRD")
#Display the first 6 rows of the data frame
head(fact_crd)
## # A tibble: 6 × 4
## A B Rep Height
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 46.5
## 2 1 1 2 55.9
## 3 1 1 3 78.7
## 4 1 2 1 49.5
## 5 1 2 2 59.5
## 6 1 2 3 78.7
with(fact_crd,fat2.crd(factor1 = A,
factor2 = B,
resp = Height,
quali = c(TRUE,TRUE),
mcomp = "tukey",
fac.names = c("Spacing", "Age")))
## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: Spacing
## FACTOR 2: Age
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Spacing 1 409.0 3 1.5207 0.241120
## Age 2 12846.3 5 23.8835 0.000066
## Spacing*Age 2 996.6 4 1.8529 0.198942
## Residuals 12 3227.2 2
## Total 17 17479.1 1
## ------------------------------------------------------------------------
## CV = 20.36 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.4529232
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## No significant interaction: analyzing the simple effect
## ------------------------------------------------------------------------
## Spacing
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 85.30000
## 2 2 75.76667
## ------------------------------------------------------------------------
## Age
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 118.0833
## b 2 65.36667
## b 1 58.15
## ------------------------------------------------------------------------
fact_crd <- fact_crd %>%
mutate(Spacing = factor(A, labels = c("10x10", "12x12")),
Age = factor(B, labels = c("6","12","24")))
head(fact_crd)
## # A tibble: 6 × 6
## A B Rep Height Spacing Age
## <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 1 1 1 46.5 10x10 6
## 2 1 1 2 55.9 10x10 6
## 3 1 1 3 78.7 10x10 6
## 4 1 2 1 49.5 10x10 12
## 5 1 2 2 59.5 10x10 12
## 6 1 2 3 78.7 10x10 12
#Create the ANOVA model
fact_mod1 <- aov(Height ~ Spacing + Age + Spacing:Age,
data = fact_crd)
#Display the ANOVA table
summary(fact_mod1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Spacing 1 409 409 1.521 0.241
## Age 2 12846 6423 23.883 6.55e-05 ***
## Spacing:Age 2 997 498 1.853 0.199
## Residuals 12 3227 269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fact_crd <- fact_crd %>%
group_by(Age) %>%
dplyr::summarize(Mean = mean(Height),
SD = sd(Height),
n = length(Height)) %>%
dplyr::mutate(SE = SD/sqrt(n))
fact_crd
## # A tibble: 3 × 5
## Age Mean SD n SE
## <fct> <dbl> <dbl> <int> <dbl>
## 1 6 58.2 12.0 6 4.89
## 2 12 65.4 10.4 6 4.24
## 3 24 118. 26.0 6 10.6
B.means <- emmeans(fact_mod1, specs="Age")
## NOTE: Results may be misleading due to involvement in interactions
B.means
## Age emmean SE df lower.CL upper.CL
## 6 58.1 6.69 12 43.6 72.7
## 12 65.4 6.69 12 50.8 80.0
## 24 118.1 6.69 12 103.5 132.7
##
## Results are averaged over the levels of: Spacing
## Confidence level used: 0.95
B.means.cld <- cld(B.means,
Letters=letters) %>%
arrange(Age) %>%
mutate(SE1 = fact_crd$SE)
B.means.cld
## Age emmean SE df lower.CL upper.CL .group SE1
## 1 6 58.15000 6.694975 12 43.56290 72.73710 a 4.888882
## 2 12 65.36667 6.694975 12 50.77957 79.95376 a 4.237269
## 3 24 118.08333 6.694975 12 103.49624 132.67043 b 10.610008
#Generate the two-way plot
ggplot(B.means.cld, aes(Age, emmean,label = .group)) +
geom_col(fill="lightgreen") +
geom_errorbar(aes(ymin=emmean - SE1, ymax = emmean + SE1),
width = 0.2, size = 0.7) +
geom_text(vjust=-4.5, size = 3) +
geom_text(label=round(B.means.cld$emmean,2),vjust=5, color = "black") +
labs(x= "Age", y= "Mean Height") +
theme_classic()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
2 Factor Factorial RCBD
fact_rbd <- read_excel("Factorial.xlsx",
sheet = "RCBD")
head(fact_rbd)
## # A tibble: 6 × 4
## Block P N Yield
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 0.9
## 2 1 1 2 1.2
## 3 1 1 3 1.3
## 4 1 1 4 1.8
## 5 1 1 5 1.1
## 6 1 2 1 0.9
with(fact_rbd,fat2.rbd(factor1 = P,
factor2 = N,
block = Block,
resp = Yield,
quali = c(TRUE,TRUE),
fac.names = c("P","N"),
mcomp="tukey"))
## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: P
## FACTOR 2: N
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Block 2 0.0640 3 1.631 0.213772
## P 2 0.1693 4 4.316 0.023244
## N 4 2.4902 6 31.732 0.000000
## P*N 8 1.0151 5 6.468 0.000090
## Residuals 28 0.5493 2
## Total 44 4.2880 1
## ------------------------------------------------------------------------
## CV = 11.12 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.1377132
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
##
##
## Significant interaction: analyzing the interaction
## ------------------------------------------------------------------------
##
## Analyzing P inside of each level of N
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## Block 2 0.06400 0.032 1.6311 0.2138
## N 4 2.49022 0.62256 31.7322 0
## P:N 1 2 0.01556 0.00778 0.3964 0.6764
## P:N 2 2 0.12667 0.06333 3.2282 0.0548
## P:N 3 2 0.01556 0.00778 0.3964 0.6764
## P:N 4 2 0.62000 0.31 15.801 0
## P:N 5 2 0.40667 0.20333 10.3641 4e-04
## Residuals 28 0.54933 0.01962
## Total 44 4.28800
## ------------------------------------------------------------------------
##
##
##
## P inside of the level 1 of N
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 0.9333333
## 2 2 0.8333333
## 3 3 0.8666667
## ------------------------------------------------------------------------
##
##
## P inside of the level 2 of N
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 1.2333333
## 2 2 0.9666667
## 3 3 1.2000000
## ------------------------------------------------------------------------
##
##
## P inside of the level 3 of N
##
## According to the F test, the means of this factor are statistical equal.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 1.400000
## 2 2 1.300000
## 3 3 1.366667
## ------------------------------------------------------------------------
##
##
## P inside of the level 4 of N
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 1.933333
## b 3 1.433333
## b 2 1.333333
## ------------------------------------------------------------------------
##
##
## P inside of the level 5 of N
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 1.666667
## b 1 1.233333
## b 3 1.2
## ------------------------------------------------------------------------
##
##
##
## Analyzing N inside of each level of P
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## Block 2 0.06400 0.032 1.6311 0.2138
## P 2 0.16933 0.08467 4.3155 0.0232
## N:P 1 4 1.63067 0.40767 20.7791 0
## N:P 2 4 1.29733 0.32433 16.5316 0
## N:P 3 4 0.57733 0.14433 7.3568 4e-04
## Residuals 28 0.54933 0.01962
## Total 44 4.28800
## ------------------------------------------------------------------------
##
##
##
## N inside of the level 1 of P
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 4 1.933333
## b 3 1.4
## bc 2 1.233333
## bc 5 1.233333
## c 1 0.9333333
## ------------------------------------------------------------------------
##
##
## N inside of the level 2 of P
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 5 1.666667
## b 4 1.333333
## b 3 1.3
## c 2 0.9666667
## c 1 0.8333333
## ------------------------------------------------------------------------
##
##
## N inside of the level 3 of P
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 4 1.433333
## a 3 1.366667
## a 2 1.2
## a 5 1.2
## b 1 0.8666667
## ------------------------------------------------------------------------
fact_rbd <- fact_rbd %>%
mutate(Block = factor(Block),
PGR = factor(P, labels = c("P1", "P2", "P3")),
Nitro = factor(N, labels = c("N1", "N2", "N3", "N4", "N5")))
head(fact_rbd)
## # A tibble: 6 × 6
## Block P N Yield PGR Nitro
## <fct> <dbl> <dbl> <dbl> <fct> <fct>
## 1 1 1 1 0.9 P1 N1
## 2 1 1 2 1.2 P1 N2
## 3 1 1 3 1.3 P1 N3
## 4 1 1 4 1.8 P1 N4
## 5 1 1 5 1.1 P1 N5
## 6 1 2 1 0.9 P2 N1
#Create the ANOVA model
fact_mod2 <- aov(Yield ~ Block + PGR + Nitro + PGR:Nitro,
data = fact_rbd)
#Display the ANOVA table
anova(fact_mod2)
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 2 0.06400 0.03200 1.6311 0.21377
## PGR 2 0.16933 0.08467 4.3155 0.02324 *
## Nitro 4 2.49022 0.62256 31.7322 4.946e-10 ***
## PGR:Nitro 8 1.01511 0.12689 6.4676 8.979e-05 ***
## Residuals 28 0.54933 0.01962
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Generate summary statistics for all treatment combinations
pn.means <- fact_rbd %>%
group_by(P,N) %>%
dplyr::summarise(Mean = mean(Yield),
SD = sd(Yield),
n = length(Yield)) %>%
dplyr::mutate(SE = SD/sqrt(n))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by P and N.
## ℹ Output is grouped by P.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(P, N))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
pn.means
## # A tibble: 15 × 6
## # Groups: P [3]
## P N Mean SD n SE
## <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 1 1 0.933 0.0577 3 0.0333
## 2 1 2 1.23 0.0577 3 0.0333
## 3 1 3 1.4 0.1 3 0.0577
## 4 1 4 1.93 0.153 3 0.0882
## 5 1 5 1.23 0.153 3 0.0882
## 6 2 1 0.833 0.0577 3 0.0333
## 7 2 2 0.967 0.115 3 0.0667
## 8 2 3 1.3 0.2 3 0.115
## 9 2 4 1.33 0.252 3 0.145
## 10 2 5 1.67 0.208 3 0.120
## 11 3 1 0.867 0.153 3 0.0882
## 12 3 2 1.2 0.2 3 0.115
## 13 3 3 1.37 0.0577 3 0.0333
## 14 3 4 1.43 0.0577 3 0.0333
## 15 3 5 1.2 0.1 3 0.0577
PN.means <- emmeans(fact_mod2, specs=c("PGR", "Nitro"))
PN.means
## PGR Nitro emmean SE df lower.CL upper.CL
## P1 N1 0.933 0.0809 28 0.768 1.099
## P2 N1 0.833 0.0809 28 0.668 0.999
## P3 N1 0.867 0.0809 28 0.701 1.032
## P1 N2 1.233 0.0809 28 1.068 1.399
## P2 N2 0.967 0.0809 28 0.801 1.132
## P3 N2 1.200 0.0809 28 1.034 1.366
## P1 N3 1.400 0.0809 28 1.234 1.566
## P2 N3 1.300 0.0809 28 1.134 1.466
## P3 N3 1.367 0.0809 28 1.201 1.532
## P1 N4 1.933 0.0809 28 1.768 2.099
## P2 N4 1.333 0.0809 28 1.168 1.499
## P3 N4 1.433 0.0809 28 1.268 1.599
## P1 N5 1.233 0.0809 28 1.068 1.399
## P2 N5 1.667 0.0809 28 1.501 1.832
## P3 N5 1.200 0.0809 28 1.034 1.366
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95
#AGenerate the letters
PN.means.cld <- cld(PN.means,
Letters=letters,
decreasing = TRUE)
PN.means.cld <- PN.means.cld %>%
dplyr::arrange(PGR,Nitro) %>%
dplyr::mutate(SE1 = pn.means$SE)
PN.means.cld
## PGR Nitro emmean SE df lower.CL upper.CL .group SE1
## 1 P1 N1 0.9333333 0.0808683 28 0.7676821 1.0989845 ef 0.03333333
## 2 P1 N2 1.2333333 0.0808683 28 1.0676821 1.3989845 cdef 0.03333333
## 3 P1 N3 1.4000000 0.0808683 28 1.2343488 1.5656512 bc 0.05773503
## 4 P1 N4 1.9333333 0.0808683 28 1.7676821 2.0989845 a 0.08819171
## 5 P1 N5 1.2333333 0.0808683 28 1.0676821 1.3989845 cdef 0.08819171
## 6 P2 N1 0.8333333 0.0808683 28 0.6676821 0.9989845 f 0.03333333
## 7 P2 N2 0.9666667 0.0808683 28 0.8010155 1.1323179 def 0.06666667
## 8 P2 N3 1.3000000 0.0808683 28 1.1343488 1.4656512 bcde 0.11547005
## 9 P2 N4 1.3333333 0.0808683 28 1.1676821 1.4989845 bcde 0.14529663
## 10 P2 N5 1.6666667 0.0808683 28 1.5010155 1.8323179 ab 0.12018504
## 11 P3 N1 0.8666667 0.0808683 28 0.7010155 1.0323179 f 0.08819171
## 12 P3 N2 1.2000000 0.0808683 28 1.0343488 1.3656512 cdef 0.11547005
## 13 P3 N3 1.3666667 0.0808683 28 1.2010155 1.5323179 bcd 0.03333333
## 14 P3 N4 1.4333333 0.0808683 28 1.2676821 1.5989845 bc 0.03333333
## 15 P3 N5 1.2000000 0.0808683 28 1.0343488 1.3656512 cdef 0.05773503
#Comparing the levels of N for each level of P
ggplot(PN.means.cld, aes(x = PGR,
y = emmean,
fill = Nitro,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-4,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "PGR", y= "Mean Yield") +
scale_fill_brewer(palette = "Blues") +
theme_classic()
#Comparing the levels of P for each level of N
ggplot(PN.means.cld, aes(x = Nitro,
y = emmean,
fill = PGR,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE1,
ymax = emmean + SE1),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-4,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "Nitro", y= "Mean Yield") +
theme_classic()
