Reptiles - Coastal sands alpha diversity boxplots and GLMs - by plot
Ron Chen
2023-04-25
## Loading required package: readxl## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following objects are masked from 'package:data.table':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionCoastal 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.
IMPORTANT: all observations with activity of “tracks”, “eggshell” and “eggs” are converted to a single observation per plot-year-species combination, with an individual count of 1, regardless of the reported count or the number of observations per such combination. This is done to reduce the inherent bias in these signs which tends to over-estimate abundance, as some of the counts of tracks are essentially activity of the same individual whereas eggs and eggshells are not a measure of abundance.
Model gma, abundance and richness
Richness done with rare species. Abundance and mean abundance done without rare species. Full models include cosine and sine of the time difference from June 21st (in radians); do not include temperature in the shade & sun as these have missing values.
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).
| richness | abundance | gma | year_ct | settlements | dunes | cos_td_rad | sin_td_rad | timediff_Jun21 | site | isCaesarea | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. :1.000 | Min. : 1.000 | Min. :1.000 | Min. :1.000 | Far :114 | semi-shifting:113 | Min. :-0.51974 | Min. :0.5234 | Length:170 | Ashdod :45 | Mode :logical | |
| 1st Qu.:3.000 | 1st Qu.: 4.000 | 1st Qu.:1.000 | 1st Qu.:2.000 | Near: 56 | shifting : 57 | 1st Qu.:-0.25119 | 1st Qu.:0.8406 | Class :difftime | Ashkelon :45 | FALSE:135 | |
| Median :4.000 | Median : 6.000 | Median :1.189 | Median :4.000 | NA | NA | Median : 0.09025 | Median :0.9243 | Mode :numeric | Caesarea :35 | TRUE :35 | |
| Mean :4.553 | Mean : 6.418 | Mean :1.222 | Mean :4.359 | NA | NA | Mean : 0.13122 | Mean :0.8873 | NA | Netiv Haasara: 9 | NA | |
| 3rd Qu.:6.000 | 3rd Qu.: 8.000 | 3rd Qu.:1.340 | 3rd Qu.:6.000 | NA | NA | 3rd Qu.: 0.54163 | 3rd Qu.:0.9813 | NA | Zikim :36 | NA | |
| Max. :9.000 | Max. :17.000 | Max. :2.466 | Max. :8.000 | NA | NA | Max. : 0.85208 | Max. :0.9999 | NA | NA | NA |
| settlements | dunes | year_ct | cos_td_rad | sin_td_rad | timediff_Jun21 | site | isCaesarea | |
|---|---|---|---|---|---|---|---|---|
| settlements | 1.0000000 | -0.4977827 | -0.0054301 | 0.0050577 | -0.0050829 | -0.0050135 | -0.0018192 | -0.0163868 |
| dunes | -0.4977827 | 1.0000000 | 0.0027030 | 0.0052568 | 0.0055945 | -0.0047070 | 0.0009055 | 0.0081571 |
| year_ct | -0.0054301 | 0.0027030 | 1.0000000 | 0.2633772 | 0.0606911 | -0.2426777 | 0.0598414 | 0.1236906 |
| cos_td_rad | 0.0050577 | 0.0052568 | 0.2633772 | 1.0000000 | -0.7177103 | -0.9984043 | 0.0334750 | -0.0354436 |
| sin_td_rad | -0.0050829 | 0.0055945 | 0.0606911 | -0.7177103 | 1.0000000 | 0.7493729 | -0.2764657 | 0.1048951 |
| timediff_Jun21 | -0.0050135 | -0.0047070 | -0.2426777 | -0.9984043 | 0.7493729 | 1.0000000 | -0.0438222 | 0.0447119 |
| site | -0.0018192 | 0.0009055 | 0.0598414 | 0.0334750 | -0.2764657 | -0.0438222 | 1.0000000 | 0.1110133 |
| isCaesarea | -0.0163868 | 0.0081571 | 0.1236906 | -0.0354436 | 0.1048951 | 0.0447119 | 0.1110133 | 1.0000000 |
cos_td_rad is in nearly perfect correlation with timediff_Jun21, but model does not converge when using timediff_Jun21 nor with td_sc instead of cos_td_rad. 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.387837016750881"Overdispersion parameter is <1 (underdispersion) and therefore Poisson is preferable to negative binomial.
## boundary (singular) fit: see help('isSingular')singular fit using the random effect isCaesarea.
The following results in a singular fit (and obviously also a model which includes both plot and site as random effects): me2 <- glmer(data = P.anal, formula = richness ~ settlements * year_ct + dunes * year_ct + cos_td_rad + sin_td_rad + (1|plot_coord), family = poisson)
the following results in an error because plot coord identifies the plot uniquely, hence incorrect to nest it inside site; if it were only numbers e.g. 1,2,3 then it would have made sense to nest it, as every site would have plots 1,2,3. me3 <- glmer(data = P.anal, formula = richness ~ settlements * year_ct + dunes * year_ct + cos_td_rad + sin_td_rad + (plot_coord|site), family = poisson)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 647.1 675.3 -314.5 629.1 161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2102 -0.3821 -0.0922 0.3466 1.9404
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0002334 0.01528
## Number of obs: 170, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.595634 0.392380 4.067 4.77e-05 ***
## settlementsNear 0.022762 0.188368 0.121 0.904
## year_ct 0.026958 0.027252 0.989 0.323
## dunesshifting -0.183530 0.194347 -0.944 0.345
## cos_td_rad 0.338533 0.145530 2.326 0.020 *
## sin_td_rad -0.264755 0.441746 -0.599 0.549
## settlementsNear:year_ct -0.007494 0.036545 -0.205 0.838
## year_ct:dunesshifting 0.024664 0.037038 0.666 0.505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf cs_td_ sn_td_ sttN:_
## settlmntsNr -0.217
## year_ct -0.004 0.603
## dunesshftng -0.212 0.487 0.586
## cos_td_rad -0.764 -0.025 -0.306 -0.033
## sin_td_rad -0.939 -0.026 -0.301 -0.023 0.806
## sttlmntsN:_ 0.189 -0.886 -0.674 -0.430 0.022 0.026
## yr_ct:dnssh 0.193 -0.438 -0.665 -0.891 0.029 0.018 0.490sin_td_rad is highly correlated with the model intercept. Center the variable to reduce this correlation.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad_c + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 647.1 675.3 -314.5 629.1 161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.21022 -0.38210 -0.09219 0.34661 1.94042
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0002331 0.01527
## Number of obs: 170, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.360760 0.136530 9.967 <2e-16 ***
## settlementsNear 0.022715 0.188368 0.121 0.904
## year_ct 0.026948 0.027252 0.989 0.323
## dunesshifting -0.183581 0.194346 -0.945 0.345
## cos_td_rad 0.338579 0.145530 2.327 0.020 *
## sin_td_rad_c -0.264602 0.441742 -0.599 0.549
## settlementsNear:year_ct -0.007484 0.036545 -0.205 0.838
## year_ct:dunesshifting 0.024674 0.037038 0.666 0.505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf cs_td_ sn_t__ sttN:_
## settlmntsNr -0.697
## year_ct -0.876 0.603
## dunesshftng -0.674 0.487 0.586
## cos_td_rad 0.118 -0.025 -0.306 -0.033
## sin_td_rd_c 0.171 -0.026 -0.301 -0.023 0.806
## sttlmntsN:_ 0.616 -0.886 -0.674 -0.430 0.022 0.026
## yr_ct:dnssh 0.606 -0.438 -0.665 -0.891 0.029 0.018 0.490year_ct is highly correlated with the model intercept. Center the variable to reduce this correlation.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richness ~ settlements * year_ct_c + dunes * year_ct_c + cos_td_rad +
## sin_td_rad_c + (1 | site)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 647.1 675.3 -314.5 629.1 161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2102 -0.3821 -0.0922 0.3466 1.9405
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0002329 0.01526
## Number of obs: 170, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.478220 0.065842 22.451 <2e-16 ***
## settlementsNear -0.009901 0.087646 -0.113 0.910
## year_ct_c 0.026959 0.027252 0.989 0.323
## dunesshifting -0.076027 0.089001 -0.854 0.393
## cos_td_rad 0.338556 0.145532 2.326 0.020 *
## sin_td_rad_c -0.264656 0.441754 -0.599 0.549
## settlementsNear:year_ct_c -0.007498 0.036545 -0.205 0.837
## year_ct_c:dunesshifting 0.024662 0.037038 0.666 0.506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yr_ct_ dnsshf cs_td_ sn_t__ stN:__
## settlmntsNr -0.654
## year_ct_c -0.013 0.071
## dunesshftng -0.643 0.487 0.074
## cos_td_rad -0.306 -0.013 -0.306 -0.020
## sin_td_rd_c -0.189 -0.009 -0.301 -0.017 0.806
## sttlmntN:__ 0.063 -0.087 -0.674 -0.051 0.022 0.026
## yr_ct_c:dns 0.057 -0.051 -0.665 -0.132 0.029 0.018 0.490perform stepwise model selection of poisson mixed model.
## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## richness ~ settlements * year_ct_c + dunes * year_ct_c + cos_td_rad +
## sin_td_rad_c + (1 | site)
## npar AIC
## <none> 647.08
## cos_td_rad 1 650.55
## sin_td_rad_c 1 645.45
## settlements:year_ct_c 1 645.13
## year_ct_c:dunes 1 645.53singular fit during backward selection. try using isCaesarea as random effect instead of site.
## boundary (singular) fit: see help('isSingular')singular fit as well. resort to glm.
## Start: AIC=645.09
## richness ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad
##
## Df Deviance AIC
## - settlements:year_ct 1 62.871 643.13
## - sin_td_rad 1 63.193 643.45
## - year_ct:dunes 1 63.274 643.53
## <none> 62.830 645.09
## - cos_td_rad 1 68.695 648.95
##
## Step: AIC=643.13
## richness ~ settlements + year_ct + dunes + cos_td_rad + sin_td_rad +
## year_ct:dunes
##
## Df Deviance AIC
## - settlements 1 62.889 641.15
## - sin_td_rad 1 63.229 641.49
## - year_ct:dunes 1 63.645 641.90
## <none> 62.871 643.13
## - cos_td_rad 1 68.760 647.02
##
## Step: AIC=641.15
## richness ~ year_ct + dunes + cos_td_rad + sin_td_rad + year_ct:dunes
##
## Df Deviance AIC
## - sin_td_rad 1 63.248 639.50
## - year_ct:dunes 1 63.661 639.92
## <none> 62.889 641.15
## - cos_td_rad 1 68.771 645.03
##
## Step: AIC=639.5
## richness ~ year_ct + dunes + cos_td_rad + year_ct:dunes
##
## Df Deviance AIC
## - year_ct:dunes 1 64.027 638.28
## <none> 63.248 639.50
## - cos_td_rad 1 85.343 659.60
##
## Step: AIC=638.28
## richness ~ year_ct + dunes + cos_td_rad
##
## Df Deviance AIC
## - dunes 1 64.712 636.97
## <none> 64.027 638.28
## - year_ct 1 67.336 639.59
## - cos_td_rad 1 86.007 658.26
##
## Step: AIC=636.97
## richness ~ year_ct + cos_td_rad
##
## Df Deviance AIC
## <none> 64.712 636.97
## - year_ct 1 68.017 638.27
## - cos_td_rad 1 86.664 656.92drop year_ct because \(\Delta AIC<2\).
## Single term deletions
##
## Model:
## richness ~ cos_td_rad
## Df Deviance AIC
## <none> 68.017 638.27
## cos_td_rad 1 95.109 663.37only time of year remains. Final model:
##
## Call:
## glm(formula = richness ~ cos_td_rad, family = poisson, data = P.anal)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.75247 -0.43108 -0.05784 0.33841 2.02961
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.44068 0.04010 35.927 < 2e-16 ***
## cos_td_rad 0.43963 0.08441 5.209 1.9e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 95.109 on 169 degrees of freedom
## Residual deviance: 68.017 on 168 degrees of freedom
## AIC: 638.27
##
## Number of Fisher Scoring iterations: 4
| Observations | 170 |
| Dependent variable | richness |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| χ²(1) | 27.092 |
| Pseudo-R² (Cragg-Uhler) | 0.150 |
| Pseudo-R² (McFadden) | 0.041 |
| AIC | 638.275 |
| BIC | 644.546 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 4.224 | 3.904 | 4.569 | 35.927 | 0.000 |
| cos_td_rad | 1.552 | 1.315 | 1.831 | 5.209 | 0.000 |
| Standard errors: MLE |
There is a statistically significant effect for time of sampling. No statistically significant effect for the other factors tested.
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. After getting a bad fit using all data points, I exclude the 3 outliers, rows 3,65,135 Fit fixed and mixed models.
Mixed model converged:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: -43.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09246 -0.68534 -0.07181 0.66976 2.67212
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.007181 0.08474
## Residual 0.034480 0.18569
## Number of obs: 167, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.960122 0.168786 5.688
## settlementsNear -0.219564 0.070924 -3.096
## year_ct -0.014937 0.010818 -1.381
## dunesshifting -0.064566 0.070920 -0.910
## cos_td_rad 0.192904 0.057062 3.381
## sin_td_rad 0.363150 0.185640 1.956
## settlementsNear:year_ct 0.022714 0.014103 1.611
## year_ct:dunesshifting 0.007608 0.014075 0.541
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf cs_td_ sn_td_ sttN:_
## settlmntsNr -0.189
## year_ct 0.004 0.596
## dunesshftng -0.189 0.518 0.595
## cos_td_rad -0.703 -0.037 -0.299 -0.035
## sin_td_rad -0.926 -0.030 -0.274 -0.030 0.773
## sttlmntsN:_ 0.160 -0.866 -0.685 -0.449 0.036 0.031
## yr_ct:dnssh 0.164 -0.450 -0.684 -0.867 0.031 0.027 0.517sin_td_rad is highly correlated with the model intercept. Center the variable to reduce this correlation.
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad_c + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: -43.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09246 -0.68534 -0.07181 0.66976 2.67212
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.007181 0.08474
## Residual 0.034480 0.18569
## Number of obs: 167, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.282304 0.064308 19.940
## settlementsNear -0.219564 0.070924 -3.096
## year_ct -0.014937 0.010818 -1.381
## dunesshifting -0.064566 0.070920 -0.910
## cos_td_rad 0.192904 0.057062 3.381
## sin_td_rad_c 0.363150 0.185640 1.956
## settlementsNear:year_ct 0.022714 0.014103 1.611
## year_ct:dunesshifting 0.007608 0.014075 0.541
##
## Correlation of Fixed Effects:
## (Intr) sttlmN yer_ct dnsshf cs_td_ sn_t__ sttN:_
## settlmntsNr -0.574
## year_ct -0.691 0.596
## dunesshftng -0.574 0.518 0.595
## cos_td_rad 0.136 -0.037 -0.299 -0.035
## sin_td_rd_c 0.131 -0.030 -0.274 -0.030 0.773
## sttlmntsN:_ 0.499 -0.866 -0.685 -0.449 0.036 0.031
## yr_ct:dnssh 0.499 -0.450 -0.684 -0.867 0.031 0.027 0.517perform stepwise model selection of Gaussian mixed model.
## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad_c + (1 | site)
## npar AIC
## <none> -66.799
## cos_td_rad 1 -57.105
## sin_td_rad_c 1 -64.797
## settlements:year_ct 1 -66.113
## year_ct:dunes 1 -68.495drop year_ct:dunes.
## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + dunes + cos_td_rad + sin_td_rad_c +
## (1 | site)
## npar AIC
## <none> -68.495
## dunes 1 -69.685
## cos_td_rad 1 -58.919
## sin_td_rad_c 1 -66.554
## settlements:year_ct 1 -67.992drop dunes.
## Single term deletions
##
## Model:
## gma ~ settlements * year_ct + cos_td_rad + sin_td_rad_c + (1 |
## site)
## npar AIC
## <none> -69.685
## cos_td_rad 1 -60.246
## sin_td_rad_c 1 -67.808
## settlements:year_ct 1 -69.196drop settlements:year because \(\Delta AIC<2\).
## Single term deletions
##
## Model:
## gma ~ settlements + year_ct + cos_td_rad + sin_td_rad_c + (1 |
## site)
## npar AIC
## <none> -69.196
## settlements 1 -59.620
## year_ct 1 -70.796
## cos_td_rad 1 -60.133
## sin_td_rad_c 1 -67.495drop year.
## Single term deletions
##
## Model:
## gma ~ settlements + cos_td_rad + sin_td_rad_c + (1 | site)
## npar AIC
## <none> -70.796
## settlements 1 -61.278
## cos_td_rad 1 -61.412
## sin_td_rad_c 1 -69.474drop sine of sampling time because \(\Delta AIC<2\) and should keep simpler model.
## Single term deletions
##
## Model:
## gma ~ settlements + cos_td_rad + (1 | site)
## npar AIC
## <none> -69.474
## settlements 1 -60.108
## cos_td_rad 1 -62.552Settlements and sampling time remain. This is the final model:
## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ settlements + cos_td_rad + (1 | site)
## Data: P.anal
##
## REML criterion at convergence: -64.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0203 -0.7660 -0.1444 0.6568 2.5593
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.007587 0.0871
## Residual 0.034894 0.1868
## Number of obs: 167, groups: site, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.21002 0.04360 27.752
## settlementsNear -0.10467 0.03063 -3.417
## cos_td_rad 0.10593 0.03555 2.980
##
## Correlation of Fixed Effects:
## (Intr) sttlmN
## settlmntsNr -0.235
## cos_td_rad -0.062 -0.005## $site
## (Intercept)
## Ashdod 0.11948511
## Ashkelon -0.03729732
## Caesarea -0.02703680
## Netiv Haasara -0.08914389
## Zikim 0.03399290
##
## with conditional variances for "site"## $site
Bad fit but couldn’t improve it, even when using different distribution families and link functions, and even after throwing out 3 outliers.
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools| Observations | 167 |
| Dependent variable | gma |
| Type | Mixed effects linear regression |
| AIC | -54.866 |
| BIC | -39.277 |
| Pseudo-R² (fixed effects) | 0.095 |
| Pseudo-R² (total) | 0.256 |
| Est. | S.E. | t val. | d.f. | p | |
|---|---|---|---|---|---|
| (Intercept) | 1.210 | 0.044 | 27.642 | 4.398 | 0.000 |
| settlementsNear | -0.105 | 0.031 | -3.417 | 160.083 | 0.001 |
| cos_td_rad | 0.106 | 0.036 | 2.949 | 163.995 | 0.004 |
| p values calculated using Kenward-Roger standard errors and d.f. |
| Group | Parameter | Std. Dev. |
|---|---|---|
| site | (Intercept) | 0.087 |
| Residual | 0.187 |
| Group | # groups | ICC |
|---|---|---|
| site | 5 | 0.179 |
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
Take results with a grain of salt due to BAD fit. Settlement proximity and time of sampling are the only significant factors.
abundance
Explore data
## [1] "ABUNDANCE WITHOUT RARE SPECIES"
PHI>1, hence choose negative binomial. 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.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00452366 (tol = 0.002, component 1)Mixed model failed to converge. This is not due to underdispersion, as theta converged on finite (and low) value. try to replace site with isCaesarea.
Model converged. 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
## 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 +
## cos_td_rad + sin_td_rad + (1 | isCaesarea)
## npar AIC
## <none> 798.03
## cos_td_rad 1 828.40
## sin_td_rad 1 799.55
## settlements:year_ct 1 797.27
## year_ct:dunes 1 798.85
## settlements:dunes 0 798.03Mixed 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.
calculate variance inflation factor (vif) to find which term to drop.
## GVIF Df GVIF^(1/(2*Df))
## settlements 5.986603 1 2.446754
## year_ct 3.031833 1 1.741216
## dunes 6.168762 1 2.483699
## cos_td_rad 3.200146 1 1.788895
## sin_td_rad 3.108638 1 1.763133
## settlements:year_ct 6.652683 1 2.579280
## year_ct:dunes 7.034432 1 2.652250
## settlements:dunes 342.937478 0 InfDrop settlements:dunes becuase vif>5
## Single term deletions
##
## Model:
## abundance ~ settlements * year_ct + dunes * year_ct + cos_td_rad +
## sin_td_rad + (1 | isCaesarea)
## npar AIC
## <none> 798.03
## cos_td_rad 1 828.40
## sin_td_rad 1 799.55
## settlements:year_ct 1 797.27
## year_ct:dunes 1 798.85drop settlements:year_ct.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0030196 (tol = 0.002, component 1)## Single term deletions
##
## Model:
## abundance ~ settlements + dunes * year_ct + cos_td_rad + sin_td_rad +
## (1 | isCaesarea)
## npar AIC
## <none> 797.27
## settlements 1 807.07
## cos_td_rad 1 827.20
## sin_td_rad 1 798.69
## dunes:year_ct 1 797.00drop dunes:year.
## Single term deletions
##
## Model:
## abundance ~ settlements + dunes + year_ct + cos_td_rad + sin_td_rad +
## (1 | isCaesarea)
## npar AIC
## <none> 796.99
## settlements 1 806.65
## dunes 1 799.41
## year_ct 1 796.25
## cos_td_rad 1 826.39
## sin_td_rad 1 798.33Drop year.
## Single term deletions
##
## Model:
## abundance ~ settlements + dunes + cos_td_rad + sin_td_rad + (1 |
## isCaesarea)
## npar AIC
## <none> 796.25
## settlements 1 805.97
## dunes 1 798.70
## cos_td_rad 1 842.89
## sin_td_rad 1 801.08settlements, dunes and time of sampling remain. The final model:
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(73.6907) ( log )
## Formula: abundance ~ settlements + dunes + cos_td_rad + sin_td_rad + (1 |
## isCaesarea)
## Data: P.anal
##
## AIC BIC logLik deviance df.resid
## 796.3 818.2 -391.1 782.3 163
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7668 -0.7381 -0.0843 0.3827 4.2061
##
## Random effects:
## Groups Name Variance Std.Dev.
## isCaesarea (Intercept) 0.01088 0.1043
## Number of obs: 170, groups: isCaesarea, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.93450 0.35457 2.636 0.008399 **
## settlementsNear -0.26915 0.07895 -3.409 0.000652 ***
## dunesshifting -0.16102 0.07640 -2.108 0.035071 *
## cos_td_rad 0.82197 0.12227 6.722 1.79e-11 ***
## sin_td_rad 0.93950 0.36264 2.591 0.009579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sttlmN dnsshf cs_td_
## settlmntsNr -0.092
## dunesshftng -0.084 0.451
## cos_td_rad -0.771 -0.009 -0.019
## sin_td_rad -0.962 -0.004 -0.016 0.772
Interpretation of abundance model:
| Observations | 170 |
| Dependent variable | abundance |
| Type | Mixed effects generalized linear model |
| Family | Negative Binomial(73.6907) |
| Link | log |
| AIC | 796.251 |
| BIC | 818.201 |
| Pseudo-R² (fixed effects) | 0.341 |
| Pseudo-R² (total) | 0.384 |
| exp(Est.) | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 2.546 | 0.355 | 2.636 | 0.008 |
| settlementsNear | 0.764 | 0.079 | -3.409 | 0.001 |
| dunesshifting | 0.851 | 0.076 | -2.108 | 0.035 |
| cos_td_rad | 2.275 | 0.122 | 6.722 | 0.000 |
| sin_td_rad | 2.559 | 0.363 | 2.591 | 0.010 |
| Group | Parameter | Std. Dev. |
|---|---|---|
| isCaesarea | (Intercept) | 0.104 |
| Group | # groups | ICC |
|---|---|---|
| isCaesarea | 2 | 0.026 |
## 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.
community analysis using package MVabund
## nb po
## 257.1134 267.6093## [1] "POISSON"
## [1] "NEGATIVE BINOMIAL"
negative binomial model is better than poisson according to residuals and AIC comparison.
## nb po nb1 nb2
## 257.1134 267.6093 237.8767 243.2391The addition of the explanatory variable isCaesarea has the most significant impact on improving the AIC of the model. stepwise selection of model:
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements * year_ct + dunes * year_ct + isCaesarea +
## cos_td_rad + sin_td_rad
## Df AIC
## <none> 2854.5
## isCaesarea 12 3085.4
## cos_td_rad 12 2902.5
## sin_td_rad 12 2868.0
## settlements:year_ct 12 2859.1
## year_ct:dunes 12 2853.0drop interaction of dunes with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements * year_ct + dunes + isCaesarea + cos_td_rad +
## sin_td_rad
## Df AIC
## <none> 2853.0
## dunes 12 2854.3
## isCaesarea 12 3081.4
## cos_td_rad 12 2899.3
## sin_td_rad 12 2864.4
## settlements:year_ct 12 2845.9drop interaction of settlements with year.
## Single term deletions
##
## Model:
## spp_no_rare ~ settlements + year_ct + dunes + isCaesarea + cos_td_rad +
## sin_td_rad
## Df AIC
## <none> 2845.9
## settlements 12 2851.9
## year_ct 12 2849.7
## dunes 12 2846.5
## isCaesarea 12 3071.8
## cos_td_rad 12 2889.0
## sin_td_rad 12 2857.6should drop dunes because \(\Delta AIC<2\), but when including it in the model it comes out as significant.
final model includes settlements, dunes, year, sampling time of year and isCaesarea.
##
## Test statistics:
## wald value Pr(>wald)
## (Intercept) 5.684 0.001 ***
## settlementsNear 5.567 0.001 ***
## dunesshifting 4.781 0.006 **
## year_ct 5.024 0.004 **
## isCaesareaTRUE 9.938 0.001 ***
## cos_td_rad 8.018 0.001 ***
## sin_td_rad 6.009 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test statistic: 16.5, 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 3 min 52 sec## Analysis of Deviance Table
##
## Model: spp_no_rare ~ settlements + dunes + year_ct + isCaesarea + cos_td_rad + sin_td_rad
##
## Multivariate test:
## Res.Df Df.diff Dev Pr(>Dev)
## (Intercept) 169
## settlements 168 1 25.63 0.007 **
## dunes 167 1 20.51 0.020 *
## year_ct 166 1 23.14 0.016 *
## isCaesarea 165 1 232.29 0.001 ***
## cos_td_rad 164 1 81.46 0.001 ***
## sin_td_rad 163 1 35.63 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)
## settlements 0.064 0.984 15.11
## dunes 3.041 0.458 0.089
## year_ct 8.153 0.032 1.649
## isCaesarea 85.068 0.001 83.576
## cos_td_rad 9.769 0.012 11.634
## sin_td_rad 0.008 1.000 25.743
## Chalcides.ocellatus Chalcides.sepsoides
## Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 0.004 2.783 0.521 0.047 0.984
## dunes 0.991 2.096 0.641 0.516 0.921
## year_ct 0.738 5.034 0.143 0.013 0.994
## isCaesarea 0.001 1.495 0.417 1.492 0.417
## cos_td_rad 0.007 7.211 0.026 4.895 0.053
## sin_td_rad 0.001 1.207 0.898 0.012 1.000
## Daboia.palaestinae Eryx.jaculus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 1.153 0.796 2.397 0.594
## dunes 1.296 0.803 0 0.991
## year_ct 0.499 0.966 0.024 0.994
## isCaesarea 3.103 0.315 2.233 0.417
## cos_td_rad 9.046 0.016 10.587 0.010
## sin_td_rad 0.008 1.000 0.75 0.916
## Lytorhynchus.diadema Macroprotodon.cucullatus
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 1.221 0.796 0.856 0.810
## dunes 0.717 0.912 0.085 0.991
## year_ct 1.338 0.781 0.258 0.966
## isCaesarea 22.906 0.001 6.647 0.049
## cos_td_rad 5.387 0.053 13.206 0.004
## sin_td_rad 0.719 0.916 1.129 0.898
## Psammophis.schokari Spalerosophis.diadema
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 1.786 0.699 0.002 0.984
## dunes 10.633 0.010 1.496 0.800
## year_ct 5.574 0.114 0.493 0.966
## isCaesarea 2.053 0.417 23.366 0.001
## cos_td_rad 5.611 0.053 0.55 0.431
## sin_td_rad 0.024 1.000 1.476 0.865
## Stenodactylus.sthenodactylus Testudo.graeca
## Dev Pr(>Dev) Dev Pr(>Dev)
## (Intercept)
## settlements 0.042 0.984 0.165 0.981
## dunes 0.033 0.991 0.505 0.921
## year_ct 0.051 0.994 0.055 0.994
## isCaesarea 0.04 0.614 0.313 0.614
## cos_td_rad 0.002 0.910 3.565 0.088
## sin_td_rad 0.002 1.000 4.556 0.253
## Arguments:
## Test statistics calculated assuming uncorrelated response (for faster computation)
## P-value calculated using 999 iterations via PIT-trap resampling.All factors (settlements, dunes, year, sampling time of year and isCaesarea) 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.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Significant temporal effect only on A. boskianus (rise). Temporal decline for A. scuttelatus (not statistically significant). plot fitted temporal trends - sampling date: Sept 15th
A. boskianus site: Caesarea
A. scuttelatus site: NOT Caesarea
