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")