8  Прогноз пополнения: от факторов до ансамбля

8.1 Введение

Начнём с ловушки, в которую легко попасть всем нам. Когда перед глазами аккуратные таблицы предикторов и графики «запас–пополнение», возникает приятное ощущение, что, выбрав правильную модель и пару красивых параметров, мы можем управлять будущим — будто поставив регулятор на приборной панели, повысим пополнение в нужный год. Это и есть иллюзия контроля: данные неполны, среда нестационарна, а часть закономерностей — просто шум, который мозг охотно принимает за сигнал. Наша задача — осознанно тормозить, проверять допущения, избегать утечки информации из будущего в прошлое и признавать, что длинные хвосты — не исключение, а рабочая среда экологии. При этом, повторюсь, скатываться в цинизм не нужно: рациональный, аккуратно спланированный анализ действительно улучшает предсказания и решения — в смысле медленного, но реального прогресса, если мы дисциплинируем энтузиазм диагностикой и воспроизводимостью. Биология редко ведёт себя как учебник: нелинейности, пороги, запаздывания и переменная «уловистость» индикаторов — всё это часть истории, и именно поэтому важно строить не одну «идеальную» модель, а цикл проверки конкурирующих объяснений, опираясь и на предметную логику, и на статистическую строгость.

Эта практическая работа — про такой цикл. Мы изучаем зависимость пополнения R3haddock (возможно это пикша) от среды и нерестового запаса, идём от подготовки данных до прогноза с доверительными интервалами, намеренно сопоставляя разные семейства моделей. Сначала приводим предикторы к численному виду, осмысленно обходимся с пропусками и сокращаем мультиколлинеарность, чтобы не строить интерпретацию на «двух зеркальных термометрах». Затем запускаем автоматический отбор признаков двумя методами — Boruta (в духе нелинейной, «лесной» логики) и LASSO (строгая сжатая линейная постановка). Их пересечение, скорректированное биологическим смыслом (например, сохранением нерестового запаса haddock68), даёт устойчивый стартовый набор факторов: температуры и океанография (T…, O…), биотические индикаторы (например, codTSB) и нерестовый запас. Далее мы намеренно не женимся на одном «красивом» уравнении: ставим рядом механистические модели «запас–пополнение» (Рикер и Бивертон–Холт) и статистические LM/GLM/GAM, чтобы увидеть, где данные действительно поддерживают горб плотностной зависимости, а где кривая честно выходит на плато; где линейная аппроксимация достаточна, а где гладкие функции выявляют оптимумы, пороги и нелинейные эффекты. Такой параллельный взгляд — это не прихоть, а способ не перепутать удобную историю с реальной динамикой.

Критический момент — валидация во времени. Обычная случайная кросс‑валидация льстит нам, подмешивая будущее в прошлое; мы избегаем этого через time‑slice с расширяющимся окном и горизонтом прогноза, а затем фиксируем результат внешним хронологическим тестом. Так мы проверяем не только «как хорошо объяснили вчера», но и «как не обманулись насчёт завтра». По итогам сравнения берём не «победителя по AIC/Р² любой ценой», а устойчивую к хронологии схему — иногда это GAM, иногда более простая GLM, а порой и ансамбль, который признаёт, что комбинирование разных источников знания (простого и гибкого) часто надёжнее любого одиночного героя. Прогноз на 2022–2024 мы даём не как цифры отлитые в граните, а как веер с 50% и 95% интервалами — потому что у природы длинные хвосты, и задача аналитика — показать диапазон правдоподобного, а не притворяться владельцем хрустального шара. И наконец — про нарратив: любой вывод модели — это маленькая история о будущем запаса. Она приемлема для управления только тогда, когда прошла проверку данными, диагностикой и попыткой опровержения альтернативами. В этом практикуме мы именно так и поступаем: даём данным говорить, а себе — сомневаться, сравнивать и проверять. Именно так рождаются прогнозы, которые выдерживают столкновение с реальностью.

И так в сухом остатке, в этой практической работе представлен цикл прикладного анализа зависимости пополнения запаса гидробионта от факторов среды (в том числе нерестового запаса): от подготовки данных и отбора предикторов до сравнения нескольких семейств моделей, выбора устойчивой к хронологии прогностической схемы и построения прогноза с доверительными интервалами. Подход ориентирован на начинающих, но использует современные приёмы: автоматический отбор признаков (Boruta, LASSO), сопоставление линейных/нелинейных моделей, time-slice валидацию и ансамблевый прогноз. Целевая переменная: R3haddock — пополнение запаса. Кандидатные предикторы: гидрометеорология (температуры T…), океанография (O…), биотические показатели (например, codTSB) и нерестовый запас (haddock68). Цель анализа: понять, какие факторы и в каких формах оказываются значимыми, отобрать рабочий набор моделей и получить прогноз на 2022–2024 с оценкой неопределенности.

Настоящий анализ разделен на несколько этапов:

1. Выбор предикторов. Скрипт можно скачать по ссылке.

2. Построение биологически мотивированных (механистических) нелинейных классических моделей «запас-пополнение» Рикера и Бивертона-Холта. Анализ их значимости и сравнение с моделями LM/GLM/GAM. Скрипт можно скачать по ссылке.

3. Построение классических статистических моделей LM/GLM/GAM. Анализ их прогностических способностей и выполнение прогноза. Скрипт можно скачать по ссылке.

4. Полный цикл ( СКРИПТ) прикладного анализа зависимости пополнения запаса гидробионта от факторов среды, включающий классические модели, машинное обучение, а также этапы:

• а) выбор предикторов;

• б) базовое сравнение различных моделей;

• в) выбор лучшей прогностической модели;

• г) ансамблевый прогноз.

Входные данные для работы скрипта: RECRUITMENT.xlsx, а также промежуточный файл с готовым набором предикторов: selected_predictors_dataset.csv.

8.2 Выбор предикторов

В процессе анализа факторов, влияющих на пополнение запасов гидробионтов, важным этапом является тщательная подготовка данных и отбор наиболее информативных предикторов, поскольку качество последующих моделей напрямую зависит от качества входных данных. Начиная с первичной обработки, мы приводим все потенциальные предикторы к числовому формату, так как большинство статистических и машинно-обучаемых моделей требуют именно такого представления данных, при этом заменяем строковые обозначения пропущенных значений «NA» на стандартные NA, что позволяет системе R корректно обрабатывать отсутствующие наблюдения. Для заполнения пропусков мы применяем медианную импутацию, которая представляет собой простой и устойчивый к выбросам метод, поскольку медиана менее чувствительна к экстремальным значениям по сравнению со средним. Хотя существуют и более сложные альтернативы, такие как множественная импутация с использованием пакета mice, KNN-импутация через recipes::step_impute_knn или даже методы, специально разработанные для временных рядов, например, фильтр Калмана или ARIMA-модели, медианная импутация остается практичным выбором для начального этапа анализа, особенно когда объем данных ограничен или временные зависимости не являются доминирующими. Следующим важным этапом является анализ корреляционной структуры данных, поскольку высокая мультиколлинеарность между предикторами может серьезно ухудшить интерпретацию моделей и завысить дисперсию оценок параметров, особенно в линейных моделях. Для автоматического выявления и устранения сильно коррелированных переменных мы используем функцию findCorrelation с пороговым значением коэффициента корреляции 0.8, что позволяет сохранить лишь один представитель из каждой группы высококоррелированных переменных. Хотя альтернативными подходами могут служить диагностика по значениям VIF или применение методов снижения размерности, таких как PLS или PCA, удаление явно коррелированных предикторов оказывается наиболее прямолинейным решением для обеспечения стабильности последующих моделей. Для автоматического отбора наиболее значимых предикторов мы применяем два дополнительных метода, которые по-разному подходят к этой задаче и тем самым обеспечивают взаимную проверку результатов. Boruta представляет собой обертку алгоритма Random Forest, которая генерирует «теневые» переменные, полученные путем случайного перемешивания исходных признаков, и сравнивает важность реальных предикторов с этими теневыми копиями, сохраняя только те переменные, чья важность статистически превосходит уровень шума. Этот метод особенно эффективен при наличии нелинейных зависимостей и взаимодействий между переменными, демонстрируя высокую устойчивость к шуму, хотя и требует больше вычислительных ресурсов и может излишне благоволить к группам коррелированных признаков. Параллельно мы применяем LASSO-регрессию из пакета glmnet, которая использует L1-регуляризацию для зануления коэффициентов слабо влияющих предикторов, тем самым выполняет отбор признаков в процессе оценки модели. При выборе оптимального значения параметра регуляризации lambda мы сознательно предпочитаем значение lambda.1se, которое соответствует более простой модели, но при этом находится в пределах одной стандартной ошибки от минимального значения ошибки, так как этот консервативный подход часто обеспечивает лучшую обобщающую способность на небольших выборках, характерных для экологических данных. Однако LASSO имеет свои ограничения: он чувствителен к масштабу переменных, что делает центрирование и стандартизацию обязательными предварительными шагами, и предполагает линейную форму зависимости между предикторами и откликом, что может не соответствовать реальной биологической природе процессов. Финальный набор предикторов формируется как объединение результатов Boruta и LASSO с учетом биологической логики, что повышает устойчивость отбора к случайным флуктуациям, присущим каждому отдельному методу, и гарантирует включение ключевых переменных, таких как нерестовый запас (haddock68), который биологически должен влиять на пополнение запаса. Для предварительной проверки значимости отобранных предикторов мы строим простую линейную модель, которая не предназначена для окончательного прогноза, но служит в качестве sanity-check, позволяя оценить порядок величины эффектов и выявить явно незначимые или противоречащие биологической логике переменные. Важно отметить несколько нюансов и потенциальных подводных камней, с которыми можно столкнуться на этом этапе: если распределение целевой переменной R3haddock сильно скошено, может потребоваться лог-трансформация или использование моделей, специально разработанных для положительных откликов, таких как Gamma GLM; корреляция между переменными не обязательно отражает причинно-следственные связи, и при удалении высококоррелированных предикторов мы можем потерять полезную информацию, поэтому в некоторых случаях лучше применять методы, сохраняющие информацию из всех переменных, например, PLS или GAM; наконец, медианная импутация, хотя и проста в применении, может быть недостаточно точной для временных рядов, где хронологически осмысленная импутация, такая как скользящая медиана или интерполяция, часто дает более реалистичные результаты, учитывающие естественную динамику экологических процессов. Таким образом, этап подготовки данных и отбора предикторов представляет собой критически важный фундамент для последующего построения качественных моделей прогнозирования пополнения рыбных запасов, где баланс между статистической строгостью и биологической интерпретируемостью определяет успех всего анализа.

Скрипт целиком можно скачать по ссылке

# ==============================================================================
# 1) ВЫБОР ПРЕДИКТОРОВ
# ------------------------------------------------------------------------------
# Цель блока: привести данные к числовому виду, обработать пропуски, сократить
# мультиколлинеарность (сильные корреляции), а затем автоматически выделить
# кандидатов-предикторов двумя методами (Boruta, LASSO). В конце сформируем
# финальный пул признаков и проверим их значимость в простой LM.
# ==============================================================================

# Установка и подключение необходимых библиотек
# Для автоматического отбора предикторов нам понадобятся дополнительные пакеты
if (!require("pacman")) install.packages("pacman")

Загрузка требуемого пакета: pacman

pacman::p_load(
  readxl, tidyverse, caret, corrplot, mgcv, randomForest, xgboost,
  Boruta,GGally, FactoMineR, glmnet, recipes, rsample  # Новые библиотеки для автоматического отбора
)

# Очистка среды и установка рабочей директории
# Совет: rm(list=ls()) очищает все объекты в памяти R; setwd задаёт папку,
# где искать/сохранять файлы. Убедитесь, что путь корректен на вашей машине.
rm(list = ls())
setwd("C:/RECRUITMENT/")

# Пакеты для расширенного отбора предикторов
# Boruta — обёртка над Random Forest для отбора признаков;
# glmnet — регуляризация (LASSO/ElasticNet) для отбора/усиления обобщающей способности;
# FactoMineR — PCA и другие многомерные методы (используем как утилиту).
library(Boruta)   # Алгоритм обертки для отбора признаков
library(glmnet)   # LASSO-регрессия
library(FactoMineR) # PCA анализ


# Загрузка и первичная обработка данных
# Шаги: фильтруем годы, приводим типы к числовому, заменяем строковые "NA" на NA.
DATA <- readxl::read_excel("RECRUITMENT.xlsx", sheet = "RECRUITMENT") %>%
  filter(YEAR > 1989 & YEAR < 2022) %>%
  # Преобразуем необходимые столбцы в числовой формат
  mutate(
    across(starts_with("T"), as.numeric),
    across(starts_with("I"), as.numeric),
    across(starts_with("O"), as.numeric),
  ) %>%
  # Обработка пропущенных значений (заменяем строку "NA" на NA)
  mutate(across(where(is.character), ~na_if(., "NA")))

# 1. Подготовка данных -------------------------------------------------------
# Выделим все возможные предикторы, включая географию и индексы трески
# Примечание: оставляем только числовые переменные, т.к. большинство моделей
# требует числовой вход без категориальных уровней.
predictors <- DATA %>% 
  select(-YEAR, -R3haddock) %>% 
  select_if(is.numeric) # Только числовые переменные

# Целевая переменная
response <- DATA$R3haddock

# В статистическом анализе мы различаем:
# - Отклик (response/target variable) - то, что мы пытаемся предсказать (в нашем случае R3haddock)
# - Предикторы (predictors/features) - переменные, которые могут объяснять изменения отклика
# Для корректного анализа важно, чтобы предикторы были числовыми или преобразованы в числовой формат.

# 2. Обработка пропусков -----------------------------------------------------
# Заполнение медианными значениями — простой и устойчивый способ справиться с NA.
# Альтернативы: множественная иммутация (mice), KNN-impute и др.
predictors_filled <- predictors %>%
  mutate(across(everything(), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))

# Заполнение медианой - простой и устойчивый метод обработки пропусков для числовых переменных.
# Медиана предпочтительнее среднего, так как менее чувствительна к выбросам.

# 3. Предварительный анализ корреляций ---------------------------------------
# Зачем: высокие корреляции затрудняют интерпретацию и могут вредить ряду моделей.
cor_matrix <- cor(predictors_filled, use = "complete.obs")
corrplot(cor_matrix, method = "circle", type = "upper", tl.cex = 0.7)

# Удаляем высокоскоррелированные предикторы (r > 0.8)
# Это механическое сокращение мультиколлинеарности до этапа отбора.
high_cor <- findCorrelation(cor_matrix, cutoff = 0.8)
predictors_filtered <- predictors_filled[, -high_cor]

# Высокая корреляция между предикторами (мультиколлинеарность) может привести к нестабильности моделей.
# Например, если два предиктора почти идентичны, модель может неустойчиво распределять их влияние на отклик.
# Удаление сильно коррелированных переменных (r > 0.8) помогает улучшить интерпретируемость и стабильность моделей.


# 4. Автоматизированный отбор Boruta (обертка Random Forest) -----------------
# Идея: определить признаки, которые важнее, чем случайный шум (shadow features).

# Визуализация результатов
plot(boruta_output, cex.axis = 0.7, las = 2)

boruta_stats <- attStats(boruta_output)
selected_vars <- getSelectedAttributes(boruta_output, withTentative = TRUE)

# Boruta - это алгоритм отбора признаков, основанный на методе случайного леса.
# Он сравнивает важность реальных переменных с "теневыми" переменными (случайными копиями),
# чтобы определить, действительно ли переменная информативна.
# Результаты Boruta показывают: 
#   - Confirmed (зеленые) - значимые предикторы
#   - Tentative (желтые) - предикторы, близкие к порогу значимости
#   - Rejected (красные) - незначимые предикторы


# 5. LASSO с более строгим критерием ------------------------------------------
# Идея: L1-регуляризация зануляет коэффициенты «слабых» предикторов.
# Выбор lambda.1se вместо lambda.min — более консервативный (простая модель).
x <- as.matrix(predictors_filtered)
y <- response

# LASSO (Least Absolute Shrinkage and Selection Operator) - метод регрессии с L1-регуляризацией,
# который одновременно выполняет отбор признаков и оценку коэффициентов. [[8]]
# Параметр lambda контролирует силу регуляризации:
#   - lambda.min дает наименьшую ошибку, но может включать шумовые переменные
#   - lambda.1se (на 1 стандартную ошибку больше) дает более простую модель с меньшим риском переобучения
# Для прогнозирования мы предпочитаем более строгий критерий (lambda.1se), чтобы модель была устойчивее. [[1]]

# Кросс-валидация
cv_fit <- cv.glmnet(x, y, alpha = 1, nfolds = 10)
plot(cv_fit)

# ИСПОЛЬЗУЕМ lambda.1se вместо lambda.min — СТРОЖЕ!
lasso_coef <- coef(cv_fit, s = "lambda.1se")  # <-- Ключевое изменение!
lasso_vars <- rownames(lasso_coef)[lasso_coef[,1] != 0][-1]  # исключаем (Intercept)


# 6. Сравнение отобранных предикторов ----------------------------------------
# Полезно видеть, какие признаки отмечают оба метода (устойчивые кандидаты).
cat("Boruta selected:", length(selected_vars), "variables\n")

Boruta selected: 3 variables

print(selected_vars)

[1] "codTSB" "T12"    "I5"    

cat("\nLASSO selected:", length(lasso_vars), "variables\n")

LASSO selected: 5 variables

print(lasso_vars)

[1] "codTSB" "T12"    "NAO3"   "NAO4"   "NAO5"  

# 7. Финальный набор предикторов (объединение результатов) -------------------
# Логика: объединяем списки, добавляем биологически важные переменные вручную.
final_vars <- union(selected_vars, lasso_vars) 

# Добавляем обязательные переменные по биологической логике
mandatory <- c("haddock68")
final_vars <- union(final_vars, mandatory) %>% unique()

# Мы объединяем результаты двух методов отбора признаков для большей надежности.
# Также добавляем переменную haddock68 (нерестовый запас), так как биологически 
# логично, что пополнение запаса напрямую зависит от численности производителей. 
# Это пример интеграции экспертных знаний в статистический анализ - важный принцип 
# при работе с данными в биологических науках.

# 8. Проверка значимости -----------------------------------------------------
# Быстрая оценка значимости с LM: не как окончательный вывод, а как sanity-check.
final_model <- lm(response ~ as.matrix(predictors_filled[, final_vars]))
summary(final_model)

Call:
lm(formula = response ~ as.matrix(predictors_filled[, final_vars]))

Residuals:
    Min      1Q  Median      3Q     Max 
-270986  -82376   -1037   98086  276129 

Coefficients:
                                                      Estimate Std. Error
(Intercept)                                         -1.082e+06  3.943e+05
as.matrix(predictors_filled[, final_vars])codTSB    -2.346e-01  5.536e-02
as.matrix(predictors_filled[, final_vars])T12        3.864e+05  7.198e+04
as.matrix(predictors_filled[, final_vars])I5        -1.825e+02  2.572e+03
as.matrix(predictors_filled[, final_vars])NAO3      -5.801e+04  3.129e+04
as.matrix(predictors_filled[, final_vars])NAO4       8.345e+04  3.035e+04
as.matrix(predictors_filled[, final_vars])NAO5      -7.278e+04  2.488e+04
as.matrix(predictors_filled[, final_vars])haddock68  1.232e-01  4.515e-01
                                                    t value Pr(>|t|)    
(Intercept)                                          -2.744 0.011305 *  
as.matrix(predictors_filled[, final_vars])codTSB     -4.238 0.000288 ***
as.matrix(predictors_filled[, final_vars])T12         5.368 1.64e-05 ***
as.matrix(predictors_filled[, final_vars])I5         -0.071 0.944028    
as.matrix(predictors_filled[, final_vars])NAO3       -1.854 0.076118 .  
as.matrix(predictors_filled[, final_vars])NAO4        2.750 0.011146 *  
as.matrix(predictors_filled[, final_vars])NAO5       -2.925 0.007412 ** 
as.matrix(predictors_filled[, final_vars])haddock68   0.273 0.787227    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 158600 on 24 degrees of freedom
Multiple R-squared:  0.7173,    Adjusted R-squared:  0.6348 
F-statistic: 8.698 on 7 and 24 DF,  p-value: 2.58e-05

# 9. Формирование финального датасета ----------------------------------------
# Собираем набор с откликом и выбранными предикторами; удалим строки с NA.
model_data <- DATA %>%
  select(R3haddock, all_of(final_vars)) %>%
  drop_na()

# Просмотр структуры финальных данных
glimpse(model_data)

Rows: 32
Columns: 8
$ R3haddock <dbl> 812363, 389416, 99474, 98946, 118812, 63028, 147657, 83270, ~
$ codTSB    <dbl> 913000, 1347064, 1687381, 2197863, 2112773, 1849957, 1697388~
$ T12       <dbl> 4.72, 4.66, 4.24, 3.90, 3.96, 4.27, 4.16, 4.07, 4.23, 5.08, ~
$ I5        <dbl> 43, 55, 26, 49, 56, 28, 52, 51, 69, 68, 41, 48, 50, 63, 40, ~
$ NAO3      <dbl> 1.46, -0.20, 0.87, 0.67, 1.26, 1.25, -0.24, 1.46, 0.87, 0.23~
$ NAO4      <dbl> 2.00, 0.29, 1.86, 0.97, 1.14, -0.85, -0.17, -1.02, -0.68, -0~
$ NAO5      <dbl> -1.53, 0.08, 2.63, -0.78, -0.57, -1.49, -1.06, -0.28, -1.32,~
$ haddock68 <dbl> 74586, 79205, 53195, 36337, 49122, 81514, 172177, 160886, 96~

# Визуализация важности переменных
# Внимание: важности от RF — относительные; сопоставляйте с предметной логикой.
var_importance <- randomForest(R3haddock ~ ., data = model_data, importance = TRUE)
varImpPlot(var_importance, main = "Важность предикторов")

# Перед окончательным выбором модели мы проверяем значимость предикторов с помощью линейной регрессии.
# Функция summary() показывает p-значения коэффициентов - если p < 0.05, переменная считается статистически значимой. 
# Визуализация важности переменных с помощью случайного леса дает дополнительную перспективу,
# показывая, какие переменные наиболее информативны для предсказания без предположений о линейности.

# ==============================================================================
#  ПОДГОТОВКА ДАННЫХ
# Создаём NAOspring, фиксируем финальный набор признаков, сохраняем CSV.
# ------------------------------------------------------------------------------
# Цель блока: стандартизировать набор признаков для дальнейшего сравнения
# моделей и обеспечить воспроизводимость (фиксированный CSV с нужными полями).
# ==============================================================================

# 1.1 Пакеты и окружение
# Примечание: блок повторяет базовую инициализацию для автономного запуска.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readxl, tidyverse, caret, corrplot)

rm(list = ls())
set.seed(123)
setwd("C:/RECRUITMENT/")

# 1.2 Загрузка исходных данных и приведение типов
DATA <- readxl::read_excel("RECRUITMENT.xlsx", sheet = "RECRUITMENT") %>%
  filter(YEAR > 1989 & YEAR < 2022) %>%
  mutate(
    across(starts_with("T"), as.numeric),
    across(starts_with("I"), as.numeric),
    across(starts_with("O"), as.numeric),
    across(where(is.character), ~na_if(., "NA"))
  )

# 1.3 Создаём NAOspring (если есть NAO3, NAO4, NAO5)
# Идея: агрегируем весенний индекс NAO как среднее за месяцы 3–5.
if (all(c("NAO3","NAO4","NAO5") %in% names(DATA))) {
  DATA <- DATA %>%
    mutate(NAOspring = rowMeans(pick(NAO3, NAO4, NAO5), na.rm = TRUE)) %>%
    select(-NAO3, -NAO4, -NAO5)
}

# NAO (North Atlantic Oscillation) - важный климатический индекс, влияющий описывающий изменения атмосферного давления
# над Северной Атлантикой. В частности, он отражает разницу в атмосферном давлении между Исландской депрессией и
# Азорским максимумом. NAO влияет на силу и направление западных ветров, а также на траектории штормов в Северной Атлантике. 
# Мы создаем NAOspring как среднее значение за весенние месяцы (марта, апреля, мая),
# так как именно в этот период происходят ключевые процессы, влияющие на нерест трески. 
# Создание составных переменных на основе экспертных знаний часто улучшает качество моделей.

# 1.4 Финальный учебный набор предикторов (фиксируем)
# Важно: проверяем присутствие нужных колонок и формируем компактный датасет.
needed <- c("codTSB", "T12", "I5", "NAOspring", "haddock68")
stopifnot(all(needed %in% names(DATA)))

# Сохраняем YEAR в CSV (ниже он будет отброшен при обучении, но нужен для графика)
model_data <- DATA %>%
  select(YEAR, all_of(needed), R3haddock) %>%
  drop_na()

write.csv(model_data, "selected_predictors_dataset.csv", row.names = FALSE)
glimpse(model_data)

Rows: 32
Columns: 7
$ YEAR      <dbl> 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, ~
$ codTSB    <dbl> 913000, 1347064, 1687381, 2197863, 2112773, 1849957, 1697388~
$ T12       <dbl> 4.72, 4.66, 4.24, 3.90, 3.96, 4.27, 4.16, 4.07, 4.23, 5.08, ~
$ I5        <dbl> 43, 55, 26, 49, 56, 28, 52, 51, 69, 68, 41, 48, 50, 63, 40, ~
$ NAOspring <dbl> 0.64333333, 0.05666667, 1.78666667, 0.28666667, 0.61000000, ~
$ haddock68 <dbl> 74586, 79205, 53195, 36337, 49122, 81514, 172177, 160886, 96~
$ R3haddock <dbl> 812363, 389416, 99474, 98946, 118812, 63028, 147657, 83270, ~

8.3 Модели «запас-пополнение» Рикера и Бивертона-Холта

Модели запас-пополнение представляют собой фундаментальный инструмент в оценке водных биоресурсов, которые гораздо больше, чем просто математические кривые, — это формализованные выражения фундаментальных биологических представлений о том, как численность родительского стада определяет успех следующего поколения. Среди классических моделей этого типа наиболее широко используются модель Рикера и модель Бивертона-Холта, каждая из которых отражает различные гипотезы о биологических процессах, происходящих в популяции. Модель Рикера, предложенная Уильямом Рикером в 1954 году и имеющая характерный горб на графике, выражается уравнением R = a*S*exp(-b*S), где R обозначает пополнение, S — нерестовый запас, а параметры a и b имеют четкую биологическую интерпретацию: a соответствует максимальной продуктивности на единицу запаса при очень низких плотностях, фактически отражая максимальное пополнение (количество рекрутов) на одного производителя, а b характеризует степень плотностной зависимости, определяющей точку, после которой начинается снижение из-за внутривидовой конкуренции. Эта модель предсказывает, что с ростом нерестового запаса пополнение сначала увеличивается, достигает максимума, а затем снижается, что отражает явление перенаселенности, когда чрезмерная плотность производителей приводит к конкуренции за ресурсы, нехватке корма для личинок, усилению каннибализма или даже эпидемиям, что в итоге снижает выход молоди — мы буквально видим, как чрезмерный успех закладывает семена будущего коллапса пополнения. В отличие от нее, модель Бивертона-Холта, разработанная в 1957 году, имеет вид R = aS/(1+b*S) и предполагает, что пополнение асимптотически приближается к предельному значению a/b, называемому Rmax, с увеличением нерестового запаса, без последующего снижения, что соответствует ситуации, когда основной лимитирующий фактор — это не внутривидовая конкуренция, а внешние условия: ограниченное количество нерестовых площадок, хищничество, которое не зависит от плотности, или просто конечная пропускная способность экосистемы для молоди. Эта модель идеально описывает сценарий, когда кривая плавно выходит на плато, символизируя насыщение, и представляет собой альтернативную логику, где главным ограничивающим фактором являются внешние, а не внутривидовые процессы.

При оценке параметров этих нелинейных моделей мы сталкиваемся с необходимостью применения специализированных методов, поскольку обычный метод наименьших квадратов не справляется с их сложной структурой; в нашем анализе мы используем улучшенный алгоритм nlsLM из пакета minpack.lm, который сочетает метод Левенберга-Марквардта с возможностью наложения ограничений на параметры, что важно для обеспечения биологической правдоподобности результатов, так как параметры a и b должны оставаться положительными. Для получения надежных начальных оценок параметров в модели Рикера мы применяем функцию srStarts из пакета FSA, которая автоматически определяет разумные стартовые значения на основе анализа данных, тогда как для модели Бивертона-Холта мы используем комбинацию автоматических и ручных подходов, оценивая a как среднее отношение R/S при низких значениях запаса и устанавливая разумные начальные значения для b с последующей защитой от некорректных значений. Однако подбор модели — это только полдела, и критически важно провести тщательную диагностику, поскольку самая большая ошибка — слепо применять эти модели, не задумываясь об их предпосылках. Мы строим график остатков, потому что любая закономерность в их распределении — это сигнал о том, что модель не уловила какой-то важный процесс в данных. Мы смотрим на доверительные интервалы параметров; если они невероятно широки, значит, наша модель перепараметризована для имеющихся данных, и её прогностическая сила будет сомнительной. Модель Рикера не будет работать, если в вашей системе нет механизма перенаселения, а модель Бивертона-Холта окажется бесполезной, если пополнение продолжает расти или, наоборот, обрушивается после достижения пика. Именно поэтому мы всегда начинаем с простого графика «запас-пополнение» — его форма сама подскажет, какая из концепций более адекватна для конкретной популяции.

Но реальный мир часто бывает сложнее этих двух идеализированных сценариев. Что если система ведёт себя по-рикеровски при высокой численности, но при низкой — работает иначе? Здесь на помощь приходит модификация — модель Рикера с порогом, или hockey-stick модель, которая сочетает в себе линейный рост при малых запасах и плато или спад при высоких, что может быть биологически более оправдано для многих запасов, находящихся под прессом промысла. И здесь мы подходим к самому главному — интеграции классики и современности. Эти модели не являются застывшими реликтами, а служат мощным инструментом для создания гипотез. Если модель Рикера плохо описывает данные, особенно в области низких значений запаса, это прямой сигнал о том, что возможно, существует какой-то дополнительный лимитирующий фактор, не учтенный в модели. Возможно, это температура воды на ключевой стадии развития икры, наличие хищников или доступность корма. Таким образом, классические модели становятся трамплином для более сложного анализа, включающего средовые предикторы. Мы можем включить параметры модели Рикера в качестве фиксированных эффектов в GAM или использовать предсказания классической модели в качестве одного из входных признаков для Random Forest. Этот синтез позволяет нам сохранить биологическую интерпретируемость классических моделей и добавить к ним гибкость и прогностическую силу машинного обучения для учета сложных, нелинейных влияний окружающей среды. В сущности, мы строим мост между глубоким, но узким знанием, заключенным в одной кривой, и широким, но зачастую “черно-ящичным” прогнозом сложного алгоритма, пытаясь получить лучшее из двух миров. Среди распространенных подводных камней при работе с моделями запас-пополнение следует отметить высокую чувствительность к начальным значениям параметров, что может приводить к сходимости к локальным минимумам, необходимость учета неоднородности дисперсии ошибок, особенно при работе с данными, охватывающими широкий диапазон значений запаса, и влияние временных лагов, поскольку пополнение в текущем году может зависеть не только от нерестового запаса в том же году, но и от условий прошлых лет. Кроме того, чистые модели запас-пополнение часто оказываются недостаточными для точного прогнозирования, так как пополнение зависит не только от размера нерестового запаса, но и от множества экологических факторов, что делает целесообразным развитие этих моделей в направлении включения дополнительных предикторов, как это продемонстрировано в последующих разделах нашего анализа. Тем не менее, классические модели Рикера и Бивертона-Холта остаются важной отправной точкой в анализе динамики рыбных популяций, предоставляя интерпретируемую основу для понимания механизмов регулирования численности и служа эталоном для оценки добавленной ценности более сложных моделей, что особенно важно в условиях ограниченных данных, характерных для многих водных экосистем.

# ==============================================================================
# ПРАКТИЧЕСКОЕ ЗАНЯТИЕ: АНАЛИЗ ФАКТОРОВ, ВЛИЯЮЩИХ НА ПОПОЛНЕНИЕ 
# (КЛАССИЧЕСКИЕ МОДЕЛИ ЗАПАС-ПОПОЛНЕНИЕ)
# Курс: "Оценка водных биоресурсов в среде R (для начинающих)"
# ==============================================================================

# Установка и подключение ТОЛЬКО необходимых библиотек
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,    # Манипуляции с данными и визуализация
  FSA,          # Начальные оценки для моделей запас-пополнение
  minpack.lm,   # Улучшенный алгоритм нелинейной регрессии (nlsLM)
  car,          # Проверка допущений моделей
  mgcv,         # Построение GAM-моделей
  investr,      # Доверительные интервалы для нелинейных моделей
   caret)       # Расчет RMSE

# Очистка среды и установка рабочей директории
rm(list = ls())
setwd("C:/RECRUITMENT/")

# 1. ЗАГРУЗКА ДАННЫХ -----------------------------------------------------------
model_data <- read.csv("selected_predictors_dataset.csv", 
                      header = TRUE, 
                      stringsAsFactors = FALSE)

# Проверка структуры данных
str(model_data)

'data.frame':   32 obs. of  7 variables:
 $ YEAR     : int  1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
 $ codTSB   : int  913000 1347064 1687381 2197863 2112773 1849957 1697388 1537459 1350918 1199169 ...
 $ T12      : num  4.72 4.66 4.24 3.9 3.96 4.27 4.16 4.07 4.23 5.08 ...
 $ I5       : int  43 55 26 49 56 28 52 51 69 68 ...
 $ NAOspring: num  0.6433 0.0567 1.7867 0.2867 0.61 ...
 $ haddock68: int  74586 79205 53195 36337 49122 81514 172177 160886 96380 37977 ...
 $ R3haddock: int  812363 389416 99474 98946 118812 63028 147657 83270 359701 386866 ...

# Проверка на пропущенные значения (должно быть 0 после предобработки)
sum(is.na(model_data))

[1] 0

# 2. ПОДГОТОВКА ДАННЫХ ДЛЯ МОДЕЛЕЙ ЗАПАС-ПОПОЛНЕНИЕ ----------------------------
rec_data <- data.frame(
  S = model_data$haddock68,  # Нерестовый запас
  R = model_data$R3haddock   # Пополнение
)

# 3. ПОДГОНКА МОДЕЛИ РИКЕРА ----------------------------------------------------
ricker_starts <- FSA::srStarts(R ~ S, data = rec_data, type = "Ricker")
ricker_model <- minpack.lm::nlsLM(
  R ~ a * S * exp(-b * S),
  data = rec_data,
  start = ricker_starts,
  lower = c(a = 0, b = 0)
)

# 4. ПОДГОНКА МОДЕЛИ БИВЕРТОНА-ХОЛТА -------------------------------------------
a_start <- mean(rec_data$R[rec_data$S < quantile(rec_data$S, 0.25)] / 
                rec_data$S[rec_data$S < quantile(rec_data$S, 0.25)], na.rm = TRUE)

if (is.na(a_start) || a_start <= 0) a_start <- 0.001

bh_model <- minpack.lm::nlsLM(
  R ~ (a * S) / (1 + b * S),
  data = rec_data,
  start = list(a = a_start, b = 0.0001),
  lower = c(a = 0.0001, b = 0.00001),
  control = nls.lm.control(maxiter = 200)
)

# 5. ОЦЕНКА КАЧЕСТВА МОДЕЛЕЙ ---------------------------------------------------

calculate_R2 <- function(model, data) {
  predicted <- predict(model, newdata = data)
  residuals <- data$R - predicted
  SSE <- sum(residuals^2)
  SST <- sum((data$R - mean(data$R))^2)
  R2 <- 1 - (SSE / SST)
  n <- nrow(data)
  p <- length(coef(model))
  adj_R2 <- 1 - ((n - 1) / (n - p - 1)) * (1 - R2)
  return(list(R2 = R2, adj_R2 = adj_R2))
}

calculate_pvalue <- function(model, data) {
  predicted <- predict(model, newdata = data)
  residuals <- data$R - predicted
  SSE <- sum(residuals^2)
  SST <- sum((data$R - mean(data$R))^2)
  SSR <- SST - SSE
  n <- nrow(data)
  p <- length(coef(model))
  F_stat <- (SSR / (p - 1)) / (SSE / (n - p))
  p_value <- pf(F_stat, df1 = p - 1, df2 = n - p, lower.tail = FALSE)
  return(p_value)
}

ricker_r2 <- calculate_R2(ricker_model, rec_data)
bh_r2 <- calculate_R2(bh_model, rec_data)

ricker_p <- calculate_pvalue(ricker_model, rec_data)
bh_p <- calculate_pvalue(bh_model, rec_data)

cat("AIC Рикера:", AIC(ricker_model), "\n")

AIC Рикера: 891.9919 

cat("AIC Бивертона-Холта:", AIC(bh_model), "\n")

AIC Бивертона-Холта: 894.2029 

# 6. ВИЗУАЛИЗАЦИЯ РЕЗУЛЬТАТОВ --------------------------------------------------
new_data <- data.frame(S = seq(min(rec_data$S), max(rec_data$S), length.out = 100))
ricker_ci <- investr::predFit(ricker_model, newdata = new_data, interval = "confidence")
bh_ci <- investr::predFit(bh_model, newdata = new_data, interval = "confidence")

plot_data <- new_data %>%
  mutate(
    ricker_pred = predict(ricker_model, newdata = .),
    ricker_lwr = ricker_ci[, "lwr"],
    ricker_upr = ricker_ci[, "upr"],
    bh_pred = predict(bh_model, newdata = .),
    bh_lwr = bh_ci[, "lwr"],
    bh_upr = bh_ci[, "upr"]
  )

ggplot() +
  geom_point(data = rec_data, aes(x = S, y = R), 
             color = "darkgray", size = 3, alpha = 0.7) +
  geom_ribbon(data = plot_data, aes(x = S, ymin = ricker_lwr, ymax = ricker_upr), 
              fill = "red", alpha = 0.2) +
  geom_ribbon(data = plot_data, aes(x = S, ymin = bh_lwr, ymax = bh_upr), 
              fill = "blue", alpha = 0.2) +
  geom_line(data = plot_data, aes(x = S, y = ricker_pred), 
            color = "red", linewidth = 1.2) +
  geom_line(data = plot_data, aes(x = S, y = bh_pred), 
            color = "blue", linewidth = 1.2, linetype = "dashed") +
  labs(
    title = "Сравнение моделей запас-пополнение",
    subtitle = paste0(
      "Рикер: R² = ", round(ricker_r2$R2, 2), ", p = ", format.pval(ricker_p, digits = 3),
      " | Бивертон-Холт: R² = ", round(bh_r2$R2, 2), ", p = ", format.pval(bh_p, digits = 3)
    ),
    x = "Нерестовый запас",
    y = "Пополнение"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none"
  )

# 7. СРАВНЕНИЕ С ДРУГИМИ ТИПАМИ МОДЕЛЕЙ ----------------------------------------

lm_model <- lm(R3haddock ~ ., data = model_data)
glm_model <- glm(R3haddock ~ ., family = Gamma(link = "log"), data = model_data)
gam_model <- mgcv::gam(R3haddock ~ s(haddock68) + s(codTSB) + s(T12) + s(I5) + s(NAOspring),
               data = model_data, method = "REML")

model_comparison <- data.frame(
  Model = c("Рикер", "Бивертон-Холт", "LM", "GLM", "GAM"),
  AIC = c(AIC(ricker_model), AIC(bh_model), AIC(lm_model), AIC(glm_model), AIC(gam_model)),
  R2 = c(
    ricker_r2$R2, 
    bh_r2$R2, 
    summary(lm_model)$r.squared,
    cor(model_data$R3haddock, predict(glm_model, type = "response"))^2,
    summary(gam_model)$r.sq
  )
)

print(model_comparison)

          Model      AIC          R2
1         Рикер 891.9919 0.072305265
2 Бивертон-Холт 894.2029 0.005938625
3            LM 877.0895 0.573972535
4           GLM 857.0346 0.566389771
5           GAM 862.9064 0.738660678

# 8. ИНТЕРПРЕТАЦИЯ РЕЗУЛЬТАТОВ -------------------------------------------------
cat("\nПараметры модели Рикера:")

Параметры модели Рикера:

cat("\na =", coef(ricker_model)[1], "- максимальная продукция потомства")

a = 7.931302 - максимальная продукция потомства

cat("\nb =", coef(ricker_model)[2], "- коэффициент плотностной зависимости")

b = 7.4007e-06 - коэффициент плотностной зависимости

cat("\n\nПараметры модели Бивертона-Холта:")

Параметры модели Бивертона-Холта:

cat("\na =", coef(bh_model)[1], "- максимальное пополнение на особь")

a = 43.77667 - максимальное пополнение на особь

cat("\nb =", coef(bh_model)[2], "- коэффициент внутривидовой конкуренции")

b = 0.0001247403 - коэффициент внутривидовой конкуренции

# ==============================================================================
# 7. СРАВНЕНИЕ С ДРУГИМИ ТИПАМИ МОДЕЛЕЙ ----------------------------------------

# Построение линейной модели LM
lm_model <- lm(R3haddock ~ ., data = model_data)

# Диагностика
par(mfrow = c(2, 2))
plot(lm_model)

vif(lm_model)  # Проверка мультиколлинеарности

     YEAR    codTSB       T12        I5 NAOspring haddock68 
 2.925549  3.439340  1.970023  1.332529  1.181897  2.522569 

# Интерпретация
summary(lm_model)

Call:
lm(formula = R3haddock ~ ., data = model_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-279162 -111056  -35757  141083  324173 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.448e+07  1.224e+07   1.183 0.248123    
YEAR        -7.863e+03  6.248e+03  -1.258 0.219869    
codTSB      -1.888e-01  7.817e-02  -2.416 0.023338 *  
T12          4.339e+05  9.738e+04   4.456 0.000153 ***
I5          -2.568e+03  3.058e+03  -0.840 0.409036    
NAOspring   -7.666e+04  5.735e+04  -1.337 0.193325    
haddock68    2.782e-01  5.485e-01   0.507 0.616444    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 190800 on 25 degrees of freedom
Multiple R-squared:  0.574, Adjusted R-squared:  0.4717 
F-statistic: 5.614 on 6 and 25 DF,  p-value: 0.0008393

# Построение обобщенной линейной модели GLM
glm_model <- glm(R3haddock ~ ., 
                family = Gamma(link = "log"), 
                data = model_data)
summary(glm_model)

Call:
glm(formula = R3haddock ~ ., family = Gamma(link = "log"), data = model_data)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.375e+01  3.667e+01   0.648   0.5231    
YEAR        -8.601e-03  1.872e-02  -0.459   0.6499    
codTSB      -5.945e-07  2.342e-07  -2.539   0.0177 *  
T12          1.411e+00  2.917e-01   4.837 5.68e-05 ***
I5           7.430e-03  9.161e-03   0.811   0.4250    
NAOspring    1.508e-02  1.718e-01   0.088   0.9307    
haddock68    9.422e-07  1.643e-06   0.573   0.5715    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.326724)

    Null deviance: 19.8880  on 31  degrees of freedom
Residual deviance:  8.4535  on 25  degrees of freedom
AIC: 857.03

Number of Fisher Scoring iterations: 9

# Построение обобщенной аддитивной модели GАM
library(mgcv)
gam_model <- gam(R3haddock ~ 
                 s(codTSB) + 
                 s(T12) + 
                 s(I5) + 
                 s(NAOspring) + 
                 s(haddock68),
               data = model_data,
               method = "REML")
summary(gam_model)

Family: gaussian 
Link function: identity 

Formula:
R3haddock ~ s(codTSB) + s(T12) + s(I5) + s(NAOspring) + s(haddock68)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   320163      23724   13.49 4.37e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value   
s(codTSB)    2.354  2.899 1.684 0.17581   
s(T12)       2.190  2.676 5.908 0.00357 **
s(I5)        4.642  5.539 1.518 0.15004   
s(NAOspring) 1.293  1.503 1.154 0.22854   
s(haddock68) 1.824  2.200 0.629 0.65124   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.739   Deviance explained = 84.2%
-REML = 361.62  Scale est. = 1.8011e+10  n = 32

plot(gam_model, pages = 1, residuals = TRUE)

# Таблица сравнения моделей
# Сравнение моделей
model_comparison <- data.frame(
  Model = c("Рикер", "Бивертон-Холт", "LM", "GLM", "GAM"),
  AIC = c(AIC(ricker_model), AIC(bh_model), AIC(lm_model), AIC(glm_model), AIC(gam_model)),
  R2 = c(ricker_r2$R2, bh_r2$R2, summary(lm_model)$r.squared, 
         cor(model_data$R3haddock, predict(glm_model))^2, 
         summary(gam_model)$r.sq),  # Используем summary(gam_model)$r.sq для R^2
  Adj_R2 = c(ricker_r2$adj_R2, bh_r2$adj_R2, summary(lm_model)$adj.r.squared, NA, 
             summary(gam_model)$r.sq),  # Используем summary(gam_model)$r.sq для Adjusted R^2
  RMSE = c(RMSE(predict(ricker_model), rec_data$R), 
           RMSE(predict(bh_model), rec_data$R),
           RMSE(predict(lm_model), model_data$R3haddock),
           RMSE(predict(glm_model, type = "response"), model_data$R3haddock),
           RMSE(predict(gam_model, type = "response"), model_data$R3haddock))
)

# Вывод таблицы
print(model_comparison)

          Model      AIC          R2       Adj_R2     RMSE
1         Рикер 891.9919 0.072305265  0.008326318 248869.6
2 Бивертон-Холт 894.2029 0.005938625 -0.062617332 257617.8
3            LM 877.0895 0.573972535  0.471725944 168650.7
4           GLM 857.0346 0.502625978           NA 170386.6
5           GAM 862.9064 0.738660678  0.738660678 102584.0

# ==============================================================================
# ВИЗУАЛИЗАЦИЯ ВСЕХ МОДЕЛЕЙ НА ОДНОМ ГРАФИКЕ 
# ==============================================================================

# Фиксируем другие предикторы на их средних значениях (исключая haddock68)
mean_values <- model_data %>%
  select(-R3haddock, -haddock68) %>%
  summarise(across(everything(), ~ mean(.x, na.rm = TRUE)))

# Расширяем new_data, добавляя средние значения других предикторов
new_data_full <- new_data %>%
  bind_cols(mean_values[rep(1, nrow(new_data)), ]) %>%
  rename(haddock68 = S)  # Переименовываем S в haddock68 для совместимости

# Получаем предсказания для всех моделей
new_data_full <- new_data_full %>%
  mutate(
    # Предсказания для моделей запаса-пополнения
    ricker_pred = predict(ricker_model, newdata = data.frame(S = haddock68)),
    bh_pred = predict(bh_model, newdata = data.frame(S = haddock68)),
    
    # Предсказания для линейной модели (LM)
    lm_pred = predict(lm_model, newdata = .),
    
    # Предсказания для обобщенной линейной модели (GLM)
    glm_pred = predict(glm_model, newdata = ., type = "response"),
    
    # Предсказания для обобщенной аддитивной модели (GAM)
    gam_pred = predict(gam_model, newdata = ., type = "response")
  )

# Создаем длинный формат данных для ggplot
plot_data <- new_data_full %>%
  select(haddock68, ricker_pred, bh_pred, lm_pred, glm_pred, gam_pred) %>%
  pivot_longer(
    cols = -haddock68,
    names_to = "model",
    values_to = "prediction"
  ) %>%
  mutate(
    model = case_when(
      model == "ricker_pred" ~ "Рикер",
      model == "bh_pred" ~ "Бивертон-Холт",
      model == "lm_pred" ~ "Линейная (LM)",
      model == "glm_pred" ~ "Обобщенная линейная (GLM)",
      model == "gam_pred" ~ "Обобщенная аддитивная (GAM)",
      TRUE ~ model
    )
  )

# Создаем палитру цветов для моделей
model_colors <- c(
  "Рикер" = "#E41A1C",          # Красный
  "Бивертон-Холт" = "#377EB8",  # Синий
  "Линейная (LM)" = "#4DAF4A",  # Зеленый
  "Обобщенная линейная (GLM)" = "#984EA3", # Фиолетовый
  "Обобщенная аддитивная (GAM)" = "#FF7F00" # Оранжевый
)

# Создаем график
ggplot() +
  # Точки исходных данных
  geom_point(data = rec_data, aes(x = S, y = R), 
             color = "darkgray", size = 2.5, alpha = 0.7) +
  
  # Линии предсказаний моделей
  geom_line(data = plot_data, 
            aes(x = haddock68, y = prediction, color = model, linetype = model),
            linewidth = 1.2) +
  
  # Настройка цветов и типов линий
  scale_color_manual(values = model_colors) +
  scale_linetype_manual(values = c(
    "Рикер" = "solid",
    "Бивертон-Холт" = "dashed",
    "Линейная (LM)" = "dotdash",
    "Обобщенная линейная (GLM)" = "longdash",
    "Обобщенная аддитивная (GAM)" = "twodash"
  )) +
  
  # Подписи и темы
  labs(
    title = "Сравнение моделей зависимости пополнения от нерестового запаса",
    subtitle = "Фиксация других предикторов на средних значениях",
    x = "Нерестовый запас (тыс. тонн)",
    y = "Пополнение (млн особей)",
    color = "Модель",
    linetype = "Модель"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray30"),
    axis.title = element_text(size = 12),
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "gray80", fill = NA, linewidth = 0.5)
  ) +
  guides(
    color = guide_legend(nrow = 2, byrow = TRUE),
    linetype = guide_legend(nrow = 2, byrow = TRUE)
  )

8.4 Статистические модели LM/GLM/GAM

Статистические модели линейной регрессии (LM), обобщенной линейной регрессии (GLM) и обобщенной аддитивной регрессии (GAM) представляют собой мощный и взаимодополняющий набор инструментов для анализа водных биоресурсов, позволяющий исследователям от простых линейных зависимостей переходить к сложным нелинейным взаимодействиям, сохраняя при этом интерпретируемость результатов. Линейная модель (LM) служит фундаментом для всего статистического анализа в гидробиологии, основываясь на предположении, что зависимость между предикторами и откликом является линейной, а остатки распределены нормально с постоянной дисперсией. Эта модель предоставляет простую интерпретацию коэффициентов как величины изменения отклика при единичном изменении предиктора, что особенно ценно при работе с такими биологическими показателями, как пополнение запаса или нерестовая биомасса. Однако при анализе водных биоресурсов мы часто сталкиваемся с данными, которые нарушают ключевые предположения LM: пополнение рыбы или беспозвоночного не может быть отрицательным, его распределение обычно сильно скошено вправо, а дисперсия часто увеличивается с ростом среднего значения. Именно здесь на помощь приходит обобщенная линейная модель (GLM), расширяющая возможности LM за счет введения двух ключевых компонентов — экспоненциального семейства распределений и связующей функции (link-function). Для данных о рыбных запасах особенно полезно Gamma-распределение с логарифмической связкой, которое учитывает положительность отклика и мультипликативную природу ошибок, характерную для биологических данных. В отличие от LM, где мы интерпретируем коэффициенты как абсолютные изменения, в GLM с лог-связкой коэффициенты отражают относительные изменения: увеличение предиктора на единицу приводит к умножению ожидаемого отклика на exp(коэффициент), что соответствует биологической реальности, где эффекты часто действуют мультипликативно, а не аддитивно. Но даже GLM сохраняет ограничение на линейность в преобразованном пространстве, что может быть недостаточным для описания сложных экологических зависимостей, таких как оптимальный диапазон температуры для нереста или пороговые эффекты средовых факторов. Здесь в игру вступают обобщенные аддитивные модели (GAM), которые заменяют линейные комбинации предикторов на гладкие функции, оцениваемые с помощью сплайнов, что позволяет моделировать практически любые нелинейные зависимости без предварительного задания их формы. GAM сохраняет интерпретируемость линейных моделей, так как каждая гладкая функция может быть визуализирована и проанализирована отдельно, показывая, как именно каждый фактор влияет на пополнение запаса, будь то монотонный рост, оптимум с максимумом или сложная колебательная зависимость. При работе с GAM особое внимание уделяется выбору степени гладкости, так как чрезмерно гибкие функции могут переобучиться на шум в данных, тогда как недостаточно гибкие не уловят реальные биологические закономерности; в пакете mgcv это решается автоматически через метод максимального правдоподобия с штрафом (REML), который балансирует качество подгонки и гладкость функций. Сравнивая эти модели с классическими моделями запас-пополнение, мы видим, что GAM может рассматриваться как их естественное обобщение: вместо фиксированной формы кривой Рикера или Бивертона-Холта GAM позволяет данным “говорить за себя”, выявляя оптимальную форму зависимости без предварительных гипотез, при этом сохраняя возможность включить нерестовый запас как один из гладких членов в модель, дополненный другими экологическими факторами. Однако при всей своей гибкости, GAM, как и LM с GLM, требует тщательной проверки предположений: мы анализируем графики остатков против предсказанных значений, чтобы убедиться в отсутствии систематических отклонений, проверяем нормальность остатков (для LM) или соответствие выбранному распределению (для GLM/GAM), и исследуем влияние влиятельных точек, которые могут исказить результаты, особенно в условиях ограниченных данных, характерных для гидробиологических исследований. Выбор между LM, GLM и GAM должен основываться не только на статистических критериях, таких как AIC или кросс-валидация, но и на биологической интерпретируемости результатов: иногда более простая модель с меньшей точностью предпочтительнее сложного “черного ящика”, особенно когда результаты должны быть понятны начальникам и менеджерам рыболовства. Практический подход, который обычно рекомендуется начинающим ихтиологам/гидробиологам, состоит в последовательном усложнении модели: начните с классической модели запас-пополнение, затем добавьте средовые факторы через LM/GLM, и только если зависимости явно нелинейны, перейдите к GAM, всегда проверяя, действительно ли усложнение модели приводит к биологически значимому улучшению понимания процесса. Важно помнить, что статистическая модель — это не самоцель, а инструмент для понимания биологических процессов, и даже самая изощренная модель бесполезна, если её результаты нельзя перевести на язык биологии и применить для устойчивого управления водными ресурсами. В конечном счете, сочетание классических представлений об экосистемах с современными статистическими методами позволяет нам строить мост между фундаментальной биологией и прикладной оценкой запасов, где каждая модель, от простой линейной регрессии до сложного GAM, вносит свой вклад в формирование целостного понимания динамики водных биоресурсов.

# ==============================================================================
# Версия: только LM / GLM(Gamma) / GAM
# Без caret/train: стандартная оценка параметров lm/glm/gam, собственная time-slice CV,
# выбор лучшей модели, прогноз 2022–2024, эмпирические интервалы и график.
# ==============================================================================

# 0) Пакеты и окружение --------------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl, tidyverse, mgcv, lmtest, car, ggplot2, corrplot
)

rm(list = ls())
set.seed(123)
setwd("C:/RECRUITMENT/")  # при необходимости измените путь


# 1) Подготовка данных ---------------------------------------------------------
# Загрузка, приведение типов, создание NAOspring, фиксируем набор признаков
DATA <- readxl::read_excel("RECRUITMENT.xlsx", sheet = "RECRUITMENT") %>%
  filter(YEAR > 1989 & YEAR < 2022) %>%
  mutate(
    across(starts_with("T"), as.numeric),
    across(starts_with("I"), as.numeric),
    across(starts_with("O"), as.numeric),
    across(where(is.character), ~na_if(., "NA"))
  )

if (all(c("NAO3","NAO4","NAO5") %in% names(DATA))) {
  DATA <- DATA %>%
    mutate(NAOspring = rowMeans(pick(NAO3, NAO4, NAO5), na.rm = TRUE)) %>%
    select(-NAO3, -NAO4, -NAO5)
}

needed <- c("codTSB", "T12", "I5", "NAOspring", "haddock68")
stopifnot(all(needed %in% names(DATA)))

model_data <- DATA %>%
  select(YEAR, all_of(needed), R3haddock) %>%
  drop_na() %>%
  arrange(YEAR)

write.csv(model_data, "selected_predictors_dataset.csv", row.names = FALSE)
glimpse(model_data)

Rows: 32
Columns: 7
$ YEAR      <dbl> 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, ~
$ codTSB    <dbl> 913000, 1347064, 1687381, 2197863, 2112773, 1849957, 1697388~
$ T12       <dbl> 4.72, 4.66, 4.24, 3.90, 3.96, 4.27, 4.16, 4.07, 4.23, 5.08, ~
$ I5        <dbl> 43, 55, 26, 49, 56, 28, 52, 51, 69, 68, 41, 48, 50, 63, 40, ~
$ NAOspring <dbl> 0.64333333, 0.05666667, 1.78666667, 0.28666667, 0.61000000, ~
$ haddock68 <dbl> 74586, 79205, 53195, 36337, 49122, 81514, 172177, 160886, 96~
$ R3haddock <dbl> 812363, 389416, 99474, 98946, 118812, 63028, 147657, 83270, ~

# 2) Формулы моделей и вспомогательные функции --------------------------------
f_lm  <- as.formula("R3haddock ~ codTSB + T12 + I5 + NAOspring + haddock68")
f_gam <- as.formula("R3haddock ~ s(codTSB,bs='tp',k=5) + s(T12,bs='tp',k=5) + s(I5,bs='tp',k=5) + s(NAOspring,bs='tp',k=5) + s(haddock68,bs='tp',k=5)")

rmse <- function(a, p) sqrt(mean((a - p)^2, na.rm = TRUE))
mae  <- function(a, p) mean(abs(a - p), na.rm = TRUE)
r2   <- function(a, p) 1 - sum((a - p)^2, na.rm = TRUE) / sum((a - mean(a))^2, na.rm = TRUE)

safe_fit <- function(expr) {
  out <- try(eval(expr), silent = TRUE)
  if (inherits(out, "try-error")) NULL else out
}


# 3) Time-slice CV (expanding window, h=3) + хронологический тест -------------
md <- model_data
md_for_fit <- md %>% select(codTSB, T12, I5, NAOspring, haddock68, R3haddock)

n <- nrow(md_for_fit)
holdout_frac <- 0.2
n_test <- max(4, ceiling(n * holdout_frac))
train_ts <- head(md_for_fit, n - n_test)
test_ts  <- tail(md_for_fit, n_test)

n_train <- nrow(train_ts)
initial_frac <- 0.6
horizon      <- 3
initialWindow <- max(10, floor(initial_frac * n_train))
if (initialWindow + horizon > n_train) initialWindow <- n_train - horizon

# Аккумулируем метрики и остатки по срезам
cv_summ <- tibble(Model = character(), RMSE = double(), MAE = double())
resids_cv <- list(LM = c(), GLM = c(), GAM = c())

slice_id <- 0
for (i in seq(initialWindow, n_train - horizon)) {
  slice_id <- slice_id + 1
  idx_tr <- 1:i
  idx_te <- (i+1):(i+horizon)
  dtr <- train_ts[idx_tr, ]
  dte <- train_ts[idx_te, ]

  # LM
  lm_fit <- safe_fit(quote(lm(f_lm, data = dtr)))
  if (!is.null(lm_fit)) {
    pr <- try(predict(lm_fit, newdata = dte), silent = TRUE)
    if (!inherits(pr, "try-error")) {
      cv_summ <- add_row(cv_summ, Model = "LM", RMSE = rmse(dte$R3haddock, pr), MAE = mae(dte$R3haddock, pr))
      resids_cv$LM <- c(resids_cv$LM, dte$R3haddock - pr)
    }
  }

  # GLM (Gamma)
  glm_fit <- safe_fit(quote(glm(f_lm, data = dtr, family = Gamma(link = "log"))))
  if (!is.null(glm_fit)) {
    pr <- try(predict(glm_fit, newdata = dte, type = "response"), silent = TRUE)
    if (!inherits(pr, "try-error")) {
      cv_summ <- add_row(cv_summ, Model = "GLM", RMSE = rmse(dte$R3haddock, pr), MAE = mae(dte$R3haddock, pr))
      resids_cv$GLM <- c(resids_cv$GLM, dte$R3haddock - pr)
    }
  }

  # GAM (Gamma log), ограничиваем сложность k для стабильности на малом n
  gam_fit <- safe_fit(quote(mgcv::gam(f_gam, data = dtr, family = Gamma(link = "log"), method = "REML", select = TRUE)))
  if (!is.null(gam_fit)) {
    pr <- try(predict(gam_fit, newdata = dte, type = "response"), silent = TRUE)
    if (!inherits(pr, "try-error")) {
      cv_summ <- add_row(cv_summ, Model = "GAM", RMSE = rmse(dte$R3haddock, pr), MAE = mae(dte$R3haddock, pr))
      resids_cv$GAM <- c(resids_cv$GAM, dte$R3haddock - pr)
    }
  }
}

# Средние метрики по моделям
cv_rank <- cv_summ %>% group_by(Model) %>% summarise(RMSE = mean(RMSE, na.rm = TRUE), MAE = mean(MAE, na.rm = TRUE), .groups = "drop") %>% arrange(RMSE, MAE)
print(cv_rank)

# A tibble: 3 x 3
  Model    RMSE     MAE
  <chr>   <dbl>   <dbl>
1 GLM   280259. 237187.
2 LM    370298. 340884.
3 GAM   504613. 432062.

best_model_name <- cv_rank$Model[1]
cat(sprintf("\nЛучшая модель по time-slice CV: %s\n", best_model_name))

Лучшая модель по time-slice CV: GLM

# Хронологический тест: обучаем на всём train_ts, прогнозируем на test_ts
fit_on <- function(model_name, data) {
  if (model_name == "LM") return(lm(f_lm, data = data))
  if (model_name == "GLM") return(glm(f_lm, data = data, family = Gamma(link = "log")))
  mgcv::gam(f_gam, data = data, family = Gamma(link = "log"), method = "REML", select = TRUE)
}

predict_on <- function(fit, newdata, model_name) {
  if (model_name == "GLM") return(predict(fit, newdata = newdata, type = "response"))
  if (inherits(fit, "gam")) return(predict(fit, newdata = newdata, type = "response"))
  predict(fit, newdata = newdata)
}

fit_train <- fit_on(best_model_name, train_ts)
pred_te   <- predict_on(fit_train, test_ts, best_model_name)
test_metrics <- tibble(
  Model = best_model_name,
  RMSE  = rmse(test_ts$R3haddock, pred_te),
  MAE   = mae (test_ts$R3haddock, pred_te),
  R2    = r2  (test_ts$R3haddock, pred_te)
)
print(test_metrics)

# A tibble: 1 x 4
  Model    RMSE     MAE    R2
  <chr>   <dbl>   <dbl> <dbl>
1 GLM   182048. 141692. 0.363

# 4) Диагностика моделей (подгонка на всех данных до 2021) -------------------
full_fit_df <- md_for_fit

lm_full  <- lm(f_lm,  data = full_fit_df)
glm_full <- glm(f_lm, data = full_fit_df, family = Gamma(link = "log"))
gam_full <- mgcv::gam(f_gam, data = full_fit_df, family = Gamma(link = "log"), method = "REML", select = TRUE)

cat("\n[LM] Сводка:\n"); print(summary(lm_full))

[LM] Сводка:

Call:
lm(formula = f_lm, data = full_fit_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-257877 -155326  -18935  101135  326940 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.189e+05  4.459e+05  -2.061 0.049455 *  
codTSB      -2.406e-01  6.722e-02  -3.579 0.001386 ** 
T12          3.679e+05  8.296e+04   4.435 0.000149 ***
I5          -1.770e+03  3.025e+03  -0.585 0.563536    
NAOspring   -5.125e+04  5.427e+04  -0.944 0.353710    
haddock68    4.385e-01  5.395e-01   0.813 0.423698    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 192900 on 26 degrees of freedom
Multiple R-squared:  0.547, Adjusted R-squared:  0.4599 
F-statistic: 6.279 on 5 and 26 DF,  p-value: 0.0006042

cat("\n[LM] VIF:\n"); print(car::vif(lm_full))

[LM] VIF:

   codTSB       T12        I5 NAOspring haddock68 
 2.487391  1.398254  1.275233  1.035349  2.386554 

cat("\n[LM] Breusch–Pagan:\n"); print(lmtest::bptest(lm_full))

[LM] Breusch–Pagan:

    studentized Breusch-Pagan test

data:  lm_full
BP = 5.1481, df = 5, p-value = 0.3981

cat("\n[LM] Durbin–Watson:\n"); print(lmtest::dwtest(lm_full))

[LM] Durbin–Watson:

    Durbin-Watson test

data:  lm_full
DW = 1.7745, p-value = 0.1264
alternative hypothesis: true autocorrelation is greater than 0

glm_resid <- residuals(glm_full, type = "pearson")
cat("\n[GLM-Gamma] Сводка:\n"); print(summary(glm_full))

[GLM-Gamma] Сводка:

Call:
glm(formula = f_lm, family = Gamma(link = "log"), data = full_fit_df)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.907e+00  1.288e+00   5.361 1.30e-05 ***
codTSB      -6.525e-07  1.942e-07  -3.359  0.00242 ** 
T12          1.341e+00  2.397e-01   5.593 7.08e-06 ***
I5           8.270e-03  8.740e-03   0.946  0.35278    
NAOspring    4.898e-02  1.568e-01   0.312  0.75725    
haddock68    1.117e-06  1.559e-06   0.717  0.48000    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.3107759)

    Null deviance: 19.8880  on 31  degrees of freedom
Residual deviance:  8.5197  on 26  degrees of freedom
AIC: 855.29

Number of Fisher Scoring iterations: 8

cat(sprintf("[GLM-Gamma] Pearson dispersion: %.3f\n", sum(glm_resid^2, na.rm = TRUE) / glm_full$df.residual))

[GLM-Gamma] Pearson dispersion: 0.311

cat("\n[GAM] Сводка:\n"); print(summary(gam_full))

[GAM] Сводка:

Family: Gamma 
Link function: log 

Formula:
R3haddock ~ s(codTSB, bs = "tp", k = 5) + s(T12, bs = "tp", k = 5) + 
    s(I5, bs = "tp", k = 5) + s(NAOspring, bs = "tp", k = 5) + 
    s(haddock68, bs = "tp", k = 5)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.48948    0.08886   140.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                   edf Ref.df     F  p-value    
s(codTSB)    1.7083407      4 4.917 7.43e-05 ***
s(T12)       0.9669993      4 7.752 3.96e-06 ***
s(I5)        0.0001638      4 0.000    0.616    
s(NAOspring) 0.0001092      4 0.000    0.980    
s(haddock68) 0.4529453      4 0.145    0.252    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.592   Deviance explained = 60.2%
-REML = 426.83  Scale est. = 0.2527    n = 32

cat("\n[GAM] Concurvity (коротко):\n")

[GAM] Concurvity (коротко):

ccv <- try(mgcv::concurvity(gam_full, full = FALSE), silent = TRUE)
if (!inherits(ccv, "try-error") && is.list(ccv)) {
  # Выведем усечённо и безопасно
  print(lapply(ccv, function(m) if (is.null(m)) NULL else round(m, 3)))
} else {
  cat("не удалось оценить concurvity\n")
}

$worst
             para s(codTSB) s(T12) s(I5) s(NAOspring) s(haddock68)
para            1     0.000  0.000 0.000        0.000        0.000
s(codTSB)       0     1.000  0.554 0.298        0.178        0.821
s(T12)          0     0.554  1.000 0.361        0.362        0.347
s(I5)           0     0.298  0.361 1.000        0.182        0.132
s(NAOspring)    0     0.178  0.362 0.182        1.000        0.216
s(haddock68)    0     0.821  0.347 0.132        0.216        1.000

$observed
             para s(codTSB) s(T12) s(I5) s(NAOspring) s(haddock68)
para            1     0.000  0.000 0.000        0.000        0.000
s(codTSB)       0     1.000  0.392 0.166        0.045        0.385
s(T12)          0     0.273  1.000 0.217        0.064        0.218
s(I5)           0     0.166  0.283 1.000        0.047        0.044
s(NAOspring)    0     0.031  0.128 0.114        1.000        0.049
s(haddock68)    0     0.441  0.204 0.113        0.116        1.000

$estimate
             para s(codTSB) s(T12) s(I5) s(NAOspring) s(haddock68)
para            1     0.000  0.000 0.000        0.000        0.000
s(codTSB)       0     1.000  0.320 0.140        0.041        0.747
s(T12)          0     0.281  1.000 0.201        0.068        0.243
s(I5)           0     0.121  0.278 1.000        0.053        0.096
s(NAOspring)    0     0.055  0.128 0.113        1.000        0.089
s(haddock68)    0     0.544  0.205 0.099        0.136        1.000

invisible(try(mgcv::gam.check(gam_full), silent = TRUE))

Method: REML   Optimizer: outer newton
full convergence after 13 iterations.
Gradient range [-4.544099e-05,0.000320414]
(score 426.834 & scale 0.2526969).
Hessian positive definite, eigenvalue range [5.946259e-06,16.95877].
Model rank =  21 / 21 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                   k'      edf k-index p-value
s(codTSB)    4.000000 1.708341    1.12    0.69
s(T12)       4.000000 0.966999    1.29    0.95
s(I5)        4.000000 0.000164    0.86    0.20
s(NAOspring) 4.000000 0.000109    0.99    0.50
s(haddock68) 4.000000 0.452945    1.10    0.73

# 5) Прогноз 2022–2024 и эмпирические интервалы ------------------------------
best_full <- switch(best_model_name,
  LM  = lm_full,
  GLM = glm_full,
  GAM = gam_full
)

# Остатки для PI: из CV выбранной модели, иначе из полного фита
resids <- if (length(resids_cv[[best_model_name]]) > 5) resids_cv[[best_model_name]] else residuals(best_full)

q025 <- as.numeric(quantile(resids, 0.025, na.rm = TRUE))
q250 <- as.numeric(quantile(resids, 0.250, na.rm = TRUE))
q750 <- as.numeric(quantile(resids, 0.750, na.rm = TRUE))
q975 <- as.numeric(quantile(resids, 0.975, na.rm = TRUE))

fc_start <- 2022
pred_cols <- c("codTSB","T12","I5","NAOspring","haddock68")
mu <- md %>% filter(YEAR > 1989 & YEAR < fc_start) %>% summarise(across(all_of(pred_cols), ~mean(.x, na.rm = TRUE))) %>% as.list()

if (!exists("user_future")) user_future <- NULL

build_future <- function(years, mu, user_df = NULL) {
  df <- tibble::tibble(YEAR = years)
  for (v in pred_cols) df[[v]] <- mu[[v]]
  if (!is.null(user_df)) {
    for (i in seq_len(nrow(user_df))) {
      yr <- user_df$YEAR[i]
      if (yr %in% years) {
        idx <- which(df$YEAR == yr)
        for (v in intersect(pred_cols, names(user_df))) {
          val <- user_df[[v]][i]
          if (!is.na(val)) df[[v]][idx] <- val
        }
      }
    }
  }
  df
}

future_years <- fc_start:2024
scenario_future <- build_future(future_years, mu, user_future)

predict_best <- function(fit, newdata, model_name) {
  if (model_name == "GLM") return(predict(fit, newdata = newdata, type = "response"))
  predict(fit, newdata = newdata)
}

pred_future <- predict_best(best_full, scenario_future, best_model_name)

forecast_tbl <- tibble::tibble(
  YEAR      = scenario_future$YEAR,
  Model     = best_model_name,
  pred_mean = as.numeric(pred_future),
  PI50_low  = pred_future + q250, PI50_high = pred_future + q750,
  PI95_low  = pred_future + q025, PI95_high = pred_future + q975
)


knitr::kable(
  forecast_tbl %>% dplyr::mutate(dplyr::across(where(is.numeric), ~round(.x, 2))),
  caption = "Holdout-метрики (округлено до 2 знаков)"
)

Holdout-метрики (округлено до 2 знаков)

YEAR	Model	pred_mean	PI50_low	PI50_high	PI95_low	PI95_high
2022	GLM	268057.6	172383	509783.3	-203196.1	1051570
2023	GLM	268057.6	172383	509783.3	-203196.1	1051570
2024	GLM	268057.6	172383	509783.3	-203196.1	1051570

# 6) Визуализация 1990–2024 ---------------------------------------------------
pred_df <- bind_rows(
  md %>% select(YEAR, all_of(pred_cols)),
  scenario_future
) %>% distinct(YEAR, .keep_all = TRUE) %>% arrange(YEAR)

pred_df$Pred      <- as.numeric(predict_best(best_full, pred_df, best_model_name))
pred_df$PI50_low  <- pred_df$Pred + q250
pred_df$PI50_high <- pred_df$Pred + q750
pred_df$PI95_low  <- pred_df$Pred + q025
pred_df$PI95_high <- pred_df$Pred + q975

hist_df <- md %>% select(YEAR, R3haddock)

ggplot() +
  geom_ribbon(data = pred_df, aes(x = YEAR, ymin = PI95_low, ymax = PI95_high), fill = "grey80", alpha = 0.25) +
  geom_ribbon(data = pred_df, aes(x = YEAR, ymin = PI50_low, ymax = PI50_high), fill = "grey60", alpha = 0.35) +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = Pred), color = "steelblue4", linewidth = 1) +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = Pred), color = "steelblue4", linewidth = 1, linetype = "dashed") +
  geom_point(data = hist_df, aes(x = YEAR, y = R3haddock), color = "black", size = 2, alpha = 0.9) +
  scale_x_continuous(expand = expansion(mult = c(0, 0))) +
  labs(title = paste0("Пополнение R3haddock: факт (1990–2021) и прогноз (2022–2024) — ", best_model_name),
       subtitle = "Прогноз — пунктир, интервалы — эмпирические из остатков",
       x = "Год", y = "R3haddock") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

# AIC-таблица (LM/GLM сопоставимы напрямую; для GAM также показываем ML)
cat("\nAIC (LM): ",  AIC(lm_full),  "\n", sep = "")

AIC (LM): 877.0549

cat("AIC (GLM): ", AIC(glm_full), "\n", sep = "")

AIC (GLM): 855.2949

gam_full_ml <- mgcv::gam(f_gam, data = full_fit_df, family = Gamma(link = "log"), method = "ML", select = TRUE)
cat("AIC (GAM, REML): ", AIC(gam_full),    "\n", sep = "")

AIC (GAM, REML): 850.8686

cat("AIC (GAM, ML):   ", AIC(gam_full_ml), "\n", sep = "")

AIC (GAM, ML):   850.7948

# ============================================================================
# Конец
# ============================================================================

8.5 Полный цикл от факторов до ансамблевого прогноза

Полный цикл анализа от идентификации ключевых факторов до создания надежного ансамблевого прогноза пополнения рыбных запасов представляет собой сложный, но систематизированный процесс, требующий как глубокого понимания биологических процессов, так и владения современными методами анализа данных. Начиная с формирования исходного набора предикторов, включающего как биологические переменные (нерестовый запас, биомасса хищников), так и комплексные океанографические показатели (температура, соленость, климатические индексы), мы проходим через строгую последовательность этапов, каждый из которых важен для конечного результата. На этапе подготовки данных мы не просто приводим информацию к числовому формату и заменяем строковые обозначения пропущенных значений «NA» на стандартные NA, но и проводим глубокий анализ корреляционной структуры, устраняя мультиколлинеарность через анализ корреляций и VIF-диагностику, что важно для корректной интерпретации последующих моделей. Для обработки пропусков мы применяем медианную импутацию, которая представляет собой простой и устойчивый к выбросам метод, хотя в некоторых случаях могут быть использованы и более сложные методы, такие как KNN-импутация или множественная импутация с использованием пакета MICE, особенно когда данные имеют сложную структуру или временные зависимости.

Затем следует этап отбора предикторов, где мы применяем два комплементарных метода: Boruta на основе Random Forest для выявления нелинейных зависимостей и LASSO-регрессию для линейного отбора с регуляризацией. Их объединение позволяет получить устойчивый набор предикторов, дополненный биологически значимыми переменными по экспертной оценке, что создает баланс между статистической значимостью и содержательной интерпретируемостью. Этот этап является мостом между классической ихтиологией и современными методами анализа, где экспертные знания биолога взаимодействуют с алгоритмической строгостью статистики, гарантируя включение ключевых факторов, таких как нерестовый запас, который должен присутствовать в модели по самой своей природе процесса пополнения.

После подготовки данных мы переходим к сравнению различных семейств моделей через единую кросс-валидационную процедуру (5-fold CV) с последующим хронологическим тестированием на отложенной выборке. Помимо линейных и обобщенных линейных моделей (LM, GLM), обобщенных аддитивных моделей (GAM), мы тестируем современные алгоритмы машинного обучения: Random Forest для улавливания сложных нелинейных зависимостей и взаимодействий между факторами, будучи при этом устойчивым к шуму и выбросам; XGBoost, с его градиентным бустингом над деревьями решений, часто дающий высочайшую точность прогноза; SVM с радиальным ядром для сложных разделяющих поверхностей; и нейронные сети для автоматического извлечения признаков. Каждая модель оценивается по комплексу метрик: RMSE, MAE, R² и MAPE, что позволяет сравнивать их прогностическую силу на разных участках данных и выявлять модели, которые лучше всего справляются с конкретными аспектами прогнозирования.

Особое внимание уделяется временным характеристикам данных, поскольку при анализе водных биоресурсов мы имеем дело с временными рядами, где случайное перемешивание данных приведет к утечке информации из будущего в прошлое, искусственно завысив качество прогноза. Для решения этой проблемы мы применяем специализированную time-slice кросс-валидацию с расширяющимся окном и горизонтом прогноза 3 года, которая имитирует реальные условия прогнозирования, обучаясь только на данных из прошлого и проверяя на последующих периодах. Это позволяет оценить устойчивость моделей к временным сдвигам и их способность к экстраполяции, что критически важно для практических задач управления рыбными запасами.

Выбор окончательной модели — это не просто вопрос максимальной точности на кросс-валидации, а сложный компромисс между точностью, интерпретируемостью и биологической правдоподобностью. Кульминацией цикла становится построение ансамблевой модели, комбинирующей сильные стороны отдельных алгоритмов. В нашем анализе оптимальный ансамбль (CUBIST + LM) строится через взвешенное усреднение предсказаний, где веса определяются на основе кросс-валидационной ошибки — например, 75% веса приходится на мощную нелинейную модель Cubist, а 25% — на простую и устойчивую линейную регрессию. Такой подход позволяет нивелировать индивидуальные недостатки моделей, сохранить интерпретируемость линейных моделей, где биолог может понять, как именно каждый фактор влияет на прогноз, и при этом использовать гибкость методов машинного обучения для захвата сложных нелинейных паттернов, которые могут ускользнуть от классических статистических методов.

Важнейшим компонентом становится оценка неопределенности через эмпирические доверительные интервалы, построенные на основе распределения остатков ансамблевой модели. Мы используем квантили остатков из кросс-валидации для построения 50% и 95% доверительных интервалов, что позволяет получить не только точечный прогноз, но и меру его надежности, важную для принятия управленческих решений. Это дает возможность визуализировать не только ожидаемое значение пополнения, но и диапазон возможных сценариев, что особенно важно в условиях высокой экологической неопределенности.

Финальная визуализация представляет собой совмещение исторических данных с прогнозом на 3 года вперед, где исторические данные отображаются сплошной линией, прогноз — пунктиром, а 50% и 95% доверительные интервалы — серыми лентами различной интенсивности. Такой график не только демонстрирует результат, но и позволяет визуально оценить точность модели на исторических данных и неопределенность будущих предсказаний, делая результаты доступными не только для статистиков, но и для управленцев и политиков, принимающих решения на основе этих прогнозов.

Представленный цикл является итеративным процессом: прогнозная точность ансамбля может быть улучшена через включение новых предикторов, изучение влияния предикторов с задержкой (лагами), тонкую настройку гиперпараметров моделей и обновление данных по мере их поступления. Этот подход представляет собой практический компромисс между статистической строгостью, вычислительной эффективностью и биологической интерпретируемостью, делая его мощным инструментом для решения прикладных задач оценки водных биоресурсов, где каждый этап, от первичной обработки данных до финального прогноза, подчинен одной цели — обеспечению устойчивого управления рыбными запасами на основе надежного научного анализа.

Скрипт лучше скачать целиком).

# ==============================================================================
# ПРАКТИЧЕСКОЕ ЗАНЯТИЕ: АНАЛИЗ ФАКТОРОВ И ПРОГНОЗ ПОПОЛНЕНИЯ ЗАПАСА
# Курс: "Оценка водных биоресурсов в среде R (для начинающих)"
# Автор: Баканев С. В. Дата:20.08.2025
# Структура:
# 1) Подготовка данных и выбор предикторов
# 2) Базовое сравнение моделей (5-fold CV + holdout)
# 3) Выбор лучшей прогностической модели (time-slice CV на 3 года + хронологический тест)
# 4) Прогноз 2022–2024 (ансамбль CUBIST+LM) и график 1990–2024 с ДИ
# ------------------------------------------------------------------------------
# Пояснения к занятию (для начинающих):
# - Мы работаем с временным рядом пополнения запаса R3haddock и набором факторов
#   среды/биомассы. Цель — построить понятные и проверяемые модели прогноза.
# - Сначала отберём информативные предикторы (Boruta и LASSO), затем сравним
#   разные модели машинного обучения на кросс-валидации (CV), после чего выберем
#   лучшую схему по time-slice CV (учитывая хронологию), и сделаем прогноз.
# ==============================================================================


# ==============================================================================
# 1) ВЫБОР ПРЕДИКТОРОВ
# ------------------------------------------------------------------------------
# Цель блока: привести данные к числовому виду, обработать пропуски, сократить
# мультиколлинеарность (сильные корреляции), а затем автоматически выделить
# кандидатов-предикторов двумя методами (Boruta, LASSO). В конце сформируем
# финальный пул признаков и проверим их значимость в простой LM.
# ==============================================================================

# Установка и подключение необходимых библиотек
# Для автоматического отбора предикторов нам понадобятся дополнительные пакеты
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl, tidyverse, caret, corrplot, mgcv, randomForest, xgboost,
  Boruta,GGally, FactoMineR, glmnet, recipes, rsample  # Новые библиотеки для автоматического отбора
)

# Очистка среды и установка рабочей директории
# Совет: rm(list=ls()) очищает все объекты в памяти R; setwd задаёт папку,
# где искать/сохранять файлы. Убедитесь, что путь корректен на вашей машине.
rm(list = ls())
setwd("C:/RECRUITMENT/")

# Пакеты для расширенного отбора предикторов
# Boruta — обёртка над Random Forest для отбора признаков;
# glmnet — регуляризация (LASSO/ElasticNet) для отбора/усиления обобщающей способности;
# FactoMineR — PCA и другие многомерные методы (используем как утилиту).
library(Boruta)   # Алгоритм обертки для отбора признаков
library(glmnet)   # LASSO-регрессия
library(FactoMineR) # PCA анализ


# Загрузка и первичная обработка данных
# Шаги: фильтруем годы, приводим типы к числовому, заменяем строковые "NA" на NA.
DATA <- readxl::read_excel("RECRUITMENT.xlsx", sheet = "RECRUITMENT") %>%
  filter(YEAR > 1989 & YEAR < 2022) %>%
  # Преобразуем необходимые столбцы в числовой формат
  mutate(
    across(starts_with("T"), as.numeric),
    across(starts_with("I"), as.numeric),
    across(starts_with("O"), as.numeric),
  ) %>%
  # Обработка пропущенных значений (заменяем строку "NA" на NA)
  mutate(across(where(is.character), ~na_if(., "NA")))

# 1. Подготовка данных -------------------------------------------------------
# Выделим все возможные предикторы, включая географию и индексы трески
# Примечание: оставляем только числовые переменные, т.к. большинство моделей
# требует числовой вход без категориальных уровней.
predictors <- DATA %>% 
  select(-YEAR, -R3haddock) %>% 
  select_if(is.numeric) # Только числовые переменные

# Целевая переменная
response <- DATA$R3haddock

# В статистическом анализе мы различаем:
# - Отклик (response/target variable) - то, что мы пытаемся предсказать (в нашем случае R3haddock)
# - Предикторы (predictors/features) - переменные, которые могут объяснять изменения отклика
# Для корректного анализа важно, чтобы предикторы были числовыми или преобразованы в числовой формат.

# 2. Обработка пропусков -----------------------------------------------------
# Заполнение медианными значениями — простой и устойчивый способ справиться с NA.
# Альтернативы: множественная иммутация (mice), KNN-impute и др.
predictors_filled <- predictors %>%
  mutate(across(everything(), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))

# Заполнение медианой - простой и устойчивый метод обработки пропусков для числовых переменных.
# Медиана предпочтительнее среднего, так как менее чувствительна к выбросам.

# 3. Предварительный анализ корреляций ---------------------------------------
# Зачем: высокие корреляции затрудняют интерпретацию и могут вредить ряду моделей.
cor_matrix <- cor(predictors_filled, use = "complete.obs")
corrplot(cor_matrix, method = "circle", type = "upper", tl.cex = 0.7)

# Удаляем высокоскоррелированные предикторы (r > 0.8)
# Это механическое сокращение мультиколлинеарности до этапа отбора.
high_cor <- findCorrelation(cor_matrix, cutoff = 0.8)
predictors_filtered <- predictors_filled[, -high_cor]

# Высокая корреляция между предикторами (мультиколлинеарность) может привести к нестабильности моделей.
# Например, если два предиктора почти идентичны, модель может неустойчиво распределять их влияние на отклик.
# Удаление сильно коррелированных переменных (r > 0.8) помогает улучшить интерпретируемость и стабильность моделей.


# 4. Автоматизированный отбор Boruta (обертка Random Forest) -----------------
# Идея: определить признаки, которые важнее, чем случайный шум (shadow features).

# Визуализация результатов
plot(boruta_output, cex.axis = 0.7, las = 2)

boruta_stats <- attStats(boruta_output)
selected_vars <- getSelectedAttributes(boruta_output, withTentative = TRUE)

# Boruta - это алгоритм отбора признаков, основанный на методе случайного леса.
# Он сравнивает важность реальных переменных с "теневыми" переменными (случайными копиями),
# чтобы определить, действительно ли переменная информативна. 
# Результаты Boruta показывают: 
#   - Confirmed (зеленые) - значимые предикторы
#   - Tentative (желтые) - предикторы, близкие к порогу значимости
#   - Rejected (красные) - незначимые предикторы


# 5. LASSO с более строгим критерием ------------------------------------------
# Идея: L1-регуляризация зануляет коэффициенты «слабых» предикторов.
# Выбор lambda.1se вместо lambda.min — более консервативный (простая модель).
x <- as.matrix(predictors_filtered)
y <- response

# LASSO (Least Absolute Shrinkage and Selection Operator) - метод регрессии с L1-регуляризацией,
# который одновременно выполняет отбор признаков и оценку коэффициентов. 
# Параметр lambda контролирует силу регуляризации:
#   - lambda.min дает наименьшую ошибку, но может включать шумовые переменные
#   - lambda.1se (на 1 стандартную ошибку больше) дает более простую модель с меньшим риском переобучения
# Для прогнозирования мы предпочитаем более строгий критерий (lambda.1se), чтобы модель была устойчивее. 

# Кросс-валидация
cv_fit <- cv.glmnet(x, y, alpha = 1, nfolds = 10)
plot(cv_fit)

# ИСПОЛЬЗУЕМ lambda.1se вместо lambda.min — СТРОЖЕ!
lasso_coef <- coef(cv_fit, s = "lambda.1se")  # <-- Ключевое изменение!
lasso_vars <- rownames(lasso_coef)[lasso_coef[,1] != 0][-1]  # исключаем (Intercept)


# 6. Сравнение отобранных предикторов ----------------------------------------
# Полезно видеть, какие признаки отмечают оба метода (устойчивые кандидаты).
cat("Boruta selected:", length(selected_vars), "variables\n")

Boruta selected: 3 variables

print(selected_vars)

[1] "codTSB" "T12"    "I5"    

cat("\nLASSO selected:", length(lasso_vars), "variables\n")

LASSO selected: 5 variables

print(lasso_vars)

[1] "codTSB" "T12"    "NAO3"   "NAO4"   "NAO5"  

# 7. Финальный набор предикторов (объединение результатов) -------------------
# Логика: объединяем списки, добавляем биологически важные переменные вручную.
final_vars <- union(selected_vars, lasso_vars) 

# Добавляем обязательные переменные по биологической логике
mandatory <- c("haddock68")
final_vars <- union(final_vars, mandatory) %>% unique()

# Мы объединяем результаты двух методов отбора признаков для большей надежности.
# Также добавляем переменную haddock68 (нерестовый запас), так как биологически 
# логично, что пополнение запаса напрямую зависит от численности производителей. 
# Это пример интеграции экспертных знаний в статистический анализ - важный принцип 
# при работе с данными в биологических науках.

# 8. Проверка значимости -----------------------------------------------------
# Быстрая оценка значимости с LM: не как окончательный вывод, а как sanity-check.
final_model <- lm(response ~ as.matrix(predictors_filled[, final_vars]))
summary(final_model)

Call:
lm(formula = response ~ as.matrix(predictors_filled[, final_vars]))

Residuals:
    Min      1Q  Median      3Q     Max 
-270986  -82376   -1037   98086  276129 

Coefficients:
                                                      Estimate Std. Error
(Intercept)                                         -1.082e+06  3.943e+05
as.matrix(predictors_filled[, final_vars])codTSB    -2.346e-01  5.536e-02
as.matrix(predictors_filled[, final_vars])T12        3.864e+05  7.198e+04
as.matrix(predictors_filled[, final_vars])I5        -1.825e+02  2.572e+03
as.matrix(predictors_filled[, final_vars])NAO3      -5.801e+04  3.129e+04
as.matrix(predictors_filled[, final_vars])NAO4       8.345e+04  3.035e+04
as.matrix(predictors_filled[, final_vars])NAO5      -7.278e+04  2.488e+04
as.matrix(predictors_filled[, final_vars])haddock68  1.232e-01  4.515e-01
                                                    t value Pr(>|t|)    
(Intercept)                                          -2.744 0.011305 *  
as.matrix(predictors_filled[, final_vars])codTSB     -4.238 0.000288 ***
as.matrix(predictors_filled[, final_vars])T12         5.368 1.64e-05 ***
as.matrix(predictors_filled[, final_vars])I5         -0.071 0.944028    
as.matrix(predictors_filled[, final_vars])NAO3       -1.854 0.076118 .  
as.matrix(predictors_filled[, final_vars])NAO4        2.750 0.011146 *  
as.matrix(predictors_filled[, final_vars])NAO5       -2.925 0.007412 ** 
as.matrix(predictors_filled[, final_vars])haddock68   0.273 0.787227    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 158600 on 24 degrees of freedom
Multiple R-squared:  0.7173,    Adjusted R-squared:  0.6348 
F-statistic: 8.698 on 7 and 24 DF,  p-value: 2.58e-05

# 9. Формирование финального датасета ----------------------------------------
# Собираем набор с откликом и выбранными предикторами; удалим строки с NA.
model_data <- DATA %>%
  select(R3haddock, all_of(final_vars)) %>%
  drop_na()

# Просмотр структуры финальных данных
glimpse(model_data)

Rows: 32
Columns: 8
$ R3haddock <dbl> 812363, 389416, 99474, 98946, 118812, 63028, 147657, 83270, ~
$ codTSB    <dbl> 913000, 1347064, 1687381, 2197863, 2112773, 1849957, 1697388~
$ T12       <dbl> 4.72, 4.66, 4.24, 3.90, 3.96, 4.27, 4.16, 4.07, 4.23, 5.08, ~
$ I5        <dbl> 43, 55, 26, 49, 56, 28, 52, 51, 69, 68, 41, 48, 50, 63, 40, ~
$ NAO3      <dbl> 1.46, -0.20, 0.87, 0.67, 1.26, 1.25, -0.24, 1.46, 0.87, 0.23~
$ NAO4      <dbl> 2.00, 0.29, 1.86, 0.97, 1.14, -0.85, -0.17, -1.02, -0.68, -0~
$ NAO5      <dbl> -1.53, 0.08, 2.63, -0.78, -0.57, -1.49, -1.06, -0.28, -1.32,~
$ haddock68 <dbl> 74586, 79205, 53195, 36337, 49122, 81514, 172177, 160886, 96~

# Визуализация важности переменных
# Внимание: важности от RF — относительные; сопоставляйте с предметной логикой.
var_importance <- randomForest(R3haddock ~ ., data = model_data, importance = TRUE)
varImpPlot(var_importance, main = "Важность предикторов")

# Перед окончательным выбором модели мы проверяем значимость предикторов с помощью линейной регрессии.
# Функция summary() показывает p-значения коэффициентов - если p < 0.05, переменная считается статистически значимой. 
# Визуализация важности переменных с помощью случайного леса дает дополнительную перспективу,
# показывая, какие переменные наиболее информативны для предсказания без предположений о линейности.

# ==============================================================================
#  ПОДГОТОВКА ДАННЫХ
# Создаём NAOspring, фиксируем финальный набор признаков, сохраняем CSV.
# ------------------------------------------------------------------------------
# Цель блока: стандартизировать набор признаков для дальнейшего сравнения
# моделей и обеспечить воспроизводимость (фиксированный CSV с нужными полями).
# ==============================================================================

# 1.1 Пакеты и окружение
# Примечание: блок повторяет базовую инициализацию для автономного запуска.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readxl, tidyverse, caret, corrplot)

rm(list = ls())
set.seed(123)
setwd("C:/RECRUITMENT/")

# 1.2 Загрузка исходных данных и приведение типов
DATA <- readxl::read_excel("RECRUITMENT.xlsx", sheet = "RECRUITMENT") %>%
  filter(YEAR > 1989 & YEAR < 2022) %>%
  mutate(
    across(starts_with("T"), as.numeric),
    across(starts_with("I"), as.numeric),
    across(starts_with("O"), as.numeric),
    across(where(is.character), ~na_if(., "NA"))
  )

# 1.3 Создаём NAOspring (если есть NAO3, NAO4, NAO5)
# Идея: агрегируем весенний индекс NAO как среднее за месяцы 3–5.
if (all(c("NAO3","NAO4","NAO5") %in% names(DATA))) {
  DATA <- DATA %>%
    mutate(NAOspring = rowMeans(pick(NAO3, NAO4, NAO5), na.rm = TRUE)) %>%
    select(-NAO3, -NAO4, -NAO5)
}

# NAO (North Atlantic Oscillation) - важный климатический индекс, влияющий описывающий изменения атмосферного давления
# над Северной Атлантикой. В частности, он отражает разницу в атмосферном давлении между Исландской депрессией и
# Азорским максимумом. NAO влияет на силу и направление западных ветров, а также на траектории штормов в Северной Атлантике. 
# Мы создаем NAOspring как среднее значение за весенние месяцы (марта, апреля, мая),
# так как именно в этот период происходят ключевые процессы, влияющие на нерест трески. 
# Создание составных переменных на основе экспертных знаний часто улучшает качество моделей.

# 1.4 Финальный учебный набор предикторов (фиксируем)
# Важно: проверяем присутствие нужных колонок и формируем компактный датасет.
needed <- c("codTSB", "T12", "I5", "NAOspring", "haddock68")
stopifnot(all(needed %in% names(DATA)))

# Сохраняем YEAR в CSV (ниже он будет отброшен при обучении, но нужен для графика)
model_data <- DATA %>%
  select(YEAR, all_of(needed), R3haddock) %>%
  drop_na()

write.csv(model_data, "selected_predictors_dataset.csv", row.names = FALSE)
glimpse(model_data)

Rows: 32
Columns: 7
$ YEAR      <dbl> 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, ~
$ codTSB    <dbl> 913000, 1347064, 1687381, 2197863, 2112773, 1849957, 1697388~
$ T12       <dbl> 4.72, 4.66, 4.24, 3.90, 3.96, 4.27, 4.16, 4.07, 4.23, 5.08, ~
$ I5        <dbl> 43, 55, 26, 49, 56, 28, 52, 51, 69, 68, 41, 48, 50, 63, 40, ~
$ NAOspring <dbl> 0.64333333, 0.05666667, 1.78666667, 0.28666667, 0.61000000, ~
$ haddock68 <dbl> 74586, 79205, 53195, 36337, 49122, 81514, 172177, 160886, 96~
$ R3haddock <dbl> 812363, 389416, 99474, 98946, 118812, 63028, 147657, 83270, ~

# (необязательно) Глянуть попарные связи и корреляции
# ggpairs может быть медленным, оставим по желанию
 ggpairs(model_data, columns = 2:7,
         lower = list(continuous = wrap("smooth", alpha = 0.3, size = 0.5)),
         upper = list(cor = wrap("cor", size = 3)))

# ==============================================================================
# 2) БАЗОВОЕ СРАВНЕНИЕ МОДЕЛЕЙ (5-FOLD CV + HOLDOUT)
# Единые фолды CV, тренировочно-тестовое разбиение, сводка метрик.
# ------------------------------------------------------------------------------
# Идея блока: быстрая «панель» сравнения разных семейств моделей на одинаковых
# условиях (одинаковые фолды CV) и внешний тест (holdout). Это помогает увидеть
# уровни ошибок и выбрать несколько лидеров для более строгой проверки далее.
# ==============================================================================

# 2.1 Пакеты и данные
pacman::p_load(mgcv, randomForest, xgboost, nnet, earth, kernlab, pls, Cubist, ranger, gbm, lattice)

model_data <- read.csv("selected_predictors_dataset.csv", header = TRUE, stringsAsFactors = FALSE)
# Если YEAR отсутствует (на всякий случай), создадим
if (!"YEAR" %in% names(model_data)) {
  model_data$YEAR <- seq(1990, by = 1, length.out = nrow(model_data))
}

# Используем только предикторы и отклик (YEAR исключаем)
model_data <- model_data %>%
  select(codTSB, T12, I5, NAOspring, haddock68, R3haddock) %>%
  na.omit()

# 2.2 Holdout и CV-контроллер
# Пропорция 80/20 обеспечивает внешний тест; внутри train — 5-fold CV для
# корректной настройки моделей и оценки средней ошибки.
train_idx <- caret::createDataPartition(model_data$R3haddock, p = 0.8, list = FALSE)
train <- model_data[train_idx, ]
test  <- model_data[-train_idx, ]

ctrl <- caret::trainControl(method = "cv", number = 5, savePredictions = "final")

# Holdout-метод: мы делим данные на обучающую (80%) и тестовую (20%) выборки.
# Кросс-валидация (5-fold CV): данные разбиваются на 5 частей, модель обучается на 4 частях и тестируется на 5-й, 
# и этот процесс повторяется 5 раз. Это дает более надежную оценку качества модели, чем одно разбиение. 


# 2.3 Кастомный GAM (mgcv) для caret (bs="tp", REML, select=TRUE)
# GAM даёт гладкие нелинейности по каждому признаку; REML стабилизирует оценку.
gam_spec <- list(
  type = "Regression", library = "mgcv", loop = NULL,
  parameters = data.frame(parameter = "none", class = "character", label = "none"),
  grid = function(x,y,len=NULL,search="grid") data.frame(none = NA),
  fit = function(x,y,...) {
    df <- x; df$R3haddock <- y
    mgcv::gam(
      R3haddock ~ s(codTSB,bs="tp") + s(T12,bs="tp") + s(I5,bs="tp") +
                  s(NAOspring,bs="tp") + s(haddock68,bs="tp"),
      data=df, method="REML", select=TRUE, ...
    )
  },
  predict = function(modelFit, newdata, submodels = NULL) {
    predict(modelFit, newdata = newdata, type = "response")
  },
  prob = NULL, sort = function(x) x
)

# 2.4 Обучение моделей
# Подсказка: разные методы по-разному чувствительны к масштабу, числу признаков
# и мультиколлинеарности. Мы применяем одинаковые фолды CV для честного сравнения.

# --- 1. Линейная регрессия (LM)
# Учебный смысл: базовая линейная модель; ориентир для сравнения.
# ПОЯСНЕНИЕ: LM предполагает линейную зависимость между предикторами и откликом.
# Это простая модель, которая служит "нижней планкой" - более сложные модели могут быть лучше LM. 
lm_model    <- caret::train(R3haddock ~ ., data = train, method = "lm", trControl = ctrl)

# --- 2. Обобщённая линейная модель (GLM: Gamma с лог-ссылкой)
# Учебный смысл: модель для положительных откликов; допускает нелинейность в шкале log.
# ПОЯСНЕНИЕ: GLM с Gamma-распределением подходит для положительных непрерывных данных 
# (как размер популяции), где дисперсия зависит от среднего значения.
glm_model   <- caret::train(R3haddock ~ ., data = train, method = "glm",
                            family = Gamma(link = "log"), trControl = ctrl)

# --- 3. Обобщённая аддитивная модель (GAM, mgcv: bs="tp", REML, select=TRUE)
# Учебный смысл: гибкие гладкие нелинейности по каждому предиктору.
# ПОЯСНЕНИЕ: GAM позволяет моделировать нелинейные зависимости с помощью гладких функций (splines),
# сохраняя интерпретируемость отдельных эффектов. Это компромисс между простотой LM и сложностью ML.
gam_model   <- caret::train(x = train[, -which(names(train)=="R3haddock")],
                            y = train$R3haddock, method = gam_spec, trControl = ctrl)

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

# --- 4. Random Forest (rf: ntree=1000, mtry=1)
# Учебный смысл: ансамбль деревьев; устойчив к шуму; нелинейности/взаимодействия "из коробки".
# ПОЯСНЕНИЕ: Random Forest строит множество деревьев решений и усредняет их результаты.
# Это мощный метод, который автоматически улавливает нелинейные зависимости и взаимодействия. 
rf_model    <- caret::train(R3haddock ~ ., data = train, method = "rf", trControl = ctrl,
                            ntree = 1000, tuneGrid = data.frame(mtry = 1), importance = TRUE)

# --- 5. XGBoost (xgbTree) 
# Учебный смысл: бустинг деревьев; сильная ML-модель, легко переобучается без валидации.
# ПОЯСНЕНИЕ: XGBoost - это градиентный бустинг над деревьями решений, который последовательно 
# строит деревья, исправляя ошибки предыдущих. Требует тщательной настройки параметров.
xgb_grid    <- expand.grid(nrounds=100, max_depth=4, eta=0.1, gamma=0,
                           colsample_bytree=0.8, min_child_weight=1, subsample=0.8)
xgb_model   <- caret::train(R3haddock ~ ., data = train, method = "xgbTree",
                            trControl = ctrl, tuneGrid = xgb_grid, verbose = 0)

# --- 6. Нейросеть (MLP, nnet: линейный выход, стандартизация)
# Учебный смысл: универсальный аппроксиматор; чувствителен к масштабу; требует регуляризации.
# ПОЯСНЕНИЕ: Нейронные сети могут моделировать сложные нелинейные отношения. 
# Используемая архитектура (1 скрытый слой) - компромисс между гибкостью и риском переобучения.
# Линейный выходной слой подходит для регрессии.
nnet_model  <- caret::train(R3haddock ~ ., data = train, method = "nnet",
                            trControl = ctrl, preProcess = c("center","scale"),
                            tuneGrid = expand.grid(size = 5, decay = 0.1),
                            linout = TRUE, trace = FALSE, MaxNWts = 5000)

# --- 7. Elastic Net (glmnet)
# Учебный смысл: регуляризация (L1/L2), борьба с мультиколлинеарностью, частичный отбор признаков.
# ПОЯСНЕНИЕ: Комбинирует L1 (лассо) и L2 (ридж) регуляризации. Автоматически отбирает признаки 
# и уменьшает влияние мультиколлинеарности. Параметр alpha балансирует между лассо и риджем.
glmnet_model<- caret::train(R3haddock ~ ., data = train, method = "glmnet",
                            trControl = ctrl, preProcess = c("center","scale"), tuneLength = 10)

# --- 8. MARS (earth)
# Учебный смысл: кусочно-линейные сплайны + простые взаимодействия; гибкая интерпретация.
# ПОЯСНЕНИЕ: Многомерные адаптивные регрессионные сплайны (MARS) строят кусочно-линейные модели 
# с автоматическим выбором точек излома. Поддерживает взаимодействия ограниченного порядка.
earth_model <- caret::train(R3haddock ~ ., data = train, method = "earth",
                            trControl = ctrl, tuneLength = 10)

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

# --- 9. SVM с радиальным ядром (svmRadial)
# Учебный смысл: ядровой метод; улавливает сложные нелинейности; важна стандартизация.
# ПОЯСНЕНИЕ: Метод опорных векторов с радиальным ядром проецирует данные в пространство 
# высокой размерности, где становится возможным линейное разделение. Параметр gamma управляет 
# гибкостью границы решения.
svm_model   <- caret::train(R3haddock ~ ., data = train, method = "svmRadial",
                            trControl = ctrl, preProcess = c("center","scale"), tuneLength = 8)

# --- 10. k-ближайших соседей (kNN)
# Учебный смысл: простая интуитивная нелинейная модель на расстояниях; чувствительна к масштабу.
# ПОЯСНЕНИЕ: Предсказание основано на усреднении значений k ближайших наблюдений. 
# Требует вычисления попарных расстояний, что может быть ресурсоемким при больших данных.
knn_model   <- caret::train(R3haddock ~ ., data = train, method = "knn",
                            trControl = ctrl, preProcess = c("center","scale"), tuneLength = 15)

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 23 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 25 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 27 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 29 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 31 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.54402860027016,
-1.04267060653844, : k = 33 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 23 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 25 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 27 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 29 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 31 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.32034464856588,
-0.795974449614684, : k = 33 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.59980040948893,
-0.134264066935858, : k = 25 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.59980040948893,
-0.134264066935858, : k = 27 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.59980040948893,
-0.134264066935858, : k = 29 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.59980040948893,
-0.134264066935858, : k = 31 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.59980040948893,
-0.134264066935858, : k = 33 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 23 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 25 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 27 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 29 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 31 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.47821802102901,
-0.982937987281781, : k = 33 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.07995722209223,
-0.0328010741655768, : k = 25 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.07995722209223,
-0.0328010741655768, : k = 27 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.07995722209223,
-0.0328010741655768, : k = 29 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.07995722209223,
-0.0328010741655768, : k = 31 exceeds number 23 of patterns

Warning in knnregTrain(train = structure(c(-1.07995722209223,
-0.0328010741655768, : k = 33 exceeds number 23 of patterns

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

# --- 11. Ranger (быстрый Random Forest)
# Учебный смысл: альтернативная/быстрая реализация леса; сравнить с randomForest.
# ПОЯСНЕНИЕ: Оптимизированная реализация Random Forest на C++. Поддерживает распараллеливание 
# и эффективную работу с категориальными переменными. Важен параметр mtry (число признаков в узле).
ranger_model<- caret::train(R3haddock ~ ., data = train, method = "ranger",
                            trControl = ctrl, tuneLength = 3, importance = "impurity")

# --- 12. GBM (классический градиентный бустинг)
# Учебный смысл: другой бустинг деревьев; полезно сравнить с XGBoost.
# ПОЯСНЕНИЕ: Градиентный бустинг строит деревья последовательно, где каждое новое дерево 
# корректирует ошибки предыдущих. Параметр shrinkage (темп обучения) контролирует скорость обучения.
gbm_model   <- caret::train(R3haddock ~ ., data = train, method = "gbm",
                            trControl = ctrl,
                            tuneGrid = expand.grid(n.trees=100, interaction.depth=1,
                                                   shrinkage=0.1, n.minobsinnode=2),
                            distribution = "gaussian", bag.fraction = 1, verbose = FALSE)

# --- 13. PLS (Partial Least Squares)
# Учебный смысл: проекция на скрытые компоненты с учетом отклика; решает мультиколлинеарность.
# ПОЯСНЕНИЕ: Частные наименьшие квадраты (PLS) проецируют предикторы в латентное пространство, 
# максимизируя ковариацию с откликом. Эффективен при высокой корреляции признаков.
pls_model   <- caret::train(R3haddock ~ ., data = train, method = "pls",
                            trControl = ctrl, preProcess = c("center","scale"), tuneLength = 10)

# --- 14. Cubist (правила + деревья)
# Учебный смысл: интерпретируемые правила с комитетами; часто силен на табличных данных.
# ПОЯСНЕНИЕ: Cubist объединяет деревья решений с линейными моделями в листьях. Генерирует 
# набор правил "если-то", что улучшает интерпретируемость. Комитеты (комитеты) уменьшают дисперсию.
cubist_model<- caret::train(R3haddock ~ ., data = train, method = "cubist",
                            trControl = ctrl, tuneLength = 5)

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

# 2.5 Метрики и оценка на тесте
# Замечание: RMSE/MAE — абсолютные ошибки; R2 — доля объяснённой вариации;
# MAPE/sMAPE — относительные ошибки (осторожно при малых значениях отклика).
rmse  <- function(a, p) sqrt(mean((a - p)^2, na.rm = TRUE))
mae   <- function(a, p) mean(abs(a - p), na.rm = TRUE)
r2    <- function(a, p) 1 - sum((a - p)^2, na.rm = TRUE) / sum((a - mean(a))^2, na.rm = TRUE)
mape  <- function(a, p) mean(abs((a - p) / a), na.rm = TRUE) * 100
smape <- function(a, p) mean(2 * abs(p - a) / (abs(a) + abs(p)), na.rm = TRUE) * 100
metrics_vec <- function(y, pred) c(RMSE=rmse(y,pred), MAE=mae(y,pred), R2=r2(y,pred),
                                   MAPE=mape(y,pred), sMAPE=smape(y,pred))

# Для оценки качества моделей мы используем несколько метрик:
#   - RMSE (Root Mean Square Error): среднеквадратичная ошибка (чувствительна к выбросам)
#   - MAE (Mean Absolute Error): средняя абсолютная ошибка (более интерпретируема)
#   - R²: коэффициент детерминации (доля объясненной дисперсии)
#   - MAPE: средняя абсолютная процентная ошибка (в процентах от фактического значения)
#   - sMAPE: симметричная MAPE (устраняет проблему деления на ноль) 

y_test <- test$R3haddock
preds_test <- list(
  LM=predict(lm_model,test), GLM=predict(glm_model,test), GAM=predict(gam_model,test),
  RF=predict(rf_model,test), XGB=predict(xgb_model,test), NNET=predict(nnet_model,test),
  ENet=predict(glmnet_model,test), MARS=predict(earth_model,test), SVM=predict(svm_model,test),
  kNN=predict(knn_model,test), RANGER=predict(ranger_model,test), GBM=predict(gbm_model,test),
  PLS=predict(pls_model,test), CUBIST=predict(cubist_model,test)
)
metrics_table <- do.call(rbind, lapply(names(preds_test), function(nm){
  data.frame(Model = nm, t(metrics_vec(y_test, preds_test[[nm]])), row.names = NULL)
})) %>% arrange(RMSE, MAE)

# Создаем копию таблицы для округления
metrics_table_rounded <- metrics_table

# Находим индексы числовых столбцов (исключая первый столбец "Model")
numeric_cols <- sapply(metrics_table_rounded, is.numeric)

# Округляем только числовые столбцы до 2 знаков
metrics_table_rounded[numeric_cols] <- round(metrics_table_rounded[numeric_cols], 2)

# Выводим округленную таблицу
print(metrics_table_rounded)

    Model      RMSE       MAE    R2   MAPE sMAPE
1     GBM  65112.49  55805.52  0.76  34.41 27.07
2    MARS  95150.47  78452.71  0.49  51.69 36.18
3      RF 125056.40  97862.47  0.12  74.27 44.45
4     PLS 130411.86 111694.80  0.04  62.51 47.69
5      LM 131592.23 113063.92  0.02  63.62 48.29
6     GAM 140393.68 121236.06 -0.11  70.56 48.99
7  CUBIST 147031.92 123320.33 -0.22  50.52 46.56
8    ENet 147470.05 124147.91 -0.23  81.27 52.96
9  RANGER 148623.73 127762.15 -0.25  87.21 54.25
10    kNN 149082.19 128657.28 -0.26  88.53 53.55
11    GLM 151391.83 129510.77 -0.30  74.40 55.01
12    SVM 167864.37 147269.48 -0.59 104.73 57.83
13    XGB 208719.19 194346.48 -1.46  88.97 65.24
14   NNET 216172.26 199594.44 -1.64  93.12 87.93

knitr::kable(
  metrics_table %>% dplyr::mutate(dplyr::across(where(is.numeric), ~round(.x, 2))),
  caption = "Holdout-метрики (округлено до 2 знаков)"
)

Holdout-метрики (округлено до 2 знаков)

Model	RMSE	MAE	R2	MAPE	sMAPE
GBM	65112.49	55805.52	0.76	34.41	27.07
MARS	95150.47	78452.71	0.49	51.69	36.18
RF	125056.40	97862.47	0.12	74.27	44.45
PLS	130411.86	111694.80	0.04	62.51	47.69
LM	131592.23	113063.92	0.02	63.62	48.29
GAM	140393.68	121236.06	-0.11	70.56	48.99
CUBIST	147031.92	123320.33	-0.22	50.52	46.56
ENet	147470.05	124147.91	-0.23	81.27	52.96
RANGER	148623.73	127762.15	-0.25	87.21	54.25
kNN	149082.19	128657.28	-0.26	88.53	53.55
GLM	151391.83	129510.77	-0.30	74.40	55.01
SVM	167864.37	147269.48	-0.59	104.73	57.83
XGB	208719.19	194346.48	-1.46	88.97	65.24
NNET	216172.26	199594.44	-1.64	93.12	87.93

# 2.6 CV-резюме
# Сводим результаты CV по всем моделям и смотрим распределения ошибок.
results <- caret::resamples(list(
  LM=lm_model, GLM=glm_model, GAM=gam_model, RF=rf_model, XGB=xgb_model, NNET=nnet_model,
  ENet=glmnet_model, MARS=earth_model, SVM=svm_model, kNN=knn_model, RANGER=ranger_model,
  GBM=gbm_model, PLS=pls_model, CUBIST=cubist_model
))
summary(results)

Call:
summary.resamples(object = results)

Models: LM, GLM, GAM, RF, XGB, NNET, ENet, MARS, SVM, kNN, RANGER, GBM, PLS, CUBIST 
Number of resamples: 5 

MAE 
            Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
LM      57204.50 173765.7 178162.9 167545.5 210890.4 217704.1    0
GLM    121162.89 123189.9 126872.5 161564.3 215648.4 220947.6    0
GAM    113581.17 186001.4 194396.8 216924.0 243663.3 346977.3    0
RF     140917.34 173192.3 183106.9 211428.5 199504.4 360421.6    0
XGB    191739.57 196545.5 204224.0 220586.8 211582.9 298842.3    0
NNET   210800.50 230696.3 231946.6 255228.9 268869.5 333831.8    0
ENet    95118.26 114161.5 186843.2 175901.3 203787.5 279596.0    0
MARS   100715.25 160470.7 226674.7 224830.7 281753.1 354539.8    0
SVM    137655.75 168734.2 186772.6 218295.0 270224.5 328087.8    0
kNN    153725.61 173448.8 204644.5 201592.1 210974.4 265167.1    0
RANGER 134039.92 173499.6 173657.3 188631.8 213348.0 248614.3    0
GBM    142143.64 174377.4 182594.4 195116.5 210961.4 265505.7    0
PLS    140180.64 169826.6 174374.4 174539.4 177221.3 211093.9    0
CUBIST  40943.57 166731.9 172415.9 151410.2 183157.2 193802.4    0

RMSE 
            Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
LM      83394.05 204460.2 213949.1 198791.4 234397.7 257756.2    0
GLM    144920.86 174777.3 199683.5 211248.8 268416.8 268445.7    0
GAM    126481.10 212799.6 231166.1 251450.9 265213.5 421594.3    0
RF     152505.15 221538.9 237863.1 264290.3 246285.2 463259.3    0
XGB    224482.10 230814.5 274803.3 282033.4 324393.6 355673.5    0
NNET   262262.48 292982.1 298735.5 322398.1 312361.7 445649.0    0
ENet    96405.72 194703.2 212260.2 215307.8 227408.1 345761.6    0
MARS   133708.92 174187.4 275486.7 264992.8 296803.1 444778.0    0
SVM    175110.30 208575.8 211026.8 261495.0 330657.7 382104.5    0
kNN    191674.11 194743.2 293783.4 261451.2 307657.4 319397.8    0
RANGER 179305.47 214964.7 252239.7 245471.8 256030.3 324818.7    0
GBM    155235.18 236112.1 265637.9 269271.4 342304.1 347067.5    0
PLS    192456.23 204052.4 206556.2 210244.7 207702.6 240455.9    0
CUBIST  51000.20 217680.6 220028.3 190972.0 226862.3 239288.7    0

Rsquared 
              Min.    1st Qu.     Median      Mean   3rd Qu.      Max. NA's
LM     0.038191154 0.34151126 0.59570117 0.5146318 0.6446897 0.9530659    0
GLM    0.036094303 0.25361489 0.52140959 0.5247907 0.8499520 0.9628829    0
GAM    0.078247285 0.20151523 0.39934478 0.3553497 0.5269644 0.5706765    0
RF     0.007930602 0.01517051 0.11630358 0.1925039 0.1207826 0.7023323    0
XGB    0.041587807 0.09841163 0.11696588 0.2060251 0.2118923 0.5612680    0
NNET   0.018389559 0.02872886 0.07737108 0.1597232 0.1157484 0.5583779    0
ENet   0.170405604 0.55004140 0.62684865 0.5443111 0.6741520 0.7001081    0
MARS   0.007289893 0.04026192 0.24769274 0.3088361 0.3755802 0.8733557    0
SVM    0.044855520 0.06565037 0.24061923 0.2637711 0.4593310 0.5083994    0
kNN    0.001714988 0.10251708 0.41997970 0.3013771 0.4201459 0.5625279    0
RANGER 0.006761073 0.20688755 0.22338943 0.2919063 0.4553028 0.5671905    0
GBM    0.056729819 0.15259735 0.27867315 0.2779509 0.3246546 0.5770994    0
PLS    0.255560646 0.39058052 0.43207618 0.5269205 0.6127945 0.9435907    0
CUBIST 0.004399571 0.20141307 0.52557866 0.5053060 0.8675585 0.9275803    0

lattice::dotplot(results, metric = "RMSE")

# ==============================================================================
# 3) ВЫБОР ЛУЧШЕЙ ПРОГНОСТИЧЕСКОЙ МОДЕЛИ (TIME-SLICE CV НА 3 ГОДА + ХРОНО-ТЕСТ)
# Делим последние годы в тест, внутри train — скользящее окно, h=3.
# ------------------------------------------------------------------------------
# Почему time-slice: временные данные нельзя случайно перемешивать, иначе мы
# «подсматриваем в будущее». Создаём серии обучающих/валидационных окон,
# увеличивая тренировочный период, и тестируем на ближайшем горизонте (3 года).
# ==============================================================================

# 3.1 Данные для time-slice (с YEAR)
model_data <- read.csv("selected_predictors_dataset.csv", header = TRUE, stringsAsFactors = FALSE)
if (!"YEAR" %in% names(model_data)) {
  model_data$YEAR <- seq(1990, by = 1, length.out = nrow(model_data))
}
# Хронологический порядок
model_data <- model_data %>% arrange(YEAR)

# Для временных рядов обычные методы кросс-валидации (случайное разбиение) неприменимы,
# так как это приведет к утечке информации из будущего в прошлое. [[1]]
# Time-slice CV (скользящее окно) имитирует реальную ситуацию прогнозирования:
#   - Мы обучаемся на данных из прошлого
#   - Прогнозируем на несколько шагов вперед
#   - Последовательно сдвигаем окно обучения вперед

# Исходные фичи (исключаем YEAR)
md_for_fit <- model_data %>% select(codTSB, T12, I5, NAOspring, haddock68, R3haddock)

# 3.2 Хронологический holdout (последние годы)
# Идея: отложим ~20% последних лет как полностью внешний тест будущего качества.
n <- nrow(md_for_fit)
holdout_frac <- 0.2
n_test <- max(4, ceiling(n * holdout_frac))
train_ts <- head(md_for_fit, n - n_test)
test_ts  <- tail(md_for_fit, n_test)

# 3.3 Time-slice CV (h=3, expanding window рекомендован: fixedWindow=FALSE)
# initialWindow — размер первого «обучающего» фрагмента; horizon — горизонт
# валидации (здесь 3 года). Далее окно расширяется.
n_train <- nrow(train_ts)
initial_frac <- 0.6
horizon      <- 3
initialWindow <- max(10, floor(initial_frac * n_train))
if (initialWindow + horizon > n_train) initialWindow <- n_train - horizon

slices <- caret::createTimeSlices(1:n_train, initialWindow = initialWindow,
                                  horizon = horizon, fixedWindow = FALSE)
ctrl_ts <- caret::trainControl(method = "cv", index = slices$train, indexOut = slices$test,
                               savePredictions = "final")

# В нашем случае:
#   - horizon = 3: прогнозируем на 3 года вперед
#   - expanding window: размер обучающей выборки увеличивается с каждым шагом
#   - initialWindow: начальный размер обучающей выборки (60% от данных)
# Этот подход наиболее реалистичен для задач прогнозирования временных рядов в гидробиологии.


# 3.4 Обучение (ядро набора, без GBM — он нестабилен на малом n в timeslice)
# Примечание: используем ту же рецептуру, что и в базовом сравнении, но с
# хронологическими срезами.

fit_ts <- function(method, form, data, ctrl, ...) {
  out <- try(caret::train(form, data = data, method = method, trControl = ctrl, ...), TRUE)
  if (inherits(out,"try-error")) NULL else out
}
lm_ts   <- fit_ts("lm",        R3haddock ~ ., train_ts, ctrl_ts)
glm_ts  <- fit_ts("glm",       R3haddock ~ ., train_ts, ctrl_ts, family = Gamma(link="log"))
gam_ts  <- caret::train(x = train_ts[, -which(names(train_ts)=="R3haddock")],
                        y = train_ts$R3haddock, method = gam_spec, trControl = ctrl_ts)

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

rf_ts   <- fit_ts("rf",        R3haddock ~ ., train_ts, ctrl_ts, ntree=1000, tuneGrid=data.frame(mtry=1))
xgb_ts  <- fit_ts("xgbTree",   R3haddock ~ ., train_ts, ctrl_ts, tuneGrid = xgb_grid, verbose = 0)
rgr_ts  <- fit_ts("ranger",    R3haddock ~ ., train_ts, ctrl_ts, tuneLength=3)
nnet_ts <- fit_ts("nnet",      R3haddock ~ ., train_ts, ctrl_ts,
                  preProcess=c("center","scale"),
                  tuneGrid=expand.grid(size=5,decay=0.1), linout=TRUE, trace=FALSE, MaxNWts=5000)

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

svm_ts  <- fit_ts("svmRadial", R3haddock ~ ., train_ts, ctrl_ts, preProcess=c("center","scale"), tuneLength=8)
knn_ts  <- fit_ts("knn",       R3haddock ~ ., train_ts, ctrl_ts, preProcess=c("center","scale"), tuneLength=15)

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 17 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 19 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 21 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 23 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 25 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 27 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 29 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 31 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.91776861098288,
-0.63635090426016, : k = 33 exceeds number 15 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 17 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 19 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 21 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 23 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 25 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 27 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 29 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 31 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-1.97509275968546,
-0.649515010512528, : k = 33 exceeds number 16 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 19 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 21 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 23 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 25 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 27 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 29 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 31 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.03591114922526,
-0.667021070399982, : k = 33 exceeds number 17 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 19 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 21 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 23 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 25 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 27 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 29 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 31 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-2.09685135650479,
-0.723233808010845, : k = 33 exceeds number 18 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 21 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 23 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 25 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 27 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 29 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 31 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.87781906605955,
-0.736313775515053, : k = 33 exceeds number 19 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 21 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 23 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 25 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 27 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 29 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 31 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.59299667707382,
-0.714710653867683, : k = 33 exceeds number 20 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 23 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 25 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 27 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 29 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 31 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.44471091878981,
-0.719611760680745, : k = 33 exceeds number 21 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 23 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 25 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 27 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 29 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 31 exceeds number 22 of patterns

Warning in knnregTrain(train = structure(c(-1.35845826709587,
-0.734816317064085, : k = 33 exceeds number 22 of patterns

Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.

enet_ts <- fit_ts("glmnet",    R3haddock ~ ., train_ts, ctrl_ts, preProcess=c("center","scale"), tuneLength=10)
mars_ts <- fit_ts("earth",     R3haddock ~ ., train_ts, ctrl_ts, tuneLength=10)
pls_ts  <- fit_ts("pls",       R3haddock ~ ., train_ts, ctrl_ts, preProcess=c("center","scale"), tuneLength=10)
cub_ts  <- fit_ts("cubist",    R3haddock ~ ., train_ts, ctrl_ts, tuneLength=5)

models_ts <- list(LM=lm_ts, GLM=glm_ts, GAM=gam_ts, RF=rf_ts, XGB=xgb_ts, RANGER=rgr_ts,
                  NNET=nnet_ts, SVM=svm_ts, kNN=knn_ts, ENet=enet_ts, MARS=mars_ts, PLS=pls_ts, CUBIST=cub_ts)
models_ts <- models_ts[!vapply(models_ts, is.null, logical(1))]

# 3.5 Ранжирование по time-slice CV и по хронологическому тесту
# Сначала ранжируем по средним ошибкам на валидационных срезах, затем — по внешнему тесту.
cv_metrics <- function(m) {
  if (is.null(m$pred) || !"Resample" %in% names(m$pred)) return(c(RMSE=NA, MAE=NA))
  by_slice <- m$pred %>% group_by(Resample) %>%
    summarise(RMSE=rmse(obs,pred), MAE=mae(obs,pred), .groups="drop")
  c(RMSE = mean(by_slice$RMSE, na.rm = TRUE), MAE = mean(by_slice$MAE, na.rm = TRUE))
}
cv_rank <- do.call(rbind, lapply(models_ts, cv_metrics)) %>% as.data.frame()
cv_rank$Model <- rownames(cv_rank)
cv_rank <- cv_rank[is.finite(cv_rank$RMSE), ] %>% relocate(Model) %>% arrange(RMSE, MAE)
cat("\nTime-slice CV (h=3), средние RMSE/MAE:\n"); print(cv_rank)

Time-slice CV (h=3), средние RMSE/MAE:

        Model     RMSE      MAE
SVM       SVM 227929.6 191921.6
kNN       kNN 234897.8 197135.0
ENet     ENet 250989.0 214248.7
XGB       XGB 277153.4 248532.6
RANGER RANGER 280255.1 249992.1
GLM       GLM 280259.5 237186.9
NNET     NNET 296856.1 264368.2
PLS       PLS 302968.3 274707.6
RF         RF 303710.0 263019.5
CUBIST CUBIST 314443.4 281437.9
LM         LM 370298.1 340883.8
MARS     MARS 427624.1 378476.5
GAM       GAM 714951.0 625520.3

preds_ts <- lapply(models_ts, function(m) try(predict(m, newdata = test_ts), TRUE))
keep <- vapply(preds_ts, function(p) is.numeric(p) && length(p)==nrow(test_ts) && all(is.finite(p)), logical(1))
preds_ts <- preds_ts[keep]
test_rank <- do.call(rbind, lapply(names(preds_ts), function(nm){
  data.frame(Model=nm, t(metrics_vec(test_ts$R3haddock, preds_ts[[nm]])), row.names = NULL)
})) %>% arrange(RMSE, MAE)
cat("\nХронологический тест (последние годы), RMSE/MAE/R2:\n"); print(test_rank)

Хронологический тест (последние годы), RMSE/MAE/R2:

    Model     RMSE      MAE         R2      MAPE     sMAPE
1  CUBIST 148248.4 107629.3  0.5774780  54.93724  38.06342
2      LM 156940.0 129604.7  0.5264822  97.30313  49.10292
3     GAM 158131.0 125038.0  0.5192677  51.33030  42.12742
4     PLS 176786.2 138249.0  0.3991499 123.21995  50.93728
5     GLM 182047.7 141692.4  0.3628531  73.71724  50.32932
6      RF 185485.8 148117.2  0.3385600  96.29497  52.15688
7     kNN 187966.9 143164.4  0.3207460  71.77551  50.80097
8    ENet 195260.5 161191.7  0.2670102 115.55953  56.27135
9  RANGER 208161.4 171717.0  0.1669526 115.93679  58.76935
10    SVM 250994.6 197801.6 -0.2111500  84.52010  67.59732
11   MARS 273186.2 208579.6 -0.4347849  85.29788  71.46314
12    XGB 288167.4 233693.4 -0.5964639 102.19238  78.99024
13   NNET 345227.3 291463.2 -1.2912878 233.62022 102.82204

# ==============================================================================
# 4) ПРОГНОЗ 2022–2024 (АНСАМБЛЬ CUBIST+LM) И ГРАФИК 1990–2024 С ДИ
# Прогнозные линии (медиана и ДИ) — пунктир; исторические — сплошные.
# Можно задать свои сценарии предикторов (user_future); по умолчанию — средние.
# ------------------------------------------------------------------------------
# Логика ансамбля: комбинируем сильную нелинейную модель (Cubist) с простой и
# устойчивой линейной (LM). Веса можно настраивать. Доверительные интервалы
# получаем эмпирически из распределения остатков (простая и наглядная эвристика).
# ==============================================================================

# 4.1 Полные модели для прогноза (на всех данных) и вес ансамбля
model_data <- read.csv("selected_predictors_dataset.csv", header = TRUE, stringsAsFactors = FALSE)
if (!"YEAR" %in% names(model_data)) {
  model_data$YEAR <- seq(1990, by = 1, length.out = nrow(model_data))
}
model_data <- model_data %>% arrange(YEAR)

cubist_full <- caret::train(R3haddock ~ codTSB + T12 + I5 + NAOspring + haddock68,
                            data = model_data, method = "cubist",
                            trControl = caret::trainControl(method="none"),
                            tuneGrid = if (exists("cubist_model")) cubist_model$bestTune else NULL,
                            tuneLength = if (exists("cubist_model")) 1 else 5)

lm_full <- caret::train(R3haddock ~ codTSB + T12 + I5 + NAOspring + haddock68,
                        data = model_data, method = "lm",
                        trControl = caret::trainControl(method="none"))

alpha_opt <- if (exists("alpha_opt")) alpha_opt else 0.75
predict_ensemble <- function(newdata, alpha = alpha_opt) {
  alpha * predict(cubist_full, newdata) + (1 - alpha) * predict(lm_full, newdata)
}

# Ансамбль моделей часто дает более точные и устойчивые прогнозы, чем отдельные модели. [[8]]
# В нашем случае:
#   - CUBIST: мощная модель, основанная на правилах, хорошо работающая с табличными данными
#   - LM: простая интерпретируемая модель, устойчивая к шуму
#   - alpha_opt = 0.75: веса ансамбля (75% CUBIST, 25% LM), оптимизированные ранее (см. скрипт "ENS_WEIGHT.R")
# Комбинирование моделей с разными сильными сторонами снижает риск систематических ошибок.

# 4.2 Остатки для ДИ (из CV, если есть; иначе — по фитам)
# Эмпирические квантилы остатков дают «практические» интервалы прогноза без
# предположения нормальности ошибок (хотя строгий PI требует аккуратности).
get_residuals_for_pi <- function() {
  if (exists("lm_model") && exists("cubist_model") &&
      !is.null(lm_model$pred) && !is.null(cubist_model$pred)) {
    pl <- lm_model$pred %>% select(Resample,rowIndex,obs,p_lm=pred)
    pc <- cubist_model$pred %>% select(Resample,rowIndex,p_cu=pred)
    inner_join(pl, pc, by=c("Resample","rowIndex")) %>%
      mutate(p_ens = alpha_opt * p_cu + (1 - alpha_opt) * p_lm,
             resid = obs - p_ens) %>%
      pull(resid) %>% .[is.finite(.)]
  } else {
    model_data$R3haddock - predict_ensemble(model_data)
  }
}
resids <- get_residuals_for_pi()
q025 <- as.numeric(quantile(resids, 0.025, na.rm = TRUE))
q250 <- as.numeric(quantile(resids, 0.250, na.rm = TRUE))
q750 <- as.numeric(quantile(resids, 0.750, na.rm = TRUE))
q975 <- as.numeric(quantile(resids, 0.975, na.rm = TRUE))

# Доверительные интервалы (ДИ) показывают неопределенность прогноза.
# Мы используем квантили остатков из кросс-валидации для построения ДИ:
#   - PI50 (50% интервал): между 25-м и 75-м процентилями
#   - PI95 (95% интервал): между 2.5-м и 97.5-м процентилями
# Это непараметрический подход, не требующий предположений о нормальности ошибок.

# 4.3 Сценарии будущего (по умолчанию — средние; можно переопределить user_future)
fc_start <- 2022
pred_cols <- c("codTSB","T12","I5","NAOspring","haddock68")
train_period <- model_data %>% filter(YEAR > 1989 & YEAR < fc_start)
mu <- train_period %>% summarise(across(all_of(pred_cols), ~mean(.x, na.rm = TRUE))) %>% as.list()

# Пример пользовательского сценария:
# user_future <- tibble::tribble(
#   ~YEAR, ~codTSB, ~T12, ~I5, ~NAOspring, ~haddock68,
#   2022, 2100000, 5.1, 48,  0.3, 120000,
#   2023, 2050000, 4.8, 50, -0.1, 115000,
#   2024, 2150000, 5.0, 47,  0.2, 118000
# )
if (!exists("user_future")) user_future <- NULL

build_future <- function(years, mu, user_df=NULL) {
  df <- tibble::tibble(YEAR = years)
  for (v in pred_cols) df[[v]] <- mu[[v]]
  if (!is.null(user_df)) {
    for (i in seq_len(nrow(user_df))) {
      yr <- user_df$YEAR[i]
      if (yr %in% years) {
        idx <- which(df$YEAR == yr)
        for (v in intersect(pred_cols, names(user_df))) {
          val <- user_df[[v]][i]
          if (!is.na(val)) df[[v]][idx] <- val
        }
      }
    }
  }
  df
}
future_years <- fc_start:2024
scenario_future <- build_future(future_years, mu, user_future)

# 4.4 Прогноз и таблица ДИ
pred_future <- predict_ensemble(scenario_future)
forecast_tbl <- tibble::tibble(
  YEAR      = scenario_future$YEAR,
  pred_mean = as.numeric(pred_future),
  PI50_low  = pred_future + q250, PI50_high = pred_future + q750,
  PI95_low  = pred_future + q025, PI95_high = pred_future + q975
)

#### Таблица прогноза 2022–2024
print(forecast_tbl)

# A tibble: 3 x 6
   YEAR pred_mean PI50_low PI50_high PI95_low PI95_high
  <int>     <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
1  2022   253815.   63865.   536150.  -46219.   668189.
2  2023   253815.   63865.   536150.  -46219.   668189.
3  2024   253815.   63865.   536150.  -46219.   668189.

# По умолчанию мы используем средние значения предикторов для прогноза.
# Однако вы можете определить собственный сценарий (user_future), указав конкретные значения
# для каждого года и каждого предиктора. Это позволяет моделировать различные экологические сценарии. 

# 4.5 Непрерывный ряд 1990–2024 и график: ленты сплошные; линии медианы/ДИ — сплошные до 2021, пунктир с 2022
pred_df <- bind_rows(
  model_data %>% select(YEAR, all_of(pred_cols)),
  scenario_future
) %>% distinct(YEAR, .keep_all = TRUE) %>% arrange(YEAR)

pred_df$Pred      <- as.numeric(predict_ensemble(pred_df))
pred_df$PI50_low  <- pred_df$Pred + q250
pred_df$PI50_high <- pred_df$Pred + q750
pred_df$PI95_low  <- pred_df$Pred + q025
pred_df$PI95_high <- pred_df$Pred + q975

hist_df <- model_data %>% select(YEAR, R3haddock)

ggplot() +
  geom_ribbon(data = pred_df, aes(x = YEAR, ymin = PI95_low, ymax = PI95_high),
              fill = "grey80", alpha = 0.25) +
  geom_ribbon(data = pred_df, aes(x = YEAR, ymin = PI50_low, ymax = PI50_high),
              fill = "grey60", alpha = 0.35) +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = PI95_low),
            color = "grey45", linewidth = 0.6) +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = PI95_high),
            color = "grey45", linewidth = 0.6) +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = PI50_low),
            color = "grey35", linewidth = 0.6) +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = PI50_high),
            color = "grey35", linewidth = 0.6) +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = PI95_low),
            color = "grey45", linewidth = 0.6, linetype = "dashed") +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = PI95_high),
            color = "grey45", linewidth = 0.6, linetype = "dashed") +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = PI50_low),
            color = "grey35", linewidth = 0.6, linetype = "dashed") +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = PI50_high),
            color = "grey35", linewidth = 0.6, linetype = "dashed") +
  geom_line(data = subset(pred_df, YEAR < fc_start), aes(x = YEAR, y = Pred),
            color = "steelblue4", linewidth = 1) +
  geom_line(data = subset(pred_df, YEAR >= fc_start-1), aes(x = YEAR, y = Pred),
            color = "steelblue4", linewidth = 1, linetype = "dashed") +
  geom_point(data = hist_df, aes(x = YEAR, y = R3haddock),
             color = "black", size = 2, alpha = 0.9) +
  scale_x_continuous(expand = expansion(mult = c(0, 0))) +
  labs(
    title = "Пополнение R3haddock: факт (1990–2021) и прогноз (2022–2024)\nАнсамбль CUBIST+LM; непрерывные ДИ, прогноз — пунктир",
    x = "Год", y = "R3haddock"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

# На графике:
#   - Черные точки: исторические данные (1990-2021)
#   - Сплошная синяя линия: прогнозные значения (1990-2021)
#   - Пунктирная синяя линия: прогноз на 2022-2024
#   - Серые ленты: 50% и 95% доверительные интервалы
# Такая визуализация позволяет легко интерпретировать как исторические данные, 
# так и будущие прогнозы с учетом неопределенности.