Mammal data - Mediterranean Maquis
Talia Gabay
2023-11-14
P <- P_byplot[grepl("Mediterranean Maquis",Unit),][order(Cycle_number,Site)][,Site:=factor(Site)]Mediterranean Maquis - Richness by camera
Monitoring started in 2014 (5 monitoring cycles), Each site has two types of plots - far or near settlements
There are 15 sites total, each site has 5 sampling points
Let’s first validate the data:
#include rare species in analysis
P.anal <- copy(P) # set a fixed variable name for analysis, if want to switch between data WITH rare #species and data WITHOUT rare species then only change once here
#Validate factors and levels
P.anal$Settlements <- as.factor(P.anal$Settlements)
print("Settlements has 2 levels")## [1] "Settlements has 2 levels"print(levels(P.anal$Settlements))## [1] "Far" "Near"P.anal$Subunit <- as.factor(P.anal$Subunit)
print("Subunit has 3 levels")## [1] "Subunit has 3 levels"print(levels(P.anal$Subunit))## [1] "Carmel" "Galilee" "Judea"P.anal$Site <- as.factor(P.anal$Site)
print("Site has 15 levels")## [1] "Site has 15 levels"print(levels(P.anal$Site))## [1] "Aderet" "Beit Oren" "Ein Yaakov" "Givat Yearim"
## [5] "Givat Yeshayahu" "Goren" "Iftach" "Kerem Maharal"
## [9] "Kfar Shamai" "Kisalon" "Margaliot" "Nehusha"
## [13] "Nir Etzion" "Ofer" "Yagur"P.anal$Transect <- as.factor(P.anal$Transect)
print("Transect has 30 levels")## [1] "Transect has 30 levels"print(levels(P.anal$Transect))## [1] "Aderet Far" "Aderet Near" "Beit Oren Far"
## [4] "Beit Oren Near" "Ein Yaakov Far" "Ein Yaakov Near"
## [7] "Givat Yearim Far" "Givat Yearim Near" "Givat Yeshayahu Far"
## [10] "Givat Yeshayahu Near" "Goren Far" "Goren Near"
## [13] "Iftach Far" "Iftach Near" "Kerem Maharal Far"
## [16] "Kerem Maharal Near" "Kfar Shamai Far" "Kfar Shamai Near"
## [19] "Kisalon Far" "Kisalon Near" "Margaliot Far"
## [22] "Margaliot Near" "Nehusha Far" "Nehusha Near"
## [25] "Nir Etzion Far" "Nir Etzion Near" "Ofer Far"
## [28] "Ofer Near" "Yagur Far" "Yagur Near"P.anal$Transect_with_date <- as.factor(P.anal$Transect_with_date)
print("Transect with date has 150 levels")## [1] "Transect with date has 150 levels"print(levels(P.anal$Transect_with_date))## [1] "Aderet Far_03_10_2017" "Aderet Far_10_08_2021"
## [3] "Aderet Far_10_10_2015" "Aderet Far_17_07_2012"
## [5] "Aderet Far_29_09_2019" "Aderet Near_03_10_2017"
## [7] "Aderet Near_10_08_2021" "Aderet Near_17_07_2012"
## [9] "Aderet Near_22_10_2015" "Aderet Near_29_09_2019"
## [11] "Beit Oren Far_03_06_2015" "Beit Oren Far_05_11_2019"
## [13] "Beit Oren Far_07_07_2017" "Beit Oren Far_07_08_2021"
## [15] "Beit Oren Far_18_06_2012" "Beit Oren Near_05_11_2019"
## [17] "Beit Oren Near_07_07_2017" "Beit Oren Near_07_08_2021"
## [19] "Beit Oren Near_18_06_2012" "Beit Oren Near_25_06_2015"
## [21] "Ein Yaakov Far_03_06_2012" "Ein Yaakov Far_06_02_2020"
## [23] "Ein Yaakov Far_06_07_2021" "Ein Yaakov Far_11_11_2015"
## [25] "Ein Yaakov Far_20_10_2017" "Ein Yaakov Near_03_06_2012"
## [27] "Ein Yaakov Near_06_02_2020" "Ein Yaakov Near_06_07_2021"
## [29] "Ein Yaakov Near_10_11_2015" "Ein Yaakov Near_20_10_2017"
## [31] "Givat Yearim Far_03_07_2015" "Givat Yearim Far_11_07_2021"
## [33] "Givat Yearim Far_15_09_2019" "Givat Yearim Far_29_07_2017"
## [35] "Givat Yearim Far_31_07_2012" "Givat Yearim Near_11_07_2021"
## [37] "Givat Yearim Near_16_09_2019" "Givat Yearim Near_21_07_2015"
## [39] "Givat Yearim Near_29_07_2017" "Givat Yearim Near_31_07_2012"
## [41] "Givat Yeshayahu Far_03_09_2019" "Givat Yeshayahu Far_16_07_2012"
## [43] "Givat Yeshayahu Far_24_08_2017" "Givat Yeshayahu Far_28_06_2021"
## [45] "Givat Yeshayahu Far_31_08_2015" "Givat Yeshayahu Near_03_09_2019"
## [47] "Givat Yeshayahu Near_16_07_2012" "Givat Yeshayahu Near_24_08_2017"
## [49] "Givat Yeshayahu Near_28_06_2021" "Givat Yeshayahu Near_29_08_2015"
## [51] "Goren Far_03_09_2015" "Goren Far_14_01_2020"
## [53] "Goren Far_17_08_2021" "Goren Far_20_10_2017"
## [55] "Goren Far_22_05_2012" "Goren Near_03_09_2015"
## [57] "Goren Near_14_01_2020" "Goren Near_17_08_2021"
## [59] "Goren Near_20_10_2017" "Goren Near_22_05_2012"
## [61] "Iftach Far_04_11_2017" "Iftach Far_05_08_2015"
## [63] "Iftach Far_09_09_2021" "Iftach Far_16_04_2020"
## [65] "Iftach Far_21_05_2012" "Iftach Near_04_11_2017"
## [67] "Iftach Near_05_08_2015" "Iftach Near_09_09_2021"
## [69] "Iftach Near_16_04_2020" "Iftach Near_21_05_2012"
## [71] "Kerem Maharal Far_01_07_2012" "Kerem Maharal Far_07_12_2019"
## [73] "Kerem Maharal Far_11_08_2017" "Kerem Maharal Far_19_05_2015"
## [75] "Kerem Maharal Far_24_07_2021" "Kerem Maharal Near_01_07_2012"
## [77] "Kerem Maharal Near_07_12_2019" "Kerem Maharal Near_11_08_2017"
## [79] "Kerem Maharal Near_19_05_2015" "Kerem Maharal Near_23_07_2021"
## [81] "Kfar Shamai Far_09_11_2021" "Kfar Shamai Far_13_07_2015"
## [83] "Kfar Shamai Far_24_05_2020" "Kfar Shamai Far_28_09_2017"
## [85] "Kfar Shamai Far_29_04_2012" "Kfar Shamai Near_09_11_2021"
## [87] "Kfar Shamai Near_13_07_2015" "Kfar Shamai Near_24_05_2020"
## [89] "Kfar Shamai Near_28_09_2017" "Kfar Shamai Near_29_04_2012"
## [91] "Kisalon Far_19_07_2015" "Kisalon Far_19_11_2021"
## [93] "Kisalon Far_21_06_2017" "Kisalon Far_22_10_2019"
## [95] "Kisalon Far_30_07_2012" "Kisalon Near_01_07_2015"
## [97] "Kisalon Near_19_11_2021" "Kisalon Near_21_06_2017"
## [99] "Kisalon Near_22_10_2019" "Kisalon Near_30_07_2012"
## [101] "Margaliot Far_04_11_2017" "Margaliot Far_16_03_2020"
## [103] "Margaliot Far_18_01_2022" "Margaliot Far_20_09_2015"
## [105] "Margaliot Far_30_04_2012" "Margaliot Near_04_11_2017"
## [107] "Margaliot Near_16_03_2020" "Margaliot Near_18_01_2022"
## [109] "Margaliot Near_28_10_2015" "Margaliot Near_30_04_2012"
## [111] "Nehusha Far_02_10_2019" "Nehusha Far_05_09_2017"
## [113] "Nehusha Far_13_08_2012" "Nehusha Far_13_09_2015"
## [115] "Nehusha Far_23_08_2021" "Nehusha Near_02_10_2019"
## [117] "Nehusha Near_05_09_2017" "Nehusha Near_13_08_2012"
## [119] "Nehusha Near_23_08_2021" "Nehusha Near_26_09_2015"
## [121] "Nir Etzion Far_02_06_2015" "Nir Etzion Far_16_06_2017"
## [123] "Nir Etzion Far_17_06_2012" "Nir Etzion Far_25_11_2019"
## [125] "Nir Etzion Far_26_06_2021" "Nir Etzion Near_02_06_2015"
## [127] "Nir Etzion Near_16_06_2017" "Nir Etzion Near_17_06_2012"
## [129] "Nir Etzion Near_25_11_2019" "Nir Etzion Near_26_06_2021"
## [131] "Ofer Far_02_07_2012" "Ofer Far_04_01_2022"
## [133] "Ofer Far_13_05_2017" "Ofer Far_16_12_2019"
## [135] "Ofer Far_18_05_2015" "Ofer Near_02_07_2012"
## [137] "Ofer Near_04_01_2022" "Ofer Near_13_05_2017"
## [139] "Ofer Near_16_12_2019" "Ofer Near_18_05_2015"
## [141] "Yagur Far_04_06_2012" "Yagur Far_10_06_2021"
## [143] "Yagur Far_23_12_2019" "Yagur Far_28_06_2015"
## [145] "Yagur Far_28_08_2017" "Yagur Near_04_06_2012"
## [147] "Yagur Near_10_06_2021" "Yagur Near_23_12_2019"
## [149] "Yagur Near_28_08_2017" "Yagur Near_29_10_2015"print("RICHNESS WITH RARE SPECIES")## [1] "RICHNESS WITH RARE SPECIES"plot_alpha_diversity(P, x_val = "Subunit", y_val = "richness", ylab_val = "richness", xlab_val = "Subunit", fill_val = "Subunit")
col_names <- c("Cycle_number", "Deployment.id_new", "Unit", "Subunit", "Site", "Settlements", "Distance_rescaled", "Lon", "Lat", "year", "rescaled_Time.Diff", "sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "Transect", "Transect_with_date")
abu_by_spp.Maquis <- all.mammal.abundances[grepl("Mediterranean Maquis",Unit),]
spp <- abu_by_spp.Maquis[,c(27:39)]
# filter out species with zero counts
print(colSums(spp))## Canis aureus Canis lupus Capra nubiana Equus hemionus Gazella dorcas
## 6971 44 0 0 0
## Gazella gazella Hyaena hyaena Hystrix indica Lepus capensis Meles meles
## 406 298 3019 10 408
## Oryx leucoryx Sus scrofa Vulpes vulpes
## 0 10124 1478spp <- spp[,.SD,.SDcols = colSums(spp)>0]
spp <- mvabund(spp)
plot(spp ~ P.anal$Subunit, overall.main="raw abundances", transformation = "no")## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR
my_cols <- c("blue", "darkgreen", "red")
dotchart(P.anal$richness, ylab = "Subunit",
groups = P.anal$Subunit, gcolor = my_cols,
color = my_cols[P.anal$Subunit],
cex = 0.9, pch = 1, xlab = "richness")
kable(summary(P.anal[,.(richness, abundance, Subunit, rescaled_Time.Diff, cosinus_Monitoring.Time.Diff, sinus_Monitoring.Time.Diff, Site, Transect, Transect_with_date, Distance_rescaled)]))| richness | abundance | Subunit | rescaled_Time.Diff | cosinus_Monitoring.Time.Diff | sinus_Monitoring.Time.Diff | Site | Transect | Transect_with_date | Distance_rescaled | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. :0.000 | Min. : 0.00 | Carmel :376 | Min. :-1.8927 | Min. :-0.99996 | Min. :-0.999990 | Goren : 87 | Goren Near : 45 | Kerem Maharal Far_19_05_2015: 10 | Min. :-0.9136 | |
| 1st Qu.:2.000 | 1st Qu.: 4.00 | Galilee:417 | 1st Qu.:-0.9172 | 1st Qu.:-0.83907 | 1st Qu.:-0.602000 | Ein Yaakov : 85 | Ein Yaakov Near : 44 | Aderet Near_10_08_2021 : 9 | 1st Qu.:-0.8371 | |
| Median :2.000 | Median : 10.00 | Judea :324 | Median :-0.1836 | Median :-0.04866 | Median : 0.000000 | Kfar Shamai: 84 | Goren Far : 42 | Beit Oren Far_07_07_2017 : 9 | Median :-0.6362 | |
| Mean :2.569 | Mean : 20.37 | NA | Mean :-0.2623 | Mean :-0.08610 | Mean : 0.008294 | Margaliot : 82 | Kfar Shamai Far : 42 | Ein Yaakov Far_06_07_2021 : 9 | Mean :-0.5033 | |
| 3rd Qu.:3.000 | 3rd Qu.: 24.00 | NA | 3rd Qu.: 0.5940 | 3rd Qu.: 0.56975 | 3rd Qu.: 0.733190 | Iftach : 79 | Kfar Shamai Near: 42 | Ein Yaakov Far_11_11_2015 : 9 | 3rd Qu.:-0.3237 | |
| Max. :7.000 | Max. :341.00 | NA | Max. : 1.1755 | Max. : 1.00000 | Max. : 0.999990 | Nir Etzion : 77 | Margaliot Near : 42 | Ein Yaakov Near_03_06_2012 : 9 | Max. : 1.6119 | |
| NA | NA | NA | NA | NA | NA | (Other) :623 | (Other) :860 | (Other) :1062 | NA |
pairs(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("Subunit","sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "rescaled_Time.Diff","Site","Transect", "Transect_with_date", "Distance_rescaled")])
kable(cor(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("Subunit","sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "rescaled_Time.Diff","Site","Transect", "Transect_with_date", "Distance_rescaled")]))| Subunit | sinus_Monitoring.Time.Diff | cosinus_Monitoring.Time.Diff | rescaled_Time.Diff | Site | Transect | Transect_with_date | Distance_rescaled | |
|---|---|---|---|---|---|---|---|---|
| Subunit | 1.0000000 | 0.0808235 | 0.1250776 | 0.0324596 | -0.3922386 | -0.3901787 | -0.3895513 | -0.1278084 |
| sinus_Monitoring.Time.Diff | 0.0808235 | 1.0000000 | 0.0926611 | 0.0071447 | 0.0360044 | 0.0345840 | 0.0366821 | -0.0592206 |
| cosinus_Monitoring.Time.Diff | 0.1250776 | 0.0926611 | 1.0000000 | 0.0589309 | -0.0782191 | -0.0774880 | -0.0767999 | -0.0119164 |
| rescaled_Time.Diff | 0.0324596 | 0.0071447 | 0.0589309 | 1.0000000 | -0.0128696 | -0.0125757 | -0.0165851 | -0.0556040 |
| Site | -0.3922386 | 0.0360044 | -0.0782191 | -0.0128696 | 1.0000000 | 0.9983334 | 0.9977772 | 0.0211274 |
| Transect | -0.3901787 | 0.0345840 | -0.0774880 | -0.0125757 | 0.9983334 | 1.0000000 | 0.9994624 | -0.0167655 |
| Transect_with_date | -0.3895513 | 0.0366821 | -0.0767999 | -0.0165851 | 0.9977772 | 0.9994624 | 1.0000000 | -0.0171086 |
| Distance_rescaled | -0.1278084 | -0.0592206 | -0.0119164 | -0.0556040 | 0.0211274 | -0.0167655 | -0.0171086 | 1.0000000 |
meanvar.plot(spp, xlab = "mean abundance of a given species across sites", ylab = "variance of the abundance of a given species across sites")
Community analysis using package MVabund
Table of abundances per species by camera ID
abundance_and_cameras <- as.data.table(abu_by_spp.Maquis[,c(5,27:39)])
dim(abundance_and_cameras)## [1] 1117 14canisau <- abundance_and_cameras[,2] > 0
sum(canisau)## [1] 678canislup <- abundance_and_cameras[,3] > 0
sum(canislup)## [1] 20capra <- abundance_and_cameras[,4] > 0
sum(capra)## [1] 0Equus <- abundance_and_cameras[,5] > 0
sum(Equus)## [1] 0Gazellador <- abundance_and_cameras[,6] > 0
sum(Gazellador)## [1] 0Gazellagaz <- abundance_and_cameras[,7] > 0
sum(Gazellagaz)## [1] 141Hyaena <- abundance_and_cameras[,8] > 0
sum(Hyaena)## [1] 153Hystrix <- abundance_and_cameras[,9] > 0
sum(Hystrix)## [1] 599Lepus <- abundance_and_cameras[,10] > 0
sum(Lepus)## [1] 4Meles <- abundance_and_cameras[,11] > 0
sum(Meles)## [1] 178Oryx <- abundance_and_cameras[,12] > 0
sum(Oryx)## [1] 0Sus <- abundance_and_cameras[,13] > 0
sum(Sus)## [1] 703Vulpus <- abundance_and_cameras[,14] > 0
sum(Vulpus)## [1] 3941117 cameras. Species with a minimum of 20 cameras - Canis aureus, Canis lupus, Gazella gazella, Hyaena hyaena, Hystrix indica, Meles meles, Sus scorfa and Vulpus vulpus
spp_no_rare <- abu_by_spp.Maquis[,27:39]
print(colSums(spp_no_rare))## Canis aureus Canis lupus Capra nubiana Equus hemionus Gazella dorcas
## 6971 44 0 0 0
## Gazella gazella Hyaena hyaena Hystrix indica Lepus capensis Meles meles
## 406 298 3019 10 408
## Oryx leucoryx Sus scrofa Vulpes vulpes
## 0 10124 1478# filter out species with less than 20 individuals (rare species)
spp_no_rare <- spp_no_rare[,.SD,.SDcols = colSums(spp_no_rare)>=20]
#spp_no_rare <- spp_no_rare[,c(-2,-6)]
spp_no_rare <- mvabund(spp_no_rare)
print(colSums(spp_no_rare))## Canis.aureus Canis.lupus Gazella.gazella Hyaena.hyaena Hystrix.indica
## 6971 44 406 298 3019
## Meles.meles Sus.scrofa Vulpes.vulpes
## 408 10124 1478#Sus scorfa abundance across subunit, sites and projectid
my_cols <- c("blue", "darkgreen", "red")
dotchart(P.anal$`Sus scrofa`, ylab = "Subunit",
groups = P.anal$Subunit, gcolor = my_cols,
color = my_cols[P.anal$Subunit],
cex = 0.9, pch = 1, xlab = "Abundance")
P.anal$Project.ID <- as.factor(P.anal$Project.ID)
my_cols <- c("blue", "darkgreen", "red", "purple", "darkorange")
dotchart(P.anal$`Sus scrofa`, ylab = "Project.ID",
groups = P.anal$Project.ID, gcolor = my_cols,
color = my_cols[P.anal$Project.ID],
cex = 0.9, pch = 1, xlab = "Abundance")
fun_color_range <- colorRampPalette(c("#1b98e0", "red"))
my_cols <- fun_color_range(15)
dotchart(P.anal$`Sus scrofa`, ylab = "Site",
groups = P.anal$Site, gcolor = my_cols,
color = my_cols[P.anal$Site],
cex = 0.9, pch = 1, xlab = "Abundance")
Species with 20 individuals or more - Canis aureus, Canis lupus, Gazella gazella, Hyaena hyaena, Hystrix indica, Meles meles, Sus scorfa and Vulpus vulpus
Left with 8 species
env_data <- abu_by_spp.Maquis[,..col_names]
mva_m0.po <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, family = "poisson", data = env_data)
mva_m0.nb <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, family = "negative.binomial", data = env_data)
aic_dt <- data.table(nb = mva_m0.nb$aic, po=mva_m0.po$aic)
colMeans(aic_dt)## nb po
## 2979.663 8038.338print("POISSON")## [1] "POISSON"plot(mva_m0.po, which=1:3)
print("NEGATIVE BINOMIAL")## [1] "NEGATIVE BINOMIAL"plot(mva_m0.nb, which=1:3)
prefer negative binomial because standardized residuals are smaller.
# mva_m0 <- mva_m0.po
# mva_m1 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + dunes*rescaled_Time.Diff + Distance_rescaled:dunes + cos_td_rad + sin_td_rad + site, family = "poisson", data = env_data)
# aic_dt$po1 <- mva_m1$aic
mva_m0 <- mva_m0.nb
mva_m1 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site + Subunit, data = env_data)
aic_dt$nb1 <- mva_m1$aic
colMeans(aic_dt)## nb po nb1
## 2979.663 8038.338 2807.215mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Transect, data = env_data)
aic_dt$nb2 <- mva_m2$aic
colMeans(aic_dt)## nb po nb1 nb2
## 2979.663 8038.338 2807.215 2762.678mva_m3 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, data = env_data)
aic_dt$nb3 <- mva_m2$aic
colMeans(aic_dt)## nb po nb1 nb2 nb3
## 2979.663 8038.338 2807.215 2762.678 2762.678mva_m4 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Subunit, data = env_data)
aic_dt$nb4 <- mva_m2$aic
colMeans(aic_dt)## nb po nb1 nb2 nb3 nb4
## 2979.663 8038.338 2807.215 2762.678 2762.678 2762.678Will proceed with subunit:
##Many GLM with subunit
mva_m0.nb <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Subunit, family = "negative.binomial", data = env_data)
drop1(mva_m0.nb)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff + Subunit
## Df AIC
## <none> 23014
## cosinus_Monitoring.Time.Diff 8 23051
## sinus_Monitoring.Time.Diff 8 23009
## Subunit 16 23837
## Distance_rescaled:rescaled_Time.Diff 8 23026Dropped sin time
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled + dunes*rescaled_Time.Diff + Distance_rescaled:dunes + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Subunit, family = "negative.binomial", data = env_data)
drop1(mva_m2)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## Subunit
## Df AIC
## <none> 23009
## cosinus_Monitoring.Time.Diff 8 23045
## Subunit 16 23828
## Distance_rescaled:rescaled_Time.Diff 8 23025# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled + dunes + rescaled_Time.Diff + Distance_rescaled:dunes + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, family = "negative.binomial", data = env_data)
# drop1(mva_m2)# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled + dunes + rescaled_Time.Diff + Distance_rescaled:dunes + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_rescaled + sinus_Monitoring.Time.Diff + Site, family = "negative.binomial", data = env_data)
# drop1(mva_m2)Final model includes Subunit, cosin time difference and the interaction between time and distance from settlement
# mva_m2 <- manyglm(formula = spp_no_rare ~ Dunes*rescaled_Time.Diff + sinus_Monitoring.Time.Diff + Transect_with_date, data = env_data)
# drop1(mva_m2)plot(mva_m2, which=1:3)
#summ_m2 <- summary(mva_m2)
#saveRDS(summ_m2, "output/summ_maquis_model_site.RDS")
summ_m2 <- readRDS("output/summ_maquis_model_site.RDS")
print(summ_m2)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 25.605 0.001 ***
## Distance_rescaled 13.896 0.001 ***
## rescaled_Time.Diff 4.126 0.106
## cosinus_Monitoring.Time.Diff 7.422 0.001 ***
## SubunitGalilee 14.087 0.001 ***
## SubunitJudea 24.436 0.001 ***
## Distance_rescaled:rescaled_Time.Diff 5.523 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 36.73, p-value: 0.001
## Arguments:
## Test statistics calculated assuming response assumed to be uncorrelated
## P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).#anov_m2.uni <- anova(mva_m2, p.uni = "adjusted")
#saveRDS(anov_m2.uni, "output/anova_maquis_model_site.RDS")
anov_m2.uni <- readRDS("output/anova_maquis_model_site.RDS")
print(anov_m2.uni)## Analysis of Deviance Table
##
## Model: spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Subunit
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 1116
## Distance_rescaled 1115 1 158.0 0.001 ***
## rescaled_Time.Diff 1114 1 33.8 0.002 **
## cosinus_Monitoring.Time.Diff 1113 1 131.1 0.001 ***
## Subunit 1111 2 844.9 0.001 ***
## Distance_rescaled:rescaled_Time.Diff 1110 1 32.3 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Canis.aureus Canis.lupus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 89.311 0.001 2.86 0.111
## rescaled_Time.Diff 2.027 0.706 1.096 0.835
## cosinus_Monitoring.Time.Diff 37.437 0.001 0.075 0.963
## Subunit 19.789 0.006 34.104 0.001
## Distance_rescaled:rescaled_Time.Diff 3.32 0.312 0.034 0.979
## Gazella.gazella Hyaena.hyaena
## Dev Pr(>Dev) Dev
## (Intercept)
## Distance_rescaled 12.227 0.016 17.54
## rescaled_Time.Diff 9.135 0.065 0.657
## cosinus_Monitoring.Time.Diff 40.549 0.001 19.628
## Subunit 91.903 0.001 71.141
## Distance_rescaled:rescaled_Time.Diff 3.771 0.299 7.084
## Hystrix.indica
## Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 0.004 4.298 0.102
## rescaled_Time.Diff 0.882 17.427 0.005
## cosinus_Monitoring.Time.Diff 0.003 0 0.999
## Subunit 0.001 20.295 0.005
## Distance_rescaled:rescaled_Time.Diff 0.122 9.791 0.057
## Meles.meles Sus.scrofa
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 8.019 0.059 6.114 0.068
## rescaled_Time.Diff 0.091 0.945 3.27 0.515
## cosinus_Monitoring.Time.Diff 0.231 0.959 31.727 0.001
## Subunit 32.585 0.002 572.242 0.001
## Distance_rescaled:rescaled_Time.Diff 0.002 0.979 7.274 0.122
## Vulpes.vulpes
## Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 17.64 0.004
## rescaled_Time.Diff 0.122 0.945
## cosinus_Monitoring.Time.Diff 1.497 0.719
## Subunit 2.821 0.434
## Distance_rescaled:rescaled_Time.Diff 1.029 0.710
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.Multivariate test shows that subunit, time, distance from settlements and timeXdistance affect species composition
coefp <- merge(data.table(t(coef(mva_m2)/log(2)),keep.rownames=TRUE), data.table(t(anov_m2.uni$uni.p), keep.rownames = TRUE), by = "rn") # coefficients are exponentiated because the link function used for poisson is log -> above 1 is positive and below 1 is negative
# add total species abundance
coefp <- merge(coefp,as.data.table(colSums(spp_no_rare),keep.rownames = TRUE), by.x = "rn", by.y = "V1")
colnames(coefp) <- c("SciName","Intercept.coef","Distance.coef","time_diff.coef","Cosin_time.coef","Galilee.coef","Judea.coef","Distance_time.coef","Intercept.p","Distance.p","time_diff.p","Cosin_time.p","Subunit.p","Distance_time.p","Species_abundance")
write.csv(coefp, "coefficients_Site_Maquis.csv")####Canis aureus####
Canis.aureus <- spp_no_rare[,"Canis.aureus"]
glm.Canis.aureus <- glm.nb(formula = Canis.aureus ~ Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
print(summary(glm.Canis.aureus))##
## Call:
## glm.nb(formula = Canis.aureus ~ Distance_rescaled + rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff + Site, data = env_data, init.theta = 0.3874865829,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7465 -1.2179 -0.6465 0.0679 3.7902
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.13230 0.21707 0.609 0.542200
## Distance_rescaled -1.03728 0.14420 -7.193 6.33e-13 ***
## rescaled_Time.Diff -0.00310 0.04860 -0.064 0.949143
## cosinus_Monitoring.Time.Diff -0.25188 0.08689 -2.899 0.003745 **
## SiteBeit Oren 0.71271 0.29790 2.392 0.016736 *
## SiteEin Yaakov 1.69092 0.28253 5.985 2.16e-09 ***
## SiteGivat Yearim 1.09014 0.31079 3.508 0.000452 ***
## SiteGivat Yeshayahu -0.76551 0.31365 -2.441 0.014660 *
## SiteGoren 1.35792 0.27592 4.921 8.60e-07 ***
## SiteIftach -0.23953 0.29588 -0.810 0.418205
## SiteKerem Maharal 0.46764 0.31274 1.495 0.134840
## SiteKfar Shamai 1.33553 0.28909 4.620 3.84e-06 ***
## SiteKisalon 0.98376 0.31268 3.146 0.001654 **
## SiteMargaliot 1.15096 0.29249 3.935 8.32e-05 ***
## SiteNehusha 0.71445 0.29414 2.429 0.015143 *
## SiteNir Etzion 0.85252 0.28637 2.977 0.002911 **
## SiteOfer 1.68926 0.28493 5.929 3.06e-09 ***
## SiteYagur 0.11416 0.29555 0.386 0.699312
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3875) family taken to be 1)
##
## Null deviance: 1440.3 on 1116 degrees of freedom
## Residual deviance: 1133.8 on 1099 degrees of freedom
## AIC: 5477.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3875
## Std. Err.: 0.0206
##
## 2 x log-likelihood: -5439.4360coef(glm.Canis.aureus)## (Intercept) Distance_rescaled
## 0.13230360 -1.03728194
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## -0.00309955 -0.25188430
## SiteBeit Oren SiteEin Yaakov
## 0.71271155 1.69091881
## SiteGivat Yearim SiteGivat Yeshayahu
## 1.09013769 -0.76551388
## SiteGoren SiteIftach
## 1.35791596 -0.23952726
## SiteKerem Maharal SiteKfar Shamai
## 0.46763507 1.33552993
## SiteKisalon SiteMargaliot
## 0.98375599 1.15096110
## SiteNehusha SiteNir Etzion
## 0.71445204 0.85252488
## SiteOfer SiteYagur
## 1.68925574 0.11415527coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719source("functions/plot_model_effect_mammals_Maquis_and_Forest.R")
plot_model_effect(P.anal = env_data, m = glm.Canis.aureus, eff2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
####Gazella gazella####
Gazella.gazella <- spp_no_rare[,"Gazella.gazella"]
glm.Gazella.gazella <- glm.nb(formula = Gazella.gazella ~ Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
print(summary(glm.Gazella.gazella))##
## Call:
## glm.nb(formula = Gazella.gazella ~ Distance_rescaled + rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff + Site, data = env_data, init.theta = 0.2833262366,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3984 -0.5595 -0.1237 0.0000 2.8167
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.369e+00 2.509e-01 5.455 4.89e-08 ***
## Distance_rescaled 1.701e+00 2.828e-01 6.016 1.79e-09 ***
## rescaled_Time.Diff 1.737e-01 1.011e-01 1.719 0.085564 .
## cosinus_Monitoring.Time.Diff 5.895e-01 1.584e-01 3.721 0.000199 ***
## SiteBeit Oren -3.690e+01 7.853e+06 0.000 0.999996
## SiteEin Yaakov -3.703e+01 7.279e+06 0.000 0.999996
## SiteGivat Yearim -1.761e+00 4.606e-01 -3.823 0.000132 ***
## SiteGivat Yeshayahu -9.362e-01 3.635e-01 -2.575 0.010013 *
## SiteGoren -3.765e+01 7.044e+06 0.000 0.999996
## SiteIftach -6.009e-01 3.775e-01 -1.592 0.111472
## SiteKerem Maharal -2.907e+00 4.909e-01 -5.921 3.19e-09 ***
## SiteKfar Shamai -2.744e+00 4.872e-01 -5.633 1.77e-08 ***
## SiteKisalon -1.716e+00 4.352e-01 -3.943 8.06e-05 ***
## SiteMargaliot -3.692e+01 7.411e+06 0.000 0.999996
## SiteNehusha -1.421e+00 3.940e-01 -3.607 0.000310 ***
## SiteNir Etzion -5.299e+00 1.070e+00 -4.952 7.34e-07 ***
## SiteOfer -4.122e+00 7.985e-01 -5.161 2.45e-07 ***
## SiteYagur -3.724e+01 7.713e+06 0.000 0.999996
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.2833) family taken to be 1)
##
## Null deviance: 747.04 on 1116 degrees of freedom
## Residual deviance: 337.70 on 1099 degrees of freedom
## AIC: 1105.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.2833
## Std. Err.: 0.0388
##
## 2 x log-likelihood: -1067.5850coef(glm.Gazella.gazella)## (Intercept) Distance_rescaled
## 1.3687213 1.7012666
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.1737436 0.5895141
## SiteBeit Oren SiteEin Yaakov
## -36.8963765 -37.0287932
## SiteGivat Yearim SiteGivat Yeshayahu
## -1.7610334 -0.9362387
## SiteGoren SiteIftach
## -37.6478873 -0.6008662
## SiteKerem Maharal SiteKfar Shamai
## -2.9068536 -2.7445015
## SiteKisalon SiteMargaliot
## -1.7157595 -36.9231722
## SiteNehusha SiteNir Etzion
## -1.4209898 -5.2993400
## SiteOfer SiteYagur
## -4.1215505 -37.2426408coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_effect(P.anal = env_data, m = glm.Gazella.gazella, eff2plot = "rescaled_Time.Diff", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
plot_model_effect(P.anal = env_data, m = glm.Gazella.gazella, eff2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
####Hyaena hyaena####
Hyaena.hyaena <- spp_no_rare[,"Hyaena.hyaena"]
glm.Hyaena.hyaena <- glm.nb(formula = Hyaena.hyaena ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
print(summary(glm.Hyaena.hyaena))##
## Call:
## glm.nb(formula = Hyaena.hyaena ~ Distance_rescaled * rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff + Site, data = env_data, init.theta = 0.3289665235,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.183 -0.657 -0.423 0.000 2.435
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.132e-01 2.906e-01 -3.143 0.00168 **
## Distance_rescaled 3.491e-01 2.460e-01 1.419 0.15590
## rescaled_Time.Diff -3.370e-01 1.178e-01 -2.862 0.00421 **
## cosinus_Monitoring.Time.Diff 4.032e-01 1.482e-01 2.720 0.00652 **
## SiteBeit Oren -4.907e-02 4.345e-01 -0.113 0.91008
## SiteEin Yaakov -3.006e+01 4.088e+05 0.000 0.99994
## SiteGivat Yearim 6.015e-01 4.092e-01 1.470 0.14156
## SiteGivat Yeshayahu -1.374e+00 5.055e-01 -2.718 0.00657 **
## SiteGoren -3.590e+00 1.072e+00 -3.349 0.00081 ***
## SiteIftach -6.941e-01 4.533e-01 -1.531 0.12571
## SiteKerem Maharal -2.505e-01 4.278e-01 -0.585 0.55823
## SiteKfar Shamai -3.005e+01 3.988e+05 0.000 0.99994
## SiteKisalon 5.001e-01 4.100e-01 1.220 0.22255
## SiteMargaliot -2.991e+01 4.159e+05 0.000 0.99994
## SiteNehusha -1.232e+00 5.173e-01 -2.382 0.01720 *
## SiteNir Etzion 2.241e-01 3.761e-01 0.596 0.55130
## SiteOfer 4.732e-01 3.846e-01 1.230 0.21863
## SiteYagur -5.129e-01 4.369e-01 -1.174 0.24043
## Distance_rescaled:rescaled_Time.Diff -4.505e-01 1.610e-01 -2.798 0.00514 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.329) family taken to be 1)
##
## Null deviance: 644.94 on 1116 degrees of freedom
## Residual deviance: 425.22 on 1098 degrees of freedom
## AIC: 1161.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3290
## Std. Err.: 0.0537
##
## 2 x log-likelihood: -1121.5250coef(glm.Hyaena.hyaena)## (Intercept) Distance_rescaled
## -0.91320224 0.34909622
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## -0.33700680 0.40315862
## SiteBeit Oren SiteEin Yaakov
## -0.04907179 -30.06320966
## SiteGivat Yearim SiteGivat Yeshayahu
## 0.60152996 -1.37381799
## SiteGoren SiteIftach
## -3.59008089 -0.69410055
## SiteKerem Maharal SiteKfar Shamai
## -0.25048355 -30.05032280
## SiteKisalon SiteMargaliot
## 0.50010434 -29.91031648
## SiteNehusha SiteNir Etzion
## -1.23236188 0.22407097
## SiteOfer SiteYagur
## 0.47315952 -0.51289805
## Distance_rescaled:rescaled_Time.Diff
## -0.45052533coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_effect(P.anal = env_data, m = glm.Hyaena.hyaena, eff2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
####Meles meles####
Meles.meles <- spp_no_rare[,"Meles.meles"]
glm.Meles.meles <- glm.nb(formula = Meles.meles ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
print(summary(glm.Meles.meles))##
## Call:
## glm.nb(formula = Meles.meles ~ Distance_rescaled * rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff + Site, data = env_data, init.theta = 0.1889359821,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9939 -0.6275 -0.5390 -0.3152 2.6116
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0005 0.4005 -4.995 5.88e-07 ***
## Distance_rescaled -0.6352 0.2623 -2.422 0.0154 *
## rescaled_Time.Diff -0.0646 0.1364 -0.474 0.6358
## cosinus_Monitoring.Time.Diff 0.3727 0.1619 2.302 0.0213 *
## SiteBeit Oren 0.5935 0.5355 1.108 0.2677
## SiteEin Yaakov -1.4415 0.6928 -2.081 0.0375 *
## SiteGivat Yearim 0.3480 0.5547 0.627 0.5304
## SiteGivat Yeshayahu 0.1425 0.5235 0.272 0.7855
## SiteGoren 0.1153 0.5109 0.226 0.8214
## SiteIftach -0.3040 0.5517 -0.551 0.5815
## SiteKerem Maharal 1.0478 0.5263 1.991 0.0465 *
## SiteKfar Shamai 1.0685 0.5222 2.046 0.0407 *
## SiteKisalon -1.3632 0.7787 -1.751 0.0800 .
## SiteMargaliot 0.3423 0.5502 0.622 0.5338
## SiteNehusha 0.4437 0.5288 0.839 0.4014
## SiteNir Etzion 0.4241 0.5070 0.837 0.4029
## SiteOfer 0.9736 0.5031 1.935 0.0530 .
## SiteYagur 2.1635 0.4850 4.460 8.19e-06 ***
## Distance_rescaled:rescaled_Time.Diff -0.1352 0.1960 -0.690 0.4904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1889) family taken to be 1)
##
## Null deviance: 580.10 on 1116 degrees of freedom
## Residual deviance: 481.81 on 1098 degrees of freedom
## AIC: 1452.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1889
## Std. Err.: 0.0241
##
## 2 x log-likelihood: -1412.5940coef(glm.Meles.meles)## (Intercept) Distance_rescaled
## -2.00052131 -0.63522354
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## -0.06459685 0.37269204
## SiteBeit Oren SiteEin Yaakov
## 0.59353778 -1.44153972
## SiteGivat Yearim SiteGivat Yeshayahu
## 0.34802774 0.14248741
## SiteGoren SiteIftach
## 0.11533686 -0.30404801
## SiteKerem Maharal SiteKfar Shamai
## 1.04783452 1.06851075
## SiteKisalon SiteMargaliot
## -1.36319446 0.34235121
## SiteNehusha SiteNir Etzion
## 0.44375014 0.42406921
## SiteOfer SiteYagur
## 0.97361906 2.16347478
## Distance_rescaled:rescaled_Time.Diff
## -0.13520454coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_effect(P.anal = env_data, m = glm.Meles.meles, eff2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
####Vulpus vulpus####
Vulpes.vulpes <- spp_no_rare[,"Vulpes.vulpes"]
glm.Vulpes.vulpes <- glm.nb(formula = Vulpes.vulpes ~ Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
print(summary(glm.Vulpes.vulpes))##
## Call:
## glm.nb(formula = Vulpes.vulpes ~ Distance_rescaled + rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff + Site, data = env_data, init.theta = 0.3068953754,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2986 -0.9602 -0.6704 -0.1543 3.5218
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.71399 0.23435 3.047 0.002314 **
## Distance_rescaled 0.33533 0.16337 2.053 0.040109 *
## rescaled_Time.Diff 0.11588 0.06130 1.890 0.058716 .
## cosinus_Monitoring.Time.Diff 0.06475 0.10915 0.593 0.553014
## SiteBeit Oren -1.51041 0.38398 -3.934 8.37e-05 ***
## SiteEin Yaakov -1.17786 0.34915 -3.373 0.000742 ***
## SiteGivat Yearim -2.13669 0.45986 -4.646 3.38e-06 ***
## SiteGivat Yeshayahu 0.76126 0.32294 2.357 0.018412 *
## SiteGoren -1.17833 0.33422 -3.526 0.000423 ***
## SiteIftach -0.01873 0.32823 -0.057 0.954484
## SiteKerem Maharal 0.30567 0.33678 0.908 0.364072
## SiteKfar Shamai 0.58302 0.32535 1.792 0.073135 .
## SiteKisalon -1.72820 0.42177 -4.098 4.18e-05 ***
## SiteMargaliot -1.60211 0.38010 -4.215 2.50e-05 ***
## SiteNehusha -0.47135 0.33974 -1.387 0.165315
## SiteNir Etzion -1.40796 0.35600 -3.955 7.66e-05 ***
## SiteOfer -0.11803 0.32703 -0.361 0.718157
## SiteYagur -0.05436 0.32966 -0.165 0.869033
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3069) family taken to be 1)
##
## Null deviance: 1020.03 on 1116 degrees of freedom
## Residual deviance: 822.95 on 1099 degrees of freedom
## AIC: 2935.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3069
## Std. Err.: 0.0237
##
## 2 x log-likelihood: -2897.2410coef(glm.Vulpes.vulpes)## (Intercept) Distance_rescaled
## 0.71398960 0.33532676
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.11587943 0.06475492
## SiteBeit Oren SiteEin Yaakov
## -1.51041470 -1.17785871
## SiteGivat Yearim SiteGivat Yeshayahu
## -2.13669084 0.76125523
## SiteGoren SiteIftach
## -1.17832816 -0.01873417
## SiteKerem Maharal SiteKfar Shamai
## 0.30567472 0.58302195
## SiteKisalon SiteMargaliot
## -1.72820136 -1.60211394
## SiteNehusha SiteNir Etzion
## -0.47135486 -1.40795969
## SiteOfer SiteYagur
## -0.11803307 -0.05435716coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_effect(P.anal = env_data, m = glm.Vulpes.vulpes, eff2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
#####Hystrix indica####
Hystrix.indica <- spp_no_rare[,"Hystrix.indica"]
glm.Hystrix.indica <- glm.nb(formula = Hystrix.indica ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, data = env_data)
coef(glm.Hystrix.indica)## (Intercept) Distance_rescaled
## 0.45860383 -0.52964293
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## -0.06413709 0.04939549
## SiteBeit Oren SiteEin Yaakov
## -0.03206429 -1.27421746
## SiteGivat Yearim SiteGivat Yeshayahu
## 0.13990662 1.00927181
## SiteGoren SiteIftach
## 0.02617104 -0.23012868
## SiteKerem Maharal SiteKfar Shamai
## 0.60189664 0.84253158
## SiteKisalon SiteMargaliot
## 0.17540887 0.61169662
## SiteNehusha SiteNir Etzion
## 0.64288936 -0.61634155
## SiteOfer SiteYagur
## 0.35859865 -0.90079811
## Distance_rescaled:rescaled_Time.Diff
## -0.35785205coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719source("functions/plot_interaction_mammals_Maquis_and_Forest.R")
plot_model_interaction(P.anal = env_data, m = glm.Hystrix.indica, eff2plot = "rescaled_Time.Diff", modvar2plot = "Distance_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")## Warning: -0.962252785766932 is outside the observed range of Distance_rescaled
sd_distance <- sd(env_data$Distance_rescaled)
print(sd_distance)## [1] 0.458914mean_distance <- mean(env_data$Distance_rescaled)
print(mean_distance)## [1] -0.5033388#calculate +/1 1 SD
small_distance <- mean_distance - sd_distance
large_distance <- mean_distance + sd_distance
#Test whether temporal trend is different from zero
mm_rescaled_time <- emtrends(glm.Hystrix.indica, specs = "Distance_rescaled", var = "rescaled_Time.Diff", type = "response", at = list(Distance_rescaled = c(small_distance,large_distance)))
print(mm_rescaled_time)## Distance_rescaled rescaled_Time.Diff.trend SE df asymp.LCL asymp.UCL
## -0.9623 0.2802 0.0659 Inf 0.151 0.4094
## -0.0444 -0.0482 0.0671 Inf -0.180 0.0833
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95test_results_mm_rescaled_time <- test(mm_rescaled_time, null = 0, adjust = "fdr")
print(test_results_mm_rescaled_time)## Distance_rescaled rescaled_Time.Diff.trend SE df z.ratio p.value
## -0.9623 0.2802 0.0659 Inf 4.250 <.0001
## -0.0444 -0.0482 0.0671 Inf -0.719 0.4722
##
## Results are averaged over the levels of: Site
## P value adjustment: fdr method for 2 tests#Testing whether abundance is higher near or far from settlements
pairwise_distance <- emmeans(object = glm.Hystrix.indica, ~ Distance_rescaled*rescaled_Time.Diff, at = list(Distance_rescaled = c(small_distance,large_distance)))
print(pairwise_distance)## Distance_rescaled rescaled_Time.Diff emmean SE df asymp.LCL asymp.UCL
## -0.9623 -0.262 0.981 0.0780 Inf 0.828 1.13
## -0.0444 -0.262 0.581 0.0815 Inf 0.421 0.74
##
## Results are averaged over the levels of: Site
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95test_results_distance <- test(pairs(pairwise_distance, by="rescaled_Time.Diff"), by=NULL, adjust="fdr")
print(test_results_distance)## contrast
## (Distance_rescaled-0.962252785766932) - (Distance_rescaled-0.0444248306413041)
## rescaled_Time.Diff estimate SE df z.ratio p.value
## -0.262 0.4 0.124 Inf 3.233 0.0012
##
## Results are averaged over the levels of: Site
## Results are given on the log (not the response) scale.##Many GLM with interactions
mva_m0.nb <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, family = "negative.binomial", data = env_data)
drop1(mva_m0.nb)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Subunit *
## Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff
## Df AIC
## <none> 22814
## cosinus_Monitoring.Time.Diff 8 22853
## sinus_Monitoring.Time.Diff 8 22813
## Distance_rescaled:rescaled_Time.Diff 8 22814
## Distance_rescaled:Subunit 16 22900
## rescaled_Time.Diff:Subunit 16 22924mva_m1 <- manyglm(formula = spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, family = "negative.binomial", data = env_data)
drop1(mva_m1)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Subunit *
## Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff
## Df AIC
## <none> 22813
## cosinus_Monitoring.Time.Diff 8 22851
## Distance_rescaled:rescaled_Time.Diff 8 22815
## Distance_rescaled:Subunit 16 22897
## rescaled_Time.Diff:Subunit 16 22920mva_m1 <- manyglm(formula = spp_no_rare ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, family = "negative.binomial", data = env_data)
drop1(mva_m1)## Single term deletions
##
## Model:
## spp_no_rare ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff +
## cosinus_Monitoring.Time.Diff
## Df AIC
## <none> 22815
## cosinus_Monitoring.Time.Diff 8 22850
## Subunit:Distance_rescaled 16 22912
## Subunit:rescaled_Time.Diff 16 22919Model includes cosin time, interaction between subunit and time and subunit and distance
plot(mva_m1, which=1:3)
#summ_m2 <- summary(mva_m1)
#saveRDS(summ_m2, "output/summ_maquis_model_subunit.RDS")
summ_m2 <- readRDS("output/summ_maquis_model_subunit.RDS")
print(summ_m2)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 25.768 0.001 ***
## SubunitGalilee 7.995 0.001 ***
## SubunitJudea 14.892 0.001 ***
## Distance_rescaled 12.740 0.001 ***
## rescaled_Time.Diff 7.723 0.001 ***
## cosinus_Monitoring.Time.Diff 7.453 0.001 ***
## SubunitGalilee:Distance_rescaled 5.768 0.002 **
## SubunitJudea:Distance_rescaled 10.044 0.001 ***
## SubunitGalilee:rescaled_Time.Diff 6.877 0.001 ***
## SubunitJudea:rescaled_Time.Diff 7.906 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 34.44, p-value: 0.001
## Arguments:
## Test statistics calculated assuming response assumed to be uncorrelated
## P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).#anov_m2.uni <- anova(mva_m1, p.uni = "adjusted")
#saveRDS(anov_m2.uni, "output/anova_maquis_model_subunit.RDS")
anov_m2.uni <- readRDS("output/anova_maquis_model_subunit.RDS")
print(anov_m2.uni)## Analysis of Deviance Table
##
## Model: spp_no_rare ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 1116
## Subunit 1114 2 863.6 0.001 ***
## Distance_rescaled 1113 1 199.4 0.001 ***
## rescaled_Time.Diff 1112 1 55.7 0.001 ***
## cosinus_Monitoring.Time.Diff 1111 1 49.2 0.001 ***
## Subunit:Distance_rescaled 1109 2 138.5 0.001 ***
## Subunit:rescaled_Time.Diff 1107 2 135.8 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Canis.aureus Canis.lupus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Subunit 36.167 0.001 33.343 0.001
## Distance_rescaled 92.681 0.001 1.134 0.517
## rescaled_Time.Diff 6.789 0.104 1.291 0.759
## cosinus_Monitoring.Time.Diff 12.928 0.028 2.367 0.617
## Subunit:Distance_rescaled 19.085 0.007 0.7 0.241
## Subunit:rescaled_Time.Diff 43.537 0.001 4.213 0.487
## Gazella.gazella Hyaena.hyaena
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Subunit 60.391 0.001 92.277 0.001
## Distance_rescaled 56.968 0.001 9.28 0.025
## rescaled_Time.Diff 10.739 0.025 0.496 0.895
## cosinus_Monitoring.Time.Diff 25.716 0.001 6.914 0.146
## Subunit:Distance_rescaled 4.176 0.194 12.171 0.012
## Subunit:rescaled_Time.Diff 0.989 0.716 2.5 0.704
## Hystrix.indica Meles.meles
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Subunit 22.772 0.001 32.007 0.001
## Distance_rescaled 0.472 0.517 8.625 0.025
## rescaled_Time.Diff 17.991 0.012 0.042 0.928
## cosinus_Monitoring.Time.Diff 0.783 0.911 0.252 0.969
## Subunit:Distance_rescaled 41.145 0.001 19.296 0.007
## Subunit:rescaled_Time.Diff 1.771 0.716 3.351 0.637
## Sus.scrofa Vulpes.vulpes
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Subunit 582.6 0.001 4.074 0.282
## Distance_rescaled 12.569 0.025 17.674 0.015
## rescaled_Time.Diff 18.183 0.012 0.123 0.928
## cosinus_Monitoring.Time.Diff 0.001 0.980 0.21 0.969
## Subunit:Distance_rescaled 23.227 0.006 18.657 0.007
## Subunit:rescaled_Time.Diff 74.121 0.001 5.291 0.439
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.coefp <- merge(data.table(t(coef(mva_m1)/log(2)),keep.rownames=TRUE), data.table(t(anov_m2.uni$uni.p), keep.rownames = TRUE), by = "rn") # coefficients are exponentiated because the link function used for poisson is log -> above 1 is positive and below 1 is negative
# add total species abundance
coefp <- merge(coefp,as.data.table(colSums(spp_no_rare),keep.rownames = TRUE), by.x = "rn", by.y = "V1")
colnames(coefp) <- c("SciName","Intercept.coef","Galilee.coef","Judea.coef","Distance.coef","Time_diff.coef","cosin_time.coef","Galilee_distance.coef","Judea_distance.coef","Galilee_time.coef","Judea_time.coef","Intercept.p","Subunit.p","Distance.p","Time_diff.p","Cosin_time.p","Subunit_distance.p","Subunit_time.p","Species_abundance")
write.csv(coefp, "coefficients_Subunit_interactions_Maquis.csv")#####Canis aureus####
Canis.aureus <- spp_no_rare[,"Canis.aureus"]
glm.Canis.aureus <- glm.nb(formula = Canis.aureus ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)
coef(glm.Canis.aureus)## (Intercept) SubunitGalilee
## 0.7379764 0.2887068
## SubunitJudea Distance_rescaled
## -1.0561383 -1.0768661
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## -0.4959726 -0.3461529
## SubunitGalilee:Distance_rescaled SubunitJudea:Distance_rescaled
## -0.4936521 -1.6603072
## SubunitGalilee:rescaled_Time.Diff SubunitJudea:rescaled_Time.Diff
## 0.7092476 0.5938064coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719source("functions/plot_interaction_mammals_Maquis_and_Forest.R")
plot_model_interaction(P.anal = env_data, m = glm.Canis.aureus, eff2plot = "Distance_rescaled", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
plot_model_interaction(P.anal = env_data, m = glm.Canis.aureus, eff2plot = "rescaled_Time.Diff", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
#Test whether temporal trend is different from zero
print("Test whether temporal trend is different from zero")## [1] "Test whether temporal trend is different from zero"mm_rescaled_time <- emtrends(glm.Canis.aureus, specs = "Subunit", var = "rescaled_Time.Diff", type = "response")
print(mm_rescaled_time)## Subunit rescaled_Time.Diff.trend SE df asymp.LCL asymp.UCL
## Carmel -0.4960 0.0876 Inf -0.6676 -0.324
## Galilee 0.2133 0.0816 Inf 0.0533 0.373
## Judea 0.0978 0.0974 Inf -0.0931 0.289
##
## Confidence level used: 0.95test_results_mm_rescaled_time <- test(mm_rescaled_time, null = 0, adjust = "fdr")
print(test_results_mm_rescaled_time)## Subunit rescaled_Time.Diff.trend SE df z.ratio p.value
## Carmel -0.4960 0.0876 Inf -5.663 <.0001
## Galilee 0.2133 0.0816 Inf 2.613 0.0134
## Judea 0.0978 0.0974 Inf 1.004 0.3153
##
## P value adjustment: fdr method for 3 tests#Test whether spatial trend is different from zero
print("Test whether spatial trend is different from zero")## [1] "Test whether spatial trend is different from zero"mm_rescaled_distance <- emtrends(glm.Canis.aureus, specs = "Subunit", var = "Distance_rescaled", type = "response")
print(mm_rescaled_distance)## Subunit Distance_rescaled.trend SE df asymp.LCL asymp.UCL
## Carmel -1.08 0.159 Inf -1.39 -0.765
## Galilee -1.57 0.257 Inf -2.08 -1.066
## Judea -2.74 0.354 Inf -3.43 -2.043
##
## Confidence level used: 0.95test_results_mm_distance <- test(mm_rescaled_distance, null = 0, adjust = "fdr")
print(test_results_mm_distance)## Subunit Distance_rescaled.trend SE df z.ratio p.value
## Carmel -1.08 0.159 Inf -6.759 <.0001
## Galilee -1.57 0.257 Inf -6.100 <.0001
## Judea -2.74 0.354 Inf -7.725 <.0001
##
## P value adjustment: fdr method for 3 tests#Compare spatial trends among subunits
print("Compare spatial trends among subunits")## [1] "Compare spatial trends among subunits"pairwise_distance <- emmeans(object = glm.Canis.aureus, ~ Subunit*Distance_rescaled)
print(pairwise_distance)## Subunit Distance_rescaled emmean SE df asymp.LCL asymp.UCL
## Carmel -0.503 1.44 0.0919 Inf 1.26 1.62
## Galilee -0.503 1.79 0.0889 Inf 1.62 1.97
## Judea -0.503 1.06 0.1042 Inf 0.86 1.27
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95test_results_distance <- test(pairs(pairwise_distance, by="Distance_rescaled"), by=NULL, adjust="fdr")
print(test_results_distance)## contrast Distance_rescaled estimate SE df z.ratio p.value
## Carmel - Galilee -0.503 -0.351 0.128 Inf -2.743 0.0067
## Carmel - Judea -0.503 0.376 0.139 Inf 2.711 0.0067
## Galilee - Judea -0.503 0.727 0.140 Inf 5.204 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests#####Hystrix indica####
Hystrix.indica <- spp_no_rare[,"Hystrix.indica"]
glm.Hystrix.indica <- glm.nb(formula = Hystrix.indica ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)
summary(glm.Hystrix.indica)##
## Call:
## glm.nb(formula = Hystrix.indica ~ Subunit * Distance_rescaled +
## Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff,
## data = env_data, init.theta = 0.4028111892, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5583 -1.2303 -0.5464 0.1358 3.1233
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.85541 0.10742 7.963 1.68e-15 ***
## SubunitGalilee -0.06300 0.20344 -0.310 0.75679
## SubunitJudea -0.69582 0.23328 -2.983 0.00286 **
## Distance_rescaled 0.32578 0.13876 2.348 0.01889 *
## rescaled_Time.Diff 0.30195 0.08675 3.480 0.00050 ***
## cosinus_Monitoring.Time.Diff -0.19146 0.08162 -2.346 0.01899 *
## SubunitGalilee:Distance_rescaled -0.46150 0.28653 -1.611 0.10726
## SubunitJudea:Distance_rescaled -2.27715 0.35563 -6.403 1.52e-10 ***
## SubunitGalilee:rescaled_Time.Diff -0.16600 0.12063 -1.376 0.16879
## SubunitJudea:rescaled_Time.Diff -0.08044 0.12467 -0.645 0.51880
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4028) family taken to be 1)
##
## Null deviance: 1166.1 on 1116 degrees of freedom
## Residual deviance: 1075.9 on 1107 degrees of freedom
## AIC: 4433.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4028
## Std. Err.: 0.0245
##
## 2 x log-likelihood: -4411.7900Interaction of subunit with time is not significant - will be dropped from H. indica model
glm.Hystrix.indica <- glm.nb(formula = Hystrix.indica ~ Subunit * Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)
summary(glm.Hystrix.indica)##
## Call:
## glm.nb(formula = Hystrix.indica ~ Subunit * Distance_rescaled +
## rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data,
## init.theta = 0.4018098414, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5446 -1.2182 -0.5649 0.1539 3.2922
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.82151 0.10519 7.810 5.73e-15 ***
## SubunitGalilee 0.00744 0.19850 0.037 0.97010
## SubunitJudea -0.65891 0.23035 -2.860 0.00423 **
## Distance_rescaled 0.29731 0.13844 2.148 0.03174 *
## rescaled_Time.Diff 0.21144 0.04902 4.314 1.61e-05 ***
## cosinus_Monitoring.Time.Diff -0.15609 0.07729 -2.020 0.04343 *
## SubunitGalilee:Distance_rescaled -0.42307 0.28597 -1.479 0.13903
## SubunitJudea:Distance_rescaled -2.22784 0.35518 -6.272 3.55e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4018) family taken to be 1)
##
## Null deviance: 1164.2 on 1116 degrees of freedom
## Residual deviance: 1076.0 on 1109 degrees of freedom
## AIC: 4431.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4018
## Std. Err.: 0.0244
##
## 2 x log-likelihood: -4413.5610coef(glm.Hystrix.indica)## (Intercept) SubunitGalilee
## 0.821514189 0.007440017
## SubunitJudea Distance_rescaled
## -0.658913418 0.297311123
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.211444812 -0.156093288
## SubunitGalilee:Distance_rescaled SubunitJudea:Distance_rescaled
## -0.423075330 -2.227842340coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_interaction(P.anal = env_data, m = glm.Hystrix.indica, eff2plot = "Distance_rescaled", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
#Test whether spatial trend is different from zero
mm_rescaled_distance <- emtrends(glm.Hystrix.indica, specs = "Subunit", var = "Distance_rescaled", type = "response")
print(mm_rescaled_distance)## Subunit Distance_rescaled.trend SE df asymp.LCL asymp.UCL
## Carmel 0.297 0.138 Inf 0.026 0.569
## Galilee -0.126 0.250 Inf -0.615 0.364
## Judea -1.931 0.326 Inf -2.570 -1.291
##
## Confidence level used: 0.95test_results_distance <- test(mm_rescaled_distance, null = 0, adjust="fdr")
print(test_results_distance)## Subunit Distance_rescaled.trend SE df z.ratio p.value
## Carmel 0.297 0.138 Inf 2.148 0.0476
## Galilee -0.126 0.250 Inf -0.504 0.6146
## Judea -1.931 0.326 Inf -5.915 <.0001
##
## P value adjustment: fdr method for 3 testsplot_model_effect(P.anal = env_data, m = glm.Hystrix.indica, eff2plot = "rescaled_Time.Diff", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/maquis/")
#####Meles meles####
Meles.meles <- spp_no_rare[,"Meles.meles"]
glm.Meles.meles <- glm.nb(formula = Meles.meles ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)
coef(glm.Meles.meles)## (Intercept) SubunitGalilee
## -0.76042055 -0.50900866
## SubunitJudea Distance_rescaled
## -2.44602759 -0.67269550
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.19389911 -0.02148272
## SubunitGalilee:Distance_rescaled SubunitJudea:Distance_rescaled
## 1.52725051 -2.15632118
## SubunitGalilee:rescaled_Time.Diff SubunitJudea:rescaled_Time.Diff
## -0.39715549 -0.22230256coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_interaction(P.anal = env_data, m = glm.Meles.meles, eff2plot = "Distance_rescaled", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
#Test whether spatial trend is different from zero
mm_distance<- emtrends(glm.Meles.meles, specs = "Subunit", var = "Distance_rescaled", type = "response")
print(mm_distance)## Subunit Distance_rescaled.trend SE df asymp.LCL asymp.UCL
## Carmel -0.673 0.261 Inf -1.1836 -0.162
## Galilee 0.855 0.459 Inf -0.0444 1.753
## Judea -2.829 0.773 Inf -4.3431 -1.315
##
## Confidence level used: 0.95test_results_distance <- test(mm_distance, null = 0, adjust = "fdr")
print(test_results_distance)## Subunit Distance_rescaled.trend SE df z.ratio p.value
## Carmel -0.673 0.261 Inf -2.581 0.0148
## Galilee 0.855 0.459 Inf 1.863 0.0624
## Judea -2.829 0.773 Inf -3.662 0.0008
##
## P value adjustment: fdr method for 3 tests####Sus scrofa####
Sus.scrofa <- spp_no_rare[,"Sus.scrofa"]
glm.Sus.scrofa <- glm.nb(formula = Sus.scrofa ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)
coef(glm.Sus.scrofa)## (Intercept) SubunitGalilee
## 1.87496631 0.83817303
## SubunitJudea Distance_rescaled
## -5.38303071 -0.26524313
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.10444555 -0.17401554
## SubunitGalilee:Distance_rescaled SubunitJudea:Distance_rescaled
## 0.09849245 -1.87063640
## SubunitGalilee:rescaled_Time.Diff SubunitJudea:rescaled_Time.Diff
## -0.07140282 1.94737676coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_interaction(P.anal = env_data, m = glm.Sus.scrofa, eff2plot = "Distance_rescaled", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
sus_distance <- interact_plot(model = glm.Sus.scrofa, data=env_data, pred = Distance_rescaled, modx = Subunit, main.title = "", interval = T,plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = 'Dark2', modxvals = NULL, line.colors = 'black', point.alpha = 0.25)
sus_distance
# as.data.table(sus_distance$data)[Subunit=="Judea",Sus.scrofa]
# 1 - 0.0005676190/0.1249574967
plot_model_interaction(P.anal = env_data, m = glm.Sus.scrofa, eff2plot = "rescaled_Time.Diff", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
interact_plot(model = glm.Sus.scrofa, data=env_data, pred = rescaled_Time.Diff, modx = Subunit, main.title = "", interval = T,plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = 'Dark2', modxvals = NULL, line.colors = 'black', point.alpha = 0.25)
#Test whether temporal trend is different from zero
mm_rescaled_time <- emtrends(glm.Sus.scrofa, specs = "Subunit", var = "rescaled_Time.Diff", type = "response")
print(mm_rescaled_time)## Subunit rescaled_Time.Diff.trend SE df asymp.LCL asymp.UCL
## Carmel 0.104 0.0695 Inf -0.0317 0.241
## Galilee 0.033 0.0653 Inf -0.0949 0.161
## Judea 2.052 0.3233 Inf 1.4182 2.685
##
## Confidence level used: 0.95test_results_mm_rescaled_time <- test(mm_rescaled_time, null = 0, adjust = "fdr")
print(test_results_mm_rescaled_time)## Subunit rescaled_Time.Diff.trend SE df z.ratio p.value
## Carmel 0.104 0.0695 Inf 1.504 0.1990
## Galilee 0.033 0.0653 Inf 0.506 0.6127
## Judea 2.052 0.3233 Inf 6.347 <.0001
##
## P value adjustment: fdr method for 3 tests#Test whether spatial trend is different from zero
mm_rescaled_distance <- emtrends(glm.Sus.scrofa, specs = "Subunit", var = "Distance_rescaled", type = "response")
print(mm_rescaled_distance)## Subunit Distance_rescaled.trend SE df asymp.LCL asymp.UCL
## Carmel -0.265 0.115 Inf -0.491 -0.0393
## Galilee -0.167 0.200 Inf -0.558 0.2249
## Judea -2.136 0.623 Inf -3.357 -0.9147
##
## Confidence level used: 0.95test_results_mm_distance <- test(mm_rescaled_distance, null = 0, adjust = "fdr")
print(test_results_mm_distance)## Subunit Distance_rescaled.trend SE df z.ratio p.value
## Carmel -0.265 0.115 Inf -2.301 0.0321
## Galilee -0.167 0.200 Inf -0.834 0.4041
## Judea -2.136 0.623 Inf -3.428 0.0018
##
## P value adjustment: fdr method for 3 testsjudea <- as.data.table(P.anal[Subunit=="Judea"])
#####Vulpes vulpes####
Vulpes.vulpes <- spp_no_rare[,"Vulpes.vulpes"]
glm.Vulpes.vulpes <- glm.nb(formula = Vulpes.vulpes ~ Subunit * Distance_rescaled + Subunit * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff, data = env_data)## Warning: glm.fit: algorithm did not convergecoef(glm.Vulpes.vulpes)## (Intercept) SubunitGalilee
## 0.56496673 -0.70024393
## SubunitJudea Distance_rescaled
## -0.63159581 0.94003503
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 0.21422379 0.01986557
## SubunitGalilee:Distance_rescaled SubunitJudea:Distance_rescaled
## -1.26906026 -1.73205012
## SubunitGalilee:rescaled_Time.Diff SubunitJudea:rescaled_Time.Diff
## -0.32932902 -0.21039373coef(mva_m2)## Canis.aureus Canis.lupus Gazella.gazella
## (Intercept) 0.7836528 -6.6568237 -2.4838399
## Distance_rescaled -1.4509961 -1.2351726 1.9030344
## rescaled_Time.Diff -0.2145653 0.4338200 0.1655719
## cosinus_Monitoring.Time.Diff -0.2993373 0.8389958 0.8634291
## SubunitGalilee 0.3145271 3.8362076 1.8341808
## SubunitJudea -0.3647187 -8.8886543 3.0389989
## Distance_rescaled:rescaled_Time.Diff -0.1939343 -0.1963558 -0.3172844
## Hyaena.hyaena Hystrix.indica Meles.meles
## (Intercept) -0.9247905 0.561179943 -0.752229612
## Distance_rescaled 0.1630534 -0.269042397 -0.587767502
## rescaled_Time.Diff -0.3160072 0.007704151 -0.018706982
## cosinus_Monitoring.Time.Diff 0.4196223 -0.065554746 0.073613614
## SubunitGalilee -2.3879468 0.232072025 -1.165114581
## SubunitJudea -0.1874807 0.606274192 -1.039321245
## Distance_rescaled:rescaled_Time.Diff -0.4406182 -0.387492686 0.007278454
## Sus.scrofa Vulpes.vulpes
## (Intercept) 1.82208786 0.37034817
## Distance_rescaled -0.39646487 0.46422533
## rescaled_Time.Diff 0.07361119 -0.04014968
## cosinus_Monitoring.Time.Diff 0.01863800 0.06875367
## SubunitGalilee 0.88245266 0.04265783
## SubunitJudea -3.32773363 0.27749662
## Distance_rescaled:rescaled_Time.Diff -0.20479750 -0.11397719plot_model_interaction(P.anal = env_data, m = glm.Vulpes.vulpes, eff2plot = "Distance_rescaled", modvar2plot = "Subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA",
fontsize=22, pdf_width=160, outpath = "output/maquis/")
#Test whether spatial trend is different from zero
mm_rescaled_distance <- emtrends(glm.Vulpes.vulpes, specs = "Subunit", var = "Distance_rescaled", type = "response")
print(mm_rescaled_distance)## Subunit Distance_rescaled.trend SE df asymp.LCL asymp.UCL
## Carmel 0.940 0.177 Inf 0.593 1.2874
## Galilee -0.329 0.337 Inf -0.990 0.3320
## Judea -0.792 0.419 Inf -1.614 0.0298
##
## Confidence level used: 0.95test_results_mm_distance <- test(mm_rescaled_distance, null = 0, adjust = "fdr")
print(test_results_mm_distance)## Subunit Distance_rescaled.trend SE df z.ratio p.value
## Carmel 0.940 0.177 Inf 5.304 <.0001
## Galilee -0.329 0.337 Inf -0.976 0.3293
## Judea -0.792 0.419 Inf -1.889 0.0883
##
## P value adjustment: fdr method for 3 testsn = 1117 (total number of cameras)
total abundance =
abundance of common species (8 species) = 22,748