############################################
#### Andrew Osborne Feb 2023 -August 2023 ####
#### Plots for Y1-Y3 dispersal ####
############################################

setwd("~/Desktop/PC_Folders/LHB/Manuscripts/BC_JIC")
list.files(pattern=".csv")
HD <-read.csv("Astley_Habitat_Dispersal_Data.csv")
names(HD)
nrow(HD)
#remove quadrats from outside any polygon
HD.P <-subset(HD,HD$Polygon!=0)
nrow(HD.P)
HD.P$Flt_Season <-as.factor(HD.P$Flt_Season)

##########################################
#### Plots  3 Years logistic curves ####
##########################################

library(ggplot2)
library(ggtext)

Flt_S.Pallette <-c("#E69F00", "#0072B2", "black" )

EvT_Seasons <-ggplot(HD.P, aes(x=Ev_tussock_count, y=C.t_PA, colour=Flt_Season)) +
  geom_smooth(method="glm",method.args=list(family="binomial"),fill="darkgrey",alpha=0.2,size=1.5)+
  coord_cartesian(xlim = c(0, 25), ylim = c(0.25, 1.0))+
  scale_colour_manual(values = c(Flt_S.Pallette),name="Flight\nSeason")+
  scale_y_continuous(breaks=seq(0.3,1.0,0.1))+
  theme(legend.position = c(.85,.2),
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        panel.grid.major = element_line(color = 'lightgrey'),
        panel.grid.minor = element_line(color = 'lightgrey'),
        text=element_text(size=14),
        axis.title.x = ggtext::element_markdown(),
        axis.title.y = ggtext::element_markdown())+
  xlab("*E.vaginatum* tussocks per 2 m quadrat")+
  ylab("Probability of *C.tullia* Presence")
EvT_Seasons

Et_Seasons <-ggplot(HD.P, aes(x=E.tetralix, y=C.t_PA, colour=Flt_Season)) +
  geom_smooth(method="glm",method.args=list(family="binomial"),fill="darkgrey",alpha=0.2,size=1.5)+
  coord_cartesian(xlim = c(0, 25), ylim = c(0.25, 1.0))+
  scale_colour_manual(values = c(Flt_S.Pallette),name="Flight\nSeason")+
  scale_y_continuous(breaks=seq(0.3,1.0,0.1))+
  theme(legend.position = c(.85,.2),
        panel.background = element_rect(fill = 'transparent', color = 'black'),
        panel.grid.major = element_line(color = 'lightgrey'),
        panel.grid.minor = element_line(color = 'lightgrey'),
        text=element_text(size=14),
        axis.title.x = ggtext::element_markdown(),
        axis.title.y = ggtext::element_markdown())+ 
  xlab("*E.tetralix* percentage cover")+
  ylab("Probability of *C.tullia* Presence")
Et_Seasons

##########################################
#### Segmented ####
##########################################

library(segmented)
library(oddsratio)
library(MASS)
nrow(HD.P)
#Ev_tussock_count
{set.seed(137)
EvT.glm <- glm(C.t_PA ~ Ev_tussock_count,  family=binomial, data=HD.P)
range(HD.P$Ev_tussock_count)
xEvT<-seq(0,22,0.01)
yEvT <- predict(EvT.glm, list(Ev_tussock_count =xEvT), type="response")

lin.modEvT <- lm(yEvT~xEvT)
segmented.mod.EvT <- segmented(lin.modEvT, seg.Z=~xEvT)

fit <- numeric(length(xEvT)) * NA
fit[complete.cases(rowSums(cbind(yEvT, xEvT)))] <- broken.line(segmented.mod.EvT)$fit
data.EvT <- data.frame(xEvT = xEvT, yEvT = yEvT, fit = fit)
}

#Stats Results (Table 5)
summary(segmented.mod.EvT)
segmented.mod.EvT$psi[,2]
data.EvT[findInterval(segmented.mod.EvT$psi[,2], data.EvT$xEvT),]
dose.p(EvT.glm, p=0.5)# @Midpoint
dose.p(EvT.glm, p=0.6947077)# @Breakpoint
summary(EvT.glm)
or_glm(data=HD.P, model=EvT.glm, incr = list(Ev_tussock_count=13.78))
#CI at Breakpoint
13.78-(3.014*1.96);13.78+(3.014*1.96)
#Plot segmented regression
EvT_Segmented <-ggplot(data.EvT, aes(x = xEvT, y = yEvT)) + 
  geom_point(size=0.25,colour="darkgrey") +
  coord_cartesian(xlim = c(0, 25), ylim = c(0.25, 1.0))+
  scale_y_continuous(breaks=seq(0.3,1.0,0.1))+
  theme(panel.background = element_rect(fill = 'transparent', color = 'black'),
        panel.grid.major = element_line(color = 'lightgrey'),
        panel.grid.minor = element_line(color = 'lightgrey'),
        text=element_text(size=14),
        axis.title.x = ggtext::element_markdown(),
        axis.title.y = ggtext::element_markdown())+
  geom_line(aes(x = xEvT, y = fit), color = "blue",size=1.5)+
  xlab("*E.vaginatum* tussocks per 2 m quadrat")+
  ylab("Probability of *C.tullia* Presence")+
  geom_segment(aes(x=13.78,y=0.25,xend=13.78,yend=0.695), size=0.5, linetype=2)+
  geom_segment(aes(x=0,y=0.5,xend=5.27,yend=0.5), size=0.5, linetype=2)+
  geom_segment(aes(x=5.27,y=0.5,xend=5.27,yend=0.25), size=0.5, linetype=2)+
  geom_segment(aes(x=0,y=0.695,xend=13.832,yend=0.695), size=0.5, linetype=2)+
  annotate(x=14.83,y=0.45,"text",label = "Breakpoint==13.78", angle = 90,parse = TRUE,size=4)+
  annotate(x=6.27,y=0.36,"text",label = "Midpoint==5.27", angle = 90,parse = TRUE,size=4)
EvT_Segmented

#E.tetralix
set.seed(137)#segmented works with iterations
E.t.glm <- glm(C.t_PA ~ E.tetralix,  family=binomial, data=HD.P)
range(HD.P$E.tetralix)
xE.t<-seq(0,23,0.01)
yE.t <- predict(E.t.glm, list(E.tetralix =xE.t), type="response")
lin.modE.t <- lm(yE.t~xE.t)
segmented.mod.E.t <- segmented(lin.modE.t, seg.Z=~xE.t)

fit <- numeric(length(xE.t)) * NA
fit[complete.cases(rowSums(cbind(yE.t, xE.t)))] <- broken.line(segmented.mod.E.t)$fit

data.E.t <- data.frame(xE.t = xE.t, yE.t = yE.t, fit = fit)


#Stats Results (Table 5)
summary(segmented.mod.E.t)
segmented.mod.E.t$psi[,2]
data.E.t[findInterval(segmented.mod.E.t$psi[,2], data.E.t$xE.t),]
dose.p(E.t.glm, p=0.5)# @Midpoint
dose.p(E.t.glm, p=0.727382)# @Breakpoint
summary(E.t.glm)
or_glm(data=HD.P, model=E.t.glm, incr = list(E.tetralix=13.355))
#CI at Breakpoint
13.35-(4.643*1.96);13.78+(4.643*1.96)
#Plot segmented regression
Et_Segmented <-ggplot(data.E.t, aes(x = xE.t, y = yE.t)) + 
  geom_point(size=0.25,colour="darkgrey") +
  coord_cartesian(xlim = c(0, 25), ylim = c(0.25, 1.0))+
  scale_y_continuous(breaks=seq(0.3,1.0,0.1))+
  theme(panel.background = element_rect(fill = 'transparent', color = 'black'),
        panel.grid.major = element_line(color = 'lightgrey'),
        panel.grid.minor = element_line(color = 'lightgrey'),
        text=element_text(size=14),
        axis.title.x = ggtext::element_markdown(),
        axis.title.y = ggtext::element_markdown())+ 
  geom_line(aes(x = xE.t, y = fit), color = "blue",size=1.5)+
  xlab("*E.tetralix* percentage cover")+
  ylab("Probability of *C.tullia* Presence")+
  geom_segment(aes(x=3.52,y=0.5,xend=3.52,yend=0.25), size=0.5, linetype=2)+
  geom_segment(aes(x=0,y=0.5,xend=3.52,yend=0.5), size=0.5, linetype=2)+
  geom_segment(aes(x=13.358,y=0.727,xend=13.358,yend=0.25), size=0.5, linetype=2)+
  geom_segment(aes(x=0,y=0.728,xend=13.358,yend=0.728), size=0.5, linetype=2)+
  annotate(x=14.36,y=0.5,"text",label = "Breakpoint==13.35", angle = 90,parse = TRUE,size=4)+
  annotate(x=4.52,y=0.36,"text",label = "Midpoint==3.52", angle = 90,parse = TRUE,size=4)
Et_Segmented


library(ggpubr)

#Arrange 4 plot types as single figure (Figure 5)
ggarrange(EvT_Seasons,ggplot() + theme_void(),EvT_Segmented,Et_Seasons,ggplot() + theme_void(),Et_Segmented,
          ncol=3, nrow=2,
          labels=c("A","","B","C","D"),
          widths=c(1,0.05,1))