---
title: "Queen_Data"
author: "Hannah Whitehead"
date: "10/18/2021"
output:
  word_document:
    toc: yes
    toc_depth: '3'
  html_document:
    toc: yes
    toc_depth: 3
  pdf_document:
    toc: yes
    toc_depth: '3'
---

#SET-UP
##Libraries
```{r Libraries,message=FALSE,warning=FALSE}
library(reshape)  
library(lattice)  
library(ggplot2)  
library(data.table) 
library(vegan)
library(MASS)
library(multcomp)
library(languageR)
library(nlme)
library(effects)
library(dplyr)
library(ggplot2)
library(cowplot)
library(car)
library(lme4)
library(lmerTest)
library(effects)
library(ResourceSelection)
library(corrplot)
library(janitor)
library(readr)
library(plyr)
library(tidyr)
library(rcompanion)
library(foreign)
library(RDocumentation)
library(betareg)
```

##Download Files
```{r Download,message=FALSE,warning=FALSE}
setwd("//cns-cafeadm.campus.ads.umass.edu/Common/VEGIPM/Hannah/7. Honey bee projects/SARE Queen Project/Queen Report/1. Final Docs")
FieldData <- read_csv("1. Queen_Field_Data.csv") #attach data that we measured in the field
LabData <- read_csv("2. Queen_Lab_Data.csv") #attach results from Tarpy Lab analysis
Queenright <- FieldData[which(FieldData$Pop_Group == "1"), ] #a sub-set of the field data, which includes only queenright colonies with complete entries
```

#Sample Sizes
```{r Summary Field Data,message=FALSE,warning=FALSE}
#FULL FIELD DATA
#Turn variables into factors so that I can look at sample sizes:
FieldData$Treatment <- as.factor(FieldData$Treatment)
FieldData$Round <- as.factor(FieldData$Round)
FieldData$Q_Status <- as.factor(FieldData$Q_Status)

summary (FieldData) #this allows me to take a look at the data overall
tabyl(FieldData, Round, Treatment) #this allows me to look at the # of queens reared per treatment and round

#QUEENRIGHT DATA 
#Turn variables into factors so that I can look at sample sizes:
Queenright$Treatment <- as.factor(Queenright$Treatment)
Queenright$Round <- as.factor(Queenright$Round)
Queenright$Worker_Pop <- as.factor(Queenright$Worker_Pop)
Queenright$Scorer <- as.factor(Queenright$Scorer)
Queenright$Q_Status <- as.factor(Queenright$Q_Status)

summary(Queenright) #this allows me to take a look at the data overall
tabyl(Queenright, Round, Treatment) #this allows me to look at the # of queens reared per treatment and round

#TARPY LAB DATA:
#Turn variables into factors so that I can look at sample sizes:
LabData$Treatment <- as.factor(LabData$Treatment)
LabData$Round <- as.factor(LabData$Round)

summary(LabData) #this allows me to take a look at the data overall
tabyl(LabData, Round, Treatment) #this allows me to look at the # of queens reared per treatment and round

```

#QUEEN SURVIVAL
##Overall numbers and figures
```{r Queen Survival Overall Numbers}
#need to re-format variables for analysis 
FieldData$Q_Status <- as.character(FieldData$Q_Status) 
FieldData$Q_Status <- as.numeric(FieldData$Q_Status)

#Looking at % queen survival by treatment and round
aggregate(Q_Status ~ Treatment, data = FieldData, mean) # % survival by treatment
aggregate(Q_Status ~ Round, data = FieldData, mean) # % survival by round

#percent surviving queens in each round AND treatment
Q_Tot_Perc <- aggregate(Q_Status ~ Round + Treatment, data = FieldData, mean) 
spread(Q_Tot_Perc, Round, Q_Status)

#number of surviving queens in each round and treatment
Q_Tot_Num <- aggregate(Q_Status ~ Treatment + Round, data = FieldData, sum) 
spread(Q_Tot_Num, Round, Q_Status)

#Graph queen survival
QTake_Mean <- aggregate(Q_Status ~ Treatment, data = FieldData, mean)
QTake_Mean$fail <- 1-QTake_Mean$Q_Status
QTake_Mean
Treatments <- c("10-Day", "48-Hour", "Walk-Away")
P_QTake<-ggplot(data=QTake_Mean, aes(x=Treatment, y=Q_Status, fill=Treatment)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Treatment",values=c("firebrick3", "steelblue3", "seagreen3"), 
                    labels=c("10-Day", "48-Hour", "Walk-Away")) +
  scale_y_continuous(limits = c(0,1), labels = scales::percent) +
  scale_x_discrete(labels= Treatments) +
  ylab("Percent Queen Take") +
  xlab("Treatment") 
P_QTake
```

##Effect of treatment on queen survival: not significant
```{r Effect of treatment on survival}
#Treatment*Round interaction
M_QSurvival<-glm(Q_Status~Treatment*Round, data=FieldData, family=binomial) 
Anova(M_QSurvival) #interaction not significant (p=0.4029)

#Treatment + round, dropped the interaction
M_QueenSurvival2<-glm(Q_Status~Treatment + Round, data=FieldData, family=binomial) 
Anova(M_QueenSurvival2) #not significant (Treatment, p=0.9448 -- Round, p=0.5262)
summary(M_QueenSurvival2) #gives you the intercepts
plot(allEffects(M_QueenSurvival2))

```

##Compare 3 WALK-AWAY GROUPS - no significant effect of treatment
```{r Comparing walk away groups}
#Queen survival
QTake_WA <- FieldData[which(FieldData$WA_Group != "0"), ] #only WA splits with a group
QTake_WA$WA_Group <- as.factor(QTake_WA$WA_Group) #format
tabyl(QTake_WA, Round, WA_Group) #looking at sample sizes

#looked for significance:
M_WA <- glm(Q_Status ~ WA_Group*Round, data=QTake_WA, family=binomial) #interaction not significant
Anova(M_WA)  

M_WA2 <- glm(Q_Status ~ WA_Group + Round, data=QTake_WA, family=binomial) # WA Group not significant
Anova(M_WA2)
plot(allEffects(M_WA2))

M_WA3 <- glm(Q_Status ~ WA_Group, data=QTake_WA, family=binomial) # WA Group not significant
Anova(M_WA3)
```


#NUC POPULATION
##Overall numbers and figures
```{r nuc population overall numbers}
#need to re-format variables for analysis. Only using queenright colonies for this analysis.
Queenright$Worker_Pop<- as.character(Queenright$Worker_Pop)  
Queenright$Worker_Pop <- as.numeric(Queenright$Worker_Pop)

#Looking at worker pop numbers by treatment and round:
aggregate(Worker_Pop ~ Treatment, data = Queenright, mean) # Avg. pop size by treatment
aggregate(Worker_Pop ~ Round, data = Queenright, mean) # Avg. pop size by round

# Avg. Size by treatment AND round:
Pop_Tot_Num <- aggregate(Worker_Pop ~ Round + Treatment, data = Queenright, mean)
spread(Pop_Tot_Num, Round, Worker_Pop)
summary(Queenright)

# looking at the data agin
# --- by treatment
PSize_Mean_se <- ddply(Queenright, c("Treatment"), summarise,
               N    = length(Worker_Pop),
               mean = mean(Worker_Pop),
               sd   = sd(Worker_Pop),
               se   = sd / sqrt(N))
PSize_Mean_se

# --- by round
PSize_Mean_se2 <- ddply(Queenright, c("Round"), summarise,
               N    = length(Worker_Pop),
               mean = mean(Worker_Pop),
               sd   = sd(Worker_Pop),
               se   = sd / sqrt(N))
PSize_Mean_se2

#Bar Plot by Treatment
#--make sure the data is properly formatted:
Queenright$Worker_Pop <- as.character(Queenright$Worker_Pop)
Queenright$Worker_Pop <- as.numeric(Queenright$Worker_Pop)

Treatment <- c("10-Day", "48-Hour", "Walk-Away")
P_PSize<-ggplot(data=PSize_Mean_se, aes(x=Treatment, y=mean, fill=Treatment)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Treatment",values=c("firebrick3", "steelblue3", "seagreen3"), 
                    labels=c("10-Day", "48-Hour", "Walk-Away")) +
  scale_y_continuous(limits = c(0,3.5)) +
  scale_x_discrete(labels= Treatments) +
  ylab("Average Population Size Rating") +
  xlab("Treatment") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1,
                 position=position_dodge(.9)) 
P_PSize

#Bar Plot by Round
Round <- c("Round 1", "Round 2")
P_PSize2<-ggplot(data=PSize_Mean_se2, aes(x=Round, y=mean, fill=Round)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Round",values=c("goldenrod2", "goldenrod4"), 
                    labels=c("Round 1", "Round 2")) +
  scale_y_continuous(limits = c(0,3.5)) +
  scale_x_discrete(labels= Round) +
  ylab("Average Population Size Rating") +
  xlab("Round") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1,
                 position=position_dodge(.9)) 
P_PSize2
```

##Effect of treatment on population - overall
```{r Effect of treatment on worker pop}
#used ordinal logistic regression since pop size has three distinct bins, that have a specific relationship to each other
Queenright$Worker_Pop <- as.factor(Queenright$Worker_Pop)
M_Pop <- polr(Worker_Pop~Treatment*Round, data=Queenright, Hess = TRUE) #NOTE: Worker_Pop MUST be a factor
Anova(M_Pop) #Interaction was not significant (p = 0.86801)

M_Pop <- polr(Worker_Pop~Treatment + Round, data=Queenright, Hess = TRUE) #tried without the interaction term
Anova(M_Pop) #Treatment not significant (p=0.07049), Round significant (p=0.04722)
summary(M_Pop)
plot(allEffects(M_Pop))

```

##Effect of scorer on population size ratings
```{r}
#How many nucs were rated per person in each round?
Queenright$Worker_Pop <- as.factor(Queenright$Worker_Pop)
Queenright$Scorer <- as.factor(Queenright$Scorer)
tabyl(Queenright, Round, Scorer) #sample sizes

#what was the average nuc population size rating per person?
Queenright$Worker_Pop <- as.character(Queenright$Worker_Pop)
Queenright$Worker_Pop <- as.numeric(Queenright$Worker_Pop)
PSize_Mean_se3 <- ddply(Queenright, c("Scorer"), summarise,
               N    = length(Worker_Pop),
               mean = mean(Worker_Pop),
               sd   = sd(Worker_Pop),
               se   = sd / sqrt(N))
PSize_Mean_se3 #average values and se

#what was the average nuc population size rating by person AND round?
Queenright$Worker_Pop<- as.character(Queenright$Worker_Pop)  
Queenright$Worker_Pop <- as.numeric(Queenright$Worker_Pop)
Pop_Tot_Num_S <- aggregate(Worker_Pop ~ Round + Scorer, data = Queenright, mean)
spread(Pop_Tot_Num_S, Scorer, Worker_Pop)

#did scorer significantly affect worker population size rating?
Queenright$Worker_Pop <- as.factor(Queenright$Worker_Pop)
M_Scorer <- polr(Worker_Pop~Scorer, data=Queenright, Hess = TRUE) 
Anova(M_Scorer) #yes, there was a significant effect
summary(M_Scorer)
plot(allEffects(M_Scorer))

#if we include scorer instead of round as a predictor, do we see an effect of treatment on pop size?
Queenright$Worker_Pop <- as.factor(Queenright$Worker_Pop)
M_Scorer2 <- polr(Worker_Pop~Treatment + Scorer, data=Queenright, Hess = TRUE) 
Anova(M_Scorer2) #no significant effect of treatment on pop size
summary(M_Scorer2)
plot(allEffects(M_Scorer2))

#bar plot to look at data visually
Scorer <- c("AR", "BK", "SC")
P_PSize3<-ggplot(data=PSize_Mean_se3, aes(x=Scorer, y=mean, fill=Scorer)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Scorer",values=c("firebrick3", "steelblue3", "seagreen3"), 
                    labels=c("AR", "BK", "SC")) +
  scale_y_continuous(limits = c(0,3.5)) +
  scale_x_discrete(labels= Scorer) +
  ylab("Average Population Size Rating") +
  xlab("Scorer") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1,
                 position=position_dodge(.9)) 
P_PSize3
```


##Effect of walk-away group on population
```{r effect of walk-away groups on population}
#How many nucs of each method per round?
Queenright2<-QTake_WA[which(QTake_WA$Pop_Group != "0"), ] 
tabyl(Queenright2, Round, WA_Group)

#Worker population
M_WApop <- polr(as.factor(Worker_Pop) ~ WA_Group*Round, data=Queenright2, Hess = TRUE) 
Anova(M_WApop) ## WA method not significant

M_WApop2 <- polr(as.factor(Worker_Pop) ~ WA_Group + Round, data=Queenright2, Hess = TRUE) 
Anova(M_WApop2) ## WA method not significant 

M_WApop3 <- polr(as.factor(Worker_Pop) ~ WA_Group, data=Queenright2, Hess = TRUE) 
Anova(M_WApop3) ## WA method not significant
```

#TARPY DATA
##-Checking sample sizes
```{r Tarpy overall numbers}
LabData$Round <- as.factor(LabData$Round)
LabData$Treatment <- as.factor(LabData$Treatment)
tabyl(LabData, Round, Treatment)
```

##-Effect of treatment and round on...
###---Weight
```{r weight}
#Format the data
LabData$Weight <- as.numeric(LabData$Weight)


#Average values by treatment
aggregate(Weight ~ Treatment, data = LabData, mean)

Weight_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Weight),
               mean = mean(Weight),
               sd   = sd(Weight),
               se   = sd / sqrt(N))
Weight_Mean_se

#Significant differences?
#W <- lm(Weight ~ Treatment*Round, data=LabData) #interaction not significant
W <- lm(Weight ~ Treatment + Round, data=LabData) 
Anova(W) 
plot(allEffects(W)) #plot data
summary(W)

#bar plot
Treatments <- c("10-Day", "48-Hour", "Walk-Away")
P_Weight<-ggplot(data=Weight_Mean_se, aes(x=Treatment, y=mean, fill=Treatment)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Treatment",values=c("firebrick3", "steelblue3", "seagreen3"), 
                    labels=c("10-Day", "48-Hour", "Walk-Away")) +
  scale_y_continuous(limits = c(0,225)) +
  scale_x_discrete(labels= Treatments) +
  ylab("Average Weight (mg)") +
  xlab("Treatment") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1,
                 position=position_dodge(.9)) 
P_Weight

```

###---Head Width
```{r Head width}
#Format the data
LabData$Head_width <- as.numeric(LabData$Head_width)

#Average values by treatment
aggregate(Head_width ~ Treatment, data = LabData, mean)

Head_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Head_width),
               mean = mean(Head_width),
               sd   = sd(Head_width),
               se   = sd / sqrt(N))
Head_Mean_se

#Significant differences?
#M_Head <- lm(Head_width ~ Treatment*Round, data=LabData) #interaction not significant
M_Head <- lm(Head_width ~ Treatment + Round, data=LabData) 
Anova(M_Head) 
plot(allEffects(M_Head)) #plot data
summary(M_Head)
```

###---Thorax Width
```{r Thorax width}
#Format the data
LabData$Thorax_width <- as.numeric(LabData$Thorax_width)

#Average values by treatment
aggregate(Thorax_width ~ Treatment, data = LabData, mean)

Thorax_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Thorax_width),
               mean = mean(Thorax_width),
               sd   = sd(Thorax_width),
               se   = sd / sqrt(N))
Thorax_Mean_se

#Significant differences?
#M_Thorax <- lm(Thorax_width ~ Treatment*Round, data=LabData) #interaction not significant
M_Thorax <- lm(Thorax_width ~ Treatment + Round, data=LabData) 
Anova(M_Thorax) 
plot(allEffects(M_Thorax)) #plot data
summary(M_Thorax)

```

###---Tot. Sperm Count
```{r Tot Sperm Count}
#Format the data
LabData$Tot_Sperm_Count <- as.numeric(LabData$Tot_Sperm_Count)
LabData$Tot_Sperm_Count_6 <- LabData$Tot_Sperm_Count/1000000 #change to 10^6

#Average values by treatment
aggregate(Tot_Sperm_Count ~ Treatment, data = LabData, mean)

SpermC_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Tot_Sperm_Count),
               mean = mean(Tot_Sperm_Count),
               sd   = sd(Tot_Sperm_Count),
               se   = sd / sqrt(N))
SpermC_Mean_se

#Average values by treatment - 10^6
aggregate(Tot_Sperm_Count_6 ~ Treatment, data = LabData, mean)

SpermC_Mean_se6 <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Tot_Sperm_Count_6),
               mean = mean(Tot_Sperm_Count_6),
               sd   = sd(Tot_Sperm_Count_6),
               se   = sd / sqrt(N))
SpermC_Mean_se6

#Significant differences?
#M_SpermC <- lm(Tot_Sperm_Count ~ Treatment*Round, data=LabData) #interaction not significant
M_SpermC <- lm(Tot_Sperm_Count ~ Treatment + Round, data=LabData) 
Anova(M_SpermC) 
plot(allEffects(M_SpermC)) #plot data
summary(M_SpermC)

#bar blot
Treatments <- c("10-Day", "48-Hour", "Walk-Away")
P_TotSperm<-ggplot(data=SpermC_Mean_se6, aes(x=Treatment, y=mean, fill=Treatment)) +
  geom_bar(stat="identity", color="black", width=0.7) + 
  theme_bw() +
  scale_fill_manual(name="Treatment",values=c("firebrick3", "steelblue3", "seagreen3"), 
                    labels=c("10-Day", "48-Hour", "Walk-Away")) +
  scale_y_continuous(limits = c(0,6.7)) +
  scale_x_discrete(labels= Treatments) +
  ylab("Total Sperm Count (x10^6)") +
  xlab("Treatment") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.1,
                 position=position_dodge(.9)) 
P_TotSperm

```

###---Viable Sperm Count
```{r Viable Sperm Count}
#Format the data
LabData$Viable_Sperm_Count <- as.numeric(LabData$Viable_Sperm_Count)
LabData$Viable_Sperm_Count_6 <- LabData$Viable_Sperm_Count/1000000  #change to 10^6

#Average values by treatment
aggregate(Viable_Sperm_Count ~ Treatment, data = LabData, mean)

Viable_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Viable_Sperm_Count),
               mean = mean(Viable_Sperm_Count),
               sd   = sd(Viable_Sperm_Count),
               se   = sd / sqrt(N))
Viable_Mean_se

#Average values by treatment - 10^6
aggregate(Viable_Sperm_Count_6 ~ Treatment, data = LabData, mean)

Viable_Mean_se6 <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Viable_Sperm_Count_6),
               mean = mean(Viable_Sperm_Count_6),
               sd   = sd(Viable_Sperm_Count_6),
               se   = sd / sqrt(N))
Viable_Mean_se6

#Significant differences?
#M_SpermV <- lm(Viable_Sperm_Count ~ Treatment*Round, data=LabData) #interaction not significant
M_SpermV <- lm(Viable_Sperm_Count ~ Treatment + Round, data=LabData) 
Anova(M_SpermV) 
plot(allEffects(M_SpermV)) #plot data
summary(M_SpermV)

```

###---Spermatheca Diameter
```{r Spermatheca Diameter}
#Format the data
LabData$Spermatheca_diameter <- as.numeric(LabData$Spermatheca_diameter)

#Average values by treatment
aggregate(Spermatheca_diameter ~ Treatment, data = LabData, mean)

Diameter_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Spermatheca_diameter),
               mean = mean(Spermatheca_diameter),
               sd   = sd(Spermatheca_diameter),
               se   = sd / sqrt(N))
Diameter_Mean_se

#Significant differences?
#M_Diameter <- lm(Spermatheca_diameter ~ Treatment*Round, data=LabData) #interaction not significant
M_Diameter <- lm(Spermatheca_diameter ~ Treatment + Round, data=LabData) 
Anova(M_Diameter) 
plot(allEffects(M_Diameter)) #plot data
summary(M_Diameter)

```

###---% Sperm Viability
```{r Sperm viability}
#Format the data
LabData$Perc_Viability <- as.numeric(LabData$Perc_Viability)
#LabData$Perc_Viability <- LabData$Perc_Viability*100

#Average values by treatment
aggregate(Perc_Viability ~ Treatment, data = LabData, mean)

Viability_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Perc_Viability),
               mean = mean(Perc_Viability),
               sd   = sd(Perc_Viability),
               se   = sd / sqrt(N))
Viability_Mean_se

#Significant differences?
#used beta regression since it's a continuous variable between 0 and 1 (e.g. percent, proportion, etc.)
#M_Viability <- betareg(Perc_Viability ~ Treatment*Round, data=LabData) #interaction not significant 
M_Viability <- betareg(Perc_Viability ~ Treatment + Round, data=LabData) 
Anova(M_Viability) 
plot(allEffects(M_Viability)) #plot data
plot(M_Viability)
summary(M_Viability)

```

###---Spermatheca % filled
```{r Spermatheca filled}
#Format the data
LabData$Spermatheca_Perc_Filled <- as.numeric(LabData$Spermatheca_Perc_Filled)

#Average values by treatment
aggregate(Spermatheca_Perc_Filled ~ Treatment, data = LabData, mean)

Filled_Mean_se <- ddply(LabData, c("Treatment"), summarise,
               N    = length(Spermatheca_Perc_Filled),
               mean = mean(Spermatheca_Perc_Filled),
               sd   = sd(Spermatheca_Perc_Filled),
               se   = sd / sqrt(N))
Filled_Mean_se

#Significant differences?
#M_Filled <- betareg(Spermatheca_Perc_Filled ~ Treatment*Round, data=LabData) #interaction not significant
M_Filled <- betareg(Spermatheca_Perc_Filled ~ Treatment + Round, data=LabData) 
Anova(M_Filled) 
plot(allEffects(M_Filled)) #plot data
summary(M_Filled)

```