Loading required packages

library(tidyverse)
library(ggridges)
library(ggstance)
library(picante)
packageVersion("picante")
## [1] '1.8.1'

Map (Figure 1)

Input Data

by_country_rich <- read.csv("by_country_richness.csv", sep=";")
mp <- NULL
mapWorld <- borders("world", colour = "#bababa", fill = "#ffffbf")

Figure1 <- by_country_rich %>% 
  ggplot(aes(x = Longitude, y = Latitude, size =richness)) + mapWorld + 
  geom_point(pch=21, fill="#313695", color="black", alpha=0.7) +
  scale_x_continuous(name = "Longitude", breaks = seq(from = -180, to = 180, by = 90), expand = c(0, 0)) + 
  scale_y_continuous(name = "Latitude", limits=c(-60, 90)) +
  scale_size_continuous(breaks=c(1, 25, 50, 100), limits=c(0, 100))+
  ggtitle(label = "Richness of medicinal mammals / country") +
  theme(panel.background = element_rect(fill = "#56B4E950", colour = "white"),
        panel.grid.minor = element_blank(), 
        panel.grid.major = element_line(colour = "#f0f0f0"),
        axis.title = element_text(face = "bold", size = 14),
        axis.text  = element_text(vjust = 0.5, size = 12),
        plot.title = element_text(hjust = 0.5, size = 16, face = "bold"))
Figure1

Figure body mass by IUCN threat category

Data input and handling

data_comp <- read.csv("full_mat_final.csv", h=TRUE)
data_comp <- data_comp[,-1]

#rescalling RI to compensate for species used, but without reference to ICD

larger_zero <- data_comp$IR[data_comp$IR>0]
data_comp$IR[data_comp$IR==0] <- min(larger_zero)/max(larger_zero)

head(data_comp)
##                  species cd1 cd2 cd3 cd4 cd5 cd6 cd7 cd8 cd9 cd10 cd11 cd12
## 1       Acinonyx_jubatus   0   0   0   0   0   0   0   0   0    0    0    0
## 2      Acomys_cineraceus   0   0   0   0   0   0   0   0   0    0    0    0
## 3     Aepyceros_melampus   0   0   0   0   0   0   0   0   0    0    0    0
## 4 Ailuropoda_melanoleuca   1   0   1   1   0   0   0   0   0    0    0    0
## 5  Alcelaphus_buselaphus   0   0   0   0   0   0   0   0   0    0    0    0
## 6            Alces_alces   0   0   0   0   0   0   0   0   0    0    0    0
##   cd13 cd14 cd15 cd16 cd17 cd18 cd19         IR Terrestrial Marine Freshwater
## 1    1    0    0    0    0    0    0 0.03000000           1      0          0
## 2    0    0    0    0    0    0    0 0.01666667           1      0          0
## 3    0    0    0    0    0    0    0 0.01666667           1      0          0
## 4    0    0    0    0    0    0    0 0.31000000           1      0          0
## 5    0    0    0    0    0    0    0 0.01666667           1      0          0
## 6    0    0    1    0    0    0    0 0.03000000           1      0          0
##   Aerial   Mass.g Log_Mass  Island.Endemicity Diet.Plant Diet.Vertebrate
## 1      0  46700.0      4.7 Occurs on mainland          0             100
## 2      0     20.3      1.3 Occurs on mainland          0               0
## 3      0  52500.1      4.7 Occurs on mainland        100               0
## 4      0 108400.0      5.0 Occurs on mainland         90              10
## 5      0 171001.5      5.2 Occurs on mainland        100               0
## 6      0 356998.0      5.6 Occurs on mainland        100               0
##   Diet.Invertebrate IUCN comu       range
## 1                 0   VU    1  441.849278
## 2               100   LC    1   86.780792
## 3                 0   LC    1  293.986996
## 4                 0   VU    1    3.290258
## 5                 0   LC    1  576.008919
## 6                 0   LC    1 3606.410850
dim(data_comp)#565 spp
## [1] 565  34

Importing phylogeny

The tree is available at the VertLife webpage

full_nexus <- read.nexus("https://github.com/n8upham/MamPhy_v1/blob/master/_DATA/MamPhy_fullPosterior_BDvr_DNAonly_4098sp_topoFree_NDexp_MCC_v2_target.tre?raw=TRUE")

full_nexus_semAnolis <- drop.tip(full_nexus, "_Anolis_carolinensis")#remove Anolis

sp_names_full <- full_nexus_semAnolis$tip.label
split_names <- sapply(strsplit(sp_names_full, split='_'), function(x) (c(x[1], x[2])))
split_names_t <- data.frame(t(split_names))
new_names <- paste(split_names_t$X1, split_names_t$X2, sep="_")

# check if it is correct

data.frame(original=full_nexus$tip.label[-1][1:10], fixed=new_names[1:10])
##                                                     original
## 1     Ornithorhynchus_anatinus_ORNITHORHYNCHIDAE_MONOTREMATA
## 2              Zaglossus_bruijnii_TACHYGLOSSIDAE_MONOTREMATA
## 3          Tachyglossus_aculeatus_TACHYGLOSSIDAE_MONOTREMATA
## 4          Rhynchocyon_petersi_MACROSCELIDIDAE_MACROSCELIDEA
## 5           Rhynchocyon_cirnei_MACROSCELIDIDAE_MACROSCELIDEA
## 6     Rhynchocyon_udzungwensis_MACROSCELIDIDAE_MACROSCELIDEA
## 7       Elephantulus_rufescens_MACROSCELIDIDAE_MACROSCELIDEA
## 8  Elephantulus_brachyrhynchus_MACROSCELIDIDAE_MACROSCELIDEA
## 9          Elephantulus_intufi_MACROSCELIDIDAE_MACROSCELIDEA
## 10      Elephantulus_rupestris_MACROSCELIDIDAE_MACROSCELIDEA
##                          fixed
## 1     Ornithorhynchus_anatinus
## 2           Zaglossus_bruijnii
## 3       Tachyglossus_aculeatus
## 4          Rhynchocyon_petersi
## 5           Rhynchocyon_cirnei
## 6     Rhynchocyon_udzungwensis
## 7       Elephantulus_rufescens
## 8  Elephantulus_brachyrhynchus
## 9          Elephantulus_intufi
## 10      Elephantulus_rupestris
full_nexus_semAnolis$tip.label <- new_names
lista <- read.table("lista_may.txt", h=TRUE)#list of species
head(lista);dim(lista)#565 spp
##                        comu
## Acinonyx_jubatus          1
## Acomys_cineraceus         1
## Aepyceros_melampus        1
## Ailuropoda_melanoleuca    1
## Alcelaphus_buselaphus     1
## Alces_alces               1
## [1] 565   1
phylo_pruna <- prune.sample(t(lista), full_nexus_semAnolis)
phylo_pruna#521 spp
## 
## Phylogenetic tree with 521 tips and 520 internal nodes.
## 
## Tip labels:
##  Dugong_dugon, Trichechus_manatus, Trichechus_inunguis, Trichechus_senegalensis, Procavia_capensis, Dendrohyrax_dorsalis, ...
## 
## Rooted; includes branch lengths.
data_comp_row <- column_to_rownames(data_comp, var = "species")
head(data_comp_row)
##                        cd1 cd2 cd3 cd4 cd5 cd6 cd7 cd8 cd9 cd10 cd11 cd12 cd13
## Acinonyx_jubatus         0   0   0   0   0   0   0   0   0    0    0    0    1
## Acomys_cineraceus        0   0   0   0   0   0   0   0   0    0    0    0    0
## Aepyceros_melampus       0   0   0   0   0   0   0   0   0    0    0    0    0
## Ailuropoda_melanoleuca   1   0   1   1   0   0   0   0   0    0    0    0    0
## Alcelaphus_buselaphus    0   0   0   0   0   0   0   0   0    0    0    0    0
## Alces_alces              0   0   0   0   0   0   0   0   0    0    0    0    0
##                        cd14 cd15 cd16 cd17 cd18 cd19         IR Terrestrial
## Acinonyx_jubatus          0    0    0    0    0    0 0.03000000           1
## Acomys_cineraceus         0    0    0    0    0    0 0.01666667           1
## Aepyceros_melampus        0    0    0    0    0    0 0.01666667           1
## Ailuropoda_melanoleuca    0    0    0    0    0    0 0.31000000           1
## Alcelaphus_buselaphus     0    0    0    0    0    0 0.01666667           1
## Alces_alces               0    1    0    0    0    0 0.03000000           1
##                        Marine Freshwater Aerial   Mass.g Log_Mass
## Acinonyx_jubatus            0          0      0  46700.0      4.7
## Acomys_cineraceus           0          0      0     20.3      1.3
## Aepyceros_melampus          0          0      0  52500.1      4.7
## Ailuropoda_melanoleuca      0          0      0 108400.0      5.0
## Alcelaphus_buselaphus       0          0      0 171001.5      5.2
## Alces_alces                 0          0      0 356998.0      5.6
##                         Island.Endemicity Diet.Plant Diet.Vertebrate
## Acinonyx_jubatus       Occurs on mainland          0             100
## Acomys_cineraceus      Occurs on mainland          0               0
## Aepyceros_melampus     Occurs on mainland        100               0
## Ailuropoda_melanoleuca Occurs on mainland         90              10
## Alcelaphus_buselaphus  Occurs on mainland        100               0
## Alces_alces            Occurs on mainland        100               0
##                        Diet.Invertebrate IUCN comu       range
## Acinonyx_jubatus                       0   VU    1  441.849278
## Acomys_cineraceus                    100   LC    1   86.780792
## Aepyceros_melampus                     0   LC    1  293.986996
## Ailuropoda_melanoleuca                 0   VU    1    3.290258
## Alcelaphus_buselaphus                  0   LC    1  576.008919
## Alces_alces                            0   LC    1 3606.410850
tidydata <- match.phylo.data(phylo_pruna, data_comp_row)
## [1] "Dropping taxa from the data because they are not present in the phylogeny:"
##  [1] "Acomys_cineraceus"       "Anomalurus_derbianus"   
##  [3] "Cabassous_centralis"     "Colobus_vellerosus"     
##  [5] "Conepatus_humboldtii"    "Crocidura_nigeriae"     
##  [7] "Dasyprocta_prymnolopha"  "Dasyprocta_variegata"   
##  [9] "Dendrohyrax_arboreus"    "Euroscaptor_micrura"    
## [11] "Funisciurus_anerythrus"  "Funisciurus_leucogenys" 
## [13] "Funisciurus_substriatus" "Graphiurus_lorraineus"  
## [15] "Heliosciurus_gambianus"  "Hemiechinus_collaris"   
## [17] "Hesperoptenus_tickelli"  "Pipistrellus_imbricatus"
## [19] "Kerivoula_africana"      "Leopardus_geoffroyi"    
## [21] "Lepus_fagani"            "Lepus_nigricollis"      
## [23] "Lycalopex_culpaeus"      "Lycalopex_gymnocercus"  
## [25] "Lycalopex_sechurae"      "Manis_crassicaudata"    
## [27] "Mazama_chunyi"           "Microtus_mexicanus"     
## [29] "Murina_fusca"            "Paraechinus_micropus"   
## [31] "Piliocolobus_preussi"    "Plecotus_ariel"         
## [33] "Plecotus_sacrimontis"    "Presbytis_hosei"        
## [35] "Presbytis_natunae"       "Presbytis_siamensis"    
## [37] "Pteropus_faunulus"       "Saguinus_fuscicollis"   
## [39] "Steatomys_jacksoni"      "Tadarida_latouchei"     
## [41] "Tamandua_mexicana"       "Tolypeutes_tricinctus"  
## [43] "Viverra_civettina"       "Vulpes_bengalensis"
tidydata$phy#521 spp
## 
## Phylogenetic tree with 521 tips and 520 internal nodes.
## 
## Tip labels:
##  Dugong_dugon, Trichechus_manatus, Trichechus_inunguis, Trichechus_senegalensis, Procavia_capensis, Dendrohyrax_dorsalis, ...
## 
## Rooted; includes branch lengths.
head(tidydata$data)
##                         cd1 cd2 cd3 cd4 cd5 cd6 cd7 cd8 cd9 cd10 cd11 cd12 cd13
## Dugong_dugon              1   0   0   0   1   0   0   0   0    0    0    1    0
## Trichechus_manatus        0   0   1   1   1   1   0   1   0    0    1    0    0
## Trichechus_inunguis       1   0   1   1   1   1   1   1   0    0    1    0    0
## Trichechus_senegalensis   1   0   0   0   0   0   0   0   0    0    0    0    0
## Procavia_capensis         0   1   0   0   0   0   0   0   0    0    0    0    0
## Dendrohyrax_dorsalis      0   0   0   0   0   0   0   0   0    0    0    1    0
##                         cd14 cd15 cd16 cd17 cd18 cd19         IR Terrestrial
## Dugong_dugon               0    0    0    0    0    0 0.29000000           0
## Trichechus_manatus         0    0    0    0    1    0 0.83000000           0
## Trichechus_inunguis        0    0    0    0    1    0 1.18000000           0
## Trichechus_senegalensis    0    0    0    0    0    0 0.03000000           0
## Procavia_capensis          0    0    0    0    0    0 0.03000000           1
## Dendrohyrax_dorsalis       0    0    0    0    0    0 0.06000000           1
##                         Marine Freshwater Aerial     Mass.g Log_Mass
## Dugong_dugon                 1          0      0   310000.0      5.5
## Trichechus_manatus           1          1      0   341000.0      5.5
## Trichechus_inunguis          1          1      0   480000.0      5.7
## Trichechus_senegalensis      1          1      0   454000.0      5.7
## Procavia_capensis            0          0      0     3030.0      3.5
## Dendrohyrax_dorsalis         0          0      0     3175.0      3.5
##                          Island.Endemicity Diet.Plant Diet.Vertebrate
## Dugong_dugon            Exclusively marine        100               0
## Trichechus_manatus      Occurs on mainland        100               0
## Trichechus_inunguis     Occurs on mainland        100               0
## Trichechus_senegalensis Occurs on mainland         95               0
## Procavia_capensis       Occurs on mainland        100               0
## Dendrohyrax_dorsalis    Occurs on mainland        100               0
##                         Diet.Invertebrate IUCN comu        range
## Dugong_dugon                            0   VU    1 5.380811e+02
## Trichechus_manatus                      0   VU <NA> 1.302891e+02
## Trichechus_inunguis                     0   VU <NA> 5.178365e+01
## Trichechus_senegalensis                 5   VU <NA> 3.829432e+01
## Procavia_capensis                       0   LC    1 1.328971e+03
## Dendrohyrax_dorsalis                    0   LC    1 2.489545e+02
tidydata$data %>% filter(IUCN=="LC" | IUCN=="EN"  | IUCN=="CR"  | IUCN=="VU"  | IUCN=="NT" ) -> iucn_dat

iucn_dat$IUCN
##   [1] "VU" "VU" "VU" "VU" "LC" "LC" "VU" "EN" "LC" "LC" "LC" "LC" "LC" "VU" "LC"
##  [16] "NT" "LC" "LC" "VU" "LC" "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
##  [31] "LC" "LC" "LC" "LC" "CR" "CR" "VU" "VU" "VU" "VU" "LC" "LC" "LC" "EN" "LC"
##  [46] "LC" "LC" "NT" "LC" "LC" "NT" "EN" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "VU"
##  [61] "VU" "EN" "LC" "LC" "LC" "VU" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
##  [76] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "VU" "LC" "LC" "NT" "LC"
##  [91] "EN" "EN" "NT" "VU" "VU" "NT" "NT" "LC" "LC" "VU" "VU" "VU" "VU" "LC" "VU"
## [106] "LC" "VU" "VU" "LC" "NT" "NT" "LC" "LC" "VU" "LC" "LC" "LC" "LC" "LC" "LC"
## [121] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "VU" "VU" "EN" "LC"
## [136] "LC" "LC" "LC" "LC" "LC" "VU" "NT" "NT" "EN" "LC" "NT" "NT" "VU" "VU" "EN"
## [151] "LC" "NT" "LC" "LC" "LC" "LC" "EN" "VU" "VU" "NT" "VU" "NT" "VU" "NT" "EN"
## [166] "EN" "CR" "NT" "CR" "CR" "VU" "EN" "EN" "VU" "EN" "LC" "LC" "CR" "VU" "LC"
## [181] "LC" "LC" "LC" "LC" "LC" "LC" "VU" "VU" "EN" "CR" "EN" "LC" "EN" "NT" "LC"
## [196] "LC" "LC" "VU" "EN" "LC" "LC" "VU" "LC" "CR" "LC" "VU" "EN" "LC" "EN" "LC"
## [211] "LC" "VU" "VU" "LC" "VU" "LC" "VU" "NT" "VU" "LC" "VU" "LC" "EN" "EN" "EN"
## [226] "VU" "LC" "LC" "LC" "LC" "LC" "LC" "NT" "EN" "CR" "LC" "EN" "CR" "LC" "LC"
## [241] "LC" "VU" "VU" "LC" "LC" "VU" "CR" "CR" "NT" "LC" "LC" "LC" "LC" "LC" "LC"
## [256] "LC" "NT" "LC" "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "NT" "NT" "VU"
## [271] "VU" "LC" "VU" "NT" "LC" "LC" "VU" "VU" "VU" "NT" "LC" "EN" "NT" "LC" "LC"
## [286] "NT" "EN" "LC" "LC" "EN" "LC" "LC" "LC" "NT" "LC" "LC" "LC" "LC" "LC" "LC"
## [301] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "EN" "NT" "VU" "LC"
## [316] "EN" "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [331] "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "NT" "LC" "LC"
## [346] "LC" "LC" "EN" "LC" "VU" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [361] "LC" "LC" "LC" "VU" "VU" "VU" "LC" "LC" "LC" "LC" "LC" "LC" "VU" "LC" "LC"
## [376] "LC" "EN" "VU" "EN" "LC" "VU" "LC" "LC" "LC" "CR" "EN" "CR" "EN" "EN" "CR"
## [391] "EN" "EN" "CR" "LC" "VU" "VU" "EN" "VU" "CR" "CR" "EN" "NT" "NT" "VU" "LC"
## [406] "CR" "EN" "EN" "CR" "EN" "EN" "NT" "EN" "VU" "EN" "VU" "VU" "NT" "LC" "LC"
## [421] "LC" "VU" "LC" "LC" "EN" "VU" "LC" "LC" "VU" "VU" "LC" "LC" "LC" "LC" "LC"
## [436] "LC" "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [451] "LC" "LC" "LC" "LC" "LC" "NT" "NT" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [466] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [481] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
## [496] "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC" "LC"
iucn_dat$IUCN <- factor(iucn_dat$IUCN, levels = c("LC", "NT", "VU", "EN", "CR"))

cd_dat <- iucn_dat %>% 
        dplyr::select(paste("cd", 1:19, sep="")) 

cd_dat2 <- mutate_all(cd_dat, function(x) as.numeric(as.character(x)))
rowSums(cd_dat2) # check those species without a specific use
##   [1]  3  7  9  1  1  1  8  7  0  1  4  2  7  6  1  2  1  8  3  1  3 14  2  2  5
##  [26]  0  2  0  1  1  0  0  0  0  4  7  1  4  5 12  3  0  0  0  1  4  6  0  4  2
##  [51]  5  1  8  1  0  1  0  0  0  0  2  1  2 NA  2  7  1 15  8  6  3  2  3  8 11
##  [76]  4  5  1  0  1  0  0  1  1  3  3  0  0  1  0  0  0  2  1  2  0  0  1  6  3
## [101]  6  4  8  6  6 10  5  3  1  6  1  5  0  0  0  1  0  0  0  0  2  1  1  0  2
## [126]  1  3  2  0  0  1  1  0  0  0  1  0  0  9  3  1  0  0  0  0  0  1  0  6  2
## [151]  4  1  1  4  0  2  8  7  4  6  3  0  1  1  1  0  3  5  6  6  7  2  0  6  2
## [176]  1  8  9  5  4 14  1  0  6  0  0  2 10  2  0 11  0  2 10  1  0  4  1  0  0
## [201]  2  1  1  0  3  1  2  2  0  5  8  6  0  0  7  6  1  4  5  1  0  2  6  6  1
## [226] 13  0  2  2  0  0  0  0  0  0  2  0  0  0  0  0  0  2  1  1  0  0  3  0  0
## [251]  0  1  0  2  0  0  1  0  0  0  1  0  3  1  0  0  0  3  1  0  0  0  0  3  0
## [276]  1  0  0  2  2  0  9  2  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
## [301]  0  0  0  0  0  0  0  0  5  4  1  0  1  0  4  0  0  1  0  0  0  0  0  0  1
## [326]  2  1  1  0  0  0  0  0  0  3  5  1  2  0  0  0  0  0  0  0  0  1  0  0  0
## [351]  0  0  0  0  0  0  0  0  0  0  0  1  3  0  1  0  0  2  0  0  1  1  0  5  1
## [376]  1  5  1  1  1  3  1  0  2  0  0  0  2  0  0  3  0  2  0  0  1  0  0  0  0
## [401]  0  0  0  1  2  0  2  0  0  0  0  1  1  0  0  0  0  4  1  0  6  0  2  2  1
## [426]  1  3  0  0  0  0  2  1  2  0  7  6  0  1  3  1  0  3  0  3  0  0  0  0  1
## [451]  3  3  0  0  1  2  2  1  2  1  4  7  4  2  0  3  9  6 12  1  4  9  3  2  1
## [476]  0  1  1  1  0  1  0  1  0  0  0  0  0  0  0  0  0  5  0  0  0  2  0  5  3
## [501]  6  0  4  4  3  7  9  3  5
iucn_dat$use <- rowSums(cd_dat2)

iucn_dat %>% 
        filter(use>0) -> iucn_dat2

names(iucn_dat2)
##  [1] "cd1"               "cd2"               "cd3"              
##  [4] "cd4"               "cd5"               "cd6"              
##  [7] "cd7"               "cd8"               "cd9"              
## [10] "cd10"              "cd11"              "cd12"             
## [13] "cd13"              "cd14"              "cd15"             
## [16] "cd16"              "cd17"              "cd18"             
## [19] "cd19"              "IR"                "Terrestrial"      
## [22] "Marine"            "Freshwater"        "Aerial"           
## [25] "Mass.g"            "Log_Mass"          "Island.Endemicity"
## [28] "Diet.Plant"        "Diet.Vertebrate"   "Diet.Invertebrate"
## [31] "IUCN"              "comu"              "range"            
## [34] "use"

Figure

iucn_dat2 %>% 
        ggplot(aes(x = as.numeric(as.character(Log_Mass)), 
                   y = IUCN, group = IUCN, fill=IUCN)) +
        theme_bw() +
        geom_density_ridges() +
        geom_boxploth(width = 0.1) +
        scale_fill_viridis_d(alpha = 0.5) +
        xlab("log Body mass (g)") +
        ylab("IUCN red list categories")+
        theme(panel.grid.minor = element_blank(),
              legend.position = "none",
              axis.title.x = element_text(size=14),
              axis.text.x  = element_text(size=12),
              axis.title.y = element_text(size=14),
              axis.text.y  = element_text(size=12))
## Picking joint bandwidth of 0.409