C_01_Validation.R

C - 01: Validation

Libraries

Code

suppressPackageStartupMessages({
  
library(here)

source(here::here("Code/00_Configuration.R"))

package_list <- 
  c(package_list, 
    "ranger", 
    "tidymodels", 
    "caret", 
    "skimr", 
    "DALEXtra",
    "DALEX",
    "ggthemes",
    "yardstick")

lapply(package_list, require, character = TRUE)

tidymodels_prefer()

})


rm(list = ls())

Read in data

Code

# Models to fit on tp 2 data
## Pooled models
res_full <-
  readRDS(here::here(
    "Data/output/B_models/",
    "B_01_list_all_results_rf.rds"
  ))

## Split models
res_split <-
  readRDS(here::here(
    "Data/output/B_models/",
    "B_02_rf_res_atlas_split.rds"
  ))

# Response to test predictions against from tp 3 data
tp3 <-
  readRDS(here::here(
    "Data/output/",
    "2_big_table_3.rds"
    )) %>%
  select(
    datasetID,
    verbatimIdentification,
    Jaccard_dissim,
    log_R3_2
  ) %>%
  mutate(
    datasetID = as.factor(datasetID)
  )

# Tp 2 data = predictors for model predictions for change towards tp3
dta2 <-
  readRDS(here::here(
    "Data/output/",
    "1_all_predictors_merged.rds"
  )) %>%
  filter(samplingPeriodID == 2) %>%
  filter(datasetID %in% c(5, 13))

Code

# Jaccard
Jaccard_validation <- 
  res_full$res[[1]] %>%
  extract_fit_parsnip() %>%
  predict(new_data = dta2) %>%
  bind_cols(dta2 %>% 
              select(verbatimIdentification, datasetID)) %>%
  left_join(tp3) %>%
  mutate(diff = Jaccard_dissim - .pred) %>%
  select(-log_R3_2)

glimpse(Jaccard_validation)

Rows: 787
Columns: 5
$ .pred                  <dbl> 0.433738362, 0.433738362, 0.065348790, 0.065348…
$ verbatimIdentification <chr> "Nucifraga caryocatactes", "Nucifraga caryocata…
$ datasetID              <fct> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
$ Jaccard_dissim         <dbl> 0.419, 0.419, 0.027, 0.027, 0.317, 0.317, 0.250…
$ diff                   <dbl> -1.473836e-02, -1.473836e-02, -3.834879e-02, -3…

Code

# Log ratio
logRatio_validation <- 
  res_full$res[[2]] %>%
  extract_fit_parsnip() %>%
  predict(new_data = dta2) %>%
  bind_cols(dta2 %>% 
              select(verbatimIdentification, datasetID)) %>%
  left_join(tp3) %>%
  mutate(diff = log_R3_2 - .pred) %>%
  select(-Jaccard_dissim)

glimpse(logRatio_validation)

Rows: 787
Columns: 5
$ .pred                  <dbl> 0.030536328, 0.030536328, 0.068346074, 0.068346…
$ verbatimIdentification <chr> "Nucifraga caryocatactes", "Nucifraga caryocata…
$ datasetID              <fct> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
$ log_R3_2               <dbl> -0.225, -0.225, 0.008, 0.008, -0.177, -0.177, -…
$ diff                   <dbl> -0.255536328, -0.255536328, -0.060346074, -0.06…

Plots

Code

hist(logRatio_validation$diff)

Code

hist(Jaccard_validation$diff)

Code

Jaccard_validation %>%
  ggplot() +
  geom_histogram(aes(x = diff, fill = datasetID), 
                 alpha = 0.5, position = "identity") +
  ggthemes::theme_par()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 31 rows containing non-finite outside the scale range
(`stat_bin()`).

Code

logRatio_validation %>%
  ggplot() +
  geom_histogram(aes(x = diff, fill = datasetID), 
                 alpha = 0.5, position = "identity") +
  ggthemes::theme_par()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 31 rows containing non-finite outside the scale range
(`stat_bin()`).

Code

p_1 <- Jaccard_validation %>%
  ggplot(aes(x = .pred, y = Jaccard_dissim, col = datasetID)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_point(alpha = 0.5) +
  geom_smooth(aes(col = datasetID), method = "lm") +
  ggthemes::theme_few() +
  scale_color_manual(values = c("#f1a340", "#998ec3"))
p_1

Code

p_2 <- logRatio_validation %>%
  ggplot(aes(x = .pred, y = log_R3_2, col = datasetID)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_point(alpha = 0.5) +
  xlim(-4, 2) +
  # geom_smooth(method = "lm")+
  ggthemes::theme_few() +
  scale_color_manual(values = c("#f1a340", "#998ec3"))


p_2

Code

ggsave("Figures/C_validation/C_01_Validation_Jaccard.svg", p_1, width = 5, height = 4)
ggsave("Figures/C_validation/C_01_Validation_logRatio.svg", p_2, width = 5, height = 4)

esquisse::ggplot_to_ppt() # saved to here("Figures/C_validation/Prediction_against_replication_3_plots.pptx")

Stats

Jaccard dissimilarity model

Code

# Split model
metrics_jaccard <- 
  Jaccard_validation %>%
  group_by(datasetID) %>%
  metrics(truth = Jaccard_dissim, 
          estimate = .pred)

Code

# Pooled model
metrics_jaccard_pooled <- 
  Jaccard_validation %>%
  metrics(truth = Jaccard_dissim, 
          estimate = .pred)

Pearson correlation for Jaccard dissimilarity

Code

# correlation between predicted and observed (split model)
cor_jaccard <- 
  cor.test(Jaccard_validation$.pred, 
                        Jaccard_validation$Jaccard_dissim, 
                        method = "pearson")

log ratio model

Code

metrics_logRatio <- 
  logRatio_validation %>%
  group_by(datasetID) %>%
  metrics(truth = log_R3_2, 
          estimate = .pred)

metrics_logRatio_pooled <- 
  logRatio_validation %>%
  metrics(truth = log_R3_2, 
          estimate = .pred)

print(metrics_jaccard)

# A tibble: 6 × 4
  datasetID .metric .estimator .estimate
  <fct>     <chr>   <chr>          <dbl>
1 5         rmse    standard      0.0885
2 13        rmse    standard      0.121 
3 5         rsq     standard      0.916 
4 13        rsq     standard      0.779 
5 5         mae     standard      0.0549
6 13        mae     standard      0.0774

Code

print(metrics_logRatio)

# A tibble: 6 × 4
  datasetID .metric .estimator .estimate
  <fct>     <chr>   <chr>          <dbl>
1 5         rmse    standard      0.428 
2 13        rmse    standard      0.638 
3 5         rsq     standard      0.167 
4 13        rsq     standard      0.0102
5 5         mae     standard      0.245 
6 13        mae     standard      0.390 

Pearson correlation for log ratio

Code

# Spearman
cor_logRatio <- 
  cor.test(logRatio_validation$.pred, 
           logRatio_validation$log_R3_2, 
           method = "spearman")

Warning in cor.test.default(logRatio_validation$.pred,
logRatio_validation$log_R3_2, : Kann exakten p-Wert bei Bindungen nicht
berechnen

Code

# Pearson
cor.test(logRatio_validation$.pred, 
         logRatio_validation$log_R3_2, 
         method = "pearson")

    Pearson's product-moment correlation

data:  logRatio_validation$.pred and logRatio_validation$log_R3_2
t = 6.2384, df = 754, p-value = 7.366e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1526506 0.2882928
sample estimates:
      cor 
0.2215431 

Code

print(cor_jaccard)

    Pearson's product-moment correlation

data:  Jaccard_validation$.pred and Jaccard_validation$Jaccard_dissim
t = 74.154, df = 754, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9285569 0.9458305
sample estimates:
      cor 
0.9377712 

Code

print(cor_logRatio)

    Spearman's rank correlation rho

data:  logRatio_validation$.pred and logRatio_validation$log_R3_2
S = 58486931, p-value = 1.968e-07
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1878328 

correctly predicted species

Code

logRatio_validation %>%
  filter(abs(diff) < 0.01) %>%
  distinct() %>%
  skimr::skim()

Data summary

Name	Piped data
Number of rows	30
Number of columns	5
_______________________
Column type frequency:
character	1
factor	1
numeric	3
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
verbatimIdentification	0	1	11	23	0	30	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
datasetID	0	1	FALSE	2	5: 23, 13: 7, 6: 0, 26: 0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
.pred	0	1	0.02	0.04	-0.03	0.00	0.01	0.02	0.16	▇▇▁▁▂
log_R3_2	0	1	0.02	0.05	-0.02	0.00	0.00	0.02	0.16	▇▂▁▁▁
diff	0	1	0.00	0.01	-0.01	-0.01	0.00	0.00	0.01	▇▆▃▇▅

jaccard histograms for each region

Code

Jaccard_validation %>%
  filter(Jaccard_dissim > 0.75) %>%
  distinct() %>%
  arrange(desc(Jaccard_dissim)) %>% 
  group_by(datasetID) %>%
  ggplot() +
  geom_histogram(aes(Jaccard_dissim)) +
  facet_wrap(~datasetID)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Inspect special-case species

Code

dta <- readRDS(here("Data/output/1_all_predictors_merged.rds"))

High Jaccard dissim species:

Code

dta %>%
  filter(Jaccard_dissim > 0.9) %>%
  filter(samplingPeriodID == 1) %>%
  select(Jaccard_dissim, log_R2_1, a, b, c, d, rel_AOO, datasetID, verbatimIdentification) %>%
  distinct() %>%
  arrange(desc(Jaccard_dissim)) %>%
  group_by(datasetID) %>%
  skimr::skim()

Data summary

Name	Piped data
Number of rows	113
Number of columns	9
_______________________
Column type frequency:
character	1
numeric	7
________________________
Group variables	datasetID

Variable type: character

skim_variable	datasetID	n_missing	complete_rate	min	max	empty	n_unique	whitespace
verbatimIdentification	5	0	1	9	25	0	20	0
verbatimIdentification	6	0	1	9	26	0	42	0
verbatimIdentification	13	0	1	9	25	0	42	0
verbatimIdentification	26	0	1	10	18	0	9	0

Variable type: numeric

skim_variable	datasetID	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Jaccard_dissim	5	0	1	0.95	0.04	0.90	0.91	0.95	1.00	1.00	▆▃▃▁▇
Jaccard_dissim	6	0	1	0.95	0.04	0.90	0.92	0.94	1.00	1.00	▇▃▅▁▇
Jaccard_dissim	13	0	1	0.95	0.03	0.90	0.93	0.94	1.00	1.00	▅▇▆▁▇
Jaccard_dissim	26	0	1	0.95	0.03	0.92	0.93	0.95	0.98	1.00	▇▁▃▂▃
log_R2_1	5	0	1	0.46	0.80	-0.80	-0.13	0.33	0.92	2.37	▅▇▆▃▁
log_R2_1	6	0	1	0.03	1.33	-1.79	-1.02	-0.24	0.70	3.75	▇▇▃▃▁
log_R2_1	13	0	1	0.86	1.72	-1.86	-0.32	0.87	1.46	7.96	▃▇▂▁▁
log_R2_1	26	0	1	-0.23	1.96	-3.17	-1.39	-0.24	0.31	3.11	▃▂▇▁▃
a	5	0	1	2.20	3.27	0.00	0.00	1.00	3.00	13.00	▇▂▁▁▁
a	6	0	1	8.98	16.50	0.00	0.00	2.50	11.25	71.00	▇▁▁▁▁
a	13	0	1	4.60	5.78	0.00	0.00	2.50	6.50	25.00	▇▁▂▁▁
a	26	0	1	3.33	3.94	0.00	1.00	1.00	3.00	11.00	▇▃▁▁▃
b	5	0	1	19.55	20.36	4.00	9.50	12.50	20.75	93.00	▇▂▁▁▁
b	6	0	1	49.98	80.98	1.00	6.50	16.00	62.75	428.00	▇▁▁▁▁
b	13	0	1	42.57	46.91	1.00	9.25	28.50	53.25	214.00	▇▂▁▁▁
b	26	0	1	30.11	26.86	0.00	2.00	25.00	56.00	67.00	▇▅▁▂▇
c	5	0	1	11.85	10.85	0.00	4.75	7.00	18.00	41.00	▇▂▃▁▂
c	6	0	1	71.50	132.38	1.00	5.25	20.00	47.25	620.00	▇▁▁▁▁
c	13	0	1	23.38	29.41	1.00	4.00	10.00	35.75	142.00	▇▃▁▁▁
c	26	0	1	25.22	23.95	0.00	11.00	22.00	25.00	81.00	▆▇▂▁▂
d	5	0	1	594.40	31.23	489.00	587.75	604.00	614.00	622.00	▁▁▁▂▇
d	6	0	1	5188.55	191.11	4532.00	5179.75	5275.00	5297.50	5317.00	▁▁▁▁▇
d	13	0	1	1113.45	67.03	905.00	1067.00	1138.00	1161.25	1181.00	▁▁▂▂▇
d	26	0	1	2762.33	41.79	2685.00	2759.00	2771.00	2795.00	2807.00	▃▁▁▇▆
rel_AOO	5	0	1	0.02	0.02	0.00	0.01	0.01	0.03	0.08	▇▂▂▁▁
rel_AOO	6	0	1	0.02	0.03	0.00	0.00	0.00	0.01	0.13	▇▁▁▁▁
rel_AOO	13	0	1	0.02	0.03	0.00	0.00	0.01	0.04	0.14	▇▂▁▁▁
rel_AOO	26	0	1	0.01	0.01	0.00	0.01	0.01	0.01	0.04	▇▆▂▁▂

High relative AOO but not zero-change species:

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  group_by(datasetID) %>%
  filter(rel_AOO >= 0.8) %>%
  select(Jaccard_dissim, log_R2_1, a, b, c, d, rel_AOO, datasetID, verbatimIdentification) %>%
  filter(log_R2_1 != 0) %>%
  distinct() %>%
  arrange(desc(Jaccard_dissim)) %>%
  group_by(datasetID) %>%
  skimr::skim()

Data summary

Name	Piped data
Number of rows	136
Number of columns	9
_______________________
Column type frequency:
character	1
numeric	7
________________________
Group variables	datasetID

Variable type: character

skim_variable	datasetID	n_missing	complete_rate	min	max	empty	n_unique	whitespace
verbatimIdentification	5	0	1	9	29	0	64	0
verbatimIdentification	6	0	1	14	23	0	32	0
verbatimIdentification	13	0	1	11	23	0	6	0
verbatimIdentification	26	0	1	9	23	0	34	0

Variable type: numeric

skim_variable	datasetID	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Jaccard_dissim	5	0	1	0.08	0.07	0.00	0.01	0.07	0.14	0.28	▇▃▃▂▁
Jaccard_dissim	6	0	1	0.14	0.07	0.03	0.08	0.12	0.20	0.26	▆▇▅▅▅
Jaccard_dissim	13	0	1	0.14	0.06	0.09	0.10	0.13	0.16	0.24	▇▁▅▁▂
Jaccard_dissim	26	0	1	0.10	0.04	0.04	0.07	0.10	0.11	0.18	▆▅▇▁▂
log_R2_1	5	0	1	0.01	0.03	-0.10	0.00	0.01	0.02	0.11	▁▁▇▂▁
log_R2_1	6	0	1	-0.01	0.04	-0.08	-0.04	0.00	0.01	0.09	▃▅▇▃▁
log_R2_1	13	0	1	0.03	0.03	-0.01	0.01	0.04	0.05	0.08	▇▁▇▃▃
log_R2_1	26	0	1	0.03	0.03	-0.04	0.02	0.03	0.04	0.09	▂▂▇▃▂
a	5	0	1	566.81	56.01	415.00	525.00	576.50	620.25	627.00	▁▁▃▃▇
a	6	0	1	4337.84	481.54	3588.00	3975.25	4307.50	4729.25	5138.00	▇▇▆▇▇
a	13	0	1	951.00	81.48	804.00	939.00	968.00	984.25	1047.00	▂▁▂▇▂
a	26	0	1	2281.94	182.59	1931.00	2150.25	2282.50	2450.75	2551.00	▃▆▆▃▇
b	5	0	1	28.72	25.42	0.00	6.00	24.00	46.50	103.00	▇▃▃▂▁
b	6	0	1	330.47	189.52	80.00	190.00	278.00	490.00	700.00	▇▇▂▅▃
b	13	0	1	100.83	31.31	44.00	92.75	111.50	118.25	131.00	▂▁▂▂▇
b	26	0	1	162.88	57.94	61.00	105.75	165.50	201.00	310.00	▆▃▇▃▁
c	5	0	1	21.09	21.32	0.00	3.00	17.50	31.50	83.00	▇▅▂▁▁
c	6	0	1	361.91	204.87	49.00	198.00	318.50	569.25	716.00	▇▇▃▆▆
c	13	0	1	57.00	42.10	14.00	26.75	57.50	64.25	130.00	▇▇▃▁▃
c	26	0	1	77.82	54.19	19.00	40.25	58.50	100.50	244.00	▇▃▂▁▁
d	5	0	1	11.38	16.72	0.00	0.00	6.50	12.25	83.00	▇▁▁▁▁
d	6	0	1	288.78	183.16	41.00	134.50	298.00	426.00	693.00	▇▂▅▃▂
d	13	0	1	75.17	34.00	32.00	50.50	76.50	98.00	119.00	▇▁▇▁▇
d	26	0	1	298.35	121.16	127.00	187.50	300.50	396.50	640.00	▇▆▆▃▁
rel_AOO	5	0	1	0.94	0.06	0.80	0.90	0.96	0.99	1.00	▁▁▃▂▇
rel_AOO	6	0	1	0.89	0.06	0.81	0.83	0.90	0.94	0.98	▇▃▃▆▆
rel_AOO	13	0	1	0.89	0.04	0.81	0.87	0.90	0.91	0.94	▃▃▁▇▇
rel_AOO	26	0	1	0.88	0.05	0.80	0.84	0.89	0.93	0.95	▆▃▅▃▇

saturated occupancy species

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  group_by(datasetID) %>%
  select(verbatimIdentification, datasetID, rel_AOO) %>%
  slice_max(rel_AOO)

# A tibble: 9 × 3
# Groups:   datasetID [4]
  verbatimIdentification datasetID rel_AOO
  <chr>                  <fct>       <dbl>
1 Turdus merula          5           1    
2 Parus major            5           1    
3 Emberiza citrinella    5           1    
4 Fringilla coelebs      5           1    
5 Hirundo rustica        5           1    
6 Phoenicurus ochruros   5           1    
7 Turdus migratorius     6           0.978
8 Cettia diphone         13          0.939
9 Motacilla alba         26          0.952

Special Case: Czechia (= small country)

Code

dta %>%
  filter(datasetID == 5) %>%
  slice_max(abs(log_R2_1))

  datasetID samplingPeriodID verbatimIdentification   scientificName
1         5                2       Egretta garzetta Egretta garzetta
2         5                1       Egretta garzetta Egretta garzetta
  GlobRangeSize_km2 Mass Migration       pd Total_area_samp Total_Ncells
1          34302524  312         2 6.799823        78308.81          671
2          34302524  312         2 6.799823        78308.81          671
  Total_Ncells_samp      AOO occ_Ncells rel_occ_Ncells rel_AOO Jaccard_dissim a
1               628 1449.940         11          0.018   0.019          0.909 1
2               628  136.033          1          0.002   0.002          0.909 1
   b c   d D_AOO_a time_span startYear1 endYear2 n_years log_R2_1
1 10 0 617   0.763         2       1985     2003      19    2.366
2 10 0 617   0.000         4       1985     2003      19    2.366
  log_R2_1_per_year joincount_delta  circNorm minDist_toBorder_centr mean_lnLac
1             0.125       0.4862301 10.511814              19280.184   1.464422
2             0.125       0.0000000  1.276133               8268.471   2.591765
   Habitat_5 Generalism Threatened
1 freshwater          0          0
2 freshwater          0          0

Zero-change but Jaccard != 0

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  group_by(datasetID) %>%
  mutate(abs_logR = abs(log_R2_1)) %>%
  filter(abs_logR <= 0.1) %>%
  select(verbatimIdentification, Jaccard_dissim, abs_logR, log_R2_1, datasetID, rel_AOO) %>% # %>% skimr::skim()
  filter(Jaccard_dissim != 0)

# A tibble: 394 × 6
# Groups:   datasetID [4]
   verbatimIdentification Jaccard_dissim abs_logR log_R2_1 datasetID rel_AOO
   <chr>                           <dbl>    <dbl>    <dbl> <fct>       <dbl>
 1 Carduelis carduelis             0.024    0.009    0.009 5           0.987
 2 Passer montanus                 0.063    0.017   -0.017 5           0.971
 3 Sylvia borin                    0.14     0.018    0.018 5           0.916
 4 Pyrrhula pyrrhula               0.117    0.001   -0.001 5           0.897
 5 Parus palustris                 0.18     0.043    0.043 5           0.891
 6 Podiceps cristatus              0.315    0.042   -0.042 5           0.617
 7 Phylloscopus collybita          0.011    0.001    0.001 5           0.994
 8 Cuculus canorus                 0.021    0.006    0.006 5           0.988
 9 Columba palumbus                0.019    0.01     0.01  5           0.987
10 Parus ater                      0.083    0.017    0.017 5           0.933
# ℹ 384 more rows

Zero Change ~ Jaccard boxplots

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  group_by(datasetID) %>%
  mutate(abs_logR = abs(log_R2_1)) %>%
  filter(abs_logR <= 0.1) %>%
  select(verbatimIdentification, Jaccard_dissim, abs_logR, log_R2_1, datasetID, rel_AOO) %>%
  ggplot(aes(y = Jaccard_dissim, x = datasetID)) +
  geom_boxplot()

Relationship absolute log ratio and Jaccard

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  mutate(
    abs_logR = abs(log_R2_1),
    abs_logR_y = abs(log_R2_1_per_year)
  ) %>%
  ggplot(aes(color = datasetID, x = abs_logR_y, y = Jaccard_dissim)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "gam") +
  xlim(0, 0.3) +
  geom_vline(xintercept = 0) +
  ylim(0, 1) +
  facet_wrap(~datasetID, scales = "free_x") +
  ggthemes::theme_par()

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

Absolute log ratio per year

Code

dta %>%
  filter(samplingPeriodID == 1) %>%
  group_by(datasetID) %>%
  mutate(
    abs_logR = abs(log_R2_1),
    abs_logR_y = abs(log_R2_1_per_year)
  ) %>%
  filter(abs_logR_y <= 0.001) %>%
  select(verbatimIdentification, abs_logR_y, Jaccard_dissim, abs_logR, log_R2_1, datasetID, rel_AOO) %>%
  filter(Jaccard_dissim <= 0.1)

# A tibble: 78 × 7
# Groups:   datasetID [4]
   verbatimIdentification abs_logR_y Jaccard_dissim abs_logR log_R2_1 datasetID
   <chr>                       <dbl>          <dbl>    <dbl>    <dbl> <fct>    
 1 Carduelis carduelis         0              0.024    0.009    0.009 5        
 2 Passer montanus             0.001          0.063    0.017   -0.017 5        
 3 Phylloscopus collybita      0              0.011    0.001    0.001 5        
 4 Cuculus canorus             0              0.021    0.006    0.006 5        
 5 Columba palumbus            0.001          0.019    0.01     0.01  5        
 6 Parus ater                  0.001          0.083    0.017    0.017 5        
 7 Sturnus vulgaris            0              0.01     0.002    0.002 5        
 8 Turdus merula               0              0.002    0.002   -0.002 5        
 9 Sitta europaea              0              0.011    0.005    0.005 5        
10 Prunella modularis          0.001          0.059    0.011    0.011 5        
# ℹ 68 more rows
# ℹ 1 more variable: rel_AOO <dbl>