Введение В задачах анализа данных типа времени жизни, когда изучаемой случайной величиной является время до наступления некоторого системного события, полученные выборки наблюдений нередко являются усеченными слева и (или) цензурированными справа [1]. Для объяснения понятия «усечение слева» рассмотрим следующий пример. Пусть имеется генеральная совокупность людей, страдающих определенной болезнью. Обозначим как функцию распределения случайной величины - времени с начала болезни до смерти. Предположим, что в момент времени началось обследование n пациентов с данной болезнью. Отметим, что люди из этой же генеральной совокупности, которые умерли до наступления момента времени , не включены в выборку. Обозначим как , , …, независимое время жизни с начала болезни до смерти пациентов, включенных в выборку. Необходимо отметить, что если рассматривать данную выборку как обычную выборку из распределения , то полученный результат будет слишком оптимистичным, поскольку чем больше продолжительность жизни пациента, тем больше у него шансов попасть в выборку, в то время как пациенты c малой продолжительностью жизни в выборку не попадают. Обозначим символом время с начала болезни i-го пациента до начала обследования . Тогда условное распределение случайной величины имеет вид , . Цензурирование возникает в случае, когда на момент окончания эксперимента системное событие для некоторых объектов исследуемой выборки еще не наступило либо объект выбыл из эксперимента до наступления системного события. Для вычисления оценок максимального правдоподобия параметров распределений по усеченным слева и цензурированным справа данным в [2-4] используется EM-алгоритм. В случае отсутствия априорной информации о виде распределения строят непараметрическую оценку Каплана - Мейера функции выживаемости [1, 5] или непараметрическую оценку функции плотности [6]. При анализе выживаемости типичной задачей является исследование зависимости вероятности дожития до некоторого момента времени от различных ковариат (факторов). В качестве ковариат могут выступать такие показатели, как возраст пациента, продолжительность болезни, схема лечения и др. Наиболее часто используемой моделью в анализе выживаемости является модель пропорциональных интенсивностей, предложенная Д. Р. Коксом [7]. Вопросам проверки гипотезы о согласии с параметрической моделью Кокса по цензурированным справа данным посвящены, в частности, работы [8-10]. В случае усеченных слева и цензурированных справа данных в [11] предложен параметрический критерий согласия для проверки гипотезы о согласии с моделью Кокса против обобщенной модели пропорциональных интенсивностей. В данной работе предлагается подход к проверке гипотезы о согласии с параметрической моделью пропорциональных интенсивностей Кокса по усеченным слева и цензурированным справа данным, основанный на формировании выборки остатков и использовании непараметрических критериев согласия Колмогорова, Крамера - Мизеса - Смирнова и Андерсона - Дарлинга. Методами компьютерного моделирования исследуются распределения статистик критериев согласия, а также распределения статистики критерия Вальда при проверке гипотезы о равенстве нулю параметров модели пропорциональных интенсивностей Кокса. Рассматривается задача анализа выживаемости некоренного населения Крайнего Севера в зависимости от наблюдаемых факторов. Модель пропорциональных интенсивностей Кокса Обозначим символом случайное время жизни пациента, которое зависит от вектора ковариат . Функция выживаемости определяется соотношением , а кумулятивная функция риска выражением , где - функция интенсивности отказов; - функция плотности распределения случайной величины . Модель пропорциональных интенсивностей Кокса имеет вид [7]: , (1) где - вектор параметров регрессии; - неотрицательная функция от ковариат; - базовая кумулятивная функция риска. Функция выживаемости для модели пропорциональных интенсивностей имеет вид . В качестве функции от ковариат обычно рассматривают логлинейную модель вида . Если в формуле (1) не вводится предположение относительно закона распределения продолжительности жизни, модель называется полупараметрической. Если же вводится параметризация, как для функции воздействий, так и для базовой кумулятивной функции риска , модель считается параметрической. Оценивание параметров модели Кокса по усеченным слева и цензурированным справа данным. Усеченную слева и цензурированную справа выборку можно представить в следующем виде: , (2) где - время наступления системного события или момент цензурирования (время завершения наблюдения за i-м индивидуумом); - время усечения; - индикатор цензурирования содержит информацию о причине прекращения наблюдения, и - значение вектора ковариат, . Если в ходе эксперимента было зафиксировано системное событие, то , и данное наблюдение называется полным. Если же неизвестно по причине окончания наблюдения в момент , то , и наблюдение называется цензурированным справа. Цензурированные справа выборки можно разделить на три основных типа и их комбинации. Если время эксперимента ограничено, т. е. наблюдение за индивидуумами ведется до заранее определенного момента времени , тогда , и полученная в результате выборка называется цензурированной I типа. В случае, если эксперимент продолжается до наступления определенного количества системных событий , и наблюдение за остальными индивидуумами прекращается в момент наступления -го системного события, то полученная в результате выборка наблюдений называется цензурированной II типа, и , где - время наступления последнего системного события. Если , - случайные величины, то выборка наблюдений называется цензурированной III типа или случайно цензурированной. Оценки неизвестных параметров модели (1) находят методом максимального правдоподобия. В случае полупараметрической модели максимизируют логарифм функции частичного правдоподобия [7]: , (3) где в качестве базовой кумулятивной функции риска используют непараметрическую оценку: . В случае параметрической модели логарифмическая функция правдоподобия имеет вид . Проверка гипотезы о незначимости регрессионных параметров модели Кокса. Одним из важнейших этапов построения регрессионной модели является отбор значимых регрессоров. В случае модели пропорциональных интенсивностей отбор значимых ковариат производится с использованием процедур включения или исключения на основе критериев Вальда или отношения правдоподобия. Рассмотрим гипотезу о незначимости вектора регрессионных параметров для модели Кокса (1), для проверки которой используется оценка параметра модели Кокса , полученная путем максимизации функции частичного правдоподобия (3). Статистика критерия Вальда для проверки данной гипотезы , (4) где ; в случае справедливости при подчиняется -распределению с степенями свободы [12]. Статистика Вальда для проверки гипотезы имеет вид , где - это i-я компонента вектора оцененных параметров ; - диагональный элемент матрицы, обратной к . При условии справедливости нулевой гипотезы в пределе при статистика подчиняется -распределению с одной степенью свободы. С использованием методики компьютерного моделирования и исследования статистических закономерностей нами было проведено исследование сходимости распределений статистики Вальда (4) к соответствующему предельному χ2-распределению в случае усеченных слева и цензурированных справа данных. На рис. 1 представлены графики зависимости расстояния между эмпирическим распределением статистики Вальда и предельным -распределением от объема выборки . Для построения эмпирических распределений генерировались усеченные слева и цензурированные (II типа) справа выборки в соответствии с заданной параметрической моделью Кокса, в которой значение регрессионного параметра было взято равным нулю; по каждой выборке вычислялось значение статистики (4). Объем выборки значений статистики . Рис. 1. Сходимость распределения статистики Вальда к предельному χ2-распределению Как видно из рис. 1, с ростом степени цензурирования распределения статистики Вальда медленнее сходятся к -распределению. При степени цензурирования 20 % предельным законом распределения можно пользоваться начиная уже с - расстояние от до χ2-распределения не превышает 0,03; при степени цензурирования 50 % - начиная с , а при степени цензурирования 80 % - с . С увеличением числа оцениваемых параметров данная зависимость сохраняется. Проверка гипотезы о согласии с параметрической моделью Кокса по усеченным слева и цензурированным справа данным Обязательным этапом построения параметрической вероятностной модели является проверка статистической гипотезы о согласии: , , . Один из подходов к проверке гипотезы о согласии с параметрической моделью пропорциональных интенсивностей Кокса основан на формировании остатков Кокса - Снелла [13], которые, в случае выборки вида (2), имеют вид , , и проверке гипотезы о принадлежности выборки остатков стандартному экспоненциальному распределению [9]. Заметим, что полученная выборка остатков является обычной цензурированной справа выборкой вида . Применение модифицированных критериев согласия Колмогорова, Крамера - Мизеса - Смирнова и Андерсона - Дарлинга для цензурированных I и II типа выборок рассмотрены в [14], а для случайно цензурированных выборок - в [15]. В модифицированном критерии Колмогорова для цензурированных выборок в качестве расстояния между эмпирическим и теоретическим законами распределения используется величина , где - оценка Каплана - Мейера [14] по выборке ; - теоретическая функция распределения, соответствующая проверяемой гипотезе, в данном случае это стандартное экспоненциальное распределение); - правая граница наблюдаемой области. В модифицированном критерии Крамера - Мизеса - Смирнова в качестве расстояния между распределениями используется величина . В модифицированном критерии Андерсона - Дарлинга в качестве меры рассматривается величина . Проверяемая гипотеза о согласии отвергается при больших значениях статистик. Аналитические выражения для распределений статистик рассматриваемых критериев неизвестны. Вследствие этого вычисление критических значений статистик (или достигнутых уровней значимости) при проверке гипотез с использованием данных критериев может опираться только на распределения статистик, полученные в результате статистического моделирования [15, 16]. Пример построения вероятностной модели выживаемости по усеченным слева и цензурированным справа данным Для исследования выживаемости некоренного населения Крайнего Севера в 1991 г. в г. Мирном были проведены сбор и анализ данных у лиц, прошедших скрининговое обследование. В 2000-2003 гг. был поведен сбор данных по выживаемости обследованных лиц. Факт смерти был зафиксирован у 20 % ранее обследованных лиц. К данному виду исхода была отнесена смерть от старости и заболеваний; в случае смерти от травм, отравлений и несчастных случаев считали, что пациент выбыл из исследования (цензурированное наблюдение). Кроме того, важно учесть, что в выборку попали только лица, живые на момент скрининга; информации об умерших до 1991 г. нет. Такие наблюдения являются усеченными слева, причем временем усечения является период от рождения до начала скрининга. Таким образом, имеем усеченную слева и цензурированную справа выборку наблюдений с ковариатами объема , степень цензурирования 80 %. В качестве ковариат рассмотрены такие факторы, как возраст приезда на Север, индекс массы тела, артериальное давление, наличие хронических патологий и др. - всего свыше 100 показателей. На основе полученной выборки построим модель пропорциональных интенсивностей Кокса (1) с логлинейной функцией от ковариат вида . Отбор статистически значимых ковариат в данном случае был проведен с использованием алгоритма включения на основе критерия Вальда. В результате отбора выявлены следующие прогностические факторы: - курение; - возраст приезда на Север; - измененные зубцы Q и QS в отведениях I, AVL, V6; - гипохолестеринемия липидов высокой плотности (ЛВП): мужчины ≤ 40 мг/дл; женщины ≤ 46 мг/дл; - перенесенный инсульт головного мозга; - наличие хронической патологии (класс 2) в соответствии с международной классификацией болезней (10-й пересмотр (МКБ-10)); - наличие хронической патологии (класс 3, МКБ-10); - наличие хронических осложнений после травм (класс 19, МКБ-10). Для построения прогноза выживаемости была введена параметризация базовой функции риска. По критерию Акаике в качестве наилучшего распределения среди рассмотренных семейств (экспоненциальное, Вейбулла, гамма-, логнормальное, логлогистическое, обобщенное распределение Вейбулла и обобщенное гамма-распределение) было выбрано обобщенное гамма-распределение с функцией плотности: . Гипотеза о согласии с полученной моделью пропорциональных интенсивностей была проверена с использованием рассмотренных в данной работе критериев. По критерию Колмогорова значение статистики , достигнутый уровень значимости ; по критерию Крамера - Мизеса - Смирнова , ; по критерию Андерсона - Дарлинга , . Вследствие того, что полученные достигнутые уровни значимости для всех рассмотренных критериев, нет оснований для отклонения гипотезы о согласии с выбранной моделью. Оценки максимального правдоподобия параметров построенной модели пропорциональных интенсивностей Кокса и результаты проверки гипотезы о незначимости регрессионных параметров представлены в таблице. Результаты оценивания параметров модели пропорциональных интенсивностей Кокса Параметр Оценка максимального правдоподобия параметра Статистика Вальда Достигнутый уровень значимости (масштаба) 71,4031 - - (формы) 6,0630 - - (формы) 47,6346 - - ( - курение) 0,8950 6,1779 0,0129 (- возраст приезда на Север) -0,0350 4,0589 0,0439 (- измененные зубцы Q и QS в отведениях I, AVL, V6) 3,8111 12,3103 0,0005 (- гипохолестеринемия ЛВП) 1,7748 5,6317 0,0176 (- перенесенный инсульт головного мозга) 1,6309 7,0379 0,008 (- наличие хронической патологии (класс 2, МКБ-10)) 3,4033 10,6136 0,0011 (- наличие хронической патологии (класс 3, МКБ-10)) 2,2176 3,875 0,049 (- наличие хронических осложнений после травм (класс 19, МКБ-10)) 2,3693 5,1028 0,0239 Согласно данным таблицы, факторы, выявленные в результате отбора, оказывают статистически значимое влияние на выживаемость исследуемой группы людей, поскольку достигнутые уровни значимости (p-value), полученные при проверке гипотез о равенстве регрессионных параметров нулю с помощью критерия Вальда, меньше заданного уровня значимости - . Наличие таких факторов риска, как курение, измененные зубцы Q и QS в отведениях I, AVL, V6, гипохолестеринемия ЛВП, перенесенный инсульт головного мозга, а также наличие хронических патологий 2 и 3 класса и хронических осложнений после травм значимо сокращает длительность жизни. Полученная отрицательная оценка параметра при переменной (возраст приезда на Север) свидетельствует о том, что чем меньше возраст приезда на Север, т. е. чем больше северный стаж, тем ниже продолжительность жизни. Рассмотрим, как влияют выявленные факторы на функцию выживаемости . На рис. 2 представлены графики функции выживаемости при различных значениях вектора ковариат . На рис. 2, а, в-з графики функции выживаемости строились при среднем возрасте приезда на Север, равном 25 годам. По рисункам можно судить о том, насколько сильно влияет каждый из выявленных факторов на функцию выживаемости при отсутствии иных факторов риска. Напомним, что значение функции выживаемости представляет собой вероятность дожития до момента времени , время измеряется в годах. а б Рис. 2. Влияние на выживаемость различных факторов: а - курение; б - возраст приезда на Север в г д е ж з Рис. 2. Влияние на выживаемость различных факторов: в - зубцы Q и QS в отведениях I, AVL, V6; г - гипохолестеринемия ЛВП (мужчины ≤ 40 мг/дл; женщины ≤ 46 мг/дл; д - перенесенный инсульт головного мозга; е - наличие онкологической патологии (класс 2); ж - наличие хронической патологии (класс 3); з - наличие хронической патологии (класс 19) Согласно данным на рис. 2, наибольшее значение для выживаемости из выявленных факторов имеет нарушение функционирования миокарда, что отражено измененными зубцами Q и QS в отведениях I, AVL, V6 (рис. 2, в) и наличием онкологической патологии класса 2 (рис. 2, е). Интересным оказалось влияние фактора «Возраст приезда на Север» - чем в более раннем возрасте человек приезжает на Крайний Север для постоянной работы на одном из промышленных предприятий, тем статистически меньше продолжительность его жизни (рис. 3). Любопытной представляется оценка того, насколько большой может быть ошибка прогноза вероятности дожития, если при исследовании не учитывать усечение выборки. Так, на рис. 3 представлены графики оценок Каплана - Мейера и теоретических функций выживаемости, соответствующих обобщенному гамма-распределению с оцененными параметрами по усеченной выборке и выборке с теми же значениями продолжительности жизни, но без усечения. Рис. 3. Теоретическая функция выживаемости и оценка Каплана - Мейера, полученные с учетом и без учета усечения данных Как видим, модель выживаемости, полученная с учетом усечения слева, оказывается менее оптимистичной. В частности, вероятность дожития до 55 лет на основе модели, построенной без учета усечения, составляет примерно 0,9, в то время как на основе модели, построенной с учетом усечения, эта вероятность равна 0,75. Модель выживаемости по усеченным слева данным оказывается более «жесткой», что позволяет выбирать большее количество участников исследования в группу риска. Заключение Таким образом, основное внимание в ходе исследований уделялось проблеме построения модели пропорциональных интенсивностей Кокса и проверке гипотезы о согласии с данной моделью по усеченным слева и цензурированным справа данным. Проведено исследование распределений статистики критерия Вальда для проверки гипотез о равенстве регрессионных параметров нулю. Получены оценки скорости сходимости распределения статистики Вальда к соответствующему предельному χ2-распределению при различной степени цензурирования. Предложен метод проверки гипотезы о согласии по усеченным данным, основанный на преобразовании исходной выборки наблюдений и вычислении остатков. В результате исследований показано, что задача проверки гипотезы по усеченным слева и цензурированным справа данным сводится к задаче проверки гипотезы по цензурированным данным. Рассмотрен пример построения вероятностной модели выживаемости, на основе полупараметрической модели пропорциональных интенсивностей выявлены прогностические факторы. Для построения прогноза введена параметризация базовой функции интенсивности, соответствующая обобщенному гамма-распределению. С использованием модифицированных критериев Колмогорова, Крамера - Мизеса - Смирнова и Андерсона - Дарлинга проверена гипотеза о согласии с полученной моделью.