---
title: "Bee Buzzes from Audacity Software"
author: "Carlos Coquinho"
date: "04/05/2020"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Peak frequencies of bee buzzes. Audacity data
```{r loading packages, results='hide'}
#install.packages("randomcoloR")
library(dplyr)
library(ggplot2)
library(randomcoloR)
```

```{r loading data}
bdat<-read.csv(file="ColonyA1cit_to_het.csv")
head(bdat)
```

In this dataframe:
variable "Bee" refers to individual bees used in the experiment; "Flower" refers to the plant taxon;
"Visit" refers to the a visit of an individual bee to a flower in which buzz-pollination occurs;
"Buzz" is the individual buzz (visually determined during analysis, see Methods);
"Frequency" is the peak frequency in Hz of each individual buzz;
"Duration" is the duration in miliseconds of each individual buzz.

Thus, each visit from a bee to a flower may comprise many buzzes with different frequencies and durations.

## Summary stats for peak frequencies of bees at Solanum citrullifolium and S. heterodoxum flowers
```{r summary stats 1}
std <- function(x) sd(x)/sqrt(length(x)) #function for std error

##number of bees for each visit (1st and 6th to S. cit)
##(7th, 12th and 13th to S. het)

visit.1 <- subset(bdat, Visit==1)
visit.1$Bee <- as.factor(visit.1$Bee)
nbees.1 <- nlevels(visit.1$Bee)

visit.6 <- subset(bdat, Visit==6)
visit.6$Bee <- as.factor(visit.6$Bee)
nbees.6 <- nlevels(visit.6$Bee)

visit.7 <- subset(bdat, Visit==7)
visit.7$Bee <- as.factor(visit.7$Bee)
nbees.7 <- nlevels(visit.7$Bee)

visit.12 <- subset(bdat, Visit==12)
visit.12$Bee <- as.factor(visit.12$Bee)
nbees.12 <- nlevels(visit.12$Bee)

visit.13 <- subset(bdat, Visit==13)
visit.13$Bee <- as.factor(visit.13$Bee)
nbees.13 <- nlevels(visit.13$Bee)
```

```{r summary stats 2}
bee_summary_stats <- bdat %>% group_by(Flower:factor(Visit)) %>%
   summarise(mean = mean(Frequency),
            se = std(Frequency),
            # median = median(Frequency),
            # iqr = IQR(Frequency),
            n = n())
bee_summary_stats
```

```{r summary stats 3}
species_summary_stats <- bdat %>% group_by(Flower) %>% 
   summarise(mean = mean(Frequency),
            se = std(Frequency),
            # median = median(Frequency),
            # iqr = IQR(Frequency),
            n = n())
species_summary_stats
```

```{r summary stats 4}
overall_bee_freq <- bdat %>% 
   summarise(mean = mean(Frequency),
            se = std(Frequency),
            # median = median(Frequency),
            # iqr = IQR(Frequency),
            n = n())
overall_bee_freq
```

```{r summary stats 5}
nbuzzes_bee <- bdat %>%
   group_by(factor(Bee):Flower:factor(Visit)) %>% 
   summarise(nbuzzes = max(Buzz))
nbuzzes_bee
```

```{r summary stats 6}
range_nbuzzes <- range(nbuzzes_bee$nbuzzes)

nbuzzes_cit <- bdat %>% subset(Flower=="Cit") %>% 
   group_by(factor(Bee)) %>% 
   summarise(mean_buzzes = mean(Buzz),
             se = std(Buzz),
             nbuzzes = n())

summary_nbuzz_cit <- nbuzzes_cit %>% 
   summarise(mean_nbuzzes = mean(nbuzzes),
             se = std(nbuzzes))
summary_nbuzz_cit
```

```{r summary stats 7}
nbuzzes_het <- bdat %>% subset(Flower=="Het") %>% 
   group_by(factor(Bee)) %>% 
   summarise(mean_buzzes = mean(Buzz),
             se = std(Buzz),
             nbuzzes = n())

summary_nbuzz_het <- nbuzzes_het %>% 
   summarise(mean_nbuzzes = mean(nbuzzes),
             se = std(nbuzzes))
summary_nbuzz_het
```


## Exploratory plot
```{r plot 1}
#mycol<-viridis(3)[c(1,2)]
n <- nlevels(factor(bdat$Bee))
my_palette <- unname(distinctColorPalette(n))

#Plot
p1<-ggplot(bdat, aes(y=Frequency, x=Visit))+
      theme_bw()+
      geom_jitter(aes(colour = factor(Bee), shape = Flower),
                  show.legend = T, alpha = .75, 
                  size = 1.5, width = .3)+
      scale_colour_brewer("Colors in Paired",
                          palette="Paired")+
      guides(colour=F)
      #scale_color_manual(values=my_palette)
p1

#+geom_point()+ scale_colour_manual(values=mycol)
```

## LME Models

```{r load packages 2, results='hide'}
library(lme4)
library(lmerTest)
```

```{r model 1}
lme2<-lmer(Frequency ~ Flower + Visit + (1|Bee), data=bdat)
summary(lme2)
anova(lme2)
```

```{r model 2}
lme2a<-lmer(Frequency ~ Visit + (1|Bee), data=subset(bdat, Flower=="Cit"))
summary(lme2a)
anova(lme2a)
```

```{r model 3}
lme2b<-lmer(Frequency ~ Visit + (1|Bee), data=subset(bdat, Flower=="Het"))
summary(lme2b)
anova(lme2b)
```