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
#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.
## 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).
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
## 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).
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
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')
We think that edges have an effect on predator abundance, particularly carabids, so adding information about edges to dataset
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
Save plot as “FigPerCover.eps” in figs folder
ggsave(plot,file="./figs/FigHarvestmen.eps",width=16,height=10,units='in')
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
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
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).
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
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
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
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
## 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).
## 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).
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
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)
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
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
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
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')
separate everything by year
Start with 2016
varespec=subset(pitfall.average,Year=="2016")[mygroups]
spe.h=decostand(varespec,'hellinger')
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')
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')
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
p=p.2016+p.2017+p.2018
ggsave(p,file="./figs/FigPCA.eps",width=16,height=5,units='in')
sentinel=read.csv("./data/raw/SentinelPrey.csv",header=TRUE)
sentinel.LONG <- gather(sentinel, waxworm, pred, WAXWORM.1:WAXWORM.2, factor_key=TRUE)
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
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')