Analysis of population and range change in European mammals from Article 17 data

Explore population trends and distribution changes of the Wildlife Comeback in Europe Report species that have records within the Article 17 data.

mammals_list = c(
"Megaptera_novaeangliae",
#"Caretta_caretta", ##Loggerheads
"Pusa_hispida",
"Lutra_lutra",
"Meles_meles",
"Martes_martes",
"Myotis_daubentonii",
"Myotis_nattereri",
"Myotis_emarginatus",
"Alces_alces",
"Bison_bonasus",
"Canis_aureus",
"Canis_lupus",
"Capra_ibex",
"Capra_pyrenaica",
"Capreolus_capreolus",
"Castor_fiber",
"Cervus_elaphus",
"Gulo_gulo",
"Halichoerus_grypus",
"Lynx_lynx",
"Lynx_pardinus",
"Phoca_vitulina",
"Rupicapra_pyrenaica",
"Rupicapra_rupicapra",
"Sus_scrofa",
"Ursus_arctos")

mammals_list = gsub("_", " ", mammals_list)

List of mammals for analysis

mammals_list
##  [1] "Megaptera novaeangliae" "Pusa hispida"           "Lutra lutra"           
##  [4] "Meles meles"            "Martes martes"          "Myotis daubentonii"    
##  [7] "Myotis nattereri"       "Myotis emarginatus"     "Alces alces"           
## [10] "Bison bonasus"          "Canis aureus"           "Canis lupus"           
## [13] "Capra ibex"             "Capra pyrenaica"        "Capreolus capreolus"   
## [16] "Castor fiber"           "Cervus elaphus"         "Gulo gulo"             
## [19] "Halichoerus grypus"     "Lynx lynx"              "Lynx pardinus"         
## [22] "Phoca vitulina"         "Rupicapra pyrenaica"    "Rupicapra rupicapra"   
## [25] "Sus scrofa"             "Ursus arctos"

How many do we match in article 17 data?

#### *** 2020 species data
#mammal_species <- fread("Article17_2020_dataset/Article17_2020_speciesEUassessment.csv", na.strings = 'NULL')
mammal_species <- fread("Article17_2020_dataset/Article17_2020_data_species_regions.csv", na.strings = 'NULL')
#write.csv(bird_species,"bird_species_check.csv")

#intersect(mammal_species$assessment_speciesname, mammals_list)
intersect(mammal_species$speciesname, mammals_list)
##  [1] "Myotis daubentonii"     "Myotis emarginatus"     "Myotis nattereri"      
##  [4] "Castor fiber"           "Canis lupus"            "Canis aureus"          
##  [7] "Ursus arctos"           "Lutra lutra"            "Martes martes"         
## [10] "Lynx lynx"              "Rupicapra rupicapra"    "Capra ibex"            
## [13] "Halichoerus grypus"     "Phoca vitulina"         "Bison bonasus"         
## [16] "Megaptera novaeangliae" "Lynx pardinus"          "Capra pyrenaica"       
## [19] "Rupicapra pyrenaica"    "Gulo gulo"              "Pusa hispida"
setdiff(mammals_list, mammal_species$speciesname)
## [1] "Meles meles"         "Alces alces"         "Capreolus capreolus"
## [4] "Cervus elaphus"      "Sus scrofa"
## 21 of 26 match****
#Meles meles
#Sus scrofa
#Alces alces
#Capreolus capreolus
#Cervus elaphus

## I think these are missing as they're not Article 17 species (check?)

# Restrict to those in list
# **
# ** Need to check matches as losing some
# **
mammal_species_short <- mammal_species[mammal_species$speciesname %in% mammals_list, ]

#names(bird_species_short)[names(bird_species_short) == "ï..specieshash"] <- "specieshash"

# Get the species codes (unique for each species)
speciescode <- unique(mammal_species$speciescode)

##Threat information

Extract and count the number of threats for each species in each country

# Threats info

#### Threats
mammal_threats <- fread("Article17_2020_dataset/Article17_2020_data_pressures_threats.csv")
mammal_threats = subset(mammal_threats, featuretype == "species")
mammal_threats$featurecode = as.numeric(mammal_threats$featurecode) # Some habitat codes become NA if we do this first

# Subset threats for species we care about
mammal_threats = mammal_threats[mammal_threats$featurecode %in% speciescode, ]


# Filter out 'no-pressure' codes 
mammal_threats = subset(mammal_threats, pressurecode != "Xxt") # Get rid of Xxt

# Get rid of all of those starting with 'X'??
mammal_threats = subset(mammal_threats, !startsWith(mammal_threats$pressurecode, "X"))
mammal_threats = subset(mammal_threats, !is.na(mammal_threats$pressurecode))

mammal_threats_both = mammal_threats %>% 
  group_by(country, region, featurename, featurecode, type) %>% 
  summarise(
    n = n_distinct(pressurecode),
    n_hl = n_distinct(substring(pressurecode, 1, 1)) #n unique codes using only first character
  )
## `summarise()` has grouped output by 'country', 'region', 'featurename', 'featurecode'. You can override using the `.groups` argument.
# Convert this new data to four columns (n_t, n_p, n_hl_t, n_hl_p)
mammal_threats_n = tidyr::pivot_wider(mammal_threats_both, 
                                      id_cols = country:featurecode, 
                                      names_from="type", 
                                      values_from = c("n", "n_hl"),
                                      values_fill = 0) # Fill missing counts with 0

# Give them better names
mammal_threats_n = rename(mammal_threats_n, 
                        n_threat = n_t, 
                        n_pressure = n_p,
                        n_hl_threat = n_hl_t,
                        n_hl_pressure = n_hl_p)
# Measures
mammal_measures <- fread("Article17_2020_dataset/Article17_2020_data_measures.csv", na.strings = 'NULL')
mammal_measures$featurecode = as.numeric(mammal_measures$featurecode) # Some habitat codes become NA
## Warning: NAs introduced by coercion
# Subset measures for species we care about
mammal_measures <- mammal_measures[mammal_measures$featurecode %in% speciescode, ]
mammal_measures = subset(mammal_measures, featuretype == "species")

mammal_measures_n = mammal_measures %>% 
  group_by(country, region, featurename, featurecode) %>% 
  summarise(
    n_measure = n_distinct(measurecode),
    n_hl_measure = n_distinct(substring(measurecode, 1, 2)) #n unique codes using only first character
  ) 
## `summarise()` has grouped output by 'country', 'region', 'featurename'. You can override using the `.groups` argument.

##Combine data

# Trend data - get columns from bird species table that we want to use
mammal_trend_data = mammal_species_short %>% select(country, region, speciescode, speciesname,
                                                  population_trend,
                                                  population_trend_long,
                                                  range_trend,
                                                  range_trend_long,
                                                  habitat_trend,
                                                  habitat_trend_long,
                                                  population_trend_long_magnitude_min, population_trend_long_magnitude_max,
                                                  population_trend_magnitude_min, population_trend_magnitude_max)

mammal_trend_data = rename(mammal_trend_data, 
                        featurecode = speciescode, 
                        featurename = speciesname)

mammal_trend_data$population_trend_long_magnitude_mid = (mammal_trend_data$population_trend_long_magnitude_min + mammal_trend_data$population_trend_long_magnitude_max)/2.0
mammal_trend_data$population_trend_magnitude_mid = (mammal_trend_data$population_trend_magnitude_min + mammal_trend_data$population_trend_magnitude_max)/2.0

# All together now!
# Combined measure data with threat data
full_dataB = merge(mammal_trend_data, mammal_threats_n, by = c('region', 'country', 'featurename', 'featurecode'), all.x =TRUE)
# Combine those with trend data
full_data = merge(full_dataB, mammal_measures_n, by=c('region', 'country', 'featurename', 'featurecode'), all.x =TRUE)


model_data = full_data

Rate of change vs Threats/Measures

Are species trends different where there are more threats/measures?

# Model data is only declines or increases (remapped to 'increasing' as 0/1)
model_data_trend = model_data %>%
  ungroup() %>%
  mutate(trend = population_trend_magnitude_mid) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(trend, n_hl_threat, n_hl_pressure, n_hl_measure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns

model_data_trend$trend = log10(model_data_trend$trend + 1)

# Are populations more likely to be declining where there are more threats/measures
m1_NULL = lmer(trend ~ (1|featurename), 
            data = model_data_trend) 
# Are populations more likely to be declining where there are more threats/measures
m1_threat = lmer(trend ~ scale(n_hl_threat) + (1|featurename), 
            data = model_data_trend) 
m1_measure = lmer(trend ~ scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m1_pressuee = lmer(trend ~scale(n_hl_pressure) + (1|featurename), 
            data = model_data_trend) 
m1A = lmer(trend ~  scale(n_hl_threat) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m1B = lmer(trend ~  scale(n_hl_threat) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m1C = lmer(trend ~  scale(n_hl_pressure) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m1D = lmer(trend ~  scale(n_hl_pressure) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
# Adding 'country' throws singularities due to sparse data
#summary(m1)

tab_model(m1_NULL, m1_threat, m1_measure, m1A, m1B, m1C, m1D, show.aicc = TRUE)
  trend trend trend trend trend trend trend
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 1.67 1.28 – 2.06 <0.001 1.68 1.28 – 2.07 <0.001 1.67 1.28 – 2.07 <0.001 1.67 1.28 – 2.07 <0.001 1.81 1.37 – 2.24 <0.001 1.69 1.29 – 2.09 <0.001 1.73 1.30 – 2.16 <0.001
n_hl_threat 0.03 -0.24 – 0.30 0.838 0.03 -0.33 – 0.40 0.867 -0.04 -0.41 – 0.34 0.850
n_hl_measure 0.02 -0.25 – 0.29 0.906 -0.00 -0.37 – 0.36 0.981 0.09 -0.30 – 0.47 0.663 -0.15 -0.44 – 0.15 0.327 -0.11 -0.44 – 0.22 0.508
n_hl_threat *
n_hl_measure 		-0.20 -0.49 – 0.09 0.173
n_hl_pressure 0.34 0.04 – 0.64 0.027 0.31 -0.02 – 0.63 0.065
n_hl_pressure *
n_hl_measure 		-0.08 -0.35 – 0.19 0.577
Random Effects
σ2 0.44 0.46 0.46 0.47 0.46 0.39 0.40
τ00 0.26 featurename 0.26 featurename 0.26 featurename 0.26 featurename 0.25 featurename 0.29 featurename 0.29 featurename
ICC 0.37 0.36 0.37 0.36 0.35 0.43 0.42
N 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename
Observations 36 36 36 36 36 36 36
Marginal R2 / Conditional R2 0.000 / 0.367 0.001 / 0.365 0.000 / 0.366 0.001 / 0.359 0.049 / 0.378 0.100 / 0.484 0.101 / 0.478
AIC 90.846 95.486 95.510 99.730 102.790 95.571 100.301 
compare_performance(m1_NULL, m1_threat, m1_measure, m1A, m1B, m1C, m1D)
## # Comparison of Model Performance Indices
## 
## Model      |    Type |   AIC |    BIC | R2_conditional | R2_marginal |  ICC | RMSE |       BF
## ---------------------------------------------------------------------------------------------
## m1_NULL    | lmerMod | 90.10 |  94.85 |           0.37 |        0.00 | 0.37 | 0.60 |         
## m1_threat  | lmerMod | 94.20 | 100.53 |           0.36 |    1.00e-03 | 0.36 | 0.60 |     0.06
## m1_measure | lmerMod | 94.22 | 100.55 |           0.37 |        0.00 | 0.37 | 0.60 |     0.06
## m1A        | lmerMod | 97.73 | 105.65 |           0.36 |    1.00e-03 | 0.36 | 0.60 | 5.00e-03
## m1B        | lmerMod | 99.89 | 109.39 |           0.38 |        0.05 | 0.35 | 0.59 | 1.00e-03
## m1C        | lmerMod | 93.57 | 101.49 |           0.48 |        0.10 | 0.43 | 0.54 |     0.04
## m1D        | lmerMod | 97.41 | 106.91 |           0.48 |        0.10 | 0.42 | 0.54 | 2.00e-03
No clear support for threats/measures helping predict change - null model preferred 
plot_model(m1_NULL, show.intercept = TRUE)

tab_model(m1_NULL)
  trend
Predictors Estimates CI p
(Intercept) 1.67 1.28 – 2.06 <0.001
Random Effects
σ2 0.44
τ00 featurename 0.26
ICC 0.37
N featurename 10
Observations 36
Marginal R2 / Conditional R2 0.000 / 0.367 
plot_model(m1_NULL, type = "pred")
## list()

Long-term change

Are species long-term trends different where there are more threats/measures?

# Model data is only declines or increases (remapped to 'increasing' as 0/1)
model_data_trend = model_data %>%
  ungroup() %>%
  mutate(trend = log10(population_trend_long_magnitude_mid)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(trend, n_hl_threat, n_hl_pressure, n_hl_measure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns

model_data_trend$trend = log10(model_data_trend$trend + 1)

# Are populations more likely to be declining where there are more threats/measures
m2_NULL = lmer(trend ~ (1|featurename), 
            data = model_data_trend) 
# Are populations more likely to be declining where there are more threats/measures
m2_threat = lmer(trend ~  scale(n_hl_threat) + (1|featurename), 
            data = model_data_trend) 
m2_measure = lmer(trend ~  scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m2_pressuee = lmer(trend ~  scale(n_hl_pressure) + (1|featurename), 
            data = model_data_trend) 
m2A = lmer(trend ~  scale(n_hl_threat) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m2B = lmer(trend ~  scale(n_hl_threat) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m2C = lmer(trend ~  scale(n_hl_pressure) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
m2D = lmer(trend ~  scale(n_hl_pressure) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_trend) 
# Adding 'country' throws singularities due to sparse data
#summary(m1)

tab_model(m2_NULL, m2_threat, m2_measure, m2A, m1B, m2C, m2D, show.aicc = TRUE)
  trend trend trend trend trend trend trend
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.51 0.43 – 0.59 <0.001 0.50 0.42 – 0.58 <0.001 0.51 0.43 – 0.59 <0.001 0.51 0.42 – 0.59 <0.001 1.81 1.37 – 2.24 <0.001 0.51 0.43 – 0.59 <0.001 0.51 0.44 – 0.59 <0.001
n_hl_threat 0.01 -0.03 – 0.05 0.535 0.01 -0.04 – 0.06 0.660 -0.04 -0.41 – 0.34 0.850
n_hl_measure 0.01 -0.03 – 0.05 0.636 0.00 -0.05 – 0.05 0.879 0.09 -0.30 – 0.47 0.663 0.01 -0.03 – 0.06 0.611 0.01 -0.03 – 0.05 0.703
n_hl_threat *
n_hl_measure 		-0.20 -0.49 – 0.09 0.173
n_hl_pressure -0.01 -0.05 – 0.04 0.812 0.00 -0.04 – 0.04 0.936
n_hl_pressure *
n_hl_measure 		-0.03 -0.05 – -0.00 0.036
Random Effects
σ2 0.00 0.00 0.00 0.00 0.46 0.00 0.00
τ00 0.01 featurename 0.01 featurename 0.01 featurename 0.01 featurename 0.25 featurename 0.01 featurename 0.01 featurename
ICC 0.80 0.80 0.81 0.79 0.35 0.80 0.79
N 8 featurename 8 featurename 8 featurename 8 featurename 10 featurename 8 featurename 8 featurename
Observations 21 21 21 21 36 21 21
Marginal R2 / Conditional R2 0.000 / 0.802 0.010 / 0.799 0.007 / 0.807 0.012 / 0.794 0.049 / 0.378 0.006 / 0.803 0.065 / 0.803
AIC -34.085 -25.378 -25.335 -16.366 102.790 -16.117 -9.348 
compare_performance(m2_NULL, m2_threat, m2_measure, m2A, m2B, m2C, m2D)
## # Comparison of Model Performance Indices
## 
## Model      |    Type |    AIC |    BIC | R2_conditional | R2_marginal |  ICC | RMSE |   BF
## ------------------------------------------------------------------------------------------
## m2_NULL    | lmerMod | -35.50 | -32.36 |           0.80 |        0.00 | 0.80 | 0.04 |     
## m2_threat  | lmerMod | -27.88 | -23.70 |           0.80 |        0.01 | 0.80 | 0.04 | 0.01
## m2_measure | lmerMod | -27.84 | -23.66 |           0.81 |    7.00e-03 | 0.81 | 0.04 | 0.01
## m2A        | lmerMod | -20.37 | -15.14 |           0.79 |        0.01 | 0.79 | 0.04 | 0.00
## m2B        | lmerMod | -16.11 |  -9.85 |           0.82 |        0.07 | 0.80 | 0.04 | 0.00
## m2C        | lmerMod | -20.12 | -14.89 |           0.80 |    6.00e-03 | 0.80 | 0.04 | 0.00
## m2D        | lmerMod | -15.35 |  -9.08 |           0.80 |        0.07 | 0.79 | 0.04 | 0.00

No clear support for threats/measures helping predict long-term increases - null model preferred

plot_model(m2_NULL, show.intercept = TRUE)

tab_model(m2_NULL)
  trend
Predictors Estimates CI p
(Intercept) 0.51 0.43 – 0.59 <0.001
Random Effects
σ2 0.00
τ00 featurename 0.01
ICC 0.80
N featurename 8
Observations 21
Marginal R2 / Conditional R2 0.000 / 0.802 
plot_model(m2_NULL, type = "pred")
## list()
sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] performance_0.5.0 tidyr_1.1.3       sjPlot_2.8.5      lme4_1.1-27.1    
##  [5] Matrix_1.2-18     ggplot2_3.3.3     data.table_1.14.0 expss_0.10.7     
##  [9] dplyr_1.0.6       plyr_1.8.6       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6         mvtnorm_1.1-1      lattice_0.20-41    zoo_1.8-8         
##  [5] assertthat_0.2.1   digest_0.6.27      utf8_1.2.1         R6_2.5.0          
##  [9] backports_1.2.1    coda_0.19-4        evaluate_0.14      highr_0.9         
## [13] pillar_1.6.0       rlang_0.4.11       multcomp_1.4-14    rstudioapi_0.13   
## [17] minqa_1.2.4        nloptr_1.2.2.2     checkmate_2.0.0    effectsize_0.3.3  
## [21] rmarkdown_2.2      ggeffects_0.16.0   splines_4.0.1      stringr_1.4.0     
## [25] foreign_0.8-80     htmlwidgets_1.5.3  munsell_0.5.0      broom_0.7.8       
## [29] modelr_0.1.8       compiler_4.0.1     xfun_0.22          pkgconfig_2.0.3   
## [33] parameters_0.8.6   htmltools_0.5.1.1  insight_0.9.6      tidyselect_1.1.1  
## [37] tibble_3.1.1       htmlTable_2.1.0    codetools_0.2-16   matrixStats_0.58.0
## [41] fansi_0.4.2        crayon_1.4.1       withr_2.4.2        MASS_7.3-51.6     
## [45] sjmisc_2.8.5       grid_4.0.1         xtable_1.8-4       nlme_3.1-148      
## [49] gtable_0.3.0       lifecycle_1.0.0    DBI_1.1.1          magrittr_2.0.1    
## [53] bayestestR_0.7.2   scales_1.1.1       estimability_1.3   stringi_1.5.3     
## [57] farver_2.1.0       ellipsis_0.3.2     generics_0.1.0     vctrs_0.3.8       
## [61] sandwich_3.0-0     boot_1.3-25        sjlabelled_1.1.7   TH.data_1.0-10    
## [65] RColorBrewer_1.1-2 tools_4.0.1        glue_1.4.2         purrr_0.3.4       
## [69] sjstats_0.18.0     emmeans_1.5.1      survival_3.1-12    parallel_4.0.1    
## [73] yaml_2.2.1         colorspace_2.0-0   knitr_1.33