Приложение 2 (рекомендуемое). Примеры построения прогнозов (краткосрочных, долгосрочных) эпизоотической активности природных очагов чумы

Приложение 2
(рекомендуемое)

Примеры построения прогнозов (краткосрочных, долгосрочных) эпизоотической активности природных очагов чумы

1. Пример проведения прогноза на основе аппроксимации исходного временного ряда уравнением линейной регрессии

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

Для реализации линейного прогноза в Excel используют следующую последовательность команд: "Сервис", "Анализ данных...", "Регрессия". В результате получают три таблички, содержащие данные об оценке надежности линии аппроксимации, о величине свободного члена и коэффициенте наклона линии регрессии, а также значения их ошибок при определенном уровне значимости. Для проведения прогноза необходимо и достаточно в полученное уравнение подставить очередное значение X.

Рассмотрим пример по изменению количества жилых нор монгольской пищухи (основной носитель, обеспечивающий циркуляцию возбудителя чумы алтайского подвида) на эндемичных территориях в Горно-Алтайском высокогорном природном очаге чумы за 1991-1999 гг. (табл. 5) и сделаем прогноз этого показателя на 2000 г. Для расчета параметров линии регрессии в программе Excel следует использовать приводимый ниже "пошаговый" алгоритм анализа временного ряда.

1. Внести данные из табл. 5 в Excel.

2. Выделить курсором графу, отражающую количество нор пищухи.

3. Активизировать команду "Мастер диаграмм" верхней панели программы.

4. Выбрать команду "график" из появившегося диалогового окна, а в ней активировать пункт "график с маркерами".

Таблица 5

Данные по изменению количества жилых нор монгольской пищухи - основного носителя возбудителя чумы алтайского подвида в Горно-Алтайском высокогорном природном очаге

Год	Число жилых нор монгольской пищухи на 1 га	Число жилых нор на 1 га, рассчитанное по линейному тренду	Квадраты остатков после удаления тренда
1991	5,8	5,6	0,035
1992	6,8	5,9	0,796
1993	6,2	6,2	0,000
1994	5,4	6,5	1,210
1995	6,7	6,8	0,009
1996	6,6	7,1	0,242
1997	6,7	7,4	0,473
1998	8,8	7,7	1,245
1999	8,2	8,0	0,048
2000 г. фактически: 9,1		2000 г. прогноз: 8,3	Сумма квадратов остатков: 4,060

5. Построить график исходного временного ряда.

6. Добавить на график линию тренда типом "линейный", а также помещаем на график уравнение аппроксимации и коэффициент детерминации - (правила выполнения этих процедур описаны выше).

В результате на рисунке временного ряда в программе Excel получим уравнение линейного тренда: у = 0,296х + 5,316 с величиной (рис. 2).

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

Выполнение всех перечисленных операций (исключая ввод информации в программу) и получение готового аналитического уравнения занимает не более 15 секунд. Программа при построении уравнений по графикам всегда производит расчеты исходя из ряда значений, различающихся на единицу: в виде 1, 2, 3,...n или 1997, 1998, 1999 ... и т.д. Следовательно, если данные по X имеют другой шаг различий, то для построения уравнения регрессии необходимо воспользоваться в Excel командой "Регрессия" из пакета "Анализ данных".

Исходя из полученного уравнения, в 2000 г. количество жилых нор монгольской пищухи должно быть равно 8,3 на 1 га, что достаточно близко к значению (9,1), фактически наблюдавшемуся в тот сезон (табл. 5). Отклонение ожидаемого значения от наблюдаемого, через которое можно охарактеризовать точность прогноза, составляет в нашем случае 5,5% фактической величины. Следует учесть, что вероятность точной реализации точечного прогноза крайне мала. Прогноз должен сопровождаться оценкой величины доверительного интервала:

, где

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

- ошибка прогноза.

В свою очередь ошибку прогноза определяют по формуле:

, где

- колеблемость ряда относительно тренда;

- поправочный коэффициент, который зависит от типа тренда и периода упреждения (предсказания).

Ошибка прогноза для уравнения линейной регрессии будет иметь следующий вид:

, где

- фактические значения наблюдений;

- расчетные значения уровней ряда по линейному тренду;

n - длина временного ряда;

m - число параметров в уравнении тренда без свободного члена;

l - период упреждения. Причем время должно быть представлено рядом натуральных чисел от 1 до n.

В нашем примере (табл. 5) при составлении прогноза на 2000 г. l = 1 и m = 1. Следовательно:

Прежде чем рассчитать доверительный интервал, необходимо решить, какой уровень ошибки прогноза мы бы считали приемлемым. Допустимо применение 10%-й ошибки. Исходя из величины ошибки 10% (т.е. р < 0,1) и числа степеней свободы исследуемого временного ряда (df = n - 2), равного 7, величина критерия Стьюдента (t) составляет 1,90. Следовательно, доверительный интервал точечного прогноза равен 1,64.

Таким образом, в 2000 г. количество жилых колоний монгольской пищухи должно с вероятностью 90% находиться внутри интервала , т.е. составлять от 6,7 до 9,9 нор на 1 га.

При проведении факторного прогноза необходимо придерживаться нескольких эмпирических правил.

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

2. В статистическую модель должно входить как можно меньше предикторов. Модель должна быть наиболее простой из возможных. На каждый предиктор должно приходиться не менее 10-15 наблюдений зависимой переменной.

3. Рекомендуется не подставлять в уравнения прогноза значения предикторов, отличающиеся больше чем на 1/3 от минимальной и максимальной величины значений, наблюдавшихся в исходном временном ряде.

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

2. Пример проведения прогноза на основе гармонического анализа

Рассмотрим особенности применения метода гармонического анализа для прогнозирования временного ряда с регулярными циклами (Фурье - аппроксимация).

С помощью ряда Фурье динамика явления представляется функцией времени, в которой слагаемые расположены по убыванию их периодов:

Величина k определяет гармонику ряда Фурье и берется как целое число, начиная с 1. Часто для аппроксимации временного ряда достаточно рассмотрения 3-4 гармоник. Параметры уравнения рассчитывают методом наименьших квадратов:

- то есть простая средняя арифметическая ряда;

и

и

и и т.д.,

или в общем виде:

и

Первому наблюдению (t) обычно присваивают значение 1 (или 0). К каждому следующему прибавляют величину , где n - длина ряда.

Ряд Фурье с одной гармоникой записывают как:

Ряд Фурье с двумя гармониками:

Ряд Фурье с тремя гармониками:

и так далее.

Для прогноза (экстраполяции) в уравнение с выбранным числом гармоник подставляют значение времени (t) необходимого порядка.

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

В качестве примера проведем гармонический анализ временного ряда, характеризующего площади буферной зоны эпизоотии , т.е. площади, где надо проводить профилактические (противоэпидемические) мероприятия, на территории Прикаспийского песчаного природного очага чумы за 1979-1990 гг. (табл. 6).

Отсчет t начнем с 1, прибавляя каждый раз к его значению величину , то есть в нашем случае 0,523. Очень важно, чтобы все расчеты были проведены в "радианах".

Таблица 6

Площади буферной зоны эпизоотий на территории Прикаспийского песчаного природного очага чумы

Год	Площадь буферной зоны эпизоотии,	Год	Площадь буферной зоны эпизоотии,
1979	1 548	1991	2611
1980	78	1992	1 004
1981	741	1993	894
1982	955	1994	1 112
1983	2 689	1995	1 971
1984	1 941	1996	495
1985	1 827	1997	0
1986	2 267	1998	331
1987	2 653	1999	313
1988	2 604	2000	258
1989	1 979	2001	1 187
1990	2 844	2002	232

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

Таблица 7

Параметры ряда Фурье, необходимые для составления уравнения прогноза с использованием 4 гармоник

Год	N п/п	У	t
1979	1	1 548	1	836,4	1 302,6	-644,2	1 407,6	-1 532,5	218,5	-1 011,8	-1 171,5
1980	2	78	1,5	3,7	77,9	-77,6	7,4	-11,1	-77,2	76,6	-14,7
1981	3	741	2,0	-339,5	658,7	-430,0	-603,5	733,4	-105,7	-242,0	700,4
1982	4	955	2,6	-803,2	516,6	396,0	-869,0	137,0	945,1	-626,5	-720,8
1983	5	2 689	3,1	-2 685,9	129,7	2 676,5	-259,1	-2 660,9	387,9	2 639,1	-515,9
1984	6	1 941	3,6	-1 726,1	-887,8	1 128,8	1 579,0	-281,6	-1 920,5	-628,0	1 836,6
1985	7	1 827	4,1	-989,6	-1 535,8	-755,0	1 663,7	1 807,5	-266,5	-1 203,0	-1 375,0
1986	8	2 267	4,7	-111,2	-2 264,3	-2 256,1	222,1	332,4	2 242,5	2 223,5	-442,0
1987	9	2 653	5,2	1 211,6	-2 360,2	-1 546,3	-2 155,8	-2 624,0	391,1	-850,4	2 513,0
1988	10	2 604	5,7	2 187,8	-1 412,2	1 072,3	-2 373,0	-385,9	-2 575,2	-1 720,8	-1 954,4
1989	11	1 979	6,2	1 976,5	-98,6	1 969,2	-197,0	1 956,9	-294,9	1 939,8	-392,0
1990	12	2 844	6,8	2 531,1	1 296,8	1 661,3	2 308,3	426,0	2 811,9	-903,1	2 696,8
Сумма		22 126	46,54	2 091,9	-4 576,5	3 194,9	730,7	-2 102,7	1 757,1	-306,8	1 160,5

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

; ; ; ; ;

; ; ;

Окончательные результаты гармонического анализа ряда, включающие 4 гармоники, приведены в табл. 8.

Таблица 8

Расчет четырех периодических составляющих и уравнений, характеризующих площади буферной зоны эпизоотий на территории Прикаспийского песчаного природного очага чумы

Число гармоник	Уравнения динамики эпизоотической активности, включающие от одной до четырех гармоник
1
2
3
4

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

Из табл. 9 следует, что уже уравнение с двумя гармониками неплохо описывает исходный ряд. Уравнение с 4 гармониками объясняет 89,8% колебаний уровней исследуемого временного ряда. Вместе с тем и уравнение, содержащее три периодических составляющих, вполне хорошо аппроксимирует исходную последовательность . Поэтому для прогноза необходимо использовать более простой ряд Фурье с 3 гармониками.

Проводится расчет ожидаемой площади буферной зоны эпизоотии на территории Прикаспийского песчаного природного очага чумы на период с 1991 по 1995 год (с n = 13 по n = 17) и сравнение ее с фактической. Для этого в соответствующее уравнение табл. 9 подставим искомые значения параметра времени: 13, 14, ... 17, преобразованные в t как это сделано в табл. 7.

Таблица 9

Сравнение исходного временного ряда буферной зоны эпизоотии на территории Прикаспийского песчаного природного очага чумы с расчетными уровнями, полученными с использованием различного числа гармоник

Год	Исходные уровни ряда	Расчетные значения уровней ряда, полученные на основе уравнений, содержащих перечисленное число гармоник
		1	2	3	4
1979	1 548	1 390,4	1 279,5	1 667,8	1 554,8
1980	78	1 098,5	579,9	339,8	253,1
1981	741	1 006,1	597,9	209,3	408,8
1982	955	1 138,0	1 248,0	1 487,5	1 375,1
1983	2 689	1 458,8	1 977,1	2 366,1	2 278,8
1984	1 941	1 882,7	2 291,4	2 052,5	2 252,1
1985	1 827	2 296,2	2 187,0	1 797,6	1 685,7
1986	2 267	2 588,6	2 070,6	2 308,9	2 221,0
1987	2 653	2 681,6	2 272,3	2 662,1	2 861,7
1988	2 604	2 550,4	2 658,7	2 421,0	2 309,7
1989	1 979	2 230,1	2 747,8	2 357,6	2 269,2
1990	2 844	1 806,3	2 216,2	2 453,3	2 652,9
Корреляция исходного ряда с расчетными (r)		0,71	0,85	0,93	0,95
Коэффициент детерминации , %		50,4	71,8	86,9	89,8

Результаты проведенных расчетов приведены в табл. 10.

Таблица 10

Прогноз буферной зоны эпизоотии на территории Прикаспийского песчаного природного очага чумы на 1991-1994 гг. на основе ряда Фурье с тремя гармониками

Год	Фактическая площадь эпизоотий (по табл. 6)	Период времени прогноза (с 1991 г., т.е. t = 6,8 + 0,523 для каждого следующего года по табл. 7)	расчетные значения по ряду Фурье с тремя гармониками (табл. 9)
1991	2 611	7,3	1 676,6
1992	1 004	7,8	347,2
1993	894	8,3	203,1
1994	1 112	8,8	1 475,3
Корреляция Пирсона уровней исходного ряда и расчетного (r)			0,73
Коэффициент детерминации для оценки точности единичных прогнозов , %			53,3

Полученный прогноз на четыре года, если обсуждать оценку тенденций изменения эпизоотической активности очага, достаточно надежный, о чем свидетельствует высокое значение (0,90) коэффициента корреляции Спирмэна. Однако отклонения численных значений ожидаемой и фактической площади буферной зоны эпизоотии на территории Прикаспийского песчаного природного очага чумы остаются высокими (величина отклонений по четырем годам колеблется от 33 до 77%).

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

3. Пример анализа и прогноза эпизоотической активности Горно-Алтайского высокогорного природного очага чумы

Каждый природный очаг чумы обладает особенностями биоценотической и пространственной структуры, характера и интенсивности эпизоотических проявлений, динамики эпизоотической активности в пространстве и времени. Поэтому при составлении эпизоотологических прогнозов для определенного очага необходимо учитывать его характерные черты. Для каждого природного очага чумы могут быть определены конкретные прогностические признаки, которые позволяют прогнозировать изменения его эпизоотической активности. На основании многолетнего изучения Горно-Алтайского высокогорного природного очага чумы установлено, что между ежегодным уровнем эпизоотической активности и уровнем численности монгольской пищухи - основного носителя возбудителя чумы алтайского подвида - наблюдается статистически существенная связь. Выявленные закономерности в изменениях численности населения монгольской пищухи и эпизоотической активности очага, обусловленной циркуляцией возбудителя чумы алтайского подвида, позволяют статистически обоснованно давать объективный краткосрочный эпизоотологический прогноз в очаге.

Предлагаемый алгоритм направлен на прогноз эпизоотической активности Горно-Алтайского высокогорного природного очага при циркуляции чумного микроба алтайского подвида. Прогнозирование интенсивности эпизоотических проявлений с участием возбудителя чумы основного подвида (выявлена его циркуляция на ряде территорий с 2012 г.) должно основываться на использовании других прогностических признаков.

Основой эпизоотологического прогноза являются многолетние показатели уровней численности монгольской пищухи и эпизоотической активности, которые представляются в виде соответствующих временных рядов. Эти характеристики определяются ежегодно при проведении эпизоотологического обследования очага. После получения показателей в текущем году их добавляют к имеющимся временным рядам в программе Microsoft Office Excel и отображают графически.

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

Динамика численности населения монгольской пищухи. Численность монгольской пищухи определяют в весенний (апрель-июнь) и осенний (август-октябрь) периоды по количеству жилых нор-колоний на 1 га. Для этого на каждом участке эпизоотологического обследования на 2-4-километровых маршрутах подсчитывают все норы-колонии в полосе 10-30 м, после чего пересчитывают их число на 1 га. Для оценки уровня численности населения зверька в очаге все учеты, проведенные в определенный сезон текущего года, усредняют. Поскольку между временными рядами весенней и осенней численности наблюдается тесная статистически значимая связь (коэффициент корреляции Пирсона около 0,7), то для прогноза используют среднегодовую численность населения монгольской пищухи на территории очага, которая рассчитывается как средняя арифметическая от средних значений весенних и осенних показателей. Определение показателя численности всего населения за год основывается на оценках нескольких независимых выборочных совокупностей. Используемый подход позволяет достаточно объективно охарактеризовать состояние населения по всей площади очага. Изменения численности монгольской пищухи в отдельных популяциях, расположенных в очаге, характеризуются определенными различиями, однако в целом циклы в разных популяциях синхронны, то есть пики и депрессии численности в них чаще всего совпадают, что дает основание рассматривать общую тенденцию в динамике численности зверька на территории очага.

Динамика среднегодовой численности населения монгольской пищухи с 1971 года по 2015 год представлена на рис. 3. Этот процесс имеет комбинированный характер.

Он складывается из направленного долговременного роста и определенных колебаний численности. Оценка центральной тенденции изменения численности с помощью регрессионного анализа показывает, что коэффициент регрессии высоко значим (, df = 43, t = 4,92, р < 0,001). Исходя из уравнения линейной регрессии среднегодовая численность в 1971 г. составляла 4,0, а в 2015 г. - 7,7 жилых нор на 1 га. Для определения закономерности процесса колебаний численности была проведена описанная выше процедура приведения временного ряда к стационарному виду. Она проводилась путем вычитания теоретических значений, полученных из уравнения линейной регрессии, от эмпирических, полученных при непосредственной оценке ежегодных показателей. Выявлено, что колебания численности населения монгольской пищухи на территории очага имеет неслучайный характер. Коэффициент автокорреляции стационарного временного ряда составил 0,408 (n = 44, р < 0,01), что свидетельствует о наличии низкочастотных циклических колебаний численности (иными словами - с большим периодом). Непосредственное обращение к данным, описывающим динамику численности (рис. 3), показывает, что период колебаний находится в пределах семи-девяти лет. Статистическая оценка величины периода циклических колебаний может быть проведена на основе построения коррелограммы. Построенная коррелограмма показывает, что период циклических колебаний численности населения монгольской пищухи в Горно-Алтайском природном очаге чумы составляет семь-восемь лет (рис. 4).

Динамика эпизоотической активности. Начиная со времени обнаружения очага в 1961 г., возбудитель чумы выделяют в нем ежегодно. При этом эпизоотическая активность очага по годам значительно варьирует.

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

Динамика показателей эпизоотической активности с 1971 года по 2015 год приведена на рис. 5. Линейный тренд временных рядов площади эпизоотий и количества эпизоотических секторов высокозначим (соответственно ; df = 43; t = 5,43; р < 0,001 и ; df = 43; t = 5,84; р < 0,001), что свидетельствует о возрастании эпизоотической активности очага во времени.

Поскольку во временных рядах площади эпизоотий и количества эпизоотических секторов наблюдается существенный тренд на увеличение, для анализа их циклической составляющей была проведена процедура преобразования числовых последовательностей в стационарные временные ряды, как это было описано выше. Коэффициент автокорреляции стационарного временного ряда первого показателя равен 0,587, второго - 0,674 (n = 44; р < 0,01). Это свидетельствует о неслучайном характере колебаний эпизоотической активности очага. Однако циклические колебания с периодом, характерным для численности населения основного носителя, с помощью построения коррелограмм не обнаруживаются.

Оценка связи между показателями эпизоотической активности и численности населения монгольской пищухи. Полученные временные ряды численности и эпизоотической активности отображаются на одном графике (рис. 5), проводится их визуальное сравнение и статистическая оценка уровня связи между ними. При этом хорошо видно, что центральные тенденции временных рядов показателей, отражающих эпизоотическую активность очага и численность населения основного носителя, совпадают. Коэффициент корреляции Пирсона между эмпирическими временными рядами численности и площадей эпизоотий составил 0,685 (df = 43, р < 0,001),