Reptiles - Coastal sands alpha diversity boxplots and GLMs - by plot
Ron Chen
2023-03-25
## Loading required package: readxlCoastal Sands
Factors are proximity to settlements and infrastructure, time, dune status. Total 5 campaigns. 4 sites with 9 plots per site, except for the first campaign (2014) where there were 3 sites.
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 Acanthodactylus.scutellatus, Chalcides.sepsoides, Stenodactylus.sthenodactylus, Acanthodactylus.boskianus, Chalcides.ocellatus, Psammophis.schokari, Spalerosophis.diadema, Daboia.palaestinae, Lytorhynchus.diadema, Macroprotodon.cucullatus, Eryx.jaculus, Testudo.graeca 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, Chalcides.sepsoides, Stenodactylus.sthenodactylus, Acanthodactylus.boskianus, Chalcides.ocellatus, Psammophis.schokari, Spalerosophis.diadema, Daboia.palaestinae, Lytorhynchus.diadema, Macroprotodon.cucullatus, Eryx.jaculus, Testudo.graeca 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.604504017722961"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.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements * year_ct + dunes * year_ct + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 695.9 718.1 -341.0 681.9 170
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52992 -0.57781 -0.04921 0.53677 2.07483
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0007826 0.02798
## Number of obs: 177, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.32565 0.12868 10.302 <2e-16 ***
## settlementsNear 0.05883 0.18026 0.326 0.7441
## year_ct 0.04388 0.02430 1.806 0.0710 .
## dunesshifting -0.34767 0.18390 -1.891 0.0587 .
## settlementsNear:year_ct -0.01203 0.03447 -0.349 0.7270
## year_ct:dunesshifting 0.04731 0.03465 1.366 0.1721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf sttN:_
## settlmntsNr -0.699
## year_ct -0.871 0.618
## dunesshftng -0.686 0.490 0.606
## sttlmntsN:_ 0.611 -0.875 -0.702 -0.427
## yr_ct:dnssh 0.608 -0.434 -0.699 -0.878 0.493perform stepwise model selection of poisson mixed model.
## Single term deletions
##
## Model:
## richness ~ settlements * year_ct + dunes * year_ct + (1 | site)
## npar AIC
## <none> 695.92
## settlements:year_ct 1 694.04
## year_ct:dunes 1 695.78drop settlements:year_ct
## Single term deletions
##
## Model:
## richness ~ settlements + dunes * year_ct + (1 | site)
## npar AIC
## <none> 694.04
## settlements 1 692.04
## dunes:year_ct 1 695.16drop settlements.
## Single term deletions
##
## Model:
## richness ~ dunes * year_ct + (1 | site)
## npar AIC
## <none> 692.04
## dunes:year_ct 1 693.16Final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ dunes * year_ct + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 692.0 707.9 -341.0 682.0 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55165 -0.56142 -0.07011 0.53679 2.05772
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0007842 0.028
## Number of obs: 177, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.35483 0.09198 14.729 <2e-16 ***
## dunesshifting -0.37686 0.16033 -2.350 0.0187 *
## year_ct 0.03793 0.01730 2.192 0.0284 *
## dunesshifting:year_ct 0.05327 0.03015 1.767 0.0773 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dnsshf yer_ct
## dunesshftng -0.551
## year_ct -0.864 0.490
## dnsshftng:_ 0.490 -0.878 -0.570
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 177 |
| Dependent variable | richness |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 692.040 |
| BIC | 707.921 |
| Pseudo-R² (fixed effects) | 0.130 |
| Pseudo-R² (total) | 0.134 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 3.876 | 0.092 | 14.729 | 0.000 |
| dunesshifting | 0.686 | 0.160 | -2.350 | 0.019 |
| year_ct | 1.039 | 0.017 | 2.192 | 0.028 |
| dunesshifting:year_ct | 1.055 | 0.030 | 1.767 | 0.077 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.028 |
| Group | # groups | ICC |
|---|---|---|
| site | 5 | 0.001 |
Significant effect for year and dunes, nearly significant interaction between the two.
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
geometric mean of abundance
Explore data.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"
Fit glm, compare gamma, gaussian (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 ~ settlements * year_ct + dunes * year_ct + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 198.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4863 -0.5580 -0.0951 0.2618 6.2647
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02946 0.1716
## Residual 0.14981 0.3870
## Number of obs: 177, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.290472 0.127316 10.136
## settlementsNear -0.261664 0.144156 -1.815
## year_ct -0.002458 0.020502 -0.120
## dunesshifting 0.060449 0.138049 0.438
## settlementsNear:year_ct 0.019748 0.028749 0.687
## year_ct:dunesshifting -0.018623 0.028039 -0.664
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf sttN:_
## settlmntsNr -0.557
## year_ct -0.678 0.597
## dunesshftng -0.580 0.514 0.626
## sttlmntsN:_ 0.483 -0.864 -0.693 -0.445
## yr_ct:dnssh 0.494 -0.437 -0.712 -0.860 0.507perform stepwise model selection of Gaussian mixed model.
## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + dunes * year_ct + (1 | site)
## npar AIC
## <none> 185.98
## settlements:year_ct 1 184.47
## year_ct:dunes 1 184.44drop interaction of dunes and year.
## Single term deletions
##
## Model:
## gma ~ dunes + settlements * year_ct + (1 | site)
## npar AIC
## <none> 184.44
## dunes 1 182.51
## settlements:year_ct 1 183.89drop dunes.
## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + (1 | site)
## npar AIC
## <none> 182.51
## settlements:year_ct 1 181.94drop settlements:year_ct.
## Single term deletions
##
## Model:
## gma ~ settlements + year_ct + (1 | site)
## npar AIC
## <none> 181.94
## settlements 1 187.09
## year_ct 1 179.98drop year.
## Single term deletions
##
## Model:
## gma ~ settlements + (1 | site)
## npar AIC
## <none> 179.98
## settlements 1 185.16Only settlements remain. This is the final model:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 179
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4021 -0.5365 -0.1303 0.2264 6.3949
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02924 0.1710
## Residual 0.14804 0.3848
## Number of obs: 177, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.27196 0.08571 14.840
## settlementsNear -0.16773 0.06219 -2.697
##
## Correlation of Fixed Effects:
## (Intr)
## settlmntsNr -0.231
Bad fit but couldn’t improve it, even when using different distribution families and link functions.
| Observations | 177 |
| Dependent variable | gma |
| Type | Mixed effects linear regression |
| AIC | 186.953 |
| BIC | 199.657 |
| Pseudo-R² (fixed effects) | 0.033 |
| Pseudo-R² (total) | 0.193 |
| Est. | S.E. | t val. | d.f. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.272 | 0.086 | 14.756 | 4.355 | 0.000 |
| settlementsNear | -0.168 | 0.062 | -2.697 | 171.063 | 0.008 |
| p values calculated using Kenward-Roger standard errors and d.f. |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.171 |
| Residual | 0.385 |
| Group | # groups | ICC |
|---|---|---|
| site | 5 | 0.165 |
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
only settlement proximity is significant.
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.
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficientMixed model is rank-deficient, meaning the IVs are not linearly independent (design matrix is singular = not invertible) or that there are more terms to estimate than observations. In this case there is high collinearity between dunes and settlements (0.5 correlation, because Near settlements is always semi-shifting), and there are only 3 observations per combination of settlements-year-dunes-site. Trying to estimate another 3 interaction terms results in a singular solution.
Perform stepwise model selection of mixed model.
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Single term deletions
##
## Model:
## abundance ~ settlements * year_ct + dunes * year_ct + settlements:dunes +
## (1 | site)
## npar AIC
## <none> 887.89
## settlements:year_ct 1 886.48
## year_ct:dunes 1 892.75
## settlements:dunes 0 887.89drop interaction of settlements and year.
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Single term deletions
##
## Model:
## abundance ~ settlements + dunes * year_ct + settlements:dunes +
## (1 | site)
## npar AIC
## <none> 886.48
## dunes:year_ct 1 890.96
## settlements:dunes 0 886.48model still rank-deficient. calculate variance inflation factor (vif) to find which term to drop.
## GVIF Df GVIF^(1/(2*Df))
## settlements 1.249319 1 1.117729
## dunes 4.699763 1 2.167894
## year_ct 1.475349 1 1.214639
## dunes:year_ct 5.007124 1 2.237660
## settlements:dunes 6.927128 0 InfDrop settlements:dunes becuase vif>5
## Single term deletions
##
## Model:
## abundance ~ settlements + dunes * year_ct + (1 | site)
## npar AIC
## <none> 886.48
## settlements 1 894.98
## dunes:year_ct 1 890.96The final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: abundance ~ settlements + dunes * year_ct + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 886.5 905.5 -437.2 874.5 171
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1437 -0.8698 -0.1938 0.6705 3.7198
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02914 0.1707
## Number of obs: 177, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.72544 0.11519 14.980 < 2e-16 ***
## settlementsNear -0.24317 0.07530 -3.229 0.001240 **
## dunesshifting -0.53099 0.14116 -3.762 0.000169 ***
## year_ct 0.03689 0.01519 2.429 0.015144 *
## dunesshifting:year_ct 0.06531 0.02565 2.547 0.010875 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN dnsshf yer_ct
## settlmntsNr -0.278
## dunesshftng -0.423 0.229
## year_ct -0.546 -0.005 0.473
## dnsshftng:_ 0.330 0.002 -0.857 -0.567
Interpretation of abundance model:
| Observations | 177 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | poisson |
| Link | log |
| AIC | 886.478 |
| BIC | 905.535 |
| Pseudo-R² (fixed effects) | 0.182 |
| Pseudo-R² (total) | 0.303 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 5.615 | 0.115 | 14.980 | 0.000 |
| settlementsNear | 0.784 | 0.075 | -3.229 | 0.001 |
| dunesshifting | 0.588 | 0.141 | -3.762 | 0.000 |
| year_ct | 1.038 | 0.015 | 2.429 | 0.015 |
| dunesshifting:year_ct | 1.067 | 0.026 | 2.547 | 0.011 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.171 |
| Group | # groups | ICC |
|---|---|---|
| site | 5 | 0.028 |
## 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.
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
## 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
## nb po
## 252.1482 265.6921negative binomial model is better than poisson according to AIC comparison.
## nb po nb1
## 252.1482 265.6921 241.1193including site as an explanatory variable improves the AIC of the model. The ordination showed that most of the difference arises from one site - Caesarea, and from two species - A. boskianus and A. scutellatus, where the former is mainly in Caesarea and the latter not found there at all but found in the other sites. Therefore I create a variable to differentiate northern sites (=Caesarea) from southern sites.
stepwise selection of model:
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements * year_ct + dunes * year_ct + site_ns
## Df AIC
## <none> 3071.0
## site_ns 13 3277.9
## settlements:year_ct 13 3063.0
## year_ct:dunes 13 3064.2drop interaction of settlements with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements + dunes * year_ct + site_ns
## Df AIC
## <none> 3063.0
## settlements 13 3061.3
## site_ns 13 3267.5
## dunes:year_ct 13 3047.2drop interaction of dunes with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements + dunes + year_ct + site_ns
## Df AIC
## <none> 3047.2
## settlements 13 3045.7
## dunes 13 3046.5
## year_ct 13 3053.8
## site_ns 13 3250.7drop settlements.
## Single term deletions
##
## Model:
## spp_no_rare ~ dunes + year_ct + site_ns
## Df AIC
## <none> 3045.7
## dunes 13 3051.5
## year_ct 13 3050.4
## site_ns 13 3245.0final model includes dunes, year and site north/south.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 7.183 0.001 ***
## dunesshifting 5.266 0.003 **
## year_ct 5.612 0.001 ***
## site_nssouth 9.499 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 12.42, 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 46 sec## Analysis of Deviance Table
##
## Model: spp_no_rare ~ dunes + year_ct + site_ns
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 176
## dunes 175 1 30.44 0.001 ***
## year_ct 174 1 28.26 0.003 **
## site_ns 173 1 225.34 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acanthodactylus.boskianus Acanthodactylus.scutellatus
## Dev Pr(>Dev) Dev
## (Intercept)
## dunes 3.656 0.371 0.855
## year_ct 8.811 0.023 2.086
## site_ns 86.673 0.001 73.412
## Chalcides.ocellatus Chalcides.sepsoides
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## dunes 0.864 5.865 0.122 1.746 0.731
## year_ct 0.643 6.071 0.096 0.358 0.941
## site_ns 0.001 1.646 0.416 1.074 0.416
## Daboia.palaestinae Eryx.jaculus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## dunes 0.007 0.929 0.952 0.864
## year_ct 0.944 0.884 0.063 0.991
## site_ns 3.406 0.322 2.305 0.416
## Lytorhynchus.diadema Macroprotodon.cucullatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## dunes 0.969 0.864 0.293 0.880
## year_ct 2.154 0.643 0.471 0.941
## site_ns 22.156 0.001 6.376 0.056
## Psammophis.schokari Spalerosophis.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## dunes 11.802 0.005 1.839 0.731
## year_ct 6.267 0.096 0.854 0.884
## site_ns 1.938 0.416 22.754 0.001
## Stenodactylus.sthenodactylus Testudo.graeca
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## dunes 0.405 0.880 0.451 0.880
## year_ct 0.077 0.991 0.006 0.991
## site_ns 0.163 0.594 0.394 0.594
## Xerotyphlops.syriacus
## Dev Pr(>Dev)
## (Intercept)
## dunes 1.604 0.765
## year_ct 0.103 0.991
## site_ns 3.045 0.357
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.All factors (site_ns, 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.
