ANALISI CARABIDI - LAST UPDATE
Upload of dataset and packages
CSS <- read.csv("Cossogno_DB.csv")
CSS <- CSS[, -6:-13] # empty coloumns
CLL <- read.csv("Colloro_DB.csv")
CLL$habitat[CLL$habitat == "Forest"] <- "Wood"
CLL$habitat[CLL$habitat == "Ferns"] <- "Terraced ferns"
CLL$habitat[CLL$habitat == "Grassland"] <- "Terraced grassland"
library(dplyr)
library(ggplot2)
library(vegan)
library(tidyr)
library(corrr)
library(ggcorrplot)
library(FactoMineR)Sex ratio
sexdLL<- read.csv("Sesso_Colloro.csv")
sexdSS<- read.csv("Sesso_Cossogno.csv")
sexRatioLL<-sum(sexdLL$m)/sum(sexdLL$f)
sexRatioSS<-sum(sexdSS$m)/sum(sexdSS$f)
print(cat("Colloro M/F sex ratio: ", sexRatioLL, " | ", "Cossogno M/F sex ratio: ", sexRatioSS))## Colloro M/F sex ratio: 0.3225806 | Cossogno M/F sex ratio: 0.5014749NULLActivity Density
Activity Density (hereafter AD)
Formula: \[\text{AD} = \frac{\text{number of individuals}}{\text{number of traps}} * \frac{\text{number of sampling}}{\text{total days of fieldwork}}\]
# Number of species by year by site
## Cossogno
sppCSS_21 <- 24
sppCSS_22 <- 22
## Colloro
sppCLL_21 <- 12
sppCLL_22 <- 11
# Days of total fieldwork by year
ggtot_2021 <- 135
ggtot_2022 <- 120
# Sampling period each trap
gg_trap <- 15
# Numb of traps by site per habitat
trapNumb_CSS <- 15
trapNumb_CLL <- 4Somma individui per habitat per anno / num trappole per sito
CLL_HabYear <- CLL%>%
group_by(habitat, anno) %>%
summarise(sum(n_sp)/trapNumb_CLL)
CSS_HabYear <- CSS%>%
group_by(habitat, anno) %>%
summarise(sum(n_sp)/trapNumb_CSS)
print(CLL_HabYear)## # A tibble: 6 × 3
## # Groups: habitat [3]
## habitat anno `sum(n_sp)/trapNumb_CLL`
## <chr> <int> <dbl>
## 1 bosco 2021 6.5
## 2 bosco 2022 5.75
## 3 felci 2021 17
## 4 felci 2022 3
## 5 prato 2021 3.75
## 6 prato 2022 4.75print(CSS_HabYear)## # A tibble: 6 × 3
## # Groups: habitat [3]
## habitat anno `sum(n_sp)/trapNumb_CSS`
## <chr> <int> <dbl>
## 1 castagneto 2021 20.3
## 2 castagneto 2022 22.8
## 3 ecotono 2021 51.6
## 4 ecotono 2022 18.5
## 5 prato 2021 35.1
## 6 prato 2022 21.3AD_LL <- c(
AD_LL_Bosco21 <- 6.5 * (gg_trap/ggtot_2021),
AD_LL_Felci21 <- 17 * (gg_trap/ggtot_2021),
AD_LL_Prato21 <- 3.75 * (gg_trap/ggtot_2021),
AD_LL_Bosco22 <- 5.75 * (gg_trap/ggtot_2022),
AD_LL_Felci22 <- 3 * (gg_trap/ggtot_2022),
AD_LL_Prato22 <- 4.75 * (gg_trap/ggtot_2022)
)
AD_SS <- c(
AD_SS_Castagneto21 <- 20.27 * (gg_trap/ggtot_2021),
AD_SS_Ecotono21 <- 51.6 * (gg_trap/ggtot_2021),
AD_SS_Prato21 <- 35.07 * (gg_trap/ggtot_2021),
AD_SS_Castagento22 <- 22.8 * (gg_trap/ggtot_2022),
AD_SS_Ecotono22 <- 18.53 * (gg_trap/ggtot_2022),
AD_SS_Prato22 <- 21.33 * (gg_trap/ggtot_2022)
)# Colloro
sum_n_spLL <- CLL %>%
group_by(anno, habitat) %>%
summarise(total_n_spLL = sum(n_sp, na.rm = TRUE))
mean_se_n_spLL <- CLL %>%
group_by(anno, habitat, specie) %>%
summarise(mean_n_sp_per_specieLL = mean(n_sp, na.rm = TRUE)) %>%
group_by(anno, habitat) %>%
summarise(mean_n_spLL = mean(mean_n_sp_per_specieLL, na.rm = TRUE),
se_n_spLL = sd(mean_n_sp_per_specieLL, na.rm = TRUE) / sqrt(n()))
combined_resultsLL <- sum_n_spLL %>%
left_join(mean_se_n_spLL, by = c("anno", "habitat"))
print(combined_resultsLL)## # A tibble: 6 × 5
## # Groups: anno [2]
## anno habitat total_n_spLL mean_n_spLL se_n_spLL
## <int> <chr> <int> <dbl> <dbl>
## 1 2021 bosco 26 0.591 0.517
## 2 2021 felci 68 1.55 0.697
## 3 2021 prato 15 0.375 0.149
## 4 2022 bosco 23 0.523 0.450
## 5 2022 felci 12 0.273 0.104
## 6 2022 prato 19 0.432 0.143## Cossogno
sum_n_spSS <- CSS %>%
group_by(anno, habitat) %>%
summarise(total_n_spSS = sum(n_sp, na.rm = TRUE))
mean_se_n_spSS <- CSS %>%
group_by(anno, habitat, specie) %>%
summarise(mean_n_sp_per_specieSS = mean(n_sp, na.rm = TRUE)) %>%
group_by(anno, habitat) %>%
summarise(mean_n_spSS = mean(mean_n_sp_per_specieSS, na.rm = TRUE),
se_n_spSS = sd(mean_n_sp_per_specieSS, na.rm = TRUE) / sqrt(n()))
combined_resultsSS <- sum_n_spSS %>%
left_join(mean_se_n_spSS, by = c("anno", "habitat"))
print(combined_resultsSS)## # A tibble: 6 × 5
## # Groups: anno [2]
## anno habitat total_n_spSS mean_n_spSS se_n_spSS
## <int> <chr> <int> <dbl> <dbl>
## 1 2021 castagneto 304 0.844 0.270
## 2 2021 ecotono 774 2.15 1.02
## 3 2021 prato 526 1.46 0.527
## 4 2022 castagneto 342 1.04 0.455
## 5 2022 ecotono 278 0.842 0.369
## 6 2022 prato 320 0.970 0.399ADdf_LL <- data.frame(
Habitat = c("Wood", "Terraced ferns", "Terraced grassland", "Wood", "Terraced ferns", "Terraced grassland"),
Year = rep(c("2021", "2022"), each = 3),
DA = AD_LL,
SE = combined_resultsLL$se_n_spLL
)
ADdf_SS <- data.frame(
Habitat = c("Chestnut grove", "Ecotone", "Grassland", "Chestnut grove", "Ecotone", "Grassland"),
Year = rep(c("2021", "2022"), each = 3),
DA = AD_SS,
SE = combined_resultsSS$se_n_spSS
)# Colloro
barplot_AD_LL<- ggplot(ADdf_LL, aes(x = Habitat, y = DA, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.75), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
geom_errorbar(aes(ymin = DA - SE, ymax = DA + SE), position = position_dodge(width = 0.75), width = 0.1) +
ylim(0, 7) + ### Y axis limit from 0 - 7
labs(title = "(A)",
x = "",
y = "AD",
color = "Year") +
scale_x_discrete(labels = c("Terraced ferns" = "Terraced\nferns", "Terraced grassland" = "Terraced\ngrassland","Wood" = "Wood")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20, color = "black"),
panel.grid = element_blank(), # togliamo le griglie
panel.background = element_rect(fill = "white") # sfondo bianco
)
# Cossogno
barplot_AD_SS <- ggplot(ADdf_SS, aes(x = Habitat, y = DA, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.75), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
geom_errorbar(aes(ymin = DA - SE, ymax = DA + SE), position = position_dodge(width = 0.75), width = 0.1) +
ylim(0, 7) +
labs(title = "(B)",
x = "",
y = "AD",
color = "Year") +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20, color = "black"),
panel.grid = element_blank(), # togliamo le griglie
panel.background = element_rect(fill = "white") # sfondo bianco
)
barplot_AD_LL
barplot_AD_SS
Average number of species per trap
dd<-read.csv("ResumeCarabids2.csv")
dd$habitat[dd$habitat == "Forest"] <- "Wood"
dd$habitat[dd$habitat == "Ferns"] <- "Terraced Ferns"
## COLLORO
ddColl<-dd[dd$sito=="colloro",]
boxColl <- ggplot(ddColl, aes(x = habitat, y = n_specie, color = as.factor(anno))) +
geom_boxplot(position = position_dodge(width = 0.85)) +
ylim(0, 9) +
labs(title = "(A)",
x = "",
y = "Average n of species",
color = "Year") +
scale_x_discrete(labels = c("Wood" = "Wood", "Terraced Ferns" = "Terraced\nferns", "Grassland" = "Terraced\ngrassland")) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) + # colori per i diversi anni
theme_classic() +
theme(text=element_text(family="Arial", size = 20),
panel.grid = element_blank(), # togliamo le griglie
panel.background = element_rect(fill = "white")) # sfondo bianco)
boxColl
## COSSOGNO
ddCoss<-dd[dd$sito=="cicogna",]
boxCoss <- ggplot(ddCoss, aes(x = habitat, y = n_specie, color = as.factor(anno))) +
geom_boxplot(position = position_dodge(width = 0.85), fill = NA) +
ylim(0, 9) +
labs(title = "(B)",
x = "",
y = "Average n of species",
color = "Year") +
scale_x_discrete(labels = c("Chestnut_grove" = "Chestnut grove", "Ecotone" = "Ecotone", "Grassland" = "Grassland")) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) + # colori per i diversi anni
theme_classic()+
theme(text=element_text(family="Arial", size = 20),
panel.grid = element_blank(), # togliamo le griglie
panel.background = element_rect(fill = "white")) # sfondo bianco)
boxCoss
Biodiversity Indexes
- Create new database with total number of individual per species, divided by habitat type.
- Substitute variable habitat (type: string) with a numeric value (type: int). The values were
| Colloro | Cossogno | ||
|---|---|---|---|
| Grassland | 1 | Chestnut grove | 1 |
| Terraced ferns | 2 | Grassland | 2 |
| Wood | 3 | Ecotone | 3 |
- Load dataset
biodivDB_CLL21 <- read.csv("BiodivDB_Colloro21.csv")
biodivDB_CLL22 <- read.csv("BiodivDB_Colloro22.csv")
biodivDB_CSS21 <- read.csv("BiodivDB_Cossogno21.csv")
biodivDB_CSS22 <- read.csv("BiodivDB_Cossogno22.csv")Shannon-Wiener’s index using Vegan package
SH_CLL21 <- diversity(biodivDB_CLL21, index = "shannon")
SH_CLL22 <- diversity(biodivDB_CLL22, index = "shannon")
SH_CSS21 <- diversity(biodivDB_CSS21, index = "shannon")
SH_CSS22 <- diversity(biodivDB_CSS22, index = "shannon")
# Colloro
habitatCLL <- c("Terraced grassland", "Terraced ferns", "Wood")
dfSH_CLL <- data.frame(Habitat = rep(habitatCLL, 2),
Year = c(rep("2021", 3), rep("2022", 3)),
Shannon = c(SH_CLL21, SH_CLL22))
SH_CLL_plot <- ggplot(dfSH_CLL, aes(x = Habitat, y = Shannon, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.65), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
ylim(0, 2.5) +
labs(title = "(A)", x = "", y = "Shannon-Wiener index (H')", color = "Year") +
scale_x_discrete(labels = c("Wood" = "Wood", "Terraced ferns" = "Terraced\nferns", "Terraced grassland" = "Terraced\ngrassland")) +
theme_classic() + # Cambiamo il tema a classic per avere uno sfondo bianco
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
)
# Cossogno
habitatCSS <- c("Chestnut grove ", "Grassland", "Ecotone")
dfSH_CSS <- data.frame(Habitat = rep(habitatCSS, 2),
Year = c(rep("2021", 3), rep("2022", 3)),
Shannon = c(SH_CSS21, SH_CSS22))
SH_CSS_plot <- ggplot(dfSH_CSS, aes(x = Habitat, y = Shannon, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.65), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
ylim(0, 2.5) +
labs(title = "(A)", x = "", y = "Shannon-Wiener index (H')", color = "Year") +
theme_classic() + # Cambiamo il tema a classic per avere uno sfondo bianco
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
)
SH_CLL_plot
SH_CSS_plot
Simpson’s index using Vegan package
SP_CLL21 <- diversity(biodivDB_CLL21, index = "simpson")
SP_CLL22 <- diversity(biodivDB_CLL22, index = "simpson")
SP_CSS21 <- diversity(biodivDB_CSS21, index = "simpson")
SP_CSS22 <- diversity(biodivDB_CSS22, index = "simpson")
# Colloro
habitatCLL <- c("Terraced grassland", "Terraced ferns", "Wood")
dfSP_CLL <- data.frame(Habitat = rep(habitatCLL, 2),
Year = c(rep("2021", 3), rep("2022", 3)),
Simpson = c(SP_CLL21, SP_CLL22))
SP_CLL_plot <- ggplot(dfSP_CLL, aes(x = Habitat, y = Simpson, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.65), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
labs(title = "(B)", x = "", y = "Simpson index (D)", color = "Year") +
scale_x_discrete(labels = c("Wood" = "Wood", "Terraced ferns" = "Terraced\nferns", "Terraced grassland" = "Terraced\ngrassland")) +
ylim(0, 0.9) +
theme_classic() + # Cambiamo il tema a classic per avere uno sfondo bianco
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
)
# Cossogno
habitatCSS <- c("Chestnut grove ", "Grassland", "Ecotone")
dfSP_CSS <- data.frame(Habitat = rep(habitatCSS, 2),
Year = c(rep("2021", 3), rep("2022", 3)),
Simpson = c(SP_CSS21, SP_CSS22))
SP_CSS_plot <- ggplot(dfSP_CSS, aes(x = Habitat, y = Simpson, color = as.factor(Year))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.65), fill = NA, width = 0.5) +
scale_color_manual(values = c("2021" = "#33cc33", "2022" = "#0033cc")) +
labs(title = "(B)", x = "", y = "Simpson index (D)", color = "Year") +
ylim(0, 0.9) +
theme_classic() + # Cambiamo il tema a classic per avere uno sfondo bianco
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
)
SP_CLL_plot
SP_CSS_plot
Resuming the values of the previous calculation:
Colloro_BiodivIndex
Cossogno_BiodivIndex
NMDS
Non-metric multidimensional scaling (euclidean distances)
# Read in your data
NMDS_Colloro_x_R <- read.csv("CLL_NMDS.csv")
NMDS_Colloro_x_R$habitat[NMDS_Colloro_x_R$habitat == "Forest"] <- "Wood"
NMDS_Colloro_x_R$habitat[NMDS_Colloro_x_R$habitat == "Ferns"] <- "Terraced ferns"
NMDS_Colloro_x_R$habitat[NMDS_Colloro_x_R$habitat == "Grassland"] <- "Terraced grassland"
# Divide the file into two datasets, one with species numbers and one with habitat and trap number
data_1_CLL <- NMDS_Colloro_x_R[,3:18]
data_2_CLL <- NMDS_Colloro_x_R[,1:2]
# Perform NMDS
NMDS_CLL <- metaMDS(data_1_CLL, distance="euclidean", k=2)## Wisconsin double standardization
## Run 0 stress 0.09721812
## Run 1 stress 0.1037548
## Run 2 stress 0.0990684
## Run 3 stress 0.1071609
## Run 4 stress 0.3359232
## Run 5 stress 0.1331572
## Run 6 stress 0.104222
## Run 7 stress 0.09906769
## Run 8 stress 0.1071602
## Run 9 stress 0.118152
## Run 10 stress 0.1086629
## Run 11 stress 0.133157
## Run 12 stress 0.1071602
## Run 13 stress 0.1071602
## Run 14 stress 0.118152
## Run 15 stress 0.1122707
## Run 16 stress 0.09906882
## Run 17 stress 0.1526491
## Run 18 stress 0.1071602
## Run 19 stress 0.09721812
## ... Procrustes: rmse 1.393687e-05 max resid 2.576055e-05
## ... Similar to previous best
## Run 20 stress 0.1181522
## *** Best solution repeated 1 times# Extract NMDS scores (the coordinates) and combine with metadata
nmds_scores_CLL <- as.data.frame(scores(NMDS_CLL, display = "sites"))
nmds_data_CLL <- cbind(data_2_CLL, nmds_scores_CLL)
nmds_data_CLL$habitat[nmds_data_CLL$habitat == "Forest"] <- "Wood"
nmds_data_CLL$habitat[nmds_data_CLL$habitat == "Ferns"] <- "Terraced ferns"
nmds_data_CLL$habitat[nmds_data_CLL$habitat == "Grassland"] <- "Terraced grassland"
# Function to calculate convex hull for each group
find_hull_CLL <- function(df_CLL) df_CLL[chull(df_CLL$NMDS1, df_CLL$NMDS2), ]
# Calculate convex hull for each habitat
hulls_CLL <- nmds_data_CLL %>% group_by(habitat) %>% do(find_hull_CLL(.))
# ggplot(data_2_CLL, aes(x = NMDS_CLL$points[,1], y = NMDS_CLL$points[,2], color = habitat, shape = habitat)) +
# geom_point(size = 2) +
# labs(title = "", x = "NMDS1", y = "NMDS2", color = "Habitat", shape = "Habitat") +
# theme_classic() +
# theme(
# text = element_text(family = "Arial", size = 20),
# panel.grid = element_blank(), # Remove grid lines
# panel.background = element_rect(fill = "white") # White background
# )
# Plot NMDS with ggplot2
ggplot(nmds_data_CLL, aes(x = NMDS1, y = NMDS2, fill = habitat, color = habitat)) +
geom_point(size = 2) +
geom_polygon(data = hulls_CLL, aes(x = NMDS1, y = NMDS2, fill = habitat, color = habitat), alpha = 0.3) +
labs(title = "", x = "NMDS1", y = "NMDS2", color = "Habitat", fill = "Habitat") +
scale_color_manual(values = c("Terraced ferns" = "#f7b358", "Wood" = "#780116", "Terraced grassland" = "#c32f27")) +
scale_fill_manual(values = c("Terraced ferns" = "#f7b358", "Wood" = "#780116", "Terraced grassland" = "#c32f27")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Remove grid lines
panel.background = element_rect(fill = "white") # White background
)
# Read in your data
NMDS_Cossogno_x_R <- read.csv("CSS_NMDS.csv")
# Divide the file into two datasets, one with species numbers and one with habitat and trap number
data_1_CSS <- NMDS_Cossogno_x_R[,3:31]
data_2_CSS <- NMDS_Cossogno_x_R[,1:2]
# Perform NMDS
NMDS_CSS <- metaMDS(data_1_CSS, distance = "euclidean", k = 2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2119595
## Run 1 stress 0.2341513
## Run 2 stress 0.2268318
## Run 3 stress 0.2344022
## Run 4 stress 0.2320901
## Run 5 stress 0.2201072
## Run 6 stress 0.2184068
## Run 7 stress 0.2175828
## Run 8 stress 0.2229201
## Run 9 stress 0.2460664
## Run 10 stress 0.2237283
## Run 11 stress 0.2360303
## Run 12 stress 0.21974
## Run 13 stress 0.2290509
## Run 14 stress 0.2229655
## Run 15 stress 0.2203989
## Run 16 stress 0.2136137
## Run 17 stress 0.2325807
## Run 18 stress 0.2207021
## Run 19 stress 0.2304669
## Run 20 stress 0.2200323
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 2: no. of iterations >= maxit
## 18: stress ratio > sratmax# Extract NMDS scores for sites (the coordinates) and combine with metadata
nmds_scores_CSS <- as.data.frame(scores(NMDS_CSS, display = "sites"))
nmds_data_CSS <- cbind(data_2_CSS, nmds_scores_CSS)
# Rename the NMDS columns for clarity
names(nmds_data_CSS)[3:4] <- c("NMDS1", "NMDS2")
# Function to calculate convex hull for each group
find_hull_CSS <- function(df) df[chull(df$NMDS1, df$NMDS2), ]
# Calculate convex hull for each habitat
hulls_CSS <- nmds_data_CSS %>% group_by(habitat) %>% do(find_hull_CSS(.))
# Plot NMDS with ggplot2 (points only)
# ggplot(nmds_data_CSS, aes(x = NMDS1, y = NMDS2, color = habitat, shape = habitat)) +
# geom_point(size = 2) +
# labs(title = "", x = "NMDS1", y = "NMDS2", color = "Habitat", shape = "Habitat") +
# theme_classic() +
# theme(
# text = element_text(family = "Arial", size = 20),
# panel.grid = element_blank(), # Remove grid lines
# panel.background = element_rect(fill = "white") # White background
# )
# Plot NMDS with ggplot2 (points and convex hull polygons)
ggplot(nmds_data_CSS, aes(x = NMDS1, y = NMDS2, fill = habitat, color = habitat)) +
geom_point(size = 2) +
geom_polygon(data = hulls_CSS, aes(x = NMDS1, y = NMDS2, fill = habitat, color = habitat), alpha = 0.3) +
labs(title = "", x = "NMDS1", y = "NMDS2", color = "Habitat", shape = "Habitat") +
scale_color_manual(values = c("Chestnut grove" = "#fe7f2d", "Ecotone" = "#fcca46", "Grassland" = "#619b8a")) +
scale_fill_manual(values = c("Chestnut grove" = "#fe7f2d", "Ecotone" = "#fcca46", "Grassland" = "#619b8a")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Remove grid lines
panel.background = element_rect(fill = "white") # White background
) +
guides(color = guide_legend("Habitat"), fill = guide_legend("Habitat"))
Community’s Dominance Structure
Colloro Community’s Dominance Structure, divided by habitat type
Colloro_DB <- read.csv("Colloro_DB.csv")
# Calcola le percentuali delle specie per habitat per l'anno 2021
data_2021_CLL <- Colloro_DB %>%
filter(anno == 2021) %>%
group_by(habitat, specie) %>%
summarise(total_n_sp = sum(n_sp), .groups = 'drop') %>%
group_by(habitat) %>%
mutate(percentage = total_n_sp / sum(total_n_sp) * 100)
# Visualizza le percentuali per il 2021
print(data_2021_CLL)## # A tibble: 34 × 4
## # Groups: habitat [3]
## habitat specie total_n_sp percentage
## <chr> <chr> <int> <dbl>
## 1 bosco "abax_contractus" 23 88.5
## 2 bosco "amara_convexior" 0 0
## 3 bosco "calathus_cinctus " 1 3.85
## 4 bosco "harpalus_rufipalpis_rufipalpis" 0 0
## 5 bosco "harpalus_serripes_serripes" 0 0
## 6 bosco "harpalus_subcylindricus" 0 0
## 7 bosco "harpalus_tardus" 0 0
## 8 bosco "metallina_lampros " 1 3.85
## 9 bosco "pseudoophonus _calceatus" 0 0
## 10 bosco "pseudoophonus_griseus" 0 0
## # ℹ 24 more rows# Calcola le percentuali delle specie per habitat per l'anno 2022
data_2022_CLL <- Colloro_DB %>%
filter(anno == 2022) %>%
group_by(habitat, specie) %>%
summarise(total_n_sp = sum(n_sp), .groups = 'drop') %>%
group_by(habitat) %>%
mutate(percentage = total_n_sp / sum(total_n_sp) * 100)
# Visualizza le percentuali per il 2022
print(data_2022_CLL)## # A tibble: 33 × 4
## # Groups: habitat [3]
## habitat specie total_n_sp percentage
## <chr> <chr> <int> <dbl>
## 1 bosco "abax_contractus" 20 87.0
## 2 bosco "amara_convexior" 0 0
## 3 bosco "amara_fulvipes" 0 0
## 4 bosco "calathus_fuscipes_graecus" 0 0
## 5 bosco "carabus_convexus_convexus " 1 4.35
## 6 bosco "harpalus_subcylindricus" 0 0
## 7 bosco "harpalus_tardus" 0 0
## 8 bosco "pseudoophonus _calceatus" 0 0
## 9 bosco "pseudoophonus_griseus" 0 0
## 10 bosco "pseudoophonus_rufipes " 0 0
## # ℹ 23 more rowsCossogno Community’s Dominance Structure, divided by habitat type.
#Carico file Cossogno_DB
Cossogno_DB <- read.csv("Cossogno_DB.csv")
# Calcola le percentuali delle specie per habitat per l'anno 2021
data_2021_CSS <- Cossogno_DB %>%
filter(anno == 2021) %>%
group_by(habitat, specie) %>%
summarise(total_n_sp = sum(n_sp), .groups = 'drop') %>%
group_by(habitat) %>%
mutate(percentage = total_n_sp / sum(total_n_sp) * 100)
# Visualizza le percentuali per il 2021
print(data_2021_CSS)## # A tibble: 72 × 4
## # Groups: habitat [3]
## habitat specie total_n_sp percentage
## <chr> <chr> <int> <dbl>
## 1 castagneto abax_baenningeri 1 0.329
## 2 castagneto abax_contractus 23 7.57
## 3 castagneto amara_convexior 0 0
## 4 castagneto calathus_fuscipes_graecus 61 20.1
## 5 castagneto calathus_rubripes 1 0.329
## 6 castagneto calosoma_sycophanta 5 1.64
## 7 castagneto carabus_ convexus_convexus 29 9.54
## 8 castagneto carabus_ granulatus_interstitialis 2 0.658
## 9 castagneto carabus_germarii_fiorii 0 0
## 10 castagneto carabus_glabratus_latior 12 3.95
## # ℹ 62 more rows# Calcola le percentuali delle specie per habitat per l'anno 2022
data_2022_CSS <- Cossogno_DB %>%
filter(anno == 2022) %>%
group_by(habitat, specie) %>%
summarise(total_n_sp = sum(n_sp), .groups = 'drop') %>%
group_by(habitat) %>%
mutate(percentage = total_n_sp / sum(total_n_sp) * 100)
# Visualizza le percentuali per il 2022
print(data_2022_CSS)## # A tibble: 66 × 4
## # Groups: habitat [3]
## habitat specie total_n_sp percentage
## <chr> <chr> <int> <dbl>
## 1 castagneto abax_baenningeri 1 0.292
## 2 castagneto abax_contractus 18 5.26
## 3 castagneto amara_bifrons 0 0
## 4 castagneto amara_convexior 0 0
## 5 castagneto calathus_fuscipes_graecus 132 38.6
## 6 castagneto calathus_neocalathus_melanocephalus 1 0.292
## 7 castagneto carabus_ convexus_convexus 1 0.292
## 8 castagneto carabus_ granulatus_interstitialis 1 0.292
## 9 castagneto carabus_glabratus_latior 11 3.22
## 10 castagneto carabus_intricatus 11 3.22
## # ℹ 56 more rowsColloro Community’s Dominance Structure, without considering the habitat type
# Carico file Colloro_DB
Colloro_DB <- read.csv("Colloro_DB.csv")
# Dividi i dati per anno
data_2021_CLL_TOT <- Colloro_DB %>% filter(anno == 2021)
data_2022_CLL_TOT <- Colloro_DB %>% filter(anno == 2022)
# Conta le occorrenze delle specie per ogni anno
species_count_2021_CLL <- data_2021_CLL_TOT %>% group_by(specie) %>% summarise(count = sum(n_sp))
species_count_2022_CLL <- data_2022_CLL_TOT %>% group_by(specie) %>% summarise(count = sum(n_sp))
# Calcola la percentuale di presenza per ogni specie per l'anno 2021
total_2021_CLL <- sum(species_count_2021_CLL$count)
species_count_2021_CLL <- species_count_2021_CLL %>%
mutate(Percentuale = count / total_2021_CLL * 100)
# Calcola la percentuale di presenza per ogni specie per l'anno 2022
total_2022_CLL <- sum(species_count_2022_CLL$count)
species_count_2022_CLL <- species_count_2022_CLL %>%
mutate(Percentuale = count / total_2022_CLL * 100)
# Visualizza i risultati
print(species_count_2021_CLL, n= Inf)## # A tibble: 12 × 3
## specie count Percentuale
## <chr> <int> <dbl>
## 1 "abax_contractus" 34 31.2
## 2 "amara_convexior" 1 0.917
## 3 "calathus_cinctus " 1 0.917
## 4 "carabus_germari_fiorii" 1 0.917
## 5 "harpalus_rufipalpis_rufipalpis" 32 29.4
## 6 "harpalus_serripes_serripes" 8 7.34
## 7 "harpalus_subcylindricus" 6 5.50
## 8 "harpalus_tardus" 17 15.6
## 9 "metallina_lampros " 1 0.917
## 10 "pseudoophonus _calceatus" 6 5.50
## 11 "pseudoophonus_griseus" 1 0.917
## 12 "pterostichus_micans" 1 0.917print(species_count_2022_CLL, n= Inf)## # A tibble: 11 × 3
## specie count Percentuale
## <chr> <int> <dbl>
## 1 "abax_contractus" 28 51.9
## 2 "amara_convexior" 1 1.85
## 3 "amara_fulvipes" 1 1.85
## 4 "calathus_fuscipes_graecus" 1 1.85
## 5 "carabus_convexus_convexus " 2 3.70
## 6 "harpalus_subcylindricus" 1 1.85
## 7 "harpalus_tardus" 3 5.56
## 8 "pseudoophonus _calceatus" 4 7.41
## 9 "pseudoophonus_griseus" 5 9.26
## 10 "pseudoophonus_rufipes " 6 11.1
## 11 "pterostichus_micans " 2 3.70Cossogno Community’s Dominance Structure, without considering the habitat type
# Dividi i dati per anno
data_2021_CSS_TOT <- Cossogno_DB %>% filter(anno == 2021)
data_2022_CSS_TOT <- Cossogno_DB %>% filter(anno == 2022)
# Conta le occorrenze delle specie per ogni anno
species_count_2021_CSS <- data_2021_CSS_TOT %>% group_by(specie) %>% summarise(count = sum(n_sp))
species_count_2022_CSS <- data_2022_CSS_TOT %>% group_by(specie) %>% summarise(count = sum(n_sp))
# Calcola la percentuale di presenza per ogni specie per l'anno 2021
total_2021_CSS <- sum(species_count_2021_CSS$count)
species_count_2021_CSS <- species_count_2021_CSS %>%
mutate(Percentuale = count / total_2021_CSS * 100)
# Calcola la percentuale di presenza per ogni specie per l'anno 2022
total_2022_CSS <- sum(species_count_2022_CSS$count)
species_count_2022_CSS <- species_count_2022_CSS %>%
mutate(Percentuale = count / total_2022_CSS * 100)
# Visualizza i risultati
print(species_count_2021_CSS, n= Inf)## # A tibble: 24 × 3
## specie count Percentuale
## <chr> <int> <dbl>
## 1 abax_baenningeri 1 0.0623
## 2 abax_contractus 44 2.74
## 3 amara_convexior 9 0.561
## 4 calathus_fuscipes_graecus 406 25.3
## 5 calathus_rubripes 3 0.187
## 6 calosoma_sycophanta 10 0.623
## 7 carabus_ convexus_convexus 54 3.37
## 8 carabus_ granulatus_interstitialis 6 0.374
## 9 carabus_germarii_fiorii 1 0.0623
## 10 carabus_glabratus_latior 342 21.3
## 11 carabus_intricatus 30 1.87
## 12 carabus_problematicus_problematicus 292 18.2
## 13 cychrus_italicus 15 0.935
## 14 cymindis_cingulata 1 0.0623
## 15 harpalus_atratus 3 0.187
## 16 harpalus_harpalus_marginellus 1 0.0623
## 17 harpalus_rufipalpis_rufipalpis 38 2.37
## 18 harpalus_tardus 24 1.50
## 19 metallina_lampros 4 0.249
## 20 parophonus_parophonus_maculicornis 1 0.0623
## 21 pseudoophonus_griseus 74 4.61
## 22 pseudoophonus_rufipes 179 11.2
## 23 pterostichus_micans 63 3.93
## 24 synthomus_truncatellus 3 0.187print(species_count_2022_CSS, n= Inf)## # A tibble: 22 × 3
## specie count Percentuale
## <chr> <int> <dbl>
## 1 abax_baenningeri 1 0.106
## 2 abax_contractus 36 3.83
## 3 amara_bifrons 1 0.106
## 4 amara_convexior 2 0.213
## 5 calathus_fuscipes_graecus 323 34.4
## 6 calathus_neocalathus_melanocephalus 2 0.213
## 7 carabus_ convexus_convexus 11 1.17
## 8 carabus_ granulatus_interstitialis 3 0.319
## 9 carabus_glabratus_latior 70 7.45
## 10 carabus_intricatus 16 1.70
## 11 carabus_problematicus_problematicus 110 11.7
## 12 cychrus_italicus 2 0.213
## 13 cymindis_cingulata 2 0.213
## 14 harpalus_atratus 1 0.106
## 15 harpalus_rufipalpis_rufipalpis 25 2.66
## 16 harpalus_tardus 13 1.38
## 17 laemostenus_janthinus_coeruleus 1 0.106
## 18 limodromus_assimilis 1 0.106
## 19 pseudoophonus_griseus 103 11.0
## 20 pseudoophonus_rufipes 197 21.0
## 21 pterostichus_micans 19 2.02
## 22 stomis_roccae_roccae 1 0.106Print the data
write.csv2(data_2021_CSS, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CSS_hab21.csv")
write.csv2(data_2022_CSS, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CSS_hab22.csv")
write.csv2(data_2021_CLL, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CLL_hab21.csv")
write.csv2(data_2022_CLL, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CLL_hab22.csv")
write.csv2(species_count_2021_CSS, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CSS_tot21.csv")
write.csv2(species_count_2022_CSS, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CSS_tot22.csv")
write.csv2(species_count_2021_CLL, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CLL_tot21.csv")
write.csv2(species_count_2022_CLL, "/home/matteo/Documents/AngeliDISAFA/RandomAnalysis/Analisi_Serena/DATI_perc/CLL_tot22.csv")Dominance tables and plots
Dominance classes: Colloro
# Calcolo le 5 classi di dominanza per il 2021 di Colloro - uso il DF species_count_2021_CLL
# Definizione delle categorie in base alla percentuale
dominance_CLL_21 <- species_count_2021_CLL %>%
mutate(Categoria = case_when(
Percentuale > 10 ~ "Eudominants",
Percentuale > 5 & Percentuale <= 10 ~ "Dominants",
Percentuale > 2 & Percentuale <= 5 ~ "Subdominants",
Percentuale > 1 & Percentuale <= 2 ~ "Recedents",
Percentuale <= 1 ~ "Subrecedents"
))
dominance_CLL_21$Categoria <- factor(dominance_CLL_21$Categoria, levels = c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"))
# Conto delle specie in ciascuna categoria
dominance_count_CLL_21 <- dominance_CLL_21 %>% count(Categoria)
# Visualizza i risultati
print(dominance_count_CLL_21)## # A tibble: 3 × 2
## Categoria n
## <fct> <int>
## 1 Eudominants 3
## 2 Dominants 3
## 3 Subrecedents 6# Calcolo le 5 classi di dominanza per il 2022 di Colloro - uso il DF species_count_2022_CLL
# Definizione delle categorie in base alla percentuale
dominance_CLL_22 <- species_count_2022_CLL %>%
mutate(Categoria = case_when(
Percentuale > 10 ~ "Eudominants",
Percentuale > 5 & Percentuale <= 10 ~ "Dominants",
Percentuale > 2 & Percentuale <= 5 ~ "Subdominants",
Percentuale > 1 & Percentuale <= 2 ~ "Recedents",
Percentuale <= 1 ~ "Subrecedents"
))
dominance_CLL_22$Categoria <- factor(dominance_CLL_22$Categoria, levels = c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"))
# Conto delle specie in ciascuna categoria
dominance_count_CLL_22 <- dominance_CLL_22 %>% count(Categoria)
# Visualizza i risultati
print(dominance_count_CLL_22)## # A tibble: 4 × 2
## Categoria n
## <fct> <int>
## 1 Eudominants 2
## 2 Dominants 3
## 3 Subdominants 2
## 4 Recedents 4Dominance classes: Colloro
# Calcolo le 5 classi di dominanza per il 2021 di Cossogno - uso il DF species_count_2021_CSS
# Definizione delle categorie in base alla percentuale
dominance_CSS_21 <- species_count_2021_CSS %>%
mutate(Categoria = case_when(
Percentuale > 10 ~ "Eudominants",
Percentuale > 5 & Percentuale <= 10 ~ "Dominants",
Percentuale > 2 & Percentuale <= 5 ~ "Subdominants",
Percentuale > 1 & Percentuale <= 2 ~ "Recedents",
Percentuale <= 1 ~ "Subrecedents"
))
dominance_CSS_21$Categoria <- factor(dominance_CSS_21$Categoria, levels = c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"))
# Conto delle specie in ciascuna categoria
dominance_count_CSS_21 <- dominance_CSS_21 %>% count(Categoria)
# Visualizza i risultati
print(dominance_count_CSS_21)## # A tibble: 4 × 2
## Categoria n
## <fct> <int>
## 1 Eudominants 4
## 2 Subdominants 5
## 3 Recedents 2
## 4 Subrecedents 13# Calcolo le 5 classi di dominanza per il 2022 di Cossogno - uso il DF species_count_2022_CSS
# Definizione delle categorie in base alla percentuale
dominance_CSS_22 <- species_count_2022_CSS %>%
mutate(Categoria = case_when(
Percentuale > 10 ~ "Eudominants",
Percentuale > 5 & Percentuale <= 10 ~ "Dominants",
Percentuale > 2 & Percentuale <= 5 ~ "Subdominants",
Percentuale > 1 & Percentuale <= 2 ~ "Recedents",
Percentuale <= 1 ~ "Subrecedents"
))
dominance_CSS_22$Categoria <- factor(dominance_CSS_22$Categoria, levels = c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"))
# Conto delle specie in ciascuna categoria
dominance_count_CSS_22 <- dominance_CSS_22 %>% count(Categoria)
# Visualizza i risultati
print(dominance_count_CSS_22)## # A tibble: 5 × 2
## Categoria n
## <fct> <int>
## 1 Eudominants 4
## 2 Dominants 1
## 3 Subdominants 3
## 4 Recedents 3
## 5 Subrecedents 11torta_CLL_21 <- dominance_count_CLL_21 %>%
mutate(Percentuale = n / sum(n) * 100,
Label = paste0(Categoria, " (", round(Percentuale, 1), "%)"))
torta_CLL_21$Label <- factor(torta_CLL_21$Label, levels = paste0(c("Eudominants", "Subrecedents"), " (", round(torta_CLL_21$Percentuale, 1), "%)"))
# Definisci i colori personalizzati
custom_colors <- c("#003366", "#006699", "#99FFFF")
# Crea il grafico a torta
ggplot(torta_CLL_21, aes(x = "", y = n, fill = Categoria)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar("y", start = 0) +
scale_fill_manual(values = custom_colors, labels = paste0(c("Eudominants", "Dominants", "Subrecedents"), " (", round(torta_CLL_21$Percentuale, 1), "%)")) +
theme_void() +
theme(legend.position = "right") +
theme(legend.position = "right",
text = element_text(family = "arial", size = 20)) +
labs(title = "(A)", fill = "Legend")
# Colloro 2022
# Uso dominance_count_CLL_22
# Aggiungi una colonna con le etichette di percentuale
torta_CLL_22 <- dominance_count_CLL_22 %>%
mutate(Percentuale = n / sum(n) * 100,
Label = paste0(Categoria, " (", round(Percentuale, 1), "%)"))
torta_CLL_22$Label <- factor(torta_CLL_22$Label, levels = paste0(c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CLL_22$Percentuale, 1), "%)"))
# Definisci i colori personalizzati
custom_colors <- c("#003366", "#006699", "#3399CC", "#66CCFF", "#99FFFF")
# Crea il grafico a torta
ggplot(torta_CLL_22, aes(x = "", y = n, fill = Categoria)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar("y", start = 0) +
scale_fill_manual(values = custom_colors, labels = paste0(c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CLL_22$Percentuale, 1), "%)")) +
theme_void() +
theme(legend.position = "right") +
theme(legend.position = "right",
text = element_text(family = "arial", size = 20)) +
labs(title = "(B)", fill = "Legend")
torta_CSS_21 <- dominance_count_CSS_21 %>%
mutate(Percentuale_CSS21 = n / sum(n) * 100,
Label = paste0(Categoria, " (", round(Percentuale_CSS21, 1), "%)"))
torta_CSS_21$Label_CSS21 <- factor(torta_CSS_21$Label, levels = paste0(c("Eudominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CSS_21$Percentuale_CSS21, 1), "%)"))
# Definisci i colori personalizzati
custom_colors <- c("#003366", "#3399CC", "#66CCFF", "#99FFFF")
# Crea il grafico a torta
ggplot(torta_CSS_21, aes(x = "", y = n, fill = Categoria)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar("y", start = 0) +
scale_fill_manual(values = custom_colors, labels = paste0(c("Eudominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CSS_21$Percentuale_CSS21, 1), "%)")) +
theme_void() +
theme(legend.position = "right") +
theme(legend.position = "right",
text = element_text(family = "arial", size = 20)) +
labs(title = "(A)", fill = "Legend")
# Cossogno 2022
# Uso dominance_count_CSS_22
# Aggiungi una colonna con le etichette di percentuale
torta_CSS_22 <- dominance_count_CSS_22 %>%
mutate(Percentuale_CSS22 = n / sum(n) * 100,
Label = paste0(Categoria, " (", round(Percentuale_CSS22, 1), "%)"))
# Ordinare la colonna 'Label' secondo l'ordine specificato
torta_CSS_22$Label_CSS22 <- factor(torta_CSS_22$Label, levels = paste0(c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CSS_22$Percentuale_CSS22, 1), "%)"))
# Definisci i colori personalizzati
custom_colors <- c("#003366", "#006699", "#3399CC", "#66CCFF", "#99FFFF")
# Crea il grafico a torta
ggplot(torta_CSS_22, aes(x = "", y = n, fill = Categoria)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar("y", start = 0) +
scale_fill_manual(values = custom_colors, labels = paste0(c("Eudominants", "Dominants", "Subdominants", "Recedents", "Subrecedents"), " (", round(torta_CSS_22$Percentuale_CSS22, 1), "%)")) +
theme_void() +
theme(legend.position = "right",
text = element_text(family = "arial", size = 20)) +
labs(title = "(B)", fill = "Legend")
# PCA COLLORO
dh1<-read.csv("CLL_NMDS.csv")
dh1 <- dh1 %>%
group_by(habitat) %>%
summarise(across(everything(), sum))
numerical_data <- dh1[,3:18]
head(numerical_data)data_normalized <- scale(numerical_data)
head(data_normalized)## abax_continuus_continuus amara_aenea calathus_cinctus
## [1,] 1.1543149 -0.5773503 1.1547005
## [2,] -0.5513146 -0.5773503 -0.5773503
## [3,] -0.6030003 1.1547005 -0.5773503
## harpalus_rufipalpis_rufipalpis harpalus_serripes_serripes
## [1,] -0.6359429 -0.5773503
## [2,] -0.5167036 -0.5773503
## [3,] 1.1526465 1.1547005
## harpalus_subcylindricus harpalus_tardus metallina_lampros
## [1,] -0.9271726 -0.6234797 1.1547005
## [2,] 1.0596259 -0.5299577 -0.5773503
## [3,] -0.1324532 1.1534374 -0.5773503
## pseudoophonus._calceatus pseudoophonus_griseus pterostichus_micans
## [1,] -1.1547005 -1.1547005 1.1547005
## [2,] 0.5773503 0.5773503 -0.5773503
## [3,] 0.5773503 0.5773503 -0.5773503
## calathus_fuscipes.graecus carabus_convexus_convexus pseudoophonus_rufipes
## [1,] -0.5773503 0.5773503 -0.7559289
## [2,] 1.1547005 0.5773503 1.1338934
## [3,] -0.5773503 -1.1547005 -0.3779645
## amara_fulvipes carabus_germari_fiorii
## [1,] -0.5773503 -0.5773503
## [2,] -0.5773503 -0.5773503
## [3,] 1.1547005 1.1547005data.pca <- prcomp(data_normalized)
summary(data.pca)## Importance of components:
## PC1 PC2 PC3
## Standard deviation 3.1719 2.4370 1.554e-16
## Proportion of Variance 0.6288 0.3712 0.000e+00
## Cumulative Proportion 0.6288 1.0000 1.000e+00data.pca$loadings[, 1:2]## NULL# install.packages("factoextra")
library("factoextra")
fviz_eig(data.pca, addlabels = TRUE)
# Graph of the variables
fviz_pca_var(data.pca, col.var = "red")
fviz_pca_var(data.pca, col.var = "cos2",
gradient.cols = c("black", "orange", "green"),
repel = TRUE)
###
numerical_data_t<-t(numerical_data)
data<-numerical_data_t
# Eseguire PCA
pca <- prcomp(data, scale. = TRUE)
# Estrarre le coordinate dei variabili (habitat)
habitat_coords <- as.data.frame(pca$rotation[, 1:2])
habitat_coords$habitat <- c("Terraced ferns", "Wood", "Grassland") # Nomi degli habitat
# Visualizzare il biplot delle specie e degli habitat
biplot <- fviz_pca_biplot(pca,
geom.ind = "point", # Mostra le specie come punti
col.ind = "blue", # Colore delle specie
pointshape = 21, # Forma dei punti
pointsize = 2, # Dimensione dei punti
fill.ind = "blue", # Riempimento dei punti
repel = TRUE, # Evita la sovrapposizione delle etichette
geom.var = "arrow", # Mostra gli habitat come frecce
col.var = habitat_coords$habitat, # Colore delle frecce basato sugli habitat
legend.title = list(color = "Habitat")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
) +
labs(title = "PCA Biplot",
x = paste0("Dim1 (", round(pca$sdev[1]^2/sum(pca$sdev^2) * 100, 1), "%)"),
y = paste0("Dim2 (", round(pca$sdev[2]^2/sum(pca$sdev^2) * 100, 1), "%)"))
# Aggiungere le etichette dei variabili (habitat)
biplot + geom_text(data = habitat_coords, aes(x = PC1, y = PC2, label = habitat),
color = "black",
vjust = -1,
hjust = -0.5) +
scale_color_manual(values = c("Terraced ferns" = "green", "Wood" = "brown", "Grassland" = "orange"))
# Assegna un numero a ogni specie
species_numbers <- data.frame(Species = colnames(numerical_data), Number = 1:ncol(numerical_data))
# Visualizza la tabella di corrispondenza
print(species_numbers)## Species Number
## 1 abax_continuus_continuus 1
## 2 amara_aenea 2
## 3 calathus_cinctus 3
## 4 harpalus_rufipalpis_rufipalpis 4
## 5 harpalus_serripes_serripes 5
## 6 harpalus_subcylindricus 6
## 7 harpalus_tardus 7
## 8 metallina_lampros 8
## 9 pseudoophonus._calceatus 9
## 10 pseudoophonus_griseus 10
## 11 pterostichus_micans 11
## 12 calathus_fuscipes.graecus 12
## 13 carabus_convexus_convexus 13
## 14 pseudoophonus_rufipes 14
## 15 amara_fulvipes 15
## 16 carabus_germari_fiorii 16# Esegui la PCA
numerical_data_t <- t(numerical_data)
data <- numerical_data_t
pca <- prcomp(data, scale. = TRUE)
# Estrarre le coordinate delle specie
species_coords <- as.data.frame(pca$x[, 1:2])
species_coords$Number <- 1:ncol(numerical_data)
# Estrarre le coordinate degli habitat
habitat_coords <- as.data.frame(pca$rotation[, 1:2])
habitat_coords$habitat <- c("Terraced ferns", "Wood", "Terraced grassland") # Nomi degli habitat
# Visualizza il biplot con i numeri delle specie
biplot <- fviz_pca_biplot(pca,
geom.ind = "point", # Mostra le specie come punti
col.ind = "blue", # Colore delle specie
pointshape = 21, # Forma dei punti
pointsize = 2, # Dimensione dei punti
fill.ind = "blue", # Riempimento dei punti
repel = TRUE, # Evita la sovrapposizione delle etichette
geom.var = "arrow", # Mostra gli habitat come frecce
col.var = habitat_coords$habitat, # Colore delle frecce basato sugli habitat
legend.title = list(color = "Habitat")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
) +
labs(title = "PCA Biplot",
x = paste0("Dim1 (", round(pca$sdev[1]^2/sum(pca$sdev^2) * 100, 1), "%)"),
y = paste0("Dim2 (", round(pca$sdev[2]^2/sum(pca$sdev^2) * 100, 1), "%)"))
# Aggiungere i numeri delle specie come etichette
biplot + geom_text(data = species_coords, aes(x = PC1, y = PC2, label = Number),
color = "black",
vjust = -1,
hjust = -0.5) +
scale_color_manual(values = c("Terraced ferns" = "#f7b358", "Wood" = "#780116", "Terraced grassland" = "#c32f27"))
write.csv(species_numbers, "species_numbersCLL.csv", row.names = FALSE)# PCA COSSOGNO
dh2<-read.csv("CSS_NMDS.csv")
dh2 <- dh2 %>%
group_by(habitat) %>%
summarise(across(everything(), sum))
numerical_data <- dh2[,3:31]
head(numerical_data)data_normalized <- scale(numerical_data)
head(data_normalized)## calathus_fuscipes_graecus carabus_problematicus_problematicus
## [1,] -0.8288558 -0.2925700
## [2,] 1.1106668 1.1136537
## [3,] -0.2818110 -0.8210836
## carabus_.convexus_convexus abax_contractus cychrus_italicus
## [1,] 0.7473516 0.8470543 0.6757374
## [2,] 0.3886228 0.2560862 0.4730162
## [3,] -1.1359744 -1.1031405 -1.1487535
## harpalus_atratus pseudoophonus_rufipes pterostichus_micans
## [1,] 1.0910895 -0.1698493 1.1542353
## [2,] -0.8728716 -0.9041979 -0.5487348
## [3,] -0.2182179 1.0740472 -0.6055005
## carabus_glabratus_latior cymindis_cingulata carabus_intricatus
## [1,] -0.7477986 1.1547005 0.7258662
## [2,] 1.1358691 -0.5773503 0.4147807
## [3,] -0.3880704 -0.5773503 -1.1406469
## carabus_.granulatus_interstitialis pseudoophonus_griseus amara_convexior
## [1,] 0 -0.4588315 -0.7758802
## [2,] 1 -0.6882472 -0.3526728
## [3,] -1 1.1470787 1.1285530
## harpalus_tardus harpalus_rufipalpis_rufipalpis syntomus_truncatellus
## [1,] -0.6450859 -0.6059653 -0.5773503
## [2,] -0.5068532 -0.5482544 -0.5773503
## [3,] 1.1519392 1.1542197 1.1547005
## metallina_lampros calathus_rubripes calosoma_sycophanta abax_baenningeri
## [1,] -0.5773503 0 1.0910895 1.1547005
## [2,] 1.1547005 1 -0.2182179 -0.5773503
## [3,] -0.5773503 -1 -0.8728716 -0.5773503
## carabus_germari_fiorii limodromus_assimilis stomis_roccae_roccae
## [1,] -0.5773503 1.1547005 1.1547005
## [2,] -0.5773503 -0.5773503 -0.5773503
## [3,] 1.1547005 -0.5773503 -0.5773503
## amara_bifrons laemostenus_janthinus_coeruleus
## [1,] -0.5773503 -0.5773503
## [2,] -0.5773503 -0.5773503
## [3,] 1.1547005 1.1547005
## calathus_neocalathus_melanocephalus parophonus_parophonus_maculicornis
## [1,] 0.5773503 -0.5773503
## [2,] -1.1547005 -0.5773503
## [3,] 0.5773503 1.1547005
## harpalus_harpalus_marginellus
## [1,] -0.5773503
## [2,] -0.5773503
## [3,] 1.1547005data.pca <- prcomp(data_normalized)
summary(data.pca)## Importance of components:
## PC1 PC2 PC3
## Standard deviation 4.3876 3.1224 3.14e-16
## Proportion of Variance 0.6638 0.3362 0.00e+00
## Cumulative Proportion 0.6638 1.0000 1.00e+00data.pca$loadings[, 1:2]## NULLfviz_eig(data.pca, addlabels = TRUE)
# Graph of the variables
fviz_pca_var(data.pca, col.var = "red")
fviz_pca_var(data.pca, col.var = "cos2",
gradient.cols = c("black", "orange", "green"),
repel = TRUE)
###
numerical_data_t<-t(numerical_data)
data<-numerical_data_t
# Eseguire PCA
pca <- prcomp(data, scale. = TRUE)
# Estrarre le coordinate dei variabili (habitat)
habitat_coords <- as.data.frame(pca$rotation[, 1:2])
habitat_coords$habitat <- c("Chersnut Grove", "Ecotone", "Grassland") # Nomi degli habitat
# Visualizzare il biplot delle specie e degli habitat
biplot <- fviz_pca_biplot(pca,
geom.ind = "point", # Mostra le specie come punti
col.ind = "blue", # Colore delle specie
pointshape = 21, # Forma dei punti
pointsize = 2, # Dimensione dei punti
fill.ind = "blue", # Riempimento dei punti
repel = TRUE, # Evita la sovrapposizione delle etichette
geom.var = "arrow", # Mostra gli habitat come frecce
col.var = habitat_coords$habitat, # Colore delle frecce basato sugli habitat
legend.title = list(color = "Habitat")) +
theme_minimal() +
labs(title = "PCA Biplot",
x = paste0("Dim1 (", round(pca$sdev[1]^2/sum(pca$sdev^2) * 100, 1), "%)"),
y = paste0("Dim2 (", round(pca$sdev[2]^2/sum(pca$sdev^2) * 100, 1), "%)"))
# Aggiungere le etichette dei variabili (habitat)
biplot + geom_text(data = habitat_coords, aes(x = PC1, y = PC2, label = habitat),
color = "black",
vjust = -1,
hjust = -0.5) +
scale_color_manual(values = c("Chersnut Grove" = "green", "Ecotone" = "brown", "Grassland" = "orange"))
# Assegna un numero a ogni specie
species_numbers <- data.frame(Species = colnames(numerical_data), Number = 1:ncol(numerical_data))
# Visualizza la tabella di corrispondenza
print(species_numbers)## Species Number
## 1 calathus_fuscipes_graecus 1
## 2 carabus_problematicus_problematicus 2
## 3 carabus_.convexus_convexus 3
## 4 abax_contractus 4
## 5 cychrus_italicus 5
## 6 harpalus_atratus 6
## 7 pseudoophonus_rufipes 7
## 8 pterostichus_micans 8
## 9 carabus_glabratus_latior 9
## 10 cymindis_cingulata 10
## 11 carabus_intricatus 11
## 12 carabus_.granulatus_interstitialis 12
## 13 pseudoophonus_griseus 13
## 14 amara_convexior 14
## 15 harpalus_tardus 15
## 16 harpalus_rufipalpis_rufipalpis 16
## 17 syntomus_truncatellus 17
## 18 metallina_lampros 18
## 19 calathus_rubripes 19
## 20 calosoma_sycophanta 20
## 21 abax_baenningeri 21
## 22 carabus_germari_fiorii 22
## 23 limodromus_assimilis 23
## 24 stomis_roccae_roccae 24
## 25 amara_bifrons 25
## 26 laemostenus_janthinus_coeruleus 26
## 27 calathus_neocalathus_melanocephalus 27
## 28 parophonus_parophonus_maculicornis 28
## 29 harpalus_harpalus_marginellus 29# Esegui la PCA
numerical_data_t <- t(numerical_data)
data <- numerical_data_t
pca <- prcomp(data, scale. = TRUE)
# Estrarre le coordinate delle specie
species_coords <- as.data.frame(pca$x[, 1:2])
species_coords$Number <- 1:ncol(numerical_data)
# Estrarre le coordinate degli habitat
habitat_coords <- as.data.frame(pca$rotation[, 1:2])
habitat_coords$habitat <- c("Chestnut Grove", "Ecotone", "Grassland") # Nomi degli habitat
# Visualizza il biplot con i numeri delle specie
biplot <- fviz_pca_biplot(pca,
geom.ind = "point", # Mostra le specie come punti
col.ind = "blue", # Colore delle specie
pointshape = 21, # Forma dei punti
pointsize = 2, # Dimensione dei punti
fill.ind = "blue", # Riempimento dei punti
repel = TRUE, # Evita la sovrapposizione delle etichette
geom.var = "arrow", # Mostra gli habitat come frecce
col.var = habitat_coords$habitat, # Colore delle frecce basato sugli habitat
legend.title = list(color = "Habitat")) +
theme_classic() +
theme(
text = element_text(family = "Arial", size = 20),
panel.grid = element_blank(), # Togliamo le griglie
panel.background = element_rect(fill = "white") # Sfondo bianco
) +
labs(title = "PCA Biplot",
x = paste0("Dim1 (", round(pca$sdev[1]^2/sum(pca$sdev^2) * 100, 1), "%)"),
y = paste0("Dim2 (", round(pca$sdev[2]^2/sum(pca$sdev^2) * 100, 1), "%)"))
# Aggiungere i numeri delle specie come etichette
biplot + geom_text(data = species_coords, aes(x = PC1, y = PC2, label = Number),
color = "black",
vjust = -1,
hjust = -0.5,
family = "Arial",
size = 5) +
scale_color_manual(values = c("Chestnut Grove" = "#fe7f2d", "Ecotone" = "#fcca46", "Grassland" = "#619b8a"))
write.csv(species_numbers, "species_numbersCSS.csv", row.names = FALSE)PCA_CSS<-read.csv("CSS_NMDS.csv")
dati_specie2 <- as.matrix(PCA_CSS[, 3:31])
# Esegui la PCA
pca_result2 <- rda(dati_specie2)
# Visualizza i risultati
summary(pca_result2)##
## Call:
## rda(X = dati_specie2)
##
## Partitioning of variance:
## Inertia Proportion
## Total 884.3 1
## Unconstrained 884.3 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Eigenvalue 424.1418 218.1065 106.9454 65.85961 40.31123 11.60781
## Proportion Explained 0.4796 0.2466 0.1209 0.07447 0.04558 0.01313
## Cumulative Proportion 0.4796 0.7262 0.8472 0.92165 0.96723 0.98036
## PC7 PC8 PC9 PC10 PC11 PC12
## Eigenvalue 5.806393 3.998147 3.230927 1.513383 1.115133 0.4414197
## Proportion Explained 0.006566 0.004521 0.003653 0.001711 0.001261 0.0004992
## Cumulative Proportion 0.986924 0.991445 0.995099 0.996810 0.998071 0.9985699
## PC13 PC14 PC15 PC16 PC17
## Eigenvalue 0.3350912 0.3119033 0.1780317 0.1050535 7.584e-02
## Proportion Explained 0.0003789 0.0003527 0.0002013 0.0001188 8.576e-05
## Cumulative Proportion 0.9989489 0.9993015 0.9995029 0.9996217 9.997e-01
## PC18 PC19 PC20 PC21 PC22
## Eigenvalue 0.0628747 5.681e-02 0.0435061 3.107e-02 1.878e-02
## Proportion Explained 0.0000711 6.424e-05 0.0000492 3.514e-05 2.124e-05
## Cumulative Proportion 0.9997785 9.998e-01 0.9998920 9.999e-01 9.999e-01
## PC23 PC24 PC25 PC26 PC27
## Eigenvalue 1.792e-02 1.446e-02 0.0091942 3.061e-03 1.060e-03
## Proportion Explained 2.027e-05 1.635e-05 0.0000104 3.461e-06 1.199e-06
## Cumulative Proportion 1.000e+00 1.000e+00 0.9999953 1.000e+00 1.000e+00# RDA
# Estrai la matrice dei dati delle specie
dati_specie2 <- as.matrix(PCA_CSS[, 3:31])
# Prepara i dati ambientali (nome degli habitat)
dati_ambientali2 <- as.factor(PCA_CSS[, 1]) # Assicurati che la prima colonna sia il nome degli habitat
# Esegui l'RDA
rda_result2 <- rda(dati_specie2 ~ dati_ambientali2)
# Visualizza i risultati
summary(rda_result2)##
## Call:
## rda(formula = dati_specie2 ~ dati_ambientali2)
##
## Partitioning of variance:
## Inertia Proportion
## Total 884.3 1
## Constrained 884.3 1
## Unconstrained 0.0 0
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 424.1418 218.1065 106.9454 65.85961 40.31123 11.60781
## Proportion Explained 0.4796 0.2466 0.1209 0.07447 0.04558 0.01313
## Cumulative Proportion 0.4796 0.7262 0.8472 0.92165 0.96723 0.98036
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## Eigenvalue 5.806393 3.998147 3.230927 1.513383 1.115133 0.4414197
## Proportion Explained 0.006566 0.004521 0.003653 0.001711 0.001261 0.0004992
## Cumulative Proportion 0.986924 0.991445 0.995099 0.996810 0.998071 0.9985699
## RDA13 RDA14 RDA15 RDA16 RDA17
## Eigenvalue 0.3350912 0.3119033 0.1780317 0.1050535 7.584e-02
## Proportion Explained 0.0003789 0.0003527 0.0002013 0.0001188 8.576e-05
## Cumulative Proportion 0.9989489 0.9993015 0.9995029 0.9996217 9.997e-01
## RDA18 RDA19 RDA20 RDA21 RDA22
## Eigenvalue 0.0628747 5.681e-02 0.0435061 3.107e-02 1.878e-02
## Proportion Explained 0.0000711 6.424e-05 0.0000492 3.514e-05 2.124e-05
## Cumulative Proportion 0.9997785 9.998e-01 0.9998920 9.999e-01 9.999e-01
## RDA23 RDA24 RDA25 RDA26 RDA27
## Eigenvalue 1.792e-02 1.446e-02 0.0091942 3.061e-03 1.060e-03
## Proportion Explained 2.027e-05 1.635e-05 0.0000104 3.461e-06 1.199e-06
## Cumulative Proportion 1.000e+00 1.000e+00 0.9999953 1.000e+00 1.000e+00
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 424.1418 218.1065 106.9454 65.85961 40.31123 11.60781
## Proportion Explained 0.4796 0.2466 0.1209 0.07447 0.04558 0.01313
## Cumulative Proportion 0.4796 0.7262 0.8472 0.92165 0.96723 0.98036
## RDA7 RDA8 RDA9 RDA10 RDA11 RDA12
## Eigenvalue 5.806393 3.998147 3.230927 1.513383 1.115133 0.4414197
## Proportion Explained 0.006566 0.004521 0.003653 0.001711 0.001261 0.0004992
## Cumulative Proportion 0.986924 0.991445 0.995099 0.996810 0.998071 0.9985699
## RDA13 RDA14 RDA15 RDA16 RDA17
## Eigenvalue 0.3350912 0.3119033 0.1780317 0.1050535 7.584e-02
## Proportion Explained 0.0003789 0.0003527 0.0002013 0.0001188 8.576e-05
## Cumulative Proportion 0.9989489 0.9993015 0.9995029 0.9996217 9.997e-01
## RDA18 RDA19 RDA20 RDA21 RDA22
## Eigenvalue 0.0628747 5.681e-02 0.0435061 3.107e-02 1.878e-02
## Proportion Explained 0.0000711 6.424e-05 0.0000492 3.514e-05 2.124e-05
## Cumulative Proportion 0.9997785 9.998e-01 0.9998920 9.999e-01 9.999e-01
## RDA23 RDA24 RDA25 RDA26 RDA27
## Eigenvalue 1.792e-02 1.446e-02 0.0091942 3.061e-03 1.060e-03
## Proportion Explained 2.027e-05 1.635e-05 0.0000104 3.461e-06 1.199e-06
## Cumulative Proportion 1.000e+00 1.000e+00 0.9999953 1.000e+00 1.000e+00# Grafico di dispersione degli habitat
plot(rda_result2, display = "sites")
# Grafico di dispersione delle specie
plot(rda_result2, display = "species")
# Estrai i risultati della PCA
pca_scores2 <- scores(pca_result2, display = "species")
pca_vars2 <- scores(pca_result2, display = "sites")
# Crea un biplot per la PCA
pca_biplot2 <- ggplot() +
geom_point(data = pca_vars2, aes(x = PC1, y = PC2), color = "blue", size = 3) +
geom_text(data = pca_vars2, aes(x = PC1, y = PC2, label = rownames(pca_vars2)), hjust = 1.1, vjust = 1.1, size = 3) +
geom_point(data = pca_scores2, aes(x = PC1, y = PC2), color = "red", size = 3) +
geom_text(data = pca_scores2, aes(x = PC1, y = PC2, label = rownames(pca_scores2)), hjust = -0.1, vjust = -0.1, size = 3) +
xlab(paste("PC1 (", round(100 * summary(pca_result2)$cont$importance[2, 1], digits = 1), "% variance)")) +
ylab(paste("PC2 (", round(100 * summary(pca_result2)$cont$importance[2, 2], digits = 1), "% variance)")) +
ggtitle("Biplot della PCA") +
theme_classic()
# Mostra il biplot della PCA
print(pca_biplot2)
# RDA
# Estrai i risultati dell'RDA
rda_scores2 <- scores(rda_result2, display = "sites")
# Crea un grafico di dispersione per l'RDA
rda_plot2 <- ggplot(data = rda_scores2, aes(x = RDA1, y = RDA2, color = dati_ambientali2)) +
geom_point(size = 3) +
xlab(paste("RDA Axis 1 (", round(100 * rda_result2$CA$eig[1]/sum(rda_result2$CA$eig), digits = 1), "% variance explained)")) +
ylab(paste("RDA Axis 2 (", round(100 * rda_result2$CA$eig[2]/sum(rda_result2$CA$eig), digits = 1), "% variance explained)")) +
ggtitle("Grafico di dispersione per l'RDA") +
theme_minimal() +
theme(legend.position = "right")
# Mostra il grafico di dispersione per l'RDA
print(rda_plot2)
- e n d -