#packages
library(DHARMa)
library(extrafont)#fonts
library(ggpubr)
library(ggplot2)
library(cowplot)
library(forcats)
library(emmeans)
library(car)
library(ggstats)
library(scales)
library(Rmisc)



# ------------------------------- #
#### DATA PHENOTYPE ####
data_phenotype<-read.csv("submissionprocb_data_phenotype.csv",header=T,sep=";",dec=",", na.string ="")
data_phenotype$mitotype <- fct_relevel(data_phenotype$mitotype, c("N", "D", "K"))
levels(data_phenotype$mitotype)


#### MALE  STATUS ####
x <- xtabs(~statut.male + mitotype, data = data_phenotype) 
x

#FIGURE 
pigmentation_pourcentage <- 
  ggplot(data_phenotype) +
  aes(x = mitotype, fill = statut.male, by = factor(mitotype), y = after_stat(prop),
      label = percent(after_stat(prop), accuracy = 1),) +
  labs(title = "", x="Mitotype", y="Percentage") +
  scale_fill_manual(labels = c("Fertile", "Sterile"), values = c("gray29", "gray"))+
  theme_test() +
  theme(legend.position = "top", legend.title = element_blank(), legend.text = element_text(size = 14, color="grey20")) +
  theme(axis.text.x = element_text(size=14), axis.title.x = element_text(size = 16, family = "Arial", color = "grey20"))+
  theme(axis.text.y = element_text(size=14), axis.title.y = element_text(size = 16, family = "Arial", color = "gray20"))+
  geom_bar(stat = "prop", position = "fill")  +
  scale_y_continuous(name="Percentage", labels = percent)
pigmentation_pourcentage

#STAT ----
pigmentation_glm <- glm(pigmentation.ponte.albinos~mitotype, data=data_phenotype, family=binomial(link=logit))
summary(pigmentation_glm)
Anova(pigmentation_glm)
emmeans(pigmentation_glm, list(pairwise ~ mitotype), adjust = "tukey")

#assumptions
testDispersion(pigmentation_glm)  
plot(simulateResiduals(fittedModel = pigmentation_glm, plot = F))




# ------------------------------- #
#### DATA PHENOTYPE STERILE/FERTILE ###
data_phenotype_sf<-read.csv("submissionprocb_data_phenotype_sf.csv",header=T,sep=";",dec=",", na.string ="")
data_phenotype_sf$mitotype <- fct_relevel(data_phenotype_sf$mitotype, c("N", "D", "K"))
levels(data_phenotype_sf$mitotype)


#### WHOLE BODY MASS ####
avg_pds_sf<- summarySE(data_phenotype_sf,measurevar="poids", na.rm=TRUE, groupvars=c("mitotype"))
avg_pds_sf

ggplot(data_phenotype_sf, aes(y=poids, x=mitotype, colour=mitotype, fill=mitotype))+
  geom_boxplot(alpha = 0.5)+
  geom_jitter(width=)+
  theme_classic() 

#STAT ----
mod_pds_sf <- lm(poids~mitotype,data=data_phenotype_sf) 
summary(mod_pds_sf)
Anova(mod_pds_sf)
emmeans(mod_pds_sf, list(pairwise ~ mitotype), adjust = "tukey")

#assumptions
par(mfrow=c(2,2))
plot(mod_pds_sf) 

ggqqplot(residuals(mod_pds_sf)) 
ggqqplot(data_phenotype_sf, "poids", facet.by = "mitotype")
shapiro.test(residuals(mod_pds_sf))#p-value>0.05 

#plot(mod_pds_sf,3)
leveneTest(residuals(mod_pds_sf)~data_phenotype_sf$mitotype, data = data_phenotype_sf) #p-value<0.05 

plot(simulateResiduals(fittedModel = mod_pds_sf, plot = F))

#FIGURE WHOLE BODY MASS----
#avg_pds_sf <- summarySE(data_phenotype_sf,measurevar="poids", na.rm=TRUE, groupvars=c("mitotype"))
avg_pds_sf

fig_poids_sf <- ggstripchart(
  data = data_phenotype_sf,
  x = "mitotype",
  y = "poids",
  color = "mitotype",
  fill = "mitotype",
  size = 3.5,
  alpha = 0.7,
  jitter = 0.4,
  xlab = "",
  ylab = "Whole body mass (g)",
  ggtheme = theme_test(),
  legend = "none",
  shape = "mitotype") +
  geom_point(data=avg_pds_sf, 
             aes (y=poids),
             size = 5, 
             color="black") +
  geom_errorbar(data=avg_pds_sf,
                aes(y=poids,
                    ymin=poids-se, ymax=poids+se),
                width=0.2,
                size = 0.9,
                color="black")+
  theme(axis.text.x = element_text(size=14), axis.title.x = element_text(size = 16, family = "Arial", color = "grey20"))+
  theme(axis.text.y = element_text(size=14), axis.title.y = element_text(size = 16, family = "Arial", color = "grey20"))+
  theme(legend.title = element_text(size=15), legend.text = element_text(size = 15))
fig_poids_sf


#### MAX LENGTH ####
avg_taille_sf <- summarySE(data_phenotype_sf,measurevar="taille", na.rm=TRUE, groupvars=c("mitotype"))
avg_taille_sf

ggplot(data_phenotype_sf, aes(y=taille, x=mitotype, colour=mitotype, fill=mitotype))+
  geom_boxplot(alpha = 0.5)+
  geom_jitter(width=)+
  theme_classic() 

#STAT ----
mod_taille_sf <- lm(taille~mitotype,data=data_phenotype_sf) 
summary(mod_taille_sf)
Anova(mod_taille_sf)
emmeans(mod_taille_sf, list(pairwise ~ mitotype), adjust = "tukey")

#assumptions
par(mfrow=c(2,2))
plot(mod_taille_sf) 

ggqqplot(residuals(mod_taille_sf)) 
ggqqplot(data_phenotype_sf, "taille", facet.by = "mitotype")
shapiro.test(residuals(mod_taille_sf))#p-value>0.05 

#plot(mod_taille_sf,3)
leveneTest(residuals(mod_taille_sf)~data_phenotype_sf$mitotype, data = data_phenotype_sf) #p-value<0.05 non-homogeneite des variances des residus

plot(simulateResiduals(fittedModel = mod_taille_sf, plot = F))

#FIGURE MAX LENGTH ----
#avg_taille_sf <- summarySE(data_phenotype_sf,measurevar="taille", na.rm=TRUE, groupvars=c("mitotype"))
avg_taille_sf

fig_taille_sf <- ggstripchart(
  data = data_phenotype,
  x = "mitotype",
  y = "taille",
  color = "mitotype",
  fill = "mitotype",
  size = 3.5,
  alpha = 0.7,
  jitter = 0.4,
  xlab = "",
  ylab = "Maximal shell length (mm)",
  ggtheme = theme_test(),
  legend = "none",
  shape = "mitotype") +
  geom_point(data=avg_taille_sf, 
             aes (y=taille),
             size = 5, 
             color="black") +
  geom_errorbar(data=avg_taille_sf,
                aes(y=taille,
                    ymin=taille-se, ymax=taille+se),
                width = 0.2,
                size = 0.9,
                color="black") +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16, family = "Arial", color = "grey20"))+
  theme(axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16, family = "Arial", color = "grey20"))+
  theme(legend.title = element_text(size=15), legend.text = element_text(size = 15))
fig_taille_sf



plot_grid(fig_poids_sf,fig_taille_sf,labels=c("A", "B"), ncol = 2, nrow = 1)