Аппарат теории цепей Маркова в профэпидемиологии. Г. В. Федорович (№1, 2013)

Г.В.Федорович, ООО «НТМ-Защита», Москва

Введение.

В предыдущих статьях [1 – 3] было показано, что анализ заболеваемости с временной утратой трудоспособности (далее – ЗВУТ) дает адекватную характеристику условий на рабочих местах. В частности, по уровню ЗВУТ можно оценить классы условий труда в различных подразделениях предприятия [2], а по нозологической структуре ЗВУТ (по данным о количестве лиц с ЗВУТ различных нозологических форм) можно оценить этиологию, т.е. виды и уровни наиболее существенных вредных производственных факторов (далее – ВПФ), действующих на работников [3]. Одновременно был выявлен ряд проблем в области эпидемиологического описания ситуации с профессионально обусловленными заболеваниями (далее – ПОЗ) на производстве.

Дело в том, что документы [4] и [5] упорядочивают анализ ЗВУТ и профессиональных заболеваний (далее – ПЗ) на индивидуальном уровне. По определению, полицевой метод, рекомендованный в этих документах, сводится к изучению только болевших индивидуумов. Например, в [4] рекомендуется собирать данные о количестве К случаев ЗВУТ и об их суммарной длительности D в обследуемом рабочем коллективе. Не определены, однако, аналогичные показатели для работников, оставшихся здоровыми, несмотря на воздействие ВПФ. Поэтому, в частности, имеется значительный произвол в формировании 2х2 таблиц сопряженности, являющихся (см. напр. [3]) эффективным инструментом для оценки степени (условной вероятности) влияния ВПФ на уровень хронических ПОЗ.

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

Общие требования, которым должны удовлетворять такие модели, были определены в работах [6 - 8]. Эти требования формулируются, если рассматривать отдельный трудовой коллектив как единичную реализацию неограниченного количества (ансамбля) идентичных коллективов. При этом внутренние (микроскопические) состояния системы могут быть любыми, совместимыми с заданными значениями внешних параметров, определяющих её макроскопическое состояние. Из всех возможных микроскопических состояний, с максимальной вероятностью реализуются те, которые имеют наибольший статистический вес. Такой подход позволил выявить форму, адекватного статистического описания эпидемиологической ситуации с ПОЗ в трудовых коллективах. Тем самым определились требования, ограничивающие структуру феноменологических моделей, которые можно использовать для создания методики формализованного статистического анализа эпидемиологических данных по ПОЗ. Одна из таких моделей рассмотрена ниже.

1.Цепи Маркова как основа анализа результатов исследований ЗВУТ

 

1.1. Наиболее характерным типом эпидемиологических исследований профессиональных заболеваний являются когортные исследования. Они могут быть как ретроспективными, так и перспективными, их результаты могут быть охарактеризованы как «продольные», т.е. они либо относятся к определенному периоду наблюдений (например, число К случаев ЗВУТ за год), либо сами описывают продолжительность наблюдаемых явлений (например, суммарная длительность D всех ЗВУТ в коллективе). Эти результаты, однако, можно использовать и для расчета «поперечных» характеристик ЗВУТ, например – среднего числа n болеющих в обследуемом коллективе. Как показано в [6], величина n связана с D формулой n = D/Y, где Y – длительность наблюдений. Кроме того, по К и D можно определить среднюю длительность болезни l = D/K , а также средний период циклов «заболевание-выздоровление» L = N*Y/K , где N = численность обследуемого коллектива. Средняя продолжительность периодов от выздоровления до следующего заболевания равна L - l . Непосредственно видно, что n / N = l / L . Обычно результаты приводятся к рекомендованной в [4] длительности наблюдений 1 год (т.е. Y = 365 дней) и численности коллектива N = 100 работников.

Величины, обратные средним длительностям a = 1 / (L-l) и β = 1 / l допускают наглядную интерпретацию как вероятность заболеть (будучи здоровым) и вероятность выздороветь (будучи больным). Обе вероятности относятся к промежутку времени Δt = 1 сутки. Величины 1 - a и 1 – β определяют вероятности остаться здоровым или больным соответственно. Одним из важных следствий такой интерпретации является то, что вероятность a определяется, в значительной степени, условиями на производстве, в то время как вероятность β, в большей степени отражает условия реабилитации – бытовые, медицинского обслуживания и т.п.

Очевидно, справедливо следующее утверждение: если работник сегодня с вероятностью р₀ здоров или с вероятностью р₁ болен (разумеется, р₀ + р₁ = 1), то завтра он будет здоров с вероятностью (1 - a)*р₀ + β*р₁ , а с вероятностью a*р₀ + (1 - β)*р₁ заболеет. Повторяя эти рассуждения, несложно рассчитать вероятности быть больным или здоровым на послезавтра и на последующие дни. Описанная последовательность представляет собой простейший вариант процесса, известного в статистике как цепь Маркова (см. напр. [9]). Развит эффективный аппарат математического анализа таких последовательностей. Им можно воспользоваться для анализа данных о ЗВУТ на различных производствах. Продемонстрируем это.

1.2. Анализ цепи Маркова позволяет оценить распределение работников по срокам ЗВУТ. Рассмотрим отдельного работника. Состояние «здоров» будем отмечать как состояние 0, соответственно «болен» - как 1. Отсчет времени начнем с момента его заболевания, т.е. считаем, что в момент t = 0 произошел переход из состояния 0 в состояние 1. Рассмотрим случайную величинуT, обозначающую время первого возвращения в состояние 0. По определению,T = i Δt , если в моменты времениt₁, t₂, … t_i-1 система остается в состоянии 1. Каждое такое событие происходит с вероятностью 1 - β , так что вероятность остаться больным в течение времени T = i Δt равна (1 – β)^i-1. В момент времениt_iпроисходит возвращение в состояние 0. Вероятность последнего события β . Результирующий закон распределения f_i случайной величиныTимеет вид:

f_i = β(1 – β)^i-1 , (i = 1,2, …) (1)

Нетрудно убедиться, что функция f_i нормирована на единицу:

S_i f_i = b / [1 – (1 - b)] = 1 (2)

Если начальное значение i принять равным нулю, то можно записать соотношение (1) в несколько ином виде:

f_i = b*exp[i*ln(1-b)] » b*exp(-i*b) » [1 – exp(-β)]*exp(-iβ) (3)

В последнем соотношении использовался факт реальной малости величины β , что позволяет записать ln(1 - β) ≈ - β и b » [1 – exp(-β)] .

Соотношение (3) можно ассоциировать с распределением длительностей ЗВУТ, полученным ранее в работе [6], при анализе ансамблевых средних в трудовых коллективах. Как отмечалось в этой работе, статистика ансамблей дает адекватную форму описания систем, состоящих из множества подобных элементов. То обстоятельство, что анализ ситуации с ЗВУТ в трудовых коллективах с использованием методов теории цепей Маркова дает ту же форму представления результатов, что и статистика ансамблей, свидетельствует об адекватности такого подхода.

1.3. Рассмотрим эволюцию начального распределения больных и здоровых работников. Представим его в виде вектора-столбца:

Вероятности переходов π₀₀ = 1 - a , π₀₁ = β , π₁₀ = a и π₁₁ = 1 - β образуют матрицу переходов π = π_jk (j,k = 0,1) . Очевидно, что случай a + β = 0 (т.е. a = 0 и β = 0) не представляет интереса, т.к. в такой системе ничего не меняется: больные остаются больными, здоровые – здоровыми. Аналогично, неинтересен и случай a + β = 2 (т.е. a = 1 и β = 1), т.к. в такой системе на каждом шаге состояние неслучайно меняется на прямо противоположное. В дальнейшем будем полагать, что |1 - a - β| < 1 .

Вектор вероятностей p в момент времени t_i = i Δt определяется матричным уравнением

p(t_i) = π p(t_i-1) (5)

Для дальнейшего представляют интерес несколько результатов теории цепей Маркова:

(1)Последовательное применение операции (5) к начальному состоянию р(t=0) приводит к распределению вероятностей состояния на i-том шаге:

 

16. (i Δt) = πⁱ p(t=0) (6)

 

2. Для i-той степени матрицы перехода π можно записать выражение:

 

(7)

 

где обозначено

р₀() = β/(+β) , р₁() = /(+β) (8)

 

Так как |1 - - | < 1 , то с ростом числа шагов в выражении (7) для матрицы πⁱ остается только первое слагаемое.

3. Подставляя выражение (7) в формулу (6), получим:

Как и в (7), последнее слагаемое в сумме (9), содержащее параметры начального распределения вероятности, убывает с ростом i . С ростом числа звеньев цепи система «забывает» начальное состояние. Независимо от него, распределение вероятности обнаружить здорового или больного работника после достаточно большого промежутка времени от начала отсчета определяется соотношениями (8). В теории цепей Маркова распределение вероятности p(¥) называется финальным. Так же, как и начальное, финальное распределение удовлетворяет условию нормировки р₀ + р₁ = 1 . Непосредственной подстановкой нетрудно убедиться, что финальное распределение удовлетворяет уравнению

p(¥) = π p(¥) (10)

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

1.4. Определение (8) финального распределения больных и здоровых через индивидуальные вероятности заболеть a и выздороветь β является нетривиальным результатом теории цепей Маркова. Последовательность соотношений (4) – (10) демонстрирует как из показателей, характерных для отдельных работников, складывается результирующая статистическая картина эпидемиологической обстановки на производстве. Определенным аргументом в пользу такого подхода можно считать отмеченное выше совпадение формы выводов теории цепей Маркова и статистики ансамблей ЗВУТ в трудовых коллективах. В этой связи следует отметить, что в работе [6] прослежены параллели между описаниями ПОЗ на языке статистики ансамблей и в эпидемиологических терминах. Это обстоятельство позволяет надеяться, что и результаты теории цепей Маркова можно интерпретировать в системе эпидемиологических понятий. Последующее изложение иллюстрирует это.

2.Способы интерпретации результатов.

2.1. Один из распространенных способов интерпретации результатов эпидемиологических исследований профессионально обусловленных заболеваний (ПОЗ) проанализирован в работе [10]. Рассматривалось использование 2х2 таблиц сопряженности для оценки степени (условной вероятности) влияния ВПФ на уровень хронических ПЗ. Проблемы возникают при попытке использовать такие таблицы для анализа ЗВУТ. Во-первых, для характеристики уровня ЗВУТ принято использовать два параметра – количество случаев К и суммарную продолжительность D. Неясно, какие именно величины следует подставлять в колонку таблицы ПОЗ=1 (а и с в обозначениях [10]). Во-вторых, обе величины относятся к заболевшим работникам. Неясно, какими величинами b и d следует характеризовать (в колонке таблицы ПОЗ=0) здоровых работников. Последнее необходимо для оценок рисков (оценочных, относительных, атрибутивных) – основных эпидемиологических выводов, к которым должно приводить использование таблиц сопряженности.

Ответы на оба вопроса дают оценки (8) финальных вероятностей обнаружить работника здоровым или больным. Если обозначить через N₁ численность работников, подвергающихся воздействию ВПФ, а через N₀ , соответственно, численность контрольной группы, то следует полагать a = N₁[a/(a+β)]₁ , b = N₁[β/(a+β)]₁ . Здесь индексом 1 отмечены величины, относящиеся к работникам, подвергающимся воздействию ВПФ. Для контрольной группы c = N₀[a/(a+β)]₀ , d = N₀[β/(a+β)]₀ . После этого все величины, характеризующие риски ЗВУТ, оцениваются по тем же правилам, что и риски хронических ПЗ.

Определяются следующие виды риска:

• Оценочный риск заболевания – отношение числа случаев заболевания в когорте к численности когорты: (a + c)/(N₀ + N₁)
• Риск возникновения заболевания в группе, не подверженной воздействию ВПФ:

R₀ = c / N₀ = [a / (a + β)]₀ (11)

• Риск возникновения заболевания в группе, подверженной воздействию ВПФ:

R₁ = = a / N₁ = [a / (a + β)]₁ (12)

• Относительный риск (RR) – величина, описывающая увеличение риска заболевания под действием ВПФ. Это отношение вероятности заболевания у работников подвергающихся воздействию ВПФ к вероятности заболевания у работников, не подвергающихся такому воздействию.

RR = R₁ / R₀ = a*N₀ / c*N₁ = [a / (a + β)]₁ / [a / (a + β)]₀ (13)

 

2.2.Продемонстрируем возможности использования таблиц сопряженности для анализа воздействия вредных производственных факторы на работников, занятых в сварочном производстве, на трех крупных предприятиях машиностроения в г. Санкт - Петербурге: ОАО «Ижорские заводы», ФГУП «Адмиралтейские верфи» и ОАО судостроительный завод «Северная верфь» [11]. Для анализа использовались данные за 2007, 2008, 2009 гг.

Из основной группы (электросварщики) и контрольной группы (токари, фрезеровщики, слесари) было подобрано 170 пар (т.е. N₀ = N₁ =170), метчированных по возрасту и стажу между основной и контрольной группами. Исходные данные по основным показателям ЗВУТ исследуемых групп приведены в табл.1

Заметим, что вероятности выздоровления b в основной и контрольной группах расходятся незначительно (менее 5 %). Можно сделать вывод, что условия реабилитации (бытовые, медицинского обслуживания) в этих группах близки. В то же время вероятности заболевания a различаются значительнее (более 15 %), что свидетельствует о заметных различиях в условиях труда.

 

Построенная методом анализа цепей Маркова таблица сопряженности имеет вид:

Анализ данных табл.2 с использованием соотношений разд.2 позволяет утверждать следующее:

• Оценочный риск ЗВУТ в обследованном коллективе составляет 0,021.
• Риск ЗВУТ в группе, подвергающейся воздействию ВПФ, составляет 0,022.
• Риск ЗВУТ в группе, не подвергающейся воздействию ВПФ, составляет 0,020.
• Относительный риск RR = 1,11 ; соответственно, этиологическая доля EF = 9,8 % Согласно Руководству [12] такие величины RR и EF свидетельствуют о малой степени причинно-следственной связи нарушений здоровья с работой.

 

3.Динамика перехода ЗВУТ в ПЗ.

3.1. Предложенная в п.1 феноменологическая модель определяет сложение индивидуальных показателей ЗВУТ в результирующую статистическую картину эпидемиологической обстановки на производстве. Ее можно модифицировать для того, чтобы описывать постепенную эволюцию состояния здоровья отдельных работников под воздействием ВПФ. Резонно полагать, что болезнь уменьшает потенциальные возможности организма сопротивляться ей. Именно, можно считать, что вероятности заболеть a и выздороветь b меняются ото дня ко дню, т.е. зависит от номера i звена в цепи Маркова. Например, если при t_i = i*Δt работник болен, то при t_i+1 = (i+1)*Δt вероятность выздороветь b уменьшается на величину d. Аналогично с вероятностью заболеть a , однако, эта величина должна возрастать со временем. Так как вероятность быть больным при t_i = i*Δt равна p₁(t_i), то

a(t_i+1) = a(t_i) + d*p₁(t_i) ; b(t_i+1) = b(t_i) - d*p₁(t_i) (14)

Нетрудно видеть, что при сделанном предположении сумма a + b остается постоянной, т.е. не зависит от номера i звена в цепи Маркова:

a(t_i) + b(t_i) = a(0) + b(0) = a₀ + b₀ (15)

Результаты будут значительно проще и нагляднее, если перейти от дискретного к непрерывному времени, т.е. от цепи Маркрва к процессу Маркова. Уравнение для непрерывной функции р(t) º p₁(t) можно вывести, пассматривая ее приращение за время Δt :

Δp = p(t + Δt) - p(t) = {a*[1 - p(t)] - b*p(t)}*Δt (16)

Здесь явно учитывается, что вероятности a и b, а также 1- a и 1 - b отнесены к единице времени. Так как интервал времени Δt , равный 1 сут., значительно меньше характерных масштабов изменения параметров системы, то из (16) можно получить дифференциальное уравнение

dp/dt + (a + b)*p = a (17)

Решение уравнения (17) показывает, что за время » 1/(a+b) » l вероятность р₁ = р(t) обнаружить работника больным приходит к своему равновесному значению a/(a+b), которое, однако, сейчас не является финальным, т.к. вероятность заболеть a меняется со временем. Для a(t) имеем уравнение, следующее из (14):

da/dt = d*p(t) = d*a(t)/(a₀+b₀) (18)

Решение этого уравнения, удовлетворяющее начальному условию a(t=0) = a₀ , описывает лавинообразный (экспоненциальный) рост вероятности заболевания:

a(t) = a₀*exp[d*t/(a₀+b₀)] (19)

Соответственно, убывает вероятность выздоровления:

b(t) = b₀ – a₀*{exp[d*t/(a₀+b₀)] – 1} (20)

Одним из последствий таких вариаций является уменьшение до нуля вероятности выздоровления b. Это происходит за время T, определяемое из соотношения:

d*Т = (a₀+b₀)*ln[(a₀+b₀)/ a₀] » b₀*ln(b₀/a₀) (21)

В последнем равенстве учтено, что реально a₀ << b₀ . Эффект зануления вероятности выздоровления соответствует переходу ЗВУТ в хроническое ПЗ с постоянной утратой трудоспособности. Соотношение (21) дает зависимость времени такого перехода от условий труда (величины среднего периода циклов «заболевание-выздоровление» L ) и быта (средней длительности болезни l ).

3.2. Рассмотренная модель, описывающая протекание ПОЗ (от ЗВУТ до хронических ПЗ с постоянной утратой трудоспособности), дает возможность выявить параметры, определяющие количественные характеристики уровней заболеваемости и суммарную эпидемиологическую обстановку на производстве. Кроме того, она позволяет установить соотношения между этими характеристиками и их эволюцию в процессе перехода ЗВУТ в ПЗ. Возможности этой модели, однако, ограничиваются ее феноменологической природой. Она, в частности, не позволяет непосредственно определить скорость d изменения вероятностей заболевания и выздоровления. Возможно только определение произведения (21) этой скорости на среднее время Т перехода ЗВУТ в ПЗ. Последняя величина, однако, определяется из данных наблюдения этого процесса на реальных производствах. Для этого можно воспользоваться, например, информацией о ПЗ из сборника [13]. В нем приведены распределения ПЗ по времени контакта работников (мужчин и женщин поотдельности) с ВПФ. Данные дифференцированы по классам условий труда. Для работников мужчин эти данные приведены на рис.1.

 

Рисунок 1

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

Видно, что за исключением КУТ 2 и КУТ 4, для которых характерен заметный «отсев» работников в начале трудового стажа (см. рис.1), среднее время контакта составляет » 26 лет с незначительными вариациями. Отсев работников, занятых в опасных условиях труда (КУТ 4) понятен – именно с этим связан известный «эффект здорового работника» [14]. Для КУТ 2 такого естественного объяснения нет.

Возвращаясь к проблеме определения скорости d изменения вероятностей заболевания и выздоровления, можно заключить, что воздействие ВПФ приводит к уменьшению защитных возможностей организма на 4 – 5 % ежегодно. Рассмотренная простая модель цепи Маркова с двумя состояниями не допускает более детального сопоставления модельных и реальных кривых распределения времени контакта работника с ВПФ до регистрации ПЗ. Здесь требуются более сложные модели и более детальные результаты натурных исследований.

Заключение.

(1)Предложено феноменологическое описание ЗВУТ на языке теории цепей Маркова. Показано, что модель простой цепи с двумя состояниями удовлетворяет требованиям, которые следуют из рассмотрения

(2)Модель цепи Маркова описывает переход от индивидуальных характеристик ЗВУТ, являющихся входными параметрами модели, к статистическому описанию в терминах профессиональных рисков () эпидемиологической Такой подход позволяет наиболее оптимально использовать исходную информацию о количестве К случаев ЗВУТ и об их суммарной длительности в обследуемом рабочем коллективе.

(3)Модель допускает наглядное расширение для описания эволюции ПОЗ: от ЗВУТ до хронического ПЗ с постоянной утратой трудоспособности. Результаты легко интерпретируются в терминах ослабления защитных возможностей организма.

(4)Для распределения времени контакта работника с ВПФ до регистрации ПЗ требуются более сложные модели и более детальные результаты натурных исследований.

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

Литература

1.Федорович Г.В. Классификация условий труда по эпидемиологическим данным // БиОТ. – 2011. - № 4. – С. 49 –52.

2.Федорович Г.В. Профессиональный риск: количественная оценка и управление // БиОТ. – 2012. - № 1 – С. 60 – 64.

3.Федорович Г.В. Этиологические структуры профзаболеваний // БиОТ. – 2012. - № 4 – (в печати).

4.Методические рекомендации по углубленному изучению заболеваемости с временной утратой трудоспособности. Минздрав СССР, № 2484-81, 1981.

5.Положение о расследовании и учете профессиональных заболеваний. Приложение к Постановлению правительства РФ № 967, 2000.

6.Федорович Г.В. Статистика ансамблей в расчетах профессиональных рисков // БиОТ – 2010 - № 4 – С. 48 – 52.

7.Федорович Г.В. Статистические ансамбли временной утраты трудоспособности // Человек и труд – 2011- № 3 - С. 57 – 61.

8. Федорович Г.В. Внутренняя структура ансамблей временной утраты трудоспособности // Человек и труд – 2011- № 4 - С. 35 – 38.

9.Тихонов В.И., Миронов В.А. Марковские процессы. - М.: Советское радио, 1977 - 488 с.

10.Федорович Г.В. Эпидемиологический анализ характеристик профессионального риска // БиОТ – 2012 - № 3 – С. 41 – 45.

11.Кусраева З.С. Современный подход к оценке профессионального риска при выполнении электродуговой сварки и резки металлов // Доклад на XI съезде гигиенистов. М., апрель 2012 г.

12.Р 2.2.1766-03 «Руководство по оценке профессионального риска для здоровья работников. Организационно-методические основы, принципы и критерии оценки». - М.: Минздрав России, 2004. – 17 с.

13. Профессиональные заболевания и их распределение по классам условий труда в РФ в 2009 г.: Информационный сборник под ред. Верещагина А.И. – М., 2010 – 108 с.
14. Эпидемиологический словарь. Дж. М. Ласт (ред) – М., 2009 – 316 с.