Birds_Maquis_analysis_by_plot
Ron Chen
2023-11-02
out <- prepare_bird_data()## Loading required package: readxl## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following objects are masked from 'package:data.table':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union## Loading required package: solartime## - Winter data excluded.
## - The counts of all passing / wintering species (not resident nor breeding in the sampled unit) were set to zero. Abreed_and_non_breed contains the non-breeders as well.var_names <- names(out)
for (i in 1:length(var_names)) {
assign(var_names[i], out[[i]])
}
rm(out)
P <- P_byplot[grepl("Maquis",unit),][order(year,site)][,site:=factor(site)][,subunit:=factor(subunit)]
P_no_rare <- P_byplot_no_rare[grepl("Maquis",unit),][order(year,site)][,subunit:=factor(subunit)]Mediterranean Maquis
Unit is divided into 3 subunits: Judea, Carmel and Galilee. Factors are proximity to settlements and time. Sampling started in spring 2012 (pilot year), but during T0 sampling was performed only in winter of 2014, and the next spring sampling was done in T1 (2015). Therefore 2012 will not be considered as pilot here but as T0. Total 5 campaigns, 3subunits per campaign, 5 sites per subunit, with 6 plots per site (total of 450 plot-campaign combinations).
Raw data Total abundance: 21855 Number of observations: 5795 Total richness: 108
Filtered data Total abundance: 13796 Number of observations: 4265 Total richness: 49
Model gma, abundance and richness
Richness done with rare species. Abundance and mean abundance done without rare species. Full models include cosine and sine of the time difference from June 21st (in radians).
richness
Explore data and plot mean-variance plot. There is a strong relationship, indicating that employing GLMs is the proper way to analyze, rather than OLS (assumption of homogeneity is violated).
# 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
print("RICHNESS WITH RARE SPECIES")## [1] "RICHNESS WITH RARE SPECIES"plot_alpha_diversity(P.anal, x_val = "settlements", y_val = "richness", ylab_val = "richness", xlab_val = "settlements")## ℹ SHA-1 hash of file is "eef55df1963ba201f57cf44513f566298756fd4f"
plot_alpha_diversity(P.anal, x_val = "subunit", y_val = "richness", ylab_val = "richness", xlab_val = "subunit")
dotchart(P.anal$richness, groups = P.anal$settlements, pch = as.numeric(P.anal$settlements), ylab = "settlements", xlab = "richness")
dotchart(P.anal$richness, groups = P.anal$subunit, pch = as.numeric(P.anal$subunit), ylab = "subunit", xlab = "richness")
IVs <- c("richness", "year_ct","site","settlements","subunit", "td_sc", "cos_td_rad", "sin_td_rad", "h_from_sunrise", "cos_hsun", "sin_hsun", "monitors_name",
"wind", "precipitation", "temperature", "clouds")
kable(summary(P.anal[,..IVs]))%>% kable_styling(font_size = 11)| richness | year_ct | site | settlements | subunit | td_sc | cos_td_rad | sin_td_rad | h_from_sunrise | cos_hsun | sin_hsun | monitors_name | wind | precipitation | temperature | clouds | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 0.000 | Min. :0.0 | Nir Etzion : 31 | Far :224 | Judean Highlands:149 | Min. :-1.9882 | Min. :-0.1671 | Min. :-1.0000 | Length:444 | Min. :-0.3289 | Min. :-0.1253 | Sassi Haham : 53 | 0 : 33 | 0 :119 | 0 : 4 | 0 : 32 | |
| 1st Qu.: 6.000 | 1st Qu.:3.0 | Ein Yaakov : 30 | Near:220 | Carmel :150 | 1st Qu.:-0.2713 | 1st Qu.: 0.5702 | 1st Qu.:-0.8215 | Class :difftime | 1st Qu.: 0.6855 | 1st Qu.: 0.2563 | Eyal Shochat: 47 | 1 : 42 | 3 : 2 | 1 : 38 | 1 : 0 | |
| Median : 8.000 | Median :5.0 | Givat Yearim : 30 | NA | Galilee :145 | Median : 0.4918 | Median : 0.8140 | Median :-0.5808 | Mode :numeric | Median : 0.8738 | Median : 0.4863 | Eran Banker : 30 | 2 : 12 | NA’s:323 | 2 : 37 | 2 : 5 | |
| Mean : 8.054 | Mean :4.8 | Givat Yeshayahu: 30 | NA | NA | Mean : 0.3205 | Mean : 0.7098 | Mean :-0.5864 | NA | Mean : 0.7997 | Mean : 0.4836 | Asaf Mayrose: 18 | 3 : 2 | NA | 3 : 18 | 3 : 13 | |
| 3rd Qu.:10.000 | 3rd Qu.:7.0 | Goren : 30 | NA | NA | 3rd Qu.: 1.1023 | 3rd Qu.: 0.9413 | 3rd Qu.:-0.3375 | NA | 3rd Qu.: 0.9666 | 3rd Qu.: 0.7281 | Ohad Sharir : 18 | NA’s:355 | NA | NA’s:347 | NA’s:394 | |
| Max. :18.000 | Max. :9.0 | Kerem Maharal : 30 | NA | NA | Max. : 1.7127 | Max. : 0.9976 | Max. :-0.0688 | NA | Max. : 1.0000 | Max. : 0.9998 | (Other) : 58 | NA | NA | NA | NA | |
| NA | NA | (Other) :263 | NA | NA | NA | NA | NA | NA | NA’s :89 | NA’s :89 | NA’s :220 | NA | NA | NA | NA |
no observation for all 4 weather variables. many NAs for sampling time of day variables.exclude from model. An extreme observation of richness>20 in Judean highlands near settlements.
pairs(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=IVs[1:11]])
kable(cor(P.anal[,lapply(X = .SD,FUN = as.numeric),.SDcols=IVs[1:11]], use = "pairwise.complete.obs")) %>% kable_styling(font_size = 11)| richness | year_ct | site | settlements | subunit | td_sc | cos_td_rad | sin_td_rad | h_from_sunrise | cos_hsun | sin_hsun | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| richness | 1.0000000 | -0.0647635 | 0.0926545 | 0.4762118 | -0.1936858 | 0.0617968 | 0.0751315 | 0.0460403 | -0.0582024 | 0.0933868 | -0.0342460 |
| year_ct | -0.0647635 | 1.0000000 | -0.0460962 | -0.0041782 | 0.0046005 | -0.6240340 | -0.5599811 | -0.6642968 | 0.0244310 | -0.0070326 | 0.0341395 |
| site | 0.0926545 | -0.0460962 | 1.0000000 | 0.0038640 | -0.1472070 | 0.0103683 | 0.0017250 | 0.0296647 | -0.0105111 | -0.0092740 | -0.0240278 |
| settlements | 0.4762118 | -0.0041782 | 0.0038640 | 1.0000000 | -0.0167082 | 0.0029695 | 0.0016285 | 0.0037459 | 0.0341209 | -0.0071595 | 0.0507160 |
| subunit | -0.1936858 | 0.0046005 | -0.1472070 | -0.0167082 | 1.0000000 | -0.0056370 | -0.0129685 | -0.0136922 | 0.1193296 | -0.1166460 | 0.1174210 |
| td_sc | 0.0617968 | -0.6240340 | 0.0103683 | 0.0029695 | -0.0056370 | 1.0000000 | 0.9788766 | 0.9693667 | -0.1975082 | 0.1811635 | -0.2043768 |
| cos_td_rad | 0.0751315 | -0.5599811 | 0.0017250 | 0.0016285 | -0.0129685 | 0.9788766 | 1.0000000 | 0.9008724 | -0.1878186 | 0.1700185 | -0.1970275 |
| sin_td_rad | 0.0460403 | -0.6642968 | 0.0296647 | 0.0037459 | -0.0136922 | 0.9693667 | 0.9008724 | 1.0000000 | -0.2005880 | 0.1868900 | -0.2040813 |
| h_from_sunrise | -0.0582024 | 0.0244310 | -0.0105111 | 0.0341209 | 0.1193296 | -0.1975082 | -0.1878186 | -0.2005880 | 1.0000000 | -0.9598144 | 0.9818950 |
| cos_hsun | 0.0933868 | -0.0070326 | -0.0092740 | -0.0071595 | -0.1166460 | 0.1811635 | 0.1700185 | 0.1868900 | -0.9598144 | 1.0000000 | -0.8930078 |
| sin_hsun | -0.0342460 | 0.0341395 | -0.0240278 | 0.0507160 | 0.1174210 | -0.2043768 | -0.1970275 | -0.2040813 | 0.9818950 | -0.8930078 | 1.0000000 |
Fit Poisson glm, check for existence of overdispersion
mdl_r.poiss.int <- glm(data = P.anal, formula = richness ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad , 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.917251305665046"Overdispersion parameter is < 1. Choose Poisson.
m0 <- mdl_r.poiss.int
me0 <- glmer(data = P.anal, formula = richness ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)vif(me0)## GVIF Df GVIF^(1/(2*Df))
## settlements 3.275497 1 1.809834
## year_ct 3.277089 1 1.810273
## subunit 1.001389 2 1.000347
## cos_td_rad 5.439321 1 2.332235
## sin_td_rad 6.692207 1 2.586930
## settlements:year_ct 4.690440 1 2.165742go with mixed model, attempt model selection yet.
summary(me0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## richness ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 2115.0 2151.8 -1048.5 2097.0 435
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.04280 -0.60906 -0.07994 0.57865 2.87628
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009067 0.09522
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.575932 0.201233 7.831 4.83e-15 ***
## settlementsNear 0.393867 0.061554 6.399 1.57e-10 ***
## year_ct -0.004072 0.009709 -0.419 0.67493
## subunitCarmel 0.082972 0.071979 1.153 0.24902
## subunitGalilee -0.198418 0.072660 -2.731 0.00632 **
## cos_td_rad 0.299159 0.148180 2.019 0.04350 *
## sin_td_rad -0.228252 0.169323 -1.348 0.17765
## settlementsNear:year_ct -0.002007 0.010854 -0.185 0.85329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct sbntCr sbntGl cs_td_ sn_td_
## settlmntsNr -0.190
## year_ct -0.027 0.552
## subunitCrml -0.195 -0.002 0.005
## subunitGall -0.199 0.003 -0.001 0.501
## cos_td_rad -0.931 0.008 -0.080 0.019 0.019
## sin_td_rad 0.866 -0.010 0.328 -0.006 -0.018 -0.853
## sttlmntsN:_ 0.158 -0.833 -0.660 0.003 0.003 -0.006 0.009perform stepwise model selection of poisson mixed model.
drop1(me0)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0035507 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## richness ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 2115.0
## subunit 2 2122.2
## cos_td_rad 1 2117.1
## sin_td_rad 1 2114.8
## settlements:year_ct 1 2113.0remove settlements X year.
me1 <- glmer(data = P.anal, formula = richness ~ settlements + year_ct + subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)
drop1(me1)## Single term deletions
##
## Model:
## richness ~ settlements + year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 2113.0
## settlements 1 2240.8
## year_ct 1 2111.5
## subunit 2 2120.3
## cos_td_rad 1 2115.1
## sin_td_rad 1 2112.8drop year
me1 <- glmer(data = P.anal, formula = richness ~ settlements + subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)
drop1(me1)## Single term deletions
##
## Model:
## richness ~ settlements + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 2111.5
## settlements 1 2239.4
## subunit 2 2118.7
## cos_td_rad 1 2113.3
## sin_td_rad 1 2110.8drop sine.
me1 <- glmer(data = P.anal, formula = richness ~ settlements + subunit + cos_td_rad + (1|site), family = poisson)
drop1(me1)## Single term deletions
##
## Model:
## richness ~ settlements + subunit + cos_td_rad + (1 | site)
## npar AIC
## <none> 2110.8
## settlements 1 2238.5
## subunit 2 2118.0
## cos_td_rad 1 2113.5Settlements and cosin time of year remain. Final model:
summary(me1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements + subunit + cos_td_rad + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 2110.8 2135.4 -1049.4 2098.8 438
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09026 -0.59675 -0.07933 0.56203 2.96581
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009619 0.09808
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.80687 0.07258 24.896 <2e-16 ***
## settlementsNear 0.38421 0.03401 11.296 <2e-16 ***
## subunitCarmel 0.08257 0.07349 1.124 0.2612
## subunitGalilee -0.20004 0.07405 -2.702 0.0069 **
## cos_td_rad 0.13535 0.06305 2.147 0.0318 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl
## settlmntsNr -0.279
## subunitCrml -0.526 0.002
## subunitGall -0.508 0.009 0.502
## cos_td_rad -0.635 -0.002 0.021 -0.003ranef(me1)## $site
## (Intercept)
## Abirim -7.321848e-02
## Aderet 6.137108e-02
## Beit Oren -2.233669e-02
## Ein Yaakov -2.970003e-02
## Givat Yearim -7.623271e-02
## Givat Yeshayahu 9.364091e-02
## Goren -7.232404e-03
## Iftach 1.773011e-01
## Kerem Maharal 9.669359e-02
## Kfar Shamai -5.892993e-02
## Margaliot -1.725282e-03
## Nehusha -4.299542e-02
## Nir Etzion -1.521589e-01
## Ofer 8.260917e-02
## Ramat Raziel -3.076398e-02
## Yagur 2.684868e-05
##
## with conditional variances for "site"dotplot(ranef(me1, condVar=TRUE))## $site
plot(me1)
qqmath(me1)
# scale location plot
plot(me1, sqrt(abs(resid(.))) ~ fitted(.),
type = c("p", "smooth"),
main="scale-location",
par.settings = list(plot.line =
list(alpha=1, col = "red",
lty = 1, lwd = 2)))
attempt to add year, see if significant:
summary(glmer(data = P.anal, formula = richness ~ year_ct + settlements + subunit + cos_td_rad + (1|site), family = poisson))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ year_ct + settlements + subunit + cos_td_rad + (1 |
## site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 2112.8 2141.5 -1049.4 2098.8 437
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0946 -0.6000 -0.0782 0.5691 2.9651
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009581 0.09788
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.8147139 0.0962891 18.847 <2e-16 ***
## year_ct -0.0008103 0.0065499 -0.124 0.9015
## settlementsNear 0.3841894 0.0340127 11.295 <2e-16 ***
## subunitCarmel 0.0824132 0.0733925 1.123 0.2615
## subunitGalilee -0.2001589 0.0739537 -2.707 0.0068 **
## cos_td_rad 0.1298876 0.0769263 1.688 0.0913 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) yer_ct sttlmN sbntCr sbntGl
## year_ct -0.658
## settlmntsNr -0.214 0.005
## subunitCrml -0.407 0.016 0.002
## subunitGall -0.391 0.013 0.009 0.502
## cos_td_rad -0.769 0.574 0.001 0.027 0.005year not significant, rightfully dropped. center time of year variable, highly correlated with intercept.
P.anal[,`:=`(cos_td_rad_c=cos_td_rad-mean(cos_td_rad),
sin_td_rad_c=sin_td_rad-mean(sin_td_rad))]
me1 <- glmer(data = P.anal, formula = richness ~ settlements + subunit + cos_td_rad_c + (1|site), family = poisson)
summary(me1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements + subunit + cos_td_rad_c + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 2110.8 2135.4 -1049.4 2098.8 438
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09026 -0.59674 -0.07933 0.56201 2.96580
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009619 0.09808
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.90293 0.05611 33.916 <2e-16 ***
## settlementsNear 0.38421 0.03401 11.296 <2e-16 ***
## subunitCarmel 0.08256 0.07349 1.124 0.2612
## subunitGalilee -0.20004 0.07404 -2.702 0.0069 **
## cos_td_rad_c 0.13534 0.06305 2.147 0.0318 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl
## settlmntsNr -0.363
## subunitCrml -0.663 0.002
## subunitGall -0.660 0.009 0.502
## cos_td_rd_c -0.023 -0.002 0.021 -0.003P.anal <- P.anal
summ(me1, exp=T, digits = 3)## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 444 |
| Dependent variable | richness |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 2110.815 |
| BIC | 2135.390 |
| Pseudo-R² (fixed effects) | 0.292 |
| Pseudo-R² (total) | 0.345 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 6.706 | 0.056 | 33.916 | 0.000 |
| settlementsNear | 1.468 | 0.034 | 11.296 | 0.000 |
| subunitCarmel | 1.086 | 0.073 | 1.124 | 0.261 |
| subunitGalilee | 0.819 | 0.074 | -2.702 | 0.007 |
| cos_td_rad_c | 1.145 | 0.063 | 2.147 | 0.032 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.098 |
| Group | # groups | ICC |
|---|---|---|
| site | 16 | 0.010 |
# effplot1 <- effect_plot(model = me1, data=P.anal, pred = settlements, interval = T, plot.points = F, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = T, colors = "Qual1", main.title = "effect of proximity to settlements (points are partial residuals)", line.colors = "black", point.alpha = 0.25)
# effplot1$layers[[2]]$geom_params$width <- 0.4
effplot1 <- plot_model_effect(P.anal = P.anal, m = me1, eff2plot = "settlements", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")## Loading required package: Cairo## Loading required package: extrafont## Registering fonts with R## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
effplot1
#Near is 46.84% higher
iv1_abs <- diff(effplot1$data$richness)
print(iv1_abs)## 2
## 3.141228iv1_rel <- iv1_abs / as.data.table(effplot1$data)[settlements=="Far",richness] * 100
print(iv1_rel)## 2
## 46.84533effplot1$data$richness## 1 2
## 6.705529 9.8467579.846757 / 6.705529 - 1## [1] 0.4684534effplot2 <- effect_plot(model = me1, data=P.anal, pred = subunit, interval = T, plot.points = F, partial.residuals = T, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, colors = "Qual1", main.title = "effect of subunit (points are partial residuals)", line.colors = "black", point.alpha = 0.25)## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.effplot2$layers[[2]]$geom_params$width <- 0.4
effplot2
#emmeans
EMM <- emmeans(object = me1, ~subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 2.10 0.0524 Inf 1.99 2.20
## Carmel 2.18 0.0518 Inf 2.08 2.28
## Galilee 1.89 0.0526 Inf 1.79 2.00
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95# Carmel vs Galilee
(exp(2.18)/exp(1.89) - 1) * 100## [1] 33.64275# Carmel vs Judea
(exp(2.18)/exp(2.1) - 1) * 100## [1] 8.328707#Judea vs Galilee
(exp(2.1)/exp(1.89) - 1) * 100## [1] 23.36781test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -0.0826 0.0735 Inf -1.124 0.2612
## Judean Highlands - Galilee 0.2000 0.0740 Inf 2.702 0.0103
## Carmel - Galilee 0.2826 0.0736 Inf 3.839 0.0004
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testsstatistically significant lower richness in Galilee and far from settlements. No significant change in richness over time.
geometric mean of abundance
Explore data. Exclude time of day because of high number of NAs.
# exclude rare species in analysis
P.anal <- copy(P_no_rare)
P.anal <- P.anal[pilot==FALSE]
print("GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES")## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"plot_alpha_diversity(P.anal, x_val = "settlements", y_val = "gma", ylab_val = "gma", xlab_val = "settlements")## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Removed 1 rows containing non-finite values (`stat_boxplot()`).
plot_alpha_diversity(P.anal, x_val = "subunit", y_val = "gma", ylab_val = "gma", xlab_val = "subunit")## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Removed 1 rows containing non-finite values (`stat_boxplot()`).
## Removed 1 rows containing non-finite values (`stat_boxplot()`).
dotchart(P.anal$gma, groups = P.anal$settlements, pch = as.numeric(P.anal$settlements), ylab = "settlements", xlab = "gma")
dotchart(P.anal$gma, groups = P.anal$subunit, pch = as.numeric(P.anal$subunit), ylab = "subunit", xlab = "gma")
IVs <- c("gma", "year_ct","site","settlements","subunit", "td_sc", "cos_td_rad", "sin_td_rad")
kable(summary(P.anal[,..IVs]))%>% kable_styling(font_size = 11)| gma | year_ct | site | settlements | subunit | td_sc | cos_td_rad | sin_td_rad | |
|---|---|---|---|---|---|---|---|---|
| Min. : 1.000 | Min. :0.0 | Nir Etzion : 31 | Far :224 | Judean Highlands:149 | Min. :-1.9882 | Min. :-0.1671 | Min. :-1.0000 | |
| 1st Qu.: 2.076 | 1st Qu.:3.0 | Ein Yaakov : 30 | Near:220 | Carmel :150 | 1st Qu.:-0.2713 | 1st Qu.: 0.5702 | 1st Qu.:-0.8215 | |
| Median : 2.560 | Median :5.0 | Givat Yearim : 30 | NA | Galilee :145 | Median : 0.4918 | Median : 0.8140 | Median :-0.5808 | |
| Mean : 2.818 | Mean :4.8 | Givat Yeshayahu: 30 | NA | NA | Mean : 0.3205 | Mean : 0.7098 | Mean :-0.5864 | |
| 3rd Qu.: 3.252 | 3rd Qu.:7.0 | Goren : 30 | NA | NA | 3rd Qu.: 1.1023 | 3rd Qu.: 0.9413 | 3rd Qu.:-0.3375 | |
| Max. :10.301 | Max. :9.0 | Kerem Maharal : 30 | NA | NA | Max. : 1.7127 | Max. : 0.9976 | Max. :-0.0688 | |
| NA’s :1 | NA | (Other) :263 | NA | NA | NA | NA | NA |
Fit glm, compare gamma, gaussian (poisson inappropriate because response is not discrete)
# mdl_g.gamma.int <- glm(data = P.anal, formula = gma ~ settlements * year_ct + subunit*year_ct + cos_td_rad + sin_td_rad, family = Gamma) # run this if there are no zeros
# plot(mdl_g.gamma.int,which = 1:3,sub.caption = "gamma")
mdl_g.gauss.int <- glm(data = P.anal, formula = gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad, family = gaussian)
plot(mdl_g.gauss.int,which = 1:3,sub.caption = "gaussian")
Remove rows 211, 376, 408. Fit fixed and mixed models.
P.anal <- P.anal[!c(211, 376, 408)]
m0 <- lm(data = P.anal, formula = gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad)
me0 <- lmer(data = P.anal, formula = gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad + (1|site))Mixed model converged
summary(me0)## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 1292.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9628 -0.6076 -0.1482 0.4344 5.8495
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.07491 0.2737
## Residual 1.02960 1.0147
## Number of obs: 440, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.99047 0.57113 1.734
## settlementsNear 0.33563 0.17729 1.893
## year_ct -0.01007 0.02642 -0.381
## subunitCarmel 0.13762 0.20953 0.657
## subunitGalilee -0.07988 0.20628 -0.387
## cos_td_rad 1.17762 0.42045 2.801
## sin_td_rad -1.48469 0.49359 -3.008
## settlementsNear:year_ct -0.02829 0.03101 -0.912
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct sbntCr sbntGl cs_td_ sn_td_
## settlmntsNr -0.157
## year_ct 0.022 0.487
## subunitCrml -0.191 -0.001 0.008
## subunitGall -0.202 0.003 0.005 0.508
## cos_td_rad -0.934 0.003 -0.106 0.012 0.015
## sin_td_rad 0.871 -0.004 0.374 0.001 -0.014 -0.856
## sttlmntsN:_ 0.132 -0.838 -0.582 0.003 0.001 -0.004 0.003perform stepwise model selection of gaussian model.
drop1(me0)## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 1290.5
## subunit 2 1287.8
## cos_td_rad 1 1296.3
## sin_td_rad 1 1297.7
## settlements:year_ct 1 1289.3drop subunit
me1 <- lmer(data = P.anal, formula = gma ~ settlements * year_ct + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 1287.8
## cos_td_rad 1 1293.6
## sin_td_rad 1 1295.0
## settlements:year_ct 1 1286.6drop year X settlements
me1 <- lmer(data = P.anal, formula = gma ~ settlements + year_ct + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements + year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 1286.6
## settlements 1 1289.0
## year_ct 1 1286.0
## cos_td_rad 1 1292.5
## sin_td_rad 1 1293.9drop year
me1 <- lmer(data = P.anal, formula = gma ~ settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements + cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 1286.0
## settlements 1 1288.2
## cos_td_rad 1 1291.1
## sin_td_rad 1 1291.9me1 <- lmer(data = P.anal, formula = gma ~ settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements + cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 1286.0
## settlements 1 1288.2
## cos_td_rad 1 1291.1
## sin_td_rad 1 1291.9This is the final model:
summary(me1)## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements + cos_td_rad + sin_td_rad + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 1281.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0616 -0.6198 -0.1484 0.4293 5.8897
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.06892 0.2625
## Residual 1.02964 1.0147
## Number of obs: 440, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.1545 0.5453 2.117
## settlementsNear 0.2002 0.0968 2.069
## cos_td_rad 1.1108 0.4160 2.670
## sin_td_rad -1.2330 0.4369 -2.823
##
## Correlation of Fixed Effects:
## (Intr) sttlmN cs_td_
## settlmntsNr -0.088
## cos_td_rad -0.964 0.000
## sin_td_rad 0.957 -0.002 -0.904ranef(me1)## $site
## (Intercept)
## Abirim 0.008982551
## Aderet 0.378982346
## Beit Oren 0.277725813
## Ein Yaakov 0.005273928
## Givat Yearim -0.367444219
## Givat Yeshayahu 0.158173367
## Goren 0.162039344
## Iftach -0.305765485
## Kerem Maharal -0.030632466
## Kfar Shamai -0.137142196
## Margaliot -0.074170019
## Nehusha -0.154001749
## Nir Etzion -0.149171842
## Ofer 0.038757806
## Ramat Raziel -0.079372801
## Yagur 0.267765624
##
## with conditional variances for "site"dotplot(ranef(me1, condVar=TRUE))## $site
plot(me1)
qqmath(me1)
# scale location plot
plot(me1, sqrt(abs(resid(.))) ~ fitted(.),
type = c("p", "smooth"),
main="scale-location",
par.settings = list(plot.line =
list(alpha=1, col = "red",
lty = 1, lwd = 2)))
Not a great fit. settlement is significant.
summ(me1, digits = 3)| Observations | 440 |
| Dependent variable | gma |
| Type | Mixed effects linear regression |
| AIC | 1293.430 |
| BIC | 1317.950 |
| Pseudo-R² (fixed effects) | 0.028 |
| Pseudo-R² (total) | 0.089 |
| Est. | S.E. | t val. | d.f. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.155 | 0.549 | 2.102 | 406.354 | 0.036 |
| settlementsNear | 0.200 | 0.097 | 2.068 | 421.708 | 0.039 |
| cos_td_rad | 1.111 | 0.419 | 2.651 | 427.380 | 0.008 |
| sin_td_rad | -1.233 | 0.440 | -2.800 | 416.685 | 0.005 |
| p values calculated using Kenward-Roger standard errors and d.f. |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.263 |
| Residual | 1.015 |
| Group | # groups | ICC |
|---|---|---|
| site | 16 | 0.063 |
plot(me1, which=1:3)
effplot1 <- plot_model_effect(P.anal = P.anal, m = me1, eff2plot = "settlements", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/", y_cutoff = 5)## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.## Warning: Removed 16 rows containing missing values (`geom_point()`).
## Removed 16 rows containing missing values (`geom_point()`).
effplot1## Warning: Removed 16 rows containing missing values (`geom_point()`).
perc_diff <- diff(as.data.table(effplot1$data)[,gma])/as.data.table(effplot1$data)[settlements=="Far",gma]*100
# EMM <- emmeans(object = me1, ~ settlements)
# print(EMM)
# test_results_settlements <- test(pairs(EMM, by=NULL, adjust="fdr"))
# print(test_results_settlements)On average, GMA in near plots is higher by 7.5128056 percent.
abundance
Explore data
# exclude rare species in analysis
P.anal <- copy(P_no_rare)
print("ABUNDANCE WITHOUT RARE SPECIES")## [1] "ABUNDANCE WITHOUT RARE SPECIES"plot_alpha_diversity(P.anal, x_val = "settlements", y_val = "abundance", ylab_val = "abundance", xlab_val = "settlements")
plot_alpha_diversity(P.anal, x_val = "subunit", y_val = "abundance", ylab_val = "abundance", xlab_val = "subunit")
dotchart(P.anal$abundance, groups = P.anal$settlements, pch = as.numeric(P.anal$settlements), ylab = "settlements", xlab = "abundance")
dotchart(P.anal$abundance, groups = P.anal$subunit, pch = as.numeric(P.anal$subunit), ylab = "subunit", xlab = "abundance")
IVs <- c("abundance", "year_ct","site","settlements","subunit", "td_sc", "cos_td_rad", "sin_td_rad")
kable(summary(P.anal[,..IVs]))%>% kable_styling(font_size = 11)| abundance | year_ct | site | settlements | subunit | td_sc | cos_td_rad | sin_td_rad | |
|---|---|---|---|---|---|---|---|---|
| Min. : 0.00 | Min. :0.0 | Nir Etzion : 31 | Far :224 | Judean Highlands:149 | Min. :-1.9882 | Min. :-0.1671 | Min. :-1.0000 | |
| 1st Qu.: 15.75 | 1st Qu.:3.0 | Ein Yaakov : 30 | Near:220 | Carmel :150 | 1st Qu.:-0.2713 | 1st Qu.: 0.5702 | 1st Qu.:-0.8215 | |
| Median : 24.00 | Median :5.0 | Givat Yearim : 30 | NA | Galilee :145 | Median : 0.4918 | Median : 0.8140 | Median :-0.5808 | |
| Mean : 30.45 | Mean :4.8 | Givat Yeshayahu: 30 | NA | NA | Mean : 0.3205 | Mean : 0.7098 | Mean :-0.5864 | |
| 3rd Qu.: 36.00 | 3rd Qu.:7.0 | Goren : 30 | NA | NA | 3rd Qu.: 1.1023 | 3rd Qu.: 0.9413 | 3rd Qu.:-0.3375 | |
| Max. :345.00 | Max. :9.0 | Kerem Maharal : 30 | NA | NA | Max. : 1.7127 | Max. : 0.9976 | Max. :-0.0688 | |
| NA | NA | (Other) :263 | NA | NA | NA | NA | NA |
Some outliers with total abundance >100. Examine:
P.anal[abundance>100,.(subunit,point_name,datetime,monitors_name,richness,gma,abundance)]## subunit point_name datetime monitors_name
## 1: Judean Highlands Givat Yeshayahu Far 2 2017-04-20 10:20:00 Eran Banker
## 2: Judean Highlands Givat Yeshayahu Near 1 2017-04-20 06:25:00 Eran Banker
## 3: Judean Highlands Givat Yeshayahu Near 2 2017-04-20 06:50:00 Eran Banker
## 4: Judean Highlands Givat Yeshayahu Near 3 2017-04-20 07:15:00 Eran Banker
## 5: Judean Highlands Aderet Near 2 2019-05-26 06:52:00 Eran Banker
## 6: Carmel Beit Oren Near 1 2019-04-19 07:49:00 Other
## 7: Carmel Beit Oren Near 3 2019-04-19 07:35:00 Other
## 8: Carmel Kerem Maharal Near 3 2021-04-17 07:41:00 Eliraz Dvir
## 9: Carmel Nir Etzion Far 11 2021-04-23 07:24:00 Eliraz Dvir
## richness gma abundance
## 1: 11 5.356635 115
## 2: 13 9.170560 140
## 3: 14 7.652753 138
## 4: 16 9.098580 221
## 5: 15 5.408017 106
## 6: 5 8.043088 129
## 7: 10 7.338970 254
## 8: 11 10.301194 345
## 9: 9 3.066833 119Exclude 3 plots with high abundance (>150) to improve model fit.
P.anal <- P.anal[abundance<150]
mdl_a.po <- glm(data = P.anal, formula = abundance ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad, family = poisson)
od <- mdl_a.po$deviance/mdl_a.po$df.residual
print("Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model.")
print(paste0('od = ',od))
plot(mdl_a.po,which = 1:2,sub.caption = "poisson")PHI>1, hence choose negative binomial. Fit fixed and mixed models. Choose mixed model if possible, otherwise choose a model with fixed-effects only.
m0 <- glm.nb(data = P.anal, formula = abundance ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad)
AIC(m0)## [1] 3689.581me0 <- glmer.nb(data = P.anal, formula = abundance ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad + (1|site))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.004264 (tol = 0.002, component 1)AIC(me0)## [1] 3681.608mixed model converged, slightly better AIC. Perform stepwise model selection of mixed model.
drop1(me0)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00218393 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## abundance ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3681.6
## subunit 2 3687.1
## cos_td_rad 1 3697.1
## sin_td_rad 1 3699.5
## settlements:year_ct 1 3679.7drop settlements X year because \(\Delta AIC<2\)
me1 <- glmer.nb(data = P.anal, formula = abundance ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## abundance ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3679.7
## settlements 1 3773.3
## subunit 2 3685.2
## year_ct 1 3678.2
## cos_td_rad 1 3695.2
## sin_td_rad 1 3697.6drop year
me1 <- glmer.nb(data = P.anal, formula = abundance ~ settlements + subunit + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## abundance ~ settlements + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3678.2
## settlements 1 3772.0
## subunit 2 3683.8
## cos_td_rad 1 3693.2
## sin_td_rad 1 3697.8Subunit, settlement and time of year remain. The final model:
summary(me1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(3.3101) ( log )
## Formula: abundance ~ settlements + subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 3678.2 3711.0 -1831.1 3662.2 436
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7208 -0.7138 -0.2446 0.3842 9.2005
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01938 0.1392
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67516 0.32509 5.153 2.57e-07 ***
## settlementsNear 0.54579 0.05556 9.823 < 2e-16 ***
## subunitCarmel 0.16956 0.11089 1.529 0.1262
## subunitGalilee -0.22714 0.11015 -2.062 0.0392 *
## cos_td_rad 1.02590 0.24294 4.223 2.41e-05 ***
## sin_td_rad -1.16744 0.24800 -4.708 2.51e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl cs_td_
## settlmntsNr -0.111
## subunitCrml -0.214 0.020
## subunitGall -0.226 0.033 0.505
## cos_td_rad -0.946 0.019 0.050 0.043
## sin_td_rad 0.938 -0.017 -0.032 -0.062 -0.903ranef(me1)## $site
## (Intercept)
## Abirim -0.11137155
## Aderet 0.16235427
## Beit Oren 0.11023799
## Ein Yaakov 0.01625407
## Givat Yearim -0.19886986
## Givat Yeshayahu 0.17665857
## Goren -0.01370126
## Iftach 0.07733150
## Kerem Maharal 0.13673950
## Kfar Shamai 0.00869344
## Margaliot 0.02457356
## Nehusha -0.07214524
## Nir Etzion -0.15864707
## Ofer -0.04190665
## Ramat Raziel -0.06677795
## Yagur -0.04538202
##
## with conditional variances for "site"dotplot(ranef(me1, condVar=TRUE))## $site
plot(me1)
qqmath(me1)
# scale location plot
plot(me1, sqrt(abs(resid(.))) ~ fitted(.),
type = c("p", "smooth"),
main="scale-location",
par.settings = list(plot.line =
list(alpha=1, col = "red",
lty = 1, lwd = 2)))
Not a great fit, high residuals. center time of year variables, highly correlated with intercept.
P.anal[,`:=`(cos_td_rad_c=cos_td_rad-mean(cos_td_rad),
sin_td_rad_c=sin_td_rad-mean(sin_td_rad))]
me1 <- glmer.nb(data = P.anal, formula = abundance ~ settlements + subunit + cos_td_rad_c + sin_td_rad_c + (1|site))
summary(me1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(3.3101) ( log )
## Formula: abundance ~ settlements + subunit + cos_td_rad_c + sin_td_rad_c +
## (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 3678.2 3711.0 -1831.1 3662.2 436
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7208 -0.7138 -0.2446 0.3842 9.2005
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01938 0.1392
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.08791 0.08411 36.711 < 2e-16 ***
## settlementsNear 0.54578 0.05556 9.823 < 2e-16 ***
## subunitCarmel 0.16956 0.11088 1.529 0.1262
## subunitGalilee -0.22714 0.11015 -2.062 0.0392 *
## cos_td_rad_c 1.02590 0.24293 4.223 2.41e-05 ***
## sin_td_rad_c -1.16743 0.24798 -4.708 2.51e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl cs_t__
## settlmntsNr -0.358
## subunitCrml -0.669 0.020
## subunitGall -0.678 0.033 0.505
## cos_td_rd_c -0.047 0.019 0.050 0.043
## sin_td_rd_c 0.044 -0.017 -0.032 -0.062 -0.903Interpretation of abundance model:
summ(me1,exp=T,digits=3)| Observations | 444 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | Negative Binomial(3.3101) |
| Link | log |
| AIC | 3678.232 |
| BIC | 3710.998 |
| Pseudo-R² (fixed effects) | 0.692 |
| Pseudo-R² (total) | 0.804 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 21.931 | 0.084 | 36.711 | 0.000 |
| settlementsNear | 1.726 | 0.056 | 9.823 | 0.000 |
| subunitCarmel | 1.185 | 0.111 | 1.529 | 0.126 |
| subunitGalilee | 0.797 | 0.110 | -2.062 | 0.039 |
| cos_td_rad_c | 2.790 | 0.243 | 4.223 | 0.000 |
| sin_td_rad_c | 0.311 | 0.248 | -4.708 | 0.000 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.139 |
| Group | # groups | ICC |
|---|---|---|
| site | 16 | 0.052 |
effplot1 <- effect_plot(model = me1, data=P.anal, pred = settlements, interval = T, plot.points = T, jitter = c(0.1,0), point.size = 3, cat.pred.point.size = 4, partial.residuals = F, colors = "Qual1", main.title = "effect of settlement proximity (points are observations)", line.colors = "black", point.alpha = 0.25) + ylim(0,150)## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.effplot1
plot_model_effect(P.anal = P.anal, m = me1, eff2plot = "settlements", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/", y_cutoff = 80)## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
#Far vs Near - near is 72.59% higher than far
print(effplot1$data$abundance)## 1 2
## 21.93111 37.85224iv1_abs <- diff(effplot1$data$abundance)
iv1_rel <- iv1_abs / as.data.table(effplot1$data)[settlements=="Far",abundance] * 100
print(iv1_rel)## 2
## 72.5961plot_model_effect(P.anal = P.anal, m = me1, eff2plot = "subunit", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/", y_cutoff = 80)## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
#emmeans
EMM <- emmeans(object = me1, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 3.36 0.0786 Inf 3.21 3.51
## Carmel 3.53 0.0782 Inf 3.38 3.68
## Galilee 3.13 0.0771 Inf 2.98 3.28
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95# Carmel vs Galilee
(exp(3.53)/exp(3.13) - 1) * 100## [1] 49.18247# Carmel vs Judea
(exp(3.53)/exp(3.36) - 1) * 100## [1] 18.53049#Judea vs Galilee
(exp(3.36)/exp(3.13) - 1) * 100## [1] 25.86test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -0.170 0.111 Inf -1.529 0.1262
## Judean Highlands - Galilee 0.227 0.110 Inf 2.062 0.0588
## Carmel - Galilee 0.397 0.110 Inf 3.607 0.0009
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testssignificantly lower abundance in Galilee subunit and far from settlements. No significant temporal trend.
community analysis using package MVabund
abu_by_spp.maquis <- abu_by_spp[grepl("Maquis",unit),]
abu_by_spp_no_rare.maquis <- abu_by_spp_no_rare[grepl("Maquis",unit),]
spp <- abu_by_spp.maquis[,.SD[,(length(col_names)+1):ncol(abu_by_spp.maquis)]]
# filter out species with zero counts (were either not observed or are rare and were excluded)
spp <- spp[,.SD,.SDcols = colSums(spp)>0]
spp <- mvabund(spp)
spp_no_rare <- abu_by_spp_no_rare.maquis[,.SD[,(length(col_names)+1):ncol(abu_by_spp_no_rare.maquis)]]
# filter out species with zero counts (were either not observed or are rare and were excluded)
spp_no_rare <- spp_no_rare[,.SD,.SDcols = colSums(spp_no_rare)>0]
spp_no_rare <- mvabund(spp_no_rare)
env_data <- abu_by_spp_no_rare.maquis[pilot==FALSE,..col_names]
env_data[,site:=factor(site)][,subunit:=factor(subunit)]
try(expr = plot(spp ~ abu_by_spp.maquis$settlements,overall.main="raw abundances", transformation = "no"))## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Curruca.melanocephala, Turdus.merula, Streptopelia.decaocto, Passer.domesticus, Pycnonotus.xanthopygos, Columba.livia, Parus.major, Garrulus.glandarius, Cinnyris.osea, Spilopelia.senegalensis, Cecropis.daurica, Chloris.chloris were included in the plot
## (the variables with highest total abundance).
try(plot(spp ~ abu_by_spp.maquis$subunit,overall.main="raw abundances", transformation = "no"))## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Curruca.melanocephala, Turdus.merula, Streptopelia.decaocto, Passer.domesticus, Pycnonotus.xanthopygos, Columba.livia, Parus.major, Garrulus.glandarius, Cinnyris.osea, Spilopelia.senegalensis, Cecropis.daurica, Chloris.chloris were included in the plot
## (the variables with highest total abundance).
plot(sort(spp_no_rare[spp_no_rare>0]))
dotchart(sort(A[grepl("Maquis",unit),count_under_250])) # these are the actual observations, before aggregation into plots
dotchart(sort(A[grepl("Maquis",unit) & count_under_250>=30,count_under_250]))
meanvar.plot(spp, xlab = "mean abundance of a given species across sites", ylab = "variance of the abundance of a given species across sites")
There are few observations with counts of >60. Examine these:
A[grepl("Maquis",unit) & count_under_250>=60,.(point_name,datetime,SciName,monitors_name,notes,rad_0_20,rad_20_100,rad_100_250,rad_over_250)]## point_name datetime SciName monitors_name
## 1: Kerem Maharal Near 3 2021-04-17 07:41:00 Chloris chloris Eliraz Dvir
## 2: Kerem Maharal Near 3 2021-04-17 07:41:00 Garrulus glandarius Eliraz Dvir
## 3: Beit Oren Near 3 2019-04-19 07:35:00 Curruca curruca Other
## 4: Beit Oren Near 1 2019-04-19 07:49:00 Curruca curruca Other
## notes rad_0_20 rad_20_100 rad_100_250 rad_over_250
## 1: 0 80 0 0
## 2: 0 0 150 0
## 3: <NA> 12 30 100 0
## 4: <NA> 14 35 50 0G. glandarius observation seems highly irregular. Remove altogether, might be a mistake. C. curruca observations are probably due to migration waves. Remove.
spp_no_rare[env_data[point_name=="Kerem Maharal Near 3" & year==2021,which = TRUE],"Garrulus.glandarius"] <- 0
spp_no_rare[env_data[point_name=="Beit Oren Near 3" & year==2019,which = TRUE],"Curruca.curruca"] <- 0
spp_no_rare[env_data[point_name=="Beit Oren Near 1" & year==2019,which = TRUE],"Curruca.curruca"] <- 0start model specification:
mva_m0.po <- manyglm(formula = spp_no_rare ~ settlements*year_ct + subunit + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
mva_m0.nb <- manyglm(formula = spp_no_rare ~ settlements*year_ct + subunit + cos_td_rad + sin_td_rad, 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
## 978.686 1390.834print("POISSON")## [1] "POISSON"plot(mva_m0.po, which=1:3)
print("NEGATIVE BINOMIAL")## [1] "NEGATIVE BINOMIAL"plot(mva_m0.nb, which=1:3)
negative binomial model is better than poisson according to residuals and AIC comparison.
mva_m0 <- mva_m0.nb
mva_m1 <- manyglm(formula = spp_no_rare ~ settlements*year_ct + subunit + cos_td_rad + sin_td_rad + site, data = env_data)
aic_dt$nb1 <- mva_m1$aic
colMeans(aic_dt)## nb po nb1
## 978.6860 1390.8344 968.2392The addition of the explanatory variable ‘site’ is somewhat improving the AIC of the model. Prefer to exclude site, for simplification. stepwise selection of model:
drop1(mva_m0)## Single term deletions
##
## Model:
## spp_no_rare ~ settlements * year_ct + subunit + cos_td_rad +
## sin_td_rad
## Df AIC
## <none> 23489
## subunit 48 23844
## cos_td_rad 24 23537
## sin_td_rad 24 23549
## settlements:year_ct 24 23481drop settlements X year.
mva_m1 <- manyglm(formula = spp_no_rare ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, family = "negative.binomial", data = env_data)
drop1(mva_m1)## Single term deletions
##
## Model:
## spp_no_rare ~ settlements + subunit + year_ct + cos_td_rad +
## sin_td_rad
## Df AIC
## <none> 23481
## settlements 24 24269
## subunit 48 23838
## year_ct 24 23502
## cos_td_rad 24 23528
## sin_td_rad 24 23542final model includes settlements, year, subunit, sampling time of year.
plot(mva_m1, which=1:3)
#summ_m1 <- summary(mva_m1)
#saveRDS(summ_m1, file = "../output/Maquis_new_MVabund_summary_output_999iter.rds")
#anov_m1.uni <- anova(mva_m1, p.uni = "adjusted", nBoot = 99, show.time = "all")
#saveRDS(anov_m1.uni, file = "../output/Maquis_new_MVabund_anova_output_99iter.rds")
summ_m1 <- readRDS(file = "../output/Maquis_new_MVabund_summary_output_999iter.rds")
anov_m1.uni <- readRDS("../output/Maquis_new_MVabund_anova_output_99iter.rds")
print(summ_m1)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 12.055 0.001 ***
## settlementsNear 28.826 0.001 ***
## subunitCarmel 14.027 0.001 ***
## subunitGalilee 15.305 0.001 ***
## year_ct 8.913 0.001 ***
## cos_td_rad 10.020 0.001 ***
## sin_td_rad 10.820 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 38.38, 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_m1.uni)## Analysis of Deviance Table
##
## Model: spp_no_rare ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 443
## settlements 442 1 770.4 0.01 **
## subunit 440 2 426.6 0.01 **
## year_ct 439 1 143.4 0.01 **
## cos_td_rad 438 1 46.9 0.04 *
## sin_td_rad 437 1 109.1 0.01 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acridotheres.tristis Alectoris.chukar
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 79.707 0.01 38.366 0.01
## subunit 39.138 0.01 10.607 0.10
## year_ct 29.844 0.01 9.298 0.07
## cos_td_rad 3.002 0.87 2.276 0.95
## sin_td_rad 5.64 0.41 0.596 1.00
## Carduelis.carduelis Cecropis.daurica
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 8.13 0.04 33.414 0.01
## subunit 21.019 0.01 1.066 0.94
## year_ct 6.527 0.29 1.257 0.96
## cos_td_rad 0.389 1.00 0.127 1.00
## sin_td_rad 1.61 0.95 3.044 0.78
## Chloris.chloris Cinnyris.osea Columba.livia
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## settlements 30.55 0.01 41.948 0.01 6.276
## subunit 0.613 0.94 17.867 0.02 19.317
## year_ct 5.072 0.46 3.939 0.62 20.916
## cos_td_rad 0.072 1.00 1.651 0.98 0.01
## sin_td_rad 0.112 1.00 2.01 0.93 11.564
## Corvus.cornix Corvus.monedula
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 0.11 22.932 0.01 0.428 0.58
## subunit 0.01 47.477 0.01 22.714 0.01
## year_ct 0.01 0.695 0.98 0.301 1.00
## cos_td_rad 1.00 1.013 1.00 0.184 1.00
## sin_td_rad 0.02 0.157 1.00 0.053 1.00
## Curruca.curruca Curruca.melanocephala
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 2.096 0.58 21.24 0.01
## subunit 27.501 0.01 41.067 0.01
## year_ct 13.288 0.02 0.747 0.98
## cos_td_rad 7.65 0.23 0.3 1.00
## sin_td_rad 4.442 0.56 0.429 1.00
## Dendrocopos.syriacus Falco.tinnunculus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 9.285 0.03 1.849 0.58
## subunit 3.648 0.66 9.122 0.11
## year_ct 0.002 1.00 0.543 0.98
## cos_td_rad 1.417 1.00 0.078 1.00
## sin_td_rad 3.35 0.71 0.641 1.00
## Garrulus.glandarius Parus.major Passer.domesticus
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## settlements 3.502 0.44 3.908 0.44 166.539
## subunit 24.625 0.01 25.606 0.01 21.305
## year_ct 4.924 0.46 0.176 1.00 1.781
## cos_td_rad 1.156 1.00 3.376 0.84 0.485
## sin_td_rad 1.081 0.99 7.701 0.17 2.962
## Prinia.gracilis Psittacula.krameri
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 0.01 48.787 0.01 15.594 0.01
## subunit 0.01 1.864 0.90 14.393 0.05
## year_ct 0.92 16.036 0.02 1.082 0.96
## cos_td_rad 1.00 0.633 1.00 0.793 1.00
## sin_td_rad 0.80 1.004 0.99 1.724 0.95
## Pycnonotus.xanthopygos Spilopelia.senegalensis
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 45.339 0.01 157.013 0.01
## subunit 25.29 0.01 27.046 0.01
## year_ct 2.87 0.78 2.224 0.88
## cos_td_rad 11.394 0.05 0.333 1.00
## sin_td_rad 0.499 1.00 0.603 1.00
## Streptopelia.decaocto Streptopelia.turtur
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 1.955 0.58 1.175 0.58
## subunit 12.168 0.09 6.566 0.33
## year_ct 0.068 1.00 18.089 0.02
## cos_td_rad 6.672 0.34 3.045 0.87
## sin_td_rad 14.008 0.01 1.847 0.93
## Troglodytes.troglodytes Turdus.merula
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 26.472 0.01 3.926 0.44
## subunit 4.411 0.62 2.142 0.90
## year_ct 0.576 0.98 3.137 0.75
## cos_td_rad 0.48 1.00 0.347 1.00
## sin_td_rad 10.786 0.03 33.19 0.01
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 99 iterations via PIT-trap resampling.Factors settlements, year, subunit, time of year have a statistically significant effect on community composition.
coefp <- merge(data.table(t(coef(mva_m1)/log(2)),keep.rownames=TRUE), data.table(t(anov_m1.uni$uni.p), keep.rownames = TRUE), by = "rn") # coefficients are on log2 scale to simplify interpretation, yet keep presenetation on log scale rather than linear. The link function used for poisson is log -> above 0 is a positive trend and below 0 is negative trend
# add total species abundance
coefp <- merge(coefp,as.data.table(colSums(spp_no_rare),keep.rownames = TRUE), by.x = "rn", by.y = "V1")
# add hebrew name and traits
sci_heb <- unique(A[,.(SciName,HebName)])
sci_heb[,SciName:=make.names(SciName)]
coefp <- merge(coefp, sci_heb, by.x = "rn", by.y = "SciName")
# set column names based on the predictor variables
colnames(coefp)[grepl("rn",colnames(coefp))] <- "SciName"
colnames(coefp)[grepl("Intercept).x",colnames(coefp))] <- "intercept.coef"
colnames(coefp)[grepl("Intercept).y",colnames(coefp))] <- "intercept.p"
colnames(coefp)[colnames(coefp)=="settlementsNear"] <- "settlementsNear.coef"
colnames(coefp)[colnames(coefp)=="settlements"] <- "settlements.p"
colnames(coefp)[grepl("cos_td_rad.x",colnames(coefp))] <- "cos_td_rad.coef"
colnames(coefp)[colnames(coefp)=="cos_td_rad.y"] <- "cos_td_rad.p"
colnames(coefp)[grepl("sin_td_rad.x",colnames(coefp))] <- "sin_td_rad.coef"
colnames(coefp)[grepl("sin_td_rad.y",colnames(coefp))] <- "sin_td_rad.p"
colnames(coefp)[colnames(coefp)=="year_ct.x"] <- "year_ct.coef"
colnames(coefp)[colnames(coefp)=="year_ct.y"] <- "year_ct.p"
colnames(coefp)[colnames(coefp)=="V2"] <- "species_abundance"
colnames(coefp)[colnames(coefp)=="dunesshifting"] <- "dunesShifting.coef"
colnames(coefp)[colnames(coefp)=="dunes"] <- "dunes.p"
colnames(coefp)[colnames(coefp)=="site"] <- "site.p"
colnames(coefp)[colnames(coefp)=="subunitCarmel"] <- "subunitCarmel.coef"
colnames(coefp)[colnames(coefp)=="subunitGalilee"] <- "subunitGalilee.coef"
colnames(coefp)[colnames(coefp)=="subunit"] <- "subunit.p"
# plot_coefs_settle <- plot_coefs_with_traits(D = A_no_rare[grepl("Maquis",unit)],
# cp = coefp[,.(SciName,settlementsNear.coef,settlements.p,species_abundance)],
# plot_title = "Effect of proximity to settlements on species abundance",
# plot_xlabel = expression(paste("coefficient of NEAR settlements, (log"[2]," scale)")),
# mark_batha = FALSE)
# plot_coefs_settle
plot_coefs_settle <- plot_coefs_with_traits_for_publication(D = A_no_rare[grepl("Maquis",unit)],
cp = coefp[,.(SciName,settlementsNear.coef,settlements.p,species_abundance)],
plot_title = "Effect of proximity to settlement on species abundance",
plot_xlabel = expression(paste("coefficient of NEAR settlement, (log"[2]," scale)")),
mark_batha = TRUE,
show_legend = FALSE,
export_plot = TRUE,
fontname = "Almoni ML v5 AAA",
effect_str = "settlements",
outpath = "../output/maquis_new/")## Loading required package: ggnewscale## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## [1] "black arrow is median of synanthrope/invasive\n\n"plot_coefs_settle
# compute mean of year coefficients (geometric mean of linear scale coefficients) to plot single year coefficient
#coefp[,mn_year.coef:=year_ct.coef+`subunitCarmel:year_ct`/3+`subunitGalilee:year_ct`/3]
# if any of the relevant terms - year_ct, subunit X year - is <0.1, set the p-value of the mean year coef to be 0.01 (arbitrary, lower than 0.1 which is the threshold)
#coefp[,mn_year.p:=apply(X = .SD, 1, FUN = \(x) 1-sum(x<0.1)*0.99), .SDcols=c("year_ct.p","subunit:year_ct")]
# plot_coefs_year <- plot_coefs_with_traits(D = A_no_rare[grepl("Maquis",unit)],
# cp = coefp[,.(SciName,mn_year.coef,mn_year.p,species_abundance)],
# plot_title = "Mean temporal effect on species abundance",
# plot_xlabel = expression(paste("coefficient of year counter, (log"[2]," scale)")),
# mark_batha = FALSE)
plot_coefs_year <- plot_coefs_with_traits_for_publication(D = A_no_rare[grepl("Maquis",unit)],
cp = coefp[,.(SciName,year_ct.coef,year_ct.p,species_abundance)],
plot_title = "Temporal effect on species abundance",
plot_xlabel = expression(paste("coefficient of year counter, (log"[2]," scale)")),
mark_batha = TRUE,
show_legend = FALSE,
export_plot = TRUE,
fontname = "Almoni ML v5 AAA",
effect_str = "year_ct",
outpath = "../output/maquis_new/")## [1] "black arrow is median of synanthrope/invasive\n\n"plot_coefs_year
plot_coefs_subunit <- plot_coefs_with_traits_for_publication(D = A_no_rare[grepl("Maquis",unit)],
cp = coefp[,.(SciName,subunitCarmel.coef,subunit.p,species_abundance)],
plot_title = "Effect of sub-unit on species abundance",
plot_xlabel = expression(paste("coefficient of Carmel sub-unit, (log"[2]," scale)")),
mark_batha = TRUE,
show_legend = FALSE,
export_plot = TRUE,
fontname = "Almoni ML v5 AAA",
effect_str = "subunit",
outpath = "../output/maquis_new/")## [1] "black arrow is median of synanthrope/invasive\n\n"plot_coefs_subunit
# plot the coefficients of year considering interaction with subunits
# coef_cutoff <- 0.25 # do not display species with all three coefs below this
# X <- coefp[,.(Judea=year_ct.coef,
# Carmel=year_ct.coef + `subunitCarmel:year_ct`,
# Galil=year_ct.coef + `subunitGalilee:year_ct`,
# species_abundance,
# SciName,HebName)]
# Xcut <- X[abs(Judea)>coef_cutoff | abs(Carmel)>coef_cutoff | abs(Galil)>coef_cutoff]
# col_scale_min <- -0.7
# p_yearsubu <- ggplot(data = Xcut, aes(x = Galil, y = Carmel, color=Judea)) +
# theme_light() +
# geom_hline(yintercept = 0, color = "darkslategrey") +
# geom_vline(xintercept = 0, color = "darkslategrey") +
# geom_point(size=3.5) +
# scale_color_gradient2(midpoint=0, low="darkred", mid="white",high="darkblue", space ="Lab",
# limits=c(col_scale_min,NA), oob=scales::squish, name=paste0("Judea (lower limit=",col_scale_min,")")) +
# geom_text_repel(aes(label = paste0("(",species_abundance,") ",HebName)), color="black", size=3.5, force_pull = 0, max.overlaps = 50,
# max.iter = 1e5, box.padding = 0.7, min.segment.length = 0) +
# ggtitle("Year coefficient for three Mediterranean Maquis subunits")
#
# p_yearsubu##########Alectoris chukar - חוגלת סלעים#########
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Alectoris chukar", ylab_val = "Alectoris chukar")
d.alechu <- cbind(env_data,data.table(Alectoris.chukar = spp_no_rare[,"Alectoris.chukar"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.alechu <- glm.nb(formula = Alectoris.chukar~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.alechu)
coef(mva_m1)[,"Alectoris.chukar"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -3.3297204 -1.8142174 1.0363737 0.3667818 0.1915418
## cos_td_rad sin_td_rad
## 1.7562492 -1.0706266coef(m.alechu)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -3.3297134 -1.8142178 1.0363750 0.3667830 0.1915409
## cos_td_rad sin_td_rad
## 1.7562446 -1.0706260alechu_time <- effect_plot(data = d.alechu, model = m.alechu, pred = year_ct, 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 = "")
alechu_time
#Temporal trend - increase of 460.61% total
time_dt <- as.data.table(alechu_time$data$Alectoris.chukar)
(max(time_dt)/min(time_dt) - 1) * 100## [1] 460.6174alechu_time$data$Alectoris.chukar## [1] 0.2333147 0.2374129 0.2415831 0.2458266 0.2501446 0.2545385 0.2590095
## [8] 0.2635591 0.2681886 0.2728994 0.2776930 0.2825708 0.2875342 0.2925848
## [15] 0.2977242 0.3029538 0.3082752 0.3136902 0.3192003 0.3248071 0.3305124
## [22] 0.3363180 0.3422255 0.3482368 0.3543537 0.3605780 0.3669117 0.3733566
## [29] 0.3799147 0.3865880 0.3933786 0.4002884 0.4073196 0.4144743 0.4217546
## [36] 0.4291629 0.4367012 0.4443720 0.4521776 0.4601202 0.4682023 0.4764265
## [43] 0.4847950 0.4933106 0.5019758 0.5107931 0.5197653 0.5288952 0.5381854
## [50] 0.5476388 0.5572582 0.5670466 0.5770070 0.5871423 0.5974556 0.6079501
## [57] 0.6186289 0.6294953 0.6405526 0.6518041 0.6632532 0.6749035 0.6867583
## [64] 0.6988215 0.7110965 0.7235871 0.7362971 0.7492304 0.7623908 0.7757825
## [71] 0.7894093 0.8032755 0.8173853 0.8317429 0.8463527 0.8612192 0.8763468
## [78] 0.8917400 0.9074037 0.9233425 0.9395613 0.9560650 0.9728586 0.9899471
## [85] 1.0073359 1.0250300 1.0430350 1.0613562 1.0799992 1.0989697 1.1182735
## [92] 1.1379163 1.1579041 1.1782431 1.1989393 1.2199990 1.2414286 1.2632347
## [99] 1.2854238 1.30800271.3080027/0.2333147 - 1## [1] 4.606174alechu_set <- effect_plot(data = d.alechu, model = m.alechu, pred = settlements, 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 = "")
#abundance far from settlements is 513.62% higher than near
set_dt <- as.data.table(alechu_set$data$Alectoris.chukar)
(max(set_dt)/min(set_dt) - 1) *100## [1] 513.62740.58505058/0.09534296 -1## [1] 5.136275#or 5.88 fold
0.6002156 / 0.1019473## [1] 5.887509####Acridotheres tristis - מיינה מצויה####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Acridotheres tristis", ylab_val = "Acridotheres tristis")
d.acrtri <- cbind(env_data,data.table(Acridotheres.tristis = spp_no_rare[,"Acridotheres.tristis"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.acrtri <- glm.nb(formula = Acridotheres.tristis ~ settlements + year_ct + subunit + cos_td_rad + sin_td_rad, data = d.acrtri)
coef(mva_m1)[,"Acridotheres.tristis"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -5.81322891 3.04949199 -0.06360271 -2.84839762 0.14854946
## cos_td_rad sin_td_rad
## 1.04586057 -2.82696864coef(m.acrtri)## (Intercept) settlementsNear year_ct subunitCarmel subunitGalilee
## -5.81322905 3.04949186 0.14854959 -0.06360287 -2.84839757
## cos_td_rad sin_td_rad
## 1.04586061 -2.82696797# plot_model_effect(P.anal = d.acrtri, m = m.acrtri, eff2plot = "year_ct", export_plot = TRUE, plot_points = FALSE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")
tristis_time_plot <- effect_plot(data = d.acrtri, model = m.acrtri, pred = year_ct, 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 trend - increase of 280.73% total
tristis_time <- as.data.table(tristis_time_plot$data$Acridotheres.tristis)
(max(tristis_time)/min(tristis_time) - 1) * 100## [1] 280.7399tristis_time_plot$data$Acridotheres.tristis## [1] 0.03293719 0.03338501 0.03383892 0.03429899 0.03476533 0.03523800
## [7] 0.03571710 0.03620271 0.03669493 0.03719384 0.03769953 0.03821209
## [13] 0.03873163 0.03925823 0.03979199 0.04033300 0.04088138 0.04143720
## [19] 0.04200059 0.04257163 0.04315044 0.04373712 0.04433177 0.04493451
## [25] 0.04554545 0.04616469 0.04679235 0.04742854 0.04807338 0.04872700
## [31] 0.04938949 0.05006100 0.05074163 0.05143152 0.05213079 0.05283956
## [37] 0.05355798 0.05428616 0.05502424 0.05577235 0.05653064 0.05729924
## [43] 0.05807828 0.05886792 0.05966830 0.06047955 0.06130184 0.06213530
## [49] 0.06298010 0.06383639 0.06470431 0.06558404 0.06647573 0.06737954
## [55] 0.06829564 0.06922419 0.07016537 0.07111935 0.07208629 0.07306639
## [61] 0.07405981 0.07506673 0.07608735 0.07712184 0.07817039 0.07923321
## [67] 0.08031047 0.08140238 0.08250913 0.08363094 0.08476799 0.08592051
## [73] 0.08708869 0.08827276 0.08947292 0.09068941 0.09192243 0.09317222
## [79] 0.09443900 0.09572300 0.09702446 0.09834361 0.09968070 0.10103597
## [85] 0.10240967 0.10380204 0.10521335 0.10664384 0.10809378 0.10956343
## [91] 0.11105307 0.11256296 0.11409338 0.11564460 0.11721692 0.11881061
## [97] 0.12042597 0.12206330 0.12372288 0.125405030.12540503/0.03293719 - 1## [1] 2.807399# or 3.16 fold increase
2.4589569/0.7779155## [1] 3.160956# tristis_set <- plot_model_effect(P.anal = d.acrtri, m = m.acrtri, eff2plot = "settlements", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")
tristis_set_plot <- effect_plot(data = d.acrtri, model = m.acrtri, pred = settlements, 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 = "")
#abundance near settlements is 2010.46% higher than far
tristis_set <- as.data.table(tristis_set_plot$data$Acridotheres.tristis)
(max(tristis_set)/min(tristis_set) - 1) * 100## [1] 2010.462tristis_set_plot$data$Acridotheres.tristis## [1] 0.0671933 1.41808891.4180889 / 0.0671933 -1## [1] 20.10462#or 21.1 fold
1.4180889 / 0.0671933## [1] 21.10462tristis_subunit <- effect_plot(data = d.acrtri, model = m.acrtri, pred = subunit, 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 = "")
#emmeans
EMM <- emmeans(object = m.acrtri, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -1.18 0.216 Inf -1.60 -0.753
## Carmel -1.24 0.220 Inf -1.67 -0.808
## Galilee -4.02 0.430 Inf -4.87 -3.180
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel - Galilee
(exp(-1.24)/exp(-4.02) - 1) * 100## [1] 1511.902#Carmel - Judea
(exp(-1.18)/exp(-1.24) - 1) * 100## [1] 6.183655#Galilee - Judea
(exp(-1.18)/exp(-4.02) - 1) * 100## [1] 1611.577test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 0.0636 0.255 Inf 0.249 0.8032
## Judean Highlands - Galilee 2.8484 0.435 Inf 6.547 <.0001
## Carmel - Galilee 2.7848 0.433 Inf 6.436 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Carduelis carduelis - חוחית####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Carduelis carduelis", ylab_val = "Carduelis carduelis")
d.carduelis <- cbind(env_data,data.table(Carduelis.carduelis = spp_no_rare[,"Carduelis.carduelis"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.carduelis <- glm.nb(formula = Carduelis.carduelis ~ settlements + year_ct + subunit + cos_td_rad + sin_td_rad, data = d.carduelis)
coef(mva_m1)[,"Carduelis.carduelis"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -6.3665442 1.2522516 2.4837754 0.9127661 -0.1859762
## cos_td_rad sin_td_rad
## 2.5842809 -2.6428473coef(m.carduelis)## (Intercept) settlementsNear year_ct subunitCarmel subunitGalilee
## -6.3665485 1.2522511 -0.1859760 2.4837743 0.9127662
## cos_td_rad sin_td_rad
## 2.5842828 -2.6428525carduelis_set <- effect_plot(data = d.carduelis, model = m.carduelis, pred = settlements, 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 = "")
carduelis_set
#abundance near settlements is 249.82% higher than far
(max(carduelis_set$data$Carduelis.carduelis)/min(carduelis_set$data$Carduelis.carduelis)-1) * 100## [1] 249.82090.07259413 / 0.02075180 -1## [1] 2.498209carduelis_subunit <- effect_plot(data = d.carduelis, model = m.carduelis, pred = subunit, 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 = "")
carduelis_subunit
#emmeans
EMM <- emmeans(object = m.carduelis, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -3.249 0.457 Inf -4.15 -2.352
## Carmel -0.765 0.259 Inf -1.27 -0.258
## Galilee -2.336 0.354 Inf -3.03 -1.642
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel - Galilee
(exp(-0.765)/exp(-2.336) - 1) * 100## [1] 381.1457#Carmel - Judea
(exp(-0.765)/exp(-3.249) - 1) * 100## [1] 1098.913#Galilee - Judea
(exp(-2.336)/exp(-3.249) - 1) * 100## [1] 149.1787test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -2.484 0.516 Inf -4.810 <.0001
## Judean Highlands - Galilee -0.913 0.561 Inf -1.627 0.1038
## Carmel - Galilee 1.571 0.432 Inf 3.635 0.0004
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests######Cecropis daurica - סנונית מערות######
d.cecropis <- cbind(env_data,data.table(Cecropis.daurica = spp_no_rare[,"Cecropis.daurica"]))
m.cecropis <- glm.nb(formula = Cecropis.daurica ~ settlements + subunit + year_ct + sin_td_rad + cos_td_rad, data = d.cecropis)
coef(mva_m1)[,"Cecropis.daurica"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 0.98095288 2.16869649 -0.45368366 -0.08874088 0.10836920
## cos_td_rad sin_td_rad
## -2.15241967 2.63508635coef(m.cecropis)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 0.98094089 2.16869257 -0.45368117 -0.08874169 0.10836879
## sin_td_rad cos_td_rad
## 2.63507576 -2.15240660cecropis_set <- effect_plot(model = m.cecropis, data=d.cecropis, pred = settlements, 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 on abundance of Passer domesticus")
#Near settlements is 774% highern than far
(max(cecropis_set$data$Cecropis.daurica)/min(cecropis_set$data$Cecropis.daurica)-1)*100## [1] 774.6841cecropis_set$data$Cecropis.daurica## [1] 0.2076556 1.81633091.8163309/0.2076556 - 1## [1] 7.746843##########Chloris chloris - ירקון#######
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Chloris chloris", ylab_val = "Chloris chloris")
d.chloris <- cbind(env_data,data.table(Chloris.chloris = spp_no_rare[,"Chloris.chloris"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.chloris <- glm.nb(formula = Chloris.chloris ~ settlements + subunit + year_ct + sin_td_rad + cos_td_rad, data = d.chloris)
coef(mva_m1)[,"Chloris.chloris"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.17270658 1.46825586 -0.42961211 -0.09156341 0.11302623
## cos_td_rad sin_td_rad
## -0.19368715 0.42817036coef(m.chloris)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.17266422 1.46823977 -0.42962756 -0.09155752 0.11302354
## sin_td_rad cos_td_rad
## 0.42819754 -0.19369059chloris_set <- effect_plot(model = m.chloris, data=d.chloris, pred = settlements, 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 on abundance of Passer domesticus")
#Near settlements is 774% highern than far
(max(chloris_set$data$Chloris.chloris)/min(chloris_set$data$Chloris.chloris)-1)*100## [1] 334.1586chloris_set$data$Chloris.chloris## [1] 0.3610366 1.56747141.5674714/0.3610366 - 1## [1] 3.341586########Cinnyris osea - צופית בוהקת#######
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Cinnyris osea", ylab_val = "Cinnyris osea")
d.cinnyris <- cbind(env_data,data.table(Cinnyris.osea = spp_no_rare[,"Cinnyris.osea"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.cinnyris <- glm.nb(formula = Cinnyris.osea~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.cinnyris)
coef(mva_m1)[,"Cinnyris.osea"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.51622489 0.87058980 -0.26514518 -0.71725657 -0.07193677
## cos_td_rad sin_td_rad
## 0.27528358 -0.82378562coef(m.cinnyris)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.51622477 0.87058991 -0.26514518 -0.71725663 -0.07193676
## cos_td_rad sin_td_rad
## 0.27528327 -0.82378562# plot_model_effect(P.anal = d.cinnyris, m = m.cinnyris, eff2plot = "settlements", export_plot = TRUE, plot_points = TRUE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")
cinnyris_set <- effect_plot(data = d.cinnyris, model = m.cinnyris, pred = settlements, 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 = "")
#abundance near settlements is 138.83% higher than far
(max(cinnyris_set$data$Cinnyris.osea)/min(cinnyris_set$data$Cinnyris.osea)-1)*100## [1] 138.8319cinnyris_set$data$Cinnyris.osea## [1] 0.832746 1.9888631.988863 / 0.832746 -1## [1] 1.388319cinnyris_subunit <- effect_plot(data = d.cinnyris, model = m.cinnyris, pred = subunit, 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 = "")
cinnyris_subunit
#emmeans
EMM <- emmeans(object = m.cinnyris, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 0.2523 0.103 Inf 0.0507 0.454
## Carmel -0.0129 0.110 Inf -0.2277 0.202
## Galilee -0.4650 0.128 Inf -0.7152 -0.215
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel - Galilee
(exp(-0.0129)/exp(-0.4650) - 1) * 100## [1] 57.16091#Judea-Carmel
(exp(0.2523)/exp(-0.0129) - 1) * 100## [1] 30.36917#Judea-Galilee
(exp(0.2523)/exp(-0.4650) - 1) * 100## [1] 104.8894test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 0.265 0.149 Inf 1.785 0.0743
## Judean Highlands - Galilee 0.717 0.162 Inf 4.420 <.0001
## Carmel - Galilee 0.452 0.167 Inf 2.715 0.0100
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests#####Columba livia - יונת סלעים#########
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Columba livia", ylab_val = "Columba livia")
d.columba <- cbind(env_data,data.table(Columba.livia = spp_no_rare[,"Columba.livia"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.columba <- glm.nb(formula = Columba.livia ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad
, data = d.columba)
coef(mva_m1)[,"Columba.livia"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -8.7541198 1.8964697 0.5402957 -1.5045681 0.1794411
## cos_td_rad sin_td_rad
## 5.1582629 -5.7721640coef(m.columba)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -8.7540662 1.8964544 0.5403010 -1.5045420 0.1794368
## cos_td_rad sin_td_rad
## 5.1582320 -5.7721438#
# plot_model_effect(P.anal = d.columba, m = m.columba, eff2plot = "year_ct", export_plot = TRUE, plot_points = FALSE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")
columba_time <- effect_plot(data = d.columba, model = m.columba, pred = year_ct, 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 trend - increase of 402.75% total
(max(columba_time$data$Columba.livia)/min(columba_time$data$Columba.livia)-1)*100## [1] 402.7542columba_time$data$Columba.livia## [1] 0.1812109 0.1841911 0.1872204 0.1902995 0.1934292 0.1966103 0.1998438
## [8] 0.2031305 0.2064712 0.2098669 0.2133184 0.2168267 0.2203927 0.2240173
## [15] 0.2277015 0.2314464 0.2352528 0.2391218 0.2430544 0.2470518 0.2511148
## [22] 0.2552447 0.2594425 0.2637094 0.2680464 0.2724548 0.2769356 0.2814901
## [29] 0.2861196 0.2908252 0.2956081 0.3004698 0.3054114 0.3104342 0.3155397
## [36] 0.3207291 0.3260039 0.3313654 0.3368151 0.3423545 0.3479849 0.3537079
## [43] 0.3595251 0.3654379 0.3714480 0.3775569 0.3837663 0.3900778 0.3964931
## [50] 0.4030139 0.4096419 0.4163790 0.4232269 0.4301873 0.4372623 0.4444536
## [57] 0.4517632 0.4591930 0.4667450 0.4744211 0.4822236 0.4901543 0.4982155
## [64] 0.5064093 0.5147378 0.5232033 0.5318080 0.5405542 0.5494443 0.5584805
## [71] 0.5676654 0.5770014 0.5864909 0.5961364 0.6059406 0.6159060 0.6260353
## [78] 0.6363313 0.6467965 0.6574338 0.6682461 0.6792363 0.6904071 0.7017617
## [85] 0.7133030 0.7250341 0.7369582 0.7490784 0.7613979 0.7739200 0.7866480
## [92] 0.7995854 0.8127356 0.8261020 0.8396882 0.8534979 0.8675347 0.8818024
## [99] 0.8963047 0.91104550.9110455/0.1812109 - 1## [1] 4.027542columba_subunit <- effect_plot(data = d.columba, model = m.columba, pred = subunit, 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 = "")
columba_subunit
#emmeans
EMM <- emmeans(object = m.columba, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 0.101 0.254 Inf -0.397 0.600
## Carmel 0.642 0.247 Inf 0.157 1.127
## Galilee -1.403 0.307 Inf -2.005 -0.801
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel - Galilee
(exp(0.642)/exp(-1.403) - 1) * 100## [1] 672.9159#Carmel-Judea
(exp(0.642)/exp(0.101) - 1) * 100## [1] 71.77237#Judea-Galilee
(exp(0.101)/exp(-1.403) - 1) * 100## [1] 349.9652test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -0.54 0.354 Inf -1.528 0.1265
## Judean Highlands - Galilee 1.50 0.394 Inf 3.815 0.0002
## Carmel - Galilee 2.04 0.390 Inf 5.243 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests########Corvus cornix - עורב אפור#################
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Corvus cornix", ylab_val = "Corvus cornix")
d.corvus <- cbind(env_data,data.table(Corvus.cornix = spp_no_rare[,"Corvus.cornix"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.corvus <- glm.nb(formula = Corvus.cornix ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.corvus)
coef(mva_m1)[,"Corvus.cornix"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -2.033847764 1.253931156 1.716084872 0.856113981 0.003008888
## cos_td_rad sin_td_rad
## 0.134731455 0.386707210coef(m.corvus)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -2.033846968 1.253930552 1.716084018 0.856113449 0.003008826
## cos_td_rad sin_td_rad
## 0.134731422 0.386706516corvus_set <- effect_plot(data = d.corvus, model = m.corvus, pred = settlements, 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 = "")
#abundance near settlements is 250.04% higher than far
(max(corvus_set$data$Corvus.cornix)/min(corvus_set$data$Corvus.cornix)-1)*100## [1] 250.4089corvus_set$data$Corvus.cornix## [1] 0.1164205 0.40794780.4079478/0.1164205 -1## [1] 2.504089#or 3.5 fold
0.4079478 / 0.1164205## [1] 3.504089corvus_subunit <- effect_plot(data = d.corvus, model = m.corvus, pred = subunit, 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 = "")
corvus_subunit
#emmeans
EMM <- emmeans(object = m.corvus, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -1.524 0.213 Inf -1.9408 -1.106
## Carmel 0.193 0.146 Inf -0.0927 0.478
## Galilee -0.667 0.173 Inf -1.0065 -0.328
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel - Galilee
(exp(0.193)/exp(-0.667) - 1) * 100## [1] 136.3161#Carmel-Judea
(exp(0.193)/exp(-1.524) - 1) * 100## [1] 456.78#Galilee-Judea
(exp(-0.667)/exp(-1.524) - 1) * 100## [1] 135.6082test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -1.716 0.255 Inf -6.724 <.0001
## Judean Highlands - Galilee -0.856 0.270 Inf -3.168 0.0015
## Carmel - Galilee 0.860 0.224 Inf 3.838 0.0002
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests#########Curruca curruca - סבכי טוחנים########
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Curruca curruca", ylab_val = "Curruca curruca")
d.curruca <- cbind(env_data,data.table(Curruca.curruca = spp_no_rare[,"Curruca.curruca"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.curruca <- glm.nb(formula = Curruca.curruca ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.curruca)
coef(mva_m1)[,"Curruca.curruca"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -19.51596678 -0.76748669 11.66738622 12.74590976 0.06434301
## cos_td_rad sin_td_rad
## 1.31234280 -6.42805839coef(m.curruca)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -41.73709897 -0.76749621 33.88839482 34.96693121 0.06434496
## cos_td_rad sin_td_rad
## 1.31237217 -6.42817587#
# plot_model_effect(P.anal = d.curruca, m = m.curruca, eff2plot = "year_ct", export_plot = TRUE, plot_points = FALSE, fontname = "Almoni ML v5 AAA", outpath = "../output/maquis_new/")
curruca_time <- effect_plot(data = d.curruca, model = m.curruca, pred = year_ct, 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 = "")## Error in qr.solve(qr.R(qr.lm(object))[p1, p1]): singular matrix 'a' in solvecurruca_time## Error in eval(expr, envir, enclos): object 'curruca_time' not found#Temporal trend - increase of 402.75% total
(max(curruca_time$data$Curruca.curruca)/min(curruca_time$data$Curruca.curruca)-1)*100## Error in eval(expr, envir, enclos): object 'curruca_time' not foundcurruca_time$data$Curruca.curruca## Error in eval(expr, envir, enclos): object 'curruca_time' not found0.9110455/0.1812109 - 1## [1] 4.027542# or 5.02 fold increase
0.9110455/0.1812109## [1] 5.027542curruca_subunit <- effect_plot(data = d.curruca, model = m.curruca, pred = subunit, 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 = "")## Error in qr.solve(qr.R(qr.lm(object))[p1, p1]): singular matrix 'a' in solvecurruca_subunit## Error in eval(expr, envir, enclos): object 'curruca_subunit' not found#emmeans
EMM <- emmeans(object = m.curruca, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -37.11 4717424 Inf -9246019 9245945
## Carmel -3.22 0 Inf -4 -2
## Galilee -2.14 0 Inf -3 -1
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Galilee-Carmel
(exp(-2.14)/exp(-3.22) - 1) * 100## [1] 194.468#Carmel-Judea
(exp(-3.22)/exp(-37.13) - 1) * 100## [1] 5.332439e+16#Galilee-Judea
(exp(-2.14)/exp(-37.13) - 1) * 100## [1] 1.570232e+17test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -33.89 4717424 Inf 0.000 1.0000
## Judean Highlands - Galilee -34.97 4717424 Inf 0.000 1.0000
## Carmel - Galilee -1.08 0 Inf -2.306 0.0633
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests########Curruca melanocephala - סבכי שחור-ראש#######
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Curruca melanocephala", ylab_val = "Curruca melanocephala")
d.melano <- cbind(env_data,data.table(Curruca.melanocephala = spp_no_rare[,"Curruca.melanocephala"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.melano <- glm.nb(formula = Curruca.melanocephala ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.melano)
coef(mva_m1)[,"Curruca.melanocephala"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 1.62651788 -0.42701053 -0.06039957 -0.66627444 -0.01213706
## cos_td_rad sin_td_rad
## 0.29890567 -0.26265011coef(m.melano)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 1.62651794 -0.42701053 -0.06039957 -0.66627442 -0.01213706
## cos_td_rad sin_td_rad
## 0.29890565 -0.26265002melano_set <- effect_plot(data = d.melano, model = m.melano, pred = settlements, 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 = "")
#abundance far from settlements is 53.26% higher than near
diff(melano_set$data$Curruca.melanocephala)/min(melano_set$data$Curruca.melanocephala)## [1] -0.5326688(6.920159/4.515104 -1) * 100## [1] 53.26688#or 1.53 fold
6.920159 / 4.515104## [1] 1.532669melano_subunit <- effect_plot(data = d.melano, model = m.melano, pred = subunit, 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 = "")
melano_subunit
#emmeans
EMM <- emmeans(object = m.melano, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 1.72 0.0713 Inf 1.581 1.86
## Carmel 1.66 0.0716 Inf 1.520 1.80
## Galilee 1.05 0.0799 Inf 0.898 1.21
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(1.66)/exp(1.05) - 1) * 100## [1] 84.04314#Judea-Carmel
(exp(1.72)/exp(1.66) - 1) * 100## [1] 6.183655#Judea-Galilee
(exp(1.72)/exp(1.05) - 1) * 100## [1] 95.42373test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 0.0604 0.101 Inf 0.598 0.5499
## Judean Highlands - Galilee 0.6663 0.107 Inf 6.229 <.0001
## Carmel - Galilee 0.6059 0.107 Inf 5.653 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests#######Dendrocopos syriacus - נקר סורי############
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Dendrocopos syriacus", ylab_val = "Dendrocopos syriacus")
d.syriacus <- cbind(env_data,data.table(Dendrocopos.syriacus = spp_no_rare[,"Dendrocopos.syriacus"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.syriacus <- glm.nb(formula = Dendrocopos.syriacus ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.syriacus)
coef(mva_m1)[,"Dendrocopos.syriacus"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -6.572022180 0.934448120 -0.650716218 -0.238886114 -0.004452597
## cos_td_rad sin_td_rad
## 3.295442948 -2.857773539coef(m.syriacus)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -6.572021986 0.934448131 -0.650716326 -0.238886302 -0.004452595
## cos_td_rad sin_td_rad
## 3.295442934 -2.857773352syriacus_set <- effect_plot(data = d.syriacus, model = m.syriacus, pred = settlements, 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 = "")
#abundance near settlements is 154.58% higher than far
(max(syriacus_set$data$Dendrocopos.syriacus)/min(syriacus_set$data$Dendrocopos.syriacus)-1)*100## [1] 154.5808syriacus_set$data$Dendrocopos.syriacus## [1] 0.07588191 0.193180780.19318078/0.07588191 -1## [1] 1.545808#or 2.54 fold
0.19318078/0.07588191## [1] 2.545808#####Passer domesticus - דרור הבית#####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Passer domesticus", ylab_val = "Passer domesticus")
d.pasdom <- cbind(env_data,data.table(Passer.domesticus = spp_no_rare[,"Passer.domesticus"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.pasdom <- glm.nb(formula = Passer.domesticus ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.pasdom)
coef(mva_m1)[,"Passer.domesticus"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -5.2174230 4.9808727 -1.2576557 0.3878957 -0.1176087
## cos_td_rad sin_td_rad
## 1.4382456 -2.3168731coef(m.pasdom)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -5.2174237 4.9808727 -1.2576557 0.3878958 -0.1176087
## cos_td_rad sin_td_rad
## 1.4382462 -2.3168737passer_set <- effect_plot(data = d.pasdom, model = m.pasdom, pred = settlements, 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 = "")
#abundance near settlements is 14,460.14% higher than far
diff(passer_set$data$Passer.domesticus)/min(passer_set$data$Passer.domesticus) * 100## [1] 14460.144.84735843 / 0.03329198 -1## [1] 144.6014#or 156.37 fold
4.84624809 / 0.03099124## [1] 156.3748passer_subunit <- effect_plot(data = d.pasdom, model = m.pasdom, pred = subunit, 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 = "")
melano_subunit
#emmeans
EMM <- emmeans(object = m.pasdom, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -0.912 0.270 Inf -1.44 -0.3838
## Carmel -2.170 0.304 Inf -2.76 -1.5744
## Galilee -0.524 0.260 Inf -1.03 -0.0151
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Galilee-Carmel
(exp(-0.524)/exp(-2.170) - 1) * 100## [1] 418.6194#Judea-Carmel
(exp(-0.912)/exp(-2.170) - 1) * 100## [1] 251.8378#Galilee-Judea
(exp(-0.524)/exp(-0.912) - 1) * 100## [1] 47.40298test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 1.258 0.313 Inf 4.022 0.0001
## Judean Highlands - Galilee -0.388 0.296 Inf -1.311 0.1899
## Carmel - Galilee -1.646 0.315 Inf -5.230 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Prinia gracilis - פשוש####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Prinia gracilis", ylab_val = "Prinia gracilis")
d.prigra <- cbind(env_data,data.table(Prinia.gracilis = spp_no_rare[,"Prinia.gracilis"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.prigra <- glm.nb(formula = Prinia.gracilis ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.prigra)
coef(mva_m1)[,"Prinia.gracilis"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.12702277 1.51717609 0.26228731 -0.06852479 -0.08504256
## cos_td_rad sin_td_rad
## -0.36574149 0.92278481coef(m.prigra)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.12702190 1.51717688 0.26228733 -0.06852537 -0.08504253
## cos_td_rad sin_td_rad
## -0.36574233 0.92278602prigra_set <- effect_plot(data = d.prigra, model = m.prigra, pred = settlements, 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 = "")
#abundance near settlements is 355.93% higher than far
diff(prigra_set$data$Prinia.gracilis)/min(prigra_set$data$Prinia.gracilis) * 100## [1] 355.9335prigra_set$data$Prinia.gracilis## [1] 0.262919 1.1987361.198736/0.262919 -1## [1] 3.559336prigra_time <- effect_plot(data = d.prigra, model = m.prigra, pred = year_ct, 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 trend - decrease of 53.48%
(max(prigra_time$data$Prinia.gracilis)-min(prigra_time$data$Prinia.gracilis))/max(prigra_time$data$Prinia.gracilis)## [1] 0.5348442prigra_time$data$Prinia.gracilis## [1] 0.3954450 0.3923996 0.3893776 0.3863788 0.3834032 0.3804505 0.3775205
## [8] 0.3746131 0.3717281 0.3688653 0.3660245 0.3632057 0.3604085 0.3576329
## [15] 0.3548786 0.3521456 0.3494336 0.3467425 0.3440721 0.3414223 0.3387929
## [22] 0.3361837 0.3335947 0.3310255 0.3284762 0.3259465 0.3234363 0.3209454
## [29] 0.3184737 0.3160210 0.3135872 0.3111722 0.3087758 0.3063978 0.3040381
## [36] 0.3016966 0.2993731 0.2970676 0.2947797 0.2925096 0.2902568 0.2880215
## [43] 0.2858033 0.2836023 0.2814181 0.2792509 0.2771002 0.2749662 0.2728486
## [50] 0.2707473 0.2686622 0.2665931 0.2645400 0.2625027 0.2604811 0.2584750
## [57] 0.2564844 0.2545092 0.2525491 0.2506041 0.2486741 0.2467590 0.2448587
## [64] 0.2429729 0.2411017 0.2392449 0.2374024 0.2355741 0.2337598 0.2319596
## [71] 0.2301732 0.2284005 0.2266416 0.2248961 0.2231641 0.2214455 0.2197400
## [78] 0.2180477 0.2163685 0.2147022 0.2130487 0.2114079 0.2097798 0.2081642
## [85] 0.2065611 0.2049703 0.2033917 0.2018253 0.2002710 0.1987287 0.1971982
## [92] 0.1956795 0.1941725 0.1926771 0.1911932 0.1897208 0.1882597 0.1868098
## [99] 0.1853712 0.18394361 - 0.1839436/0.3954450 ## [1] 0.534844####Psittacula krameri - דררה#####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Psittacula krameri", ylab_val = "Psittacula krameri")
d.psitta <- cbind(env_data,data.table(Psittacula.krameri = spp_no_rare[,"Psittacula.krameri"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.psitta <- glm.nb(formula = Psittacula.krameri ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.psitta)
coef(mva_m1)[,"Psittacula.krameri"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -5.269088823 1.668902826 1.190144622 -0.452394245 -0.009682361
## cos_td_rad sin_td_rad
## 0.995954283 -2.204367675coef(m.psitta)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -5.269099503 1.668904042 1.190146788 -0.452393130 -0.009682248
## cos_td_rad sin_td_rad
## 0.995960536 -2.204373874psitta_set_plot <- effect_plot(data = d.psitta, model = m.psitta, pred = settlements, 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 = "")
#abundance near settlements is 430.63% higher than far
pristta_set <- as.data.table(psitta_set_plot$data$Psittacula.krameri)
(max(pristta_set)/min(pristta_set) - 1) * 100## [1] 430.6349psitta_set_plot$data$Psittacula.krameri## [1] 0.03629751 0.192607260.19260726/0.03629751 -1## [1] 4.306349#or 5.3 fold
0.19260726/0.03629751## [1] 5.306349psitta_subunit <- effect_plot(data = d.psitta, model = m.psitta, pred = subunit, 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 = "")
psitta_subunit
#emmeans
EMM <- emmeans(object = m.psitta, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -2.48 0.347 Inf -3.16 -1.801
## Carmel -1.29 0.263 Inf -1.81 -0.777
## Galilee -2.93 0.404 Inf -3.72 -2.143
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(-1.29)/exp(-2.93) - 1) * 100## [1] 415.517#Carmel-Judea
(exp(-1.29)/exp(-2.48) - 1) * 100## [1] 228.7081#Judea-Galilee
(exp(-2.48)/exp(-2.93) - 1) * 100## [1] 56.83122test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -1.190 0.419 Inf -2.839 0.0068
## Judean Highlands - Galilee 0.452 0.510 Inf 0.887 0.3750
## Carmel - Galilee 1.643 0.465 Inf 3.533 0.0012
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests#####Pycnonotus.xanthopygos - בולבול ממושקף#####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Pycnonotus xanthopygos", ylab_val = "Pycnonotus xanthopygos")
d.pycno <- cbind(env_data,data.table(Pycnonotus.xanthopygos = spp_no_rare[,"Pycnonotus.xanthopygos"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.pycno <- glm.nb(formula = Pycnonotus.xanthopygos ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.pycno)
coef(mva_m1)[,"Pycnonotus.xanthopygos"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.887826774 1.006461173 -0.586201462 -0.823030493 -0.002510203
## cos_td_rad sin_td_rad
## 1.409560683 -0.456734212coef(m.pycno)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -0.887826608 1.006461068 -0.586201536 -0.823030497 -0.002510177
## cos_td_rad sin_td_rad
## 1.409560652 -0.456733882pycno_set <- effect_plot(data = d.pycno, model = m.pycno, pred = settlements, 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 = "")
#abundance near settlements is 173.59% higher than far
(max(pycno_set$data$Pycnonotus.xanthopygos)/min(pycno_set$data$Pycnonotus.xanthopygos)-1)*100## [1] 173.5902pycno_set$data$Pycnonotus.xanthopygos## [1] 1.445443 3.9545893.954589/1.445443 -1## [1] 1.735901#or 1.73 fold
3.954589 / 1.445443## [1] 2.735901pycno_subunit <- effect_plot(data = d.pycno, model = m.pycno, pred = subunit, 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 = "")
pycno_subunit
#emmeans
EMM <- emmeans(object = m.pycno, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 0.8716 0.111 Inf 0.6546 1.089
## Carmel 0.2854 0.120 Inf 0.0497 0.521
## Galilee 0.0486 0.128 Inf -0.2018 0.299
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(0.2854 )/exp(0.0486) - 1) * 100## [1] 26.71877#Judea-Carmel
(exp(0.8716)/exp(0.2854) - 1) * 100## [1] 79.71463#Judea-Galilee
(exp(0.8716)/exp(0.0486) - 1) * 100## [1] 127.7322test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 0.586 0.163 Inf 3.605 0.0005
## Judean Highlands - Galilee 0.823 0.168 Inf 4.899 <.0001
## Carmel - Galilee 0.237 0.174 Inf 1.357 0.1746
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Spilopelia senegalensis - צוצלת####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Spilopelia senegalensis", ylab_val = "Spilopelia senegalensis")
d.spilop <- cbind(env_data,data.table(Spilopelia.senegalensis = spp_no_rare[,"Spilopelia.senegalensis"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.spilop <- glm.nb(formula = Spilopelia.senegalensis ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.spilop)
coef(mva_m1)[,"Spilopelia.senegalensis"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -2.48717563 2.85687260 -0.13160712 -1.14301540 -0.04214909
## cos_td_rad sin_td_rad
## 0.63738395 -0.60406669coef(m.spilop)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -2.48716893 2.85687304 -0.13160658 -1.14301634 -0.04214947
## cos_td_rad sin_td_rad
## 0.63737835 -0.60406475spilo_set <- effect_plot(data = d.spilop, model = m.spilop, pred = settlements, 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 = "")
#abundance near settlements is 1640.7% higher than far
(max(spilo_set$data$Spilopelia.senegalensis)/min(spilo_set$data$Spilopelia.senegalensis)-1)*100## [1] 1640.701spilo_set$data$Spilopelia.senegalensis## [1] 0.152156 2.6485802.64858/0.152156 - 1## [1] 16.407#or 17.4 fold
2.64858/0.152156## [1] 17.407spilop_subunit <- effect_plot(data = d.spilop, model = m.spilop, pred = subunit, 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 = "")
spilop_subunit
#emmeans
EMM <- emmeans(object = m.spilop, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -0.454 0.154 Inf -0.756 -0.153
## Carmel -0.586 0.157 Inf -0.895 -0.277
## Galilee -1.597 0.198 Inf -1.985 -1.210
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(-0.586)/exp(-1.597) - 1) * 100## [1] 174.8348#Judea-Carmel
(exp(-0.454)/exp(-0.586) - 1) * 100## [1] 14.11083#Judea-Galilee
(exp(-0.454)/exp(-1.597) - 1) * 100## [1] 213.6163test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 0.132 0.190 Inf 0.692 0.4887
## Judean Highlands - Galilee 1.143 0.220 Inf 5.195 <.0001
## Carmel - Galilee 1.011 0.222 Inf 4.560 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Streptopelia decaocto - תור צווארון####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Streptopelia decaocto", ylab_val = "Streptopelia decaocto")
d.decaocto <- cbind(env_data,data.table(Streptopelia.decaocto=spp_no_rare[,"Streptopelia.decaocto"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.decaocto <- glm.nb(formula = Streptopelia.decaocto ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.decaocto)
coef(mva_m1)[,"Streptopelia.decaocto"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.436931728 0.157512825 0.290118464 -0.165355750 -0.005583385
## cos_td_rad sin_td_rad
## 2.033709987 -1.869831164coef(m.decaocto)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.43693148 0.15751286 0.29011842 -0.16535581 -0.00558339
## cos_td_rad sin_td_rad
## 2.03370977 -1.86983108decaocto_subunit <- effect_plot(data = d.decaocto, model = m.decaocto, pred = subunit, 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 = "")
decaocto_subunit
#emmeans
EMM <- emmeans(object = m.decaocto, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 1.15 0.0887 Inf 0.981 1.33
## Carmel 1.45 0.0852 Inf 1.278 1.61
## Galilee 0.99 0.0918 Inf 0.810 1.17
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(1.45)/exp(0.99) - 1) * 100## [1] 58.4074#Carmel-Judea
(exp(1.45)/exp(1.15) - 1) * 100## [1] 34.98588#Judea-Galilee
(exp(1.15)/exp(0.99) - 1) * 100## [1] 17.35109test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -0.290 0.123 Inf -2.359 0.0275
## Judean Highlands - Galilee 0.165 0.127 Inf 1.297 0.1945
## Carmel - Galilee 0.455 0.125 Inf 3.637 0.0008
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Streptopelia turtur - תור מצוי####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Streptopelia turtur", ylab_val = "Streptopelia turtur")
d.strtur <- cbind(env_data,data.table(Streptopelia.turtur = spp_no_rare[,"Streptopelia.turtur"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.strtur <- glm.nb(formula = Streptopelia.turtur ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.strtur)
coef(mva_m1)[,"Streptopelia.turtur"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -4.0568827 -0.3138816 0.8126151 0.2758993 -0.1499697
## cos_td_rad sin_td_rad
## 2.8254657 -1.6813797coef(m.strtur)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -4.0568824 -0.3138816 0.8126152 0.2758993 -0.1499697
## cos_td_rad sin_td_rad
## 2.8254655 -1.6813795turtur_time <- effect_plot(data = d.strtur, model = m.strtur, pred = year_ct, 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 trend - decrese of 285.63% total
(max(turtur_time$data$Streptopelia.turtur)-min(turtur_time$data$Streptopelia.turtur))/max(turtur_time$data$Streptopelia.turtur) * 100## [1] 74.0689turtur_time$data$Streptopelia.turtur## [1] 0.34455932 0.33989362 0.33529109 0.33075088 0.32627215 0.32185407
## [7] 0.31749581 0.31319658 0.30895555 0.30477196 0.30064501 0.29657395
## [13] 0.29255802 0.28859646 0.28468855 0.28083356 0.27703077 0.27327947
## [19] 0.26957896 0.26592857 0.26232761 0.25877541 0.25527131 0.25181465
## [25] 0.24840481 0.24504114 0.24172301 0.23844982 0.23522095 0.23203580
## [31] 0.22889378 0.22579431 0.22273681 0.21972071 0.21674545 0.21381048
## [37] 0.21091526 0.20805923 0.20524188 0.20246269 0.19972112 0.19701668
## [43] 0.19434886 0.19171716 0.18912110 0.18656020 0.18403397 0.18154195
## [49] 0.17908367 0.17665868 0.17426653 0.17190677 0.16957897 0.16728268
## [55] 0.16501749 0.16278298 0.16057872 0.15840431 0.15625934 0.15414342
## [61] 0.15205615 0.14999714 0.14796602 0.14596240 0.14398591 0.14203618
## [67] 0.14011286 0.13821557 0.13634398 0.13449774 0.13267649 0.13087991
## [73] 0.12910765 0.12735939 0.12563481 0.12393358 0.12225538 0.12059991
## [79] 0.11896685 0.11735591 0.11576679 0.11419918 0.11265280 0.11112736
## [85] 0.10962257 0.10813816 0.10667385 0.10522937 0.10380445 0.10239883
## [91] 0.10101223 0.09964442 0.09829512 0.09696410 0.09565110 0.09435588
## [97] 0.09307820 0.09181782 0.09057451 0.089348031 - 0.08934803/0.34455932## [1] 0.740689# or 3.85 fold decrese
0.34455932/0.08934803## [1] 3.856373####Troglodytes troglodytes####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Troglodytes troglodytes", ylab_val = "Troglodytes troglodytes")
d.troglo <- cbind(env_data,data.table(Troglodytes.troglodytes = spp_no_rare[,"Troglodytes.troglodytes"]))
# this is not the same model as the manyglm, to show that the trend is mainly near settlements.
m.troglo <- glm.nb(formula = Troglodytes.troglodytes ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.troglo)
coef(mva_m1)[,"Troglodytes.troglodytes"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -3.84236334 -1.19509372 -0.04556621 0.45798127 -0.10261925
## cos_td_rad sin_td_rad
## 2.44609621 -3.55529291coef(m.troglo)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -3.84238142 -1.19509312 -0.04556577 0.45798207 -0.10261927
## cos_td_rad sin_td_rad
## 2.44610919 -3.55530696troglo_set <- effect_plot(data = d.troglo, model = m.troglo, pred = settlements, 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 = "")
troglo_set$data$Troglodytes.troglodytes## [1] 0.5981707 0.1810518#abundance far from settlements is 230.38% higher than near
(max(troglo_set$data$Troglodytes.troglodytes)/min(troglo_set$data$Troglodytes.troglodytes)-1)*100## [1] 230.3865troglo_set$data$Troglodytes.troglodytes## [1] 0.5981707 0.18105180.5981707/0.1810518 - 1## [1] 2.303865#or 3.3 fold
0.5981707/0.1810518## [1] 3.303865####Garrulus glandarius - עורבני####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Garrulus glandarius", ylab_val = "Garrulus glandarius")
d.glandarius <- cbind(env_data,data.table(Garrulus.glandarius=spp_no_rare[,"Garrulus.glandarius"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.glandarius <- glm.nb(formula = Garrulus.glandarius ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.glandarius)
coef(mva_m1)[,"Garrulus.glandarius"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.11756667 0.25786121 0.33158628 -0.69131910 -0.05805483
## cos_td_rad sin_td_rad
## 1.00690364 -0.79489729coef(m.decaocto)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.43693148 0.15751286 0.29011842 -0.16535581 -0.00558339
## cos_td_rad sin_td_rad
## 2.03370977 -1.86983108glandarius_subunit <- effect_plot(data = d.glandarius, model = m.glandarius, pred = subunit, 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 = "")
glandarius_subunit
#emmeans
EMM <- emmeans(object = m.glandarius, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -0.0865 0.135 Inf -0.35092 0.178
## Carmel 0.2451 0.126 Inf -0.00262 0.493
## Galilee -0.7778 0.161 Inf -1.09430 -0.461
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(0.2451)/exp(-0.7778) - 1) * 100## [1] 178.1249#Carmel-Judea
(exp(0.2451)/exp(-0.0865) - 1) * 100## [1] 39.31955#Judea-Galilee
(exp(-0.0865)/exp(-0.7778) - 1) * 100## [1] 99.6309test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -0.332 0.185 Inf -1.796 0.0725
## Judean Highlands - Galilee 0.691 0.210 Inf 3.294 0.0015
## Carmel - Galilee 1.023 0.205 Inf 4.993 <.0001
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Parus major - ירגזי מצוי####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Parus major", ylab_val = "Parus major")
d.parus <- cbind(env_data,data.table(Parus.major=spp_no_rare[,"Parus.major"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.parus <- glm.nb(formula = Parus.major ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.parus)
coef(mva_m1)[,"Parus.major"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.94131820 0.21506433 1.02233232 0.59434990 -0.06096473
## cos_td_rad sin_td_rad
## 0.90791148 -2.05889505coef(m.parus)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## -1.94131734 0.21506408 1.02233220 0.59435010 -0.06096475
## cos_td_rad sin_td_rad
## 0.90791116 -2.05889431parus_subunit <- effect_plot(data = d.parus, model = m.parus, pred = subunit, 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 = "")
parus_subunit
#emmeans
EMM <- emmeans(object = m.parus, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands -0.275 0.145 Inf -0.5580 0.00871
## Carmel 0.748 0.123 Inf 0.5071 0.98826
## Galilee 0.320 0.132 Inf 0.0618 0.57767
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(0.748)/exp(0.320) - 1) * 100## [1] 53.41861#Carmel-Judea
(exp(0.748)/exp(-0.275) - 1) * 100## [1] 178.1527#Galilee-Judea
(exp(0.320)/exp(-0.275) - 1) * 100## [1] 81.30309test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel -1.022 0.189 Inf -5.397 <.0001
## Judean Highlands - Galilee -0.594 0.195 Inf -3.043 0.0035
## Carmel - Galilee 0.428 0.180 Inf 2.380 0.0173
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests####Corvus.monedula - קאק####
plot_alpha_diversity(P = abu_by_spp.maquis, x_val = "settlements", y_val = "Corvus monedula", ylab_val = "Corvus monedula")
d.monedula <- cbind(env_data,data.table(Corvus.monedula=spp_no_rare[,"Corvus.monedula"]))
# this is not the same model as the manyglm, to show the temporal trend better.
m.monedula <- glm.nb(formula = Corvus.monedula ~ settlements + subunit + year_ct + cos_td_rad + sin_td_rad, data = d.monedula)
coef(mva_m1)[,"Corvus.monedula"]## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 0.82864990 -0.01170426 -2.01179194 -3.12892756 -0.02734277
## cos_td_rad sin_td_rad
## -0.01871584 0.59352267coef(m.monedula)## (Intercept) settlementsNear subunitCarmel subunitGalilee year_ct
## 0.82895652 -0.01174288 -2.01172173 -3.12898154 -0.02737127
## cos_td_rad sin_td_rad
## -0.01884941 0.59362987monedula_subunit <- effect_plot(data = d.monedula, model = m.monedula, pred = subunit, 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 = "")
monedula_subunit
#emmeans
EMM <- emmeans(object = m.monedula, ~ subunit)
print(EMM)## subunit emmean SE df asymp.LCL asymp.UCL
## Judean Highlands 0.33 0.377 Inf -0.409 1.070
## Carmel -1.68 0.416 Inf -2.497 -0.866
## Galilee -2.80 0.506 Inf -3.790 -1.808
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Carmel-Galilee
(exp(-1.68)/exp(-2.80) - 1) * 100## [1] 206.4854#Judea-Carmel
(exp(0.33 )/exp(-1.68) - 1) * 100## [1] 646.3317#Judea-Galilee
(exp(0.33 )/exp(-2.80) - 1) * 100## [1] 2187.398test_results_subunit <- test(pairs(EMM, by=NULL, adjust="fdr"))
print(test_results_subunit)## contrast estimate SE df z.ratio p.value
## Judean Highlands - Carmel 2.01 0.562 Inf 3.581 0.0005
## Judean Highlands - Galilee 3.13 0.630 Inf 4.963 <.0001
## Carmel - Galilee 1.12 0.654 Inf 1.707 0.0878
##
## Results are averaged over the levels of: settlements
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 tests# Y <- as.data.table(spp_no_rare)
# rept_abu_settle <- stackedsdm(y = Y, formula_X = ~year_ct*subunit + cos_td_rad + sin_td_rad, data = env_data)
# rept_abu_subunit <- stackedsdm(y = Y, formula_X = ~settlements + year_ct + cos_td_rad + sin_td_rad , data = env_data)
# rept_lv_settle <- cord(rept_abu_settle)
# rept_lv_subunit <- cord(rept_abu_subunit)
#
# plot(rept_lv_settle, biplot = TRUE, site.col = ifelse(as.numeric(env_data$settlements)==2, "green", "brown"), main = "settlements (Near in green, Far in brown)")
#
# mycolors <- c('#98FB98','#FFC1C1','#87CEFA')
# plot(rept_lv_subunit, biplot = TRUE, site.col = mycolors[as.numeric(env_data$subunit)], main = "Subunits (Judea green, Carmel pink, Galilee blue)")