Mammal data - Arid Desert
Talia Gabay & Ron Chen
2024-04-09
P <- P_byplot[grepl("Arid Desert",Unit),][order(Cycle_number,Site)][,Site:=factor(Site)]Arid Desert - Richness by camera
Monitoring started in 2014 (5 monitoring cycles), Each site as two types of plots - near or far from Agriculture. There are 5 sites total - Ein Yahav, Lotan, Paran, Yotvata, Zofar Each has 5 sampling points except for Paran (near) which is missing samples from the first year (2014)
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$Agriculture <- as.factor(P.anal$Agriculture)
print("Agriculture has 2 levels")## [1] "Agriculture has 2 levels"print(levels(P.anal$Agriculture))## [1] "Far" "Near"P.anal$Site <- as.factor(P.anal$Site)
print("Site has 5 levels")## [1] "Site has 5 levels"print(levels(P.anal$Site))## [1] "Ein Yahav" "Lotan" "Paran" "Yotvata" "Zofar"P.anal$Transect <- as.factor(P.anal$Transect)
print("Transect has 10 levels")## [1] "Transect has 10 levels"print(levels(P.anal$Transect))## [1] "Ein Yahav Far" "Ein Yahav Near" "Lotan Far" "Lotan Near"
## [5] "Paran Far" "Paran Near" "Yotvata Far" "Yotvata Near"
## [9] "Zofar Far" "Zofar Near"P.anal$Transect_with_date <- as.factor(P.anal$Transect_with_date)
print("Transect with date has 49 levels - Missing data for Sde Boker near in 2014")## [1] "Transect with date has 49 levels - Missing data for Sde Boker near in 2014"print(levels(P.anal$Transect_with_date))## [1] "Ein Yahav Far_06_10_2016" "Ein Yahav Far_07_10_2018"
## [3] "Ein Yahav Far_12_06_2014" "Ein Yahav Far_13_09_2020"
## [5] "Ein Yahav Far_16_07_2022" "Ein Yahav Near_06_10_2016"
## [7] "Ein Yahav Near_07_10_2018" "Ein Yahav Near_12_06_2014"
## [9] "Ein Yahav Near_13_09_2020" "Ein Yahav Near_16_07_2022"
## [11] "Lotan Far_12_09_2016" "Lotan Far_21_02_2019"
## [13] "Lotan Far_22_05_2014" "Lotan Far_24_08_2022"
## [15] "Lotan Far_25_09_2020" "Lotan Near_12_09_2016"
## [17] "Lotan Near_21_02_2019" "Lotan Near_22_05_2014"
## [19] "Lotan Near_24_08_2022" "Lotan Near_25_09_2020"
## [21] "Paran Far_09_08_2022" "Paran Far_11_09_2018"
## [23] "Paran Far_16_11_2016" "Paran Far_24_09_2020"
## [25] "Paran Far_26_05_2014" "Paran Near_09_08_2022"
## [27] "Paran Near_11_09_2018" "Paran Near_16_11_2016"
## [29] "Paran Near_24_09_2020" "Yotvata Far_05_10_2020"
## [31] "Yotvata Far_06_09_2022" "Yotvata Far_06_10_2018"
## [33] "Yotvata Far_11_12_2016" "Yotvata Far_22_05_2014"
## [35] "Yotvata Near_05_10_2020" "Yotvata Near_06_09_2022"
## [37] "Yotvata Near_06_10_2018" "Yotvata Near_11_12_2016"
## [39] "Yotvata Near_16_05_2014" "Zofar Far_11_09_2018"
## [41] "Zofar Far_12_06_2014" "Zofar Far_14_09_2020"
## [43] "Zofar Far_16_08_2016" "Zofar Far_28_07_2022"
## [45] "Zofar Near_11_09_2018" "Zofar Near_12_06_2014"
## [47] "Zofar Near_14_09_2020" "Zofar Near_16_08_2016"
## [49] "Zofar Near_28_07_2022"#distance.rescaled <- scale(all.mammal.abundances$Distance)
#write.csv(distance.rescaled, file = "Distance_agri_rescaled.csv")print("RICHNESS WITH RARE SPECIES")## [1] "RICHNESS WITH RARE SPECIES"plot_alpha_diversity(P, x_val = "Agriculture", y_val = "richness", ylab_val = "richness", xlab_val = "Agriculture", fill_val = "Agriculture")
col_names <- c("Cycle_number", "Deployment.id_new", "Unit", "Subunit", "Site", "Distance_agri_rescaled","Lon", "Lat", "year", "rescaled_Time.Diff", "sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "Transect", "Transect_with_date")
abu_by_spp.AridDesert <- all.mammal.abundances[grepl("Arid Desert",Unit),]
spp <- abu_by_spp.AridDesert[,27:40]
#spp <- abu_by_spp.AridDesert[,25:38]
# filter out species with zero counts
print(colSums(spp))## Canis aureus Canis lupus Capra nubiana Dama mesopotamica
## 264 237 394 0
## Equus hemionus Gazella dorcas Gazella gazella Hyaena hyaena
## 100 733 0 96
## Hystrix indica Lepus capensis Meles meles Oryx leucoryx
## 30 662 0 72
## Sus scrofa Vulpes vulpes
## 0 447spp <- spp[,.SD,.SDcols = colSums(spp)>0]
spp <- mvabund(spp)
plot(spp ~ P.anal$Agriculture, 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")
dotchart(P.anal$richness, ylab = "Agriculture",
groups = P.anal$Agriculture, gcolor = my_cols,
color = my_cols[P.anal$Agriculture],
cex = 0.9, pch = 1, xlab = "richness")
my_cols <- c("blue", "darkgreen", "red", "orange", "purple")
dotchart(P.anal$richness, ylab = "Site",
groups = P.anal$Site, gcolor = my_cols,
color = my_cols[P.anal$Site],
cex = 0.9, pch = 1, xlab = "richness")
kable(summary(P.anal[,.(richness, abundance, rescaled_Time.Diff, cosinus_Monitoring.Time.Diff, sinus_Monitoring.Time.Diff, Site, Transect, Transect_with_date, Distance_agri_rescaled)]))| richness | abundance | rescaled_Time.Diff | cosinus_Monitoring.Time.Diff | sinus_Monitoring.Time.Diff | Site | Transect | Transect_with_date | Distance_agri_rescaled | |
|---|---|---|---|---|---|---|---|---|---|
| Min. :1.000 | Min. : 1.00 | Min. :-1.2473 | Min. :-0.9844 | Min. :-0.9538 | Ein Yahav:59 | Lotan Far : 38 | Ein Yahav Far_16_07_2022 : 9 | Min. :-0.68190 | |
| 1st Qu.:1.000 | 1st Qu.: 2.00 | 1st Qu.:-0.4567 | 1st Qu.:-0.1804 | 1st Qu.:-0.4121 | Lotan :64 | Yotvata Near : 34 | Ein Yahav Near_16_07_2022: 9 | 1st Qu.:-0.47903 | |
| Median :2.000 | Median : 5.00 | Median : 0.1395 | Median : 0.3755 | Median :-0.1324 | Paran :51 | Zofar Far : 34 | Lotan Far_24_08_2022 : 9 | Median : 0.03654 | |
| Mean :1.884 | Mean : 10.02 | Mean : 0.2774 | Mean : 0.2386 | Mean : 0.1085 | Yotvata :67 | Yotvata Far : 33 | Paran Near_16_11_2016 : 9 | Mean : 0.19731 | |
| 3rd Qu.:2.000 | 3rd Qu.: 11.00 | 3rd Qu.: 1.3302 | 3rd Qu.: 0.8532 | 3rd Qu.: 0.9200 | Zofar :62 | Ein Yahav Far: 31 | Yotvata Far_06_10_2018 : 9 | 3rd Qu.: 0.79490 | |
| Max. :6.000 | Max. :118.00 | Max. : 1.3751 | Max. : 0.9912 | Max. : 0.9995 | NA | Paran Far : 29 | Yotvata Near_06_09_2022 : 9 | Max. : 2.37740 | |
| NA | NA | NA | NA | NA | NA | (Other) :104 | (Other) :249 | NA |
pairs(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "rescaled_Time.Diff","Site","Transect", "Transect_with_date", "Distance_agri_rescaled")])
kable(cor(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("sinus_Monitoring.Time.Diff", "cosinus_Monitoring.Time.Diff", "rescaled_Time.Diff","Site","Transect", "Transect_with_date", "Distance_agri_rescaled")]))| sinus_Monitoring.Time.Diff | cosinus_Monitoring.Time.Diff | rescaled_Time.Diff | Site | Transect | Transect_with_date | Distance_agri_rescaled | |
|---|---|---|---|---|---|---|---|
| sinus_Monitoring.Time.Diff | 1.0000000 | -0.0086829 | -0.0362440 | -0.3317552 | -0.3358655 | -0.3280949 | 0.0415841 |
| cosinus_Monitoring.Time.Diff | -0.0086829 | 1.0000000 | 0.2815373 | -0.1389841 | -0.1416740 | -0.1498081 | -0.0534359 |
| rescaled_Time.Diff | -0.0362440 | 0.2815373 | 1.0000000 | 0.0201873 | 0.0375285 | 0.0539297 | -0.0584075 |
| Site | -0.3317552 | -0.1389841 | 0.0201873 | 1.0000000 | 0.9851579 | 0.9795745 | 0.0882773 |
| Transect | -0.3358655 | -0.1416740 | 0.0375285 | 0.9851579 | 1.0000000 | 0.9950592 | -0.0436044 |
| Transect_with_date | -0.3280949 | -0.1498081 | 0.0539297 | 0.9795745 | 0.9950592 | 1.0000000 | -0.0566318 |
| Distance_agri_rescaled | 0.0415841 | -0.0534359 | -0.0584075 | 0.0882773 | -0.0436044 | -0.0566318 | 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")
mdl_r.poiss.int <- glm(data = P.anal, formula = richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site + Transect_with_date, family = poisson)
od <- mdl_r.poiss.int$deviance/mdl_r.poiss.int$df.residual
print("Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model.")## [1] "Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model."print(paste0('od = ',od))## [1] "od = 0.324325977068647"Overdispersion parameter is <1 (underdispersion) and therefore will use Poisson
m0 <- glm(data = P.anal, formula = richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, family = poisson)
me0transectwithdate <- glmer(data = P.anal, formula = richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Transect_with_date), family = poisson, verbose = 1)## start par. = 1 fn = 948.8433
## At return
## eval: 15 fn: 871.06563 par: 0.00000
## (NM) 20: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 40: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 60: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 80: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 100: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 120: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 140: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 160: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 180: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 200: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 220: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788## boundary (singular) fit: see help('isSingular')me0transect <- glmer(data = P.anal, formula = richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Transect), family = poisson, verbose = 1)## start par. = 1 fn = 904.4763
## At return
## eval: 14 fn: 871.06563 par: 0.00000
## (NM) 20: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 40: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 60: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 80: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 100: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 120: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 140: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 160: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 180: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788## boundary (singular) fit: see help('isSingular')me0site <- glmer(data = P.anal, formula = richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Site), family = poisson, verbose = 1)## start par. = 1 fn = 891.3006
## At return
## eval: 14 fn: 871.06563 par: 0.00000
## (NM) 20: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 40: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 60: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 80: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 100: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 120: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 140: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 160: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788
## (NM) 180: f = 871.066 at 0 0.51983 0.168859 0.182429 0.0502167 -0.0250879 -0.0165788## boundary (singular) fit: see help('isSingular')All GLMERs result in singular fit, will use regular GLM:
m1 <- step(m0)## Start: AIC=883.07
## richness ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - Distance_agri_rescaled:rescaled_Time.Diff 1 128.69 881.13
## - sinus_Monitoring.Time.Diff 1 128.78 881.22
## - cosinus_Monitoring.Time.Diff 1 129.26 881.70
## <none> 128.62 883.07
##
## Step: AIC=881.13
## richness ~ Distance_agri_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - sinus_Monitoring.Time.Diff 1 128.82 879.26
## - cosinus_Monitoring.Time.Diff 1 129.35 879.79
## <none> 128.69 881.13
## - Distance_agri_rescaled 1 138.09 888.53
## - rescaled_Time.Diff 1 140.82 891.26
##
## Step: AIC=879.26
## richness ~ Distance_agri_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - cosinus_Monitoring.Time.Diff 1 129.51 877.95
## <none> 128.82 879.26
## - Distance_agri_rescaled 1 138.19 886.64
## - rescaled_Time.Diff 1 141.14 889.58
##
## Step: AIC=877.95
## richness ~ Distance_agri_rescaled + rescaled_Time.Diff
##
## Df Deviance AIC
## <none> 129.51 877.95
## - Distance_agri_rescaled 1 138.64 885.08
## - rescaled_Time.Diff 1 144.35 890.80Final model includes distance from agriculture and time difference
summ(m1, exp = TRUE, digits=3)| Observations | 303 |
| Dependent variable | richness |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| χ²(2) | 22.896 |
| Pseudo-R² (Cragg-Uhler) | 0.077 |
| Pseudo-R² (McFadden) | 0.026 |
| AIC | 877.952 |
| BIC | 889.093 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.699 | 1.543 | 1.871 | 10.796 | 0.000 |
| Distance_agri_rescaled | 1.171 | 1.059 | 1.296 | 3.065 | 0.002 |
| rescaled_Time.Diff | 1.207 | 1.096 | 1.330 | 3.813 | 0.000 |
| Standard errors: MLE |
effect_plot(model = m1, data=P.anal, pred = rescaled_Time.Diff, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1", main.title = "Temporal effect (points are partial residuals)")
effect_plot(model = m1, data=P.anal, pred = Distance_agri_rescaled, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1", main.title = "Effect of distance from agriculture")
Time difference and distance from agriculture significantly affect species richness
community analysis using package MVabund
Table of abundances per species by camera ID
abundance_and_cameras <- as.data.table(abu_by_spp.AridDesert[,c(5,27:40)])
dim(abundance_and_cameras)## [1] 303 15canisau <- abundance_and_cameras[,2] > 0
sum(canisau)## [1] 37canislup <- abundance_and_cameras[,3] > 0
sum(canislup)## [1] 52capra <- abundance_and_cameras[,4] > 0
sum(capra)## [1] 37Doma <- abundance_and_cameras[,5] > 0
sum(Doma)## [1] 0Equus <- abundance_and_cameras[,6] > 0
sum(Equus)## [1] 18Gazellador <- abundance_and_cameras[,7] > 0
sum(Gazellador)## [1] 106Gazellagaz <- abundance_and_cameras[,8] > 0
sum(Gazellagaz)## [1] 0Hyaena <- abundance_and_cameras[,9] > 0
sum(Hyaena)## [1] 69Hystrix <- abundance_and_cameras[,10] > 0
sum(Hystrix)## [1] 10Lepus <- abundance_and_cameras[,11] > 0
sum(Lepus)## [1] 109Meles <- abundance_and_cameras[,12] > 0
sum(Meles)## [1] 0Oryx <- abundance_and_cameras[,13] > 0
sum(Oryx)## [1] 8Sus <- abundance_and_cameras[,14] > 0
sum(Sus)## [1] 0Vulpus <- abundance_and_cameras[,15] > 0
sum(Vulpus)## [1] 125303 cameras total. Species with a minimum of 20 cameras - Canis aureus, Canis lupus, Capra nubiana, Gazella dorcas, Hyaena hyaena, Lepus capensis and vulpus vulpus
spp_no_rare <- abu_by_spp.AridDesert[,27:40]
print(colSums(spp_no_rare))## Canis aureus Canis lupus Capra nubiana Dama mesopotamica
## 264 237 394 0
## Equus hemionus Gazella dorcas Gazella gazella Hyaena hyaena
## 100 733 0 96
## Hystrix indica Lepus capensis Meles meles Oryx leucoryx
## 30 662 0 72
## Sus scrofa Vulpes vulpes
## 0 447# 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 <- mvabund(spp_no_rare[,c(-4,-7,-9)])
print(colSums(spp_no_rare))## Canis.aureus Canis.lupus Capra.nubiana Gazella.dorcas Hyaena.hyaena
## 264 237 394 733 96
## Lepus.capensis Vulpes.vulpes
## 662 447Left with 7 species
env_data <- abu_by_spp.AridDesert[,..col_names]
env_data$Site <- factor(env_data$Site, levels=c("Yotvata","Lotan","Paran", "Zofar", "Ein Yahav"))
mva_m0.po <- manyglm(formula = spp_no_rare ~ Distance_agri_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_agri_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
## 630.101 1517.274print("POISSON")## [1] "POISSON"plot(mva_m0.po, which=1:3)## Warning in default.plot.manyglm(x, res.type = res.type, which = which, caption
## = caption, : Only the first 7 colors will be used for plotting.
print("NEGATIVE BINOMIAL")## [1] "NEGATIVE BINOMIAL"plot(mva_m0.nb, which=1:3)## Warning in default.plot.manyglm(x, res.type = res.type, which = which, caption
## = caption, : Only the first 7 colors will be used for plotting.
prefer negative binomial because standardized residuals are smaller.
# mva_m0 <- mva_m0.po
# mva_m1 <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + dunes*rescaled_Time.Diff + Distance_agri_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_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Transect_with_date, data = env_data)
aic_dt$nb1 <- mva_m1$aic
colMeans(aic_dt)## nb po nb1
## 630.1010 1517.2736 577.3289mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_agri_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
## 630.1010 1517.2736 577.3289 575.0241mva_m3 <- manyglm(formula = spp_no_rare ~ Distance_agri_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
## 630.1010 1517.2736 577.3289 575.0241 575.0241Transect and Site improve the model.
Many GLM with Site:
mva_m0.nb <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, family = "negative.binomial", data = env_data)
drop1(mva_m0.nb)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff + Site
## Df AIC
## <none> 4130.7
## cosinus_Monitoring.Time.Diff 7 4160.9
## sinus_Monitoring.Time.Diff 7 4168.7
## Site 28 4410.7
## Distance_agri_rescaled:rescaled_Time.Diff 7 4141.8Final model includes all factors.
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled + dunes*rescaled_Time.Diff + Distance_agri_rescaled:dunes + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + Site, family = "negative.binomial", data = env_data)
# drop1(mva_m2)# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled + dunes + rescaled_Time.Diff + Distance_agri_rescaled:dunes + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
# mva_m2 <- manyglm(formula = spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + Site, family = "negative.binomial", data = env_data)
# drop1(mva_m2)plot(mva_m0.nb, which=1:3)## Warning in default.plot.manyglm(x, res.type = res.type, which = which, caption
## = caption, : Only the first 7 colors will be used for plotting.
#summ_m2 <- summary(mva_m0.nb)
#saveRDS(summ_m2, "output/summ_model_arid_site.RDS")
summ_m2 <- readRDS("output/summ_model_arid_site.RDS")
#anov_m2.uni <- anova(mva_m0.nb, p.uni = "adjusted")
#saveRDS(anov_m2.uni, "output/anova_arid_site.RDS")
anov_m2.uni <- readRDS("output/anova_arid_site.RDS")
print(summ_m2)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 13.525 0.001 ***
## Distance_agri_rescaled 13.810 0.001 ***
## rescaled_Time.Diff 8.096 0.001 ***
## cosinus_Monitoring.Time.Diff 6.887 0.001 ***
## sinus_Monitoring.Time.Diff 6.383 0.001 ***
## SiteLotan 10.473 0.001 ***
## SiteParan 9.145 0.001 ***
## SiteZofar 13.015 0.001 ***
## SiteEin Yahav 12.449 0.001 ***
## Distance_agri_rescaled:rescaled_Time.Diff 5.569 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 24.74, 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).print(anov_m2.uni)## Analysis of Deviance Table
##
## Model: spp_no_rare ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 302
## Distance_agri_rescaled 301 1 82.2 0.001 ***
## rescaled_Time.Diff 300 1 77.9 0.002 **
## cosinus_Monitoring.Time.Diff 299 1 30.0 0.003 **
## sinus_Monitoring.Time.Diff 298 1 94.2 0.001 ***
## Site 294 4 332.2 0.001 ***
## Distance_agri_rescaled:rescaled_Time.Diff 293 1 25.1 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Canis.aureus Canis.lupus
## Dev Pr(>Dev) Dev
## (Intercept)
## Distance_agri_rescaled 12.305 0.037 2.561
## rescaled_Time.Diff 10.235 0.056 29.11
## cosinus_Monitoring.Time.Diff 0.459 0.559 4.057
## sinus_Monitoring.Time.Diff 0.755 0.780 2.012
## Site 72.196 0.001 8.393
## Distance_agri_rescaled:rescaled_Time.Diff 4.896 0.165 0
## Capra.nubiana
## Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_agri_rescaled 0.545 43.606 0.001
## rescaled_Time.Diff 0.002 12.357 0.024
## cosinus_Monitoring.Time.Diff 0.357 1.504 0.432
## sinus_Monitoring.Time.Diff 0.638 69.131 0.001
## Site 0.153 47.828 0.001
## Distance_agri_rescaled:rescaled_Time.Diff 0.994 7.384 0.057
## Gazella.dorcas Hyaena.hyaena
## Dev Pr(>Dev) Dev
## (Intercept)
## Distance_agri_rescaled 15.31 0.007 2.458
## rescaled_Time.Diff 24.287 0.003 0.03
## cosinus_Monitoring.Time.Diff 5.104 0.307 3.235
## sinus_Monitoring.Time.Diff 14.232 0.009 0.971
## Site 37.488 0.001 9.441
## Distance_agri_rescaled:rescaled_Time.Diff 8.543 0.038 0.457
## Lepus.capensis
## Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_agri_rescaled 0.545 0.121 0.818
## rescaled_Time.Diff 0.856 1.552 0.651
## cosinus_Monitoring.Time.Diff 0.425 3.041 0.425
## sinus_Monitoring.Time.Diff 0.780 0.33 0.780
## Site 0.153 78.365 0.001
## Distance_agri_rescaled:rescaled_Time.Diff 0.772 1.383 0.598
## Vulpes.vulpes
## Dev Pr(>Dev)
## (Intercept)
## Distance_agri_rescaled 5.813 0.210
## rescaled_Time.Diff 0.306 0.851
## cosinus_Monitoring.Time.Diff 12.64 0.023
## sinus_Monitoring.Time.Diff 6.817 0.096
## Site 78.519 0.001
## Distance_agri_rescaled:rescaled_Time.Diff 2.396 0.452
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.Multivariate test shows that all factors (Distance from agriculture, time difference and the interaction between the two, sin and cosin time difference and Site) significantly affect species composition
coefp <- merge(data.table(t(coef(mva_m0.nb)/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_agri.coef","Time_diff.coef","cosin_time.coef","sin_time.coef","Site_Lotan.coef","Site_Paran.coef","Site_Zofar.coef","Site_Ein Yahav.coef","Distance_agri_time.coef","Intercept.p","Distance_agri.p","Time_diff.p","cosin_time.p","sin_time.p","Site.p","Distance_agri_time.p","species_abundance")
write.csv(coefp, "coefficients_Arid_Desert.csv")Per-species models - for Canis aureus only used one Site (Yotvata) and for Capra nubiana only two sites (Yotvata and Lotan) because of many zeros
####Capra nubiana####
only_Capra_data <- P.anal[,c(1:26,29)]
Capra.nubiana.two.sites <- only_Capra_data[Site %in% c("Yotvata", "Lotan")]
glm.Capra.nubiana <- glm.nb(`Capra nubiana` ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, data = Capra.nubiana.two.sites, maxit = 999)
coef(glm.Capra.nubiana)## (Intercept)
## -5.3221834
## Distance_agri_rescaled
## 4.6121649
## rescaled_Time.Diff
## 0.1922945
## cosinus_Monitoring.Time.Diff
## 2.4648134
## sinus_Monitoring.Time.Diff
## 3.6125795
## SiteYotvata
## -3.1391983
## Distance_agri_rescaled:rescaled_Time.Diff
## -1.2799111coef(mva_m0.nb)## Canis.aureus Canis.lupus
## (Intercept) 0.69252453 -1.267224310
## Distance_agri_rescaled -1.18759646 0.502207397
## rescaled_Time.Diff 0.56347272 1.601455556
## cosinus_Monitoring.Time.Diff -0.31595816 0.556084286
## sinus_Monitoring.Time.Diff 0.06535263 0.908012936
## SiteLotan -3.29327281 -1.100984348
## SiteParan -5.11279350 0.095262226
## SiteZofar -4.01381142 -1.688544946
## SiteEin Yahav -5.43535222 -1.260133169
## Distance_agri_rescaled:rescaled_Time.Diff -0.83830144 0.004557944
## Capra.nubiana Gazella.dorcas
## (Intercept) -8.460181 -3.7674126
## Distance_agri_rescaled 4.611579 2.8143264
## rescaled_Time.Diff 0.192041 0.9843444
## cosinus_Monitoring.Time.Diff 2.464587 0.7685791
## sinus_Monitoring.Time.Diff 3.612071 -0.6974734
## SiteLotan 3.138868 3.0096643
## SiteParan -8.705339 3.5213536
## SiteZofar -6.446088 3.8564988
## SiteEin Yahav -11.167216 2.5766650
## Distance_agri_rescaled:rescaled_Time.Diff -1.279655 -0.9174263
## Hyaena.hyaena Lepus.capensis
## (Intercept) -1.3226630 -1.47294721
## Distance_agri_rescaled 0.2034286 -0.01452709
## rescaled_Time.Diff -0.1322795 0.46253783
## cosinus_Monitoring.Time.Diff 0.3764462 -0.06510093
## sinus_Monitoring.Time.Diff 0.1265160 -0.09094098
## SiteLotan -0.4513221 0.47870879
## SiteParan -0.2005448 1.04514367
## SiteZofar 0.6285472 2.04366525
## SiteEin Yahav -0.2965664 3.47451516
## Distance_agri_rescaled:rescaled_Time.Diff 0.1286649 0.27140401
## Vulpes.vulpes
## (Intercept) 1.4261256
## Distance_agri_rescaled -0.8380766
## rescaled_Time.Diff 0.1301859
## cosinus_Monitoring.Time.Diff -0.3238206
## sinus_Monitoring.Time.Diff 0.2567386
## SiteLotan -1.1567613
## SiteParan -1.2343164
## SiteZofar -3.2675189
## SiteEin Yahav -1.9164041
## Distance_agri_rescaled:rescaled_Time.Diff 0.2614349source("functions/plot_interaction_mammals.R")
plot_model_interaction(P.anal = Capra.nubiana.two.sites, m = glm.Capra.nubiana, eff2plot = "rescaled_Time.Diff", modvar2plot = "Distance_agri_rescaled", plot_points = FALSE, export_plot = TRUE, outpath = "output/arid/", fontname = "Almoni ML v5 AAA")
# capra_interact <- interact_plot(model = glm.Capra.nubiana, data=Capra.nubiana.two.sites, pred = rescaled_Time.Diff, modx = Distance_agri_rescaled, partial.residuals = F, modx.values = "plus-minus", jitter = c(0.2), point.size = 3, rug = F, point.shape = F, interval = T, colors = "Dark2") + scale_x_continuous(breaks=c(-1.233466888,-0.601840099,0.028922631,0.66054942,1.291312151), labels=c("June 2014", "June 2016", "June 2018", "June 2020", "June 2022"), name = "Time")
#Testing whether trends and means are different from zero
#Calculate mean and sd of distance from agriculture
sd_distance_agri <- sd(env_data$Distance_agri_rescaled)
print(sd_distance_agri)## [1] 0.7725932mean_distance_agri <- mean(env_data$Distance_agri_rescaled)
print(mean_distance_agri)## [1] 0.1973055#calculate +/1 1 SD
small_distance <- mean_distance_agri - sd_distance_agri
large_distance <- mean_distance_agri + sd_distance_agri
#Test whether temporal trend is different from zero
mm_rescaled_time <- emtrends(glm.Capra.nubiana, specs = "Distance_agri_rescaled", var = "rescaled_Time.Diff", type = "response", at = list(Distance_agri_rescaled = c(small_distance,large_distance)))
print(mm_rescaled_time)## Distance_agri_rescaled rescaled_Time.Diff.trend SE df asymp.LCL asymp.UCL
## -0.575 0.929 0.723 Inf -0.488 2.345
## 0.970 -1.049 0.187 Inf -1.416 -0.683
##
## 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_agri_rescaled rescaled_Time.Diff.trend SE df z.ratio p.value
## -0.575 0.929 0.723 Inf 1.285 0.1988
## 0.970 -1.049 0.187 Inf -5.609 <.0001
##
## 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 agriculture
pairwise_distance <- emmeans(object = glm.Capra.nubiana, ~Distance_agri_rescaled*rescaled_Time.Diff, at = list(Distance_agri_rescaled = c(small_distance,large_distance)))
print(pairwise_distance)## Distance_agri_rescaled rescaled_Time.Diff emmean SE df asymp.LCL asymp.UCL
## -0.575 0.281 -7.212 1.291 Inf -9.74 -4.6812
## 0.970 0.281 -0.642 0.355 Inf -1.34 0.0541
##
## 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_agri_rescaled-0.575287691303192) - Distance_agri_rescaled0.96989861358702
## rescaled_Time.Diff estimate SE df z.ratio p.value
## 0.281 -6.57 1.05 Inf -6.255 <.0001
##
## Results are averaged over the levels of: Site
## Results are given on the log (not the response) scale.####Gazella dorcas####
Gazella.dorcas <- spp_no_rare[,"Gazella.dorcas"]
glm.Gazella.dorcas <- glm.nb(Gazella.dorcas ~ Distance_agri_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, data = env_data, maxit = 999)
coef(glm.Gazella.dorcas)## (Intercept)
## -3.7672758
## Distance_agri_rescaled
## 2.8142456
## rescaled_Time.Diff
## 0.9843487
## cosinus_Monitoring.Time.Diff
## 0.7685811
## sinus_Monitoring.Time.Diff
## -0.6974723
## SiteLotan
## 3.0095339
## SiteParan
## 3.5212432
## SiteZofar
## 3.8563009
## SiteEin Yahav
## 2.5765897
## Distance_agri_rescaled:rescaled_Time.Diff
## -0.9174520coef(mva_m0.nb)## Canis.aureus Canis.lupus
## (Intercept) 0.69252453 -1.267224310
## Distance_agri_rescaled -1.18759646 0.502207397
## rescaled_Time.Diff 0.56347272 1.601455556
## cosinus_Monitoring.Time.Diff -0.31595816 0.556084286
## sinus_Monitoring.Time.Diff 0.06535263 0.908012936
## SiteLotan -3.29327281 -1.100984348
## SiteParan -5.11279350 0.095262226
## SiteZofar -4.01381142 -1.688544946
## SiteEin Yahav -5.43535222 -1.260133169
## Distance_agri_rescaled:rescaled_Time.Diff -0.83830144 0.004557944
## Capra.nubiana Gazella.dorcas
## (Intercept) -8.460181 -3.7674126
## Distance_agri_rescaled 4.611579 2.8143264
## rescaled_Time.Diff 0.192041 0.9843444
## cosinus_Monitoring.Time.Diff 2.464587 0.7685791
## sinus_Monitoring.Time.Diff 3.612071 -0.6974734
## SiteLotan 3.138868 3.0096643
## SiteParan -8.705339 3.5213536
## SiteZofar -6.446088 3.8564988
## SiteEin Yahav -11.167216 2.5766650
## Distance_agri_rescaled:rescaled_Time.Diff -1.279655 -0.9174263
## Hyaena.hyaena Lepus.capensis
## (Intercept) -1.3226630 -1.47294721
## Distance_agri_rescaled 0.2034286 -0.01452709
## rescaled_Time.Diff -0.1322795 0.46253783
## cosinus_Monitoring.Time.Diff 0.3764462 -0.06510093
## sinus_Monitoring.Time.Diff 0.1265160 -0.09094098
## SiteLotan -0.4513221 0.47870879
## SiteParan -0.2005448 1.04514367
## SiteZofar 0.6285472 2.04366525
## SiteEin Yahav -0.2965664 3.47451516
## Distance_agri_rescaled:rescaled_Time.Diff 0.1286649 0.27140401
## Vulpes.vulpes
## (Intercept) 1.4261256
## Distance_agri_rescaled -0.8380766
## rescaled_Time.Diff 0.1301859
## cosinus_Monitoring.Time.Diff -0.3238206
## sinus_Monitoring.Time.Diff 0.2567386
## SiteLotan -1.1567613
## SiteParan -1.2343164
## SiteZofar -3.2675189
## SiteEin Yahav -1.9164041
## Distance_agri_rescaled:rescaled_Time.Diff 0.2614349plot_model_interaction(P.anal = env_data, m = glm.Gazella.dorcas, eff2plot = "rescaled_Time.Diff", modvar2plot = "Distance_agri_rescaled", plot_points = FALSE, export_plot = TRUE, outpath = "output/arid/", fontname = "Almoni ML v5 AAA", legend_position = "bottom")
# interact_plot(model = glm.Gazella.dorcas, data=env_data, pred = rescaled_Time.Diff, modx = Distance_agri_rescaled, partial.residuals = F, modx.values = "plus-minus", jitter = c(0.2), point.size = 3, rug = F, point.shape = F, interval = T, colors = "Dark2") + scale_x_continuous(breaks=c(-1.233466888,-0.601840099,0.028922631,0.66054942,1.291312151), labels=c("June 2014", "June 2016", "June 2018", "June 2020", "June 2022"), name = "Time")
#Test whether temporal trend is different from zero
mm_rescaled_time <- emtrends(glm.Gazella.dorcas, specs = "Distance_agri_rescaled", var = "rescaled_Time.Diff", type = "response", at = list(Distance_agri_rescaled = c(small_distance,large_distance)))
print(mm_rescaled_time)## Distance_agri_rescaled rescaled_Time.Diff.trend SE df asymp.LCL asymp.UCL
## -0.575 1.5121 0.316 Inf 0.892 2.132
## 0.970 0.0945 0.214 Inf -0.325 0.514
##
## 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_agri_rescaled rescaled_Time.Diff.trend SE df z.ratio p.value
## -0.575 1.5121 0.316 Inf 4.780 <.0001
## 0.970 0.0945 0.214 Inf 0.441 0.6591
##
## 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 agriculture
pairwise_distance <- emmeans(object = glm.Gazella.dorcas, ~Distance_agri_rescaled*rescaled_Time.Diff, at = list(Distance_agri_rescaled = c(small_distance,large_distance)))
print(pairwise_distance)## Distance_agri_rescaled rescaled_Time.Diff emmean SE df asymp.LCL asymp.UCL
## -0.575 0.277 -2.27 0.314 Inf -2.88 -1.65
## 0.970 0.277 1.69 0.172 Inf 1.35 2.03
##
## 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_agri_rescaled-0.575287691303192) - Distance_agri_rescaled0.96989861358702
## rescaled_Time.Diff estimate SE df z.ratio p.value
## 0.277 -3.96 0.374 Inf -10.570 <.0001
##
## Results are averaged over the levels of: Site
## Results are given on the log (not the response) scale.Capra:
Temporal trend is significantly different from zero in large distances from agriculture (p < 0.0001) and not significantly different from zero in small distances (p = 0.19)
On average over time, ibex abundance is significantly higher near agriculture (p < 0.0001)
Gazella dorcas:
Temporal trend is significantly different from zero in small distances from agriculture (p < 0.0001) and not significantly different from zero in large distances (p = 0.65)
On average over time, Gazelle abundance is significantly higher near agriculture (p < 0.0001)
####Canis aureus####
only_canis_aur_data <- P.anal[,c(1:27)]
Canis.aureus.only.yotvata <- only_canis_aur_data[Site %in% c("Yotvata")]
Canis.aureus <- spp_no_rare[,"Canis.aureus"]
glm.Canis.aureus <- glm.nb(`Canis aureus` ~ Distance_agri_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, data = Canis.aureus.only.yotvata, maxit = 999)
#coef(glm.Canis.aureus)
#coef(mva_m0.nb)
source("functions/plot_model_effect_mammals.R")
plot_model_effect(P.anal = Canis.aureus.only.yotvata, m = glm.Canis.aureus, 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/arid/")
# effect_plot(model = glm.Canis.aureus, data=Canis.aureus.only.yotvata, pred = rescaled_Time.Diff, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1") + ylim(0,5) + scale_x_continuous(breaks=c(-1.233466888,-0.601840099,0.028922631,0.66054942,1.291312151), labels=c("June 2014", "June 2016", "June 2018", "June 2020", "June 2022"), name = "Time")
plot_model_effect(P.anal = Canis.aureus.only.yotvata, m = glm.Canis.aureus, eff2plot = "Distance_agri_rescaled", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "output/arid/")
# effect_plot(model = glm.Canis.aureus, data=Canis.aureus.only.yotvata, pred = Distance_agri_rescaled, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1") + scale_x_continuous(breaks=c(-0.589201946,0.245068396,1.172035444,2.099002491), labels=c("100", "1000", "2000", "3000"), name = "Distance (m)")
#
####Canis lupus####
Canis.lupus <- spp_no_rare[,"Canis.lupus"]
glm.Canis.lupus <- glm.nb(Canis.lupus ~ Distance_agri_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + Site, data = env_data, maxit = 999)
coef(glm.Canis.lupus)## (Intercept) Distance_agri_rescaled
## -1.26773275 0.50367984
## rescaled_Time.Diff cosinus_Monitoring.Time.Diff
## 1.60219090 0.55545107
## sinus_Monitoring.Time.Diff SiteLotan
## 0.90840035 -1.10096642
## SiteParan SiteZofar
## 0.09479091 -1.68726381
## SiteEin Yahav
## -1.26068664coef(mva_m0.nb)## Canis.aureus Canis.lupus
## (Intercept) 0.69252453 -1.267224310
## Distance_agri_rescaled -1.18759646 0.502207397
## rescaled_Time.Diff 0.56347272 1.601455556
## cosinus_Monitoring.Time.Diff -0.31595816 0.556084286
## sinus_Monitoring.Time.Diff 0.06535263 0.908012936
## SiteLotan -3.29327281 -1.100984348
## SiteParan -5.11279350 0.095262226
## SiteZofar -4.01381142 -1.688544946
## SiteEin Yahav -5.43535222 -1.260133169
## Distance_agri_rescaled:rescaled_Time.Diff -0.83830144 0.004557944
## Capra.nubiana Gazella.dorcas
## (Intercept) -8.460181 -3.7674126
## Distance_agri_rescaled 4.611579 2.8143264
## rescaled_Time.Diff 0.192041 0.9843444
## cosinus_Monitoring.Time.Diff 2.464587 0.7685791
## sinus_Monitoring.Time.Diff 3.612071 -0.6974734
## SiteLotan 3.138868 3.0096643
## SiteParan -8.705339 3.5213536
## SiteZofar -6.446088 3.8564988
## SiteEin Yahav -11.167216 2.5766650
## Distance_agri_rescaled:rescaled_Time.Diff -1.279655 -0.9174263
## Hyaena.hyaena Lepus.capensis
## (Intercept) -1.3226630 -1.47294721
## Distance_agri_rescaled 0.2034286 -0.01452709
## rescaled_Time.Diff -0.1322795 0.46253783
## cosinus_Monitoring.Time.Diff 0.3764462 -0.06510093
## sinus_Monitoring.Time.Diff 0.1265160 -0.09094098
## SiteLotan -0.4513221 0.47870879
## SiteParan -0.2005448 1.04514367
## SiteZofar 0.6285472 2.04366525
## SiteEin Yahav -0.2965664 3.47451516
## Distance_agri_rescaled:rescaled_Time.Diff 0.1286649 0.27140401
## Vulpes.vulpes
## (Intercept) 1.4261256
## Distance_agri_rescaled -0.8380766
## rescaled_Time.Diff 0.1301859
## cosinus_Monitoring.Time.Diff -0.3238206
## sinus_Monitoring.Time.Diff 0.2567386
## SiteLotan -1.1567613
## SiteParan -1.2343164
## SiteZofar -3.2675189
## SiteEin Yahav -1.9164041
## Distance_agri_rescaled:rescaled_Time.Diff 0.2614349plot(glm.Canis.lupus)
plot_model_effect(P.anal = env_data, m = glm.Canis.lupus, 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/arid/")
# effect_plot(model = glm.Canis.lupus, data=env_data, pred = rescaled_Time.Diff, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1") + scale_x_continuous(breaks=c(-1.233466888,-0.601840099,0.028922631,0.66054942,1.291312151), labels=c("June 2014", "June 2016", "June 2018", "June 2020", "June 2022"), name = "Time")n = 303 (total number of cameras)
total abundance = 3035 (10 species)
total abundance without rare species = 2933 (7 species)
Session information
session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.3 (2023-03-15 ucrt)
## os Windows 10 x64 (build 22631)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate Hebrew_Israel.utf8
## ctype Hebrew_Israel.utf8
## tz Asia/Jerusalem
## date 2024-04-09
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-7 2017-09-03 [1] R-Forge (R 4.2.3)
## betareg 3.2-0 2021-02-09 [1] R-Forge (R 4.2.3)
## boot 1.3-28.1 2022-11-22 [1] CRAN (R 4.2.3)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.3)
## cachem 1.0.7 2023-02-24 [1] CRAN (R 4.2.3)
## Cairo * 1.6-0 2022-07-05 [1] CRAN (R 4.2.2)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.3)
## car * 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
## carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.3)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.2.3)
## chron * 2.3-61 2023-05-02 [1] CRAN (R 4.2.3)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3)
## cluster 2.1.4 2022-08-22 [1] CRAN (R 4.2.3)
## coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.3)
## codetools 0.2-19 2023-02-01 [1] CRAN (R 4.2.2)
## colorspace 2.1-1 2023-03-08 [1] R-Forge (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.3)
## curl 5.0.0 2023-01-12 [1] CRAN (R 4.2.3)
## data.table * 1.14.8 2023-02-17 [1] CRAN (R 4.2.3)
## devtools * 2.4.5 2022-10-11 [1] CRAN (R 4.2.3)
## digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.3)
## doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.2.3)
## dplyr * 1.1.1 2023-03-22 [1] CRAN (R 4.2.3)
## ecoCopula * 1.0.2 2022-03-02 [1] CRAN (R 4.2.3)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.3)
## emmeans * 1.8.6 2023-05-11 [1] CRAN (R 4.2.3)
## estimability 1.4.1 2022-08-05 [1] CRAN (R 4.2.1)
## evaluate 0.20 2023-01-17 [1] CRAN (R 4.2.3)
## extrafont * 0.19 2023-01-18 [1] CRAN (R 4.2.2)
## extrafontdb 1.0 2012-06-11 [1] CRAN (R 4.2.0)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.3)
## farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.3)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.3)
## flexmix 2.3-19 2023-03-16 [1] CRAN (R 4.2.3)
## foreach 1.5.2 2022-02-02 [1] CRAN (R 4.2.3)
## Formula 1.2-6 2023-02-25 [1] R-Forge (R 4.2.2)
## fs 1.6.1 2023-02-06 [1] CRAN (R 4.2.3)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.3)
## ggplot2 * 3.5.0 2024-02-23 [1] CRAN (R 4.2.3)
## glm2 1.2.1 2018-08-11 [1] CRAN (R 4.2.0)
## glmmTMB * 1.1.8 2023-10-07 [1] CRAN (R 4.2.3)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.3)
## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.3)
## highr 0.10 2022-12-22 [1] CRAN (R 4.2.3)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.3)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.3)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.3)
## httr 1.4.5 2023-02-24 [1] CRAN (R 4.2.3)
## insight 0.19.9 2024-03-15 [1] CRAN (R 4.2.3)
## interactions * 1.1.5 2021-07-02 [1] CRAN (R 4.2.3)
## iterators 1.0.14 2022-02-05 [1] CRAN (R 4.2.3)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.3)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.3)
## jtools * 2.2.1 2022-12-02 [1] CRAN (R 4.2.3)
## kableExtra * 1.4.0 2024-01-24 [1] CRAN (R 4.2.3)
## knitr 1.42 2023-01-25 [1] CRAN (R 4.2.3)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.2.0)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.3)
## lattice * 0.21-8 2023-04-05 [1] CRAN (R 4.2.3)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.3)
## lme4 * 1.1-32 2023-03-14 [1] CRAN (R 4.2.3)
## lmtest 0.9-40 2022-03-21 [1] CRAN (R 4.2.3)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.3)
## MASS * 7.3-58.3 2023-03-07 [1] CRAN (R 4.2.3)
## Matrix * 1.5-5 2023-04-05 [1] R-Forge (R 4.2.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.3)
## mgcv 1.8-42 2023-03-02 [1] CRAN (R 4.2.3)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.3)
## minqa 1.2.5 2022-10-19 [1] CRAN (R 4.2.3)
## modeltools 0.2-23 2020-03-05 [1] CRAN (R 4.2.0)
## multcomp 1.4-23 2023-03-09 [1] CRAN (R 4.2.3)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.3)
## mvabund * 4.2.1 2022-02-16 [1] CRAN (R 4.2.3)
## mvtnorm 1.2-0 2023-04-05 [1] R-Forge (R 4.2.3)
## nlme 3.1-162 2023-01-31 [1] CRAN (R 4.2.3)
## nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.2.3)
## nnet 7.3-18 2022-09-28 [1] CRAN (R 4.2.3)
## numDeriv 2022.9-1 2022-09-27 [1] R-Forge (R 4.2.1)
## ordinal 2022.11-16 2022-11-16 [1] CRAN (R 4.2.3)
## pander 0.6.5 2022-03-18 [1] CRAN (R 4.2.3)
## performance * 0.10.3 2023-04-07 [1] CRAN (R 4.2.3)
## permute * 0.9-7 2022-01-27 [1] CRAN (R 4.2.3)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
## pkgbuild 1.4.2.9000 2023-07-11 [1] Github (r-lib/pkgbuild@7048654)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.3)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.3)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.3)
## processx 3.8.0 2022-10-26 [1] CRAN (R 4.2.3)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.3)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.3)
## ps 1.7.4 2023-04-02 [1] CRAN (R 4.2.3)
## purrr 1.0.1 2023-01-10 [1] CRAN (R 4.2.3)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.3)
## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.2.0)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.3)
## readxl * 1.4.2 2023-02-09 [1] CRAN (R 4.2.3)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.3)
## rlang * 1.1.0 2023-03-14 [1] CRAN (R 4.2.3)
## rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.3)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.3)
## Rttf2pt1 1.3.12 2023-01-22 [1] CRAN (R 4.2.2)
## sandwich 3.1-0 2023-04-04 [1] R-Forge (R 4.2.3)
## sass 0.4.5 2023-01-24 [1] CRAN (R 4.2.3)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.2.3)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.3)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.3)
## statmod 1.5.0 2023-01-06 [1] CRAN (R 4.2.3)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.2)
## stringr 1.5.0 2022-12-02 [1] CRAN (R 4.2.3)
## survival 3.5-5 2023-03-12 [1] CRAN (R 4.2.3)
## svglite 2.1.1 2023-01-10 [1] CRAN (R 4.2.3)
## systemfonts 1.0.4 2022-02-11 [1] CRAN (R 4.2.3)
## TH.data 1.1-2 2022-11-07 [1] R-Forge (R 4.2.3)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.3)
## TMB 1.9.6 2023-08-11 [1] CRAN (R 4.2.3)
## tweedie 2.3.5 2022-08-17 [1] CRAN (R 4.2.3)
## ucminf 1.1-4.1 2022-09-29 [1] CRAN (R 4.2.1)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.3)
## usethis * 2.1.6 2022-05-25 [1] CRAN (R 4.2.3)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.3)
## vctrs 0.6.1 2023-03-22 [1] CRAN (R 4.2.3)
## vegan * 2.6-4 2022-10-11 [1] CRAN (R 4.2.3)
## viridisLite 0.4.1 2022-08-22 [1] CRAN (R 4.2.3)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.3)
## xfun 0.38 2023-03-24 [1] CRAN (R 4.2.3)
## xml2 1.3.3 2021-11-30 [1] CRAN (R 4.2.3)
## xtable 1.8-6 2020-06-19 [1] R-Forge (R 4.2.3)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.3)
## zoo 1.8-11 2022-09-17 [1] CRAN (R 4.2.3)
##
## [1] C:/Users/Ron Chen/AppData/Local/R/win-library/4.2.3
## [2] C:/Program Files/R/R-4.2.3/library
##
## ──────────────────────────────────────────────────────────────────────────────