Reptiles - Loess alpha diversity boxplots and GLMs
Ron Chen
2023-02-21
Loess Covered Areas in the Northern Negev
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, Laudakia.vulgaris, Hemidactylus.turcicus, 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.56110041446309"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)
## df AIC
## mdl_r.poiss.int 6 145.0835
## mdl_r.gauss.int 7 138.4883
Seems that the gaussian fits the data better. Fit fixed and mixed models. Find the most parsimonious model using AIC (lowest AIC is the most probable model, preventing loss of information)
## Warning in glmer(data = P.anal, formula = richness ~ offset(effort) + habitat *
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly## [1] "which model is a better fit?"## df AIC
## m0 7 138.4883
## me0 8 138.7756## [1] "mixed model and fixed-effect model have essentially the same AIC. Keep the model with fewest parameters, according to the analysis workflow protocol."## Start: AIC=138.49
## richness ~ offset(effort) + habitat * year_ct
##
## Df Deviance AIC
## <none> 55.573 138.49
## - habitat:year_ct 2 64.971 140.58##
## Call: glm(formula = richness ~ offset(effort) + habitat * year_ct,
## family = gaussian, data = P.anal)
##
## Coefficients:
## (Intercept) habitatkkl plantings
## -0.23864 0.33864
## habitatloess year_ct
## 2.33864 0.05682
## habitatkkl plantings:year_ct habitatloess:year_ct
## 0.04318 -0.45682
##
## Degrees of Freedom: 38 Total (i.e. Null); 33 Residual
## Null Deviance: 68.31
## Residual Deviance: 55.57 AIC: 138.5## [1] "full model is the best model."
## [1] "summary of m0: full glm model with interactions between habitat and year:"| Observations | 39 |
| Dependent variable | richness |
| Type | Linear regression |
| χ²(5) | 12.735 |
| Pseudo-R² (Cragg-Uhler) | 0.193 |
| Pseudo-R² (McFadden) | 0.061 |
| AIC | 138.488 |
| BIC | 150.133 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | -0.239 | 0.731 | -0.326 | 0.746 |
| habitatkkl plantings | 0.339 | 1.060 | 0.319 | 0.751 |
| habitatloess | 2.339 | 1.060 | 2.206 | 0.034 |
| year_ct | 0.057 | 0.155 | 0.367 | 0.716 |
| habitatkkl plantings:year_ct | 0.043 | 0.228 | 0.189 | 0.851 |
| habitatloess:year_ct | -0.457 | 0.228 | -2.004 | 0.053 |
| Standard errors: MLE |
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools## Loading required namespace: broom.mixed
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.## Outcome is based on a total of 1 exposures
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures.
## Outcome is based on a total of 1 exposures
Significant positive effect for richness in Loess habitat, and nearly significant negative temporal trend in Loess, losing nearly one species every four years. No other significant effects found.
geometric mean of abundance
Explore data.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
## Removed 2 rows containing non-finite values (`stat_boxplot()`).
## Removed 2 rows containing non-finite values (`stat_boxplot()`).
Fit glmm, compare gamma, gaussian and poisson
## df AIC
## mdl_g.gamma.int 7 124.3535
## mdl_g.gauss.int 7 117.0376
## mdl_g.poiss.int 6 Inf
Seems that the gaussian fits the data better. Fit fixed and mixed models. Find the most parsimonious model using AIC (lowest AIC is the most probable model, preventing loss of information)
## Warning in glmer(data = P.anal, formula = gma ~ offset(effort) + habitat * :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly## [1] "which model is a better fit?"## df AIC
## m0 7 117.0376
## me0 8 126.9430## [1] "fixed-effect model is better than mixed model by AIC."## Start: AIC=117.04
## gma ~ offset(effort) + habitat * year_ct
##
## Df Deviance AIC
## - habitat:year_ct 2 34.650 116.06
## <none> 32.062 117.04
##
## Step: AIC=116.06
## gma ~ habitat + year_ct + offset(effort)
##
## Df Deviance AIC
## - year_ct 1 34.848 114.29
## <none> 34.650 116.06
## - habitat 2 41.457 119.06
##
## Step: AIC=114.29
## gma ~ habitat + offset(effort)
##
## Df Deviance AIC
## <none> 34.848 114.29
## - habitat 2 41.607 117.20##
## Call: glm(formula = gma ~ habitat + offset(effort), family = gaussian,
## data = P.anal)
##
## Coefficients:
## (Intercept) habitatkkl plantings habitatloess
## -0.64440 -0.88987 0.02658
##
## Degrees of Freedom: 38 Total (i.e. Null); 36 Residual
## Null Deviance: 41.61
## Residual Deviance: 34.85 AIC: 114.3## [1] "best model includes only habitat."## [1] "diagnostic plots of full model:"
## [1] "diagnostic plots of habitat only model:"
## [1] "best model shows a slightly worse fit than full model, with positive trend for the standardized Pearson residuals. Therefore prefer the full model."## [1] "summary of m0: full glm model with interactions between habitat and year:"| Observations | 39 |
| Dependent variable | gma |
| Type | Linear regression |
| χ²(5) | 9.545 |
| Pseudo-R² (Cragg-Uhler) | 0.243 |
| Pseudo-R² (McFadden) | 0.090 |
| AIC | 117.038 |
| BIC | 128.683 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | 0.129 | 0.555 | 0.233 | 0.818 |
| habitatkkl plantings | -1.855 | 0.805 | -2.304 | 0.028 |
| habitatloess | -1.010 | 0.805 | -1.255 | 0.218 |
| year_ct | -0.184 | 0.117 | -1.568 | 0.126 |
| habitatkkl plantings:year_ct | 0.232 | 0.173 | 1.340 | 0.190 |
| habitatloess:year_ct | 0.250 | 0.173 | 1.444 | 0.158 |
| Standard errors: MLE |
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures
## Using data P.anal from global environment. This could cause incorrect
## results if P.anal has been altered since the model was fit. You can
## manually provide the data to the "data =" argument.
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures.
## Outcome is based on a total of 1 exposures
gma: only significant effects found for KKL plantings habitat, having lower gma. Negative temporal trend for bedouin agriculture is near significance.
abundance
Explore data
## [1] "ABUNDANCE WITHOUT RARE SPECIES"
Fit glmm, compare gamma, gaussian and poisson
## df AIC
## mdl_a.gamma.int 7 209.3605
## mdl_a.gauss.int 7 229.3747
## mdl_a.poiss.int 6 215.8985
Seems that the Gamma fits the data better. Fit fixed and mixed models. Find the most parsimonious model using AIC (lowest AIC is the most probable model, preventing loss of information)
## [1] "which model is a better fit?"## df AIC
## m0 7 209.3605
## me0 8 210.9597## [1] "fixed-effect model is better than mixed model by AIC."## Start: AIC=209.36
## abundance ~ offset(effort) + habitat * year_ct
##
## Df Deviance AIC
## - habitat:year_ct 2 15.258 209.24
## <none> 13.771 209.36
##
## Step: AIC=209.61
## abundance ~ habitat + year_ct + offset(effort)##
## Call: glm(formula = abundance ~ habitat + year_ct + offset(effort),
## family = Gamma(link = "log"), data = P.anal)
##
## Coefficients:
## (Intercept) habitatkkl plantings habitatloess
## -0.63185 -0.33144 0.02717
## year_ct
## -0.08540
##
## Degrees of Freedom: 38 Total (i.e. Null); 35 Residual
## Null Deviance: 17.93
## Residual Deviance: 15.26 AIC: 209.6## [1] "best model is the full model."## [1] "diagnostic plots of full model:"
## [1] "summary of m0: full glm model with interactions between habitat and year:"| Observations | 39 |
| Dependent variable | abundance |
| Type | Generalized linear model |
| Family | Gamma |
| Link | log |
| χ²(5) | 4.160 |
| Pseudo-R² (Cragg-Uhler) | 0.247 |
| Pseudo-R² (McFadden) | 0.053 |
| AIC | 209.360 |
| BIC | 221.005 |
| Est. | S.E. | t val. | p | |
|---|---|---|---|---|
| (Intercept) | -0.460 | 0.349 | -1.318 | 0.196 |
| habitatkkl plantings | -1.168 | 0.506 | -2.309 | 0.027 |
| habitatloess | 0.094 | 0.506 | 0.186 | 0.853 |
| year_ct | -0.127 | 0.074 | -1.725 | 0.094 |
| habitatkkl plantings:year_ct | 0.197 | 0.109 | 1.807 | 0.080 |
| habitatloess:year_ct | -0.020 | 0.109 | -0.187 | 0.853 |
| Standard errors: MLE |
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures.
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures
## Outcome is based on a total of 1 exposures
Abundance data: significant effects found for KKL plantings habitat being less abundant than the other two units. There is a nearly significant positive temporal trend for KKL plantings, and negative for bedouin agriculture and loess.
community analysis using package MVabund
## nb po
## 101.2271 100.1315## [1] "poisson model is best."## nb po po1
## 101.2271 100.1315 105.9190## [1] "including site as an explanatory variable improves the AIC of the model. However as this cannot be included as a random variable, the model without 'site' is chosen."
## Time elapsed: 0 hr 0 min 11 sec## Warning in anova.manyglm(mva_m1, cor.type = "shrink", p.uni = "adjusted"): The
## likelihood ratio test can only be used if correlation matrix of the abundances
## is is assumed to be the Identity matrix. The Wald Test will be used.## Time elapsed: 0 hr 0 min 6 sec##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 9.949 0.001 ***
## habitatkkl plantings 4.327 0.007 **
## habitatloess 5.002 0.002 **
## year_ct 6.058 0.001 ***
## habitatkkl plantings:year_ct 2.515 0.278
## habitatloess:year_ct 1.248 0.777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 11.51, 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 ~ habitat * year_ct
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 96
## habitat 94 2 242.98 0.001 ***
## year_ct 93 1 70.44 0.001 ***
## habitat:year_ct 91 2 15.49 0.544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Acanthodactylus.beershebensis Chalcides.ocellatus
## Dev Pr(>Dev) Dev
## (Intercept)
## habitat 79.846 0.001 11.045
## year_ct 20.226 0.010 1.879
## habitat:year_ct 0 0.996 0.806
## Chamaeleo.chamaeleon Hemidactylus.turcicus
## Pr(>Dev) Dev Pr(>Dev) Dev
## (Intercept)
## habitat 0.188 16.249 0.055 4.951
## year_ct 0.779 0.179 0.992 2.258
## habitat:year_ct 0.996 0.092 0.996 4.536
## Heremites.vittatus Laudakia.vulgaris
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 0.409 7.747 0.296 11.341 0.188
## year_ct 0.751 0.009 0.992 9.748 0.045
## habitat:year_ct 0.687 0.76 0.996 7 0.342
## Malpolon.insignitus Mesalina.bahaeldini
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 3.454 0.409 8.83 0.265
## year_ct 0.2 0.992 0.711 0.947
## habitat:year_ct 0.186 0.996 0.113 0.996
## Ophisops.elegans Ptyodactylus.guttatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## habitat 95.034 0.001 4.487 0.409
## year_ct 35.186 0.001 0.042 0.992
## habitat:year_ct 1.604 0.980 0.388 0.996
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.
## [1] "Acanthodactylus beershebensis - שנונית באר שבע"## (Intercept) habitat year_ct habitat:year_ct
## NA 0.001 0.010 0.996
## [1] "Ophisops elegans - עינחש"## (Intercept) habitat year_ct habitat:year_ct
## NA 0.001 0.001 0.980
## [1] "Chamaeleo chamaeleon - זיקית"## (Intercept) habitat year_ct habitat:year_ct
## NA 0.055 0.992 0.996
## [1] "Laudakia vulgaris - חרדון מצוי"## (Intercept) habitat year_ct habitat:year_ct
## NA 0.188 0.045 0.342