library(readxl)
library(tidyverse)
everything<-read_excel("2022-PeerJ-MetzTarpy-revised-data.xlsx",sheet="Everything",skip=1)
source.rename<-tibble(Source=c(14,16,32,61,71,"DLD"),n.Source=c("A","B","C","D","E","F"))
everything<-left_join(everything,source.rename)
everything$Source<-NULL
names(everything)[names(everything)=="n.Source"]<-"Source"
#This ÒunweightsÓ the data for the dead drones
everything<-everything[rep(seq.int(1,nrow(everything)),everything$weight),]

#collection numbers
everything<-everything%>%mutate(Tagged=ifelse(is.na(em.Body)=="FALSE","Tagged","NT"))
everything<-everything%>%filter(c.Location!= "Introduction")
everything<-everything%>%filter(is.na(Source)!= "TRUE")

everything<-everything%>%filter(Tagged=="Tagged",c.Location=="Trap")%>%mutate(c.Condition="Dead")%>%full_join(.,filter(everything,Tagged!="Tagged"|c.Location!="Trap"))

everything<-everything%>%filter(Tagged=="Tagged",c.Location=="Lost")%>%mutate(c.Condition="Lost")%>%full_join(.,filter(everything,Tagged!="Tagged"|c.Location!="Lost"))

#Binary flight (censoring)
everything<- everything%>%mutate(c.Flight=as.logical(ifelse(c.Condition=="Live","TRUE","FALSE")))

#collection table(s)
everything%>%filter(Tagged=="Tagged")%>%group_by(Source,c.Condition,)%>%summarise(n=n())%>%pivot_wider(names_from="c.Condition",values_from="n")%>%write_csv("collection_tagged.csv")

everything%>%filter(Tagged=="NT")%>%group_by(Source,Foster,c.Condition, Dissected)%>%summarise(n=n())%>%mutate(Condition=ifelse(c.Condition=="Dead","Dead",ifelse(Dissected=="TRUE","Dissected","Tallied")))%>%subset(select=-c(3,4))%>%pivot_wider(names_from="Condition",values_from="n")%>%write_csv("collection_marked.csv")

everything%>%group_by(Source,Foster, c.Flight,c.Location)%>%summarise(n=n())%>%mutate(Flew=ifelse(c.Location!="Cage","Exclude",ifelse(c.Flight=="TRUE","Flew","Censored")))%>%subset(select=-c(c.Flight,c.Location))%>%group_by(Source,Foster,Flew)%>%summarise(n=sum(n))%>%pivot_wider(names_from="Flew",values_from="n")%>%subset(select=c(1,2,5,3,4))%>%write_csv("collection_flew.csv")

#survival analysis
library(survival)
library(survminer)
library(gtsummary)

#model setup
flew<-everything%>%filter(c.Location=="Cage")
flight<-Surv(time=flew$c.Age,event=flew$c.Flight)

#testing
survfit(flight~1)%>%quantile()
coxph(flight~Source*Foster,data=flew)%>%summary()

#because there was a significant interaction, we split foster colony out.
AM<-flew%>%filter(Foster=="AM")
PM<-flew%>%filter(Foster=="PM")
am.flight<-Surv(time=AM$c.Age,event=AM$c.Flight)
pm.flight<-Surv(time=PM$c.Age,event=PM$c.Flight)

survfit(am.flight~1)%>%quantile()
coxph(am.flight~Source,data=AM)%>%summary()

p2<-coxph(am.flight~Source,data=AM)%>%ggforest(main="Hazard Ratio (AM)") +geom_text(aes(x=.1,y=.9,label="c)"),size=10)

survfit(pm.flight~1)%>%quantile()
coxph(pm.flight~Source,data=PM)%>%summary()

p3<-coxph(pm.flight~Source,data=PM)%>%ggforest(main="Hazard Ratio (PM)")+geom_text(aes(x=.1,y=.9,label="d)"),size=10)

#capture age figure
medians<-flew%>%group_by(Source,Foster)%>%summarise(Age=median(c.Age,na.rm=TRUE),L=quantile(c.Age,probs=.25,na.rm=TRUE),H=quantile(c.Age,probs=.75,na.rm=TRUE))%>%as_tibble()
medians<-medians%>%mutate(Diff=H-Age,y=ifelse(Diff>0,Age+.5,Age-.5))

p1letter=tibble(label=c("a)","b)"),x=c(1,1),y=c(35,35),Foster=c("AM","PM"))
p1<- ggplot(data=flew,aes(x=Source,y=c.Age))+ geom_point(aes(color=Source),position=position_jitterdodge())+geom_boxplot(aes(fill=Source),outlier.alpha=0,alpha=.5)+geom_text(data=medians,aes(x=Source,y=y+10,label=Age))+geom_text(data=p1letter,aes(x=x,y=y,label=label),color="black",size=10)+facet_wrap(facets=vars(Foster))+theme(legend.position="none",axis.title=element_text(face="bold",size=15),strip.text=element_text(face="bold",size=15),axis.text=element_text(size=13))+labs(y="Capture Age (d)",x="Colony Source")

#figure 1
library(gridExtra) 
lay=rbind(c(1),c(2,3))
figure1<-grid.arrange(p1,p2,p3,layout_matrix=lay)


#emergence properties
#summary statistics
everything%>%filter(em.Body!= "NA")%>%summarise(mass=mean(em.Body),sd1=sd(em.Body),head=mean(el.Head,na.rm=TRUE),sd2=sd(el.Head,na.rm=TRUE),thorax=mean(el.Thorax,na.rm=TRUE),sd3=sd(el.Thorax,na.rm=TRUE))

#ANOVAs
Emerged<-everything%>%filter(em.Body!="NA")%>%subset(select=c("Source","em.Body","el.Head","el.Thorax"))%>%pivot_longer(cols = c(2:4),names_to = "parameter")
Parameters<-data.frame(Source=c("F","A","B","D","E","C"),cell=c(13.93,12.88,13.33,13.42,12.98,13.27),mites=c(16,21,3,13,4,9),n=c(23,33,0,112,34,12))
Emerged<-merge(Emerged,Parameters)
library(broom)
#this gives the Source-wise ANOVAs for each emergence characteristic
lms.Emerged<-Emerged%>%group_by(parameter)%>%nest()%>%mutate(models=lapply(data,function(df)anova(lm(value~Source,data=df))),results=map(models,tidy))%>%unnest(c(results))%>%subset(select=-c(2,3,6,7))%>%pivot_wider(names_from = "term",values_from = c("df","statistic","p.value"))%>%subset(select=-c(5,7))

#This will give the TukeyÕs significance groups
library(agricolae)
tukey.Emerged<-Emerged%>%group_by(parameter)%>%nest()%>%mutate(models=(lapply(data, function(df) HSD.test(aov(value~Source,data=df),"Source"))))
tukey.Emerged<-rbind(tukey.Emerged[[3]][[1]][["groups"]]%>%subset(select=2)%>%t()%>%as.data.frame()%>%'rownames<-'(tukey.Emerged$parameter[1]), tukey.Emerged[[3]][[2]][["groups"]]%>%subset(select=2)%>%t()%>%as.data.frame()%>%'rownames<-'(tukey.Emerged$parameter[2]), tukey.Emerged[[3]][[3]][["groups"]]%>%subset(select=2)%>%t()%>%as.data.frame()%>%'rownames<-'(tukey.Emerged$parameter[3]))

#building the tukey groups for an annotated figure
tukey.Groups<-as.data.frame(t(tukey.Emerged))
tukey.Groups<-tukey.Groups%>%mutate(Source=row.names(tukey.Groups))%>%pivot_longer(cols=c(1:3),names_to = "parameter",values_to = "group")
groupheight<-Emerged%>%group_by(Source,parameter)%>%summarise(Y=quantile(value,.73,na.rm=TRUE))
tukey.Groups<-left_join(tukey.Groups,groupheight)

letters<-tibble(y=c(3.77,4.87,173),x=c(1,1,1),parameter=c("el.Head","el.Thorax","em.Body"),label=c("a)","b)","c)"))
p2labs=c("Head width (mm)","Thorax width (mm)","Body mass (mg)")
names(p2labs)<-c("el.Head","el.Thorax","em.Body")

p2ac<- ggplot(Emerged,aes(x=Source,y=value))+ geom_point(aes(color=Source),position=position_jitterdodge())+geom_boxplot(aes(fill=Source),alpha=0.5,outlier.alpha=0)+geom_text(data=tukey.Groups,aes(x=Source,y=Y,label=group),vjust=1.1,size=3)+geom_text(data=letters,aes(x=x,y=y,label=label), size=6)+facet_wrap(facets = vars(parameter),scales = "free",labeller=labeller(parameter=p2labs))+labs(y=NULL)+theme(strip.text = element_text(face = "bold"))+expand_limits(x=c(0,7))+scale_fill_manual(values=c("#F8766D","#00BA38","#00BFC4","#619CFF","#F564E3"))+scale_color_manual(values=c("#F8766D","#00BA38","#00BFC4","#619CFF","#F564E3"))


#emergence survival analysis

emerged<-everything%>%filter(em.Body!="NA")%>%filter(c.Location=="Cage")
eflight<-Surv(time=emerged$c.Age,event=emerged$c.Flight)
survfit(eflight~1)%>%quantile()
coxph(eflight~Source,data=emerged)%>%summary()

coxph(eflight~em.Body+el.Head+el.Thorax+Source,data=emerged)%>%summary()

emmedians<-emerged%>%group_by(Source)%>%summarise(Age=median(c.Age,na.rm=TRUE),L=quantile(c.Age,probs=.25,na.rm=TRUE),H=quantile(c.Age,probs=.75,na.rm=TRUE))%>%as_tibble()
emmedians<-emmedians%>%mutate(Diff=H-Age,y=ifelse(Diff>0,Age+.5,Age-.5))

p2dlabs<-tibble(x=.75,y=24,label="d) ")

p2d<-ggplot(data=emerged,aes(x=Source,y=c.Age))+ geom_point(aes(color=Source),position=position_jitterdodge())+geom_boxplot(aes(fill=Source),outlier.alpha=0,alpha=.5)+geom_text(data=emmedians,aes(x=Source,y=y,label=Age))+geom_text(data=p2dlabs,aes(x=x,y=y,label=label),size=7)+theme(legend.position="none",axis.title=element_text(face="bold",size=15),strip.text=element_text(face="bold",size=15),axis.text=element_text(size=13))+labs(y="Capture Age (d)",x="Colony Source")+scale_fill_manual(values=c("#F8766D","#00BA38","#00BFC4","#619CFF","#F564E3"))+scale_color_manual(values=c("#F8766D","#00BA38","#00BFC4","#619CFF","#F564E3"))

p2e<- coxph(eflight~em.Body+el.Head+el.Thorax+Source,data=emerged)%>% ggforest() +geom_text(aes(x=.1,y=.9,label="e)"),size=7)


#figure 2
library(gridExtra) 
lay=rbind(c(1),c(2,3))
figure2<-grid.arrange(p2ac,p2d,p2e,layout_matrix=lay)

#Causes & correlates of drone variation in emergence characters
#mite count and cell size
lm(em.Body~Mite.Count,data=everything)%>%summary()
lm(el.Thorax~Mite.Count,data=everything)%>%summary()
lm(el.Head ~Mite.Count,data=everything)%>%summary()
lm(em.Body~Cell.Size,data=everything)%>%summary()
lm(el.Thorax~Cell.Size,data=everything)%>%summary()
lm(el.Head ~Cell.Size,data=everything)%>%summary()

#changes from emergence to flight
everything<-everything%>%mutate(dm.Body=em.Body-cm.Body,dl.Head=el.Head-cl.Head,dl.Thorax=el.Thorax-cl.Thorax)
everything%>%filter(dm.Body!="NA")%>%summarise(delta.Body=mean(dm.Body,na.rm=TRUE),se1=sd(dm.Body,na.rm=TRUE)/sqrt(n()))
everything%>%filter(dl.Head!="NA")%>%summarise(delta.Head=mean(dl.Head,na.rm=TRUE),se1=sd(dl.Head,na.rm=TRUE)/sqrt(n()))
everything%>%filter(dl.Thorax!="NA")%>%summarise(delta.Thorax=mean(dl.Thorax,na.rm=TRUE),se1=sd(dl.Thorax,na.rm=TRUE)/sqrt(n()))

t.test(everything$dm.Body)
t.test(everything$dl.Thorax)
t.test(everything$dl.Head)

#percentages
everything%>%subset(select=c("QD","e.Date","c.Date","c.Age","em.Body","el.Head","el.Thorax","cm.Body","cl.Head","cl.Thorax","dm.Body","dl.Head","dl.Thorax"))%>%na.omit()%>%mutate(pdm.Body=dm.Body/em.Body,pdl.Head=dl.Head/el.Head,pdl.Thorax=dl.Thorax/el.Thorax)%>%summarise(pdm.Body=mean(pdm.Body),pdl.Head=mean(pdl.Head),pld.Thorax=mean(pdl.Thorax))

#delta measures correlations
library(Hmisc)
everything%>%subset(select=c("c.Age","em.Body","cm.Body","dm.Body","el.Head","cl.Head","dl.Head","el.Thorax","cl.Thorax","dl.Thorax"))%>%filter(em.Body!="NA",cm.Body!="NA")%>%as.matrix()%>%rcorr(type="spearman")

#Age of capture, source, and emergence characteristics effects on reproductive params
everything%>%filter(cc.Sperm>0)%>%lm(log(cc.Sperm)~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>%aov()%>%summary()
everything%>%filter(cc.Sperm>0)%>%lm(log(cc.Sperm)~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>%summary()


everything%>%filter(cv.Sperm>0)%>%lm(asin(sqrt((cv.Sperm)))~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>%aov()%>%summary()
everything%>%filter(cv.Sperm>0)%>%lm(asin(sqrt((cv.Sperm)))~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>% summary()

everything%>%filter(cl.Mucus>0)%>%lm(cl.Mucus~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>%aov()%>%summary()
everything%>%filter(cl.Mucus>0)%>%lm(cl.Mucus~c.Age+Source+em.Body+el.Head+el.Thorax,data=.)%>% summary()

#mass effect figure
everything<-everything%>%mutate(bem.Body=everything$em.Body%>%cut(.,breaks=6))
p4<-everything%>%filter(em.Body!="NA")%>%filter(c.Location=="Cage")%>%ggplot(aes(x=c.Age,y=cc.Sperm))+geom_smooth(method="lm",aes(linetype=Source),color="grey50",fill=NA)+geom_point(aes(color=bem.Body,size=bem.Body))+scale_color_manual(values=c("#f1eef6","#d0d1e6","#a6bddb","#74a9cf","#2b8cbe","#045a8d"))+scale_linetype_manual(values=c("twodash","solid","dotted","longdash","dotdash"))

p4let<- ggplot_build(p4)$data[[1]]%>%group_by(linetype)%>%summarise(xmin=min(x),xmax=max(x),ymin=min(y),ymax=max(y))%>%mutate(Colony=c("F","D","E","C","A"))

p4a<-p4+geom_text(data=p4let,aes(y=ymin,x=xmin,label=Colony),fontface="bold",color="grey50",size=5)+geom_text (data=p4let,aes(y=ymax,x=xmax,label=Colony), fontface="bold",color="grey50",size=5) +theme(panel.background=element_rect(fill="white"),panel.grid=element_line(color="grey90"),legend.box="horizontal",legend.key=element_rect(fill="white"),legend.position=c(.26,.87),legend.background=element_rect(fill=NA))+labs(x="Capture Age (d)",y="Total sperm count (millions)",color="Emerged Body Mass (mg) (Min,Max]",size="Emerged Body Mass (mg) (Min,Max]")+scale_y_continuous(limits=c(0,18))+scale_x_continuous(breaks=c(5,10,15,20,25))+guides(linetype="none",color=guide_legend(nrow=3,byrow=TRUE),size=guide_legend(nrow=3,byrow=TRUE))

p42<-everything%>%filter(em.Body!="NA")%>%filter(c.Location=="Cage")%>%ggplot(aes(x=c.Age,y=cv.Sperm))+geom_smooth(method="lm",aes(linetype=Source),color="grey50",fill=NA)+geom_point(aes(color=bem.Body,size=bem.Body))+scale_color_manual(values=c("#fee5d9","#fcbba1","#fc9272","#fb6a4a","#de2d26","#a50f15"))+scale_linetype_manual(values=c("twodash","solid","dotted","longdash","dotdash"))

p42let<- ggplot_build(p42)$data[[1]]%>%group_by(linetype)%>%summarise(xmin=min(x),xmax=max(x),ymin=min(y),ymax=max(y))%>%mutate(Colony=c("F","D","E","C","A"))

p4b<-p42+geom_text(data=p42let,aes(y=ymin,x=xmin,label=Colony),fontface="bold",color="grey50",size=5)+geom_text (data=p42let,aes(y=ymax,x=xmax,label=Colony), fontface="bold",color="grey50",size=5) +theme(panel.background=element_rect(fill="white"),panel.grid=element_line(color="grey90"),legend.box="horizontal",legend.key=element_rect(fill="white"),legend.position=c(.77,.26),legend.background=element_rect(fill=NA))+labs(x="Capture Age (d)",y="Sperm viability  (Live/Total)",color="Emerged Body Mass (mg) (Min,Max]",size="Emerged Body Mass (mg) (Min,Max]")+scale_y_continuous(limits=c(0,1))+scale_x_continuous(breaks=c(5,10,15,20,25))+guides(linetype="none",color=guide_legend(nrow=3,byrow=TRUE),size=guide_legend(nrow=3,byrow=TRUE))

library(cowplot)
plot_grid(p4a,p4b,labels=c("a)","b)"),label_size=20)