README.Rmd

---
title: "Penang fisheries"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE, fig.width = 10, fig.height = 7, fig.showtext = TRUE, warning = F)
library(ggplot2)
library(magrittr)
library(ggdist)
library(ggforce)
library(brms)
library(showtext)
library(ggrepel)
font_paths("fonts")
font_add(family = "RobotoC", regular = "fonts/Roboto_Condensed/RobotoCondensed-Regular.ttf", 
         bold = "fonts/Roboto_Condensed/RobotoCondensed-Bold.ttf")
showtext_auto()
# theme_set(theme_minimal(base_size = 16, base_family = "Roboto Condensed"))
```

```{r load-data}

drake::loadd(trips)
drake::loadd(points)
drake::loadd(landings_clean)
drake::loadd(landing_model)
drake::loadd(vessel_model)
drake::loadd(species_clean)
drake::loadd(landing_sites)

trips <- trips %>%
  dplyr::filter(date < as.Date("2021-09-01"))
```

```{r trip-stats}
full_info_trips <- dplyr::filter(trips, !is.na(trip_id_landing), !is.na(trip_id_track))

landing_trips <- dplyr::filter(trips, !is.na(trip_id_landing))
tracking_trips <- dplyr::filter(trips, !is.na(trip_id_track))

landing_fisher_tabyl <- landing_trips %>%
  janitor::tabyl(fisher)

main_contributor <- trips %>%
  dplyr::group_by(fisher) %>%
  dplyr::summarise(n_trips_landing_fisher = dplyr::n_distinct(trip_id_landing, na.rm = T), 
                   n_trips = dplyr::n(), 
                   prop = n_trips_landing_fisher/n_trips) %>% 
  dplyr::filter(n_trips_landing_fisher == max(n_trips_landing_fisher)) 
```

We report data collected between `r format(min(trips$date), "%d %B %Y")` and `r format(max(trips$date), "%d %B %Y")`. 
Overall, during these years we recorded information about `r nrow(trips)` fishing trips in Penang. 
The information about these trips was contributed by a total of `r dplyr::n_distinct(trips$fisher)` local fishers and is a combination of tracking and landings. 

Up to `r format(max(trips$date), "%d %B %Y")`, WorldFish installed solar-powered GPS trackers sourced from Pelagic Data Systems Inc. in 24 boats at four key landing sites of Teluk Bahang, Balik Pulau, Kuala Binjai, Kampung Binjai, Kuala Mudah and Sungai Batu.
From those boats, `r dplyr::n_distinct(points$boat_name)` transmitted data successfully.
Most of the trackers were installed in July or August 2020, but three of them were installed in September or November 2019 which allow us to get a better picture of their activity across 2020. 
Until `r format(max(trips$date), "%d %B %Y")` these units tracked a total of `r nrow(tracking_trips)` trips.  

In addition to the tracking data, we also collected landings data. 
Data collection was performed informally through a WhatsApp group set up with participating fishers. 
In total, we collected landings information for `r nrow(landing_trips)` trips by `r dplyr::n_distinct(landing_trips$fisher)` fishers. 
Landings data is available between `r format(min(landing_trips$date), "%d %B %Y")` and `r format(max(landing_trips$date), "%d %B %Y")`.

Three of the fishers that contributed landings data also had a tracker installed on their boats. 
Consequently we have both tracking and landing data for `r nrow(full_info_trips)` trips.
Most of the landing data, however, comes from the single most active fisher in the group.
This fisher alone contributed landing data for `r max(landing_fisher_tabyl$n)` trips which corresponds to about `r scales::percent(main_contributor$prop[1])` of all their fishing trips over the reporting period. 

```{r tracks-figure, fig.height=10}
points %>%
  dplyr::filter(time < lubridate::ymd("2021-09-01")) %>%
  dplyr::left_join(trips, by = c("trip" = "trip_id_track")) %>% 
  dplyr::group_by(trip) %>%
  dplyr::slice_sample(prop = 0.1) %>%
  dplyr::ungroup() %>%
  dplyr::arrange(time) %>%
  dplyr::mutate(trip = forcats::fct_reorder(as.character(trip), !is.na(trip_id_landing))) %>%
  ggplot(aes(x = lng, y = lat, group = trip, colour = !is.na(trip_id_landing))) +
  geom_path(size = 0.3, alpha = 0.3) +
  geom_point(data = landing_sites, colour = "black", group = "") +
  geom_text_repel(data = landing_sites, 
                  aes(group = name, label = name), colour = "black", 
                  box.padding = 2, size = 3.7, segment.size = 0.25, family = "RobotoC") + 
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  scale_color_manual(#values = c("grey70", "grey30"), 
                     values = c("grey70", "#1a4985"),
                     labels = c("trips without landing data", "trips with landing data")) +
  coord_quickmap(xlim = c(99.65, 100.75)) +
  theme_minimal() +
  theme(legend.position = c(1,1),
        legend.justification = c(1,1),
        legend.title = element_blank(), 
        axis.title = element_blank(), 
        axis.text = element_blank(), 
        panel.grid = element_blank(), 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold")) +
  labs(title = paste("A VISUAL SAMPLE OF", nrow(tracking_trips), "FISHING TRIPS IN PENANG"), 
       caption = "* Trips recorded between 01 September 2019 and 31 August 2021", 
       x = "Longitude", y = "Latitude")
```

## Catch weight and income

We are interested in calculating the total catch from artisanal fisheries in Penang and the income it generates. 
We need four key components to perform these calculations. 
First, we need to estimate the expected catch weight from each trip.
Second, we estimate the expected income that the catch provides to artisanal fishers. 
Third, we need to obtain an estimate of the vessel activity throughout the year. 
Fourth, we need an estimate of the number of fishers in Penang. 

We use hierarchical Bayesian models as the framework to estimate these figures for four main reasons. 
First, Bayesian models perform particularly well with relatively small amount of data. 
Second, a hierarchical model allow us to minimise the biases that a arise from unbalances in the contributed data; both across fishers and across time. 
Third, using models, instead of simple averages allow us to obtain estimates of the catch even for periods for which we have no comprehensive sampling. 
Lastly, and most importantly, a Bayesian framework allow us, not only to obtain an estimate of the numbers we are interested on, but also allow us to better understand the uncertainties involved in the calculation. 

```{r model-stats}
models <- list(
  catch = landing_model[[1]], 
  income = landing_model[[2]], 
  vessel = vessel_model)

predict_overall_mean <- function(model){
  model %>%
    posterior_epred(newdata = data.frame(dummy = 1), re_formula = NA) %>%
    as.data.frame() %>%
    tibble::as_tibble() %>%
    dplyr::rename(Estimate = V1)
}

est_sum <- models %>%
  purrr::map(. %>%
               predict_overall_mean() %>%
               posterior_summary(probs = c(0.05, 0.95)) %>%
               purrr::array_tree(margin = 2))
```

First, we look at the expected catch weight from a single trip. 
There was a lot of variability in the catch weight with some trips fishing as little as `r round(min(landing_trips$weight_kg), 1)`kg and some as much as `r round(max(landing_trips$weight_kg), 0)`kg. 
Nevertheless, using ten thousand Monte Carlo simulations in our Bayesian models, we found that the average catch was very likely between `r round(est_sum$catch$Q5, 0)`kg and `r round(est_sum$catch$Q95, 0)`kg (90% credible intervals, mean `r round(est_sum$catch$Estimate, 1)`kg).
The average catch showed important differences among fishers, with some consistently landing more catch by weight than the average. 
More data from a larger number of fishers will allow us to better understand what drives these differences. 

Second, we look at the expected income obtained from a fishing trip. 
We found that each fishing trip was very likely to provide an income between RM`r signif(est_sum$income$Q5, 2)` and RM`r signif(est_sum$income$Q95, 2)` (mean RM`r round(est_sum$income$Estimate)`) to artisanal fishers. 
Similar as the catch weight, average income levels showed a large variation. 
Some trips provided only RM`r round(min(landing_trips$total_price))` as income for the fishers while some trips provided as much as RM`r round(max(landing_trips$total_price))`. 

Third, we investigate the vessel activity coefficient during the study period. 
We found that this coefficient was very likely to be between `r round(est_sum$vessel$Q5, 2)` and `r round(est_sum$vessel$Q95, 2)` (mean `r round(est_sum$vessel$Estimate, 2)`). 
Mathematically, this coefficient can be interpreted in two (equivalent) ways. 
When applied to a single boat, it can be interpreted as the probability that this boat would to a fishing trip in a given day. 
When applied to a group of boats, for example the Penang artisanal fishing fleet, it can be interpreted as the proportion of boats that can be expected to be going on a trip in a given day. 
Although some differences between fishers exist, the largest variations were across time. 
Vessel activity appeared to be lessened during Fridays and Sundays and in some particular weeks of the year, presumably in response to weahter events. 

Lastly, we obtained the information about Penang's fishing fleet based on the number of active and registered fishers. 
In 2018 there were 1,658 fishers registered. 
This number might differ from the number of fishers during the reporting period, but we expect the updated number to be within 50 fishers from the 2018 number. 
We incorporate this uncertainty in all estimates we present.

```{r estimates-period}
estimates_period <- list(
  min = as.Date("2019-09-01"), 
  max = as.Date("2020-08-31")
)
```


Because we only have landings information up until `r format(max(landing_trips$date), "%d %B %Y")`, we present estimates for only the first year of our reporting period between `r format(estimates_period$min, "%d %B %Y")` and `r format(estimates_period$max, "%d %B %Y")`.

```{r yearly-estimates}
predict_yearly_error <- function(model){
  pred_data <- tibble::tibble(date = seq(from = estimates_period$min, 
                                       to = estimates_period$max,
                                       by = 1)) %>%
  dplyr::mutate(month = lubridate::month(date), 
                week = lubridate::isoweek(date), 
                fort = as.numeric(week) %/% 2,
                wday = lubridate::wday(date, label = T))
  
  posterior_epred(
    model, 
    newdata = pred_data, nsamples = 10000, 
    re_formula = ~ (1 | wday) + (1 | week), allow_new_levels = TRUE) %>% 
    set_colnames(as.character(pred_data$date)) %>%
    as.data.frame.table() %>%
    dplyr::mutate(date = as.Date(Var2)) %>%
    dplyr::rename(Estimate = Freq, sample = Var1) %>%
    dplyr::inner_join(pred_data, by = "date") %>%
    tibble::as_tibble() %>%
    dplyr::select(-Var2)
}

posterior_draws <- models %>%
  purrr::map_dfr(predict_yearly_error, .id = "model") %>%
  tidyr::pivot_wider(id_cols = c("date", "sample", "month", "week", "wday", "fort"), 
                     names_from = "model", 
                     values_from = "Estimate") %>%
  dplyr::mutate(n_fishers = round(rnorm(dplyr::n(), 1658, 20)), 
                income = income * vessel * n_fishers,
                catch = catch * vessel * n_fishers, 
                vessel = vessel) %>% 
  dplyr::group_by(month, sample) %>%
  dplyr::summarise(date = min(date), 
                   income = sum(income), 
                   catch = sum(catch), 
                   vessel = mean(vessel))

annual_totals <- posterior_draws %>%
  dplyr::group_by(sample) %>%
  dplyr::summarise(dplyr::across(tidyselect:::where(is.numeric), sum))

ann_sum <- annual_totals %>%
  dplyr::select_if(is.numeric) %>%
  posterior_summary(probs = c(0.05, 0.95)) %>%
  set_names(NULL) %>%
  purrr::array_tree()
```

Using all these pieces of information we estimate that between the period `r format(estimates_period$min, "%d %B %Y")` and `r format(estimates_period$max, "%d %B %Y")` **artisanal fishers were very likely (90% probability) to catch between `r round(ann_sum$catch$Q5/1000000, 1)` and `r round(ann_sum$catch$Q95/1000000, 1)` thousand tonnes in Penang waters, which provided income between RM`r round(ann_sum$income$Q5/1000000, 0)` and RM`r round(ann_sum$income$Q95/1000000, 0)` million to the local communities**. 

```{r annual-figure, fig.height=5}

annual_catch_dist_plot <- annual_totals %>%
  dplyr::slice_sample(n = 1000, replace = F) %>%
  ggplot(aes(x = catch)) +
  stat_dotsinterval(.width = c(.90, .66), fill = "grey90", slab_colour = "grey90") +
  geom_text(aes(x = mean(catch), 
                label = stat(paste0("mean: ", round(mean(x)/1000000, 1), " thousand tonnes")), 
                y = 0, group = 1), stat = "unique", vjust = -1, family = "RobotoC") +
  scale_x_continuous(labels = scales::label_number(scale = 1/1000000, suffix = " thousand\ntonnes", accuracy = 1)) +
  theme_minimal() +
  theme(axis.title = element_blank(), 
        axis.text.y = element_blank(), 
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(), 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold")) +
  labs(title = toupper("Estimated annual catch by artisanal fishers in Penang"), 
       subtitle = paste("Weight"))

annual_income_dist_plot <- annual_totals %>%
  dplyr::slice_sample(n = 1000, replace = F) %>%
  ggplot(aes(x = income)) +
  stat_dotsinterval( .width = c(.90, .66), fill = "grey90", slab_colour = "grey90") +
  geom_text(aes(x = mean(income), 
                label = stat(paste0("mean: RM", round(mean(x)/1000000), " million")), 
                y = 0, group = 1), stat = "unique", vjust = -1, family = "RobotoC") +
  scale_x_continuous(labels = scales::label_dollar(prefix = "RM", scale = 1/1000000, suffix = "\nmillion")) +
  theme_minimal() +
  theme(axis.title = element_blank(), 
        axis.text.y = element_blank(), 
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(), 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold")) +
  labs(subtitle = "Income", 
       caption = "* Each dot corresponds to a simulation in our model. Horizontal lines indicate the 66 and 90% credible intervals")

cowplot::plot_grid(annual_catch_dist_plot, 
                   annual_income_dist_plot, 
                   ncol = 1,
                   align = "v")
```

Although these wide ranges suggest a large uncertainty, it is actually a remarkable achievement considering the limited amount of data used to generate them.
During the year we observe that the catch weight and income stayed within the credible intervals of our estimates throughout the years. 

```{r timeseries-figure}
set.seed(7)

label_date <- function(x){
  format(x, "%B") %>%
    stringr::str_replace("January", paste0("January\n2020"))
}

label_income <- function(x){
  x <- scales::label_dollar(prefix = "RM", scale = 1/1000000)(x)
  x[length(x):(length(x)-1)] <- paste0(x[length(x):(length(x)-1)], "\nmillion")
  x
}

label_weight <- function(x){
  x <- scales::label_number(scale = 1/1000, suffix = "")(x)
  x[length(x):(length(x)-1)] <- paste0(x[length(x):(length(x)-1)], "\ntonnes")
  x
}

plot_ts <- function(x, col){
  x %>% 
  ggplot(aes_string(x = "date", y = col, group = 1)) +
  stat_lineribbon(size = 0.5,
                  alpha = 1/2,
                  .width = c(.90, .66)) +
  scale_fill_brewer(palette = "Greys") +
  expand_limits(y = 0) +
  scale_x_date(date_breaks = "month", minor_breaks = NULL, labels = label_date) +
  theme_minimal() +
    theme(axis.title = element_blank(), legend.position = "none", 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold"))
}

catch_ts_plot <-
  posterior_draws %>%
  plot_ts("catch") +
    scale_y_continuous(labels = label_weight) +
    labs(title = toupper("Monthly catch totals by artisanal fishers in Penang"), 
       subtitle = "Weight")
  
income_ts_plot  <- 
  posterior_draws %>%
  plot_ts("income") +
    scale_y_continuous(labels = label_income) +
  labs(subtitle = "Income",
       caption = "* Shaded areas indicate the 66 and 90% credible intervals")

vessel_ts_plot <- 
  posterior_draws %>%
  plot_ts("vessel")  +
  scale_y_continuous(labels = scales::label_percent()) +
  labs(subtitle = "Vessel activity")

cowplot::plot_grid(#title_plot, 
                   catch_ts_plot, 
                   income_ts_plot, 
                   rel_heights = c(1.05, 1),
                   # vessel_ts_plot,
                   ncol = 1, align = "v")
```

A larger sample size, both for landings and tracking (which we use to determine vessel activity) will allow us to detect smaller variations across time and answer other questions to improve the management of the fisheries and improve the livelihoods of artisanal fishers. 

## Catch composition

We also looked at the catch composition in the trips where we recorded data.
Although composition data is, so far, biased and should not be taken as indicative of Penang's catch, we can already distinguish some interesting patterns. 
*Text to be written once species/grade data is of sufficient quality*.

```{r catch-composition-figure, fig.height = 3.5}

get_root_name <- function(x){
  x %>%
    stringr::str_extract("[a-z\b]+") %>%
    stringr::str_sub(start = 1, end = -2)
}

species_summary <- landings_clean %>%
  dplyr::filter(!is.na(price_kg)) %>%
  dplyr::left_join(species_clean) %>%
  dplyr::mutate(species = malay_common_name, 
                price = weight_kg * price_kg) %>%
  dplyr::group_by(species, common_malay) %>%
  dplyr::summarise(
    frequency = dplyr::n_distinct(trip_id), 
    weight = sum(weight_kg, na.rm = TRUE), 
    price_kg = mean(price_kg, na.rm = TRUE), 
    price = sum(price)) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(species_p = forcats::fct_lump_prop(species, 
                                              w = price,
                                              prop = 0.02, 
                                              other_level = NA),
                species_w = forcats::fct_lump_prop(species, 
                                                   w = weight,
                                                   prop = 0.02, 
                                                   other_level = NA)) %>%
  dplyr::mutate(dplyr::across(dplyr::starts_with("species"), as.character),
                species = dplyr::coalesce(species_w, species_p), 
                species = tidyr::replace_na(species, "Other"), 
                dplyr::across(tidyselect:::where(is.numeric), ~ . / sum(.))) 

p1 <- species_summary %>%
  dplyr::mutate(dplyr::across(c("species", "common_malay"), 
                              ~ forcats::fct_reorder(., price, sum)), 
                species_p = forcats::fct_reorder(species_p, price, sum, .desc = T),
                focal = as.numeric(species_p) <= 3,
                focal = tidyr::replace_na(focal, FALSE),
                species = forcats::fct_relevel(species, "Other")) %>%
  ggplot(aes(x = species, y = price)) +
  geom_col(aes(fill = focal), alpha = 0.7) +
  # geom_text(aes(label = scales::percent(price, accuracy = 1))) +
  coord_flip() +
  theme_minimal() +
  scale_y_continuous(expand = expansion(mult = c(0,0)), 
                     labels = scales::label_percent()) +
  scale_fill_manual(values = c("grey70", "#1a4985"))+
  theme(legend.position = "none", axis.title = element_blank(), 
        # plot.margin = unit(c(1,4,1,0), units = "mm"), 
        plot.subtitle = element_text(hjust = 0.5), 
        axis.text.y = element_text(hjust = 0.5), 
        panel.grid.major.y = element_blank(), 
        plot.background = element_blank(), 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold")) +
  labs(subtitle = "Income", 
       caption = "* Species that accounted for less than 1% of the total catch or weight were grouped together")

p2  <- species_summary %>%
  dplyr::mutate(dplyr::across(c("species", "common_malay"), 
                              ~ forcats::fct_reorder(., weight, sum)),
                species_p = forcats::fct_reorder(species_p, price, sum, .desc = T),
                focal = as.numeric(species_p) <= 3,
                focal = tidyr::replace_na(focal, FALSE),
                  species = forcats::fct_relevel(species, "Other")) %>%
  ggplot(aes(x = species, y = weight)) +
  geom_col(aes(group = common_malay, fill = focal), alpha = 0.7) +
  coord_flip()+
  scale_y_reverse(expand = expansion(mult = 0), 
                     labels = scales::label_percent()) +
  scale_x_discrete(position = "top") +
  scale_fill_manual(values = c("grey70", "#1a4985"))+
  theme_minimal() +
  theme(legend.position = "none", axis.title = element_blank(), 
        # plot.margin = unit(c(1,0,1, 4), units = "mm"), 
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.subtitle = element_text(hjust = 0.5), 
        panel.grid.major.y = element_blank(), 
        plot.background = element_blank(), 
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(face = "bold")) +
  labs(subtitle = "Weight")

plink <- species_summary %>%
  dplyr::mutate(species_p = forcats::fct_reorder(species, price, sum),
                species_w = forcats::fct_reorder(species, weight, sum), 
                species_p2 = forcats::fct_reorder(species_p, price, sum, 
                                                  .desc = T),
                focal = as.numeric(species_p2) <= 3,
                focal = tidyr::replace_na(focal, FALSE),
                dplyr::across(dplyr::starts_with("species"), 
                              ~ forcats::fct_relevel(., "Other")), 
                dplyr::across(dplyr::starts_with("species"), 
                              as.numeric)) %>% 
  dplyr::distinct(species_p, species_w, focal) %>%
  ggplot(aes(colour = focal)) +
  geom_diagonal(aes(x = 1, xend = 2, y = species_w, yend = species_p), strength = 0.5, alpha = 0.7) +
  geom_point(aes(x = 1, y = species_w), alpha = 0.7) +
  geom_point(aes(x = 2, y = species_p), alpha = 0.7) +
  scale_color_manual(values = c("grey70", "#1a4985"))+
  theme_no_axes(theme_minimal())+
  theme(legend.position = "none", axis.title = element_blank(), 
        plot.margin = unit(c(1,0,1, 0), units = "mm"),
        text = element_text(family = "RobotoC"), 
        plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = toupper("Contribution of different species to reported landings"))
  

cowplot::plot_grid(p2, plink, p1, 
                   nrow = 1, align ="h", rel_widths = c(1,0.2, 1))
```

## Other unfinished plots

(not to be included in final document)

```{r}
predict_fisher_differences <- function(model){
  model %>%
    posterior_epred(newdata = data.frame(fisher = unique(model$data$fisher)), re_formula = ~ (1 | fisher)) %>%
    set_colnames(unique(model$data$fisher)) %>%
    as.data.frame.table() %>%
    dplyr::rename(fisher = Var2) %>%
    dplyr::rename(Estimate = Freq, sample = Var1)
}

models %>%
  purrr::map_dfr(.f = predict_fisher_differences, .id = "model") %>%
  dplyr::mutate(
    # Randomise names to protect identity
    fisher = factor(
      x = fisher, 
      labels =  randomNames::randomNames(
        n = dplyr::n_distinct(fisher), 
        ethnicity = 6, 
        which.names = "first", 
        sample.with.replacement = FALSE)),
    # Order fishers by catch weight
    fisher = forcats::fct_reorder2(
      .f = fisher, 
      .x = Estimate, 
      .y = as.numeric(model == "catch"), 
      .fun = weighted.mean, 
      .desc = TRUE), 
    # Reverse order so that fishers without catch show at the bottom of the plot
    fisher = forcats::fct_rev(f = fisher),
    model = forcats::fct_relevel(model, "catch")) %>%
  ggplot(aes(x = Estimate, y = fisher)) +
  stat_pointinterval(point_size = 1.5) +
  facet_wrap("model", scales = "free_x") +
  theme_minimal() +
  labs(
    title = "There were important differences among fishers",
    caption = "* Fisher names have been changed to protect their identity")

predict_wday_differences <- function(model){
  model %>%
    posterior_epred(newdata = data.frame(wday = unique(model$data$wday)), re_formula = ~ (1 | wday)) %>%
    set_colnames(unique(model$data$wday)) %>%
    as.data.frame.table() %>%
    dplyr::rename(wday = Var2) %>%
    dplyr::rename(Estimate = Freq, sample = Var1)
}

models %>%
  purrr::map_dfr(.f = predict_wday_differences, .id = "model")%>%
  ggplot(aes(x = Estimate, y = wday)) +
  stat_pointinterval() +
  facet_wrap("model", scales = "free_x")

predict_week_differences <- function(model){
  model %>%
    posterior_epred(newdata = data.frame(week = unique(model$data$week)), re_formula = ~ (1 | week)) %>%
    set_colnames(unique(model$data$week)) %>%
    as.data.frame.table() %>%
    dplyr::rename(week = Var2) %>%
    dplyr::rename(Estimate = Freq, sample = Var1)
}

models %>%
  purrr::map_dfr(.f = predict_week_differences, .id = "model") %>%
  dplyr::mutate(week = forcats::fct_reorder(week, as.numeric(as.character(week)))) %>%
  ggplot(aes(x = Estimate, y = week)) +
  stat_pointinterval() +
  facet_wrap("model", scales = "free_x")
```


```{r, message=FALSE}
trip_summary <- landings_clean %>%
  dplyr::group_by(trip_id) %>%
  dplyr::summarise(
    date = dplyr::first(date), 
    fisher = dplyr::first(fisher), 
    total_price = sum(total_price), 
    weight_kg = sum(weight_kg),
    n_species = dplyr::n_distinct(common_malay))

fisher_summary <- landings_clean %>%
  dplyr::group_by(fisher) %>%
  dplyr::summarise(
    imei = dplyr::first(imei),
    first_trip = min(date), 
    last_trip = max(date),
    n_trips = dplyr::n_distinct(trip_id),
    total_price = sum(total_price), 
    weight_kg = sum(weight_kg),
    n_species = dplyr::n_distinct(common_malay)) %>%
  dplyr::arrange(-n_trips)
```


```{r}
trip_summary %>%
  ggplot(aes(x = date)) +
  geom_histogram(binwidth = 14, boundary = as.Date("2020-01-01")) +
  theme_minimal()
```

```{r}
fisher_summary %>%
  dplyr::mutate(fisher = forcats::fct_reorder(fisher, first_trip, .desc = TRUE)) %>%
  ggplot(aes(y = fisher)) +
  geom_segment(aes(x = first_trip, 
                   xend = last_trip, yend = fisher, colour = is.na(imei))) +
  geom_point(data = trip_summary, 
             aes(x = date), size = 1, shape = 3) +
  theme_minimal()
```
Between `r format(min(trip_summary$date), "%d %B %Y")` and `r format(max(trip_summary$date), "%d %B %Y")` we recorded landings for `r nrow(trip_summary)` trips by `r dplyr::n_distinct(landings_clean$fisher)` fishers. 



During this time, GPS trackers were installed in `r length(na.omit(fisher_summary$imei))` boats. 

```{r}
trip_summary %>%
  ggplot(aes(x = date, y = weight_kg)) +
  geom_point(aes(colour = fisher))

trip_summary %>%
  ggplot(aes(x = date, y = total_price)) +
  geom_point(aes(colour = fisher))
```