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'))
 
##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)

#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)
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)
##Models including only sex on accell
Bakka_spawn<-Acc_Temp_All_Mat_imputedb_osen_bakka%>%filter(Main_Area!="osen"&Month%in%c(2:3))
Bakka_feed<-Acc_Temp_All_Mat_imputedb%>%filter(Main_Area=="Bakkasund"&Month%in%c(5:12))
Osen_spawn<-Acc_Temp_All_Mat_imputedb_osen_bakka%>%filter(Main_Area=="osen"&Month%in%c(2:3))
Osen_feed<-Acc_Temp_All_Mat_imputedb%>%filter(Main_Area=="osen"&Month%in%c(5:12))

##Binned by 1 m and remove data points deeper than 60 m
Spawn_All<-rbind(Bakka_spawn,Osen_spawn)
Spawn_All_sub<-Spawn_All%>%filter(Year==TagYear&Depth_round<61)
Feed_All<-rbind(Bakka_feed,Osen_feed)
Feed_All_sub<-Feed_All%>%filter(Year==TagYear&Depth_round<61)

##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_JAE.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)%>%
  summarise(Acc_mean=mean(imputed_data),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_JAE.doc")

##Feeding
Mod_BakkaF<-lmer(Acc_mean~Sex*Length+(1|Year),weights = log(Det),data=Test_File_GenF_Bakka,na.action = "na.fail")
summary(Mod_BakkaF)
M2<-dredge(Mod_BakkaF)
M2
#Chosen model
Mod_BakkaF_fin<-lmer(Acc_mean~1+(1|Year),weights = log(Det),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_Table1_JAE.doc")

##More complex models with depth as categorical variable
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_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_osen<-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<-Test_File_Cat%>%mutate(Depth_CatF=ifelse(Depth_Cat==0,"<20m",">20m"))



library(glmmTMB)
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),
                   na.action = "na.fail")
summary(Bakka_Cat)
M2<-dredge(Bakka_Cat)
M2#Full model chosen


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),
                                 na.action = "na.fail")
summary(osen_Cat)
M2<-dredge(osen_Cat)
M2
#Full model chosen
#Saved model output for Table 5 in MS
tab_model(Bakka_Cat,osen_Cat,file="Table5_JAE.doc",transform = NULL)



##Figures for manuscript
Test_File_Cat_osen<-Test_File_Cat_osen%>%mutate(Depth_CatF=ifelse(Depth_Cat==0,"<20m",">20m"))
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_new<-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_new
ggsave(Fig.4_new,file="Fig4_JAE.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
  
names(Test_File_Cat)
names(Test_File_Cat_osen)
Test_File_Cat$Main_Area<-"Bakkasund"
Test_File_Cat_osen$Main_Area<-"Osen"

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()+
  facet_grid(Period~Main_Area)+theme_JE_2+coord_cartesian(ylim=c(0,2.3),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_JAE.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_JAE.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_1JAE.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)

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_JAE.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")
Net_pen_imp_data_sum<-Net_pen_imp_data%>%filter(Month!=4)%>%group_by(serial,sex,year,length)%>%
  summarise(Mean_Acc=mean(imputed_data),Numb=n())


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_JAE.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()
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_JAE.tiff",dpi = 150,width = 7.04,height = 6.62,bg="white")