---
title: "Population Analysis by Lineage - FS vs bottleneck Lineages. Related to Fig S8D"
author: "Danielle Ireland"
date: "2025-06-24"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Load libraries

```{r cars}
library(tidyverse)
library(dplyr)
library(readxl)
library(ggplot2)
library(rstatix)

set.seed(228)
```

## Read in division data for non-bottleneck lineages and only keep data for FS families.


```{r }
data_NB_original<-as.data.frame(read_excel("Dataset_S1.xlsx", sheet="Dataset S1A"))
data_NB<-data_NB_original %>%
  filter((Family==1) | (Family==4 ))
#create group variables (1=FS, 2=Bottleneck)
data_NB$Group<-1
```

# Read in bottleneck lineage data. Add Recode lineage and group names and isolate important variables
```{r }
data_B<-read_excel("Dataset_S1.xlsx", sheet="Dataset S1C")
data_B<-data_B %>% mutate(FamilyName=recode(Family,
                                            "0"="B1",
                                             "1"="B2"))
data_B$Group<-2
```


## Bin data by Birth Month to see if differences are persistent over time
```{r}
All_data<-rbind(data_NB, data_B)
#Bin birth month by every 20% of the data. Since all data last for 31 months, each bin is ~6 months
All_data <-All_data %>%
  mutate_at(c("FamilyName","Group"), factor)%>%
  group_by(FamilyName) %>%
  mutate(Birth.Bin=findInterval(Birth.Month, c(-Inf,quantile(Birth.Month,c(0.2, 0.4, 0.6, 0.8)),Inf)))%>%
  mutate_at("Birth.Bin", factor)

#group by FamilyName and find median and sample size
All_data_grouped<-split(All_data, All_data$FamilyName)
RWT_summary_stats_binned<-All_data %>% 
  group_by(FamilyName,Birth.Bin) %>% 
  summarise(median=median(TtD), n=n())

#Plot RWT per Birth bin (Fig S7D)
F_binned_plot<-ggplot(All_data_grouped$F, aes(x=Birth.Bin,y=TtD))+
  geom_violin(color='black') +
  stat_summary(fun=median, geom="point", size=2, color='black') +
  coord_flip()+
  xlab("Birth month bin (~6 months)")+
  ylab("Reproductive waiting time (days)")+
  ylim(0,300)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
S_binned_plot<-ggplot(All_data_grouped$S, aes(x=Birth.Bin,y=TtD))+
  geom_violin(color='black') +
  stat_summary(fun=median, geom="point", size=2, color='black') +
  coord_flip()+
  xlab("Birth month bin (~6 months)")+
  ylab("Reproductive waiting time (days)")+
  ylim(0,300)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
B1_binned_plot<-ggplot(All_data_grouped$B1, aes(x=Birth.Bin,y=TtD))+
  geom_violin(color='black') +
  stat_summary(fun=median, geom="point", size=2, color='black') +
  coord_flip()+
  xlab("Birth month bin (~6 months)")+
  ylab("Reproductive waiting time (days)")+
  ylim(0,300)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
B2_binned_plot<-ggplot(All_data_grouped$B2, aes(x=Birth.Bin,y=TtD))+
  geom_violin(color='black') +
  stat_summary(fun=median, geom="point", size=2, color='black') +
  coord_flip()+
  xlab("Birth month bin (~6 months)")+
  ylab("Reproductive waiting time (days)")+
  ylim(0,300)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))
library(gridExtra)
grid.arrange(F_binned_plot, S_binned_plot, B1_binned_plot,B2_binned_plot, ncol=2)


#statistical analysis: Kruskal Wallis test on each family separately 
res.kruskal_F<-All_data_grouped$F %>% kruskal_test (TtD~Birth.Bin)
res.kruskal_S<-All_data_grouped$S %>% kruskal_test (TtD~Birth.Bin)
res.kruskal_B1<-All_data_grouped$B1 %>% kruskal_test (TtD~Birth.Bin)
res.kruskal_B2<-All_data_grouped$B2 %>% kruskal_test (TtD~Birth.Bin)
KW_binned<-rbind(res.kruskal_F,res.kruskal_S,res.kruskal_B1,res.kruskal_B2)
KW_binned
```