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

Моя цель - предложение широкого ассортимента товаров и услуг на постоянно высоком качестве обслуживания по самым выгодным ценам.

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-части и потом считать её независимой оценкой.

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

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

Признак

Смысл

population

Численность населения

density

Плотность населения

vaccination_rate

Доля вакцинированного населения

avg_temp

Средняя температура за неделю

cases_abs

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

cases_per_100k

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

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

[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.

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

Метрика

Значение

Recall

0,048

F1

0,069

Этот эксперимент не доказывает, что большая модель хуже маленькой. На синтетических данных сложно отделить влияние архитектуры от влияния генератора, разметки, случайного 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",
    )

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

Метрика

Значение

Recall

0,286

Precision

0,047

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

Однако ещё важнее другое: по этому эксперименту нельзя делать вывод, что 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)]

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

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

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

Вероятность

Статус

Действие

Меньше 0,30

Источник: https://habr.com/ru/articles/1089990/


Интересные статьи

Интересные статьи

Привет! Мы — партнеры Битрикс24 и создаем для бизнеса приложения с генеративным ИИ, интегрируя их с корпоративными порталами.В этой статье мы расскажем о созданном нами решении, которое позволяет любы...
Меня зовут Михаил Кузнецов, я product owner в команде, которая развивает внутреннюю платформу разработки Mindbox. В этой статье я расскажу, как мы отказались от легаси-фреймворка, который пронизывал в...
Так уж повелось на российском рынке маркетинга (а может и на мировом тоже), что каждое агентство утверждает, что делает анализ конкурентов перед запуском рекламы.Утверждает каждое, а делают единицы. ...
В данной статье мы рассмотрим систему аутентификации пользователей и внешних систем в личном кабинете через сервер аутентификации Blitz Identity Provider. Согласно требованиям проекта, который мы...
В 1С-Битрикс: Управление сайтом (как и в Битрикс24) десятки, если не сотни настраиваемых типов данных (или сущностей): инфоблоки, пользователи, заказы, склады, форумы, блоги и т.д. Стр...