library(ggpubr)
library(lme4)
library(lmerTest)
library(emmeans)
library(dplyr)

#-----------------------#
## linear correlations ##
#-----------------------#

#Linear correlations between the number of flowers visited and the time spent on flowers within a single whorl and on the whole plant

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

behav$sub_sex <- as.factor(behav$sub_sex)

behav_cuckoo <- subset(behav, subgenus == "C")
behav_trueM <- subset(behav, sub_sex == "mT")
behav_trueW <- subset(behav, sub_sex == "wT")

#whorl_tot <- cor.test(behav$n_flowers_per_whorl, behav$timevisit_whorl, method=c("pearson"))
#plant_tot <- cor.test(behav$n_flowers_total, behav$timevisit_plant, method=c("pearson"))

whorl_cuckoo <- cor.test(behav_cuckoo$n_flowers_per_whorl, behav_cuckoo$timevisit_whorl, method=c("pearson"))
plant_cuckoo <- cor.test(behav_cuckoo$n_flowers_total, behav_cuckoo$timevisit_plant, method=c("pearson"))

whorl_trueM <- cor.test(behav_trueM$n_flowers_per_whorl, behav_trueM$timevisit_whorl, method=c("pearson"))
plant_trueM <- cor.test(behav_trueM$n_flowers_total, behav_trueM$timevisit_plant, method=c("pearson"))

whorl_trueW <- cor.test(behav_trueW$n_flowers_per_whorl, behav_trueW$timevisit_whorl, method=c("pearson"))
plant_trueW <- cor.test(behav_trueW$n_flowers_total, behav_trueW$timevisit_plant, method=c("pearson"))


##plots
behav$log_whorl <- log(behav$timevisit_whorl)#add log-transformed time visit per whorl
behav$log_plant <- log(behav$timevisit_plant)#add log-transformed time visit per plant

w <- ggscatter(behav, x = "n_flowers_per_whorl", y = "log_whorl", 
               add = "reg.line", conf.int = TRUE, 
               color = "sub_sex", palette = c("#666666", "#1B9E77", "#D95F02"), shape = "sub_sex",
               size = 3, alpha = 0.6,
               xlab = "N. of flowers visited", ylab = "log2(time spent within whorl)"); w

p <- ggscatter(behav, x = "n_flowers_per_whorl", y = "log_plant", 
               add = "reg.line", conf.int = TRUE, 
               color = "sub_sex", palette = c("#666666", "#1B9E77", "#D95F02"), shape = "sub_sex",
               size = 3, alpha = 0.6,
               xlab = "N. of flowers visited", ylab = "log2(time spent on plant)"); p

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

#----------#
## MODELS ##
#----------#

#Linear mixed-effects models (LMM) with data log-transformed (gaussian distribution, link = identity)
#Error term coded as 12 data points each consisting of an interval within day within year
#The intercept was excluded to obtain the precise values for each variable (i.e. bumblebee type)


##Time spent within each whorl
tver <- lmer(log(timevisit_whorl) ~ -1 + sub_sex + (1 | int_tot), behav, REML=F, na.action = na.pass)
anova(tver)
summary(tver)
plot(density(resid(tver))) #plot density of residuals

emmeans(tver, pairwise~sub_sex) #pairwise comparisons


##Time spent on plants - without accounting for number of whorls
behav_singleID <- behav %>% distinct(insect_ID, .keep_all = TRUE) #exclude multiple IDs (for plant level analyses)

tplant <- lmer(log(timevisit_plant) ~ -1 + sub_sex + (1 | int_tot), behav_singleID, REML=F)
anova(tplant)
summary(tplant)
plot(density(resid(tplant))) #plot density of residuals


##Time spent on plants accounting for number of whorls - without interaction
##This is the best model after model selection (see below)
tplver_noint <- lmer(log(timevisit_plant) ~ -1 + sub_sex + tot_whorls + (1 | int_tot), behav_singleID, REML=F)
anova(tplver_noint)
summary(tplver_noint)
plot(density(resid(tplver_noint))) #plot density of residuals

emmeans(tplver_noint, pairwise~sub_sex) #pairwise comparisons


##Time spent on plants accounting for number of whorls - with interaction between bumblebee type and whorls
tplver <- lmer(log(timevisit_plant) ~ -1 + sub_sex + tot_whorls + sub_sex:tot_whorls + (1 | int_tot), behav_singleID, REML=F)
anova(tplver)
summary(tplver)
plot(density(resid(tplver))) #plot density of residuals


##Model selection table based on AICc values
library(MuMIn)
model.sel(tplant, tplver, tplver_noint)



##Number of flowers visited within whorl
flwhorl <- lmer(log(n_flowers_per_whorl) ~ -1 + sub_sex + (1 | int_tot), behav, REML=F)
anova(flwhorl)
summary(flwhorl)
plot(density(resid(flwhorl))) #plot density of residuals

emmeans(flwhorl, pairwise~sub_sex)

##Number of flowers visited per plant
flplant <- lmer(log(n_flowers_total) ~ -1 + sub_sex + (1 | int_tot), behav_singleID, REML=F)
anova(flplant)
summary(flplant)
plot(density(resid(flplant))) #plot density of residuals

emmeans(flplant, pairwise~sub_sex)

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

#----------------#
## VIOLIN PLOTS ##
#----------------#

#Significance symbols above the plots were changed from those automatically generated 
#to be coherent with the results obtained from the LMMs

behav$log_timevisit_whorl <- log(behav$timevisit_whorl) #add log_visit column
behav_singleID$log_timevisit_plant <- log(behav_singleID$timevisit_plant) #add log_visit column

my_comparisons_sex <- list( c("mC", "mT"), c("mT", "wT"), c("mC", "wT") ) #add significance levels comparing groups

##Time spent within each whorl
ggviolin(behav, x = "sub_sex", y = "log_timevisit_whorl", fill = "sub_sex",
         palette = c("#666666", "#1B9E77", "#D95F02"), alpha = .5,
         add = "boxplot", add.params = list(fill = "white")) +
  labs(x = "\nBumblebee type", y= "log2(time spent within whorl)\n") +
  scale_x_discrete(labels=c("mC" = "Cuckoo males", "mT" = "True males", "wT" = "True workers")) +
  theme(legend.position = "none") +
  stat_compare_means(comparisons = my_comparisons_sex, label = "p.signif")


##Time spent on plant
ggviolin(behav_singleID, x = "sub_sex", y = "log_timevisit_plant", fill = "sub_sex",
         palette = c("#666666", "#1B9E77", "#D95F02"), alpha = .5,
         add = "boxplot", add.params = list(fill = "white")) +
  labs(x = "\nBumblebee type", y= "log2(time spent on plants)\n") +
  scale_x_discrete(labels=c("mC" = "Cuckoo males", "mT" = "True males", "wT" = "True workers")) +
  theme(legend.position = "none") +
  stat_compare_means(comparisons = my_comparisons_sex, label = "p.signif")


##Number of flowers visited within each whorl
ggviolin(behav, x = "sub_sex", y = "n_flowers_per_whorl", fill = "sub_sex",
         palette = c("#666666", "#1B9E77", "#D95F02"), alpha = .5,
         add = "boxplot", add.params = list(fill = "white")) +
  labs(x = "\nBumblebee type", y= "N. of flowers visited per whorl\n") +
  scale_x_discrete(labels=c("mC" = "Cuckoo males", "mT" = "True males", "wT" = "True workers")) +
  theme(legend.position = "none") +
  stat_compare_means(comparisons = my_comparisons_sex, label = "p.signif")


##Number of flowers visited on plant
ggviolin(behav_singleID, x = "sub_sex", y = "n_flowers_total", fill = "sub_sex",
         palette = c("#666666", "#1B9E77", "#D95F02"), alpha = .5,
         add = "boxplot", add.params = list(fill = "white")) +
  labs(x = "\nBumblebee type", y= "N. of flowers visited per plant\n") +
  scale_x_discrete(labels=c("mC" = "Cuckoo males", "mT" = "True males", "wT" = "True workers")) +
  theme(legend.position = "none") +
  stat_compare_means(comparisons = my_comparisons_sex, label = "p.signif")