pacman::p_load(pacman, tidyverse)
# Read the dataset
data <- read_csv("Aridity_Gradient_Israel.csv")
## Rows: 420 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Species
## dbl (3): MAP, Individual, FA_index
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Print the first few rows of the dataset
head(data)
## # A tibble: 6 × 4
## Species MAP Individual FA_index
## <chr> <dbl> <dbl> <dbl>
## 1 Anagallis arvensis 780 1 0.02
## 2 Anagallis arvensis 780 2 0.031
## 3 Anagallis arvensis 780 3 0.078
## 4 Anagallis arvensis 780 4 0.032
## 5 Anagallis arvensis 780 5 0.01
## 6 Anagallis arvensis 780 6 0.006
# Check for missing values and handle them (e.g., removing rows with NA in FA_index)
data_clean <- data %>%
filter(!is.na(FA_index))
# Summary statistics for FA_index across different species and MAP levels
summary_stats <- data_clean %>%
group_by(Species, MAP) %>%
summarize(
mean_FA = mean(FA_index, na.rm = TRUE),
sd_FA = sd(FA_index, na.rm = TRUE),
n = n()
)
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.
# Display the complete summary statistics
options(tibble.print_max = Inf, tibble.width = Inf)
print(summary_stats)
## # A tibble: 42 × 5
## # Groups: Species [12]
## Species MAP mean_FA sd_FA n
## <chr> <dbl> <dbl> <dbl> <int>
## 1 Anagallis arvensis 300 0.0194 0.0128 10
## 2 Anagallis arvensis 371 0.0251 0.0245 10
## 3 Anagallis arvensis 540 0.0224 0.0143 10
## 4 Anagallis arvensis 780 0.0281 0.0222 8
## 5 Erodium malacoides 300 0.0294 0.00974 10
## 6 Erodium malacoides 371 0.0455 0.0451 10
## 7 Erodium malacoides 540 0.0414 0.0460 9
## 8 Hedypnois rhagadioloides 300 0.0255 0.0165 10
## 9 Hedypnois rhagadioloides 371 0.0234 0.0130 10
## 10 Hedypnois rhagadioloides 540 0.0217 0.0115 10
## 11 Hedypnois rhagadioloides 780 0.0115 0.00859 10
## 12 Helianthemum salicifolium 300 0.0298 0.0286 10
## 13 Helianthemum salicifolium 371 0.0198 0.0189 10
## 14 Helianthemum salicifolium 540 0.0329 0.0305 10
## 15 Hymenocarpos circinnatus 300 0.0364 0.0256 10
## 16 Hymenocarpos circinnatus 371 0.0515 0.0504 10
## 17 Hymenocarpos circinnatus 540 0.0449 0.0337 10
## 18 Hymenocarpos circinnatus 780 0.0392 0.0194 10
## 19 Lotus peregrinus 371 0.0424 0.0315 10
## 20 Lotus peregrinus 540 0.0434 0.0277 10
## 21 Lotus peregrinus 780 0.0356 0.0233 10
## 22 Picris galilaea 300 0.0198 0.0147 10
## 23 Picris galilaea 371 0.0236 0.0171 10
## 24 Picris galilaea 540 0.0189 0.00840 10
## 25 Picris galilaea 780 0.0254 0.0132 10
## 26 Pterocephalus brevis 300 0.0223 0.0184 10
## 27 Pterocephalus brevis 371 0.0211 0.0179 10
## 28 Pterocephalus brevis 540 0.0248 0.0196 10
## 29 Pterocephalus brevis 780 0.0234 0.0185 9
## 30 Scorpiurus muricatus 300 0.0236 0.0215 10
## 31 Scorpiurus muricatus 540 0.0323 0.0268 10
## 32 Scorpiurus muricatus 780 0.0248 0.0139 10
## 33 Theligonum cynocrambe 371 0.0301 0.0204 10
## 34 Theligonum cynocrambe 540 0.0309 0.0202 10
## 35 Theligonum cynocrambe 780 0.0279 0.0217 10
## 36 Trifolium stellatum 371 0.0371 0.0363 10
## 37 Trifolium stellatum 540 0.0174 0.0138 10
## 38 Trifolium stellatum 780 0.0171 0.0139 10
## 39 Urospermum picroides 300 0.0188 0.0146 10
## 40 Urospermum picroides 371 0.0207 0.0218 10
## 41 Urospermum picroides 540 0.0186 0.0136 10
## 42 Urospermum picroides 780 0.0199 0.0195 10
# Visualise FA_index across different MAP levels for each species
ggplot(data_clean, aes(x = as.factor(MAP), y = FA_index, fill = Species)) +
geom_boxplot() +
facet_wrap(~ Species, scales = "free_y") +
labs(
title = "Fluctuating Asymmetry Index Across Different MAP Levels",
x = "Mean Annual Precipitation (MAP) in mm",
y = "Fluctuating Asymmetry Index (FA_index)"
) +
theme_minimal()
# Scatter plot to show the relationship between MAP and FA_index for each species
ggplot(data_clean, aes(x = MAP, y = FA_index, color = Species)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~ Species, scales = "free") +
labs(
title = "Relationship Between MAP and FA_index for Each Species",
x = "Mean Annual Precipitation (MAP) in mm",
y = "Fluctuating Asymmetry Index (FA_index)"
) +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
# Linear models to analyse the effect of MAP on FA_index for each species
lm_results <- data_clean %>%
group_by(Species) %>%
do(model = lm(FA_index ~ MAP, data = .)) %>%
summarize(
Species,
intercept = coef(model)[1],
slope = coef(model)[2],
p_value = summary(model)$coefficients[2, 4]
)
# Display the complete linear model results
options(tibble.print_max = Inf, tibble.width = Inf)
print(lm_results)
## # A tibble: 12 × 4
## Species intercept slope p_value
## <chr> <dbl> <dbl> <dbl>
## 1 Anagallis arvensis 0.0171 0.0000134 0.433
## 2 Erodium malacoides 0.0236 0.0000379 0.590
## 3 Hedypnois rhagadioloides 0.0345 -0.0000282 0.0119
## 4 Helianthemum salicifolium 0.0176 0.0000246 0.612
## 5 Hymenocarpos circinnatus 0.0458 -0.00000560 0.848
## 6 Lotus peregrinus 0.0504 -0.0000177 0.556
## 7 Picris galilaea 0.0181 0.00000777 0.508
## 8 Pterocephalus brevis 0.0207 0.00000442 0.784
## 9 Scorpiurus muricatus 0.0255 0.00000256 0.898
## 10 Theligonum cynocrambe 0.0329 -0.00000590 0.792
## 11 Trifolium stellatum 0.0497 -0.0000458 0.0909
## 12 Urospermum picroides 0.0193 0.000000477 0.974