Statistical Analysis of Experiments Using R and RStudio
NE Milla, Jr.
2024-10-09
Loading libraries
Analysis of Experiments in RCBD
Importing the data (Make sure to set the appropriate working directory)
agrondata <- read_excel("agrondata.xlsx")
head(agrondata) #displays the 1st 6 rows of the data frame## # A tibble: 6 × 10
## treatment block DTE HT15 HT30 HT45 HT60 DTF DTM HERBAGE
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 4 56.2 102. 152. 210. 44 60 158000
## 2 1 2 4 53.9 92.9 140 194. 44 60 192000
## 3 1 3 4 53.4 90.4 134. 189. 44 60 180000
## 4 1 4 4 52.2 89 129. 183. 44 60 155000
## 5 1 5 4 54.6 102. 150. 207. 44 60 185000
## 6 2 1 7 47.3 73.5 131. 173. 76 90 255000ANOVA for RCBD experiment with post hoc
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 2 2416.67 1208.33 232.588 0.000000
## Block 4 74.27 18.57 3.574 0.059085
## Residuals 8 41.56 5.20
## Total 14 2532.50
## ------------------------------------------------------------------------
## CV = 5.7 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.1688107
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## ------------------------------------------------------------------------
## Homogeneity of variances test
## p-value: 0.07919401
## According to the test of oneillmathews at 5% of significance, the variances can be considered homocedastic.
## ------------------------------------------------------------------------
##
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 54.06
## b 2 42.6
## c 3 23.3
## ------------------------------------------------------------------------Presenting the data (tabular form)
## # A tibble: 3 × 3
## treatment Mean SD
## <dbl> <dbl> <dbl>
## 1 1 54.1 1.48
## 2 2 42.6 4.90
## 3 3 23.3 1.67Presenting the results visually
Compact letter display
#Run ANOVA the traditional way
Trt <- as.factor(agrondata$treatment)
Blk <- as.factor(agrondata$block)
mod1 <- aov(HT15 ~ Trt + Blk,
data = agrondata)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## Trt 2 2416.7 1208.3 232.588 8.17e-08 ***
## Blk 4 74.3 18.6 3.574 0.0591 .
## Residuals 8 41.6 5.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Trt emmean SE df lower.CL upper.CL
## 1 54.1 1.02 8 51.7 56.4
## 2 42.6 1.02 8 40.2 45.0
## 3 23.3 1.02 8 20.9 25.7
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95## Trt emmean SE df lower.CL upper.CL .group
## 3 23.3 1.02 8 20.9 25.7 a
## 2 42.6 1.02 8 40.2 45.0 b
## 1 54.1 1.02 8 51.7 56.4 c
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#Create graphs with SE and letters
ggplot(trt.means.cld, aes(Trt, emmean,label = .group)) +
geom_col(fill = "lightblue", col="black") +
geom_errorbar(aes(ymin=emmean - SE, ymax = emmean + SE),
width = 0.2, size = 0.7) +
geom_text(vjust=-1.5) +
labs(x= "Treatment", y= "Mean Plant Height (15 days)") +
lims(y = c(0,60)) +
theme_classic()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Analysis of Factorial Experiments in CRD
Native tree experiment
- Factor A: spacing (10x10 m2, 12x12 m2)
- Factor B: age at planting (6, 12, 24 mo.)
- Response: height (ft)
Import and explore data
#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.7with(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
## ------------------------------------------------------------------------Graphical display
#Create factors from integer codes
fact_crd$Spacing = factor(fact_crd$A)
fact_crd$Age = factor(fact_crd$B)
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 1 1
## 2 1 1 2 55.9 1 1
## 3 1 1 3 78.7 1 1
## 4 1 2 1 49.5 1 2
## 5 1 2 2 59.5 1 2
## 6 1 2 3 78.7 1 2#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#Generate the summary statistics for the levels of Factor B
B.means <- emmeans(fact_mod1, specs="Age")
B.means#Print the table of means## Age emmean SE df lower.CL upper.CL
## 1 58.1 6.69 12 43.6 72.7
## 2 65.4 6.69 12 50.8 80.0
## 3 118.1 6.69 12 103.5 132.7
##
## Results are averaged over the levels of: Spacing
## Confidence level used: 0.95#Assign the letter designations to the means of the Factor B
B.means.cld <- cld(B.means,
Letters=letters,
decreasing = TRUE)
B.means.cld #Print the table of means## Age emmean SE df lower.CL upper.CL .group
## 3 118.1 6.69 12 103.5 132.7 a
## 2 65.4 6.69 12 50.8 80.0 b
## 1 58.1 6.69 12 43.6 72.7 b
##
## Results are averaged over the levels of: Spacing
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#Generate the two-way plot
ggplot(B.means.cld, aes(Age, emmean,label = .group)) +
geom_col(fill="darkgray") +
geom_errorbar(aes(ymin=emmean - SE, ymax = emmean + SE),
width = 0.2, size = 0.7) +
geom_text(vjust=-3.2, size = 3) +
labs(x= "Age", y= "Mean Height") +
theme_classic()
Analysis of Factorial Experiments in RCBD
A factorial experiment was carried out in order to determine the response of one variety of a crop to three levels of a new plant growth regulator (PGR), denoted by P1, P2, P3, and five levels of nitrogen fertilizer (N1, N2, N3, N4, N5)
A randomized block design with three replications was adopted and the 15 treatments were randomly assigned to the plots within each block
One of the parameters that were recorded is plot yield (kg/ha)
Importing the data into R
## # 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.9Factorial ANOVA (RCBD)
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
## ------------------------------------------------------------------------Graphical display
#Converting the integer codes to factors
fact_rbd$Block = factor(fact_rbd$Block)
fact_rbd$P = factor(fact_rbd$P)
fact_rbd$N = factor(fact_rbd$N)
head(fact_rbd) ## # A tibble: 6 × 4
## Block P N Yield
## <fct> <fct> <fct> <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#Create the ANOVA model
fact_mod2 <- aov(Yield ~ Block + P + N + P:N,
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
## P 2 0.16933 0.08467 4.3155 0.02324 *
## N 4 2.49022 0.62256 31.7322 4.946e-10 ***
## P:N 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 <- emmeans(fact_mod2, specs=c("P", "N"))
PN.means#Print the table of means## P N emmean SE df lower.CL upper.CL
## 1 1 0.933 0.0809 28 0.768 1.099
## 2 1 0.833 0.0809 28 0.668 0.999
## 3 1 0.867 0.0809 28 0.701 1.032
## 1 2 1.233 0.0809 28 1.068 1.399
## 2 2 0.967 0.0809 28 0.801 1.132
## 3 2 1.200 0.0809 28 1.034 1.366
## 1 3 1.400 0.0809 28 1.234 1.566
## 2 3 1.300 0.0809 28 1.134 1.466
## 3 3 1.367 0.0809 28 1.201 1.532
## 1 4 1.933 0.0809 28 1.768 2.099
## 2 4 1.333 0.0809 28 1.168 1.499
## 3 4 1.433 0.0809 28 1.268 1.599
## 1 5 1.233 0.0809 28 1.068 1.399
## 2 5 1.667 0.0809 28 1.501 1.832
## 3 5 1.200 0.0809 28 1.034 1.366
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95#Assign the letters to the treatment means
PN.means.cld <- cld(PN.means,
Letters=letters,
decreasing = TRUE)
PN.means.cld#Print the table of means## P N emmean SE df lower.CL upper.CL .group
## 1 4 1.933 0.0809 28 1.768 2.099 a
## 2 5 1.667 0.0809 28 1.501 1.832 ab
## 3 4 1.433 0.0809 28 1.268 1.599 bc
## 1 3 1.400 0.0809 28 1.234 1.566 bc
## 3 3 1.367 0.0809 28 1.201 1.532 bcd
## 2 4 1.333 0.0809 28 1.168 1.499 bcde
## 2 3 1.300 0.0809 28 1.134 1.466 bcde
## 1 5 1.233 0.0809 28 1.068 1.399 cdef
## 1 2 1.233 0.0809 28 1.068 1.399 cdef
## 3 5 1.200 0.0809 28 1.034 1.366 cdef
## 3 2 1.200 0.0809 28 1.034 1.366 cdef
## 2 2 0.967 0.0809 28 0.801 1.132 def
## 1 1 0.933 0.0809 28 0.768 1.099 ef
## 3 1 0.867 0.0809 28 0.701 1.032 f
## 2 1 0.833 0.0809 28 0.668 0.999 f
##
## Results are averaged over the levels of: Block
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 15 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#Comparing the levels of N for each level of P
ggplot(PN.means.cld, aes(x = P,
y = emmean,
fill = N,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE,
ymax = emmean + SE),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.2,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "P", 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 = N,
y = emmean,
fill = P,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE,
ymax = emmean + SE),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.2,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "N", y= "Mean Yield") +
theme_classic()
Analysis of Experiments with Nested Factors
A study was conducted to assess the diversity of an orchid species in forest reserves located in three municipalities in Leyte. In each municipality, four (4) sites were randomly identified and three (3) 10m x 10m quadrats were laid out in each selected site. Diversity is recorded per quadrat.
Reading the data into R
#Reads the data and converts integer codes into nominal (factor) codes
biod <- read_xlsx("biodiversity.xlsx")
biod$mun <- as.factor(biod$Mun)
biod$site <- as.factor(biod$Site)#Runs nested ANOVA and store the results in m1
m1 <- aov(Diversity~mun/site, data= biod)
#Displays the ANOVA table
anova(m1)## Analysis of Variance Table
##
## Response: Diversity
## Df Sum Sq Mean Sq F value Pr(>F)
## mun 2 4.50 2.25 0.5000 0.61271
## mun:site 9 128.25 14.25 3.1667 0.01156 *
## Residuals 24 108.00 4.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## NOTE: A nesting structure was detected in the fitted model:
## site %in% mun## $emmeans
## mun = 1:
## site emmean SE df lower.CL upper.CL
## 1 11 1.22 24 8.47 13.53
## 2 10 1.22 24 7.47 12.53
## 3 10 1.22 24 7.47 12.53
## 4 12 1.22 24 9.47 14.53
##
## mun = 2:
## site emmean SE df lower.CL upper.CL
## 1 11 1.22 24 8.47 13.53
## 2 11 1.22 24 8.47 13.53
## 3 9 1.22 24 6.47 11.53
## 4 9 1.22 24 6.47 11.53
##
## mun = 3:
## site emmean SE df lower.CL upper.CL
## 1 13 1.22 24 10.47 15.53
## 2 13 1.22 24 10.47 15.53
## 3 7 1.22 24 4.47 9.53
## 4 7 1.22 24 4.47 9.53
##
## Confidence level used: 0.95
##
## $contrasts
## mun = 1:
## contrast estimate SE df t.ratio p.value
## site1 - site2 1 1.73 24 0.577 1.0000
## site1 - site3 1 1.73 24 0.577 1.0000
## site1 - site4 -1 1.73 24 -0.577 1.0000
## site2 - site3 0 1.73 24 0.000 1.0000
## site2 - site4 -2 1.73 24 -1.155 1.0000
## site3 - site4 -2 1.73 24 -1.155 1.0000
##
## mun = 2:
## contrast estimate SE df t.ratio p.value
## site1 - site2 0 1.73 24 0.000 1.0000
## site1 - site3 2 1.73 24 1.155 1.0000
## site1 - site4 2 1.73 24 1.155 1.0000
## site2 - site3 2 1.73 24 1.155 1.0000
## site2 - site4 2 1.73 24 1.155 1.0000
## site3 - site4 0 1.73 24 0.000 1.0000
##
## mun = 3:
## contrast estimate SE df t.ratio p.value
## site1 - site2 0 1.73 24 0.000 1.0000
## site1 - site3 6 1.73 24 3.464 0.0121
## site1 - site4 6 1.73 24 3.464 0.0121
## site2 - site3 6 1.73 24 3.464 0.0121
## site2 - site4 6 1.73 24 3.464 0.0121
## site3 - site4 0 1.73 24 0.000 1.0000
##
## P value adjustment: bonferroni method for 6 tests#Generates the bar plot with error bars for Mun=1
pw1 <- emmeans(m1,
pairwise ~ site,
at = list(mun="1"),
adjust="bonferroni")## NOTE: Results may be misleading due to involvement in interactionspw1.cld <- cld(pw1,
Letters=letters,
decreasing = TRUE)
g1 <- ggplot(pw1.cld, aes(site, emmean,label = "")) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE, ymax = emmean + SE),
width = 0.2, size = 0.7) +
geom_text(vjust=-2.2) +
lims(y = c(0,15)) +
labs(x= "Sites",
y= "Mean Biodiversity Index",
title = "Mun = 1") +
theme_classic()
g1
#Generates the bar plot with error bars for Mun=2
pw2 <- emmeans(m1,
pairwise ~ site,
at = list(mun="2"),
adjust="bonferroni")## NOTE: Results may be misleading due to involvement in interactionspw2.cld <- cld(pw2,
Letters=letters,
decreasing = TRUE)
g2 <- ggplot(pw2.cld, aes(site, emmean,label = "")) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE, ymax = emmean + SE),
width = 0.2, size = 0.7) +
geom_text(vjust=-2.2) +
lims(y = c(0,15)) +
labs(x= "Sites",
y= "Mean Biodiversity Index",
title = "Mun = 2") +
theme_classic() +
theme(
axis.text.y = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()
)
g2
#Generates the bar plot with error bars for Mun=2
pw3 <- emmeans(m1,
pairwise ~ site,
at = list(mun="3"),
adjust="bonferroni")## NOTE: Results may be misleading due to involvement in interactionspw3.cld <- cld(pw3,
Letters=letters,
decreasing = TRUE)
g3 <- ggplot(pw3.cld, aes(site, emmean,label = .group)) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE, ymax = emmean + SE),
width = 0.2, size = 0.7) +
geom_text(vjust=-1.0,
hjust=1.1) +
lims(y = c(0,15)) +
labs(x= "Sites",
y= "Mean Biodiversity Index",
title = "Mun=3") +
theme_classic() +
theme(
axis.text.y = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()
)
g3
Analysis of Split-Plot Experiments in RCBD
An experiment was carried out to compare the effects of three concentrations of a chemical seed dressing on the yield of oats. Three varieties were used in the trial because it is suspected that the response to the seed treatment would depend on the variety used. A split design was laid out in five randomized blocks. Varieties were assigned at random to the main plots within each block, and the seed treatment and control were assigned at random to the subplots within each main plot.
Importing the data
## # A tibble: 6 × 4
## Block Variety Concentration Yield
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 42.9
## 2 1 1 4 44.4
## 3 1 1 3 49.5
## 4 1 1 2 53.8
## 5 1 2 3 59.8
## 6 1 2 1 53.3ANOVA for split plot experiments
with(spplot, split2.rbd(factor1 = Variety,
factor2 = Concentration,
block = Block,
resp = Yield,
quali = c(TRUE, TRUE),
fac.names = c("Variety", "Concentration"),
mcomp = "tukey",
unfold = "1"))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1 (plot): Variety
## FACTOR 2 (split-plot): Concentration
## ------------------------------------------------------------------------
##
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr(>Fc)
## Variety 2 3631.5 1815.77 60.775 1.5e-05 ***
## Block 4 2931.3 732.82 24.528 0.000152 ***
## Error a 8 239.0 29.88
## Concentration 3 350.2 116.73 8.316 0.000248 ***
## Variety*Concentration 6 368.3 61.38 4.372 0.002049 **
## Error b 36 505.4 14.04
## Total 59 8025.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ------------------------------------------------------------------------
## CV 1 = 10.6146 %
## CV 2 = 7.275885 %
##
## No significant interaction: analyzing the simple effects
## ------------------------------------------------------------------------
## Variety
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.255
## b 2 52.88
## c 1 41.35
## ------------------------------------------------------------------------
##
## Concentration
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 55.45333
## b 3 51.22
## b 4 50.29333
## b 1 49.01333
## ------------------------------------------------------------------------
##
##
##
##
## Significant interaction: analyzing the interaction
## ------------------------------------------------------------------------
##
## Analyzing Variety inside of each level of Concentration
## ------------------------------------------------------------------------
## DF SS MS Fc p.value
## Variety : Concentration 1 2.00000 1573.2013 786.60067 43.70568 0.000000
## Variety : Concentration 2 2.00000 495.6253 247.81267 13.76915 0.000048
## Variety : Concentration 3 2.00000 422.6680 211.33400 11.74229 0.000148
## Variety : Concentration 4 2.00000 1508.3253 754.16267 41.90333 0.000000
## Pooled Error 32.22141 579.9106 17.99768 NA NA
## ------------------------------------------------------------------------
##
##
## Variety inside of Concentration 1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.84
## b 2 50.34
## c 1 35.86
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 62.78
## b 2 54.84
## b 1 48.74
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 3
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 56.96
## a 2 52.54
## b 1 44.16
## ------------------------------------------------------------------------
##
## Variety inside of Concentration 4
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 60.44
## b 2 53.8
## c 1 36.64
## ------------------------------------------------------------------------
##
##
## Analyzing Concentration inside of each level of Variety
## ------------------------------------------------------------------------
## DF SS MS Fc p.value
## Concentration : Variety 1 3 574.1620 191.38733 13.633626 0.000004
## Concentration : Variety 2 3 56.2760 18.75867 1.336288 0.277818
## Concentration : Variety 3 3 88.0455 29.34850 2.090663 0.118618
## Error b 36 505.3640 14.03789 NA NA
## ------------------------------------------------------------------------
##
##
## Concentration inside of Variety 1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 48.74
## a 3 44.16
## b 4 36.64
## b 1 35.86
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
##
##
## Concentration inside of Variety 2
## ------------------------------------------------------------------------
## According to F test, the means of this factor are not different.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 50.34
## 2 2 54.84
## 3 3 52.54
## 4 4 53.80
## ------------------------------------------------------------------------
##
## Concentration inside of Variety 3
## ------------------------------------------------------------------------
## According to F test, the means of this factor are not different.
## ------------------------------------------------------------------------
## Levels Means
## 1 1 60.84
## 2 2 62.78
## 3 3 56.96
## 4 4 60.44
## ------------------------------------------------------------------------Another way to generate ANOVA for split plot experiments
#Converts integer codes to nominal codes
spplot <- spplot |>
mutate(Blk = as.factor(Block),
Var = as.factor(Variety),
Conc = as.factor(Concentration))#Generates ANOVA
spplot.out <- aov(Yield ~ Blk + Var + Blk:Var + Conc + Var:Conc,
data = spplot)
anova(spplot.out)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Blk 4 2931.3 732.82 52.2028 1.691e-14 ***
## Var 2 3631.5 1815.77 129.3477 < 2.2e-16 ***
## Conc 3 350.2 116.73 8.3156 0.0002483 ***
## Blk:Var 8 239.0 29.88 2.1283 0.0583461 .
## Var:Conc 6 368.3 61.38 4.3725 0.0020493 **
## Residuals 36 505.4 14.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We need to correct the F ratio and p-value for Blk and Var
#Extract the Mean Square for Blk:Var and use this as MSE(a)
MSE_a <- anova(spplot.out)["Mean Sq"][4,1]
MS_Blk <- anova(spplot.out)["Mean Sq"][1,1]
MS_Var <- anova(spplot.out)["Mean Sq"][2,1]
print(c(MSE_a,MS_Blk, MS_Var))## [1] 29.87704 732.81692 1815.76850The corrected F ratio and p-value for Blk is:
\[ \begin{aligned} F_{Block} &= \frac{MS_{Block}}{MSE(a)} \notag \\ &= 24.528\notag \\ p-value &= 1.5179\times 10^{-4} \end{aligned} \]
The corrected F ratio and p-value for Var is:
\[ \begin{aligned} F_{Var} &= \frac{MS_{Var}}{MSE(a)} \notag \\ &= 60.775\notag \\ p-value &= 5.02\times 10^{-6} \end{aligned} \]
Generate the visualization of the interaction effect
#Generate summary statistics for all treatment combinations
VC.means <- emmeans(spplot.out, specs=c("Var", "Conc"))
VC.means#Print the table of means## Var Conc emmean SE df lower.CL upper.CL
## 1 1 35.9 1.68 36 32.5 39.3
## 2 1 50.3 1.68 36 46.9 53.7
## 3 1 60.8 1.68 36 57.4 64.2
## 1 2 48.7 1.68 36 45.3 52.1
## 2 2 54.8 1.68 36 51.4 58.2
## 3 2 62.8 1.68 36 59.4 66.2
## 1 3 44.2 1.68 36 40.8 47.6
## 2 3 52.5 1.68 36 49.1 55.9
## 3 3 57.0 1.68 36 53.6 60.4
## 1 4 36.6 1.68 36 33.2 40.0
## 2 4 53.8 1.68 36 50.4 57.2
## 3 4 60.4 1.68 36 57.0 63.8
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95#Assign the letters to the treatment means
VC.means.cld <- cld(VC.means,
Letters=letters,
decreasing = TRUE)
VC.means.cld#Print the table of means## Var Conc emmean SE df lower.CL upper.CL .group
## 3 2 62.8 1.68 36 59.4 66.2 a
## 3 1 60.8 1.68 36 57.4 64.2 ab
## 3 4 60.4 1.68 36 57.0 63.8 abc
## 3 3 57.0 1.68 36 53.6 60.4 abcd
## 2 2 54.8 1.68 36 51.4 58.2 abcd
## 2 4 53.8 1.68 36 50.4 57.2 bcd
## 2 3 52.5 1.68 36 49.1 55.9 cd
## 2 1 50.3 1.68 36 46.9 53.7 de
## 1 2 48.7 1.68 36 45.3 52.1 de
## 1 3 44.2 1.68 36 40.8 47.6 ef
## 1 4 36.6 1.68 36 33.2 40.0 fg
## 1 1 35.9 1.68 36 32.5 39.3 g
##
## Results are averaged over the levels of: Blk
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 12 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.#Comparing the levels of Concentration for each level of Variety
ggplot(VC.means.cld, aes(x = Var,
y = emmean,
fill = Conc,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE,
ymax = emmean + SE),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.2,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "Variety", y= "Mean Yield") +
scale_fill_brewer(palette = "Greens") +
theme_classic()
#Comparing the levels of Variety for each level of Concentration
ggplot(VC.means.cld, aes(x = Conc,
y = emmean,
fill = Var,
label = .group)) +
geom_col(position = "dodge") +
geom_errorbar(aes(ymin = emmean - SE,
ymax = emmean + SE),
width = 0.2,
size = 0.7,
position = position_dodge(0.9)) +
geom_text(vjust=-3.2,
position = position_dodge(0.9),
size = 2.5) +
labs(x= "Concentration", y= "Mean Yield") +
scale_fill_brewer(palette = "Greens") +
theme_classic()