pacman::p_load(pacman, ggplot2, dplyr, tidyr)
# Load the data
data <- read.csv("Greenhouse_experiment.csv")
# Display the first few rows of the data to check the structure
head(data)
## Species Watering_treatments Locality Individual
## 1 Hedypnois rhagadioloides 30 Mediterranean 1
## 2 Hedypnois rhagadioloides 30 Mediterranean 2
## 3 Hedypnois rhagadioloides 30 Mediterranean 3
## 4 Hedypnois rhagadioloides 30 Mediterranean 4
## 5 Hedypnois rhagadioloides 30 Semi-arid 1
## 6 Hedypnois rhagadioloides 30 Semi-arid 2
## FA_Index
## 1 0.086509812
## 2 0.030443128
## 3 0.029132707
## 4 0.071766978
## 5 0.025756231
## 6 0.009048129
# Convert relevant columns to factors
data$Species <- as.factor(data$Species)
data$Watering_treatments <- as.factor(data$Watering_treatments)
data$Locality <- as.factor(data$Locality)
# Summary statistics
summary_stats <- data %>%
group_by(Species, Watering_treatments, Locality) %>%
summarise(
count = n(),
mean_FA_Index = mean(FA_Index, na.rm = TRUE),
sd_FA_Index = sd(FA_Index, na.rm = TRUE)
)
## `summarise()` has grouped output by 'Species', 'Watering_treatments'. You can
## override using the `.groups` argument.
print(summary_stats)
## # A tibble: 12 × 6
## # Groups: Species, Watering_treatments [6]
## Species Watering_treatments Locality count mean_FA_Index sd_FA_Index
## <fct> <fct> <fct> <int> <dbl> <dbl>
## 1 Hedypnois rhaga… 30 Mediter… 4 0.0545 0.0291
## 2 Hedypnois rhaga… 30 Semi-ar… 6 0.0324 0.0323
## 3 Hedypnois rhaga… 50 Mediter… 10 0.0336 0.0211
## 4 Hedypnois rhaga… 50 Semi-ar… 10 0.0331 0.0115
## 5 Hedypnois rhaga… 90 Mediter… 10 0.0366 0.0548
## 6 Hedypnois rhaga… 90 Semi-ar… 9 0.0197 0.0223
## 7 Hymenocarpos ci… 30 Mediter… 10 0.0253 0.0149
## 8 Hymenocarpos ci… 30 Semi-ar… 9 0.0385 0.0247
## 9 Hymenocarpos ci… 50 Mediter… 10 0.0455 0.0430
## 10 Hymenocarpos ci… 50 Semi-ar… 10 0.0302 0.0223
## 11 Hymenocarpos ci… 90 Mediter… 10 0.0319 0.0247
## 12 Hymenocarpos ci… 90 Semi-ar… 10 0.0250 0.0270
# Perform a three-way ANOVA
anova_result <- aov(FA_Index ~ Species * Watering_treatments * Locality, data = data)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 1 0.00001 0.0000069 0.008 0.929
## Watering_treatments 2 0.00116 0.0005822 0.665 0.517
## Locality 1 0.00116 0.0011561 1.320 0.253
## Species:Watering_treatments 2 0.00091 0.0004562 0.521 0.596
## Species:Locality 1 0.00042 0.0004153 0.474 0.493
## Watering_treatments:Locality 2 0.00060 0.0003009 0.344 0.710
## Species:Watering_treatments:Locality 2 0.00247 0.0012349 1.410 0.249
## Residuals 96 0.08405 0.0008756
# Boxplot for FA_Index by Species and Watering_treatments
ggplot(data, aes(x = Watering_treatments, y = FA_Index, fill = Species)) +
geom_boxplot() +
facet_wrap(~ Locality) +
labs(title = "FA Index by Watering Treatments, Species, and Locality",
x = "Watering Treatments in ml",
y = "FA Index") +
theme_minimal()

# Interaction plots to see the interaction effects
interaction.plot(data$Watering_treatments, data$Species, data$FA_Index,
col = c("red", "blue"),
xlab = "Watering Treatments in ml",
ylab = "Mean FA Index",
trace.label = "Species")

interaction.plot(data$Watering_treatments, data$Locality, data$FA_Index,
col = c("darkgreen", "orange"),
xlab = "Watering Treatments in ml",
ylab = "Mean FA Index",
trace.label = "Locality")

interaction.plot(data$Species, data$Locality, data$FA_Index,
col = c("purple", "brown"),
xlab = "Species",
ylab = "Mean FA Index",
trace.label = "Locality")
