Birds_Maquis_analysis_by_plot_with_subunit_interactions
Ron Chen
2023-12-28
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.
## - Species less likely to be interacting with sampling location were EXCLUDEDvar_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 * year_ct + settlements * 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.882425207773876"Overdispersion parameter is < 1. Choose Poisson.
m0 <- mdl_r.poiss.int
me0 <- glmer(data = P.anal, formula = richness ~ settlements * year_ct + subunit * year_ct + settlements * subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0171254 (tol = 0.002, component 1)vif(me0)## GVIF Df GVIF^(1/(2*Df))
## settlements 5.291552 1 2.300337
## year_ct 4.773684 1 2.184876
## subunit 5.531759 2 1.533613
## cos_td_rad 5.682659 1 2.383833
## sin_td_rad 7.443852 1 2.728342
## settlements:year_ct 4.716119 1 2.171663
## year_ct:subunit 7.312211 2 1.644418
## settlements:subunit 5.020712 2 1.496895go 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 * year_ct + settlements *
## subunit + cos_td_rad + sin_td_rad + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 2106.3 2159.5 -1040.1 2080.3 431
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0347 -0.5637 -0.0940 0.5555 2.7091
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.008162 0.09034
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.384683 0.219874 6.298 3.02e-10 ***
## settlementsNear 0.516578 0.078388 6.590 4.40e-11 ***
## year_ct 0.003225 0.011715 0.275 0.783123
## subunitCarmel 0.161534 0.105491 1.531 0.125707
## subunitGalilee 0.097938 0.110020 0.890 0.373367
## cos_td_rad 0.348390 0.151344 2.302 0.021337 *
## sin_td_rad -0.314279 0.178202 -1.764 0.077797 .
## settlementsNear:year_ct -0.003232 0.010891 -0.297 0.766631
## year_ct:subunitCarmel -0.006346 0.012992 -0.488 0.625254
## year_ct:subunitGalilee -0.025262 0.014608 -1.729 0.083754 .
## settlementsNear:subunitCarmel -0.078696 0.080807 -0.974 0.330118
## settlementsNear:subunitGalilee -0.311468 0.086784 -3.589 0.000332 ***
## ---
## 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_ sttN:_ yr_c:C
## settlmntsNr -0.229
## year_ct -0.231 0.384
## subunitCrml -0.357 0.266 0.295
## subunitGall -0.395 0.284 0.297 0.519
## cos_td_rad -0.911 0.005 0.043 0.102 0.139
## sin_td_rad 0.871 -0.008 0.096 -0.145 -0.200 -0.857
## sttlmntsN:_ 0.157 -0.664 -0.579 -0.008 -0.052 -0.009 0.012
## yr_ct:sbntC 0.323 -0.012 -0.495 -0.587 -0.330 -0.161 0.251 0.017
## yr_ct:sbntG 0.358 -0.047 -0.463 -0.301 -0.621 -0.203 0.298 0.071 0.516
## sttlmntsN:C 0.114 -0.543 -0.001 -0.474 -0.241 0.007 -0.008 0.000 0.004
## sttlmntsN:G 0.120 -0.520 -0.010 -0.237 -0.462 -0.002 0.008 0.022 0.003
## yr_c:G sttN:C
## settlmntsNr
## year_ct
## subunitCrml
## subunitGall
## cos_td_rad
## sin_td_rad
## sttlmntsN:_
## yr_ct:sbntC
## yr_ct:sbntG
## sttlmntsN:C -0.002
## sttlmntsN:G 0.014 0.490
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0171254 (tol = 0.002, component 1)perform 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.00297726 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00500952 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00464399 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00212765 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## richness ~ settlements * year_ct + subunit * year_ct + settlements *
## subunit + cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 2106.3
## cos_td_rad 1 2109.7
## sin_td_rad 1 2107.4
## settlements:year_ct 1 2104.3
## year_ct:subunit 2 2105.5
## settlements:subunit 2 2115.9remove settlements X year.
me1 <- glmer(data = P.anal, formula = richness ~ subunit * year_ct + settlements * subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00576349 (tol = 0.002, component 1)drop1(me1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00494382 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## richness ~ subunit * year_ct + settlements * subunit + cos_td_rad +
## sin_td_rad + (1 | site)
## npar AIC
## <none> 2104.3
## cos_td_rad 1 2107.7
## sin_td_rad 1 2105.4
## subunit:year_ct 2 2103.5
## subunit:settlements 2 2114.0drop subunit X year
me1 <- glmer(data = P.anal, formula = richness ~ year_ct + settlements * subunit + cos_td_rad + sin_td_rad + (1|site), family = poisson)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00490422 (tol = 0.002, component 1)drop1(me1)## Single term deletions
##
## Model:
## richness ~ year_ct + settlements * subunit + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 2103.5
## year_ct 1 2102.1
## cos_td_rad 1 2105.7
## sin_td_rad 1 2103.4
## settlements:subunit 2 2113.0drop 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> 2102.1
## cos_td_rad 1 2104.0
## sin_td_rad 1 2101.4
## settlements:subunit 2 2111.5me1 <- 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> 2101.4
## cos_td_rad 1 2104.1
## settlements:subunit 2 2110.8interaction between settlements and subunit 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
## 2101.4 2134.2 -1042.7 2085.4 436
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.11475 -0.57133 -0.08616 0.56975 2.86205
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009297 0.09642
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.73590 0.07833 22.161 < 2e-16 ***
## settlementsNear 0.50123 0.05865 8.547 < 2e-16 ***
## subunitCarmel 0.13210 0.08801 1.501 0.133368
## subunitGalilee -0.02205 0.08867 -0.249 0.803605
## cos_td_rad 0.13449 0.06301 2.134 0.032811 *
## settlementsNear:subunitCarmel -0.08012 0.08080 -0.992 0.321372
## settlementsNear:subunitGalilee -0.30862 0.08675 -3.558 0.000374 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl cs_td_ sttN:C
## settlmntsNr -0.468
## subunitCrml -0.594 0.417
## subunitGall -0.575 0.414 0.515
## cos_td_rad -0.588 -0.001 0.020 -0.004
## sttlmntsN:C 0.342 -0.726 -0.565 -0.301 -0.003
## sttlmntsN:G 0.315 -0.676 -0.282 -0.562 0.003 0.491ranef(me1)## $site
## (Intercept)
## Abirim -7.039279e-02
## Aderet 5.965112e-02
## Beit Oren -2.170240e-02
## Ein Yaakov -2.735210e-02
## Givat Yearim -7.520626e-02
## Givat Yeshayahu 9.290074e-02
## Goren -9.374680e-03
## Iftach 1.704414e-01
## Kerem Maharal 9.581824e-02
## Kfar Shamai -5.626600e-02
## Margaliot -7.456596e-04
## Nehusha -4.231655e-02
## Nir Etzion -1.509844e-01
## Ofer 8.169541e-02
## Ramat Raziel -3.012056e-02
## Yagur -8.290748e-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
## 2103.4 2140.3 -1042.7 2085.4 435
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.12043 -0.57650 -0.08122 0.55700 2.85914
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009248 0.09617
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.746047 0.100486 17.376 < 2e-16 ***
## year_ct -0.001054 0.006548 -0.161 0.872183
## settlementsNear 0.501269 0.058646 8.547 < 2e-16 ***
## subunitCarmel 0.131946 0.087902 1.501 0.133343
## subunitGalilee -0.022122 0.088572 -0.250 0.802767
## cos_td_rad 0.127410 0.076851 1.658 0.097342 .
## settlementsNear:subunitCarmel -0.080185 0.080801 -0.992 0.321011
## settlementsNear:subunitGalilee -0.308758 0.086756 -3.559 0.000372 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) yer_ct sttlmN sbntCr sbntGl cs_td_ sttN:C
## year_ct -0.627
## settlmntsNr -0.362 -0.004
## subunitCrml -0.469 0.011 0.418
## subunitGall -0.451 0.005 0.415 0.515
## cos_td_rad -0.735 0.573 -0.003 0.022 -0.001
## sttlmntsN:C 0.263 0.005 -0.726 -0.566 -0.301 0.000
## sttlmntsN:G 0.239 0.010 -0.676 -0.282 -0.563 0.008 0.491year 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
## 2101.4 2134.2 -1042.7 2085.4 436
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.11475 -0.57133 -0.08615 0.56974 2.86206
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.009297 0.09642
## Number of obs: 444, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.83136 0.06338 28.896 < 2e-16 ***
## settlementsNear 0.50124 0.05865 8.547 < 2e-16 ***
## subunitCarmel 0.13209 0.08801 1.501 0.133380
## subunitGalilee -0.02206 0.08867 -0.249 0.803561
## cos_td_rad_c 0.13450 0.06301 2.134 0.032808 *
## settlementsNear:subunitCarmel -0.08012 0.08080 -0.992 0.321373
## settlementsNear:subunitGalilee -0.30862 0.08675 -3.558 0.000374 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN sbntCr sbntGl cs_t__ sttN:C
## settlmntsNr -0.579
## subunitCrml -0.720 0.417
## subunitGall -0.714 0.414 0.515
## cos_td_rd_c -0.021 -0.001 0.020 -0.004
## sttlmntsN:C 0.420 -0.726 -0.565 -0.301 -0.003
## sttlmntsN:G 0.391 -0.676 -0.282 -0.562 0.003 0.491P.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 | 2101.431 |
| BIC | 2134.197 |
| Pseudo-R² (fixed effects) | 0.293 |
| Pseudo-R² (total) | 0.344 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 6.242 | 0.063 | 28.896 | 0.000 |
| settlementsNear | 1.651 | 0.059 | 8.547 | 0.000 |
| subunitCarmel | 1.141 | 0.088 | 1.501 | 0.133 |
| subunitGalilee | 0.978 | 0.089 | -0.249 | 0.804 |
| cos_td_rad_c | 1.144 | 0.063 | 2.134 | 0.033 |
| settlementsNear:subunitCarmel | 0.923 | 0.081 | -0.992 | 0.321 |
| settlementsNear:subunitGalilee | 0.734 | 0.087 | -3.558 | 0.000 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.096 |
| Group | # groups | ICC |
|---|---|---|
| site | 16 | 0.009 |
cat_plot(me1, pred = settlements, modx = subunit, interval = TRUE)
plot_model_interaction_cat (P.anal = P.anal, m = me1, eff2plot = "settlements", modvar2plot = "subunit",
export_plot = T, plot_points = T, fontname = "Almoni ML v5 AAA",
outpath = "../output/maquis/")## Loading required package: Cairo## Loading required package: extrafont## Registering fonts with R## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
#emmeans
EMM <- emmeans(object = me1, ~subunit*settlements)
print(EMM)## subunit settlements emmean SE df asymp.LCL asymp.UCL
## Judean Highlands Far 1.83 0.0634 Inf 1.71 1.96
## Carmel Far 1.96 0.0611 Inf 1.84 2.08
## Galilee Far 1.81 0.0621 Inf 1.69 1.93
## Judean Highlands Near 2.33 0.0561 Inf 2.22 2.44
## Carmel Near 2.38 0.0556 Inf 2.28 2.49
## Galilee Near 2.00 0.0600 Inf 1.88 2.12
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Judea near vs far
(exp(2.33)/exp(1.83) -1) * 100## [1] 64.87213#Carmel near vs far
(exp(2.38)/exp(1.96) -1) * 100## [1] 52.19616#Galilee near vs far
(exp(2.00)/exp(1.81) -1) * 100## [1] 20.92496test_results_settlements <- test(pairs(EMM, by="subunit"), by=NULL, adjust="fdr")
print(test_results_settlements)## contrast subunit estimate SE df z.ratio p.value
## Far - Near Judean Highlands -0.501 0.0586 Inf -8.547 <.0001
## Far - Near Carmel -0.421 0.0556 Inf -7.576 <.0001
## Far - Near Galilee -0.193 0.0639 Inf -3.013 0.0026
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 teststest_results_subunit <- test(pairs(EMM, by="settlements"), by=NULL, adjust="fdr")
print(test_results_subunit)## contrast settlements estimate SE df z.ratio p.value
## Judean Highlands - Carmel Far -0.1321 0.0880 Inf -1.501 0.2001
## Judean Highlands - Galilee Far 0.0221 0.0887 Inf 0.249 0.8036
## Carmel - Galilee Far 0.1541 0.0871 Inf 1.771 0.1532
## Judean Highlands - Carmel Near -0.0520 0.0790 Inf -0.658 0.6126
## Judean Highlands - Galilee Near 0.3307 0.0821 Inf 4.028 0.0002
## Carmel - Galilee Near 0.3826 0.0817 Inf 4.681 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 6 tests# iv1_abs <- diff(effplot1$data$richness)
# iv1_rel <- iv1_abs / as.data.table(effplot1$data)[settlements=="Far",richness] * 100
#
# iv2_abs <- as.data.table(effplot2$data)[subunit!="Galilee"][,mean(richness)]-as.data.table(effplot2$data)[subunit=="Galilee",richness]
# iv2_rel <- iv2_abs / as.data.table(effplot2$data)[subunit!="Galilee"][,mean(richness)] * 100
#
# m_coef <- fixef(me1)
# X <- unique(data.table(samp_date = as_date(as.duration(P.anal$timediff_Jun21+172)), contrib = exp(m_coef["cos_td_rad_c"]*P.anal$cos_td_rad_c))[order(samp_date)])
# plot(X$samp_date, X$contrib,xlab = "time of year", ylab = "Effect on richness (multiplier)", type = "b")statistically 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 * year_ct + subunit * settlements + 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 * year_ct + subunit * settlements + cos_td_rad + sin_td_rad)
me0 <- lmer(data = P.anal, formula = gma ~ settlements * year_ct + subunit * year_ct + subunit * settlements + 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 * year_ct + subunit * settlements +
## cos_td_rad + sin_td_rad + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 1300.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0787 -0.6105 -0.1632 0.4070 5.8272
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.07482 0.2735
## Residual 1.03131 1.0155
## Number of obs: 440, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.685947 0.630081 1.089
## settlementsNear 0.301688 0.223451 1.350
## year_ct 0.015827 0.032149 0.492
## subunitCarmel 0.458922 0.304589 1.507
## subunitGalilee -0.002097 0.305651 -0.007
## cos_td_rad 1.291527 0.432120 2.989
## sin_td_rad -1.686855 0.522970 -3.226
## settlementsNear:year_ct -0.028772 0.031040 -0.927
## year_ct:subunitCarmel -0.063263 0.039188 -1.614
## year_ct:subunitGalilee -0.031299 0.040874 -0.766
## settlementsNear:subunitCarmel -0.039179 0.236155 -0.166
## settlementsNear:subunitGalilee 0.148822 0.237935 0.625
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct sbntCr sbntGl cs_td_ sn_td_ sttN:_ yr_c:C
## settlmntsNr -0.183
## year_ct -0.212 0.317
## subunitCrml -0.361 0.206 0.318
## subunitGall -0.397 0.206 0.306 0.519
## cos_td_rad -0.918 0.005 0.041 0.116 0.148
## sin_td_rad 0.880 -0.008 0.108 -0.158 -0.208 -0.862
## sttlmntsN:_ 0.123 -0.664 -0.484 -0.003 -0.003 -0.006 0.006
## yr_ct:sbntC 0.336 -0.006 -0.513 -0.615 -0.341 -0.176 0.259 0.010
## yr_ct:sbntG 0.371 -0.004 -0.481 -0.324 -0.630 -0.215 0.308 0.007 0.530
## sttlmntsN:C 0.092 -0.527 0.003 -0.386 -0.192 0.001 -0.001 -0.002 0.002
## sttlmntsN:G 0.100 -0.524 0.006 -0.193 -0.385 -0.006 0.011 -0.001 0.001
## yr_c:G sttN:C
## settlmntsNr
## year_ct
## subunitCrml
## subunitGall
## cos_td_rad
## sin_td_rad
## sttlmntsN:_
## yr_ct:sbntC
## yr_ct:sbntG
## sttlmntsN:C -0.002
## sttlmntsN:G 0.001 0.496perform stepwise model selection of gaussian model.
drop1(me0)## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + subunit * year_ct + subunit * settlements +
## cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 1295.1
## cos_td_rad 1 1302.1
## sin_td_rad 1 1303.8
## settlements:year_ct 1 1294.0
## year_ct:subunit 2 1293.8
## settlements:subunit 2 1291.8drop subunit X year
me1 <- lmer(data = P.anal, formula = gma ~ settlements * year_ct + subunit * settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + subunit * settlements + cos_td_rad +
## sin_td_rad + (1 | site)
## npar AIC
## <none> 1293.8
## cos_td_rad 1 1299.6
## sin_td_rad 1 1301.0
## settlements:year_ct 1 1292.6
## settlements:subunit 2 1290.5drop subunit X settlements
me1 <- lmer(data = P.anal, formula = gma ~ settlements * year_ct + subunit + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## 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 settlements X year
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 + sin_td_rad + cos_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ settlements + sin_td_rad + cos_td_rad + (1 | site)
## npar AIC
## <none> 1286.0
## settlements 1 1288.2
## sin_td_rad 1 1291.9
## cos_td_rad 1 1291.1drop settlements
me1 <- lmer(data = P.anal, formula = gma ~ sin_td_rad + cos_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## gma ~ sin_td_rad + cos_td_rad + (1 | site)
## npar AIC
## <none> 1288.2
## sin_td_rad 1 1294.1
## cos_td_rad 1 1293.3only time of year remains in the model.
This is the final model:
summary(me1)## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ sin_td_rad + cos_td_rad + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 1282.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9564 -0.6075 -0.1802 0.4110 5.9690
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.06916 0.263
## Residual 1.03746 1.019
## Number of obs: 440, groups: site, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.2541 0.5452 2.300
## sin_td_rad -1.2317 0.4385 -2.809
## cos_td_rad 1.1108 0.4175 2.660
##
## Correlation of Fixed Effects:
## (Intr) sn_td_
## sin_td_rad 0.960
## cos_td_rad -0.968 -0.904plot(me1, which=1:3)
# summ(me1)
#
# 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/")
# effplot1abundance
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 * year_ct + subunit * settlements + 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 * year_ct + subunit * settlements + cos_td_rad + sin_td_rad)
AIC(m0)## [1] 3687.163me0 <- glmer.nb(data = P.anal, formula = abundance ~ settlements * year_ct + subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad + (1|site))
AIC(me0)## [1] 3679.845mixed 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.00331907 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## abundance ~ settlements * year_ct + subunit * year_ct + subunit *
## settlements + cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 3679.8
## cos_td_rad 1 3697.8
## sin_td_rad 1 3700.4
## settlements:year_ct 1 3677.9
## year_ct:subunit 2 3682.9
## settlements:subunit 2 3678.7drop subunit X year
me1 <- glmer.nb(data = P.anal, formula = abundance ~ settlements * year_ct + subunit * settlements + 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.00242417 (tol = 0.002, component 1)drop1(me1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00250601 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## abundance ~ settlements * year_ct + subunit * settlements + cos_td_rad +
## sin_td_rad + (1 | site)
## npar AIC
## <none> 3682.9
## cos_td_rad 1 3698.3
## sin_td_rad 1 3700.7
## settlements:year_ct 1 3681.0
## settlements:subunit 2 3681.6drop settlements X year
me1 <- glmer.nb(data = P.anal, formula = abundance ~ year_ct + subunit * settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## abundance ~ year_ct + subunit * settlements + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3681.0
## year_ct 1 3679.5
## cos_td_rad 1 3696.4
## sin_td_rad 1 3698.8
## subunit:settlements 2 3679.7drop year
me1 <- glmer.nb(data = P.anal, formula = abundance ~ subunit * settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## abundance ~ subunit * settlements + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3679.5
## cos_td_rad 1 3694.4
## sin_td_rad 1 3699.0
## subunit:settlements 2 3678.2drop settlements X subunit
me1 <- glmer.nb(data = P.anal, formula = abundance ~ subunit + settlements + cos_td_rad + sin_td_rad + (1|site))
drop1(me1)## Single term deletions
##
## Model:
## abundance ~ subunit + settlements + cos_td_rad + sin_td_rad +
## (1 | site)
## npar AIC
## <none> 3678.2
## subunit 2 3683.8
## settlements 1 3772.0
## cos_td_rad 1 3693.2
## sin_td_rad 1 3697.8Subunit, settlement and time of year remain, same as no interaction analysis.
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*year_ct + subunit*settlements + cos_td_rad + sin_td_rad, family = "poisson", data = env_data)
mva_m0.nb <- manyglm(formula = spp_no_rare ~ settlements*year_ct + subunit*year_ct + subunit*settlements + 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
## 976.7852 1371.0086print("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*year_ct + subunit*settlements + cos_td_rad + sin_td_rad + site, data = env_data)
aic_dt$nb1 <- mva_m1$aic
colMeans(aic_dt)## nb po nb1
## 976.7852 1371.0086 966.7749The 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 * year_ct + subunit *
## settlements + cos_td_rad + sin_td_rad
## Df AIC
## <none> 23443
## cos_td_rad 24 23501
## sin_td_rad 24 23515
## settlements:year_ct 24 23426
## year_ct:subunit 48 23471
## settlements:subunit 48 23458drop settlements X year.
mva_m1 <- manyglm(formula = spp_no_rare ~ subunit*year_ct + subunit*settlements + cos_td_rad + sin_td_rad, family = "negative.binomial", data = env_data)
drop1(mva_m1)## Single term deletions
##
## Model:
## spp_no_rare ~ subunit * year_ct + subunit * settlements + cos_td_rad +
## sin_td_rad
## Df AIC
## <none> 23426
## cos_td_rad 24 23484
## sin_td_rad 24 23501
## subunit:year_ct 48 23460
## subunit:settlements 48 23439final model includes subunit x settlements, subunit x year_ct and sampling time of year.
plot(mva_m1, which=1:3)
#summ_m1 <- summary(mva_m1)
#saveRDS(summ_m1, file = "../output/Maquis_new_MVabund_interactions_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_interactions_anova_output_99iter.rds")
summ_m1 <- readRDS(file = "../output/20230831_Maquis_new_MVabund_interactions_summary_output_999iter.rds")
anov_m1.uni <- readRDS("../output/20230831_Maquis_new_MVabund_interactions_anova_output_99iter.rds")
print(summ_m1)##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 12.183 0.001 ***
## subunitCarmel 8.923 0.001 ***
## subunitGalilee 6.796 0.010 **
## year_ct 6.584 0.036 *
## settlementsNear 18.449 0.001 ***
## cos_td_rad 10.665 0.001 ***
## sin_td_rad 11.531 0.001 ***
## subunitCarmel:year_ct 7.738 0.006 **
## subunitGalilee:year_ct 7.594 0.004 **
## subunitCarmel:settlementsNear 4.572 0.538
## subunitGalilee:settlementsNear 7.429 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 39.66, 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 ~ subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 443
## subunit 441 2 393.0 0.01 **
## year_ct 440 1 125.5 0.01 **
## settlements 439 1 821.9 0.01 **
## cos_td_rad 438 1 46.9 0.03 *
## sin_td_rad 437 1 109.1 0.01 **
## subunit:year_ct 435 2 137.4 0.01 **
## subunit:settlements 433 2 109.8 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)
## subunit 44.029 0.01 9.54 0.24
## year_ct 24.181 0.01 12.413 0.01
## settlements 80.478 0.01 36.319 0.01
## cos_td_rad 3.002 0.93 2.276 0.97
## sin_td_rad 5.64 0.35 0.596 0.99
## subunit:year_ct 8.217 0.39 2.436 0.90
## subunit:settlements 12.215 0.13 0.135 1.00
## Carduelis.carduelis Cecropis.daurica
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 18.526 0.03 0.804 0.91
## year_ct 7.364 0.23 0.498 1.00
## settlements 9.786 0.04 34.435 0.01
## cos_td_rad 0.389 1.00 0.127 1.00
## sin_td_rad 1.61 0.95 3.044 0.78
## subunit:year_ct 12.307 0.11 0.41 1.00
## subunit:settlements 0.982 1.00 0.481 1.00
## Chloris.chloris Cinnyris.osea
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 0.606 0.91 18.184 0.03
## year_ct 7.455 0.23 4.452 0.48
## settlements 28.175 0.01 41.118 0.01
## cos_td_rad 0.072 1.00 1.651 0.98
## sin_td_rad 0.112 0.99 2.01 0.91
## subunit:year_ct 20.167 0.02 12.677 0.11
## subunit:settlements 3.829 0.96 0.91 1.00
## Columba.livia Corvus.cornix
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 11.932 0.11 34.187 0.01
## year_ct 7.544 0.22 0.269 1.00
## settlements 27.033 0.01 36.648 0.01
## cos_td_rad 0.01 1.00 1.013 0.99
## sin_td_rad 11.564 0.03 0.157 0.99
## subunit:year_ct 2.764 0.90 7.109 0.45
## subunit:settlements 8.924 0.34 21.928 0.01
## Corvus.monedula Curruca.curruca
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 23.137 0.01 26.636 0.01
## year_ct 0.306 1.00 12.975 0.01
## settlements 0.001 0.99 3.273 0.52
## cos_td_rad 0.184 1.00 7.65 0.28
## sin_td_rad 0.053 0.99 4.442 0.55
## subunit:year_ct 5.316 0.78 0.001 1.00
## subunit:settlements 9.568 0.27 1.179 1.00
## Curruca.melanocephala Dendrocopos.syriacus
## Dev Pr(>Dev) Dev
## (Intercept)
## subunit 38.349 0.01 3.624
## year_ct 0.745 0.98 0.002
## settlements 23.96 0.01 9.308
## cos_td_rad 0.3 1.00 1.417
## sin_td_rad 0.429 0.99 3.35
## subunit:year_ct 9.877 0.27 0.232
## subunit:settlements 0.509 1.00 5.453
## Falco.tinnunculus Garrulus.glandarius
## Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## subunit 0.79 8.726 0.26 26.08
## year_ct 1.00 0.424 1.00 4.802
## settlements 0.04 2.364 0.56 2.168
## cos_td_rad 0.99 0.078 1.00 1.156
## sin_td_rad 0.77 0.641 0.98 1.081
## subunit:year_ct 1.00 5.248 0.78 5.242
## subunit:settlements 0.78 0.526 1.00 8.69
## Parus.major Passer.domesticus
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 0.01 26.708 0.01 18.116 0.03
## year_ct 0.43 0.357 1.00 0.957 0.97
## settlements 0.56 2.626 0.55 170.553 0.01
## cos_td_rad 0.99 3.376 0.91 0.485 1.00
## sin_td_rad 0.97 7.701 0.19 2.962 0.78
## subunit:year_ct 0.78 6.652 0.57 2.283 0.90
## subunit:settlements 0.34 4.806 0.91 4.667 0.91
## Prinia.gracilis Psittacula.krameri
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## subunit 1.649 0.87 12.826 0.10
## year_ct 10.434 0.08 1.303 0.92
## settlements 54.604 0.01 16.939 0.01
## cos_td_rad 0.633 1.00 0.793 1.00
## sin_td_rad 1.004 0.97 1.724 0.95
## subunit:year_ct 13.996 0.06 3.863 0.85
## subunit:settlements 7.053 0.46 0.813 1.00
## Pycnonotus.xanthopygos Spilopelia.senegalensis
## Dev Pr(>Dev) Dev
## (Intercept)
## subunit 22.998 0.01 20.899
## year_ct 4.112 0.53 2.986
## settlements 46.389 0.01 162.398
## cos_td_rad 11.394 0.10 0.333
## sin_td_rad 0.499 0.99 0.603
## subunit:year_ct 3.75 0.85 4.692
## subunit:settlements 2.425 1.00 0.038
## Streptopelia.decaocto Streptopelia.turtur
## Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## subunit 0.01 12.15 0.11 6.199
## year_ct 0.64 0.079 1.00 18.03
## settlements 0.01 1.962 0.56 1.6
## cos_td_rad 1.00 6.672 0.38 3.045
## sin_td_rad 0.99 14.008 0.02 1.847
## subunit:year_ct 0.80 1.683 0.91 1.887
## subunit:settlements 1.00 7.741 0.42 2.078
## Troglodytes.troglodytes Turdus.merula
## Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## subunit 0.44 4.791 0.63 2.302
## year_ct 0.01 0.88 0.97 2.975
## settlements 0.56 25.789 0.01 3.928
## cos_td_rad 0.93 0.48 1.00 0.347
## sin_td_rad 0.95 10.786 0.03 33.19
## subunit:year_ct 0.91 5.003 0.78 1.539
## subunit:settlements 1.00 2.959 1.00 1.928
##
## Pr(>Dev)
## (Intercept)
## subunit 0.84
## year_ct 0.64
## settlements 0.39
## cos_td_rad 1.00
## sin_td_rad 0.01
## subunit:year_ct 0.91
## subunit:settlements 1.00
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 99 iterations via PIT-trap resampling.Interactions are significant only for three species: subunit X year: prinia gracilis and chloris chloris subunit X settlements: corvus cornix Therefore present only the analysis EXCLUDING subunit interactions with time / settlements
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/")
# 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/")
#
# plot_coefs_year
# 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_yearsubuChloris chloris
#####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~subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad, data = d.chloris)
coef(mva_m1)[,"Chloris.chloris"]## (Intercept) subunitCarmel
## -1.38223485 -3.31810481
## subunitGalilee year_ct
## 0.28573227 0.04493925
## settlementsNear cos_td_rad
## 1.31113634 0.45112987
## sin_td_rad subunitCarmel:year_ct
## 0.19135232 0.37580003
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.06881828 1.04485733
## subunitGalilee:settlementsNear
## -0.19537232coef(m.chloris)## (Intercept) subunitCarmel
## -1.38223300 -3.31810335
## subunitGalilee year_ct
## 0.28573282 0.04493950
## settlementsNear cos_td_rad
## 1.31113632 0.45112850
## sin_td_rad subunitCarmel:year_ct
## 0.19135596 0.37579996
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.06881835 1.04485589
## subunitGalilee:settlementsNear
## -0.19537267plot_model_interaction(P.anal = d.chloris, m = m.chloris, eff2plot = "year_ct", modvar2plot = "subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "../output/maquis/")
chloris_interact <- interact_plot(model = m.chloris, data=d.chloris, pred = year_ct, modx=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 = "Temporal effect on abundance of Prinia gracilis")
chloris_interact
emm_year_ct <- emtrends(m.chloris, specs = "subunit", var = "year_ct", type = "response")
print(emm_year_ct)## subunit year_ct.trend SE df asymp.LCL asymp.UCL
## Judean Highlands 0.0449 0.0676 Inf -0.0876 0.177
## Carmel 0.4207 0.1022 Inf 0.2205 0.621
## Galilee -0.0239 0.0844 Inf -0.1893 0.142
##
## Results are averaged over the levels of: settlements
## Confidence level used: 0.95test_year_ct <- test(emm_year_ct, null=0, adjust="fdr")
print(test_year_ct)## subunit year_ct.trend SE df z.ratio p.value
## Judean Highlands 0.0449 0.0676 Inf 0.665 0.7595
## Carmel 0.4207 0.1022 Inf 4.118 0.0001
## Galilee -0.0239 0.0844 Inf -0.283 0.7773
##
## Results are averaged over the levels of: settlements
## P value adjustment: fdr method for 3 testsPrinia gracilis
#####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~subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad, data = d.prigra)
coef(mva_m1)[,"Prinia.gracilis"]## (Intercept) subunitCarmel
## -2.26251893 1.38322101
## subunitGalilee year_ct
## 1.93593465 0.01501659
## settlementsNear cos_td_rad
## 2.38276491 0.26398831
## sin_td_rad subunitCarmel:year_ct
## 0.01373954 -0.07295477
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.25436101 -1.00105072
## subunitGalilee:settlementsNear
## -1.31478216coef(m.prigra)## (Intercept) subunitCarmel
## -2.26251907 1.38322118
## subunitGalilee year_ct
## 1.93593773 0.01501663
## settlementsNear cos_td_rad
## 2.38276513 0.26398826
## sin_td_rad subunitCarmel:year_ct
## 0.01373969 -0.07295484
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.25436150 -1.00105054
## subunitGalilee:settlementsNear
## -1.31478370plot_model_interaction(P.anal = d.prigra, m = m.prigra, eff2plot = "year_ct", modvar2plot = "subunit", plot_points=FALSE, plot_residuals=FALSE, export_plot=TRUE, ylabel = NULL, fontname = "Almoni ML v5 AAA", fontsize=22, pdf_width=160, outpath = "../output/maquis/")
prigra_interact <- interact_plot(model = m.prigra, data=d.prigra, pred = year_ct, modx=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 = "Temporal effect on abundance of Prinia gracilis")
prigra_interact
#Decrease of 40.63%
carmel_time <- as.data.table(prigra_interact$data)[subunit=="Carmel",Prinia.gracilis]
print(carmel_time)## [1] 0.4965926 0.4939838 0.4913888 0.4888074 0.4862396 0.4836852 0.4811443
## [8] 0.4786167 0.4761024 0.4736013 0.4711134 0.4686385 0.4661766 0.4637277
## [15] 0.4612916 0.4588683 0.4564578 0.4540599 0.4516746 0.4493018 0.4469415
## [22] 0.4445936 0.4422580 0.4399347 0.4376236 0.4353247 0.4330378 0.4307629
## [29] 0.4285000 0.4262490 0.4240098 0.4217824 0.4195666 0.4173626 0.4151700
## [36] 0.4129890 0.4108195 0.4086614 0.4065146 0.4043790 0.4022547 0.4001416
## [43] 0.3980395 0.3959485 0.3938685 0.3917994 0.3897412 0.3876938 0.3856571
## [50] 0.3836311 0.3816158 0.3796111 0.3776169 0.3756332 0.3736599 0.3716969
## [57] 0.3697443 0.3678020 0.3658698 0.3639478 0.3620359 0.3601340 0.3582421
## [64] 0.3563602 0.3544881 0.3526259 0.3507735 0.3489308 0.3470977 0.3452743
## [71] 0.3434605 0.3416562 0.3398614 0.3380761 0.3363001 0.3345334 0.3327760
## [78] 0.3310278 0.3292888 0.3275590 0.3258383 0.3241265 0.3224238 0.3207300
## [85] 0.3190452 0.3173691 0.3157019 0.3140434 0.3123937 0.3107526 0.3091201
## [92] 0.3074963 0.3058809 0.3042740 0.3026756 0.3010856 0.2995039 0.2979305
## [99] 0.2963654 0.29480851 - 0.2948085/0.4965926 ## [1] 0.4063373(max(carmel_time) - min(carmel_time))/max(carmel_time)## [1] 0.4063373#Decreace of 88.39%
galilee_time <- as.data.table(prigra_interact$data)[subunit=="Galilee",Prinia.gracilis]
(max(galilee_time) - min(galilee_time))/max(galilee_time)## [1] 0.8839929#Increase of 14.47%
judea_time <- as.data.table(prigra_interact$data)[subunit=="Judean Highlands",Prinia.gracilis]
(max(judea_time) - min(judea_time))/min(judea_time)## [1] 0.1447081emm_year_ct <- emtrends(m.prigra, specs = "subunit", var = "year_ct", type = "response")
print(emm_year_ct)## subunit year_ct.trend SE df asymp.LCL asymp.UCL
## Judean Highlands 0.0150 0.0537 Inf -0.0903 0.1203
## Carmel -0.0579 0.0553 Inf -0.1663 0.0505
## Galilee -0.2393 0.0667 Inf -0.3701 -0.1086
##
## Results are averaged over the levels of: settlements
## Confidence level used: 0.95test_year_ct <- test(emm_year_ct, null=0, adjust="fdr")
print(test_year_ct)## subunit year_ct.trend SE df z.ratio p.value
## Judean Highlands 0.0150 0.0537 Inf 0.279 0.7799
## Carmel -0.0579 0.0553 Inf -1.048 0.4422
## Galilee -0.2393 0.0667 Inf -3.587 0.0010
##
## Results are averaged over the levels of: settlements
## P value adjustment: fdr method for 3 testsCorvus cornix
####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 ~ subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad, data = d.corvus)
coef(mva_m1)[,"Corvus.cornix"]## (Intercept) subunitCarmel
## -15.41907632 15.69624138
## subunitGalilee year_ct
## 13.87109757 0.12635685
## settlementsNear cos_td_rad
## 14.23172291 0.17633434
## sin_td_rad subunitCarmel:year_ct
## 0.21810658 -0.22276179
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.07136465 -13.53439460
## subunitGalilee:settlementsNear
## -13.01077692coef(m.corvus)## (Intercept) subunitCarmel
## -33.95095922 34.22813307
## subunitGalilee year_ct
## 32.40298781 0.12636112
## settlementsNear cos_td_rad
## 32.76359461 0.17632863
## sin_td_rad subunitCarmel:year_ct
## 0.21811556 -0.22276607
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.07136854 -32.06626545
## subunitGalilee:settlementsNear
## -31.54264820judea_far <- as.data.table(d.corvus)[subunit=="Judea"]
judea_far <- as.data.table(judea_far)[settlements=="Far"]
sum(judea_far$Corvus.cornix)## [1] 0corcor_interact <- cat_plot(m.corvus, pred = settlements, modx = subunit, interval = TRUE)
corcor_interact
EMM <- emmeans(object = m.corvus, ~ subunit*settlements)
print(EMM)## subunit settlements emmean SE df asymp.LCL asymp.UCL
## Judean Highlands Far -33.3472 1.19e+06 Inf -2.34e+06 2336490.5
## Carmel Far -0.1883 2.00e-01 Inf -6.00e-01 0.2
## Galilee Far -1.2868 3.00e-01 Inf -1.80e+00 -0.8
## Judean Highlands Near -0.5836 2.00e-01 Inf -1.00e+00 -0.1
## Carmel Near 0.5091 2.00e-01 Inf 2.00e-01 0.9
## Galilee Near -0.0658 2.00e-01 Inf -5.00e-01 0.3
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Judea near vs far
(exp(-0.5836)/exp(-33.3498) -1) *100## [1] 1.698947e+16#Carmel near vs far
(exp(0.5091)/exp(-0.1883) -1) *100## [1] 100.8524#Galilee near vs far
(exp(-0.0658)/exp(-1.2868) -1) *100## [1] 239.0577test_results_subunit <- test(pairs(EMM, by="subunit"), by=NULL, adjust="fdr")
print(test_results_subunit)## contrast subunit estimate SE df z.ratio p.value
## Far - Near Judean Highlands -32.764 1.19e+06 Inf 0.000 1.0000
## Far - Near Carmel -0.697 3.00e-01 Inf -2.523 0.0175
## Far - Near Galilee -1.221 3.00e-01 Inf -3.567 0.0011
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testsSpelopelia senegalensis
#####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 ~ subunit * year_ct + subunit * settlements + cos_td_rad + sin_td_rad, data = d.spilop)
coef(mva_m1)[,"Spilopelia.senegalensis"]## (Intercept) subunitCarmel
## -3.08146657 -0.20315471
## subunitGalilee year_ct
## -0.47620439 -0.02488994
## settlementsNear cos_td_rad
## 2.89634113 0.93835577
## sin_td_rad subunitCarmel:year_ct
## -1.05321858 0.02072174
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.12883003 -0.02949221
## subunitGalilee:settlementsNear
## -0.12900968coef(m.spilop)## (Intercept) subunitCarmel
## -3.08146699 -0.20315375
## subunitGalilee year_ct
## -0.47620426 -0.02488990
## settlementsNear cos_td_rad
## 2.89634111 0.93835602
## sin_td_rad subunitCarmel:year_ct
## -1.05321870 0.02072158
## subunitGalilee:year_ct subunitCarmel:settlementsNear
## -0.12883006 -0.02949293
## subunitGalilee:settlementsNear
## -0.12900959cat_plot(m.spilop, pred = settlements, modx = subunit, interval = TRUE)
EMM <- emmeans(object = m.spilop, ~ subunit*settlements)
print(EMM)## subunit settlements emmean SE df asymp.LCL asymp.UCL
## Judean Highlands Far -1.917 0.327 Inf -2.559 -1.276
## Carmel Far -2.021 0.341 Inf -2.688 -1.354
## Galilee Far -3.012 0.522 Inf -4.035 -1.989
## Judean Highlands Near 0.979 0.144 Inf 0.697 1.261
## Carmel Near 0.846 0.147 Inf 0.559 1.133
## Galilee Near -0.244 0.191 Inf -0.618 0.129
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95#Judea near vs far
(exp(0.979)/exp(-1.917) -1) *100## [1] 1710.159#Carmel near vs far
(exp(0.846 )/exp(-2.021) -1) *100## [1] 1658.419#Galilee near vs far
(exp(-0.244)/exp(-3.012) -1) *100## [1] 1492.675test_results_subunit <- test(pairs(EMM, by="subunit"), by=NULL, adjust="fdr")
print(test_results_subunit)## contrast subunit estimate SE df z.ratio p.value
## Far - Near Judean Highlands -2.90 0.357 Inf -8.102 <.0001
## Far - Near Carmel -2.87 0.371 Inf -7.733 <.0001
## Far - Near Galilee -2.77 0.551 Inf -5.020 <.0001
##
## 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)")