## 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
## Loading required package: solartime
## - Winter data excluded.
## - The counts of all passing / wintering species (not resident nor breeding in the sampled unit) were set to zero. Abreed_and_non_breed contains the non-breeders as well.

Planted Conifer Forest

Unit is divided into 3 subunits: Judea, Carmel and Galilee. Only factors is time. All plots are located far from settlements. Sampling started in spring 2014. Total 5 campaigns, 3 subunits per campaign, 5 sites per subunit, with 3 plots per site. Total of 45 plots per campaign, and 225 plot-campaign combinations.

Raw data Total abundance: 6230 Number of observations: 2161 Total richness: 86

Filtered data Total abundance: 4150 Number of observations: 1631 Total richness: 43

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

richness

Explore data.

## [1] "RICHNESS WITH RARE SPECIES"
## ℹ SHA-1 hash of file is "eef55df1963ba201f57cf44513f566298756fd4f"
richness year_ct site subunit td_sc cos_td_rad sin_td_rad h_from_sunrise cos_hsun sin_hsun monitors_name wind precipitation temperature clouds
Min. : 1.000 Min. :2.0 Aderet : 15 Judean Highlands:75 Min. :-2.3697 Min. :-0.3335 Min. :-0.9959 Length:225 Min. :-0.2639 Min. :-0.1201 Eyal Shochat: 35 0 : 24 0 : 56 0 : 0 0 : 9
1st Qu.: 5.000 1st Qu.:3.0 Amatzia : 15 Carmel :75 1st Qu.:-0.2713 1st Qu.: 0.5702 1st Qu.:-0.8215 Class :difftime 1st Qu.: 0.6155 1st Qu.: 0.2784 Asaf Mayrose: 23 1 : 13 3 : 3 1 : 20 1 : 3
Median : 6.000 Median :5.0 Bat Shlomo: 15 Galilee :75 Median : 0.3010 Median : 0.7611 Median :-0.6486 Mode :numeric Median : 0.8098 Median : 0.5867 Sassi Haham : 18 2 : 7 NA’s:166 2 : 26 2 : 1
Mean : 6.022 Mean :5.2 Eitanim : 15 NA Mean : 0.2127 Mean : 0.6933 Mean :-0.6326 NA Mean : 0.7332 Mean : 0.5336 Eliraz Dvir : 12 3 : 1 NA 3 : 1 3 : 2
3rd Qu.: 7.000 3rd Qu.:7.0 Elyakim : 15 NA 3rd Qu.: 0.6444 3rd Qu.: 0.8521 3rd Qu.:-0.5234 NA 3rd Qu.: 0.9605 3rd Qu.: 0.7881 Eran Banker : 9 NA’s:180 NA NA’s:178 NA’s:210
Max. :13.000 Max. :9.0 Eshtaol : 15 NA Max. : 1.6364 Max. : 0.9947 Max. :-0.1031 NA Max. : 1.0000 Max. : 1.0000 (Other) : 20 NA NA NA NA
NA NA (Other) :135 NA NA NA NA NA NA’s :45 NA’s :45 NA’s :108 NA NA NA NA

no observation for all 4 weather variables. many NAs for sampling time of day variables.exclude from model. An extreme observation of richness>20 in Judean highlands near settlements.

richness year_ct site subunit td_sc cos_td_rad sin_td_rad h_from_sunrise cos_hsun sin_hsun
richness 1.0000000 -0.0869667 -0.1526824 -0.2098492 0.2468122 0.2620866 0.2032932 -0.1291245 0.1729942 -0.0748826
year_ct -0.0869667 1.0000000 0.0000000 0.0000000 -0.4304647 -0.4084156 -0.4336689 -0.0881196 0.0985736 -0.0618585
site -0.1526824 0.0000000 1.0000000 0.7181325 -0.1237969 -0.0708077 -0.2124888 0.2856522 -0.2977022 0.2576585
subunit -0.2098492 0.0000000 0.7181325 1.0000000 -0.1180592 -0.0606624 -0.2128418 0.2944203 -0.3327003 0.2303051
td_sc 0.2468122 -0.4304647 -0.1237969 -0.1180592 1.0000000 0.9753453 0.9307073 0.0619157 -0.0348851 0.0841439
cos_td_rad 0.2620866 -0.4084156 -0.0708077 -0.0606624 0.9753453 1.0000000 0.8324627 0.1128374 -0.0873088 0.1309387
sin_td_rad 0.2032932 -0.4336689 -0.2124888 -0.2128418 0.9307073 0.8324627 1.0000000 -0.0274966 0.0538205 -0.0012408
h_from_sunrise -0.1291245 -0.0881196 0.2856522 0.2944203 0.0619157 0.1128374 -0.0274966 1.0000000 -0.9603558 0.9618286
cos_hsun 0.1729942 0.0985736 -0.2977022 -0.3327003 -0.0348851 -0.0873088 0.0538205 -0.9603558 1.0000000 -0.8509778
sin_hsun -0.0748826 -0.0618585 0.2576585 0.2303051 0.0841439 0.1309387 -0.0012408 0.9618286 -0.8509778 1.0000000

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.658637694737327"

Overdispersion parameter is < 1. Choose Poisson.

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

go with mixed model with subunits as a fixed effect (want to make statement about subunits), attempt model selection yet.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ subunit * year_ct + cos_td_rad + sin_td_rad + (1 |  
##     site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    965.2    995.9   -473.6    947.2      216 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7752 -0.5322 -0.0193  0.4039  2.7190 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01239  0.1113  
## Number of obs: 225, groups:  site, 15
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)   
## (Intercept)             1.04905    0.33058   3.173  0.00151 **
## subunitCarmel           0.13102    0.17236   0.760  0.44718   
## subunitGalilee         -0.07910    0.17679  -0.447  0.65459   
## year_ct                 0.01929    0.02006   0.961  0.33638   
## cos_td_rad              0.67536    0.24778   2.726  0.00642 **
## sin_td_rad             -0.39325    0.28670  -1.372  0.17017   
## subunitCarmel:year_ct  -0.03749    0.02601  -1.441  0.14950   
## subunitGalilee:year_ct -0.02178    0.02704  -0.805  0.42054   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sbntCr sbntGl yer_ct cs_td_ sn_td_ sbnC:_
## subunitCrml -0.029                                          
## subunitGall -0.026  0.548                                   
## year_ct     -0.148  0.622  0.603                            
## cos_td_rad  -0.933 -0.180 -0.179 -0.050                     
## sin_td_rad   0.841  0.317  0.307  0.300 -0.844              
## sbntCrml:y_  0.131 -0.808 -0.426 -0.708  0.034 -0.158       
## sbntGll:yr_  0.045 -0.439 -0.823 -0.688  0.120 -0.227  0.499
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00332887 (tol = 0.002, component 1)

perform stepwise model selection of poisson mixed model.

## Single term deletions
## 
## Model:
## richness ~ subunit * year_ct + cos_td_rad + sin_td_rad + (1 | 
##     site)
##                 npar    AIC
## <none>               965.18
## cos_td_rad         1 971.36
## sin_td_rad         1 965.10
## subunit:year_ct    2 963.27

remove subunit X year.

## Single term deletions
## 
## Model:
## richness ~ subunit + year_ct + cos_td_rad + sin_td_rad + (1 | 
##     site)
##            npar    AIC
## <none>          963.27
## subunit       2 962.88
## year_ct       1 961.30
## cos_td_rad    1 969.76
## sin_td_rad    1 964.06

drop year.

## Single term deletions
## 
## Model:
## richness ~ subunit + cos_td_rad + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          961.30
## subunit       2 960.89
## cos_td_rad    1 967.83
## sin_td_rad    1 962.09

drop subunit.

## Single term deletions
## 
## Model:
## richness ~ cos_td_rad + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          960.89
## cos_td_rad    1 967.34
## sin_td_rad    1 961.41

drop sine.

## Single term deletions
## 
## Model:
## richness ~ cos_td_rad + (1 | site)
##            npar    AIC
## <none>          961.41
## cos_td_rad    1 968.85

only time of year remains. Final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ cos_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    961.4    971.7   -477.7    955.4      222 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.81116 -0.55209 -0.05948  0.41685  2.46546 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01481  0.1217  
## Number of obs: 225, groups:  site, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.54330    0.09217  16.745  < 2e-16 ***
## cos_td_rad   0.34666    0.11534   3.005  0.00265 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## cos_td_rad -0.892
## $site
##                  (Intercept)
## Aderet           0.041958260
## Amatzia          0.194200215
## Bat Shlomo       0.050085708
## Eitanim         -0.057342406
## Elyakim          0.017471976
## Eshtaol         -0.142787807
## Givat Yeshayahu  0.134084761
## Kabri           -0.117575421
## Kerem Maharal    0.066797701
## Manara          -0.094902416
## Meron           -0.110849885
## Ofer            -0.009267953
## Ramat Hashofet   0.008894432
## Ramot Naftali    0.055138130
## Zuriel          -0.008776173
## 
## with conditional variances for "site"
## $site

attempt to add year, see if significant:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: richness ~ year_ct + cos_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    963.4    977.0   -477.7    955.4      221 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.80811 -0.54823 -0.05958  0.40392  2.42124 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.01476  0.1215  
## Number of obs: 225, groups:  site, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.522036   0.133873  11.369  < 2e-16 ***
## year_ct     0.002545   0.011607   0.219  0.82644    
## cos_td_rad  0.358225   0.126935   2.822  0.00477 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) yer_ct
## year_ct    -0.725       
## cos_td_rad -0.860  0.416

year not significant, rightfully dropped.

No statistically significant effects of subunit or time found for richness.

geometric mean of abundance

Explore data. Exclude time of day because of high number of NAs.

## [1] "GEOMETRIC MEAN ABUNDANCE WITHOUT RARE SPECIES"
gma year_ct site settlements subunit td_sc cos_td_rad sin_td_rad
Min. :1.000 Min. :2.0 Aderet : 15 Far : 9 Judean Highlands:75 Min. :-2.3697 Min. :-0.3335 Min. :-0.9959
1st Qu.:1.861 1st Qu.:3.0 Amatzia : 15 Near: 0 Carmel :75 1st Qu.:-0.2713 1st Qu.: 0.5702 1st Qu.:-0.8215
Median :2.330 Median :5.0 Bat Shlomo: 15 NA’s:216 Galilee :75 Median : 0.3010 Median : 0.7611 Median :-0.6486
Mean :2.574 Mean :5.2 Eitanim : 15 NA NA Mean : 0.2127 Mean : 0.6933 Mean :-0.6326
3rd Qu.:3.061 3rd Qu.:7.0 Elyakim : 15 NA NA 3rd Qu.: 0.6444 3rd Qu.: 0.8521 3rd Qu.:-0.5234
Max. :6.327 Max. :9.0 Eshtaol : 15 NA NA Max. : 1.6364 Max. : 0.9947 Max. :-0.1031
NA NA (Other) :135 NA NA NA NA NA

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

Remove rows 67, 112.

Fit fixed and mixed models.

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

Mixed model did not converge, attempt model selection.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: gma ~ subunit * year_ct + cos_td_rad + sin_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    564.6    598.7   -272.3    544.6      213 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6655 -0.6553 -0.1655  0.5333  3.0425 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept) 0.0003395 0.01843 
##  Residual             0.1162229 0.34091 
## Number of obs: 223, groups:  site, 15
## 
## Fixed effects:
##                         Estimate Std. Error t value Pr(>|z|)    
## (Intercept)             0.539774   0.111683   4.833 1.34e-06 ***
## subunitCarmel          -0.016579   0.050921  -0.326  0.74474    
## subunitGalilee          0.157705   0.059452   2.653  0.00799 ** 
## year_ct                 0.011174   0.006922   1.614  0.10649    
## cos_td_rad             -0.146400   0.083722  -1.749  0.08035 .  
## sin_td_rad              0.159644   0.094525   1.689  0.09124 .  
## subunitCarmel:year_ct  -0.006663   0.008195  -0.813  0.41616    
## subunitGalilee:year_ct -0.016063   0.009618  -1.670  0.09489 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sbntCr sbntGl yer_ct cs_td_ sn_td_ sbnC:_
## subunitCrml -0.024                                          
## subunitGall -0.013  0.546                                   
## year_ct     -0.195  0.676  0.582                            
## cos_td_rad  -0.942 -0.174 -0.157  0.006                     
## sin_td_rad   0.850  0.329  0.291  0.245 -0.845              
## sbntCrml:y_  0.176 -0.830 -0.450 -0.780 -0.002 -0.135       
## sbntGll:yr_  0.056 -0.469 -0.863 -0.671  0.096 -0.206  0.535
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0245391 (tol = 0.002, component 1)

perform stepwise model selection of gaussian model.

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00783015 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00313316 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00368728 (tol = 0.002, component 1)
## Single term deletions
## 
## Model:
## gma ~ subunit * year_ct + cos_td_rad + sin_td_rad + (1 | site)
##                 npar    AIC
## <none>               564.62
## cos_td_rad         1 565.95
## sin_td_rad         1 565.58
## subunit:year_ct    2 563.40

drop year : subunit

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00368728 (tol = 0.002, component 1)
## Single term deletions
## 
## Model:
## gma ~ subunit + year_ct + cos_td_rad + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          563.40
## subunit       2 569.61
## year_ct       1 562.38
## cos_td_rad    1 564.21
## sin_td_rad    1 563.39

drop year

## Single term deletions
## 
## Model:
## gma ~ subunit + cos_td_rad + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          562.38
## subunit       2 568.55
## cos_td_rad    1 563.38
## sin_td_rad    1 561.91

drop sine

## Single term deletions
## 
## Model:
## gma ~ subunit + cos_td_rad + (1 | site)
##            npar    AIC
## <none>          561.91
## subunit       2 569.17
## cos_td_rad    1 561.68

drop cosine.

## Single term deletions
## 
## Model:
## gma ~ subunit + (1 | site)
##         npar    AIC
## <none>       561.68
## subunit    2 569.51

Only subunit remains. This is the final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: gma ~ subunit + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##    561.7    578.7   -275.8    551.7      218 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6633 -0.6407 -0.1850  0.5666  3.2713 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  site     (Intercept) 0.0002287 0.01512 
##  Residual             0.1224116 0.34987 
## Number of obs: 223, groups:  site, 15
## 
## Fixed effects:
##                Estimate Std. Error t value Pr(>|z|)    
## (Intercept)     0.40298    0.01839  21.915  < 2e-16 ***
## subunitCarmel  -0.06942    0.02443  -2.841  0.00449 ** 
## subunitGalilee  0.06403    0.02772   2.310  0.02088 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sbntCr
## subunitCrml -0.748       
## subunitGall -0.662  0.497
## $site
##                   (Intercept)
## Aderet           0.0014332238
## Amatzia          0.0052545903
## Bat Shlomo       0.0053741289
## Eitanim          0.0011619522
## Elyakim          0.0030863418
## Eshtaol         -0.0063275188
## Givat Yeshayahu -0.0019589870
## Kabri           -0.0087587560
## Kerem Maharal   -0.0092614993
## Manara           0.0043994590
## Meron            0.0011939807
## Ofer            -0.0017370411
## Ramat Hashofet   0.0017625159
## Ramot Naftali    0.0027354100
## Zuriel           0.0002541877
## 
## with conditional variances for "site"
## $site

## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
Observations 223
Dependent variable gma
Type Mixed effects generalized linear model
Family Gamma
Link inverse
AIC 561.68
BIC 578.72
Pseudo-R² (fixed effects) NA
Pseudo-R² (total) NA
Fixed Effects
Est. S.E. t val. p
(Intercept) 0.40 0.02 21.92 0.00
subunitCarmel -0.07 0.02 -2.84 0.00
subunitGalilee 0.06 0.03 2.31 0.02
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.02
Residual 0.35
Grouping Variables
Group # groups ICC
site 15 0.00 
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
contrast estimate SE df z.ratio p.value
Judean Highlands - Carmel 0.0694207 0.0244337 Inf 2.841183 0.0067420
Judean Highlands - Galilee -0.0640337 0.0277188 Inf -2.310115 0.0208818
Carmel - Galilee -0.1334544 0.0263182 Inf -5.070793 0.0000012

GMA is significantly different among all three subunits, with the highest being Carmel and the lowest Galilee (post-hoc z-test for pairwise contrasts, fdr correction for 3 tests). No significant temporal trend.

abundance

Explore data

## [1] "ABUNDANCE WITHOUT RARE SPECIES"
abundance year_ct site settlements subunit td_sc cos_td_rad sin_td_rad
Min. : 1.00 Min. :2.0 Aderet : 15 Far : 9 Judean Highlands:75 Min. :-2.3697 Min. :-0.3335 Min. :-0.9959
1st Qu.:10.00 1st Qu.:3.0 Amatzia : 15 Near: 0 Carmel :75 1st Qu.:-0.2713 1st Qu.: 0.5702 1st Qu.:-0.8215
Median :16.00 Median :5.0 Bat Shlomo: 15 NA’s:216 Galilee :75 Median : 0.3010 Median : 0.7611 Median :-0.6486
Mean :17.14 Mean :5.2 Eitanim : 15 NA NA Mean : 0.2127 Mean : 0.6933 Mean :-0.6326
3rd Qu.:21.00 3rd Qu.:7.0 Elyakim : 15 NA NA 3rd Qu.: 0.6444 3rd Qu.: 0.8521 3rd Qu.:-0.5234
Max. :69.00 Max. :9.0 Eshtaol : 15 NA NA Max. : 1.6364 Max. : 0.9947 Max. :-0.1031
NA NA (Other) :135 NA NA NA NA NA

Two outliers with total abundance >50. Examine:

##             subunit              point_name            datetime monitors_name
## 1: Judean Highlands  Aderet KKL Plantings 3 2017-05-25 07:45:00          <NA>
## 2:           Carmel Elyakim KKL Plantings 1 2019-04-26 06:07:00         Other
##    richness      gma abundance
## 1:       11 3.684917        54
## 2:        6 5.614404        69

Keep all plots for now.

PHI>1, hence choose negative binomial. Fit fixed and mixed models. Choose mixed model if possible, otherwise choose a model with fixed-effects only.

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

mixed model failed to converge, attempt model selection of mixed model and see if eventually it converges.

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00412153 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00794915 (tol = 0.002, component 1)
## Single term deletions
## 
## Model:
## abundance ~ subunit * year_ct + cos_td_rad + sin_td_rad + (1 | 
##     site)
##                 npar    AIC
## <none>               1582.7
## cos_td_rad         1 1588.9
## sin_td_rad         1 1584.8
## subunit:year_ct    2 1579.3

drop subunit X year

## Single term deletions
## 
## Model:
## abundance ~ subunit + year_ct + cos_td_rad + sin_td_rad + (1 | 
##     site)
##            npar    AIC
## <none>          1579.3
## subunit       2 1592.1
## year_ct       1 1577.7
## cos_td_rad    1 1585.1
## sin_td_rad    1 1581.2

drop year

## Single term deletions
## 
## Model:
## abundance ~ subunit + cos_td_rad + sin_td_rad + (1 | site)
##            npar    AIC
## <none>          1577.7
## subunit       2 1590.6
## cos_td_rad    1 1583.8
## sin_td_rad    1 1579.3

drop sine because \(\Delta AIC<2\)

## Single term deletions
## 
## Model:
## abundance ~ subunit + cos_td_rad + (1 | site)
##            npar    AIC
## <none>          1579.2
## subunit       2 1595.5
## cos_td_rad    1 1582.8

subunit and time of year remain. The final model:

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(5.0582)  ( log )
## Formula: abundance ~ subunit + cos_td_rad + (1 | site)
##    Data: P.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##   1579.2   1599.7   -783.6   1567.2      219 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7006 -0.6340 -0.2512  0.4708  4.8965 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.001695 0.04117 
## Number of obs: 225, groups:  site, 15
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     2.66798    0.11700  22.803  < 2e-16 ***
## subunitCarmel   0.14408    0.08567   1.682   0.0926 .  
## subunitGalilee -0.41092    0.08844  -4.646 3.38e-06 ***
## cos_td_rad      0.33290    0.13893   2.396   0.0166 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sbntCr sbntGl
## subunitCrml -0.390              
## subunitGall -0.399  0.491       
## cos_td_rad  -0.852  0.021  0.047
## $site
##                  (Intercept)
## Aderet           0.002221368
## Amatzia          0.003288417
## Bat Shlomo      -0.009378040
## Eitanim         -0.016330075
## Elyakim          0.002662661
## Eshtaol         -0.013511033
## Givat Yeshayahu  0.024408095
## Kabri            0.007108884
## Kerem Maharal    0.012060266
## Manara          -0.010333935
## Meron           -0.018844743
## Ofer             0.006270886
## Ramat Hashofet  -0.011544826
## Ramot Naftali    0.015570814
## Zuriel           0.006594637
## 
## with conditional variances for "site"
## $site

Interpretation of abundance model:

Observations 225
Dependent variable abundance
Type Mixed effects generalized linear model
Family Negative Binomial(5.0582)
Link log
AIC 1579.156
BIC 1599.653
Pseudo-R² (fixed effects) 0.524
Pseudo-R² (total) 0.537
Fixed Effects
exp(Est.) S.E. z val. p
(Intercept) 14.411 0.117 22.803 0.000
subunitCarmel 1.155 0.086 1.682 0.093
subunitGalilee 0.663 0.088 -4.646 0.000
cos_td_rad 1.395 0.139 2.396 0.017
Random Effects
Group Parameter Std. Dev.
site (Intercept) 0.041
Grouping Variables
Group # groups ICC
site 15 0.006 
## Confidence intervals for merMod models is an experimental feature. The
## intervals reflect only the variance of the fixed effects, not the random
## effects.
contrast estimate SE df z.ratio p.value
Judean Highlands - Carmel -0.1440769 0.0856729 Inf -1.681710 0.0926251
Judean Highlands - Galilee 0.4109161 0.0884443 Inf 4.646046 0.0000051
Carmel - Galilee 0.5549930 0.0878459 Inf 6.317801 0.0000000

significantly lower abundance in Galilee subunit compared to other two subunits (z-test post-hoc with fdr correction for 3 tests). No significant temporal trend.

Carmel and Judean highland plots have on average 8.929437 and 6.1163948 individuals more than Galilee plots, which is 74.1928811 and 6.1163948 percent higher, respectively.

community analysis using package MVabund

## Overlapping points were shifted along the y-axis to make them visible.
## 
##  PIPING TO 2nd MVFACTOR
## Only the variables Curruca.melanocephala, Streptopelia.decaocto, Turdus.merula, Parus.major, Garrulus.glandarius, Troglodytes.troglodytes, Chloris.chloris, Corvus.cornix, Alectoris.chukar, Carduelis.carduelis, Cinnyris.osea, Streptopelia.turtur 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 Curruca.melanocephala, Streptopelia.decaocto, Turdus.merula, Parus.major, Garrulus.glandarius, Troglodytes.troglodytes, Chloris.chloris, Corvus.cornix, Alectoris.chukar, Carduelis.carduelis, Cinnyris.osea, Streptopelia.turtur were included in the plot 
## (the variables with highest total abundance).

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). start model specification:

##       nb       po 
## 544.1508 621.8585
## [1] "POISSON"

## [1] "NEGATIVE BINOMIAL"

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

##       nb       po      nb1 
## 544.1508 621.8585 528.0333

The addition of the explanatory variable ‘site’ is improving the AIC of the model. stepwise selection of model:

## Single term deletions
## 
## Model:
## spp_no_rare ~ subunit * year_ct + cos_td_rad + sin_td_rad + site
##                  Df    AIC
## <none>              7392.5
## cos_td_rad       14 7424.5
## sin_td_rad       14 7409.6
## site            196 7618.1
## subunit:year_ct  28 7428.0

final model includes year, subunit, subunit X year, sampling time of year.

## 
## Coefficients: (2 not defined because of singularities)
##                        wald value Pr(>wald)   
## (Intercept)                 4.673     0.029 * 
## subunitCarmel               0.128     0.341   
## subunitGalilee              8.317     0.203   
## year_ct                     9.078     0.042 * 
## cos_td_rad                  4.890     0.048 * 
## sin_td_rad                  4.764     0.072 . 
## siteAmatzia               216.021     0.002 **
## siteBat Shlomo              0.433     0.228   
## siteEitanim                 3.908     0.154   
## siteElyakim                 0.361     0.252   
## siteEshtaol                 5.035     0.166   
## siteGivat Yeshayahu         3.083     0.594   
## siteKabri                 183.169     0.005 **
## siteKerem Maharal           0.788     0.158   
## siteManara                  3.559     0.560   
## siteMeron                  13.655     0.122   
## siteOfer                    0.320     0.269   
## siteRamat Hashofet            Inf     0.267   
## siteRamot Naftali         155.849     0.029 * 
## siteZuriel                 58.601     0.360   
## subunitCarmel:year_ct     868.986     0.082 . 
## subunitGalilee:year_ct    430.799     0.467   
## --- 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Test statistic:    NaN, 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 ~ subunit * year_ct + cos_td_rad + sin_td_rad + site
## 
## Multivariate test:
##                 Res.Df Df.diff   Dev Pr(>Dev)   
## (Intercept)        224                          
## subunit            222       2 225.4     0.01 **
## year_ct            221       1  27.3     0.03 * 
## cos_td_rad         220       1  25.8     0.04 * 
## sin_td_rad         219       1  26.7     0.04 * 
## site               205      14 609.4     0.01 **
## subunit:year_ct    205       2  91.5     0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Univariate Tests:
##                 Alectoris.chukar          Carduelis.carduelis         
##                              Dev Pr(>Dev)                 Dev Pr(>Dev)
## (Intercept)                                                           
## subunit                    8.156     0.13               1.022     0.68
## year_ct                    0.531     0.99               0.369     0.99
## cos_td_rad                  0.33     0.99               0.889     0.94
## sin_td_rad                  1.35     0.89               0.853     0.91
## site                       22.31     0.16              32.744     0.05
## subunit:year_ct            1.025     0.98               0.821     0.98
##                 Chloris.chloris          Cinnyris.osea          Corvus.cornix
##                             Dev Pr(>Dev)           Dev Pr(>Dev)           Dev
## (Intercept)                                                                  
## subunit                   7.569     0.15         2.521     0.50        20.134
## year_ct                   0.002     1.00         0.729     0.98         0.381
## cos_td_rad                1.148     0.94         0.491     0.98         0.388
## sin_td_rad                0.242     0.98         3.904     0.44         0.214
## site                      32.49     0.05        43.967     0.01        68.936
## subunit:year_ct          26.268     0.01         0.329     0.98         0.346
##                          Curruca.melanocephala          Garrulus.glandarius
##                 Pr(>Dev)                   Dev Pr(>Dev)                 Dev
## (Intercept)                                                                
## subunit             0.01                55.321     0.01              21.225
## year_ct             0.99                 0.253     0.99               2.582
## cos_td_rad          0.99                 0.141     0.99               0.015
## sin_td_rad          0.98                 1.283     0.89               1.559
## site                0.01               112.609     0.01               36.36
## subunit:year_ct     0.98                 6.213     0.54               6.028
##                          Parus.major          Prinia.gracilis         
##                 Pr(>Dev)         Dev Pr(>Dev)             Dev Pr(>Dev)
## (Intercept)                                                           
## subunit             0.01       7.029     0.21          13.377     0.02
## year_ct             0.73        0.79     0.98           2.431     0.75
## cos_td_rad          0.99       0.388     0.99           6.876     0.12
## sin_td_rad          0.89       0.071     0.98           0.063     0.98
## site                0.02      25.601     0.09          43.647     0.01
## subunit:year_ct     0.54       3.497     0.83          10.608     0.12
##                 Pycnonotus.xanthopygos          Streptopelia.decaocto         
##                                    Dev Pr(>Dev)                   Dev Pr(>Dev)
## (Intercept)                                                                   
## subunit                         28.636     0.01                22.355     0.01
## year_ct                          4.077     0.39                 5.457     0.26
## cos_td_rad                       6.935     0.11                 1.756     0.93
## sin_td_rad                       2.053     0.80                 2.518     0.71
## site                            28.335     0.06                44.621     0.01
## subunit:year_ct                  4.756     0.67                 2.621     0.95
##                 Streptopelia.turtur          Troglodytes.troglodytes         
##                                 Dev Pr(>Dev)                     Dev Pr(>Dev)
## (Intercept)                                                                  
## subunit                       6.122     0.28                   8.512     0.10
## year_ct                       7.971     0.06                   1.738     0.88
## cos_td_rad                    1.918     0.89                   0.423     0.99
## sin_td_rad                    0.339     0.98                    1.53     0.89
## site                         37.084     0.02                  66.751     0.01
## subunit:year_ct              19.455     0.02                   1.953     0.97
##                 Turdus.merula         
##                           Dev Pr(>Dev)
## (Intercept)                           
## subunit                23.453     0.01
## year_ct                 0.014     1.00
## cos_td_rad              4.108     0.46
## sin_td_rad             10.699     0.03
## site                   13.943     0.48
## subunit:year_ct         7.626     0.38
## Arguments:
##  Test statistics calculated assuming uncorrelated response (for faster computation) 
## P-value calculated using 99 iterations via PIT-trap resampling.

subunit factor is not significant, only one interaction with year is mildly significant. re-run model without subunit:year interaction and without site for simplification of the interpretation of results

## Single term deletions
## 
## Model:
## spp_no_rare ~ subunit + year_ct + cos_td_rad + sin_td_rad
##            Df    AIC
## <none>        7645.4
## subunit    28 7793.5
## year_ct    14 7645.1
## cos_td_rad 14 7653.7
## sin_td_rad 14 7644.1

drop sine

## Single term deletions
## 
## Model:
## spp_no_rare ~ subunit + year_ct + cos_td_rad
##            Df    AIC
## <none>        7644.1
## subunit    28 7803.4
## year_ct    14 7645.5
## cos_td_rad 14 7641.9

drop cosine

## Single term deletions
## 
## Model:
## spp_no_rare ~ subunit + year_ct
##         Df    AIC
## <none>     7641.9
## subunit 28 7807.1
## year_ct 14 7641.2

drop year

## Single term deletions
## 
## Model:
## spp_no_rare ~ subunit
##         Df    AIC
## <none>     7641.2
## subunit 28 7810.6

final model includes subunit only.

## 
## Test statistics:
##                wald value Pr(>wald)    
## (Intercept)        25.922     0.001 ***
## subunitCarmel       9.582     0.001 ***
## subunitGalilee     11.418     0.001 ***
## --- 
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Test statistic:  15.21, 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 ~ subunit
## 
## Multivariate test:
##             Res.Df Df.diff   Dev Pr(>Dev)   
## (Intercept)    224                          
## subunit        222       2 225.4     0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Univariate Tests:
##             Alectoris.chukar          Carduelis.carduelis         
##                          Dev Pr(>Dev)                 Dev Pr(>Dev)
## (Intercept)                                                       
## subunit                8.156     0.13               1.022     0.69
##             Chloris.chloris          Cinnyris.osea          Corvus.cornix
##                         Dev Pr(>Dev)           Dev Pr(>Dev)           Dev
## (Intercept)                                                              
## subunit               7.569     0.13         2.521     0.48        20.134
##                      Curruca.melanocephala          Garrulus.glandarius
##             Pr(>Dev)                   Dev Pr(>Dev)                 Dev
## (Intercept)                                                            
## subunit         0.01                55.321     0.01              21.225
##                      Parus.major          Prinia.gracilis         
##             Pr(>Dev)         Dev Pr(>Dev)             Dev Pr(>Dev)
## (Intercept)                                                       
## subunit         0.01       7.029     0.13          13.377     0.01
##             Pycnonotus.xanthopygos          Streptopelia.decaocto         
##                                Dev Pr(>Dev)                   Dev Pr(>Dev)
## (Intercept)                                                               
## subunit                     28.636     0.01                22.355     0.01
##             Streptopelia.turtur          Troglodytes.troglodytes         
##                             Dev Pr(>Dev)                     Dev Pr(>Dev)
## (Intercept)                                                              
## subunit                   6.122     0.20                   8.512     0.13
##             Turdus.merula         
##                       Dev Pr(>Dev)
## (Intercept)                       
## subunit            23.453     0.01
## Arguments:
##  Test statistics calculated assuming uncorrelated response (for faster computation) 
## P-value calculated using 99 iterations via PIT-trap resampling.

Factor subunit 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.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [[1]]

## 
## [[2]]