-------------------------------------------------------------------------------------
1. Differences between temperature treatments (Table S7)

##repeated measurement anova
library(tidyr)
library(lme4)
library(emmeans)

#all data
## prepare temperature ##
temp22 <- read.table("temp_2022.txt", sep = "\t", header = T) 

temp22 <- subset(temp22, Day_Night == "Day") #subset sunrise to sunset

temp22$DOY <- date.to.DOY(temp22$Date, format = "dd/mm/yyyy") #convert date to day of year

temp22$trtmt_arena <- paste(temp22$Treatment, temp22$Arena, sep="_") #concatenate treatment and arena

#summarize average temperature by DOY, trtmt and arena
temp22_avg <- temp22 %>%
  group_by(DOY, trtmt_arena) %>%
  summarise(temp_avg = mean(Temp))

unique(temp22_avg$DOY) #check


temp <- separate(temp22_avg, "trtmt_arena", into = c("trtmt", "arena"), sep = "_") 
m1 <- lmer(temp_avg ~ trtmt + (1|DOY), data = temp)
summary(m1)

emmeans(m1, pairwise ~ trtmt) #pairwise comparisons


--------------------------------------------------------------------------------------
2. Volumetric water content (Table S5, Figure S2)

library(dplyr)
library(TDPanalysis)

## prepare temperature ##
temp22_avg <- read.table("temp2022_avg.txt", sep = "\t", header = T) 


## prepare water content ##
VWC <- read.table("VWC_soil.txt", sep = "\t", header = T)

VWC$DOY <- date.to.DOY(VWC$date, format = "mm/dd/yyyy") #convert date to day of year

VWC$arena <- "TR" #add column with arena ID

VWC$trtmt_arena <- paste(VWC$trtmt, VWC$arena, sep="_") #concatenate treatment and arena

unique(VWC$DOY)

## join temp and VWC per arena ##
temp_VWC <- merge(temp22_avg, VWC)
unique(temp_VWC$DOY)

temp_VWC$prop_VWC <- temp_VWC$VWC/100
temp_VWC$species <- as.factor(temp_VWC$species)


library(glmmTMB)
VW_tot <- glmmTMB(prop_VWC ~ temp_avg + (1|DOY), family = beta_family(link = "logit"), data = temp_VWC) #general model
summary(VW_tot)

VW_CH <- glmmTMB(prop_VWC ~ temp_avg + (1|DOY), family = beta_family(link = "logit"), data = subset(temp_VWC, species == "CH")) #model for Collinsia
summary(VW_CH)

VW_NM <- glmmTMB(prop_VWC ~ temp_avg + (1|DOY), family = beta_family(link = "logit"), data = subset(temp_VWC, species == "NM")) #model for Nemophila
summary(VW_NM)

VW_PC <- glmmTMB(prop_VWC ~ temp_avg + (1|DOY), family = beta_family(link = "logit"), data = subset(temp_VWC, species == "PC")) #model for Phacelia
summary(VW_PC)


#plot
library(ggplot2)
library(ggeffects)
pr3 <- predict_response(VW_tot, terms = "temp_avg", type = "fixed", bias_correction = TRUE)

plt3 <- plot(pr3, show_data = TRUE)


xx <- plt3$plot_env$rawdat$x #use raw data for colored points
yy <- (plt3$plot_env$rawdat$response)*100
predx <- plt3$plot_env$plot_data$x #use predicted data for line
predy <- (plt3$plot_env$plot_data$predicted)*100
CIhigh <- (plt3$plot_env$plot_data$conf.high)*100
CIlow <- (plt3$plot_env$plot_data$conf.low)*100
species <- c("CH", "NM", "PC")

A <- temp_VWC[,c(3,7,9)]
B <- data.frame(temp_avg = predx, VWC = predy, species, CIlow = CIlow, CIhigh = CIhigh)



ggplot(A) +
  geom_point(data = A, aes(x = temp_avg, y = VWC, color = species, shape = species), size = 3, alpha = 0.4) +
  geom_line(data = B, aes(x = temp_avg, y = VWC), color='black', linewidth = 1) +
  geom_ribbon(data = B, aes(x = temp_avg, ymin = CIlow, ymax = CIhigh), alpha = 0.2) +
  #geom_line(data = B, aes(x = temp_avg, y = CIlow), color='darkgrey', linewidth = 1, linetype = 2) +
  #geom_line(data = B, aes(x = temp_avg, y = CIhigh), color='darkgrey', linewidth = 1, linetype = 2) +
  scale_color_manual(values = c("#B12A90FF", "#31688EFF", "orange")) +
  labs(title = "", x = "Temperature (°C)", y = "Volumetric water content (%)") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title = element_text(size=38)) +
  theme(legend.title = element_text(size = 18), legend.text = element_text(size = 14))


------------------------------------------------------------------------------------------------------------
3. Phenology analyses

--------------------------------------------------------
3a. Phenological distributions (Table S1, Figures 2, S1)
library(dplyr)
library(tidyr)
library(moments)
library(ggplot2)
library(lme4)
library(lmerTest)

##prepare data
#flowers
flwr_count <- read.table("flwrcount_temp2022.txt", sep = "\t", header = T)

flwr_count <- separate(flwr_count, "trtmt_arena", into = c("trtmt","arena"), sep = "_") #separate trtmt and arena

flwr_sum <- subset(flwr_count, arena != "TR") %>% 
  group_by(DOY, trtmt, arena) %>% #grouped species to look at plant community effect
  summarise(flowers = sum(flwr_count),
            tmpavg = mean(temp_avg))

#bees
bee_count <- read.table("bee_abund_2022.txt", sep = "\t", header = T)

bee_count$arena <- as.character(bee_count$arena)

flwr_avg2 <- subset(flwr_count, arena != "TR") %>% #extract average temperature per arena
  group_by(DOY, trtmt, arena) %>% 
  summarise(tmpavg = mean(temp_avg))

bees <- left_join(bee_count, flwr_avg2)

bee_sum <- bees %>% 
  group_by(DOY, trtmt, arena) %>% 
  summarise(bees = sum(bees),
            tmpavg = mean(tmpavg))


#join flower and bee df, and format
flwr_bees <- flwr_sum %>% left_join(bee_sum)
flwr_bees[is.na(flwr_bees)] <- 0 #replace NAs with 0s for DOYs in bees
flwr_bees <- flwr_bees[, c(1,2,3,5,4,6)] #reorder cols
flwr_bees_long <- flwr_bees %>% pivot_longer(cols = flowers:bees, names_to = "type", values_to = "abund") #wide to long format


##calculate mean temperature for arena within treatment
arena_temp <- flwr_bees_long %>% 
  group_by(trtmt, arena) %>% 
  summarise(meantmp = mean(tmpavg))


##subset data by arena within treatment
#flowers
A1_flwr <- subset(flwr_sum, trtmt == "A" & arena == "1")
A2_flwr <- subset(flwr_sum, trtmt == "A" & arena == "2")
A3_flwr <- subset(flwr_sum, trtmt == "A" & arena == "3")
A4_flwr <- subset(flwr_sum, trtmt == "A" & arena == "4")
A5_flwr <- subset(flwr_sum, trtmt == "A" & arena == "5")

W2_1_flwr <- subset(flwr_sum, trtmt == "W2" & arena == "1")
W2_2_flwr <- subset(flwr_sum, trtmt == "W2" & arena == "2")
W2_3_flwr <- subset(flwr_sum, trtmt == "W2" & arena == "3")
W2_4_flwr <- subset(flwr_sum, trtmt == "W2" & arena == "4")
W2_5_flwr <- subset(flwr_sum, trtmt == "W2" & arena == "5")

W4_1_flwr <- subset(flwr_sum, trtmt == "W4" & arena == "1")
W4_2_flwr <- subset(flwr_sum, trtmt == "W4" & arena == "2")
W4_3_flwr <- subset(flwr_sum, trtmt == "W4" & arena == "3")
W4_4_flwr <- subset(flwr_sum, trtmt == "W4" & arena == "4")
W4_5_flwr <- subset(flwr_sum, trtmt == "W4" & arena == "5")

W6_1_flwr <- subset(flwr_sum, trtmt == "W6" & arena == "1")
W6_2_flwr <- subset(flwr_sum, trtmt == "W6" & arena == "2")
W6_3_flwr <- subset(flwr_sum, trtmt == "W6" & arena == "3")
W6_4_flwr <- subset(flwr_sum, trtmt == "W6" & arena == "4")
W6_5_flwr <- subset(flwr_sum, trtmt == "W6" & arena == "5")

#bees
A1_bee <- subset(bee_sum, trtmt == "A" & arena == "1")
A2_bee <- subset(bee_sum, trtmt == "A" & arena == "2")
A3_bee <- subset(bee_sum, trtmt == "A" & arena == "3")
A4_bee <- subset(bee_sum, trtmt == "A" & arena == "4")
A5_bee <- subset(bee_sum, trtmt == "A" & arena == "5")

W2_1_bee <- subset(bee_sum, trtmt == "W2" & arena == "1")
W2_2_bee <- subset(bee_sum, trtmt == "W2" & arena == "2")
W2_3_bee <- subset(bee_sum, trtmt == "W2" & arena == "3")
W2_4_bee <- subset(bee_sum, trtmt == "W2" & arena == "4")
W2_5_bee <- subset(bee_sum, trtmt == "W2" & arena == "5")

W4_1_bee <- subset(bee_sum, trtmt == "W4" & arena == "1")
W4_2_bee <- subset(bee_sum, trtmt == "W4" & arena == "2")
W4_3_bee <- subset(bee_sum, trtmt == "W4" & arena == "3")
W4_4_bee <- subset(bee_sum, trtmt == "W4" & arena == "4")
W4_5_bee <- subset(bee_sum, trtmt == "W4" & arena == "5")

W6_1_bee <- subset(bee_sum, trtmt == "W6" & arena == "1")
W6_2_bee <- subset(bee_sum, trtmt == "W6" & arena == "2")
W6_3_bee <- subset(bee_sum, trtmt == "W6" & arena == "3")
W6_4_bee <- subset(bee_sum, trtmt == "W6" & arena == "4")
W6_5_bee <- subset(bee_sum, trtmt == "W6" & arena == "5")


###############
## Skeweness ##
###############

## flowers by arena ##
A1_sk <- skewness(A1_flwr$flowers)
A1_agt <- agostino.test(A1_flwr$flowers)
A2_sk <- skewness(A2_flwr$flowers)
A2_agt <- agostino.test(A2_flwr$flowers)
A3_sk <- skewness(A3_flwr$flowers)
A3_agt <- agostino.test(A3_flwr$flowers)
A4_sk <- skewness(A4_flwr$flowers)
A4_agt <- agostino.test(A4_flwr$flowers)
A5_sk <- skewness(A5_flwr$flowers)
A5_agt <- agostino.test(A5_flwr$flowers)

W2_1_sk <- skewness(W2_1_flwr$flowers)
W2_1_agt <- agostino.test(W2_1_flwr$flowers)
W2_2_sk <- skewness(W2_2_flwr$flowers)
W2_2_agt <- agostino.test(W2_2_flwr$flowers)
W2_3_sk <- skewness(W2_3_flwr$flowers)
W2_3_agt <- agostino.test(W2_3_flwr$flowers)
W2_4_sk <- skewness(W2_4_flwr$flowers)
W2_4_agt <- agostino.test(W2_4_flwr$flowers)
W2_5_sk <- skewness(W2_5_flwr$flowers)
W2_5_agt <- agostino.test(W2_5_flwr$flowers)

W4_1_sk <- skewness(W4_1_flwr$flowers)
W4_1_agt <- agostino.test(W4_1_flwr$flowers)
W4_2_sk <- skewness(W4_2_flwr$flowers)
W4_2_agt <- agostino.test(W4_2_flwr$flowers)
W4_3_sk <- skewness(W4_3_flwr$flowers)
W4_3_agt <- agostino.test(W4_3_flwr$flowers)
W4_4_sk <- skewness(W4_4_flwr$flowers)
W4_4_agt <- agostino.test(W4_4_flwr$flowers)
W4_5_sk <- skewness(W4_5_flwr$flowers)
W4_5_agt <- agostino.test(W4_5_flwr$flowers)

W6_1_sk <- skewness(W6_1_flwr$flowers)
W6_1_agt <- agostino.test(W6_1_flwr$flowers)
W6_2_sk <- skewness(W6_2_flwr$flowers)
W6_2_agt <- agostino.test(W6_2_flwr$flowers)
W6_3_sk <- skewness(W6_3_flwr$flowers)
W6_3_agt <- agostino.test(W6_3_flwr$flowers)
W6_4_sk <- skewness(W6_4_flwr$flowers)
W6_4_agt <- agostino.test(W6_4_flwr$flowers)
W6_5_sk <- skewness(W6_5_flwr$flowers)
W6_5_agt <- agostino.test(W6_5_flwr$flowers)

mm <- data.frame(t(data.frame(c(A1_agt[1], A2_agt[1], A3_agt[1], A4_agt[1], A5_agt[1],
                                W2_1_agt[1], W2_2_agt[1], W2_3_agt[1], W2_4_agt[1], W2_5_agt[1],
                                W4_1_agt[1], W4_2_agt[1], W4_3_agt[1], W4_4_agt[1], W4_5_agt[1],
                                W6_1_agt[1], W6_2_agt[1], W6_3_agt[1], W6_4_agt[1], W6_5_agt[1]))))

ww <- data.frame(t(data.frame(c(A1_agt[2], A2_agt[2], A3_agt[2], A4_agt[2], A5_agt[2],
                                W2_1_agt[2], W2_2_agt[2], W2_3_agt[2], W2_4_agt[2], W2_5_agt[2],
                                W4_1_agt[2], W4_2_agt[2], W4_3_agt[2], W4_4_agt[2], W4_5_agt[2],
                                W6_1_agt[2], W6_2_agt[2], W6_3_agt[2], W6_4_agt[2], W6_5_agt[2]))))

dd <- cbind(mm, ww)

colnames(dd) <- c("skewness","z","p-val") #rename columns
rownames(dd) <- NULL #remove rownames

sk_flwr <- cbind(arena_temp, dd)


## bees by arena ##
A1_sk_bee <- skewness(A1_bee$bees)
A1_agt_bee <- agostino.test(A1_bee$bees, alternative = "less")
A2_sk_bee <- skewness(A2_bee$bees)
A2_agt_bee <- agostino.test(A2_bee$bees, alternative = "less")
A3_sk_bee <- skewness(A3_bee$bees)
A3_agt_bee <- agostino.test(A3_bee$bees, alternative = "less")
A4_sk_bee <- skewness(A4_bee$bees)
A4_agt_bee <- agostino.test(A4_bee$bees, alternative = "less")
A5_sk_bee <- skewness(A5_bee$bees)
A5_agt_bee <- agostino.test(A5_bee$bees, alternative = "less")

W2_1_sk_bee <- skewness(W2_1_bee$bees)
W2_1_agt_bee <- agostino.test(W2_1_bee$bees, alternative = "less")
W2_2_sk_bee <- skewness(W2_2_bee$bees)
W2_2_agt_bee <- agostino.test(W2_2_bee$bees, alternative = "less")
W2_3_sk_bee <- skewness(W2_3_bee$bees)
W2_3_agt_bee <- agostino.test(W2_3_bee$bees, alternative = "less")
W2_4_sk_bee <- skewness(W2_4_bee$bees)
W2_4_agt_bee <- agostino.test(W2_4_bee$bees, alternative = "less")
W2_5_sk_bee <- skewness(W2_5_bee$bees)
W2_5_agt_bee <- agostino.test(W2_5_bee$bees, alternative = "less")

W4_1_sk_bee <- skewness(W4_1_bee$bees)
W4_1_agt_bee <- agostino.test(W4_1_bee$bees, alternative = "less")
W4_2_sk_bee <- skewness(W4_2_bee$bees)
W4_2_agt_bee <- agostino.test(W4_2_bee$bees, alternative = "less")
W4_3_sk_bee <- skewness(W4_3_bee$bees)
W4_3_agt_bee <- agostino.test(W4_3_bee$bees, alternative = "less")
W4_4_sk_bee <- skewness(W4_4_bee$bees)
W4_4_agt_bee <- agostino.test(W4_4_bee$bees, alternative = "less")
W4_5_sk_bee <- skewness(W4_5_bee$bees)
W4_5_agt_bee <- agostino.test(W4_5_bee$bees, alternative = "less")

W6_1_sk_bee <- skewness(W6_1_bee$bees)
W6_1_agt_bee <- agostino.test(W6_1_bee$bees, alternative = "less")
W6_2_sk_bee <- skewness(W6_2_bee$bees)
W6_2_agt_bee <- agostino.test(W6_2_bee$bees, alternative = "less")
W6_3_sk_bee <- skewness(W6_3_bee$bees)
W6_3_agt_bee <- agostino.test(W6_3_bee$bees, alternative = "less")
W6_4_sk_bee <- skewness(W6_4_bee$bees)
W6_4_agt_bee <- agostino.test(W6_4_bee$bees, alternative = "less")
W6_5_sk_bee <- skewness(W6_5_bee$bees)
W6_5_agt_bee <- agostino.test(W6_5_bee$bees, alternative = "less")

nn <- data.frame(t(data.frame(c(A1_agt_bee[1], A2_agt_bee[1], A3_agt_bee[1], A4_agt_bee[1], A5_agt_bee[1],
                                W2_1_agt_bee[1], W2_2_agt_bee[1], W2_3_agt_bee[1], W2_4_agt_bee[1], W2_5_agt_bee[1],
                                W4_1_agt_bee[1], W4_2_agt_bee[1], W4_3_agt_bee[1], W4_4_agt_bee[1], W4_5_agt_bee[1],
                                W6_1_agt_bee[1], W6_2_agt_bee[1], W6_3_agt_bee[1], W6_4_agt_bee[1], W6_5_agt_bee[1]))))


qq <- data.frame(t(data.frame(c(A1_agt_bee[2], A2_agt_bee[2], A3_agt_bee[2], A4_agt_bee[2], A5_agt_bee[2],
                                W2_1_agt_bee[2], W2_2_agt_bee[2], W2_3_agt_bee[2], W2_4_agt_bee[2], W2_5_agt_bee[2],
                                W4_1_agt_bee[2], W4_2_agt_bee[2], W4_3_agt_bee[2], W4_4_agt_bee[2], W4_5_agt_bee[2],
                                W6_1_agt_bee[2], W6_2_agt_bee[2], W6_3_agt_bee[2], W6_4_agt_bee[2], W6_5_agt_bee[2]))))

ff <- cbind(nn, qq)

colnames(ff) <- c("skewness","z","p-val") #rename columns
rownames(ff) <- NULL #remove rownames

sk_bees <- cbind(arena_temp, ff)


#bind bee and flower skewness
sk_flwr$type <- "flwr"
sk_bees$type <- "bee"
sk_tot <- rbind(sk_flwr, sk_bees)
sk_tot$type = factor(sk_tot$type, levels = c("flwr","bee"))


###########################################################################################


##############
## Kurtosis ##
##############

## flowers by arena ##
A1_ku <- kurtosis(A1_flwr$flowers)
A1_agl <- anscombe.test(A1_flwr$flowers)
A2_ku <- kurtosis(A2_flwr$flowers)
A2_agl <- anscombe.test(A2_flwr$flowers)
A3_ku <- kurtosis(A3_flwr$flowers)
A3_agl <- anscombe.test(A3_flwr$flowers)
A4_ku <- kurtosis(A4_flwr$flowers)
A4_agl <- anscombe.test(A4_flwr$flowers)
A5_ku <- kurtosis(A5_flwr$flowers)
A5_agl <- anscombe.test(A5_flwr$flowers)

W2_1_ku <- kurtosis(W2_1_flwr$flowers)
W2_1_agl <- anscombe.test(W2_1_flwr$flowers)
W2_2_ku <- kurtosis(W2_2_flwr$flowers)
W2_2_agl <- anscombe.test(W2_2_flwr$flowers)
W2_3_ku <- kurtosis(W2_3_flwr$flowers)
W2_3_agl <- anscombe.test(W2_3_flwr$flowers)
W2_4_ku <- kurtosis(W2_4_flwr$flowers)
W2_4_agl <- anscombe.test(W2_4_flwr$flowers)
W2_5_ku <- kurtosis(W2_5_flwr$flowers)
W2_5_agl <- anscombe.test(W2_5_flwr$flowers)

W4_1_ku <- kurtosis(W4_1_flwr$flowers)
W4_1_agl <- anscombe.test(W4_1_flwr$flowers)
W4_2_ku <- kurtosis(W4_2_flwr$flowers)
W4_2_agl <- anscombe.test(W4_2_flwr$flowers)
W4_3_ku <- kurtosis(W4_3_flwr$flowers)
W4_3_agl <- anscombe.test(W4_3_flwr$flowers)
W4_4_ku <- kurtosis(W4_4_flwr$flowers)
W4_4_agl <- anscombe.test(W4_4_flwr$flowers)
W4_5_ku <- kurtosis(W4_5_flwr$flowers)
W4_5_agl <- anscombe.test(W4_5_flwr$flowers)

W6_1_ku <- kurtosis(W6_1_flwr$flowers)
W6_1_agl <- anscombe.test(W6_1_flwr$flowers)
W6_2_ku <- kurtosis(W6_2_flwr$flowers)
W6_2_agl <- anscombe.test(W6_2_flwr$flowers)
W6_3_ku <- kurtosis(W6_3_flwr$flowers)
W6_3_agl <- anscombe.test(W6_3_flwr$flowers)
W6_4_ku <- kurtosis(W6_4_flwr$flowers)
W6_4_agl <- anscombe.test(W6_4_flwr$flowers)
W6_5_ku <- kurtosis(W6_5_flwr$flowers)
W6_5_agl <- anscombe.test(W6_5_flwr$flowers)


bb <- data.frame(t(data.frame(c(A1_agl[1], A2_agl[1], A3_agl[1], A4_agl[1], A5_agl[1],
                                W2_1_agl[1], W2_2_agl[1], W2_3_agl[1], W2_4_agl[1], W2_5_agl[1],
                                W4_1_agl[1], W4_2_agl[1], W4_3_agl[1], W4_4_agl[1], W4_5_agl[1],
                                W6_1_agl[1], W6_2_agl[1], W6_3_agl[1], W6_4_agl[1], W6_5_agl[1]))))

rr <- data.frame(t(data.frame(c(A1_agl[2], A2_agl[2], A3_agl[2], A4_agl[2], A5_agl[2],
                                W2_1_agl[2], W2_2_agl[2], W2_3_agl[2], W2_4_agl[2], W2_5_agl[2],
                                W4_1_agl[2], W4_2_agl[2], W4_3_agl[2], W4_4_agl[2], W4_5_agl[2],
                                W6_1_agl[2], W6_2_agl[2], W6_3_agl[2], W6_4_agl[2], W6_5_agl[2]))))

ss <- cbind(bb, rr)

colnames(ss) <- c("kurtosis","z","p-val") #rename columns
rownames(ss) <- NULL #remove rownames

ku_flwr <- cbind(arena_temp, ss)


## bees by arena ##
A1_ku_bee <- kurtosis(A1_bee$bees)
A1_agl_bee <- anscombe.test(A1_bee$bees, alternative = "less")
A2_ku_bee <- kurtosis(A2_bee$bees)
A2_agl_bee <- anscombe.test(A2_bee$bees, alternative = "less")
A3_ku_bee <- kurtosis(A3_bee$bees)
A3_agl_bee <- anscombe.test(A3_bee$bees, alternative = "less")
A4_ku_bee <- kurtosis(A4_bee$bees)
A4_agl_bee <- anscombe.test(A4_bee$bees, alternative = "less")
A5_ku_bee <- kurtosis(A5_bee$bees)
A5_agl_bee <- anscombe.test(A5_bee$bees, alternative = "less")

W2_1_ku_bee <- kurtosis(W2_1_bee$bees)
W2_1_agl_bee <- anscombe.test(W2_1_bee$bees, alternative = "less")
W2_2_ku_bee <- kurtosis(W2_2_bee$bees)
W2_2_agl_bee <- anscombe.test(W2_2_bee$bees, alternative = "less")
W2_3_ku_bee <- kurtosis(W2_3_bee$bees)
W2_3_agl_bee <- anscombe.test(W2_3_bee$bees, alternative = "less")
W2_4_ku_bee <- kurtosis(W2_4_bee$bees)
W2_4_agl_bee <- anscombe.test(W2_4_bee$bees, alternative = "less")
W2_5_ku_bee <- kurtosis(W2_5_bee$bees)
W2_5_agl_bee <- anscombe.test(W2_5_bee$bees, alternative = "less")

W4_1_ku_bee <- kurtosis(W4_1_bee$bees)
W4_1_agl_bee <- anscombe.test(W4_1_bee$bees, alternative = "less")
W4_2_ku_bee <- kurtosis(W4_2_bee$bees)
W4_2_agl_bee <- anscombe.test(W4_2_bee$bees, alternative = "less")
W4_3_ku_bee <- kurtosis(W4_3_bee$bees)
W4_3_agl_bee <- anscombe.test(W4_3_bee$bees, alternative = "less")
W4_4_ku_bee <- kurtosis(W4_4_bee$bees)
W4_4_agl_bee <- anscombe.test(W4_4_bee$bees, alternative = "less")
W4_5_ku_bee <- kurtosis(W4_5_bee$bees)
W4_5_agl_bee <- anscombe.test(W4_5_bee$bees, alternative = "less")

W6_1_ku_bee <- kurtosis(W6_1_bee$bees)
W6_1_agl_bee <- anscombe.test(W6_1_bee$bees, alternative = "less")
W6_2_ku_bee <- kurtosis(W6_2_bee$bees)
W6_2_agl_bee <- anscombe.test(W6_2_bee$bees, alternative = "less")
W6_3_ku_bee <- kurtosis(W6_3_bee$bees)
W6_3_agl_bee <- anscombe.test(W6_3_bee$bees, alternative = "less")
W6_4_ku_bee <- kurtosis(W6_4_bee$bees)
W6_4_agl_bee <- anscombe.test(W6_4_bee$bees, alternative = "less")
W6_5_ku_bee <- kurtosis(W6_5_bee$bees)
W6_5_agl_bee <- anscombe.test(W6_5_bee$bees, alternative = "less")


vv <- data.frame(t(data.frame(c(A1_agl_bee[1], A2_agl_bee[1], A3_agl_bee[1], A4_agl_bee[1], A5_agl_bee[1],
                                W2_1_agl_bee[1], W2_2_agl_bee[1], W2_3_agl_bee[1], W2_4_agl_bee[1], W2_5_agl_bee[1],
                                W4_1_agl_bee[1], W4_2_agl_bee[1], W4_3_agl_bee[1], W4_4_agl_bee[1], W4_5_agl_bee[1],
                                W6_1_agl_bee[1], W6_2_agl_bee[1], W6_3_agl_bee[1], W6_4_agl_bee[1], W6_5_agl_bee[1]))))


tt <- data.frame(t(data.frame(c(A1_agl_bee[2], A2_agl_bee[2], A3_agl_bee[2], A4_agl_bee[2], A5_agl_bee[2],
                                W2_1_agl_bee[2], W2_2_agl_bee[2], W2_3_agl_bee[2], W2_4_agl_bee[2], W2_5_agl_bee[2],
                                W4_1_agl_bee[2], W4_2_agl_bee[2], W4_3_agl_bee[2], W4_4_agl_bee[2], W4_5_agl_bee[2],
                                W6_1_agl_bee[2], W6_2_agl_bee[2], W6_3_agl_bee[2], W6_4_agl_bee[2], W6_5_agl_bee[2]))))

gg <- cbind(vv, tt)

colnames(gg) <- c("kurtosis","z","p-val") #rename columns
rownames(gg) <- NULL #remove rownames

ku_bees <- cbind(arena_temp, gg)


#bind bee and flower kurtosis
ku_flwr$type <- "flwr"
ku_bees$type <- "bee"
ku_tot <- rbind(ku_flwr, ku_bees)
ku_tot$type = factor(ku_tot$type, levels = c("flwr","bee"))

###########################################################################################


#######################
## Distribution mean ##
#######################

## flowers by arena ##
A1_mn <- mean(A1_flwr$flowers)
A2_mn <- mean(A2_flwr$flowers)
A3_mn <- mean(A3_flwr$flowers)
A4_mn <- mean(A4_flwr$flowers)
A5_mn <- mean(A5_flwr$flowers)

W2_1_mn <- mean(W2_1_flwr$flowers)
W2_2_mn <- mean(W2_2_flwr$flowers)
W2_3_mn <- mean(W2_3_flwr$flowers)
W2_4_mn <- mean(W2_4_flwr$flowers)
W2_5_mn <- mean(W2_5_flwr$flowers)

W4_1_mn <- mean(W4_1_flwr$flowers)
W4_2_mn <- mean(W4_2_flwr$flowers)
W4_3_mn <- mean(W4_3_flwr$flowers)
W4_4_mn <- mean(W4_4_flwr$flowers)
W4_5_mn <- mean(W4_5_flwr$flowers)

W6_1_mn <- mean(W6_1_flwr$flowers)
W6_2_mn <- mean(W6_2_flwr$flowers)
W6_3_mn <- mean(W6_3_flwr$flowers)
W6_4_mn <- mean(W6_4_flwr$flowers)
W6_5_mn <- mean(W6_5_flwr$flowers)

oo <- data.frame(c(A1_mn, A2_mn, A3_mn, A4_mn, A5_mn,
                   W2_1_mn, W2_2_mn, W2_3_mn, W2_4_mn, W2_5_mn,
                   W4_1_mn, W4_2_mn, W4_3_mn, W4_4_mn, W4_5_mn,
                   W6_1_mn, W6_2_mn, W6_3_mn, W6_4_mn, W6_5_mn))


colnames(oo) <- "mean_distr" #rename columns
rownames(oo) <- NULL #remove rownames

mean_flwr <- cbind(arena_temp, oo)


## bees by arena ##
A1_mn_bees <- mean(A1_bee$bees)
A2_mn_bees <- mean(A2_bee$bees)
A3_mn_bees <- mean(A3_bee$bees)
A4_mn_bees <- mean(A4_bee$bees)
A5_mn_bees <- mean(A5_bee$bees)

W2_1_mn_bees <- mean(W2_1_bee$bees)
W2_2_mn_bees <- mean(W2_2_bee$bees)
W2_3_mn_bees <- mean(W2_3_bee$bees)
W2_4_mn_bees <- mean(W2_4_bee$bees)
W2_5_mn_bees <- mean(W2_5_bee$bees)

W4_1_mn_bees <- mean(W4_1_bee$bees)
W4_2_mn_bees <- mean(W4_2_bee$bees)
W4_3_mn_bees <- mean(W4_3_bee$bees)
W4_4_mn_bees <- mean(W4_4_bee$bees)
W4_5_mn_bees <- mean(W4_5_bee$bees)

W6_1_mn_bees <- mean(W6_1_bee$bees)
W6_2_mn_bees <- mean(W6_2_bee$bees)
W6_3_mn_bees <- mean(W6_3_bee$bees)
W6_4_mn_bees <- mean(W6_4_bee$bees)
W6_5_mn_bees <- mean(W6_5_bee$bees)

pp <- data.frame(c(A1_mn_bees, A2_mn_bees, A3_mn_bees, A4_mn_bees, A5_mn_bees,
                   W2_1_mn_bees, W2_2_mn_bees, W2_3_mn_bees, W2_4_mn_bees, W2_5_mn_bees,
                   W4_1_mn_bees, W4_2_mn_bees, W4_3_mn_bees, W4_4_mn_bees, W4_5_mn_bees,
                   W6_1_mn_bees, W6_2_mn_bees, W6_3_mn_bees, W6_4_mn_bees, W6_5_mn_bees))


colnames(pp) <- "mean_distr" #rename columns
rownames(pp) <- NULL #remove rownames

mean_bees <- cbind(arena_temp, pp)


#bind bee and flower mean
mean_flwr$type <- "flwr"
mean_bees$type <- "bee"
mean_tot <- rbind(mean_flwr, mean_bees)
mean_tot$type = factor(mean_tot$type, levels = c("flwr","bee"))


##############################################################################
######################
# interaction models #
######################

lm_mean <- lm(mean_distr ~ type * meantmp, data = mean_tot) 
summary(lm_mean)

lm_sk <- lm(skewness ~ type * meantmp, data = sk_tot) 
summary(lm_sk)

lm_ku <- lm(kurtosis ~ type * meantmp, data = ku_tot) 
summary(lm_ku)



################################################################################
#########
# plots #
#########
library(ggplot2)
library(patchwork)

##mean
mf <- ggplot(mean_flwr, aes(x = meantmp, y = mean_distr)) +
  geom_point(size = 3, alpha = 0.6, color = "#1b9e77") +
  stat_smooth(method = "lm", color = "#1b9e77", linewidth = 2) +
  #ylim(-0.2,1.6) +
  labs(title = "", x = "", y = "Mean phenology") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title = element_text(size=38), axis.text.x=element_text(vjust=0.5), axis.title.x=element_text(vjust=0.5))
mf

mb <- ggplot(mean_bee, aes(x = meantmp, y = mean_distr)) +
  geom_point(size = 3, alpha = 0.6, color = "#d95f02") +
  stat_smooth(method = "lm", color = "#d95f02", linewidth = 2) +
  #ylim(-0.2,1.6) +
  labs(title = "", x = "Temperature (°C)", y = "Mean phenology") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title = element_text(size=38), axis.text.x=element_text(vjust=0.5), axis.title.x=element_text(vjust=0.5))
mb

mf / mb +
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 28))


##skewness
sk <- ggplot(sk_tot, aes(x = meantmp, y = skewness, color = type)) +
  geom_point(size = 3, alpha = 0.6) +
  stat_smooth(method = "lm", linewidth = 2) +
  scale_color_manual(values = c("#1b9e77", "#d95f02")) +
  #ylim(-0.2,1.6) +
  labs(title = "", x = "Temperature (°C)", y = "Skewness") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title = element_text(size=38), axis.text.x=element_text(vjust=0.5), axis.title.x=element_text(vjust=0.5))
sk


##breadth
ku <- ggplot(ku_tot, aes(x = meantmp, y = kurtosis, color = type)) +
  geom_point(size = 3, alpha = 0.6) +
  stat_smooth(method = "lm", linewidth = 2, linetype = "dashed") +
  scale_color_manual(values = c("#1b9e77", "#d95f02")) +
  #ylim(-0.2,1.6) +
  labs(title = "", x = "", y = "Breadth") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title = element_text(size=38), axis.text.x=element_text(vjust=0.5), axis.title.x=element_text(vjust=0.5))
ku

ku / sk 


-----------------------------------------------
3b. Phenological shifts (Table S2, Figure 2)

library(tidyr)
library(dplyr)
library(glmmTMB)
library(interactions)
library(ggplot2)

#############
## flowers ##
#############
flwr_count <- read.table("flwrcount_temp2022.txt", sep = "\t", header = T)

flwr_count <- separate(flwr_count, "trtmt_arena", into = c("trtmt","arena"), sep = "_") #separate trtmt and arena


flwrs <- subset(flwr_count, arena != "TR") %>% 
  group_by(DOY, trtmt, arena) %>% 
  summarise(flowers = sum(flwr_count),
            tmpavg = mean(temp_avg))

flwrs <- flwrs %>% mutate(               #"un-nest" arenas from treatment
  new_arena = case_when(
    trtmt == "A" & arena == 1 ~ 1,
    trtmt == "A" & arena == 2 ~ 2,
    trtmt == "A" & arena == 3 ~ 3,
    trtmt == "A" & arena == 4 ~ 4,
    trtmt == "A" & arena == 5 ~ 5,
    trtmt == "W2" & arena == 1 ~ 6,
    trtmt == "W2" & arena == 2 ~ 7,
    trtmt == "W2" & arena == 3 ~ 8,
    trtmt == "W2" & arena == 4 ~ 9,
    trtmt == "W2" & arena == 5 ~ 10,
    trtmt == "W4" & arena == 1 ~ 11,
    trtmt == "W4" & arena == 2 ~ 12,
    trtmt == "W4" & arena == 3 ~ 13,
    trtmt == "W4" & arena == 4 ~ 14,
    trtmt == "W4" & arena == 5 ~ 15,
    trtmt == "W6" & arena == 1 ~ 16,
    trtmt == "W6" & arena == 2 ~ 17,
    trtmt == "W6" & arena == 3 ~ 18,
    trtmt == "W6" & arena == 4 ~ 19,
    trtmt == "W6" & arena == 5 ~ 20)
)

#model
glm_spl <- glmmTMB(flowers ~ tmpavg*splines::bs(DOY, degree = 2) + (1|new_arena) + (1|DOY), family = nbinom1(link = 'log'), data = flwrs) #using splines was the only method that allowed to get predictions by temperature
summary(glm_spl)

#plot
ifl <- interact_plot(glm_spl, pred = DOY, modx = tmpavg, modx.values = c(24.4, 25.8, 26.9, 30.4),
                     data = flwrs, interval = FALSE,
                     x.label = "", y.label = "Flower abundance") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title=element_text(size=38),
        legend.title=element_text(size=0), 
        legend.text=element_text(size=24)) #plot interaction effects

ifl

##########
## bees ##
##########
bee_count <- read.table("bee_abund_2022.txt", sep = "\t", header = T)

flwr_avg2 <- subset(flwr_count, arena != "TR") %>% 
  group_by(DOY, trtmt, arena) %>% 
  summarise(tmpavg = mean(temp_avg))

flwr_avg2$arena <- as.numeric(flwr_avg2$arena)


bees <- left_join(flwr_avg2, bee_count)

bees[is.na(bees)] <- 0 #replace NAs with 0s for DOYs in bees to include whole observation period (i.e., flowers open but no active bees)

bees <- bees %>% mutate(               #"un-nest" arenas from treatment
  new_arena = case_when(
    trtmt == "A" & arena == 1 ~ 1,
    trtmt == "A" & arena == 2 ~ 2,
    trtmt == "A" & arena == 3 ~ 3,
    trtmt == "A" & arena == 4 ~ 4,
    trtmt == "A" & arena == 5 ~ 5,
    trtmt == "W2" & arena == 1 ~ 6,
    trtmt == "W2" & arena == 2 ~ 7,
    trtmt == "W2" & arena == 3 ~ 8,
    trtmt == "W2" & arena == 4 ~ 9,
    trtmt == "W2" & arena == 5 ~ 10,
    trtmt == "W4" & arena == 1 ~ 11,
    trtmt == "W4" & arena == 2 ~ 12,
    trtmt == "W4" & arena == 3 ~ 13,
    trtmt == "W4" & arena == 4 ~ 14,
    trtmt == "W4" & arena == 5 ~ 15,
    trtmt == "W6" & arena == 1 ~ 16,
    trtmt == "W6" & arena == 2 ~ 17,
    trtmt == "W6" & arena == 3 ~ 18,
    trtmt == "W6" & arena == 4 ~ 19,
    trtmt == "W6" & arena == 5 ~ 20)
)

#model
glm_spl_bee <- glmmTMB(bees ~ tmpavg*splines::bs(DOY, degree = 2) + (1|new_arena) + (1|DOY), family = poisson(link = 'log'), data = bees) #using splines was the only method that allowed to get predictions by temperature
summary(glm_spl_bee)

#plot
ib <- interact_plot(glm_spl_bee, pred = DOY, modx = tmpavg, modx.values = c(24.4, 25.8, 26.9, 30.4),
                    data = bees, interval = FALSE,
                    x.label = "Day of year", y.label = "Bee abundance") +
  theme_classic() +
  theme(axis.text=element_text(size=28), axis.title=element_text(size=38),
        legend.title=element_text(size=0), 
        legend.text=element_text(size=24)) #plot interaction effects

ib


---------------------------------------
3c. Phenology overlap (Table S3)

library(tidyr)
library(dplyr)

#prepare data
#flowers
flwr_count <- read.table("flwrcount_temp2022.txt", sep = "\t", header = T)

flwr_count <- separate(flwr_count, "trtmt_arena", into = c("trtmt","arena"), sep = "_") #separate trtmt and arena

flwr_sum <- subset(flwr_count, arena != "TR") %>% 
  group_by(DOY, trtmt, arena) %>% #clumped species to look at plant community effect
  summarise(flowers = sum(flwr_count),
            tmpavg = mean(temp_avg))


#bees
bee_count <- read.table("bee_abund_2022.txt", sep = "\t", header = T)

bee_count$arena <- as.character(bee_count$arena)

flwr_avg2 <- subset(flwr_count, arena != "TR") %>% #extract average temperature per arena
  group_by(DOY, trtmt, arena) %>% 
  summarise(tmpavg = mean(temp_avg))

bees <- left_join(bee_count, flwr_avg2)

bee_sum <- bees %>% 
  group_by(DOY, trtmt, arena) %>% 
  summarise(bees = sum(bees),
            tmpavg = mean(tmpavg))


#check DOYs for flowers and bees
unique(flwr_sum$DOY)
unique(bee_sum$DOY)


#join flower and bee df, and format
flwr_bees <- flwr_sum %>% left_join(bee_sum)

flwr_bees[is.na(flwr_bees)] <- 0 #replace NAs with 0s for DOYs in bees

flwr_bees <- flwr_bees[, c(1,2,3,5,4,6)] #reorder cols

flwr_bees_long <- flwr_bees %>% pivot_longer(cols = flowers:bees, names_to = "type", values_to = "abund") #wide to long format


##calculate mean temperature for arena within treatment
arena_temp <- flwr_bees_long %>% 
  group_by(trtmt, arena) %>% 
  summarise(meantmp = mean(tmpavg))


#############
## OVERLAP ##
#############
library(overlap)

#Use estimator Dhat1 if smaller of two samples has size <50, otherwise use Dhat4 (see Ridout & Linkie 2009 - Estimating overlap of daily activity patterns from camera trap data, Meredith & Ridout 2018 - Overview of the overlap package, Meredith et al. 2024 - R-package overlap)

flwr_bees_long$rad <- ((flwr_bees_long$DOY/365) * 360) * pi / 180 #convert doy to radians

## Ambient
#Arena 1
A1_flwr <- subset(flwr_bees_long, trtmt == "A" & arena == "1" & type == "flowers")
A1_flwr_rad <- A1_flwr[, c(6,7)]
A1_flwr_exp <- A1_flwr_rad[rep(seq_len(dim(A1_flwr_rad)[1]), A1_flwr_rad$abund),] #expand radians in rows based on abundance
densityPlot(A1_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A1_bees <- subset(flwr_bees_long, trtmt == "A" & arena == "1" & type == "bees")
A1_bees_rad <- A1_bees[, c(6,7)]
A1_bees_exp <- A1_bees_rad[rep(seq_len(dim(A1_bees_rad)[1]), A1_bees_rad$abund),] #expand radians in rows based on abundance  
densityPlot(A1_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A1_ovrlp <- overlapEst(A1_flwr_exp$rad, A1_bees_exp$rad, type="Dhat4")
overlapPlot(A1_flwr_exp$rad, A1_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="Ambient - arena 1", xlab = "doyrad")

#Arena 2
A2_flwr <- subset(flwr_bees_long, trtmt == "A" & arena == "2" & type == "flowers")
A2_flwr_rad <- A2_flwr[, c(6,7)]
A2_flwr_exp <- A2_flwr_rad[rep(seq_len(dim(A2_flwr_rad)[1]), A2_flwr_rad$abund),] 
densityPlot(A2_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A2_bees <- subset(flwr_bees_long, trtmt == "A" & arena == "2" & type == "bees")
A2_bees_rad <- A2_bees[, c(6,7)]
A2_bees_exp <- A2_bees_rad[rep(seq_len(dim(A2_bees_rad)[1]), A2_bees_rad$abund),]  
densityPlot(A2_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A2_ovrlp <- overlapEst(A2_flwr_exp$rad, A2_bees_exp$rad, type="Dhat1")
overlapPlot(A2_flwr_exp$rad, A2_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="Ambient - arena 2", xlab = "doyrad")

#Arena 3
A3_flwr <- subset(flwr_bees_long, trtmt == "A" & arena == "3" & type == "flowers")
A3_flwr_rad <- A3_flwr[, c(6,7)]
A3_flwr_exp <- A3_flwr_rad[rep(seq_len(dim(A3_flwr_rad)[1]), A3_flwr_rad$abund),] 
densityPlot(A3_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A3_bees <- subset(flwr_bees_long, trtmt == "A" & arena == "3" & type == "bees")
A3_bees_rad <- A3_bees[, c(6,7)]
A3_bees_exp <- A3_bees_rad[rep(seq_len(dim(A3_bees_rad)[1]), A3_bees_rad$abund),]  
densityPlot(A3_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A3_ovrlp <- overlapEst(A3_flwr_exp$rad, A3_bees_exp$rad, type="Dhat4")
overlapPlot(A3_flwr_exp$rad, A3_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="Ambient - arena 3", xlab = "doyrad")

#Arena 4
A4_flwr <- subset(flwr_bees_long, trtmt == "A" & arena == "4" & type == "flowers")
A4_flwr_rad <- A4_flwr[, c(6,7)]
A4_flwr_exp <- A4_flwr_rad[rep(seq_len(dim(A4_flwr_rad)[1]), A4_flwr_rad$abund),] 
densityPlot(A4_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A4_bees <- subset(flwr_bees_long, trtmt == "A" & arena == "4" & type == "bees")
A4_bees_rad <- A4_bees[, c(6,7)]
A4_bees_exp <- A4_bees_rad[rep(seq_len(dim(A4_bees_rad)[1]), A4_bees_rad$abund),]  
densityPlot(A4_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A4_ovrlp <- overlapEst(A4_flwr_exp$rad, A4_bees_exp$rad, type="Dhat4")
overlapPlot(A4_flwr_exp$rad, A4_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="Ambient - arena 4", xlab = "doyrad")

#Arena 5
A5_flwr <- subset(flwr_bees_long, trtmt == "A" & arena == "5" & type == "flowers")
A5_flwr_rad <- A5_flwr[, c(6,7)]
A5_flwr_exp <- A5_flwr_rad[rep(seq_len(dim(A5_flwr_rad)[1]), A5_flwr_rad$abund),] 
densityPlot(A5_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A5_bees <- subset(flwr_bees_long, trtmt == "A" & arena == "5" & type == "bees")
A5_bees_rad <- A5_bees[, c(6,7)]
A5_bees_exp <- A5_bees_rad[rep(seq_len(dim(A5_bees_rad)[1]), A5_bees_rad$abund),]  
densityPlot(A5_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

A5_ovrlp <- overlapEst(A5_flwr_exp$rad, A5_bees_exp$rad, type="Dhat4")
overlapPlot(A5_flwr_exp$rad, A5_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="Ambient - arena 5", xlab = "doyrad")


## W2
#Arena 1
W2_1_flwr <- subset(flwr_bees_long, trtmt == "W2" & arena == "1" & type == "flowers")
W2_1_flwr_rad <- W2_1_flwr[, c(6,7)]
W2_1_flwr_exp <- W2_1_flwr_rad[rep(seq_len(dim(W2_1_flwr_rad)[1]), W2_1_flwr_rad$abund),] #expand radians in rows based on abundance
densityPlot(W2_1_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_1_bees <- subset(flwr_bees_long, trtmt == "W2" & arena == "1" & type == "bees")
W2_1_bees_rad <- W2_1_bees[, c(6,7)]
W2_1_bees_exp <- W2_1_bees_rad[rep(seq_len(dim(W2_1_bees_rad)[1]), W2_1_bees_rad$abund),] #expand radians in rows based on abundance  
densityPlot(W2_1_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_1_ovrlp <- overlapEst(W2_1_flwr_exp$rad, W2_1_bees_exp$rad, type="Dhat1")
overlapPlot(W2_1_flwr_exp$rad, W2_1_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W2 - arena 1", xlab = "doyrad")

#Arena 2
W2_2_flwr <- subset(flwr_bees_long, trtmt == "W2" & arena == "2" & type == "flowers")
W2_2_flwr_rad <- W2_2_flwr[, c(6,7)]
W2_2_flwr_exp <- W2_2_flwr_rad[rep(seq_len(dim(W2_2_flwr_rad)[1]), W2_2_flwr_rad$abund),] 
densityPlot(W2_2_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_2_bees <- subset(flwr_bees_long, trtmt == "W2" & arena == "2" & type == "bees")
W2_2_bees_rad <- W2_2_bees[, c(6,7)]
W2_2_bees_exp <- W2_2_bees_rad[rep(seq_len(dim(W2_2_bees_rad)[1]), W2_2_bees_rad$abund),]   
densityPlot(W2_2_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_2_ovrlp <- overlapEst(W2_2_flwr_exp$rad, W2_2_bees_exp$rad, type="Dhat1")
overlapPlot(W2_2_flwr_exp$rad, W2_2_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W2 - arena 2", xlab = "doyrad")

#Arena 3
W2_3_flwr <- subset(flwr_bees_long, trtmt == "W2" & arena == "3" & type == "flowers")
W2_3_flwr_rad <- W2_3_flwr[, c(6,7)]
W2_3_flwr_exp <- W2_3_flwr_rad[rep(seq_len(dim(W2_3_flwr_rad)[1]), W2_3_flwr_rad$abund),] 
densityPlot(W2_3_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_3_bees <- subset(flwr_bees_long, trtmt == "W2" & arena == "3" & type == "bees")
W2_3_bees_rad <- W2_3_bees[, c(6,7)]
W2_3_bees_exp <- W2_3_bees_rad[rep(seq_len(dim(W2_3_bees_rad)[1]), W2_3_bees_rad$abund),]   
densityPlot(W2_3_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_3_ovrlp <- overlapEst(W2_3_flwr_exp$rad, W2_3_bees_exp$rad, type="Dhat1")
overlapPlot(W2_3_flwr_exp$rad, W2_3_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W2 - arena 3", xlab = "doyrad")

#Arena 4
W2_4_flwr <- subset(flwr_bees_long, trtmt == "W2" & arena == "4" & type == "flowers")
W2_4_flwr_rad <- W2_4_flwr[, c(6,7)]
W2_4_flwr_exp <- W2_4_flwr_rad[rep(seq_len(dim(W2_4_flwr_rad)[1]), W2_4_flwr_rad$abund),] 
densityPlot(W2_4_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_4_bees <- subset(flwr_bees_long, trtmt == "W2" & arena == "4" & type == "bees")
W2_4_bees_rad <- W2_4_bees[, c(6,7)]
W2_4_bees_exp <- W2_4_bees_rad[rep(seq_len(dim(W2_4_bees_rad)[1]), W2_4_bees_rad$abund),]   
densityPlot(W2_4_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_4_ovrlp <- overlapEst(W2_4_flwr_exp$rad, W2_4_bees_exp$rad, type="Dhat1")
overlapPlot(W2_4_flwr_exp$rad, W2_4_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W2 - arena 4", xlab = "doyrad")

#Arena 5
W2_5_flwr <- subset(flwr_bees_long, trtmt == "W2" & arena == "5" & type == "flowers")
W2_5_flwr_rad <- W2_5_flwr[, c(6,7)]
W2_5_flwr_exp <- W2_5_flwr_rad[rep(seq_len(dim(W2_5_flwr_rad)[1]), W2_5_flwr_rad$abund),] 
densityPlot(W2_5_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_5_bees <- subset(flwr_bees_long, trtmt == "W2" & arena == "5" & type == "bees")
W2_5_bees_rad <- W2_5_bees[, c(6,7)]
W2_5_bees_exp <- W2_5_bees_rad[rep(seq_len(dim(W2_5_bees_rad)[1]), W2_5_bees_rad$abund),]   
densityPlot(W2_5_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W2_5_ovrlp <- overlapEst(W2_5_flwr_exp$rad, W2_5_bees_exp$rad, type="Dhat1")
overlapPlot(W2_5_flwr_exp$rad, W2_5_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W2 - arena 5", xlab = "doyrad")


## W4
#Arena 1
W4_1_flwr <- subset(flwr_bees_long, trtmt == "W4" & arena == "1" & type == "flowers")
W4_1_flwr_rad <- W4_1_flwr[, c(6,7)]
W4_1_flwr_exp <- W4_1_flwr_rad[rep(seq_len(dim(W4_1_flwr_rad)[1]), W4_1_flwr_rad$abund),] #expand radians in rows based on abundance
densityPlot(W4_1_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_1_bees <- subset(flwr_bees_long, trtmt == "W4" & arena == "1" & type == "bees")
W4_1_bees_rad <- W4_1_bees[, c(6,7)]
W4_1_bees_exp <- W4_1_bees_rad[rep(seq_len(dim(W4_1_bees_rad)[1]), W4_1_bees_rad$abund),] #expand radians in rows based on abundance  
densityPlot(W4_1_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_1_ovrlp <- overlapEst(W4_1_flwr_exp$rad, W4_1_bees_exp$rad, type="Dhat1")
overlapPlot(W4_1_flwr_exp$rad, W4_1_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W4 - arena 1", xlab = "doyrad")

#Arena 2
W4_2_flwr <- subset(flwr_bees_long, trtmt == "W4" & arena == "2" & type == "flowers")
W4_2_flwr_rad <- W4_2_flwr[, c(6,7)]
W4_2_flwr_exp <- W4_2_flwr_rad[rep(seq_len(dim(W4_2_flwr_rad)[1]), W4_2_flwr_rad$abund),]
densityPlot(W4_2_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_2_bees <- subset(flwr_bees_long, trtmt == "W4" & arena == "2" & type == "bees")
W4_2_bees_rad <- W4_2_bees[, c(6,7)]
W4_2_bees_exp <- W4_2_bees_rad[rep(seq_len(dim(W4_2_bees_rad)[1]), W4_2_bees_rad$abund),]  
densityPlot(W4_2_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_2_ovrlp <- overlapEst(W4_2_flwr_exp$rad, W4_2_bees_exp$rad, type="Dhat1")
overlapPlot(W4_2_flwr_exp$rad, W4_2_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W4 - arena 2", xlab = "doyrad")

#Arena 3
W4_3_flwr <- subset(flwr_bees_long, trtmt == "W4" & arena == "3" & type == "flowers")
W4_3_flwr_rad <- W4_3_flwr[, c(6,7)]
W4_3_flwr_exp <- W4_3_flwr_rad[rep(seq_len(dim(W4_3_flwr_rad)[1]), W4_3_flwr_rad$abund),]
densityPlot(W4_3_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_3_bees <- subset(flwr_bees_long, trtmt == "W4" & arena == "3" & type == "bees")
W4_3_bees_rad <- W4_3_bees[, c(6,7)]
W4_3_bees_exp <- W4_3_bees_rad[rep(seq_len(dim(W4_3_bees_rad)[1]), W4_3_bees_rad$abund),]  
densityPlot(W4_3_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_3_ovrlp <- overlapEst(W4_3_flwr_exp$rad, W4_3_bees_exp$rad, type="Dhat1")
overlapPlot(W4_3_flwr_exp$rad, W4_3_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W4 - arena 3", xlab = "doyrad")

#Arena 4
W4_4_flwr <- subset(flwr_bees_long, trtmt == "W4" & arena == "4" & type == "flowers")
W4_4_flwr_rad <- W4_4_flwr[, c(6,7)]
W4_4_flwr_exp <- W4_4_flwr_rad[rep(seq_len(dim(W4_4_flwr_rad)[1]), W4_4_flwr_rad$abund),]
densityPlot(W4_4_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_4_bees <- subset(flwr_bees_long, trtmt == "W4" & arena == "4" & type == "bees")
W4_4_bees_rad <- W4_4_bees[, c(6,7)]
W4_4_bees_exp <- W4_4_bees_rad[rep(seq_len(dim(W4_4_bees_rad)[1]), W4_4_bees_rad$abund),]  
densityPlot(W4_4_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_4_ovrlp <- overlapEst(W4_4_flwr_exp$rad, W4_4_bees_exp$rad, type="Dhat1")
overlapPlot(W4_4_flwr_exp$rad, W4_4_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W4 - arena 4", xlab = "doyrad")

#Arena 5
W4_5_flwr <- subset(flwr_bees_long, trtmt == "W4" & arena == "5" & type == "flowers")
W4_5_flwr_rad <- W4_5_flwr[, c(6,7)]
W4_5_flwr_exp <- W4_5_flwr_rad[rep(seq_len(dim(W4_5_flwr_rad)[1]), W4_5_flwr_rad$abund),]
densityPlot(W4_5_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_5_bees <- subset(flwr_bees_long, trtmt == "W4" & arena == "5" & type == "bees")
W4_5_bees_rad <- W4_5_bees[, c(6,7)]
W4_5_bees_exp <- W4_5_bees_rad[rep(seq_len(dim(W4_5_bees_rad)[1]), W4_5_bees_rad$abund),]  
densityPlot(W4_5_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W4_5_ovrlp <- overlapEst(W4_5_flwr_exp$rad, W4_5_bees_exp$rad, type="Dhat1")
overlapPlot(W4_5_flwr_exp$rad, W4_5_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W4 - arena 5", xlab = "doyrad")


## W6
#Arena 1
W6_1_flwr <- subset(flwr_bees_long, trtmt == "W6" & arena == "1" & type == "flowers")
W6_1_flwr_rad <- W6_1_flwr[, c(6,7)]
W6_1_flwr_exp <- W6_1_flwr_rad[rep(seq_len(dim(W6_1_flwr_rad)[1]), W6_1_flwr_rad$abund),] #expand radians in rows based on abundance
densityPlot(W6_1_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_1_bees <- subset(flwr_bees_long, trtmt == "W6" & arena == "1" & type == "bees")
W6_1_bees_rad <- W6_1_bees[, c(6,7)]
W6_1_bees_exp <- W6_1_bees_rad[rep(seq_len(dim(W6_1_bees_rad)[1]), W6_1_bees_rad$abund),] #expand radians in rows based on abundance  
densityPlot(W6_1_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_1_ovrlp <- overlapEst(W6_1_flwr_exp$rad, W6_1_bees_exp$rad, type="Dhat1")
overlapPlot(W6_1_flwr_exp$rad, W6_1_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W6 - arena 1", xlab = "doyrad")

#Arena 2
W6_2_flwr <- subset(flwr_bees_long, trtmt == "W6" & arena == "2" & type == "flowers")
W6_2_flwr_rad <- W6_2_flwr[, c(6,7)]
W6_2_flwr_exp <- W6_2_flwr_rad[rep(seq_len(dim(W6_2_flwr_rad)[1]), W6_2_flwr_rad$abund),]
densityPlot(W6_2_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_2_bees <- subset(flwr_bees_long, trtmt == "W6" & arena == "2" & type == "bees")
W6_2_bees_rad <- W6_2_bees[, c(6,7)]
W6_2_bees_exp <- W6_2_bees_rad[rep(seq_len(dim(W6_2_bees_rad)[1]), W6_2_bees_rad$abund),]  
densityPlot(W6_2_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_2_ovrlp <- overlapEst(W6_2_flwr_exp$rad, W6_2_bees_exp$rad, type="Dhat1")
overlapPlot(W6_2_flwr_exp$rad, W6_2_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W6 - arena 2", xlab = "doyrad")

#Arena 3
W6_3_flwr <- subset(flwr_bees_long, trtmt == "W6" & arena == "3" & type == "flowers")
W6_3_flwr_rad <- W6_3_flwr[, c(6,7)]
W6_3_flwr_exp <- W6_3_flwr_rad[rep(seq_len(dim(W6_3_flwr_rad)[1]), W6_3_flwr_rad$abund),] 
densityPlot(W6_3_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_3_bees <- subset(flwr_bees_long, trtmt == "W6" & arena == "3" & type == "bees")
W6_3_bees_rad <- W6_3_bees[, c(6,7)]
W6_3_bees_exp <- W6_3_bees_rad[rep(seq_len(dim(W6_3_bees_rad)[1]), W6_3_bees_rad$abund),]
densityPlot(W6_3_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_3_ovrlp <- overlapEst(W6_3_flwr_exp$rad, W6_3_bees_exp$rad, type="Dhat1")
overlapPlot(W6_3_flwr_exp$rad, W6_3_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W6 - arena 3", xlab = "doyrad")

#Arena 4
W6_4_flwr <- subset(flwr_bees_long, trtmt == "W6" & arena == "4" & type == "flowers")
W6_4_flwr_rad <- W6_4_flwr[, c(6,7)]
W6_4_flwr_exp <- W6_4_flwr_rad[rep(seq_len(dim(W6_4_flwr_rad)[1]), W6_4_flwr_rad$abund),] 
densityPlot(W6_4_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_4_bees <- subset(flwr_bees_long, trtmt == "W6" & arena == "4" & type == "bees")
W6_4_bees_rad <- W6_4_bees[, c(6,7)]
W6_4_bees_exp <- W6_4_bees_rad[rep(seq_len(dim(W6_4_bees_rad)[1]), W6_4_bees_rad$abund),]   
densityPlot(W6_4_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_4_ovrlp <- overlapEst(W6_4_flwr_exp$rad, W6_4_bees_exp$rad, type="Dhat1")
overlapPlot(W6_4_flwr_exp$rad, W6_4_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W6 - arena 4", xlab = "doyrad")

#Arena 5
W6_5_flwr <- subset(flwr_bees_long, trtmt == "W6" & arena == "5" & type == "flowers")
W6_5_flwr_rad <- W6_5_flwr[, c(6,7)]
W6_5_flwr_exp <- W6_5_flwr_rad[rep(seq_len(dim(W6_5_flwr_rad)[1]), W6_5_flwr_rad$abund),] 
densityPlot(W6_5_flwr_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_5_bees <- subset(flwr_bees_long, trtmt == "W6" & arena == "5" & type == "bees")
W6_5_bees_rad <- W6_5_bees[, c(6,7)]
W6_5_bees_exp <- W6_5_bees_rad[rep(seq_len(dim(W6_5_bees_rad)[1]), W6_5_bees_rad$abund),]   
densityPlot(W6_5_bees_exp$rad, xscale = NA, rug=TRUE, extend = NULL, adjust = 2)

W6_5_ovrlp <- overlapEst(W6_5_flwr_exp$rad, W6_5_bees_exp$rad, type="Dhat1")
overlapPlot(W6_5_flwr_exp$rad, W6_5_bees_exp$rad, rug = F, xscale = NA, adjust = 2, main="W6 - arena 5", xlab = "doyrad")


tot_ovrlp <- data.frame(c(A1_ovrlp[[1]], A2_ovrlp[[1]], A3_ovrlp[[1]], A4_ovrlp[[1]], A5_ovrlp[[1]],  #combine all overlap values (without names)
               W2_1_ovrlp[[1]], W2_2_ovrlp[[1]], W2_3_ovrlp[[1]], W2_4_ovrlp[[1]], W2_5_ovrlp[[1]],
               W4_1_ovrlp[[1]], W4_2_ovrlp[[1]], W4_3_ovrlp[[1]], W4_4_ovrlp[[1]], W4_5_ovrlp[[1]],
               W6_1_ovrlp[[1]], W6_2_ovrlp[[1]], W6_3_ovrlp[[1]], W6_4_ovrlp[[1]], W6_5_ovrlp[[1]]))

colnames(tot_ovrlp) <- "overlap" #rename column

temp_ovrlp <- cbind(arena_temp, tot_ovrlp) #bind overlap to temperatures


##model
library(glmmTMB)
temp_ovrlp$meantmp_sq <- (temp_ovrlp$meantmp)^2 #add square term

beta_glm <- glmmTMB(overlap ~ meantmp + meantmp_sq, family = beta_family(link = "logit"), data = temp_ovrlp)
summary(beta_glm)


--------------------------------------------------------------------------------------------------------------------
4. Data imputations

library(tidyr)
library(dplyr)
library(Amelia)

vars_tot <- read.table("vars_tot2022.txt", sep = "\t", header = T) 
vars_tot <- vars_tot[, -8]


## imputations ##
T_imp <- amelia(vars_tot, m = 100, ts = "DOY", cs = "trtmt_sp", polytime = 2, intercs = TRUE, p2s = 2, empri = .01 * nrow(vars_tot), logs = c("flwr_count", "visits", "hand_freq", "width1", "width2", "vol", "conc", "bees"))
summary(T_imp)
plot(T_imp) 

#save output list
save(T_imp, file = "imputations_temporal_tot.RData")

## average imputed datasets ##

#Tot
#extract and bind columns for each variable
col12 <- as.data.frame(T_imp$imputations[[1]])
info = as.data.frame(col12[,1:2])

flwr_count = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
visits = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
hand_freq = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
width1 = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
width2 = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
vol = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
conc = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)
bees = data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp)

for (i in 1:length(T_imp$imputations)) {
  df <- as.data.frame(T_imp$imputations[[i]])
  
  flower_count_i = as.data.frame(df[,3])
  colnames(flower_count_i) = 'flower_count'
  
  visits_i = as.data.frame(df[,4])
  colnames(visits_i) = 'visits'
  
  hand_freq_i = as.data.frame(df[,5])
  colnames(hand_freq_i) = 'hand_freq'
  
  width1_i = as.data.frame(df[,6])
  colnames(width1_i) = 'width1'
  
  width2_i = as.data.frame(df[,7])
  colnames(width2_i) = 'width2'
  
  vol_i = as.data.frame(df[,8])
  colnames(vol_i) = 'vol'
  
  conc_i = as.data.frame(df[,9])
  colnames(conc_i) = 'conc'
  
  bees_i = as.data.frame(df[,10])
  colnames(bees_i) = 'bees'
  
  flwr_count = cbind(flwr_count, flower_count_i)
  visits = cbind(visits, visits_i)
  hand_freq = cbind(hand_freq, hand_freq_i)
  width1 = cbind(width1, width1_i)
  width2 = cbind(width2, width2_i)
  vol = cbind(vol, vol_i)
  conc = cbind(conc, conc_i)
  bees = cbind(bees, bees_i)
  
  print(i)
  
}

#remove from the environment unnecessary files
rm(df, flower_count_i, visits_i, hand_freq_i, width1_i, width2_i, vol_i, conc_i, bees_i)


#average variables and bind
imput_vars_T <- data.frame(doy = info$DOY, trtmt_sp = info$trtmt_sp, 
                             flwrcount_mean = rowMeans(flwr_count[,3:101]),
                             visits_mean = rowMeans(visits[,3:101]),
                             handling_mean = rowMeans(hand_freq[,3:101]),
                             width1 = rowMeans(width1[,3:101]),
                             width2 = rowMeans(width2[,3:101]),
                             vol_mean = rowMeans(vol[,3:101]),
                             conc_mean = rowMeans(conc[,3:101]),
                             bees_mean = rowMeans(bees[,3:101])
)


#round flower count and bees
imput_vars_T <- imput_vars_T %>% mutate_at(3, round, 0)
imput_vars_T <- imput_vars_T %>% mutate_at(10, round, 0)

#save table
write.table(imput_vars_T, file = 'imputed_temporal_vars.txt', sep = "\t")


## merge temperature ##
imputed_T <- read.table("imputed_temporal_vars.txt", sep = "\t", header = T)

temp <- read.table("vars_tot_temp.txt", sep = "\t", header = T)

imputed_T$doy_trtmt_sp <- paste(imputed_T$doy, imputed_T$trtmt_sp, sep="_")

temp$doy_trtmt_sp <- paste(temp$doy, temp$trtmt_sp, sep="_")

imputed_T <- imputed_T[,c(11,3:8,10)]

temp <- temp[,c(4,3)]

it_T <- cbind(imputed_T, temp)
it_T <- it_T[,c(1,10,2:8)]
it_T <- separate(it_T, "doy_trtmt_sp", into = c("doy","trtmt", "sp"), sep = "_")



##calculate flower landing surface area
it_T <- it_T %>% mutate(
  surf = case_when(sp == "NM" | sp == "PC" ~ (width1/2)*(width2/2)*pi,
                   sp == "CH" ~ (width1/2)*(width2/2)*pi*2))

it_T <- it_T[, c(1:7,12,10,11)] #remove unwanted columns

#rename columns
it_T <- it_T %>% rename(temp = temp_avg, flwrcount = flwrcount_mean, visits = visits_mean, handling = handling_mean, vol = vol_mean, bees = bees_mean)

#save table
write.table(it_T, file = 'imputed_temporal_vars_T.txt', sep = "\t")



--------------------------------------------------------------------------------------------------------------------------------------------
5. Dynamic pSEM (Table S4, Figure 3)

library(glmmTMB)
library(piecewiseSEM)

imputed_T <- read.table("imputed_temporal_vars_T.txt", sep = "\t", header = T)

#transform data
imputed_T$visits_log <- log(imputed_T$visits + 1)
imputed_T$handling_log <- log(imputed_T$handling + 1)
imputed_T$surf_log <- log(imputed_T$surf)
imputed_T$vol_log <- log(imputed_T$vol+1)


##pSEM on temporal data
pSEM_t <- psem(
  
  glmmTMB(flwrcount ~ temp + (1|doy) + (1|sp) + (1|arena), family = poisson(link = 'log'), data = imputed_T),
  
  glmmTMB(bees ~ temp + (1|doy) + (1|arena), family = poisson(link = 'log'), data = imputed_T),
  
  glmmTMB(visits_log ~ flwrcount + bees + surf_log + vol_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T),
  
  glmmTMB(handling_log ~ flwrcount + bees + visits_log + surf_log + vol_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T),
  
  glmmTMB(surf_log ~ temp + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T),
  
  glmmTMB(vol_log ~ temp + surf_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T)
  
)

summary(pSEM_t)


#update model
pSEM_t2 <- update(pSEM_t, surf_log %~~% bees, flwrcount %~~% bees, visits_log %~~% temp, handling_log %~~% temp, handling_log %~~% surf_log, vol_log %~~% bees, surf_log %~~% flwrcount)
summary(pSEM_t2)


##########################################
#Regression plots (panels B-G in Figure 3)
library(ggeffects)
library(ggplot2)

bees_poi_tmb <- glmmTMB(bees ~ temp + (1|doy) + (1|arena), family = poisson(link = 'log'), data = imputed_T)
vsts_tmb <- glmmTMB(visits_log ~ flwrcount + bees + surf_log + vol_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T)
hand_tmb <- glmmTMB(handling_log ~ flwrcount + bees + visits_log + surf_log + vol_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T)
surf_tmb <- glmmTMB(surf_log ~ temp + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T)
vol_tmb <- glmmTMB(vol_log ~ temp + surf_log + (1|doy) + (1|sp) + (1|arena), family = gaussian(link = 'identity'), data = imputed_T)


pred_beeabnd <- predict_response(bees_poi_tmb, "temp[all]", bias_correction = TRUE)
plot(pred_beeabnd, show_data = F, color = "#f5a134") +
  geom_line(linewidth = 3, color = "#f5a134") +
  labs(x = "Temperature (°C)", y = "Bee abundance", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))
  
  
pred_vis <- predict_response(vsts_tmb, "bees[all]", bias_correction = TRUE)
ggplot(pred_vis, aes(x = x, y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#488cdf") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#488cdf", alpha = 0.2) +
  labs(x = "Bee abundance", y = "No. flower visits", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_hand <- predict_response(hand_tmb, "visits_log[all]", bias_correction = TRUE)
ggplot(pred_hand, aes(x = exp(x), y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#488cdf") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#488cdf", alpha = 0.2) +
  labs(x = "No. flower visits", y = "Flower handling time", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_surf <- predict_response(surf_tmb, "temp", bias_correction = TRUE)
ggplot(pred_surf, aes(x = x, y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Temperature (°C)", y = "Flower size", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_vol <- predict_response(vol_tmb, "temp", bias_correction = TRUE)
ggplot(pred_vol, aes(x = x, y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Temperature (°C)", y = "Nectar volume", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_vis_vol <- predict_response(vsts_tmb, "vol_log[0:4]", bias_correction = TRUE)
ggplot(pred_vis_vol, aes(x = x, y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Nectar volume", y = "No. flower visits", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))




----------------------------------------------------------------------------------------------------------
5. pSEM for plant fitness (Table S6, Figure 4)

library(tidyr)
library(glmmTMB)
library(piecewiseSEM)

seeds_germ_pSEM <- read.table("data_seeds_germ_pSEM.txt", sep = "\t", header = T)
seeds_germ_pSEM <- separate(seeds_germ_pSEM, "plant_flwr_ID", into = c("plant_ID","flwr_ID"), sep = "_")
seeds_germ_pSEM$mort <- seeds_germ_pSEM$germ - seeds_germ_pSEM$cotyl #seed mortality

#transform data
seeds_germ_pSEM$visits_avg_log <- log(seeds_germ_pSEM$visits_avg + 1)
seeds_germ_pSEM$handling_avg_log <- log(seeds_germ_pSEM$hand_avg + 1)
seeds_germ_pSEM$surf_avg_log <- log(seeds_germ_pSEM$surf)
seeds_germ_pSEM$mass_seed_log <- log(seeds_germ_pSEM$mass_seed)


## Piecewise SEM on plant fitness ##
pSEM_b <- psem(
  
  glmmTMB(flwr_avg ~ temp_avg + (1|sp), family = poisson(link = "log"), data = seeds_germ_pSEM),
  
  glmmTMB(bee_abund ~ temp_avg + (1|arena), family = poisson(link = "log"), data = seeds_germ_pSEM),
  
  glmmTMB(visits_avg_log ~ bee_abund + flwr_avg + (1|sp), family = gaussian(link = 'identity'), data = seeds_germ_pSEM),
  
  glmmTMB(ovules_tot ~ temp_avg + (1|sp) , family = poisson(link = "log"), data = seeds_germ_pSEM),
  
  glmmTMB(seeds ~ temp_avg + ovules_tot + visits_avg_log + (1|sp), family = poisson(link = "log"), data = seeds_germ_pSEM),
  
  glmmTMB(germ ~ temp_avg + visits_avg_log + seeds + (1|arena), family = poisson(link = "log"), data = seeds_germ_pSEM),
  
  glmmTMB(mort ~ temp_avg + visits_avg_log + seeds + germ + (1|arena), family = poisson(link = "log"), data = seeds_germ_pSEM)
)
  
summary(pSEM_b)

#update model
pSEM_b2 <- update(pSEM_b, bee_abund %~~% flwr_avg, bee_abund %~~% ovules_tot, bee_abund %~~% seeds, ovules_tot %~~% flwr_avg, seeds %~~% flwr_avg, visits_avg_log %~~% ovules_tot)
summary(pSEM_b2)


##########################################
#Regression plots (panels B-G in Figure 4)
library(ggeffects)
library(ggplot2)

sset <- glmmTMB(seeds ~ temp_avg + ovules_tot + visits_avg_log + (1|sp), family = poisson(link = "log"), data = seeds_germ_pSEM)
grm_poi <- glmmTMB(germ ~ temp_avg + visits_avg_log + seeds + (1|arena), family = poisson(link = "log"), data = seeds_germ_pSEM)
mor_poi <- glmmTMB(mort ~ temp_avg + visits_avg_log + seeds + germ + (1|arena), family = poisson(link = "log"), data = seeds_germ_pSEM)


pred_sset <- predict_response(sset, "temp_avg")
plot(pred_sset, show_data = F, color = "#f5a134") +
  geom_line(linewidth = 3, color = "#f5a134") +
  labs(x = "Temperature (°C)", y = "Nb seeds", title = "") +
  scale_y_continuous(breaks = seq(10, 20, by = 5)) +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))

pred_sset_ov <- predict_response(sset, "ovules_tot")
plot(pred_sset_ov, show_data = F, color = "#488cdf") +
  geom_line(linewidth = 3, color = "#488cdf") +
  labs(x = "Nb ovules", y = "Nb seeds", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_vstgerm <- predict_response(grm_poi, "visits_avg_log[all]")
ggplot(pred_vstgerm, aes(x = exp(x), y = predicted)) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Nb flower visits", y = "Seed germination", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_vismort <- predict_response(mor_poi, "visits_avg_log[all]")
ggplot(pred_vismort, aes(x = exp(x), y = predicted)) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Nb flower visits", y = "Seedling mortality", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_tmpgerm <- predict_response(grm_poi, "temp_avg")
plot(pred_tmpgerm, show_data = F, color = "#488cdf") +
  geom_line(linewidth = 3, color = "#488cdf") +
  labs(x = "Temperature (°C)", y = "Seed germination", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_tmpmort <- predict_response(mor_poi, "temp_avg")
plot(pred_tmpmort, show_data = F, color = "#488cdf") +
  geom_line(linewidth = 3, color = "#488cdf") +
  labs(x = "Temperature (°C)", y = "Seedling mortality", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


------------------------------------------------------------------------------------------------------------------
6. pSEM for bee fitness (Table S7, Figure 5)

library(tidyr)
library(glmmTMB)
library(piecewiseSEM)

bee_nest_pSEM <- read.table("data_bee_nest_pSEM.txt", sep = "\t", header = T)
bee_nest_pSEM <- separate(bee_nest_pSEM, "plant_flwr_ID", into = c("plant_ID","flwr_ID"), sep = "_")

#transform data
bee_nest_pSEM$visits_avg_log <- log(bee_nest_pSEM$visits_avg + 1)
bee_nest_pSEM$handling_avg_log <- log(bee_nest_pSEM$hand_avg + 1)
bee_nest_pSEM$vol_log <- log(bee_nest_pSEM$vol + 1)


## Piecewise SEM on bee fitness ##
pSEM_c <- psem(
  
  glmmTMB(flwr_avg ~ temp_avg + (1|sp), family = poisson(link = "log"), data = bee_nest_pSEM),
  
  glmmTMB(bee_abund ~ temp_avg + (1|arena), family = poisson(link = "log"), data = bee_nest_pSEM),
  
  glmmTMB(visits_avg_log ~ bee_abund + flwr_avg + vol + (1|sp), family = gaussian(link = 'identity'), data = bee_nest_pSEM),
  
  glmmTMB(vol ~ temp_avg + (1|sp) , family = gaussian(link = 'identity'), data = bee_nest_pSEM), 
  
  glmmTMB(handling_avg_log ~ vol + visits_avg_log + (1|sp), family = gaussian(link = 'identity'), data = bee_nest_pSEM),
  
  glmmTMB(cells ~ bee_abund + visits_avg_log + handling_avg_log + (1|arena), family = poisson(link = "log"), data = bee_nest_pSEM),
  
  glmmTMB(cocoons ~ temp_avg + visits_avg_log + handling_avg_log + cells + (1|arena), family = poisson(link = "log"), data = bee_nest_pSEM)

  )

summary(pSEM_c)


#update model
pSEM_c2 <- update(pSEM_c, bee_abund %~~% flwr_avg, cells %~~% flwr_avg, cells %~~% vol)
summary(pSEM_c2)


##########################################
#Regression plots (panels B-G in Figure 5)
library(ggeffects)
library(ggplot2)

hnd <- glmmTMB(handling_avg_log ~ vol + visits_avg_log + (1|sp), family = gaussian(link = 'identity'), data = bee_nest_pSEM)
cll <- glmmTMB(cells ~ bee_abund + visits_avg_log + handling_avg_log + (1|arena), family = poisson(link = "log"), data = bee_nest_pSEM)


pred_hand <- predict_response(hnd, "vol")
ggplot(pred_hand, aes(x = x, y = exp(predicted))) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Nectar volume", y = "Flower handling time", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_cll_hand <- predict_response(cll, "handling_avg_log[all]")
ggplot(pred_cll_hand, aes(x = exp(x), y = predicted)) +   
  geom_line(linewidth = 3, color = "#f5a134") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "#f5a134", alpha = 0.2) +
  labs(x = "Flower handling time", y = "Nb nest cells", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_cll_bee <- predict_response(cll, "bee_abund[all]")
plot(pred_cll_bee, show_data = F, color = "#488cdf") +
  geom_line(linewidth = 3, color = "#488cdf") +
  labs(x = "Bee abundance", y = "Nb nest cells", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


pred_cll_vis <- predict_response(cll, "visits_avg_log[all]")
ggplot(pred_cll_vis, aes(x = exp(x), y = predicted)) +   
  geom_line(linewidth = 3, color = "#488cdf") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "#488cdf", alpha = 0.2) +
  labs(x = "Nb flower visits", y = "Nb nest cells", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


ggplot(bee_nest_pSEM, aes(x = cells, y = cocoons)) +
  stat_smooth(method = "glm", formula = y ~ x, method.args = list(family = "poisson"), se = TRUE, level = 0.99, color = "#488cdf", fill = "#488cdf", linewidth = 3, alpha = 0.2) +
  labs(x = "Nb nest cells", y = "Nb cocoons", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))


ggplot(bee_nest_pSEM, aes(x = temp_avg, y = cocoons)) +
  stat_smooth(method = "glm", formula = y ~ x, method.args = list(family = "poisson"), se = TRUE, level = 0.99, color = "#f5a134", fill = "#f5a134", linewidth = 3, alpha = 0.2) +
  labs(x = "Temperature (°C)", y = "Nb cocoons", title = "") +
  theme_classic() +
  theme(axis.text.x = element_text(size=40, face = "bold"),
        axis.text.y = element_text(size=40, face = "bold")) +
  theme(axis.title.x = element_text(size=46, face = "bold"),
        axis.title.y = element_text(size=46, face = "bold"))