Reptiles - Inland sands alpha diversity boxplots and GLMs - by plot
Ron Chen
2023-03-24
## Loading required package: readxlInland Sands
Beer Milka has Far and Near plots, both in shifting and semi-shifting dunes. Secher and Shunra East both have only far plots, in shifting and semi-shifting dunes.
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"
## [1] "only Beer Milka data:"
## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Acanthodactylus.scutellatus, Acanthodactylus.aegyptius, Chalcides.sepsoides, Stenodactylus.petrii, Lytorhynchus.diadema, Varanus.griseus, Scincus.scincus, Cerastes.vipera, Cerastes.cerastes, Chamaeleo.chamaeleon, Spalerosophis.diadema, Macroprotodon.cucullatus were included in the plot
## (the variables with highest total abundance).
## Overlapping points were shifted along the y-axis to make them visible.##
## PIPING TO 2nd MVFACTOR## Only the variables Acanthodactylus.scutellatus, Acanthodactylus.aegyptius, Chalcides.sepsoides, Stenodactylus.petrii, Lytorhynchus.diadema, Varanus.griseus, Scincus.scincus, Cerastes.vipera, Cerastes.cerastes, Chamaeleo.chamaeleon, Spalerosophis.diadema, Macroprotodon.cucullatus 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.367450089554735"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')## [1] "mixed model generates a singular fit. Hence choose regular glm."## Start: AIC=292.46
## richness ~ agriculture * year_ct + dunes * year_ct
##
## Df Deviance AIC
## - year_ct:dunes 1 24.394 290.49
## - agriculture:year_ct 1 25.814 291.91
## <none> 24.358 292.45
##
## Step: AIC=290.49
## richness ~ agriculture + year_ct + dunes + agriculture:year_ct
##
## Df Deviance AIC
## - dunes 1 24.396 288.49
## - agriculture:year_ct 1 25.850 289.95
## <none> 24.394 290.49
##
## Step: AIC=288.49
## richness ~ agriculture + year_ct + agriculture:year_ct
##
## Df Deviance AIC
## - agriculture:year_ct 1 25.852 287.95
## <none> 24.396 288.49
##
## Step: AIC=287.95
## richness ~ agriculture + year_ct
##
## Df Deviance AIC
## - agriculture 1 26.282 286.38
## <none> 25.852 287.95
## - year_ct 1 29.890 289.99
##
## Step: AIC=286.38
## richness ~ year_ct
##
## Df Deviance AIC
## <none> 26.282 286.38
## - year_ct 1 30.320 288.42year effect remains.
| Observations | 72 |
| Dependent variable | richness |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| χ²(1) | 4.04 |
| Pseudo-R² (Cragg-Uhler) | 0.06 |
| Pseudo-R² (McFadden) | 0.01 |
| AIC | 286.38 |
| BIC | 290.93 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 8.08 | 5.62 | 11.62 | 11.29 | 0.00 |
| year_ct | 0.94 | 0.89 | 1.00 | -2.01 | 0.04 |
| Standard errors: MLE |
geometric mean of abundance
Explore data.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"
## [1] "only Beer Milka data:"
Fit glm, compare gamma, gaussian (poisson inappropriate because response is not discrete)
Seems that the gamma fits the data better. Fit fixed and mixed models.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00324284 (tol = 0.002, component 1)Gamma mixed model doesn’t converge. Try Gaussian mixed model.
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ agriculture * year_ct + dunes * year_ct + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 14.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0486 -0.5761 -0.2102 0.3937 3.1642
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01435 0.1198
## Residual 0.04904 0.2215
## Number of obs: 72, groups: site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.51248 0.16682 9.066
## agricultureNear 0.12870 0.23282 0.553
## year_ct -0.04239 0.02441 -1.736
## dunesshifting -0.27878 0.19877 -1.403
## agricultureNear:year_ct -0.02542 0.03691 -0.689
## year_ct:dunesshifting 0.04466 0.03197 1.397
##
## Correlation of Fixed Effects:
## (Intr) agrclN yer_ct dnsshf agrN:_
## agricultrNr -0.339
## year_ct -0.878 0.360
## dunesshftng -0.596 0.000 0.632
## agrcltrNr:_ 0.332 -0.951 -0.378 0.000
## yr_ct:dnssh 0.575 0.000 -0.655 -0.965 0.000perform stepwise model selection of Gaussian mixed model.
## Single term deletions
##
## Model:
## gma ~ agriculture * year_ct + dunes * year_ct + (1 | site)
## npar AIC
## <none> 2.9166
## agriculture:year_ct 1 1.4252
## year_ct:dunes 1 2.9858drop interaction of agriculture and year.
## Single term deletions
##
## Model:
## gma ~ agriculture + dunes * year_ct + (1 | site)
## npar AIC
## <none> 1.42522
## agriculture 1 -0.52184
## dunes:year_ct 1 1.47941drop agriculture.
## Single term deletions
##
## Model:
## gma ~ dunes * year_ct + (1 | site)
## npar AIC
## <none> -0.52184
## dunes:year_ct 1 -0.47377This is the final model:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ dunes * year_ct + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 6.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9485 -0.5666 -0.1716 0.4145 3.3133
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01278 0.1131
## Residual 0.04814 0.2194
## Number of obs: 72, groups: site, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.54663 0.15400 10.043
## dunesshifting -0.27878 0.19692 -1.416
## year_ct -0.04874 0.02239 -2.177
## dunesshifting:year_ct 0.04466 0.03167 1.410
##
## Correlation of Fixed Effects:
## (Intr) dnsshf yer_ct
## dunesshftng -0.639
## year_ct -0.872 0.682
## dnsshftng:_ 0.617 -0.965 -0.707
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 72 |
| Dependent variable | gma |
| Type | Mixed effects linear regression |
| AIC | 18.834 |
| BIC | 32.494 |
| Pseudo-R² (fixed effects) | 0.051 |
| Pseudo-R² (total) | 0.250 |
| Est. | S.E. | t val. | d.f. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.547 | 0.154 | 10.039 | 31.528 | 0.000 |
| dunesshifting | -0.279 | 0.197 | -1.416 | 66.002 | 0.162 |
| year_ct | -0.049 | 0.022 | -2.177 | 66.002 | 0.033 |
| dunesshifting:year_ct | 0.045 | 0.032 | 1.410 | 66.002 | 0.163 |
| p values calculated using Kenward-Roger standard errors and d.f. |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.113 |
| Residual | 0.219 |
| Group | # groups | ICC |
|---|---|---|
| site | 3 | 0.210 |
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
only year is significant.
abundance
Explore data
## [1] "ABUNDANCE WITHOUT RARE SPECIES"
## [1] "only Beer Milka data:"
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 agriculture based on boxplots.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## abundance ~ agriculture * year_ct + dunes * year_ct + agriculture:dunes +
## (1 | site)
## Data: P.anal
## AIC BIC logLik deviance df.resid
## 365.0418 383.2552 -174.5209 349.0418 64
## Random effects:
## Groups Name Std.Dev.
## site (Intercept) 0.1831
## Number of obs: 72, groups: site, 3
## Fixed Effects:
## (Intercept) agricultureNear
## 2.88524 -0.74804
## year_ct dunesshifting
## -0.12723 -0.33486
## agricultureNear:year_ct year_ct:dunesshifting
## 0.08782 0.03716
## agricultureNear:dunesshifting
## 0.36216Perform stepwise model selection of mixed model.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00294099 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## abundance ~ agriculture * year_ct + dunes * year_ct + agriculture:dunes +
## (1 | site)
## npar AIC
## <none> 365.04
## agriculture:year_ct 1 365.44
## year_ct:dunes 1 363.58
## agriculture:dunes 1 366.88drop interaction of dunes and year.
## Single term deletions
##
## Model:
## abundance ~ agriculture * year_ct + dunes + agriculture:dunes +
## (1 | site)
## npar AIC
## <none> 363.58
## agriculture:year_ct 1 364.19
## agriculture:dunes 1 365.63The final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: abundance ~ agriculture * year_ct + dunes + agriculture:dunes +
## (1 | site)
## Data: P.anal
## AIC BIC logLik deviance df.resid
## 363.5834 379.5201 -174.7917 349.5834 65
## Random effects:
## Groups Name Std.Dev.
## site (Intercept) 0.1831
## Number of obs: 72, groups: site, 3
## Fixed Effects:
## (Intercept) agricultureNear
## 2.78564 -0.77215
## year_ct dunesshifting
## -0.10972 -0.12260
## agricultureNear:year_ct agricultureNear:dunesshifting
## 0.09119 0.37078
Interpretation of abundance model:
| Observations | 72 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 363.583 |
| BIC | 379.520 |
| Pseudo-R² (fixed effects) | 0.173 |
| Pseudo-R² (total) | 0.356 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 16.210 | 0.209 | 13.300 | 0.000 |
| agricultureNear | 0.462 | 0.364 | -2.124 | 0.034 |
| year_ct | 0.896 | 0.029 | -3.721 | 0.000 |
| dunesshifting | 0.885 | 0.095 | -1.286 | 0.198 |
| agricultureNear:year_ct | 1.095 | 0.056 | 1.616 | 0.106 |
| agricultureNear:dunesshifting | 1.449 | 0.185 | 2.007 | 0.045 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.183 |
| Group | # groups | ICC |
|---|---|---|
| site | 3 | 0.032 |
## Loading required namespace: broom.mixed
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
Abundance data: significant effects found for shifting dunes being less abundant than semi-shifting FAR from agriculture, and vice versa NEAR agriculture. There is a significant negative overall temporal trend.
community analysis using package MVabund
## nb po
## 111.3625 112.2653## [1] "negative binomial model is best."## nb po nb1
## 111.3625 112.2653 115.1957## [1] "including site as an explanatory variable does not improve the AIC of the model."## [1] "stepwise selection of model:"## Single term deletions
##
## Model:
## spp_no_rare ~ agriculture * year_ct + dunes * year_ct
## Df AIC
## <none> 1781.8
## agriculture:year_ct 16 1759.0
## year_ct:dunes 16 1761.5drop interaction of agriculture with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ agriculture + dunes * year_ct
## Df AIC
## <none> 1759
## agriculture 16 1762
## dunes:year_ct 16 1739drop interaction of dunes with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ agriculture + dunes + year_ct
## Df AIC
## <none> 1739.0
## agriculture 16 1742.0
## dunes 16 1767.9
## year_ct 16 1745.3final model includes agriculture, dunes and year.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 11.113 0.001 ***
## agricultureNear 4.954 0.003 **
## dunesshifting 7.011 0.001 ***
## year_ct 4.725 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 9.626, 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).## Time elapsed: 0 hr 0 min 52 sec## Analysis of Deviance Table
##
## Model: spp_no_rare ~ agriculture + dunes + year_ct
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 71
## agriculture 70 1 29.22 0.012 *
## dunes 69 1 59.97 0.001 ***
## year_ct 68 1 38.37 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acanthodactylus.aegyptius Acanthodactylus.scutellatus
## Dev Pr(>Dev) Dev
## (Intercept)
## agriculture 4.902 0.193 2.001
## dunes 14.815 0.002 18.398
## year_ct 1.317 0.894 3.879
## Cerastes.cerastes Cerastes.vipera
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.755 0.521 0.977 1.913 0.755
## dunes 0.002 0.067 1.000 0.931 0.962
## year_ct 0.380 0 0.995 0.893 0.917
## Chalcides.ocellatus Chalcides.sepsoides
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.07 0.977 0.318 0.977
## dunes 6.931 0.072 0.301 0.997
## year_ct 0 0.995 0.018 0.995
## Chamaeleo.chamaeleon Lytorhynchus.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.128 0.977 0.189 0.977
## dunes 1.646 0.851 0.209 0.997
## year_ct 0.15 0.990 1.26 0.894
## Macroprotodon.cucullatus Psammophis.schokari
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.208 0.977 3.016 0.527
## dunes 0 1.000 0 1.000
## year_ct 7.479 0.058 12.207 0.009
## Scincus.scincus Spalerosophis.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.938 0.931 5.754 0.118
## dunes 8.878 0.016 0.403 0.997
## year_ct 4.854 0.245 3.948 0.380
## Stenodactylus.petrii Testudo.kleinmanni
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.246 0.977 2.301 0.728
## dunes 0.385 0.997 0 1.000
## year_ct 1.481 0.869 0 0.995
## Trapelus.savignii Varanus.griseus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.541 0.977 6.173 0.116
## dunes 6.931 0.072 0.078 1.000
## year_ct 0.302 0.988 0.577 0.956
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.All factors (agriculture, dunes and year) have 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.
| year | Psammophis schokari | Macroprotodon cucullatus | Spalerosophis diadema | Scincus scincus | Stenodactylus petrii | Acanthodactylus scutellatus |
|---|---|---|---|---|---|---|
| 2017 | 6 | 5 | 5 | 20 | 24 | 60 |
| 2019 | 0 | 1 | 5 | 9 | 25 | 66 |
| 2021 | 0 | 0 | 0 | 9 | 16 | 30 |