Reptiles - Loess covered areas in the Northern Negev - alpha diversity boxplots and GLMs - by plot
Ron Chen
2023-03-26
## Loading required package: readxlLoess covered areas in the Northern Negev
Factors are landscape use (=habitat), time. habitats are natural loess, extensive bedouin agriculture, KKL plantings. Total 4 campaigns. 7 sites, total of 9 plots per each of 3 habitats. total 27 plots.
Model gma, abundance and richness
Richness done with rare species. Abundance and mean abundance done without rare species.
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).
## [1] "RICHNESS WITH RARE SPECIES"
## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Ophisops.elegans, Chalcides.ocellatus, Acanthodactylus.beershebensis, Heremites.vittatus, Mesalina.bahaeldini, Hemidactylus.turcicus, Laudakia.vulgaris, Chamaeleo.chamaeleon, Ptyodactylus.guttatus, Acanthodactylus.boskianus, Malpolon.insignitus, Eumeces.schneiderii were included in the plot
## (the variables with highest total abundance).
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.343922713188086"Overdispersion parameter is <1 (underdispersion) and therefore Poisson is preferable to negative binomial. Check to see if gamma / gaussian are a better fit (essentially OLS regression) Poisson family for glm / glme.
## boundary (singular) fit: see help('isSingular')mixed model is singular. try to remove interaction.
## boundary (singular) fit: see help('isSingular')use fixed model. perform stepwise model selection of poisson fixed model.
## Start: AIC=277.57
## richness ~ habitat * year_ct
##
## Df Deviance AIC
## - habitat:year_ct 2 33.292 274.88
## <none> 31.985 277.57
##
## Step: AIC=274.88
## richness ~ habitat + year_ct
##
## Df Deviance AIC
## - habitat 2 34.443 272.03
## - year_ct 1 33.988 273.57
## <none> 33.292 274.88
##
## Step: AIC=272.03
## richness ~ year_ct
##
## Df Deviance AIC
## - year_ct 1 35.134 270.72
## <none> 34.443 272.03
##
## Step: AIC=270.72
## richness ~ 1All terms were removed in stepwise selection.
geometric mean of abundance
Explore data.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"
Fit glm. cannot use Gamma because there is a single zero (poisson inappropriate because response is not discrete)
Gamma and gaussian seem similar. Use gaussian. Fit fixed and mixed models.
Mixed model converged:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ habitat * year_ct + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 351.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1704 -0.3175 -0.2078 0.2234 7.3963
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01189 0.109
## Residual 1.94393 1.394
## Number of obs: 99, groups: site, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.1475 0.5023 6.266
## habitatkkl plantings -1.6829 0.7287 -2.309
## habitatloess -0.3071 0.6973 -0.440
## year_ct -0.3503 0.1128 -3.105
## habitatkkl plantings:year_ct 0.3035 0.1620 1.874
## habitatloess:year_ct 0.1473 0.1558 0.945
##
## Correlation of Fixed Effects:
## (Intr) hbttkp hbttls yer_ct hpln:_
## hbttkklplnt -0.686
## habitatloss -0.720 0.500
## year_ct -0.872 0.602 0.629
## hplntngs:y_ 0.608 -0.877 -0.438 -0.696
## hbttlss:yr_ 0.632 -0.436 -0.864 -0.724 0.504perform stepwise model selection of Gaussian mixed model.
## Single term deletions
##
## Model:
## gma ~ habitat * year_ct + (1 | site)
## npar AIC
## <none> 356.99
## habitat:year_ct 2 356.67Improvement by removing the habitat*year interaction is quite minor. tried to keep the interaction but the interaction of year with both levels of habitat was not significant. Hence remove interaction.
## Single term deletions
##
## Model:
## gma ~ habitat + year_ct + (1 | site)
## npar AIC
## <none> 356.67
## habitat 2 357.21
## year_ct 1 364.32This is also the final model.
Interpretation of gma model:
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 99 |
| Dependent variable | gma |
| Type | Mixed effects linear regression |
| AIC | 363.272 |
| BIC | 378.843 |
| Pseudo-R² (fixed effects) | 0.130 |
| Pseudo-R² (total) | 0.138 |
| Est. | S.E. | t val. | d.f. | p | |
|---|---|---|---|---|---|
| (Intercept) | 2.576 | 0.367 | 7.016 | 36.161 | 0.000 |
| habitatkkl plantings | -0.478 | 0.382 | -1.253 | 23.528 | 0.223 |
| habitatloess | 0.262 | 0.393 | 0.666 | 10.710 | 0.519 |
| year_ct | -0.203 | 0.065 | -3.097 | 93.398 | 0.003 |
| p values calculated using Kenward-Roger standard errors and d.f. |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.136 |
| Residual | 1.404 |
| Group | # groups | ICC |
|---|---|---|
| site | 7 | 0.009 |
## Loading required namespace: broom.mixed
Significant decrease in gma over time, with an interaction with natural loess. Significantly lower gma in KKL plantings compared to the other two habitats.
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
abundance
Explore data
## [1] "ABUNDANCE WITHOUT RARE SPECIES"
Choose poisson based on the type of data (abundance). Fit fixed and mixed models. Choose mixed model if possible, otherwise choose a model with fixed-effects only. add interaction between dunes and settlements based on boxplots.
Perform stepwise model selection of mixed model.
## Single term deletions
##
## Model:
## abundance ~ habitat * year_ct + (1 | site)
## npar AIC
## <none> 385.61
## habitat:year_ct 2 393.10This is already the final model.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: abundance ~ habitat * year_ct + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 385.6 403.8 -185.8 371.6 92
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4739 -0.7389 -0.1036 0.5033 4.1774
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.00508 0.07127
## Number of obs: 99, groups: site, 7
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.84874 0.19007 9.727 < 2e-16 ***
## habitatkkl plantings -1.17692 0.32739 -3.595 0.000325 ***
## habitatloess -0.31017 0.26684 -1.162 0.245091
## year_ct -0.25162 0.05438 -4.627 3.71e-06 ***
## habitatkkl plantings:year_ct 0.26649 0.07951 3.352 0.000804 ***
## habitatloess:year_ct 0.12772 0.07137 1.790 0.073534 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hbttkp hbttls yer_ct hpln:_
## hbttkklplnt -0.571
## habitatloss -0.716 0.420
## year_ct -0.813 0.472 0.578
## hplntngs:y_ 0.558 -0.847 -0.401 -0.683
## hbttlss:yr_ 0.619 -0.360 -0.802 -0.762 0.520
Interpretation of abundance model:
| Observations | 99 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 385.613 |
| BIC | 403.779 |
| Pseudo-R² (fixed effects) | 0.306 |
| Pseudo-R² (total) | 0.316 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 6.352 | 0.190 | 9.727 | 0.000 |
| habitatkkl plantings | 0.308 | 0.327 | -3.595 | 0.000 |
| habitatloess | 0.733 | 0.267 | -1.162 | 0.245 |
| year_ct | 0.778 | 0.054 | -4.627 | 0.000 |
| habitatkkl plantings:year_ct | 1.305 | 0.080 | 3.352 | 0.001 |
| habitatloess:year_ct | 1.136 | 0.071 | 1.790 | 0.074 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.071 |
| Group | # groups | ICC |
|---|---|---|
| site | 7 | 0.005 |
Significant decrese in abundance over time with interactions with habitat. Lower abundance in KKL plantings.
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
community analysis using package MVabund
Let’s look for visual patterns. Also, test to see if sites are distinct as far as community composition is concerned.
fit community model.
## nb po
## 103.3515 110.2898negative binomial model is better than poisson according to AIC comparison.
## nb po nb1
## 103.3515 110.2898 103.8338including site as an explanatory variable does not improves the AIC of the model, therefore do not include site in the fixed model. stepwise selection of model:
## Single term deletions
##
## Model:
## spp_no_rare ~ habitat * year_ct
## Df AIC
## <none> 1033.5
## habitat:year_ct 20 1007.2drop interaction of habitat with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ habitat + year_ct
## Df AIC
## <none> 1007.2
## habitat 20 1095.7
## year_ct 10 1015.1drop habitat.
## Single term deletions
##
## Model:
## spp_no_rare ~ year_ct
## Df AIC
## <none> 1095.7
## year_ct 10 1095.3final model includes nothing… strange considering that there is a loess-specific species and that there is a decrease in abundance and gma. try to fit poisson.
## Single term deletions
##
## Model:
## spp_no_rare ~ habitat * year_ct
## Df AIC
## <none> 1102.9
## habitat:year_ct 20 1079.6drop interaction of habitat with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ habitat + year_ct
## Df AIC
## <none> 1079.6
## habitat 20 1291.1
## year_ct 10 1125.0Final model includes habitat and year.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 11.564 0.055 .
## habitatkkl plantings 5.644 0.001 ***
## habitatloess 7.076 0.293
## year_ct 7.328 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 11.47, p-value: 0.327
## 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).## Time elapsed: 0 hr 0 min 6 sec## Analysis of Deviance Table
##
## Model: spp_no_rare ~ habitat + year_ct
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 98
## habitat 96 2 253.69 0.001 ***
## year_ct 95 1 65.35 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acanthodactylus.beershebensis Chalcides.ocellatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 79.1 0.001 11.114 0.243
## year_ct 17.74 0.023 2.392 0.721
## Chamaeleo.chamaeleon Hemidactylus.turcicus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 18.083 0.043 6.215 0.367
## year_ct 0.162 0.998 0 0.998
## Heremites.vittatus Laudakia.vulgaris
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 7.967 0.367 11.531 0.243
## year_ct 0.058 0.998 9.151 0.089
## Malpolon.insignitus Mesalina.bahaeldini
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 3.367 0.380 8.616 0.367
## year_ct 0.277 0.997 0.477 0.997
## Ophisops.elegans Ptyodactylus.guttatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 103.319 0.001 4.379 0.380
## year_ct 35.004 0.002 0.087 0.998
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.Proximity to settlements has a statistically significant effect on community composition.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
coefficients are exponentiated because the link function used for negative binomial is log -> above 1 is positive and below 1 is negative.
| habitat | Ophisops elegans | Chamaeleo chamaeleon | Acanthodactylus beershebensis |
|---|---|---|---|
| bedouin agriculture | 63 | 0 | 0 |
| kkl plantings | 2 | 10 | 0 |
| loess | 3 | 1 | 36 |
| year | Ophisops elegans | Chamaeleo chamaeleon | Acanthodactylus beershebensis |
|---|---|---|---|
| 2014 | 31 | 2 | 23 |
| 2016 | 30 | 5 | 7 |
| 2018 | 3 | 2 | 1 |
| 2020 | 4 | 2 | 5 |