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)

Which species have population trend data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$population_trend)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           2  0  3  0  0   1
##   Canis aureus            0  0 11  4  0   5
##   Canis lupus             5  1 23 13  2   1
##   Capra ibex              2  1  3  1  0   0
##   Capra pyrenaica         0  0  6  0  0   0
##   Castor fiber            0  1 31  4  0   0
##   Gulo gulo               0  0  3  1  0   0
##   Halichoerus grypus      1  0 13  2  1   2
##   Lutra lutra             1  0 25 21  1   6
##   Lynx lynx               2  4  5 13  1   8
##   Lynx pardinus           0  0  2  0  0   0
##   Martes martes           0  0  8 22  3  20
##   Megaptera novaeangliae  2  0  0  0  1  10
##   Myotis daubentonii      0  4  4 24  7  11
##   Myotis emarginatus      1  5 12 16  1   7
##   Myotis nattereri        0  0  9 18  9  14
##   Phoca vitulina          1  0  9  1  1   3
##   Pusa hispida            1  0  0  0  0   1
##   Rupicapra pyrenaica     0  0  3  3  0   0
##   Rupicapra rupicapra     3  0  1 10  0   3
##   Ursus arctos            0  3 12 10  1   0

Which species have long-term population trend data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$population_trend_long)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           4  0  2  0  0   0
##   Canis aureus           16  0  3  1  0   0
##   Canis lupus            32  1  9  2  1   0
##   Capra ibex              4  0  3  0  0   0
##   Capra pyrenaica         4  0  2  0  0   0
##   Castor fiber           21  0 15  0  0   0
##   Gulo gulo               4  0  0  0  0   0
##   Halichoerus grypus     12  0  7  0  0   0
##   Lutra lutra            36  1 16  1  0   0
##   Lynx lynx              29  1  1  2  0   0
##   Lynx pardinus           1  0  1  0  0   0
##   Martes martes          40  0  6  6  0   1
##   Megaptera novaeangliae  9  0  0  0  0   4
##   Myotis daubentonii     32  1  0  7  2   8
##   Myotis emarginatus     30  2  6  1  0   3
##   Myotis nattereri       35  1  4  5  0   5
##   Phoca vitulina          8  0  6  0  1   0
##   Pusa hispida            2  0  0  0  0   0
##   Rupicapra pyrenaica     2  0  2  2  0   0
##   Rupicapra rupicapra    13  0  3  1  0   0
##   Ursus arctos           19  0  6  1  0   0

Which species have regional range change data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$range_trend)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           0  0  3  0  0   3
##   Canis aureus            0  0 14  4  0   2
##   Canis lupus             4  0 23 15  2   1
##   Capra ibex              2  0  1  4  0   0
##   Capra pyrenaica         0  0  2  4  0   0
##   Castor fiber            0  0 32  4  0   0
##   Gulo gulo               0  0  2  2  0   0
##   Halichoerus grypus      3  0  1 15  0   0
##   Lutra lutra             1  0 25 25  0   3
##   Lynx lynx               2  3  6 15  2   5
##   Lynx pardinus           0  0  2  0  0   0
##   Martes martes           0  0  8 37  0   8
##   Megaptera novaeangliae  4  0  0  2  0   7
##   Myotis daubentonii      0  0  3 42  0   5
##   Myotis emarginatus      1  1  2 32  0   6
##   Myotis nattereri        0  1  4 36  3   6
##   Phoca vitulina          2  0  3  9  0   1
##   Pusa hispida            2  0  0  0  0   0
##   Rupicapra pyrenaica     0  0  2  4  0   0
##   Rupicapra rupicapra     3  0  1 12  0   1
##   Ursus arctos            0  2 10 14  0   0

Which species have regional long-term range change data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$range_trend_long)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           4  0  2  0  0   0
##   Canis aureus           16  0  4  0  0   0
##   Canis lupus            30  0 11  3  0   1
##   Capra ibex              4  0  1  2  0   0
##   Capra pyrenaica         3  0  2  1  0   0
##   Castor fiber           21  0 15  0  0   0
##   Gulo gulo               3  0  1  0  0   0
##   Halichoerus grypus     14  0  3  1  1   0
##   Lutra lutra            35  1 16  2  0   0
##   Lynx lynx              28  1  2  2  0   0
##   Lynx pardinus           1  0  1  0  0   0
##   Martes martes          40  0  4  9  0   0
##   Megaptera novaeangliae  7  0  0  0  1   5
##   Myotis daubentonii     34  0  0 12  0   4
##   Myotis emarginatus     29  1  4  3  0   5
##   Myotis nattereri       37  0  2  7  0   4
##   Phoca vitulina         10  0  1  3  1   0
##   Pusa hispida            2  0  0  0  0   0
##   Rupicapra pyrenaica     3  0  1  2  0   0
##   Rupicapra rupicapra    12  0  2  3  0   0
##   Ursus arctos           19  1  5  1  0   0

Which species have regional habitat change data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$habitat_trend)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           0  0  1  4  0   1
##   Canis aureus            2  0  6 12  0   0
##   Canis lupus             8  1 10 26  0   0
##   Capra ibex              2  0  0  5  0   0
##   Capra pyrenaica         0  0  0  6  0   0
##   Castor fiber            0  0 14 20  1   1
##   Gulo gulo               0  0  0  4  0   0
##   Halichoerus grypus      2  0  0 13  1   3
##   Lutra lutra             1  1  4 41  1   6
##   Lynx lynx               2  1  1 19  1   9
##   Lynx pardinus           0  0  0  1  1   0
##   Martes martes           0  1  2 29  5  16
##   Megaptera novaeangliae  4  0  0  1  0   8
##   Myotis daubentonii      0  5  0 29  0  16
##   Myotis emarginatus      1 12  0 20  1   8
##   Myotis nattereri        0  2  2 29  3  14
##   Phoca vitulina          2  0  1  8  0   4
##   Pusa hispida            1  0  0  0  0   1
##   Rupicapra pyrenaica     0  0  0  3  0   3
##   Rupicapra rupicapra     3  0  0  8  0   6
##   Ursus arctos            1  0  2 22  0   1

Which species have regional long-term habitat change data available?

# Population trend vs species 
table(mammal_species_short$speciesname, mammal_species_short$habitat_trend_long)
##                         
##                              D  I  S  U Unk
##   Bison bonasus           4  0  1  1  0   0
##   Canis aureus           18  0  0  2  0   0
##   Canis lupus            35  0  4  6  0   0
##   Capra ibex              5  0  0  1  1   0
##   Capra pyrenaica         3  1  2  0  0   0
##   Castor fiber           26  0  2  4  0   4
##   Gulo gulo               4  0  0  0  0   0
##   Halichoerus grypus     15  0  0  3  1   0
##   Lutra lutra            39  1  6  6  1   1
##   Lynx lynx              30  0  1  2  0   0
##   Lynx pardinus           2  0  0  0  0   0
##   Martes martes          41  0  3  9  0   0
##   Megaptera novaeangliae  9  0  0  0  1   3
##   Myotis daubentonii     34  1  0  8  1   6
##   Myotis emarginatus     32  3  0  4  0   3
##   Myotis nattereri       38  1  0  5  0   6
##   Phoca vitulina          8  0  1  5  1   0
##   Pusa hispida            2  0  0  0  0   0
##   Rupicapra pyrenaica     3  0  0  2  0   1
##   Rupicapra rupicapra    13  0  1  0  0   3
##   Ursus arctos           19  0  0  6  0   1

##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)


mammal_threats_n
## # A tibble: 7,311 x 8
## # Groups:   country, region, featurename, featurecode [7,311]
##    country region featurename      featurecode n_pressure n_threat n_hl_pressure
##    <chr>   <chr>  <chr>                  <dbl>      <int>    <int>         <int>
##  1 AT      ALP    Acipenser ruthe…        2487          2        2             2
##  2 AT      ALP    Adenophora lili…        4068          6        6             5
##  3 AT      ALP    Alburnus mento          5289          4        4             3
##  4 AT      ALP    Anisus vorticul…        4056          3        4             3
##  5 AT      ALP    Apium repens            1614          2        2             1
##  6 AT      ALP    Aquilegia alpina        1480          2        2             2
##  7 AT      ALP    Arnica montana          1762          6        6             4
##  8 AT      ALP    Artemisia genipi        1764          2        2             2
##  9 AT      ALP    Aspius aspius           1130          2        2             2
## 10 AT      ALP    Asplenium adult…        4066          3        3             3
## # … with 7,301 more rows, and 1 more variable: n_hl_threat <int>
# Couple of quick plots to look at relationship between n_threat and n_hl_threat 
ggplot(mammal_threats_n, aes(x = n_threat, y = n_hl_threat)) + 
  geom_point(alpha = 0.05, size = 5) + 
  ggtitle("Number of unique threats vs. number of unique 'high-level' threats") + 
  xlab("Number of unique threats") + 
  ylab("Number of unique high-level threats") + 
  ggpubr::theme_classic2()
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

ggplot(mammal_threats_n, aes(x = n_pressure, y = n_hl_pressure)) + 
  geom_point(alpha = 0.05, size = 5) + 
  ggtitle("Number of unique pressures vs. number of unique 'high-level' pressures") + 
  xlab("Number of unique pressures") + 
  ylab("Number of unique high-level pressures") + 
  ggpubr::theme_classic2()

# 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")


# ** Do we need to filter out 'no-measure' codes? 
#mammal_threats = subset(mammal_threats, measurecode != "Xxt") ?

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.
mammal_measures_n
## # A tibble: 4,671 x 6
## # Groups:   country, region, featurename [4,671]
##    country region featurename                 featurecode n_measure n_hl_measure
##    <chr>   <chr>  <chr>                             <dbl>     <int>        <int>
##  1 AT      ALP    Adenophora lilifolia               4068         1            1
##  2 AT      ALP    Alburnus mento                     5289         2            2
##  3 AT      ALP    Anisus vorticulus                  4056         1            1
##  4 AT      ALP    Apium repens                       1614         1            1
##  5 AT      ALP    Aspius aspius                      1130         2            2
##  6 AT      ALP    Astacus astacus                    1091         7            5
##  7 AT      ALP    Austropotamobius pallipes          1092         1            1
##  8 AT      ALP    Austropotamobius torrentium        1093         6            4
##  9 AT      ALP    Barbastella barbastellus           1308         2            1
## 10 AT      ALP    Barbus barbus                      5085         3            3
## # … with 4,661 more rows
# Couple of quick plots to look at relationship between n_threat and n_hl_threat 
ggplot(mammal_measures_n, aes(x = n_measure, y = n_hl_measure)) + 
  geom_point(alpha = 0.05, size = 5) + 
  ggtitle("Number of unique measures vs. number of unique 'high-level' measures (mammals)") + 
  xlab("Number of unique threats") + 
  ylab("Number of unique high-level threats") + 
  ggpubr::theme_classic2()

##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)

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

# 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)

How many species do we end up with?

length(unique(full_data$featurename))
## [1] 21

How many country records (threats/measures) do we have for each species?

table(full_data$featurename)
## 
##          Bison bonasus           Canis aureus            Canis lupus 
##                      6                     20                     45 
##             Capra ibex        Capra pyrenaica           Castor fiber 
##                      7                      6                     36 
##              Gulo gulo     Halichoerus grypus            Lutra lutra 
##                      4                     19                     54 
##              Lynx lynx          Lynx pardinus          Martes martes 
##                     33                      2                     53 
## Megaptera novaeangliae     Myotis daubentonii     Myotis emarginatus 
##                     13                     50                     42 
##       Myotis nattereri         Phoca vitulina           Pusa hispida 
##                     50                     15                      2 
##    Rupicapra pyrenaica    Rupicapra rupicapra           Ursus arctos 
##                      6                     17                     26

Combine countries to regions

Population/habitat change data is given per regions, so aggregate n_threat/n_measure to region

# Average threats over countries to regions as that's where we have response data
# model_data = full_data %>%
#   group_by(region, featurename, featurecode, population_trend, habitat_trend) %>% 
#   summarise(
#     mean_n_pressure = mean(n_pressure),
#     mean_n_threat = mean(n_threat),
#     mean_n_hl_pressure = mean(n_hl_pressure),
#     mean_n_hl_threat = mean(n_hl_threat),
#     mean_n_measure = mean(n_measure),
#     mean_n_hl_measure = mean(n_hl_measure),
#   )
model_data = full_data

How many species are there per region?

# How many do we have
addmargins(table(model_data$featurename, model_data$region))
##                         
##                          ALP ATL BLS BOR CON MATL MBAL MED MMAC MMED PAN STE
##   Bison bonasus            3   0   0   0   3    0    0   0    0    0   0   0
##   Canis aureus             5   0   2   0   7    0    0   2    0    0   3   1
##   Canis lupus             12   6   1   5  13    0    0   6    0    0   2   0
##   Capra ibex               5   0   0   0   1    0    0   1    0    0   0   0
##   Capra pyrenaica          2   2   0   0   0    0    0   2    0    0   0   0
##   Castor fiber             7   6   0   5  11    0    0   2    0    0   4   1
##   Gulo gulo                2   0   0   2   0    0    0   0    0    0   0   0
##   Halichoerus grypus       0   0   0   0   0   10    8   0    1    0   0   0
##   Lutra lutra             13   9   2   5  14    0    0   6    0    0   4   1
##   Lynx lynx               11   1   1   5  10    0    0   3    0    0   2   0
##   Lynx pardinus            0   0   0   0   0    0    0   2    0    0   0   0
##   Martes martes           13   9   1   5  14    0    0   6    0    0   4   1
##   Megaptera novaeangliae   0   0   0   0   0    6    0   0    2    5   0   0
##   Myotis daubentonii      11   9   1   5  14    0    0   6    0    0   3   1
##   Myotis emarginatus      11   7   1   0  12    0    0   7    0    0   3   1
##   Myotis nattereri        11   9   1   5  14    0    0   7    0    0   3   0
##   Phoca vitulina           0   0   0   0   0   10    4   0    1    0   0   0
##   Pusa hispida             0   0   0   0   0    2    0   0    0    0   0   0
##   Rupicapra pyrenaica      2   2   0   0   0    0    0   2    0    0   0   0
##   Rupicapra rupicapra      7   0   0   0   7    0    0   3    0    0   0   0
##   Ursus arctos            12   1   0   4   5    0    0   4    0    0   0   0
##   Sum                    127  61  10  41 125   28   12  59    4    5  28   6
##                         
##                          Sum
##   Bison bonasus            6
##   Canis aureus            20
##   Canis lupus             45
##   Capra ibex               7
##   Capra pyrenaica          6
##   Castor fiber            36
##   Gulo gulo                4
##   Halichoerus grypus      19
##   Lutra lutra             54
##   Lynx lynx               33
##   Lynx pardinus            2
##   Martes martes           53
##   Megaptera novaeangliae  13
##   Myotis daubentonii      50
##   Myotis emarginatus      42
##   Myotis nattereri        50
##   Phoca vitulina          15
##   Pusa hispida             2
##   Rupicapra pyrenaica      6
##   Rupicapra rupicapra     17
##   Ursus arctos            26
##   Sum                    506

First simple model.. do species appear to have more conservation measures in place when there are more threats reported?

model_data_measure = model_data %>%
  mutate(population_trend = factor(population_trend)) %>%
  filter(population_trend != "") %>%  # Subset to only increases or declines
  select(region, n_hl_threat, n_hl_measure, population_trend_long, featurename) %>% # Columns for modelling
  drop_na() %>% # Get rid of NAs in any of these columns
  droplevels()


#*** Including species ('featurename') as a random effect is tricky (throw singularities) as there's only a couple
# of records for some species. 
m1_measures_NULL = lmer(n_hl_measure ~ (1|region), 
                        data = model_data_measure) 
m1_measuresA = lmer(n_hl_measure ~ population_trend_long + n_hl_threat + (1|region), 
                    data = model_data_measure) 
m1_measuresB = lmer(n_hl_measure ~ population_trend_long * n_hl_threat + (1|region), 
                    data = model_data_measure) 

# Adding 'country' throws singularities due to sparse data
#summary(m1_measures)
tab_model(m1_measures_NULL, m1_measuresA, m1_measuresB, show.aicc = TRUE)
  n_hl_measure n_hl_measure n_hl_measure
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 3.27 2.68 – 3.87 <0.001 1.52 0.95 – 2.09 <0.001 1.59 0.91 – 2.27 <0.001
population_trend_long [D] -0.20 -1.65 – 1.25 0.788 -2.46 -5.83 – 0.90 0.152
population_trend_long [I] -0.26 -0.70 – 0.17 0.228 -0.36 -1.38 – 0.66 0.493
population_trend_long [S] -0.46 -1.69 – 0.77 0.466 -1.65 -5.02 – 1.72 0.338
population_trend_long [U] -1.10 -3.14 – 0.95 0.294 3.05 -2.17 – 8.28 0.252
population_trend_long
[Unk]
-0.37 -1.21 – 0.47 0.394 -0.28 -2.59 – 2.02 0.808
n_hl_threat 0.46 0.36 – 0.57 <0.001 0.44 0.31 – 0.58 <0.001
population_trend_long [D]
* n_hl_threat
0.64 -0.22 – 1.51 0.146
population_trend_long [I]
* n_hl_threat
0.02 -0.19 – 0.24 0.817
population_trend_long [S]
* n_hl_threat
0.22 -0.34 – 0.78 0.448
population_trend_long [U]
* n_hl_threat
-1.18 -2.54 – 0.18 0.089
population_trend_long
[Unk] * n_hl_threat
-0.04 -0.72 – 0.65 0.920
Random Effects
σ2 2.64 2.10 2.09
τ00 0.80 region 0.25 region 0.27 region
ICC 0.23 0.10 0.11
N 12 region 12 region 12 region
Observations 265 265 265
Marginal R2 / Conditional R2 0.000 / 0.233 0.244 / 0.323 0.252 / 0.337
AIC 1035.007 975.498 983.129
compare_performance(m1_measures_NULL, m1_measuresA, m1_measuresB)
## # Comparison of Model Performance Indices
## 
## Name             |   Model |      AIC |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
## --------------------------------------------------------------------------------------------------
## m1_measures_NULL | lmerMod | 1034.915 | 1045.654 |      0.233 |      0.000 | 0.233 | 1.598 | 1.625
## m1_measuresA     | lmerMod |  974.792 | 1007.010 |      0.323 |      0.244 | 0.105 | 1.413 | 1.449
## m1_measuresB     | lmerMod |  981.449 | 1031.565 |      0.337 |      0.252 | 0.114 | 1.395 | 1.446
# 
# library(rstanarm)
# sm1_measures_NULL = stan_glmer(n_hl_measure ~ (1|featurename), 
#                         data = model_data_measure, diagnostic_file = file.path(tempdir(), "df0.csv"))
# sm1_measuresA = stan_glmer(n_hl_measure ~ population_trend + n_hl_threat + (1|featurename), 
#                     data = model_data_measure, diagnostic_file = file.path(tempdir(), "df1.csv"))
# sm1_measuresB = stan_glmer(n_hl_measure ~ population_trend * n_hl_threat + (1|featurename), 
#                     data = model_data_measure, diagnostic_file = file.path(tempdir(), "df2.csv"))
# 
# 
# library(bridgesampling)
# bridge_0 <- bridge_sampler(sm1_measures_NULL)
# bridge_1 <- bridge_sampler(sm1_measuresA)
# bridge_2 <- bridge_sampler(sm1_measuresB)
# 
# bf(bridge_1, bridge_0)
# bf(bridge_2, bridge_1)
# 
# plot(sm1_measuresA)

In general, popualtions with more threats have significantly more measures, but this doesn’t appear to vary with the type of population trend (e.g. declines, increases, etc)

plot_model(m1_measuresA, show.intercept = TRUE)

tab_model(m1_measuresA) 
  n_hl_measure
Predictors Estimates CI p
(Intercept) 1.52 0.95 – 2.09 <0.001
population_trend_long [D] -0.20 -1.65 – 1.25 0.788
population_trend_long [I] -0.26 -0.70 – 0.17 0.228
population_trend_long [S] -0.46 -1.69 – 0.77 0.466
population_trend_long [U] -1.10 -3.14 – 0.95 0.294
population_trend_long
[Unk]
-0.37 -1.21 – 0.47 0.394
n_hl_threat 0.46 0.36 – 0.57 <0.001
Random Effects
σ2 2.10
τ00 region 0.25
ICC 0.10
N region 12
Observations 265
Marginal R2 / Conditional R2 0.244 / 0.323
plot_model(m1_measuresA, type = "pred")
## $population_trend_long

## 
## $n_hl_threat

#plot_model(m1_measuresB, type = "pred", terms = c("mean_n_hl_threat", "population_trend")) # Hard to see significant effects on orediction plot because of broad CIs of D pops

Declines vs Threats/Measures

Are species more likely to be declining where there are more threats/measures?

# Model data is only declines or increases (remapped to 'increasing' as 0/1)
model_data_binom = model_data %>%
  ungroup() %>%
  mutate(population_trend = factor(population_trend)) %>%
  filter(population_trend == "D" | population_trend == "I" | population_trend == "S") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(population_trend == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(increasing, n_hl_threat, n_hl_pressure, n_hl_measure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns


# Are populations more likely to be declining where there are more threats/measures
m1_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
# Are populations more likely to be declining where there are more threats/measures
m1_threat = glmer(increasing ~  scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1_measure = glmer(increasing ~  scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1_pressuee = glmer(increasing ~  scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1A = glmer(increasing ~  scale(n_hl_threat) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1B = glmer(increasing ~  scale(n_hl_threat) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1C = glmer(increasing ~  scale(n_hl_pressure) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1D = glmer(increasing ~  scale(n_hl_pressure) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
# 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)
  increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 1.49 0.63 – 3.53 0.368 1.51 0.64 – 3.59 0.348 1.46 0.61 – 3.46 0.393 1.51 0.64 – 3.57 0.346 1.51 0.63 – 3.63 0.358 1.50 0.64 – 3.50 0.351 1.49 0.63 – 3.55 0.368
n_hl_threat 1.12 0.82 – 1.55 0.473 1.43 0.96 – 2.13 0.077 1.43 0.95 – 2.17 0.087
n_hl_measure 0.81 0.59 – 1.10 0.179 0.67 0.45 – 0.98 0.037 0.66 0.45 – 0.98 0.038 0.67 0.46 – 0.98 0.040 0.67 0.46 – 0.98 0.041
n_hl_threat *
n_hl_measure
1.00 0.73 – 1.37 0.981
n_hl_pressure 1.42 0.95 – 2.13 0.088 1.42 0.94 – 2.16 0.095
n_hl_pressure *
n_hl_measure
1.01 0.74 – 1.38 0.961
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 2.35 featurename 2.34 featurename 2.36 featurename 2.30 featurename 2.30 featurename 2.24 featurename 2.23 featurename
ICC 0.42 0.42 0.42 0.41 0.41 0.40 0.40
N 18 featurename 18 featurename 18 featurename 18 featurename 18 featurename 18 featurename 18 featurename
Observations 229 229 229 229 229 229 229
Marginal R2 / Conditional R2 0.000 / 0.417 0.002 / 0.417 0.008 / 0.422 0.022 / 0.424 0.022 / 0.424 0.021 / 0.417 0.021 / 0.417
AIC 286.813 288.361 287.068 285.988 288.078 286.221 288.309
compare_performance(m1_NULL, m1_threat, m1_measure, m1A, m1B, m1C, m1D)
## # Comparison of Model Performance Indices
## 
## Name       |    Model |     AIC |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## ------------------------------------------------------------------------------------------------------------------------------------
## m1_NULL    | glmerMod | 286.760 | 293.627 |      0.417 |      0.000 | 0.417 | 0.421 | 1.000 |    0.533 |  -112.100 |           0.007
## m1_threat  | glmerMod | 288.254 | 298.555 |      0.417 |      0.002 | 0.416 | 0.421 | 1.000 |    0.532 |  -112.337 |           0.009
## m1_measure | glmerMod | 286.961 | 297.262 |      0.422 |      0.008 | 0.417 | 0.419 | 1.000 |    0.529 |  -112.968 |           0.005
## m1A        | glmerMod | 285.809 | 299.544 |      0.424 |      0.022 | 0.411 | 0.417 | 1.000 |    0.523 |  -114.723 |           0.007
## m1B        | glmerMod | 287.809 | 304.977 |      0.424 |      0.022 | 0.411 | 0.417 | 1.000 |    0.523 |  -114.720 |           0.007
## m1C        | glmerMod | 286.042 | 299.777 |      0.417 |      0.021 | 0.405 | 0.418 | 1.000 |    0.525 |  -114.417 |           0.008
## m1D        | glmerMod | 288.040 | 305.208 |      0.417 |      0.021 | 0.404 | 0.418 | 1.000 |    0.525 |  -114.404 |           0.007
No clear support for threats/measures helping predict increases - null model preferred

Support for m1A (including both threats and measures, but no interaction)

plot_model(m1A, show.intercept = TRUE)

tab_model(m1A) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 1.51 0.64 – 3.57 0.346
n_hl_threat 1.43 0.96 – 2.13 0.077
n_hl_measure 0.67 0.45 – 0.98 0.037
Random Effects
σ2 3.29
τ00 featurename 2.30
ICC 0.41
N featurename 18
Observations 229
Marginal R2 / Conditional R2 0.022 / 0.424
plot_model(m1A, type = "pred")
## $n_hl_threat

## 
## $n_hl_measure

Long-term change

Are species more likely to be declining where there are more threats/measures?

# Model data is only declines or increases (remapped to 'increasing' as 0/1)
model_data_binom_long = model_data %>%
  ungroup() %>%
  mutate(population_trend_long = factor(population_trend_long)) %>%
  filter(population_trend_long == "D" | population_trend_long == "I") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(population_trend_long == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(increasing, n_hl_threat, n_hl_measure, n_hl_pressure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns


# Are populations more likely to be declining where there are more threats/measures
m2_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
# Are populations more likely to be declining where there are more threats/measures
m2_threat = glmer(increasing ~  scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2_measure = glmer(increasing ~  scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2_pressure = glmer(increasing ~  scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2A = glmer(increasing ~  scale(n_hl_threat) + scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2B = glmer(increasing ~  scale(n_hl_threat) * scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2C = glmer(increasing ~  scale(n_hl_threat) + scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2D = glmer(increasing ~  scale(n_hl_threat) * scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_long, family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
## boundary (singular) fit: see ?isSingular
# Adding 'country' throws singularities due to sparse data
#summary(m1)

tab_model(m2_NULL, m2_threat, m2_measure, m2A, m2B, m2C, m2D, show.aicc = TRUE)
  increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 21.10 3.21 – 138.50 0.001 20.82 3.39 – 127.92 0.001 20.97 3.22 – 136.31 0.001 21.75 3.31 – 143.07 0.001 26.58 2.57 – 275.34 0.006 20.36 3.26 – 127.04 0.001 11.35 3.08 – 41.78 <0.001
n_hl_threat 1.56 0.51 – 4.75 0.433 1.98 0.43 – 9.22 0.382 1.92 0.37 – 10.01 0.441 1.87 0.11 – 30.67 0.660 2.27 0.16 – 31.58 0.541
n_hl_measure 1.06 0.35 – 3.15 0.922 0.71 0.18 – 2.86 0.628 0.69 0.15 – 3.03 0.618
n_hl_threat *
n_hl_measure
0.81 0.22 – 2.97 0.745
n_hl_pressure 0.82 0.05 – 12.90 0.889 1.11 0.10 – 12.12 0.930
n_hl_threat *
n_hl_pressure
2.15 0.38 – 12.10 0.384
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 0.75 featurename 0.46 featurename 0.73 featurename 0.46 featurename 0.66 featurename 0.39 featurename 0.00 featurename
ICC 0.19 0.12 0.18 0.12 0.17 0.11  
N 12 featurename 12 featurename 12 featurename 12 featurename 12 featurename 12 featurename 12 featurename
Observations 68 68 68 68 68 68 68
Marginal R2 / Conditional R2 0.000 / 0.186 0.050 / 0.166 0.001 / 0.182 0.078 / 0.192 0.090 / 0.241 0.053 / 0.153 0.274 / NA
AIC 34.449 36.027 36.630 38.051 40.276 38.268 39.637
compare_performance(m2_NULL, m2_threat, m2_measure, m2A, m2B, m2C, m2D)
## Random effect variances not available. Returned R2 does not account for random effects.
## # Comparison of Model Performance Indices
## 
## Name       |    Model |    AIC |    BIC | R2 (cond.) | R2 (marg.) |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   ICC
## ----------------------------------------------------------------------------------------------------------------------------------
## m2_NULL    | glmerMod | 34.265 | 38.704 |      0.186 |      0.000 | 0.228 | 1.000 |    0.193 |  -189.585 |           0.036 | 0.186
## m2_threat  | glmerMod | 35.652 | 42.310 |      0.166 |      0.050 | 0.228 | 1.000 |    0.198 |  -190.175 |           0.015 | 0.122
## m2_measure | glmerMod | 36.255 | 42.914 |      0.182 |  7.336e-04 | 0.228 | 1.000 |    0.194 |  -189.298 |           0.017 | 0.181
## m2A        | glmerMod | 37.416 | 46.294 |      0.192 |      0.078 | 0.229 | 1.000 |    0.197 |  -192.934 |           0.015 | 0.123
## m2B        | glmerMod | 39.308 | 50.406 |      0.241 |      0.090 | 0.227 | 1.000 |    0.190 |  -197.320 |           0.015 | 0.166
## m2C        | glmerMod | 37.633 | 46.511 |      0.153 |      0.053 | 0.229 | 1.000 |    0.201 |  -189.137 |           0.015 | 0.106
## m2D        | glmerMod | 38.669 | 49.767 |            |      0.274 | 0.232 | 1.000 |    0.211 |  -197.271 |           0.015 |

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) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 21.10 3.21 – 138.50 0.001
Random Effects
σ2 3.29
τ00 featurename 0.75
ICC 0.19
N featurename 12
Observations 68
Marginal R2 / Conditional R2 0.000 / 0.186
plot_model(m2_NULL, type = "pred")
## list()

Range-change vs Threats/Measures

Are species ranges more likely to be shrinking where there are more threats/measures?

model_data_binom_range = model_data %>%
  ungroup() %>%
  mutate(range_trend = factor(range_trend)) %>%
  filter(range_trend == "D" | range_trend == "I" | range_trend == "S") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(range_trend == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(region,increasing, n_hl_measure, n_hl_threat, n_hl_pressure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns


# Are species ranges more likely to be shrinking where there are more threats/measures
m1r_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom_range, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1r_threat = glmer(increasing ~ scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom_range, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1r_measure = glmer(increasing ~ scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_range, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1r_pressure = glmer(increasing ~ scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_range, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1rA = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_threat) + (1|featurename), 
           data = model_data_binom_range, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1rB = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_threat) + (1|featurename), 
          data = model_data_binom_range, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  
m1rC = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_pressure) + (1|featurename), 
           data = model_data_binom_range, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1rD = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_pressure) + (1|featurename), 
          data = model_data_binom_range, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  

tab_model(m1r_NULL, m1r_threat, m1r_measure, m1rA, m1rB, m1rC, m1rD, show.aicc = TRUE)
  increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 0.39 0.17 – 0.91 0.029 0.41 0.18 – 0.92 0.031 0.39 0.17 – 0.91 0.029 0.40 0.18 – 0.92 0.031 0.35 0.15 – 0.80 0.013 0.40 0.18 – 0.91 0.029 0.35 0.15 – 0.80 0.013
n_hl_threat 1.24 0.90 – 1.71 0.193 1.33 0.91 – 1.92 0.137 1.50 1.00 – 2.25 0.052
n_hl_measure 1.00 0.74 – 1.37 0.976 0.88 0.61 – 1.25 0.470 0.81 0.56 – 1.18 0.278 0.88 0.61 – 1.26 0.472 0.82 0.56 – 1.19 0.300
n_hl_measure *
n_hl_threat
1.32 0.96 – 1.82 0.092
n_hl_pressure 1.32 0.91 – 1.92 0.146 1.47 0.98 – 2.21 0.064
n_hl_measure *
n_hl_pressure
1.28 0.94 – 1.76 0.121
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 2.35 featurename 2.19 featurename 2.35 featurename 2.23 featurename 2.15 featurename 2.18 featurename 2.09 featurename
ICC 0.42 0.40 0.42 0.40 0.39 0.40 0.39
N 19 featurename 19 featurename 19 featurename 19 featurename 19 featurename 19 featurename 19 featurename
Observations 253 253 253 253 253 253 253
Marginal R2 / Conditional R2 0.000 / 0.417 0.008 / 0.404 0.000 / 0.417 0.010 / 0.410 0.024 / 0.409 0.010 / 0.405 0.022 / 0.403
AIC 280.443 280.833 282.491 282.387 281.639 282.488 282.181
compare_performance(m1r_NULL, m1r_threat, m1r_measure, m1rA, m1rB, m1rC, m1rD)
## # Comparison of Model Performance Indices
## 
## Name        |    Model |     AIC |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## -------------------------------------------------------------------------------------------------------------------------------------
## m1r_NULL    | glmerMod | 280.395 | 287.462 |      0.417 |      0.000 | 0.417 | 0.393 | 1.000 |    0.468 |   -48.904 |           0.023
## m1r_threat  | glmerMod | 280.737 | 291.337 |      0.404 |      0.008 | 0.400 | 0.393 | 1.000 |    0.467 |   -49.219 |           0.019
## m1r_measure | glmerMod | 282.395 | 292.995 |      0.417 |  3.996e-06 | 0.417 | 0.393 | 1.000 |    0.468 |   -48.901 |           0.023
## m1rA        | glmerMod | 282.226 | 296.360 |      0.410 |      0.010 | 0.404 | 0.392 | 1.000 |    0.466 |   -49.406 |           0.019
## m1rB        | glmerMod | 281.396 | 299.063 |      0.409 |      0.024 | 0.395 | 0.389 | 1.000 |    0.461 |   -50.100 |           0.013
## m1rC        | glmerMod | 282.327 | 296.460 |      0.405 |      0.010 | 0.399 | 0.392 | 1.000 |    0.467 |   -49.335 |           0.020
## m1rD        | glmerMod | 281.938 | 299.605 |      0.403 |      0.022 | 0.389 | 0.390 | 1.000 |    0.463 |   -49.940 |           0.014

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

plot_model(m1r_NULL, show.intercept = TRUE)

tab_model(m1r_NULL) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 0.39 0.17 – 0.91 0.029
Random Effects
σ2 3.29
τ00 featurename 2.35
ICC 0.42
N featurename 19
Observations 253
Marginal R2 / Conditional R2 0.000 / 0.417
plot_model(m1r_NULL, type = "pred")
## list()

Are species long-term ranges more likely to be shrinking where there are more threats/measures?

model_data_binom_range_long = model_data %>%
  ungroup() %>%
  mutate(range_trend_long = factor(range_trend_long)) %>%
  filter(range_trend_long == "D" | range_trend_long == "I" | range_trend_long == "S") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(range_trend_long == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(region, increasing, n_hl_measure, n_hl_threat, n_hl_pressure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns


# Are species ranges more likely to be shrinking where there are more threats/measures
m2r_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom_range_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2r_threat = glmer(increasing ~ scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom_range_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2r_measure = glmer(increasing ~ scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_range_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2r_pressure = glmer(increasing ~ scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_range_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2rA = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_threat) + (1|featurename), 
           data = model_data_binom_range_long, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2rB = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_threat) + (1|featurename), 
          data = model_data_binom_range_long, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  
m2rC = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_pressure) + (1|featurename), 
           data = model_data_binom_range_long, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2rD = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_pressure) + (1|featurename), 
          data = model_data_binom_range_long, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  

tab_model(m2r_NULL, m2r_threat, m2r_measure, m2r_pressure, m2rA, m2rB, m2rC, m2rD, show.aicc = TRUE)
  increasing increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 2.97 1.15 – 7.71 0.025 3.09 1.24 – 7.73 0.016 2.98 1.13 – 7.85 0.027 3.06 1.21 – 7.72 0.018 3.27 1.25 – 8.52 0.016 2.76 1.01 – 7.50 0.047 3.14 1.21 – 8.11 0.018 2.40 0.94 – 6.13 0.067
n_hl_threat 1.43 0.72 – 2.85 0.308 2.03 0.78 – 5.32 0.148 2.28 0.86 – 6.05 0.099
n_hl_measure 0.90 0.46 – 1.76 0.751 0.60 0.24 – 1.46 0.260 0.62 0.26 – 1.49 0.283 0.70 0.30 – 1.62 0.405 0.72 0.32 – 1.63 0.427
n_hl_pressure 1.29 0.64 – 2.57 0.474 1.60 0.66 – 3.85 0.295 2.09 0.83 – 5.28 0.118
n_hl_measure *
n_hl_threat
1.41 0.60 – 3.30 0.433
n_hl_measure *
n_hl_pressure
1.87 0.79 – 4.44 0.157
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 1.43 featurename 1.13 featurename 1.52 featurename 1.21 featurename 1.23 featurename 1.08 featurename 1.26 featurename 0.97 featurename
ICC 0.30 0.26 0.32 0.27 0.27 0.25 0.28 0.23
N 14 featurename 14 featurename 14 featurename 14 featurename 14 featurename 14 featurename 14 featurename 14 featurename
Observations 68 68 68 68 68 68 68 68
Marginal R2 / Conditional R2 0.000 / 0.303 0.028 / 0.278 0.002 / 0.317 0.014 / 0.278 0.070 / 0.324 0.079 / 0.308 0.034 / 0.302 0.089 / 0.296
AIC 75.374 76.554 77.468 77.070 77.485 79.193 78.620 78.831
compare_performance(m2r_NULL, m2r_threat, m2r_measure, m2r_pressure, m2rA, m2rB, m2rC, m2rD)
## # Comparison of Model Performance Indices
## 
## Name         |    Model |    AIC |    BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## ------------------------------------------------------------------------------------------------------------------------------------
## m2r_NULL     | glmerMod | 75.189 | 79.628 |      0.303 |      0.000 | 0.303 | 0.359 | 1.000 |    0.410 |   -78.782 |           0.027
## m2r_threat   | glmerMod | 76.179 | 82.838 |      0.278 |      0.028 | 0.257 | 0.361 | 1.000 |    0.418 |   -79.701 |           0.021
## m2r_measure  | glmerMod | 77.093 | 83.751 |      0.317 |      0.002 | 0.315 | 0.357 | 1.000 |    0.406 |   -79.215 |           0.029
## m2r_pressure | glmerMod | 76.695 | 83.353 |      0.278 |      0.014 | 0.268 | 0.362 | 1.000 |    0.417 |   -79.076 |           0.015
## m2rA         | glmerMod | 76.850 | 85.728 |      0.324 |      0.070 | 0.273 | 0.352 | 1.000 |    0.404 |   -83.118 |           0.015
## m2rB         | glmerMod | 78.225 | 89.322 |      0.308 |      0.079 | 0.248 | 0.356 | 1.000 |    0.408 |   -83.012 |           0.015
## m2rC         | glmerMod | 77.986 | 86.864 |      0.302 |      0.034 | 0.278 | 0.358 | 1.000 |    0.410 |   -80.732 |           0.016
## m2rD         | glmerMod | 77.864 | 88.961 |      0.296 |      0.089 | 0.228 | 0.361 | 1.000 |    0.411 |   -82.678 |           0.015

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

plot_model(m2r_NULL, show.intercept = TRUE)

tab_model(m2r_NULL) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 2.97 1.15 – 7.71 0.025
Random Effects
σ2 3.29
τ00 featurename 1.43
ICC 0.30
N featurename 14
Observations 68
Marginal R2 / Conditional R2 0.000 / 0.303
plot_model(m2r_NULL, type = "pred")
## list()

Habitat-change vs Threats/Measures

Are species habitats more likely to be shrinking where there are more threats/measures?

model_data_binom_habitat = model_data %>%
  ungroup() %>%
  mutate(habitat_trend = factor(habitat_trend)) %>%
  filter(habitat_trend == "D" | habitat_trend == "I") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(habitat_trend == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(region, increasing, n_hl_measure, n_hl_threat, n_hl_pressure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns

# Are species ranges more likely to be shrinking where there are more threats/measures
m1h_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom_habitat, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1h_threat = glmer(increasing ~ scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom_habitat, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1h_measure = glmer(increasing ~ scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_habitat, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1h_pressure = glmer(increasing ~ scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_habitat, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1hA = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_threat) + (1|featurename), 
           data = model_data_binom_habitat, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1hB = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_threat) + (1|featurename), 
          data = model_data_binom_habitat, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  
m1hC = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_pressure) + (1|featurename), 
           data = model_data_binom_habitat, 
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m1hD = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_pressure) + (1|featurename), 
          data = model_data_binom_habitat, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  

tab_model(m1h_NULL, m1h_threat, m1h_measure, m1h_pressure, m1hA, m1hB, m1hC, m1hD, show.aicc = TRUE)
  increasing increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 18.57 0.35 – 979.19 0.149 84.89 0.02 – 299367.00 0.287 257.20 0.10 – 682828.69 0.168 34.81 0.18 – 6793.30 0.187 291.09 0.06 – 1522604.31 0.194 1099.33 0.15 – 7808897.08 0.122 1359.38 0.01 – 246625716.90 0.243 461.71 0.04 – 5261651.27 0.198
n_hl_threat 0.25 0.01 – 5.36 0.377 0.86 0.02 – 47.92 0.941 0.12 0.00 – 67.64 0.510
n_hl_measure 0.08 0.00 – 4.10 0.209 0.08 0.00 – 5.71 0.250 0.34 0.01 – 21.46 0.607 0.02 0.00 – 78.35 0.340 0.04 0.00 – 24.49 0.331
n_hl_pressure 0.41 0.06 – 2.65 0.350 3.48 0.03 – 397.31 0.606 2.85 0.03 – 244.72 0.645
n_hl_measure *
n_hl_threat
0.05 0.00 – 4.91 0.199
n_hl_measure *
n_hl_pressure
0.45 0.03 – 6.17 0.553
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 10.54 featurename 26.53 featurename 42.55 featurename 16.50 featurename 45.35 featurename 45.60 featurename 69.16 featurename 42.07 featurename
ICC 0.76 0.89 0.93 0.83 0.93 0.93 0.95 0.93
N 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename 10 featurename
Observations 44 44 44 44 44 44 44 44
Marginal R2 / Conditional R2 0.000 / 0.762 0.060 / 0.896 0.121 / 0.937 0.038 / 0.840 0.121 / 0.941 0.176 / 0.945 0.144 / 0.961 0.150 / 0.938
AIC 35.164 35.806 32.855 36.354 35.272 34.510 34.936 37.202
compare_performance(m1h_NULL, m1h_threat, m1h_measure, m1h_pressure, m1hA, m1hB, m1hC, m1hD)
## # Comparison of Model Performance Indices
## 
## Name         |    Model |    AIC |    BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## ------------------------------------------------------------------------------------------------------------------------------------
## m1h_NULL     | glmerMod | 34.872 | 38.440 |      0.762 |      0.000 | 0.762 | 0.225 | 1.000 |    0.179 |   -69.197 |           0.072
## m1h_threat   | glmerMod | 35.206 | 40.559 |      0.896 |      0.060 | 0.890 | 0.198 | 1.000 |    0.133 |   -88.502 |           0.061
## m1h_measure  | glmerMod | 32.255 | 37.608 |      0.937 |      0.121 | 0.928 | 0.175 | 1.000 |    0.098 |  -118.506 |           0.061
## m1h_pressure | glmerMod | 35.754 | 41.107 |      0.840 |      0.038 | 0.834 | 0.214 | 1.000 |    0.152 |   -77.985 |           0.062
## m1hA         | glmerMod | 34.247 | 41.383 |      0.941 |      0.121 | 0.932 | 0.173 | 1.000 |    0.097 |  -120.294 |           0.061
## m1hB         | glmerMod | 32.931 | 41.852 |      0.945 |      0.176 | 0.933 | 0.157 | 1.000 |    0.079 |  -136.775 |           0.065
## m1hC         | glmerMod | 33.911 | 41.047 |      0.961 |      0.144 | 0.955 | 0.167 | 1.000 |    0.088 |  -148.091 |           0.064
## m1hD         | glmerMod | 35.623 | 44.544 |      0.938 |      0.150 | 0.927 | 0.176 | 1.000 |    0.098 |  -125.941 |           0.062

Some (weak) evidence that populations with more measures in place have declining habitat trends.

plot_model(m1h_measure, show.intercept = TRUE)

tab_model(m1h_measure) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 257.20 0.10 – 682828.69 0.168
n_hl_measure 0.08 0.00 – 4.10 0.209
Random Effects
σ2 3.29
τ00 featurename 42.55
ICC 0.93
N featurename 10
Observations 44
Marginal R2 / Conditional R2 0.121 / 0.937
plot_model(m1h_measure, type = "pred")
## $n_hl_measure

Are species long-term habitats more likely to be shrinking where there are more threats/measures?

model_data_binom_habitat_long = model_data %>%
  ungroup() %>%
  mutate(habitat_trend_long = factor(habitat_trend_long)) %>%
  filter(habitat_trend_long == "D" | habitat_trend_long == "I" | habitat_trend_long == "S") %>%  # Subset to only increases or declines
  mutate(increasing = ifelse(habitat_trend_long == "I", 1, 0)) %>% # Response variable with be 'increasing' or otherwise 'decreasing'
  select(region, increasing, n_hl_measure, n_hl_threat, n_hl_pressure, featurename) %>% # Columns for modelling
  drop_na() # Get rid of NAs in any of these columns

# Are species ranges more likely to be shrinking where there are more threats/measures
m2h_NULL = glmer(increasing ~ (1|featurename), 
            data = model_data_binom_habitat_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2h_threat = glmer(increasing ~ scale(n_hl_threat) + (1|featurename), 
            data = model_data_binom_habitat_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2h_measure = glmer(increasing ~ scale(n_hl_measure) + (1|featurename), 
            data = model_data_binom_habitat_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2h_pressure = glmer(increasing ~ scale(n_hl_pressure) + (1|featurename), 
            data = model_data_binom_habitat_long, 
            family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2hA = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_threat) + (1|featurename), 
           data = model_data_binom_habitat_long,
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2hB = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_threat) + (1|featurename), 
          data = model_data_binom_habitat_long, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  
m2hC = glmer(increasing ~ scale(n_hl_measure) + scale(n_hl_pressure) + (1|featurename), 
           data = model_data_binom_habitat_long,
           family = "binomial", control = glmerControl(optimizer = "bobyqa")) 
m2hD = glmer(increasing ~ scale(n_hl_measure) * scale(n_hl_pressure) + (1|featurename), 
          data = model_data_binom_habitat_long, 
          family = "binomial", control = glmerControl(optimizer = "bobyqa"))  

tab_model(m2h_NULL, m2h_threat, m2h_measure, m2h_pressure, m2hA, m2hB, m2hC, m2hD, show.aicc = TRUE)
  increasing increasing increasing increasing increasing increasing increasing increasing
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 0.46 0.20 – 1.05 0.065 0.45 0.19 – 1.05 0.063 0.47 0.21 – 1.03 0.060 0.46 0.20 – 1.04 0.062 0.45 0.19 – 1.05 0.065 0.49 0.20 – 1.21 0.124 0.46 0.20 – 1.04 0.061 0.53 0.22 – 1.26 0.152
n_hl_threat 1.48 0.75 – 2.93 0.263 1.50 0.65 – 3.48 0.343 1.41 0.60 – 3.35 0.430
n_hl_measure 1.22 0.66 – 2.25 0.534 0.98 0.45 – 2.11 0.951 1.04 0.47 – 2.29 0.925 1.05 0.51 – 2.16 0.885 1.20 0.56 – 2.56 0.636
n_hl_pressure 1.38 0.72 – 2.66 0.335 1.34 0.63 – 2.86 0.445 1.20 0.54 – 2.68 0.657
n_hl_measure *
n_hl_threat
0.82 0.41 – 1.67 0.594
n_hl_measure *
n_hl_pressure
0.68 0.33 – 1.38 0.288
Random Effects
σ2 3.29 3.29 3.29 3.29 3.29 3.29 3.29 3.29
τ00 0.43 featurename 0.46 featurename 0.35 featurename 0.40 featurename 0.47 featurename 0.44 featurename 0.39 featurename 0.42 featurename
ICC 0.12 0.12 0.10 0.11 0.13 0.12 0.11 0.11
N 11 featurename 11 featurename 11 featurename 11 featurename 11 featurename 11 featurename 11 featurename 11 featurename
Observations 50 50 50 50 50 50 50 50
Marginal R2 / Conditional R2 0.000 / 0.116 0.039 / 0.158 0.010 / 0.107 0.027 / 0.134 0.039 / 0.160 0.052 / 0.165 0.028 / 0.130 0.068 / 0.173
AIC 67.774 68.707 69.656 69.078 71.070 73.259 71.424 72.705
compare_performance(m2h_NULL, m2h_threat, m2h_measure, m2h_pressure, m2hA, m2hB, m2hC, m2hD)
## # Comparison of Model Performance Indices
## 
## Name         |    Model |    AIC |    BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
## ------------------------------------------------------------------------------------------------------------------------------------
## m2h_NULL     | glmerMod | 67.518 | 71.343 |      0.116 |      0.000 | 0.116 | 0.441 | 1.000 |    0.571 |    -6.203 |           0.087
## m2h_threat   | glmerMod | 68.185 | 73.921 |      0.158 |      0.039 | 0.123 | 0.433 | 1.000 |    0.555 |    -6.316 |           0.061
## m2h_measure  | glmerMod | 69.134 | 74.870 |      0.107 |      0.010 | 0.097 | 0.443 | 1.000 |    0.576 |    -6.250 |           0.085
## m2h_pressure | glmerMod | 68.556 | 74.292 |      0.134 |      0.027 | 0.109 | 0.437 | 1.000 |    0.565 |    -6.288 |           0.068
## m2hA         | glmerMod | 70.182 | 77.830 |      0.160 |      0.039 | 0.126 | 0.432 | 1.000 |    0.554 |    -6.315 |           0.061
## m2hB         | glmerMod | 71.895 | 81.456 |      0.165 |      0.052 | 0.119 | 0.432 | 1.000 |    0.554 |    -6.330 |           0.052
## m2hC         | glmerMod | 70.535 | 78.183 |      0.130 |      0.028 | 0.105 | 0.438 | 1.000 |    0.567 |    -6.293 |           0.069
## m2hD         | glmerMod | 71.341 | 80.901 |      0.173 |      0.068 | 0.113 | 0.430 | 1.000 |    0.552 |    -6.380 |           0.054

This doesn’t appear to be the case for the long-term trends, no sig. relationship between threats/measures and habitat change

plot_model(m2h_NULL, show.intercept = TRUE)

tab_model(m2h_NULL) 
  increasing
Predictors Odds Ratios CI p
(Intercept) 0.46 0.20 – 1.05 0.065
Random Effects
σ2 3.29
τ00 featurename 0.43
ICC 0.12
N featurename 11
Observations 50
Marginal R2 / Conditional R2 0.000 / 0.116
plot_model(m2h_NULL, type = "pred")
## list()

Notes

  • Missing species? (these aren’t in Article 12 I think )
  • Do we have enough information to weight these species by pop size?
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 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.7.0 tidyr_1.1.2       sjPlot_2.8.7      lme4_1.1-26      
##  [5] Matrix_1.2-18     ggplot2_3.3.3     data.table_1.13.6 expss_0.10.7     
##  [9] dplyr_1.0.7       plyr_1.8.6       
## 
## loaded via a namespace (and not attached):
##  [1] splines_4.0.3      carData_3.0-4      modelr_0.1.8       assertthat_0.2.1  
##  [5] statmod_1.4.35     highr_0.9          cellranger_1.1.0   yaml_2.2.1        
##  [9] bayestestR_0.8.2   pillar_1.6.1       backports_1.2.1    lattice_0.20-41   
## [13] glue_1.4.2         digest_0.6.27      RColorBrewer_1.1-2 ggsignif_0.6.1    
## [17] checkmate_2.0.0    minqa_1.2.4        colorspace_2.0-0   htmltools_0.5.1.1 
## [21] pkgconfig_2.0.3    broom_0.7.5        haven_2.3.1        purrr_0.3.4       
## [25] xtable_1.8-4       mvtnorm_1.1-1      scales_1.1.1       openxlsx_4.2.3    
## [29] rio_0.5.26         emmeans_1.5.4      htmlTable_2.1.0    tibble_3.1.2      
## [33] farver_2.1.0       generics_0.1.0     car_3.0-10         sjlabelled_1.1.7  
## [37] ellipsis_0.3.2     ggpubr_0.4.0       withr_2.4.2        cli_3.0.0         
## [41] readxl_1.3.1       magrittr_2.0.1     crayon_1.4.1       effectsize_0.4.3  
## [45] estimability_1.3   evaluate_0.14      fansi_0.5.0        nlme_3.1-149      
## [49] MASS_7.3-53        rstatix_0.7.0      forcats_0.5.1      foreign_0.8-81    
## [53] tools_4.0.3        hms_1.0.0          lifecycle_1.0.0    matrixStats_0.58.0
## [57] stringr_1.4.0      munsell_0.5.0      zip_2.2.0          ggeffects_1.0.1   
## [61] compiler_4.0.3     rlang_0.4.11       grid_4.0.3         nloptr_1.2.2.2    
## [65] parameters_0.11.0  rstudioapi_0.13    htmlwidgets_1.5.3  labeling_0.4.2    
## [69] rmarkdown_2.6      boot_1.3-25        gtable_0.3.0       abind_1.4-5       
## [73] sjstats_0.18.1     DBI_1.1.1          curl_4.3.2         sjmisc_2.8.6      
## [77] R6_2.5.0           knitr_1.33         utf8_1.2.1         insight_0.12.0    
## [81] stringi_1.6.2      Rcpp_1.0.7         vctrs_0.3.8        tidyselect_1.1.1  
## [85] xfun_0.24          coda_0.19-4