Birds_Loess_analysis_by_plot
Ron Chen
2023-11-26
## 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 EXCLUDEDLoess covered areas in the northern Negev
Factors are land use (KKL plantings, Bedouine agriculture and natural loess) and time. Total 5 campaigns, one of which is pilot in 2012 which does not contain the bedouin agriculture land use level. 3 levels of land use: KKL plantings: 5 campaigns 2012-2020, 9 plots in each campaign. Bedouin agriculture: 4 campaigns 2014-2020, plot numbers: 9, 11, 12 ,12 respective to years. Loess: 5 campaigns 2012-2020, plot numbers: 15,9,9,9,9 respective to years. Year 2012 had 3 sites that were not repeated: Goral (3 plots), Hatzerim Base (6 plots), Nevatim Base (3 plots). If need to throw out data, best is to throw out Nevatim and the 3 plots named Hatzerim Loess 1,2,3 (have low / outlying counts of biodiversity metrics). 5 sites with 9 plots per site (3 for each land use X proximity combination).
Raw data Total abundance: 6040 Number of observations: 1447 Total richness: 92
Filtered data Total abundance: 2448 Number of observations: 780 Total richness: 41
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.
## [1] "RICHNESS WITH RARE SPECIES"## [1] "data with 2012"## ℹ SHA-1 hash of file is "eef55df1963ba201f57cf44513f566298756fd4f"
## [1] "data without 2012"
## [1] "data with 2012"
## [1] "data without 2012"
| richness | year_ct | site | land_use | td_sc | cos_td_rad | sin_td_rad | h_from_sunrise | cos_hsun | sin_hsun | monitors_name | wind | precipitation | temperature | clouds | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 1.000 | Min. :0.000 | Sayeret Shaked :30 | Loess :51 | Min. :-1.5653 | Min. :0.02295 | Min. :-0.9997 | Length:140 | Min. :0.1951 | Min. :-0.1642 | Eyal Shochat : 30 | 0 : 0 | 0 : 0 | 0 : 0 | 0 : 0 | |
| 1st Qu.: 3.000 | 1st Qu.:2.000 | Nahal Ashan :24 | Bedouin Agriculture:44 | 1st Qu.:-0.7244 | 1st Qu.:0.34157 | 1st Qu.:-0.9201 | Class :difftime | 1st Qu.:0.8607 | 1st Qu.: 0.1331 | Adi Domer : 0 | 1 : 0 | 3 : 0 | 1 : 0 | 1 : 0 | |
| Median : 4.000 | Median :4.000 | Mishmar Hanegev:23 | KKL Plantings :45 | Median : 0.6826 | Median :0.85208 | Median :-0.5087 | Mode :numeric | Median :0.9542 | Median : 0.2990 | Asaf Mayrose : 0 | 2 : 0 | NA’s:140 | 2 : 0 | 2 : 0 | |
| Mean : 4.593 | Mean :4.214 | Givot Goral :15 | NA | Mean : 0.6387 | Mean :0.67279 | Mean :-0.4482 | NA | Mean :0.8921 | Mean : 0.3286 | Eliraz Dvir : 0 | 3 : 0 | NA | 3 : 0 | 3 : 0 | |
| 3rd Qu.: 6.000 | 3rd Qu.:6.000 | Eshel Hanasi :12 | NA | 3rd Qu.: 1.5220 | 3rd Qu.:0.90882 | 3rd Qu.:-0.1543 | NA | 3rd Qu.:0.9893 | 3rd Qu.: 0.5090 | Eliraz Dvir and Yoav Barak: 0 | NA’s:140 | NA | NA’s:140 | NA’s:140 | |
| Max. :15.000 | Max. :8.000 | Givot Bar :12 | NA | Max. : 4.9558 | Max. :1.00000 | Max. : 0.9845 | NA | Max. :1.0000 | Max. : 0.9808 | (Other) : 0 | NA | NA | NA | NA | |
| NA | NA | (Other) :24 | NA | NA | NA | NA | NA | NA’s :51 | NA’s :51 | NA’s :110 | NA | NA | NA | NA |
no observation for all 4 weather variables. many NAs for sampling time of day variables.exclude from model.
## Warning in cor(P.anal[, lapply(X = .SD, FUN = as.numeric), .SDcols =
## IVs[1:11]], : the standard deviation is zero| richness | year_ct | site | land_use | td_sc | cos_td_rad | sin_td_rad | h_from_sunrise | cos_hsun | sin_hsun | monitors_name | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| richness | 1.0000000 | 0.0358986 | -0.2439373 | 0.3623180 | -0.0031931 | 0.1716892 | -0.0409555 | -0.1092363 | 0.1844295 | -0.0750554 | NA |
| year_ct | 0.0358986 | 1.0000000 | -0.0131923 | 0.0784319 | -0.6891876 | -0.4527528 | -0.6804369 | -0.1097558 | 0.1168999 | -0.1052239 | NA |
| site | -0.2439373 | -0.0131923 | 1.0000000 | -0.1194452 | -0.1113288 | -0.1171673 | -0.0980130 | 0.1422600 | -0.0542126 | 0.1561513 | NA |
| land_use | 0.3623180 | 0.0784319 | -0.1194452 | 1.0000000 | -0.0996863 | 0.0851086 | -0.1008691 | -0.0772346 | 0.1350010 | -0.0453437 | NA |
| td_sc | -0.0031931 | -0.6891876 | -0.1113288 | -0.0996863 | 1.0000000 | 0.5324060 | 0.9841321 | 0.0203523 | -0.1285853 | -0.0237796 | NA |
| cos_td_rad | 0.1716892 | -0.4527528 | -0.1171673 | 0.0851086 | 0.5324060 | 1.0000000 | 0.4312002 | 0.0061055 | -0.1126261 | -0.0359190 | NA |
| sin_td_rad | -0.0409555 | -0.6804369 | -0.0980130 | -0.1008691 | 0.9841321 | 0.4312002 | 1.0000000 | 0.0492202 | -0.1557296 | 0.0031361 | NA |
| h_from_sunrise | -0.1092363 | -0.1097558 | 0.1422600 | -0.0772346 | 0.0203523 | 0.0061055 | 0.0492202 | 1.0000000 | -0.9303548 | 0.9897728 | NA |
| cos_hsun | 0.1844295 | 0.1168999 | -0.0542126 | 0.1350010 | -0.1285853 | -0.1126261 | -0.1557296 | -0.9303548 | 1.0000000 | -0.8708307 | NA |
| sin_hsun | -0.0750554 | -0.1052239 | 0.1561513 | -0.0453437 | -0.0237796 | -0.0359190 | 0.0031361 | 0.9897728 | -0.8708307 | 1.0000000 | NA |
| monitors_name | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
Fit Poisson glm, check for existence of overdispersion
## [1] "Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model."## [1] "od = 0.760798125335199"Overdispersion parameter is < 1. Poisson more appropriate. Compare Poisson and negative binomial.
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached## df AIC
## m0.po 16 586.8881
## m0.nb 17 588.8897## [1] "poisson"
## [1] "neg bin"
negative binomial did not converge because there is underdispersion and poisson is better.
go with mixed model.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 588.5 615.0 -285.2 570.5 131
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6859 -0.6177 -0.1608 0.5767 3.6774
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01245 0.1116
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.06782 0.17914 5.961 2.51e-09 ***
## land_useBedouin Agriculture -0.26481 0.23580 -1.123 0.2614
## land_useKKL Plantings 0.32321 0.16197 1.995 0.0460 *
## year_ct -0.01841 0.03094 -0.595 0.5518
## cos_td_rad 0.37660 0.16078 2.342 0.0192 *
## sin_td_rad -0.09158 0.13362 -0.685 0.4931
## land_useBedouin Agriculture:year_ct 0.10169 0.04355 2.335 0.0195 *
## land_useKKL Plantings:year_ct 0.01537 0.03572 0.430 0.6670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ sn_td_ l_BA:_
## lnd_sBdnAgr -0.222
## lnd_sKKLPln -0.297 0.383
## year_ct -0.492 0.269 0.454
## cos_td_rad -0.730 -0.160 -0.208 0.083
## sin_td_rad 0.147 0.013 0.144 0.403 -0.267
## lnd_sBAgr:_ 0.202 -0.805 -0.321 -0.534 0.123 0.043
## lnd_sKKLP:_ 0.233 -0.259 -0.735 -0.659 0.147 -0.061 0.492perform stepwise model selection of poisson mixed model.
## Single term deletions
##
## Model:
## richness ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 588.49
## cos_td_rad 1 592.12
## sin_td_rad 1 586.97
## land_use:year_ct 2 590.59remove sine of time of year.
## Single term deletions
##
## Model:
## richness ~ land_use * year_ct + cos_td_rad + (1 | site)
## npar AIC
## <none> 586.97
## cos_td_rad 1 590.14
## land_use:year_ct 2 589.43land_use * year and cosine of time of year remain. Final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ land_use * year_ct + cos_td_rad + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 587.0 610.5 -285.5 571.0 132
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6438 -0.6020 -0.1577 0.5573 3.6664
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01149 0.1072
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.083714 0.174602 6.207 5.41e-10 ***
## land_useBedouin Agriculture -0.263209 0.234868 -1.121 0.2624
## land_useKKL Plantings 0.339773 0.160606 2.116 0.0344 *
## year_ct -0.009524 0.027935 -0.341 0.7332
## cos_td_rad 0.348352 0.154097 2.261 0.0238 *
## land_useBedouin Agriculture:year_ct 0.102878 0.043278 2.377 0.0174 *
## land_useKKL Plantings:year_ct 0.013666 0.035496 0.385 0.7002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ l_BA:_
## lnd_sBdnAgr -0.223
## lnd_sKKLPln -0.321 0.389
## year_ct -0.608 0.294 0.438
## cos_td_rad -0.720 -0.173 -0.191 0.203
## lnd_sBAgr:_ 0.196 -0.811 -0.332 -0.602 0.148
## lnd_sKKLP:_ 0.242 -0.265 -0.737 -0.693 0.147 0.494## $site
## (Intercept)
## Eshel Hanasi -0.045137191
## Givot Bar 0.002827638
## Givot Goral 0.106788644
## Hatzerim Base -0.001968742
## Mishmar Hanegev 0.075006970
## Nahal Ashan 0.065199853
## Nevatim Base 0.005947562
## Park Loess -0.093317558
## Sayeret Shaked -0.103992836
##
## with conditional variances for "site"## $site
Model with interaction of land use with time, to test for overall temporal trend in richness:
## Single term deletions
##
## Model:
## richness ~ land_use + year_ct + cos_td_rad + (1 | site)
## npar AIC
## <none> 589.43
## land_use 2 596.72
## year_ct 1 588.65
## cos_td_rad 1 591.67drop year.
## Single term deletions
##
## Model:
## richness ~ land_use + cos_td_rad + (1 | site)
## npar AIC
## <none> 588.65
## land_use 2 597.06
## cos_td_rad 1 589.67drop cosine.
## Single term deletions
##
## Model:
## richness ~ land_use + (1 | site)
## npar AIC
## <none> 589.67
## land_use 2 598.91After removing the interaction, overall temporal trend is not significant. only land use remains in the model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ land_use + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 589.7 601.4 -290.8 581.7 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8990 -0.6971 -0.1091 0.5584 3.9218
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01637 0.128
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.29743 0.09395 13.809 < 2e-16 ***
## land_useBedouin Agriculture 0.23985 0.13314 1.801 0.071628 .
## land_useKKL Plantings 0.39823 0.10824 3.679 0.000234 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA
## lnd_sBdnAgr -0.695
## lnd_sKKLPln -0.654 0.588## $site
## (Intercept)
## Eshel Hanasi -0.05918680
## Givot Bar 0.02700116
## Givot Goral 0.15312505
## Hatzerim Base 0.01695068
## Mishmar Hanegev 0.10556380
## Nahal Ashan 0.03871069
## Nevatim Base -0.01360889
## Park Loess -0.11623294
## Sayeret Shaked -0.13582493
##
## with conditional variances for "site"## $site
model interpretation
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 140 |
| Dependent variable | richness |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 586.971 |
| BIC | 610.504 |
| Pseudo-R² (fixed effects) | 0.187 |
| Pseudo-R² (total) | 0.231 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 2.956 | 0.175 | 6.207 | 0.000 |
| land_useBedouin Agriculture | 0.769 | 0.235 | -1.121 | 0.262 |
| land_useKKL Plantings | 1.405 | 0.161 | 2.116 | 0.034 |
| year_ct | 0.991 | 0.028 | -0.341 | 0.733 |
| cos_td_rad | 1.417 | 0.154 | 2.261 | 0.024 |
| land_useBedouin Agriculture:year_ct | 1.108 | 0.043 | 2.377 | 0.017 |
| land_useKKL Plantings:year_ct | 1.014 | 0.035 | 0.385 | 0.700 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.107 |
| Group | # groups | ICC |
|---|---|---|
| site | 9 | 0.011 |
## Warning: Removed 25 rows containing missing values (`geom_path()`).
Test whether temporal trends are significantly different from zero, using emmeans package, FDR correction for multiple comparisons.
## land_use year_ct.trend SE df z.ratio p.value
## Loess -0.00952 0.0279 Inf -0.341 0.8725
## Bedouin Agriculture 0.09335 0.0346 Inf 2.698 0.0209
## KKL Plantings 0.00414 0.0258 Inf 0.161 0.8725
##
## P value adjustment: fdr method for 3 testsTest for pairwise differences among levels of land use, apply FDR correction for multiple comparisons.
## contrast year_ct estimate SE df z.ratio p.value
## Loess - Bedouin Agriculture 4.21 -0.170 0.138 Inf -1.239 0.2154
## Loess - KKL Plantings 4.21 -0.397 0.113 Inf -3.521 0.0013
## Bedouin Agriculture - KKL Plantings 4.21 -0.227 0.116 Inf -1.959 0.0751
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testsNo significant overall temporal trend in richness in the Loess unit.
Significant increase of 9.7850508% in richness per year on average for Bedouin agriculture (total increase of 75.0891538% across six years 2014-2020), no significant temporal trend for Loess and KKL Plantings.
on average over time, Loess richness (3.589255 species) is significantly lower than KKL (5.340458 species; P=0.001); Bedouin agriculture (4.2558382 species) is not significantly different from both of them.
geometric mean of abundance
Explore data. Exclude time of day because of high number of NAs.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"## [1] "data with 2012"
## [1] "data without 2012"
## [1] "data with 2012"
## [1] "data without 2012"
| gma | year_ct | site | land_use | td_sc | cos_td_rad | sin_td_rad | h_from_sunrise | cos_hsun | sin_hsun | monitors_name | wind | precipitation | temperature | clouds | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 1.260 | Min. :0.000 | Sayeret Shaked :30 | Loess :51 | Min. :-1.5653 | Min. :0.02295 | Min. :-0.9997 | Length:140 | Min. :0.1951 | Min. :-0.1642 | Eyal Shochat : 30 | 0 : 0 | 0 : 0 | 0 : 0 | 0 : 0 | |
| 1st Qu.: 2.000 | 1st Qu.:2.000 | Nahal Ashan :24 | Bedouin Agriculture:44 | 1st Qu.:-0.7244 | 1st Qu.:0.34157 | 1st Qu.:-0.9201 | Class :difftime | 1st Qu.:0.8607 | 1st Qu.: 0.1331 | Adi Domer : 0 | 1 : 0 | 3 : 0 | 1 : 0 | 1 : 0 | |
| Median : 2.406 | Median :4.000 | Mishmar Hanegev:23 | KKL Plantings :45 | Median : 0.6826 | Median :0.85208 | Median :-0.5087 | Mode :numeric | Median :0.9542 | Median : 0.2990 | Asaf Mayrose : 0 | 2 : 0 | NA’s:140 | 2 : 0 | 2 : 0 | |
| Mean : 2.926 | Mean :4.214 | Givot Goral :15 | NA | Mean : 0.6387 | Mean :0.67279 | Mean :-0.4482 | NA | Mean :0.8921 | Mean : 0.3286 | Eliraz Dvir : 0 | 3 : 0 | NA | 3 : 0 | 3 : 0 | |
| 3rd Qu.: 3.207 | 3rd Qu.:6.000 | Eshel Hanasi :12 | NA | 3rd Qu.: 1.5220 | 3rd Qu.:0.90882 | 3rd Qu.:-0.1543 | NA | 3rd Qu.:0.9893 | 3rd Qu.: 0.5090 | Eliraz Dvir and Yoav Barak: 0 | NA’s:140 | NA | NA’s:140 | NA’s:140 | |
| Max. :10.371 | Max. :8.000 | Givot Bar :12 | NA | Max. : 4.9558 | Max. :1.00000 | Max. : 0.9845 | NA | Max. :1.0000 | Max. : 0.9808 | (Other) : 0 | NA | NA | NA | NA | |
| NA | NA | (Other) :24 | NA | NA | NA | NA | NA | NA’s :51 | NA’s :51 | NA’s :110 | NA | NA | NA | NA |
Fit glm, compare gamma, gaussian (poisson inappropriate because response is not discrete)
Gamma seems better than gaussian. Fit fixed and mixed models.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.012328 (tol = 0.002, component 1)Mixed model did not converge, use glm:
##
## Call:
## glm(formula = gma ~ land_use * year_ct + cos_td_rad + sin_td_rad +
## site, family = Gamma, data = P.anal)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.95553 -0.28198 -0.07873 0.15650 1.16860
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3840219 0.1001764 3.833 0.0002 ***
## land_useBedouin Agriculture -0.0842592 0.0824937 -1.021 0.3091
## land_useKKL Plantings 0.0001899 0.0583558 0.003 0.9974
## year_ct -0.0089596 0.0088011 -1.018 0.3107
## cos_td_rad -0.0206125 0.0613463 -0.336 0.7374
## sin_td_rad 0.0954936 0.0531921 1.795 0.0750 .
## siteGivot Bar -0.0780948 0.0546019 -1.430 0.1552
## siteGivot Goral 0.0075556 0.0627164 0.120 0.9043
## siteHatzerim Base -0.0337372 0.0918198 -0.367 0.7139
## siteMishmar Hanegev -0.0260323 0.0584789 -0.445 0.6570
## siteNahal Ashan 0.0940335 0.0822483 1.143 0.2551
## siteNevatim Base 0.0066576 0.1716479 0.039 0.9691
## sitePark Loess 0.0131366 0.0860731 0.153 0.8789
## siteSayeret Shaked 0.0386733 0.0779441 0.496 0.6207
## land_useBedouin Agriculture:year_ct 0.0276513 0.0126212 2.191 0.0303 *
## land_useKKL Plantings:year_ct 0.0285655 0.0124291 2.298 0.0232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1956952)
##
## Null deviance: 28.575 on 139 degrees of freedom
## Residual deviance: 21.703 on 124 degrees of freedom
## AIC: 445.98
##
## Number of Fisher Scoring iterations: 5perform stepwise model selection of Gaussian model.
## Start: AIC=445.98
## gma ~ land_use * year_ct + cos_td_rad + sin_td_rad + site
##
## Df Deviance AIC
## - site 8 23.307 438.18
## - cos_td_rad 1 21.725 444.09
## <none> 21.703 445.98
## - sin_td_rad 1 22.370 447.39
## - land_use:year_ct 2 23.147 449.36
##
## Step: AIC=440.23
## gma ~ land_use + year_ct + cos_td_rad + sin_td_rad + land_use:year_ct
##
## Df Deviance AIC
## - cos_td_rad 1 23.479 439.08
## <none> 23.307 440.23
## - sin_td_rad 1 24.148 442.36
## - land_use:year_ct 2 24.632 442.74
##
## Step: AIC=439.29
## gma ~ land_use + year_ct + sin_td_rad + land_use:year_ct
##
## Df Deviance AIC
## <none> 23.479 439.29
## - sin_td_rad 1 24.167 440.64
## - land_use:year_ct 2 24.952 442.45keep removing terms, because \(\Delta AIC<2\). remove sine.
## Single term deletions
##
## Model:
## gma ~ land_use * year_ct
## Df Deviance AIC
## <none> 24.167 441.45
## land_use:year_ct 2 25.657 444.51year X land use remain. This is the final model:
##
## Call:
## glm(formula = gma ~ land_use * year_ct, family = Gamma, data = P.anal)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9292 -0.3539 -0.1036 0.1449 1.2127
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.380616 0.036636 10.389 <2e-16 ***
## land_useBedouin Agriculture -0.137185 0.062042 -2.211 0.0287 *
## land_useKKL Plantings -0.013240 0.058533 -0.226 0.8214
## year_ct -0.013438 0.007054 -1.905 0.0589 .
## land_useBedouin Agriculture:year_ct 0.024062 0.011716 2.054 0.0419 *
## land_useKKL Plantings:year_ct 0.028566 0.012488 2.288 0.0237 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2111039)
##
## Null deviance: 28.575 on 139 degrees of freedom
## Residual deviance: 24.167 on 134 degrees of freedom
## AIC: 441.45
##
## Number of Fisher Scoring iterations: 5
| Observations | 140 |
| Dependent variable | gma |
| Type | Generalized linear model |
| Family | Gamma |
| Link | inverse |
| χ²(5) | 4.408 |
| Pseudo-R² (Cragg-Uhler) | 0.165 |
| Pseudo-R² (McFadden) | 0.054 |
| AIC | 441.450 |
| BIC | 462.041 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | 0.381 | 0.037 | 10.389 | 0.000 |
| land_useBedouin Agriculture | -0.137 | 0.062 | -2.211 | 0.029 |
| land_useKKL Plantings | -0.013 | 0.059 | -0.226 | 0.821 |
| year_ct | -0.013 | 0.007 | -1.905 | 0.059 |
| land_useBedouin Agriculture:year_ct | 0.024 | 0.012 | 2.054 | 0.042 |
| land_useKKL Plantings:year_ct | 0.029 | 0.012 | 2.288 | 0.024 |
| Standard errors: MLE |
According to the model, no significant temporal trend for Loess, but interaction of year and land use is significant. Test whether temporal trends are significantly different from zero, using emmeans package, FDR correction for multiple comparisons.
## land_use year_ct.trend SE df t.ratio p.value
## Loess -0.0134 0.00705 134 -1.905 0.1768
## Bedouin Agriculture 0.0106 0.00935 134 1.136 0.2581
## KKL Plantings 0.0151 0.01030 134 1.468 0.2166
##
## P value adjustment: fdr method for 3 testsNo significant temporal trend detected for all three land uses. In view of this, re-fit model without interaction term.
## Start: AIC=446.07
## gma ~ land_use + year_ct
##
## Df Deviance AIC
## - year_ct 1 25.657 444.07
## <none> 25.657 446.07
## - land_use 2 28.526 454.73
##
## Step: AIC=444.07
## gma ~ land_use
##
## Df Deviance AIC
## <none> 25.657 444.07
## - land_use 2 28.575 453.04only land use remains.
##
## Call:
## glm(formula = gma ~ land_use, family = Gamma, data = P.anal)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8448 -0.3748 -0.1082 0.1509 1.3829
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.32831 0.02181 15.056 <2e-16 ***
## land_useBedouin Agriculture -0.03117 0.03045 -1.024 0.3078
## land_useKKL Plantings 0.09527 0.03705 2.571 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2249982)
##
## Null deviance: 28.575 on 139 degrees of freedom
## Residual deviance: 25.657 on 137 degrees of freedom
## AIC: 444.07
##
## Number of Fisher Scoring iterations: 5
| Observations | 140 |
| Dependent variable | gma |
| Type | Generalized linear model |
| Family | Gamma |
| Link | inverse |
| χ²(2) | 2.918 |
| Pseudo-R² (Cragg-Uhler) | 0.110 |
| Pseudo-R² (McFadden) | 0.034 |
| AIC | 444.071 |
| BIC | 455.837 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | 0.328 | 0.022 | 15.056 | 0.000 |
| land_useBedouin Agriculture | -0.031 | 0.030 | -1.024 | 0.308 |
| land_useKKL Plantings | 0.095 | 0.037 | 2.571 | 0.011 |
| Standard errors: MLE |
Test for pairwise differences among levels of land use, apply FDR correction for multiple comparisons.
## contrast estimate SE df t.ratio p.value
## Loess - Bedouin Agriculture 0.0312 0.0304 137 1.024 0.3078
## Loess - KKL Plantings -0.0953 0.0370 137 -2.571 0.0168
## Bedouin Agriculture - KKL Plantings -0.1264 0.0367 137 -3.443 0.0023
##
## Note: contrasts are still on the inverse scale
## P value adjustment: fdr method for 3 testsNo significant difference in GMA between Loess (3.0458883) and Bedouin agriculture (3.3653697). Compared to KKL Plantings (2.3608274), natural Loess is 0.2901783% higher (P=0.017) and Bedouin agriculture is 0.4255044% higher (P=0.002).
No significant temporal trend detected overall in loess unit as well as for all three land uses separately.
abundance
Explore data
## [1] "ABUNDANCE WITHOUT RARE SPECIES"## [1] "data with 2012"
## [1] "data without 2012"
## [1] "data with 2012"
## [1] "data without 2012"
| abundance | year_ct | site | land_use | td_sc | cos_td_rad | sin_td_rad | |
|---|---|---|---|---|---|---|---|
| Min. : 2.00 | Min. :0.000 | Sayeret Shaked :30 | Loess :51 | Min. :-1.5653 | Min. :0.02295 | Min. :-0.9997 | |
| 1st Qu.: 7.00 | 1st Qu.:2.000 | Nahal Ashan :24 | Bedouin Agriculture:44 | 1st Qu.:-0.7244 | 1st Qu.:0.34157 | 1st Qu.:-0.9201 | |
| Median : 10.00 | Median :4.000 | Mishmar Hanegev:23 | KKL Plantings :45 | Median : 0.6826 | Median :0.85208 | Median :-0.5087 | |
| Mean : 16.34 | Mean :4.214 | Givot Goral :15 | NA | Mean : 0.6387 | Mean :0.67279 | Mean :-0.4482 | |
| 3rd Qu.: 17.25 | 3rd Qu.:6.000 | Eshel Hanasi :12 | NA | 3rd Qu.: 1.5220 | 3rd Qu.:0.90882 | 3rd Qu.:-0.1543 | |
| Max. :125.00 | Max. :8.000 | Givot Bar :12 | NA | Max. : 4.9558 | Max. :1.00000 | Max. : 0.9845 | |
| NA | NA | (Other) :24 | NA | NA | NA | NA |
outliers in loess with total abundance >30, in bedouin agriculture with total abundance >120, in KKL with abundance>60. Examine:
## unit subunit site year year_ct
## 1: Loess Covered Areas in the Northern Negev <NA> Nahal Ashan 2020 8
## settlements agriculture habitat dunes land_use point_name
## 1: <NA> <NA> <NA> <NA> Loess Nahal Ashan Loess 2
## date datetime td_sc td_rad cos_td_rad
## 1: 2020-03-25 02:00:00 2020-03-25 09:25:40 -1.48902 -1.513416 0.0573492
## sin_td_rad timediff_Jun21 monitors_name wind precipitation temperature
## 1: -0.9983542 -87.91667 days Eyal Shochat <NA> <NA> <NA>
## clouds h_from_sunrise cos_hsun sin_hsun pilot richness abundance gma
## 1: <NA> 1.7 hours 0.9025853 0.4305111 FALSE 5 38 4.967254## unit subunit site year year_ct
## 1: Loess Covered Areas in the Northern Negev <NA> Givot Goral 2012 0
## settlements agriculture habitat dunes land_use point_name
## 1: <NA> <NA> <NA> <NA> KKL Plantings Goral KKL Plantings 3
## date datetime td_sc td_rad cos_td_rad sin_td_rad timediff_Jun21
## 1: 2012-06-12 <NA> 1.521963 -0.1549279 0.9880227 -0.1543088 -9 days
## monitors_name wind precipitation temperature clouds h_from_sunrise cos_hsun
## 1: <NA> <NA> <NA> <NA> <NA> NA hours NA
## sin_hsun pilot richness abundance gma
## 1: NA TRUE 7 66 2.830611## unit subunit site year year_ct
## 1: Loess Covered Areas in the Northern Negev <NA> Givot Bar 2018 6
## settlements agriculture habitat dunes land_use
## 1: <NA> <NA> <NA> <NA> Bedouin Agriculture
## point_name date datetime td_sc
## 1: Givot Bar Bedouin Agriculture 2 2018-05-22 2018-05-22 05:21:00 0.6825762
## td_rad cos_td_rad sin_td_rad timediff_Jun21 monitors_name wind
## 1: -0.5336404 0.860961 -0.5086709 -31 days <NA> <NA>
## precipitation temperature clouds h_from_sunrise cos_hsun sin_hsun pilot
## 1: <NA> <NA> <NA> -0.42 hours 0.993961 -0.1097343 FALSE
## richness abundance gma
## 1: 7 125 5.011216Exclude all three plots.
## [1] "Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model."## [1] "od = 5.38981891399102"
PHI>1, hence choose negative binomial. Fit fixed and mixed models. Choose mixed model if possible, otherwise choose a model with fixed-effects only.
mixed model converged.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(4.1441) ( log )
## Formula: abundance ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 935.2 964.4 -457.6 915.2 127
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6521 -0.6367 -0.1579 0.3512 4.0056
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.05731 0.2394
## Number of obs: 137, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.591604 0.233268 6.823 8.91e-12 ***
## land_useBedouin Agriculture 0.417051 0.297780 1.401 0.161353
## land_useKKL Plantings 0.725205 0.210081 3.452 0.000556 ***
## year_ct 0.008704 0.036647 0.238 0.812257
## cos_td_rad 0.641133 0.201439 3.183 0.001459 **
## sin_td_rad -0.391850 0.168180 -2.330 0.019809 *
## land_useBedouin Agriculture:year_ct 0.054927 0.050512 1.087 0.276857
## land_useKKL Plantings:year_ct -0.093647 0.044448 -2.107 0.035128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ sn_td_ l_BA:_
## lnd_sBdnAgr -0.301
## lnd_sKKLPln -0.275 0.388
## year_ct -0.471 0.281 0.441
## cos_td_rad -0.694 -0.118 -0.200 0.086
## sin_td_rad 0.104 0.073 0.155 0.399 -0.255
## lnd_sBAgr:_ 0.199 -0.744 -0.314 -0.534 0.126 0.042
## lnd_sKKLP:_ 0.165 -0.210 -0.747 -0.597 0.178 -0.057 0.450## $site
## (Intercept)
## Eshel Hanasi -0.18031807
## Givot Bar 0.31083900
## Givot Goral -0.02950704
## Hatzerim Base -0.04610457
## Mishmar Hanegev 0.34287916
## Nahal Ashan -0.07148836
## Nevatim Base 0.08196214
## Park Loess -0.13717054
## Sayeret Shaked -0.25752018
##
## with conditional variances for "site"## $site
Perform stepwise model selection of mixed model.
## Single term deletions
##
## Model:
## abundance ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 935.23
## cos_td_rad 1 943.18
## sin_td_rad 1 938.70
## land_use:year_ct 2 940.97land use * year and sampling time of year remain. The final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(4.1441) ( log )
## Formula: abundance ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 935.2 964.4 -457.6 915.2 127
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6521 -0.6367 -0.1579 0.3512 4.0056
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.05731 0.2394
## Number of obs: 137, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.591604 0.233268 6.823 8.91e-12 ***
## land_useBedouin Agriculture 0.417051 0.297780 1.401 0.161353
## land_useKKL Plantings 0.725205 0.210081 3.452 0.000556 ***
## year_ct 0.008704 0.036647 0.238 0.812257
## cos_td_rad 0.641133 0.201439 3.183 0.001459 **
## sin_td_rad -0.391850 0.168180 -2.330 0.019809 *
## land_useBedouin Agriculture:year_ct 0.054927 0.050512 1.087 0.276857
## land_useKKL Plantings:year_ct -0.093647 0.044448 -2.107 0.035128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ sn_td_ l_BA:_
## lnd_sBdnAgr -0.301
## lnd_sKKLPln -0.275 0.388
## year_ct -0.471 0.281 0.441
## cos_td_rad -0.694 -0.118 -0.200 0.086
## sin_td_rad 0.104 0.073 0.155 0.399 -0.255
## lnd_sBAgr:_ 0.199 -0.744 -0.314 -0.534 0.126 0.042
## lnd_sKKLP:_ 0.165 -0.210 -0.747 -0.597 0.178 -0.057 0.450## $site
## (Intercept)
## Eshel Hanasi -0.18031807
## Givot Bar 0.31083900
## Givot Goral -0.02950704
## Hatzerim Base -0.04610457
## Mishmar Hanegev 0.34287916
## Nahal Ashan -0.07148836
## Nevatim Base 0.08196214
## Park Loess -0.13717054
## Sayeret Shaked -0.25752018
##
## with conditional variances for "site"## $site
Not the greatest fit in the world. Take into account when reporting the results. It seems that temporal trends in the different land uses are weak and probably not significant. Interaction of year and KKL plantings is marginally significant, but not sure if it is different from zero. Test for difference of temporal trend from zero for each land use.
## land_use year_ct.trend SE df z.ratio p.value
## Loess 0.0087 0.0366 Inf 0.238 0.8123
## Bedouin Agriculture 0.0636 0.0438 Inf 1.453 0.2193
## KKL Plantings -0.0849 0.0371 Inf -2.292 0.0658
##
## P value adjustment: fdr method for 3 testsAll three temporal trends are not significantly different from zero. Re-fit model without interaction of time with subunit
## Single term deletions
##
## Model:
## abundance ~ land_use + year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 940.53
## land_use 2 949.75
## year_ct 1 938.71
## cos_td_rad 1 949.89
## sin_td_rad 1 944.78drop year.
## Single term deletions
##
## Model:
## abundance ~ land_use + cos_td_rad + sin_td_rad + (1 | site)
## npar AIC
## <none> 938.71
## land_use 2 947.74
## cos_td_rad 1 949.88
## sin_td_rad 1 943.87Final model includes land use and sampling time of year.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(3.7354) ( log )
## Formula: abundance ~ land_use + cos_td_rad + sin_td_rad + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 938.7 959.2 -462.4 924.7 130
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4920 -0.5606 -0.1735 0.2896 4.7000
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0414 0.2035
## Number of obs: 137, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5549 0.2002 7.766 8.12e-15 ***
## land_useBedouin Agriculture 0.7848 0.1958 4.009 6.10e-05 ***
## land_useKKL Plantings 0.3965 0.1429 2.775 0.005526 **
## cos_td_rad 0.7135 0.1938 3.682 0.000232 ***
## sin_td_rad -0.3887 0.1466 -2.651 0.008032 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA l_KKLP cs_td_
## lnd_sBdnAgr -0.343
## lnd_sKKLPln -0.260 0.579
## cos_td_rad -0.756 -0.046 -0.135
## sin_td_rad 0.443 0.340 0.204 -0.469## $site
## (Intercept)
## Eshel Hanasi -0.16462507
## Givot Bar 0.21226536
## Givot Goral 0.03544298
## Hatzerim Base 0.01315824
## Mishmar Hanegev 0.26595661
## Nahal Ashan -0.08509449
## Nevatim Base 0.06613307
## Park Loess -0.10126718
## Sayeret Shaked -0.23221676
##
## with conditional variances for "site"## $site
Still not the greatest fit in the world. Take into account when reporting the results. Interpretation of abundance model:
| Observations | 137 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | Negative Binomial(3.7354) |
| Link | log |
| AIC | 938.711 |
| BIC | 959.151 |
| Pseudo-R² (fixed effects) | 0.612 |
| Pseudo-R² (total) | 0.750 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 1.555 | 0.200 | 7.766 | 0.000 |
| land_useBedouin Agriculture | 0.785 | 0.196 | 4.009 | 0.000 |
| land_useKKL Plantings | 0.397 | 0.143 | 2.775 | 0.006 |
| cos_td_rad | 0.713 | 0.194 | 3.682 | 0.000 |
| sin_td_rad | -0.389 | 0.147 | -2.651 | 0.008 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.203 |
| Group | # groups | ICC |
|---|---|---|
| site | 9 | 0.080 |
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
Test for pairwise differences in average abundance among levels of land use, apply FDR correction for multiple comparisons.
## contrast estimate SE df z.ratio p.value
## Loess - Bedouin Agriculture -0.785 0.196 Inf -4.009 0.0002
## Loess - KKL Plantings -0.397 0.143 Inf -2.775 0.0083
## Bedouin Agriculture - KKL Plantings 0.388 0.162 Inf 2.392 0.0168
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testsSignificant difference in average abundance among all three land uses - Loess (9.1048207), Bedouin agriculture (19.9565935) and KKL Plantings (13.5357733).
Total abundance of a KKL plantings plot is, on average, 48.6660062% higher compared to a natural loess plot; Bedouin agriculture is 119.187112% higher compared to natural loess.
community analysis using package MVabund
Explore single species observations, including mean-variance plot. There is a strong mean-var relationship, indicating that employing GLMs is the proper way to analyze, rather than OLS (assumption of homogeneity is violated).
## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Galerida.cristata, Passer.domesticus, Columba.livia, Streptopelia.decaocto, Prinia.gracilis, Lanius.excubitor, Corvus.cornix, Burhinus.oedicnemus, Chloris.chloris, Spilopelia.senegalensis, Acridotheres.tristis, Falco.tinnunculus were included in the plot
## (the variables with highest total abundance).
## [1] "abundance observations, aggregated into plots and species (i.e., several observations from the same species in the same plot are aggregated):"
## [1] "these are the actual observations, before aggregation into plots:"
## [1] "zoom in on high abundance observations:"
start model specification:
## nb po
## 227.5317 318.7130## [1] "POISSON"
## [1] "NEGATIVE BINOMIAL"
negative binomial model is better than poisson according to residuals and AIC comparison.
## nb po nb1
## 227.5317 318.7130 227.6228The addition of the explanatory variable ‘site’ is not improving the AIC of the model. stepwise selection of model:
## Single term deletions
##
## Model:
## spp_no_rare ~ land_use * year_ct + cos_td_rad + sin_td_rad
## Df AIC
## <none> 3868.0
## cos_td_rad 17 3913.9
## sin_td_rad 17 3871.3
## land_use:year_ct 34 3892.2close call for sine. when included in the final model comes out borderline, P=0.05. remove.
## Single term deletions
##
## Model:
## spp_no_rare ~ land_use * year_ct + cos_td_rad
## Df AIC
## <none> 3871.3
## cos_td_rad 17 3897.0
## land_use:year_ct 34 3891.6final model includes land use * year, sampling time of year.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 11.318 0.001 ***
## land_useBedouin Agriculture 5.697 0.001 ***
## land_useKKL Plantings 8.128 0.001 ***
## year_ct 6.616 0.001 ***
## cos_td_rad 7.148 0.001 ***
## land_useBedouin Agriculture:year_ct 3.231 0.365
## land_useKKL Plantings:year_ct 6.670 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 16.33, 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).## Analysis of Deviance Table
##
## Model: spp_no_rare ~ land_use * year_ct + cos_td_rad
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 139
## land_use 137 2 246.61 0.01 **
## year_ct 136 1 48.13 0.01 **
## cos_td_rad 135 1 47.55 0.01 **
## land_use:year_ct 133 2 88.35 0.01 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acridotheres.tristis Burhinus.oedicnemus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 5.752 0.33 6.596 0.23
## year_ct 12.957 0.01 0.219 1.00
## cos_td_rad 0 1.00 15.536 0.02
## land_use:year_ct 2.76 0.90 14.101 0.03
## Chloris.chloris Columba.livia Corvus.cornix
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## land_use 4.174 0.44 13.207 0.03 5.238
## year_ct 0.414 1.00 0.354 1.00 0.331
## cos_td_rad 0.219 1.00 0.078 1.00 0.113
## land_use:year_ct 2.824 0.90 15.761 0.02 1.046
## Falco.tinnunculus Galerida.cristata
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 0.33 1.984 0.72 12.595 0.03
## year_ct 1.00 3.545 0.69 2.399 0.87
## cos_td_rad 1.00 3.665 0.60 4.986 0.39
## land_use:year_ct 0.97 0.695 0.97 3.698 0.76
## Iduna.pallida Lanius.excubitor
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 27.694 0.01 34.644 0.01
## year_ct 0.014 1.00 0.015 1.00
## cos_td_rad 1.292 0.98 0.283 1.00
## land_use:year_ct 0 0.97 3.911 0.71
## Passer.domesticus Prinia.gracilis
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 29.549 0.01 1.698 0.72
## year_ct 2.472 0.87 1.002 0.99
## cos_td_rad 3.581 0.62 0.35 1.00
## land_use:year_ct 17.406 0.02 4.83 0.62
## Pycnonotus.xanthopygos Spilopelia.senegalensis
## Dev Pr(>Dev) Dev
## (Intercept)
## land_use 13.186 0.03 20.893
## year_ct 0.091 1.00 0.113
## cos_td_rad 0.619 1.00 3.265
## land_use:year_ct 7.731 0.26 0.214
## Streptopelia.decaocto Streptopelia.turtur
## Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## land_use 0.01 40.218 0.01 0.511
## year_ct 1.00 13.932 0.01 6.625
## cos_td_rad 0.69 0.583 1.00 5.012
## land_use:year_ct 0.97 2.149 0.90 5.806
## Upupa.epops Vanellus.spinosus
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 0.84 10.642 0.03 18.028 0.01
## year_ct 0.17 1.411 0.97 2.235 0.87
## cos_td_rad 0.39 0.279 1.00 7.688 0.18
## land_use:year_ct 0.54 1.122 0.96 4.293 0.68
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 99 iterations via PIT-trap resampling.All factors (land_use, year, land_use:year, time of year) have a statistically significant effect on community composition.
The interaction land use * year is significantly different from Loess only for KKL plantings, not for Bedouin agriculture. In addition, this term is only significant for 3 species: stone curlew, rock pigeon and house sparrow. Therefore, plot species - level model for these 3 species which includes the interaction, and re-run manyglm without interaction to simplify interpretation.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 11.726 0.001 ***
## land_useBedouin Agriculture 9.586 0.001 ***
## land_useKKL Plantings 10.209 0.001 ***
## year_ct 7.001 0.001 ***
## cos_td_rad 6.012 0.024 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 15.03, 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).## Analysis of Deviance Table
##
## Model: spp_no_rare ~ land_use + year_ct + cos_td_rad
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 139
## land_use 137 2 246.61 0.01 **
## year_ct 136 1 48.13 0.01 **
## cos_td_rad 135 1 47.55 0.01 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acridotheres.tristis Burhinus.oedicnemus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 5.752 0.32 6.596 0.22
## year_ct 12.957 0.01 0.219 1.00
## cos_td_rad 0 1.00 15.536 0.01
## Chloris.chloris Columba.livia Corvus.cornix
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## land_use 4.174 0.40 13.207 0.03 5.238
## year_ct 0.414 0.99 0.354 1.00 0.331
## cos_td_rad 0.219 1.00 0.078 1.00 0.113
## Falco.tinnunculus Galerida.cristata
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 0.33 1.984 0.70 12.595 0.03
## year_ct 1.00 3.545 0.63 2.399 0.83
## cos_td_rad 1.00 3.665 0.76 4.986 0.50
## Iduna.pallida Lanius.excubitor Passer.domesticus
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## land_use 27.694 0.01 34.644 0.01 29.549
## year_ct 0.014 1.00 0.015 1.00 2.472
## cos_td_rad 1.292 0.98 0.283 1.00 3.581
## Prinia.gracilis Pycnonotus.xanthopygos
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 0.01 1.698 0.70 13.186 0.03
## year_ct 0.83 1.002 0.98 0.091 1.00
## cos_td_rad 0.76 0.35 1.00 0.619 1.00
## Spilopelia.senegalensis Streptopelia.decaocto
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## land_use 20.893 0.01 40.218 0.01
## year_ct 0.113 1.00 13.932 0.01
## cos_td_rad 3.265 0.76 0.583 1.00
## Streptopelia.turtur Upupa.epops Vanellus.spinosus
## Dev Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## land_use 0.511 0.79 10.642 0.03 18.028
## year_ct 6.625 0.18 1.411 0.94 2.235
## cos_td_rad 5.012 0.49 0.279 1.00 7.688
##
## Pr(>Dev)
## (Intercept)
## land_use 0.01
## year_ct 0.83
## cos_td_rad 0.16
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 99 iterations via PIT-trap resampling.## 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.## [[1]]
##
## [[2]]
##
## [[3]]
## [[1]]
##
## [[2]]
Significant negative temporal effect on house sparrow and positive effect on common myna. Significant effect land use: most species do not prefer natural loess over the other two habitats. Nonetheless, batha specialist species are not detected in our survey: is it possible they have disappeared?
check general trend in species that were filtered due to rareness in dataset.
this text refers to endangered species being vulnerable and up. Analysis was updated to include NT category. Of the 24 species that did not make the rare species cutoff, 9 (38%) appeared only in the natural loess vs. only 3 (12%) in KKL plantings and 2 (8%) in the Bedouin agriculture. Of these 9 species, 9 are defined as batha specialists, vs. 1 in KKL plantings and 0 in Bedouin agriculture, respectively. Moreover, of the 9 species unique to the natural loess, 7 are endangered vs. 1 in KKL plantings and 0 in Bedouin agriculture, respectively.
Test whether there is a temporal trend within batha specialists, without filtering rare species
## [1] "RICHNESS OF ONLY BATHA SPECIALISTS, WITH RARE SPECIES"
## [1] "ABUNDANCE OF ONLY BATHA SPECIALISTS, WITH RARE SPECIES"
## [1] "Estimating overdispersion parameter phi: (Res. Dev.)/(n-p) where n=number of observations; p=number of parameters in the model."## [1] "od = 1.1336296386876"
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
Prefer poisson because phi ~ 1 and neg bin model did not converge.
mixed model converged.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: abundance ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 627.9 654.3 -304.9 609.9 131
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8460 -0.7303 -0.0729 0.5858 3.8127
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01736 0.1318
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.60431 0.15943 10.063 < 2e-16 ***
## land_useBedouin Agriculture -0.19108 0.22635 -0.844 0.399
## land_useKKL Plantings -0.75252 0.17780 -4.232 2.31e-05 ***
## year_ct -0.00278 0.02513 -0.111 0.912
## cos_td_rad 0.21577 0.14848 1.453 0.146
## sin_td_rad -0.15301 0.13005 -1.177 0.239
## land_useBedouin Agriculture:year_ct -0.04693 0.04253 -1.103 0.270
## land_useKKL Plantings:year_ct 0.04248 0.03604 1.179 0.238
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ sn_td_ l_BA:_
## lnd_sBdnAgr -0.202
## lnd_sKKLPln -0.174 0.243
## year_ct -0.495 0.210 0.330
## cos_td_rad -0.751 -0.114 -0.169 0.151
## sin_td_rad 0.153 -0.044 0.073 0.368 -0.240
## lnd_sBAgr:_ 0.184 -0.816 -0.207 -0.413 0.053 0.147
## lnd_sKKLP:_ 0.144 -0.169 -0.791 -0.479 0.122 -0.005 0.311## $site
## (Intercept)
## Eshel Hanasi -0.102802228
## Givot Bar 0.042769996
## Givot Goral 0.053940917
## Hatzerim Base -0.022533918
## Mishmar Hanegev 0.067533671
## Nahal Ashan -0.192021477
## Nevatim Base 0.148233206
## Park Loess 0.019603381
## Sayeret Shaked 0.003242354
##
## with conditional variances for "site"## $site
model selection
## Single term deletions
##
## Model:
## abundance ~ land_use * year_ct + cos_td_rad + sin_td_rad + (1 |
## site)
## npar AIC
## <none> 627.85
## cos_td_rad 1 628.00
## sin_td_rad 1 627.34
## land_use:year_ct 2 627.65drop sine
## Single term deletions
##
## Model:
## abundance ~ land_use * year_ct + cos_td_rad + (1 | site)
## npar AIC
## <none> 627.34
## cos_td_rad 1 626.94
## land_use:year_ct 2 626.63drop interaction land use * year
## Single term deletions
##
## Model:
## abundance ~ land_use + year_ct + cos_td_rad + (1 | site)
## npar AIC
## <none> 626.63
## land_use 2 652.65
## year_ct 1 625.14
## cos_td_rad 1 626.18drop year.
## Single term deletions
##
## Model:
## abundance ~ land_use + cos_td_rad + (1 | site)
## npar AIC
## <none> 625.14
## land_use 2 650.65
## cos_td_rad 1 624.21drop cosine.
## Single term deletions
##
## Model:
## abundance ~ land_use + (1 | site)
## npar AIC
## <none> 624.21
## land_use 2 649.15final model includes only land use:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: abundance ~ land_use + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 624.2 636.0 -308.1 616.2 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8662 -0.8424 -0.0583 0.6095 3.8781
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01181 0.1087
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.77746 0.07508 23.674 < 2e-16 ***
## land_useBedouin Agriculture -0.39532 0.11637 -3.397 0.000681 ***
## land_useKKL Plantings -0.55276 0.10501 -5.264 1.41e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA
## lnd_sBdnAgr -0.612
## lnd_sKKLPln -0.487 0.375## $site
## (Intercept)
## Eshel Hanasi -0.08274936
## Givot Bar 0.05382039
## Givot Goral 0.03473399
## Hatzerim Base -0.05185206
## Mishmar Hanegev 0.07160145
## Nahal Ashan -0.14196387
## Nevatim Base 0.06076581
## Park Loess 0.04465879
## Sayeret Shaked 0.02313821
##
## with conditional variances for "site"## $site
| Observations | 140 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 624.210 |
| BIC | 635.976 |
| Pseudo-R² (fixed effects) | 0.212 |
| Pseudo-R² (total) | 0.257 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 1.777 | 0.075 | 23.674 | 0.000 |
| land_useBedouin Agriculture | -0.395 | 0.116 | -3.397 | 0.001 |
| land_useKKL Plantings | -0.553 | 0.105 | -5.264 | 0.000 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.109 |
| Group | # groups | ICC |
|---|---|---|
| site | 9 | 0.012 |
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
plot observations on ternary plot
## land_use year number_of_sites number_of_plots
## 1: Bedouin Agriculture 2014 3 9
## 2: Bedouin Agriculture 2016 4 11
## 3: Bedouin Agriculture 2018 4 12
## 4: Bedouin Agriculture 2020 4 12
## 5: KKL Plantings 2012 3 9
## 6: KKL Plantings 2014 3 9
## 7: KKL Plantings 2016 3 9
## 8: KKL Plantings 2018 3 9
## 9: KKL Plantings 2020 3 9
## 10: Loess 2012 4 12
## 11: Loess 2014 3 9
## 12: Loess 2016 3 9
## 13: Loess 2018 3 9
## 14: Loess 2020 3 9randomly choose 9 plots from bedouin agriculture for each year in the range 2016-2020
## land_use tot_plots
## 1: Bedouin Agriculture 44
## 2: KKL Plantings 36
## 3: Loess 36## land_use SciName abundance
## 1: Bedouin Agriculture Acridotheres tristis 18
## 2: Bedouin Agriculture Alectoris chukar 18
## 3: Bedouin Agriculture Anthus cervinus 0
## 4: Bedouin Agriculture Anthus trivialis 0
## 5: Bedouin Agriculture Apus apus 0
## ---
## 150: Loess Pterocles alchata 11
## 151: Loess Pterocles orientalis 0
## 152: Loess Streptopelia decaocto 19
## 153: Loess Streptopelia turtur 2
## 154: Loess Upupa epops 3Test differences among land uses for abundance of batha specialists:
## contrast estimate SE df z.ratio p.value
## Loess - Bedouin Agriculture 0.395 0.116 Inf 3.397 0.0010
## Loess - KKL Plantings 0.553 0.105 Inf 5.264 <.0001
## Bedouin Agriculture - KKL Plantings 0.157 0.124 Inf 1.269 0.2045
##
## Results are given on the log (not the response) scale.
## P value adjustment: fdr method for 3 testsSignificant difference in total abundance of BATHA SPECIALISTS per plot among all three land uses - Loess (5.9148344), Bedouin agriculture (3.9834139) and KKL Plantings (3.4031412).
Total abundance of a Bedouin agriculture plot is, on average, 17.0510912% higher compared to a KKL plantings plot (P=0.205); natural loess is 73.8051407% higher compared to KKL plantings (P<0.001); natural loess is 48.4865616% higher compared to Bedouin agriculture (P=0.001).
No general temporal trend found, nor such a separate trend for the different land uses
Test whether there is a temporal trend for the proportion of endangered species, without filtering rare species
## [1] "PROPORTION OF ENDANGERED SPECIES, INCLUDING RARE SPECIES"
## boundary (singular) fit: see help('isSingular')mixed model did not converge.
##
## Call:
## glm(formula = prop_rich_endang ~ land_use * year_ct + cos_td_rad +
## sin_td_rad + site, family = binomial, data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5476 -0.5967 -0.3685 -0.0554 1.7772
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3931 1.8927 -1.264 0.20609
## land_useBedouin Agriculture 0.2334 2.1701 0.108 0.91435
## land_useKKL Plantings -1.7451 0.5466 -3.193 0.00141 **
## year_ct -0.3958 0.1360 -2.911 0.00360 **
## cos_td_rad 1.1976 0.9067 1.321 0.18655
## sin_td_rad -0.5972 0.7853 -0.760 0.44699
## siteGivot Bar -0.6794 1.5294 -0.444 0.65688
## siteGivot Goral 0.3524 1.6227 0.217 0.82807
## siteHatzerim Base -0.1831 1.7286 -0.106 0.91562
## siteMishmar Hanegev -0.7806 1.5539 -0.502 0.61544
## siteNahal Ashan 0.6022 1.6517 0.365 0.71542
## siteNevatim Base 1.9235 2.2099 0.870 0.38409
## sitePark Loess 1.3679 1.7399 0.786 0.43174
## siteSayeret Shaked 0.3648 1.6497 0.221 0.82498
## land_useBedouin Agriculture:year_ct -0.2871 0.4451 -0.645 0.51883
## land_useKKL Plantings:year_ct 0.3594 0.1546 2.325 0.02009 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 131.854 on 139 degrees of freedom
## Residual deviance: 72.297 on 124 degrees of freedom
## AIC: 171.15
##
## Number of Fisher Scoring iterations: 7
##
## Call:
## glm(formula = prop_rich_endang ~ land_use * year_ct + cos_td_rad +
## sin_td_rad + site, family = quasibinomial, data = P.anal,
## weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5476 -0.5967 -0.3685 -0.0554 1.7772
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3931 1.5609 -1.533 0.127796
## land_useBedouin Agriculture 0.2334 1.7898 0.130 0.896448
## land_useKKL Plantings -1.7451 0.4508 -3.871 0.000174 ***
## year_ct -0.3958 0.1121 -3.530 0.000585 ***
## cos_td_rad 1.1976 0.7478 1.602 0.111796
## sin_td_rad -0.5972 0.6476 -0.922 0.358292
## siteGivot Bar -0.6794 1.2613 -0.539 0.591103
## siteGivot Goral 0.3524 1.3382 0.263 0.792736
## siteHatzerim Base -0.1831 1.4256 -0.128 0.897990
## siteMishmar Hanegev -0.7806 1.2815 -0.609 0.543584
## siteNahal Ashan 0.6022 1.3622 0.442 0.659209
## siteNevatim Base 1.9235 1.8225 1.055 0.293309
## sitePark Loess 1.3679 1.4349 0.953 0.342288
## siteSayeret Shaked 0.3648 1.3606 0.268 0.789038
## land_useBedouin Agriculture:year_ct -0.2871 0.3671 -0.782 0.435553
## land_useKKL Plantings:year_ct 0.3594 0.1275 2.819 0.005614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.6801662)
##
## Null deviance: 131.854 on 139 degrees of freedom
## Residual deviance: 72.297 on 124 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7
dispersion parameter is different than 1, hence use quasibinomial. BUT: AIC undefined for quasibinomial, use binomial.
## Start: AIC=171.15
## prop_rich_endang ~ land_use * year_ct + cos_td_rad + sin_td_rad +
## site
##
## Df Deviance AIC
## - site 8 81.158 164.01
## - sin_td_rad 1 72.919 169.77
## - cos_td_rad 1 74.240 171.09
## <none> 72.297 171.15
## - land_use:year_ct 2 79.158 174.01
##
## Step: AIC=164.01
## prop_rich_endang ~ land_use + year_ct + cos_td_rad + sin_td_rad +
## land_use:year_ct
##
## Df Deviance AIC
## - cos_td_rad 1 81.752 162.61
## - sin_td_rad 1 82.318 163.17
## <none> 81.158 164.01
## - land_use:year_ct 2 85.557 164.41
##
## Step: AIC=162.61
## prop_rich_endang ~ land_use + year_ct + sin_td_rad + land_use:year_ct
##
## Df Deviance AIC
## - sin_td_rad 1 83.429 162.28
## <none> 81.752 162.61
## - land_use:year_ct 2 86.079 162.93
##
## Step: AIC=162.28
## prop_rich_endang ~ land_use + year_ct + land_use:year_ct
##
## Df Deviance AIC
## - land_use:year_ct 2 86.951 161.80
## <none> 83.429 162.28
##
## Step: AIC=161.8
## prop_rich_endang ~ land_use + year_ct
##
## Df Deviance AIC
## <none> 86.951 161.80
## - year_ct 1 99.771 172.62
## - land_use 2 105.881 176.74final model includes no interaction, only year and land use.
##
## Call:
## glm(formula = prop_rich_endang ~ land_use + year_ct, family = binomial,
## data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.62244 -0.63499 -0.39624 0.01404 2.64293
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.11484 0.24868 -4.483 7.36e-06 ***
## land_useBedouin Agriculture -2.06201 0.63969 -3.223 0.001267 **
## land_useKKL Plantings -1.12434 0.35648 -3.154 0.001611 **
## year_ct -0.21354 0.06359 -3.358 0.000786 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 131.854 on 139 degrees of freedom
## Residual deviance: 86.951 on 136 degrees of freedom
## AIC: 161.8
##
## Number of Fisher Scoring iterations: 6
remove observation 100, re-fit
##
## Call:
## glm(formula = prop_rich_endang ~ land_use + year_ct, family = binomial,
## data = P.anal2, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5664 -0.5511 -0.3595 -0.0645 1.8973
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.02459 0.24951 -4.106 4.02e-05 ***
## land_useBedouin Agriculture -1.94585 0.64556 -3.014 0.002577 **
## land_useKKL Plantings -1.36528 0.38922 -3.508 0.000452 ***
## year_ct -0.25916 0.07026 -3.688 0.000226 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 128.974 on 138 degrees of freedom
## Residual deviance: 78.923 on 135 degrees of freedom
## AIC: 151.01
##
## Number of Fisher Scoring iterations: 6
Test for pairwise differences among levels of land use, apply FDR correction for multiple comparisons.
## contrast year_ct estimate SE df z.ratio p.value
## Loess - Bedouin Agriculture 4.2 1.946 0.646 Inf 3.014 0.0039
## Loess - KKL Plantings 4.2 1.365 0.389 Inf 3.508 0.0014
## Bedouin Agriculture - KKL Plantings 4.2 -0.581 0.691 Inf -0.840 0.4007
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: fdr method for 3 testsTest whether there is a temporal trend for the proportion of batha specialist species, without filtering rare species
## [1] "PROPORTION OF BATHA SPECIALIST SPECIES, INCLUDING RARE SPECIES"
all models converged.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: prop_rich_batha ~ land_use * year_ct + cos_td_rad + sin_td_rad +
## (1 | site)
## Data: P.anal
## Weights: richness
##
## AIC BIC logLik deviance df.resid
## 314.2 340.7 -148.1 296.2 131
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0495 -0.3775 0.1247 0.5130 2.3382
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.004138 0.06433
## Number of obs: 140, groups: site, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22685 0.35595 0.637 0.524
## land_useBedouin Agriculture -0.57734 0.52504 -1.100 0.272
## land_useKKL Plantings -1.49484 0.37966 -3.937 8.24e-05 ***
## year_ct -0.10910 0.06727 -1.622 0.105
## cos_td_rad -0.44633 0.34923 -1.278 0.201
## sin_td_rad -0.06857 0.28534 -0.240 0.810
## land_useBedouin Agriculture:year_ct -0.02940 0.09720 -0.302 0.762
## land_useKKL Plantings:year_ct 0.11181 0.07964 1.404 0.160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lnd_BA ln_KKLP yer_ct cs_td_ sn_td_ l_BA:_
## lnd_sBdnAgr -0.100
## lnd_sKKLPln -0.232 0.299
## year_ct -0.546 0.201 0.401
## cos_td_rad -0.763 -0.229 -0.208 0.159
## sin_td_rad -0.023 0.006 0.189 0.572 -0.100
## lnd_sBAgr:_ 0.132 -0.871 -0.267 -0.419 0.162 -0.034
## lnd_sKKLP:_ 0.224 -0.226 -0.782 -0.580 0.119 -0.168 0.371## $site
## (Intercept)
## Eshel Hanasi 0.006934564
## Givot Bar -0.001013631
## Givot Goral 0.001136048
## Hatzerim Base -0.008359415
## Mishmar Hanegev -0.007469215
## Nahal Ashan -0.030336018
## Nevatim Base 0.003267121
## Park Loess 0.013597723
## Sayeret Shaked 0.022627877
##
## with conditional variances for "site"## $site
##
## Call:
## glm(formula = prop_rich_batha ~ land_use * year_ct + cos_td_rad +
## sin_td_rad + site, family = binomial, data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.56301 -0.36220 0.07949 0.49957 2.19630
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.38516 0.71065 0.542 0.587829
## land_useBedouin Agriculture -0.56451 0.66098 -0.854 0.393075
## land_useKKL Plantings -1.44329 0.39666 -3.639 0.000274 ***
## year_ct -0.09413 0.07317 -1.286 0.198308
## cos_td_rad -0.42804 0.42986 -0.996 0.319370
## sin_td_rad -0.01908 0.38374 -0.050 0.960344
## siteGivot Bar -0.22195 0.48034 -0.462 0.644032
## siteGivot Goral -0.17574 0.47654 -0.369 0.712287
## siteHatzerim Base -0.46526 0.67657 -0.688 0.491655
## siteMishmar Hanegev -0.37229 0.45173 -0.824 0.409858
## siteNahal Ashan -0.66266 0.57188 -1.159 0.246561
## siteNevatim Base 0.11422 1.05237 0.109 0.913567
## sitePark Loess 0.23407 0.66354 0.353 0.724271
## siteSayeret Shaked 0.03297 0.54441 0.061 0.951711
## land_useBedouin Agriculture:year_ct -0.03619 0.10937 -0.331 0.740739
## land_useKKL Plantings:year_ct 0.12154 0.08366 1.453 0.146262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 111.756 on 139 degrees of freedom
## Residual deviance: 66.864 on 124 degrees of freedom
## AIC: 320.26
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = prop_rich_batha ~ land_use * year_ct + cos_td_rad +
## sin_td_rad, family = binomial, data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4230 -0.3929 0.1248 0.5003 2.3105
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22939 0.35334 0.649 0.5162
## land_useBedouin Agriculture -0.58005 0.52244 -1.110 0.2669
## land_useKKL Plantings -1.49627 0.37909 -3.947 7.91e-05 ***
## year_ct -0.11108 0.06546 -1.697 0.0897 .
## cos_td_rad -0.44529 0.34840 -1.278 0.2012
## sin_td_rad -0.07442 0.28048 -0.265 0.7908
## land_useBedouin Agriculture:year_ct -0.02804 0.09634 -0.291 0.7710
## land_useKKL Plantings:year_ct 0.11183 0.07954 1.406 0.1597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 111.756 on 139 degrees of freedom
## Residual deviance: 74.838 on 132 degrees of freedom
## AIC: 312.23
##
## Number of Fisher Scoring iterations: 4
##
## Call:
## glm(formula = prop_rich_batha ~ land_use * year_ct + cos_td_rad +
## sin_td_rad + site, family = quasibinomial, data = P.anal,
## weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.56301 -0.36220 0.07949 0.49957 2.19630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.38516 0.50196 0.767 0.4444
## land_useBedouin Agriculture -0.56451 0.46687 -1.209 0.2289
## land_useKKL Plantings -1.44329 0.28018 -5.151 9.88e-07 ***
## year_ct -0.09413 0.05169 -1.821 0.0710 .
## cos_td_rad -0.42804 0.30363 -1.410 0.1611
## sin_td_rad -0.01908 0.27105 -0.070 0.9440
## siteGivot Bar -0.22195 0.33928 -0.654 0.5142
## siteGivot Goral -0.17574 0.33660 -0.522 0.6025
## siteHatzerim Base -0.46526 0.47789 -0.974 0.3322
## siteMishmar Hanegev -0.37229 0.31907 -1.167 0.2455
## siteNahal Ashan -0.66266 0.40394 -1.640 0.1034
## siteNevatim Base 0.11422 0.74333 0.154 0.8781
## sitePark Loess 0.23407 0.46868 0.499 0.6184
## siteSayeret Shaked 0.03297 0.38454 0.086 0.9318
## land_useBedouin Agriculture:year_ct -0.03619 0.07725 -0.468 0.6403
## land_useKKL Plantings:year_ct 0.12154 0.05909 2.057 0.0418 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.498918)
##
## Null deviance: 111.756 on 139 degrees of freedom
## Residual deviance: 66.864 on 124 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
## df AIC
## me0b 9 314.2154
## m0b 16 320.2589
## m1b 8 312.2327
## m0qb 16 NAmodel excluding site has lower AIC compared to regular binomial glm and binomial mixed model, hence use model without site. binomial models have better QQ plots (range of residuals) compared to quasibinomial.
## Start: AIC=312.23
## prop_rich_batha ~ land_use * year_ct + cos_td_rad + sin_td_rad
##
## Df Deviance AIC
## - sin_td_rad 1 74.909 310.30
## - land_use:year_ct 2 77.604 311.00
## - cos_td_rad 1 76.465 311.86
## <none> 74.838 312.23
##
## Step: AIC=310.3
## prop_rich_batha ~ land_use + year_ct + cos_td_rad + land_use:year_ct
##
## Df Deviance AIC
## - land_use:year_ct 2 77.604 309.00
## - cos_td_rad 1 76.616 310.01
## <none> 74.909 310.30
##
## Step: AIC=309
## prop_rich_batha ~ land_use + year_ct + cos_td_rad
##
## Df Deviance AIC
## - cos_td_rad 1 79.408 308.80
## <none> 77.604 309.00
## - year_ct 1 80.521 309.92
## - land_use 2 100.840 328.23
##
## Step: AIC=308.8
## prop_rich_batha ~ land_use + year_ct
##
## Df Deviance AIC
## - year_ct 1 80.859 308.25
## <none> 79.408 308.80
## - land_use 2 107.244 332.64
##
## Step: AIC=308.25
## prop_rich_batha ~ land_use
##
## Df Deviance AIC
## <none> 80.859 308.25
## - land_use 2 111.756 335.15final model includes only land use.
##
## Call:
## glm(formula = prop_rich_batha ~ land_use, family = binomial,
## data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5116 -0.3869 0.1542 0.5237 2.3376
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3962 0.1524 -2.599 0.00935 **
## land_useBedouin Agriculture -0.9671 0.2275 -4.250 2.13e-05 ***
## land_useKKL Plantings -1.1656 0.2266 -5.143 2.70e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 111.756 on 139 degrees of freedom
## Residual deviance: 80.859 on 137 degrees of freedom
## AIC: 308.25
##
## Number of Fisher Scoring iterations: 4
Test for pairwise differences among levels of land use, apply FDR correction for multiple comparisons.
## contrast estimate SE df z.ratio p.value
## Loess - Bedouin Agriculture 0.967 0.228 Inf 4.250 <.0001
## Loess - KKL Plantings 1.166 0.227 Inf 5.143 <.0001
## Bedouin Agriculture - KKL Plantings 0.199 0.238 Inf 0.834 0.4044
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: fdr method for 3 testsLoess land use has significantly higher proportion of batha specialist species compared to KKL plantings and Bedouin agriculture (the latter two do not differ)
Test whether there is a temporal trend for the proportion of batha specialist endangered species, without filtering rare species
## [1] "PROPORTION OF BATHA SPECIALIST SPECIES WHICH ARE ALSO ENDAGERED, INCLUDING RARE SPECIES"
It appears there is quasi-complete separation in the data: the land use variable predicts almost perfectly the proportion of endangered and batha specialist species: - In bedouin agriculture and KKL plantings the proportion is always 0 (except for a single plot in KKL plantings) - In natural loess the proportion is 0-0.5 As it is not possible to analyze a category with zero observations with standard logistic glm, the categories of KKL plantings and Bedouin agriculture will be unified.
## df AIC
## me0b 7 90.13603
## m0b 14 93.71395
## m1b 6 88.24418
## m0qb 14 NAcompare binomial and quasibinomial models without site variable
##
## Call:
## glm(formula = prop_rich_batha_endang ~ is_loess * year_ct + cos_td_rad +
## sin_td_rad, family = binomial, data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.41950 -0.23791 -0.14923 -0.09499 2.77008
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -0.3382 0.8255 -0.410
## is_loessBedouin agriculture / KKL plantings -4.1596 1.7665 -2.355
## year_ct -0.5025 0.1962 -2.561
## cos_td_rad -1.3609 0.9389 -1.449
## sin_td_rad -0.8680 0.6185 -1.403
## is_loessBedouin agriculture / KKL plantings:year_ct 0.2338 0.3736 0.626
## Pr(>|z|)
## (Intercept) 0.6820
## is_loessBedouin agriculture / KKL plantings 0.0185 *
## year_ct 0.0104 *
## cos_td_rad 0.1472
## sin_td_rad 0.1605
## is_loessBedouin agriculture / KKL plantings:year_ct 0.5316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 98.401 on 139 degrees of freedom
## Residual deviance: 46.729 on 134 degrees of freedom
## AIC: 88.244
##
## Number of Fisher Scoring iterations: 8
##
## Call:
## glm(formula = prop_rich_batha_endang ~ is_loess * year_ct + cos_td_rad +
## sin_td_rad, family = quasibinomial, data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.41950 -0.23791 -0.14923 -0.09499 2.77008
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.3382 0.9190 -0.368
## is_loessBedouin agriculture / KKL plantings -4.1596 1.9666 -2.115
## year_ct -0.5025 0.2184 -2.301
## cos_td_rad -1.3609 1.0453 -1.302
## sin_td_rad -0.8680 0.6885 -1.261
## is_loessBedouin agriculture / KKL plantings:year_ct 0.2338 0.4159 0.562
## Pr(>|t|)
## (Intercept) 0.7134
## is_loessBedouin agriculture / KKL plantings 0.0363 *
## year_ct 0.0229 *
## cos_td_rad 0.1952
## sin_td_rad 0.2096
## is_loessBedouin agriculture / KKL plantings:year_ct 0.5751
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.239325)
##
## Null deviance: 98.401 on 139 degrees of freedom
## Residual deviance: 46.729 on 134 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 8
Regular binomial glm and quasibinomial glm yield similar diagnostic plots, qb has somewhat lower residuals, however cannot perform AIC based model selection using qb, therefore stick with binomial glm.
## Start: AIC=88.24
## prop_rich_batha_endang ~ is_loess * year_ct + cos_td_rad + sin_td_rad
##
## Df Deviance AIC
## - is_loess:year_ct 1 47.111 86.625
## <none> 46.729 88.244
## - sin_td_rad 1 48.758 88.273
## - cos_td_rad 1 48.935 88.450
##
## Step: AIC=86.63
## prop_rich_batha_endang ~ is_loess + year_ct + cos_td_rad + sin_td_rad
##
## Df Deviance AIC
## - sin_td_rad 1 48.976 86.491
## <none> 47.111 86.625
## - cos_td_rad 1 49.285 86.800
## - year_ct 1 56.262 93.777
## - is_loess 1 71.175 108.690
##
## Step: AIC=86.49
## prop_rich_batha_endang ~ is_loess + year_ct + cos_td_rad
##
## Df Deviance AIC
## - cos_td_rad 1 49.799 85.314
## <none> 48.976 86.491
## - year_ct 1 58.239 93.754
## - is_loess 1 73.409 108.923
##
## Step: AIC=85.31
## prop_rich_batha_endang ~ is_loess + year_ct
##
## Df Deviance AIC
## <none> 49.799 85.314
## - year_ct 1 58.239 91.754
## - is_loess 1 79.956 113.471final model includes only land use.
##
## Call:
## glm(formula = prop_rich_batha_endang ~ is_loess + year_ct, family = binomial,
## data = P.anal, weights = richness)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5293 -0.2939 -0.1409 -0.0826 2.8267
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -1.5362 0.3071 -5.002
## is_loessBedouin agriculture / KKL plantings -3.6164 1.0387 -3.482
## year_ct -0.2710 0.1045 -2.592
## Pr(>|z|)
## (Intercept) 5.68e-07 ***
## is_loessBedouin agriculture / KKL plantings 0.000499 ***
## year_ct 0.009533 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 98.401 on 139 degrees of freedom
## Residual deviance: 49.799 on 137 degrees of freedom
## AIC: 85.314
##
## Number of Fisher Scoring iterations: 8
## Loading required package: Cairo## Loading required package: extrafont## Registering fonts with R
Species defined as both endangered and batha specialist are seen almost entirely in the natural loess plots. In addition, there is a significant decline in the probability of observation of such species of 86.4333228% over the 8 years of monitoring.
