C_03_Zero_trend_species.R

C - 03: Zero change species

Libraries

Code

library(dplyr)

Attache Paket: 'dplyr'

Die folgenden Objekte sind maskiert von 'package:stats':

    filter, lag

Die folgenden Objekte sind maskiert von 'package:base':

    intersect, setdiff, setequal, union

Code

library(here)

here() starts at C:/Users/wolke/OneDrive - CZU v Praze/Frieda_PhD_files/02_StaticPatterns/Git

Code

library(ggplot2)

Warning: Paket 'ggplot2' wurde unter R Version 4.4.3 erstellt

Code

library(ggthemes)
library(rstatix)

Attache Paket: 'rstatix'

Das folgende Objekt ist maskiert 'package:stats':

    filter

Small helper function

Code

get_mode_from_hist <- function(x, bins = 100) {
  h <- hist(x, breaks = bins, plot = TRUE)
  mode_bin_center <- h$mids[which.max(h$counts)]
  return(mode_bin_center)
}

Data for zero-trend species

log ratio >= log(0.9) & log(1.1) a) for sampling period 1 and 2 b) for sampling period 3

Code

# a) Sampling periods 1 and 2
zero1_2 <- 
  readRDS(here::here("Data/output/1_all_predictors_merged.rds")) %>%
  filter(log_R2_1 >= log(0.9) & log_R2_1 <= log(1.1))

# b) Sampling period 3
zero3 <- 
  readRDS(here::here("Data/output/2_big_table_3.rds")) %>%
  filter(log_R3_2 >= log(0.9) & log_R3_2 <= log(1.1)) %>%
  mutate(
    datasetID = as.factor(as.character(datasetID)),
    samplingPeriodID = as.factor(as.character(samplingPeriodID))
  )

Plots

Code

log_ratios <- 
  zero1_2 %>%
  mutate(log_ratio = log_R2_1) %>%
  full_join(zero3 %>% mutate(log_ratio = log_R3_2)) %>%
  select(-log_R2_1, -log_R3_2)

Joining with `by = join_by(datasetID, samplingPeriodID, verbatimIdentification,
Total_area_samp, Total_Ncells, Total_Ncells_samp, AOO, occ_Ncells,
rel_occ_Ncells, rel_AOO, Jaccard_dissim, a, b, c, d, D_AOO_a, log_ratio)`

Code

p <- 
  ggplot(log_ratios) +
  aes(x = log_ratio, 
      y = Jaccard_dissim, 
      fill = datasetID, 
      color = datasetID) +
  geom_smooth(method = "gam", se = TRUE, alpha = 0.2) +
  scale_fill_manual(values = c(
    "5" = "#E66101",
    "6" = "#FDB863",
    "13" = "#B2ABD2",
    "26" = "#5E3C99"
  )) +
  scale_color_manual(values = c(
    "5" = "#E66101",
    "6" = "#FDB863",
    "13" = "#B2ABD2",
    "26" = "#5E3C99"
  )) +
  labs(
    x = "log ratio AOO",
    y = "J",
    title = "Zero-change species",
    subtitle = "Jaccard~log ratio",
    fill = "Dataset",
    color = "Dataset"
  ) +
  ggthemes::theme_par() +
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(face = "italic"),
    axis.title.x = element_text(face = "italic")
  ) +
  facet_wrap(vars(samplingPeriodID))

p

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

Code

ggplot(log_ratios) +
  aes(x = log_ratio, 
      y = Jaccard_dissim, 
      fill = samplingPeriodID, 
      color = samplingPeriodID) +
  geom_smooth(method = "gam", se = TRUE, alpha = 0.2) +
  scale_fill_manual(values = c(
    "1" = "#E66101",
    "2" = "#FDB863",
    "3" = "#B2ABD2"
    )) +
  scale_color_manual(values = c(
    "1" = "#E66101",
    "2" = "#FDB863",
    "3" = "#B2ABD2"
  )) +
  labs(
    x = "log ratio AOO",
    y = "J",
    title = "Zero-change species",
    fill = "Sampling Period",
    color = "Sampling Period"
  ) +
  ggthemes::theme_par() +
  theme(
    legend.position = "bottom",
    axis.title.y = element_text(face = "italic"),
    axis.title.x = element_text(face = "italic")
  ) +
  facet_wrap(vars(datasetID),
             labeller = labeller(datasetID = c("5" = "Czechia",
                                               "6" = "New York",
                                               "13" = "Japan",
                                               "26" = "Europe")))

`geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

Summary stats

Code

log_ratios %>% 
  filter(log_ratio == 0) %>% 
  group_by(datasetID, samplingPeriodID) %>%
  summarize(
    n_sp = n_distinct(verbatimIdentification),
    n_sp_BL = n_distinct(scientificName))

`summarise()` has grouped output by 'datasetID'. You can override using the
`.groups` argument.

# A tibble: 10 × 4
# Groups:   datasetID [4]
   datasetID samplingPeriodID  n_sp n_sp_BL
   <fct>     <fct>            <int>   <int>
 1 5         1                    5       5
 2 5         2                    8       6
 3 5         3                    4       1
 4 6         1                    6       6
 5 6         2                    6       6
 6 13        1                    2       2
 7 13        2                    7       3
 8 13        3                    6       1
 9 26        1                    4       4
10 26        2                    4       4

Code

zero1_2 %>% 
  get_summary_stats(type = "common")

# A tibble: 26 × 10
   variable      n    min    max  median     iqr    mean      sd      se      ci
   <fct>     <dbl>  <dbl>  <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
 1 GlobRang…   800 8.18e2 1.18e8 1.17e+7 1.32e+7 1.48e+7 1.38e+7 4.89e+5 9.60e+5
 2 Mass        800 5.54e0 1.07e4 4.44e+1 2.07e+2 2.78e+2 8.45e+2 2.99e+1 5.87e+1
 3 pd          800 1.47e0 5.11e1 6.23e+0 5.05e+0 7.89e+0 5.50e+0 1.94e-1 3.82e-1
 4 Total_ar…   804 7.83e4 5.91e6 3.68e+5 5.83e+6 2.68e+6 2.86e+6 1.01e+5 1.98e+5
 5 Total_Nc…   804 6.71e2 5.34e3 5.08e+3 4.41e+3 3.58e+3 2.09e+3 7.36e+1 1.45e+2
 6 Total_Nc…   804 6.28e2 5.32e3 2.82e+3 2.19e+3 2.56e+3 1.64e+3 5.79e+1 1.14e+2
 7 AOO         804 4.98e0 5.77e6 1.14e+5 1.83e+6 1.22e+6 1.82e+6 6.44e+4 1.26e+5
 8 occ_Ncel…   804 1   e0 5.23e3 6.25e+2 1.59e+3 1.27e+3 1.35e+3 4.75e+1 9.33e+1
 9 rel_occ_…   804 0      1   e0 6.03e-1 6.73e-1 5.37e-1 3.45e-1 1.2 e-2 2.4 e-2
10 rel_AOO     804 0      1   e0 6.41e-1 6.87e-1 5.54e-1 3.52e-1 1.2 e-2 2.4 e-2
# ℹ 16 more rows

Code

mode_estimate <- get_mode_from_hist(log_ratios$log_ratio)

Code

mode_estimate

[1] -0.001