library (ggplot2)
library(readxl)
library(forcats)
library (glmmTMB)
library (dplyr)
library (DHARMa)
library (multcomp)
library(emmeans)
library (eoffice)
library(gridExtra)
library(car)
data <- read_excel("municipal_pollinator_data_merged_bumblebees.xlsx",
sheet = "global_database_goodnames")
ants <- subset (data, insect_group == "Ants")
flies <- subset (data, insect_group == "Other flies")
honeybees <- subset (data, insect_group == "Honey bees")
wasps <- subset (data, insect_group == "Wasps")
wildbees <- subset (data, insect_group == "Wild bees")
bees_and_wildbees <- rbind(honeybees, wildbees)
bees_wildbees_flies_wasps <- rbind(honeybees, wildbees, flies, wasps)
First some summary statistics: mean visitation rate for each cultivar
mean_insectsM2 <- data %>%
group_by(cultivar) %>%
summarise(mean_insects_m2 = mean(insects_m2_full_coverage, na.rm = TRUE))
mean_insectsM2
## # A tibble: 35 × 2
## cultivar mean_insects_m2
## <chr> <dbl>
## 1 "Achillea filipendulina" 0.160
## 2 "Alcaltheae spp." 0.529
## 3 "Anemone japonica" 0.914
## 4 "Aster ageratoides" 1.61
## 5 "Begonia semperflorens \"Eureka Scarlet\"" 0.377
## 6 "Begonia semperflorens \"Eureka White\"" 0.318
## 7 "Begonia semperflorens \"red\"" 0.735
## 8 "Begonia semperflorens \"white\"" 0.334
## 9 "Bistorta affinis" 0.857
## 10 "Campanula persicifolia" 0.828
## # ℹ 25 more rows
now some descriptive plots
# mean bar plot
all_insects_by_cultivar_plot <- ggplot(data, aes(x = reorder(cultivar, insects_m2_full_coverage), y = insects_m2_full_coverage)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(ylim = c(0, 6)) +
stat_summary(fun.y = "mean", geom = "bar", size = 3) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 0.3),
axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5)) + # Set angle to 0 for the Y-axis text
labs(x = "Insects per m2", y = "cultivar")
# Rotate the plot by 90 degrees
all_insects_by_cultivar_plot <- all_insects_by_cultivar_plot + coord_flip()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
plot (all_insects_by_cultivar_plot)
#topptx (all_insects_by_cultivar_plot, "fig2_all_insects_by_cultivar_plot.pptx", width = 6, height = 8) #this function is used to export the plot to powerpoint for further editiing and then including in the manuscript.Commented out for the markdown file
Now individual plots of cultivar attractiveness forthe top 4 most common insect groups (honeybees wildbees flies ants wasps)
insectplot_facet_common <- ggplot (bees_wildbees_flies_wasps, aes(x = reorder(cultivar,insects_m2_full_coverage), y = insects_m2_full_coverage))+
facet_wrap(~ insect_group, nrow = 4, ncol = 1)+
scale_y_continuous(expand = c(0, 0)) + # forces X axis to 0, but in this case is overriden by ribbon
coord_cartesian(ylim = c(0, 26)) +
stat_summary(fun.y = "mean", geom = "bar", size = 3, fill = "red") +
ylab("Average insects per M2") +
xlab("Cultivar") +
theme_bw(15) +
theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust = 0.3))
plot (insectplot_facet_common)
First some summary statistics
data_summary <- data %>%
group_by(insect_group) %>%
summarise(total_count = sum(count)) %>%
mutate(percentage = (total_count / sum(total_count)) * 100) %>%
arrange(desc(percentage))
data_summary
## # A tibble: 10 × 3
## insect_group total_count percentage
## <chr> <dbl> <dbl>
## 1 Honey bees 305 36.2
## 2 Wild bees 189 22.4
## 3 Other flies 124 14.7
## 4 Ants 90 10.7
## 5 Wasps 56 6.64
## 6 Bumblebees 30 3.56
## 7 Hover flies 22 2.61
## 8 True bugs 15 1.78
## 9 Beetles 6 0.712
## 10 Butterflies 6 0.712
first just one bar chart showing contribution of all insect groups
# Calculate the total count for each insect group
total_counts <- data %>%
group_by(insect_group) %>%
summarise(total_count = sum(count)) %>%
arrange(desc(total_count)) # Order by total count in descending order
# Create a bar chart
insect_contribution_plot <- ggplot(total_counts, aes(x = reorder(insect_group, -total_count), y = total_count)) +
geom_bar(stat = "identity") +
labs( x = "Insect Group",
y = "Total Count") +
theme_minimal()
plot (insect_contribution_plot)
#topptx (insect_contribution_plot, "insect_contribution_plot.pptx", width = 6, height = 6)
This helps us see who the visitors are to each cultivar
# Calculate the percentage contribution of 'insects_m2_full_coverage' for each 'insect_group' within each 'cultivar'
data_percent <- data %>%
group_by(cultivar, insect_group) %>%
summarise(total_insects_m2 = sum(insects_m2_full_coverage),
original_values = mean(insects_m2_full_coverage)) %>%
group_by(cultivar) %>%
mutate(percentage = total_insects_m2 / sum(total_insects_m2) * 100)
## `summarise()` has grouped output by 'cultivar'. You can override using the
## `.groups` argument.
# Create a stacked bar chart showing the percentage contribution
filledstackedplot_cultivar <- ggplot(data_percent, aes(x = reorder(cultivar, original_values), y = percentage, fill = insect_group)) +
geom_bar(stat = "identity") +
labs(y = "Cultivar", x = "Percentage Contribution (%)") + # Swap x and y axis labels
scale_fill_brewer(palette = "Set3") +
theme_minimal() +
theme(axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.3)) +
coord_flip()
# Display the plot
plot(filledstackedplot_cultivar)
#topptx (filledstackedplot_cultivar, "fig3_filledstackedplot_cultivar.pptx", width = 6, height = 8)
Various publications claim that there is no relationship between planted area size and insect visitation rates per unit area. I will quickly confirm visually.
ggplot (data, aes(x = total_area_cultivars, y = insects_m2_full_coverage)) +
geom_jitter(width = 10, height = 1, alpha = 0.3) +
geom_smooth(method = "lm") +
theme_bw(15) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
labs( x = "total area of planted plants in location", y = "Insects per m2")
## `geom_smooth()` using formula = 'y ~ x'
Confirmed. No apparent effect. Naturally, this is also heavily skewed by
the zero inflation, but certainly nothing obvious here. A quick
statistical analysis to be found in the stats part.
To reduce dimensionality and persue clear questions, here we are focussing in on just wild bees and honeybees. Is there an overlap between wild bees and honeybees? What plants do wild bees like but honeybees don’t?
fig3_honey_wildbees <- ggplot (bees_and_wildbees, aes(x = reorder(cultivar,insects_m2_full_coverage), y = insects_m2_full_coverage, fill = insect_group))+
scale_y_continuous(expand = c(0, 0)) + # forces X axis to 0, but in this case is overriden by ribbon
coord_cartesian(ylim = c(0, 20)) +
stat_summary(fun.y = "mean", geom = "bar", shape = 23, size = 3, alpha = 0.5) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.3)) + # hjust and vjust nudge the x-axis label positions
labs( x = "cultivar", y = "Insect visits per m2") +
coord_flip()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
plot (fig3_honey_wildbees)
#topptx (fig3_honey_wildbees, "fig3_wildvshoney.pptx", width = 6, height = 6)
some descriptive states here:
sum of count
# Calculate the sum of "count" for each combination of "cultivar" and "insect_group"
bees_and_wildbees %>%
group_by(cultivar, insect_group) %>%
summarise(total_count = sum(count))
## `summarise()` has grouped output by 'cultivar'. You can override using the
## `.groups` argument.
## # A tibble: 70 × 3
## # Groups: cultivar [35]
## cultivar insect_group total_count
## <chr> <chr> <dbl>
## 1 "Achillea filipendulina" Honey bees 0
## 2 "Achillea filipendulina" Wild bees 0
## 3 "Alcaltheae spp." Honey bees 2
## 4 "Alcaltheae spp." Wild bees 1
## 5 "Anemone japonica" Honey bees 27
## 6 "Anemone japonica" Wild bees 12
## 7 "Aster ageratoides" Honey bees 7
## 8 "Aster ageratoides" Wild bees 13
## 9 "Begonia semperflorens \"Eureka Scarlet\"" Honey bees 4
## 10 "Begonia semperflorens \"Eureka Scarlet\"" Wild bees 0
## # ℹ 60 more rows
bees_and_wildbees %>%
group_by(insect_group) %>%
summarise(total_count = sum(count))
## # A tibble: 2 × 2
## insect_group total_count
## <chr> <dbl>
## 1 Honey bees 305
## 2 Wild bees 189
Questions we want to ask: - do cultivars differ in attractiveness overall? - do cultivears differ in attractiveness for specific target groups? - does total patch area affect insect visitation rates - then do pairwise comparisons - which cultivars has more wild bees than honeybees?
m1 <- glmmTMB(insects_m2_full_coverage ~ cultivar
+ (1 | location),
ziformula=~1,
family=nbinom1,
data = data)
looking at results, testing model fit
# testing model fit
simulateResiduals(m1, n = 500, plot = T) #dharma
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
## Object of Class DHARMa with simulated residuals based on 500 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.7908908 0.966 0.988 0.4696725 0.3231412 0.5229729 0.06865622 0.1879041 0.8175542 0.8212179 0.8565376 0.4296221 0.964 0.1516091 0.01738105 0.6381312 0.2055175 0.4215931 0.4491866 0.2059352 ...
# top level summary
Anova (m1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: insects_m2_full_coverage
## Chisq Df Pr(>Chisq)
## cultivar 416.29 34 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#adjust for multiple comparisons
adjX <- glht(m1) # multcomp
summary(adjX, test = adjusted("BH"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmmTMB(formula = insects_m2_full_coverage ~ cultivar + (1 |
## location), data = data, family = nbinom1, ziformula = ~1,
## dispformula = ~1)
##
## Linear Hypotheses:
## Estimate Std. Error z value
## (Intercept) == 0 1.69469 0.60692 2.792
## cultivarAlcaltheae spp. == 0 0.14181 0.64412 0.220
## cultivarAnemone japonica == 0 0.59901 0.61460 0.975
## cultivarAster ageratoides == 0 0.40294 0.62316 0.647
## cultivarBegonia semperflorens "Eureka Scarlet" == 0 0.53668 0.65169 0.824
## cultivarBegonia semperflorens "Eureka White" == 0 -0.04651 0.62876 -0.074
## cultivarBegonia semperflorens "red" == 0 0.23150 0.72878 0.318
## cultivarBegonia semperflorens "white" == 0 0.25822 0.68265 0.378
## cultivarBistorta affinis == 0 0.37178 0.62214 0.598
## cultivarCampanula persicifolia == 0 0.84611 0.66669 1.269
## cultivarCanna indica == 0 0.16360 0.64022 0.256
## cultivarCleome spinosa == 0 -0.01248 0.62178 -0.020
## cultivarEchinacea purpurea "Alba" == 0 0.45830 0.62413 0.734
## cultivarEchinacea purpurea "pink" == 0 0.50941 0.61761 0.825
## cultivarErigeron x cultorum == 0 0.75820 0.63069 1.202
## cultivarErynginum planum "Blauer Zwerg" == 0 2.40088 0.61694 3.892
## cultivarEryngium planum "Flachblatt-Mannstreu" == 0 1.59038 0.64522 2.465
## cultivarGalinsoga parviflora == 0 0.26557 0.73107 0.363
## cultivarGeranium x magnificum == 0 0.40775 0.60995 0.669
## cultivarHibiscus syriacus == 0 0.32815 0.63489 0.517
## cultivarHosta x cultorum "Funkie" == 0 0.82405 0.64045 1.287
## cultivarHydrangea macrophylla == 0 -3.18657 1.22587 -2.599
## cultivarNepeta grandiflora == 0 0.22392 0.62556 0.358
## cultivarPetunia cultivars "Nightsky" == 0 -3.23855 0.97183 -3.332
## cultivarPetunia cultivars "Queen of Hearts" == 0 0.46119 0.66103 0.698
## cultivarPhlox amplifolia == 0 0.56841 0.69933 0.813
## cultivarRudbeckia fulgida == 0 0.39484 0.60951 0.648
## cultivarSalvia farianacea == 0 0.31851 0.61009 0.522
## cultivarSalvia splendens "purple" == 0 0.70736 0.62108 1.139
## cultivarSalvia splendens "Red" == 0 0.45030 0.61740 0.729
## cultivarSalvia splendens "white" == 0 0.33364 0.65074 0.513
## cultivarSalvia yangii == 0 1.00350 0.61436 1.633
## cultivarSedum telephium == 0 0.43070 0.61380 0.702
## cultivarTagetes patula "orange" == 0 0.53862 0.63136 0.853
## cultivarTagetes patula "yellow" == 0 0.46713 0.62104 0.752
## Pr(>|z|)
## (Intercept) == 0 0.06106 .
## cultivarAlcaltheae spp. == 0 0.87580
## cultivarAnemone japonica == 0 0.78809
## cultivarAster ageratoides == 0 0.78809
## cultivarBegonia semperflorens "Eureka Scarlet" == 0 0.78809
## cultivarBegonia semperflorens "Eureka White" == 0 0.96872
## cultivarBegonia semperflorens "red" == 0 0.84762
## cultivarBegonia semperflorens "white" == 0 0.84044
## cultivarBistorta affinis == 0 0.78835
## cultivarCampanula persicifolia == 0 0.78809
## cultivarCanna indica == 0 0.87314
## cultivarCleome spinosa == 0 0.98398
## cultivarEchinacea purpurea "Alba" == 0 0.78809
## cultivarEchinacea purpurea "pink" == 0 0.78809
## cultivarErigeron x cultorum == 0 0.78809
## cultivarErynginum planum "Blauer Zwerg" == 0 0.00349 **
## cultivarEryngium planum "Flachblatt-Mannstreu" == 0 0.09595 .
## cultivarGalinsoga parviflora == 0 0.84044
## cultivarGeranium x magnificum == 0 0.78809
## cultivarHibiscus syriacus == 0 0.78835
## cultivarHosta x cultorum "Funkie" == 0 0.78809
## cultivarHydrangea macrophylla == 0 0.08171 .
## cultivarNepeta grandiflora == 0 0.84044
## cultivarPetunia cultivars "Nightsky" == 0 0.01507 *
## cultivarPetunia cultivars "Queen of Hearts" == 0 0.78809
## cultivarPhlox amplifolia == 0 0.78809
## cultivarRudbeckia fulgida == 0 0.78809
## cultivarSalvia farianacea == 0 0.78835
## cultivarSalvia splendens "purple" == 0 0.78809
## cultivarSalvia splendens "Red" == 0 0.78809
## cultivarSalvia splendens "white" == 0 0.78835
## cultivarSalvia yangii == 0 0.59723
## cultivarSedum telephium == 0 0.78809
## cultivarTagetes patula "orange" == 0 0.78809
## cultivarTagetes patula "yellow" == 0 0.78809
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- BH method)
#emmeans
meanie <- emmeans(m1, pairwise ~ cultivar)
# print (meanie) # no point plotting all pairwise comparisons for all pairs - we'd have a heck of a lot of pairs. Just look at the plot.
plot (meanie)
Going to use the raw count, because we’re going to compare only within cultivar. Means we can use poisson distribution too
Going to only look at top 6 most attractive cultivars:
unique(data$cultivar)
# List of cultivars to filter
selected_cultivars <- c("Salvia yangii", "Aster ageratoides", "Erynginum planum \"Blauer Zwerg\"", "Echinacea purpurea \"Alba\"", "Campanula persicifolia", "Anemone japonica")
# Creating a subset dataset containing only the specified cultivars
top6cultivars <- bees_and_wildbees[bees_and_wildbees$cultivar %in% selected_cultivars, ]
make a figure, check all is in order
ggplot (top6cultivars , aes(x = cultivar, y = count, fill = insect_group))+
scale_y_continuous(expand = c(0, 0)) + # forces X axis to 0, but in this case is overriden by ribbon
coord_cartesian(ylim = c(0, 45)) +
stat_summary(fun.y = "sum", geom = "bar", shape = 23, size = 3, alpha = 0.5) +
theme_bw(12) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.3)) + # hjust and vjust nudge the x-axis label positions +
labs( x = "cultivar", y = "Total insects counted", title ="Sum count honey and wildbees, top 6 cutlivars")
## Warning in stat_summary(fun.y = "sum", geom = "bar", shape = 23, size = 3, :
## Ignoring unknown parameters: `shape`
ok now lets see if we get do sensible modelling and pairwise comparisons with this
##NOTE IM USING THE RAW COUNT HERE, NOT THE M2 VARIABLE.
#model
m2beestop6 <- glmmTMB(count ~ as.factor(cultivar) * insect_group
+ (1 | location),
ziformula=~1,
family=poisson,
data = top6cultivars)
Anova (m2beestop6)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: count
## Chisq Df Pr(>Chisq)
## as.factor(cultivar) 9.0092 5 0.108698
## insect_group 6.8307 1 0.008960 **
## as.factor(cultivar):insect_group 19.5148 5 0.001541 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Dharma
simulateResiduals(m2beestop6, n = 500, plot = T) #dharma
## Object of Class DHARMa with simulated residuals based on 500 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.7831745 0.9960342 0.3899917 0.8681415 0.3261663 0.8148226 0.3691393 0.4257348 0.949657 0.1081472 0.4299255 0.954141 0.9653428 0.9459946 0.7635606 0.8870011 0.2395449 0.5948557 0.04640933 0.8098306 ...
adjX <- glht(m2beestop6) # multcomp
summary(adjX, test = adjusted("BH"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmmTMB(formula = count ~ as.factor(cultivar) * insect_group +
## (1 | location), data = top6cultivars, family = poisson, ziformula = ~1,
## dispformula = ~1)
##
## Linear Hypotheses:
## Estimate
## (Intercept) == 0 -9.271e-02
## as.factor(cultivar)Aster ageratoides == 0 8.551e-03
## as.factor(cultivar)Campanula persicifolia == 0 -1.264e+00
## as.factor(cultivar)Echinacea purpurea "Alba" == 0 5.316e-01
## as.factor(cultivar)Erynginum planum "Blauer Zwerg" == 0 -1.957e+01
## as.factor(cultivar)Salvia yangii == 0 9.413e-01
## insect_groupWild bees == 0 -7.861e-01
## as.factor(cultivar)Aster ageratoides:insect_groupWild bees == 0 1.389e+00
## as.factor(cultivar)Campanula persicifolia:insect_groupWild bees == 0 1.451e+00
## as.factor(cultivar)Echinacea purpurea "Alba":insect_groupWild bees == 0 -1.266e+00
## as.factor(cultivar)Erynginum planum "Blauer Zwerg":insect_groupWild bees == 0 1.992e+01
## as.factor(cultivar)Salvia yangii:insect_groupWild bees == 0 -2.897e+00
## Std. Error
## (Intercept) == 0 2.922e-01
## as.factor(cultivar)Aster ageratoides == 0 5.530e-01
## as.factor(cultivar)Campanula persicifolia == 0 1.038e+00
## as.factor(cultivar)Echinacea purpurea "Alba" == 0 3.520e-01
## as.factor(cultivar)Erynginum planum "Blauer Zwerg" == 0 8.003e+03
## as.factor(cultivar)Salvia yangii == 0 3.252e-01
## insect_groupWild bees == 0 3.624e-01
## as.factor(cultivar)Aster ageratoides:insect_groupWild bees == 0 6.188e-01
## as.factor(cultivar)Campanula persicifolia:insect_groupWild bees == 0 1.300e+00
## as.factor(cultivar)Echinacea purpurea "Alba":insect_groupWild bees == 0 8.461e-01
## as.factor(cultivar)Erynginum planum "Blauer Zwerg":insect_groupWild bees == 0 8.003e+03
## as.factor(cultivar)Salvia yangii:insect_groupWild bees == 0 1.080e+00
## z value
## (Intercept) == 0 -0.317
## as.factor(cultivar)Aster ageratoides == 0 0.015
## as.factor(cultivar)Campanula persicifolia == 0 -1.217
## as.factor(cultivar)Echinacea purpurea "Alba" == 0 1.510
## as.factor(cultivar)Erynginum planum "Blauer Zwerg" == 0 -0.002
## as.factor(cultivar)Salvia yangii == 0 2.894
## insect_groupWild bees == 0 -2.169
## as.factor(cultivar)Aster ageratoides:insect_groupWild bees == 0 2.244
## as.factor(cultivar)Campanula persicifolia:insect_groupWild bees == 0 1.116
## as.factor(cultivar)Echinacea purpurea "Alba":insect_groupWild bees == 0 -1.496
## as.factor(cultivar)Erynginum planum "Blauer Zwerg":insect_groupWild bees == 0 0.002
## as.factor(cultivar)Salvia yangii:insect_groupWild bees == 0 -2.682
## Pr(>|z|)
## (Intercept) == 0 0.9980
## as.factor(cultivar)Aster ageratoides == 0 0.9980
## as.factor(cultivar)Campanula persicifolia == 0 0.3831
## as.factor(cultivar)Echinacea purpurea "Alba" == 0 0.2694
## as.factor(cultivar)Erynginum planum "Blauer Zwerg" == 0 0.9980
## as.factor(cultivar)Salvia yangii == 0 0.0439
## insect_groupWild bees == 0 0.0902
## as.factor(cultivar)Aster ageratoides:insect_groupWild bees == 0 0.0902
## as.factor(cultivar)Campanula persicifolia:insect_groupWild bees == 0 0.3964
## as.factor(cultivar)Echinacea purpurea "Alba":insect_groupWild bees == 0 0.2694
## as.factor(cultivar)Erynginum planum "Blauer Zwerg":insect_groupWild bees == 0 0.9980
## as.factor(cultivar)Salvia yangii:insect_groupWild bees == 0 0.0439
##
## (Intercept) == 0
## as.factor(cultivar)Aster ageratoides == 0
## as.factor(cultivar)Campanula persicifolia == 0
## as.factor(cultivar)Echinacea purpurea "Alba" == 0
## as.factor(cultivar)Erynginum planum "Blauer Zwerg" == 0
## as.factor(cultivar)Salvia yangii == 0 *
## insect_groupWild bees == 0 .
## as.factor(cultivar)Aster ageratoides:insect_groupWild bees == 0 .
## as.factor(cultivar)Campanula persicifolia:insect_groupWild bees == 0
## as.factor(cultivar)Echinacea purpurea "Alba":insect_groupWild bees == 0
## as.factor(cultivar)Erynginum planum "Blauer Zwerg":insect_groupWild bees == 0
## as.factor(cultivar)Salvia yangii:insect_groupWild bees == 0 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- BH method)
#pairwise
meanie <- emmeans(m2beestop6, pairwise ~ insect_group|cultivar)
summary (meanie)
## $emmeans
## cultivar = Anemone japonica:
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees -0.0927 0.292 Inf -6.65e-01 0.480
## Wild bees -0.8788 0.371 Inf -1.61e+00 -0.152
##
## cultivar = Aster ageratoides:
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees -0.0842 0.551 Inf -1.16e+00 0.995
## Wild bees 0.5184 0.462 Inf -3.87e-01 1.424
##
## cultivar = Campanula persicifolia:
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees -1.3566 1.045 Inf -3.40e+00 0.692
## Wild bees -0.6912 0.760 Inf -2.18e+00 0.799
##
## cultivar = Echinacea purpurea "Alba":
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees 0.4389 0.342 Inf -2.32e-01 1.109
## Wild bees -1.6128 0.758 Inf -3.10e+00 -0.128
##
## cultivar = Erynginum planum "Blauer Zwerg":
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees -19.6665 8003.425 Inf -1.57e+04 15666.758
## Wild bees -0.5284 0.699 Inf -1.90e+00 0.841
##
## cultivar = Salvia yangii:
## insect_group emmean SE df asymp.LCL asymp.UCL
## Honey bees 0.8486 0.269 Inf 3.22e-01 1.375
## Wild bees -2.8349 1.031 Inf -4.86e+00 -0.813
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## cultivar = Anemone japonica:
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees 0.786 0.362 Inf 2.169 0.0301
##
## cultivar = Aster ageratoides:
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees -0.603 0.499 Inf -1.207 0.2276
##
## cultivar = Campanula persicifolia:
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees -0.665 1.248 Inf -0.533 0.5940
##
## cultivar = Echinacea purpurea "Alba":
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees 2.052 0.766 Inf 2.679 0.0074
##
## cultivar = Erynginum planum "Blauer Zwerg":
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees -19.138 8003.425 Inf -0.002 0.9981
##
## cultivar = Salvia yangii:
## contrast estimate SE df z.ratio p.value
## Honey bees - Wild bees 3.684 1.020 Inf 3.612 0.0003
##
## Results are given on the log (not the response) scale.