library(Rmisc)
library(ggplot2)
library(tidyverse)
library(tidytext)
library(dplyr)
library(stringr)
library (MASS)
library(performance)
library(DHARMa)

#Read in data
PupSurvivalandRecruitmentEstimates_2019to2025 <-read.csv("PupSurvivalandRecruitmentEstimates_2019to2025.csv")

##Annual averages of survival and recruitment metrics 
annual_averages <- PupSurvivalandRecruitmentEstimates_2019to2025 %>%
  group_by(Year) %>%
  summarise(total_pups_born = sum(Litter_Size[!is.na(Litter_Size)], na.rm = TRUE), total_pups_survived = sum(Recruitment[!is.na(Litter_Size)], na.rm = TRUE), annual_survival = total_pups_survived / total_pups_born,
    mean_recruitment = mean(Recruitment, na.rm = TRUE),
    n_packs_total = n(),
    n_packs_survival = sum(!is.na(Litter_Size)))
#View annual results
print(annual_averages)

##Overall averages across all years (mean of annual metrics)
overall_averages <- annual_averages %>%
  summarise(mean_annual_survival = mean(annual_survival, na.rm = TRUE),mean_annual_recruitment = mean(mean_recruitment, na.rm = TRUE), n_years = n())
#View overall results
print(overall_averages)


##Read in new data set for regression tests
LitterSizeRegressionData <- read.csv("LitterSizeRegressionData.csv")

##Binomial Logistic Regression for comparison of Litter Size and Pup Survival
#Add failures column
LitterSizeRegressionData$failures <- LitterSizeRegressionData$Littersize - LitterSizeRegressionData$Pups_Recruited

# Fit Logistic regression
model <- glm(cbind(Pups_Recruited, failures) ~ Littersize,
             data = LitterSizeRegressionData, family = binomial)
summary(model)

#Model fit check
DHARMa::simulateResiduals(model) |> DHARMa::plotResiduals()
performance::check_model(model)
performance::model_performance(model)

# Calculate p-value for annotation
p_val <- signif(coef(summary(model))["Littersize", "Pr(>|z|)"], 3)

# Prediction for plotting
new_data <- data.frame(Littersize = seq(min(LitterSizeRegressionData$Littersize),
                                        max(LitterSizeRegressionData$Littersize), by=0.1))
pred <- predict(model, newdata=new_data, type="link", se.fit=TRUE)
new_data$fit <- plogis(pred$fit)
new_data$lower <- plogis(pred$fit - 1.96 * pred$se.fit)
new_data$upper <- plogis(pred$fit + 1.96 * pred$se.fit)

# Plot regression
ggplot(LitterSizeRegressionData, aes(Littersize, Pups_Recruited / Littersize)) +
  geom_point(size=2) +
  geom_line(data=new_data, aes(Littersize, fit), color="navy", size=1) +
  geom_ribbon(data=new_data, aes(ymin=lower, ymax=upper), fill="blue", alpha=0.2) +
  scale_y_continuous(limits=c(-0.05, 1), breaks=seq(0,1,0.2)) +
  labs(x="Litter Size", y="Survival Rate") +
  annotate("text", x=max(LitterSizeRegressionData$Littersize)-0.5, y=0.95,
           label=paste0("p = ", round(p_val, 3)), hjust=1, vjust=1, size=5)


##Linear Regression for comparison of Litter Size and Pup Recruitment##
#Regression Model Y ~ X
fit <- lm(Littersize ~ Pups_Recruited, data = LitterSizeRegressionData)  #create a model called "fit"
summary(fit)

sessionInfo()