Reptiles - Forest conifer forest alpha diversity boxplots and GLMs - by plot
Ron Chen
2023-03-26
## Loading required package: readxlForest conifer forest
Factor is time. Total 5 campaigns. 3 sub-units, 5 sites in each sub-unit with 3 plots per site far from settlements.
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 Ablepharus.rueppellii, Mediodactylus.orientalis, Hemidactylus.turcicus, Heremites.vittatus, Phoenicolacerta.laevis, Laudakia.vulgaris, Chalcides.ocellatus, Ophisops.elegans, Chalcides.guentheri, Chamaeleo.chamaeleon, Platyceps.collaris, 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.378494787083707"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')## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ year_ct + (1 | subunit)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 497.4 506.9 -245.7 491.4 173
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6157 -0.5595 0.1228 0.2241 1.7381
##
## Random effects:
## Groups Name Variance Std.Dev.
## subunit (Intercept) 0.002309 0.04806
## Number of obs: 176, groups: subunit, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.521766 0.120108 4.344 1.4e-05 ***
## year_ct 0.007706 0.022753 0.339 0.735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## year_ct -0.840perform stepwise model selection of poisson mixed model.
## Single term deletions
##
## Model:
## richness ~ year_ct + (1 | subunit)
## npar AIC
## <none> 497.36
## year_ct 1 495.47drop year. No fixed effects remain. Final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ (1 | subunit)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 495.5 501.8 -245.7 491.5 174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5885 -0.5689 0.1596 0.2366 1.6968
##
## Random effects:
## Groups Name Variance Std.Dev.
## subunit (Intercept) 0.002491 0.04991
## Number of obs: 176, groups: subunit, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.55577 0.06554 8.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1geometric mean of abundance
Explore data.
## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"
Fit glm. cannot use Gamma because there are zeros (poisson inappropriate because response is not discrete) Use gaussian. Fit fixed and mixed models.
Mixed model converged:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ year_ct + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: 483.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5353 -0.4978 -0.1824 0.0554 4.7454
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.1496 0.3868
## Residual 0.7996 0.8942
## Number of obs: 176, groups: site, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.422711 0.170987 8.321
## year_ct -0.009593 0.027229 -0.352
##
## Correlation of Fixed Effects:
## (Intr)
## year_ct -0.706perform stepwise model selection of Gaussian mixed model.
## Single term deletions
##
## Model:
## gma ~ year_ct + (1 | site)
## npar AIC
## <none> 483.32
## year_ct 1 481.45drop year. None of the fixed factors remain in the model - none of them contribute enough.
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 subunit based on boxplots.
Perform stepwise model selection of mixed model.
## Single term deletions
##
## Model:
## abundance ~ year_ct + (1 | site)
## npar AIC
## <none> 713.49
## year_ct 1 711.61drop year. None of the fixed factors remain in the model - none of them contribute enough.
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
## 189.2303 212.3053negative binomial model is better than poisson according to AIC comparison.
## nb po nb11 nb1
## 189.2303 212.3053 191.1502 187.8103including site as an explanatory variable does not improve the model explanatory power. including subunit as an explanatory variable improves the AIC of the model; however not significantly and the ordination above does not show subunits as distinct. Therefore do not include subunit nor site in the fixed model. stepwise selection of model:
## Single term deletions
##
## Model:
## spp_no_rare ~ year_ct
## Df AIC
## <none> 1703.1
## year_ct 9 1711.7final model includes year.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 14.664 0.001 ***
## year_ct 4.844 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 4.844, p-value: 0.002
## 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 22 sec## Analysis of Deviance Table
##
## Model: spp_no_rare ~ year_ct
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 175
## year_ct 174 1 26.59 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Tests:
## Ablepharus.rueppellii Chalcides.guentheri
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## year_ct 0.559 0.917 7.262 0.043
## Chalcides.ocellatus Hemidactylus.turcicus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## year_ct 0.062 0.996 3.591 0.307
## Heremites.vittatus Laudakia.vulgaris
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## year_ct 10.408 0.007 0.046 0.996
## Mediodactylus.orientalis Ophisops.elegans
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## year_ct 2.934 0.409 0.003 0.996
## Phoenicolacerta.laevis
## Dev Pr(>Dev)
## (Intercept)
## year_ct 1.726 0.648
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.There is a statistically significant temporal effect on community composition.
## Warning in geom_point(aes(fill = factor(year_ct.p < 0.1)), color = "black", :
## Ignoring unknown parameters: `linewidth`## 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.