Mammal data - Negev Highlands
Talia Gabay & Ron Chen
2024-04-09
P <- P_byplot[grepl("Negev Highlands",Unit),][order(Cycle_number,Site)][,Site:=factor(Site)]Negev Highlands - Richness by camera
Monitoring started in 2014 (5 monitoring cycles), Each site as two types of plots - near or far from Settlements. There are 5 sites total - Bislach, Ezuz, Merhav Am, Sde Boker, Yeruham Each has 5 sampling points except for Sde Boker(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$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$Site <- as.factor(P.anal$Site)
print("Site has 5 levels")## [1] "Site has 5 levels"print(levels(P.anal$Site))## [1] "Bislach" "Ezuz" "Merhav Am" "Sde Boker" "Yeruham"P.anal$Transect <- as.factor(P.anal$Transect)
print("Transect has 10 levels")## [1] "Transect has 10 levels"print(levels(P.anal$Transect))## [1] "Bislach Far" "Bislach Near" "Ezuz Far" "Ezuz Near"
## [5] "Merhav Am Far" "Merhav Am Near" "Sde Boker Far" "Sde Boker Near"
## [9] "Yeruham Far" "Yeruham 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] "Bislach Far_02_07_2014" "Bislach Far_04_01_2019"
## [3] "Bislach Far_12_10_2022" "Bislach Far_16_10_2020"
## [5] "Bislach Far_19_10_2016" "Bislach Near_02_07_2014"
## [7] "Bislach Near_04_01_2019" "Bislach Near_12_10_2022"
## [9] "Bislach Near_16_10_2020" "Bislach Near_19_10_2016"
## [11] "Ezuz Far_01_02_2019" "Ezuz Far_06_11_2022"
## [13] "Ezuz Far_09_10_2020" "Ezuz Far_20_09_2016"
## [15] "Ezuz Far_26_06_2014" "Ezuz Near_01_02_2019"
## [17] "Ezuz Near_06_11_2022" "Ezuz Near_09_10_2020"
## [19] "Ezuz Near_20_09_2016" "Ezuz Near_26_06_2014"
## [21] "Merhav Am Far_05_10_2020" "Merhav Am Far_19_06_2014"
## [23] "Merhav Am Far_20_01_2019" "Merhav Am Far_23_06_2016"
## [25] "Merhav Am Far_29_09_2022" "Merhav Am Near_05_10_2020"
## [27] "Merhav Am Near_19_06_2014" "Merhav Am Near_20_01_2019"
## [29] "Merhav Am Near_23_06_2016" "Merhav Am Near_29_09_2022"
## [31] "Sde Boker Far_11_08_2016" "Sde Boker Far_18_09_2022"
## [33] "Sde Boker Far_20_10_2020" "Sde Boker Far_29_06_2014"
## [35] "Sde Boker Far_30_12_2018" "Sde Boker Near_11_08_2016"
## [37] "Sde Boker Near_18_09_2022" "Sde Boker Near_20_10_2020"
## [39] "Sde Boker Near_30_12_2018" "Yeruham Far_11_12_2018"
## [41] "Yeruham Far_14_07_2014" "Yeruham Far_17_09_2020"
## [43] "Yeruham Far_24_08_2016" "Yeruham Far_31_10_2022"
## [45] "Yeruham Near_11_12_2018" "Yeruham Near_14_07_2014"
## [47] "Yeruham Near_17_09_2020" "Yeruham Near_24_08_2016"
## [49] "Yeruham Near_31_10_2022"#distance.rescaled <- scale(all.mammal.abundances$Distance)
#write.csv(distance.rescaled, file = "Distance_rescaled.csv")print("RICHNESS WITH RARE SPECIES")## [1] "RICHNESS WITH RARE SPECIES"plot_alpha_diversity(P, x_val = "Settlements", y_val = "richness", ylab_val = "richness", xlab_val = "Settlements", fill_val = "Settlements")
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.NegevHighlands <- all.mammal.abundances[grepl("Negev Highlands",Unit),]
spp <- abu_by_spp.NegevHighlands[,27:40]
#spp <- abu_by_spp.NegevHighlands[,25:38]
# filter out species with zero counts
print(colSums(spp))## Canis aureus Canis lupus Capra nubiana Dama mesopotamica
## 973 4 11 0
## Equus hemionus Gazella dorcas Gazella gazella Hyaena hyaena
## 16 190 0 128
## Hystrix indica Lepus capensis Meles meles Oryx leucoryx
## 607 177 0 0
## Sus scrofa Vulpes vulpes
## 0 465spp <- spp[,.SD,.SDcols = colSums(spp)>0]
spp <- mvabund(spp)
plot(spp ~ P.anal$Settlements, 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 = "Settlements",
groups = P.anal$Settlements, gcolor = my_cols,
color = my_cols[P.anal$Settlements],
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, Distance_rescaled, rescaled_Time.Diff, cosinus_Monitoring.Time.Diff, sinus_Monitoring.Time.Diff, Site, Transect, Transect_with_date)]))| richness | abundance | Distance_rescaled | rescaled_Time.Diff | cosinus_Monitoring.Time.Diff | sinus_Monitoring.Time.Diff | Site | Transect | Transect_with_date | |
|---|---|---|---|---|---|---|---|---|---|
| Min. :1.000 | Min. : 1.000 | Min. :-0.91363 | Min. :-1.2179 | Min. :-0.9990 | Min. :-0.99999 | Bislach :53 | Yeruham Near : 37 | Bislach Near_12_10_2022 : 9 | |
| 1st Qu.:1.000 | 1st Qu.: 2.000 | 1st Qu.:-0.63636 | 1st Qu.:-0.4809 | 1st Qu.:-0.4161 | 1st Qu.:-0.50637 | Ezuz :54 | Sde Boker Near: 34 | Sde Boker Near_18_09_2022: 9 | |
| Median :2.000 | Median : 5.000 | Median : 0.06859 | Median : 0.7539 | Median : 0.3590 | Median :-0.09718 | Merhav Am:59 | Yeruham Far : 32 | Sde Boker Near_20_10_2020: 9 | |
| Mean :2.235 | Mean : 8.628 | Mean :-0.01456 | Mean : 0.4169 | Mean : 0.1602 | Mean :-0.01436 | Sde Boker:63 | Bislach Far : 30 | Sde Boker Near_30_12_2018: 9 | |
| 3rd Qu.:3.000 | 3rd Qu.:11.000 | 3rd Qu.: 0.49839 | 3rd Qu.: 1.3855 | 3rd Qu.: 0.8623 | 3rd Qu.: 0.67023 | Yeruham :69 | Merhav Am Near: 30 | Yeruham Far_31_10_2022 : 9 | |
| Max. :6.000 | Max. :70.000 | Max. : 1.12569 | Max. : 1.4278 | Max. : 0.9994 | Max. : 0.99882 | NA | Ezuz Far : 29 | Yeruham Near_11_12_2018 : 9 | |
| NA | NA | NA | NA | NA | NA | NA | (Other) :106 | (Other) :244 |
pairs(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("Distance_rescaled","rescaled_Time.Diff","cosinus_Monitoring.Time.Diff","sinus_Monitoring.Time.Diff","Site","Transect","Transect_with_date")])
kable(cor(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=c("Distance_rescaled","rescaled_Time.Diff","cosinus_Monitoring.Time.Diff","sinus_Monitoring.Time.Diff","Site","Transect","Transect_with_date")]))| Distance_rescaled | rescaled_Time.Diff | cosinus_Monitoring.Time.Diff | sinus_Monitoring.Time.Diff | Site | Transect | Transect_with_date | |
|---|---|---|---|---|---|---|---|
| Distance_rescaled | 1.0000000 | 0.0760419 | -0.0276814 | 0.0574501 | -0.1913315 | -0.3192192 | -0.3191631 |
| rescaled_Time.Diff | 0.0760419 | 1.0000000 | 0.4169221 | -0.0144697 | 0.0116821 | 0.0229156 | 0.0334149 |
| cosinus_Monitoring.Time.Diff | -0.0276814 | 0.4169221 | 1.0000000 | 0.0923914 | 0.0500147 | 0.0434014 | 0.0495119 |
| sinus_Monitoring.Time.Diff | 0.0574501 | -0.0144697 | 0.0923914 | 1.0000000 | 0.2655165 | 0.2569541 | 0.2716732 |
| Site | -0.1913315 | 0.0116821 | 0.0500147 | 0.2655165 | 1.0000000 | 0.9853042 | 0.9806464 |
| Transect | -0.3192192 | 0.0229156 | 0.0434014 | 0.2569541 | 0.9853042 | 1.0000000 | 0.9952240 |
| Transect_with_date | -0.3191631 | 0.0334149 | 0.0495119 | 0.2716732 | 0.9806464 | 0.9952240 | 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_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.451806633794825"Overdispersion parameter is <1 (underdispersion) and therefore will use Poisson
m0 <- glm(data = P.anal, formula = richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff, family = poisson)
me0transectwithdate <- glmer(data = P.anal, formula = richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Transect_with_date), family = poisson, verbose = 1)## start par. = 1 fn = 1009.829
## At return
## eval: 15 fn: 922.62887 par: 0.00000
## (NM) 20: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 40: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 60: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 80: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 100: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 120: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 140: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 160: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 180: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 200: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 220: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238## boundary (singular) fit: see help('isSingular')me0transect <- glmer(data = P.anal, formula = richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Transect), family = poisson, verbose = 1)## start par. = 1 fn = 959.4909
## At return
## eval: 14 fn: 922.62887 par: 0.00000
## (NM) 20: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 40: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 60: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 80: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 100: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 120: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 140: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 160: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 180: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238## boundary (singular) fit: see help('isSingular')me0site <- glmer(data = P.anal, formula = richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + sinus_Monitoring.Time.Diff + (1|Site), family = poisson, verbose = 1)## start par. = 1 fn = 945.2112
## At return
## eval: 14 fn: 922.62887 par: 0.00000
## (NM) 20: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 40: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 60: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 80: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 100: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 120: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 140: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 160: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238
## (NM) 180: f = 922.629 at 0 0.659477 -0.0301886 0.262144 0.0522413 0.00395857 -0.0145238## boundary (singular) fit: see help('isSingular')All GLMERs result in singular fit, will use regular GLM:
m1 <- step(m0)## Start: AIC=934.63
## richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - sinus_Monitoring.Time.Diff 1 149.68 932.63
## - Distance_rescaled:rescaled_Time.Diff 1 149.71 932.66
## - cosinus_Monitoring.Time.Diff 1 150.41 933.36
## <none> 149.67 934.63
##
## Step: AIC=932.63
## richness ~ Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## Distance_rescaled:rescaled_Time.Diff
##
## Df Deviance AIC
## - Distance_rescaled:rescaled_Time.Diff 1 149.71 930.67
## - cosinus_Monitoring.Time.Diff 1 150.43 931.38
## <none> 149.68 932.63
##
## Step: AIC=930.67
## richness ~ Distance_rescaled + rescaled_Time.Diff + cosinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - Distance_rescaled 1 150.11 929.06
## - cosinus_Monitoring.Time.Diff 1 150.45 929.41
## <none> 149.71 930.67
## - rescaled_Time.Diff 1 172.77 951.72
##
## Step: AIC=929.06
## richness ~ rescaled_Time.Diff + cosinus_Monitoring.Time.Diff
##
## Df Deviance AIC
## - cosinus_Monitoring.Time.Diff 1 150.93 927.89
## <none> 150.11 929.06
## - rescaled_Time.Diff 1 172.80 949.75
##
## Step: AIC=927.89
## richness ~ rescaled_Time.Diff
##
## Df Deviance AIC
## <none> 150.93 927.89
## - rescaled_Time.Diff 1 184.96 959.91# m1 <- glmer(data = P.anal, formula = richness ~ Distance_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff + (1|Transect_with_date), family = poisson)
# drop1(m1)Final model includes time difference
summ(m1, exp = TRUE, digits=3)| Observations | 298 |
| Dependent variable | richness |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| χ²(1) | 34.027 |
| Pseudo-R² (Cragg-Uhler) | 0.112 |
| Pseudo-R² (McFadden) | 0.036 |
| AIC | 927.886 |
| BIC | 935.280 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.933 | 1.756 | 2.128 | 13.469 | 0.000 |
| rescaled_Time.Diff | 1.328 | 1.205 | 1.464 | 5.702 | 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)")
#m_coef <- coef(m2)
#P.anal$Cycle_number <- as.numeric(P.anal$Cycle_number)
#plot(P.anal$Cycle_number,
# exp(m_coef["rescaled_Time.Diff"]*P.anal$rescaled_Time.Diff),
#xlab = "Cycle number", ylab = "contribution to richness")Time difference significantly affects species richness
community analysis using package MVabund
Table of abundances per species by camera ID
abundance_and_cameras <- as.data.table(abu_by_spp.NegevHighlands[,c(5,27:40)])
dim(abundance_and_cameras)## [1] 298 15canisau <- abundance_and_cameras[,2] > 0
sum(canisau)## [1] 163canislup <- abundance_and_cameras[,3] > 0
sum(canislup)## [1] 3capra <- abundance_and_cameras[,4] > 0
sum(capra)## [1] 4Doma <- abundance_and_cameras[,5] > 0
sum(Doma)## [1] 0Equus <- abundance_and_cameras[,6] > 0
sum(Equus)## [1] 8Gazellador <- abundance_and_cameras[,7] > 0
sum(Gazellador)## [1] 56Gazellagaz <- abundance_and_cameras[,8] > 0
sum(Gazellagaz)## [1] 0Hyaena <- abundance_and_cameras[,9] > 0
sum(Hyaena)## [1] 74Hystrix <- abundance_and_cameras[,10] > 0
sum(Hystrix)## [1] 138Lepus <- abundance_and_cameras[,11] > 0
sum(Lepus)## [1] 72Meles <- abundance_and_cameras[,12] > 0
sum(Meles)## [1] 0Oryx <- abundance_and_cameras[,13] > 0
sum(Oryx)## [1] 0Sus <- abundance_and_cameras[,14] > 0
sum(Sus)## [1] 0Vulpus <- abundance_and_cameras[,15] > 0
sum(Vulpus)## [1] 148298 cameras
Species with a minimum of 20 cameras - Canis aureus, Gazella Dorcas, Hyaena hyaena, Hystrix indica, Lepus capensis, Vulpus vulpus
spp_no_rare <- abu_by_spp.NegevHighlands[,27:40]
print(colSums(spp_no_rare))## Canis aureus Canis lupus Capra nubiana Dama mesopotamica
## 973 4 11 0
## Equus hemionus Gazella dorcas Gazella gazella Hyaena hyaena
## 16 190 0 128
## Hystrix indica Lepus capensis Meles meles Oryx leucoryx
## 607 177 0 0
## Sus scrofa Vulpes vulpes
## 0 465# 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)
print(colSums(spp_no_rare))## Canis.aureus Gazella.dorcas Hyaena.hyaena Hystrix.indica Lepus.capensis
## 973 190 128 607 177
## Vulpes.vulpes
## 465Left with 6 species after filtering
env_data <- abu_by_spp.NegevHighlands[,..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
## 791.4671 1237.0394print("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 6 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 6 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_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 + Transect_with_date, data = env_data)
aic_dt$nb1 <- mva_m1$aic
colMeans(aic_dt)## nb po nb1
## 791.4671 1237.0394 761.0316mva_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
## 791.4671 1237.0394 761.0316 765.2318mva_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
## 791.4671 1237.0394 761.0316 765.2318 765.2318Transect + date, Transect and Site improve the model.
Many GLM with Site:
mva_m0.nb <- manyglm(formula = spp_no_rare ~ Distance_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_rescaled * rescaled_Time.Diff + cosinus_Monitoring.Time.Diff +
## sinus_Monitoring.Time.Diff + Site
## Df AIC
## <none> 4627.8
## cosinus_Monitoring.Time.Diff 6 4625.3
## sinus_Monitoring.Time.Diff 6 4622.1
## Site 24 4748.8
## Distance_rescaled:rescaled_Time.Diff 6 4637.2Dropped 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 + Site, 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 +
## Site
## Df AIC
## <none> 4622.1
## cosinus_Monitoring.Time.Diff 6 4619.1
## Site 24 4754.1
## Distance_rescaled:rescaled_Time.Diff 6 4630.4Dropped cosin monitoring 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 + Site, family = "negative.binomial", data = env_data)
drop1(mva_m2)## Single term deletions
##
## Model:
## spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Site
## Df AIC
## <none> 4619.1
## Site 24 4758.6
## Distance_rescaled:rescaled_Time.Diff 6 4627.3Final model includes Site and the interaction between distance from settlement and time difference
# 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 + Site, family = "negative.binomial", data = env_data)
# drop1(mva_m2)plot(mva_m2, which=1:3)## Warning in default.plot.manyglm(x, res.type = res.type, which = which, caption
## = caption, : Only the first 6 colors will be used for plotting.
# summ_m2 <- summary(mva_m2)
# saveRDS(summ_m2, file = "output/summ_negev_highlands_999iter.rds")
summ_m2 <- readRDS(file = "output/summ_negev_highlands_999iter.rds")
print(summ_m2)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 5.829 0.001 ***
## Distance_rescaled 6.848 0.001 ***
## rescaled_Time.Diff 8.470 0.001 ***
## SiteEzuz 5.615 0.001 ***
## SiteMerhav Am 5.633 0.002 **
## SiteSde Boker 5.499 0.001 ***
## SiteYeruham 5.695 0.001 ***
## Distance_rescaled:rescaled_Time.Diff 4.522 0.010 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 18.35, 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, file = "output/anova_negev_highlands_999iter.rds")
anov_m2.uni <- readRDS(file = "output/anova_negev_highlands_999iter.rds")
print(anov_m2.uni)## Analysis of Deviance Table
##
## Model: spp_no_rare ~ Distance_rescaled * rescaled_Time.Diff + Site
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 297
## Distance_rescaled 296 1 97.14 0.001 ***
## rescaled_Time.Diff 295 1 56.62 0.001 ***
## Site 291 4 187.52 0.001 ***
## Distance_rescaled:rescaled_Time.Diff 290 1 20.29 0.009 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Canis.aureus Gazella.dorcas
## Dev Pr(>Dev) Dev
## (Intercept)
## Distance_rescaled 21.032 0.001 5.396
## rescaled_Time.Diff 28.206 0.001 9.241
## Site 13.747 0.087 41.683
## Distance_rescaled:rescaled_Time.Diff 8.429 0.038 4.08
## Hyaena.hyaena
## Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 0.040 2.31 0.150
## rescaled_Time.Diff 0.011 7.259 0.030
## Site 0.001 22.716 0.015
## Distance_rescaled:rescaled_Time.Diff 0.220 1.285 0.619
## Hystrix.indica Lepus.capensis
## Dev Pr(>Dev) Dev
## (Intercept)
## Distance_rescaled 40.08 0.001 19.714
## rescaled_Time.Diff 6.751 0.030 0.978
## Site 30.577 0.003 25.846
## Distance_rescaled:rescaled_Time.Diff 0.384 0.661 5.381
## Vulpes.vulpes
## Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## Distance_rescaled 0.001 8.61 0.010
## rescaled_Time.Diff 0.323 4.185 0.090
## Site 0.007 52.954 0.001
## Distance_rescaled:rescaled_Time.Diff 0.137 0.733 0.661
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.Multivariate test shows that site, distance from settlements, time difference and the interaction between distance from settlement and time difference statistically 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","Site_Ezuz.coef","Site_Merhav_Am.coef","Site_Sde_Boker.coef","Site_Yeruham.coef","Distance_time.coef","Intercept.p","Distance.p","Time_diff.p","Site.p","Distance_time.p","species_abundance")
write.csv(coefp, "coefficients_Negev_Highland.csv")
# plot.Time_diff <- ggplot(coefp[order(Time_diff.coef)], aes(Time_diff.coef,as.numeric(rownames(coefp)))) +
# xlim (-2, 2) +
# geom_point(aes(fill=factor(Time_diff.p<0.1)),color = "black", size=3, shape=21) +
# geom_text(aes(label = paste0("(",species_abundance,") ",SciName)), size=3, nudge_y = -0.3) +
# geom_vline(xintercept = 0, color="black", size=1.5, linetype = "dotted") +
# labs(fill = "") +
# xlab("Effect of time difference on species abundances") +
# ylab("")+
# scale_fill_discrete(breaks = factor(TRUE), labels = "P<0.1", type = c("grey","black"))
#
#
# plot.Distance <- ggplot(coefp[order(Distance.coef)],
# aes(Distance.coef,as.numeric(rownames(coefp)))) +
# xlim (-3, 3) +
# geom_point(aes(fill=factor(Distance.p<0.1)),color = "black", size=3, shape=21) +
# geom_text(aes(label = paste0("(",species_abundance,") ",SciName)), size=3, nudge_y = -0.3) +
# geom_vline(xintercept = 0, color="black", size=1.5, linetype = "dotted") +
# labs(fill = "") +
# xlab("Effect of distance from settlement on species abundances") +
# ylab("")+
# scale_fill_discrete(breaks = factor(TRUE), labels = "P<0.1", type = c("black","grey"))
# plot.Distance.Agri <- ggplot(coefp[order(Distance_agri.coef)],
# aes(Distance_agri.coef,as.numeric(rownames(coefp)))) +
# xlim (-3, 3) +
# geom_point(aes(fill=factor(Distance_agri.p<0.1)),color = "black", size=3, shape=21) +
# geom_text(aes(label = paste0("(",species_abundance,") ",SciName)), size=3, nudge_y = -0.3) +
# geom_vline(xintercept = 0, color="black", size=1.5, linetype = "dotted") +
# labs(fill = "") +
# xlab("Effect of distance from agriculture on species abundances") +
# ylab("")+
# scale_fill_discrete(breaks = factor(TRUE), labels = "P<0.1", type = c("black","grey"))
# plot.Time_diff
# plot.Distance
# plot.Distance.Agri####Canis aureus####
Canis.aureus <- spp_no_rare[,"Canis.aureus"]
glm.Canis.aureus <- glm.nb(formula = Canis.aureus ~ Distance_rescaled * rescaled_Time.Diff + Site, data = env_data)
coef(glm.Canis.aureus)## (Intercept) Distance_rescaled
## 0.202438466 -0.288343650
## rescaled_Time.Diff SiteEzuz
## 0.667695727 0.694395358
## SiteMerhav Am SiteSde Boker
## 0.671853164 -0.002126329
## SiteYeruham Distance_rescaled:rescaled_Time.Diff
## 0.418181180 -0.614087730coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621Canis_inetract <- interact_plot(model = glm.Canis.aureus, data=env_data, pred = rescaled_Time.Diff, modx = Distance_rescaled, modx.values = "plus-minus", partial.residuals = F, 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")
# Canis.interaction.near <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = -0.6386647)
#
# predict(glm.Canis.aureus, newdata = Canis.interaction.near,
# type = "response")
#
# #Increase of 1352.67% over time near settlements
#
# (6.1172275 /0.4211005 ) - 1
#
# #or 14.52 fold
# (6.1172275/0.4211005)
#
#
# Canis.interaction.far <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = 0.6095490)
#
# predict(glm.Canis.aureus, newdata = Canis.interaction.far,
# type = "response")
#
# #Decrease of 109.74% over time far from settlements
#
# (1.5234634 /0.7263442 ) - 1
#
# #or 2.09 fold
# (1.5234634/0.7263442)
canis_aur_interact_plot <- plot_model_interaction(P.anal = env_data,
m = glm.Canis.aureus,
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/negev_highlands/",
legend_position = "bottom")## Loading required package: extrafont## Registering fonts with R
#Testing whether trends and means are different from zero
#Calculate mean and sd of distance
sd_distance <- sd(env_data$Distance_rescaled)
print(sd_distance)## [1] 0.6241069mean_distance <- mean(env_data$Distance_rescaled)
print(mean_distance)## [1] -0.01455781#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.Canis.aureus, 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.639 1.060 0.169 Inf 0.7280 1.392
## 0.610 0.293 0.177 Inf -0.0544 0.641
##
## 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.639 1.060 0.169 Inf 6.259 <.0001
## 0.610 0.293 0.177 Inf 1.653 0.0983
##
## 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.Canis.aureus, ~Distance_rescaled*rescaled_Time.Diff, at = list(Distance_rescaled = c(small_distance,large_distance)))
test_results_distance <- test(pairs(pairwise_distance, by="rescaled_Time.Diff"), by=NULL, adjust="fdr")
print(test_results_distance)## contrast
## (Distance_rescaled-0.638664675609571) - Distance_rescaled0.609549046300846
## rescaled_Time.Diff estimate SE df z.ratio p.value
## 0.417 0.679 0.205 Inf 3.309 0.0009
##
## Results are averaged over the levels of: Site
## Results are given on the log (not the response) scale.temporal trend is significantly different from zero in small distances from settlements (p < 0.0001) and not significantly different from zero in large distances (p = 0.098)
On average over time, Jackal abundance is significantly higher near settlements (p = 0.0009)
####Gazella dorcas####
Gazella.dorcas <- spp_no_rare[,"Gazella.dorcas"]
glm.Gazella.dorcas <- glm.nb(formula = Gazella.dorcas ~ Distance_rescaled + rescaled_Time.Diff + Site, data = env_data)
print(summary(glm.Gazella.dorcas))##
## Call:
## glm.nb(formula = Gazella.dorcas ~ Distance_rescaled + rescaled_Time.Diff +
## Site, data = env_data, init.theta = 0.2088578803, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1017 -0.7419 -0.4797 0.0000 2.2826
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.264e+00 4.337e-01 -2.915 0.003559 **
## Distance_rescaled 1.164e+00 3.724e-01 3.126 0.001772 **
## rescaled_Time.Diff 7.440e-01 2.231e-01 3.335 0.000853 ***
## SiteEzuz 1.179e+00 5.201e-01 2.267 0.023396 *
## SiteMerhav Am 7.277e-02 5.845e-01 0.125 0.900915
## SiteSde Boker 3.752e-01 5.184e-01 0.724 0.469245
## SiteYeruham -3.606e+01 7.395e+06 0.000 0.999996
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.2089) family taken to be 1)
##
## Null deviance: 199.38 on 297 degrees of freedom
## Residual deviance: 130.99 on 291 degrees of freedom
## AIC: 480.61
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.2089
## Std. Err.: 0.0423
##
## 2 x log-likelihood: -464.6090coef(glm.Gazella.dorcas)## (Intercept) Distance_rescaled rescaled_Time.Diff SiteEzuz
## -1.26403563 1.16427008 0.74402049 1.17903270
## SiteMerhav Am SiteSde Boker SiteYeruham
## 0.07276738 0.37518689 -36.05655008coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621# Time graph - need to change y axis scale in plot_model_effect_mammals (uncomment lines 128 - 129 in function)
source("functions/plot_model_effect_mammals.R")
plot_model_effect(P.anal = env_data, m = glm.Gazella.dorcas, 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/negev_highlands/") 
dorcas_time <- effect_plot(model = glm.Gazella.dorcas, 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")
dorcas_time
# #Time references - 01/08/2014 and 01/08/2022
# dorcase.time <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = -0.01455781)
#
# predict(glm.Gazella.dorcas, newdata = dorcase.time,
# type = "response")
#
# #Increase of 554.35% during monitoring year or 6.54 fold
# (0.7550233 /0.1153835 ) - 1
# (0.7550233/0.1153835)
#
# sum(P.anal$`Gazella dorcas`[P.anal$year=="2014"])
# sum(P.anal$`Gazella dorcas`[P.anal$year=="2022"])
#
# #Increase of 1500% during monitoring years or 16 fold
#
# (96/6 - 1) *100
# 96/6
plot_model_effect(P.anal = env_data, m = glm.Gazella.dorcas, 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/negev_highlands/")
dorcas_distance <- effect_plot(model = glm.Gazella.dorcas, data=env_data, pred = Distance_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.820445737,0.018246686,0.950127155,1.882007625), labels=c("100", "1000", "2000", "3000"), name = "Distance (m)")
# dorcase.distance <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = 0.4169025,
# Distance_rescaled = c(-0.820445737,-0.447693549,
# 0.018246686,0.48418692,0.950127155))
#
# predict(glm.Gazella.dorcas, newdata = dorcase.distance,
# type = "response")
#
# #Decrease from 100 to 2000 m of 685.71% or 7.85 fold
# (1.1645712/0.1482179) - 1
# (1.1645712/0.1482179)
#
# # Decrease of 172.02% every 500 m
#
# 0.3935262 / 0.2287586
# 0.6769707 / 0.3935262
####Hyaena hyaena####
Hyaena.hyaena <- spp_no_rare[,"Hyaena.hyaena"]
glm.Hyaena.hyaena <- glm.nb(formula = Hyaena.hyaena ~ Distance_rescaled + rescaled_Time.Diff + Site, data = env_data)
print(summary(glm.Hyaena.hyaena))##
## Call:
## glm.nb(formula = Hyaena.hyaena ~ Distance_rescaled + rescaled_Time.Diff +
## Site, data = env_data, init.theta = 0.6204385678, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1327 -0.8026 -0.5875 -0.2873 3.4130
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.59934 0.26764 -2.239 0.025132 *
## Distance_rescaled 0.09614 0.19570 0.491 0.623243
## rescaled_Time.Diff 0.43844 0.15626 2.806 0.005018 **
## SiteEzuz -1.99830 0.54309 -3.680 0.000234 ***
## SiteMerhav Am -0.73115 0.39828 -1.836 0.066393 .
## SiteSde Boker -0.90650 0.38835 -2.334 0.019584 *
## SiteYeruham -0.04114 0.33606 -0.122 0.902572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.6204) family taken to be 1)
##
## Null deviance: 233.72 on 297 degrees of freedom
## Residual deviance: 197.48 on 291 degrees of freedom
## AIC: 491.64
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.620
## Std. Err.: 0.179
##
## 2 x log-likelihood: -475.642coef(glm.Hyaena.hyaena)## (Intercept) Distance_rescaled rescaled_Time.Diff SiteEzuz
## -0.59933736 0.09614149 0.43844023 -1.99829860
## SiteMerhav Am SiteSde Boker SiteYeruham
## -0.73115347 -0.90649712 -0.04113771coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621# Time graph - need to change y axis scale in plot_model_effect_mammals (uncomment lines 128 - 129 in function)
plot_model_effect(P.anal = env_data, m = glm.Hyaena.hyaena, 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/negev_highlands/", y_cutoff = 2, y_breaks = seq(0,2, by = 0.5))
hyaena_time <- effect_plot(model = glm.Hyaena.hyaena, 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")
hyaena_time
#Time references - 01/08/2014 and 01/08/2022
# hyaena.time <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = -0.01455781)
#
# predict(glm.Hyaena.hyaena, newdata = hyaena.time,
# type = "response")
#
# #Increase of 202.51% during monitoring year or 3.02 fold
# (0.9886008 /0.3267927 ) - 1
# (0.9886008/0.3267927)
#
# sum(P.anal$`Hyaena hyaena`[P.anal$year=="2014"])
# sum(P.anal$`Hyaena hyaena`[P.anal$year=="2022"])
#Increase of 1966.66% during monitoring years or 20.66 fold
# (62/3 - 1) *100
# 62/3
#### Hystrix indica####
Hystrix.indica <- spp_no_rare[,"Hystrix.indica"]
glm.Hystrix.indica <- glm.nb(formula = Hystrix.indica ~ Distance_rescaled + rescaled_Time.Diff + Site, data = env_data)
print(summary(glm.Hystrix.indica))##
## Call:
## glm.nb(formula = Hystrix.indica ~ Distance_rescaled + rescaled_Time.Diff +
## Site, data = env_data, init.theta = 0.472641149, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.72740 -1.04949 -0.73835 0.09607 2.05553
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3108 0.2742 -1.133 0.257020
## Distance_rescaled -0.8373 0.1694 -4.942 7.71e-07 ***
## rescaled_Time.Diff 0.4672 0.1262 3.702 0.000214 ***
## SiteEzuz -0.1266 0.3710 -0.341 0.732819
## SiteMerhav Am 0.3407 0.3594 0.948 0.343077
## SiteSde Boker 0.0197 0.3588 0.055 0.956218
## SiteYeruham 1.3618 0.3338 4.080 4.51e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.4726) family taken to be 1)
##
## Null deviance: 363.71 on 297 degrees of freedom
## Residual deviance: 265.39 on 291 degrees of freedom
## AIC: 963.22
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.4726
## Std. Err.: 0.0659
##
## 2 x log-likelihood: -947.2160coef(glm.Hystrix.indica)## (Intercept) Distance_rescaled rescaled_Time.Diff SiteEzuz
## -0.31077024 -0.83728054 0.46723666 -0.12664469
## SiteMerhav Am SiteSde Boker SiteYeruham
## 0.34075174 0.01969547 1.36184299coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621# Time graph - need to change y axis scale in plot_model_effect_mammals (uncomment lines 128 - 129 in function)
plot_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/negev_highlands/", y_cutoff = 2.7, y_breaks = seq(0,2.5, by = 0.5))
hystrix_time <- effect_plot(model = glm.Hystrix.indica, 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")
hystrix_time
#Time references - 01/08/2014 and 01/08/2022
# hystrix.time <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = -0.01455781)
#
# predict(glm.Hystrix.indica, newdata = hystrix.time,
# type = "response")
#
# #Increase of 225.32% during monitoring year or 3.25 fold
# (1.3901255 /0.4272973 ) - 1
# (1.3901255/0.4272973)
#
# sum(P.anal$`Hystrix indica`[P.anal$year=="2014"])
# sum(P.anal$`Hystrix indica`[P.anal$year=="2022"])
#Increase of 216% during monitoring years or 3.16 fold
# (158/50 - 1) *100
# 158/50
plot_model_effect(P.anal = env_data, m = glm.Hystrix.indica, 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/negev_highlands/")
hystrix_distance<- effect_plot(model = glm.Hystrix.indica, data=env_data, pred = Distance_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.820445737,0.018246686,0.950127155,1.882007625), labels=c("100", "1000", "2000", "3000"), name = "Distance (m)")
hystrix_distance
# hystrix.distance <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = 0.4169025,
# Distance_rescaled = c(-0.820445737,-0.447693549,
# 0.018246686,0.48418692,0.950127155))
#
# predict(glm.Hystrix.indica, newdata = hystrix.distance,
# type = "response")
#
# #Decrease from 100 to 2000 m of 340.37% or 4.40 fold
# (1.7699745/0.4019205) - 1
# (1.7699745/0.4019205)
#
# # Decrease of 67.69% every 500 m
#
# 0.5937016 / 0.8769935
# 0.8769935 / 1.2954615
####Lepus capensis####
Lepus.capensis <- spp_no_rare[,"Lepus.capensis"]
glm.Lepus.capensis <- glm.nb(formula = Lepus.capensis ~ Distance_rescaled + rescaled_Time.Diff + Site, data = env_data)
coef(glm.Lepus.capensis)## (Intercept) Distance_rescaled rescaled_Time.Diff SiteEzuz
## -0.9426391 1.2967504 0.1351331 1.1455811
## SiteMerhav Am SiteSde Boker SiteYeruham
## -0.4684696 -0.7146170 -0.2738756coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621plot_model_effect(P.anal = env_data, m = glm.Lepus.capensis, 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/negev_highlands/")
lepus_distance <- effect_plot(model = glm.Lepus.capensis, data=env_data, pred = Distance_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.820445737,0.018246686,0.950127155,1.882007625), labels=c("100", "1000", "2000", "3000"), name = "Distance (m)")
lepus_distance
# lepus.distance <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = 0.4169025,
# Distance_rescaled = c(-0.820445737,-0.447693549,
# 0.018246686,0.48418692,0.950127155))
#
# predict(glm.Lepus.capensis, newdata = lepus.distance,
# type = "response")
#
# #Increase from 100 to 2000 m of 893.42% or 9.93 fold
# (1.4130860/0.1422435) - 1
# (1.4130860/0.1422435)
#
# # Increase of 182.98% every 500 m
#
# 0.4220462 / 0.2306512
# 0.772261 / 0.4220462
####Vulpus vulpus####
Vulpes.vulpes <- spp_no_rare[,"Vulpes.vulpes"]
glm.Vulpes.vulpes <- glm.nb(formula = Vulpes.vulpes ~ Distance_rescaled + rescaled_Time.Diff + Site, data = env_data)
print(summary(glm.Vulpes.vulpes))##
## Call:
## glm.nb(formula = Vulpes.vulpes ~ Distance_rescaled + rescaled_Time.Diff +
## Site, data = env_data, init.theta = 0.7237539999, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6226 -0.9097 -0.6285 0.1305 3.6909
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8257 0.2705 -3.052 0.002271 **
## Distance_rescaled -0.2356 0.1441 -1.635 0.102124
## rescaled_Time.Diff 0.3220 0.1107 2.910 0.003615 **
## SiteEzuz -0.3374 0.3846 -0.877 0.380356
## SiteMerhav Am 1.5615 0.3275 4.768 1.86e-06 ***
## SiteSde Boker 1.5141 0.3177 4.766 1.88e-06 ***
## SiteYeruham 1.1820 0.3185 3.711 0.000206 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.7238) family taken to be 1)
##
## Null deviance: 352.32 on 297 degrees of freedom
## Residual deviance: 276.91 on 291 degrees of freedom
## AIC: 942.57
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.724
## Std. Err.: 0.111
##
## 2 x log-likelihood: -926.574coef(glm.Vulpes.vulpes)## (Intercept) Distance_rescaled rescaled_Time.Diff SiteEzuz
## -0.8257402 -0.2356295 0.3220456 -0.3373799
## SiteMerhav Am SiteSde Boker SiteYeruham
## 1.5614529 1.5141002 1.1820256coef(mva_m2)## Canis.aureus Gazella.dorcas Hyaena.hyaena
## (Intercept) 0.202471624 -1.0331143 -0.673103688
## Distance_rescaled -0.288347160 1.5811782 -0.063669196
## rescaled_Time.Diff 0.667674257 0.8084377 0.447480193
## SiteEzuz 0.694390321 0.9537089 -1.931047366
## SiteMerhav Am 0.671816264 -0.4328389 -0.660501490
## SiteSde Boker -0.002144141 -0.1411651 -0.830281862
## SiteYeruham 0.418133246 -13.9665941 -0.004746017
## Distance_rescaled:rescaled_Time.Diff -0.614116326 -0.9795058 0.285629373
## Hystrix.indica Lepus.capensis
## (Intercept) -0.26649537 -1.08313007
## Distance_rescaled -0.77978978 1.07360882
## rescaled_Time.Diff 0.44329402 0.04304816
## SiteEzuz -0.16920989 1.37732752
## SiteMerhav Am 0.30876418 -0.36155749
## SiteSde Boker -0.03221625 -0.57325986
## SiteYeruham 1.35263904 -0.22208773
## Distance_rescaled:rescaled_Time.Diff -0.13312087 0.68802607
## Vulpes.vulpes
## (Intercept) -0.8808563
## Distance_rescaled -0.3101542
## rescaled_Time.Diff 0.3323375
## SiteEzuz -0.2880808
## SiteMerhav Am 1.6010889
## SiteSde Boker 1.5908155
## SiteYeruham 1.2124690
## Distance_rescaled:rescaled_Time.Diff 0.1554621plot_model_effect(P.anal = env_data, m = glm.Vulpes.vulpes, 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/negev_highlands/")
vulpes_time <- effect_plot(model = glm.Vulpes.vulpes, 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")
vulpes_time
#Time references - 01/08/2014 and 01/08/2022
# vulpes.time <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = c(-1.180759317,1.344019721),
# Distance_rescaled = -0.01455781)
#
# predict(glm.Vulpes.vulpes, newdata = vulpes.time,
# type = "response")
#
# #Increase of 125.45% during monitoring year or 2.25 fold
# (0.6774129 /0.3004216 ) - 1
# (0.6774129/0.3004216)
# sum(P.anal$`Vulpes vulpes`[P.anal$year=="2016"])
# sum(P.anal$`Vulpes vulpes`[P.anal$year=="2022"])
# range(P.anal$Distance[P.anal$Settlements=="Near"])
# range(P.anal$Distance[P.anal$Settlements=="Far"])
#
# #Increase of 39.6% during monitoring years (here only from 2016 because of no observations in 2014) or 1.39 fold
#
# (141/101 - 1) *100
# 141/101
vulpes_distance <- plot_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/negev_highlands/")
vulpes_distance <- effect_plot(model = glm.Vulpes.vulpes, data=env_data, pred = Distance_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.820445737,0.018246686,0.950127155,1.882007625), labels=c("100", "1000", "2000", "3000"), name = "Distance (m)")
vulpes_distance
# vulpes.distance <- data.frame(Site = "Bislach",
# rescaled_Time.Diff = 0.4169025,
# Distance_rescaled = c(-0.820445737,-0.447693549,
# 0.018246686,0.48418692,0.950127155))
#
# predict(glm.Vulpes.vulpes, newdata = vulpes.distance,
# type = "response")
#
# #Decrase from 100 to 2000 m of 51.77% or 1.51 fold
# (0.6076493 /0.4003738 ) - 1
# (0.6076493 /0.4003738 )
#
# # Decrease of 89.6% every 500 m
#
# 0.4468344 / 0.4986863
# 0.4986863 / 0.5565553 n = 298 (total number of cameras)
total abundance (all species) = 2540, abundance without rare species is the same as total abundance
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)
## 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)
## 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
##
## ──────────────────────────────────────────────────────────────────────────────