#packages used
library(tidyverse)
library(lme4)
library(emmeans)
library(ggeffects)
library(nlme)
library(ggbeeswarm)
library(MASS)
library(glmmTMB)
library(DHARMa)

#reading in data
dat<-read.csv(file="PupRecruitmentEstimates_2019-2022.csv") #entire dataset
dat3<-read.csv(file="PupRecruitmentEstimates_2019-2022_minus5litters.csv") #data minus the 5 litters of min/max estimates

#changing variable type
dat$recruited<-as.numeric(dat$recruited)
dat$packSize<-as.numeric(dat$packSize)
dat$year<-as.factor(dat$year)
dat$max<-as.numeric(dat$max)

#creating new dataset
dat2<-dat

#average pack size and recruitment by year
dat2 %>% group_by(year) %>% summarize(avg=mean(recruited),
                                      ps=mean(packSize))

#sample size by year and treatment
dat2 %>% group_by(year, type) %>% summarize(count=n())



#poisson model using entire dataset to assess impact of treatment type and year on recruitment
dat2$year <- relevel(as.factor(dat2$year), ref = "2021")
pupPOIS<-glmer(recruited~type+year+(1|pack), data=dat2, family = poisson(link = "log")) #with random effect
summary(pupPOIS)
confint(pupPOIS)

#poisson with max estimates of recruitment
dat2$year <- relevel(as.factor(dat2$year), ref = "2021")
pupMAX<-glmer(max~type+year+(1|pack), data=dat2, family = poisson(link = "log")) #with random effect
summary(pupMAX)
confint(pupMAX)


#poisson with data on 5 litters with min/max estimates removed
dat3$year <- relevel(as.factor(dat3$year), ref = "2021")
pupMinus<-glmer(recruited~type+year+(1|pack), data=dat3, family = poisson(link = "log")) #with random effect; random effect was tiny so getting singular fit.
summary(pupMinus)
confint(pupMinus)




#getting model estimates
handling<-ggeffect(pupPOIS, terms=c("type"), type="fixed")

#plotting effect of treatment type on pup recruitment
xy<-ggplot(handling, aes(x=x,y=predicted, color=x))+geom_point(size=8,position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin=conf.low, ymax=conf.high),size=2, position = position_dodge(width = 0.4))+theme_bw()+
  coord_cartesian(ylim=c(-0.05,5.1), expand=c(1,1))+
  geom_beeswarm(dat2, mapping=aes(x=type,y=recruited, color=type),size=3,cex=4,alpha=0.7, method="center")+
  ylab("Recruitment (pups/pack)")+theme_classic()+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))

#saving figure
ggsave("figures/Recruitment_DisturbanceFig.png",plot=xy, height=6, width=6 , dpi = 300)



#getting estimates from the poisson model
ann<-ggeffect(pupPOIS, terms=c("year"), type="fixed")

#getting average recruitment data by year
rawREC<-dat2 %>% group_by(year,type) %>% summarize(avg=mean(recruited))
rawREC$type<-as.factor(rawREC$type)

#renaming factor for plotting
dat2$year<-factor(dat2$year, levels=c('2020','2021','2022','2023'))

#plotting annual trends in recruitment by treatment
RR<-ggplot(dat2, aes(x=x,y=predicted))+
  geom_beeswarm(dat2, mapping=aes(x=year,y=recruited, color=type),size=3,cex=4,alpha=0.7, method="center")+
  geom_line(rawREC,mapping=aes(x=year,y=avg, group=type, color=type), size=1)+
  geom_point(rawREC,mapping=aes(x=year,y=avg, group=type, color=type), size=4)+
  ylab("Recruitment (pups/pack)")+theme_classic()+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position=c(0.2,0.9),
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))

ggsave("figures/AnnualRecruitment_RAWData_Fig.png",plot=RR, height=6, width=6 , dpi = 300)

#renaming factor for plotting
ann$x<-factor(ann$x, levels=c('2020','2021','2022','2023'))

#plotting annual trend in recruitment overall
rt<-ggplot(ann, aes(x=x,y=predicted))+geom_point(size=4)+
  geom_beeswarm(dat2, mapping=aes(x=year,y=recruited, color=type),size=3,cex=4,alpha=0.7, method="center")+
  geom_linerange(aes(ymin=conf.low, ymax=conf.high))+coord_cartesian(ylim=c(0,5))+
  ylab("Recruitment (pups/pack)")+theme_classic()+
  geom_point(size=4)+geom_line(aes(x=as.numeric(x),y=predicted), size=1)+
  geom_line(aes(x=as.numeric(x),y=predicted), size=1)+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position=c(0.2,0.9),
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))


ggsave("figures/AnnualRecruitment_Fig.png",plot=rt, height=6, width=6 , dpi = 300)





#pack size vs. marking of pu
psDat<-dat

#truncated poisson model to examine impact of marking and year on pack size
psDat$year <- relevel(as.factor(psDat$year), ref = "2021")
packPOIS<-glmmTMB(packSize~type+year+(1|pack), data=psDat, family = truncated_compois(link = "log"))
summary(packPOIS)
confint(packPOIS)
PR<-ggeffect(packPOIS, terms=c("type","year"), type="fixed")

#verifying model is well specified
simulationOutput <- simulateResiduals(fittedModel = pupPOIS)
plot(simulationOutput)

#getting model predictions
dv<-ggeffect(packPOIS, terms=c("type"), type="fixed")

#Effects of disturbance type on pack size
yy<-ggplot(dv, aes(x=x,y=predicted, color=x))+geom_point(size=8, position = position_dodge(width = 0.4))+
  geom_linerange(aes(ymin=conf.low, ymax=conf.high),size=2, position = position_dodge(width = 0.4))+
  coord_cartesian(ylim=c(1.8,8.2))+
  geom_beeswarm(dat2, mapping=aes(x=type,y=packSize, color=type),size=3,cex=4,alpha=0.7, method="center")+
  ylab("Pack size (wolves/pack)")+theme_classic()+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))

ggsave("figures/PackSize_DisturbanceFig.png",plot=yy, height=6, width=6 , dpi = 300)


#annual pack size estimates
ps<-ggeffect(packPOIS, terms=c("year"), type="fixed")

#getting data together for plotting
rawPS<-dat2 %>% group_by(year,type) %>% summarize(avg=mean(packSize))
rawPS$type<-as.factor(rawPS$type)
rawPS$year<-factor(rawPS$year, levels=c('2020','2021','2022','2023'))

#annual pack size trends by treatment type
RPS<-ggplot(ps, aes(x=x,y=predicted))+
  geom_beeswarm(dat2, mapping=aes(x=year,y=packSize, color=type),size=3,cex=4,alpha=0.7, method="center")+
  geom_line(rawPS,mapping=aes(x=year,y=avg, group=type, color=type), size=1)+
  geom_point(rawPS,mapping=aes(x=year,y=avg, group=type, color=type), size=4)+
  ylab("Pack size (wolves/pack)")+theme_classic()+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position="none",
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))

ggsave("figures/AnnualPackSize_RAWData_Fig.png",plot=RPS, height=6, width=6 , dpi = 300)

#overall trend in pack size over time
pa<-ggplot(ps, aes(x=x,y=predicted))+geom_point(size=4)+geom_line(aes(x=as.numeric(x),y=predicted), size=1)+
  geom_beeswarm(dat2, mapping=aes(x=year,y=packSize, color=type),size=3,cex=4,alpha=0.7, method="center")+
  geom_linerange(aes(ymin=conf.low, ymax=conf.high))+coord_cartesian(ylim=c(1.6,8.3))+
  ylab("Pack size (wolves/pack)")+theme_classic()+
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(size=24, family = "Gill Sans"), #editing font size of x-axis title
        axis.text.x = element_text(family = "Gill Sans", size=16),#editing the angle, placement and size of x-axis text
        axis.text.y = element_text(family = "Gill Sans",size=16),
        legend.position="none",
        legend.title = element_blank(),
        legend.text = element_text(size=16, family="Gill Sans"))+
  scale_colour_manual(values=c("Undisturbed"="#7570B3","Disturbed"="#E6AB02"))

ggsave("figures/AnnualPackSize_Fig.png",plot=pa, height=6, width=6 , dpi = 300)