#R code for the article in Oecologia : Solitary foundation or colony fission in ants: an intraspecific study shows that worker presence and number increase colony foundation success
#Author: Basile Finand

#R Version: 4.2.2


setwd("C:/Users/bafinand/OneDrive - University of Helsinki/PhD/Article expérimentations/Donnees")

#The script has been done in french:
Sys.setlocale("LC_TIME", "French")


####Packages####

library(stringr)
library(survival)
library(survminer)
library(ggplot2)
library(lme4)
library(car)
library(emmeans)


#####Data formatting####

donnees=read.table(file = "Data.txt", h=T)

#Creating a variable representing the name of the colony

donnees=cbind(donnees, Couvain = (donnees$Oeufs+donnees$Larves+donnees$Nymphes))
nouv_ind=donnees$ouvrieres-donnees$Nombre_ouvrieres
donnees=cbind(donnees, nouv_ind=nouv_ind, croissance = nouv_ind+donnees$Couvain)
colonie=c()
for (i in 1:length(donnees$Reference)){
  coloniei= str_split(donnees$Reference, "A")[[i]][1]
  colonie=c(colonie,coloniei)}
donnees=cbind(donnees, colonie=colonie)


#Put the date as a date synchronized per week

donnees$Date_accouplement=as.Date(x = as.character(donnees$Date_accouplement), format = "%d/%m/%Y")
donnees$Date=as.Date(x = as.character(donnees$Date), format = "%d/%m/%Y")
Jours_hiv = donnees$Date - as.Date(x = as.character("06/04/2021"), format = "%d/%m/%Y")
date_hiv = c("05/04/2021","06/04/2021","07/04/2021","08/04/2021","09/04/2021")
date_hiv = as.Date(x = as.character(date_hiv), format = "%d/%m/%Y")
semaine_hiv=c()

for (i in 1:length(donnees$Reference)){
  for (j in 1:length(date_hiv)){
    if (as.numeric(date_hiv[j]-donnees$Date_accouplement[i])%%7 == 0){
      semaine_hiv=c(semaine_hiv,(as.numeric(donnees$Date[i])-as.numeric(date_hiv[j]))/7)
      break
    }
  }
}



Jours_semaine = weekdays(donnees$Date)
donnees=cbind(donnees, Jours_semaine=Jours_semaine)

donnees=cbind(donnees,Semaine_lundi=as.Date(c(rep(NA,length(donnees$Reference)))))
for (i in 1:length(donnees[,1])){
  if(is.na(donnees$Jours_semaine[i]) == FALSE & donnees$Jours_semaine[i] == "lundi"){donnees$Semaine_lundi[i] = donnees$Date[i]}
  if(is.na(donnees$Jours_semaine[i]) == FALSE & donnees$Jours_semaine[i] == "mardi"){donnees$Semaine_lundi[i] = donnees$Date[i]-1}
  if(is.na(donnees$Jours_semaine[i]) == FALSE & donnees$Jours_semaine[i] == "mercredi"){donnees$Semaine_lundi[i] = donnees$Date[i]-2}
  if(is.na(donnees$Jours_semaine[i]) == FALSE & donnees$Jours_semaine[i] == "jeudi"){donnees$Semaine_lundi[i] = donnees$Date[i]-3}
  if(is.na(donnees$Jours_semaine[i]) == FALSE & donnees$Jours_semaine[i] == "vendredi"){donnees$Semaine_lundi[i] = donnees$Date[i]-4}
  
}

donnees=cbind(donnees,Jours_hiv=Jours_hiv,Semaine_hiv=semaine_hiv)
donnees=droplevels(donnees)
donnees$Reference = as.factor(donnees$Reference)


#Function to summarize the data

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}

donnees_completes = donnees



###############################
#####Virgin queens analysis####
###############################


donnees = donnees_completes

nb_col0_vivante_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "0"])))
nb_col0_vivante_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "0"])))
nb_col0_morte_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "0"])))
nb_col0_morte_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "0"])))
nb_col2_vivante_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "2"])))
nb_col2_vivante_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "2"])))
nb_col2_morte_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "2"])))
nb_col2_morte_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "2"])))
nb_col4_vivante_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "4"])))
nb_col4_vivante_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "VIVANTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "4"])))
nb_col4_morte_sppleine = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "PLEINE" & donnees$Nombre_ouvrieres == "4"])))
nb_col4_morte_spvide = length(levels(droplevels(donnees$Reference[donnees$Statut == "MORTE" & donnees$Spermatheque == "VIDE" & donnees$Nombre_ouvrieres == "4"])))


M <- as.table(rbind(c(nb_col0_vivante_sppleine, nb_col2_vivante_sppleine,nb_col4_vivante_sppleine), c(nb_col0_vivante_spvide,nb_col2_vivante_spvide ,nb_col4_vivante_spvide))) # création d'une table 2 lignes/3 colonnes
dimnames(M) <- list(Spermatheque=c("Pleine","Vide"), Nb_ouv=c("0ouv","2ouv", "4ouv"))# entête colonne et ligne
(test <- chisq.test(M)) 


###########################
#####Survival analyzis#####
###########################


donnees_survie = read.table("donnees_survie.txt", h=T)
donnees_survie$Nombre_ouvrieres = as.factor(donnees_survie$Nombre_ouvrieres)
surv_diff <- survdiff(Surv(Temps, Statut) ~ Nombre_ouvrieres, data = donnees_survie)
surv_diff
(t=pairwise_survdiff(Surv(Temps, Statut) ~ Nombre_ouvrieres, data = donnees_survie, p.adjust.method = "BH"))

data_courbe_surv = data.frame(Nombre_ouvrieres=c(0,2,4), Temps=c(0,0,0), pourc_surv=c(1,1,1))
for(j in levels(donnees_survie$Nombre_ouvrieres)){ 
  nb_ind = length(donnees_survie$Reference[donnees_survie$Nombre_ouvrieres == j])
  for (i in 1:max(donnees_survie$Temps)){
    pourc_surv = (length(donnees_survie$Reference[donnees_survie$Nombre_ouvrieres == j]) - length(donnees_survie$Reference[donnees_survie$Nombre_ouvrieres == j & donnees_survie$Statut == 1 & donnees_survie$Temps <= i]))/nb_ind
    data_courbe_surv = rbind(data_courbe_surv,c(j,i,pourc_surv))
  }
}

data_courbe_surv$Nombre_ouvrieres = as.factor(data_courbe_surv$Nombre_ouvrieres)
data_courbe_surv$Temps = as.numeric(data_courbe_surv$Temps)
data_courbe_surv$pourc_surv = as.numeric(data_courbe_surv$pourc_surv)

#Figure 1

ggplot(data = data_courbe_surv, aes(x = Temps, y = pourc_surv, colour = Nombre_ouvrieres))+
  geom_vline(xintercept = 4, linetype="dotted")+
  geom_vline(xintercept = 11, linetype="dotted")+
  geom_rect(aes(xmin = 4, xmax = 11, ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_vline(xintercept = 23, linetype="dotted")+
  geom_vline(xintercept = 33, linetype="dotted")+
  geom_rect(aes(xmin = 23, xmax = 33, ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = 11, xmax = 23, ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_line(size=2)+
  scale_color_manual(values=c("#009E73","#E69F00","#D55E00"))+
  xlab("Number of weeks after matting")+
  ylab("Percentage of survival")+
  theme_classic()+ 
  theme(legend.position="bottom")

#########################
#####Growth analyzis#####
#########################


donnees=donnees_completes[donnees_completes$Statut == "VIVANTE",]
donnees=droplevels(donnees)
(nb_col_vide0=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==0 & donnees$Spermatheque=="VIDE"]))))
(nb_col_vide2=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==2 & donnees$Spermatheque=="VIDE"]))))
(nb_col_vide4=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==4 & donnees$Spermatheque=="VIDE"]))))
(nb_col_pleine0=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==0 & donnees$Spermatheque=="PLEINE"]))))
(nb_col_pleine2=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==2 & donnees$Spermatheque=="PLEINE"]))))
(nb_col_pleine4=length(levels(droplevels(donnees$Reference[donnees$Nombre_ouvrieres==4 & donnees$Spermatheque=="PLEINE"]))))

M <- as.table(rbind(c(nb_col_vide0,nb_col_vide2, nb_col_vide4), c(nb_col_pleine0,nb_col_pleine2,nb_col_pleine4))) # création d'une table 2 lignes/3 colonnes
dimnames(M) <- list(Statut=c("Vide","Pleine"), Condition=c("0ouv","2ouv", "4ouv"))# entête colonne et ligne
(test <- chisq.test(M, correct=FALSE))

donnees=donnees[donnees$Spermatheque == "PLEINE",]
donnees=droplevels(donnees)

(nb_col_0=length(donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))][donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))]=="0"])) # Nombre de colonie 0
(nb_col_2=length(donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))][donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))]=="2"]))  # Nombre de colonie 2
(nb_col_4=length(donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))][donnees$Nombre_ouvrieres[1:length(levels(donnees$Reference))]=="4"]))  # Nombre de colonie 4


#Figure 3

donneessummo=data_summary(donnees, varname="Oeufs", groupnames=c("Semaine_lundi","Nombre_ouvrieres"))
donneessumml=data_summary(donnees, varname="Larves", groupnames=c("Semaine_lundi","Nombre_ouvrieres"))
donneessummn=data_summary(donnees, varname="Nymphes", groupnames=c("Semaine_lundi","Nombre_ouvrieres"))
donneessummi=data_summary(donnees, varname="nouv_ind", groupnames=c("Semaine_lundi","Nombre_ouvrieres"))

donneessummtot=cbind(donneessummo[,1:3],Larves = donneessumml$Larves,Nymphes = donneessummn$Nymphes,Ouvrieres = donneessummi$nouv_ind)

colors = rev(c("skyblue","olivedrab2","gold1","slateblue3"))

graph1=ggplot(donneessummtot[donneessummtot$Nombre_ouvrieres==0,], aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres))+
  geom_rect(aes(xmin = as.Date("2020-11-05"), xmax = as.Date("2020-11-26"), ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = as.Date("2020-11-26"), xmax = as.Date("2021-04-06"), ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_vline(xintercept = as.Date("2020-11-26"), linetype="dotted")+
  geom_area( aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres),fill=colors[1])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), fill=colors[2])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves),fill=colors[3])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs),fill=colors[4])+
  geom_line(col=colors[1],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), col=colors[2],size=1)+
  geom_line(aes(x =Semaine_lundi, y = Oeufs+Larves), col=colors[3],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs), col=colors[4], size=1)+
  ylim(0,max(na.omit(donneessummtot$Oeufs+donneessummtot$Larves+donneessummtot$Nymphes+donneessummtot$Ouvrieres)))+
  theme_classic()
graph1

graph2=ggplot(donneessummtot[donneessummtot$Nombre_ouvrieres==2,], aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres))+
  geom_rect(aes(xmin = as.Date("2020-11-05"), xmax = as.Date("2020-11-26"), ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = as.Date("2020-11-26"), xmax = as.Date("2021-04-06"), ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_vline(xintercept = as.Date("2020-11-26"), linetype="dotted")+
  geom_area( aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres),fill=colors[1])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), fill=colors[2])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves),fill=colors[3])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs),fill=colors[4])+
  geom_line(col=colors[1],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), col=colors[2],size=1)+
  geom_line(aes(x =Semaine_lundi, y = Oeufs+Larves), col=colors[3],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs), col=colors[4], size=1)+
  ylim(0,max(na.omit(donneessummtot$Oeufs+donneessummtot$Larves+donneessummtot$Nymphes+donneessummtot$Ouvrieres)))+
  theme_classic()
graph2

graph3=ggplot(donneessummtot[donneessummtot$Nombre_ouvrieres==4,], aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres))+
  geom_rect(aes(xmin = as.Date("2020-11-05"), xmax = as.Date("2020-11-26"), ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = as.Date("2020-11-26"), xmax = as.Date("2021-04-06"), ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_vline(xintercept = as.Date("2020-11-26"), linetype="dotted")+
  geom_area( aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes+Ouvrieres),fill=colors[1])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), fill=colors[2])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs+Larves),fill=colors[3])+
  geom_area(aes(x = Semaine_lundi, y = Oeufs),fill=colors[4])+
  geom_line(col=colors[1],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs+Larves+Nymphes), col=colors[2],size=1)+
  geom_line(aes(x =Semaine_lundi, y = Oeufs+Larves), col=colors[3],size=1)+
  geom_line(aes(x = Semaine_lundi, y = Oeufs), col=colors[4], size=1)+
  ylim(0,max(na.omit(donneessummtot$Oeufs+donneessummtot$Larves+donneessummtot$Nymphes+donneessummtot$Ouvrieres)))+
  theme_classic()
graph3


graph7=ggplot(donneessumm, aes(x=Semaine_lundi, y=croissance, group=as.factor(Nombre_ouvrieres), color=as.factor(Nombre_ouvrieres), ymin=croissance-sd, ymax=croissance+sd ))+
  geom_rect(aes(xmin = as.Date("2020-11-05"), xmax = as.Date("2020-11-26"), ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = as.Date("2020-11-26"), xmax = as.Date("2021-04-06"), ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_vline(xintercept = as.Date("2020-11-26"), linetype="dotted")+
  geom_ribbon(data=donneessumm[donneessumm$Nombre_ouvrieres == 0,], alpha = 0.2, mapping = aes(color = NULL), fill="#009E73" ) +
  geom_ribbon(data=donneessumm[donneessumm$Nombre_ouvrieres == 2,], alpha = 0.2, mapping = aes(color = NULL), fill="#E69F00" ) +
  geom_ribbon(data=donneessumm[donneessumm$Nombre_ouvrieres == 4,], alpha = 0.2, mapping = aes(color = NULL), fill="#D55E00" ) +
  scale_color_manual(values=c("#009E73","#E69F00","#D55E00"))+
  xlab("Temps")+
  ylab("Croissance")+
  geom_point(size=2)+ 
  geom_line(size=1)+
  theme_classic()+
  theme(legend.position="none")+
  labs(color="Nombre d'ouvrières")+
  guides(colour = guide_legend(override.aes = list(fill=c("red", "blue", "green"), alpha=1)))
graph7


#Figure 4

donnees_dpt =subset(donnees,donnees$Semaine_hiv == 27)
donnees_dpt$Nombre_ouvrieres = as.factor(donnees_dpt$Nombre_ouvrieres)
ggplot(donnees_dpt, aes(x = as.factor(Nombre_ouvrieres), y = croissance, fill=as.factor(Nombre_ouvrieres)))+
  scale_fill_manual(values=c("#009E73","#E69F00","#D55E00"))+
  geom_boxplot(alpha=0.8,outlier.alpha=0)+
  geom_jitter(shape=16, position=position_jitter(width=0.1, height=0))+
  theme_classic()


#Statistical analyzis of the growth

glm=glmer(croissance ~ (1|colonie) + Nombre_ouvrieres, data=donnees_dpt, family=poisson())
Anova(glm)
emmeans(glm, list(pairwise ~ Nombre_ouvrieres), adjust = "tukey")
shapiro.test(residuals(glm))

overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun(glm)


#Pic of growth before hibernation


donnees_pic1 =subset(donnees,donnees$Semaine_lundi == "2020-10-12")
donnees_pic1$Nombre_ouvrieres = as.factor(donnees_pic1$Nombre_ouvrieres)
donnees_pic1 = donnees_pic1[donnees_pic1$Statut == "VIVANTE" & donnees_pic1$Spermatheque == "PLEINE",]
mean(donnees_pic1$croissance[donnees_pic1$Nombre_ouvrieres==4])
sd(donnees_pic1$croissance[donnees_pic1$Nombre_ouvrieres==4])
mean(donnees_pic1$croissance[donnees_pic1$Nombre_ouvrieres==2])
sd(donnees_pic1$croissance[donnees_pic1$Nombre_ouvrieres==2])

#Pic of growth after hibernation


donnees_pic20 = subset(donnees,donnees$Semaine_lundi == "2021-06-28" & donnees$Nombre_ouvrieres == 0)
donnees_pic22 = subset(donnees,donnees$Semaine_lundi == "2021-06-14" & donnees$Nombre_ouvrieres == 2)
donnees_pic24 = subset(donnees,donnees$Semaine_lundi == "2021-06-14" & donnees$Nombre_ouvrieres == 4)
donnees_pic2 = rbind(donnees_pic20,donnees_pic22,donnees_pic24)
donnees_pic2$Nombre_ouvrieres = as.factor(donnees_pic2$Nombre_ouvrieres)
donnees_pic2 = donnees_pic2[donnees_pic2$Statut == "VIVANTE" & donnees_pic2$Spermatheque == "PLEINE",]

mean(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==4])
sd(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==4])
mean(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==2])
sd(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==2])
mean(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==0])
sd(donnees_pic2$croissance[donnees_pic2$Nombre_ouvrieres==0])


######################################
####Presence of queens in the nest####
######################################

data_nid = read.table("donnees_presence_nid.txt", h=T)
data_nid$Semaine = as.Date(x = as.character(data_nid$Semaine), format = "%d/%m/%Y")

#Figure 2

ggplot(data=data_nid, aes(x=Semaine, y=Presence_nid, group=as.factor(Conditions), color=as.factor(Conditions) ))+
  geom_rect(aes(xmin = as.Date("2020-11-05"), xmax = as.Date("2020-11-26"), ymin = -Inf, ymax = Inf), fill="grey", alpha=0.01, color="white" )+
  geom_rect(aes(xmin = as.Date("2020-11-26"), xmax = as.Date("2021-04-06"), ymin = -Inf, ymax = Inf), fill="darkgrey", alpha=0.01, color="white" )+
  geom_vline(xintercept = as.Date("2020-11-26"), linetype="dotted")+
  geom_line(size=2)+
  xlab("Semaines")+
  ylab("Proportion de reines à l'intérieur du nid")+
  scale_color_manual(values=c("#009E73","#E69F00","#D55E00"))+
  labs(color="Nombre d'ouvrières")+
  guides(colour = guide_legend(override.aes = list(fill=c("#009E73","#E69F00","#D55E00"), alpha=1)))+
  theme_classic()+
  theme(legend.position="bottom")