Analysis of Factorial Experiment Using R
NE Milla, Jr.
2026-02-25
Loading libraries
library(readxl)
library(ExpDes)
library(agricolae)
library(tidyverse)
library(ggstatsplot)
library(car)
library(patchwork)
library(emmeans)
library(multcomp)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 <- 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 ' ' 1Graphical display
#Generate the summary statistics for the levels of Factor B (Age)
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## 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#Assign letters
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()
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 <- 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))
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## 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()