Прогнозирование вспышек инфекций

MVP мониторинга ОРВИ на синтетических данных: пайплайн работает, а выводов о модели пока нет

Я собрал MVP системы поддержки принятия решений для эпидемиологического мониторинга. Цель была технической: проверить полный путь от временного ряда по территории до расчёта риска и выдачи результата в диалоговом интерфейсе.

Важное ограничение: эпидемиологическая динамика в этом MVP синтетическая. Реальными были названия муниципалитетов, площадь, численность населения, плотность и температура из Яндекс Погоды. Показатели вакцинации, число случаев, заболеваемость и целевая метка генерировались программно.

Поэтому приведённые метрики не говорят ничего о способности модели прогнозировать реальные вспышки ОРВИ. На таких данных нельзя делать вывод, что одна архитектура лучше другой: модель в первую очередь учится закономерностям, заложенным в генератор синтетического ряда.

Тем не менее стенд оказался полезен. Он позволил проверить построение временных окон, хронологическое разбиение, нормализацию, обучение, инференс, пороги уведомлений и обработку редкого положительного класса.

Что я хотел проверить

В будущем система должна помогать эпидемиологам и сотрудникам МИАЦ замечать потенциально опасную динамику раньше, чем она станет очевидной в отчётах и дашбордах.

На уровне MVP задача была сформулирована так:

По истории показателей за восемь недель оценить вероятность того, что через две недели территория превысит заданный модельный эпидемический порог.

Система не должна ставить диагнозы, автоматически вводить ограничения или принимать решения вместо специалиста. Её задача, показать расчётный уровень риска и помочь обратить внимание на территорию, которую нужно проверить.

В будущем реальная версия должна получать агрегированные показатели из источников Роспотребнадзора. На этапе MVP я проверял только инженерную схему работы.

Данные: реальные параметры и синтетическая динамика

Для MVP мне был нужен стенд, на котором можно проверить полный ML-пайплайн до подключения реальных источников статистики.

В качестве географической основы я использовал муниципалитеты Адыгеи, Краснодарского края и Ростовской области. Это компактный демонстрационный фрагмент Юга России.

Реальными были:

• названия муниципалитетов;

• площадь территории;

• численность населения;

• рассчитанная плотность населения;

• средняя температура из Яндекс Погоды.

Остальные показатели генерировались программно:

• охват вакцинацией;

• абсолютное число случаев;

• заболеваемость на 100 тыс. населения;

• целевая метка события.

Для каждого муниципалитета генерировался ряд длиной 260 недель. Температура менялась по сезонной синусоиде с добавлением случайного шума. Число случаев строилось через распределение Пуассона: его ожидаемое значение зависело от населения, плотности, условного охвата вакцинацией и сезонного коэффициента.

Упрощённо генератор выглядел так:

seasonality_factor = (
    1.0
    + 0.8 * np.sin((week - 10) / 52 * 2 * np.pi)
)

density_multiplier = (
    np.log10(max(current_density, 1)) / 2
)

susceptibility = (
    (1 - vaccination_rate)
    * density_multiplier
)

lambda_cases = (
    (population / 10000)
    * susceptibility
    * seasonality_factor
    * 10
)

cases_abs = np.random.poisson(
    max(lambda_cases, 1)
)

В периоды высокой сезонности генератор с небольшой вероятностью дополнительно увеличивал число случаев в 2,5-4 раза. Так в ряду появлялись редкие пики, на которых можно было проверить разметку и обучение классификатора.

Положительный класс был редким. В обучающей части было порядка 200 положительных окон при примерно 9 200 отрицательных. Это стало главным источником дисбаланса в задаче.

Такой датасет подходит для отладки инженерного контура. Он не подходит для оценки медицинской, эпидемиологической или прогностической эффективности модели.

Как формировалась целевая метка

Положительный класс не брался из реальных эпидемиологических расследований. Он строился программно по эвристике, согласованной с сотрудником МИАЦ.

Территория считалась перешедшей модельный эпидемический порог, если выполнялись два условия:

1. Заболеваемость на 100 тыс. населения превышала восьминедельное скользящее среднее более чем на два стандартных отклонения.

2. Сам показатель был выше 50 случаев на 100 тыс. населения.

В коде это выглядит так:

def label_outbreaks(group):
    rolling_mean = (
        group["cases_per_100k"]
        .rolling(window=8, min_periods=1)
        .mean()
    )

    rolling_std = (
        group["cases_per_100k"]
        .rolling(window=8, min_periods=1)
        .std()
        .fillna(0)
    )

    epidemic_threshold = (
        rolling_mean
        + 2 * rolling_std
    )

    group["is_outbreak"] = (
        (group["cases_per_100k"] > epidemic_threshold)
        & (group["cases_per_100k"] > 50)
    ).astype(int)

    return group

Есть важное методологическое ограничение. Целевая метка строится по cases_per_100k, и этот же признак входит во входное окно модели. Поэтому задача частично сводится к изучению правила, использованного в генераторе и разметке.

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

Хронологический split

В итоговом наборе получилось:

• 9 405 строк для обучения;

• 2 295 строк для validation.

Я не перемешивал данные случайным образом. Разделение было хронологическим:

cutoff_date = pd.to_datetime("2022-12-31")

train_df = df[df["date"] <= cutoff_date]
val_df = df[df["date"] > cutoff_date]

Для временного ряда это важнее случайного split. Если перемешать недели, модель может увидеть в обучении фрагменты будущей динамики той же территории. Тогда validation будет выглядеть лучше, чем сценарий реального прогноза.

Отдельного test-набора в MVP пока нет. Все цифры ниже относятся к validation-выборке. Это ещё одно ограничение: нельзя долго настраивать архитектуру на одной validation-части и потом считать её независимой оценкой.

Окно истории и горизонт прогноза

На вход модель получает историю из восьми недель по шести признакам:

После подготовки один объект имеет форму:

[1, 8, 6]

где:

• 1 это размер batch;

• 8 это число недель в истории;

• 6 это число признаков.

Прогноз строится на две недели вперёд. Это задаётся параметром horizon=2:

for i in range(
    len(group) - lookback - horizon + 1
):
    self.x_data.append(
        features[i : i + lookback]
    )

    self.y_data.append(
        targets[i + lookback + horizon - 1]
    )

То есть модель получает значения за восемь недель и пытается определить, сработает ли модельная метка через две недели.

Нормализация и сохранение артефактов

Численность населения измеряется сотнями тысяч или миллионами, плотность может быть десятками или тысячами человек на квадратный километр, температура измеряется градусами, а охват вакцинацией находится в диапазоне от 0 до 1.

Без нормализации признаки имеют несопоставимый масштаб. Поэтому я использовал StandardScaler.

if is_train:
    self.scaler = StandardScaler()

    df[self.feature_cols] = (
        self.scaler.fit_transform(
            df[self.feature_cols]
        )
    )
else:
    self.scaler = scaler

    df[self.feature_cols] = (
        self.scaler.transform(
            df[self.feature_cols]
        )
    )

После обучения scaler сохраняется отдельно:

joblib.dump(
    dataset_train.scaler,
    "scaler.pkl",
)

В рабочем режиме нельзя создавать новый scaler для каждого расчёта. Нужно применять тот же объект, который был обучен на train-выборке. Иначе входные данные попадут в другую шкалу, и результат модели потеряет смысл. StandardScaler рассчитывает среднее и стандартное отклонение по данным обучения, а затем использует их для преобразования новых наблюдений.[scikit-learn]

Почему я выбрал Transformer

Задача выглядит как временная последовательность: есть восемь недель наблюдений, изменение температуры, рост случаев, условный охват вакцинацией и параметры территории.

Поэтому первым вариантом стал компактный Temporal Transformer:

• d_model=64;

• n_heads=4;

• два encoder layer;

• dropout=0.15;

• классификация по embedding последней недели.

Transformer не знает порядок элементов последовательности сам по себе. Поэтому к embedding каждой недели добавляется positional encoding.

import math
import torch
import torch.nn as nn


class PositionalEncoding(nn.Module):
    def __init__(self, d_model, max_len=50):
        super().__init__()

        pe = torch.zeros(max_len, d_model)

        position = torch.arange(
            0,
            max_len,
            dtype=torch.float,
        ).unsqueeze(1)

        div_term = torch.exp(
            torch.arange(
                0,
                d_model,
                2,
            ).float()
            * (-math.log(10000.0) / d_model)
        )

        pe[:, 0::2] = torch.sin(
            position * div_term
        )

        pe[:, 1::2] = torch.cos(
            position * div_term
        )

        self.register_buffer(
            "pe",
            pe.unsqueeze(0),
        )

    def forward(self, x):
        return x + self.pe[:, :x.size(1), :]

Сама модель:

class EpidemicTransformerFast(nn.Module):
    def __init__(
        self,
        num_features=6,
        d_model=64,
        n_heads=4,
        num_layers=2,
        dropout=0.15,
    ):
        super().__init__()

        self.input_projection = nn.Linear(
            num_features,
            d_model,
        )

        self.pos_encoder = PositionalEncoding(
            d_model
        )

        encoder_layer = nn.TransformerEncoderLayer(
            d_model=d_model,
            nhead=n_heads,
            dim_feedforward=d_model * 2,
            dropout=dropout,
            batch_first=True,
        )

        self.transformer_encoder = (
            nn.TransformerEncoder(
                encoder_layer,
                num_layers,
            )
        )

        self.classifier = nn.Sequential(
            nn.Linear(d_model, d_model // 2),
            nn.GELU(),
            nn.Dropout(dropout),
            nn.Linear(d_model // 2, 1),
        )

    def forward(self, x):
        x = self.input_projection(x)
        x = self.pos_encoder(x)

        out = self.transformer_encoder(x)

        last_week = out[:, -1, :]

        return self.classifier(
            last_week
        ).squeeze(-1)

Модель возвращает логит. Вероятность рассчитывается через sigmoid:

probability = torch.sigmoid(logit).item()

Дисбаланс: модель либо поднимает тревоги, либо молчит

Положительный класс был редким. Если не учитывать дисбаланс, модели выгодно почти всегда предсказывать нулевой класс.

Первой попыткой был BCEWithLogitsLoss с высоким весом положительного класса:

criterion = nn.BCEWithLogitsLoss(
    pos_weight=torch.tensor([45.0]).to(device)
)

Логика была понятной: пропуск положительного события должен стоить модели дороже.

На практике модель начала чаще поднимать тревоги. Для будущего интерфейса это плохой сценарий. Если система слишком часто предупреждает пользователя, он перестаёт воспринимать уведомления всерьёз.

Я также попробовал более ёмкую конфигурацию Transformer:

• d_model=128;

• n_heads=8;

• dropout=0.3;

• pos_weight=40.

Лучший сохранённый артефакт показал:

Этот эксперимент не доказывает, что большая модель хуже маленькой. На синтетических данных сложно отделить влияние архитектуры от влияния генератора, разметки, случайного seed и порога классификации.

Эксперимент с Focal Loss

Следующей попыткой был Focal Loss. Он уменьшает вклад простых примеров и сильнее фокусирует обучение на сложных.

class FocalLoss(nn.Module):
    def __init__(self, alpha=0.85, gamma=2.0):
        super().__init__()

        self.alpha = alpha
        self.gamma = gamma

    def forward(self, inputs, targets):
        bce = F.binary_cross_entropy_with_logits(
            inputs,
            targets,
            reduction="none",
        )

        pt = torch.exp(-bce)

        loss = (
            self.alpha
            * (1 - pt) ** self.gamma
            * bce
        )

        return loss.mean()

В этом эксперименте Focal Loss не помог. На validation модель перестала выделять положительный класс, а recall упал до 0,000.

Скорее всего, модель подстроилась под структуру синтетических данных. Точную причину этого поведения я не выяснял: для такого вывода понадобились бы повторные запуски с разными seed, анализ распределений логитов, PR-кривые и отдельные эксперименты с параметрами alpha, gamma и порогом классификации.

Поэтому этот результат я интерпретирую не как свойство Focal Loss вообще, а как наблюдение для конкретной конфигурации MVP.

Финальная конфигурация MVP

В последней итерации я вернулся к меньшей архитектуре и снизил вес положительного класса до 15:

criterion = nn.BCEWithLogitsLoss(
    pos_weight=torch.tensor([15.0]).to(device)
)

Порог классификации на validation был установлен на 0,30:

preds = (probs > 0.30).float()

Лучшие веса сохранялись по F1:

if val_f1 > best_val_f1:
    best_val_f1 = val_f1

    torch.save(
        model.state_dict(),
        "best_epidemic_transformer_final.pth",
    )

Итоговые метрики этой конфигурации:

Это слабые показатели, и для рабочего мониторинга их недостаточно.

Однако ещё важнее другое: по этому эксперименту нельзя делать вывод, что Transformer не подходит для задачи. Паттерн целевого события, сезонность и редкие пики определяет генератор синтетических данных. Модель может учиться структуре генератора, а не закономерностям реальной эпидемиологической динамики.

MVP позволил проверить технический контур и увидеть чувствительность обучения к pos_weight, порогу классификации и функции потерь. Но сравнение архитектур имеет смысл только после появления реальных временных рядов и независимой разметки событий.

Что сравнивать на реальных данных

CatBoost я пока не запускал, поэтому говорить о его преимуществе над Transformer рано.

Когда появятся реальные агрегированные данные, сравнение стоит строить от простого к сложному:

1. Базовое правило по динамике заболеваемости.

2. Логистическая регрессия.

3. Градиентный бустинг, например CatBoost.

4. Последовательностная модель, если временная структура действительно добавляет сигнал.

CatBoost выглядит разумным кандидатом для первого сильного baseline, поскольку текущие признаки табличные, а длина истории небольшая. Но это гипотеза, а не результат проведённого сравнения.

Следующая итерация должна быть посвящена не увеличению числа attention-heads, а сбору данных и формированию признаков.

Например, можно проверить пространственные лаги:

• динамику случаев в соседних муниципалитетах;

• темпы роста заболеваемости у соседей;

• связь территорий через общие границы;

• транспортные потоки;

• различия в плотности населения и сезонности.

Диалоговый интерфейс MVP

Для демонстрации сценария я сделал простой чат-интерфейс. Он последовательно запрашивает данные по территории:

1. название территории;

2. численность населения;

3. плотность населения;

4. охват вакцинацией;

5. среднюю температуру;

6. число зарегистрированных случаев.

После каждого шага выполняется валидация. Например, население не может быть отрицательным:

try:
    population = int(message.text)

    if population <= 0:
        raise ValueError

except ValueError:
    return await message.answer(
        "Введите положительное целое число."
    )

Для охвата вакцинацией проверяется диапазон от 0 до 1:

try:
    vaccination = float(
        message.text.replace(",", ".")
    )

    if not 0.0 <= vaccination <= 1.0:
        raise ValueError

except ValueError:
    return await message.answer(
        "Введите значение от 0.0 до 1.0."
    )

В демонстрационном режиме ручной ввод не заменяет полноценную историю. После ввода одной недели она повторяется восемь раз, чтобы сформировать тензор нужной формы:

full_history = [current_week for _ in range(8)]

Это заглушка для демонстрации идеи интерфейса. В рабочей версии восемь недель должны автоматически собираться из внешнего источника статистики, а не восстанавливаться из одного введённого значения.

Уровни риска и ложные тревоги

Результат модели нельзя показывать пользователю как окончательное утверждение. Поэтому я добавил три зоны интерпретации:

В коде:

if probability >= 0.70:
    tier = "🔴 ВЫСОКИЙ РИСК"

    action = (
        "Проверить исходные показатели и "
        "передать информацию эпидемиологу."
    )

elif probability >= 0.30:
    tier = "🟡 ПОВЫШЕННЫЙ РИСК"

    action = (
        "Проверить динамику случаев и "
        "готовность организационных ресурсов."
    )

else:
    tier = "🟢 ШТАТНЫЙ РЕЖИМ"

    action = (
        "Продолжить штатный мониторинг."
    )

Пороги 0,30 и 0,70 пока не утверждены. Это продуктовая гипотеза, которую нужно обсуждать со специалистами и проверять на реальных данных.

Проблема здесь не только в метриках. Слишком частые жёлтые и красные сигналы создают Alert Fatigue. Если уведомления регулярно оказываются ложными, пользователи перестают реагировать даже на важные сообщения.

Что ещё не хватает в оценке

В MVP я сохранял recall, precision и F1 при фиксированном пороге. Для редкого класса этого недостаточно.

В следующей итерации нужно добавить:

• PR-AUC;

• PR-кривую;

• матрицу ошибок;

• распределение вероятностей по классам;

• калибровку вероятностей;

• несколько запусков с разными random seed;

• независимый test-набор.

Расчёт PR-AUC и матрицы ошибок можно добавить после validation:

from sklearn.metrics import (
    average_precision_score,
    confusion_matrix,
)

val_targets = []
val_probabilities = []

model.eval()

with torch.no_grad():
    for batch_x, batch_y in val_loader:
        batch_x = batch_x.to(device)

        logits = model(batch_x)
        probabilities = torch.sigmoid(logits)

        val_targets.extend(
            batch_y.numpy()
        )

        val_probabilities.extend(
            probabilities.cpu().numpy()
        )

threshold = 0.30

val_predictions = (
    np.array(val_probabilities) > threshold
).astype(int)

pr_auc = average_precision_score(
    val_targets,
    val_probabilities,
)

tn, fp, fn, tp = confusion_matrix(
    val_targets,
    val_predictions,
).ravel()

print(f"PR-AUC: {pr_auc:.4f}")
print(f"TN={tn}, FP={fp}, FN={fn}, TP={tp}")

Добавлять конкретные значения PR-AUC и матрицы ошибок в статью можно только после повторного запуска расчёта.

Что нужно до работы с реальными данными

MVP позволил проверить цепочку от данных до результата, но до реального пилота ещё далеко.

Минимальный список следующей итерации:

• получить согласованный источник реальных агрегированных данных;

• накопить историю по территориям;

• отдельно описать правила разметки событий;

• сформировать независимый test-набор;

• сравнить простые baseline-модели, CatBoost и Transformer;

• проверить устойчивость к пропускам и задержкам статистики;

• добавить пространственные признаки;

• откалибровать вероятности;

• согласовать пороги риска со специалистами;

• реализовать защищённый backend;

• провести правовую и регуляторную оценку будущей системы.

Модель не должна автоматически вводить ограничения, менять маршрутизацию или принимать решения о ресурсах. Она может только подсказать, где стоит проверить данные. Решение остаётся за эпидемиологом.

Итоги

MVP позволил собрать и проверить полный технический контур: генерацию рядов, формирование восьминедельных окон, хронологическое разбиение, нормализацию, обучение, сохранение артефактов, инференс и вывод результата через диалоговый интерфейс.

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

При этом стенд дал несколько практических наблюдений:

• большой pos_weight быстро превращает модель в источник ложных тревог;

• смена функции потерь не гарантирует улучшения, в одном эксперименте Focal Loss привёл к нулевому recall;

• порог классификации существенно меняет поведение системы;

• даже демонстрационный интерфейс требует валидации ввода и явного разделения логики интерфейса и ML-инференса;

• для будущей системы критичнее получить реальные данные и независимую метку события, чем усложнять архитектуру.

Следующий этап, получить реальные агрегированные ряды, подготовить независимый test-набор и сравнить простые baseline-модели с CatBoost и Transformer. Только после этого можно будет обсуждать качество прогноза и пригодность системы для пилота.

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.