NoMowMay Summary of Results

No Mow May 2020:

Who Participated?

• 435 Appleton Homes in 2020, >2000 in 2021, 81 Cities in 2023 across the USA
• Bee diversity and abundance data collected at 20 homes and 15 mowed urban green spaces in 2020
• we can extrapolate to the area and say that we protected between 19 and 24 acres of pollinator habiat across the city in 2020 but unable to reliably estimate for 2021

##   NoMowCities State Year Participants       Lon      Lat
## 1    Appleton    WI 2020           NA -88.40725 44.26180
## 2     Oshkosh    WI 2021           NA -88.54001 44.02080
## 3      Wausau    WI 2021           NA -89.62026 44.92619

What is the average size of a NoMow Lawn Participant?

• The average participant had a yard of 195 sq meters or 2100 sq. ft!

Were there more floral resourcess in NoMow Lawns?

• Mowed Areas had 35% fewer species than NoMow lawns
• Mowed areas had 64% less area covered by angiosperms than No Mow areas.

#Floral Richness
FR<- ggplot (NoMow, aes (x=Mow_NoMow, y=Floral_Richness, fill=Mow_NoMow)) + 
  geom_boxplot(alpha=0.5, color="black") + theme_bw() + ylim(0,15) +
  scale_fill_manual (values=c("black", "olivedrab3"))  +
  ylab("Floral Richness") + xlab("") + 
  annotate ("text", x=1, y=13, label="A", size=16) + theme (legend.position = "none") + 
  scale_x_discrete(labels=c("Mow" = "Regular Mowing", "NoMow" = "Not Mowed"))
ggplotly (FR)

#Median value for No Mow was 7, median value for mowed was 4.5. Percent increase calculation:
# ((7-4.5)/4.5)*100 = 55.5% so a 55% increase compared from mowed to No Mow

# These results are signigicant at the p<0.001. 
kruskal.test(Floral_Richness ~ Mow_NoMow, data=NoMow)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Floral_Richness by Mow_NoMow
## Kruskal-Wallis chi-squared = 13.095, df = 1, p-value = 0.000296

#Floral Density

FD<- ggplot (NoMow, aes (x=Mow_NoMow, y=Mean_Floral_Density, fill=Mow_NoMow)) +
  geom_boxplot(alpha=0.5, color="black") + theme_bw() +
  scale_fill_manual (values=c("black", "olivedrab3"))  +
  ylab("Floral Percent Cover") + xlab("") + 
  annotate ("text", x=1, y=80, label="B", size=16) + theme (legend.position = "none") + 
  scale_x_discrete(labels=c("Mow" = "Regular Mowing", "NoMow" = "Not Mowed"))
ggplotly (FD)

# Median value 18.7 for mowed, 52 for no mowed. 

# These results are signigicant at the p<0.0001. 
kruskal.test(Mean_Floral_Density ~ Mow_NoMow, data=NoMow)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Mean_Floral_Density by Mow_NoMow
## Kruskal-Wallis chi-squared = 16.828, df = 1, p-value = 4.092e-05

Are there more bees in NoMow homes compared to mowed greenspaces?

• Yes. There are 70% more bees in non-mowed lawns than in the mowed urban green spaces

summary(NoMow)

##      Date            AddressCode        lon              lat       
##  Length:35          Min.   : 1.0   Min.   :-88.42   Min.   :44.25  
##  Class :character   1st Qu.: 9.5   1st Qu.:-88.40   1st Qu.:44.26  
##  Mode  :character   Median :18.0   Median :-88.39   Median :44.27  
##                     Mean   :18.0   Mean   :-88.39   Mean   :44.27  
##                     3rd Qu.:26.5   3rd Qu.:-88.38   3rd Qu.:44.27  
##                     Max.   :35.0   Max.   :-88.37   Max.   :44.29  
##    NoMowArea      Mow_NoMow             Time              Duration     
##  Min.   : 91.0   Length:35          Length:35          Min.   : 27.00  
##  1st Qu.:150.0   Class :character   Class :character   1st Qu.: 45.00  
##  Median :153.0   Mode  :character   Mode  :character   Median : 45.00  
##  Mean   :188.2                                         Mean   : 54.94  
##  3rd Qu.:213.0                                         3rd Qu.: 64.00  
##  Max.   :446.0                                         Max.   :107.00  
##  Floral_Density_1 Floral_Density_2 Floral_Density_3 Floral_Density_4
##  Min.   : 0.00    Min.   :  0.00   Min.   :  0.00   Min.   : 0.00   
##  1st Qu.: 5.00    1st Qu.: 25.00   1st Qu.: 15.00   1st Qu.: 7.50   
##  Median :40.00    Median : 50.00   Median : 40.00   Median :40.00   
##  Mean   :34.77    Mean   : 47.03   Mean   : 40.69   Mean   :39.66   
##  3rd Qu.:55.00    3rd Qu.: 77.50   3rd Qu.: 60.00   3rd Qu.:70.00   
##  Max.   :95.00    Max.   :100.00   Max.   :100.00   Max.   :90.00   
##  Floral_Density_5 Mean_Floral_Density Floral_Richness  Bee_Abundance   
##  Min.   :  0.00   Min.   : 1.00       Min.   : 2.000   Min.   : 1.000  
##  1st Qu.: 15.00   1st Qu.:22.00       1st Qu.: 5.000   1st Qu.: 3.000  
##  Median : 35.00   Median :45.00       Median : 7.000   Median : 8.000  
##  Mean   : 41.63   Mean   :40.80       Mean   : 6.571   Mean   : 8.543  
##  3rd Qu.: 67.50   3rd Qu.:56.50       3rd Qu.: 8.000   3rd Qu.:11.000  
##  Max.   :100.00   Max.   :88.75       Max.   :12.000   Max.   :24.000  
##   Bee_Richness  
##  Min.   :1.000  
##  1st Qu.:2.000  
##  Median :4.000  
##  Mean   :3.886  
##  3rd Qu.:5.500  
##  Max.   :7.000

bee.abund <- ggplot (NoMow, aes (x=Mow_NoMow, y=Bee_Abundance, fill=Mow_NoMow)) +
    geom_boxplot (alpha=0.5, color="black") +  theme_bw() +xlab("") + ylab ("Bee Abundance") + 
    scale_fill_manual (values=c("black", "#FFCC00")) + 
    annotate ("text", x=.75, y=22, label="C", size=16) + theme (legend.position = "none") + 
    scale_x_discrete(labels=c("Mow" = "Regular Mowing", "NoMow" = "Not Mowed"))

ggplotly (bee.abund)

#median value 3 for mowed, 10 for no mowed

# These results are signigicant at the p<0.0001. 
kruskal.test(Bee_Abundance ~ Mow_NoMow, data=NoMow)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Bee_Abundance by Mow_NoMow
## Kruskal-Wallis chi-squared = 19.135, df = 1, p-value = 1.218e-05

Are there more genera of bees in NoMow homes compared to mowed greenspaces?

• There are 60% more bee genera in No Mow spaces compared to mowed spaces

bee.rich <- ggplot (NoMow, aes (x=Mow_NoMow, y=Bee_Richness, fill=Mow_NoMow)) +
  geom_boxplot (alpha=0.5, color="black") +  theme_bw() +xlab("") + ylab ("Bee Genera Richness") + ylim(0,10)+
  scale_fill_manual (values=c("black", "#FFCC00")) + 
  annotate ("text", x=.75, y=9, label="D", size=16) + theme (legend.position = "none") + 
  scale_x_discrete(labels=c("Mow" = "Regular Mowing", "NoMow" = "Not Mowed"))
ggplotly (bee.rich)

#median value for mowed is 2, no mow is 5.
# These results are signigicant at the p<0.001. 
kruskal.test(Bee_Richness ~ Mow_NoMow, data=NoMow)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Bee_Richness by Mow_NoMow
## Kruskal-Wallis chi-squared = 14.458, df = 1, p-value = 0.0001434

What are the variables that best predict increases in abundance of bees?

• Area unmowed is the best predictor of abundances of bees
• Floral species richness is also a good predictor of bee abundance in unmowed (city parks)
  

A quick check for multi-colinearity

library (corrgram)
library (tidyverse)
NoMowPreds<-NoMow %>% dplyr::select (NoMowArea, Mean_Floral_Density, Floral_Richness); head (NoMowPreds)

##   NoMowArea Mean_Floral_Density Floral_Richness
## 1    157.32               88.75               8
## 2    145.30               19.60               5
## 3    190.58               55.00               7
## 4    153.00               58.00               7
## 5     93.00               51.00               8
## 6    150.00               15.20               4

corrgram (NoMowPreds, upper.panel = panel.cor)

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

#note all variables are scaled to make them comparable 
NoMow_GLM_abund<-glm (Bee_Abundance ~ Mow_NoMow + scale (NoMowArea) + scale(Mean_Floral_Density) + scale (Floral_Richness) , family="poisson", data=NoMow)
#The AIC of the global model is 194.56 with the classification of lawn type considered as well as the total area.
summary (NoMow_GLM_abund)

## 
## Call:
## glm(formula = Bee_Abundance ~ Mow_NoMow + scale(NoMowArea) + 
##     scale(Mean_Floral_Density) + scale(Floral_Richness), family = "poisson", 
##     data = NoMow)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 1.16912    0.18889   6.189 6.04e-10 ***
## Mow_NoMowNoMow              1.20521    0.22502   5.356 8.51e-08 ***
## scale(NoMowArea)            0.20716    0.05348   3.873 0.000107 ***
## scale(Mean_Floral_Density)  0.02726    0.08818   0.309 0.757183    
## scale(Floral_Richness)     -0.07819    0.07912  -0.988 0.322987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 146.770  on 34  degrees of freedom
## Residual deviance:  54.228  on 30  degrees of freedom
## AIC: 194.56
## 
## Number of Fisher Scoring iterations: 4

#reducing the model shows that area sampled might be the best predictor of bee abundance 
stepAIC (NoMow_GLM_abund)

## Start:  AIC=194.56
## Bee_Abundance ~ Mow_NoMow + scale(NoMowArea) + scale(Mean_Floral_Density) + 
##     scale(Floral_Richness)
## 
##                              Df Deviance    AIC
## - scale(Mean_Floral_Density)  1   54.323 192.66
## - scale(Floral_Richness)      1   55.217 193.55
## <none>                            54.228 194.56
## - scale(NoMowArea)            1   67.962 206.30
## - Mow_NoMow                   1   84.649 222.99
## 
## Step:  AIC=192.66
## Bee_Abundance ~ Mow_NoMow + scale(NoMowArea) + scale(Floral_Richness)
## 
##                          Df Deviance    AIC
## - scale(Floral_Richness)  1   55.223 191.56
## <none>                        54.323 192.66
## - scale(NoMowArea)        1   68.697 205.03
## - Mow_NoMow               1   98.548 234.88
## 
## Step:  AIC=191.56
## Bee_Abundance ~ Mow_NoMow + scale(NoMowArea)
## 
##                    Df Deviance    AIC
## <none>                  55.223 191.56
## - scale(NoMowArea)  1   68.754 203.09
## - Mow_NoMow         1  104.075 238.41

## 
## Call:  glm(formula = Bee_Abundance ~ Mow_NoMow + scale(NoMowArea), family = "poisson", 
##     data = NoMow)
## 
## Coefficients:
##      (Intercept)    Mow_NoMowNoMow  scale(NoMowArea)  
##           1.1972            1.1641            0.1876  
## 
## Degrees of Freedom: 34 Total (i.e. Null);  32 Residual
## Null Deviance:       146.8 
## Residual Deviance: 55.22     AIC: 191.6

area.abund<- ggplot (NoMow, aes (x=NoMowArea, y=Bee_Abundance, color=Mow_NoMow)) + geom_point() + 
  geom_smooth(method="lm") + theme_bw() + scale_color_manual(values=c("grey", "black")) + 
   xlab("Size of Now Mow Area in sq. meters") + ylab ("Bee Abundance")

ggplotly (area.abund)

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

floral.abund<-ggplot (NoMow, aes (x=Floral_Richness, y=Bee_Abundance, color=Mow_NoMow)) + 
  geom_point() +  geom_smooth(method="loess") + theme_bw() + 
  scale_color_manual(values=c("grey","black")) + xlab("Number of Flowering Species in Lawn") + 
  ylab ("Bee Abundance")

ggplotly (floral.abund)

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

What are the variables that best predict increases in species richness of bee genera?

• Area sampled is the best predictor of bee genera richness.

NoMow_GLM_rich<-glm (Bee_Richness ~ Mow_NoMow + scale (NoMowArea) + scale(Mean_Floral_Density) + 
                       scale (Floral_Richness), family="poisson", data=NoMow) 

summary (NoMow_GLM_rich); stepAIC (NoMow_GLM_rich)

## 
## Call:
## glm(formula = Bee_Richness ~ Mow_NoMow + scale(NoMowArea) + scale(Mean_Floral_Density) + 
##     scale(Floral_Richness), family = "poisson", data = NoMow)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 1.07359    0.23435   4.581 4.63e-06 ***
## Mow_NoMowNoMow              0.31369    0.31038   1.011   0.3122    
## scale(NoMowArea)            0.15246    0.08568   1.779   0.0752 .  
## scale(Mean_Floral_Density)  0.09998    0.13581   0.736   0.4616    
## scale(Floral_Richness)      0.11553    0.11407   1.013   0.3112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 33.348  on 34  degrees of freedom
## Residual deviance: 12.928  on 30  degrees of freedom
## AIC: 131.83
## 
## Number of Fisher Scoring iterations: 4

## Start:  AIC=131.83
## Bee_Richness ~ Mow_NoMow + scale(NoMowArea) + scale(Mean_Floral_Density) + 
##     scale(Floral_Richness)
## 
##                              Df Deviance    AIC
## - scale(Mean_Floral_Density)  1   13.471 130.37
## - scale(Floral_Richness)      1   13.940 130.84
## - Mow_NoMow                   1   13.949 130.85
## <none>                            12.928 131.83
## - scale(NoMowArea)            1   15.883 132.79
## 
## Step:  AIC=130.37
## Bee_Richness ~ Mow_NoMow + scale(NoMowArea) + scale(Floral_Richness)
## 
##                          Df Deviance    AIC
## - scale(Floral_Richness)  1   15.290 130.19
## <none>                        13.471 130.37
## - scale(NoMowArea)        1   15.938 130.84
## - Mow_NoMow               1   16.422 131.32
## 
## Step:  AIC=130.19
## Bee_Richness ~ Mow_NoMow + scale(NoMowArea)
## 
##                    Df Deviance    AIC
## <none>                  15.290 130.19
## - scale(NoMowArea)  1   19.569 132.47
## - Mow_NoMow         1   22.513 135.42

## 
## Call:  glm(formula = Bee_Richness ~ Mow_NoMow + scale(NoMowArea), family = "poisson", 
##     data = NoMow)
## 
## Coefficients:
##      (Intercept)    Mow_NoMowNoMow  scale(NoMowArea)  
##           0.8978            0.5943            0.1653  
## 
## Degrees of Freedom: 34 Total (i.e. Null);  32 Residual
## Null Deviance:       33.35 
## Residual Deviance: 15.29     AIC: 130.2

#best fit reduced model
glm(formula = Bee_Richness ~ Mow_NoMow + NoMowArea, family = "poisson", 
    data = NoMow)

## 
## Call:  glm(formula = Bee_Richness ~ Mow_NoMow + NoMowArea, family = "poisson", 
##     data = NoMow)
## 
## Coefficients:
##    (Intercept)  Mow_NoMowNoMow       NoMowArea  
##       0.469417        0.594296        0.002277  
## 
## Degrees of Freedom: 34 Total (i.e. Null);  32 Residual
## Null Deviance:       33.35 
## Residual Deviance: 15.29     AIC: 130.2

area.rich<- ggplot (NoMow, aes (x=NoMowArea, y=Bee_Richness, color=Mow_NoMow)) + geom_point() + 
  geom_smooth(method="lm") + theme_bw() + scale_color_manual(values=c("grey", "black")) + 
   xlab("Size of Now Mow Area in sq. meters") + ylab ("Bee Genera Richness")

ggplotly (area.rich)

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

No Mow May 2021 Data Analyses

Are there more flowers in No Mow Yards during the four weeks of sampling?

• Appleton and Oshkosh peak in floral abundance in plots during week 2 of NoMow May. Wausau continously increases from weeks 2-4 though with a lot of variablity.
• when data are pooled between cities there are consistently more flowers in no mow plots relative to mowed plots

# How many participants? 39 unique paired sites across three cities
# 
 col= data.summary<-data %>% group_by(City, Week, Treatment) %>%
  summarize (mean_flowers= mean(`Number of flowers`), se_flowers= sd(`Number of flowers`)/sqrt(n())) 

## `summarise()` has grouped output by 'City', 'Week'. You can override using the
## `.groups` argument.

# 
 data.summary2<-data %>% group_by(City, Week,Treatment) %>% 
   summarize (mean_flowers= mean(`Number of flowers`), se_flowers= sd(`Number of flowers`)/sqrt(n())) 

## `summarise()` has grouped output by 'City', 'Week'. You can override using the
## `.groups` argument.

# 
 plot1<- ggplot (data.summary, aes(color=City)) + geom_point(aes (x=Week, y=mean_flowers)) + 
   geom_line(aes (x=Week, y=data.summary$mean_flowers)) + theme_bw() + 
   geom_errorbar(aes(x=Week, ymin=mean_flowers - se_flowers, ymax= mean_flowers + se_flowers, width=.1)) +
   facet_grid(.~Treatment) + ylab("Mean number of flowers per plot") +  scale_color_brewer(palette = "Accent")
 #plot1
 
 ggplotly(plot1)

#treatment/week Boxplot 
plot2<-ggplot (data, aes(x=Treatment, y=log (`Number of flowers`))) + geom_boxplot() +
  facet_grid(.~Week) + theme_bw() + ylab("Number of Flowers in Plot (log)")

ggplotly (plot2)

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

#wilcoxon test by week 
# Week 1 test 
week1<- data %>% filter (Week==1)
wilcox.test (week1$`Number of flowers`~ week1$Treatment)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  week1$`Number of flowers` by week1$Treatment
## W = 1321.5, p-value = 1.151e-06
## alternative hypothesis: true location shift is not equal to 0

# Week 2 test 
week2<- data %>% filter (Week==2)
wilcox.test (week2$`Number of flowers`~ week2$Treatment)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  week2$`Number of flowers` by week2$Treatment
## W = 1752.5, p-value = 1.764e-07
## alternative hypothesis: true location shift is not equal to 0

# Week 3 test 
week3<- data %>% filter (Week==3)
wilcox.test (week3$`Number of flowers`~ week3$Treatment)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  week3$`Number of flowers` by week3$Treatment
## W = 1619.5, p-value = 0.0001725
## alternative hypothesis: true location shift is not equal to 0

# Week 4 test
week4<- data %>% filter (Week==4)
wilcox.test (week4$`Number of flowers`~ week4$Treatment)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  week4$`Number of flowers` by week4$Treatment
## W = 635, p-value = 0.001493
## alternative hypothesis: true location shift is not equal to 0

#more accuratley a mixed model repeated measures can take into account differences between weeks 
lmeModel = lmer(`Number of flowers`  ~ Treatment * Week + (1|City), data=data)
anova(lmeModel)

## Analysis of Variance Table
##                npar  Sum Sq Mean Sq F value
## Treatment         1 12927.4 12927.4 23.8455
## Week              1   417.6   417.6  0.7704
## Treatment:Week    1    57.7    57.7  0.1064

#most significant differences between Treatments than between weeks given the error in cities. 

Is percent coverage of flowering plants greater in No Mow Plots?

cover<-ggplot (data, aes(x=Treatment, y=data$`Percent Cover Angiosperms`, 
                         color=City, fill= City)) + geom_boxplot (alpha=0.5) + 
  theme_bw() +ylab("% Cover of Angiosperms") +facet_grid(City~Week) +   
  scale_color_brewer(palette = "Accent") + scale_fill_brewer(palette = "Accent") 

#cover

# GLM 
#cover_glm<- glm (`Percent Cover Angiosperms`~ Treatment + City + Week,  data=data)
#summary(cover_glm)
#stepAIC(cover_glm)

#repeated measures by week considering error within city 
anovaModelRM = aov(`Percent Cover Angiosperms` ~ Treatment * Week + Error(City), data = data)
summary(anovaModelRM)

## 
## Error: City
##           Df Sum Sq Mean Sq
## Treatment  1  16476   16476
## Week       1   2736    2736
## 
## Error: Within
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## Treatment        1  47949   47949   61.09 2.95e-14 ***
## Week             1  14019   14019   17.86 2.80e-05 ***
## Treatment:Week   1   1547    1547    1.97    0.161    
## Residuals      531 416811     785                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Mixed Model- best estimator for nonlinear trend
lmeModel = lmer(`Percent Cover Angiosperms`  ~ Treatment*Week + (1|City), data=data)
anova(lmeModel)

## Analysis of Variance Table
##                npar Sum Sq Mean Sq F value
## Treatment         1  49626   49626 63.2165
## Week              1  13705   13705 17.4582
## Treatment:Week    1   1558    1558  1.9851

ggplotly (cover)

## Warning: Use of `` data$`Percent Cover Angiosperms` `` is discouraged.
## ℹ Use `Percent Cover Angiosperms` instead.

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

#Pooled abundance differences in two major pollinator groups.

#abundance (count) of insects by group Hymenoptera, Diptera, Other bar plot 
#pol.sub<- pollinators %>% dplyr::select (Treatment, City, Hymenoptera.Abund, Diptera.Abund) %>%gather (Treatment, key=City, value="Abundance", 3:4)
pol.long<-read.csv("abund21.csv")
ggplot (pol.long, aes (x=Treatment, y=Abundance, fill=Taxon)) + 
  geom_bar(stat = "identity", alpha=0.75) + 
  facet_grid(.~City) + theme_bw() + scale_fill_brewer(palette = "Set2")

## Warning: Removed 1 rows containing missing values (`position_stack()`).

pol.long2<-read.csv("rich21.csv")
ggplot (pol.long2, aes (x=Treatment, y=Richness, fill=Taxon)) + 
  geom_boxplot(alpha=0.75) +
  facet_grid(.~City) + theme_bw() + scale_fill_brewer(palette = "Set2")

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

comparison of pollinator abundance vs treatment

h.abund<- ggplot (pollinators, aes (x=Treatment, y=Hymenoptera.Abund, fill=City)) + geom_boxplot(alpha=0.5) + 
    theme_bw() +ylab("Hymenoptera Abundance (Number of Individuals) ") + theme (legend.position = "none") +
      geom_point(aes(x=Treatment, y=Hymenoptera.Abund, color=City),  size=.25) + 
      facet_grid(.~City) + geom_line(aes (group=Comp_site), size=0.25)  + scale_fill_brewer(palette = "Accent") 

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#h.abund
ggplotly(h.abund)

t.test (Hymenoptera.Abund ~ Treatment, paired=T, data=pollinators) #7X more hymenoptera in No Mow Yards 

## 
##  Paired t-test
## 
## data:  Hymenoptera.Abund by Treatment
## t = -7.408, df = 79, p-value = 1.228e-10
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -8.80153 -5.07347
## sample estimates:
## mean difference 
##         -6.9375

# Mixed model 
lmeModel.abund = lmer(Hymenoptera.Abund  ~ Treatment + (1|Comp_site) + (1|City), data=pollinators)

## boundary (singular) fit: see help('isSingular')

anova(lmeModel.abund)

## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## Treatment    1 1925.2  1925.2  54.989

lmeModel.abund

## Linear mixed model fit by REML ['lmerMod']
## Formula: Hymenoptera.Abund ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## REML criterion at convergence: 1018.939
## Random effects:
##  Groups    Name        Std.Dev.
##  Comp_site (Intercept) 0.000   
##  City      (Intercept) 0.000   
##  Residual              5.917   
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## Fixed Effects:
##    (Intercept)  TreatmentNoMow  
##          0.650           6.938  
## optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

summary (lmeModel.abund)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Hymenoptera.Abund ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## 
## REML criterion at convergence: 1018.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2823 -0.3528 -0.1099  0.0697  4.1259 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Comp_site (Intercept)  0.00    0.000   
##  City      (Intercept)  0.00    0.000   
##  Residual              35.01    5.917   
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)      0.6500     0.6615   0.983
## TreatmentNoMow   6.9375     0.9355   7.415
## 
## Correlation of Fixed Effects:
##             (Intr)
## TreatmntNMw -0.707
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

dip.abund<- ggplot (pollinators, aes (x=Treatment, y=Diptera.Abund, fill=City)) + geom_boxplot(alpha=0.5) + 
    theme_bw() +ylab("Diptera Abundance (Number of Individuals) ") + theme (legend.position = "none") +
      geom_point(aes(x=Treatment, y=Diptera.Abund, color=City),  size=.25) + 
      facet_grid(.~City) + geom_line(aes (group=Comp_site), size=0.25)  + scale_fill_brewer(palette = "Accent")

#dip.abund
ggplotly(dip.abund) 

# Mixed model 
lmeModel.abund2 = lmer(Diptera.Abund  ~ Treatment + (1|Comp_site) + (1|City), data=pollinators)
anova(lmeModel.abund2)

## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## Treatment    1 126.03  126.03  40.312

lmeModel.abund2

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diptera.Abund ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## REML criterion at convergence: 643.2084
## Random effects:
##  Groups    Name        Std.Dev.
##  Comp_site (Intercept) 0.2836  
##  City      (Intercept) 0.3649  
##  Residual              1.7681  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## Fixed Effects:
##    (Intercept)  TreatmentNoMow  
##         0.2309          1.7750

summary (lmeModel.abund2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diptera.Abund ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## 
## REML criterion at convergence: 643.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1326 -0.3857 -0.1288  0.0665  5.4743 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Comp_site (Intercept) 0.08044  0.2836  
##  City      (Intercept) 0.13312  0.3649  
##  Residual              3.12626  1.7681  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)      0.2309     0.2998   0.770
## TreatmentNoMow   1.7750     0.2796   6.349
## 
## Correlation of Fixed Effects:
##             (Intr)
## TreatmntNMw -0.466

Abund<-grid.arrange(h.abund,dip.abund )

ggsave(filename = "abundance.png", plot=Abund, width = 20, height = 20, units = "cm")

Comparison of Hymenoptera & Diptera Richness

h.rich<- ggplot (pollinators, aes (x=Treatment, y=Hymenoptera.Rich, fill=City)) + 
  geom_boxplot(alpha=0.5) + theme_bw() + ylab("Hymenoptera Richness (Number  of Species) ") + 
    xlab("") + theme (legend.position = "none") +
    geom_point(aes(x=Treatment, y=Hymenoptera.Rich, color=City),  size=.25) + 
    facet_grid(.~City) + geom_line(aes (group=Comp_site), size=0.25)  +
    scale_fill_brewer(palette = "Accent") 
#h.rich

ggplotly(h.rich)

t.test (Hymenoptera.Rich ~ Treatment, paired=T, data=pollinators) #80 sites, paired observations of richness. Richness is 1.8X higher in No Mow plots 

## 
##  Paired t-test
## 
## data:  Hymenoptera.Rich by Treatment
## t = -7.722, df = 79, p-value = 3.028e-11
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -2.279696 -1.345304
## sample estimates:
## mean difference 
##         -1.8125

hy.summary<- pollinators %>% group_by(Treatment) %>% summarize (mean.rich=  mean (Hymenoptera.Rich ))
hy.summary

## # A tibble: 2 × 2
##   Treatment mean.rich
##   <chr>         <dbl>
## 1 Mow           0.175
## 2 NoMow         1.99

#Consider a mixed model with city and site as random factors 
lmeModel = lmer(Hymenoptera.Rich  ~ Treatment + (1|Comp_site) + (1|City), data=pollinators)
lmeModel

## Linear mixed model fit by REML ['lmerMod']
## Formula: Hymenoptera.Rich ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## REML criterion at convergence: 588.0341
## Random effects:
##  Groups    Name        Std.Dev.
##  Comp_site (Intercept) 0.2507  
##  City      (Intercept) 0.2649  
##  Residual              1.4845  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## Fixed Effects:
##    (Intercept)  TreatmentNoMow  
##         0.2413          1.8125

anova(lmeModel) #highly sig. effect of treatment on hymenoptera richness

## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## Treatment    1 131.41  131.41  59.629

summary (lmeModel)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Hymenoptera.Rich ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## 
## REML criterion at convergence: 588
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5012 -0.4807 -0.0744  0.0639  5.0715 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Comp_site (Intercept) 0.06284  0.2507  
##  City      (Intercept) 0.07019  0.2649  
##  Residual              2.20375  1.4845  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)      0.2413     0.2346   1.029
## TreatmentNoMow   1.8125     0.2347   7.722
## 
## Correlation of Fixed Effects:
##             (Intr)
## TreatmntNMw -0.500

# Now Diptera
d.rich<- ggplot (pollinators, aes (x=Treatment, y=Diptera.Rich, fill=City)) + geom_boxplot(alpha=0.5) + 
    theme_bw() +ylab("Diptera Richness (Number of Species) ") + theme (legend.position = "none") +
      geom_point(aes(x=Treatment, y=Diptera.Rich, color=City),  size=.25) + 
      facet_grid(.~City) + geom_line(aes (group=Comp_site), size=0.25) + scale_fill_brewer(palette = "Accent") 

#d.rich

ggplotly(d.rich)

t.test (Diptera.Rich ~ Treatment, paired=T, data=pollinators) #80 sites, paired observations of richness. Richness is 1.8X higher in No Mow plots 

## 
##  Paired t-test
## 
## data:  Diptera.Rich by Treatment
## t = -5.8625, df = 79, p-value = 1.007e-07
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.8204591 -0.4045409
## sample estimates:
## mean difference 
##         -0.6125

#Consider a mixed model with city and site as random factors 
lmeModel.dip = lmer(Diptera.Rich  ~ Treatment + (1|Comp_site) + (1|City), data=pollinators)
lmeModel.dip

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diptera.Rich ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## REML criterion at convergence: 336.54
## Random effects:
##  Groups    Name        Std.Dev.
##  Comp_site (Intercept) 0.1437  
##  City      (Intercept) 0.2120  
##  Residual              0.6608  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## Fixed Effects:
##    (Intercept)  TreatmentNoMow  
##         0.1368          0.6125

anova(lmeModel.dip) #highly sig. effect of treatment on diptera richness

## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## Treatment    1 15.006  15.006  34.368

summary (lmeModel.dip)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diptera.Rich ~ Treatment + (1 | Comp_site) + (1 | City)
##    Data: pollinators
## 
## REML criterion at convergence: 336.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0541 -0.8725 -0.0618  0.0611  4.8031 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Comp_site (Intercept) 0.02066  0.1437  
##  City      (Intercept) 0.04494  0.2120  
##  Residual              0.43663  0.6608  
## Number of obs: 160, groups:  Comp_site, 80; City, 3
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)      0.1368     0.1476   0.927
## TreatmentNoMow   0.6125     0.1045   5.862
## 
## Correlation of Fixed Effects:
##             (Intr)
## TreatmntNMw -0.354

dip.summary<- pollinators %>% group_by(Treatment) %>% summarize (mean.rich=  mean (Diptera.Rich ))
dip.summary

## # A tibble: 2 × 2
##   Treatment mean.rich
##   <chr>         <dbl>
## 1 Mow           0.075
## 2 NoMow         0.688

summary (pollinators)

##       Site          Comp_site      Treatment             City          
##  Min.   :  2.00   Min.   : 2.00   Length:160         Length:160        
##  1st Qu.: 43.75   1st Qu.:21.75   Class :character   Class :character  
##  Median :150.50   Median :43.50   Mode  :character   Mode  :character  
##  Mean   :147.18   Mean   :47.17                                        
##  3rd Qu.:243.25   3rd Qu.:73.25                                        
##  Max.   :299.00   Max.   :99.00                                        
##  Hymenoptera.Abund Hymenoptera.Rich Diptera.Abund  Diptera.Rich   
##  Min.   : 0.000    Min.   : 0.000   Min.   : 0    Min.   :0.0000  
##  1st Qu.: 0.000    1st Qu.: 0.000   1st Qu.: 0    1st Qu.:0.0000  
##  Median : 0.000    Median : 0.000   Median : 0    Median :0.0000  
##  Mean   : 4.119    Mean   : 1.081   Mean   : 1    Mean   :0.3812  
##  3rd Qu.: 6.000    3rd Qu.: 2.000   3rd Qu.: 1    3rd Qu.:1.0000  
##  Max.   :32.000    Max.   :10.000   Max.   :12    Max.   :4.0000  
##  Other.Insect.Abund Other.Insect.Rich  Total_Abund       Total_Rich    
##  Min.   : 0.0000    Min.   :0.0000    Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 0.0000    1st Qu.:0.0000    1st Qu.: 0.000   1st Qu.: 0.000  
##  Median : 0.0000    Median :0.0000    Median : 0.000   Median : 0.000  
##  Mean   : 0.7688    Mean   :0.2062    Mean   : 5.888   Mean   : 1.669  
##  3rd Qu.: 0.0000    3rd Qu.:0.0000    3rd Qu.: 9.000   3rd Qu.: 2.250  
##  Max.   :15.0000    Max.   :4.0000    Max.   :48.000   Max.   :11.000

# Rich<-grid.arrange(h.rich, d.rich); Rich
# ggsave(filename = "richness.png", plot=Rich, width = 20, height = 20, units = "cm")