Celem poniższego projektu jest analiza zbioru danych dotyczących pomiarów śledzi atlantyckich i warunków w jakich żyją, zbieranych w ciągu ostatnich 60 lat, i próba znalezienia przyczyn stopniowego spadku ich długości w ostatnich latach. W procesie tej analizy zbadano ewentualne zależności pomiędzy atrybutami zbioru i utworzono klasyfikator dokonujący predykcji rozmiaru śledzi na podstawie pozostałych pomiarów. W wyniku poniższych kroków sformułowano trzy główne czynniki wpływające na spadek długości śledzi w ostatnich latach.

Poniższy kod ładuje wymagane biblioteki i zapewnia powtarzalność wyników:

library(ggplot2)
library(plotly)
library(VIM)
library(knitr)
library(ggcorrplot)
library(caret)       
library(randomForest)
set.seed(42)

Zbiór danych składa się z pomiarów śledzi i warunków w jakich żyją, zebranych przez ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Kolejne kolumny w zbiorze danych to:

• length --- długość złowionego śledzia [cm];
• cfin1 --- dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1];
• cfin2 --- dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2];
• chel1 --- dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1];
• chel2 --- dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2];
• lcop1 --- dostępność planktonu [zagęszczenie widłonogów gat. 1];
• lcop2 --- dostępność planktonu [zagęszczenie widłonogów gat. 2];
• fbar --- natężenie połowów w regionie [ułamek pozostawionego narybku];
• recr --- roczny narybek [liczba śledzi];
• cumf --- łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku];
• totaln --- łączna liczba ryb złowionych w ramach połowu [liczba śledzi];
• sst --- temperatura przy powierzchni wody [°C];
• sal --- poziom zasolenia wody [Knudsen ppt];
• xmonth --- miesiąc połowu [numer miesiąca];
• nao --- oscylacja północnoatlantycka [mb].

Wczytywanie oryginalnego zbioru danych z wyspecyfikowaniem typu danych w kolumnach:

sample <- read.csv("sledzie.csv", nrows = 100)
df_classes <- sapply(sample, class) 
df_classes <- replace(df_classes, df_classes=="factor", "numeric")
all_df <- read.csv("sledzie.csv", colClasses = df_classes, fill=TRUE, na.string=c("NA", "?"))

Próbka danych ze zbioru:

Oryginalny zbiór składa się z 52582 pomiarów, z czego 10094 (około 19%) zawiera brakujące wartości:

##   n_samples n_missing_values percent_missing_values
## 1     52582            10094              0.1919668

Poniżej zostały zaprezentowane szczegółowe informacje na temat brakujących wartości pomiarów:

## 
##  Variables sorted by number of missings: 
##  Variable      Count
##     lcop1 0.03143661
##     lcop2 0.03025750
##       sst 0.03012438
##     cfin1 0.03006732
##     chel2 0.02959188
##     chel1 0.02957286
##     cfin2 0.02921152
##         X 0.00000000
##    length 0.00000000
##      fbar 0.00000000
##      recr 0.00000000
##      cumf 0.00000000
##    totaln 0.00000000
##       sal 0.00000000
##    xmonth 0.00000000
##       nao 0.00000000

Można zauważyć, że braki w pomiarach dotyczą 7 atrybutów - dostępności poszczególnych gatunków planktonu oraz temperatury przy powierzchni wody.

Możliwym rozwiązaniem tego problemu byłoby użycie metody k-nn, aby uzupełniać brakujące dane w pomiarze na podstawie tych znajdujących się w najpodobniejszych do niego pomiarach.

Po dokładniejszym przyjrzeniu się danym, okazuje się, że w przypadku większości takich wierszy można łatwo uzupełnić te braki, ponieważ otaczające je pomiary, pochodzące z tego samego połowu, są bardzo podobne lub identyczne, mogą więc być wykorzystane do wypełnienia braków, bez uszczerbku na ogólnej jakości zbioru.

Ponieważ dane są uporządkowane chronologicznie można założyć, że jeśli dwa sąsiednie pomiary mają taką samą wartość atrybutu xmonth to mogą posłużyć do uzupełnienia braków sąsiada. Aby uprościć proces, w przypadku, jeśli sąsiednie pomiary pochodzą już z innego miesiąca lub mają takie same braki jak pomiar aktualnie uzupełniany to wartość brakująca jest nadpisywana przez średnią dla danego atrybutu w całym zbiorze.

Brakujące dane uzupełnia następujący kod:

missing_vals <- apply(all_df, 1, anyNA)

if(missing_vals[1]){
  for (x in 1:ncol(all_df)) {
    if (is.na(all_df[1, x])) {
      all_df[1, x] <- all_df[2, x]
    }
  }
}

for(i in 2:nrow(all_df)-1){
  if(missing_vals[i]){
    for (x in 1:ncol(all_df)) {
      if (is.na(all_df[i, x])) {
        if(!is.na(all_df[i - 1, x]) && (all_df$xmonth[i] == all_df$xmonth[i - 1])){
          all_df[i, x] <- all_df[i - 1, x]
        }else if(!is.na(all_df[i + 1, x]) && (all_df$xmonth[i] == all_df$xmonth[i + 1])){
          all_df[i, x] <- all_df[i + 1, x]
        }else{
          all_df[i, x] <- mean(all_df[,x], na.rm = TRUE)
        }
      }
    }
  }
}

if(missing_vals[nrow(all_df)]){
  for (x in 1:ncol(all_df)) {
    if (is.na(all_df[nrow(all_df), x])) {
      all_df[nrow(all_df), x] <- all_df[nrow(all_df)-1, x]
    }
  }
}

Niestety atrybut xmonth nie wystarczy by jednoznacznie oddzielić kolejne lata pomiarów, ponieważ jego wartość zmienia się częściej niż by to wynikało z rzeczywistego okresu zbierania pomiarów:

months <- 0
curr_month <-0
for(i in 1:nrow(all_df)){
  if(curr_month != all_df$xmonth[i]){
    months <- months + 1
    curr_month <- all_df$xmonth[i]
  }
}
months/12

## [1] 317.5833

W związku z tym zostały utworzone sztuczne atrybuty year oraz years_group, dzielące zbiór odpowiednio na 60 i 10 równej wielkości grup. Jest to podział dokonany na potrzeby analizy i wizualizacji danych.

Przebieg poniższego wykresu obrazuje stopniowy spadek rozmiaru śledzia z czasem:

Na wykresach można zaobserwować zmniejszenie się dostępności wszystkich gatunków planktonu w ostatnich latach, co pokrywa się z mniejszymi rozmiarami osiąganymi przez śledzie. Zależność wydaje się być szczególnie widoczna w przypadku chel1 i lcop1. Potwierdza to intuicyjny wniosek, że dostępność pokarmu może mieć znaczący wpływ na rozmiary ryb.

Powyższe wykresy nie obrazują wyraźnych zależności. Szczególnie mała instensywność połowów przypada zarówno na okres osiągania przez śledzie największych jak i najmniejszych odnotowanych rozmiarów. Dodatkowo ze względu na problem z wyraźnym podziałem pomiarów na lata w zbiorze danych atrybuty "roczne" stają się mniej wiarygodne.

Chociaż poziom zasolenia nie wydaje się mieć wielkiego wpływu na rozmiary śledzi, to wykresy sugerują znaczący wpływ temperatury wody i oscylacji północnoatlantyckiej w tej kwestii.

Badanie wykazało istnienie korelacji między długością śledzi a dostępnością niektórych gatunków planktonu (chel1 i lcop1) oraz ułamkiem pozostawionego narybku (fbar), a także odwrotnej korelacji z temperaturą przy powierzchni wody (sst) oraz oscylacją północnoatlantycką (nao). Można więc wnioskować, że na zmniejszenie się rozmiarów osiąganych przez śledzie wpłynęły zmniejszona dostępność pokarmu, niekorzystne warunki środowiska oraz zwiększone natężenie połowów. Dodatkowo można także zauważyć zależności między oscylacją a wspomnianymi miarami dostępności planktonu oraz temperaturą wody.

W celu predykcji długości wyłowionej ryby została wykorzystana metoda RandomForest, a przyjętymi miarami oceny jakości klasyfikacji $R^2$ oraz $RMSE$. Poniżej zaprezentowany został kod tworzący i uczący model oraz wyniki jego działania.

rf_df <- all_df[,c(2:14,16)]

inTraining <- createDataPartition(y = rf_df$length, p = 0.8, list = FALSE)

training <- rf_df[inTraining, ]
testing <- rf_df[-inTraining, ]

Mtries <- tuneRF(training, training$length, ntree = 30, plot = FALSE)

## mtry = 4  OOB error = 0.07235467 
## Searching left ...
## mtry = 2 	OOB error = 0.2741899 
## -2.789527 0.05 
## Searching right ...
## mtry = 8 	OOB error = 0.003078526 
## 0.9574523 0.05 
## mtry = 14 	OOB error = 0.0002513633 
## 0.9183495 0.05

rfGrid <- expand.grid(mtry = Mtries)
gridCtrl <- trainControl(method = "repeatedcv", number = 2, repeats = 5)

fit <- train(length ~ ., data = training, method = "rf", metric = "RMSE", preProc = c("center", "scale"), trControl = gridCtrl, tuneGrid = rfGrid, ntree = 30, importance=TRUE)

rfClasses <- predict(fit, newdata = testing)
postResample(pred = rfClasses, obs = testing$length)

##      RMSE  Rsquared       MAE 
## 1.1937860 0.4712525 0.9453771

Wykres poniżej wskazuje jako najważniejsze atrybuty związane z dostępnością planktonu, temperaturę blisko powierzchni wody oraz łączną liczbę ryb wyłowionych w ramach połowu:

Na podstawie przeprowadzonej analizy można wyróżnić trzy czynniki, mające znaczny wpływ na spadającą długość śledzi w ostatnich latach:

• dostępność pokarmu --- plankton stanowi główne źródło pokarmu śledzia atlantyckiego, więc jego dostępność stanowi dobrą podstawę do prób predykcji rozmiarów osiąganych przez ryby
• warunki środowiskowe --- oscylacje północnoatlantyckie mogą wpływać na temperaturę wody, ta z kolei razem z zasoleniem, mogą stwarzać niekorzystne warunki dla śledzi, mogą też bezpośrednio wpływać na zmniejszenie dostępności planktonu
• wpływ człowieka --- analiza wykazała, że intensywność połowów także ma wpływ na długość śledzi, która osiągała średnio większe wartości w okresie, gdy pozostawiano więcej narybku.

Warto jednak podkreślić, że śledź atlantycki jest częścią większego ekosystemu i na jego rozwój wpływać będą liczne czynniki, które nie zostały uwzględnione w zbiorze danych. Być może na podstawie szerszego zbioru atrybutów udałoby się znaleźć więcej takich zależności.