## Loading required package: readxl

Desert-Med transition zone

Factors are proximity to settlements, time. Total 4 campaigns. 5 sites with 6 plots per site.

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 Heremites.vittatus, Ophisops.elegans, Ptyodactylus.guttatus, Chalcides.ocellatus, Laudakia.vulgaris, Chalcides.guentheri, Eumeces.schneiderii, Testudo.graeca, Hemidactylus.turcicus, Ablepharus.rueppellii, Letheobia.simoni, Platyceps.collaris 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.739281100916263"

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 + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    408.0    421.8   -199.0    398.0      112 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2571 -0.6365 -0.1173  0.5837  3.0235 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.02953  0.1718  
## Number of obs: 117, groups:  site, 5
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.69107    0.20162   3.428 0.000609 ***
## settlementsNear          0.09775    0.24949   0.392 0.695210    
## year_ct                  0.03922    0.03857   1.017 0.309185    
## settlementsNear:year_ct  0.01431    0.05174   0.276 0.782174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sttlmN yer_ct
## settlmntsNr -0.688              
## year_ct     -0.821  0.663       
## sttlmntsN:_  0.613 -0.887 -0.746

perform stepwise model selection of poisson mixed model.

## Single term deletions
## 
## Model:
## richness ~ settlements * year_ct + (1 | site)
##                     npar    AIC
## <none>                   407.99
## settlements:year_ct    1 406.07

drop settlements:year_ct

## Single term deletions
## 
## Model:
## richness ~ settlements + year_ct + (1 | site)
##             npar    AIC
## <none>           406.07
## settlements    1 405.98
## year_ct        1 407.45

drop settlements.

## Single term deletions
## 
## Model:
## richness ~ year_ct + (1 | site)
##         npar    AIC
## <none>       405.98
## year_ct    1 407.25

Final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ year_ct + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    406.0    414.3   -200.0    400.0      114 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.20618 -0.68083 -0.09647  0.59266  2.74668 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.02905  0.1704  
## Number of obs: 117, groups:  site, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.74282    0.14607   5.085 3.67e-07 ***
## year_ct      0.04651    0.02571   1.809   0.0705 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## year_ct -0.753

## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
Observations 117
Dependent variable richness
Type Mixed effects generalized linear model
Family poisson
Link log
AIC 405.977
BIC 414.263
Pseudo-R² (fixed effects) 0.030
Pseudo-R² (total) 0.108
Fixed Effects
exp(Est.) S.E. z val. p
(Intercept) 2.102 0.146 5.085 0.000
year_ct 1.048 0.026 1.809 0.070
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.170
Grouping Variables
Group # groups ICC
site 5 0.028

Nearly significant effect for year.

## 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. cannot use Gamma because there is a single zero (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 + (1 | site)
##    Data: P.anal
## 
## REML criterion at convergence: 220.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2918 -0.4621 -0.3628  0.1707  6.1692 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.003033 0.05507 
##  Residual             0.343378 0.58598 
## Number of obs: 117, groups:  site, 5
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              1.267946   0.164582   7.704
## settlementsNear          0.116834   0.223738   0.522
## year_ct                 -0.007364   0.034939  -0.211
## settlementsNear:year_ct -0.009983   0.048370  -0.206
## 
## Correlation of Fixed Effects:
##             (Intr) sttlmN yer_ct
## settlmntsNr -0.719              
## year_ct     -0.871  0.641       
## sttlmntsN:_  0.629 -0.875 -0.722

perform stepwise model selection of Gaussian mixed model.

## Single term deletions
## 
## Model:
## gma ~ settlements * year_ct + (1 | site)
##                     npar    AIC
## <none>                   215.75
## settlements:year_ct    1 213.78

drop interaction of settlements and year.

## Single term deletions
## 
## Model:
## gma ~ settlements + year_ct + (1 | site)
##             npar    AIC
## <none>           213.78
## settlements    1 212.28
## year_ct        1 212.07

drop year.

## Single term deletions
## 
## Model:
## gma ~ settlements + (1 | site)
##             npar    AIC
## <none>           212.07
## settlements    1 210.58

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 settlements based on boxplots.

Perform stepwise model selection of mixed model.

## Single term deletions
## 
## Model:
## abundance ~ settlements * year_ct + (1 | site)
##                     npar    AIC
## <none>                   506.48
## settlements:year_ct    1 505.15

drop interaction of settlements and year.

## Single term deletions
## 
## Model:
## abundance ~ settlements + year_ct + (1 | site)
##             npar    AIC
## <none>           505.15
## settlements    1 505.13
## year_ct        1 504.59

should drop year by AIC model selection, but will test in model because of visual trend.

The final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: abundance ~ settlements + year_ct + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    505.1    516.2   -248.6    497.1      113 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8727 -0.8157 -0.1938  0.4321  4.5751 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.06335  0.2517  
## Number of obs: 117, groups:  site, 5
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.95212    0.16718   5.695 1.23e-08 ***
## settlementsNear  0.14529    0.10299   1.411    0.158    
## year_ct          0.02758    0.02291   1.204    0.229    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sttlmN
## settlmntsNr -0.342       
## year_ct     -0.578  0.016

Interpretation of abundance model:

Observations 117
Dependent variable abundance
Type Mixed effects generalized linear model
Family poisson
Link log
AIC 505.147
BIC 516.196
Pseudo-R² (fixed effects) 0.026
Pseudo-R² (total) 0.207
Fixed Effects
exp(Est.) S.E. z val. p
(Intercept) 2.591 0.167 5.695 0.000
settlementsNear 1.156 0.103 1.411 0.158
year_ct 1.028 0.023 1.204 0.229
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.252
Grouping Variables
Group # groups ICC
site 5 0.060 
## Loading required namespace: broom.mixed

None of the effects is significant.

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 
## 136.1871 140.9204

negative binomial model is better than poisson according to AIC comparison.

##       nb       po      nb1 
## 136.1871 140.9204 134.6786

including site as an explanatory variable improves the AIC of the model; however not significantly and the ordination above does not show sites as distinct. Therefore do not include site in the fixed model. stepwise selection of model:

## Single term deletions
## 
## Model:
## spp_no_rare ~ settlements * year_ct
##                     Df    AIC
## <none>                 1770.4
## settlements:year_ct 13 1759.5

drop interaction of settlements with year.

## Single term deletions
## 
## Model:
## spp_no_rare ~ settlements + year_ct
##             Df    AIC
## <none>         1759.5
## settlements 13 1765.3
## year_ct     13 1751.7

drop year.

## Single term deletions
## 
## Model:
## spp_no_rare ~ settlements
##             Df    AIC
## <none>         1751.7
## settlements 13 1757.6

final model includes only settlements.

## 
## Test statistics:
##                 wald value Pr(>wald)    
## (Intercept)         15.611     0.001 ***
## settlementsNear      5.323     0.005 ** 
## --- 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Test statistic:  5.323, p-value: 0.005 
## 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 20 sec
## Analysis of Deviance Table
## 
## Model: spp_no_rare ~ settlements
## 
## Multivariate test:
##             Res.Df Df.diff   Dev Pr(>Dev)  
## (Intercept)    116                         
## settlements    115       1 31.89    0.014 *
## ---
## 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)                                                            
## settlements                 5.484    0.276               9.879    0.030
##             Chalcides.ocellatus          Dolichophis.jugularis         
##                             Dev Pr(>Dev)                   Dev Pr(>Dev)
## (Intercept)                                                            
## settlements               1.313    0.853                 1.438    0.853
##             Eumeces.schneiderii          Hemidactylus.turcicus         
##                             Dev Pr(>Dev)                   Dev Pr(>Dev)
## (Intercept)                                                            
## settlements               1.624    0.853                 1.228    0.853
##             Heremites.vittatus          Laudakia.vulgaris         
##                            Dev Pr(>Dev)               Dev Pr(>Dev)
## (Intercept)                                                       
## settlements              0.287    0.853             2.931    0.691
##             Letheobia.simoni          Ophisops.elegans         
##                          Dev Pr(>Dev)              Dev Pr(>Dev)
## (Intercept)                                                    
## settlements            0.391    0.853            1.139    0.853
##             Platyceps.collaris          Ptyodactylus.guttatus         
##                            Dev Pr(>Dev)                   Dev Pr(>Dev)
## (Intercept)                                                           
## settlements              2.843    0.691                 2.147    0.810
##             Testudo.graeca         
##                        Dev Pr(>Dev)
## (Intercept)                        
## settlements          1.181    0.853
## Arguments:
##  Test statistics calculated assuming uncorrelated response (for faster computation) 
## P-value calculated using 999 iterations via PIT-trap resampling.

Proximity to settlements has 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.