Libraries

install.packages("ggplot2")
install.packages("MASS")
install.packages("nlme")
install.packages("lme4")
install.packages("lmerTest")
install.packages("plyr")
install.packages("car")
install.packages("emmeans")
install.packages("openxlsx")
install.packages("data.table")
install.packages("glmmTMB")
install.packages("vegan")
install.packages("ggfortify")
install.packages("patchwork")
install.packages("tidyr")

Install libraries

## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## 
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area

Set colors

Figure 2: Percent cover in May of 2017 and 2018 in response to treatments

1) Import data

2) Run data analysis for 2017 and 2018 separately

#2017
model.4=lm(log(Mean_per_cover)~Fertilizer+Cover, data=subset(cover.dat, Year=="2017"))
summary(model.4)
## 
## Call:
## lm(formula = log(Mean_per_cover) ~ Fertilizer + Cover, data = subset(cover.dat, 
##     Year == "2017"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47485 -0.14471 -0.06535  0.12063  0.43881 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.2113     0.1182  27.180 3.79e-12 ***
## FertilizerManure   0.2582     0.1398   1.847   0.0895 .  
## CoverCover         0.3544     0.1398   2.535   0.0262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2694 on 12 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4681, Adjusted R-squared:  0.3795 
## F-statistic:  5.28 on 2 and 12 DF,  p-value: 0.02264
Anova(model.4,type="III")
## Anova Table (Type III tests)
## 
## Response: log(Mean_per_cover)
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 53.626  1 738.7301 3.789e-12 ***
## Fertilizer   0.248  1   3.4114   0.08954 .  
## Cover        0.467  1   6.4267   0.02617 *  
## Residuals    0.871 12                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(model.4$residuals)

hist(asin(cover.dat$Mean_per_cover/100))

hist(log(cover.dat$Mean_per_cover))

emmeans(model.4, specs=pairwise~Fertilizer:Cover)
## $emmeans
##  Fertilizer Cover  emmean    SE df lower.CL upper.CL
##  Inorganic  Bare     3.21 0.118 12     2.95     3.47
##  Manure     Bare     3.47 0.129 12     3.19     3.75
##  Inorganic  Cover    3.57 0.118 12     3.31     3.82
##  Manure     Cover    3.82 0.118 12     3.57     4.08
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                          estimate    SE df t.ratio p.value
##  Inorganic,Bare - Manure,Bare       -0.2582 0.140 12 -1.847  0.2999 
##  Inorganic,Bare - Inorganic,Cover   -0.3544 0.140 12 -2.535  0.1042 
##  Inorganic,Bare - Manure,Cover      -0.6126 0.191 12 -3.216  0.0326 
##  Manure,Bare - Inorganic,Cover      -0.0962 0.205 12 -0.470  0.9642 
##  Manure,Bare - Manure,Cover         -0.3544 0.140 12 -2.535  0.1042 
##  Inorganic,Cover  - Manure,Cover    -0.2582 0.140 12 -1.847  0.2999 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
#2018
model.1=lm(log(Mean_per_cover)~Fertilizer*Cover, data=subset(cover.dat, Year=="2018"))
summary(model.1)
## 
## Call:
## lm(formula = log(Mean_per_cover) ~ Fertilizer * Cover, data = subset(cover.dat, 
##     Year == "2018"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3263 -0.3308  0.1141  0.4940  0.7355 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.0333     0.3399   8.923 1.21e-06 ***
## FertilizerManure               0.3181     0.4808   0.662    0.521    
## CoverCover                    -0.1514     0.4808  -0.315    0.758    
## FertilizerManure:CoverCover    0.5488     0.6799   0.807    0.435    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6799 on 12 degrees of freedom
## Multiple R-squared:  0.2415, Adjusted R-squared:  0.05185 
## F-statistic: 1.273 on 3 and 12 DF,  p-value: 0.3277
Anova(model.1,type="III")
## Anova Table (Type III tests)
## 
## Response: log(Mean_per_cover)
##                  Sum Sq Df F value    Pr(>F)    
## (Intercept)      36.802  1 79.6151 1.209e-06 ***
## Fertilizer        0.202  1  0.4377    0.5207    
## Cover             0.046  1  0.0992    0.7583    
## Fertilizer:Cover  0.301  1  0.6516    0.4352    
## Residuals         5.547 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(model.1$residuals)

emmeans(model.1, specs=pairwise~Fertilizer:Cover)
## $emmeans
##  Fertilizer Cover  emmean   SE df lower.CL upper.CL
##  Inorganic  Bare     3.03 0.34 12     2.29     3.77
##  Manure     Bare     3.35 0.34 12     2.61     4.09
##  Inorganic  Cover    2.88 0.34 12     2.14     3.62
##  Manure     Cover    3.75 0.34 12     3.01     4.49
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                          estimate    SE df t.ratio p.value
##  Inorganic,Bare - Manure,Bare        -0.318 0.481 12 -0.662  0.9094 
##  Inorganic,Bare - Inorganic,Cover     0.151 0.481 12  0.315  0.9886 
##  Inorganic,Bare - Manure,Cover       -0.716 0.481 12 -1.488  0.4735 
##  Manure,Bare - Inorganic,Cover        0.469 0.481 12  0.977  0.7651 
##  Manure,Bare - Manure,Cover          -0.397 0.481 12 -0.827  0.8408 
##  Inorganic,Cover  - Manure,Cover     -0.867 0.481 12 -1.803  0.3184 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates
model.5=lm(log(Mean_per_cover)~Fertilizer, data=subset(cover.dat, Year=="2018"))
summary(model.5)
## 
## Call:
## lm(formula = log(Mean_per_cover) ~ Fertilizer, data = subset(cover.dat, 
##     Year == "2018"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.52507 -0.17247  0.05321  0.50715  0.63063 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.9576     0.2297  12.876 3.77e-09 ***
## FertilizerManure   0.5925     0.3248   1.824   0.0896 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6497 on 14 degrees of freedom
## Multiple R-squared:  0.192,  Adjusted R-squared:  0.1343 
## F-statistic: 3.327 on 1 and 14 DF,  p-value: 0.08956
Anova(model.5,type="III")
## Anova Table (Type III tests)
## 
## Response: log(Mean_per_cover)
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 69.977  1 165.800 3.766e-09 ***
## Fertilizer   1.404  1   3.327   0.08956 .  
## Residuals    5.909 14                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hist(model.5$residuals)

emmeans(model.5, specs=pairwise~Fertilizer)
## $emmeans
##  Fertilizer emmean   SE df lower.CL upper.CL
##  Inorganic    2.96 0.23 14     2.46     3.45
##  Manure       3.55 0.23 14     3.06     4.04
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.592 0.325 14 -1.824  0.0896 
## 
## Results are given on the log (not the response) scale.

Cover is significant in 2017 (manure is marginally significant), Manure is significant in 2018, likely because we fertilized with manure earlier in the year.

3) Plot effect of manure/cover for years

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

Save plot as “FigPerCover.eps” in figs folder

ggsave(plot,file="./figs/FigPerCover.eps",width=13,height=8,units='in')
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

Supplemental figure S1, Residue biomass on the soil surface in 2018

1) Import data

2) Run statistics

Here we are running residue weight (multiplying by 2 so that it is per meter, not per half meter)

m=lmer(Residue.weight..g.*2~Fertilizer*Cover+(1|Block), data=residue.data)
summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Residue.weight..g. * 2 ~ Fertilizer * Cover + (1 | Block)
##    Data: residue.data
## 
## REML criterion at convergence: 130.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9406 -0.6678 -0.2818  0.4553  1.4665 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 3322     57.63   
##  Residual             3497     59.14   
## Number of obs: 15, groups:  Block, 4
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                  343.305     41.287   6.424   8.315 0.000114 ***
## FertilizerManure              44.725     41.815   7.761   1.070 0.316937    
## CoverCover                   -54.173     45.954   7.953  -1.179 0.272529    
## FertilizerManure:CoverCover   34.123     62.131   7.865   0.549 0.598109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FrtlzM CvrCvr
## FertilzrMnr -0.506              
## CoverCover  -0.461  0.455       
## FrtlzrMn:CC  0.341 -0.673 -0.740
Anova(m,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Residue.weight..g. * 2
##                    Chisq Df Pr(>Chisq)    
## (Intercept)      69.1404  1     <2e-16 ***
## Fertilizer        1.1440  1     0.2848    
## Cover             1.3896  1     0.2385    
## Fertilizer:Cover  0.3016  1     0.5829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(m)

emmeans(m,specs=pairwise~Fertilizer:Cover)
## $emmeans
##  Fertilizer Cover emmean   SE   df lower.CL upper.CL
##  Inorganic  Bare     343 41.3 6.71      245      442
##  Manure     Bare     388 41.3 6.71      290      487
##  Inorganic  Cover    289 45.9 8.25      184      395
##  Manure     Cover    368 41.3 6.71      269      466
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                         estimate   SE   df t.ratio p.value
##  Inorganic,Bare - Manure,Bare        -44.7 41.8 8.00 -1.070  0.7162 
##  Inorganic,Bare - Inorganic,Cover     54.2 46.4 8.19  1.167  0.6617 
##  Inorganic,Bare - Manure,Cover       -24.7 41.8 8.00 -0.590  0.9323 
##  Manure,Bare - Inorganic,Cover        98.9 46.4 8.19  2.131  0.2208 
##  Manure,Bare - Manure,Cover           20.1 41.8 8.00  0.479  0.9615 
##  Inorganic,Cover - Manure,Cover      -78.8 46.4 8.19 -1.699  0.3824 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates
m.null=lmer(Residue.weight..g.*2~1+(1|Block), data=residue.data)
anova(m.null,m)
## refitting model(s) with ML (instead of REML)
## Data: residue.data
## Models:
## m.null: Residue.weight..g. * 2 ~ 1 + (1 | Block)
## m: Residue.weight..g. * 2 ~ Fertilizer * Cover + (1 | Block)
##        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m.null  3 178.01 180.13 -86.003   172.01                         
## m       6 178.36 182.61 -83.182   166.36 5.6425      3     0.1304

3) Plot effect of manure/cover on residue

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Save plot as “FigResidue.eps” in figs folder

ggsave(plot,file="./figs/FigResidue.eps",width=8,height=8,units='in')
## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Mesofauna analysis

1) Import data

2) Run statistics

Stats for collembola - multiply by .47 to make collembola density per liter cubed. We are separating collembola stats by month (June and August 2018)

June stats for collembola

m=glm.nb(collembola~Fertilizer*Cover, data=subset(meso.data,Month=="June"))
summary(m)
## 
## Call:
## glm.nb(formula = collembola ~ Fertilizer * Cover, data = subset(meso.data, 
##     Month == "June"), init.theta = 0.5979425775, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7653  -1.1655  -0.5290   0.1025   1.2219  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                   0.8109     0.7275   1.115    0.265
## FertilizerManure              1.2040     0.9903   1.216    0.224
## CoverCover                    0.2007     1.0189   0.197    0.844
## FertilizerManure:CoverCover  -1.0369     1.4088  -0.736    0.462
## 
## (Dispersion parameter for Negative Binomial(0.5979) family taken to be 1)
## 
##     Null deviance: 18.840  on 15  degrees of freedom
## Residual deviance: 16.899  on 12  degrees of freedom
## AIC: 85.18
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.598 
##           Std. Err.:  0.257 
## 
##  2 x log-likelihood:  -75.180
Anova(m,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: collembola
##                  LR Chisq Df Pr(>Chisq)
## Fertilizer        1.43044  1     0.2317
## Cover             0.03876  1     0.8439
## Fertilizer:Cover  0.53852  1     0.4630
emmeans(m,specs=pairwise~Fertilizer:Cover,type='response')
## $emmeans
##  Fertilizer Cover response   SE  df asymp.LCL asymp.UCL
##  Inorganic  Bare      2.25 1.64 Inf     0.541      9.36
##  Manure     Bare      7.50 5.04 Inf     2.010     27.99
##  Inorganic  Cover     2.75 1.96 Inf     0.679     11.13
##  Manure     Cover     3.25 2.29 Inf     0.818     12.91
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df z.ratio p.value
##  Inorganic,Bare / Manure,Bare     0.300 0.297 Inf -1.216  0.6168 
##  Inorganic,Bare / Inorganic,Cover 0.818 0.834 Inf -0.197  0.9973 
##  Inorganic,Bare / Manure,Cover    0.692 0.701 Inf -0.363  0.9836 
##  Manure,Bare / Inorganic,Cover    2.727 2.673 Inf  1.024  0.7355 
##  Manure,Bare / Manure,Cover       2.308 2.245 Inf  0.860  0.8257 
##  Inorganic,Cover / Manure,Cover   0.846 0.848 Inf -0.167  0.9984 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
m.null=glm.nb(collembola~1, data=subset(meso.data,Month=="June"))
anova(m.null,m,test='chisq')
## Warning in anova.negbin(m.null, m, test = "chisq"): only Chi-squared LR tests
## are implemented
## Likelihood ratio tests of Negative Binomial Models
## 
## Response: collembola
##                Model     theta Resid. df    2 x log-lik.   Test    df LR stat.
## 1                  1 0.5188986        15       -77.00139                      
## 2 Fertilizer * Cover 0.5979426        12       -75.17967 1 vs 2     3 1.821716
##     Pr(Chi)
## 1          
## 2 0.6102208

August stats for collembola

m=glm.nb(collembola~Fertilizer*Cover, data=subset(meso.data,Month=="August"))
summary(m)
## 
## Call:
## glm.nb(formula = collembola ~ Fertilizer * Cover, data = subset(meso.data, 
##     Month == "August"), init.theta = 1.848839148, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6726  -1.0775  -0.6502   0.6208   1.2446  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  2.35138    0.39879   5.896 3.72e-09 ***
## FertilizerManure            -0.18232    0.56817  -0.321    0.748    
## CoverCover                  -0.15415    0.56747  -0.272    0.786    
## FertilizerManure:CoverCover  0.02198    0.84237   0.026    0.979    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8488) family taken to be 1)
## 
##     Null deviance: 15.851  on 14  degrees of freedom
## Residual deviance: 15.581  on 11  degrees of freedom
## AIC: 105.52
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.849 
##           Std. Err.:  0.783 
## 
##  2 x log-likelihood:  -95.524
Anova(m,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: collembola
##                  LR Chisq Df Pr(>Chisq)
## Fertilizer       0.102901  1     0.7484
## Cover            0.073754  1     0.7859
## Fertilizer:Cover 0.000681  1     0.9792
emmeans(m,specs=pairwise~Fertilizer:Cover,type='response')
## $emmeans
##  Fertilizer Cover response   SE  df asymp.LCL asymp.UCL
##  Inorganic  Bare     10.50 4.19 Inf      4.81      22.9
##  Manure     Bare      8.75 3.54 Inf      3.96      19.3
##  Inorganic  Cover     9.00 3.63 Inf      4.08      19.9
##  Manure     Cover     7.67 3.63 Inf      3.03      19.4
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df z.ratio p.value
##  Inorganic,Bare / Manure,Bare     1.200 0.682 Inf  0.321  0.9886 
##  Inorganic,Bare / Inorganic,Cover 1.167 0.662 Inf  0.272  0.9930 
##  Inorganic,Bare / Manure,Cover    1.370 0.847 Inf  0.508  0.9571 
##  Manure,Bare / Inorganic,Cover    0.972 0.556 Inf -0.049  1.0000 
##  Manure,Bare / Manure,Cover       1.141 0.711 Inf  0.212  0.9966 
##  Inorganic,Cover / Manure,Cover   1.174 0.730 Inf  0.258  0.9940 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
m.null=glm.nb(collembola~1, data=subset(meso.data,Month=="August"))
anova(m.null,m,test='chisq')
## Warning in anova.negbin(m.null, m, test = "chisq"): only Chi-squared LR tests
## are implemented
## Likelihood ratio tests of Negative Binomial Models
## 
## Response: collembola
##                Model    theta Resid. df    2 x log-lik.   Test    df LR stat.
## 1                  1 1.812569        14       -95.79103                      
## 2 Fertilizer * Cover 1.848839        11       -95.52366 1 vs 2     3 0.267366
##     Pr(Chi)
## 1          
## 2 0.9660445

June stats for mites - use poisson distribution because mean = 22, sd = 18 (close enough that a poisson will fit). However the AIC value of model with poisson is higher than with negative binomial.

m=glm(total.mites~Fertilizer*Cover, family="poisson",data=subset(meso.data,Month=="June"))
summary(m)
## 
## Call:
## glm(formula = total.mites ~ Fertilizer * Cover, family = "poisson", 
##     data = subset(meso.data, Month == "June"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.7386  -2.7334  -0.8217   1.5690   7.7009  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   2.6210     0.1348  19.438  < 2e-16 ***
## FertilizerManure              0.7376     0.1639   4.499 6.82e-06 ***
## CoverCover                    0.3365     0.1765   1.906   0.0567 .  
## FertilizerManure:CoverCover  -0.3993     0.2216  -1.801   0.0716 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 204.0  on 15  degrees of freedom
## Residual deviance: 176.5  on 12  degrees of freedom
## AIC: 257.69
## 
## Number of Fisher Scoring iterations: 5
Anova(m,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: total.mites
##                  LR Chisq Df Pr(>Chisq)    
## Fertilizer        21.6396  1   3.29e-06 ***
## Cover              3.6838  1    0.05494 .  
## Fertilizer:Cover   3.2696  1    0.07057 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(m,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate   SE  df asymp.LCL asymp.UCL
##  Bare  Inorganic  13.8 1.85 Inf      10.6      17.9
##  Cover Inorganic  19.2 2.19 Inf      15.4      24.1
##  Bare  Manure     28.8 2.68 Inf      23.9      34.5
##  Cover Manure     27.0 2.60 Inf      22.4      32.6
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio     SE  df z.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.714 0.1261 Inf -1.906  0.2255 
##  Bare,Inorganic / Bare,Manure     0.478 0.0784 Inf -4.499  <.0001 
##  Bare,Inorganic / Cover,Manure    0.509 0.0844 Inf -4.074  0.0003 
##  Cover,Inorganic / Bare,Manure    0.670 0.0986 Inf -2.724  0.0327 
##  Cover,Inorganic / Cover,Manure   0.713 0.1063 Inf -2.268  0.1055 
##  Bare,Manure / Cover,Manure       1.065 0.1427 Inf  0.469  0.9659 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
m.null=glm(total.mites~1, family="poisson",data=subset(meso.data,Month=="June"))
anova(m.null,m,test='Chisq')
## Analysis of Deviance Table
## 
## Model 1: total.mites ~ 1
## Model 2: total.mites ~ Fertilizer * Cover
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        15      204.0                          
## 2        12      176.5  3   27.493 4.641e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

August stats for mites - use poisson distribution because mean = 52, sd = 12 (close to normal actually). However the AIC value of model with poisson is higher than with a normal distribution.

m=glm(total.mites~Fertilizer*Cover, family="poisson",data=subset(meso.data,Month=="August"))
summary(m)
## 
## Call:
## glm(formula = total.mites ~ Fertilizer * Cover, family = "poisson", 
##     data = subset(meso.data, Month == "August"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.7385  -0.6744  -0.0360   0.8082   3.1557  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  3.87640    0.07198  53.853   <2e-16 ***
## FertilizerManure             0.02051    0.10128   0.203    0.839    
## CoverCover                   0.12638    0.09873   1.280    0.201    
## FertilizerManure:CoverCover  0.07659    0.14265   0.537    0.591    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 47.476  on 14  degrees of freedom
## Residual deviance: 41.523  on 11  degrees of freedom
## AIC: 136.06
## 
## Number of Fisher Scoring iterations: 4
Anova(m,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: total.mites
##                  LR Chisq Df Pr(>Chisq)
## Fertilizer        0.04103  1     0.8395
## Cover             1.64187  1     0.2001
## Fertilizer:Cover  0.28810  1     0.5914
emmeans(m,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate   SE  df asymp.LCL asymp.UCL
##  Bare  Inorganic  48.2 3.47 Inf      41.9      55.6
##  Cover Inorganic  54.8 3.70 Inf      48.0      62.5
##  Bare  Manure     49.2 3.51 Inf      42.8      56.6
##  Cover Manure     60.3 4.48 Inf      52.2      69.8
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio     SE  df z.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.881 0.0870 Inf -1.280  0.5757 
##  Bare,Inorganic / Bare,Manure     0.980 0.0992 Inf -0.203  0.9971 
##  Bare,Inorganic / Cover,Manure    0.800 0.0827 Inf -2.160  0.1346 
##  Cover,Inorganic / Bare,Manure    1.112 0.1092 Inf  1.078  0.7029 
##  Cover,Inorganic / Cover,Manure   0.907 0.0912 Inf -0.967  0.7684 
##  Bare,Manure / Cover,Manure       0.816 0.0840 Inf -1.971  0.1988 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
m.null=glm(total.mites~1, family="poisson",data=subset(meso.data,Month=="August"))
anova(m.null,m,test='Chisq')
## Analysis of Deviance Table
## 
## Model 1: total.mites ~ 1
## Model 2: total.mites ~ Fertilizer * Cover
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        14     47.476                     
## 2        11     41.523  3   5.9531   0.1139

3) Plot effect of manure and cover on collembola and mites in June and august

Start by making data tidy (so we can use facet grid)

long.data=melt(meso.data,id.vars=c("Year","Plot","Month","Month.Order","sample.date","other","Cover","Fertilizer"),variable.name="meso_group",value.name="abundance")

Use these tidy data to graph Save plot as “FigPerCover.eps” in figs folder

ggsave(plot,file="./figs/FigMesofauna.eps",width=13,height=10,units='in')

Pitfall analysis

1) Import data

We think that edges have an effect on predator abundance, particularly carabids, so adding information about edges to dataset

Harvestmen activity-density

2) Harvestmen statistics

Using a zero-inflated poisson distribution, running stats for each year separately

All years togehter harvestmen decrease overall in manure ferilized plots

harvestmen.model=glmmTMB(harvestmen~Fertilizer+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(harvestmen.model)
##  Family: poisson  ( log )
## Formula:          
## harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 |      Year)
## Zero inflation:              ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##    836.2    858.5   -412.1    824.2      299 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 2.650e-01 0.5148130
##  Plot        (Intercept) 1.923e-01 0.4385199
##  Year        (Intercept) 1.624e-09 0.0000403
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        0.6172     0.3058   2.018  0.04359 * 
## FertilizerManure  -0.6290     0.2100  -2.995  0.00274 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6924     0.2101  -3.296 0.000981 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(harvestmen.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: harvestmen
##              Chisq Df Pr(>Chisq)   
## (Intercept) 4.0725  1   0.043586 * 
## Fertilizer  8.9722  1   0.002741 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
harvestmen.null=glmmTMB(harvestmen~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(harvestmen.model,harvestmen.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## harvestmen.null: harvestmen ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## harvestmen.model: harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 | , zi=~1, disp=~1
## harvestmen.model:     Year), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## harvestmen.null   5 843.08 861.68 -416.54   833.08                            
## harvestmen.model  6 836.16 858.48 -412.08   824.16 8.9197      1   0.002821 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(harvestmen.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE  df lower.CL upper.CL
##  Inorganic  1.854 0.567 299    1.015     3.38
##  Manure     0.988 0.316 299    0.527     1.85
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure  1.88 0.394 299 2.995   0.0030 
## 
## Tests are performed on the log scale

2016

harvestmen.model=glmmTMB(harvestmen~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(harvestmen.model)
##  Family: poisson  ( log )
## Formula:          harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    241.4    254.0   -115.7    231.4       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 3.796e-03 6.161e-02
##  Plot        (Intercept) 9.205e-10 3.034e-05
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.3646     0.1787   2.041   0.0413 *
## FertilizerManure  -0.6285     0.2586  -2.430   0.0151 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -1.3148     0.6126  -2.146   0.0319 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(harvestmen.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: harvestmen
##              Chisq Df Pr(>Chisq)  
## (Intercept) 4.1638  1    0.04130 *
## Fertilizer  5.9065  1    0.01509 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
harvestmen.null=glmmTMB(harvestmen~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(harvestmen.model,harvestmen.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## harvestmen.null: harvestmen ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## harvestmen.model: harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## harvestmen.null   4 245.05 255.13 -118.52   237.05                           
## harvestmen.model  5 241.36 253.97 -115.68   231.36 5.6877      1    0.01708 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(harvestmen.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE df lower.CL upper.CL
##  Inorganic  1.440 0.257 87    1.010     2.05
##  Manure     0.768 0.189 87    0.471     1.25
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure  1.87 0.485 87 2.430   0.0171 
## 
## Tests are performed on the log scale

2017

harvestmen.model=glmmTMB(harvestmen~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(harvestmen.model)
##  Family: poisson  ( log )
## Formula:          
## harvestmen ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    210.2    227.4    -98.1    196.2       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.25906  0.509   
##  Plot        (Intercept) 0.02527  0.159   
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                              Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  0.004635   0.450929   0.010    0.992
## FertilizerManure            -0.509865   0.449073  -1.135    0.256
## CoverCover                   0.296186   0.370583   0.799    0.424
## FertilizerManure:CoverCover -0.221555   0.630221  -0.352    0.725
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.7699     0.5549  -1.387    0.165
      Anova(harvestmen.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: harvestmen
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.0001  1     0.9918
## Fertilizer       1.2891  1     0.2562
## Cover            0.6388  1     0.4241
## Fertilizer:Cover 0.1236  1     0.7252
harvestmen.null=glmmTMB(harvestmen~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(harvestmen.model,harvestmen.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## harvestmen.null: harvestmen ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## harvestmen.model: harvestmen ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## harvestmen.null   4 207.85 217.66 -99.923   199.85                         
## harvestmen.model  7 210.17 227.35 -98.085   196.17 3.6751      3     0.2988
emmeans(harvestmen.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE df lower.CL upper.CL
##  Bare  Inorganic  1.005 0.453 79    0.409     2.46
##  Cover Inorganic  1.351 0.604 79    0.555     3.29
##  Bare  Manure     0.603 0.303 79    0.222     1.64
##  Cover Manure     0.650 0.336 79    0.232     1.82
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.744 0.276 79 -0.799  0.8546 
##  Bare,Inorganic / Bare,Manure     1.665 0.748 79  1.135  0.6689 
##  Bare,Inorganic / Cover,Manure    1.545 0.712 79  0.945  0.7810 
##  Cover,Inorganic / Bare,Manure    2.239 0.972 79  1.856  0.2552 
##  Cover,Inorganic / Cover,Manure   2.078 0.930 79  1.634  0.3657 
##  Bare,Manure / Cover,Manure       0.928 0.472 79 -0.147  0.9989 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
harvestmen.model=glmmTMB(harvestmen~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(harvestmen.model)
##  Family: poisson  ( log )
## Formula:          harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    206.8    219.0    -98.4    196.8       81 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.23627  0.4861  
##  Plot        (Intercept) 0.06579  0.2565  
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.1371     0.4099   0.334   0.7380  
## FertilizerManure  -0.6286     0.3350  -1.877   0.0606 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.8180     0.5855  -1.397    0.162
      Anova(harvestmen.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: harvestmen
##              Chisq Df Pr(>Chisq)  
## (Intercept) 0.1119  1    0.73799  
## Fertilizer  3.5216  1    0.06057 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
harvestmen.null=glmmTMB(harvestmen~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(harvestmen.model,harvestmen.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## harvestmen.null: harvestmen ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## harvestmen.model: harvestmen ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)  
## harvestmen.null   4 207.85 217.66 -99.923   199.85                          
## harvestmen.model  5 206.76 219.03 -98.378   196.76 3.089      1    0.07882 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(harvestmen.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate   SE df lower.CL upper.CL
##  Inorganic  1.147 0.47 81    0.507     2.59
##  Manure     0.612 0.27 81    0.254     1.47
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure  1.88 0.628 81 1.877   0.0642 
## 
## Tests are performed on the log scale

2018

harvestmen.model=glmmTMB(harvestmen~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(harvestmen.model)
##  Family: poisson  ( log )
## Formula:          
## harvestmen ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    369.7    389.6   -177.9    355.7      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.7159   0.8461  
##  Plot        (Intercept) 0.3419   0.5847  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  0.34513    0.57052   0.605    0.545
## FertilizerManure            -0.39699    0.48237  -0.823    0.411
## CoverCover                  -0.47939    0.48567  -0.987    0.324
## FertilizerManure:CoverCover -0.09644    0.70624  -0.137    0.891
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.4481     0.5469  -2.648  0.00809 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(harvestmen.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: harvestmen
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.3659  1     0.5452
## Fertilizer       0.6773  1     0.4105
## Cover            0.9743  1     0.3236
## Fertilizer:Cover 0.0186  1     0.8914
harvestmen.null=glmmTMB(harvestmen~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(harvestmen.model,harvestmen.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## harvestmen.null: harvestmen ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## harvestmen.model: harvestmen ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## harvestmen.null   4 366.88 378.26 -179.44   358.88                         
## harvestmen.model  7 369.72 389.63 -177.86   355.72 3.1623      3     0.3673
emmeans(harvestmen.model,specs=pairwise~Fertilizer:Cover,type='response')
## $emmeans
##  Fertilizer Cover  rate    SE  df lower.CL upper.CL
##  Inorganic  Bare  1.412 0.806 120    0.456     4.37
##  Manure     Bare  0.949 0.533 120    0.313     2.88
##  Inorganic  Cover 0.874 0.494 120    0.285     2.68
##  Manure     Cover 0.534 0.314 120    0.167     1.71
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Inorganic,Bare / Manure,Bare      1.49 0.717 120 0.823   0.8435 
##  Inorganic,Bare / Inorganic,Cover  1.62 0.784 120 0.987   0.7571 
##  Inorganic,Bare / Manure,Cover     2.65 1.355 120 1.899   0.2341 
##  Manure,Bare / Inorganic,Cover     1.09 0.527 120 0.170   0.9982 
##  Manure,Bare / Manure,Cover        1.78 0.911 120 1.124   0.6753 
##  Inorganic,Cover / Manure,Cover    1.64 0.844 120 0.957   0.7738 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Make figures for Harvestmen - looking to adjust 2016 to be half the size of 2017 and 2018

Save plot as “FigPerCover.eps” in figs folder

ggsave(plot,file="./figs/FigHarvestmen.eps",width=16,height=10,units='in')

Spider activity-density

2) Wolf-spdider statistics

Overall wolf-spider model

wolf.model=glmmTMB(wolf.spiders~Fertilizer+Cover+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(wolf.model)
##  Family: poisson  ( log )
## Formula:          wolf.spiders ~ Fertilizer + Cover + (1 | Month.Order) + (1 |  
##     Plot) + (1 | Year)
## Zero inflation:                ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##     1144     1170     -565     1130      298 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 1.59233  1.2619  
##  Plot        (Intercept) 0.09151  0.3025  
##  Year        (Intercept) 0.28920  0.5378  
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       0.32538    0.71728   0.454   0.6501  
## FertilizerManure -0.23048    0.11930  -1.932   0.0534 .
## CoverCover       -0.05005    0.14202  -0.352   0.7245  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.6810     0.3215  -5.228 1.71e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(wolf.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: wolf.spiders
##              Chisq Df Pr(>Chisq)  
## (Intercept) 0.2058  1    0.65009  
## Fertilizer  3.7325  1    0.05336 .
## Cover       0.1242  1    0.72452  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wolf.null=glmmTMB(wolf.spiders~1+(1|Month.Order)+(1|Plot)+(1|Year),,data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(wolf.model,wolf.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## wolf.null: wolf.spiders ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## wolf.model: wolf.spiders ~ Fertilizer + Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## wolf.model:     Plot) + (1 | Year), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wolf.null   5 1144.1 1162.7 -567.06   1134.1                         
## wolf.model  7 1144.0 1170.0 -564.98   1130.0 4.1648      2     0.1246
emmeans(wolf.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  1.38 0.993 298    0.337     5.68
##  Cover Inorganic  1.32 0.951 298    0.318     5.46
##  Bare  Manure     1.10 0.789 298    0.268     4.51
##  Cover Manure     1.05 0.756 298    0.252     4.34
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic  1.05 0.149 298 0.352   0.9850 
##  Bare,Inorganic / Bare,Manure      1.26 0.150 298 1.932   0.2169 
##  Bare,Inorganic / Cover,Manure     1.32 0.249 298 1.494   0.4423 
##  Cover,Inorganic / Bare,Manure     1.20 0.219 298 0.985   0.7581 
##  Cover,Inorganic / Cover,Manure    1.26 0.150 298 1.932   0.2169 
##  Bare,Manure / Cover,Manure        1.05 0.149 298 0.352   0.9850 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Using a zero-inflated poisson distribution, running stats for each year separately

2016

wolf.model=glmmTMB(wolf.spiders~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(wolf.model)
##  Family: poisson  ( log )
## Formula:          wolf.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    440.5    453.1   -215.2    430.5       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.90266  0.9501  
##  Plot        (Intercept) 0.03424  0.1850  
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.9689     0.4973   1.948   0.0514 .
## FertilizerManure  -0.1535     0.1776  -0.864   0.3874  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.4601     0.5742  -4.284 1.83e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(wolf.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: wolf.spiders
##              Chisq Df Pr(>Chisq)  
## (Intercept) 3.7964  1    0.05136 .
## Fertilizer  0.7471  1    0.38740  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wolf.null=glmmTMB(wolf.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(wolf.model,wolf.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## wolf.null: wolf.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## wolf.model: wolf.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wolf.null   4 439.16 449.25 -215.58   431.16                         
## wolf.model  5 440.46 453.07 -215.23   430.46 0.6988      1     0.4032
emmeans(wolf.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer rate   SE df lower.CL upper.CL
##  Inorganic  2.63 1.31 87    0.981     7.08
##  Manure     2.26 1.12 87    0.842     6.06
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure  1.17 0.207 87 0.864   0.3898 
## 
## Tests are performed on the log scale

2017

wolf.model=glmmTMB(wolf.spiders~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(wolf.model)
##  Family: poisson  ( log )
## Formula:          wolf.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 |  
##     Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    187.1    204.3    -86.6    173.1       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 4.180e+00 2.045e+00
##  Plot        (Intercept) 8.318e-10 2.884e-05
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -0.8507     1.2617  -0.674  0.50015    
## FertilizerManure             -0.8455     0.2565  -3.296  0.00098 ***
## CoverCover                   -0.6096     0.3248  -1.877  0.06056 .  
## FertilizerManure:CoverCover   0.7336     0.4366   1.680  0.09294 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.3457     0.4681  -2.875  0.00404 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(wolf.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: wolf.spiders
##                    Chisq Df Pr(>Chisq)    
## (Intercept)       0.4546  1  0.5001531    
## Fertilizer       10.8655  1  0.0009797 ***
## Cover             3.5220  1  0.0605610 .  
## Fertilizer:Cover  2.8227  1  0.0929406 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wolf.null=glmmTMB(wolf.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(wolf.model,wolf.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## wolf.null: wolf.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## wolf.model: wolf.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## wolf.model:     Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## wolf.null   4 190.50 200.31 -91.247   182.50                           
## wolf.model  7 187.11 204.29 -86.557   173.11 9.3806      3    0.02464 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(wolf.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE df lower.CL upper.CL
##  Bare  Inorganic  0.427 0.539 79   0.0347     5.26
##  Cover Inorganic  0.232 0.297 79   0.0182     2.96
##  Bare  Manure     0.183 0.233 79   0.0147     2.29
##  Cover Manure     0.208 0.264 79   0.0165     2.60
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 1.840 0.598 79  1.877  0.2464 
##  Bare,Inorganic / Bare,Manure     2.329 0.597 79  3.296  0.0079 
##  Bare,Inorganic / Cover,Manure    2.058 0.536 79  2.768  0.0347 
##  Cover,Inorganic / Bare,Manure    1.266 0.444 79  0.672  0.9073 
##  Cover,Inorganic / Cover,Manure   1.118 0.395 79  0.317  0.9889 
##  Bare,Manure / Cover,Manure       0.883 0.258 79 -0.424  0.9742 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2018

wolf.model=glmmTMB(wolf.spiders~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(wolf.model)
##  Family: poisson  ( log )
## Formula:          wolf.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 |  
##     Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    427.8    447.8   -206.9    413.8      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 2.7178   1.6486  
##  Plot        (Intercept) 0.1094   0.3307  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                 -0.23606    0.86044  -0.274    0.784
## FertilizerManure             0.21971    0.27541   0.798    0.425
## CoverCover                   0.09842    0.27848   0.353    0.724
## FertilizerManure:CoverCover -0.34947    0.39201  -0.892    0.373
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -20.77    7161.69  -0.003    0.998
      Anova(wolf.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: wolf.spiders
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.0753  1     0.7838
## Fertilizer       0.6364  1     0.4250
## Cover            0.1249  1     0.7238
## Fertilizer:Cover 0.7948  1     0.3727
wolf.null=glmmTMB(wolf.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(wolf.model,wolf.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## wolf.null: wolf.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## wolf.model: wolf.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## wolf.model:     Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wolf.null   4 422.83 434.21 -207.42   414.83                         
## wolf.model  7 427.85 447.76 -206.93   413.85 0.9849      3     0.8049
emmeans(wolf.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  0.790 0.680 120    0.144     4.34
##  Cover Inorganic  0.871 0.750 120    0.159     4.79
##  Bare  Manure     0.984 0.845 120    0.180     5.39
##  Cover Manure     0.765 0.659 120    0.139     4.21
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.906 0.252 120 -0.353  0.9848 
##  Bare,Inorganic / Bare,Manure     0.803 0.221 120 -0.798  0.8553 
##  Bare,Inorganic / Cover,Manure    1.032 0.290 120  0.112  0.9995 
##  Cover,Inorganic / Bare,Manure    0.886 0.243 120 -0.442  0.9710 
##  Cover,Inorganic / Cover,Manure   1.139 0.318 120  0.465  0.9665 
##  Bare,Manure / Cover,Manure       1.285 0.356 120  0.907  0.8012 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2) Other spiders statistics

Using a zero-inflated poisson distribution, running stats for each year separately

Overall model

other.spid.model=glmmTMB(other.spiders~Fertilizer+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(other.spid.model)
##  Family: poisson  ( log )
## Formula:          
## other.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot) +      (1 | Year)
## Zero inflation:                 ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##    767.4    789.8   -377.7    755.4      299 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 1.82737  1.3518  
##  Plot        (Intercept) 0.04301  0.2074  
##  Year        (Intercept) 0.22105  0.4702  
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      -0.13953    0.74724  -0.187    0.852
## FertilizerManure -0.08735    0.15904  -0.549    0.583
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.4763     0.2050  -2.323   0.0202 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(other.spid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: other.spiders
##              Chisq Df Pr(>Chisq)
## (Intercept) 0.0349  1     0.8519
## Fertilizer  0.3016  1     0.5829
other.null=glmmTMB(other.spiders~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(other.spid.model,other.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## other.null: other.spiders ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## other.spid.model: other.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + , zi=~1, disp=~1
## other.spid.model:     (1 | Year), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## other.null        5 765.75 784.35 -377.87   755.75                         
## other.spid.model  6 767.44 789.76 -377.72   755.44 0.3082      1     0.5788
emmeans(other.spid.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE  df lower.CL upper.CL
##  Inorganic  0.870 0.650 299    0.200     3.78
##  Manure     0.797 0.598 299    0.182     3.49
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure  1.09 0.174 299 0.549   0.5833 
## 
## Tests are performed on the log scale

2016

other.spid.model=glmmTMB(other.spiders~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(other.spid.model)
##  Family: poisson  ( log )
## Formula:          other.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:                 ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    170.2    182.8    -80.1    160.2       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 8.003e-01 8.946e-01
##  Plot        (Intercept) 4.902e-09 7.001e-05
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       -0.9978     0.5545   -1.80   0.0719 .
## FertilizerManure   0.5914     0.3305    1.79   0.0735 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.8766     0.5926  -1.479    0.139
      Anova(other.spid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: other.spiders
##              Chisq Df Pr(>Chisq)  
## (Intercept) 3.2383  1    0.07194 .
## Fertilizer  3.2023  1    0.07353 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
other.null=glmmTMB(other.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(other.spid.model,other.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## other.null: other.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## other.spid.model: other.spiders ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## other.null        4 170.86 180.94 -81.428   162.86                         
## other.spid.model  5 170.21 182.81 -80.103   160.21 2.6508      1     0.1035
emmeans(other.spid.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE df lower.CL upper.CL
##  Inorganic  0.369 0.204 87    0.122     1.11
##  Manure     0.666 0.357 87    0.230     1.93
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure 0.554 0.183 87 -1.789  0.0770 
## 
## Tests are performed on the log scale

2017

other.spid.model=glmmTMB(other.spiders~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(other.spid.model)
##  Family: poisson  ( log )
## Formula:          
## other.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 |      Plot)
## Zero inflation:                 ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    138.0    155.2    -62.0    124.0       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 6.1659   2.4831  
##  Plot        (Intercept) 0.3979   0.6308  
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                  -1.1600     1.7402  -0.667  0.50504   
## FertilizerManure             -2.9403     1.0322  -2.849  0.00439 **
## CoverCover                    0.1208     0.6756   0.179  0.85812   
## FertilizerManure:CoverCover  -0.6665     1.5782  -0.422  0.67279   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.001046   0.469427  -0.002    0.998
      Anova(other.spid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: other.spiders
##                   Chisq Df Pr(>Chisq)   
## (Intercept)      0.4443  1    0.50504   
## Fertilizer       8.1148  1    0.00439 **
## Cover            0.0320  1    0.85812   
## Fertilizer:Cover 0.1784  1    0.67279   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
other.null=glmmTMB(other.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(other.spid.model,other.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## other.null: other.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## other.spid.model: other.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## other.spid.model:     Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## other.null        4 144.14 153.95 -68.068   136.14                            
## other.spid.model  7 138.01 155.19 -62.004   124.01 12.129      3   0.006955 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(other.spid.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer   rate     SE df lower.CL upper.CL
##  Bare  Inorganic  0.3135 0.5455 79 0.009816   10.012
##  Cover Inorganic  0.3537 0.6177 79 0.010945   11.433
##  Bare  Manure     0.0166 0.0321 79 0.000351    0.783
##  Cover Manure     0.0096 0.0199 79 0.000155    0.595
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                          ratio     SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic  0.886  0.599 79 -0.179  0.9980 
##  Bare,Inorganic / Bare,Manure     18.922 19.531 79  2.849  0.0281 
##  Bare,Inorganic / Cover,Manure    32.656 41.184 79  2.764  0.0351 
##  Cover,Inorganic / Bare,Manure    21.351 20.947 79  3.120  0.0132 
##  Cover,Inorganic / Cover,Manure   36.848 45.078 79  2.948  0.0214 
##  Bare,Manure / Cover,Manure        1.726  2.456 79  0.383  0.9807 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2018

other.spid.model=glmmTMB(other.spiders~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(other.spid.model)
##  Family: poisson  ( log )
## Formula:          
## other.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 |      Plot)
## Zero inflation:                 ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    369.4    389.3   -177.7    355.4      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 1.7988   1.3412  
##  Plot        (Intercept) 0.1425   0.3775  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  -0.4005     0.7346  -0.545    0.586
## FertilizerManure              0.3208     0.3576   0.897    0.370
## CoverCover                    0.1904     0.3580   0.532    0.595
## FertilizerManure:CoverCover  -0.3394     0.5015  -0.677    0.498
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.0629     0.5135  -4.018 5.88e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(other.spid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: other.spiders
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.2972  1     0.5856
## Fertilizer       0.8049  1     0.3696
## Cover            0.2829  1     0.5948
## Fertilizer:Cover 0.4582  1     0.4985
other.null=glmmTMB(other.spiders~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(other.spid.model,other.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## other.null: other.spiders ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## other.spid.model: other.spiders ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## other.spid.model:     Plot), zi=~1, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## other.null        4 364.14 375.51 -178.07   356.14                         
## other.spid.model  7 369.35 389.26 -177.68   355.35 0.7858      3     0.8529
emmeans(other.spid.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  0.670 0.492 120    0.156     2.87
##  Cover Inorganic  0.811 0.590 120    0.192     3.43
##  Bare  Manure     0.923 0.672 120    0.219     3.90
##  Cover Manure     0.796 0.584 120    0.186     3.40
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.827 0.296 120 -0.532  0.9512 
##  Bare,Inorganic / Bare,Manure     0.726 0.259 120 -0.897  0.8063 
##  Bare,Inorganic / Cover,Manure    0.842 0.308 120 -0.469  0.9657 
##  Cover,Inorganic / Bare,Manure    0.878 0.302 120 -0.379  0.9814 
##  Cover,Inorganic / Cover,Manure   1.019 0.359 120  0.053  0.9999 
##  Bare,Manure / Cover,Manure       1.161 0.409 120  0.423  0.9745 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Make figures for all spiders

looking to adjust 2016 to be half the size of 2017 and 2018

Start by making data tidy (so we can use facet grid)

long.pit.data=melt(pitfall.data,id.vars=c("Year","Month","Month.Order","start.date","end.date","Plot","PLOT.ID","row","Cover","Fertilizer","notes"),variable.name="pitfall_group",value.name="activityDensity")
## Warning in melt.data.table(pitfall.data, id.vars = c("Year", "Month",
## "Month.Order", : 'measure.vars' [carabid.beetles, carabid.juveniles,
## rove.beetles, other.beetles, ...] are not all of the same type. By order of
## hierarchy, the molten data value column will be of type 'character'. All measure
## variables not of type 'character' will be coerced too. Check DETAILS in ?
## melt.data.table for more on coercion.
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_summary).

Save plot as “FigPerCover.eps” in figs folder

ggsave(plot,file="./figs/FigSpiders.eps",width=16,height=12,units='in')
## Warning: Removed 1 rows containing non-finite values (stat_ydensity).
## Warning: Removed 1 rows containing non-finite values (stat_summary).

Ant activity-density

2) Ant statistics

Overall model

ant.model=glmmTMB(ants~Fertilizer*Cover+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(ant.model)
##  Family: poisson  ( log )
## Formula:          ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot) +  
##     (1 | Year)
## Zero inflation:        ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##   1242.5   1272.3   -613.3   1226.5      297 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 6.453e-02 2.540e-01
##  Plot        (Intercept) 9.328e-02 3.054e-01
##  Year        (Intercept) 5.262e-09 7.254e-05
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  0.92802    0.18011   5.152 2.57e-07 ***
## FertilizerManure            -0.01448    0.14400  -0.101    0.920    
## CoverCover                  -0.22507    0.20257  -1.111    0.267    
## FertilizerManure:CoverCover  0.21129    0.30132   0.701    0.483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.2242     0.1811  -6.761 1.37e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(ant.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: ants
##                    Chisq Df Pr(>Chisq)    
## (Intercept)      26.5470  1  2.572e-07 ***
## Fertilizer        0.0101  1     0.9199    
## Cover             1.2344  1     0.2665    
## Fertilizer:Cover  0.4917  1     0.4832    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ant.null=glmmTMB(ants~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(ant.model,ant.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## ant.null: ants ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## ant.model: ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot) + , zi=~1, disp=~1
## ant.model:     (1 | Year), zi=~1, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ant.null   5 1237.8 1256.4 -613.89   1227.8                         
## ant.model  8 1242.5 1272.3 -613.25   1226.5 1.2759      3     0.7349
emmeans(ant.model,specs=pairwise~Fertilizer,type='response')
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer rate    SE  df lower.CL upper.CL
##  Inorganic  2.26 0.389 297     1.61     3.17
##  Manure     2.48 0.441 297     1.74     3.52
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure 0.913 0.133 297 -0.626  0.5319 
## 
## Results are averaged over the levels of: Cover 
## Tests are performed on the log scale

Using a zero-inflated poisson distribution, running stats for each year separately

2016

ant.model=glmmTMB(ants~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(ant.model)
##  Family: poisson  ( log )
## Formula:          ants ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:        ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    377.7    390.4   -183.9    367.7       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.01308  0.1144  
##  Plot        (Intercept) 0.03670  0.1916  
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       0.91973    0.16908    5.44 5.34e-08 ***
## FertilizerManure  0.06262    0.21564    0.29    0.772    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9569     0.2805  -3.412 0.000645 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(ant.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: ants
##               Chisq Df Pr(>Chisq)    
## (Intercept) 29.5908  1  5.336e-08 ***
## Fertilizer   0.0843  1     0.7715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ant.null=glmmTMB(ants~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(ant.model,ant.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## ant.null: ants ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## ant.model: ants ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ant.null   4 375.83 385.92 -183.91   367.83                         
## ant.model  5 377.75 390.35 -183.87   367.75 0.0844      1     0.7714
emmeans(ant.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer rate    SE df lower.CL upper.CL
##  Inorganic  2.51 0.424 87     1.79     3.51
##  Manure     2.67 0.445 87     1.92     3.72
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure 0.939 0.203 87 -0.290  0.7722 
## 
## Tests are performed on the log scale

2017

ant.model=glmmTMB(ants~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(ant.model)
##  Family: poisson  ( log )
## Formula:          ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:        ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    349.0    366.2   -167.5    335.0       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.4134   0.6430  
##  Plot        (Intercept) 0.3629   0.6024  
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  0.29803    0.53213   0.560    0.575
## FertilizerManure            -0.03225    0.51056  -0.063    0.950
## CoverCover                  -0.13553    0.51425  -0.264    0.792
## FertilizerManure:CoverCover  0.78837    0.69555   1.133    0.257
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.5760     0.4273  -3.688 0.000226 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(ant.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: ants
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.3137  1     0.5754
## Fertilizer       0.0040  1     0.9496
## Cover            0.0695  1     0.7921
## Fertilizer:Cover 1.2847  1     0.2570
ant.null=glmmTMB(ants~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(ant.model,ant.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## ant.null: ants ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## ant.model: ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##           Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## ant.null   4 345.16 354.97 -168.58   337.16                        
## ant.model  7 348.97 366.15 -167.49   334.97 2.185      3     0.5349
emmeans(ant.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate    SE df lower.CL upper.CL
##  Bare  Inorganic  1.35 0.717 79    0.467     3.89
##  Cover Inorganic  1.18 0.618 79    0.414     3.35
##  Bare  Manure     1.30 0.683 79    0.460     3.70
##  Cover Manure     2.51 1.200 79    0.966     6.50
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 1.145 0.589 79  0.264  0.9935 
##  Bare,Inorganic / Bare,Manure     1.033 0.527 79  0.063  0.9999 
##  Bare,Inorganic / Cover,Manure    0.538 0.268 79 -1.245  0.6006 
##  Cover,Inorganic / Bare,Manure    0.902 0.460 79 -0.202  0.9970 
##  Cover,Inorganic / Cover,Manure   0.469 0.227 79 -1.567  0.4031 
##  Bare,Manure / Cover,Manure       0.521 0.252 79 -1.346  0.5365 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2018

ant.model=glmmTMB(ants~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(ant.model)
##  Family: poisson  ( log )
## Formula:          ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:        ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    471.3    491.2   -228.7    457.3      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.039916 0.1998  
##  Plot        (Intercept) 0.004449 0.0667  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   0.9367     0.1632   5.741 9.42e-09 ***
## FertilizerManure              0.1295     0.1723   0.751   0.4524    
## CoverCover                   -0.1371     0.1923  -0.713   0.4759    
## FertilizerManure:CoverCover  -0.7016     0.2954  -2.375   0.0175 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.9335     0.4396  -4.398 1.09e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(ant.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: ants
##                    Chisq Df Pr(>Chisq)    
## (Intercept)      32.9576  1  9.419e-09 ***
## Fertilizer        0.5646  1    0.45241    
## Cover             0.5083  1    0.47587    
## Fertilizer:Cover  5.6419  1    0.01754 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ant.null=glmmTMB(ants~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(ant.model,ant.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## ant.null: ants ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## ant.model: ants ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## ant.null   4 477.14 488.51 -234.57   469.14                            
## ant.model  7 471.34 491.25 -228.67   457.34 11.798      3   0.008107 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(ant.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  2.55 0.416 120    1.847     3.52
##  Cover Inorganic  2.22 0.407 120    1.549     3.19
##  Bare  Manure     2.90 0.462 120    2.119     3.98
##  Cover Manure     1.26 0.281 120    0.806     1.96
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 1.147 0.221 120  0.713  0.8918 
##  Bare,Inorganic / Bare,Manure     0.879 0.151 120 -0.751  0.8759 
##  Bare,Inorganic / Cover,Manure    2.032 0.466 120  3.095  0.0129 
##  Cover,Inorganic / Bare,Manure    0.766 0.145 120 -1.411  0.4950 
##  Cover,Inorganic / Cover,Manure   1.772 0.423 120  2.396  0.0833 
##  Bare,Manure / Cover,Manure       2.313 0.526 120  3.686  0.0019 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Manure + Cover significantly lower ant count in 2018

Rove beetle activity-density

2) Rove statistics

Overall model: using zero-inflated poisson

rove.model=glmmTMB(rove.beetles~Fertilizer+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(rove.model)
##  Family: poisson  ( log )
## Formula:          rove.beetles ~ Fertilizer + (1 | Month.Order) + (1 | Plot) +  
##     (1 | Year)
## Zero inflation:                ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##    734.8    757.1   -361.4    722.8      299 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.3264   0.5713  
##  Plot        (Intercept) 0.2070   0.4550  
##  Year        (Intercept) 0.1909   0.4369  
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       -0.4105     0.4433  -0.926   0.3545  
## FertilizerManure   0.3986     0.2225   1.792   0.0732 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.3683     0.2371  -1.553     0.12
      Anova(rove.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: rove.beetles
##              Chisq Df Pr(>Chisq)  
## (Intercept) 0.8574  1    0.35446  
## Fertilizer  3.2105  1    0.07317 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rove.null=glmmTMB(rove.beetles~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(rove.model,rove.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## rove.null: rove.beetles ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## rove.model: rove.beetles ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + , zi=~1, disp=~1
## rove.model:     (1 | Year), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## rove.null   5 735.88 754.48 -362.94   725.88                           
## rove.model  6 734.80 757.13 -361.40   722.80 3.0788      1    0.07932 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(rove.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE  df lower.CL upper.CL
##  Inorganic  0.663 0.294 299    0.277     1.59
##  Manure     0.988 0.428 299    0.422     2.32
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure 0.671 0.149 299 -1.792  0.0742 
## 
## Tests are performed on the log scale

Using a zero-inflated poisson distribution, running stats for each year separately

2016

rove.model=glmmTMB(rove.beetles~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(rove.model)
##  Family: poisson  ( log )
## Formula:          rove.beetles ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    188.0    200.6    -89.0    178.0       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 7.155e-01 0.8458746
##  Plot        (Intercept) 3.211e-08 0.0001792
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      -0.38874    0.52159  -0.745    0.456
## FertilizerManure  0.06459    0.31632   0.204    0.838
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.5466     0.4512  -1.212    0.226
      Anova(rove.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: rove.beetles
##              Chisq Df Pr(>Chisq)
## (Intercept) 0.5555  1     0.4561
## Fertilizer  0.0417  1     0.8382
rove.null=glmmTMB(rove.beetles~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(rove.model,rove.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## rove.null: rove.beetles ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## rove.model: rove.beetles ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## rove.null   4 186.00 196.08 -88.998   178.00                         
## rove.model  5 187.95 200.56 -88.977   177.95 0.0417      1     0.8382
emmeans(rove.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE df lower.CL upper.CL
##  Inorganic  0.678 0.354 87    0.240     1.91
##  Manure     0.723 0.364 87    0.266     1.97
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure 0.937 0.297 87 -0.204  0.8387 
## 
## Tests are performed on the log scale

2017

rove.model=glmmTMB(rove.beetles~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(rove.model)
##  Family: poisson  ( log )
## Formula:          rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 |  
##     Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    114.2    131.3    -50.1    100.2       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 3.3417   1.8280  
##  Plot        (Intercept) 0.3986   0.6313  
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  -1.9946     1.4679  -1.359    0.174
## FertilizerManure              0.6365     0.8211   0.775    0.438
## CoverCover                   -1.8644     1.2769  -1.460    0.144
## FertilizerManure:CoverCover   2.0317     1.4616   1.390    0.165
## 
## Zero-inflation model:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.007282   0.548181  -0.013    0.989
      Anova(rove.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: rove.beetles
##                   Chisq Df Pr(>Chisq)
## (Intercept)      1.8464  1     0.1742
## Fertilizer       0.6009  1     0.4382
## Cover            2.1319  1     0.1443
## Fertilizer:Cover 1.9321  1     0.1645
rove.null=glmmTMB(rove.beetles~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(rove.model,rove.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## rove.null: rove.beetles ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## rove.model: rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## rove.model:     Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## rove.null   4 114.42 124.24 -53.211   106.42                           
## rove.model  7 114.16 131.34 -50.082   100.16 6.2577      3    0.09972 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(rove.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer   rate     SE df lower.CL upper.CL
##  Bare  Inorganic  0.1361 0.1997 79 0.007326    2.527
##  Cover Inorganic  0.0211 0.0367 79 0.000659    0.675
##  Bare  Manure     0.2572 0.3589 79 0.015984    4.137
##  Cover Manure     0.3040 0.4149 79 0.020086    4.600
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                          ratio     SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 6.4519 8.2383 79  1.460  0.4663 
##  Bare,Inorganic / Bare,Manure     0.5291 0.4345 79 -0.775  0.8654 
##  Bare,Inorganic / Cover,Manure    0.4476 0.3557 79 -1.012  0.7432 
##  Cover,Inorganic / Bare,Manure    0.0820 0.1010 79 -2.030  0.1859 
##  Cover,Inorganic / Cover,Manure   0.0694 0.0846 79 -2.188  0.1354 
##  Bare,Manure / Cover,Manure       0.8460 0.5949 79 -0.238  0.9952 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2018

rove.model=glmmTMB(rove.beetles~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(rove.model)
##  Family: poisson  ( log )
## Formula:          rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 |  
##     Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    403.5    423.4   -194.7    389.5      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.3721   0.6100  
##  Plot        (Intercept) 0.2295   0.4791  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                 -0.15678    0.48174  -0.325   0.7448  
## FertilizerManure             0.76114    0.45495   1.673   0.0943 .
## CoverCover                   0.09751    0.48224   0.202   0.8398  
## FertilizerManure:CoverCover -0.70941    0.65201  -1.088   0.2766  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.9470     0.4133  -2.291   0.0219 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(rove.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: rove.beetles
##                   Chisq Df Pr(>Chisq)  
## (Intercept)      0.1059  1    0.74485  
## Fertilizer       2.7990  1    0.09432 .
## Cover            0.0409  1    0.83976  
## Fertilizer:Cover 1.1839  1    0.27657  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rove.null=glmmTMB(rove.beetles~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(rove.model,rove.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## rove.null: rove.beetles ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## rove.model: rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## rove.model:     Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## rove.null   4 400.44 411.82 -196.22   392.44                        
## rove.model  7 403.48 423.39 -194.74   389.48 2.959      3      0.398
emmeans(rove.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  0.855 0.412 120    0.329     2.22
##  Cover Inorganic  0.942 0.452 120    0.364     2.44
##  Bare  Manure     1.830 0.816 120    0.757     4.42
##  Cover Manure     0.992 0.459 120    0.398     2.48
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.907 0.437 120 -0.202  0.9971 
##  Bare,Inorganic / Bare,Manure     0.467 0.213 120 -1.673  0.3424 
##  Bare,Inorganic / Cover,Manure    0.861 0.408 120 -0.315  0.9891 
##  Cover,Inorganic / Bare,Manure    0.515 0.236 120 -1.451  0.4706 
##  Cover,Inorganic / Cover,Manure   0.950 0.451 120 -0.109  0.9995 
##  Bare,Manure / Cover,Manure       1.844 0.813 120  1.389  0.5090 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Centipedes activity-density

2) centipede statistics

Overall model

centi.model=glmmTMB(centipedes~Fertilizer+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
      summary(centi.model)
##  Family: poisson  ( log )
## Formula:          
## centipedes ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 |      Year)
## Zero inflation:              ~1
## Data: subset(pitfall.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##    342.7    365.0   -165.3    330.7      299 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 3.150e-01 5.613e-01
##  Plot        (Intercept) 2.961e-09 5.441e-05
##  Year        (Intercept) 4.724e-01 6.873e-01
## Number of obs: 305, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -1.6208     0.5899  -2.748    0.006 **
## FertilizerManure   0.7228     0.2844   2.541    0.011 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.1535     0.4310  -0.356    0.722
      Anova(centi.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: centipedes
##              Chisq Df Pr(>Chisq)   
## (Intercept) 7.5503  1    0.00600 **
## Fertilizer  6.4584  1    0.01104 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
centi.null=glmmTMB(centipedes~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(pitfall.data),ziformula=~1,family=poisson)
anova(centi.model,centi.null,test="Chisq")
## Data: subset(pitfall.data)
## Models:
## centi.null: centipedes ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## centi.model: centipedes ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 | , zi=~1, disp=~1
## centi.model:     Year), zi=~1, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## centi.null   5 346.83 365.43 -168.41   336.83                           
## centi.model  6 342.69 365.01 -165.34   330.69 6.1428      1     0.0132 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(centi.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE  df lower.CL upper.CL
##  Inorganic  0.198 0.117 299   0.0619    0.631
##  Manure     0.407 0.230 299   0.1338    1.240
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure 0.485 0.138 299 -2.541  0.0115 
## 
## Tests are performed on the log scale

Using a zero-inflated poisson distribution, running stats for each year separately

2016

centi.model=glmmTMB(centipedes~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(centi.model)
##  Family: poisson  ( log )
## Formula:          centipedes ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    167.5    180.1    -78.8    157.5       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.46250  0.6801  
##  Plot        (Intercept) 0.02929  0.1711  
## Number of obs: 92, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       -0.9525     0.5575  -1.708   0.0876 .
## FertilizerManure   0.9768     0.4176   2.339   0.0193 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.1093     0.4515  -0.242    0.809
      Anova(centi.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: centipedes
##              Chisq Df Pr(>Chisq)  
## (Intercept) 2.9187  1    0.08756 .
## Fertilizer  5.4706  1    0.01934 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
centi.null=glmmTMB(centipedes~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2016"),ziformula=~1,family=poisson)
anova(centi.model,centi.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2016")
## Models:
## centi.null: centipedes ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## centi.model: centipedes ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## centi.null   4 169.93 180.01 -80.962   161.93                           
## centi.model  5 167.54 180.15 -78.769   157.54 4.3867      1    0.03622 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(centi.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer  rate    SE df lower.CL upper.CL
##  Inorganic  0.386 0.215 87    0.127     1.17
##  Manure     1.025 0.498 87    0.390     2.69
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure 0.377 0.157 87 -2.339  0.0216 
## 
## Tests are performed on the log scale

2017

centi.model=glmmTMB(centipedes~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(centi.model)
##  Family: poisson  ( log )
## Formula:          
## centipedes ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:              ~1
## Data: subset(pitfall.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##     70.7     87.9    -28.4     56.7       79 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance  Std.Dev. 
##  Month.Order (Intercept) 2.500e-01 5.000e-01
##  Plot        (Intercept) 9.074e-12 3.012e-06
## Number of obs: 86, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                  -2.3952     0.7857  -3.049   0.0023 **
## FertilizerManure             -0.2134     1.0009  -0.213   0.8312   
## CoverCover                    0.4797     0.8671   0.553   0.5801   
## FertilizerManure:CoverCover  -0.8329     1.5045  -0.554   0.5799   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -19.66   19736.55  -0.001    0.999
      Anova(centi.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: centipedes
##                   Chisq Df Pr(>Chisq)   
## (Intercept)      9.2932  1     0.0023 **
## Fertilizer       0.0455  1     0.8312   
## Cover            0.3061  1     0.5801   
## Fertilizer:Cover 0.3065  1     0.5799   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
centi.null=glmmTMB(centipedes~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2017"),ziformula=~1,family=poisson)
anova(centi.model,centi.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2017")
## Models:
## centi.null: centipedes ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## centi.model: centipedes ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## centi.null   4 65.976 75.793 -28.988   57.976                         
## centi.model  7 70.710 87.891 -28.355   56.710 1.2654      3     0.7374
emmeans(centi.model,specs=pairwise~Cover*Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer   rate     SE df lower.CL upper.CL
##  Bare  Inorganic  0.0912 0.0716 79  0.01908    0.435
##  Cover Inorganic  0.1473 0.0910 79  0.04306    0.504
##  Bare  Manure     0.0736 0.0585 79  0.01514    0.358
##  Cover Manure     0.0517 0.0545 79  0.00636    0.421
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.619 0.537 79 -0.553  0.9454 
##  Bare,Inorganic / Bare,Manure     1.238 1.239 79  0.213  0.9965 
##  Bare,Inorganic / Cover,Manure    1.762 2.159 79  0.462  0.9670 
##  Cover,Inorganic / Bare,Manure    2.000 1.732 79  0.800  0.8540 
##  Cover,Inorganic / Cover,Manure   2.847 3.190 79  0.934  0.7868 
##  Bare,Manure / Cover,Manure       1.423 1.746 79  0.288  0.9916 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

2018

rove.model=glmmTMB(rove.beetles~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(rove.model)
##  Family: poisson  ( log )
## Formula:          rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 |  
##     Plot)
## Zero inflation:                ~1
## Data: subset(pitfall.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    403.5    423.4   -194.7    389.5      120 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.3721   0.6100  
##  Plot        (Intercept) 0.2295   0.4791  
## Number of obs: 127, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                 -0.15678    0.48174  -0.325   0.7448  
## FertilizerManure             0.76114    0.45495   1.673   0.0943 .
## CoverCover                   0.09751    0.48224   0.202   0.8398  
## FertilizerManure:CoverCover -0.70941    0.65201  -1.088   0.2766  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.9470     0.4133  -2.291   0.0219 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(rove.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: rove.beetles
##                   Chisq Df Pr(>Chisq)  
## (Intercept)      0.1059  1    0.74485  
## Fertilizer       2.7990  1    0.09432 .
## Cover            0.0409  1    0.83976  
## Fertilizer:Cover 1.1839  1    0.27657  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rove.null=glmmTMB(rove.beetles~1+(1|Month.Order)+(1|Plot),data=subset(pitfall.data,Year=="2018"),ziformula=~1,family=poisson)
anova(rove.model,rove.null,test="Chisq")
## Data: subset(pitfall.data, Year == "2018")
## Models:
## rove.null: rove.beetles ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## rove.model: rove.beetles ~ Fertilizer * Cover + (1 | Month.Order) + (1 | , zi=~1, disp=~1
## rove.model:     Plot), zi=~1, disp=~1
##            Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## rove.null   4 400.44 411.82 -196.22   392.44                        
## rove.model  7 403.48 423.39 -194.74   389.48 2.959      3      0.398
emmeans(rove.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  0.855 0.412 120    0.329     2.22
##  Cover Inorganic  0.942 0.452 120    0.364     2.44
##  Bare  Manure     1.830 0.816 120    0.757     4.42
##  Cover Manure     0.992 0.459 120    0.398     2.48
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.907 0.437 120 -0.202  0.9971 
##  Bare,Inorganic / Bare,Manure     0.467 0.213 120 -1.673  0.3424 
##  Bare,Inorganic / Cover,Manure    0.861 0.408 120 -0.315  0.9891 
##  Cover,Inorganic / Bare,Manure    0.515 0.236 120 -1.451  0.4706 
##  Cover,Inorganic / Cover,Manure   0.950 0.451 120 -0.109  0.9995 
##  Bare,Manure / Cover,Manure       1.844 0.813 120  1.389  0.5090 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale

Carabid beetles

1) Import data

Characterize the beetles

Identify plots with edges

Analyze totals for each predator group and graph

Overall model with cover and fertilizer

carabid.model=glmmTMB(Total~Fertilizer+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(carabid.data),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          
## Total ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 | Year)
## Zero inflation:         ~1
## Data: subset(carabid.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##   1115.2   1137.2   -551.6   1103.2      282 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.1947   0.4412  
##  Plot        (Intercept) 0.3524   0.5936  
##  Year        (Intercept) 0.1633   0.4041  
## Number of obs: 288, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.8743     0.3706   2.359   0.0183 *
## FertilizerManure  -0.1230     0.1483  -0.829   0.4070  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6524     0.1564   -4.17 3.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##              Chisq Df Pr(>Chisq)  
## (Intercept) 5.5641  1    0.01833 *
## Fertilizer  0.6875  1    0.40701  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(carabid.data),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data)
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## carabid.model: Total ~ Fertilizer + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## carabid.null   5 1113.9 1132.2 -551.95   1103.9                         
## carabid.model  6 1115.2 1137.2 -551.60   1103.2 0.6854      1     0.4077
emmeans(carabid.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer rate    SE  df lower.CL upper.CL
##  Inorganic  2.40 0.888 282     1.16     4.97
##  Manure     2.12 0.784 282     1.02     4.39
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE  df t.ratio p.value
##  Inorganic / Manure  1.13 0.168 282 0.829   0.4077 
## 
## Tests are performed on the log scale

Edge effect overall

# Edge effects
carabid.model=glmmTMB(Total~Edge+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(carabid.data),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Edge + (1 | Month.Order) + (1 | Plot) + (1 | Year)
## Zero inflation:         ~1
## Data: subset(carabid.data)
## 
##      AIC      BIC   logLik deviance df.resid 
##   1111.7   1137.4   -548.9   1097.7      281 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.1912   0.4372  
##  Plot        (Intercept) 0.2378   0.4877  
##  Year        (Intercept) 0.1200   0.3464  
## Number of obs: 288, groups:  Month.Order, 4; Plot, 16; Year, 3
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.53233    0.42213   1.261   0.2073  
## Edge1        0.08442    0.31396   0.269   0.7880  
## Edge2        0.67344    0.36337   1.853   0.0638 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6676     0.1574  -4.241 2.23e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##              Chisq Df Pr(>Chisq)  
## (Intercept) 1.5903  1    0.20729  
## Edge        7.1412  2    0.02814 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot)+(1|Year),data=subset(carabid.data),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data)
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
## carabid.model: Total ~ Edge + (1 | Month.Order) + (1 | Plot) + (1 | Year), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## carabid.null   5 1113.9 1132.2 -551.95   1103.9                           
## carabid.model  7 1111.7 1137.3 -548.86   1097.7 6.1812      2    0.04547 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(carabid.model,specs=pairwise~Edge,type='response')
## $emmeans
##  Edge rate    SE  df lower.CL upper.CL
##  0    1.70 0.719 281    0.742     3.91
##  1    1.85 0.643 281    0.936     3.67
##  2    3.34 1.203 281    1.643     6.79
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast ratio    SE  df t.ratio p.value
##  0 / 1    0.919 0.289 281 -0.269  0.9609 
##  0 / 2    0.510 0.185 281 -1.853  0.1544 
##  1 / 2    0.555 0.125 281 -2.605  0.0261 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

2016

carabid.model=glmmTMB(Total~Fertilizer+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Fertilizer + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    305.1    315.9   -147.5    295.1       59 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.02594  0.1611  
##  Plot        (Intercept) 0.33304  0.5771  
## Number of obs: 64, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        0.9978     0.3712   2.688  0.00719 **
## FertilizerManure   0.3861     0.4723   0.817  0.41371   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.7778     0.3242  -2.399   0.0164 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##              Chisq Df Pr(>Chisq)   
## (Intercept) 7.2245  1   0.007191 **
## Fertilizer  0.6681  1   0.413712   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2016"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2016")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Fertilizer + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## carabid.null   4 303.79 312.43 -147.90   295.79                         
## carabid.model  5 305.09 315.89 -147.55   295.09 0.7005      1     0.4026
emmeans(carabid.model,specs=pairwise~Fertilizer,type='response')
## $emmeans
##  Fertilizer rate   SE df lower.CL upper.CL
##  Inorganic  2.71 1.01 59     1.29      5.7
##  Manure     3.99 1.28 59     2.10      7.6
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio    SE df t.ratio p.value
##  Inorganic / Manure  0.68 0.321 59 -0.817  0.4170 
## 
## Tests are performed on the log scale
#edge effects
carabid.model=glmmTMB(Total~Edge+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2016"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Edge + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2016")
## 
##      AIC      BIC   logLik deviance df.resid 
##    305.3    316.1   -147.7    295.3       59 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.0262   0.1619  
##  Plot        (Intercept) 0.3462   0.5884  
## Number of obs: 64, groups:  Month.Order, 4; Plot, 8
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.0344     0.3775   2.740  0.00614 **
## Edge2         0.3157     0.4826   0.654  0.51303   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.7768     0.3247  -2.392   0.0168 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##              Chisq Df Pr(>Chisq)   
## (Intercept) 7.5074  1   0.006145 **
## Edge        0.4279  1   0.513034   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2016"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2016")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Edge + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## carabid.null   4 303.79 312.43 -147.90   295.79                        
## carabid.model  5 305.34 316.14 -147.67   295.34 0.452      1     0.5014
emmeans(carabid.model,specs=pairwise~Edge,type='response')
## $emmeans
##  Edge rate   SE df lower.CL upper.CL
##  1    2.81 1.06 59     1.32     5.99
##  2    3.86 1.27 59     2.00     7.44
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast ratio    SE df t.ratio p.value
##  1 / 2    0.729 0.352 59 -0.654  0.5156 
## 
## Tests are performed on the log scale

2017

carabid.model=glmmTMB(Total~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    427.6    445.5   -206.8    413.6       89 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.3592   0.5993  
##  Plot        (Intercept) 0.4583   0.6770  
## Number of obs: 96, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                  1.22231    0.58244   2.099   0.0359 *
## FertilizerManure            -0.48026    0.61174  -0.785   0.4324  
## CoverCover                  -0.15739    0.59808  -0.263   0.7924  
## FertilizerManure:CoverCover -0.07022    0.82999  -0.085   0.9326  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.5486     0.2549  -2.152   0.0314 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##                   Chisq Df Pr(>Chisq)  
## (Intercept)      4.4042  1    0.03585 *
## Fertilizer       0.6163  1    0.43242  
## Cover            0.0692  1    0.79243  
## Fertilizer:Cover 0.0072  1    0.93257  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2017"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2017")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## carabid.null   4 422.83 433.08 -207.41   414.83                         
## carabid.model  7 427.58 445.53 -206.79   413.58 1.2464      3     0.7419
emmeans(carabid.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer rate    SE df lower.CL upper.CL
##  Bare  Inorganic  3.40 1.977 89    1.067    10.80
##  Cover Inorganic  2.90 1.492 89    1.044     8.06
##  Bare  Manure     2.10 1.112 89    0.734     6.01
##  Cover Manure     1.67 0.916 89    0.564     4.96
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic  1.17 0.700 89 0.263   0.9936 
##  Bare,Inorganic / Bare,Manure      1.62 0.989 89 0.785   0.8611 
##  Bare,Inorganic / Cover,Manure     2.03 1.272 89 1.129   0.6725 
##  Cover,Inorganic / Bare,Manure     1.38 0.753 89 0.592   0.9342 
##  Cover,Inorganic / Cover,Manure    1.73 0.975 89 0.979   0.7615 
##  Bare,Manure / Cover,Manure        1.26 0.723 89 0.395   0.9789 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
# Edge effects
carabid.model=glmmTMB(Total~Edge+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2017"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Edge + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    419.5    434.9   -203.7    407.5       90 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.3412   0.5841  
##  Plot        (Intercept) 0.2647   0.5145  
## Number of obs: 96, groups:  Month.Order, 3; Plot, 16
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.8145     0.4857   1.677   0.0936 .
## Edge1        -0.2730     0.4132  -0.661   0.5087  
## Edge2         0.9388     0.4505   2.084   0.0372 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.5285     0.2481  -2.131   0.0331 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##               Chisq Df Pr(>Chisq)   
## (Intercept)  2.8121  1   0.093554 . 
## Edge        10.8035  2   0.004509 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2017"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2017")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Edge + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## carabid.null   4 422.83 433.08 -207.41   414.83                           
## carabid.model  6 419.49 434.88 -203.75   407.49 7.3322      2    0.02558 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(carabid.model,specs=pairwise~Edge,type='response')
## $emmeans
##  Edge rate    SE df lower.CL upper.CL
##  0    2.26 1.097 90    0.860     5.93
##  1    1.72 0.711 90    0.756     3.91
##  2    5.77 2.611 90    2.352    14.18
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast ratio    SE df t.ratio p.value
##  0 / 1    1.314 0.543 90  0.661  0.7867 
##  0 / 2    0.391 0.176 90 -2.084  0.0989 
##  1 / 2    0.298 0.111 90 -3.249  0.0046 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

2018

carabid.model=glmmTMB(Total~Fertilizer*Cover+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    371.8    391.8   -178.9    357.8      121 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.44743  0.6689  
##  Plot        (Intercept) 0.07061  0.2657  
## Number of obs: 128, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##                              Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  0.266778   0.409943   0.651    0.515
## FertilizerManure            -0.393665   0.346551  -1.136    0.256
## CoverCover                   0.182825   0.313069   0.584    0.559
## FertilizerManure:CoverCover  0.004506   0.477724   0.009    0.992
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -1.122      0.371  -3.025  0.00248 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##                   Chisq Df Pr(>Chisq)
## (Intercept)      0.4235  1     0.5152
## Fertilizer       1.2904  1     0.2560
## Cover            0.3410  1     0.5592
## Fertilizer:Cover 0.0001  1     0.9925
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2018"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2018")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Fertilizer * Cover + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## carabid.null   4 368.53 379.94 -180.26   360.53                         
## carabid.model  7 371.83 391.79 -178.91   357.83 2.6996      3     0.4403
emmeans(carabid.model,specs=pairwise~Cover:Fertilizer,type='response')
## $emmeans
##  Cover Fertilizer  rate    SE  df lower.CL upper.CL
##  Bare  Inorganic  1.306 0.535 121    0.580     2.94
##  Cover Inorganic  1.568 0.665 121    0.676     3.63
##  Bare  Manure     0.881 0.389 121    0.368     2.11
##  Cover Manure     1.062 0.442 121    0.466     2.42
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                         ratio    SE  df t.ratio p.value
##  Bare,Inorganic / Cover,Inorganic 0.833 0.261 121 -0.584  0.9367 
##  Bare,Inorganic / Bare,Manure     1.482 0.514 121  1.136  0.6681 
##  Bare,Inorganic / Cover,Manure    1.229 0.386 121  0.657  0.9127 
##  Cover,Inorganic / Bare,Manure    1.780 0.626 121  1.639  0.3606 
##  Cover,Inorganic / Cover,Manure   1.476 0.482 121  1.192  0.6330 
##  Bare,Manure / Cover,Manure       0.829 0.294 121 -0.528  0.9520 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log scale
#Edge effects
carabid.model=glmmTMB(Total~Edge+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2018"),ziformula=~1,family=poisson)
      summary(carabid.model)
##  Family: poisson  ( log )
## Formula:          Total ~ Edge + (1 | Month.Order) + (1 | Plot)
## Zero inflation:         ~1
## Data: subset(carabid.data, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    362.9    380.1   -175.5    350.9      122 
## 
## Random effects:
## 
## Conditional model:
##  Groups      Name        Variance Std.Dev.
##  Month.Order (Intercept) 0.42745  0.6538  
##  Plot        (Intercept) 0.02194  0.1481  
## Number of obs: 128, groups:  Month.Order, 4; Plot, 16
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.4993     0.4492  -1.111 0.266377    
## Edge1         0.5604     0.3157   1.775 0.075850 .  
## Edge2         1.0969     0.3297   3.327 0.000879 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.4154     0.4633  -3.055  0.00225 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      Anova(carabid.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Total
##               Chisq Df Pr(>Chisq)   
## (Intercept)  1.2353  1   0.266377   
## Edge        12.7334  2   0.001718 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
carabid.null=glmmTMB(Total~1+(1|Month.Order)+(1|Plot),data=subset(carabid.data,Year=="2018"),ziformula=~1,family=poisson)
anova(carabid.model,carabid.null,test="Chisq")
## Data: subset(carabid.data, Year == "2018")
## Models:
## carabid.null: Total ~ 1 + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
## carabid.model: Total ~ Edge + (1 | Month.Order) + (1 | Plot), zi=~1, disp=~1
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## carabid.null   4 368.53 379.94 -180.26   360.53                            
## carabid.model  6 362.94 380.05 -175.47   350.94 9.5873      2   0.008282 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(carabid.model,specs=pairwise~Edge,type='response')
## $emmeans
##  Edge  rate    SE  df lower.CL upper.CL
##  0    0.607 0.273 122    0.249     1.48
##  1    1.063 0.394 122    0.511     2.21
##  2    1.818 0.679 122    0.867     3.81
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast ratio    SE  df t.ratio p.value
##  0 / 1    0.571 0.180 122 -1.775  0.1822 
##  0 / 2    0.334 0.110 122 -3.327  0.0033 
##  1 / 2    0.585 0.126 122 -2.488  0.0375 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

3) Plot Carabids

Make figures for Carabid- looking to adjust 2016 to be half the size of 2017 and 2018

## Warning: Removed 2 rows containing non-finite values (stat_ydensity).
## Warning: Removed 2 rows containing non-finite values (stat_summary).

Save plot as “FigPerCarabids.eps” in figs folder

ggsave(plot,file="./figs/FigCarbid.eps",width=16,height=10,units='in')
## Warning: Removed 2 rows containing non-finite values (stat_ydensity).
## Warning: Removed 2 rows containing non-finite values (stat_summary).

Edge effects on carabids

## Warning: Removed 2 rows containing non-finite values (stat_ydensity).
## Warning: Removed 2 rows containing non-finite values (stat_summary).

Save plot as “FigCarabidsEdge.eps” in figs folder

ggsave(plot,file="./figs/FigCarbidEdge.eps",width=16,height=10,units='in')
## Warning: Removed 2 rows containing non-finite values (stat_ydensity).
## Warning: Removed 2 rows containing non-finite values (stat_summary).

Make figures for Carabid- Looking at cover crop effect

does not exist

## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 60 rows containing missing values (geom_point).

Save plot as “FigPerCarabids.eps” in figs folder

ggsave(plot,file="./figs/FigCarbid.eps",width=16,height=10,units='in')
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 60 rows containing missing values (geom_point).

I want to explore what pitfall trap capture looks like in the field using heat maps

Manure Carabid diversity analysis

Make data set with total sum of carabid by species over the season

macro.gr=ddply(carabid.data,.(Year, Plot, Fertilizer, Cover), summarize, Total=sum(Total), Pterasticus.melanarius = sum(Pterasticus.melanarius),
               Pterasticus.mutus=sum(Pterasticus.mutus), Harpalus.pensylvanicus=sum(Harpalus.pensylvanicus),Harpalus.affinis=sum(Harpalus.affinis),
               Harpalus.herbivagus=sum(Harpalus.herbivagus),Patrobus.longicornus=sum(Patrobus.longicornus),Chlaenius.tricolor=sum(Chlaenius.tricolor),
               Poecilus.chalcites=sum(Poecilus.chalcites),Poecilus.lucublandus=sum(Poecilus.lucublandus),Notiobia.sayi=sum(Notiobia.sayi),
               Bembidion.quadrimaculatum=sum(Bembidion.quadrimaculatum),Amara.neoscotica=sum(Amara.neoscotica),Anisodactylus.sanctaecrucis=sum(Anisodactylus.sanctaecrucis),
               Anisodactylus.harrisii=sum(Anisodactylus.harrisii),Dyschrius.globulosus=sum(Dyschrius.globulosus),Colliuris.pensylvanica=sum(Colliuris.pensylvanica),
               Scarites.quadriceps=sum(Scarites.quadriceps),Notiophilus.biguttaus=sum(Notiophilus.biguttaus),Agonum.placidum=sum(Agonum.placidum),
               Cicindela.punctulata=sum(Cicindela.punctulata),
               Stenolophus.conjunctus=sum(Stenolophus.conjunctus),Dicaelus.elongatus=sum(Dicaelus.elongatus))
length(macro.gr[,1])
## [1] 40
shan=diversity(macro.gr[,6:27])
S <- specnumber(macro.gr[,6:27])
J <- shan/log(S)
carabid.div=cbind(macro.gr,shan,S,J)

Shannon-Wiener index

All years

diverse.model=lmer(shan~Fertilizer+(1|Year), data=subset(carabid.div))
summary(diverse.model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: shan ~ Fertilizer + (1 | Year)
##    Data: subset(carabid.div)
## 
## REML criterion at convergence: 38
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86998 -0.78393  0.05583  0.70455  1.92285 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Year     (Intercept) 0.1660   0.4075  
##  Residual             0.1166   0.3415  
## Number of obs: 40, groups:  Year, 3
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)        0.8212     0.2479  2.2607   3.312   0.0679 .
## FertilizerManure   0.2023     0.1080 36.0529   1.873   0.0692 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## FertilzrMnr -0.218
Anova(diverse.model,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: shan
##               Chisq Df Pr(>Chisq)    
## (Intercept) 10.9685  1  0.0009267 ***
## Fertilizer   3.5088  1  0.0610421 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model)

diverse.model.null=lm(shan~1, data=subset(carabid.div))
anova(diverse.model,diverse.model.null)
## refitting model(s) with ML (instead of REML)
## Data: subset(carabid.div)
## Models:
## diverse.model.null: shan ~ 1
## diverse.model: shan ~ Fertilizer + (1 | Year)
##                    Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## diverse.model.null  2 62.344 65.722 -29.172   58.344                        
## diverse.model       4 42.185 48.940 -17.092   34.185 24.16      2  5.673e-06
##                       
## diverse.model.null    
## diverse.model      ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(diverse.model,specs=pairwise~Fertilizer)
## $emmeans
##  Fertilizer emmean    SE  df lower.CL upper.CL
##  Inorganic   0.821 0.248 2.2  -0.1576      1.8
##  Manure      1.023 0.248 2.2   0.0446      2.0
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.202 0.108 36 -1.873  0.0692 
## 
## Degrees-of-freedom method: kenward-roger

2016

diverse.model.2016=lm(shan~Fertilizer, data=subset(carabid.div,Year=="2016"))
summary(diverse.model.2016)
## 
## Call:
## lm(formula = shan ~ Fertilizer, data = subset(carabid.div, Year == 
##     "2016"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6004 -0.1298  0.0364  0.2081  0.4185 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        0.6497     0.1739   3.736  0.00967 **
## FertilizerManure   0.4743     0.2459   1.929  0.10204   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3478 on 6 degrees of freedom
## Multiple R-squared:  0.3827, Adjusted R-squared:  0.2798 
## F-statistic: 3.719 on 1 and 6 DF,  p-value: 0.102
Anova(diverse.model.2016,type="III")
## Anova Table (Type III tests)
## 
## Response: shan
##              Sum Sq Df F value  Pr(>F)   
## (Intercept) 1.68828  1 13.9566 0.00967 **
## Fertilizer  0.44993  1  3.7195 0.10204   
## Residuals   0.72580  6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2016)

diverse.model.null=lm(shan~1, data=subset(carabid.div,Year=="2016"))
anova(diverse.model.2016,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: shan ~ Fertilizer
## Model 2: shan ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      6 0.7258                           
## 2      7 1.1757 -1  -0.44993 3.7195  0.102
emmeans(diverse.model.2016,specs=pairwise~Fertilizer)
## $emmeans
##  Fertilizer emmean    SE df lower.CL upper.CL
##  Inorganic    0.65 0.174  6    0.224     1.08
##  Manure       1.12 0.174  6    0.698     1.55
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.474 0.246  6 -1.929  0.1020

2017

diverse.model.2017=lm(shan~Fertilizer*Cover, data=subset(carabid.div,Year=="2017"))
summary(diverse.model.2017)
## 
## Call:
## lm(formula = shan ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2017"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31737 -0.19164 -0.01163  0.15092  0.47636 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                   0.3174     0.1337   2.374   0.0351 *
## FertilizerManure              0.3601     0.1891   1.904   0.0811 .
## CoverCover                    0.1228     0.1891   0.649   0.5283  
## FertilizerManure:CoverCover  -0.1275     0.2674  -0.477   0.6420  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2674 on 12 degrees of freedom
## Multiple R-squared:  0.3078, Adjusted R-squared:  0.1347 
## F-statistic: 1.778 on 3 and 12 DF,  p-value: 0.2047
Anova(diverse.model.2017,type="III")
## Anova Table (Type III tests)
## 
## Response: shan
##                   Sum Sq Df F value  Pr(>F)  
## (Intercept)      0.40290  1  5.6360 0.03514 *
## Fertilizer       0.25928  1  3.6270 0.08110 .
## Cover            0.03015  1  0.4217 0.52833  
## Fertilizer:Cover 0.01626  1  0.2274 0.64204  
## Residuals        0.85785 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2017)

diverse.model.null=lm(shan~1, data=subset(carabid.div,Year=="2017"))
anova(diverse.model.2017,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: shan ~ Fertilizer * Cover
## Model 2: shan ~ 1
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1     12 0.85785                           
## 2     15 1.23924 -3  -0.38139 1.7783 0.2047
emmeans(diverse.model.2017,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean     SE df lower.CL upper.CL
##  Inorganic   0.379 0.0945 12    0.173    0.585
##  Manure      0.675 0.0945 12    0.469    0.881
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.296 0.134 12 -2.216  0.0467 
## 
## Results are averaged over the levels of: Cover

2018

diverse.model.2018=lm(shan~Fertilizer*Cover, data=subset(carabid.div,Year=="2018"))
summary(diverse.model.2018)
## 
## Call:
## lm(formula = shan ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2018"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6283 -0.2691  0.1525  0.2876  0.5092 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   1.5079     0.2026   7.443  7.8e-06 ***
## FertilizerManure             -0.2567     0.2865  -0.896    0.388    
## CoverCover                   -0.2844     0.2865  -0.993    0.340    
## FertilizerManure:CoverCover   0.4578     0.4052   1.130    0.281    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4052 on 12 degrees of freedom
## Multiple R-squared:  0.1025, Adjusted R-squared:  -0.1219 
## F-statistic: 0.4569 on 3 and 12 DF,  p-value: 0.7174
Anova(diverse.model.2018,type="III")
## Anova Table (Type III tests)
## 
## Response: shan
##                  Sum Sq Df F value    Pr(>F)    
## (Intercept)      9.0946  1 55.4036 7.805e-06 ***
## Fertilizer       0.1318  1  0.8029    0.3879    
## Cover            0.1618  1  0.9855    0.3404    
## Fertilizer:Cover 0.2096  1  1.2769    0.2806    
## Residuals        1.9698 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2018)

diverse.model.null=lm(shan~1, data=subset(carabid.div,Year=="2018"))
anova(diverse.model.2018,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: shan ~ Fertilizer * Cover
## Model 2: shan ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     12 1.9698                           
## 2     15 2.1948 -3  -0.22501 0.4569 0.7174
emmeans(diverse.model.2018,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean    SE df lower.CL upper.CL
##  Inorganic    1.37 0.143 12     1.05     1.68
##  Manure       1.34 0.143 12     1.03     1.65
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   0.0278 0.203 12 0.137   0.8932 
## 
## Results are averaged over the levels of: Cover

Richness

2016

diverse.model.2016=lm(S~Fertilizer, data=subset(carabid.div,Year=="2016"))
summary(diverse.model.2016)
## 
## Call:
## lm(formula = S ~ Fertilizer, data = subset(carabid.div, Year == 
##     "2016"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##     -2     -1      0      1      2 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)        3.0000     0.7071   4.243  0.00542 **
## FertilizerManure   4.0000     1.0000   4.000  0.00712 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.414 on 6 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.6818 
## F-statistic:    16 on 1 and 6 DF,  p-value: 0.007119
Anova(diverse.model.2016,type="III")
## Anova Table (Type III tests)
## 
## Response: S
##             Sum Sq Df F value   Pr(>F)   
## (Intercept)     36  1      18 0.005424 **
## Fertilizer      32  1      16 0.007119 **
## Residuals       12  6                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2016)

diverse.model.null=lm(S~1, data=subset(carabid.div,Year=="2016"))
anova(diverse.model.2016,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: S ~ Fertilizer
## Model 2: S ~ 1
##   Res.Df RSS Df Sum of Sq  F   Pr(>F)   
## 1      6  12                            
## 2      7  44 -1       -32 16 0.007119 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(diverse.model.2016,specs=pairwise~Fertilizer)
## $emmeans
##  Fertilizer emmean    SE df lower.CL upper.CL
##  Inorganic       3 0.707  6     1.27     4.73
##  Manure          7 0.707  6     5.27     8.73
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate SE df t.ratio p.value
##  Inorganic - Manure       -4  1  6 -4.000  0.0071

2017

diverse.model.2017=lm(S~Fertilizer*Cover, data=subset(carabid.div,Year=="2017"))
summary(diverse.model.2017)
## 
## Call:
## lm(formula = S ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2017"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.250 -0.500 -0.375  0.500  1.750 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   2.2500     0.4507   4.992 0.000313 ***
## FertilizerManure              0.2500     0.6374   0.392 0.701764    
## CoverCover                    0.2500     0.6374   0.392 0.701764    
## FertilizerManure:CoverCover  -0.2500     0.9014  -0.277 0.786231    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9014 on 12 degrees of freedom
## Multiple R-squared:  0.01887,    Adjusted R-squared:  -0.2264 
## F-statistic: 0.07692 on 3 and 12 DF,  p-value: 0.9713
Anova(diverse.model.2017,type="III")
## Anova Table (Type III tests)
## 
## Response: S
##                   Sum Sq Df F value    Pr(>F)    
## (Intercept)      20.2500  1 24.9231 0.0003133 ***
## Fertilizer        0.1250  1  0.1538 0.7017639    
## Cover             0.1250  1  0.1538 0.7017639    
## Fertilizer:Cover  0.0625  1  0.0769 0.7862314    
## Residuals         9.7500 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2017)

diverse.model.null=lm(S~1, data=subset(carabid.div,Year=="2017"))
anova(diverse.model.2017,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: S ~ Fertilizer * Cover
## Model 2: S ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     12 9.7500                           
## 2     15 9.9375 -3   -0.1875 0.0769 0.9713
emmeans(diverse.model.2017,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean    SE df lower.CL upper.CL
##  Inorganic    2.38 0.319 12     1.68     3.07
##  Manure       2.50 0.319 12     1.81     3.19
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.125 0.451 12 -0.277  0.7862 
## 
## Results are averaged over the levels of: Cover

2018

diverse.model.2018=lm(S~Fertilizer*Cover, data=subset(carabid.div,Year=="2018"))
summary(diverse.model.2018)
## 
## Call:
## lm(formula = S ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2018"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.750 -1.125  0.250  1.000  2.250 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    5.500      0.863   6.373 3.54e-05 ***
## FertilizerManure              -1.500      1.220  -1.229    0.243    
## CoverCover                    -0.750      1.220  -0.615    0.550    
## FertilizerManure:CoverCover    1.750      1.726   1.014    0.331    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.726 on 12 degrees of freedom
## Multiple R-squared:  0.1159, Adjusted R-squared:  -0.1051 
## F-statistic: 0.5245 on 3 and 12 DF,  p-value: 0.6737
Anova(diverse.model.2018,type="III")
## Anova Table (Type III tests)
## 
## Response: S
##                   Sum Sq Df F value    Pr(>F)    
## (Intercept)      121.000  1 40.6154 3.541e-05 ***
## Fertilizer         4.500  1  1.5105    0.2426    
## Cover              1.125  1  0.3776    0.5504    
## Fertilizer:Cover   3.062  1  1.0280    0.3306    
## Residuals         35.750 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2018)

diverse.model.null=lm(S~1, data=subset(carabid.div,Year=="2018"))
anova(diverse.model.2018,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: S ~ Fertilizer * Cover
## Model 2: S ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     12 35.750                           
## 2     15 40.438 -3   -4.6875 0.5245 0.6737
emmeans(diverse.model.2018,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean   SE df lower.CL upper.CL
##  Inorganic    5.12 0.61 12     3.80     6.45
##  Manure       4.50 0.61 12     3.17     5.83
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure    0.625 0.863 12 0.724   0.4828 
## 
## Results are averaged over the levels of: Cover

Evenness

2016

diverse.model.2016=lm(J~Fertilizer, data=subset(carabid.div,Year=="2016"))
summary(diverse.model.2016)
## 
## Call:
## lm(formula = J ~ Fertilizer, data = subset(carabid.div, Year == 
##     "2016"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28268 -0.08504  0.02167  0.10854  0.28178 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.63652    0.10179   6.253 0.000776 ***
## FertilizerManure -0.06752    0.14395  -0.469 0.655578    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2036 on 6 degrees of freedom
## Multiple R-squared:  0.03537,    Adjusted R-squared:  -0.1254 
## F-statistic:  0.22 on 1 and 6 DF,  p-value: 0.6556
Anova(diverse.model.2016,type="III")
## Anova Table (Type III tests)
## 
## Response: J
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 1.62063  1  39.103 0.0007756 ***
## Fertilizer  0.00912  1   0.220 0.6555777    
## Residuals   0.24867  6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2016)

diverse.model.null=lm(J~1, data=subset(carabid.div,Year=="2016"))
anova(diverse.model.2016,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: J ~ Fertilizer
## Model 2: J ~ 1
##   Res.Df     RSS Df  Sum of Sq    F Pr(>F)
## 1      6 0.24867                          
## 2      7 0.25779 -1 -0.0091192 0.22 0.6556
emmeans(diverse.model.2016,specs=pairwise~Fertilizer)
## $emmeans
##  Fertilizer emmean    SE df lower.CL upper.CL
##  Inorganic   0.637 0.102  6    0.387    0.886
##  Manure      0.569 0.102  6    0.320    0.818
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   0.0675 0.144  6 0.469   0.6556

2017

diverse.model.2017=lm(J~Fertilizer*Cover, data=subset(carabid.div,Year=="2017"))
summary(diverse.model.2017)
## 
## Call:
## lm(formula = J ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2017"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28147 -0.18557  0.03183  0.13509  0.40666 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                  0.41965    0.13496   3.110  0.00993 **
## FertilizerManure             0.29729    0.17853   1.665  0.12407   
## CoverCover                   0.14464    0.17853   0.810  0.43501   
## FertilizerManure:CoverCover -0.08213    0.24329  -0.338  0.74204   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2337 on 11 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3095, Adjusted R-squared:  0.1212 
## F-statistic: 1.644 on 3 and 11 DF,  p-value: 0.236
Anova(diverse.model.2017,type="III")
## Anova Table (Type III tests)
## 
## Response: J
##                   Sum Sq Df F value   Pr(>F)   
## (Intercept)      0.52831  1  9.6691 0.009934 **
## Fertilizer       0.15151  1  2.7729 0.124066   
## Cover            0.03587  1  0.6564 0.435006   
## Fertilizer:Cover 0.00623  1  0.1140 0.742041   
## Residuals        0.60103 11                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2017)

diverse.model.null=lm(J~1, data=subset(carabid.div,Year=="2017"))
anova(diverse.model.2017,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: J ~ Fertilizer * Cover
## Model 2: J ~ 1
##   Res.Df     RSS Df Sum of Sq      F Pr(>F)
## 1     11 0.60103                           
## 2     14 0.87043 -3  -0.26941 1.6436  0.236
emmeans(diverse.model.2017,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean     SE df lower.CL upper.CL
##  Inorganic   0.492 0.0893 11    0.295    0.688
##  Manure      0.748 0.0826 11    0.566    0.930
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE df t.ratio p.value
##  Inorganic - Manure   -0.256 0.122 11 -2.106  0.0590 
## 
## Results are averaged over the levels of: Cover

2018

diverse.model.2018=lm(J~Fertilizer*Cover, data=subset(carabid.div,Year=="2018"))
summary(diverse.model.2018)
## 
## Call:
## lm(formula = J ~ Fertilizer * Cover, data = subset(carabid.div, 
##     Year == "2018"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22275 -0.02227  0.01520  0.08018  0.11471 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  0.88976    0.05437  16.363 1.43e-09 ***
## FertilizerManure             0.03136    0.07690   0.408    0.691    
## CoverCover                  -0.00446    0.07690  -0.058    0.955    
## FertilizerManure:CoverCover -0.02681    0.10875  -0.247    0.809    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1087 on 12 degrees of freedom
## Multiple R-squared:  0.02263,    Adjusted R-squared:  -0.2217 
## F-statistic: 0.09261 on 3 and 12 DF,  p-value: 0.9627
Anova(diverse.model.2018,type="III")
## Anova Table (Type III tests)
## 
## Response: J
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)      3.1667  1 267.7630 1.432e-09 ***
## Fertilizer       0.0020  1   0.1663    0.6906    
## Cover            0.0000  1   0.0034    0.9547    
## Fertilizer:Cover 0.0007  1   0.0608    0.8094    
## Residuals        0.1419 12                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(diverse.model.2018)

diverse.model.null=lm(J~1, data=subset(carabid.div,Year=="2018"))
anova(diverse.model.2018,diverse.model.null)
## Analysis of Variance Table
## 
## Model 1: J ~ Fertilizer * Cover
## Model 2: J ~ 1
##   Res.Df     RSS Df  Sum of Sq      F Pr(>F)
## 1     12 0.14192                            
## 2     15 0.14520 -3 -0.0032856 0.0926 0.9627
emmeans(diverse.model.2018,specs=pairwise~Fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  Fertilizer emmean     SE df lower.CL upper.CL
##  Inorganic   0.888 0.0384 12    0.804    0.971
##  Manure      0.905 0.0384 12    0.822    0.989
## 
## Results are averaged over the levels of: Cover 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate     SE df t.ratio p.value
##  Inorganic - Manure   -0.018 0.0544 12 -0.330  0.7469 
## 
## Results are averaged over the levels of: Cover

Community Analysis PCA Predator community

make data set

 pitfall.average=ddply(pitfall.data,.(Year,Plot,Fertilizer,Cover),summarize,carabids=mean(carabid.beetles,na.rm=TRUE),
                            rove=mean(rove.beetles,na.rm=TRUE),
                            wolf.spiders=mean(wolf.spiders,na.rm=TRUE),
                            other.spiders=mean(other.spiders,na.rm=TRUE),
                            harvestmen=mean(harvestmen,na.rm=TRUE),
                            centipedes=mean(centipedes,na.rm=TRUE),
                            ants=mean(ants,na.rm=TRUE))
      mygroups=c("carabids","rove","wolf.spiders","other.spiders","harvestmen",'centipedes','ants')

Transform dataset: using hellinger transformation (From Numerial ecology with R, p 140)

separate everything by year

Start with 2016

varespec=subset(pitfall.average,Year=="2016")[mygroups]
spe.h=decostand(varespec,'hellinger')

Make PCA

PCs=prcomp(spe.h,scale.=TRUE)
PCs$rotation
##                       PC1         PC2         PC3         PC4         PC5
## carabids       0.56230120 -0.06417411  0.03244215 -0.24176593 -0.04172368
## rove          -0.03077411 -0.41613608 -0.61551820  0.62572971  0.09699622
## wolf.spiders  -0.55151053 -0.15248383  0.13819252 -0.09653732 -0.31650993
## other.spiders  0.03596392  0.54949068  0.31685288  0.66843884 -0.32791762
## harvestmen    -0.41413147  0.34692972  0.01150996 -0.01133036  0.75508376
## centipedes     0.31368923 -0.33611180  0.55792289  0.29958189  0.44555275
## ants           0.32788296  0.51399392 -0.43494999 -0.06327672  0.11168541
##                      PC6         PC7
## carabids      -0.5554375 -0.55672623
## rove          -0.1208941 -0.17736896
## wolf.spiders   0.2341756 -0.69939362
## other.spiders -0.1837393 -0.09094105
## harvestmen    -0.3169356 -0.19306581
## centipedes     0.3988140 -0.17329401
## ants           0.5733888 -0.30638085

create variables for graphing with ggplot

PC1=as.numeric(as.character(PCs$x[,1]))
PC2=as.numeric(as.character(PCs$x[,2]))
PC3=as.numeric(as.character(PCs$x[,3]))
plot.data=data.frame(PC1,PC2,PC3,Fertilizer=subset(pitfall.average,Year=="2016")$Fertilizer)

Do a PERMANOVA

adonis(spe.h ~ Fertilizer, data = plot.data, method='eu')
## 
## Call:
## adonis(formula = spe.h ~ Fertilizer, data = plot.data, method = "eu") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs  MeanSqs F.Model     R2 Pr(>F)
## Fertilizer  1  0.082647 0.082647  2.2611 0.2737  0.176
## Residuals   6  0.219310 0.036552         0.7263       
## Total       7  0.301957                  1.0000

Build plot 1

p.2016=autoplot(PCs, data=plot.data, colour='Fertilizer',shape="Fertilizer",size=3,loadings=TRUE,loadings.colour = 'black',
                 loadings.label = TRUE, loadings.label.size = 5,loadings.label.colour='black',loadings.label.repel=T,frame=T,frame.color='Fertilizer',frame.alpha=1)
p.2016=p.2016+theme_classic()+scale_color_manual(values=fert.colors)+scale_fill_manual(values=fert.colors)+theme(legend.position='none')
p.2016=p.2016+theme(axis.text.x=element_text(size=20,angle=0,vjust=0.5),
                axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),
                axis.text.y=element_text(size=20))+ggtitle('2016')
p.2016 = p.2016 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_rect(colour = "black", size=1))
p.2016

do 2017

varespec=subset(pitfall.average,Year=="2017")[mygroups]
spe.h=decostand(varespec,'hellinger')

Make PCA

PCs=prcomp(spe.h,scale.=TRUE)
PCs$rotation
##                       PC1         PC2         PC3        PC4         PC5
## carabids       0.55990034 -0.07351944  0.02979476 -0.3780767  0.61018342
## rove          -0.49996566  0.40548150 -0.24655598  0.1349323  0.20255616
## wolf.spiders  -0.22750033  0.18014248  0.56541779 -0.6145584 -0.35260361
## other.spiders  0.16732180 -0.42498622 -0.51764396 -0.2106331 -0.56330977
## harvestmen    -0.38530601 -0.54046096  0.05600392  0.1941280  0.21692624
## centipedes     0.06230622 -0.42072718  0.57381292  0.3558423 -0.05392307
## ants           0.45216167  0.38471501  0.13517459  0.5024874 -0.30827798
##                       PC6       PC7
## carabids       0.05367983 0.4025810
## rove          -0.23881080 0.6391893
## wolf.spiders   0.19210680 0.2391206
## other.spiders -0.07851040 0.3944449
## harvestmen     0.65726383 0.1989130
## centipedes    -0.53058300 0.2807055
## ants           0.42837474 0.3134297
summary(PCs)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.4646 1.2854 1.2507 0.81185 0.74125 0.56732 0.32879
## Proportion of Variance 0.3064 0.2360 0.2235 0.09416 0.07849 0.04598 0.01544
## Cumulative Proportion  0.3064 0.5425 0.7659 0.86009 0.93858 0.98456 1.00000

create variables for graphing with ggplot

PC1=as.numeric(as.character(PCs$x[,1]))
PC2=as.numeric(as.character(PCs$x[,2]))
PC3=as.numeric(as.character(PCs$x[,3]))
treatment=paste(Fertilizer=subset(pitfall.average,Year=="2017")$Fertilizer,Cover=subset(pitfall.average,Year=="2017")$Cover)
plot.data=data.frame(PC1,PC2,PC3,Fertilizer=subset(pitfall.average,Year=="2017")$Fertilizer,Cover=subset(pitfall.average,Year=="2017")$Cover,treatment)

Do a PERMANOVA

adonis(spe.h ~ Fertilizer*Cover, data = plot.data, method='eu')
## 
## Call:
## adonis(formula = spe.h ~ Fertilizer * Cover, data = plot.data,      method = "eu") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                  Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Fertilizer        1    0.6224 0.62238 2.99656 0.18085  0.010 **
## Cover             1    0.0923 0.09225 0.44417 0.02681  0.896   
## Fertilizer:Cover  1    0.2345 0.23445 1.12880 0.06813  0.343   
## Residuals        12    2.4924 0.20770         0.72422          
## Total            15    3.4415                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Build plot 2

p.2017=autoplot(PCs, data=plot.data, colour='Fertilizer',shape="treatment",size=3,loadings=TRUE,loadings.colour = 'black',
                 loadings.label = TRUE, loadings.label.size = 5,loadings.label.colour='black',loadings.label.repel=T,frame=T,frame.color='Fertilizer',frame.alpha=1)
p.2017=p.2017+theme_classic()+scale_color_manual(values=fert.colors)+scale_fill_manual(values=fert.colors)+theme(legend.position='none')
p.2017=p.2017+theme(axis.text.x=element_text(size=20,angle=0,vjust=0.5),
                axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),
                axis.text.y=element_text(size=20))+ggtitle('2017')
p.2017 = p.2017 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_rect(colour = "black", size=1))
p.2017

do 2018

varespec=subset(pitfall.average,Year=="2018")[mygroups]
spe.h=decostand(varespec,'hellinger')

Make PCA

PCs=prcomp(spe.h,scale.=TRUE)
PCs$rotation
##                       PC1          PC2         PC3         PC4         PC5
## carabids      -0.46345911 -0.006338239  0.46398859  0.14466028 -0.64954778
## rove          -0.01682015 -0.642438342  0.30897752  0.17948780  0.42220596
## wolf.spiders   0.60646544  0.093909879 -0.11930263  0.41919276 -0.08510473
## other.spiders  0.38636080 -0.130539541  0.05890843 -0.80122834 -0.28336367
## harvestmen    -0.20932745  0.584059761  0.31080485 -0.25121234  0.53495359
## centipedes    -0.31805197 -0.453558298 -0.27352057 -0.25585912  0.14476457
## ants          -0.35053202  0.120549705 -0.70719941  0.02442485 -0.07183238
##                       PC6       PC7
## carabids      -0.09113014 0.3445815
## rove           0.37932566 0.3703529
## wolf.spiders  -0.35525411 0.5477012
## other.spiders  0.16540577 0.2838902
## harvestmen    -0.13061421 0.3900173
## centipedes    -0.71897644 0.1223377
## ants           0.40031937 0.4432291

create variables for graphing with ggplot

PC1=as.numeric(as.character(PCs$x[,1]))
PC2=as.numeric(as.character(PCs$x[,2]))
PC3=as.numeric(as.character(PCs$x[,3]))
treatment=paste(Fertilizer=subset(pitfall.average,Year=="2018")$Fertilizer,Cover=subset(pitfall.average,Year=="2018")$Cover)
plot.data=data.frame(PC1,PC2,PC3,Fertilizer=subset(pitfall.average,Year=="2018")$Fertilizer,Cover=subset(pitfall.average,Year=="2018")$Cover,treatment)

Do a PERMANOVA

adonis(spe.h ~ Fertilizer*Cover, data = plot.data, method='eu')
## 
## Call:
## adonis(formula = spe.h ~ Fertilizer * Cover, data = plot.data,      method = "eu") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                  Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## Fertilizer        1   0.07440 0.074398 1.18607 0.08133  0.333
## Cover             1   0.06011 0.060108 0.95826 0.06571  0.461
## Fertilizer:Cover  1   0.02757 0.027571 0.43955 0.03014  0.867
## Residuals        12   0.75272 0.062727         0.82283       
## Total            15   0.91480                  1.00000

Build plot 2

p.2018=autoplot(PCs, data=plot.data, colour='Fertilizer',shape="treatment",size=3,loadings=TRUE,loadings.colour = 'black',
                 loadings.label = TRUE, loadings.label.size = 5,loadings.label.colour='black',loadings.label.repel=T,frame=T,frame.color='Fertilizer',frame.alpha=1)
p.2018=p.2018+theme_classic()+scale_color_manual(values=fert.colors)+scale_fill_manual(values=fert.colors)
p.2018=p.2018+theme(axis.text.x=element_text(size=20,angle=0,vjust=0.5),
                axis.title.x=element_text(size=20),axis.title.y=element_text(size=20),
                axis.text.y=element_text(size=20))+ggtitle('2018')
p.2018= p.2018 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_rect(colour = "black", size=1))
p.2018

Put all the plots together

p=p.2016+p.2017+p.2018
ggsave(p,file="./figs/FigPCA.eps",width=16,height=5,units='in')

Sentinel Predation

1) Import data

sentinel=read.csv("./data/raw/SentinelPrey.csv",header=TRUE)
sentinel.LONG <- gather(sentinel, waxworm, pred, WAXWORM.1:WAXWORM.2, factor_key=TRUE)

2) Run stats

Use binomial distribution for all of 2017 together

model0=glmmTMB(pred~1+(1|Month.Order)+(1|Plot/Location),family='binomial',data=subset(sentinel.LONG,Year=="2017"))
model2017=glmmTMB(pred~fertilizer*cover+(1|Month.Order)+(1|Plot/Location),family='binomial',data=subset(sentinel.LONG,Year=="2017"))
summary(model2017)
##  Family: binomial  ( logit )
## Formula:          
## pred ~ fertilizer * cover + (1 | Month.Order) + (1 | Plot/Location)
## Data: subset(sentinel.LONG, Year == "2017")
## 
##      AIC      BIC   logLik deviance df.resid 
##    273.2    298.7   -129.6    259.2      278 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Month.Order   (Intercept) 6.790e-01 8.240e-01
##  Location:Plot (Intercept) 1.086e+00 1.042e+00
##  Plot          (Intercept) 8.841e-09 9.402e-05
## Number of obs: 285, groups:  Month.Order, 3; Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        1.3581     0.6534   2.079   0.0377 *
## fertilizermanure                   0.4692     0.6334   0.741   0.4589  
## covercover crop                    0.8308     0.6594   1.260   0.2077  
## fertilizermanure:covercover crop  -0.6137     0.9322  -0.658   0.5104  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model2017,type='III')
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)  
## (Intercept)      4.3203  1    0.03766 *
## fertilizer       0.5487  1    0.45886  
## cover            1.5873  1    0.20771  
## fertilizer:cover 0.4333  1    0.51035  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model0,model2017,test='Chisq')
## Data: subset(sentinel.LONG, Year == "2017")
## Models:
## model0: pred ~ 1 + (1 | Month.Order) + (1 | Plot/Location), zi=~0, disp=~1
## model2017: pred ~ fertilizer * cover + (1 | Month.Order) + (1 | Plot/Location), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model0     4 269.09 283.70 -130.55   261.09                         
## model2017  7 273.17 298.74 -129.59   259.17 1.9171      3     0.5898
emmeans(model2017,specs=pairwise~fertilizer)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  fertilizer emmean    SE  df lower.CL upper.CL
##  inorganic    1.77 0.599 278    0.594     2.95
##  manure       1.94 0.599 278    0.758     3.11
## 
## Results are averaged over the levels of: cover 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast           estimate    SE  df t.ratio p.value
##  inorganic - manure   -0.162 0.464 278 -0.350  0.7265 
## 
## Results are averaged over the levels of: cover 
## Results are given on the log odds ratio (not the response) scale.

JULY 2017

july2017=subset(sentinel.LONG,Month=="July"&Year=="2017")

model.jul.2017=glmmTMB(pred~fertilizer*cover+(1|Plot/Location),family=binomial,data=july2017)
summary(model.jul.2017)
##  Family: binomial  ( logit )
## Formula:          pred ~ fertilizer * cover + (1 | Plot/Location)
## Data: july2017
## 
##      AIC      BIC   logLik deviance df.resid 
##    101.9    117.3    -45.0     89.9       89 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Location:Plot (Intercept) 2.600e+02 16.124373
##  Plot          (Intercept) 2.227e-06  0.001492
## Number of obs: 95, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        6.8608     3.0942   2.217   0.0266 *
## fertilizermanure                   1.5293     3.1898   0.479   0.6316  
## covercover crop                    1.5291     3.1897   0.479   0.6317  
## fertilizermanure:covercover crop  -0.9162     4.3734  -0.210   0.8341  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.jul.2017,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)  
## (Intercept)      4.9165  1     0.0266 *
## fertilizer       0.2299  1     0.6316  
## cover            0.2298  1     0.6317  
## fertilizer:cover 0.0439  1     0.8341  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.jul.2017.null=glmmTMB(pred~1+(1|Plot/Location),family=binomial,data=july2017)
anova(model.jul.2017,model.jul.2017.null,test='Chisq')
## Data: july2017
## Models:
## model.jul.2017.null: pred ~ 1 + (1 | Plot/Location), zi=~0, disp=~1
## model.jul.2017: pred ~ fertilizer * cover + (1 | Plot/Location), zi=~0, disp=~1
##                     Df     AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model.jul.2017.null  3  96.431 104.09 -45.216   90.431                         
## model.jul.2017       6 101.932 117.25 -44.966   89.932 0.4994      3      0.919
emmeans(model.jul.2017,specs=pairwise~fertilizer:cover,type='response')
## $emmeans
##  fertilizer cover        prob        SE df lower.CL upper.CL
##  inorganic  bare       0.9990 0.0032361 89   0.6710        1
##  manure     bare       0.9998 0.0005731 89   0.9669        1
##  inorganic  cover crop 0.9998 0.0005732 89   0.9669        1
##  manure     cover crop 0.9999 0.0003161 89   0.9801        1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio    SE df t.ratio p.value
##  inorganic,bare / manure,bare                  0.217 0.691 89 -0.479  0.9635 
##  inorganic,bare / inorganic,cover crop         0.217 0.691 89 -0.479  0.9635 
##  inorganic,bare / manure,cover crop            0.117 0.387 89 -0.649  0.9156 
##  manure,bare / inorganic,cover crop            1.000 2.914 89  0.000  1.0000 
##  manure,bare / manure,cover crop               0.542 1.632 89 -0.203  0.9970 
##  inorganic,cover crop / manure,cover crop      0.542 1.632 89 -0.203  0.9970 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

AUGUST 2017

aug2017=subset(sentinel.LONG,Month=="August"&Year=="2017")

model.aug.2017=glmmTMB(pred~cover*fertilizer+(1|Plot/Location),family=binomial,data=aug2017)
summary(model.aug.2017)
##  Family: binomial  ( logit )
## Formula:          pred ~ cover * fertilizer + (1 | Plot/Location)
## Data: aug2017
## 
##      AIC      BIC   logLik deviance df.resid 
##     91.7    107.0    -39.9     79.7       88 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Location:Plot (Intercept) 1.147e+00 1.071e+00
##  Plot          (Intercept) 4.780e-09 6.913e-05
## Number of obs: 94, groups:  Location:Plot, 47; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        1.6390     0.7479   2.192   0.0284 *
## covercover crop                    1.2539     1.1060   1.134   0.2569  
## fertilizermanure                   0.5930     0.9830   0.603   0.5464  
## covercover crop:fertilizermanure  -1.8246     1.4970  -1.219   0.2229  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.aug.2017,Type='III')
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)
## cover            0.1060  1     0.7447
## fertilizer       0.0951  1     0.7578
## cover:fertilizer 1.4856  1     0.2229
model.aug.2017.null=glmmTMB(pred~1+(1|Plot/Location),family=binomial,data=aug2017)
anova(model.aug.2017,model.aug.2017.null,test="Chisq")
## Data: aug2017
## Models:
## model.aug.2017.null: pred ~ 1 + (1 | Plot/Location), zi=~0, disp=~1
## model.aug.2017: pred ~ cover * fertilizer + (1 | Plot/Location), zi=~0, disp=~1
##                     Df    AIC     BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model.aug.2017.null  3 87.691  95.321 -40.845   81.691                         
## model.aug.2017       6 91.735 106.995 -39.867   79.735 1.9559      3     0.5816
emmeans(model.aug.2017,specs=pairwise~fertilizer:cover,type='response')
## $emmeans
##  fertilizer cover       prob     SE df lower.CL upper.CL
##  inorganic  bare       0.837 0.1018 88    0.538    0.958
##  manure     bare       0.903 0.0796 88    0.604    0.983
##  inorganic  cover crop 0.947 0.0553 88    0.665    0.994
##  manure     cover crop 0.840 0.1039 88    0.530    0.961
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio    SE df t.ratio p.value
##  inorganic,bare / manure,bare                  0.553 0.543 88 -0.603  0.9308 
##  inorganic,bare / inorganic,cover crop         0.285 0.316 88 -1.134  0.6697 
##  inorganic,bare / manure,cover crop            0.978 0.879 88 -0.025  1.0000 
##  manure,bare / inorganic,cover crop            0.516 0.591 88 -0.577  0.9386 
##  manure,bare / manure,cover crop               1.770 1.737 88  0.581  0.9375 
##  inorganic,cover crop / manure,cover crop      3.427 3.762 88  1.122  0.6770 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

SEPTEMBER 2017

sept2017=subset(sentinel.LONG,Month=="September"&Year=="2017")

model.sept.2017=glmmTMB(pred~cover*fertilizer+(1|Plot/Location),family=binomial,data=sept2017)
summary(model.sept.2017)
##  Family: binomial  ( logit )
## Formula:          pred ~ cover * fertilizer + (1 | Plot/Location)
## Data: sept2017
## 
##      AIC      BIC   logLik deviance df.resid 
##     61.0     76.4    -24.5     49.0       90 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev.
##  Location:Plot (Intercept) 4.027e+01 6.346136
##  Plot          (Intercept) 6.811e-08 0.000261
## Number of obs: 96, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                        7.7642     2.9348   2.646  0.00816 **
## covercover crop                   -0.4575     3.0629  -0.149  0.88127   
## fertilizermanure                  -1.3438     2.8239  -0.476  0.63418   
## covercover crop:fertilizermanure   1.0142     3.8981   0.260  0.79472   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.sept.2017)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)
## cover            0.0079  1     0.9291
## fertilizer       0.1736  1     0.6769
## cover:fertilizer 0.0677  1     0.7947
emmeans(model.sept.2017,specs=pairwise~fertilizer:cover,type='response')
## $emmeans
##  fertilizer cover        prob       SE df lower.CL upper.CL
##  inorganic  bare       0.9996 0.001245 90   0.8737        1
##  manure     bare       0.9984 0.004042 90   0.8133        1
##  inorganic  cover crop 0.9993 0.001783 90   0.8830        1
##  manure     cover crop 0.9991 0.002396 90   0.8658        1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio    SE df t.ratio p.value
##  inorganic,bare / manure,bare                  3.833 10.83 90  0.476  0.9642 
##  inorganic,bare / inorganic,cover crop         1.580  4.84 90  0.149  0.9988 
##  inorganic,bare / manure,cover crop            2.197  6.45 90  0.268  0.9932 
##  manure,bare / inorganic,cover crop            0.412  1.06 90 -0.345  0.9858 
##  manure,bare / manure,cover crop               0.573  1.38 90 -0.231  0.9956 
##  inorganic,cover crop / manure,cover crop      1.390  3.74 90  0.122  0.9993 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

Use binomial distribution for all of 2018 together

model0=glmer(pred~(1|Plot),family='binomial',data=subset(sentinel.LONG,Year=="2018"))
model2018=glmer(pred~fertilizer*Month+(1|Plot/Location),family='binomial',data=subset(sentinel.LONG,Year=="2018"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0143922 (tol = 0.001, component 1)
summary(model2018)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: pred ~ fertilizer * Month + (1 | Plot/Location)
##    Data: subset(sentinel.LONG, Year == "2018")
## 
##      AIC      BIC   logLik deviance df.resid 
##    433.8    473.3   -206.9    413.8      374 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2471 -0.8341  0.3167  0.6627  1.3329 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Location:Plot (Intercept) 0.4103   0.6406  
##  Plot          (Intercept) 0.0787   0.2805  
## Number of obs: 384, groups:  Location:Plot, 48; Plot, 16
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       2.9077     0.6318   4.602 4.18e-06 ***
## fertilizermanure                 -1.2908     0.7506  -1.720 0.085502 .  
## MonthJuly                        -2.6295     0.6808  -3.862 0.000112 ***
## MonthJune                        -2.5346     0.6808  -3.723 0.000197 ***
## MonthSeptember                   -0.5738     0.7676  -0.747 0.454783    
## fertilizermanure:MonthJuly        1.0150     0.8343   1.217 0.223747    
## fertilizermanure:MonthJune        1.5912     0.8394   1.896 0.058010 .  
## fertilizermanure:MonthSeptember  -0.1624     0.9168  -0.177 0.859390    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                  (Intr) frtlzr MnthJl MnthJn MnthSp frtlzrmnr:MnthJl
## fertilzrmnr      -0.827                                             
## MonthJuly        -0.862  0.714                                      
## MonthJune        -0.861  0.713  0.796                               
## MonthSptmbr      -0.744  0.625  0.690  0.690                        
## frtlzrmnr:MnthJl  0.691 -0.813 -0.806 -0.640 -0.562                 
## frtlzrmnr:MnthJn  0.691 -0.808 -0.640 -0.806 -0.559  0.726          
## frtlzrmn:MS       0.618 -0.733 -0.574 -0.574 -0.837  0.660          
##                  frtlzrmnr:MnthJn
## fertilzrmnr                      
## MonthJuly                        
## MonthJune                        
## MonthSptmbr                      
## frtlzrmnr:MnthJl                 
## frtlzrmnr:MnthJn                 
## frtlzrmn:MS       0.655          
## convergence code: 0
## Model failed to converge with max|grad| = 0.0143922 (tol = 0.001, component 1)
Anova(model2018)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                    Chisq Df Pr(>Chisq)    
## fertilizer        1.5842  1     0.2082    
## Month            28.9252  3  2.322e-06 ***
## fertilizer:Month  7.3369  3     0.0619 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model0,model2018,test='Chisq')
## Data: subset(sentinel.LONG, Year == "2018")
## Models:
## model0: pred ~ (1 | Plot)
## model2018: pred ~ fertilizer * Month + (1 | Plot/Location)
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## model0     2 467.65 475.55 -231.82   463.65                             
## model2018 10 433.83 473.34 -206.91   413.83 49.816      8  4.432e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(model2018,specs=pairwise~fertilizer|Month)
## $emmeans
## Month = August:
##  fertilizer  emmean    SE  df asymp.LCL asymp.UCL
##  inorganic  2.90771 0.632 Inf    1.6694     4.146
##  manure     1.61693 0.422 Inf    0.7896     2.444
## 
## Month = July:
##  fertilizer  emmean    SE  df asymp.LCL asymp.UCL
##  inorganic  0.27825 0.348 Inf   -0.4041     0.961
##  manure     0.00249 0.347 Inf   -0.6773     0.682
## 
## Month = June:
##  fertilizer  emmean    SE  df asymp.LCL asymp.UCL
##  inorganic  0.37311 0.350 Inf   -0.3125     1.059
##  manure     0.67351 0.359 Inf   -0.0301     1.377
## 
## Month = September:
##  fertilizer  emmean    SE  df asymp.LCL asymp.UCL
##  inorganic  2.33394 0.517 Inf    1.3214     3.346
##  manure     0.88075 0.368 Inf    0.1595     1.602
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
## Month = August:
##  contrast           estimate    SE  df z.ratio p.value
##  inorganic - manure    1.291 0.751 Inf  1.720  0.0855 
## 
## Month = July:
##  contrast           estimate    SE  df z.ratio p.value
##  inorganic - manure    0.276 0.491 Inf  0.561  0.5747 
## 
## Month = June:
##  contrast           estimate    SE  df z.ratio p.value
##  inorganic - manure   -0.300 0.500 Inf -0.601  0.5480 
## 
## Month = September:
##  contrast           estimate    SE  df z.ratio p.value
##  inorganic - manure    1.453 0.628 Inf  2.313  0.0207 
## 
## Results are given on the log odds ratio (not the response) scale.

JUNE 2018

june2018=subset(sentinel.LONG,Month=="June"&Year=="2018")

model.jun.2018=glmmTMB(pred~cover*fertilizer+(1|Plot/Location),family=binomial,data=june2018)
summary(model.jun.2018)
##  Family: binomial  ( logit )
## Formula:          pred ~ cover * fertilizer + (1 | Plot/Location)
## Data: june2018
## 
##      AIC      BIC   logLik deviance df.resid 
##    137.0    152.3    -62.5    125.0       90 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Location:Plot (Intercept) 1.062e-08 1.031e-04
##  Plot          (Intercept) 1.623e-09 4.029e-05
## Number of obs: 96, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)                       1.671e-01  4.097e-01   0.408    0.683
## covercover crop                   3.438e-01  5.879e-01   0.585    0.559
## fertilizermanure                 -1.188e-08  5.794e-01   0.000    1.000
## covercover crop:fertilizermanure  5.878e-01  8.577e-01   0.685    0.493
Anova(model.jun.2018)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)
## cover            2.0971  1     0.1476
## fertilizer       0.3941  1     0.5302
## cover:fertilizer 0.4696  1     0.4932
emmeans(model.jun.2018,specs=pairwise~fertilizer:cover,type='response')
## $emmeans
##  fertilizer cover       prob     SE df lower.CL upper.CL
##  inorganic  bare       0.542 0.1017 90    0.344    0.727
##  manure     bare       0.542 0.1017 90    0.344    0.727
##  inorganic  cover crop 0.625 0.0988 90    0.419    0.794
##  manure     cover crop 0.750 0.0884 90    0.540    0.884
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio    SE df t.ratio p.value
##  inorganic,bare / manure,bare                  1.000 0.579 90  0.000  1.0000 
##  inorganic,bare / inorganic,cover crop         0.709 0.417 90 -0.585  0.9364 
##  inorganic,bare / manure,cover crop            0.394 0.246 90 -1.492  0.4468 
##  manure,bare / inorganic,cover crop            0.709 0.417 90 -0.585  0.9364 
##  manure,bare / manure,cover crop               0.394 0.246 90 -1.492  0.4468 
##  inorganic,cover crop / manure,cover crop      0.556 0.351 90 -0.929  0.7892 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

JULY 2018 = cover by manure interaction

july2018=subset(sentinel.LONG,Month=="July"&Year=="2018")

model.jul.2018=glmmTMB(pred~cover*fertilizer+(1|Plot/Location),family=binomial,data=july2018)
summary(model.jul.2018)
##  Family: binomial  ( logit )
## Formula:          pred ~ cover * fertilizer + (1 | Plot/Location)
## Data: july2018
## 
##      AIC      BIC   logLik deviance df.resid 
##    136.5    151.9    -62.3    124.5       90 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Location:Plot (Intercept) 9.121e-01 9.550e-01
##  Plot          (Intercept) 7.240e-09 8.509e-05
## Number of obs: 96, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        0.6268     0.5564   1.127   0.2599  
## covercover crop                   -0.6268     0.7690  -0.815   0.4150  
## fertilizermanure                  -1.4624     0.8152  -1.794   0.0728 .
## covercover crop:fertilizermanure   2.3059     1.1635   1.982   0.0475 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.jul.2018)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)  
## cover            0.6421  1    0.42296  
## fertilizer       0.2316  1    0.63037  
## cover:fertilizer 3.9279  1    0.04749 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(model.jul.2018,specs=pairwise~cover:fertilizer,type='response')
## $emmeans
##  cover      fertilizer  prob    SE df lower.CL upper.CL
##  bare       inorganic  0.652 0.126 90    0.383    0.850
##  cover crop inorganic  0.500 0.133 90    0.258    0.742
##  bare       manure     0.302 0.120 90    0.123    0.572
##  cover crop manure     0.699 0.120 90    0.427    0.879
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio    SE df t.ratio p.value
##  bare,inorganic / cover crop,inorganic         1.872 1.439 90  0.815  0.8472 
##  bare,inorganic / bare,manure                  4.316 3.519 90  1.794  0.2830 
##  bare,inorganic / cover crop,manure            0.805 0.624 90 -0.280  0.9923 
##  cover crop,inorganic / bare,manure            2.306 1.791 90  1.076  0.7051 
##  cover crop,inorganic / cover crop,manure      0.430 0.336 90 -1.080  0.7024 
##  bare,manure / cover crop,manure               0.187 0.155 90 -2.018  0.1891 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

AUGUST 2018

aug2018=subset(sentinel.LONG,Month=="August"&Year=="2018")

model.aug.2018=glmmTMB(pred~cover*fertilizer+(1|Plot/Location),family=binomial,data=aug2018)
summary(model.aug.2018)
##  Family: binomial  ( logit )
## Formula:          pred ~ cover * fertilizer + (1 | Plot/Location)
## Data: aug2018
## 
##      AIC      BIC   logLik deviance df.resid 
##     69.0     84.4    -28.5     57.0       90 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance Std.Dev.
##  Location:Plot (Intercept) 63.049   7.940   
##  Plot          (Intercept)  4.997   2.235   
## Number of obs: 96, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                        8.7827     3.3101   2.653  0.00797 **
## covercover crop                   -0.6913     3.5577  -0.194  0.84594   
## fertilizermanure                  -2.1593     3.3208  -0.650  0.51553   
## covercover crop:fertilizermanure   1.8760     4.4333   0.423  0.67217   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.aug.2018)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##                   Chisq Df Pr(>Chisq)
## cover            0.0575  1     0.8105
## fertilizer       0.2529  1     0.6150
## cover:fertilizer 0.1791  1     0.6722
emmeans(model.aug.2018,specs=pairwise~fertilizer:cover,type='response')
## $emmeans
##  fertilizer cover        prob        SE df lower.CL upper.CL
##  inorganic  bare       0.9998 0.0005075 90   0.9008        1
##  manure     bare       0.9987 0.0034568 90   0.8088        1
##  inorganic  cover crop 0.9997 0.0008737 90   0.9182        1
##  manure     cover crop 0.9996 0.0011134 90   0.9138        1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                 odds.ratio     SE df t.ratio p.value
##  inorganic,bare / manure,bare                  8.665 28.775 90  0.650  0.9152 
##  inorganic,bare / inorganic,cover crop         1.996  7.102 90  0.194  0.9974 
##  inorganic,bare / manure,cover crop            2.650  9.181 90  0.281  0.9922 
##  manure,bare / inorganic,cover crop            0.230  0.642 90 -0.527  0.9524 
##  manure,bare / manure,cover crop               0.306  0.814 90 -0.445  0.9704 
##  inorganic,cover crop / manure,cover crop      1.328  3.903 90  0.096  0.9997 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

SEPTEMBER 2018

sept2018=subset(sentinel.LONG,Month=="September"&Year=="2018")

model.sept.2018=glmmTMB(pred~fertilizer+(1|Plot/Location),family=binomial,data=sept2018)
summary(model.sept.2018)
##  Family: binomial  ( logit )
## Formula:          pred ~ fertilizer + (1 | Plot/Location)
## Data: sept2018
## 
##      AIC      BIC   logLik deviance df.resid 
##     95.8    106.0    -43.9     87.8       92 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev. 
##  Location:Plot (Intercept) 3.941e+00 1.9852726
##  Plot          (Intercept) 1.425e-08 0.0001194
## Number of obs: 96, groups:  Location:Plot, 48; Plot, 16
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)         3.408      1.273   2.677  0.00744 **
## fertilizermanure   -2.043      1.074  -1.903  0.05706 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(model.sept.2018)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: pred
##             Chisq Df Pr(>Chisq)  
## fertilizer 3.6207  1    0.05706 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.sept.2018.null=glmmTMB(pred~1+(1|Plot/Location),family=binomial,data=sept2018)
anova(model.sept.2018,model.sept.2018.null)
## Data: sept2018
## Models:
## model.sept.2018.null: pred ~ 1 + (1 | Plot/Location), zi=~0, disp=~1
## model.sept.2018: pred ~ fertilizer + (1 | Plot/Location), zi=~0, disp=~1
##                      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model.sept.2018.null  3 98.792 106.48 -46.396   92.792                         
## model.sept.2018       4 95.773 106.03 -43.887   87.773 5.0182      1    0.02508
##                       
## model.sept.2018.null  
## model.sept.2018      *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(model.sept.2018,specs=pairwise~fertilizer,type='response')
## $emmeans
##  fertilizer  prob     SE df lower.CL upper.CL
##  inorganic  0.968 0.0395 92    0.707    0.997
##  manure     0.796 0.1258 92    0.456    0.948
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast           odds.ratio   SE df t.ratio p.value
##  inorganic / manure       7.72 8.29 92 1.903   0.0602 
## 
## Tests are performed on the log odds ratio scale

3) Plot sentinel prey

Summarize predation rates

prey.means = ddply(sentinel.LONG,.(Year,Month,Month.Order,fertilizer,cover,Plot),summarize, mean=mean(pred,na.rm=TRUE), 
                   N=length(pred),
                   sd=sd(pred,na.rm=TRUE),
                   se=sd/sqrt(N))

Plot

plot=ggplot(prey.means,aes(x=reorder(Month,Month.Order),y=mean,fill=paste(fertilizer,cover)))
plot=plot+geom_boxplot()#+geom_jitter(width=.1)
plot=plot+scale_fill_manual(values=manure.color)
#plot=plot + stat_summary(fun = mean, fun.min = mean, fun.max = mean, geom = "crossbar",  fatten=2, width=.75,color="black")
plot=plot+facet_grid(.~Year)
plot=plot+theme_bw()+ylim(0,1.2)
plot = plot+ylab("% Eaten")

plot=plot+theme(axis.text.x=element_text(size=20,angle=0,vjust=0.5),axis.title.x=element_blank(),axis.title.y=element_text(size=20),axis.text.y=element_text(size=20)) 
plot=plot+theme(strip.placement = "outside",strip.background = element_blank(),strip.text.x=element_text(size=20),strip.text.y=element_text(size=20))
plot = plot + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_rect(color = "black", size=1))

plot

# export plot

ggsave(plot,file="./figs/FigSentinelPrey.eps",width=16,height=8,units='in')