---
title: "obs_data_final"
author: "Sonja Glasser"
date: "2026-02-18"
output: html_document
---

We assessed B. impatiens visitation to Asteraceae flowers by observing the 3-4 most abundant cut flower and border flower species during the 2nd and 3rd visit to each farm. Observations occurred during floral surveys or the following day. We conducted three rounds of 15-minute observations per cut flower species. For each observation, we haphazardly selected a 1 m row section, counted the number of flowers and tallied number of bumble bee visits and pollen-collecting behavior (yes/no). Floral species varied between sites resulting in uneven observation time per species. Therefore, we were unable to analyze bumble bee preference but could verify if bumble bees were collecting pollen from floral species. 

So here we summarize the data to show the total number of visits / flowers observed during 15 min period. Then take the average of the per floral species. We will also look at what proportion of those visits resulting in pollen foraging. 

```{r setup, include=FALSE}
require(here)#telling R where to look for a file
require(lme4) # generalized mixed models 
library(performance) # check model performance
library(tidyverse) # language & ggplot
library(car) # Anova
library(emmeans) #extract the estimated marginal means from the model
library(effects)#allEffects()
library(MuMIn)#dredge () model selection
library(AICcmodavg) #compare models 
library(GGally) # pairwise plots 
library(lattice) # another graphing package
library(MASS)#glm.nb
library("glmmTMB")
library("bbmle") 
library(influence.ME)#check for outlier in GLMM
library(DHARMa) #model fit
library(ggpubr)#ggarrange
library(pscl)
library(emmeans)
```

## Read in data 

read in raw data and fix variables to correct format

```{r}
data <- read.csv(here("raw_data", "flwr_obs_submit.csv"), header=TRUE) # originally named flwr_obs_raw_2 

data$visit <- factor(data$visit, 
                         levels = c("1", "2", "3"))
#separate julian.date column
data$date<-as.Date(data$date, "%m/%d/%y") # read as date, but makes them 2020 instead of 2021.

# other variables are just normal factors
cols <- c("sp", "genus", "farm", "pollen_col", "gender")
data<-data %>%
  mutate_each_(funs(factor(.)),cols)

# make a column of yes o no bee visits. 

data$bee_visit<-ifelse(data$sp==0,"N","Y")
```


```{r}
#capitalize genus
data <- data %>% 
  mutate(
    genus_cap = str_to_title(genus)   # Capitalize genus
  ) 

#group flowers based on "border vs. cut"
data <- data %>% 
  mutate(floral_group = ifelse(genus == "eutrochium" | genus == "silphium", "Border Flowers", "Cut Flowers"))

# Reorder panels
data$floral_group <- factor(data$floral_group,
                               levels = c("Cut Flowers", "Border Flowers"))

```



Make data frame with just B. impatiens and summarizing the total number of visits then finding the rate of visit per flower. 
```{r}
visit_imp<-data%>%
  filter(sp == c("imp","0"))%>%
  group_by(genus_cap, farm, n_inflor, bee_visit, floral_group)%>%
  count(bee_visit)%>%
  filter(bee_visit!="NA")%>%
  mutate(visit_per_flwr=n/n_inflor)


```

final plot of visitation rate. 

```{r}

v <- ggplot(visit_imp, aes(x = genus_cap, y = visit_per_flwr)) +
  
  # Raw data FIRST (behind box)
  geom_jitter(width = 0.10,
              size = 1.6,
              alpha = 0.5,
              color = "#0C7BDC") +
  
  # Boxplot on top
  geom_boxplot(width = 0.6,
               outlier.shape = NA,
               linewidth = 1,
               fill = "white",
               color = "black") +
  
  facet_wrap(~ floral_group, nrow = 1, scales = "free_x", space = "free_x") +
  
  theme_classic(base_size = 14) +
  labs(
    x = "Asteraceae genus",
    y = "Visits per flower (15 min)"
  ) +
  theme(
    axis.text.x = element_text(face = "italic", angle = 45, hjust = 1, vjust = 1),
    axis.text.y = element_text(color = "black"),
    axis.title = element_text(face = "bold"),
    axis.line = element_line(linewidth = 0.8),
    
    # Facet polish
    strip.background = element_rect(fill = "white", color = "black", linewidth = 0.8),
    strip.text = element_text(face = "bold")
  )

v

ggsave("visit_rate_boxplot.tiff",
       v,
       width = 8,
       height = 4,
       units = "in",
       dpi = 600,
       compression = "lzw")

```


```{r}
data <- data %>% 
  mutate(pollen_col = ifelse(pollen_col == "Y", 1, 0))

vis_ave<-data%>%
  filter(sp == c("imp","0"))%>%
  group_by(genus_cap, farm, n_inflor, bee_visit, floral_group)%>%
  count(bee_visit)%>%
  filter(bee_visit!="NA")%>%
  mutate(visit_per_flwr=n/n_inflor)%>%
  group_by(genus_cap, floral_group) %>%
  summarise(
    mean_visit = mean(visit_per_flwr, na.rm = TRUE),
    sd_visit   = sd(visit_per_flwr, na.rm = TRUE),
    n             = n(),
    se_visit  = sd_visit / sqrt(n),
 )

pol_ave<-data%>%
  filter(sp == c("imp","0"), pollen_col != "NA")%>%
  group_by(genus_cap, farm, n_inflor, bee_visit, pollen_col, floral_group)%>%
  count(bee_visit)%>%
  filter(bee_visit!="NA")%>%
  mutate(visit_per_flwr=n/n_inflor)%>%
  group_by(genus_cap, floral_group) %>%
  summarise(
    mean_pol = mean(pollen_col, na.rm = TRUE),
    sd_pol   = sd(pollen_col, na.rm = TRUE),
    n             = n(),
    se_pol  = sd_pol / sqrt(n)
  )

visit_ave_pol_ave <- data%>%
  filter(sp == c("imp","0"), pollen_col != "NA")%>%
  group_by(genus_cap, farm, n_inflor, bee_visit, pollen_col, floral_group)%>%
  count(bee_visit)%>%
  filter(bee_visit!="NA")%>%
  mutate(visit_per_flwr=n/n_inflor)%>%
  group_by(genus_cap, floral_group) %>%
  summarise(
    mean_visit = mean(visit_per_flwr, na.rm = TRUE),
    sd_visit   = sd(visit_per_flwr, na.rm = TRUE),
    n             = n(),
    se_visit  = sd_visit / sqrt(n),
     mean_pol = mean(pollen_col, na.rm = TRUE),
    sd_pol   = sd(pollen_col, na.rm = TRUE),
    n             = n(),
    se_pol  = sd_pol / sqrt(n)
 )
```
```{r}
p1 <- ggplot(vis_ave, aes(x = genus_cap, y = mean_visit)) +
  
  geom_point(size = 3) +
  
  geom_errorbar(aes(ymin = mean_visit - se_visit,
                    ymax = mean_visit + se_visit),
                width = 0.15,
                linewidth = 0.7) +
  
  facet_wrap(~ floral_group, nrow = 1, scales = "free_x", space = "free_x") +
  
  labs(
    x = NULL,
    y = "Visits per flower (15 min)"
  ) +
  
  theme_classic(base_size = 14) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    strip.background = element_rect(fill = "white", color = "black"),
    strip.text = element_text(face = "bold")
  )

p2 <- ggplot(pol_ave, aes(x = genus_cap, y = mean_pol)) +
  
  geom_point(size = 3) +
  
  geom_errorbar(aes(ymin = mean_pol - se_pol,
                    ymax = mean_pol + se_pol),
                width = 0.15,
                linewidth = 0.7) +
  
  facet_wrap(~ floral_group, nrow = 1, scales = "free_x", space = "free_x") +
  
  scale_y_continuous(limits = c(0, 1)) +
  
  labs(
    x = "Flower genus",
    y = "Pollen collection proportion"
  ) +
  
  theme_classic(base_size = 14) +
  theme(
    axis.text.x = element_text(face = "italic", angle = 45, hjust = 1),
    strip.background = element_blank(),
    strip.text = element_blank()
  )

library(patchwork)

combined_plot <- p1 / p2 +
  plot_layout(heights = c(1, 1))

combined_plot

```


```{r}
library(viridis)
 visit_ave_pol_ave_fig <- ggplot(visit_ave_pol_ave, aes(x = genus_cap, y = mean_visit)) +
  
  geom_point(aes(fill = mean_pol, size = n),
             shape = 21,
             color = "black", 
             alpha = 0.7) +
  
  geom_errorbar(aes(ymin = mean_visit - se_visit,
                    ymax = mean_visit + se_visit),
                width = 0.15,
                linewidth = 0.7) +
  
  facet_wrap(~ floral_group, nrow = 1, scales = "free_x", space = "free_x") +
  
  scale_fill_gradientn(colours = viridis(256, option = "D"),
                      name = "Pollen collection\nproportion") +
  
  scale_size(range = c(4, 7), guide = "none") +
  
  theme_classic(base_size = 14) +
  theme(
    axis.text.x = element_text(face = "italic", angle = 45, hjust = 1),
    strip.background = element_rect(fill = "white", color = "black"),
    strip.text = element_text(face = "bold")
  ) +
  
  labs(
    x = "Flower genus",
    y = "Visits per flower (15 min)"
  )  

ggsave("visit_ave_pol_ave_fig.tiff",
       visit_ave_pol_ave_fig,
       width = 8,
       height = 4,
       units = "in",
       dpi = 600,
       compression = "lzw")
```

This is the figure with raw visits, not rate of visits. 
```{r}
pol_sum <- pol_imp %>% 
  mutate(
    genus_cap = str_to_title(genus)   # Capitalize genus
  ) %>%
  group_by(genus_cap, floral_group, pollen_col) %>%
  summarise(n_visits = n(), .groups = "drop") %>%
  group_by(genus_cap, floral_group) %>%
  mutate(total_visits = sum(n_visits)) %>%
  ungroup()

pol_sum <- pol_sum %>% 
  group_by(floral_group) %>%
  arrange(desc(total_visits)) %>%
  mutate(genus_cap = factor(genus_cap, levels = unique(genus_cap))) %>%
  ungroup()

# Reorder panels
pol_sum$floral_group <- factor(pol_sum$floral_group,
                               levels = c("Cut Flowers", "Border Flowers"))

# Recode pollen_col → No/Yes
pol_sum <- pol_sum %>% 
  mutate(pollen_cat = ifelse(pollen_col == 1, "Yes", "No"))

# Updated colors
cb_cols <- c("No" = "#E69F00",   # orange
             "Yes" = "#56B4E9") # blue

# Final plot
flwr_obs_pol_col<- ggplot(pol_sum,
       aes(x = genus_cap,
           y = n_visits,
           fill = pollen_cat)) +
  
  geom_bar(stat = "identity") +

  facet_wrap(~ floral_group, nrow = 1, scales = "free_x") +
  
  labs(
    x = "Asteraceae Genus",
    y = "Number of Visits",
    fill = "Pollen Collected"
  ) +
  
  scale_fill_manual(values = cb_cols) +
  
  theme_classic(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "italic", size = 14.5),
    strip.text = element_text(size = 16, face = "bold"),
    panel.spacing.x = unit(1, "lines"),
    legend.position = "right"
  ) 

flwr_obs_pol_col

ggsave("flwr_obs_pol_col.png", flwr_obs_pol_col, width = 7.5, height = 4.5, dpi = 150, units = "in", device='png')


```