── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks terra::extract()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(viridis)
Warning: package 'viridis' was built under R version 4.3.3
Biome data includes WWF biome classification, and ecoregion type. The layer was clipped to only include biomes 6 (Boreal) and 11 (Tundra) in arcgis pro. Then the shapefiles are exported into R for analysis.
The data for MSA is derived from the ‘MSA tool’ (Abrosio et al., 2024). The output was a shapefile for the whole world, and in Arcgis pro, only the areas that fall under the biomes defined above are included. This is raster data in the form of a tif file, which is then imported into R.
Warning: The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
ggsave("msa_overtime.svg", width =14, height =6)
Warning: The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
#pressure filepressure_data <-read_csv("SSP2/pressure and region/gpotveg_pressure_impact_plants.csv")
Rows: 152 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): region, pressure
dbl (3): 2020, 2050, 2100
ℹ 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.
pressure_data <- pressure_data %>%mutate(biome =case_when( region %in%c("ice", "tundra") ~"Tundra", region %in%c("wooded tundra", "Boreal forest", "Cool conifer forest") ~"Boreal",TRUE~NA_character_ ) ) %>%filter(!is.na(biome))pressure_long <- pressure_data %>%pivot_longer(cols =c(`2020`, `2050`, `2100`),names_to ="year",values_to ="msa_loss" )#loss in each yeardriver_summary_pct <- pressure_long %>%group_by(biome, year, pressure) %>%summarise(total_loss =sum(msa_loss, na.rm =TRUE), .groups ="drop") %>%group_by(biome, year) %>%mutate(percent = total_loss /sum(total_loss, na.rm =TRUE) *100,label =paste0(round(percent, 1), "%") ) %>%filter(percent >0.09)#top 8 affected ecoregions in each yeartop_pressures <- driver_summary_pct %>%group_by(biome, year) %>%slice_max(total_loss, n =8) %>%ungroup()#boreal forest top driversboreal_plot <-ggplot(top_pressures %>%filter(biome =="Boreal"),aes(x =reorder(pressure, total_loss), y = total_loss, fill = pressure)) +geom_col(show.legend =FALSE) +geom_text(aes(label = label), hjust =-0.1, size =3) +facet_wrap(~ year, scales ="free_y") +coord_flip() +labs(title ="Top Drivers of Biodiversity Loss in Boreal Biome",x ="Driver",y ="Total MSA Loss" ) +theme_minimal()# top drivers in tundratundra_plot <-ggplot(top_pressures %>%filter(biome =="Tundra"),aes(x =reorder(pressure, total_loss), y = total_loss, fill = pressure)) +geom_col(show.legend =FALSE) +geom_text(aes(label = label), hjust =-0.1, size =3) +facet_wrap(~ year, scales ="free_y") +coord_flip() +labs(title ="Top Drivers of Biodiversity Loss in Tundra Biome",x ="Driver",y ="Total MSA Loss" ) +theme_minimal()print(boreal_plot)
print(tundra_plot)
# Convert to long formatpressure_long <- pressure_data %>%pivot_longer(cols =c(`2020`, `2050`, `2100`),names_to ="year",values_to ="msa_loss") %>%filter(!is.na(biome)) %>%mutate(year =as.integer(year)) %>%# Drop very small or negative values (optional)filter(msa_loss >1e-4)# Aggregate total loss per pressure/biome/yearpressure_summary <- pressure_long %>%group_by(biome, year, pressure) %>%summarise(total_loss =sum(msa_loss, na.rm =TRUE), .groups ="drop")#stacked bar chart ggplot(pressure_summary, aes(x =factor(year), y = total_loss, fill = pressure)) +geom_bar(stat ="identity") +facet_wrap(~ biome, scales ="free_y") +labs(title ="Total MSA Loss by Pressure in Boreal and Tundra Biomes",x ="Year",y ="Total MSA Loss",fill ="Pressure" ) +theme_minimal(base_size =13) +theme(plot.title =element_text(face ="bold", size =16, hjust =0.5),legend.position ="bottom" )
to test if there’s a significant interaction between biome and pressure
A two-way ANOVA revealed significant main effects of biome (p = 0.005) and pressure (p < 0.001) on total biodiversity loss, as well as a significant interaction (p < 0.001) between biome and pressure. This indicates that the effect of a given pressure varies depending on the biome, confirming that certain drivers disproportionately impact either Boreal or Tundra ecosystems.
