Sys.setenv(LANGUAGE="en")
Sys.setlocale("LC_TIME", "English")
Sys.setenv(TZ='GMT')

library(ggplot2)
library(tidyr)
library(dplyr)
library(tidyverse)
library(lubridate)
library(suncalc)
library(MuMIn)
library(MASS)
library(sjPlot)
library(lubridate)
library(lme4)
rm(list=ls())
path <- file.path("C:", "Spawnseis", "Acell_calc", fsep="\\")
path
setwd(path)

 theme_JE<-theme_bw()+theme(plot.background = element_blank(),
                            panel.grid.major = element_blank(),panel.grid.minor = element_blank())+theme(axis.line = element_line(color = 'black'))
 theme_JE_2<-theme_bw(base_size = 14)+theme(plot.background = element_blank(),
                                            panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
   theme(axis.line = element_line(color = 'black'))
##Sound exposure dates Bakkasund
SE_Bakka<-tibble(EXP_Start_Date=c("10.02.2020","15.02.2021","13.02.2022",
                        "28.02.2024","30.06.2025",
                        "27.05.2023"),EXP_End_Date=c("14.02.2020","20.02.2021","17.02.2022",
                                                     "16.06.2023","01.03.2024","03.07.2025")) 
SE_Bakka<-SE_Bakka%>%mutate(EXP_Start_Date=as.Date(EXP_Start_Date,format="%d.%m.%Y",tz="UTC"),
                            EXP_End_Date=as.Date(EXP_End_Date,format="%d.%m.%Y",tz="UTC" ))
SE_Bakka<-SE_Bakka%>%mutate(Year=year(EXP_Start_Date),Month=month(EXP_Start_Date))
names(SE_Bakka) 
##FIELD DATA
###Data files for the field data, one for each spawning ground. 
Accel_Temp_Data_Osen<-read_csv("Accel_Data_Osen_Sub.csv")
Accel_Temp_Data_Bakkasund<-read_csv("Accel_Data_Bakka_Sub.csv")
names(Accel_Temp_Data_Osen)
names(Accel_Temp_Data_Bakkasund)

##Supplementary figure 3
Fem_1581336<-Accel_Temp_Data_Bakkasund%>%filter(Serial==1581336&Month%in%c(2:3)&Date>"2024-02-13"&Date<"2024-03-02")
names(Fem_1581336)
Mal_1581341<-Accel_Temp_Data_Bakkasund%>%filter(Serial==1581341&Month%in%c(2:3)&Date>"2024-02-13"&Date<"2024-03-02")
names(Mal_1581341)
Sup3a<-ggplot(Fem_1581336,aes(DateTime,Depth2))+geom_point(aes(col="red",alpha=0.3))+
  geom_point(data=Mal_1581341,aes(DateTime,Depth2),col="blue",alpha=0.3)+
  theme_JE_2+ylim(50,0)+ylab("Depth (m)")+theme(legend.position = "none")+
  scale_x_datetime(date_breaks = "2 days", date_labels =  "%d.%m")+xlab("Date") 
Sup3b<-ggplot(Fem_1581336,aes(DateTime,Acc))+geom_point(aes(col="red",alpha=0.3))+
  geom_point(data=Mal_1581341,aes(DateTime,Acc),col="blue",alpha=0.3)+
  theme_JE_2+ylim(0,5)+labs(y ="Activity levels"~~(m/s^2))+theme(legend.position = "none")+
  scale_x_datetime(date_breaks = "2 days", date_labels =  "%d.%m")+xlab("Date") 
require(gridExtra)  
Sup3<-grid.arrange(Sup3a,Sup3b,nrow=2)
ggsave(Sup3,file="Sup3.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
#Make imputed datasets to be used in later analyses
library(censlm)
Acc_Temp_All<-rbind(Accel_Temp_Data_Bakkasund,Accel_Temp_Data_Osen)
Acc_Temp_All_Mat<-Acc_Temp_All%>%filter(MatStageComb=="Mat"&str_detect(Sensor.Unit,pattern="s"))#1498210
Acc_sub<-Acc_Temp_All_Mat%>%filter(Acc>1)
Acc_rest<-Acc_Temp_All_Mat%>%filter(Acc<=1)#Make decay function work
Acc_sub_1<-Acc_sub%>%filter(TagYear!=2024)
Acc_sub_2<-Acc_sub%>%anti_join(Acc_sub_1)
obs <- replace(Acc_sub_1$Acc, Acc_sub_1$Acc> 3.46494, 3.46494)
obs2<- replace(Acc_sub_2$Acc, Acc_sub_2$Acc > 4.9011, 4.9011)
fit <- clm(log(I(obs)) ~ 1, right = log(3.46494))
summary(fit)
imputed_data <- exp(imputed(fit))
imputed_data_frame<-as.data.frame(imputed_data)

##Percentage of imputed data less than 4.9011, i.e. less than the maximum value in 
imp_sum<-imputed_data_frame%>%filter(imputed_data>3.46494)
imp_sum2<-imputed_data_frame%>%filter(imputed_data>4.9011)
PerCent<-1-(length(imp_sum2$imputed_data)/length(imp_sum$imputed_data))

Acc_Temp_All_Mat_imputed_1<-cbind(Acc_sub_1,imputed_data_frame)
fit2 <- clm(log(I(obs2)) ~ 1, right = log(4.9011))
summary(fit2)
imputed_data <- exp(imputed(fit2))
imputed_data_frame2<-as.data.frame(imputed_data)
Acc_Temp_All_Mat_imputed_2<-cbind(Acc_sub_2,imputed_data_frame2)
Acc_rest<-Acc_rest%>%mutate(imputed_data=Acc)
Acc_Temp_All_Mat_imputedb<-rbind(Acc_Temp_All_Mat_imputed_1,Acc_Temp_All_Mat_imputed_2,Acc_rest)
Acc_Temp_All_Mat_imputedb<-Acc_Temp_All_Mat_imputedb%>%mutate(censored=ifelse(Acc==max(Acc),1,0))
Acc_Temp_All_Mat_imputedb<-Acc_Temp_All_Mat_imputedb%>%filter(is.na(Depth2)!=TRUE)
Acc_Temp_All_Mat_imputedb<-Acc_Temp_All_Mat_imputedb%>%mutate(imp2=ifelse(Acc<3.46494&TagYear<2024,Acc,
                                                                   ifelse(Acc<4.9011&TagYear==2024,Acc,imputed_data)))

Acc_Temp_All_Mat_imputedb_osen_bakka<-Acc_Temp_All_Mat_imputedb
Acc_Temp_All_Mat_imputedb_osen_bakka<-Acc_Temp_All_Mat_imputedb_osen_bakka%>%mutate(TagLocation=ifelse(TagLocation=="Osen","Osen","Bakkasund"))


library(nlme)
library(lme4)
library(pscl)
library(lmerTest)
names(Acc_Temp_All_Mat_imputedb_osen_bakka)
Acc_Temp_All_Mat_imputedb_osen_bakka_24<-Acc_Temp_All_Mat_imputedb_osen_bakka%>%
  filter(TagYear==2024)
Acc_Temp_All_Mat_imputedb_osen_bakka_24<-Acc_Temp_All_Mat_imputedb_osen_bakka_24%>%arrange(Serial,DateTime)%>%
  mutate(LineN=1:577797)
ML<-tibble(LineN=c(seq(1,577797,5)))
Acc_Temp_All_Mat_imputedb_osen_bakka_24<-Acc_Temp_All_Mat_imputedb_osen_bakka_24%>%inner_join(ML)
Acc_Temp_All_Mat_imputedb_osen_bakka_rest<-Acc_Temp_All_Mat_imputedb_osen_bakka%>%
  filter(TagYear!=2024)
Acc_Temp_All_Mat_imputedb_osen_bakka_24<-Acc_Temp_All_Mat_imputedb_osen_bakka_24%>%dplyr::select(-LineN)

Acc_Temp_All_Mat_imputedb_osen_bakka_thinned<-rbind(Acc_Temp_All_Mat_imputedb_osen_bakka_24,
                                                    Acc_Temp_All_Mat_imputedb_osen_bakka_rest)
##Models including only sex on accell
Bakka_spawn<-Acc_Temp_All_Mat_imputedb_osen_bakka_thinned%>%filter(Main_Area!="osen"&Month%in%c(2:3))
Bakka_feed<-Acc_Temp_All_Mat_imputedb_osen_bakka_thinned%>%filter(Main_Area=="Bakkasund"&Month%in%c(5:12))
Osen_spawn<-Acc_Temp_All_Mat_imputedb_osen_bakka_thinned%>%filter(Main_Area=="osen"&Month%in%c(2:3))
Osen_feed<-Acc_Temp_All_Mat_imputedb_osen_bakka_thinned%>%filter(Main_Area=="osen"&Month%in%c(5:12))
##Filtering out data past exposure dates in spawning and feeding period.
names(Bakka_spawn)
SE_Bakka_Spawn<-SE_Bakka%>%filter(Month%in%c(2:3))%>%dplyr::select(EXP_Start_Date,EXP_End_Date,Year)
Bakka_spawn<-Bakka_spawn%>%left_join(SE_Bakka_Spawn,by=c("Year"="Year"))
unique(Bakka_spawn$Year)
#Pre_filt 686103
Bakka_spawn_filt<-Bakka_spawn%>%filter(Year!=2019|Year!=2023)%>%filter(Date<EXP_Start_Date|Date>EXP_End_Date)
Bakka_spawn_filt2<-Bakka_spawn%>%filter(Year==2019|Year==2023)
Bakka_spawn_2<-rbind(Bakka_spawn_filt,Bakka_spawn_filt2)
names(Bakka_spawn_2)
Bakka_spawn_2<-Bakka_spawn_2%>%dplyr::select(-EXP_Start_Date,-EXP_End_Date)
##659610. Binned by 1 m and remove data points deeper than 60 m
Spawn_All<-rbind(Bakka_spawn_2,Osen_spawn)
Spawn_All_sub<-Spawn_All%>%filter(Year==TagYear&Depth_round<61)
Spawn_All_subb<-Spawn_All%>%filter(Year==TagYear)

Spawn_sum_Area<-Spawn_All_sub%>%group_by(Main_Area)%>%summarise(Det=n())

##Feeding pre-filt
SE_Bakka_Feed<-SE_Bakka%>%filter(Month%in%c(5:12))%>%dplyr::select(EXP_Start_Date,EXP_End_Date,Year)
Bakka_feed<-Bakka_feed%>%left_join(SE_Bakka_Feed,by=c("Year"="Year"))
Bakka_feed_filt<-Bakka_feed%>%filter(Year==2023)%>%filter(Date<EXP_Start_Date|Date>EXP_End_Date)
Bakka_feed_filt2<-Bakka_feed%>%filter(Year!=2023)
Bakka_feed_2<-rbind(Bakka_feed_filt,Bakka_feed_filt2)
Bakka_feed_2<-Bakka_feed_2%>%dplyr::select(-EXP_Start_Date,-EXP_End_Date)
Feed_All<-rbind(Bakka_feed_2,Osen_feed)
Feed_All_sub<-Feed_All%>%filter(Year==TagYear&Depth_round<61)
Feed_All_subb<-Feed_All%>%filter(Year==TagYear)

Feed_sum_Area<-Feed_All_sub%>%group_by(Main_Area)%>%summarise(Det=n())



##Create input for overview table of tagged fish in the field, i.e. Table 2 in MS 
Table_Field<-rbind(Spawn_All_sub,Feed_All_sub)
names(Table_Field)
Table_Field<-Table_Field%>%mutate(TagLocation=ifelse(TagLocation=="Osen","Osen","Bakkasund"))%>%
  group_by(TagLocation,TagYear,Sex,Serial)%>%summarise(MLength=mean(Length),Det=n())%>%
  group_by(TagLocation,TagYear,Sex)%>%summarise(MLength=mean(MLength/10),Numb=n())
write.csv(Table_Field,"Table_2_EE_revised_thinned_new.csv",row.names = FALSE)


##Create files with aggregated data for analyses
#Spawning
Test_File_Gen<-Spawn_All_sub%>%
group_by(Serial,Sex,Length,Year,Main_Area)%>%mutate(HV=ifelse(imputed_data>2,1,0))%>%
  summarise(Acc_mean=mean(imputed_data),P_HV=(sum(HV)/n()),Det=n(),Period="Spawning")

Test_File_Gen_Bakka<-Test_File_Gen%>%filter(Main_Area=="Bakkasund")
Test_File_Gen_Osen<-Test_File_Gen%>%filter(Main_Area=="osen")


#Feeding
Test_File_GenF<-Feed_All_sub%>%
  group_by(Serial,Sex,Length,Year,Main_Area)%>%
  summarise(Acc_mean=mean(imputed_data),Det=n(),Period="Feeding")
Test_File_GenF_Bakka<-Test_File_GenF%>%filter(Main_Area=="Bakkasund")
Test_File_GenF_Osen<-Test_File_GenF%>%filter(Main_Area=="osen")

#Starting models for generic test on accell differences between sexes in the field 
#Spawning periods
#Bakkasund
Mod_Bakka<-lmer(Acc_mean~Sex*Length+(1|Year),weights = log(Det),data=Test_File_Gen_Bakka,na.action = "na.fail")
summary(Mod_Bakka)
#Dredge of initial model
M2<-dredge(Mod_Bakka)
M2
#Chosen model
Mod_Bakka_fin<-lmer(Acc_mean~Sex+(1|Year),weights = log(Det),subset(Test_File_Gen,Main_Area=="Bakkasund"))
summary(Mod_Bakka_fin)


#Osen
Mod_osen_lm<-lm(Acc_mean~Sex*Length,weights = log(Det),data=Test_File_Gen_Osen,na.action = "na.fail")
summary(Mod_osen_lm)
M2<-dredge(Mod_osen_lm)
M2
#Chosen model
Mod_osen_lm_fin<-lm(Acc_mean~Sex,weights = log(Det),data=Test_File_Gen_Osen,na.action = "na.fail")
summary(Mod_osen_lm_fin)


#Final generic models  
tab_model(Mod_Bakka_fin,Mod_osen_lm_fin,file="Table_4_EE_thinned_new.doc")

##Feeding
Mod_BakkaF<-lmer(Acc_mean~Sex*Length+(1|Year),weights = log(Det+1),data=Test_File_GenF_Bakka,na.action = "na.fail")
summary(Mod_BakkaF)

AICc(Mod_BakkaF)
M2<-dredge(Mod_BakkaF)
M2

#Chosen model
Mod_BakkaF_fin<-lmer(Acc_mean~1+(1|Year),weights = log(Det+1),data=Test_File_GenF_Bakka,na.action = "na.fail")
summary(Mod_BakkaF_fin)

Mod_osenF_lm<-lm(Acc_mean~Sex*Length,weights = log(Det),data=Test_File_GenF_Osen,na.action = "na.fail")
summary(Mod_osenF_lm)
M2<-dredge(Mod_osenF_lm)
M2

#Chosen model
Mod_osenF_lm_fin<-lm(Acc_mean~Length,weights = log(Det),data=Test_File_GenF_Osen,na.action = "na.fail")
summary(Mod_osenF_lm_fin)
##Feeding results-supp table 1
tab_model(Mod_BakkaF_fin,Mod_osenF_lm_fin,file="Sup_Table3_EE_thinned_new.doc")

################################More complex models with depth as categorical variable
Test_File_Cat<-Spawn_All_sub%>%mutate(Depth_Cat=ifelse(Depth2<=20,0,1))%>%filter(Main_Area=="Bakkasund")%>%
  group_by(Serial,Sex,Length,Depth_Cat,Year,Main_Area)%>%summarise(Acc_mean=mean(imputed_data),
                                                                               Acc_median=median(imputed_data),Det=n())%>%ungroup()

#################################
Test_File_Cat<-Spawn_All_sub%>%mutate(Depth_Cat=ifelse(Depth2<=20,0,1),HV=ifelse(imputed_data>2,1,0))%>%filter(Main_Area=="Bakkasund")%>%
  group_by(Serial,Sex,Length,Depth_Cat,Year,Main_Area)%>%summarise(Acc_mean=mean(imputed_data),
                                                                   Acc_median=median(imputed_data),P_HV=sum(HV)/n(),
                                                                   Det=n())%>%ungroup()
Test_File_Cat<-Test_File_Cat%>%mutate(Depth_CatF=ifelse(Depth_Cat==0,"<20m",">20m"))


Test_File_Cat_osen<-Spawn_All_sub%>%mutate(Depth_Cat=ifelse(Depth2<=20,0,1),
                                           HV=ifelse(imputed_data>2,1,0))%>%filter(Main_Area!="Bakkasund")%>%
  group_by(Serial,Sex,Length,Depth_Cat,Year,Main_Area)%>%summarise(Acc_mean=mean(imputed_data),
                                                                             Acc_median=median(imputed_data),
                                                                   P_HV=sum(HV)/n(),Det=n())%>%ungroup()
Test_File_Cat_osen<-Test_File_Cat_osen%>%mutate(Depth_CatF=ifelse(Depth_Cat==0,"<20m",">20m"))


library(glmmTMB)
##Bakkasund
Bakka_Cat<- glmmTMB(Acc_mean ~ Sex*factor(Depth_Cat)+
                     (1|Serial)+(1|Year),family = Gamma(link = "log"),                      
                   control=glmmTMBControl(optCtrl = list(maxit=200)),
                   data = Test_File_Cat,weights = log(Det+1),
                   na.action = "na.fail")
summary(Bakka_Cat)
M2<-dredge(Bakka_Cat)
M2#Full model chosen

names(Test_File_Cat)


##Propotion High Accell Analyses
transform01 <- function(x) {
  (x * (length(x) - 1) + 0.5) / (length(x))
}
Test_File_Cat<-Test_File_Cat%>%
  mutate(P_HV_Scaled=transform01(P_HV))
Bakka_Cat_HV<- glmmTMB(P_HV_Scaled ~ Sex*factor(Depth_Cat)+
                         (1|Serial)+(1|Year),family=beta_family(link="logit"),                      
                       control=glmmTMBControl(optCtrl = list(maxit=200)),
                       data = Test_File_Cat,weights = log(Det+1),
                       na.action = "na.fail")

summary(Bakka_Cat_HV)
M2<-dredge(Bakka_Cat_HV)
M2#Full model chosen

##Osen
osen_Cat<- glmmTMB(Acc_mean ~ Sex*factor(Depth_Cat)+
                                  (1|Serial),family = Gamma(link = "log"),                      
                                 control=glmmTMBControl(optCtrl = list(maxit=200)),
                                 data = Test_File_Cat_osen,weights = log(Det+1),
                                 na.action = "na.fail")
summary(osen_Cat)
M2<-dredge(osen_Cat)
M2#Full model chosen

Test_File_Cat_osen<-Test_File_Cat_osen%>%mutate(P_HV_Scaled=transform01(P_HV))
Osen_Cat_HV<- glmmTMB(P_HV_Scaled ~ Sex*factor(Depth_Cat)+
                         (1|Serial),family=beta_family(link="logit"),                      
                       control=glmmTMBControl(optCtrl = list(maxit=200)),
                       data = Test_File_Cat_osen,weights = log(Det),
                       na.action = "na.fail")

summary(Osen_Cat_HV)

M2<-dredge(Osen_Cat_HV)
M2#Full model chosen
#Saved model output for Table 5 in MS
tab_model(Bakka_Cat,osen_Cat,file="Table5_EE_thinned_new.doc",transform = NULL)
tab_model(Bakka_Cat_HV,Osen_Cat_HV,file="Table6_EE_thinned_new.doc",transform = NULL)
##Figures for manuscript
names(Test_File_Cat)
names(Test_File_Cat_osen)
Test_File_Cat_Tot<-rbind(Test_File_Cat,Test_File_Cat_osen)
Test_File_Cat_Tot<-Test_File_Cat_Tot%>%
  mutate(Main_Area=ifelse(Main_Area=="Bakkasund","Bakkasund","Osen"))
Test_File_Cat_Tot<-rbind(Test_File_Cat,Test_File_Cat_osen)
Test_File_Cat_Tot<-Test_File_Cat_Tot%>%mutate(Main_Area=ifelse(Main_Area=="osen","Osen",Main_Area))

Fig.4_newa<-ggplot(Test_File_Cat_Tot,aes(Depth_CatF,Acc_mean,col=Sex))+
  geom_boxplot(outlier.shape = NA)+facet_wrap(~Main_Area)+
  theme_JE_2+theme(legend.position = "none")+
  scale_color_manual(breaks = c("F","M"),values = c("red","blue"))+
  labs(y ="Activity levels"~~(m/s^2))+xlab("Depth")
Fig.4_newa
Fig.4_newb<-ggplot(Test_File_Cat_Tot,aes(Depth_CatF,P_HV,col=Sex))+
  geom_boxplot(outlier.shape = NA)+facet_wrap(~Main_Area)+
  theme_JE_2+theme(legend.position = "none")+coord_cartesian(ylim=c(0,0.3))+
  scale_color_manual(breaks = c("F","M"),values = c("red","blue"))+
  labs(y ="Prop_HV")+xlab("Depth")
Fig.4_newb


library(cowplot)
library(egg)
Fig4_EE<-ggarrange(Fig.4_newa+
                            theme(axis.title.x = element_blank(),
                                  axis.text.x = element_blank()),
                          Fig.4_newb,nrow=2)
Fig4_EE
#ggsave(Fig.4_new,file="Fig4_JAE.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
#ggsave(Fig.4_newb,file="Fig5_EE.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
ggsave(Fig4_EE,file="Fig4_EE_thinned_new.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")  
Test_File_All<-Spawn_All_sub%>%mutate(Depth_Round=round(Depth2))%>%
  group_by(Serial,Sex,Length,Depth_round,Year,Main_Area,Day_Night)%>%summarise(Acc_mean=mean(imputed_data),
                                                                               Acc_median=median(imputed_data),Det=n())%>%ungroup()

Test_File_All_Feed<-Feed_All_sub%>%mutate(Depth_Round=round(Depth2))%>%
  group_by(Serial,Sex,Length,Depth_round,Year,Main_Area,Day_Night)%>%summarise(Acc_mean=mean(imputed_data),
                                                                               Acc_median=median(imputed_data),Det=n())%>%ungroup()
Test_File_All_sub<-Test_File_All%>%filter(Depth_round>-1)
Test_File_All_Feed<-Test_File_All_Feed%>%filter(Depth_round>-1)
names(Test_File_All_sub)
Test_File_All_sub<-Test_File_All_sub%>%mutate(Period="Spawning")
Test_File_All_Feed<-Test_File_All_Feed%>%mutate(Period="Feeding")
Test_File_All_Comb<-rbind(Test_File_All_sub,Test_File_All_Feed)
Test_File_All_Comb<-Test_File_All_Comb%>%mutate(Main_Area=ifelse(Main_Area=="Bakkasund","Bakkasund","Osen"))
Test_File_All_Comb<-Test_File_All_Comb%>%mutate(Period=fct_relevel(Period,"Spawning"))

#Figure 2 in MS
Fig.2<-ggplot(Test_File_All_Comb,aes(Depth_round,Acc_mean,group=Sex,col=Sex))+
  stat_summary(size=0.4)+
  facet_grid(Period~Main_Area)+theme_JE_2+coord_cartesian(ylim=c(0,2.4),xlim=c(0,61))+
  theme(legend.position = c(0.85,0.90))+geom_vline(xintercept = c(20),lty=2)+
  scale_color_manual(breaks=c("F","M"),values = c("red","blue"))+
  labs(y ="Activity levels"~~(m/s^2))+xlab("Depth (m)")
Fig.2
ggsave(Fig.2,file="Fig2_EE_thinned_new.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")

##Make split violin plots, i.e.figure 5
##Code below copied from net. Makes violin plots split comparing akes and females side by side
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

##Now new geom can be used on data from the locations overall distrubtion across stations withn area
names(Acc_Temp_All_Mat)
Acc_Temp_All_Mat<-Acc_Temp_All_Mat%>%filter(Month!=1&Month!=4&Sensor.Unit!="m")
Acc_Temp_All_Mat<-Acc_Temp_All_Mat%>%filter(Main_Area=="Bakkasund"|Main_Area=="osen")
Acc_Temp_All_Mat<-Acc_Temp_All_Mat%>%mutate(Day_NightF=ifelse(Day_Night==0,"Day","Night"),
                                            Main_Area=ifelse(Main_Area=="osen","Osen",Main_Area))

##Figure 5
Fig.5<-ggplot(subset(Acc_Temp_All_Mat,Period==0),aes(Main_Area,Depth2,fill=factor(Sex)))+geom_split_violin()+theme_JE_2+
  theme(legend.position = c(0.9,0.2))+xlab("Spawning Ground")+ylab("Depth (m)")+scale_y_reverse()+#facet_wrap(~Main_Area)+
  scale_fill_manual(values = c("red","blue"))+labs(fill="")
Fig.5

ggsave(Fig.5,file="Fig5_EE_thinned_new.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")



##Net_Pen data
Net_Pen_Accell<-read_csv("Net_Pen_Accell_sub.csv")
names(Net_Pen_Accell)
Net_Pen_Accell$Date<-as.Date(Net_Pen_Accell$datetime)
Net_Pen_Accell<-Net_Pen_Accell%>%mutate(Month=month(datetime))

Table_NetPen<-Net_Pen_Accell%>%group_by(serial,year,bag,sex,length)%>%summarise(Det=n())%>%
  group_by(year,bag,sex)%>%summarise(MLength=mean(length/10),Numb=n())
write.csv(Table_NetPen,file="Table_1EE_new.csv",row.names = FALSE)

##Do data imputation before further analyses
Net_Pen_Accell_F<-Net_Pen_Accell
Imp_sub<-Net_Pen_Accell_F%>%filter(value>1)
Imp_rest<-Net_Pen_Accell_F%>%filter(value<=1)#Make decay function work
obs <- replace(Imp_sub$value, Imp_sub$value> 3.46494, 3.46494)
#clm.control(maxModIter=10)
fit <- clm(log(I(obs)) ~ 1, right = log(3.46494))
summary(fit)
imputed_data <- exp(imputed(fit))
imputed_data_frame<-as.data.frame(imputed_data)
max(imputed_data)

Netpen_imputed<-cbind(Imp_sub,imputed_data_frame)
Imp_rest<-Imp_rest%>%mutate(imputed_data=value)
Net_pen_imp_data<-rbind(Netpen_imputed,Imp_rest)
names(Net_pen_imp_data)
Sup1a<-ggplot(Net_pen_imp_data,aes(value))+geom_histogram(fill="red",col="black",binwidth = 0.07)+
  labs(x ="Activity levels"~~(m/s^2))+coord_cartesian(xlim = c(0,10),ylim = c(0,3000))+
  theme_JE_2
Sup1a
Sup1b<-ggplot(Net_pen_imp_data,aes(imputed_data))+geom_histogram(fill="red",col="black",binwidth = 0.07)+
  labs(x ="Activity levels"~~(m/s^2))+coord_cartesian(xlim = c(0,10),ylim = c(0,3000))+
  theme_JE_2
Sup1b
require(gridExtra)
Fig_Sup<-grid.arrange(Sup1a,Sup1b,nrow=2)
ggsave(Fig_Sup,file="FigS1_EE_new.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
names(Net_pen_imp_data)
Net_sum<-Net_pen_imp_data%>%group_by(bag,year)%>%summarise(Det=n())
Net_Pen_Exposure<-tibble(bag=c("exp1","exp2",1:3),EXP_Start_Date=c("01.04.2019","01.04.2019","09.03.2020",
    "09.03.2020","09.03.2020"),EXP_End_Date=c("06.04.2019","06.04.2019","11.03.2020","11.03.2020","11.03.2020"))
Net_Pen_Exposure<-Net_Pen_Exposure%>%mutate(EXP_Start_Date=as.Date(EXP_Start_Date,format="%d.%m.%Y",tz="UTC"),
                         EXP_End_Date=as.Date(EXP_End_Date,format="%d.%m.%Y",tz="UTC" ))
Net_Pen_Exposure<-Net_Pen_Exposure%>%mutate(year=year(EXP_Start_Date))
Net_pen_imp_data_sum_start<-Net_pen_imp_data%>%left_join(Net_Pen_Exposure,by=c("year"="year","bag"="bag"))%>%
  filter(Month!=4)
Net_pen_imp_data_suma<-Net_pen_imp_data_sum_start%>%filter(bag!="con")%>%
  filter(Date<EXP_Start_Date|Date>EXP_End_Date)%>%mutate(HV=ifelse(imputed_data>2,1,0))%>%
  group_by(serial,sex,year,length)%>%
  summarise(Mean_Acc=mean(imputed_data),P_HV=sum(HV)/n(),Numb=n())
Net_pen_imp_data_sumb<-Net_pen_imp_data_sum_start%>%filter(bag=="con")%>%
 mutate(HV=ifelse(imputed_data>2,1,0))%>%
  group_by(serial,sex,year,length)%>%
  summarise(Mean_Acc=mean(imputed_data),P_HV=sum(HV)/n(),Numb=n())

Net_pen_imp_data_sum<-rbind(Net_pen_imp_data_suma,Net_pen_imp_data_sumb)

sum(Net_pen_imp_data_sum$Numb)

names(Net_pen_imp_data)
lm_general_netpen <- lm(Mean_Acc ~ sex*length, data=Net_pen_imp_data_sum,weights = log(Numb),na.action = "na.fail")#,weights = Numb)
summary(lm_general_netpen)
M2<-dredge(lm_general_netpen)
M2
lm_general_netpen_fin <- lm(Mean_Acc ~ sex+length, data=Net_pen_imp_data_sum,weights = log(Numb),na.action = "na.fail")#,weights = Numb)
summary(lm_general_netpen_fin)
##Model output for Table 2 in MS
tab_model(lm_general_netpen_fin,file="Table3_EE_new.doc")
##Make plot showing overall values for field and net - pen data, i.e. Figure 3 in MS 
Plot_Test_All<-rbind(Test_File_Gen,Test_File_GenF)
Plot_Test_All<-Plot_Test_All%>%mutate(Sp_Period=paste(Main_Area,Period,sep = "_"))
Plot_Test_All<-Plot_Test_All%>%mutate(Period=fct_relevel(Period,"Spawning"))


names(Plot_Test_All)
names(Net_pen_imp_data_sum)
Plot_Test_All_m<-Plot_Test_All%>%dplyr::select(Sex,Acc_mean,Sp_Period)
Net_pen_imp_data_sum<-Net_pen_imp_data_sum%>%mutate(Sex=sex,Acc_mean=Mean_Acc,Sp_Period="Net_Pen")
Net_pen_imp_data_sum_m<-Net_pen_imp_data_sum%>%ungroup%>%dplyr::select(Sex,Acc_mean,Sp_Period,Acc_mean)%>%
  ungroup()
Net_pen_imp_data_sum_m<-Net_pen_imp_data_sum%>%ungroup%>%dplyr::select(Sex,Acc_mean,Sp_Period,Acc_mean)%>%
  ungroup()
Plot_Test_All_m<-Plot_Test_All_m%>%mutate(Sp_Period=ifelse(Sp_Period=="Bakkasund_Spawning","Bakka_Spawn",
                                         ifelse(Sp_Period=="Bakkasund_Feeding","Bakka_Feed",
                                         ifelse(Sp_Period=="osen_Feeding","Osen_Feed","Osen_Spawn"))))
Net_Fiel<-rbind(Plot_Test_All_m,Net_pen_imp_data_sum_m)
Net_Fiel$Sp_Period<-fct_relevel(Net_Fiel$Sp_Period,"Net_Pen","Bakka_Spawn",
                                 "Osen_Spawn","Bakka_Feed","Osen_Feed")
Fig.3<-ggplot(Net_Fiel,aes(Sp_Period,Acc_mean,col=Sex))+stat_summary(position = position_dodge(width = 0.1))+theme_JE+
  scale_color_manual(breaks = c("F","M"),values = c("red","blue"))+theme_JE_2+
  theme(legend.position = c(0.8,0.8))+xlab("")+ylab("Activity levels"~~(m/s^2))
Fig.3
ggsave(Fig.3,file="Fig3_EE_thinned_new.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")