Введение в анализ выживаемости
Методы анализа выживаемости интенсивно применяются в медицине, биологии, страховании и промышленности.
Одной из важных характеристик, описывающих течение болезни, является продолжительность жизни пациентов с момента поступления в клинику или после проведения операции.
В принципе, для описания средних времен жизни и сравнения новой методики со старой можно использовать стандартные статистические методы.
Однако рассматриваемые данные имеют специфику, которую следует учитывать. Дело в том, что в медицинской практике мы часто имеем дело с неполными данными.
Это связано с тем, что трудно наблюдать все время жизни пациента после операции, так как пациент мог быть выписан или переведен в другую клинику и связь с ним была утеряна. При этом мы располагаем не полной информацией о времени жизни пациента, а лишь частичной.
Естественное желание исследователя использовать все данные, т. е. анализировать как полные времена жизни, так и неполные, и не терять с трудом собранную информацию.
Для этого и предназначены методы анализа выживаемости, которые позволяют изучать неполные или цензурированные данные.
Наблюдения, которые содержат неполную информацию, называются неполными или цензурированными (например, «пациент А был жив по крайней мере 4 месяца после того, как был переведен в другую клинику и контакт с ним был потерян»). Это пример цензурированного наблюдения: информация о том, что пациент был жив 4 месяца, важна и может быть использована для построения оценок.
Наблюдения от момента операции до летального исхода называется полными.
Итак, в анализе выживаемости различают полные (по-английски complete) и неполные, или цензурированные, наблюдения (по-английски censored).
Конечно, можно было использовать только полные времена жизни, но тогда мы имели бы в своем распоряжении очень мало наблюдений и соответственно неточные оценки.
Использование, наряду с полными наблюдениями, неполных или цензурированных наблюдений является главной особенностью методов анализа выживаемости.
Прежде всего, постараемся оценить вероятность того, что пациент прожил больше t дней после операции. Это важный показатель, называемый функцией выживания.
Наиболее естественный способ описания функции выживаемости состоит в построении Таблиц времен жизни.
Это один из старейших приемов анализа данных о выживаемости и традиционно используется, например, в страховании, где такие таблицы называются таблицами дожития.
Организация данных
Исходный файл данных имеет вид:
Организация файла следующая.
Пациенты располагаются в строках. В столбцах записаны даты операции и даты завершения пребывания в больнице. Например, из • первой строки видно, что пациенту была сделана операция 6 января 1968 (первые три клетки), выписался 21 января 1968 года (вторая тройка клеток). Далее связь с ним была утеряна, таким образом, это неполное наблюдение (значение переменной номер 7 — censored).
Восьмая переменная AGE содержит возраст пациентов.
Переменные 9, 10 содержат специальную медицинскую информацию об особенностях операции.
Значение переменной 11 — название госпиталя, где сделана операция. Ниже показана таблица жизни для этого файла данных.
Конечно, подобную таблицу жизни можно рассматривать как «расцшренную» таблицу частот. Однако обычная таблица частот строится по полным наблюдениям. В таблице жизни учтены как полные, так и неполные наблюдения.
Идея таблиц жизни, или дожития, в терминологии страхования, проста. Нам нужно вычислить простейшие статистики, чтобы описать время выживания пациентов.
Для этого временная ось разбивается на некоторое число интервалов. В приведенной выше таблице это число равно 12. В системе STATISTICA количество интервалов на временной оси пользователь может выбрать по своему усмотрению.
Для каждого интервала вычисляется число объектов, которые в начале рассматриваемого интервала были «живы» (см. соответствующий столбец в электронной таблице — переменная ЧИСЛО В НАЧАЛЕ), и число объектов, которые «умерли» в данном интервале (переменная ЧИСЛО УМЕРШИХ).
Также вычисляется число цензурированных или изъятых из наблюдения объектов на каждом интервале — переменная ЧИСЛО ИЗЪЯТЫХ (в таблицах жизни употребляют термин изъятые — withdrawn для цензурированных наблюдений, в данном примере это выписанные больные). Вычисляются доли этих объектов.
Для понимания таблиц полезно помнить, что на данном временном интервале наблюдение может быть либо цензурировано (больной выписан или переведен в другую клинику), либо наблюдается фатальный исход.
Рассмотрим более формально переменные в электронной таблице жизни.
Число в начале
Это число объектов, которые были «живы» в начале рассматриваемого временного интервала.
Число изъятых
Это число цензурированных на данном интервале объектов (объектов, изъятых из
наблюдения). Эти объекты имеют метку цензурированые (censored).
Число изучаемых
Это число объектов, которые были «живы» в начале рассматриваемого временного интервала, минус половина от числа изъятых.
Число умерших
Это число объектов, умерших на данном интервале. Умершие или отказавшие объекты обычно имеют метку complete.
Доля умерших
Эта отношение числа объектов, умерших в соответствующем интервале, к числу объектов, изучаемых на этом интервале.
Прокрутим электронную таблицу вправо и рассмотрим оставшиеся переменные таблицы.
Доля выживших
Эта доля равна единице минус доля умерших.
Кумулятивная доля выживших объектов, или функция выживания Это — оценка функции выживания, то есть вероятность того, что пациент переживет данный интервал. Она равна произведению долей выживших объектов по всем предыдущим интервалам. Если посмотреть на столбец КУМ.ДОЛЯ ВЫЖИВШ. приведенной выше таблицы, то можно увидеть, например, что 0,582759= 0,672414 × 0,866667, 0,569514= 0,582759×0,977273 и т. д.
Плотность вероятности
Это плотность вероятности смерти на данном интервале, когда из функции выживания на данном интервале вычитается функция выживания на следующем интервале и делится на длину интервала, показанную во втором столбце таблицы.
Например, (1 - 0,672414)/161,3636= 0,00203.
На графике оценки плотности можно видно, что вероятность смерти в первые 160 дней после операции максимальна. Далее она резко падает.
Большие вероятности смерти расположены также в интервалах от 161 до 322, от 968 до 1129 и др.
Интенсивность отказов, или функция мгновенного риска
Это также одна из важных характеристик, описывающих течение болезни. Функция мгновенного риска является важной прогностической характеристикой, описывающей течение болезни.
Формально функция мгновенного риска равна вероятности того, что пациент умрет в следующем интервале наблюдения, при условии, что вначале он был жив.
График функции риска достаточно наглядно показывает, что в первые дни риск смерти очень велик, затем он падает и спустя некоторое время вновь начинает возрастать. Заметим, что именно функция риска используется для прогностических целей.
Позвольте сделать отступление. Один из лейтмотивов нашей книги является непредвзятость и критическое отношение к полученным результатам. Такая критичность особенно важна в медицине. Мы доверяем результатам, полученным с помощью компьютера, однако всесторонне их проверяем.
Итак, нас интересует функция риска, однако реально мы получаем лишь оценку риска. Поэтому важна точность полученных оценок. Из простых соображений следует, что мы не доверяем оценкам с большой погрешностью. Например, мы не будем доверять оценкам, погрешность которых имеет тот же порядок, что и сами оценки. Поэтому внимательно просмотрите построенную таблицу и выбросите из нее плохие оценки (оценки с большой погрешностью). Это чрезвычайно важный принцип анализа данных!
Известно, что для получения надежных оценок параметров и ошибок в таблицах жизни требуется как минимум 30 наблюдений.
Взгляните на таблицу. Заметьте, в ней наряду с оценками приведены стандартные ошибки полученных оценок.
Медиана ожидаемого времени жизни
Это моменты времени, в которых функция выживания равна ?. Например, из первой строчки таблицы следует, что пациент с вероятностью ? проживет больше 809 дней после операции.
Если пациент пережил первый интервал, то медиана его времени жизни равна 1036 и т. д.
В общем случае таблица времен жизни дает хорошее представление о распределении отказов или смертей, если наблюдений достаточно много.
Однако для прогноза часто необходимо знать форму функции выживания. Для этой цели используются различные семейства распределений.
Наиболее важны следующие семейства распределений: экспоненциальное, Вей-булла и распределение Гомперца.
Эти распределения имеют неизвестные параметры, которые программа оценивает. Процедура оценивания параметров основана на методе наименьших квадратов. Для проведения оценивания применима модель линейной регрессии, поскольку все перечисленные семейства распределений могут быть «сведены к линейным» (относительно параметров) с помощью подходящих преобразований. Такие преобразования приводят иногда к тому, что дисперсия остатков зависит от интервалов (то есть дисперсия различна на разных интервалах). Чтобы учесть это, в алгоритмах подгонки дополнительно используются оценки взвешенных наименьших квадратов двух типов.
Напомним, что одна из задач в анализе выживаемости состоит в том, чтобы оценить функцию выживания, то есть вероятность того, что пациент проживет определенное время после операции.
Оказывается, что для цензурированных наблюдений функцию выживания можно оценить непосредственно, не используя таблицу времен жизни. Такой метод впервые предложили Каплан и Мейер в 1958 году.
Представьте, что вы имеете файл, в котором записаны в хронологическом порядке отдельные события. Тогда имеет место следующая оценка функции выживания:
S(t) = П[(n-j)/(n-j+l)õ(j)]
В этом выражении S(t) — оценка функции выживания, n — общее число событий (объем выборки), j — порядковый (хронологически) номер отдельного события, õ(j) равно 1, если j-е событие означает отказ (смерть) и õ(j) равно 0, если j-е событие означает потерю наблюдения (индикатор цензурирования), П означает произведение по всем наблюдениям j, завершившимся к моменту t.
Данная оценка функции выживания состоит из произведения нескольких сомножителей, поэтому она также называется множительной оценкой.
Рассмотрим тот же файл данных; что и для таблиц времен жизни. Оценка Каплана—Мейера функции выживания, построенная по этим данным, показана в следующей таблице:
Из таблицы видно, например, что вероятность того, что пациент проживет больше 25 дней, равна 0,966, вероятность того, что пациент проживет больше 39 дней, равна 0,9299 и т. д.
В первом столбце таблицы показаны номера наблюдений, для которых в данный момент времени произошло некоторое событие, знак + означает, что пациент цензурирован (был выписан).
Прокрутите электронную таблицу с результатами вниз по временной оси:
Обратите внимание на ошибки оценок. Стандартная ошибка функция выживания достаточно мала (сравните с ошибками для таблиц времен жизни). Ниже показан график функции выживания.
Отметим, что для удобства интерпретации на графике полные наблюдения помечены точками, неполные наблюдения отмечены крестиками.
Преимущество метода Каплана—Мейера (по сравнению с методом таблиц жизни) состоит в том, что оценки не зависят от разбиения времен жизни на интервалы.
Таким образом, нам не нужно разбивать временную ось на интервалы. Оценки Каплана—Мейера строятся в STATISTICA одним щелчком мыши.
Сравнение выживаемости в группах
Интересно сравнить времена жизни пациентов в различных группах, например, в группах: мужчины и женщины. В STATISTICA имеются специальные процедуры для сравнения выживаемости в группах.
Если количество групп две, то используется диалог Сравнение двух выборок.
Если количество групп больше двух, то используется диалог Сравнение нескольких выборок.
Для сравнения выживаемости в группах имеется несколько критериев: вариант известного непараметрического критерия Вилкоксона, предложенный для неполных наблюдений Геханом и Пето, а также F-критерий Кокса и логарифмический ранговый критерий.
Большинство этих критериев приводят соответствующие z-значения (нормального приближения), которые могут быть использованы для статистической проверки различий между группами.
Однако критерии дают надежные результаты лишь при достаточно больших объемах выборок. При малых объемах выборок эти критерии не столь надежны. В любом случае всегда полезны визуальные методы.
Эти графики позволяют увидеть различие между группами.
Кроме этого STATISTICA содержит программу на STATISTICA BASIC (файл Ma.nthaen.stb), вычисляющую критерий Ментела-Хенцела для сравнения двух групп данных (см. Lee E. Т. (1992) Statistical methods for survival data analysis). Этот критерий может быть полезен во многих клинических и эпидемиологических работах для того, чтобы контролировать эффект смешивающих переменных.
Критерий основан на анализе таблиц 2×2 (например, Группировка 1/2 и Выживаемость), стратифицированных или расслоенных с помощью категориальной переменной (смешанной переменной; например, Положением). Критерий позволяет проверить, являются две переменные в таблицах 2×2, например, переменная Группировка и Выживаемость, зависимыми или нет.
Не существует твердо установленных рекомендаций по применению определенных критериев.
Известно, что F-критерий Кокса обычно мощнее, чем критерий Вилкоксона— Гехана, если:
В работе Lee, Desu, and Gehan (1975) A Monte-Carlo study of the power of some two-sample tests, Biometrika, 62, p. 425-532, критерий Гехана сравнивался с некоторыми другими критериями. Показано, например, что критерий Кокса—Ментела и логарифмический ранговый критерий являются более мощными, если выборки имеет определенное распределение, например, экспоненциальное или Вейбулла. При этих условиях между критерием Кокса—Ментела и логарифмическим ранговым критерием почти нет различия.
В работе Ли (Lee Е. Т. (1980) Statistical methods for survival data analysis. Belmont, CA: Lifitime Learning) обсуждается мощность различных критериев более детально. Если вас затрудняет выбор определенного критерия, рекомендуем обратиться к этим работам.
Если сравниваются две или более группы, важно проверить доли цензурированных наблюдений в каждой. В частности, в медицинских исследованиях степень цензурирования может зависеть, например, от различий в методике лечения: пациенты, которым стало много лучше или стало хуже, с большой вероятностью теряются из наблюдения. Различие в степени цензурирования может привести к смещению в статистических выводах.
Это очень важный момент. Чтобы подогнать результат, недобросовестный исследователь может искусственно исключить из исследования тяжелых больных. Поэтому при проведении сравнения различных методик нужно руководствоваться здравым смыслом. Ясно, что если в одной группе доля цензурированных наблюдений существенно больше, чем в другой, нужно принять естественные меры предосторожности, по крайней мере, точно указать проблему.
Регрессионные модели в анализе выживаемости
В предыдущих разделах мы кратко обсуждали задачу оценивания функции выживания на основе реальных данных.
Более трудной задачей является оценка функции мгновенного риска, которая представляет собой вероятность летального исхода в малый промежуток времени при условии, что в начале исследуемого промежутка пациент был жив. Это важная характеристика прогноза развития болезни.
Непосредственная оценка функции мгновенного риска может потребовать большого количества наблюдений, поэтому применяются специальные модели, одна из которых — это модель Кокса пропорциональных рисков, или, на языке теории надежности, пропорциональных интенсивностей.
Большая проблема медицинских и биологических исследований состоит в выяснении того, являются ли некоторые переменные связанными с наблюдаемыми временами жизни. Если зависимость есть, то ее нужно оценить численно.
Существуют две главные причины, по которым в таких исследованиях нельзя непосредственно использовать классическую регрессию. Во-первых, времена жизни обычно не являются простыми линейными функциями от соответствующих регрессоров, поэтому анализ методами Множественной регрессии может привести к ошибочным выводам, например, не позволит обнаружить важных регрессоров. Во-вторых, вновь возникает проблема неполных наблюдений, т. к. некоторые наблюдения могут быть незавершенными.
Анализ выживаемости предлагает пять общих регрессионных моделей для неполных данных:
1) Модель пропорциональных интенсивностей Кокса (Сох (1972) Regression models and life tables, Journal of the Royal Statistical Sociaty, 34, p. 187-220);
2) Модель Кокса с зависящими от времени ковариатами;
3) Экспоненциальную регрессионную модель (см. книги Prentice (1973) Exponential survivals with censoring and explanatory variables, Biometrika, 60, p. 279-288);
4) Нормальную линейную регрессионную модель (см. например, Wolynetz (1979) Maximum likelihood estimation in a linear model from confined and censored normal data, Applied Statistics, 28, p. 185-206);
5) Логнормальную линейную регрессионную модель (являющуюся модификацией нормальной модели).
Для каждой из этих моделей STATISTICA позволяет вычислить оценки максимального правдоподобия (Maximum likelihood estimations).
Модель пропорциональных интенсивностей или пропорциональных рисков Кокса — наиболее общая регрессионная модель, в которой предполагается, что функция интенсивности имеет вид: h(t) = h0(t) y(z1,...,zm). Множитель h0(t) называется базовой функцией интенсивности.
Модель может быть параметризована, например, в виде:
h[(t),(z1, z2,..., zm)] = h0(t)× exp(b1 x z1 +...+ bm × zm)
Заметьте, в правой части стоит произведение двух функций, причем каждая из них зависит от своего множества переменных.
Функция интенсивности h0(t) может рассматриваться как функция интенсивности при равенстве нулю всех ковариат. Она не зависит от переменных z (называемых ковариатами). Второй сомножитель зависит от переменных z, которые, возможно, зависят от t.
Приведем пример такой модели.
Пусть изучается воздействие некоторого препарата на состояние больного, a z — категориальная переменная со значениями 1 для больных, принимавших новое лекарство, и 0 — для больных, не принимавших это лекарство. Тогда функцию риска можно записать в виде:
h(t,z) = h0(t) × exp{b1× z+b2× [z × log(t)-100]}
Обратите внимание, что функция интенсивности в момент t (левая часть формулы) есть функция: 1) функции интенсивности h0, 2) ковариаты z и 3) z, умноженной на логарифм времени.
Умножение ковариаты z на логарифм времени позволяет учесть, например, фактор времени при приеме нового лекарства.
Константа 100 в этом примере использована просто как нормировка, т. к. среднее логарифма времени жизни для этого множества данных равно 100.
Зная оценки параметров b1,b2 и функцию интенсивности h0, можно оценить функцию мгновенного риска через время t после операции.
Самое замечательное, что такие модели позволяют учитывать интуицию медицинских исследователей. Построение и оценка адекватности модели в конкретных исследованиях — отдельная нетривиальная задача.
Другой пример, h(t,s,x)- риск коронарной смерти для пациента возраста t лет при условии, что в возрасте s его систолическое артериальное давление было х (см. Meshalkin L. D., Kagan А. В. (1972) A contribution to the discussion upon the paper «Regression models and life tables» by D. R. Cox, J. R. Statist. Soc. Ser. B, № 2).
Итак, функция мгновенного риска в модели Кокса представлена в виде произведения двух сомножителей, один из которых характеризует объект, другой — базовую функцию мгновенного риска.
Предикторы определяются постановкой задачи, например, пол пациента, возраст, наличие определенных сопутствующих заболеваний или прием нового лекарства. Выбор предикторов определяется интуицией исследователя. Врач может попытаться предсказать на основе определенного набора предикторов степень риска на ближайшие несколько дней. Имея прогноз, он может изменить методику лечения.
Займемся некоторой математической кухней. Модель Кокса можно линеаризовать, поделив обе части соотношения на h0(t) и взяв натуральный логарифм от обеих частей:
log{h[(t),(z...)]/h0(t)} = b1× 2,+...+ bm× zm
Таким образом, мы получили линейную модель.
Итак, еще раз отметим, в основе модели Кокса лежат два предположения. Во-первых, зависимость между функцией интенсивности и логлинейной функцией ковариат является мультипликативной. Это предположение называется гипотезой пропорциональности. Реально оно означает, что для двух заданных наблюдений с различными значениями независимых переменных отношение их функций интенсивности не зависит от времени (чтобы ослабить это предположение, используются ковариаты, зависящие от времени; см. ниже). Второе предположение состоит в логлинейной зависимости функции интенсивности и регрессоров.
Предположение пропорциональности рисков часто подвергается сомнению. Например, рассмотрим гипотетическое исследование, в котором ковариатой является категориальная переменная, а именно индикатор того, подвергнут пациент хирургической операции или нет. Пусть пациент 1 подвергнут операции, в то время как пациент 2 — нет.
Согласно предположению пропорциональности, отношение функций интенсивностей для обоих пациентов не зависит от времени и означает, что риск для прооперированного пациента постоянно более высокий (или более низкий), чем риск пациента, не подвергнутого операции (при условии, что оба дожили до рассматриваемого момента).
Реалистичней другая модель, когда сразу после операции риск прооперированного пациента выше, но при благоприятном исходе операции с течением времени убывает и становится меньше риска не оперированного пациента. В этом случае используются регрессоры, зависящие от времени.
Можно привести много других примеров, где предположение о пропорциональности неприемлемо. Так, при изучении физического здоровья возраст является одним из факторов выживаемости после хирургической операции. Ясно, что возраст — более важный предиктор для риска сразу после операции, чем по прошествии некоторого времени после операции (например, вслед за первыми признаками выздоровления).
В случае категориальных ковариат, например, учитывающих, был или не был пациент подвергнут хирургической операции, рекомендуется обратиться к стратифицированному анализу выживаемости, в котором, исходя из априорных знаний, исследователь разбивает пациентов на однородные по фактору риска группы.
Можно провести подгонку модели пропорциональных интенсивностей отдельно для каждой группы наблюдений. Таким образом, можно явно представить функцию интенсивности для каждой группы. Иногда предположение пропорциональности не выполняется. В таком случае можно явно определить ковариаты как функции времени.
В главе Подгонка вероятностных распределений показано, как с помощью критерия хи-квадрат проверяется выполнимость предположений модели Кокса в системе STATISTICS
Заметим, что арифметические выражения, которые определяют ковариаты, не должны содержать ссылок на длительности жизни. Однако допускается, чтобы некоторые ковариаты были функциями двух или большего числа других ковариат. Это, например, удобно в моделях многофакторных экспериментов. Для каждого фактора можно создать переменную в файле данных, чтобы установить желаемые контрасты. Логика и выбор априорных значений коэффициентов контрастов те же, что и в дисперсионном анализе. Если специфицируются ковариаты для регрессионной модели пропорциональных интенсивностей, то можно также определить взаимодействия факторов.
Например, предположим, что фактор А имеет 2 уровня. Всем субъектам, отнесенным к первому уровню этого фактора, мы приписываем -1 как значение соответствующей переменной (переменной А) в файле данных. Аналогично всем субъектам, отнесенным ко второму уровню, приписываем значение +1. Второй фактор, также с двумя уровнями, будет закодирован тем же способом (переменная В). После того как переменные АиВ определены как ковариаты, выражение А *В есть третья ковариата для проверки взаимодействия между этими двумя факторами.
Для задания зависящих от времени ковариат можно использовать тот же самый синтаксис, который используется в формулах электронной таблицы.
В некоторых случаях есть основание предполагать, что влияние одной или нескольких ковариат на функцию интенсивности не является непрерывным по времени. Например, риск для пациента после операции может зависеть от времени, прошедшего после операции в течение первых двух дней, и, во вторую очередь, от некоторых других факторов. В таком случае можно использовать некоторые логические операции, которые также поддерживаются при вводе формул электронных таблиц.
Например, можно определить зависящую от времени ковариату с помощью следующего выражения:
Age × (T_<2)
Логическое выражение Т< 2 равно 0 (ложь), если после операции прошло больше 2 дней, и равно 1 (истина), если меньше. Таким образом, здесь явно учтен эффект первых двух послеоперационных дней.
Эта модель записывается в виде:
S(z) = ехр(а + b1 × Z1 + b2 × z2 + ... + bm × zm)
S(z) обозначает время жизни, а — неизвестная константа, bi — параметры регрессии.
Вновь можно использовать критерий согласия хи-квадрат, чтобы оценить адекватность модели.
Статистика хи-квадрат может быть вычислена как функция логарифма правдоподобия для модели со всеми оцененными параметрами (L1)и логарифма правдоподобия модели, в которой все ковариаты обращаются в 0 (L0).
Если значение хи-квадрат значимо, отвергаем нулевую гипотезу и принимаем, что независимые переменные значимо влияют на время жизни.
Один из способов проверить адекватность экспоненциальной модели — построить остатки времен жизни и сравнить их со значениями стандартных экспоненциальных порядковых статистик.
Если предположение о том, что данные имеют экспоненциальное распределение, справедливо, то все точки на графике хорошо ложатся на прямую линию.
Нормальная и логнормальная регрессия
В этой модели предполагается, что времена жизни (или их логарифмы) имеют нормальное распределение. Модель совпадает с обычной моделью множественной регрессии и может быть записана следующим образом:
t = a + b1×z1 + b1×z1 + ... + bm × zm ,
где t — время жизни.
Если принимается модель логнормальной регрессии, то t заменяется ln t.
Модель нормальной регрессии особенно полезна, поскольку часто данные можно преобразовать в приблизительно нормальные с помощью подходящего преобразования.
Таким образом, в некотором смысле это наиболее общая параметрическая модель (в противоположность модели пропорциональных интенсивностей Кокса, которая является непараметрической).
Для всех регрессионных моделей в системе STATISTICA доступен стратифицированный анализ, который открывается в окне Результаты.
Цель стратифицированного анализа — проверить гипотезу о том, что одна и та же регрессионная кривая подходит для разных групп данных. Итак, стандартным образом мы разбиваем данные на несколько однородных групп.
Затем строятся регрессионные модели отдельно для каждой группы. Сумма логарифмов правдоподобия для разных моделей представляет собой логарифм правдоподобия модели с разными коэффициентами регрессии (и свободными членами, если требуется) в разных группах.
Далее ко всем данным обычным образом подгоняется регрессионная модель, не учитывая разбиение на группы, и вычисляется общий логарифм правдоподобия. По разности двух логарифмов правдоподобия проверяется значимость различия между группами.
В стратифицированном анализе на основе априорных соображений исследователь разбивает объекты на однородные группы риска, которые называются стратами, и проводит регрессионный анализ внутри каждой группы (см. например, книгу Кокрен У. (1976) «Методы выборочного исследования», где всесторонне обсуждаются методы построения групп). Во многих ситуациях риск-группы заранее известны, технически их можно получить, введя группирующие переменные.
Для модели пропорциональных интенсивностей Кокса система STATISTICA предлагает опцию подгонки к стратифицированным данным модели с общими коэффициентами для разных групп, но с разными базовыми функциями интенсивности. В результате наблюдения в отдельной группе удовлетворяют предположению пропорциональности, но это предположение не обязательно выполняется для наблюдений объединенных групп.
STATISTICA позволяет исследовать модель Кокса с ковариатами, зависящими от времени, а также сравнить модель с зависимыми от времени ковариатами и постоянными ковариатами.
Подробное введение в анализ выживаемости можно найти, например, в работах Bain (1978), Barlow and Proschan (1975) — русский перевод: Барлоу Р., Прошан Ф. Статистическая теория надежности и испытаний на безотказность. М. Наука, 1984, Сох and Oakes (1984) — русский перевод: Кокс Д. Р., Дукс Д. Анализ данных типа времени жизни. М., Финансы и статистика, 1988, Elandt-Johnson and Johnson (1980), Gross and Clark (1975), Lawless (1982), Lee (1980, 1992), Miller (1981), and Nelson (1982). Инженерные приложения этой техники обсуждены у Hahn and Shapiro (1967) — русский перевод: Хан Г., Шапиро С. Статистические модели в инженерных задачах. М., Мир, 1969.
На этом мы закончим общий обзор методов анализа выживаемости и перейдем к их реализации в системе STATISTICA, а также к примерам.
Модуль Анализ выживаемости системы STATISTICA предназначен для анализа цензурированных или неполных данных о выживаемости и отказах.
Модуль содержит процедуры для описания времен жизни и оценивания функций выживания, интенсивности и плотности вероятности, для подгонки теоретических распределений выживаемости к данным и для сравнения выживаемости в двух и более выборках. Модуль Анализ выживаемости содержит также регрессионные процедуры для подгонки объясняющих моделей к цензурированным данным (модель пропорциональных интенсивностей Кокса, в том числе с зависящими от времени ковариатами, экспоненциальная регрессия, нормальная и логнормальная регрессия).
Все процедуры в модуле Анализ выживаемости автоматически преобразуют данные в числовой формат. Таким образом, чтобы получить интересующие данные, пользователь может записать даты начала и даты окончания наблюдений, связанные с отказами или цензурированием (потерями объектов).
Таблицы времен жизни могут быть построены по исходным данным. Однако можно анализировать и готовые таблицы времен жизни.
Для всех регрессионных моделей доступны оценки максимального правдоподобия. При вычислении этих оценок для моделей пропорциональных интенсивностей и экспоненциальной регрессионной модели используется процедура безусловной максимизации. Для нормальной и логнормальной регрессионных моделей оценки параметров проводятся с помощью EМ-алгоритма. Этот алгоритм был впервые предложен в работе Dempster, Laird, and Rubin (1977) Maximum likelihood from incomplete data via the EM algorithm, Journal of the Royal Statistical Sociaty, 39, p. 1-38, и обсуждается в книге Сох and Oakes (1984) Analysis of survival data, New York: Chapman&Hall.
Общая значимость регрессионной модели может быть оценена с помощью критерия хи-квадрат, вычисляемого на основе логарифмов правдоподобия для подогнанной и нулевой моделей.
Для оценки адекватности подогнанной модели предоставляется большой выбор графических опций. В случае моделей пропорциональных интенсивностей пользователь может построить функции выживания для различных значений независимых переменных. Для экспоненциальной регрессионной модели есть возможность построения графиков зависимости остатков и экспоненциальной порядковой статистики, остатков и предсказанных с помощью регрессионного уравнения времен жизни, остатков и логарифмов наблюдаемых времен жизни. Для нормальной и логнормальной линейной регрессионной модели пользователь может воспроизвести на экране график зависимости наблюдаемых и подогнанных времена жизни, подогнанных времен жизни и остатков подгонки, а также нормальный вероятностный график остатков.
Альтернативные процедуры возможны для нецензурированных данных.
Если данные о продолжительности жизни (безотказной работы) нецензурированы, то применимо большинство непараметрических статистик. Для нецензурированных данных можно также использовать нелинейное оценивание, чтобы подогнать определенную регрессионную модель (включая пробит, логит и экспоненциальную модели) к данным.
Если продолжительность жизни или безотказной работы описывается бинарной переменной, то могут быть применены логит или пробит регрессионные модели.
Другой общий метод сравнения выживаемости в различных группах реализуется с помощью таблиц частот. Если времена жизни, или наработки до момента отказа, распределены по нескольким временным интервалам, может быть использована общая логлинейная модель.
Пример 1. Таблицы времен жизни
В этом примере мы рассчитаем таблицу времен жизни, оценим функцию выживания, плотность вероятности и функцию интенсивности для различных временных интервалов, а также найдем теоретическое распределение, наилучшим образом согласующееся с данными. Данные основаны на работе Crowley, J., & Ни, М., (1977) Covariance analysis of heart transplant survival data, Journal of the American Statistical Association, 72, p. 27-36.
В модуле Анализ выживаемости откройте файл Heart.sta.
Далее выберите Таблицы и распределения времен жизни из стартовой панели Анализ выживаемости и времен отказов.
Можно анализировать как исходный файл данных, так и сгруппированные данные. В данном случае мы анализируем исходные данные.
Нажмите кнопку Переменные и выберите 6 переменных в первом списке.
Первые три переменные — дата начала (например, дата операции), оставшиеся три переменные — дата наступления события.
Программа интерпретирует первую и четвертую переменные как месяцы, вторую и пятую — как дни, а третью и шестую — как год.
Заметим, что можно сразу ввести времена жизни (одна переменная в файле данных или даты в другом формате — 2 переменные).
Далее необходимо определить переменную Censored как индикатор цензурирования во втором списке.
Диалоговое окно Таблицы и распределения времен жизни будет теперь выглядеть так:
Поскольку были использованы коды по умолчанию для индикатора цензурирования (0 -полное, 1 — неполное), STATISTICA автоматически отображает Код для завершенных наблюдений и Код для неполных или цензурированных наблюдений.
Дополнительно можно определить для таблицы времен жизни Число интервалов или Ширину интервалов.
Процедура подгонки теоретического распределения к данным невозможна при наличии интервалов, не содержащих ни смертей (отказов), ни изъятых наблюдений.
Если вы хотите сделать подгонку, установите флажок Исправить интервалы, не содержащие смертей/отказов.
Если таблица времен жизни используется только в описательных целях и не предполагается подгонка распределения, то корректировку интервалов делать не нужно.
Оставив опции по умолчанию, нажмите ОК. После того как все наблюдения обработаны, откроется диалоговое окно Результаты для таблиц и распределений времен жизни.
Нажмите на кнопку Таблица времен жизни, чтобы отобразить на экране полную таблицу результатов времен жизни.
На рисунке показана часть полной таблицы жизни.
Можно подгонять к данным основные семейства распределений, используя обычный метод наименьших квадратов или две модификации метода взвешенных наименьших квадратов.
Чтобы выбрать наиболее подходящее семейство распределений, сначала рассмотрим модель с экспоненциальным распределением (выбрав позицию Экспоненциальная в поле Модель).
Оценка согласия проводится с помощью критерия хи-квадарт.
Нажмите кнопку Оценки параметров, чтобы посмотреть оценки для данного семейства распределений, а также значение критерия хи-квадрат.
Если критерий значим, делается заключение, что подогнанное распределение значимо расходится с наблюдаемыми данными. Поэтому мы отвергаем это семейство распределений и говорим, что оно не согласуется с данными.
Из таблицы результатов следует, что ни один метод подгонки не дает экспоненциального распределения удовлетворительного согласия. Тот же результат хорошо виден на графиках.
Нажмите кнопку График функции выживания. На приведенных ниже графиках ни одна из экспонент также не аппроксимирует наблюдаемую функцию выживания удовлетворительно. Видно, что оцененная функция выживания сильно отклоняется от аппроксимирующих функций выживания.
Можно просмотреть оценки параметров для различных семейств распределений. Вначале выберите соответствующее семейство из поля списка Модель, а затем нажмите кнопку Оценки параметров. Если проанализировать все эти семейства, можно сделать вывод, что только для семейства Вейбулла (см. главу Вероятностные распределения) нет значимого отличия от наблюдаемых значений при оценивании параметров по минимуму суммы взвешенных квадратов.
Ниже показаны графики функции выживания из семейства Вейбулла, подогнанные тремя разными способами.
Для третьего набора параметров (соответствующего Weight 3) имеется удовлетворительное согласие с данными. Хи-квадрат — критерий для этой ситуации — не дает значимого отклонения (р=0,56). Следовательно, можно сделать вывод, что распределение Вейбулла с этим набором параметров удовлетворительно описывает наблюдаемые времена жизни.
В заключение заметим, что модуль Анализ выживаемости STATISTICA позволяет анализировать также табулированные данные (для этого нужно выбрать опцию Таблица времен жизни в поле списка Входные данные).
Файл с табулированными данными должен содержать 3 переменные со следующей информацией:
1) нижняя граница временных интервалов;
2) число цензурированных или неполных наблюдений;
3) число отказов (число умерших в каждом временном интервале).
После выбора Таблиц времен жизни откроется диалоговое окно Таблицы и распределения времен жизни, в котором можно выбрать эти переменные.
Пример 2. Регрессионная модель Кокса
Файл данных Heart.sta содержит дополнительные переменные: возраст пациента во время трансплантации (переменная Возраст — Age) и медицинские характеристики: мера антигенной несовместимости (переменная Антиген — Antigen) и мера тканевой несовместимости (переменная Несовместимость — Mismatch).
Представляет интерес зависимость между переменными Возраст — Age, Антиген — Antigen и Несовместимость — Mismatch и временами жизни. Наиболее общей регрессионной моделью, не накладывающей ограничения на форму функции выживания, является модель пропорциональных интенсивностей Кокса. Рассмотрим, как можно оценить коэффициенты регрессии для этих трех независимых переменных для того, чтобы предсказать времена жизни с помощью модели пропорциональных интенсивностей Кокса.
Нажмите опцию Регрессионные модели на Стартовой панели, чтобы открыть диалоговое окно Регрессионные модели для цензурированных данных.
Чтобы выбрать переменные для анализа, нажмите кнопку переменные и задайте все времена жизни и цензурирующую переменную, как это было сделано ранее.
Необходимо также выбрать независимые переменные или регрессоры (Возраст — Age, Антиген — Antigen, Несовместимость — Mismatch).
Группирующую переменную в данном примере мы не отмечаем.
Теперь выберите коды для цензурирующей переменной. С помощью этих кодов STATISTICA разобьет данные на 2 группы: полные и неполные. По умолчанию STATISTICA использует следующий код: 0 = завершенное наблюдение, 1 = ценэури-рованное.
Если вы используете другой код, дважды щелкните по полю ввода Коды завершенного наблюдения и Коды цензурированного наблюдения и выберите коды из списка.
Диалоговое окно Регрессионные методы для цензурированных данных появится на экране:
Выберите в списке Модель позицию Регрессионная модель Кокса. Нажмите ОК и откройте диалоговое окно Регрессионная модель, оценивание.
Это диалоговое окно позволяет задать параметры процедуры оценивания.
Процедура оценивания максимизирует логарифм правдоподобия регрессионной модели с помощью метода Ньютона — Рафсона.
Алгоритм оценивание параметров является итеративным и начинается с некоторых начальных значений параметров (кнопка Начальные значения). Далее программа делает несколько итераций, последовательно приближаясь к оценкам неизвестных параметров. Разность между текущими оценками и оценками, полученными на предыдущем шаге, называется невязкой. Если невязка удовлетворяет критерию сходимости (см. поле Критерий сходимости), то процесс приближения завершается. Максимальное число итераций и критерий сходимости указываются в соответствующих полях.
Значения, предлагаемые программой по умолчанию, обычно приемлемы, поэтому просто нажмите ОК и начните процедуру оценивания.
С помощью этого диалогового окна можно наглядно проследить, как происходит процесс оценивания. В столбцах Параметры показаны оценки параметров на каждом таге.
После того как критерий сходимости будет выполнен, процедура оценивания останавливается.
Обычно процедура поиска быстро сходится, если приближения за заданное число итераций неудовлетворительны, программа запросит дополнительно некоторое количество итераций. Вы можете изменить начальные значения, используя, например, оценки параметров, полученные на предыдущем экспериментальном материале.
В данном примере наилучшие оценки параметров найдены, итеративная процедура сходится, поэтому предлагается нажать ОК, чтобы перейти в диалоговое окно Результаты.
Это диалоговое окно позволяет просмотреть результаты. Значение статистики критерия хи-квадрат для данной модели высокозначимо, поэтому можно заключить, что, по крайней мере, некоторые независимые переменные значимо действительно связаны с выживаемостью.
Нажмите кнопку Оценки параметров, чтобы увидеть оценки параметров и их стандартные ошибки.
Стандартные ошибки вычисляются как часть процедуры оценивания и по своей природе являются асимптотическими. Они вычисляются на основе частных производных второго порядка от логарифма функции правдоподобия. Это означает, что t-значения тоже должны рассматриваться только как приближенные. Обычно любая оценка параметра (регрессионной модели), которая по крайней мере в два раза превосходит свою стандартную ошибку (t>2,0), может рассматриваться как статистически значимая (на уровне р<0,05).
Электронная таблица С результатами также содержит статистику критерия Валъда для каждого коэффициента (см. книгу Рао С. Р. Линейные статистические методы и их применения). Из приведенной таблицы следует, что возраст пациента и тканевая несовместимость — наиболее важные предикторы для функции мгновенного риска.
Итак, значимые переменные в модели — AGE и MISMATCH. Рассмотрим графики функции выживания как функции независимых переменных. Пусть все независимые переменные равны своим средним значениям, тогда график функции выживания имеет вид (нажмите кнопку График выживаемости для средних):
Средние значения независимых переменных и стандартные ошибки можно посмотреть в таблице:
Зададим определенные значения предикторов. Мы имеем значимые переменные: AGE — возраст и MISMATCH — тканевая несовместимость. Увеличим возраст больного до 55 лет.
График функции выживания изменится и будет иметь вид:
В заключение заметим, что с помощью кнопки Редактор данных графика можно представить функцию выживания в численном виде:
Таким образом проводится регрессионный анализ в модуле Анализ выживаемости.
|
|
 |
