Statistical Analysis of Experiments with Nested Factors using R
NE Milla, Jr.
2026-02-26
Loading libraries
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-by-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 <- biod %>%
mutate(Mun = factor(Mun),
Site = factor(Site))
biod## # A tibble: 36 × 4
## Mun Site Quadrat Diversity
## <fct> <fct> <dbl> <dbl>
## 1 1 1 1 10
## 2 1 1 2 14
## 3 1 1 3 9
## 4 1 2 1 12
## 5 1 2 2 8
## 6 1 2 3 10
## 7 1 3 1 8
## 8 1 3 2 10
## 9 1 3 3 12
## 10 1 4 1 13
## # ℹ 26 more rows#Runs nested ANOVA and store the results in m1
mod1 <- aov(Diversity~Mun/Site, data= biod)
#Displays the ANOVA table
anova(mod1)## 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 ' ' 1TRT.means <- biod %>%
group_by(Mun, Site) %>%
dplyr::summarize(Mean=mean(Diversity),
SD = sd(Diversity),
n = length(Diversity)) %>%
dplyr::mutate(SE = SD/sqrt(n))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Mun and Site.
## ℹ Output is grouped by Mun.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Mun, Site))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.## # A tibble: 12 × 6
## # Groups: Mun [3]
## Mun Site Mean SD n SE
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 1 1 11 2.65 3 1.53
## 2 1 2 10 2 3 1.15
## 3 1 3 10 2 3 1.15
## 4 1 4 12 1 3 0.577
## 5 2 1 11 3 3 1.73
## 6 2 2 11 2 3 1.15
## 7 2 3 9 1 3 0.577
## 8 2 4 9 1 3 0.577
## 9 3 1 13 3 3 1.73
## 10 3 2 13 1 3 0.577
## 11 3 3 7 2 3 1.15
## 12 3 4 7 3 3 1.73#Calculates summary statistics for Sites in Mun 1
Mun1.means <- biod %>%
filter(Mun==1) %>%
group_by(Site) %>%
dplyr::summarize(Mean=mean(Diversity),
SD = sd(Diversity),
n = length(Diversity)) %>%
dplyr::mutate(SE = SD/sqrt(n))
pw1 <- emmeans(mod1,
pairwise ~ Site,
at = list(Mun="1"),
adjust="bonferroni")
#Assigns letters
pw1.cld <- cld(pw1,
Letters=letters,
decreasing = TRUE)
pw1.cld$SE1=Mun1.means$SE
pw1.cld## Site emmean SE df lower.CL upper.CL .group SE1
## 4 12 1.22 24 9.47 14.5 a 1.528
## 1 11 1.22 24 8.47 13.5 a 1.155
## 3 10 1.22 24 7.47 12.5 a 1.155
## 2 10 1.22 24 7.47 12.5 a 0.577
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 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.#Generates the bar plot with error bars for Mun=1
g1 <- ggplot(pw1.cld, aes(Site, emmean,label = "")) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE1, ymax = emmean + SE1),
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
#Calculates summary statistics for Sites in Mun 2
Mun2.means <- biod %>%
filter(Mun==2) %>%
group_by(Site) %>%
dplyr::summarize(Mean=mean(Diversity),
SD = sd(Diversity),
n = length(Diversity)) %>%
dplyr::mutate(SE = SD/sqrt(n))
pw2 <- emmeans(mod1,
pairwise ~ Site,
at = list(Mun="2"),
adjust="bonferroni")
#Assigns letters
pw2.cld <- cld(pw2,
Letters=letters,
decreasing = TRUE)
pw2.cld$SE1=Mun2.means$SE
pw2.cld## Site emmean SE df lower.CL upper.CL .group SE1
## 2 11 1.22 24 8.47 13.5 a 1.732
## 1 11 1.22 24 8.47 13.5 a 1.155
## 4 9 1.22 24 6.47 11.5 a 0.577
## 3 9 1.22 24 6.47 11.5 a 0.577
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 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.#Generates the bar plot with error bars for Mun=2
g2 <- ggplot(pw2.cld, aes(Site, emmean,label = "")) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE1, ymax = emmean + SE1),
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()
g2
#Calculates summary statistics for Sites in Mun 3
Mun3.means <- biod %>%
filter(Mun==3) %>%
group_by(Site) %>%
dplyr::summarize(Mean=mean(Diversity),
SD = sd(Diversity),
n = length(Diversity)) %>%
dplyr::mutate(SE = SD/sqrt(n))
pw3 <- emmeans(mod1,
pairwise ~ Site,
at = list(Mun="3"),
adjust="bonferroni")
#Assigns letters
pw3.cld <- cld(pw3,
Letters=letters,
decreasing = TRUE)
pw3.cld$SE1=Mun3.means$SE
pw3.cld## Site emmean SE df lower.CL upper.CL .group SE1
## 2 13 1.22 24 10.47 15.53 a 1.732
## 1 13 1.22 24 10.47 15.53 a 0.577
## 3 7 1.22 24 4.47 9.53 b 1.155
## 4 7 1.22 24 4.47 9.53 b 1.732
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 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.#Generates the bar plot with error bars for Mun=2
#Generates the bar plot with error bars for Mun=2
g3 <- ggplot(pw3.cld, aes(Site, emmean,label = .group)) +
geom_col(fill="gray") +
geom_errorbar(aes(ymin=emmean - SE1, ymax = emmean + SE1),
width = 0.2, size = 0.7) +
geom_text(vjust=-0.8,
hjust=1.5) +
lims(y = c(0,15)) +
labs(x= "Sites",
y= "Mean Biodiversity Index",
title = "Mun=3") +
theme_classic()
g3
#Removes y-axis
g2 <- g2 +
theme(
axis.text.y = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()
)
g3 <- g3 +
theme(
axis.text.y = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank()
)
#Combines the plots
g1 + g2 + g3 +
plot_layout(axis_titles = "collect") 