Обоснование статистической модели для описания и прогноза рисков населения подвергнуться атакам иксодовых клещей

Клещевой вирусный энцефалит (КВЭ) является эндемичным заболеванием в Центральной, Восточной и некоторых частях Северной Европы, а также в Северной и Центральной Азии, где ежегодно регистрируется 10–12 тыс. случаев [1].

Уральский федеральный округ (УФО) является высокоэндемичным по клещевому энцефалиту, при этом наибольшее количество пострадавших от укусов клещей и заболеваемость приходятся на Свердловскую и Курганскую области [2, 3]. Регистрируе-

Мищенко Владимир Алексеевич – научный сотрудник лаборатории трансмиссивных вирусных инфекций и клещевого энцефалита, инженер-исследователь лаборатории эволюционной экологии (e-mail: , ; тел.: 8 (343) 261-99-47, 8 (965) 515-34-89; ORCID: .

Кшнясев Иван Александрович – кандидат биологических наук, старший научный сотрудник лаборатории популяционной экологии и моделирования (e-mail: ; тел.: 8 (343) 210-38-58; ORCID: .

Вялых Иван Владимирович – кандидат ветеринарных наук, ведущий научный сотрудник, заведующий лабораторией трансмиссивных вирусных инфекций и клещевого энцефалита (e-mail: ; тел.: 8 (343) 261-99-47; ORCID: .

мая заболеваемость КВЭ коррелирует с обращаемостью в связи с укусами клещей. Также для КВЭ характерны долгосрочные периодические изменения, наблюдаемые на протяжении многих десятилетий и связанные с такими факторами, как демография человеческой популяции, изменения в землепользовании и соответствующей плотности организмов (популяций) в условиях дикой природы, рекреационная деятельность людей. Необходимо отметить важное влияние изменений в климатической системе как возможной движущей силе для цикличности процессов, связанных с КВЭ [4].

В Западной Европе вирус клещевого энцефалита (ВКЭ) передается главным образом клещами Ixodes ricinus , тогда как переносчиком для сибирских и дальневосточных подтипов вируса является I. persulcatus [5]. Для иксодид характерно наличие нескольких стадий развития в жизненном цикле. Каждая стадия жизненного цикла клеща длится от нескольких месяцев до одного года, поэтому весь цикл обычно завершается за два или три года, хотя, в зависимости от географического положения и актов питания, может варьироваться от двух до шести лет. Стадийность развития иксодид непосредственно связана со сменой их основных хозяев – прокормителей для личинок и нимф, а также и (пусть, возможно, и в меньшей степени) для имаго – мелких млекопитающих (ММ) [6]. Выраженные резкие изменения плотности популяций мелких млекопитающих с амплитудой 2–3-го порядка и квазипериодом 3–5 лет, известны как пример так называемых популяционных циклов [7, 8]. Предполагается, что обилие мелких млекопитающих положительно сказывается (с некоторыми лагами) на численности генераций последовательных стадий (в основном личинок и нимф) иксодовых клещей (ИК). Через контур численности голодных и свободных ИК – не нашедших обычного прокормителя – существует прямая связь с рисками для населения пострадать от укусов клещей [9].

Цель исследования – анализ и синтез адекватной формализованной / параметризованной статистической модели для описания и прогноза рисков населения подвергнуться атакам иксодовых клещей.

Материалы и методы. Исходные данные (количество пострадавших от укусов клещей) в Свердловской области получены из формы федерального статистического наблюдения (ФФСН) № 2 «Сведения об инфекционных и паразитарных заболеваниях» (раздел 1) (за 2007–2021 гг.) и из материалов ФБУЗ «Центр гигиены и эпидемиологии в Свердловской области» (за 1992–2006 гг.).

Учет ММ проводили на стационарных участках методом ловушко-линий с помощью ловушек-плашек и деревянных трапиковых живоловок [10] в весенний и осенний периоды с 1992 по 2021 г.

С целью количественной оценки влияния плотности ММ на обращаемость пострадавших в связи с укусами клещей рассчитывали шансы атак. Шансы ( odds ) – отношение числа пострадавших от укусов клещей ( N 1 ) к населению Свердловской области, за исключением пострадавших ( N 0 ).

От шансов возможно перейти к интуитивно понятной шкале вероятностей или рисков (1):

exp ( LogOdds ) _ Odds

P ( X" )             /             \                 "

1 + exp ( LogOdds ) 1 + Odds

Использовали стандартный аппарат теории обобщенных линейных моделей (GLM)1 – логит-регрессию (2):

{                \

Ln I -2l I = b0 + £biXi.

Оценивали эффекты следующих (экзогенных) предикторов ( X i ): плотность ММ: весной и осенью предыдущего года (SMspr ( t –1) и SMfall ( t –1)) и весной текущего года (SMspr ( t )). А для учёта явной нестационарности (тренда) временного ряда зависимой переменной – логарифма шансов атак ИК (Log odds ticks attacks – LOTA) использовали несколько (различной степени строгости или, наоборот, эвристичности) подходов:

• 1)    включение в модель члена авторегрессии первого порядка: AR (1) X – где дополнительным предиктором служило наблюдаемое значение LOTA ( t –1) в предыдущем году;

• 2)    предварительное сглаживание (и последующее включение в модель как дополнительного предиктора) временного ряда LOTA с помощью локальной регрессии ( loess ); оптимальное значение параметра сглаживания ( span = 0,427) установлено с помощью 10-кратной перекрестной проверки (crossvalidation);

• 3)    предварительная STL (Seasonal and Trend decomposition using Loess) – декомпозиция временного ряда LOTA, и выделение тренда с помощью сглаживания методом локальной регрессии.

Исходный временной ряд и варианты сглаживания представлены на рис. 1.

Отношения шансов ( OR ) и их доверительные интервалы (95 % ДИ) приведены после преобразования: OR = exp ( b_i ) или OR ^-1 = 1/exp ( b_i ), где b_i – параметры логит-регрессии (логарифмы отношения шансов). Отношения шансов для редких событий (частота < 10^-1) приближаются к значениям относительного риска, что упрощает интерпретацию.

Для сравнения и ранжирования моделей логит-регрессии использовали информационный критерий Акаике ( Akaike information criterion – AIC ), определяющий оптимальность как компромисс между

Год

Рис. 1. Динамика логарифма шансов атак иксодовых клещей и варианты сглаживания временного ряда для выделения тренда. Свердловская область, 1992–2021 гг.

точностью и сложностью модели. Наименьшей величине AIC соответствует и статистически более адекватная модель. Сравнение моделей выполнено на основе модификации исходного AIC – состоятельного критерия Акаике (CAIC), рассчитанного по формуле (3):

CAIC = - 2 LL + k [ 1 + In ( m ) ] , (3)

где LL – логарифм максимума функции правдоподобия, k – число параметров, m – число наблюдений. В данной модификации, в сравнении с AIC , назначает более жесткий «штраф» за дополнительные параметры².

«Вес» (относительное правдоподобие) каждой модели, рассчитанный по формуле (4), использовали для ранжирования и сравнения конкурирующих моделей:

exp ( - 0,5 A CAIC ) w = Z exp ( - 0,5 A CAIC ) ’

Представленный «вес» может быть интерпретирован как апостериорная вероятность того, что -я модель является лучшей при исследованном множестве прочих моделей-претендентов. Если «вес» отличается менее чем на 10 % от w max , считали, что эти модели идентичны по качеству наилучшей [11].

Статистическая обработка результатов и их визуализация проведены с использованием пакета прикладных программ Statistica v. 10.0 (StatSoft, Inс) и статистической среды R (v. 3.4.4) [12].

Результаты и их обсуждение. Наблюдаемая многолетняя динамика шансов для населения пострадать от укусов ИК представляет собой очевидно нестационарный (ни по среднему, ни по амплитуде, ни по периоду) временной ряд (рис. 1), в котором можно выделить:

• –    тренд – долговременный дрейф средних значений, относительно медленные изменения во времени;

• –    квазипериодическую относительно быструю повторяющуюся компоненту;

• –    остатки – нерегулярный компонент ряда, относительно высокочастотный шум.

Тренд, вероятно, является «эхом»: 1) как длинного циркадианного погодного-климатического цикла, связанного с крупномасштабными феноменами (формированием и продолжительностью / устойчивостью форм атмосферной циркуляции и таких феноменов, как антициклоны, блокирующие западный перенос), вызывающими характерные многолетние колебания температуры воздуха и количества осадков в регионе [13, 14]; 2) так и особенностей социально-экономической ситуации в 90-е гг. ХХ в. в Российской Федерации и последующим её улучшением [15].

Квазитрехлетние циклы ( T = 2–3 года) могут быть связаны и с параметрами жизненного цикла ИК, и с динамикой численности ММ [16], для которых характерны такие периоды [17]. Популяции ММ и ИК связаны – годы с высокой активностью клещей (оцениваемой здесь по числу пострадавших) следуют за годами с максимальной плотностью грызунов [18]. Личинки и нимфы питаются в основном на мелких млекопитающих, что приводит к естественной передаче ВКЭ.

Динамика логарифма шансов атак иксодовых клещей оптимально ( w ≈ 1, CAIC = 13005588) описывается моделью с трендом, полученным путем STL-декомпозиции. Модели авторегрессии и локального сглаживания значительно «проигрывают» вышеуказанной: w ~ 0, CAIC = 13009459 и w ~ 0, CAIC = 13072412 соответственно.

В оптимальную модель включено пять параметров: тренд, выделенный из временного ряда LOTA с помощью STL-декомпозиции (trend-STL); плотность ММ: весной текущего, весной предыдущего и осенью предыдущего годов; свободный член b₀ (таблица).

Включение в список предикторов модели ло-гит-регрессии дополнительной компоненты – оценки тренда позволяет учесть, пусть и эвристически, нестационарность ряда и оценить интересующие нас эффекты (рис. 2, a ).

Рост плотности ММ весной предыдущего года приведет к возрастанию ( p < 0,01) шансов атак ИК (относительно населения Свердловской области) в 1,04 раза (больше на 3,53 %, рис. 2, б ). Вероятно, за счет обилия прокормителей весной прошлого года в текущем году ожидается увеличение численности иксодид различных стадий (нимфы, имаго), вышедших из диапаузы.

Оценки параметров оптимальной модели для описания многолетней динамики логарифма шансов атак иксодовых клещей. Свердловская область, 1992–2021 гг. («лучшая» модель логит-регрессии: LR (4) = 99827; p < 0,0001; R ² = 0,66)

Параметр	b	SE ( b )	Z -Вальда	p -значение	Отношение шансов
					OR	95 % ДИ
b 0	1,10	0,02	47,72	< 0,0001	–	–	–
trend-STL	1,24	0,004	272,10	< 0,0001	3,45	3,42	3,48
SMspr (t–1)	0,03	0,001	41,04	< 0,0001	1,04	1,03	1,04
SMspr (t)	– 0,06	0,001	– 80,81	< 0,0001	1,06-1	1,057-1	1,061-1
SMfall (t–1)	0,02	0,001	20,83	< 0,0001	1,02	1,017	1,021

Примечание: b₀ – свободный член; OR – odds ratio – отношение шансов ; LR (df) – тест отношения правдоподобий с количеством степеней свободы; trend-STL – тренд, выделенный из временного ряда LOTA с помощью STL-декомпозиции; плотность ММ – SM: spr (t–1) – весной прошлого года, spr (t) – весной текущего года, fall (t–1) – осенью прошлого года; -1 – OR ^-1 = 1/exp ( b_i ).

Рис. 2. Зависимость вероятности атак иксодовых клещей от плотности мелких млекопитающих: a – зависимость от логарифма шансов атак иксодовых клещей; зависимость от плотности мелких млекопитающих: б – весной прошлого года, в – весной текущего года; г – осенью прошлого года

Увеличение плотности ММ весной текущего года приведет к статистически значимому снижению шансов атак ИК в 1,06 раза (меньше на 5,87 %, рис. 2, в ). Таким образом, популяции ММ в фазе роста, являются «перехватчиками» нимф и имаго ИК – достаточное количество прокормителей способствует меньшей вероятности атак на крупных млекопитающих и людей.

Увеличение плотности ММ осенью предшествующего года приведет к статистически значимому росту шансов атак ИК в 1,02 раза (больше на 1,98 %, рис. 2, г ). Состояние популяций ММ в осенний сезон является проекцией весенне-летней плотности

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

На рис. 3 показаны наблюдаемые и ожидаемые (предсказанные) моделью риски (вероятности) атак ИК для населения Свердловской области.

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

Рис. 3. Наблюдаемые и ожидаемые значения вероятности атак иксодовых клещей для населения Свердловской области, 1992–2021 гг.

корреляции Спирмена – r S = 0,70; p < 0,0001). Коэффициент детерминации достаточно высок: R ² = 0,66, следовательно, модель в состоянии объяснить 66 % дисперсии динамики логарифма риска (или вероятности) атак ИК. Ожидается, что включение в модель дополнительных предикторов, связанных с динамикой популяций ИК и прокормителей, климатическими изменениями и социально-экономическими факторами, может улучшить предсказательную способность модели.

Выводы:

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

• 2.    С помощью модели логит-регрессии показано, что численность мелких млекопитающих в предшествующем году и начале текущего эпидсезона активности ИК может служить упреждающим предиктором риска для населения пострадать от укусов иксодовых клещей: высокая плотность мелких млекопитающих в весенний период данного клещевого сезона приводит к выявленному эффекту «перехвата» и количественно оцененному снижению вероятности населения пострадать от укусов клещей. Однако если плотность мелких млекопитающих была высокой в весенний и осенний сезоны года, предшествующего текущему клещевому сезону, то следует ожидать увеличения количества перелинявших и перезимовавших иксодид различных стадий (нимфы, имаго) и роста риска для населения подвергнуться атакам клещей.

• 3.    Прогнозные значения полученной статистической модели адекватно описывают исходный временной ряд шансов / вероятностей атак иксодовых клещей. Включение дополнительных параметров возможно улучшит прогнозные свойства модели.

Финансирование. Исследование проводилось в рамках НИР «Комплексный анализ популяционного иммунитета к вирусу клещевого энцефалита в субъектах Уральского федерального округа» (№ регистрации: 121040500096-4).