DengAI: прогнозирование вспышек лихорадки денге

Постановка задачи

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

Наша цель — прогнозировать количество случаев денге, которые будут возникать в Сан-Хуане и Икитосе каждую неделю на основе переменных окружающей среды, таких как изменения температуры, осадки и растительность.

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

Исследовательский анализ данных

Набор данных включает еженедельные данные за 18 лет, собранные в Сан-Хуане, и еженедельные данные за 10 лет, собранные в Икитосе. Он содержит 20 функций, описывающих погодные переменные, такие как количество осадков, процент влажности, температура воздуха, температура точки росы и температурные диапазоны. Целевой переменной является число случаев лихорадки денге в неделю.

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

На визуализации мы видим, что масштабы вспышек лихорадки денге в Сан-Хуане со временем уменьшились. Это может быть связано с рядом факторов, таких как городское развитие и соответствующее сокращение зеленых насаждений, где могут размножаться комары, экзогенными факторами, такими как модернизация жилья и более широкое использование противомоскитных сеток, или другими неизвестными переменными. Одной из характеристик в нашем наборе данных является NDVI (нормализованный разностный индекс растительности), который дает оценку плотности растительности в районе на основе спутниковых изображений. График этого показателя во времени потенциально может дать представление о том, могло ли сокращение растительности способствовать уменьшению вспышек лихорадки денге.

Мы видим, что NDVI в Сан-Хуане демонстрирует особенно значительную тенденцию к снижению, что может свидетельствовать о том, что наше предположение о корреляции между NDVI и случаями лихорадки денге верно. Кроме того, мы видим, что Сан-Хуан кажется более городским, чем Икитос, о чем свидетельствует его более низкий общий NDVI. Это потенциально может повлиять на поведение вируса денге в обоих местах, что является одной из причин, по которой нам потребуется обучать отдельные модели для каждого города.

Построив целевую переменную по неделям, мы можем лучше понять сезонность наших данных.

Визуально видно, что в наших данных присутствует сезонный аспект. В Сан-Хуане обычно выпадает больше осадков во время сезона ураганов, который длится с конца августа по октябрь, в то время как в Икитосе обычно выпадает больше осадков в зимние месяцы. Учитывая, что два города расположены в разных полушариях, сезонные всплески заболеваемости не совпадают друг с другом.

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

Как видно из визуализации, заболеваемость лихорадкой денге в Сан-Хуане, как правило, снижается в более засушливые месяцы с февраля по июнь, а затем снова возрастает с началом сезона ураганов. Пик заболеваемости в Сан-Хуане обычно приходится на октябрь, ближе к концу сезона ураганов. С другой стороны, в Икитосе мы наблюдаем рост числа случаев с октября по март, что примерно соответствует сезону дождей в регионе.

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

Большинство функций в нашем наборе данных говорят сами за себя. Однако термин «суточный температурный диапазон» может быть незнаком неметеорологам. Эта характеристика описывает диапазон между минимальной и максимальной температурой воздуха, наблюдаемой за 24-часовой период. Более высокий диапазон создает худшие условия для выживания комаров, и поэтому эта особенность отрицательно коррелирует со случаями лихорадки денге.

Чтобы подготовить наши функции к моделированию, мы выполнили несколько шагов предварительной обработки. Во-первых, мы очистили набор данных, объединив повторяющиеся измерения, записанные из разных источников (например, четыре разных показателя осадков), в единые признаки. Мы работали с нулевыми значениями с помощью линейной интерполяции, которая заменяет отсутствующие значения путем проведения прямой линии между точкой перед отсутствующим значением и точкой после отсутствующего значения. Этот подход помогает сохранить основные тенденции в наборе данных. Для масштабирования наших данных мы использовали MinMaxScaler, который преобразует все наши значения, чтобы они попадали в диапазон от 0 до 1. Кроме того, поскольку наши данные основаны на времени, мы внедрили временные задержки, чтобы включить информацию за прошлые недели в наши прогнозы. Мы обнаружили, что использование временных задержек в четыре недели дало лучшую производительность наших моделей, а использование большего количества задержек привело только к переоснащению. Ниже приведена наша полная функция предварительной обработки.

def preprocess(data,train=True):
    
    #Reset index to prevent issues when merging
    reindexed_data = data.reset_index(drop = True)
    
    #Fill nas with interpolation
    features = reindexed_data.interpolate(method='linear')
    
    new_columns = pd.DataFrame()
    #Scale and average total precipitation features
    precipitation_to_avg = features.loc[:,['reanalysis_sat_precip_amt_mm','precipitation_amt_mm','station_precip_mm','reanalysis_precip_amt_kg_per_m2']]
    scaled_precip = pd.DataFrame(MinMaxScaler().fit_transform(precipitation_to_avg), 
                                columns=precipitation_to_avg.columns)
    new_columns.loc[:,'avg_total_precipitation'] = scaled_precip.mean(axis=1)
    
    #Scale and average temperature features
    temps_to_avg = features.loc[:,['reanalysis_air_temp_k','reanalysis_avg_temp_k','station_avg_temp_c']]
    scaled_temps = pd.DataFrame(MinMaxScaler().fit_transform(temps_to_avg), 
                                columns=temps_to_avg.columns)
    new_columns.loc[:,'avg_temp'] = scaled_temps.mean(axis=1)
    
    #Scale and average max temperature features
    max_temps_to_avg = features.loc[:,['station_max_temp_c','reanalysis_max_air_temp_k']]
    scaled_max_temps = pd.DataFrame(MinMaxScaler().fit_transform(max_temps_to_avg), 
                                columns=max_temps_to_avg.columns)
    new_columns.loc[:,'avg_max_temp'] = scaled_max_temps.mean(axis=1)
    
    #Scale and average min temperature features
    min_temps_to_avg = features.loc[:,['station_min_temp_c','reanalysis_min_air_temp_k']]
    scaled_min_temps = pd.DataFrame(MinMaxScaler().fit_transform(min_temps_to_avg), 
                                columns=min_temps_to_avg.columns)
    new_columns.loc[:,'avg_min_temp'] = scaled_min_temps.mean(axis=1)
    
    #Scale and average diurnal temperature features
    diurnal_temps_to_avg = features.loc[:,['station_diur_temp_rng_c','reanalysis_tdtr_k']]
    scaled_diurnal_temps = pd.DataFrame(MinMaxScaler().fit_transform(diurnal_temps_to_avg), 
                                columns=diurnal_temps_to_avg.columns)
    new_columns.loc[:,'avg_diurnal_temp'] = scaled_diurnal_temps.mean(axis=1)
    
    #Scale humidity and dew point features
    remaining_feats_to_scale = features.loc[:,['reanalysis_specific_humidity_g_per_kg','reanalysis_relative_humidity_percent','reanalysis_dew_point_temp_k']]
    scaled_feats = pd.DataFrame(MinMaxScaler().fit_transform(remaining_feats_to_scale), 
                                columns=remaining_feats_to_scale.columns)
    
    df = new_columns.join(scaled_feats)
    
    #Add lagged features for 4 weeks
    to_shift = ['avg_max_temp','reanalysis_relative_humidity_percent',
       'reanalysis_specific_humidity_g_per_kg','avg_total_precipitation','reanalysis_dew_point_temp_k','avg_min_temp'] 
    for i in to_shift:
        df[i+'_1lag'] = df[i].shift(+1)
        df[i+'_2lag'] = df[i].shift(+2)
        df[i+'_3lag'] = df[i].shift(+3)
        df[i+'_4lag'] = df[i].shift(+4)
    df = df.fillna(method='bfill')
    
    #Merge together

    if train:
        non_scaled_feats = features[['ndvi_ne', 'ndvi_nw',
           'ndvi_se', 'ndvi_sw','total_cases']]
    
    else:
        non_scaled_feats = features[['ndvi_ne', 'ndvi_nw',
           'ndvi_se', 'ndvi_sw']]
    
    final_df = df.join(non_scaled_feats)
    
    return final_df

После предварительной обработки наших данных мы были готовы приступить к этапу моделирования.

Моделирование

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

Чтобы установить базовый балл для сравнения, мы использовали среднее значение общего количества еженедельных случаев для каждого города в наших данных в качестве нашего прогнозируемого значения. Эта базовая модель дала MAE 35 случаев для объединенных двух городов. Это послужило эталоном, по которому мы оценивали производительность наших более сложных моделей. Поэкспериментировав с различными алгоритмами, мы обнаружили, что регрессор случайного леса лучше всего работает с нашим набором данных. Ниже приведен пример прогнозов, сгенерированных этой моделью.

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

def peak_boost(preds, scalar, mixed=None):
    '''A function to boost peaks of the predictions'''
    #Take derivative of predictions
    deriv = np.gradient(preds)
    #Scale derivative
    new_deriv = deriv * scalar
    #Integrate the new derivative starting from 0
    integral = cumtrapz(new_deriv, initial=0)
    length = range(len(preds))
    #Add intercept to match initial range
    boosted = [x+preds[0] for x in integral]
    
    #Options to put floors on the output, preventing negative case counts
    if mixed == 'max':
        out = [max(boosted[i], preds[i]) for i in length]
    elif mixed == 'pos':
        out = [max(boosted[i], 0) for i in length]
    else:
        out = boosted
    
    #Doubling values above a 0.99 threshold in order to better estimate peaks
    Ser = pd.Series(out)
    
    quantile_value = Ser.quantile(q=0.99)
    print(f'Values are doubled over this threshold: {quantile_value}')
    
    for k,v in Ser.items():
        if (v > quantile_value):
            Ser[k] = int(Ser[k]*2)
    
    return np.array([int(x) for x in Ser])

Вот как выглядели наши прогнозы случайного леса после применения функции повышения пика.

Ускоренные прогнозы оказались гораздо более точными. С применением функции повышения пика наши модели получили вот такие оценки.

Наша лучшая модель достигла MAE в 23 случая для двух объединенных городов. Это была модель, которую мы решили запустить в производство и использовать для нашего приложения.

Внешний интерфейс

Чтобы прогнозировать случаи лихорадки денге в будущем, нам нужно было использовать API прогнозирования погоды, который включал бы как можно больше функций, на которых была обучена наша модель. Мы нашли бесплатный API с открытым исходным кодом под названием Open-Meteo, который предоставлял данные прогноза за две недели и включал большинство необходимых нам функций.

Используя модуль даты и времени Python, наше приложение ежедневно обновляет данные, полученные из API, и генерирует новые прогнозы. Мы также присваиваем оценку риска нашим прогнозируемым случаям, как показано на снимке экрана ниже.

Оценка риска рассчитывается относительно исторических случаев, наблюдаемых в каждом городе. Если прогнозируемое количество случаев превышает 70-й процентиль прошлых случаев, мы присваиваем классификацию «Высокий». Если он ниже 70-го процентиля, но выше 40-го процентиля, мы присваиваем классификацию «Средний», а если он ниже 40-го процентиля, мы присваиваем ему классификацию «Низкий».

Будущая работа

Существуют ограничения точности наших прогнозов из-за ограниченности наших данных. Например, мы предположили, что динамика передачи лихорадки денге не изменилась с 2010 года, тогда как на самом деле численность населения, городская структура и другие факторы могли изменить способ распространения болезни. Кроме того, наши прогнозы настолько хороши, насколько хорош прогноз погоды, предоставляемый используемым нами API.

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

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

Полную структуру кода и обновления нашего прогресса можно найти в нашем репозитории GitHub.

Ниже запись того, как я представляю этот проект.