Reptiles - Inland sands alpha diversity boxplots and GLMs
Ron Chen
2023-03-23
## 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.215326328105353"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=115.77
## richness ~ agriculture * dunes * year_ct
##
## Df Deviance AIC
## - agriculture:dunes:year_ct 1 3.4458 113.77
## <none> 3.4452 115.77
##
## Step: AIC=113.77
## richness ~ agriculture + dunes + year_ct + agriculture:dunes +
## agriculture:year_ct + dunes:year_ct
##
## Df Deviance AIC
## - dunes:year_ct 1 3.4473 111.77
## - agriculture:dunes 1 3.5058 111.83
## - agriculture:year_ct 1 3.8034 112.13
## <none> 3.4458 113.77
##
## Step: AIC=111.77
## richness ~ agriculture + dunes + year_ct + agriculture:dunes +
## agriculture:year_ct
##
## Df Deviance AIC
## - agriculture:dunes 1 3.5067 109.83
## - agriculture:year_ct 1 3.8042 110.13
## <none> 3.4473 111.77
##
## Step: AIC=109.83
## richness ~ agriculture + dunes + year_ct + agriculture:year_ct
##
## Df Deviance AIC
## - dunes 1 3.6781 108.00
## - agriculture:year_ct 1 3.8636 108.19
## <none> 3.5067 109.83
##
## Step: AIC=108
## richness ~ agriculture + year_ct + agriculture:year_ct
##
## Df Deviance AIC
## - agriculture:year_ct 1 4.0351 106.36
## <none> 3.6781 108.00
##
## Step: AIC=106.36
## richness ~ agriculture + year_ct
##
## Df Deviance AIC
## - agriculture 1 4.0919 104.42
## - year_ct 1 5.0643 105.39
## <none> 4.0351 106.36
##
## Step: AIC=104.42
## richness ~ year_ct
##
## Df Deviance AIC
## - year_ct 1 5.1211 103.44
## <none> 4.0919 104.42
##
## Step: AIC=103.45
## richness ~ 1There is no significant trend for any of the independent variables tested - year, proximity to agriculture and dune status (shifting vs. semi-shifting). The best model is a constant model.
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.0266882 (tol = 0.002, component 1)## [1] "Mixed model generates singular fit, choose regular glm."## Start: AIC=39.2
## gma ~ agriculture * dunes * year_ct
##
## Df Deviance AIC
## - agriculture:dunes:year_ct 1 0.78968 37.839
## <none> 0.75904 39.203
##
## Step: AIC=38.16
## gma ~ agriculture + dunes + year_ct + agriculture:dunes + agriculture:year_ct +
## dunes:year_ct
##
## Df Deviance AIC
## - dunes:year_ct 1 0.78971 36.159
## - agriculture:year_ct 1 0.82465 36.906
## - agriculture:dunes 1 0.87209 37.920
## <none> 0.78968 38.158
##
## Step: AIC=36.16
## gma ~ agriculture + dunes + year_ct + agriculture:dunes + agriculture:year_ct
##
## Df Deviance AIC
## - agriculture:year_ct 1 0.82549 34.969
## - agriculture:dunes 1 0.87280 36.041
## <none> 0.78971 36.159
##
## Step: AIC=35.23
## gma ~ agriculture + dunes + year_ct + agriculture:dunes
##
## Df Deviance AIC
## - year_ct 1 0.88215 34.549
## - agriculture:dunes 1 0.90814 35.155
## <none> 0.82549 35.229
##
## Step: AIC=34.83
## gma ~ agriculture + dunes + agriculture:dunes
##
## Df Deviance AIC
## - agriculture:dunes 1 0.96504 34.741
## <none> 0.88215 34.831
##
## Step: AIC=35
## gma ~ agriculture + dunes## [1] "best model includes only agriculture."## [1] "diagnostic plots of full model:"
## [1] "diagnostic plots of agriculture only model:"
## [1] "summary of m1"| Observations | 24 |
| Dependent variable | gma |
| Type | Generalized linear model |
| Family | Gamma |
| Link | inverse |
| χ²(2) | 0.123 |
| Pseudo-R² (Cragg-Uhler) | 0.160 |
| Pseudo-R² (McFadden) | 0.097 |
| AIC | 35.000 |
| BIC | 39.713 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | 0.502 | 0.033 | 15.119 | 0.000 |
| agricultureNear | -0.050 | 0.045 | -1.123 | 0.274 |
| dunesshifting | -0.049 | 0.041 | -1.183 | 0.250 |
| Standard errors: MLE |
gma: model is overall not statistically significant, and so is the only factor included (agriculture).
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.
the mixed model shows a singular fit - this arises because the combination of dunes and year is singular in two of the three sites. Choose regular model (fixed effects).
## Start: AIC=172.41
## abundance ~ agriculture * dunes * year_ct
##
## Df Deviance AIC
## - agriculture:dunes:year_ct 1 36.559 170.97
## <none> 36.005 172.41
##
## Step: AIC=170.97
## abundance ~ agriculture + dunes + year_ct + agriculture:dunes +
## agriculture:year_ct + dunes:year_ct
##
## Df Deviance AIC
## - dunes:year_ct 1 37.101 169.51
## <none> 36.559 170.97
## - agriculture:year_ct 1 38.959 171.37
## - agriculture:dunes 1 40.402 172.81
##
## Step: AIC=169.51
## abundance ~ agriculture + dunes + year_ct + agriculture:dunes +
## agriculture:year_ct
##
## Df Deviance AIC
## <none> 37.101 169.51
## - agriculture:year_ct 1 39.706 170.11
## - agriculture:dunes 1 41.149 171.56## [1] "remaining terms are agriculture, dunes, year, an interaction between agriculture and year and between agriculture and dunes."
## [1] "summary of m1:"| Observations | 24 |
| Dependent variable | abundance |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| χ²(5) | 19.367 |
| Pseudo-R² (Cragg-Uhler) | 0.554 |
| Pseudo-R² (McFadden) | 0.109 |
| AIC | 169.508 |
| BIC | 176.576 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 49.422 | 34.680 | 70.432 | 21.580 | 0.000 |
| agricultureNear | 0.535 | 0.265 | 1.081 | -1.743 | 0.081 |
| dunesshifting | 0.885 | 0.734 | 1.067 | -1.285 | 0.199 |
| year_ct | 0.896 | 0.846 | 0.949 | -3.718 | 0.000 |
| agricultureNear:dunesshifting | 1.449 | 1.008 | 2.082 | 2.006 | 0.045 |
| agricultureNear:year_ct | 1.095 | 0.981 | 1.224 | 1.615 | 0.106 |
| Standard errors: MLE |
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools## Loading required namespace: broom.mixed
## [1] "Only year and dunes are significant."
Abundance data: significant effects found for shifting dunes being less abundant than semi-shifting. 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.001 ***
## dunesshifting 7.011 0.001 ***
## year_ct 4.725 0.017 *
## ---
## 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 1 min 8 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.010 **
## 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.196 2.001
## dunes 14.815 0.001 18.398
## year_ct 1.317 0.886 3.879
## Cerastes.cerastes Cerastes.vipera
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.746 0.521 0.978 1.913 0.746
## dunes 0.001 0.067 0.999 0.931 0.957
## year_ct 0.376 0 0.997 0.893 0.937
## Chalcides.ocellatus Chalcides.sepsoides
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.07 0.981 0.318 0.978
## dunes 6.931 0.075 0.301 0.996
## year_ct 0 0.997 0.018 0.997
## Chamaeleo.chamaeleon Lytorhynchus.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.128 0.981 0.189 0.981
## dunes 1.646 0.821 0.209 0.997
## year_ct 0.15 0.995 1.26 0.886
## Macroprotodon.cucullatus Psammophis.schokari
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.208 0.981 3.016 0.511
## dunes 0 0.999 0 0.999
## year_ct 7.479 0.070 12.207 0.008
## Scincus.scincus Spalerosophis.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.938 0.947 5.754 0.116
## dunes 8.878 0.025 0.403 0.996
## year_ct 4.854 0.234 3.948 0.376
## Stenodactylus.petrii Testudo.kleinmanni
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.246 0.978 2.301 0.712
## dunes 0.385 0.996 0 0.999
## year_ct 1.481 0.865 0 0.997
## Trapelus.savignii Varanus.griseus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## agriculture 0.541 0.978 6.173 0.112
## dunes 6.931 0.075 0.078 0.999
## year_ct 0.302 0.985 0.577 0.962
## 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.## Warning: Using alpha for a discrete variable is not advised.
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Using alpha for a discrete variable is not advised.
| 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 |