## 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, union

Inland Sands

Beer Milka has Far and Near plots, both in shifting and semi-shifting dunes. Secher and Shunra East both have only far plots, in shifting and semi-shifting dunes.

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"

## [1] "only Beer Milka data:"

## Overlapping points were shifted along the y-axis to make them visible.
## 
##  PIPING TO 2nd MVFACTOR
## Only the variables Acanthodactylus.scutellatus, Acanthodactylus.aegyptius, Chalcides.sepsoides, Stenodactylus.petrii, Lytorhynchus.diadema, Scincus.scincus, Varanus.griseus, Cerastes.vipera, Cerastes.cerastes, Chamaeleo.chamaeleon, Spalerosophis.diadema, Macroprotodon.cucullatus 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, Acanthodactylus.aegyptius, Chalcides.sepsoides, Stenodactylus.petrii, Lytorhynchus.diadema, Scincus.scincus, Varanus.griseus, Cerastes.vipera, Cerastes.cerastes, Chamaeleo.chamaeleon, Spalerosophis.diadema, Macroprotodon.cucullatus were included in the plot 
## (the variables with highest total abundance).
richness abundance gma year_ct agriculture dunes cos_td_rad sin_td_rad timediff_Jun21 site
Min. : 2.000 Min. : 2.00 Min. :1.000 Min. :4 Far :54 semi-shifting:36 Min. :0.5982 Min. :-0.1373 Length:72 Beer Milka :36
1st Qu.: 5.000 1st Qu.: 6.00 1st Qu.:1.000 1st Qu.:4 Near:18 shifting :36 1st Qu.:0.8429 1st Qu.: 0.1329 Class :difftime Secher :18
Median : 6.000 Median : 8.00 Median :1.210 Median :6 NA NA Median :0.9226 Median : 0.3857 Mode :numeric Shunra East:18
Mean : 5.625 Mean : 8.25 Mean :1.256 Mean :6 NA NA Mean :0.8821 Mean : 0.3464 NA NA
3rd Qu.: 6.000 3rd Qu.:10.00 3rd Qu.:1.407 3rd Qu.:8 NA NA 3rd Qu.:0.9887 3rd Qu.: 0.5380 NA NA
Max. :10.000 Max. :17.00 Max. :2.154 Max. :8 NA NA Max. :0.9987 Max. : 0.8014 NA NA
agriculture dunes year_ct cos_td_rad sin_td_rad timediff_Jun21 site
agriculture 1.0000000 0.0000000 0.0000000 0.1724342 -0.1658646 -0.1712629 -0.5222330
dunes 0.0000000 1.0000000 0.0000000 -0.0141948 -0.0206225 -0.0153811 0.0000000
year_ct 0.0000000 0.0000000 1.0000000 -0.6058285 0.7336831 0.7239156 0.0000000
cos_td_rad 0.1724342 -0.0141948 -0.6058285 1.0000000 -0.9025273 -0.9307838 0.0461771
sin_td_rad -0.1658646 -0.0206225 0.7336831 -0.9025273 1.0000000 0.9973550 -0.0379353
timediff_Jun21 -0.1712629 -0.0153811 0.7239156 -0.9307838 0.9973550 1.0000000 -0.0377632
site -0.5222330 0.0000000 0.0000000 0.0461771 -0.0379353 -0.0377632 1.0000000

timediff_Jun21 is nearly completely correlated with sin_td_rad. sin in high correlation with cos. 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.300846684351762"

Overdispersion parameter is <1 (underdispersion) and therefore Poisson is preferable to negative binomial. Poisson family for glm / glme.

m0 <- glm(data = P.anal, formula = richness ~ agriculture * year_ct + dunes  * year_ct + cos_td_rad + sin_td_rad, family = poisson)
me0 <- glmer(data = P.anal, formula =  richness ~ agriculture * year_ct + dunes  * year_ct + cos_td_rad + sin_td_rad + (1|site), family = poisson)
## boundary (singular) fit: see help('isSingular')

Mixed model generates a singular fit. Try removing agriculture Near (exists only in Beer Milka):

me0 <- glmer(data = P.anal[agriculture=="Far"], formula =  richness ~ dunes  * year_ct + cos_td_rad + sin_td_rad + (1|site), family = poisson)
## boundary (singular) fit: see help('isSingular')

Still singular fit. Hence choose regular glm.

## Start:  AIC=291.35
## richness ~ agriculture * year_ct + dunes * year_ct + cos_td_rad + 
##     sin_td_rad
## 
##                       Df Deviance    AIC
## - year_ct:dunes        1   19.317 289.41
## - agriculture:year_ct  1   19.645 289.74
## <none>                     19.254 291.35
## - cos_td_rad           1   23.850 293.95
## - sin_td_rad           1   24.018 294.12
## 
## Step:  AIC=289.41
## richness ~ agriculture + year_ct + dunes + cos_td_rad + sin_td_rad + 
##     agriculture:year_ct
## 
##                       Df Deviance    AIC
## - dunes                1   19.341 287.44
## - agriculture:year_ct  1   19.718 287.81
## <none>                     19.317 289.41
## - cos_td_rad           1   23.896 291.99
## - sin_td_rad           1   24.049 292.15
## 
## Step:  AIC=287.44
## richness ~ agriculture + year_ct + cos_td_rad + sin_td_rad + 
##     agriculture:year_ct
## 
##                       Df Deviance    AIC
## - agriculture:year_ct  1   19.748 285.85
## <none>                     19.341 287.44
## - cos_td_rad           1   23.897 289.99
## - sin_td_rad           1   24.051 290.15
## 
## Step:  AIC=285.84
## richness ~ agriculture + year_ct + cos_td_rad + sin_td_rad
## 
##               Df Deviance    AIC
## - agriculture  1   20.310 284.41
## <none>             19.748 285.85
## - year_ct      1   24.655 288.75
## - sin_td_rad   1   25.237 289.33
## - cos_td_rad   1   25.404 289.50
## 
## Step:  AIC=284.41
## richness ~ year_ct + cos_td_rad + sin_td_rad
## 
##              Df Deviance    AIC
## <none>            20.310 284.41
## - year_ct     1   24.783 286.88
## - sin_td_rad  1   25.463 287.56
## - cos_td_rad  1   26.049 288.15

year effect remains with the time of year (both sine and cosine).

Observations 72
Dependent variable richness
Type Generalized linear model
Family poisson
Link log
χ²(3) 10.010
Pseudo-R² (Cragg-Uhler) 0.132
Pseudo-R² (McFadden) 0.035
AIC 284.407
BIC 293.514
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.957 0.146 6.276 -0.045 0.964
year_ct 0.908 0.829 0.994 -2.083 0.037
cos_td_rad 9.262 1.447 59.297 2.350 0.019
sin_td_rad 2.951 1.148 7.590 2.246 0.025
Standard errors: MLE

On average, there is a decrease of 9% every year in richness. However this effect is not statistically strong (though significant at \(\alpha=0.05\)).

geometric mean of abundance

Explore data.

## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"

## [1] "only Beer Milka data:"

Fit glm, compare gamma, gaussian (poisson inappropriate because response is not discrete)

Seems that the gamma fits the data better. Fit fixed and mixed models.

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0186238 (tol = 0.002, component 1)

Gamma mixed model fails to converge during backward stepwise selection. Try Gaussian mixed model.

## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ agriculture * year_ct + dunes * year_ct + cos_td_rad +  
##     sin_td_rad + (1 | site)
##    Data: P.anal
## 
## REML criterion at convergence: 11.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2969 -0.4887 -0.1479  0.4890  3.1882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.03597  0.1897  
##  Residual             0.04622  0.2150  
## Number of obs: 72, groups:  site, 3
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)              1.71554    0.54523   3.146
## agricultureNear          0.10076    0.23455   0.430
## year_ct                 -0.08396    0.03150  -2.665
## dunesshifting           -0.26917    0.19334  -1.392
## cos_td_rad              -0.03629    0.49027  -0.074
## sin_td_rad               0.30692    0.25311   1.213
## agricultureNear:year_ct -0.03138    0.03725  -0.842
## year_ct:dunesshifting    0.03613    0.03106   1.163
## 
## Correlation of Fixed Effects:
##             (Intr) agrclN yer_ct dnsshf cs_td_ sn_td_ agrN:_
## agricultrNr -0.334                                          
## year_ct     -0.315  0.245                                   
## dunesshftng -0.205  0.012  0.442                            
## cos_td_rad  -0.935  0.259  0.046  0.037                     
## sin_td_rad  -0.644  0.225 -0.386  0.061  0.767              
## agrcltrNr:_  0.352 -0.951 -0.299 -0.009 -0.273 -0.197       
## yr_ct:dnssh  0.181 -0.006 -0.465 -0.965 -0.016 -0.040  0.004
## $site
##             (Intercept)
## Beer Milka   0.12916852
## Secher      -0.20711889
## Shunra East  0.07795038
## 
## with conditional variances for "site"
## $site

cos_td_rad is highly correlated with the model intercept. Center the variable to reduce this correlation:

perform stepwise model selection of Gaussian mixed model.

## Single term deletions
## 
## Model:
## gma ~ agriculture * year_ct + dunes * year_ct + cos_td_rad_c + 
##     sin_td_rad + (1 | site)
##                     npar    AIC
## <none>                   3.1117
## cos_td_rad_c           1 1.1122
## sin_td_rad             1 2.6513
## agriculture:year_ct    1 1.9251
## year_ct:dunes          1 2.6066

drop cosine of time difference phase

## Single term deletions
## 
## Model:
## gma ~ agriculture * year_ct + dunes * year_ct + sin_td_rad + 
##     (1 | site)
##                     npar     AIC
## <none>                   1.11218
## sin_td_rad             1 2.60193
## agriculture:year_ct    1 0.00261
## year_ct:dunes          1 0.60658

drop interaction of agriculture and year.

## Single term deletions
## 
## Model:
## gma ~ agriculture + dunes * year_ct + sin_td_rad + (1 | site)
##               npar      AIC
## <none>              0.00261
## agriculture      1 -0.52154
## sin_td_rad       1  1.49997
## dunes:year_ct    1 -0.52329

drop agriculture.

## Single term deletions
## 
## Model:
## gma ~ dunes * year_ct + sin_td_rad + (1 | site)
##               npar      AIC
## <none>             -0.52154
## sin_td_rad       1  1.47983
## dunes:year_ct    1 -1.10351

Drop interaction of dunes and year.

## Single term deletions
## 
## Model:
## gma ~ dunes + year_ct + sin_td_rad + (1 | site)
##            npar     AIC
## <none>          -1.1035
## dunes         1 -2.0192
## year_ct       1  4.6232
## sin_td_rad    1  0.9842

drop dunes.

## Single term deletions
## 
## Model:
## gma ~ year_ct + sin_td_rad + (1 | site)
##            npar     AIC
## <none>          -2.0192
## year_ct       1  3.7374
## sin_td_rad    1  0.1386

This is the final model:

## Linear mixed model fit by REML ['lmerMod']
## Formula: gma ~ year_ct + sin_td_rad + (1 | site)
##    Data: P.anal
## 
## REML criterion at convergence: -0.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8885 -0.6476 -0.2118  0.4196  3.4895 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 0.03259  0.1805  
##  Residual             0.04608  0.2147  
## Number of obs: 72, groups:  site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  1.59851    0.15723  10.167
## year_ct     -0.08049    0.02584  -3.115
## sin_td_rad   0.37229    0.15842   2.350
## 
## Correlation of Fixed Effects:
##            (Intr) yer_ct
## year_ct    -0.696       
## sin_td_rad  0.427 -0.800
## $site
##             (Intercept)
## Beer Milka   0.12916852
## Secher      -0.20711889
## Shunra East  0.07795038
## 
## with conditional variances for "site"
## $site

## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
Observations 72
Dependent variable gma
Type Mixed effects linear regression
AIC 9.235
BIC 20.618
Pseudo-R² (fixed effects) 0.095
Pseudo-R² (total) 0.470
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 1.599 0.158 10.098 8.250 0.000
year_ct -0.080 0.026 -3.038 68.975 0.003
sin_td_rad 0.372 0.165 2.261 68.690 0.027
p values calculated using Kenward-Roger standard errors and d.f.
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.181
Residual 0.215
Grouping Variables
Group # groups ICC
site 3 0.414 
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.

year and time of year are significant. There is a significant decrease in mean abundance over time.

abundance

Explore data

## [1] "ABUNDANCE WITHOUT RARE SPECIES"

## [1] "only Beer Milka data:"

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 = 1.21897720119714"

Overdispersion parameter is >1 (overdispersion) and therefore negative binomial is preferable to poisson. Fit fixed and mixed models. Choose mixed model if possible, otherwise choose a model with fixed-effects only. add interaction between dunes and agriculture based on boxplots.

## theta.ml: iter 0 'theta = 10.967700'
## theta.ml: iter1 theta =16.7431
## theta.ml: iter2 theta =25.0887
## theta.ml: iter3 theta =37.1591
## theta.ml: iter4 theta =54.7096
## theta.ml: iter5 theta =80.4105
## theta.ml: iter6 theta =118.315
## theta.ml: iter7 theta =174.551
## theta.ml: iter8 theta =258.352
## theta.ml: iter9 theta =383.591
## theta.ml: iter10 theta =571.084
## theta.ml: iter11 theta =852.049
## theta.ml: iter12 theta =1273.29
## theta.ml: iter13 theta =1905.02
## theta.ml: iter14 theta =2852.51
## theta.ml: iter15 theta =4273.68
## theta.ml: iter16 theta =6405.39
## theta.ml: iter17 theta =9602.92
## theta.ml: iter18 theta =14399.2
## theta.ml: iter19 theta =21593.6
## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached
## th := est_theta(glmer(..)) = 21593.59
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00306768 (tol = 0.002, component 1)
##  --> dev.= -2*logLik(.) = 336.2834 
##  1: th=    10635.44843, dev=  336.28918812, beta[1]=    1.67154751
##  2: th=    43842.36020, dev=  336.28061437, beta[1]=    1.67153304
##  3: th=    105213.1908, dev=  336.27901436, beta[1]=    1.67161098
##  4: th=    111494.2445, dev=  336.27894999, beta[1]=    1.67164257
##  5: th=    187324.9118, dev=  336.27851348, beta[1]=    1.67164971
##  6: th=    356495.0586, dev=  336.27820897, beta[1]=    1.67166007
##  7: th=    278811.1036, dev=  336.27830288, beta[1]=    1.67165922
##  8: th=    324549.6499, dev=  336.27824212, beta[1]=    1.67163688
##  9: th=    384219.8098, dev=  336.27818467, beta[1]=    1.67166185
## 10: th=    402422.3319, dev=  336.27817050, beta[1]=    1.67168007
## 11: th=    414100.7250, dev=  336.27816201, beta[1]=    1.67166419
## 12: th=    421487.2038, dev=  336.27815688, beta[1]=    1.67161608
## 13: th=    426118.0278, dev=  336.27815387, beta[1]=    1.67164422
## 14: th=    429005.4381, dev=  336.27815196, beta[1]=    1.67163422
## 15: th=    430799.7301, dev=  336.27815077, beta[1]=    1.67164320
## 16: th=    431912.4137, dev=  336.27815006, beta[1]=    1.67161195
## 17: th=    432601.5265, dev=  336.27814958, beta[1]=    1.67162341
## 18: th=    433027.9713, dev=  336.27814938, beta[1]=    1.67165261
## 19: th=    433291.7388, dev=  336.27814924, beta[1]=    1.67161828
## 20: th=    433454.8364, dev=  336.27814908, beta[1]=    1.67162529
## 21: th=    433555.6670, dev=  336.27814889, beta[1]=    1.67165462
## 22: th=    433617.9954, dev=  336.27814900, beta[1]=    1.67163033
## 23: th=    433547.8648, dev=  336.27814900, beta[1]=    1.67168292
## 24: th=    433582.9616, dev=  336.27814902, beta[1]=    1.67171308
## 25: th=    433566.0924, dev=  336.27814900, beta[1]=    1.67170130
## 26: th=    433555.6670, dev=  336.27814890, beta[1]=    1.67163179
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(433555.7)  ( log )
## Formula: 
## abundance ~ agriculture * year_ct + dunes * year_ct + agriculture:dunes +  
##     cos_td_rad + sin_td_rad + (1 | site)
##    Data: P.anal
##       AIC       BIC    logLik  deviance  df.resid 
##  358.2781  383.3215 -168.1391  336.2781        61 
## Random effects:
##  Groups Name        Std.Dev.
##  site   (Intercept) 0.283   
## Number of obs: 72, groups:  site, 3
## Fixed Effects:
##                   (Intercept)                agricultureNear  
##                       1.67163                       -0.43620  
##                       year_ct                  dunesshifting  
##                      -0.21923                       -0.27299  
##                    cos_td_rad                     sin_td_rad  
##                       1.38201                        1.32464  
##       agricultureNear:year_ct          year_ct:dunesshifting  
##                       0.03262                        0.02924  
## agricultureNear:dunesshifting  
##                       0.38837

theta is very large and the model failed to converge –> second term in variance formula tends to zero –> probably model is underdispersed and poisson is better.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## abundance ~ agriculture * year_ct + dunes * year_ct + agriculture:dunes +  
##     cos_td_rad + sin_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    356.3    379.0   -168.1    336.3       62 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9599 -0.5864 -0.1080  0.3654  1.9584 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.08008  0.283   
## Number of obs: 72, groups:  site, 3
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.67159    0.92723   1.803  0.07142 .  
## agricultureNear               -0.43619    0.39561  -1.103  0.27021    
## year_ct                       -0.21922    0.05580  -3.929 8.53e-05 ***
## dunesshifting                 -0.27295    0.32725  -0.834  0.40424    
## cos_td_rad                     1.38203    0.85038   1.625  0.10412    
## sin_td_rad                     1.32464    0.41478   3.194  0.00141 ** 
## agricultureNear:year_ct        0.03262    0.06338   0.515  0.60682    
## year_ct:dunesshifting          0.02924    0.05431   0.538  0.59036    
## agricultureNear:dunesshifting  0.38838    0.19480   1.994  0.04618 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) agrclN yer_ct dnsshf cs_td_ sn_td_ agrN:_ yr_ct:
## agricultrNr -0.245                                                 
## year_ct     -0.315  0.233                                          
## dunesshftng -0.163 -0.022  0.414                                   
## cos_td_rad  -0.940  0.165  0.048  0.002                            
## sin_td_rad  -0.606  0.128 -0.427  0.026  0.735                     
## agrcltrNr:_  0.287 -0.922 -0.275  0.056 -0.210 -0.144              
## yr_ct:dnssh  0.156  0.057 -0.437 -0.954 -0.002 -0.027 -0.055       
## agrcltrNr:d -0.103 -0.195 -0.012 -0.114  0.136  0.141 -0.067 -0.038

cos_td_rad is highly correlated with intercept. use centered variable. same for year_ct.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## abundance ~ agriculture * year_ct_c + dunes * year_ct_c + agriculture:dunes +  
##     cos_td_rad_c + sin_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    356.3    379.0   -168.1    336.3       62 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9600 -0.5864 -0.1080  0.3654  1.9584 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.08009  0.283   
## Number of obs: 72, groups:  site, 3
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.57535    0.23423   6.726 1.75e-11 ***
## agricultureNear               -0.24049    0.15399  -1.562  0.11836    
## year_ct_c                     -0.21923    0.05580  -3.929 8.52e-05 ***
## dunesshifting                 -0.09753    0.09856  -0.990  0.32238    
## cos_td_rad_c                   1.38201    0.85038   1.625  0.10413    
## sin_td_rad                     1.32465    0.41478   3.194  0.00141 ** 
## agricultureNear:year_ct_c      0.03262    0.06338   0.515  0.60679    
## year_ct_c:dunesshifting        0.02924    0.05431   0.538  0.59027    
## agricultureNear:dunesshifting  0.38838    0.19480   1.994  0.04618 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) agrclN yr_ct_ dnsshf cs_t__ sn_td_ agN:__ yr_c_:
## agricultrNr -0.105                                                 
## year_ct_c    0.336 -0.081                                          
## dunesshftng -0.198  0.300 -0.073                                   
## cos_td_rd_c -0.452 -0.095  0.048  0.000                            
## sin_td_rad  -0.652 -0.026 -0.427 -0.004  0.735                     
## agrcltrN:__  0.070  0.101 -0.275  0.004 -0.210 -0.144              
## yr_ct_c:dns -0.013  0.009 -0.437  0.137 -0.002 -0.027 -0.055       
## agrcltrNr:d  0.010 -0.668 -0.012 -0.502  0.136  0.141 -0.067 -0.038

Perform stepwise model selection of mixed model.

## Single term deletions
## 
## Model:
## abundance ~ agriculture * year_ct_c + dunes * year_ct_c + agriculture:dunes + 
##     cos_td_rad_c + sin_td_rad + (1 | site)
##                       npar    AIC
## <none>                     356.28
## cos_td_rad_c             1 356.97
## sin_td_rad               1 364.53
## agriculture:year_ct_c    1 354.54
## year_ct_c:dunes          1 354.57
## agriculture:dunes        1 358.26

drop interaction of agriculture and year.

## Single term deletions
## 
## Model:
## abundance ~ agriculture + dunes * year_ct_c + agriculture:dunes + 
##     cos_td_rad_c + sin_td_rad + (1 | site)
##                   npar    AIC
## <none>                 354.54
## cos_td_rad_c         1 355.75
## sin_td_rad           1 363.59
## dunes:year_ct_c      1 352.86
## agriculture:dunes    1 356.70

drop interaction of dunes and year.

## Single term deletions
## 
## Model:
## abundance ~ agriculture + dunes + year_ct_c + agriculture:dunes + 
##     cos_td_rad_c + sin_td_rad + (1 | site)
##                   npar    AIC
## <none>                 352.86
## year_ct_c            1 369.30
## cos_td_rad_c         1 354.10
## sin_td_rad           1 362.12
## agriculture:dunes    1 355.07

drop cosine of sampling time (because \(\Delta AIC<2\) and should prefer simpler model in such a case)

## Single term deletions
## 
## Model:
## abundance ~ agriculture + dunes + year_ct_c + agriculture:dunes + 
##     sin_td_rad + (1 | site)
##                   npar    AIC
## <none>                 354.10
## year_ct_c            1 369.97
## sin_td_rad           1 360.91
## agriculture:dunes    1 355.41

drop interaction of agri and dunes (because \(\Delta AIC<2\) and should prefer simpler model in such a case)

## Single term deletions
## 
## Model:
## abundance ~ agriculture + dunes + year_ct_c + sin_td_rad + (1 | 
##     site)
##             npar    AIC
## <none>           355.41
## agriculture    1 353.60
## dunes          1 353.44
## year_ct_c      1 370.44
## sin_td_rad     1 361.51

drop dunes.

## Single term deletions
## 
## Model:
## abundance ~ agriculture + year_ct_c + sin_td_rad + (1 | site)
##             npar    AIC
## <none>           353.44
## agriculture    1 351.63
## year_ct_c      1 368.54
## sin_td_rad     1 359.59

drop agri.

## Single term deletions
## 
## Model:
## abundance ~ year_ct_c + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          351.63
## year_ct_c     1 367.45
## sin_td_rad    1 358.36

The final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: abundance ~ year_ct_c + sin_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    351.6    360.7   -171.8    343.6       68 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6573 -0.6922 -0.2195  0.3943  2.6207 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.09452  0.3074  
## Number of obs: 72, groups:  site, 3
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.69462    0.21291   7.959 1.73e-15 ***
## year_ct_c   -0.19423    0.04748  -4.091 4.29e-05 ***
## sin_td_rad   0.82584    0.28005   2.949  0.00319 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) yr_ct_
## year_ct_c   0.433       
## sin_td_rad -0.506 -0.826
## $site
##             (Intercept)
## Beer Milka    0.1947242
## Secher       -0.4162581
## Shunra East   0.2299138
## 
## with conditional variances for "site"
## $site

Interpretation of abundance model:

Observations 72
Dependent variable abundance
Type Mixed effects generalized linear model
Family poisson
Link log
AIC 351.632
BIC 360.739
Pseudo-R² (fixed effects) 0.177
Pseudo-R² (total) 0.531
Fixed Effects
exp(Est.) S.E. z val. p
(Intercept) 5.445 0.213 7.959 0.000
year_ct_c 0.823 0.047 -4.091 0.000
sin_td_rad 2.284 0.280 2.949 0.003
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.307
Grouping Variables
Group # groups ICC
site 3 0.086 
## Loading required namespace: broom.mixed

Year and time of year are significant.

## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.

Abundance data: There is a highly significant negative overall temporal trend (decrease of 18% per year). Significant effect also for time of year.

community analysis using package MVabund

##       nb       po 
## 142.4615 143.0549
## [1] "POISSON"

## [1] "NEGATIVE BINOMIAL"

prefer negative binomial because standardized residuals are smaller.

##       nb       po      nb1 
## 142.4615 143.0549 145.6183

including site as an explanatory variable does not improve the AIC of the model. stepwise selection of model:

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture * year_ct + dunes * year_ct + agriculture:dunes + 
##     cos_td_rad + sin_td_rad
##                     Df    AIC
## <none>                 1567.1
## cos_td_rad          11 1569.1
## sin_td_rad          11 1565.2
## agriculture:year_ct 11 1555.5
## year_ct:dunes       11 1556.0
## agriculture:dunes   11 1562.8

drop interaction of agriculture with year.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes * year_ct + agriculture:dunes + 
##     cos_td_rad + sin_td_rad
##                   Df    AIC
## <none>               1555.5
## cos_td_rad        11 1555.5
## sin_td_rad        11 1551.7
## dunes:year_ct     11 1543.1
## agriculture:dunes 11 1549.6

drop interaction of dunes with year.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes + year_ct + agriculture:dunes + 
##     cos_td_rad + sin_td_rad
##                   Df    AIC
## <none>               1543.1
## year_ct           11 1535.2
## cos_td_rad        11 1541.7
## sin_td_rad        11 1538.0
## agriculture:dunes 11 1536.9

drop year.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes + agriculture:dunes + cos_td_rad + 
##     sin_td_rad
##                   Df    AIC
## <none>               1535.2
## cos_td_rad        11 1531.6
## sin_td_rad        11 1528.6
## agriculture:dunes 11 1527.6

drop interaction of agri and dunes.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes + cos_td_rad + sin_td_rad
##             Df    AIC
## <none>         1527.6
## agriculture 11 1528.7
## dunes       11 1556.8
## cos_td_rad  11 1523.6
## sin_td_rad  11 1521.2

drop sine of radian time difference from June 21st.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes + cos_td_rad
##             Df    AIC
## <none>         1521.2
## agriculture 11 1521.1
## dunes       11 1547.5
## cos_td_rad  11 1511.3

drop cosine.

## Single term deletions
## 
## Model:
## spp_no_rare ~ agriculture + dunes
##             Df    AIC
## <none>         1511.3
## agriculture 11 1511.5
## dunes       11 1538.2

final model includes agriculture and dunes (keep agriculture even though \(\Delta AIC<2\) because agriculture comes out significant and is also significant for single species).

## 
## Test statistics:
##                 wald value Pr(>wald)    
## (Intercept)         11.590     0.001 ***
## agricultureNear      4.117     0.009 ** 
## dunesshifting        7.127     0.001 ***
## --- 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Test statistic:  8.103, 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 0 min 22 sec
## Analysis of Deviance Table
## 
## Model: spp_no_rare ~ agriculture + dunes
## 
## Multivariate test:
##             Res.Df Df.diff   Dev Pr(>Dev)    
## (Intercept)     71                           
## agriculture     70       1 17.45    0.022 *  
## dunes           69       1 48.87    0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Univariate Tests:
##             Acanthodactylus.aegyptius          Acanthodactylus.scutellatus
##                                   Dev Pr(>Dev)                         Dev
## (Intercept)                                                               
## agriculture                     4.902    0.138                       2.001
## dunes                          14.815    0.002                      18.398
##                      Cerastes.cerastes          Cerastes.vipera         
##             Pr(>Dev)               Dev Pr(>Dev)             Dev Pr(>Dev)
## (Intercept)                                                             
## agriculture    0.562             0.521    0.880           1.913    0.562
## dunes          0.001             0.067    0.866           0.931    0.688
##             Chalcides.sepsoides          Chamaeleo.chamaeleon         
##                             Dev Pr(>Dev)                  Dev Pr(>Dev)
## (Intercept)                                                           
## agriculture               0.318    0.880                0.128    0.880
## dunes                     0.301    0.866                1.646    0.547
##             Lytorhynchus.diadema          Scincus.scincus         
##                              Dev Pr(>Dev)             Dev Pr(>Dev)
## (Intercept)                                                       
## agriculture                0.189    0.880           0.938    0.779
## dunes                      0.209    0.866           8.878    0.011
##             Spalerosophis.diadema          Stenodactylus.petrii         
##                               Dev Pr(>Dev)                  Dev Pr(>Dev)
## (Intercept)                                                             
## agriculture                 5.754    0.089                0.246    0.880
## dunes                       0.403    0.866                0.385    0.866
##             Varanus.griseus         
##                         Dev Pr(>Dev)
## (Intercept)                         
## agriculture           0.537    0.880
## dunes                  2.84    0.307
## Arguments:
##  Test statistics calculated assuming uncorrelated response (for faster computation) 
## P-value calculated using 999 iterations via PIT-trap resampling.

Dunes factor has a statistically significant effect on community composition. no significant temporal trend.

## 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.

coefficients are exponentiated because the link function used for negative binomial is log -> above 1 is positive and below 1 is negative.