For now, only females are included in all analyses
Show the code
## find honeybee counts <- df_sites %>% group_by (park_name) %>% summarise (total_hb = sum (total_honey_bees, na.rm = TRUE ))# simple categories for flowers <- df_flower_raw %>% mutate (floral_type = case_when (== "Native" ~ "native" ,== "Non-Native" & cultivar_weedy %in% c ("Weedy" , "Weedy/Cultivar" ) ~ "weedy" ,== "Non-Native" & cultivar_weedy %in% c ("Cultivar" ) ~ "cultivar" ,TRUE ~ "cultivar" # find floral abundance by species at each site <- df_flower_clean %>% mutate (total_abundance = sum (transect_1_quantity + transect_2_quantity + transect_3_quantity, na.rm = TRUE )) %>% group_by (park_name, latin_name) %>% summarise (sum_abundance = sum (total_abundance), plant_type = first (floral_type))# summarize floral data across all visits and losing plant identity <- df_flower_by_species %>% group_by (park_name) %>% summarise (species_richness = n (), total_abundance = log (sum (sum_abundance)), count_native = sum (plant_type == "native" ), count_weedy = sum (plant_type == "weedy" ), count_cultivar = sum (plant_type == "cultivar" ), prop_native = count_native/ species_richness, prop_weedy = count_weedy/ species_richness, prop_cultivar = count_cultivar/ species_richness)# join nancy's bee measurements with full bee data so it can be joined to various metrics <- df_bombus %>% :: select (id, ave_of_measurements) %>% mutate (genus_species = "Bombus fervidus" ,mean_itd = ave_of_measurements) %>% :: select (id, genus_species, mean_itd)<- df_halictus %>% :: select (id, ave_of_measurements) %>% mutate (genus_species = "Halictus ligatus" ,mean_itd = ave_of_measurements) %>% :: select (id, genus_species, mean_itd)<- df_megachile %>% :: select (id, ave_of_measurements) %>% mutate (genus_species = "Megachile rotundata" ,mean_itd = ave_of_measurements) %>% :: select (id, genus_species, mean_itd)<- bind_rows (df_flt_bombus, df_flt_halictus, df_flt_megachile)<- inner_join (df_full_bees, df_body_sizes, by = c ("bee_id" = "id" ))<- full_join (df_joined_bees, df_hb_count, by = c ("park_name" )) %>% full_join (., df_soceco, by = c ("park_name" )) %>% full_join (., df_flower_abun_rich, by = c ("park_name" )) %>% filter (sex == "F" )
Plot of body size as a function of honeybee abundance
Show the code
%>% ggplot (., aes (x = total_hb, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.3 ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 23 ) + labs (x = "Honeybee Abundance" , y = "Body Size (ITD)" , color = "Bee Species" ) + scale_color_manual (values = c ("goldenrod" , "#F7879A" , "#74A089" )) + theme (legend.position = "none" )
Show the code
<- df_all_joined %>% ungroup () %>% group_by (genus_species) %>% nest () %>% mutate (fit = purrr:: map (data, ~ lm (mean_itd ~ total_hb, data = .x)),tidied = purrr:: map (fit, tidy),glanced = purrr:: map (fit, glance)):: unnest (mod_hb, tidied) %>% :: select_if (sapply (., class) != "list" ) %>% mutate (across (where (is.numeric), round, 4 ))
# A tibble: 6 × 6
# Groups: genus_species [3]
genus_species term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Bombus fervidus (Intercept) 3.26 0.0257 127. 0
2 Bombus fervidus total_hb 0.0001 0 1.94 0.0538
3 Halictus ligatus (Intercept) 1.57 0.0115 137. 0
4 Halictus ligatus total_hb 0 0 -1.70 0.0911
5 Megachile rotundata (Intercept) 2.16 0.0156 139. 0
6 Megachile rotundata total_hb 0 0 -0.706 0.481
Plot of body size as a function of median household income within 1km
Show the code
%>% ggplot (., aes (x = x1km_householdincome_medhinc_cy, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.5 ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 22 ) + labs (x = "Median Household \n Income (1 km)" , y = "Body Size (ITD)" , color = "Bee Species" ) + scale_color_manual (values = c ("goldenrod" , "#F7879A" , "#74A089" )) + theme (legend.position = "none" , axis.text.x = element_text (size = 15 ))
Show the code
<- df_all_joined %>% ungroup () %>% group_by (genus_species) %>% nest () %>% mutate (fit = purrr:: map (data, ~ lm (mean_itd ~ x1km_householdincome_medhinc_cy, data = .x)),tidied = purrr:: map (fit, tidy),glanced = purrr:: map (fit, glance)):: unnest (mod_income, tidied) %>% :: select_if (sapply (., class) != "list" ) %>% mutate (across (where (is.numeric), round, 4 ))
# A tibble: 6 × 6
# Groups: genus_species [3]
genus_species term estim…¹ std.e…² stati…³ p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Bombus fervidus (Intercept) 3.33 0.0412 80.9 0
2 Bombus fervidus x1km_householdincome_medh… 0 0 -0.704 0.482
3 Halictus ligatus (Intercept) 1.56 0.0235 66.3 0
4 Halictus ligatus x1km_householdincome_medh… 0 0 -0.0725 0.942
5 Megachile rotundata (Intercept) 2.15 0.0259 82.8 0
6 Megachile rotundata x1km_householdincome_medh… 0 0 0.211 0.833
# … with abbreviated variable names ¹estimate, ²std.error, ³statistic
Plot of body size as a function of proportion of park that is native plants
Show the code
%>% ggplot (., aes (x = prop_native, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.5 ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 22 ) + labs (x = "Proportion Native Plants" , y = "Body Size (ITD)" , color = "Bee Species" ) + scale_color_manual (values = c ("goldenrod" , "#F7879A" , "#74A089" )) + theme (legend.position = "none" )
Show the code
%>% filter (genus_species == "Halictus ligatus" ) %>% ggplot (., aes (x = prop_native, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.5 , color = "#F7879A" ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 22 ) + labs (x = "Proportion Native Plants" , y = "Body Size (ITD)" , color = "Bee Species" ) + # scale_color_manual(values = c("goldenrod", "#F7879A", "#74A089")) + theme (legend.position = "none" )
Show the code
<- df_all_joined %>% ungroup () %>% group_by (genus_species) %>% nest () %>% mutate (fit = purrr:: map (data, ~ lm (mean_itd ~ prop_native, data = .x)),tidied = purrr:: map (fit, tidy),glanced = purrr:: map (fit, glance)):: unnest (mod_native, tidied) %>% :: select_if (sapply (., class) != "list" ) %>% mutate (across (where (is.numeric), round, 7 ))
# A tibble: 6 × 6
# Groups: genus_species [3]
genus_species term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Bombus fervidus (Intercept) 3.28 0.0422 77.7 0
2 Bombus fervidus prop_native 0.0810 0.133 0.609 0.543
3 Halictus ligatus (Intercept) 1.44 0.0281 51.2 0
4 Halictus ligatus prop_native 0.303 0.0700 4.32 0.0000215
5 Megachile rotundata (Intercept) 2.14 0.0328 65.1 0
6 Megachile rotundata prop_native 0.0519 0.0969 0.535 0.593
Plot of body size as a function of park floral abundance
Show the code
%>% ggplot (., aes (x = total_abundance, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.5 ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 21.8 ) + labs (x = "Floral Abundance (log10)" , y = "Body Size (ITD)" , color = "Bee Species" ) + scale_color_manual (values = c ("goldenrod" , "#F7879A" , "#74A089" )) + theme (legend.position = "none" )
Show the code
<- df_all_joined %>% ungroup () %>% group_by (genus_species) %>% nest () %>% mutate (fit = purrr:: map (data, ~ lm (mean_itd ~ total_abundance, data = .x)),tidied = purrr:: map (fit, tidy),glanced = purrr:: map (fit, glance)):: unnest (mod_pabun, tidied) %>% :: select_if (sapply (., class) != "list" ) %>% mutate (across (where (is.numeric), round, 4 ))
# A tibble: 6 × 6
# Groups: genus_species [3]
genus_species term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Bombus fervidus (Intercept) 3.05 0.481 6.35 0
2 Bombus fervidus total_abundance 0.0154 0.0299 0.515 0.607
3 Halictus ligatus (Intercept) 2.55 0.302 8.44 0
4 Halictus ligatus total_abundance -0.0622 0.0189 -3.30 0.0011
5 Megachile rotundata (Intercept) 2.08 0.278 7.45 0
6 Megachile rotundata total_abundance 0.0048 0.0174 0.276 0.783
Plot of body size as a function of park floral richness
Show the code
%>% ggplot (., aes (x = species_richness, y = mean_itd, color = genus_species)) + geom_jitter (size = 3 , alpha = 0.5 ) + geom_smooth (method = "lm" ) + theme_classic (base_size = 23 ) + labs (x = "Floral Richness" , y = "Body Size (ITD)" , color = "Bee Species" ) + scale_color_manual (values = c ("goldenrod" , "#F7879A" , "#74A089" )) + theme (legend.position = "none" )
Show the code
<- df_all_joined %>% ungroup () %>% group_by (genus_species) %>% nest () %>% mutate (fit = purrr:: map (data, ~ lm (mean_itd ~ species_richness, data = .x)),tidied = purrr:: map (fit, tidy),glanced = purrr:: map (fit, glance)):: unnest (mod_prich, tidied) %>% :: select_if (sapply (., class) != "list" ) %>% mutate (across (where (is.numeric), round, 4 ))
# A tibble: 6 × 6
# Groups: genus_species [3]
genus_species term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Bombus fervidus (Intercept) 3.28 0.0316 104. 0
2 Bombus fervidus species_richness 0.0005 0.0006 0.827 0.409
3 Halictus ligatus (Intercept) 1.60 0.0187 85.5 0
4 Halictus ligatus species_richness -0.0009 0.0004 -2.55 0.0113
5 Megachile rotundata (Intercept) 2.16 0.0183 117. 0
6 Megachile rotundata species_richness -0.0001 0.0003 -0.205 0.838
Generally, sites “add” floral species richness by having more horticultural cultivars - rather than increased numbers of native plants
Show the code
%>% ggplot (., aes (x = prop_cultivar, y = species_richness)) + geom_point ()
Not actually predicted by abundance of plants used by Halictus, though! Perhaps suggest something about the management style of the park more generally? (i.e. more suitable nesting habitat at more “native” parks?)
Show the code
<- df_all_joined %>% filter (genus == "Halictus" , species == "ligatus" ) %>% distinct (floral_association_genus, floral_association_species) %>% mutate (latin_name = paste (floral_association_genus, floral_association_species)) %>% pull (latin_name)%>% filter (genus == "Halictus" , species == "ligatus" ) %>% mutate (latin_name = paste (floral_association_genus, floral_association_species)) %>% group_by (latin_name) %>% tally () %>% arrange (desc (n))
# A tibble: 45 × 2
latin_name n
<chr> <int>
1 Ratibida columnifera 55
2 Rudbeckia hirta 40
3 Erigeron divergens 36
4 Achillea millefolium 33
5 Heterotheca villosa 15
6 Helianthus annuus 12
7 Achillea tomentosa 8
8 Erigeron speciosus 8
9 Coreopsis verticillata 6
10 Taraxacum officinale 6
# … with 35 more rows
Show the code
<- df_flower_by_species %>% filter (latin_name %in% v_halictus_flws) %>% group_by (park_name) %>% summarise (halictus_fl_abun = log (sum (sum_abundance)))%>% filter (genus_species == "Halictus ligatus" ) %>% full_join (., df_halictus_flowers, by = c ("park_name" )) %>% ggplot (., aes (x = halictus_fl_abun, y = mean_itd)) + geom_jitter (size = 3 , alpha = 0.5 , color = "#F7879A" ) + geom_smooth (color = "#F7879A" , method = "lm" ) + theme_classic (base_size = 23 ) + labs (x = "H. ligatus host flower \n abundance (log10)" , y = "Body Size (ITD)" ) + theme (axis.title.x = element_text (size = 20 ))
Show the code
%>% filter (genus_species == "Halictus ligatus" ) %>% full_join (., df_halictus_flowers, by = c ("park_name" )) %>% lm (data = ., mean_itd ~ halictus_fl_abun) %>% summary ()
Call:
lm(formula = mean_itd ~ halictus_fl_abun, data = .)
Residuals:
Min 1Q Median 3Q Max
-0.44625 -0.09496 -0.01051 0.08349 0.39583
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.28880 0.45441 5.037 8.33e-07 ***
halictus_fl_abun -0.04820 0.02988 -1.613 0.108
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1458 on 290 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.008896, Adjusted R-squared: 0.005479
F-statistic: 2.603 on 1 and 290 DF, p-value: 0.1077
