#Poka results
library(ggplot2)
setwd


# F0 days to 50% emergence ------------------------------------------------

#F0 median time to 50% emergence

s<- read.csv("F0 emergence 2.csv")
head(s)

F0_50E_graph <- ggplot(s, aes(Treatment, Days_since_beginning, fill=Sex)) +
  geom_boxplot(notch = FALSE) +
  labs(x = 'Dung type', y = 'Days to 50% emergence') +
  theme_bw() +
  scale_fill_discrete(labels = c("Female", "Male")) +
  theme(
    legend.position = 'right',
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 10)
  )
F0_50E_graph


g<-read.csv("F0 emergence 3.csv")
#Poisson regression females
(p3 <- glm(F_Days_since_beginning ~ Treatment, family="poisson", data=g))
summary(p3)
#RD greater than DF -> use quasipoisson
(p3 <- glm(F_Days_since_beginning ~ Treatment, family="quasipoisson", data=g))
summary(p3)

anova(p3, test = "Chisq")

#males
(p4 <- glm(M_Days_since_beginning ~ Treatment, family="poisson", data=g))
summary(p4)
#RD greater than DF -> use quasipoisson
(p4 <- glm(M_Days_since_beginning ~ Treatment, family="quasipoisson", data=g))
summary(p4)

anova(p4, test = "Chisq")
# F1 total emergence ------------------------------------------------------

df<-read.csv("F1 total emergence 2.csv")
head(df)

#F1 Total offspring emergence
# The code below will generate the graph:
F1_TE_graph <- ggplot(df, aes(x = Species, y = Total)) +
  geom_point(aes(colour=Species), size=3,stat = 'summary', fun.y = 'mean') +
  geom_errorbar(aes(colour=Species), stat = 'summary', fun.data = 'mean_se', 
                width=0.1, fun.args = list(mult = 1)) +
  labs(x = 'Dung type', y = 'Mean offspring')+
  labs(colour = 'Legend')+ 
  labs(colour = "Treatment") +
  theme_bw() +
  theme(
    legend.position = 'none',
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11)
  )
F1_TE_graph

#Poisson regression
(p1 <- glm(Total ~ Species, family="poisson", data=df))
summary(p1)
#RD greater than 1.5 x DF -> use quasipoisson
(p1 <- glm(Total ~ Species, family="quasipoisson", data=df))
summary(p1)

a3 <- anova(p1, test = "Chisq")
summary(a3)


# F1 day of emergence -----------------------------------------------------


k<- read.csv("F1 day of emergence 3.csv")
head(k)

F1_DOE_graph <- ggplot(k, aes(Treatment, Median_DTE, fill=Treatment)) +
  geom_boxplot(notch = FALSE) +
  labs(x = 'Dung type', y = 'Days to offspring emergence') +
  theme_bw() +
  theme(
    legend.position = 'none',
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11)
  )
F1_DOE_graph

#Poisson regression
(p2 <- glm(Median_DTE ~ Treatment, family="poisson", data=k))
summary(p2)

anova(p2, test = "Chisq")





# F1 pronotal width -------------------------------------------------------

w<- read.csv("F1 fitness 2.csv")
head(w)

#pronotal width x treatment + sex
# The code below will generate the graph:
F1_P_graph <- ggplot(w, aes(Sex, Mean_pronotum, fill = Species)) +
  geom_boxplot(notch = FALSE) +
  labs(x = 'Sex', y = 'Pronotal width (mm)') +
  scale_x_discrete(labels=c("F" = "Female", "Ma" = "Major male",
                            "Mi" = "Minor male")) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 11)
  )
F1_P_graph

#pronotal width comparison for each sex
i<- read.csv("Pronotal_mar31.csv")
head(i)

#Female
a_F<- aov(F_pronotum ~ Species, data=i)
summary(a_F)

#Major male
a_Ma<- aov(Ma_pronotum ~ Species, data=i)
summary(a_Ma)

#Minor male
a_Mi<- aov(Mi_pronotum ~ Species, data=i)
summary(a_Mi)
# F1 offspring count and male ratio ---------------------------------------------

library(cowplot)

#proportion females to males
t<- read.csv("ratio_fm.csv")
head(t)

F1_MC_graph3 <- ggplot(t, aes(y=Percent_female, x=Species, fill=Species)) + 
  geom_boxplot(notch = FALSE) +
  labs(x = 'Dung type', y = 'Proportion of females to total offspring') +
  theme_bw() +
  theme(
    legend.position = 'none',
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background=element_rect(colour="white",
                                  fill="white")
  ) 

F1_MC_graph3


#major/total male ratio
o<- read.csv("F1 male ratio.csv")
head(o)

F1_MC_graph5 <- ggplot(o, aes(y=Count, x=Species, fill=Species)) + 
  geom_boxplot(notch = FALSE) +
  labs(x = 'Dung type', y = 'Proportion of major males to total males') +
  theme_bw() +
  theme(
    legend.position = 'none',
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) 

F1_MC_graph5

#merge graphs
dual_plot<-plot_grid(F1_MC_graph3+theme(legend.position="none"),
                     F1_MC_graph5+theme(legend.position="none"), labels = c('A', 'B'), label_size = 12)
              
dual_plot

#analysis for proportion males to total offspring
h<- read.csv("F1 proportion.csv")
head(h)
a3 <- aov(Percent_male ~ Species, data = h)
summary(a3)

#analysis for proportion major to total males
a2 <- aov(Count ~ Species, data = o)
summary(a2)

tukey.test.a2 <- TukeyHSD(a2)
tukey.test.a2