identification of metapopulation

based on code developed during Tomas Niarfex intership

exploration/NEA_sensitivity_analysis/metapopulationIdentification.R

+library(tidyverse)
+
+computeAutochtonousRate = function(data){
+  data %>% filter(migrationBasin == originBasin) %>% 
+    select(-originBasin) %>% 
+    inner_join(data %>% group_by(year, migrationBasin) %>% 
+                 summarise(totalRun = sum(effective), .groups = 'drop'), by = c('migrationBasin', 'year')) %>% 
+    mutate(autochtonousRate = effective / totalRun) %>% 
+    select(year, migrationBasin, autochtonousRate)
+}
+
+computeHomingRate = function(data){
+  data %>% filter(migrationBasin == originBasin) %>% 
+    select(-migrationBasin) %>% 
+    inner_join(data %>% group_by(year, originBasin) %>% 
+                 summarise(production = sum(effective), .groups = 'drop'), by = c('originBasin', 'year')) %>% 
+    mutate(homingRate = effective / production) %>% 
+    select(year, originBasin, homingRate)
+}
+
+# comparison between autochtonous rate and homing rate
+# computeAutochtonousRate(exchangesData) %>% 
+#   inner_join(computeHomingRate(data = exchangesData), by =  c('migrationBasin' = 'originBasin', 'year'))
+
+
+
+identyMetapopulation = function(exchangesData, threshold = .95, verbose = FALSE) {
+  exchangesDataUpdated = exchangesData
+  
+  # initialise the connection table between basin and metapopulation
+  metapopulation = exchangesDataUpdated %>% 
+    distinct(year, basin = migrationBasin) %>% 
+    mutate(metapop = basin)
+  
+  iteration = 1
+  
+  # find the basin with the minimum autochtonous rate for the first time,
+  basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+    group_by(year) %>%
+    slice(which.min(autochtonousRate))
+  
+  # loop while autochtonousRate is still < 0.95
+  
+  while (basinWithMinAutochtonousRate$autochtonousRate < threshold ) {
+    #  while (iteration <= 4 ) {
+    if (verbose) cat(iteration, ": ", basinWithMinAutochtonousRate$autochtonousRate, '\n' )
+    
+    # basinWithMinAutochtonousRate will be merged with origin basin sending the maximum number of fish in this basin
+    metapopsToBeMerged <-   exchangesDataUpdated %>% 
+      inner_join(basinWithMinAutochtonousRate %>% select(-autochtonousRate), 
+                 by = c('year', 'migrationBasin' )) %>%
+      filter(migrationBasin != originBasin) %>% # to avoid self merging
+      group_by(year) %>% 
+      slice(which.max(effective)) %>% 
+      ungroup() %>% 
+      select(year, migrationBasin, originBasin)
+    
+    # update the connection table between basin and metapopulation by merging metapopulation
+    for (i in 1:nrow(metapopsToBeMerged)) {
+      mergedName = paste0("M", iteration, '_',i)
+      metapopsToBeMerged_i <- metapopsToBeMerged %>% slice(i) 
+      if (verbose) cat("\t", metapopsToBeMerged_i$migrationBasin, ' + ', metapopsToBeMerged_i$originBasin, ' = ', mergedName, "\n")
+      
+      metapopulation <- metapopulation %>% 
+        mutate(metapop = if_else(year == metapopsToBeMerged_i$year &
+                                   metapop %in%  (metapopsToBeMerged_i %>% select(ends_with("Basin")) %>% unlist(use.names = FALSE)), 
+                                 mergedName,
+                                 metapop)) 
+    }
+    # sum by migration and origin basins according to updated metapopulations
+    # computed with initial exchangesData
+    exchangesDataUpdated <- exchangesData %>% 
+      #sum up on migration basins
+      inner_join(metapopulation, 
+                 by = c('year', 'migrationBasin' = 'basin')) %>% 
+      group_by(year, metapop, originBasin) %>% 
+      summarise(effective = sum(effective), .groups = 'drop') %>% 
+      rename(migrationBasin = metapop) %>% 
+      # sum on origine basins
+      inner_join(metapopulation, 
+                 by = c('year', 'originBasin' = 'basin')) %>% 
+      group_by(year, metapop, migrationBasin) %>% 
+      summarise(effective = sum(effective), .groups = 'drop') %>% 
+      rename(originBasin = metapop)
+    
+    # minimum of autochtonous rate
+    basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+      group_by(year) %>%
+      slice(which.min(autochtonousRate))
+    
+    iteration = iteration + 1
+  }
+  return(list(metapopulation = metapopulation, exchangesData = exchangesDataUpdated))
+}
+
+# upload data =====
+# longer table of effective of different origin basins in each migration basin 
+exchangesData <- read_csv("../../data/output/northeastamerica/effectiveFluxes_1-observed.csv")
+
+exchangesData %>% group_by(year, basin = originBasin) %>% 
+  summarise(production = sum(effective), .groups = 'drop') %>% 
+  inner_join(exchangesData %>% group_by(year, basin = migrationBasin) %>% 
+               summarise(totalRun = sum(effective), .groups = 'drop'), by = c('year', 'basin')) %>% 
+  arrange(totalRun, production) %>% print(n = Inf)
+
+#TODO explain how it is possible to have a production >0  and spawners run = 0: PB in JAVA or 
+# no reproductive success some years
+exchangesData %>% filter(migrationBasin == 'Musquodoboit' ) %>% print(n = Inf)
+exchangesData %>% filter(originBasin == 'Musquodoboit') %>% print(n = Inf)
+
+# remove basin with no run of spawner
+selectedBasins <- exchangesData %>% group_by(year, basin = originBasin) %>% 
+  summarise(production = sum(effective), .groups = 'drop') %>% 
+  inner_join(exchangesData %>% group_by(year, basin = migrationBasin) %>% 
+               summarise(totalRun = sum(effective), .groups = 'drop'), by = c('year', 'basin')) %>% 
+  filter(totalRun > 0) %>%  select(year, basin)
+
+selectedExchangesData <- exchangesData %>% 
+  inner_join(selectedBasins, by = c('year', 'migrationBasin' = 'basin')) %>% 
+  inner_join(selectedBasins, by = c('year', 'originBasin' = 'basin')) 
+
+# results ====
+metapopulations = identyMetapopulation(selectedExchangesData, threshold = .95)
+
+metapopulations$metapopulation %>% group_by(metapop) %>% summarise(n_basin = n(), .groups = 'drop') %>% arrange(desc(n_basin)) %>% 
+  mutate(nb_basin = sum(n_basin)) %>% 
+  print(n = Inf)
+
+computeAutochtonousRate(exchangesData)
+computeAutochtonousRate(metapopulations$exchangesData) %>% print(n = Inf)
+
+# =================================
+# idées abandonnées
+# =================================
+# # compute ratio strayer-from over production
+# # the basin with the minimum of autochonous (maximum of allotochnous) in the spawner run is merged with the origin basin that sends 
+# # the maximum strayers accornding to its production
+#   metapopsToBeMerged  <- exchangesData %>% 
+#   inner_join(basinWithMinAutochtonousRate, by = c('year', 'migrationBasin')) %>% 
+#     #joint with production in each catchment
+#   inner_join(exchangesData %>% group_by(year, originBasin) %>% 
+#                summarise(production = sum(effective)) 
+#              , by = c('year', 'originBasin')) %>% 
+#   filter(migrationBasin != originBasin) %>% 
+#   mutate(ratioStrayerProdution =  effective / production) %>% 
+#   group_by(year) %>% 
+#   slice(which.max(ratioStrayerProdution)) %>% 
+#   ungroup() %>% 
+#   mutate(mergedName = paste0("M", iteration, '_',row_number())) 
\ No newline at end of file