---
title: "2_biomass_analysis"
author: "Bram Parmentier"
output: html_document
---

# Biomass INLA analysis for Parmentier et al (2025): "Small fish biomass in the North Sea is far greater than previously estimated"

Objective:
- Create biomass predictions for:
  - all small species combined
  - 7 most abundant species individual + all remaining species
  - 8-31g individuals (Greenstreet 2007)
  - Seal prey for DFS and BTS sampling areas (Aarts et al 2019)

load libraries
```{r, warning=F}
library(xlsx);library(ggplot2)
library(sf)
library(INLA);library(mapdata)
library(stringr); library(fields)
library(gstat); library(lattice);library(latticeExtra)
library(dplyr); library(VGAM);library(raster); library(sp)
library(patchwork)
library(mgcv)
```

load data
```{r}
# Set locale to English
Sys.setlocale("LC_TIME", "C")

#set path
path <- ""

trawldata <- read.csv(file.path(path,"/trawldata_depth_mud.csv"))
fishdata <- read.csv(file.path(path,"/fishbiomass.csv"))

fishdata <- fishdata[!is.na(fishdata$StationID),]
fishdata <- fishdata[fishdata$StationID %in% trawldata$StationID,]
trawldata <- trawldata[trawldata$StationID %in%  fishdata$StationID,] 

#environmental data
depth <- read.csv(file.path(path,"/other data/depth_NorthSea.csv"))
depth$depth <- round(depth$depth*-1,0)

mud <- read.csv(file.path(path,"/other data/mud.csv"))
mud$mud_perc <- round(mud$mud_perc,2)

#round environmental data
trawldata$depth <- round(trawldata$depth,0)
trawldata$mud_perc <- round(trawldata$mud_perc,2)

#BTS and DFS shapefiles
load(file.path(path,"shapefiles/DFS_BTS_areas.rdata")) 
```

Flat or round fish larger than 30cm are discarded. The used gear (Triple-D) caught them only very occasionally but they contribute largely to the observed biomass. However large (>=30cm) individuals of Ammodytidae, Syngnathiformes, and Myxine glutinosa are not removed. 
```{r}
large_individuals <- fishdata[fishdata$length_cm >= 30 & !is.na(fishdata$length_cm),]
stayin <- c("Ammodytidae", "Syngnathiformes", "Myxine glutinosa")
large_individuals <- large_individuals[!large_individuals$Species_reported %in% stayin,]
table(large_individuals$Species_reported)

#large fish in sample with smaller fish
Sample <- large_individuals[large_individuals$Weight_type == "Sample" & !is.na(large_individuals$Weight_type),]
Sample <- Sample[Sample$Species_reported != "Amblyraja radiata",] #these ray are both larger than 30cm

#update WW_g by removing the expected weight of the large individual
fishdata$StationID_fish <- paste0(fishdata$StationID, fishdata$Species_reported) #will be removed 
Sample$StationID_fish <- paste0(Sample$StationID, Sample$Species_reported) 

for ( i in unique(Sample$StationID_fish)) {
  fishdata[fishdata$StationID_fish == i,]$WW_g[1] <- sum(fishdata[fishdata$StationID_fish == i,]$WW_g, na.rm=T) - fishdata[fishdata$StationID_fish == i & fishdata$length_cm >= 30,]$expected_WW_g
}

#update fishdata
fishdata <- fishdata[!fishdata$unique_code %in% large_individuals$unique_code,!names(fishdata) %in% ("StationID_fish")]
```

Not all records have a WW_g (either 'Sample' or 'Entry'), a new column is made which also includes the predicted weight of fish. This predicted weight (WW_g_pred_group) consist of weights derived from the length-weight species specific relationship (when a length was available), based on the length distribution of measured fish, or broken estimated weight. 
```{r}
fishdata$Weight <- fishdata$WW_g
weight_type <- c("Entry", NA)
fishdata[is.na(fishdata$Weight) & fishdata$Weight_type %in% weight_type,]$Weight <- fishdata[is.na(fishdata$Weight) & fishdata$Weight_type %in% weight_type,]$WW_g_pred_group
```

In some cases the catch is subsampled and also the fished distance differs between samples. In the models this is dealt with by applying an offset.
```{r}
fishdata$Fraction <- as.numeric(fishdata$Fraction)
table(fishdata$Fraction)
table(trawldata$Dist)
```

However, in some cases there are multiple fractions within one sample station. This happened, for example, when there were many gobies, and only in half of the catch is processed for gobies while other species the whole fraction was sorted.
```{r}
#which stations do have multiple 
multipleFractions <- as.data.frame(table(fishdata$StationID, fishdata$Fraction))
colnames(multipleFractions) <- c("StationID","Fraction", "Freq")

multipleFractions[multipleFractions$Freq > 0,]$Freq <- 1 #to create only present value
multipleFractions <- aggregate(Freq~StationID, multipleFractions, FUN="sum")
multipleFractions <- multipleFractions[multipleFractions$Freq > 1,]
```

We solved this by dividing the weight (count is not analysed) by fraction in the case there are multiple fractions
```{r}
fishdata_fraction_cor <- fishdata
fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Weight <- fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Weight/fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Fraction

fishdata_fraction_cor[fishdata_fraction_cor$StationID %in% multipleFractions$StationID,]$Fraction <- 1
```

combine trawl and fishdata_fraction_cor
```{r}
combined <- merge(fishdata_fraction_cor, trawldata, by="StationID", all.x = T) 
```

Weight per station is determined. This will be the response variable in the models.
```{r}
#all fish
combined[is.na(combined$Weight),]$Weight <- 0 
totperstation <- aggregate(Weight ~ StationID, combined, sum, drop = FALSE) 

#add fraction and Date
totperstation <- merge(totperstation, combined[!duplicated(combined$StationID),c("StationID", "Fraction"),], by= "StationID", all.x = T )

#create utm
trawldata_sf <- st_as_sf(trawldata, coords = c("Lon_mid", "Lat_mid"), crs = 4326)
trawldata_utm <- st_transform(trawldata_sf, crs = 32631)
utm_coords <- st_coordinates(trawldata_utm)
trawldata$Xkm <- utm_coords[,1]/1000
trawldata$Ykm <- utm_coords[,2]/1000

#include trawl information with utm
totperstation <- merge(totperstation, trawldata, by= "StationID", all.x = T)
totperstation <- totperstation[!is.na(totperstation$Lat_mid),]

#exclude Survey which goes to the north part of north sea
totperstation <- totperstation[totperstation$Lat_mid < 56,]
totperstation <- totperstation[totperstation$Fraction > 0.1,] #too low fraction not reliable
```

Observed weight per station, but corrected for fraction and fished distance
```{r}
totperstation$Weight_cor <- totperstation$Weight/totperstation$Fraction
totperstation$Weight_cor <- totperstation$Weight_cor/totperstation$Dist*5 #width of dregde is 0.2m, so times 5 to get to m2

summary(totperstation$Weight_cor)

#Fig S6
#pdf(file = file.path(path, "figures/supplement/S6_observedbiomass_hist.pdf"), width = 170 / 25.4, height = 170 / 25.4)  # mm to inches
hist(totperstation$Weight_cor, breaks=100, xlab= expression("Wet weight (g m"^{-2}*")"), main="", ylab="Triple-D samples", 
     xaxt='n')
     axis(1, at=seq(0, max(totperstation$Weight_cor +10 , na.rm = TRUE), by=10), las=1)
     abline(v=c(quantile(totperstation$Weight_cor, probs = c(0.05, 0.95))), col="red", lty=2)
     abline(v=quantile(totperstation$Weight_cor, probs = 0.5), col="blue", lty=2)
#dev.off()
```

#INLA spatial mesh
first load in some spatial objects
```{r}
#land areas
land <- st_read(file.path(path,"shapefiles/WGS84/WGS84geoKust500.shp"))
land <- as(land, "Spatial")
projection(land) <- CRS("+proj=longlat +datum=WGS84")
land_sf <- st_as_sf(land)

#coast from informatiehuis marien
coast <- st_read(file.path(path,"shapefiles/AU_baseline.shp")) 
coast <- as(coast, "Spatial")
proj_string <- "+proj=longlat +datum=WGS84 +no_defs"
projection(coast) <- CRS(proj_string)

#creating shape files for the areas of interest
#Dogger bank
dog_2more <- c("64PE438#141", "64PE438#142", "64PE438#151")
Loc_dog <- cbind(totperstation[totperstation$StationID %in% dog_2more,]$Lon_mid, totperstation[totperstation$StationID %in% dog_2more,]$Lat_mid)
Loc_dog1 <- cbind(totperstation[totperstation$Survey == "ACTN23",]$Lon_mid, totperstation[totperstation$Survey == "ACTN23" ,]$Lat_mid)
Loc_dog <- rbind(Loc_dog, Loc_dog1)
ConvHull_dog <- inla.nonconvex.hull(points=Loc_dog, convex=-0.1, resolution=50)

#polygon
doggerbank <- as.data.frame(cbind(ConvHull_dog$loc[,1], ConvHull_dog$loc[,2]))
colnames(doggerbank) <- c("Long", "Lat")
coords <- cbind(doggerbank$Long, doggerbank$Lat)
doggerbank_sf <- st_as_sf(SpatialPolygons(list(Polygons(list(Polygon(coords)), ID = "doggerbank"))))
st_crs(doggerbank_sf) <- proj_string

#NCP
Loc_NCP <- cbind(totperstation[totperstation$Survey != "ACTN23" & !totperstation$StationID %in% dog_2more,]$Lon_mid, totperstation[totperstation$Survey != "ACTN23" & !totperstation$StationID %in% dog_2more,]$Lat_mid)
ConvHull_NCP <- inla.nonconvex.hull(points=Loc_NCP, convex=-0.03, resolution=100)

#NCP polygon is overlapping coastline
plot(ConvHull_NCP)
lines(land)

#creating polygon
NCP <- as.data.frame(cbind(ConvHull_NCP$loc[,1], ConvHull_NCP$loc[,2]))
colnames(NCP) <- c("Lon", "Lat")
coords <- cbind(NCP$Lon, NCP$Lat)
NCP_polygon <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID = "NCP")))

#Create coast spatial pologyon instead of a line 
temp <- cbind(7,51)
coast@lines[[1]]@Lines[[1]]@coords <- rbind(coast@lines[[1]]@Lines[[1]]@coords, temp)

#close polygon
closed <- sapply(coast@lines, function(x) all(x@Lines[[1]]@coords[1, ] == x@Lines[[1]]@coords[nrow(x@Lines[[1]]@coords), ]))
if (!all(closed)) {
  for (i in which(!closed)) {
    coast@lines[[i]]@Lines[[1]]@coords <- rbind(coast@lines[[i]]@Lines[[1]]@coords, coast@lines[[i]]@Lines[[1]]@coords[1, ])
  }
}

#creating polygon
coords <- coast@lines[[i]]@Lines[[1]]@coords
Coast_polygon <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID = "NCP")))

#Convert SpatialPolygons objects to sf objects
NCP_polygon_sf <- st_as_sf(NCP_polygon)
Coast_polygon_sf <- st_as_sf(Coast_polygon)
NCP_polygon_modified_sf <- st_difference(NCP_polygon_sf, Coast_polygon_sf)
NCP_polygon_modified_sf <- st_difference(NCP_polygon_sf, Coast_polygon_sf)
st_crs(NCP_polygon_modified_sf ) <- proj_string 

#doggerbank 40m depth line
dog40a <- read.csv(file.path(path,"shapefiles/iso40mDog_a.csv"))
dog40b <- read.csv(file.path(path,"shapefiles/iso40mDog_b.csv"))
dog40a_sf <- st_as_sf(dog40a, coords = c("Long.40m.dog.a", "Lat.40m.dog.a"), crs = proj_string)
dog40a_line <- dog40a_sf %>%
  st_combine() %>%   
  st_cast("LINESTRING") 
dog40b_sf <- st_as_sf(dog40b, coords = c("Long.40m.dog.b", "Lat.40m.dog.b"), crs = proj_string)
dog40b_line <- dog40b_sf %>%
  st_combine() %>%   
  st_cast("LINESTRING")
```

Location of samples will be used to generate spatial meshes for the data.
```{r, eval=T, echo=TRUE}
Loc <- cbind(totperstation$Xkm, totperstation$Ykm)
colnames(Loc) <- c("Longitude", "Latitude")
D <- dist(Loc)
hist(D, breaks = 4000)
plot(x=sort(D), y=(1:length(D)/length(D)), type="l", xlab= "Distance between haul locations", ylab= "Cumulative proportion")

#in longitude and latitude
Loc <- cbind(totperstation$Lon_mid, totperstation$Lat_mid)
colnames(Loc) <- c("Longitude", "Latitude")
plot(Loc, asp=1)
lines(land)
```

Creating INLA mesh, with Dogger bank area and the Dutch EEZ
```{r}
ConvHull <- inla.nonconvex.hull(points=Loc, convex=-0.063, resolution=100) 
mesh <- inla.mesh.2d(boundary=ConvHull, max.edge=c(0.09), cutoff = 0.05) 
plot(mesh)
points(Loc)
mesh$n

#plot in ggplot
# Extract nodes (vertices)
nodes_df <- as.data.frame(mesh$loc)
colnames(nodes_df) <- c("Longitude", "Latitude")
proj_string <- "+proj=longlat +datum=WGS84 +no_defs"
nodes_sf <- st_as_sf(nodes_df, coords = c("Longitude", "Latitude"), crs = proj_string)

# Extract edges from the triangles
edges_indices <- mesh$graph$tv
edges_list <- lapply(seq_len(nrow(edges_indices)), function(i) {
  # Define the coordinates for each edge of the triangle
  edge_coords <- matrix(c(
    nodes_df[edges_indices[i, 1], "Longitude"], nodes_df[edges_indices[i, 1], "Latitude"],
    nodes_df[edges_indices[i, 2], "Longitude"], nodes_df[edges_indices[i, 2], "Latitude"],
    nodes_df[edges_indices[i, 2], "Longitude"], nodes_df[edges_indices[i, 2], "Latitude"],
    nodes_df[edges_indices[i, 3], "Longitude"], nodes_df[edges_indices[i, 3], "Latitude"],
    nodes_df[edges_indices[i, 3], "Longitude"], nodes_df[edges_indices[i, 3], "Latitude"],
    nodes_df[edges_indices[i, 1], "Longitude"], nodes_df[edges_indices[i, 1], "Latitude"]
  ), ncol = 2, byrow = TRUE)
  st_sfc(st_linestring(edge_coords), crs = proj_string)
})

# Create sf object for edges
edges_sf <- st_sf(geometry = do.call(c, edges_list))

xlim <- c(0,7)
ylim <- c(52,56)

loc_sf <- st_as_sf(data.frame(Loc), coords = c("Longitude", "Latitude"), crs = proj_string)

#Fig. S4
ggplot() +
  geom_sf(data = nodes_sf, color = "grey", size = 0.5, shape = 16) +  
  geom_sf(data = edges_sf, color = "lightblue", linewidth=0.005) +         
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  geom_sf(data = loc_sf, color = "brown", size = 0.5, shape = 16) + 
  geom_sf(data = doggerbank_sf, fill = "lightblue", color = "blue", alpha = 0.3) + 
  geom_sf(data = NCP_polygon_modified_sf , fill = "lightblue", color = "blue", alpha = 0.3) +
  coord_sf() + coord_sf(xlim = xlim, ylim = ylim) +
  labs(x = "Longitude", y = "Latitude") +
  theme_bw() + theme(text = element_text(size = 7.34)) 
#ggsave(file = file.path(path, "figures/supplement/S4_INLAmesh.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```
create prediction dataframe (per 1 km2)
```{r}
rangex <- range(totperstation$Xkm)
rangey <- range(totperstation$Ykm)
gridcellsize <- 1 

## sf object unique grid_pred points
grid_pred <- expand.grid(XkmRounded = seq(rangex[1]-100, rangex[2]+100, gridcellsize),
                    YkmRounded = seq(rangey[1]-100, rangey[2]+100, gridcellsize))

original_crs <- "+proj=utm +zone=31 +datum=WGS84 +units=km +no_defs"
grid_pred_sf <- st_as_sf(grid_pred, coords = c("XkmRounded", "YkmRounded"), crs = original_crs )
grid_pred_sf_longlat  <- st_transform(grid_pred_sf, crs = proj_string)

#Extract coordinates (longitude and latitude) from the sf object
coordinates <- st_coordinates(grid_pred_sf_longlat)
grid_pred_longlat <- data.frame(Lon_mid = coordinates[, 1], 
                           Lat_mid = coordinates[, 2], 
                           grid_pred_sf_longlat)

#only in selected areas
points_within_DOG <- st_within(grid_pred_sf_longlat, doggerbank_sf, sparse = FALSE)
points_within_NCP <- st_within(grid_pred_sf_longlat, NCP_polygon_modified_sf, sparse = FALSE)

grid_pred_longlat$inside_DOG <- points_within_DOG[, 1]
grid_pred_longlat$inside_NCP <- points_within_NCP[, 1]

grid_pred_longlat <- grid_pred_longlat[grid_pred_longlat$inside_DOG == T | grid_pred_longlat$inside_NCP == T,]
grid_pred_longlat <- grid_pred_longlat[,!(names(grid_pred_longlat) %in% c("geometry"))] 

plot(grid_pred_longlat$Lon_mid, grid_pred_longlat$Lat_mid, asp=1.5, xlab="Longitude", ylab="Latitude")
```

Set predictions in predication dataframe for a haul distiance of 50m, and a sorting fraction of 1. In the EEZ interest area the Survey effect of "2019feb_187" was set as baseline, likewise the DOY (Day of Year since first of October) was set. In the UK Dogger Bank area we used Survey "2023Jun_65" and DOY = 265 as baselines. Furthermore, depth and mud content values falling outside the observed range at the sampling locations were capped at the respective minimum and maximum values recorded at those locations to prevent extrapolation. 
```{r}
#prediction defined
grid_pred_longlat$Fraction <- 1
grid_pred_longlat$Dist <- 50

#EEZ
grid_pred_longlat$Survey2 <- "2019Feb_187"
grid_pred_longlat$DOY_adjusted <- 143

#UK Dogger Bank
grid_pred_longlat[grid_pred_longlat$inside_DOG == T,]$Survey2 <- "2023Jun_65"
grid_pred_longlat[grid_pred_longlat$inside_DOG == T,]$DOY_adjusted <- 265 

#add nearest depth
grid_pred_longlat_sf <- st_as_sf(grid_pred_longlat, coords = c("Lon_mid", "Lat_mid"), crs = 4326)
depth_sf <- st_as_sf(depth, coords = c("long", "lat"), crs = 4326)
nearest_indices <- st_nearest_feature(grid_pred_longlat_sf, depth_sf)
nearest_depth <- depth_sf[nearest_indices, ]

grid_pred_longlat <- cbind(
  grid_pred_longlat,
  depth = nearest_depth$depth
)

#add nearest mud
mud_sf <- st_as_sf(mud, coords = c("Longitude", "Latitude"), crs = 4326)
nearest_indices <- st_nearest_feature(grid_pred_longlat_sf, mud_sf)
nearest_mud <- mud_sf[nearest_indices, ]

grid_pred_longlat <- cbind(
  grid_pred_longlat,
  mud_perc = nearest_mud$mud_perc#,
  #mud_logit = nearest_mud$mud_logit
)

#prevent extrapolation
#depth
grid_pred_longlat[grid_pred_longlat$depth < min(totperstation[!is.na(totperstation$StationID),]$depth),]$depth <- min(totperstation[!is.na(totperstation$StationID),]$depth)
grid_pred_longlat[grid_pred_longlat$depth > max(totperstation[!is.na(totperstation$StationID),]$depth),]$depth <- max(totperstation[!is.na(totperstation$StationID),]$depth)

#mud
grid_pred_longlat[grid_pred_longlat$mud_perc < min(totperstation[!is.na(totperstation$StationID),]$mud_perc),]$mud_perc <- min(totperstation[!is.na(totperstation$StationID),]$mud_perc)
grid_pred_longlat[grid_pred_longlat$mud_perc > max(totperstation[!is.na(totperstation$StationID),]$mud_perc),]$mud_perc <- max(totperstation[!is.na(totperstation$StationID),]$mud_perc)

#combine catch data with prediction data
data_comb <- bind_rows(totperstation, grid_pred_longlat)
```

create projector matrix and inform INLA about the SPDE approach + define spatial field.
```{r}
A <- inla.spde.make.A(mesh, cbind(data_comb$Lon_mid, data_comb$Lat_mid))

dim(A)
table(A[1,]) 
sum(A[1,]) #confirm that the sum is 1

spde <- inla.spde2.matern(mesh, alpha = 2)

w.index <- inla.spde.make.index(
  name = 'w',
  n.spde = spde$n.spde)
```
Define mesh, response variables and covariates
```{r}
data_comb$fSurvey <- factor(data_comb$Survey2)

Xmatrix <- model.matrix(~ Dist + Fraction , data=data_comb)
X <- as.data.frame(Xmatrix[,-1])
names(X) <- c(gsub("[:]",".",names(X)))
head(X)

N <- nrow(data_comb)
Stk_depth_mud_DOY <- inla.stack(tag = "Fit",
                   data = list(y = round(data_comb$Weight)),
                   A = list(1,1,1,1,1,1,A),
                   effects = list(Intercept=rep(1,N),
                     X = X,
                     depth = data_comb$depth,
                     mud = data_comb$mud_perc,
                     DOY = data_comb$DOY_adjusted,
                     fSurvey = data_comb$fSurvey,
                     w = w.index))

#DOY == rw2, to have a more smooth seasonal pattern
f1 <- y ~ -1 + Intercept + offset(log(Dist/Fraction)) + f(w, model = spde) + f(fSurvey, model ="iid") + f(depth, model="rw1") + f(mud, model="rw1") + f(DOY, model="rw2")

M_allfish_depth_mud_DOY <- inla(eval(f1),
           family = "nbinomial",
           data = inla.stack.data(Stk_depth_mud_DOY),
           control.compute = list(dic = TRUE,
                                  waic = TRUE,
                                  config = TRUE,
                                  cpo=T),
           control.predictor = list(A = inla.stack.A(Stk_depth_mud_DOY)), verbose = F)
```

Survey effect
```{r}
#Survey effect
knitr::kable(M_allfish_depth_mud_DOY$summary.random$fSurvey[, c(1:4, 6)], "simple", digits = 2)
```

Random Walk smoothers for depth, mud and DOY
```{r}
#depth
plot_data <- data.frame(
  depth = M_allfish_depth_mud_DOY$summary.random$depth$ID,
  median = M_allfish_depth_mud_DOY$summary.random$depth$`0.5quant`,
  quant0.025 = M_allfish_depth_mud_DOY$summary.random$depth$`0.025quant`,
  quant0.975 = M_allfish_depth_mud_DOY$summary.random$depth$`0.975quant`
)

plot_depth <- ggplot(plot_data, aes(x = depth, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=depth), sides = "b", inherit.aes = F, size=0.1) +
  theme_bw()

#mud
plot_data <- data.frame(
  mud = M_allfish_depth_mud_DOY$summary.random$mud$ID,
  median = M_allfish_depth_mud_DOY$summary.random$mud$`0.5quant`,
  quant0.025 = M_allfish_depth_mud_DOY$summary.random$mud$`0.025quant`,
  quant0.975 = M_allfish_depth_mud_DOY$summary.random$mud$`0.975quant`
)

plot_mud <- ggplot(plot_data, aes(x = mud, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=mud_perc), sides = "b", inherit.aes = F, size=0.1) +
  theme_bw() +
  labs(x = "mud (%)")

#Day of Year
plot_data <- data.frame(
  DOY = M_allfish_depth_mud_DOY$summary.random$DOY$ID,
  median = M_allfish_depth_mud_DOY$summary.random$DOY$`0.5quant`,
  quant0.025 = M_allfish_depth_mud_DOY$summary.random$DOY$`0.025quant`,
  quant0.975 = M_allfish_depth_mud_DOY$summary.random$DOY$`0.975quant`
)

plot_DOY <- ggplot(plot_data, aes(x = DOY, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=DOY_adjusted), sides = "b", inherit.aes = F, size=0.1) +
  theme_bw() +
  labs(x = "day (since October 1st)")

final_plot <- (plot_depth / plot_mud / plot_DOY)

#Fig. S7
print(final_plot)
ggsave(file = file.path(path, "figures/supplement/S7_allspecies_smoothers.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```

spatial field
```{r}
w.pm1 <- M_allfish_depth_mud_DOY$summary.random$w$mean
wproj <- inla.mesh.projector(mesh,  dims = c(1000, 1000))
w.pm100_100 <- inla.mesh.project(wproj, w.pm1)

Grid <- expand.grid(Lon = wproj$x,
                    Lat = wproj$y)
Grid$w.pm <- as.vector(w.pm100_100)
Grid <- Grid[!is.na(Grid$w.pm),]

ggplot() +
  geom_tile(data = Grid, aes(x = Lon, y = Lat, fill = w.pm), show.legend = TRUE) +
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = "Weight") +
  geom_sf(data = doggerbank_sf, color = alpha("white",0.8), size = 3, linetype = "solid", fill=NA) +
  geom_sf(data = NCP_polygon_modified_sf, color = alpha("white",0.8), size = 3, linetype = "solid", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) 
```

#Biomass 
prediction based on the cells from the prediction dataframe
```{r}
#weight observed
totperstation$Weight_m  <- totperstation$Weight/totperstation$Fraction/totperstation$Dist*5 #x5 to get to 1m2 

#weight predicted
pred <- M_allfish_depth_mud_DOY$summary.fitted.values[c((nrow(totperstation)+1):nrow(data_comb)),] 
pred <- cbind(grid_pred_longlat, pred)

pred$g_m <- exp(pred$'0.5quant')/50*5 #dividing by 50 (=0.2m*1m) *5 to 1m2 

#Fig. 3
ggplot() +
  geom_tile(data = pred, aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + 
  scale_fill_gradientn(colors = c("white", "khaki", "green", "yellow", "orange", "red", "darkred", "purple"),
                        values = c(0, 0.1, 0.15, 0.25, 0.45, 0.6, 0.75, 0.9, 1)) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", linewidth = 0.05) +
  labs(x = "Longitude", y = "Latitude", fill = expression(Biomass~(g~m^{-2}))) +
  geom_sf(data = dog40a_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  geom_point(data = totperstation[totperstation$Weight_m > 10, ],
             aes(x = Lon_mid, y = Lat_mid, size = Weight_m),
             color = alpha("black", .2), shape = 1) +
  geom_point(data = totperstation[totperstation$Weight_m < 10, ],
             aes(x = Lon_mid, y = Lat_mid),
             size = 0.6, color = alpha("black", .2), shape = 1) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15)) + theme(text = element_text(size = 7.34)) 

#ggsave(file = file.path(path, "figures/main/F3_allspecies_distribution_alt.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```

Predicted biomass per area of interest. Mean of predictions cells. 95% interval is determined by taking the mean of the 0.025 and 0.975 quantiles. 
```{r}
#DOG
mean(pred[pred$inside_DOG == T,]$g_m)
mean(exp(pred[pred$inside_DOG == T,]$'0.025quant')/10)
mean(exp(pred[pred$inside_DOG == T,]$'0.975quant')/10)
dog_allfish <- mean(pred[pred$inside_DOG == T,]$g_m)

#EEZ
mean(pred[pred$inside_NCP == T,]$g_m)
mean(exp(pred[pred$inside_NCP == T,]$'0.025quant')/10)
mean(exp(pred[pred$inside_NCP == T,]$'0.975quant')/10)
EEZ_allfish <- mean(pred[pred$inside_NCP == T,]$g_m)
```

#Biomass in the DFS and BTS areas 
For the seal prey comparison, but with no species selection.
```{r}
#BTS
BTS_sf <- st_as_sf(BTS, coords = c("Lon", "Lat"))
st_crs(BTS_sf) <- proj_string

pred_sf <- st_as_sf(pred, coords = c("Lon_mid", "Lat_mid"), crs = proj_string)

# Check if points in Grid_sf are within the BTS polygon
points_within_BTS <- st_within(pred_sf, BTS_sf, sparse = FALSE)
pred$BTS_area <- points_within_BTS[, 1]

ggplot() +
  geom_tile(data = pred[pred$BTS_area == T,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = "Biomass (g m⁻²)") + #, title="Small fish (<30cm) biomass"
  geom_sf(data = dog40a_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  # geom_point(data = totperstation[totperstation$Weight_m > 10, ],
  #            aes(x = Lon_mid, y = Lat_mid, size = Weight_m),
  #            color = alpha("black", .5), shape = 1) +
  # geom_point(data = totperstation[totperstation$Weight_m < 10, ],
  #            aes(x = Lon_mid, y = Lat_mid),
  #            size = 0.6, color = alpha("black", .5), shape = 1) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))

mean(pred[pred$BTS_area == T,]$g_m)
mean(exp(pred[pred$BTS_area == T,]$'0.025quant')/10)
mean(exp(pred[pred$BTS_area == T,]$'0.975quant')/10)

#DFS
DFS_sf <- st_as_sf(DFS, coords = c("Lon", "Lat"))
st_crs(DFS_sf) <- proj_string

# Check if points in Grid_sf are within the DFS polygon
points_within_DFS <- st_within(pred_sf, DFS_sf, sparse = FALSE)
pred$DFS_area <- points_within_DFS[, 1]

ggplot() +
  geom_tile(data = pred[pred$DFS_area == T,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = "Biomass (g m⁻²)") + #, title="Small fish (<30cm) biomass"
  geom_sf(data = dog40a_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  geom_point(data = totperstation[totperstation$Weight_m > 10, ],
             aes(x = Lon_mid, y = Lat_mid, size = Weight_m),
             color = alpha("black", .5), shape = 1) +
  geom_point(data = totperstation[totperstation$Weight_m < 10, ],
             aes(x = Lon_mid, y = Lat_mid),
             size = 0.6, color = alpha("black", .5), shape = 1) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))

mean(pred[pred$DFS_area == T,]$g_m)
mean(exp(pred[pred$DFS_area == T,]$'0.025quant')/10)
mean(exp(pred[pred$DFS_area == T,]$'0.975quant')/10)
```

#Residuals
Following Zuur 2017: "Spatial, Temporal and Spatial-Temporal Data Analysis with R-Inla"
```{r}
mu <- M_allfish_depth_mud_DOY$summary.fitted.values[1:1307, "mean"]
E1 <- (totperstation$Weight - mu)/ sqrt(mu)

#Pearson residuals vs fitted
plot(mu, E1, xlab="Fitted values", ylab="Pearson residuals") # not too bad

#per coveriate
plot(totperstation$depth, mu, ylab="Fitted values", xlab="depth") # not too bad
T1 <- gam(E1 ~s(totperstation$depth))
plot(T1)
abline(h=0, lty=2)

plot(totperstation$mud_perc, mu, ylab="Fitted values", xlab="mud") # not too bad
T1 <- gam(E1 ~s(totperstation$mud_perc))
plot(T1)
abline(h=0, lty=2)

#sample variogram
mydata <- data.frame(E1 = E1,
                     Xkm = totperstation$Xkm,
                     Ykm = totperstation$Ykm)

coordinates(mydata) <- c("Xkm", "Ykm")
V1 <- variogram(E1 ~1, mydata, cressie = TRUE)
plot(V1)
```

#Simulations
#Predictive Check
Simulation with the posterior parameters
```{r}
NSim <- 100
SimData <- inla.posterior.sample(n = NSim, result = M_allfish_depth_mud_DOY)
```

```{r}
rownames_before_colon <- sub(":.*", "", rownames(SimData[[1]]$latent))
table(rownames_before_colon)

#APredictor 
predictions <- inla.posterior.sample.eval("APredictor", SimData) #length trawldata?
predictions <- predictions[1:nrow(totperstation),] #sample locations
predictions <- exp(predictions)

#Fig s8
pdf(file = file.path(path, "figures/supplement/S8_posterior_distribution.pdf"), width = 170 / 25.4, height = 170 / 25.4)  # mm to inches
hist(data_comb$Weight, xlim = c(0,2100),ylim =c(0, 0.005),  breaks=  100, freq = F, col=alpha("blue", 0.3), main="", xlab="Biomass (g)") #/10m²
hist(predictions, xlim = c(0,2100), breaks= 100, freq = F, col=alpha("red", 0.4), add=T)
dev.off()
```

fitted values
```{r}
fitted <- M_allfish_depth_mud_DOY$summary.fitted.values[c(1:nrow(totperstation)),"0.5quant"]
combined_fitted <- cbind(data_comb[c(1:nrow(totperstation)),],fitted)

ggplot(combined_fitted , aes(x = Weight, y = fitted)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  labs(title = "Observed vs. Fitted Values", x = "observed values", y = "fitted values") +
  theme_bw()
```

#per species
Biomass predictions are made for the 7 most common species. The method is similar to the previous combined total biomass.

Create English name of the 8 most common fish species. 
```{r}
fishdata$Species_English <- "Remaining"
fishdata[fishdata$Species_reported == "Gobiidae",]$Species_English <- "Goby spp."
fishdata[fishdata$Species_reported == "Ammodytidae",]$Species_English <- "Sandeel spp."
fishdata[fishdata$Species_reported == "Buglossidium luteum",]$Species_English <- "Solenette"
fishdata[fishdata$Species_reported == "Limanda limanda",]$Species_English <- "Dab"
fishdata[fishdata$Species_reported == "Arnoglossus laterna",]$Species_English <- "Scaldfish"
fishdata[fishdata$Species_reported == "Callionymus sp.",]$Species_English <- "Dragonet spp."
fishdata[fishdata$Species_reported == "Pleuronectes platessa",]$Species_English <- "Plaice"
```

Weight per station
```{r}
interested_species <- unique(fishdata$Species_English) 

for (i in interested_species) {
print(i)
fishdata_specific <- fishdata[fishdata$Species_English == i,]
stations_notinfishdata_specific <- unique(fishdata[!fishdata$StationID %in% fishdata_specific$StationID,]$StationID)

#multiple fractions within a species
multipleFractions <- as.data.frame(table(fishdata_specific$StationID, fishdata_specific$Fraction))
colnames(multipleFractions) <- c("StationID","Fraction", "Freq")

multipleFractions[multipleFractions$Freq > 0,]$Freq <- 1 #to create only present value
multipleFractions <- aggregate(Freq~StationID, multipleFractions, FUN="sum")
multipleFractions <- multipleFractions[multipleFractions$Freq > 1,]

for (j in unique(multipleFractions$StationID)) {
min_frac <- min(fishdata_specific[fishdata_specific$StationID == j,]$Fraction, na.rm=T)
max_frac <- max(fishdata_specific[fishdata_specific$StationID == j,]$Fraction, na.rm=T)  
multiplier <- max_frac/min_frac

fishdata_specific[fishdata_specific$StationID == j & fishdata_specific$Fraction == min_frac,]$Weight <- fishdata_specific[fishdata_specific$StationID == j & fishdata_specific$Fraction == min_frac,]$Weight * multiplier 

fishdata_specific[fishdata_specific$StationID == j & fishdata_specific$Fraction == min_frac,]$Fraction <- fishdata_specific[fishdata_specific$StationID == j & fishdata_specific$Fraction == min_frac,]$Fraction * multiplier 
}

#aggregate
totperstation <- aggregate(Weight ~ StationID + Species_English, fishdata_specific, sum, drop = FALSE)

#add fraction information
totperstation <- merge(totperstation, fishdata_specific[!duplicated(fishdata_specific$StationID),c("StationID", "Fraction")], by= "StationID", all.x = T )

#the samples without interested species
totperstation1 <- fishdata[fishdata$StationID %in% stations_notinfishdata_specific,c("StationID", "Fraction")]
totperstation1 <- aggregate(Fraction~StationID, totperstation1, max)
totperstation1$Species_English <- i
totperstation1$Weight <- 0
totperstation <- rbind(totperstation, totperstation1)

#add trawl information
totperstation <- merge(totperstation, trawldata[,c("StationID", "Survey2", "Dist", "Lon_mid", "Lat_mid","depth","mud_perc","DOY_adjusted")], by= "StationID", all.x = T )

#exclude Survey which goes to the north part of north sea 
totperstation <- totperstation[totperstation$Lat_mid < 56,]
totperstation <- totperstation[totperstation$Fraction > 0.1,]

#weight observed
totperstation$Weight_m  <- totperstation$Weight/totperstation$Fraction/totperstation$Dist*5 #x5 to get to 1m2 
assign(paste0("observed_", i),totperstation)

#model
totperstation$fSurvey <- factor(totperstation$Survey2)

#prediction field
data_comb <- bind_rows(totperstation, grid_pred_longlat[,!names(grid_pred_longlat) %in% c("inside_NCP","inside_DOG")])

Xmatrix <- model.matrix(~ Dist + Fraction , data=data_comb)
X <- as.data.frame(Xmatrix[,-1])
names(X) <- c(gsub("[:]",".",names(X)))
head(X)

N <- nrow(data_comb)

Stk_depth_mud_DOY <- inla.stack(tag = "Fit",
                   data = list(y = round(data_comb$Weight)),
                   A = list(1,1,1,1,1,1,A),
                   effects = list(Intercept=rep(1,N),
                     X = X,
                     depth = data_comb$depth,
                     mud = data_comb$mud_perc,
                     fSurvey = data_comb$fSurvey,
                     DOY = data_comb$DOY_adjusted,
                     w = w.index))

f1 <- y ~ -1 + Intercept + offset(log(Dist/Fraction)) + f(fSurvey, model ="iid") + f(w, model = spde) + f(depth, model="rw1") + f(mud, model="rw1") + f(DOY, model="rw2")
model_name <- gsub(" ", "_", i)
 assign(paste0("M_", model_name, "_depth_mud_DOY"),
         inla(eval(f1),
              family = "nbinomial",
              data = inla.stack.data(Stk_depth_mud_DOY),
              control.compute = list(dic = TRUE, waic = TRUE, cpo=T, config = TRUE),
              control.predictor = list(A = inla.stack.A(Stk_depth_mud_DOY)), verbose = F))
  print(paste("Model depth_mud for", i, "done"))
}
```

simulations
```{r}
interested_species <- sort(unique(fishdata$Species_English))

predictions_list <- list()

# Loop over each species
for (species in interested_species) {
  # Replace spaces with underscores in species names
  species_name <- gsub(" ", "_", species)

  model_name <- paste0("M_", species_name, "_Remaining_depth_mud_DOY")
  model_depth_mud <- get(model_name)
    
  # Set the number of simulations
  NSim <- 100
    
  # Run posterior samples
  SimData <- inla.posterior.sample(n = NSim, result = model_depth_mud)
    
  predictions <- inla.posterior.sample.eval("APredictor", SimData)
  predictions <- predictions[1:nrow(totperstation), ]
  predictions <- exp(predictions)
    
  # Store predictions in the list with the species name as the key
  predictions_list[[species]] <- predictions
}
```

Simulated versus observed
```{r}
#dab
hist(observed_Dab$Weight, xlim = c(0,100), ylim=c(0,800), breaks=  seq(0, max(observed_Dab$Weight+1, na.rm = TRUE), by = 1), freq = T, col=alpha("blue", 0.3), main="dab", xlab="Biomass (g)")
hist(predictions_list$Dab, xlim = c(0,100), breaks= seq(0, max(predictions_list$Dab+1, na.rm = TRUE), by = 1), freq = T, col=alpha("red", 0.4), add=T)

#dragonet
hist(`observed_Dragonet spp.`$Weight, xlim = c(0,100), ylim=c(0,800), breaks=  seq(0, max(`observed_Dragonet spp.`$Weight+1, na.rm = TRUE), by = 1), freq = T, col=alpha("blue", 0.3), main="dragonet", xlab="Biomass (g)")
hist(predictions_list$`Dragonet spp.`, xlim = c(0,100), breaks= seq(0, max(predictions_list$`Dragonet spp.`+1, na.rm = TRUE), by = 1), freq = T, col=alpha("red", 0.4), add=T)

#goby
hist(`observed_Goby spp.`$Weight, xlim = c(0,100), ylim=c(0,450), breaks=  seq(0, max(`observed_Goby spp.`$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="goby", xlab="Biomass (g)")
hist(predictions_list$`Goby spp.`, xlim = c(0,100), breaks= seq(0, max(predictions_list$`Goby spp.`, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)

#plaice
hist(observed_Plaice$Weight, xlim = c(0,100), ylim=c(0,800), breaks=seq(0, max(observed_Plaice$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="plaice", xlab="Biomass (g)")
hist(predictions_list$Plaice, xlim = c(0,100), breaks= seq(0, max(predictions_list$Plaice, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)

#sandeel
hist(`observed_Sandeel spp.`$Weight, xlim = c(0,100), ylim=c(0,800), breaks=  seq(0, max(`observed_Sandeel spp.`$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="sandeel", xlab="Biomass (g)")
hist(predictions_list$Sandeel, xlim = c(0,100), breaks= seq(0, max(predictions_list$`Sandeel spp.`, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)

#scaldfish
hist(observed_Scaldfish$Weight, xlim = c(0,50), ylim=c(0,800), breaks=  seq(0, max(observed_Scaldfish$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="scaldfish", xlab="Biomass (g)")
hist(predictions_list$Scaldfish, xlim = c(0,50), breaks= seq(0, max(predictions_list$Scaldfish, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)

#solenette
hist(observed_Solenette$Weight, xlim = c(0,100), ylim=c(0,800), breaks=  seq(0, max(observed_Solenette$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="solenette", xlab="Biomass (g)")
hist(predictions_list$Solenette, xlim = c(0,100), breaks= seq(0, max(predictions_list$Solenette, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)

#whiting
hist(observed_Whiting$Weight, xlim = c(0,100), ylim=c(0,1100), breaks=  seq(0, max(observed_Whiting$Weight, na.rm = TRUE)+1,1), freq = T, col=alpha("blue", 0.3), main="whiting", xlab="Biomass (g)")
hist(predictions_list$Whiting, xlim = c(0,100), breaks= seq(0, max(predictions_list$Whiting, na.rm = TRUE)+1,1), freq = T, col=alpha("red", 0.4), add=T)
```

```{r}
plots <- list()  

for (i in c("Dab", "Dragonet spp.", "Goby spp.", "Plaice", "Sandeel spp.", "Scaldfish", "Solenette", "Remaining")) {
print(i)
    species_name <- gsub(" ", "_", i)

  # Depth plot
  model_depth_mud <- get(paste0("M_", species_name, "_depth_mud_DOY")) 
  plot_data_depth_mud <- data.frame(
    depth = model_depth_mud$summary.random$depth$ID,
    median = model_depth_mud$summary.random$depth$`0.5quant`,
    quant0.025 = model_depth_mud$summary.random$depth$`0.025quant`,
    quant0.975 = model_depth_mud$summary.random$depth$`0.975quant`
  )

  p_depth <- ggplot(plot_data_depth_mud, aes(x = depth, y = median, ymin = quant0.025, ymax = quant0.975)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(data = trawldata, aes(x = depth), sides = "b", inherit.aes = F, size= 0.1) +
    theme_bw() + geom_hline(yintercept = 0, linetype = "dashed")

  p_depth <- p_depth + ggtitle(i) +
    theme(plot.title = element_text(hjust = 0, vjust = 1, face = "bold", size = 8)) + theme(text = element_text(size = 7.34)) 

  if (!i  %in% c("Remaining")) { 
    p_depth <- p_depth + theme(
        axis.title.x = element_blank(),  # Remove x-axis title
        axis.text.x = element_blank(),   # Remove x-axis tick labels
        axis.ticks.x = element_blank()   # Remove x-axis ticks
    )
}
  
  # Mud plot
  plot_data_mud <- data.frame(
    mud = model_depth_mud$summary.random$mud$ID,
    median = model_depth_mud$summary.random$mud$`0.5quant`,
    quant0.025 = model_depth_mud$summary.random$mud$`0.025quant`,
    quant0.975 = model_depth_mud$summary.random$mud$`0.975quant`
  )

  p_mud <- ggplot(plot_data_mud, aes(x = mud, y = median, ymin = quant0.025, ymax = quant0.975)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(data = trawldata, aes(x = mud_perc), sides = "b", inherit.aes = F, size= 0.1) +
    theme_bw() + geom_hline(yintercept = 0, linetype = "dashed") + theme(text = element_text(size = 7.34)) 
   
  if (!i %in% c("Remaining")){ 
    p_mud <- p_mud + theme(
        axis.title.x = element_blank(),  # Remove x-axis title
        axis.text.x = element_blank(),   # Remove x-axis tick labels
        axis.ticks.x = element_blank()   # Remove x-axis ticks
    )
  }
  
  # DOY plot
  plot_data_DOY <- data.frame(
    DOY = model_depth_mud$summary.random$DOY$ID,
    median = model_depth_mud$summary.random$DOY$`0.5quant`,
    quant0.025 = model_depth_mud$summary.random$DOY$`0.025quant`,
    quant0.975 = model_depth_mud$summary.random$DOY$`0.975quant`
  )

  p_DOY <- ggplot(plot_data_DOY, aes(x = DOY, y = median, ymin = quant0.025, ymax = quant0.975)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(data = trawldata, aes(x = DOY_adjusted), sides = "b", inherit.aes = F, size= 0.1) + #Day adjusted
    theme_bw() + geom_hline(yintercept = 0, linetype = "dashed") + labs(x="day since 1st of October")  + theme(text = element_text(size = 7.34)) 
   
  if (i != "Remaining") {
    p_DOY <- p_DOY + theme(
        axis.title.x = element_blank(),  # Remove x-axis title
        axis.text.x = element_blank(),   # Remove x-axis tick labels
        axis.ticks.x = element_blank()   # Remove x-axis ticks
    )
}
  
p_depth <- p_depth +  scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
p_mud <- p_mud + scale_y_continuous(labels = scales::number_format(accuracy = 0.1))
p_DOY <- p_DOY + scale_y_continuous(labels = scales::number_format(accuracy = 0.1))

  
  # Store plots
  plots[[length(plots) + 1]] <- p_depth
  plots[[length(plots) + 1]] <- p_mud
  plots[[length(plots) + 1]] <- p_DOY
}

# Arrange all four plots in a 2x2 layout with shared titles
final_plot <- (plots[[1]] | plots[[2]] | plots[[3]]) / (plots[[4]] | plots[[5]] | plots[[6]]) / (plots[[7]] | plots[[8]] | plots[[9]])/ (plots[[10]] | plots[[11]] | plots[[12]])/ (plots[[13]] | plots[[14]] | plots[[15]])/ (plots[[16]] | plots[[17]] | plots[[18]])/ (plots[[19]] | plots[[20]] | plots[[21]])/ (plots[[22]] | plots[[23]] | plots[[24]])

#Fig S9
print(final_plot)
#ggsave(file = file.path(path, "figures/supplement/S9_perspecies_smoothers.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```


```{r}
interested_species <- unique(fishdata$Species_English) 
for (i in interested_species) {

model <- get(paste0("M_", gsub(" ", "_", i)  ,"_depth_mud_DOY"))  
pred <- model$summary.fitted.values[c((nrow(totperstation)+1):nrow(data_comb)),] 
pred <- cbind(grid_pred_longlat, pred)

pred$g_m <- exp(pred$'0.5quant')/50*5 #dividing by 50 (=0.2m*1m) *5 to 1m2 

#observations per species
totperstation_specific <- get(paste0("observed_", i))

p <- ggplot() +
  geom_tile(data = pred[pred$g_m < 100,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = "Biomass (g m⁻²)", title = i) +
  #geom_sf(data = dog40a_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  #geom_sf(data = dog40b_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  geom_point(data = totperstation_specific[totperstation_specific$Weight_m > 1, ],
             aes(x = Lon_mid, y = Lat_mid, size = Weight_m),
             color = alpha("black", .5), shape = 1) +
  geom_point(data = totperstation_specific[totperstation_specific$Weight_m < 1, ],
             aes(x = Lon_mid, y = Lat_mid),
             size = 0.6, color = alpha("black", .5), shape = 1) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))

plot(p)

#separate dataframes for pred
pred$Species_English <- i
assign(paste0("pred_", i), pred)
}

#combine dataframes
combined_pred <- data.frame()
for (i in interested_species) {
  df <- get(paste0("pred_", i))

  if (nrow(combined_pred) == 0) {
    combined_pred <- df
  } else {
    combined_pred <- rbind(combined_pred, df)
  }
}

observed_allspecies <- data.frame()
for (i in interested_species) {
  df <- get(paste0("observed_", i))

  if (nrow(observed_allspecies) == 0) {
    observed_allspecies <- df
  } else {
    observed_allspecies <- rbind(observed_allspecies, df)
  }
}
```


```{r}
# Define the desired breaks and their corresponding log-transformed values
desired_breaks <- c(0.001, 1, 4, 9, 16, 25, 36, 49, 64, 81)
desired_breaks_lab <- c(0, 1, 4, 9, 16, 25, 36, 49, 64, 81)

desired_order <- c("Dab", "Dragonet spp.", "Goby spp.", "Plaice", "Sandeel spp.", "Scaldfish", "Solenette", "Remaining")

# Apply the factor levels in this order
combined_pred$Species_English <- factor(combined_pred$Species_English, levels = desired_order)
observed_allspecies$Species_English <- factor(observed_allspecies$Species_English, levels = desired_order)

#Fig 4
p <- ggplot() +
  geom_tile(data = combined_pred, aes(x = Lon_mid, y = Lat_mid, fill = g_m, width = 0.05, height = 0.05)) + #
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  # Facet by species
  facet_wrap(~ Species_English, nrow = 4) +
  scale_fill_gradientn(
    colors = c("white", "beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple"),
    values = c(0, 0.03, 0.08, 0.15, 0.25, 0.4, 0.7, 0.9, 1.0),   
    name = expression("Biomass (g m"^{-2}*")"), 
    breaks = desired_breaks,    
    labels = desired_breaks_lab, 
    trans = 'sqrt'
  ) +
  #Add species-specific observations
  geom_point(data = observed_allspecies[observed_allspecies$Weight_m > 0, ],
             aes(x = Lon_mid, y = Lat_mid, size = log(Weight_m  / 15 + 0.6)),
             color = alpha("black", .1), shape = 1, show.legend = F) +
  geom_point(data = observed_allspecies[observed_allspecies$Weight_m  == 0, ],
             aes(x = Lon_mid, y = Lat_mid),
             size = 0.6, color = alpha("black", .2), shape = 16, show.legend = F) +
  coord_sf(xlim = xlim, ylim = ylim) +  
  labs(x = "Longitude", y = "Latitude") +
  scale_size_continuous(range = c(0.5, 8)) + # Adjust the range to make size differences more noticeable
  theme_bw() + theme(text = element_text(size = 7.34)) 

plot(p)
#ggsave(file = file.path(path, "figures/main/F4_perspecies_distribution.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```

mean biomass per area per species
```{r}
mean_results <- data.frame(
  species = character(),
  area = character(),
  mean_g_m = numeric(),
  mean_025quant = numeric(),
  mean_975quant = numeric(),
  stringsAsFactors = FALSE
)

for (i in unique(combined_pred$Species_English)) {

#DOG
mean_g_m_DOG <- mean(combined_pred[combined_pred$Species_English == i & combined_pred$inside_DOG == TRUE, ]$g_m, na.rm = TRUE)
mean_025quant_DOG <- mean(exp(combined_pred[combined_pred$Species_English == i & combined_pred$inside_DOG == TRUE, ]$'0.025quant') / 10, na.rm = TRUE)
mean_975quant_DOG <- mean(exp(combined_pred[combined_pred$Species_English == i & combined_pred$inside_DOG == TRUE, ]$'0.975quant') / 10, na.rm = TRUE)
  
# Append DOG results
mean_results <- rbind(mean_results, data.frame(
  species = i,
  area = "DOG",
  mean_g_m = mean_g_m_DOG,
  mean_025quant = mean_025quant_DOG,
  mean_975quant = mean_975quant_DOG
  ))

# EEZ Area
mean_g_m_NCP <- mean(combined_pred[combined_pred$Species_English == i & combined_pred$inside_NCP == TRUE, ]$g_m, na.rm = TRUE)
mean_025quant_NCP <- mean(exp(combined_pred[combined_pred$Species_English == i & combined_pred$inside_NCP == TRUE, ]$'0.025quant') / 10, na.rm = TRUE)
mean_975quant_NCP <- mean(exp(combined_pred[combined_pred$Species_English == i & combined_pred$inside_NCP == TRUE, ]$'0.975quant') / 10, na.rm = TRUE)
  
# Append NCP results
mean_results <- rbind(mean_results, data.frame(
  species = i,
  area = "NCP",
  mean_g_m = mean_g_m_NCP,
  mean_025quant = mean_025quant_NCP,
  mean_975quant = mean_975quant_NCP
))
  
}

mean_results[, sapply(mean_results, is.numeric)] <- round(mean_results[, sapply(mean_results, is.numeric)], 1)

#contribution to all fish biomass
mean_results$contribution <- NA
mean_results[mean_results$area == "DOG",]$contribution <- round(mean_results[mean_results$area == "DOG",]$mean_g_m/dog_allfish*100,1)
mean_results[mean_results$area == "NCP",]$contribution <- round(mean_results[mean_results$area == "NCP",]$mean_g_m/EEZ_allfish*100,1)
```


###############################################################################################################################
#Dutch waddencoast, seal prey
In Aarts et al. (2019), prey seal biomass was estimated based on the BTS and DFS Surveys. For the same area and same species estimates are made based on the Triple-D data 

We grouped in species (e.g dragonets), however in the study of Aarts 2019 some specific (prey)species were selected. Therefore those species are again put into species level.  
```{r}
subset_Callionymuslyra <- fishdata[grepl('Callionymus lyra' , fishdata$unique_code), ]
fishdata[fishdata$unique_code %in% subset_Callionymuslyra$unique_code,]$Species_reported <- "Callionymus lyra"

subset_Ciliatamustel <- fishdata[grepl('Ciliata mustela' , fishdata$unique_code), ]
fishdata[fishdata$unique_code %in% subset_Ciliatamustel$unique_code,]$Species_reported <- "Ciliata mustela"

subset_Myoxocephalusscorpius <- fishdata[grepl('Myoxocephalus scorpius' , fishdata$unique_code), ]
fishdata[fishdata$unique_code %in% subset_Myoxocephalusscorpius$unique_code,]$Species_reported <- "Myoxocephalus scorpius"

subset_ammodytes <- fishdata[!grepl('Hyperoplus lanceolatus' , fishdata$unique_code) & fishdata$Species_reported == "Ammodytidae", ]
fishdata[fishdata$unique_code %in% subset_ammodytes$unique_code,]$Species_reported <- "Ammodytes sp."
```

Weight per station in which only the 9 selected prey species are included
```{r}
#choose area
#in the DFS Survey, sandeel is included
#Survey <- "DFS"
#Survey <- "DFS_without_Sandeel"
Survey <- "BTS"

species_seal <- c( "Callionymus lyra", "Ciliata mustela", "Gadus morhua", "Limanda limanda", "Merlangius merlangus", "Myoxocephalus scorpius", "Platichthys flesus",  "Pleuronectes platessa", "Solea solea")

if(Survey == "DFS"){
  species_seal <- c("Ammodytes sp.", "Callionymus lyra", "Ciliata mustela", "Gadus morhua", "Limanda limanda", "Merlangius merlangus", "Myoxocephalus scorpius", "Platichthys flesus",  "Pleuronectes platessa", "Solea solea")
}

if(Survey == "DFS_without_Sandeel"){
  species_seal <- c("Callionymus lyra", "Ciliata mustela", "Gadus morhua", "Limanda limanda", "Merlangius merlangus", "Myoxocephalus scorpius", "Platichthys flesus",  "Pleuronectes platessa", "Solea solea")
}

#species selection
preydata <- fishdata[fishdata$Species_reported %in% species_seal,]
stations_notinpreyhdata <- unique(fishdata[!fishdata$StationID %in% preydata$StationID,]$StationID)

#multiple fractions within a species
multipleFractions <- as.data.frame(table(preydata$StationID, preydata$Fraction))
colnames(multipleFractions) <- c("StationID","Fraction", "Freq")

multipleFractions[multipleFractions$Freq > 0,]$Freq <- 1 #to create only present value
multipleFractions <- aggregate(Freq~StationID, multipleFractions, FUN="sum")
multipleFractions <- multipleFractions[multipleFractions$Freq > 1,]

for (j in unique(multipleFractions$StationID)) {
min_frac <- min(preydata[preydata$StationID == j,]$Fraction, na.rm=T)
max_frac <- max(preydata[preydata$StationID == j,]$Fraction, na.rm=T)  
multiplier <- max_frac/min_frac

preydata[preydata$StationID == j & preydata$Fraction == min_frac,]$Weight <- preydata[preydata$StationID == j & preydata$Fraction == min_frac,]$Weight * multiplier 

preydata[preydata$StationID == j & preydata$Fraction == min_frac,]$Fraction <- preydata[preydata$StationID == j & preydata$Fraction == min_frac,]$Fraction * multiplier 
}

totperstation <- aggregate(Weight ~ StationID, preydata, sum, drop = FALSE)

#add fraction information
totperstation <- merge(totperstation, preydata[!duplicated(preydata$StationID),c("StationID", "Fraction")], by= "StationID", all.x = T )

#the samples without interested species
totperstation1 <- fishdata[fishdata$StationID %in% stations_notinpreyhdata,c("StationID", "Fraction")]
totperstation1 <- aggregate(Fraction~StationID, totperstation1, max)
totperstation1$Weight <- 0
totperstation <- rbind(totperstation, totperstation1)

#add trawl information
totperstation <- merge(totperstation, trawldata[,c("StationID", "Survey2", "Dist", "Lon_mid", "Lat_mid", "depth", "mud_perc", "DOY_adjusted")], by= "StationID", all.x = T )

#exclude Survey which goes to the north part of north sea 
totperstation <- totperstation[totperstation$Lat_mid < 56,]

#prediction datafield
data_comb <- bind_rows(totperstation, grid_pred_longlat[,!names(grid_pred_longlat) %in% c("inside_NCP","inside_DOG")])
```

define response variables and covariates. Fit model.
```{r}
data_comb$fSurvey <- factor(data_comb$Survey2)

Xmatrix <- model.matrix(~ Dist + Fraction , data=data_comb) 
X <- as.data.frame(Xmatrix[,-1])
names(X) <- c(gsub("[:]",".",names(X)))
head(X)

N <- nrow(data_comb)
Stk_depth_mud_DOY <- inla.stack(tag = "Fit",
                   data = list(y = round(data_comb$Weight)),
                   A = list(1,1,1,1,1,1,A),
                   effects = list(Intercept=rep(1,N),
                     X = X,
                     depth = data_comb$depth,
                     mud = data_comb$mud_perc,
                     DOY = data_comb$DOY_adjusted,
                     fSurvey = data_comb$fSurvey,
                     w = w.index))

#DOY == rw2, to have a more smooth annual pattern
f1 <- y ~ -1 + Intercept + offset(log(Dist/Fraction)) + f(w, model = spde) + f(fSurvey, model ="iid") + f(depth, model="rw1") + f(mud, model="rw1") + f(DOY, model="rw2") 

M_seals <- inla(eval(f1),
           family = "nbinomial",
           data = inla.stack.data(Stk_depth_mud_DOY),
           control.compute = list(dic = TRUE,
                                  waic = TRUE,
                                  config = TRUE,
                                  cpo=T),
           control.predictor = list(A = inla.stack.A(Stk_depth_mud_DOY)), verbose = F)

knitr::kable(M_seals$summary.random$fSurvey[, c(1:4, 6)], "simple", digits = 2)
```

prediction
```{r}
pred <- M_seals$summary.fitted.values[c((nrow(totperstation)+1):nrow(data_comb)),] 
pred <- cbind(grid_pred_longlat, pred)

pred$g_m <- exp(pred$'0.5quant')/50*5 #dividing by 50 (=0.2m*1m) *5 to 1m2 
pred_sf <- st_as_sf(pred, coords = c("Lon_mid", "Lat_mid"), crs = proj_string)

if (Survey == "BTS"){
#BTS

  # Check if points in Grid_sf are within the BTS polygon
points_within_BTS <- st_within(pred_sf, BTS_sf, sparse = FALSE)
pred$BTS_area <- points_within_BTS[, 1]

p <- ggplot() +
  geom_tile(data = pred[pred$BTS_area == T,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +  theme(text = element_text(size = 7.34)) +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = expression(Biomass~(g~m^{-2}))) + 
  geom_sf(data = dog40a_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))

print(p)
#Fig S11_BTS
ggsave(file = file.path(path, "figures/supplement/S11_BTS.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)

print(mean(pred[pred$BTS_area == T,]$g_m))
print(mean(exp(pred[pred$BTS_area == T,]$'0.025quant')/10))
print(mean(exp(pred[pred$BTS_area == T,]$'0.975quant')/10))}

if (Survey == "DFS"){
#DFS
DFS_sf <- st_as_sf(DFS, coords = c("Lon", "Lat"))
st_crs(DFS_sf) <- proj_string

# Check if points in Grid_sf are within the DFS polygon
points_within_DFS <- st_within(pred_sf, DFS_sf, sparse = FALSE)
pred$DFS_area <- points_within_DFS[, 1]

p <- ggplot() +
  geom_tile(data = pred[pred$DFS_area == T,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +  theme(text = element_text(size = 7.34)) +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = expression(Biomass~(g~m^{-2}))) + 
  geom_sf(data = dog40a_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))
print(p)

#Fig S11_DFS
ggsave(file = file.path(path, "figures/supplement/S11_DFS.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)

print(mean(pred[pred$DFS_area == T,]$g_m))
print(mean(exp(pred[pred$DFS_area == T,]$'0.025quant')/10))
print(mean(exp(pred[pred$DFS_area == T,]$'0.975quant')/10))}

if (Survey == "DFS_without_Sandeel"){
#DFS
# Check if points in Grid_sf are within the DFS polygon
points_within_DFS <- st_within(pred_sf, DFS_sf, sparse = FALSE)
pred$DFS_area <- points_within_DFS[, 1]

p <- ggplot() +
  geom_tile(data = pred[pred$DFS_area == T,], aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = "Biomass (g m⁻²)") + #, title="Small fish (<30cm) biomass"
  geom_sf(data = dog40a_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), size = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  guides(size = "none") +
  scale_size_continuous(range = c(0.5, 15))
print(p)

print(mean(pred[pred$DFS_area == T,]$g_m))
print(mean(exp(pred[pred$DFS_area == T,]$'0.025quant')/10))
print(mean(exp(pred[pred$DFS_area == T,]$'0.975quant')/10))}
```

###################################################################################################################################################
#Greenstreet small fish (8-31g)
Greenstreet et al.  (2007) estimated the small bottom fish (8-31g) biomass based on IBTS data from 1998 to 2004. Here we compare it with the biomass estimate based on the Triple-D. Sandeel is excluded in their study (considered as a pelagic species), so we also exclude it also from our estimate. 

An individual weight is needed to be able to select for fish which are 8-31 gram.
```{r}
fishdata$Weight_individual <- NA

#Entry
fishdata[fishdata$Weight_type == "Entry" & !is.na(fishdata$Weight_type) & fishdata$is_Partial_WW == 0,]$Weight_individual <- fishdata[fishdata$Weight_type == "Entry" & !is.na(fishdata$Weight_type)& fishdata$is_Partial_WW == 0,]$WW_g/fishdata[fishdata$Weight_type == "Entry" & !is.na(fishdata$Weight_type)  & fishdata$is_Partial_WW == 0,]$Count
summary(fishdata$Weight_individual)# 24813 NA's

#Sample with an expected_WW_g
fishdata[is.na(fishdata$Weight_individual),]$Weight_individual  <- fishdata[is.na(fishdata$Weight_individual),]$expected_WW_g/fishdata[is.na(fishdata$Weight_individual),]$Count
summary(fishdata$Weight_individual)# 9732 NA's

#Predicted individual
fishdata[is.na(fishdata$Weight_individual) & fishdata$is_Partial_WW == 0,]$Weight_individual <- fishdata[is.na(fishdata$Weight_individual) & fishdata$is_Partial_WW == 0,]$WW_g_pred_individual
summary(fishdata$Weight_individual)# 858 NA's
```

damaged fish which are within the selection range are included on base of their length (if length is available)
```{r}
#length-weight conversion
conversion <- read.csv(file.path(path, "/otherdata/bioconversion.csv"))

#what's left
left <- fishdata[is.na(fishdata$Weight_individual),]
left$weight_selection_in <- NA

exluded <- c("Hippoglossoides platessoides", "Myxine glutinosa", "Pholis gunnellus", "Liparis liparis" ) #outside selection range
for( i in unique(left[!is.na(left$length_cm) & !left$Species_reported %in% exluded ,]$Species_reported)){
print(i)
if (i == "Pleuronectiformes"){ i <- "Pleuronectes platessa"}
temp <- data.frame("length_cm" = sort(unique(left[left$Species_reported == i & !is.na(left$length_cm),]$length_cm)))
temp$weigth_g <- round(conversion[conversion$Soort == i,]$A_factor * (10*temp$length_cm)^conversion[conversion$Soort == i,]$B_exponent,1)
temp <- temp[temp$weigth_g >= 8 & temp$weigth_g <= 31,]
left$weight_selection_in[left$Species_reported == i & left$length_cm %in% temp$length_cm] <- "YES"
}
left <- left[left$weight_selection_in == "YES" & !is.na(left$weight_selection_in),]
```

selection
```{r}
smallfish <- fishdata[fishdata$Weight_individual >= 8 & fishdata$Weight_individual <= 31 & !is.na(fishdata$Weight_individual),]
smallfish <- smallfish[!is.na(smallfish$StationID),]

#include damaged but weight within selection range
smallfishtest <- rbind(smallfish, left[,!names(left)  %in% "weight_selection_in"])

withoutsandeel <- "No"
if(withoutsandeel == "Yes"){
smallfish <- smallfish[smallfish$Species_reported != "Ammodytidae",]
}
table(smallfish$is_Partial_WW)

#we need individual weights of the Sample data
smallfish[smallfish$Weight_type == "Sample" & !is.na(smallfish$Weight_type),]$Weight <- smallfish[smallfish$Weight_type == "Sample" & !is.na(smallfish$Weight_type),]$expected_WW_g
```

```{r}
stations_notinsmallfishdata <- unique(fishdata[!fishdata$StationID %in% smallfish$StationID,]$StationID)

#multiple fractions within a species
multipleFractions <- as.data.frame(table(smallfish$StationID, smallfish$Fraction))
colnames(multipleFractions) <- c("StationID","Fraction", "Freq")

multipleFractions[multipleFractions$Freq > 0,]$Freq <- 1 #to create only present value
multipleFractions <- aggregate(Freq~StationID, multipleFractions, FUN="sum")
multipleFractions <- multipleFractions[multipleFractions$Freq > 1,]

for (j in unique(multipleFractions$StationID)) {
min_frac <- min(smallfish[smallfish$StationID == j,]$Fraction, na.rm=T)
max_frac <- max(smallfish[smallfish$StationID == j,]$Fraction, na.rm=T)  
multiplier <- max_frac/min_frac

smallfish[smallfish$StationID == j & smallfish$Fraction == min_frac,]$Weight <- smallfish[smallfish$StationID == j & smallfish$Fraction == min_frac,]$Weight * multiplier 

smallfish[smallfish$StationID == j & smallfish$Fraction == min_frac,]$Fraction <- smallfish[smallfish$StationID == j & smallfish$Fraction == min_frac,]$Fraction * multiplier 
}

totperstation <- aggregate(Weight ~ StationID, smallfish, sum, drop = FALSE)

#add fraction information
totperstation <- merge(totperstation, smallfish[!duplicated(smallfish$StationID),c("StationID", "Fraction")], by= "StationID", all.x = T )

#the samples without interested species
totperstation1 <- fishdata[fishdata$StationID %in% stations_notinsmallfishdata,c("StationID", "Fraction")]
totperstation1 <- aggregate(Fraction~StationID, totperstation1, max)
totperstation1$Weight <- 0
totperstation <- rbind(totperstation, totperstation1)

#add trawl information
totperstation <- merge(totperstation, trawldata[,c("StationID", "Survey2", "Dist", "Lon_mid", "Lat_mid", "depth", "mud_perc", "DOY_adjusted")], by= "StationID", all.x = T )

#exclude Survey which goes to the north part of north sea 
totperstation <- totperstation[totperstation$Lat_mid < 56,]

#prediction datafield
data_comb <- bind_rows(totperstation, grid_pred_longlat[,!names(grid_pred_longlat) %in% c("inside_NCP","inside_DOG")])
```

define response variables and covariates for INLA
```{r}
data_comb$fSurvey <- factor(data_comb$Survey2)

Xmatrix <- model.matrix(~ Dist + Fraction , data=data_comb) 
X <- as.data.frame(Xmatrix[,-1])
names(X) <- c(gsub("[:]",".",names(X)))
head(X)

N <- nrow(data_comb)
Stk_depth_mud_DOY<- inla.stack(tag = "Fit",
                   data = list(y = round(data_comb$Weight)),
                   A = list(1,1,1,1,1,1,A),
                   effects = list(Intercept=rep(1,N),
                     X = X,
                     depth = data_comb$depth,
                     mud = data_comb$mud_perc,
                     DOY = data_comb$DOY_adjusted,
                     fSurvey = data_comb$fSurvey,
                     w = w.index))

#DOY == rw2, to have a more smooth annual pattern
f1 <- y ~ -1 + Intercept + offset(log(Dist/Fraction)) + f(w, model = spde) + f(fSurvey, model ="iid") + f(depth, model="rw1") + f(mud, model="rw1") + f(DOY, model="rw2") 

M_smallfish_depth_mud_DOY <- inla(eval(f1),
           family = "nbinomial",
           data = inla.stack.data(Stk_depth_mud_DOY),
           control.compute = list(dic = TRUE,
                                  waic = TRUE,
                                  config = TRUE,
                                  cpo=T),
           control.predictor = list(A = inla.stack.A(Stk_depth_mud_DOY)), verbose = F)

knitr::kable(M_smallfish_depth_mud_DOY$summary.random$fSurvey[, c(1:4, 6)], "simple", digits = 2)
```

smoothers
```{r}
#depth + mud +  DOY

#depth
plot_data <- data.frame(
  depth = M_smallfish_depth_mud_DOY$summary.random$depth$ID,
  median = M_smallfish_depth_mud_DOY$summary.random$depth$`0.5quant`,
  quant0.025 = M_smallfish_depth_mud_DOY$summary.random$depth$`0.025quant`,
  quant0.975 = M_smallfish_depth_mud_DOY$summary.random$depth$`0.975quant`
)

plot_depth <- ggplot(plot_data, aes(x = depth, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=depth), sides = "b", inherit.aes = F) +
  theme_bw()

#mud
plot_data <- data.frame(
  mud = M_smallfish_depth_mud_DOY$summary.random$mud$ID,
  median = M_smallfish_depth_mud_DOY$summary.random$mud$`0.5quant`,
  quant0.025 = M_smallfish_depth_mud_DOY$summary.random$mud$`0.025quant`,
  quant0.975 = M_smallfish_depth_mud_DOY$summary.random$mud$`0.975quant`
)

plot_mud <- ggplot(plot_data, aes(x = mud, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=mud_perc), sides = "b", inherit.aes = F) +
  theme_bw() +
  labs(x = "mud (%)")

#Day of Year
plot_data <- data.frame(
  DOY = M_smallfish_depth_mud_DOY$summary.random$DOY$ID,
  median = M_smallfish_depth_mud_DOY$summary.random$DOY$`0.5quant`,
  quant0.025 = M_smallfish_depth_mud_DOY$summary.random$DOY$`0.025quant`,
  quant0.975 = M_smallfish_depth_mud_DOY$summary.random$DOY$`0.975quant`
)

plot_DOY <- ggplot(plot_data, aes(x = DOY, y = median, ymin = quant0.025, ymax = quant0.975)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(data = trawldata, aes(x=DOY_adjusted), sides = "b", inherit.aes = F) +
  theme_bw() +
  labs(x = "day (since October 1st)")


final_plot <- (plot_depth | plot_mud | plot_DOY) #width 14, height 4
final_plot <- (plot_depth / plot_mud / plot_DOY)

# Print the final figure
print(final_plot)
#ggsave("C:/Users/bparmentier/OneDrive - NIOZ/Documents/NIOZ/code/visbiomassa/figures/smoother_Greenstreet.png", width = 6, height = 10, dpi = 300, units = "in")
```

```{r}
totperstation$Weight_m <- (totperstation$Weight/totperstation$Dist*5)/totperstation$Fraction

pred <- M_smallfish_depth_mud_DOY$summary.fitted.values[c((nrow(totperstation)+1):nrow(data_comb)),] 
pred <- cbind(grid_pred_longlat, pred)

pred$g_m <- exp(pred$'0.5quant')/50*5 #dividing by 50 (=0.2m*1m) *5 to 1m2 

ggplot() +
  geom_tile(data = pred, aes(x = Lon_mid, y = Lat_mid, fill = g_m), show.legend = TRUE, width = 0.05, height = 0.05) + #, width = 0.5, height = 0.
  scale_fill_gradientn(colors = c("beige", "khaki", "green", "yellow", "orange", "red", "darkred", "purple")) +
  theme_bw() +
  geom_sf(data = land_sf, fill = "grey90", color = "black", size = 0.5) +
  labs(x = "Longitude", y = "Latitude", fill = expression(Biomass~(g~m^{-2}))) + 
  geom_sf(data = dog40a_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  geom_sf(data = dog40b_line, color = alpha("black",0.5), linewidth = 0.3, linetype = "dashed", fill=NA) +
  coord_sf(xlim = xlim, ylim = ylim) +
  geom_point(data = totperstation[totperstation$Weight_m > 0, ],
             aes(x = Lon_mid, y = Lat_mid, size = Weight_m),
             color = alpha("black", .2), shape = 1) +
  geom_point(data = totperstation[totperstation$Weight_m ==0 , ],
             aes(x = Lon_mid, y = Lat_mid),
             size = 0.6, color = alpha("black", .6), shape = 16) +
  guides(size = "none") + scale_size_continuous(range = c(0.5, 15)) + theme(text = element_text(size = 7.34))
ggsave(file = file.path(path, "figures/supplement/S10_8-31selection.pdf"), units = "mm", width = 170, height = 255, dpi = 1000)
```

95% interval is determined by taking the mean of the 0.025 and 0.975 quantiles. 
```{r}
#DOG
mean(pred[pred$inside_DOG == T,]$g_m)
mean(exp(pred[pred$inside_DOG == T,]$'0.025quant')/10)
mean(exp(pred[pred$inside_DOG == T,]$'0.975quant')/10)

#NCP
mean(pred[pred$inside_NCP == T,]$g_m)
mean(exp(pred[pred$inside_NCP == T,]$'0.025quant')/10)
mean(exp(pred[pred$inside_NCP == T,]$'0.975quant')/10)
```