Abundance

sum(pit$abun)

## [1] 18542

sum(abundata$abun)

## [1] 16640

m1 <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = poisson(link = "log"), data = abundata)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m1)

## # Overdispersion test
## 
##        dispersion ratio =    50.874
##   Pearson's Chi-Squared = 31338.579
##                 p-value =   < 0.001

## Overdispersion detected.

m1 <- glmmTMB(abun ~ bnll +   (1|Site), family = "nbinom2", data = abundata)
summary(m1)

##  Family: nbinom2  ( log )
## Formula:          abun ~ bnll + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##   5208.8   5226.6  -2600.4   5200.8      621 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3341   0.578   
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 0.964 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.07465    0.29592  10.390   <2e-16 ***
## bnllpresent  0.09674    0.39687   0.244    0.807    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m1)

## [1] 5208.824

m1.nb <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = abundata)
summary(m1.nb)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##   5176.2   5220.6  -2578.1   5156.2      615 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05496  0.2344  
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.01 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      2.8624240  0.6026925   4.749 2.04e-06 ***
## bnllpresent     -0.4961149  0.2122171  -2.338   0.0194 *  
## arid             0.6222207  0.1271684   4.893 9.94e-07 ***
## Micrositeshrub   0.4388926  0.1049844   4.181 2.91e-05 ***
## mean.height     -0.0328309  0.0148368  -2.213   0.0269 *  
## mean.cover       0.0043418  0.0035863   1.211   0.2260    
## dry.veg.percent  0.0017713  0.0013437   1.318   0.1874    
## NDVI            -0.0002470  0.0002086  -1.184   0.2364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1.nb, type = 2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: abun
##                   Chisq Df Pr(>Chisq)    
## bnll             5.4652  1    0.01940 *  
## arid            23.9403  1  9.937e-07 ***
## Microsite       17.4770  1  2.908e-05 ***
## mean.height      4.8965  1    0.02691 *  
## mean.cover       1.4657  1    0.22603    
## dry.veg.percent  1.7377  1    0.18743    
## NDVI             1.4019  1    0.23641    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#which variable to keep

msmall <- glmmTMB(abun ~ bnll+mean.height +Microsite + arid + NDVI  +(1|Site), family = "nbinom2", data = abundata)
summary(msmall)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + mean.height + Microsite + arid + NDVI + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##   5176.3   5211.8  -2580.1   5160.3      617 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.04125  0.2031  
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 0.998 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     3.3930242  0.5292287   6.411 1.44e-10 ***
## bnllpresent    -0.5359589  0.1909696  -2.807  0.00501 ** 
## mean.height    -0.0347232  0.0139836  -2.483  0.01302 *  
## Micrositeshrub  0.4917439  0.0998475   4.925 8.44e-07 ***
## arid            0.6273386  0.1155885   5.427 5.72e-08 ***
## NDVI           -0.0003035  0.0002089  -1.453  0.14626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(msmall, type = 2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: abun
##               Chisq Df Pr(>Chisq)    
## bnll         7.8765  1   0.005008 ** 
## mean.height  6.1660  1   0.013023 *  
## Microsite   24.2551  1  8.438e-07 ***
## arid        29.4561  1  5.720e-08 ***
## NDVI         2.1108  1   0.146263    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.nbboth <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ site.shrub + mean.cover + dry.veg.percent + NDVI + (1|Site), family = "nbinom2", data = abundata)
summary(m1.nbboth)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + site.shrub + mean.cover +  
##     dry.veg.percent + NDVI + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##   5178.1   5226.9  -2578.0   5156.1      614 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05274  0.2296  
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.01 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.8860626  0.6063467   4.760 1.94e-06 ***
## bnllpresent          -0.5376023  0.2432119  -2.210   0.0271 *  
## arid                  0.6420374  0.1388907   4.623 3.79e-06 ***
## Micrositeshrub        0.4303985  0.1080659   3.983 6.81e-05 ***
## mean.height          -0.0350902  0.0163375  -2.148   0.0317 *  
## site.shrubunshrubbed -0.0836598  0.2508122  -0.334   0.7387    
## mean.cover            0.0041661  0.0036150   1.152   0.2491    
## dry.veg.percent       0.0018003  0.0013467   1.337   0.1813    
## NDVI                 -0.0002304  0.0002145  -1.074   0.2828    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.nbs <- glmmTMB(abun ~ bnll + arid + Microsite+ site.shrub + mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = abundata)
summary(m1.nbs)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + site.shrub + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##   5180.3   5224.7  -2580.2   5160.3      615 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.09327  0.3054  
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.01 
## 
## Conditional model:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.5090914  0.6035754   4.157 3.22e-05 ***
## bnllpresent          -0.4664656  0.3046587  -1.531 0.125743    
## arid                  0.6388620  0.1754252   3.642 0.000271 ***
## Micrositeshrub        0.4456415  0.1080009   4.126 3.69e-05 ***
## site.shrubunshrubbed  0.1353728  0.2870828   0.472 0.637251    
## mean.cover            0.0062807  0.0036157   1.737 0.082375 .  
## dry.veg.percent       0.0015115  0.0013341   1.133 0.257231    
## NDVI                 -0.0003895  0.0002145  -1.816 0.069376 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.nbbarid <- glmmTMB(abun ~ bnll  + Microsite+ mean.height + mean.cover + dry.veg.percent + NDVI + (1|Site), family = "nbinom2", data = abundata)
summary(m1.nbbarid)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: abundata
## 
##      AIC      BIC   logLik deviance df.resid 
##     5186     5226    -2584     5168      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2691   0.5187  
## Number of obs: 625, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.01 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.2864781  0.5971811   5.503 3.73e-08 ***
## bnllpresent      0.0557141  0.3586369   0.155 0.876546    
## Micrositeshrub   0.4114203  0.1073036   3.834 0.000126 ***
## mean.height     -0.0272202  0.0176414  -1.543 0.122838    
## mean.cover       0.0063686  0.0037181   1.713 0.086742 .  
## dry.veg.percent  0.0017586  0.0013463   1.306 0.191458    
## NDVI            -0.0001827  0.0002253  -0.811 0.417462    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m1.nb, m1.nbs, m1.nbboth, m1.nbbarid)

##            df      AIC
## m1.nb      10 5176.173
## m1.nbs     10 5180.331
## m1.nbboth  11 5178.064
## m1.nbbarid  9 5186.043

t.test(abun ~ bnll, data = abundata)

## 
##  Welch Two Sample t-test
## 
## data:  abun by bnll
## t = -1.4783, df = 622.92, p-value = 0.1398
## alternative hypothesis: true difference in means between group absent and group present is not equal to 0
## 95 percent confidence interval:
##  -13.028985   1.837676
## sample estimates:
##  mean in group absent mean in group present 
##              23.47253              29.06818

#mean height only

check_collinearity(m1.nb)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##             Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##             bnll 1.43 [1.30, 1.59]         1.19      0.70     [0.63, 0.77]
##             arid 1.46 [1.33, 1.63]         1.21      0.69     [0.61, 0.75]
##        Microsite 1.13 [1.06, 1.28]         1.06      0.88     [0.78, 0.94]
##      mean.height 1.21 [1.13, 1.36]         1.10      0.83     [0.74, 0.89]
##       mean.cover 1.12 [1.05, 1.26]         1.06      0.89     [0.79, 0.95]
##  dry.veg.percent 1.18 [1.10, 1.32]         1.09      0.85     [0.76, 0.91]
##             NDVI 1.12 [1.05, 1.26]         1.06      0.90     [0.79, 0.95]

s <- as.data.frame(round(summary(m1.nb)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/ground_all_abun.csv")
car::Anova(m1.nb, type = 2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: abun
##                   Chisq Df Pr(>Chisq)    
## bnll             5.4652  1    0.01940 *  
## arid            23.9403  1  9.937e-07 ***
## Microsite       17.4770  1  2.908e-05 ***
## mean.height      4.8965  1    0.02691 *  
## mean.cover       1.4657  1    0.22603    
## dry.veg.percent  1.7377  1    0.18743    
## NDVI             1.4019  1    0.23641    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplot(abundata, aes(bnll, abun)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Pitfall Arthropod Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank())

mean(abundata$abun)

## [1] 26.624

sd(abundata$abun)

## [1] 48.41677

s <- abundata %>% group_by(Site) %>% summarise(s = sum(abun))
mean(s$s)

## [1] 1848.889

sd(s$s)

## [1] 1136.405

Richness

m1 <- glmmTMB(Species ~  bnll  +  (1|month) +  (1|Site), family = poisson(link = "log"), data = pit)
summary(m1)

##  Family: poisson  ( log )
## Formula:          Species ~ bnll + (1 | month) + (1 | Site)
## Data: pit
## 
##      AIC      BIC   logLik deviance df.resid 
##   2893.0   2910.8  -1442.5   2885.0      622 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  month  (Intercept) 0.02304  0.1518  
##  Site   (Intercept) 0.06169  0.2484  
## Number of obs: 626, groups:  month, 3; Site, 9
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.64097    0.15434  10.632   <2e-16 ***
## bnllpresent -0.07088    0.17054  -0.416    0.678    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mh <- glmmTMB(H ~  bnll  +  (1|month) +  (1|Site), data = pit)
summary(mh)

##  Family: gaussian  ( identity )
## Formula:          H ~ bnll + (1 | month) + (1 | Site)
## Data: pit
## 
##      AIC      BIC   logLik deviance df.resid 
##   1033.9   1056.1   -511.9   1023.9      621 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  month    (Intercept) 0.009732 0.09865 
##  Site     (Intercept) 0.018773 0.13701 
##  Residual             0.290919 0.53937 
## Number of obs: 626, groups:  month, 3; Site, 9
## 
## Dispersion estimate for gaussian family (sigma^2): 0.291 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.14800    0.09494  12.092   <2e-16 ***
## bnllpresent -0.18174    0.10173  -1.787    0.074 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shapiro.test(resid(mh))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mh)
## W = 0.99298, p-value = 0.004956

m2 <- glmmTMB(Species ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = poisson(link = "log"), data = pit)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

#m2.nb <- glmmTMB(Species ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = pit)

AIC(m2)

## [1] 2842.149

#<2 difference AIC
#results the same too
check_overdispersion(m2)

## # Overdispersion test
## 
##        dispersion ratio =   1.081
##   Pearson's Chi-Squared = 667.279
##                 p-value =   0.079

## No overdispersion detected.

check_collinearity(m2)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##             Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##             bnll 1.42 [1.30, 1.58]         1.19      0.71     [0.63, 0.77]
##             arid 1.41 [1.29, 1.57]         1.19      0.71     [0.63, 0.78]
##        Microsite 1.13 [1.06, 1.27]         1.06      0.89     [0.79, 0.94]
##      mean.height 1.24 [1.15, 1.39]         1.11      0.81     [0.72, 0.87]
##       mean.cover 1.14 [1.07, 1.28]         1.07      0.88     [0.78, 0.94]
##  dry.veg.percent 1.18 [1.10, 1.33]         1.09      0.84     [0.75, 0.91]
##             NDVI 1.16 [1.09, 1.30]         1.08      0.86     [0.77, 0.92]

diagnose(m2)

## Unusually large Z-statistics (|x|>5):
## 
##    mean.height theta_1|Site.1 
##      -8.375056      -6.375236 
## 
## Large Z-statistics (estimate/std err) suggest a *possible* failure of
## the Wald approximation - often also associated with parameters that are
## at or near the edge of their range (e.g. random-effects standard
## deviations approaching 0).  (Alternately, they may simply represent
## very well-estimated parameters; intercepts of non-centered models may
## fall in this category.) While the Wald p-values and standard errors
## listed in summary() may be unreliable, profile confidence intervals
## (see ?confint.glmmTMB) and likelihood ratio test p-values derived by
## comparing models (e.g. ?drop1) are probably still OK.  (Note that the
## LRT is conservative when the null value is on the boundary, e.g. a
## variance or zero-inflation value of 0 (Self and Liang 1987; Stram and
## Lee 1994; Goldman and Whelan 2000); in simple cases these p-values are
## approximately twice as large as they should be.)

summary(m2)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll + arid + Microsite + mean.height + mean.cover +  
##     dry.veg.percent + NDVI + (1 | Site)
## Data: pit
## 
##      AIC      BIC   logLik deviance df.resid 
##   2842.1   2882.1  -1412.1   2824.1      617 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02698  0.1643  
## Number of obs: 626, groups:  Site, 9
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.1326091  0.2802254   4.042 5.30e-05 ***
## bnllpresent     -0.2211261  0.1380661  -1.602    0.109    
## arid             0.1883350  0.0803187   2.345    0.019 *  
## Micrositeshrub   0.1924537  0.0476377   4.040 5.35e-05 ***
## mean.height     -0.0564753  0.0067433  -8.375  < 2e-16 ***
## mean.cover       0.0023991  0.0015465   1.551    0.121    
## dry.veg.percent -0.0000872  0.0005830  -0.150    0.881    
## NDVI             0.0003915  0.0001000   3.915 9.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m2, type = 2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Species
##                   Chisq Df Pr(>Chisq)    
## bnll             2.5651  1    0.10924    
## arid             5.4983  1    0.01903 *  
## Microsite       16.3212  1  5.346e-05 ***
## mean.height     70.1416  1  < 2.2e-16 ***
## mean.cover       2.4067  1    0.12081    
## dry.veg.percent  0.0224  1    0.88111    
## NDVI            15.3256  1  9.048e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m2)[["coefficients"]][["cond"]] , 3))
write.csv(s,"modelresultscsv/ground_all_rich.csv")


ggplot(pit, aes(bnll, Species, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Pitfall Species Richness") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank())

mean(pit$Species)

## [1] 5.142173

sd(pit$Species)

## [1] 2.915434

Beetle abundance

#abundance
m4 <- glmmTMB(abun ~ bnll+ arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +(1|Site), family = "poisson", data = col)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m4)

## # Overdispersion test
## 
##        dispersion ratio =   1.395
##   Pearson's Chi-Squared = 860.593
##                 p-value = < 0.001

## Overdispersion detected.

m4 <- glmmTMB(abun ~ bnll +   (1|Site), family = "nbinom2", data = col)

summary(m4)

##  Family: nbinom2  ( log )
## Formula:          abun ~ bnll + (1 | Site)
## Data: col
## 
##      AIC      BIC   logLik deviance df.resid 
##   1666.0   1683.7   -829.0   1658.0      622 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2637   0.5135  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 2.66 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.1854     0.2654   0.698    0.485
## bnllpresent  -0.5661     0.3596  -1.574    0.115

m4.nb <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = col)
summary(m4.nb)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: col
## 
##      AIC      BIC   logLik deviance df.resid 
##   1656.1   1700.5   -818.1   1636.1      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.004007 0.0633  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family ():  2.8 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -6.629e-02  5.856e-01  -0.113  0.90987    
## bnllpresent     -1.134e+00  1.683e-01  -6.737 1.62e-11 ***
## arid             6.320e-01  9.937e-02   6.360 2.02e-10 ***
## Micrositeshrub   9.437e-03  1.267e-01   0.074  0.94062    
## mean.height     -3.589e-02  1.305e-02  -2.750  0.00596 ** 
## mean.cover      -1.284e-03  3.857e-03  -0.333  0.73917    
## dry.veg.percent -1.081e-03  1.504e-03  -0.718  0.47256    
## NDVI             6.952e-05  2.141e-04   0.325  0.74535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#partial models
m4.bnl <- glmmTMB(abun ~ bnll+ (1|Site), family = "nbinom2", data = col)
m4.nbb <- glmmTMB(abun ~ bnll+ arid +(1|Site), family = "nbinom2", data = col)
m4.arid <- glmmTMB(abun ~  arid +(1|Site), family = "nbinom2", data = col)
AIC(m4.bnl, m4.nbb, m4.arid)

##         df      AIC
## m4.bnl   4 1665.983
## m4.nbb   5 1655.089
## m4.arid  4 1667.652

m4.nbbnl <- glmmTMB(abun ~ bnll  + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = col)

m4.nbar <- glmmTMB(abun ~ arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI +  (1|Site), family = "nbinom2", data = col)

AIC(m4.nb, m4.nbbnl, m4.nbar)

##          df      AIC
## m4.nb    10 1656.138
## m4.nbbnl  9 1667.590
## m4.nbar   9 1669.631

#model with both a big improvement
knitr::kable(car::Anova(m4.nb, type = 2))

s <- as.data.frame(round(summary(m4.nb)[["coefficients"]][["cond"]], 3))

write.csv(s,"modelresultscsv/ground_beet_abun.csv")

ggplot(col, aes(arid, abun, color =bnll)) + geom_smooth(method = "lm")

## `geom_smooth()` using formula = 'y ~ x'

ggplot(col, aes(bnll, abun)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Teneb Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

mean(col$abun)

## [1] 1.019169

sd(col$abun)

## [1] 1.333729

s <- col %>% group_by(Site) %>% summarise(sum(abun))
mean(s$`sum(abun)`)

## [1] 70.88889

sd(s$`sum(abun)`)

## [1] 38.1623

sum(col$abun)

## [1] 638

mean(col$Species)

## [1] 0.8083067

sd(col$Species)

## [1] 0.9320917

Abundance

# Abundance

m2 <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "poisson", data = noants)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m2)

## # Overdispersion test
## 
##        dispersion ratio =    2.146
##   Pearson's Chi-Squared = 1324.090
##                 p-value =  < 0.001

## Overdispersion detected.

m2 <- glmmTMB(abun ~ bnll +  (1|Site), family = "nbinom2", data = noants)

summary(m2)

##  Family: nbinom2  ( log )
## Formula:          abun ~ bnll + (1 | Site)
## Data: noants
## 
##      AIC      BIC   logLik deviance df.resid 
##   3144.5   3162.2  -1568.2   3136.5      622 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05585  0.2363  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 3.58 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.6148     0.1255  12.872   <2e-16 ***
## bnllpresent  -0.1238     0.1684  -0.735    0.462    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2.nb <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "nbinom2", data = noants)
AIC(m2, m2.nb)

##       df      AIC
## m2     4 3144.492
## m2.nb 10 3073.857

summary(m2.nb)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: noants
## 
##      AIC      BIC   logLik deviance df.resid 
##   3073.9   3118.3  -1526.9   3053.9      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.03458  0.186   
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 4.68 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.9606500  0.3876191   2.478 0.013200 *  
## bnllpresent     -0.2332713  0.1640728  -1.422 0.155097    
## arid             0.1256104  0.0989658   1.269 0.204359    
## Micrositeshrub   0.3057563  0.0695673   4.395 1.11e-05 ***
## mean.height     -0.0660522  0.0097044  -6.806 1.00e-11 ***
## mean.cover       0.0025249  0.0021814   1.158 0.247064    
## dry.veg.percent  0.0010124  0.0008630   1.173 0.240753    
## NDVI             0.0004925  0.0001473   3.344 0.000825 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m2.nb, type = 2))

s <- summary(m2.nb)

s <- round(s[["coefficients"]][["cond"]], 3)

write.csv(s,"modelresultscsv/ground_noants_abun.csv")

ggplot(noants, aes(bnll, abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Pitfall Abundance without ants") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

Abundance

# Abundance

m2 <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "poisson", data = spids)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m2)

## # Overdispersion test
## 
##        dispersion ratio =   1.612
##   Pearson's Chi-Squared = 994.631
##                 p-value = < 0.001

## Overdispersion detected.

m2 <- glmmTMB(abun ~ bnll + (1|Site), family = "nbinom2", data = spids)
summary(m2)

##  Family: nbinom2  ( log )
## Formula:          abun ~ bnll + (1 | Site)
## Data: spids
## 
##      AIC      BIC   logLik deviance df.resid 
##   1892.1   1909.9   -942.1   1884.1      622 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.1527   0.3908  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.36 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0006394  0.2114433   0.003    0.998
## bnllpresent 0.2553805  0.2823705   0.904    0.366

m2.nb <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "nbinom2", data = spids)
AIC(m2, m2.nb)

##       df      AIC
## m2     4 1892.108
## m2.nb 10 1830.009

summary(m2.nb)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: spids
## 
##      AIC      BIC   logLik deviance df.resid 
##   1830.0   1874.4   -905.0   1810.0      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2338   0.4836  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 2.07 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.1747983  0.7962215  -2.731  0.00631 ** 
## bnllpresent      0.0857740  0.3984915   0.215  0.82957    
## arid             0.1681242  0.2312897   0.727  0.46729    
## Micrositeshrub   0.7369311  0.1255636   5.869 4.38e-09 ***
## mean.height     -0.0674879  0.0185749  -3.633  0.00028 ***
## mean.cover       0.0043735  0.0040121   1.090  0.27568    
## dry.veg.percent  0.0036545  0.0016111   2.268  0.02331 *  
## NDVI             0.0009229  0.0002892   3.191  0.00142 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m2.nb, type = 2))

s <- summary(m2.nb)

s <- round(s[["coefficients"]][["cond"]], 3)

write.csv(s,"modelresultscsv/ground_spids_abun.csv")

ggplot(spids, aes(bnll, abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Pitfall Abundance without ants") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

## Warning: Removed 8 rows containing missing values (`stat_boxplot()`).

mean(spids$abun)

## [1] 1.242902

sd(spids$abun)

## [1] 1.714446

s <- spids %>% group_by(Site) %>% summarise(s = sum(abun))
mean(s$s)

## [1] 78.8

sd(s$s)

## [1] 42.07876

sum(spids$abun)

## [1] 788

mean(spids$Species)

## [1] 0.9258675

sd(spids$Species)

## [1] 1.02382

Species richness

# Species richness
m1 <- glmmTMB(Species ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "poisson", data = nosings)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m1)

## # Overdispersion test
## 
##        dispersion ratio =   1.070
##   Pearson's Chi-Squared = 660.048
##                 p-value =   0.112

## No overdispersion detected.

summary(m1)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll + arid + Microsite + mean.height + mean.cover +  
##     dry.veg.percent + NDVI + (1 | Site)
## Data: nosings
## 
##      AIC      BIC   logLik deviance df.resid 
##   2826.9   2866.8  -1404.4   2808.9      617 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02755  0.166   
## Number of obs: 626, groups:  Site, 9
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.103e+00  2.825e-01   3.906 9.40e-05 ***
## bnllpresent     -2.229e-01  1.395e-01  -1.598   0.1100    
## arid             1.849e-01  8.120e-02   2.277   0.0228 *  
## Micrositeshrub   1.878e-01  4.800e-02   3.913 9.13e-05 ***
## mean.height     -5.684e-02  6.815e-03  -8.340  < 2e-16 ***
## mean.cover       2.420e-03  1.557e-03   1.554   0.1202    
## dry.veg.percent -4.895e-05  5.872e-04  -0.083   0.9336    
## NDVI             4.023e-04  1.009e-04   3.988 6.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(m1, type = 2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Species
##                   Chisq Df Pr(>Chisq)    
## bnll             2.5539  1    0.11002    
## arid             5.1859  1    0.02277 *  
## Microsite       15.3086  1  9.130e-05 ***
## mean.height     69.5530  1  < 2.2e-16 ***
## mean.cover       2.4152  1    0.12016    
## dry.veg.percent  0.0069  1    0.93357    
## NDVI            15.9067  1  6.654e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m1)[["coefficients"]][["cond"]],3))
write.csv(s,"modelresultscsv/ground_nosingles_rich.csv")

Orthoptera abundance

# Abundance

m8.bn <- glmmTMB(abun ~ bnll + (1|Site), family = "poisson", data = hops)
check_overdispersion(m8.bn)

## # Overdispersion test
## 
##        dispersion ratio =    1.679
##   Pearson's Chi-Squared = 1046.137
##                 p-value =  < 0.001

## Overdispersion detected.

m8.bnnb <- glmmTMB(abun ~ bnll + (1|Site), family = "nbinom2", data = hops)

summary(m8.bnnb)

##  Family: nbinom2  ( log )
## Formula:          abun ~ bnll + (1 | Site)
## Data: hops
## 
##      AIC      BIC   logLik deviance df.resid 
##    816.1    833.9   -404.1    808.1      622 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.8189   0.9049  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 0.511 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.3092     0.4764  -2.748  0.00599 **
## bnllpresent  -0.4333     0.6576  -0.659  0.51000   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m8 <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site),  family = "poisson", data = hops)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m8)

## # Overdispersion test
## 
##        dispersion ratio =   1.612
##   Pearson's Chi-Squared = 994.729
##                 p-value = < 0.001

## Overdispersion detected.

m8.nb <- glmmTMB(abun ~ bnll + arid + Microsite+ mean.height+ mean.cover + dry.veg.percent + NDVI + (1|Site), family = "nbinom2", data = hops)

AIC(m8, m8.nb)

##       df      AIC
## m8     9 877.1816
## m8.nb 10 819.9716

summary(m8.nb)

##  Family: nbinom2  ( log )
## Formula:          
## abun ~ bnll + arid + Microsite + mean.height + mean.cover + dry.veg.percent +  
##     NDVI + (1 | Site)
## Data: hops
## 
##      AIC      BIC   logLik deviance df.resid 
##    820.0    864.4   -400.0    800.0      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.7026   0.8382  
## Number of obs: 626, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 0.554 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     -1.8542529  1.4485410  -1.280   0.2005  
## bnllpresent     -0.2672522  0.7477473  -0.357   0.7208  
## arid            -0.0228044  0.4387062  -0.052   0.9585  
## Micrositeshrub  -0.6485160  0.2601943  -2.492   0.0127 *
## mean.height     -0.0277925  0.0376771  -0.738   0.4607  
## mean.cover       0.0030307  0.0081438   0.372   0.7098  
## dry.veg.percent  0.0030655  0.0031566   0.971   0.3315  
## NDVI             0.0003086  0.0005803   0.532   0.5949  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m8.nb)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/ground_hoppers_abun.csv")

ggplot(hops, aes(bnll, abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Pitfall Orthoptera Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

sum(hops$abun)

## [1] 180

mean(hops$abun)

## [1] 0.2875399

sd(hops$abun)

## [1] 0.7799926

s <- hops %>% group_by(Site) %>% summarise(s = sum(abun))
mean(s$s)

## [1] 20

sd(s$s)

## [1] 14.10674

mean(hops$Species)

## [1] 0.201278

sd(hops$Species)

## [1] 0.4282784

Canopy richness

m6 <- glmmTMB(Species ~ bnll.x+ arid +  shrub.x + mean.cover+ mean.height  + NDVI +   (1|Site),  family = "poisson", data = active.shrub)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m6)

## # Overdispersion test
## 
##        dispersion ratio =   1.089
##   Pearson's Chi-Squared = 394.078
##                 p-value =   0.118

## No overdispersion detected.

AIC(m6)

## [1] 1047.972

summary(m6)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll.x + arid + shrub.x + mean.cover + mean.height +  
##     NDVI + (1 | Site)
## Data: active.shrub
## 
##      AIC      BIC   logLik deviance df.resid 
##   1048.0   1079.3   -516.0   1032.0      362 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4269   0.6534  
## Number of obs: 370, groups:  Site, 6
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.3579412  0.8978045  -1.513  0.13040    
## bnll.xpresent -0.5359217  0.6282545  -0.853  0.39364    
## arid          -0.1057266  0.3836087  -0.276  0.78285    
## shrub.x        0.0023844  0.0005020   4.750 2.04e-06 ***
## mean.cover    -0.0049395  0.0037698  -1.310  0.19009    
## mean.height   -0.0603603  0.0183524  -3.289  0.00101 ** 
## NDVI           0.0011055  0.0002447   4.518 6.24e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m6, type = 3))

check_collinearity(m6)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##       bnll.x 1.19 [1.10, 1.39]         1.09      0.84     [0.72, 0.91]
##         arid 1.21 [1.11, 1.41]         1.10      0.83     [0.71, 0.90]
##      shrub.x 1.04 [1.00, 1.62]         1.02      0.96     [0.62, 1.00]
##   mean.cover 1.08 [1.02, 1.34]         1.04      0.93     [0.75, 0.98]
##  mean.height 1.30 [1.18, 1.51]         1.14      0.77     [0.66, 0.85]
##         NDVI 1.21 [1.11, 1.41]         1.10      0.83     [0.71, 0.90]

s <- as.data.frame(round(summary(m6)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/canopy_rich.csv")


m6bn <- glmmTMB(Species ~ bnll.x + (1|Site),  family = "poisson", data = active.shrub)

summary(m6bn)

##  Family: poisson  ( log )
## Formula:          Species ~ bnll.x + (1 | Site)
## Data: active.shrub
## 
##      AIC      BIC   logLik deviance df.resid 
##   1091.6   1103.4   -542.8   1085.6      367 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.1899   0.4358  
## Number of obs: 370, groups:  Site, 6
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)     0.4010     0.3174   1.263    0.206
## bnll.xpresent  -0.3877     0.3907  -0.992    0.321

#ggplot(active.shrub, aes(bnll, Species, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Canopy Species Richness") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

#open.wide <- read.csv("Clean Data/")
shrub.long <- read.csv("Clean Data/active.shrub.long.csv")
grp_shrub <- shrub.long %>% group_by(highest.rtu) %>% summarise(sum = sum(quantity))
sum(grp_shrub$sum)

## [1] 764

No singles

shrub_nosingles <- read.csv("Clean Data/active.shrub.cleaned.nosings.csv")
shrub_nosingles$Site <- gsub("LoK", "Lok", shrub_nosingles$Site)
shrub_nosingles <- left_join(shrub_nosingles, sites, by = c("Site", "month"))
shrub_nosingles <- left_join(shrub_nosingles, sitesrem, by = c("Site", "month"))
shrub_nosingles <- rename(shrub_nosingles, day = day.x)

shrub_nosingles <- left_join(shrub_nosingles, temp, by = c("Site", "month", "day"))

m7 <- glmmTMB(Species ~ bnll+ arid +  shrub.x + mean.cover+ mean.height + NDVI +   (1|Site),  family = "poisson", data = shrub_nosingles)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m7)

## # Overdispersion test
## 
##        dispersion ratio =   1.032
##   Pearson's Chi-Squared = 373.563
##                 p-value =   0.326

## No overdispersion detected.

AIC(m7)

## [1] 1004.822

summary(m7)

##  Family: poisson  ( log )
## Formula:          Species ~ bnll + arid + shrub.x + mean.cover + mean.height +  
##     NDVI + (1 | Site)
## Data: shrub_nosingles
## 
##      AIC      BIC   logLik deviance df.resid 
##   1004.8   1036.1   -494.4    988.8      362 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4484   0.6696  
## Number of obs: 370, groups:  Site, 6
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.4069373  0.9255937  -1.520  0.12850    
## bnllpresent -0.4926953  0.6440348  -0.765  0.44426    
## arid        -0.1987783  0.3959372  -0.502  0.61564    
## shrub.x      0.0024690  0.0005198   4.750 2.03e-06 ***
## mean.cover  -0.0044168  0.0038630  -1.143  0.25289    
## mean.height -0.0607460  0.0190838  -3.183  0.00146 ** 
## NDVI         0.0011250  0.0002541   4.428 9.51e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(shrub_nosingles$act.abun.no)

## [1] 734

m7.nb <- glmmTMB(Species ~ bnll+ arid +  shrub.x + mean.cover+   NDVI +  (1|Site),  family = "nbinom2", data = shrub_nosingles)
summary(m7.nb)

##  Family: nbinom2  ( log )
## Formula:          
## Species ~ bnll + arid + shrub.x + mean.cover + NDVI + (1 | Site)
## Data: shrub_nosingles
## 
##      AIC      BIC   logLik deviance df.resid 
##   1015.3   1046.6   -499.7    999.3      362 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3194   0.5651  
## Number of obs: 370, groups:  Site, 6
## 
## Dispersion parameter for nbinom2 family ():  235 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.8499124  0.8303871  -2.228 0.025896 *  
## bnllpresent -0.3315383  0.5467452  -0.606 0.544259    
## arid        -0.2964241  0.3454346  -0.858 0.390827    
## shrub.x      0.0027173  0.0005146   5.281 1.29e-07 ***
## mean.cover  -0.0018072  0.0037146  -0.487 0.626608    
## NDVI         0.0008123  0.0002261   3.593 0.000327 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m6, type = 3))

check_collinearity(m6)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##       bnll.x 1.19 [1.10, 1.39]         1.09      0.84     [0.72, 0.91]
##         arid 1.21 [1.11, 1.41]         1.10      0.83     [0.71, 0.90]
##      shrub.x 1.04 [1.00, 1.62]         1.02      0.96     [0.62, 1.00]
##   mean.cover 1.08 [1.02, 1.34]         1.04      0.93     [0.75, 0.98]
##  mean.height 1.30 [1.18, 1.51]         1.14      0.77     [0.66, 0.85]
##         NDVI 1.21 [1.11, 1.41]         1.10      0.83     [0.71, 0.90]

s <- as.data.frame(round(summary(m7.nb)[["coefficients"]][["cond"]],3))
write.csv(s,"modelresultscsv/shrub_nosingles_rich.csv")

Abundance

active.open$month <- factor(active.open$month, c("July", "Aug", "Sept"))

m3 <- glmmTMB(wide.sweep.abun ~ arid + bnll.x+ mean.cover+ mean.height + NDVI + (1|Site), family = "poisson", data = active.open)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

m3b <- glmmTMB(wide.sweep.abun ~ bnll.x+  (1|Site),family = "nbinom2", data = active.open)

summary(m3b)

##  Family: nbinom2  ( log )
## Formula:          wide.sweep.abun ~ bnll.x + (1 | Site)
## Data: active.open
## 
##      AIC      BIC   logLik deviance df.resid 
##   1630.6   1647.9   -811.3   1622.6      556 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.7647   0.8745  
## Number of obs: 560, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.27 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -0.03963    0.44763  -0.089    0.929
## bnll.xpresent -0.04232    0.60078  -0.070    0.944

check_overdispersion(m3)

## # Overdispersion test
## 
##        dispersion ratio =    2.317
##   Pearson's Chi-Squared = 1281.396
##                 p-value =  < 0.001

## Overdispersion detected.

m3.nb <- glmmTMB(wide.sweep.abun ~ bnll.x+ arid +  mean.height+ mean.cover + NDVI + (1|Site),family = "nbinom2", data = active.open)
AIC(m3, m3.nb)

##       df      AIC
## m3     7 1851.797
## m3.nb  8 1622.716

summary(m3.nb)

##  Family: nbinom2  ( log )
## Formula:          wide.sweep.abun ~ bnll.x + arid + mean.height + mean.cover +  
##     NDVI + (1 | Site)
## Data: active.open
## 
##      AIC      BIC   logLik deviance df.resid 
##   1622.7   1657.3   -803.4   1606.7      552 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.6514   0.8071  
## Number of obs: 560, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.33 
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -1.409e+00  1.100e+00  -1.281   0.2003  
## bnll.xpresent -6.403e-01  6.523e-01  -0.982   0.3263  
## arid           7.982e-01  3.737e-01   2.136   0.0327 *
## mean.height   -2.721e-02  2.443e-02  -1.114   0.2654  
## mean.cover     1.193e-02  4.802e-03   2.484   0.0130 *
## NDVI           3.658e-06  3.684e-04   0.010   0.9921  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m3.nb, type = 3))

check_collinearity(m3.nb)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##       bnll.x 1.37 [1.25, 1.55]         1.17      0.73     [0.65, 0.80]
##         arid 1.38 [1.26, 1.56]         1.18      0.72     [0.64, 0.79]
##  mean.height 1.19 [1.11, 1.34]         1.09      0.84     [0.74, 0.90]
##   mean.cover 1.12 [1.05, 1.28]         1.06      0.90     [0.78, 0.95]
##         NDVI 1.08 [1.03, 1.26]         1.04      0.92     [0.79, 0.97]

s <- as.data.frame(round(summary(m3.nb)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/open_abun.csv")

ggplot(active.open, aes(bnll.x, wide.sweep.abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Open Arthropod Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

s <- active.open %>% group_by(Site) %>% summarise(s = sum(wide.sweep.abun))
mean(s$s)

## [1] 85.44444

sd(s$s)

## [1] 84.50608

sum(active.open$wide.sweep.abun)

## [1] 769

mean(active.open$wide.sweep.abun)

## [1] 1.373214

sd(active.open$wide.sweep.abun)

## [1] 2.58764

mean(active.open$Species)

## [1] 0.85

sd(active.open$Species)

## [1] 1.047529

Richness

m3 <- glmmTMB(Species~ bnll.x+ arid +mean.height+ mean.cover+  NDVI  + (1|Site), family = "poisson", data = active.open)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m3)

## # Overdispersion test
## 
##        dispersion ratio =   1.054
##   Pearson's Chi-Squared = 582.608
##                 p-value =   0.185

## No overdispersion detected.

summary(m3)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll.x + arid + mean.height + mean.cover + NDVI + (1 |      Site)
## Data: active.open
## 
##      AIC      BIC   logLik deviance df.resid 
##   1313.7   1344.0   -649.9   1299.7      553 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.276    0.5253  
## Number of obs: 560, groups:  Site, 9
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -1.407e+00  8.780e-01  -1.602   0.1091  
## bnll.xpresent -2.628e-01  4.347e-01  -0.605   0.5454  
## arid           4.827e-01  2.448e-01   1.972   0.0487 *
## mean.height   -1.618e-02  1.873e-02  -0.864   0.3877  
## mean.cover     8.870e-03  3.971e-03   2.234   0.0255 *
## NDVI           3.527e-05  3.177e-04   0.111   0.9116  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m3)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/open_rich.csv")

m3b <- glmmTMB(Species~ bnll.x  + (1|Site), family = "poisson", data = active.open)
summary(m3b)

##  Family: poisson  ( log )
## Formula:          Species ~ bnll.x + (1 | Site)
## Data: active.open
## 
##      AIC      BIC   logLik deviance df.resid 
##   1317.3   1330.3   -655.6   1311.3      557 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3847   0.6202  
## Number of obs: 560, groups:  Site, 9
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -0.39595    0.32107  -1.233    0.217
## bnll.xpresent  0.08481    0.43068   0.197    0.844

#ggplot(active.open, aes(bnll, Species, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Open Arthropod Species Richness") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))


#open.wide <- read.csv("Clean Data/")
open.long <- read.csv("Clean Data/active_sweeps_long.csv")
grp_open <- open.long %>% group_by(highest.rtu) %>% summarise(sum = sum(abun))
sum(grp_open$sum)

## [1] 769

sum(active.open$wide.sweep.abun)

## [1] 769

Richness no singles

openno <- read.csv("Clean Data/active_sweeps_nosingles.csv")
openno <- rename(openno, Site = site)
openno <- left_join(openno, sites, by = c("Site", "month"))
openno <- left_join(openno, sitesrem, by = c("Site", "month"))
openno <- rename(openno, day = day.x)

openno <- left_join(openno, temp, by = c("Site", "month", "day"))


m4 <- glmmTMB(Species~ bnll+ arid +mean.height+ mean.cover+  NDVI  + (1|Site), family = "poisson", data = openno)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m4)

## # Overdispersion test
## 
##        dispersion ratio =   0.988
##   Pearson's Chi-Squared = 546.311
##                 p-value =   0.572

## No overdispersion detected.

summary(m4)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll + arid + mean.height + mean.cover + NDVI + (1 |      Site)
## Data: openno
## 
##      AIC      BIC   logLik deviance df.resid 
##   1266.1   1296.4   -626.1   1252.1      553 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2888   0.5374  
## Number of obs: 560, groups:  Site, 9
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.484e+00  9.114e-01  -1.628   0.1036  
## bnllpresent -2.539e-01  4.451e-01  -0.570   0.5684  
## arid         4.858e-01  2.515e-01   1.932   0.0534 .
## mean.height -2.097e-02  1.920e-02  -1.092   0.2749  
## mean.cover   8.945e-03  4.100e-03   2.182   0.0291 *
## NDVI         6.583e-05  3.297e-04   0.200   0.8417  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m4)[["coefficients"]][["cond"]],3))
write.csv(s,"modelresultscsv/open_nosingles_rich.csv")

Sweep net grasshopper abundance

active.hoppers$month <- factor(active.hoppers$month, c("July", "Aug", "Sept"))
# Abundance 
m13 <- glmmTMB(wide.hop.abun ~arid + bnll.x+ mean.cover+ mean.height + NDVI +(1|Site), family = "poisson", data = active.hoppers)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m13)

## # Overdispersion test
## 
##        dispersion ratio =    2.055
##   Pearson's Chi-Squared = 1136.521
##                 p-value =  < 0.001

## Overdispersion detected.

m13.nb <- glmmTMB(wide.hop.abun ~ bnll.x + arid + mean.height+ mean.cover + NDVI  + (1|Site), family = "nbinom2", data = active.hoppers)

AIC(m13, m13.nb)

##        df      AIC
## m13     7 1173.465
## m13.nb  8 1068.364

summary(m13.nb)

##  Family: nbinom2  ( log )
## Formula:          
## wide.hop.abun ~ bnll.x + arid + mean.height + mean.cover + NDVI +  
##     (1 | Site)
## Data: active.hoppers
## 
##      AIC      BIC   logLik deviance df.resid 
##   1068.4   1103.0   -526.2   1052.4      552 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.9782   0.989   
## Number of obs: 560, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family ():  1.2 
## 
## Conditional model:
##                 Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -3.236e-01  1.486e+00  -0.218   0.8276  
## bnll.xpresent -2.888e-01  8.035e-01  -0.359   0.7193  
## arid           2.451e-01  4.587e-01   0.534   0.5931  
## mean.height   -6.484e-02  3.256e-02  -1.991   0.0465 *
## mean.cover     2.478e-03  6.202e-03   0.400   0.6895  
## NDVI          -2.871e-05  5.188e-04  -0.055   0.9559  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m13.bnnl <- glmmTMB(wide.hop.abun ~ bnll.x +  (1|Site), family = "nbinom2", data = active.hoppers)
summary(m13.bnnl)

##  Family: nbinom2  ( log )
## Formula:          wide.hop.abun ~ bnll.x + (1 | Site)
## Data: active.hoppers
## 
##      AIC      BIC   logLik deviance df.resid 
##   1066.6   1083.9   -529.3   1058.6      556 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.9552   0.9773  
## Number of obs: 560, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 1.14 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -0.8711     0.5107  -1.706   0.0881 .
## bnll.xpresent  -0.1230     0.6812  -0.180   0.8567  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplot(active.hoppers, aes(bnll.x, wide.hop.abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Open Grasshopper Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

sum(active.hoppers$wide.hop.abun)

## [1] 377

Abundance

mal$month <- factor(mal$month, c("July", "Aug", "Sept"))

m3 <- glmmTMB(wide.mal.abun~ bnll+ arid +  mean.height+ mean.cover + NDVI + (1|Site), family = "poisson", data = mal)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m3)

## # Overdispersion test
## 
##        dispersion ratio =  22.626
##   Pearson's Chi-Squared = 452.522
##                 p-value = < 0.001

## Overdispersion detected.

m3.nb <- glmmTMB(wide.mal.abun ~ bnll+ arid +  mean.height+ mean.cover + NDVI + (1|Site), family = "nbinom2", data = mal)
AIC(m3, m3.nb)

##       df      AIC
## m3     7 656.1214
## m3.nb  8 297.9362

summary(m3.nb)

##  Family: nbinom2  ( log )
## Formula:          
## wide.mal.abun ~ bnll + arid + mean.height + mean.cover + NDVI +      (1 | Site)
## Data: mal
## 
##      AIC      BIC   logLik deviance df.resid 
##    297.9    308.3   -141.0    281.9       19 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.09969  0.3157  
## Number of obs: 27, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 5.17 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.6707203  1.0696336   4.367 1.26e-05 ***
## bnllpresent  0.2001823  0.3378356   0.593    0.553    
## arid         0.2166418  0.2062801   1.050    0.294    
## mean.height -0.1208385  0.0268920  -4.493 7.01e-06 ***
## mean.cover  -0.0006959  0.0069713  -0.100    0.920    
## NDVI         0.0005237  0.0004267   1.227    0.220    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m3.nb, type = 3))

check_collinearity(m3.nb)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##         bnll 1.48 [1.17,  2.35]         1.22      0.67     [0.43, 0.85]
##         arid 1.58 [1.23,  2.48]         1.26      0.63     [0.40, 0.82]
##  mean.height 1.30 [1.08,  2.15]         1.14      0.77     [0.46, 0.93]
##   mean.cover 1.05 [1.00, 10.59]         1.03      0.95     [0.09, 1.00]
##         NDVI 1.15 [1.02,  2.32]         1.07      0.87     [0.43, 0.98]

m3.bnl <- glmmTMB(wide.mal.abun ~ bnll+  (1|Site), family = "nbinom2", data = mal)

summary(m3.bnl)

##  Family: nbinom2  ( log )
## Formula:          wide.mal.abun ~ bnll + (1 | Site)
## Data: mal
## 
##      AIC      BIC   logLik deviance df.resid 
##    308.0    313.2   -150.0    300.0       23 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3944   0.628   
## Number of obs: 27, groups:  Site, 9
## 
## Dispersion parameter for nbinom2 family (): 3.35 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.3621     0.3545  12.303   <2e-16 ***
## bnllpresent   0.3790     0.4760   0.796    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

s <- as.data.frame(round(summary(m3.nb)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/mal_abun.csv")





ggplot(mal, aes(bnll, wide.mal.abun, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Malaise Abundance") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

s <- mal %>% group_by(Site) %>% summarise(s = sum(wide.mal.abun))
mean(s$s)

## [1] 380.5556

sd(s$s)

## [1] 407.9489

mean(mal$wide.mal.abun)

## [1] 126.8519

sd(mal$wide.mal.abun)

## [1] 170.8323

mean(mal$Species)

## [1] 21.55556

sd(mal$Species)

## [1] 6.302218

Richness

m4 <- glmmTMB(Species~ bnll+ arid +  mean.height+ mean.cover + NDVI + (1|Site), family = "poisson", data = mal)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m4)

## # Overdispersion test
## 
##        dispersion ratio =  0.978
##   Pearson's Chi-Squared = 19.566
##                 p-value =  0.485

## No overdispersion detected.

summary(m4)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll + arid + mean.height + mean.cover + NDVI + (1 |      Site)
## Data: mal
## 
##      AIC      BIC   logLik deviance df.resid 
##    170.4    179.5    -78.2    156.4       20 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.006041 0.07772 
## Number of obs: 27, groups:  Site, 9
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.8278599  0.4886875   7.833 4.77e-15 ***
## bnllpresent -0.1802356  0.1278134  -1.410  0.15850    
## arid         0.1258564  0.0738943   1.703  0.08853 .  
## mean.height -0.0294452  0.0113862  -2.586  0.00971 ** 
## mean.cover  -0.0022127  0.0031415  -0.704  0.48121    
## NDVI        -0.0001364  0.0001843  -0.740  0.45920    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m4, type = 3))

check_collinearity(m4)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##         bnll 1.67 [1.26,  2.69]         1.29      0.60     [0.37, 0.79]
##         arid 1.65 [1.25,  2.66]         1.28      0.61     [0.38, 0.80]
##  mean.height 1.17 [1.02,  2.39]         1.08      0.86     [0.42, 0.98]
##   mean.cover 1.11 [1.01,  3.14]         1.05      0.90     [0.32, 0.99]
##         NDVI 1.04 [1.00, 43.88]         1.02      0.96     [0.02, 1.00]

s <- as.data.frame(round(summary(m4)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/mal_rich.csv")
m4bnl <- glmmTMB(Species~ bnll+ (1|Site), family = "poisson", data = mal)
summary(m4bnl)

##  Family: poisson  ( log )
## Formula:          Species ~ bnll + (1 | Site)
## Data: mal
## 
##      AIC      BIC   logLik deviance df.resid 
##    176.7    180.6    -85.3    170.7       24 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02748  0.1658  
## Number of obs: 27, groups:  Site, 9
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.09218    0.10321  29.960   <2e-16 ***
## bnllpresent -0.06447    0.13938  -0.463    0.644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cor.test(sites$mean.height, sites$arid)

## 
##  Pearson's product-moment correlation
## 
## data:  sites$mean.height and sites$arid
## t = -1.1475, df = 25, p-value = 0.262
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5564005  0.1708552
## sample estimates:
##        cor 
## -0.2236818

ggplot(mal, aes(bnll, Species, fill = month)) + geom_boxplot() + scale_fill_brewer(palette="Paired") + ylab("Malaise Richness") + xlab("BNLL Status") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank()) + labs(fill = element_blank()) + scale_x_discrete(limits = c("present", "absent"))

### No singles

malno <- read.csv("Clean Data/mal_nosingles.csv")
malno$Site <- gsub("LoK", "Lok", malno$Site)
malno$Site <- gsub("PaPL", "PaPl", malno$Site)
malno <- left_join(malno, sites, by = c("Site", "month"))
malno <- left_join(malno, sitesrem, by = c("Site", "month"))



m5 <- glmmTMB(Species~ bnll+ arid +  mean.height+ mean.cover + NDVI + (1|Site), family = "poisson", data = malno)

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

check_overdispersion(m5)

## # Overdispersion test
## 
##        dispersion ratio =  1.035
##   Pearson's Chi-Squared = 20.697
##                 p-value =  0.415

## No overdispersion detected.

summary(m5)

##  Family: poisson  ( log )
## Formula:          
## Species ~ bnll + arid + mean.height + mean.cover + NDVI + (1 |      Site)
## Data: malno
## 
##      AIC      BIC   logLik deviance df.resid 
##    164.4    173.5    -75.2    150.4       20 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.
##  Site   (Intercept) 0.0002094 0.01447 
## Number of obs: 27, groups:  Site, 9
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.8249427  0.4852532   7.882 3.21e-15 ***
## bnllpresent -0.1877633  0.1171949  -1.602   0.1091    
## arid         0.1441664  0.0671540   2.147   0.0318 *  
## mean.height -0.0246801  0.0109567  -2.253   0.0243 *  
## mean.cover  -0.0036065  0.0031929  -1.130   0.2587    
## NDVI        -0.0001737  0.0001834  -0.947   0.3436    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

knitr::kable(car::Anova(m5, type = 3))

check_collinearity(m5)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##         Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##         bnll 1.78 [1.33,     2.86]         1.33      0.56     [0.35, 0.75]
##         arid 1.75 [1.31,     2.81]         1.32      0.57     [0.36, 0.76]
##  mean.height 1.19 [1.03,     2.32]         1.09      0.84     [0.43, 0.97]
##   mean.cover 1.15 [1.02,     2.47]         1.07      0.87     [0.40, 0.98]
##         NDVI 1.02 [1.00, 40501.06]         1.01      0.98     [0.00, 1.00]

s <- as.data.frame(round(summary(m5)[["coefficients"]][["cond"]], 3))
write.csv(s,"modelresultscsv/mal_richnosings.csv")