Анализ выживаемости в системе “Statistica”

МИНИСТЕРСТВО НАУКИ И ОБРАЗОВАНИЯ  РФ

 

МОСКОВСКИЙ  ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ЭКОНОМИКИ, СТАТИСТИКИ И ИНФОРМАТИКИ

 

Кафедра математической статистики и  эконометрики

 

 

 

 

 

 

 

 

 

КУРСОВАЯ  РАБОТА

на тему:

Анализ  выживаемости в системе “Statistica”

 

 

 

 

 

 

 

 

Выполнила:

Студентка группы ДЭС-401

Вакуленко Екатерина

Вариант 3

 

Преподаватель:

Звездина Н.В.

 

 

 

 

 

г. Москва, 
2011 г.

Содержание:

 

Введение 3

1.1 Аппроксимация  эмпирических данных теоретическим  распределением. 5

1.2 Оценки  Каплана–Майера 15

1.3 Сравнение  выживаемости в группах 18

1.4 Регрессионная  модель Кокса 20

Заключение 22

Список использованной литературы 23

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Введение

Методы анализа выживаемости интенсивно применяются в медицине, биологии, страховании и промышленности.

Одной из важных характеристик, описывающих  течение болезни, является продолжительность  жизни пациентов с момента  поступления в клинику или  после проведения операции. В принципе, для описания средних времен жизни и сравнения новой методики со старой можно использовать стандартные статистические методы. Однако рассматриваемые данные имеют специфику, которую следует учитывать. Дело в том, что в медицинской практике мы часто имеем дело с неполными данными.

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

Естественное желание исследователя  использовать все данные, т. е. анализировать  как полные времена жизни, так  и неполные, и не терять с трудом собранную информацию.

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

Целью данной курсовой работы является оценить функцию выживания после проведения операции на сердце в трех клиниках.

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

Объектом исследования является совокупность больных в трех клиниках.

Предметом исследования является продолжительность жизни больных, после того как им была сделана операция на сердце.

Информационной базой исследования являются данные  о выживаемости пациентов, перенесших операцию по трансплантации сердца в трех клиниках. Обработка статистической информации осуществлена с использованием программного обеспечения «Statistica».

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

 

 

 

 

 

 

 

1.Таблицы времён жизни

1.1. Исходные данные

Таблица 1

Данные  о пациентах, перенесших операцию на сердце

	MONTH_1	DAY_1	YEAR_1	MONTH_2	DAY_2	YEAR_2	CENSORED	AGE	ANTIGEN	MISMATCH	HOSPITAL
1	JANUARY	6	71	JANUARY	21	71	CENSORED	54	0	1,11	HILLVIEW
2	MAY	2	71	MAY	5	71	CENSORED	40	0	1,66	HILLVIEW
3	AUGUST	31	71	MAY	17	73	COMPLETE	51	0	1,32	HILLVIEW
4	AUGUST	22	71	OCTOBER	7	71	COMPLETE	42	0	0,61	ST_AND
5	SEPTEMBR	9	71	JANUARY	14	72	CENSORED	48	0	0,36	ST_AND
6	OCTOBER	5	71	DECEMBER	8	72	COMPLETE	54	0	1,89	ST_AND
7	OCTOBER	26	71	JULY	7	75	COMPLETE	54	0	0,87	BINER
8	NOVEMBER	22	71	AUGUST	29	72	COMPLETE	49	0	1,12	BINER
9	NOVEMBER	20	71	DECEMBER	13	71	CENSORED	56	0	2,05	HILLVIEW
10	FEBRUARY	15	72	FEBRUARY	25	72	COMPLETE	55	1	2,76	HILLVIEW
11	FEBRUARY	8	72	NOVEMBER	29	74	COMPLETE	43	0	1,13	BINER
12	MARCH	29	72	MAY	7	72	COMPLETE	42	0	1,38	HILLVIEW
13	APRIL	13	72	APRIL	13	74	COMPLETE	58	0	0,96	ST_AND
14	JULY	16	72	NOVEMBER	29	72	COMPLETE	52	1	1,62	ST_AND
15	MAY	22	72	APRIL	1	77	CENSORED	33	0	1,06	ST_AND
16	AUGUST	16	72	AUGUST	17	72	CENSORED	54	0	0,47	BINER
17	SEPTEMBR	3	72	DECEMBER	18	74	COMPLETE	44	0	1,58	BINER
18	SEPTEMBR	14	72	NOVEMBER	13	72	COMPLETE	64	0	0,69	HILLVIEW
19	JANUARY	16	73	APRIL	1	77	CENSORED	49	0	0,91	BINER
20	JANUARY	3	73	APRIL	1	77	CENSORED	40	0	0,38	HILLVIEW
21	MAY	19	73	JULY	12	73	COMPLETE	49	0	2,09	HILLVIEW
22	MAY	13	73	JUNE	29	73	COMPLETE	61	1	0,87	ST_AND
23	MAY	9	73	MAY	9	73	CENSORED	41	0	0,87	ST_AND
24	JULY	4	73	APRIL	1	77	CENSORED	48	0	0,75	BINER
25	OCTOBER	15	73	APRIL	1	77	CENSORED	45	0	0,98	BINER
26	JANUARY	5	74	FEBRUARY	18	74	CENSORED	36	0	0,00	ST_AND
27	JANUARY	11	74	OCTOBER	1	76	COMPLETE	48	0	0,81	BINER
28	FEBRUARY	22	74	APRIL	14	74	COMPLETE	47	0	1,38	HILLVIEW
29	MARCH	22	74	APRIL	1	77	CENSORED	36	0	1,35	HILLVIEW
30	APRIL	24	74	JANUARY	2	75	COMPLETE	48	1	1,08	HILLVIEW
31	AUGUST	18	74	OCTOBER	8	74	COMPLETE	52	0	1,51	ST_AND
32	NOVEMBER	8	74	APRIL	1	77	CENSORED	38	0	0,98	ST_AND
33	OCTOBER	13	74	AUGUST	30	75	COMPLETE	48	1	1,82	ST_AND
34	DECEMBER	15	74	APRIL	1	77	CENSORED	41	0	0,19	BINER
35	NOVEMBER	20	74	JANUARY	9	75	COMPLETE	49	0	0,66	BINER
36	JANUARY	7	75	APRIL	1	77	CENSORED	32	1	1,93	BINER
37	MARCH	4	75	SEPTEMBR	6	76	CENSORED	48	0	0,12	HILLVIEW
38	MARCH	17	75	MAY	22	75	COMPLETE	51	0	1,12	HILLVIEW
39	MAY	18	75	JANUARY	1	76	CENSORED	19	0	1,02	HILLVIEW
40	APRIL	9	75	JUNE	13	75	COMPLETE	45	1	1,68	ST_AND
41	JUNE	10	75	APRIL	1	77	CENSORED	48	0	1,20	ST_AND
42	JUNE	21	75	JULY	16	75	COMPLETE	53	1	1,68	ST_AND
43	AUGUST	20	75	APRIL	1	77	CENSORED	47	0	0,97	BINER
44	AUGUST	17	75	APRIL	1	77	CENSORED	26	1	1,46	BINER
45	OCTOBER	7	75	DECEMBER	9	75	COMPLETE	56	1	2,16	BINER
46	SEPTEMBR	22	75	OCTOBER	4	75	CENSORED	29	0	0,61	HILLVIEW
47	NOVEMBER	18	75	APRIL	1	77	CENSORED	52	1	1,70	HILLVIEW
48	MAY	31	76	APRIL	1	77	CENSORED	49	0	0,81	HILLVIEW
49	FEBRUARY	4	76	MARCH	5	76	COMPLETE	54	0	1,08	ST_AND
50	DECEMBER	31	75	APRIL	1	77	CENSORED	46	0	1,41	ST_AND
51	JANUARY	17	76	APRIL	1	77	CENSORED	52	1	1,94	ST_AND
52	FEBRUARY	24	76	APRIL	13	76	CENSORED	53	0	3,05	BINER
53	MARCH	7	76	DECEMBER	29	76	COMPLETE	42	0	0,60	BINER
54	MARCH	8	76	APRIL	1	77	CENSORED	48	1	1,44	BINER
55	MAY	19	76	JULY	8	76	COMPLETE	46	0	2,25	HILLVIEW
56	APRIL	27	76	APRIL	1	77	CENSORED	54	0	0,68	HILLVIEW
57	AUGUST	21	76	OCTOBER	28	76	COMPLETE	51	1	1,33	HILLVIEW
58	SEPTEMBR	12	76	OCTOBER	8	76	CENSORED	52	1	0,82	ST_AND
59	MARCH	2	77	APRIL	1	77	CENSORED	45	0	0,16	ST_AND
60	AUGUST	7	76	APRIL	1	77	CENSORED	47	0	0,33	ST_AND
61	SEPTEMBR	17	76	FEBRUARY	25	77	COMPLETE	43	0	1,20	BINER
62	OCTOBER	16	76	APRIL	1	77	CENSORED	26	0	0,46	BINER
63	DECEMBER	12	76	APRIL	1	77	CENSORED	23	1	1,78	BINER
64	MARCH	19	77	APRIL	1	77	CENSORED	28	1	0,77	HILLVIEW
65	MARCH	31	77	APRIL	1	77	CENSORED	35	0	0,67	ST_AND

 

В строках располагаются данные о каждом из прооперированных пациентов. В столбцах указаны даты начала наблюдения за пациентом (дата поступления в клинику/дата операции) – первые три переменные, даты окончания наблюдения (пациент выписался, и связь с ним была потеряна или умер) – последние три переменные. Программа интерпретирует первую и четвёртую переменные как месяцы, вторую и пятую – как дни, а третью и шестую – как год.

Так, например, из пятой строки видно, что пациенту под номером 5 была сделана  операция 9 сентября 1971, а выписался он 14 января 1972 года. Так как далее связь с этим пациентом была утеряна, то имеем неполное (цензурированное) наблюдение. Ему соответствует значение стоящей в седьмом столбце переменной – censored (цензурирован).

Следующая за ней переменная в столбце 8 (AGE) характеризует возраст пациентов.

Переменные в 9-м и 10-м столбцах содержат специальную медицинскую информацию об особенностях операции (ANTIGEN, MISMATCH).

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

Файл исходных данных содержит 65 наблюдений, т.е. данные о 65 пациентах трех клиник.

1.2. Построение таблиц  времени жизни

 

На основе данных таблиц времен жизни (таблиц смертности - в терминологии страхования) определяется ряд элементарных статистик, необходимых  для описания времени жизни пациентов (клиентов - в страховании).

В некоторых случаях времена  отказов (failure time) представляются в виде сгруппированных данных. Это объясняется тем, что во многих реальных исследованиях сложно оценить время отказов с достаточной точностью, однако можно определить, сколько отказов произошло или сколько наблюдений было цензурировано в течение определенного интервала времени. Такого рода данные называются таблицами времен жизни.

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

Количество интервалов на временной  оси пользователь может задать самостоятельно. В приведенной ниже таблице это число равно 12 (с учетом того, что стандартный период наблюдения за пациентом составляет обычно 1 год).

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

В модуле «Анализ выживаемости» предусмотрена возможность, обрабатывать как непосредственно файл первичных данных, так и сгруппированные данные. Ниже приведена таблица времен жизни, полученная в результате обработки исходной информации:

Таблица 2

Таблица времен жизни

 

Обратимся к интерпретации переменных, составляющих содержание полученной электронной  таблицы времен жизни (по столбцам):

• Номер интервала (Interval/Intno=Interval Number) для сгруппированных данных.
• Нижняя граница интервала (Interval Start)
• Середина интервала (Mid Point)
• Ширина интервала (Interval Width)
• Число в начале (Number Entering)

Число пациентов, которые были живы в начале рассматриваемого временного интервала.

• Число изъятых (Number Withdrwn) объектов

Число пациентов, связь с которыми была утеряна (т.е. изъятых из дальнейшего  рассмотрения после того, как они  выписались/перевелись из данной клиники). Эти объекты имеют метку цензурированные (censored) в файле исходных данных.

• Число изучаемых (Number Exposed) объектов

Число пациентов, которые были живы в начале рассматриваемого временного интервала, за вычетом половины от числа изъятых (цензурированных).

• Число умерших (Number Dying)

Число пациентов, умерших на данном отрезке времени (интервалe). Умершие объекты имеют метку complete.

• Доля умерших (Proportn Dead)

Отношение числа объектов, умерших в соответствующем  интервале, к общему числу объектов, попавших в этот интервал.

Рассмотрим  остальные столбцы построенной  электронной таблицы.

• Доля выживших (Proportn Survivng)

Получается  как разность между единицей и  долей умерших. Например, значение доли выживших в третьей строке получено как:

1 – 0,043478 = 0,956522.

• Кумулятивная доля выживших объектов  или функция выживания (Cum. Prop Survivng)

Это кумулятивная доля выживших к началу соответствующего временного интервала. Полученная доля, как функция от времени, представляет собой оценку функции выживания, то есть вероятность того, что пациент переживет данный период времени. Поскольку вероятности выживания считаются независимыми на разных интервалах, эта доля равна произведению долей выживших объектов по всем предыдущим интервалам. Если посмотреть на столбец (Cum. Prop Survivng), приведенной выше таблицы, то можно убедиться, что:

 и т. д. 

• Плотность вероятности (Problty Density)

Это оценка вероятности смерти (отказа) на соответствующем интервале. Получается в результате вычитания из значения функции выживания на данном интервале значения функции выживания на следующем интервале с последующим делением на ширину соответствующего интервала:

где - оценка вероятности смерти (отказа) в i-м интервале, - кумулятивная доля выживших объектов (функция выживания) к началу i-го интервала, - ширина i-го интервала.

Например, значение второй строки столбца Problty Density рассчитывается следующим образом:

.

На графике  оценки плотности вероятности видно (рис.1), что вероятность смерти в первые 160 дней после операции максимальна. Далее она резко падает.

Большие вероятности смерти расположены  также в интервалах от 161 до 322, от 968 до 1129, от 1614 до 1775 и др. (см. Таблицу 2).

 

 

Рис. 1. Функция плотности вероятности смерти

 

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

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

Рис. 2. Функция мгновенного риска

 

График функции мгновенного  риска наглядно свидетельствует  о том, что в первые дни после  операции на сердце риск смерти очень  велик, затем он падает, а спустя некоторое время  вновь начинает возрастать. Заметим, что именно функция риска используется исследователем в дальнейшем для прогностических целей.

Итак, исследователя интересует функция  риска, однако реально возможно получить лишь оценку функции риска. Поэтому  важна точность получаемых оценок. Понятно, что нельзя доверять оценкам, имеющим большую погрешность (например, если погрешность имеет тот же порядок, что и сами оценки). Поэтому следует внимательно просмотреть построенную таблицу и, если позволяет объем выборки, удалить из неё все «плохие» оценки, т.е. оценки с большой погрешностью. Это чрезвычайно важный принцип анализа данных!

С этой целью в таблице наряду с оценками приведены их стандартные ошибки для каждой из трех описанных выше функций (Std. Err. Cum. Proportion Surviving, Probability Density, Hazard Rate).

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

Медиана ожидаемого времени жизни (Median Life Exp).

По определению, медиана соответствует точке на временной оси, в которой кумулятивная функция выживания принимает значение 0,5. Например, из первой строчки таблицы столбца Median Life Exp видно, что пациент с вероятностью 0,5 будет жить 826 дней после операции. Если пациент пережил первый временной интервал (161 день после операции на сердце), то с вероятностью 0,5 он проживет еще 1008 дней, что соответствует второй строке таблицы и т.д. Другие процентили (например, 25-й и 75-й процентили или квартили) кумулятивной функции выживания вычисляются по такому же принципу. Следует иметь ввиду, что 50-й процентиль (медиана) кумулятивной функции выживания обычно не совпадает с точкой выживания 50% наблюдений данной выборки! Такое совпадение возможно только тогда, когда в течение прошедшего отрезка времени не было цензурированных наблюдений.  

1.3. Аппроксимация  эмпирических данных теоретическим  распределением.

 

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

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

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

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

Однако, такие преобразования приводят иногда к тому, что дисперсия остатков зависит от интервалов (то есть дисперсия различна на разных интервалах). Чтобы учесть это, в алгоритмах подгонки дополнительно используются оценки метода взвешенных наименьших квадратов двух типов. Программа по умолчанию сама выбирает те из них, которые производят лучшую аппроксимацию (на основе критерия c²). На практике оба подхода приводят к очень близким значениям оценок параметров. Возможно также для оценки параметров сгруппированных данных применение метода максимального правдоподобия.

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

Таблица 3

Процедура оценки параметров экспоненциального  распределения

 

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

Из приведенной таблицы  видно, что ни один из представленных методов оценивания (подгонки) не даёт для экспоненциального распределения  удовлетворительного согласия.

 

Рис. 3. Графическое представление эмпирической функции выживания и теоретических кривых экспоненциального распределения

Такую же картину можно  наблюдать на приведенном выше графике  эмпирической функции выживания  и кривых  экспоненциального распределения: ни одна из трех экспонент (соответствующих  трем различным алгоритмам оценивания) не аппроксимирует наблюдаемую функцию  выживания удовлетворительно. Эмпирическая функция выживания сильно отклоняется  от второй аппроксимирующей функции (Weight 2); согласованность с двумя другими теоретическими кривыми (Weight 1, Weight 3) несколько лучше, но при этом сохраняется значимое их отличие от «волнообразного» характера поведения рассматриваемой эмпирической функции. Поэтому необходимо продолжить поиск лучшей аппроксимации.

Теперь  рассмотрим модель с линейной интенсивностью (Linear Hazard).

Таблица 4

Процедура оценки параметров линейного распределения

 

 

Рис. 4. Графическое представление эмпирической функции выживания и теоретических кривых линейного распределения

Эмпирическая функция  выживания сильно отклоняется от второй аппроксимирующей функции (Weight 2); согласованность с двумя другими теоретическими кривыми (Weight 1, Weight 3) несколько лучше, но при этом сохраняется значимое их отличие от «волнообразного» характера поведения рассматриваемой эмпирической функции. Поэтому необходимо продолжить поиск лучшей аппроксимации.

Теперь  рассмотрим модель Гомпертца (Gompertz).

Таблица 5

Процедура оценки параметров распределения Гомпертца

 

Рис. 5. Графическое представление эмпирической функции выживания и теоретических кривых распределения Гомпертца

Эмпирическая функция  выживания сильно отклоняется от первой аппроксимирующей функции (Weight 1); согласованность с двумя другими теоретическими кривыми (Weight 2, Weight 3) лучше, но всё же необходимо продолжить поиск лучшей аппроксимации.

Наконец, рассмотрим модель  Вейбулла (Weibull).

Таблица 6

Процедура оценки параметров распределения Вейбулла

 

Сравнив оценки параметров для  остальных семейств распределений, предлагаемых системой «Statistica», можно сделать вывод, что только для распределения Вейбулла (при оценивании по минимуму суммы взвешенных квадратов, т.е. по третьему алгоритму Weight 3) отсутствует значимое отличие от наблюдаемых значений: c²-критерий не даёт значимого отклонения (p=0,557).  Следовательно, распределение Вейбулла с таким набором параметров описывает наблюдаемые времена жизни наилучшим образом. Однако стоит заметить, что исследователь ограничен в выборе  лишь из трех представленных наборов параметров.

Ниже представлены графики  функции выживания для семейства  распределений Вейбулла, подогнанные на основе трех алгоритмов (Weight1, Weight2, Weight3).

 

Рис. 6. Графическое представление эмпирической функции выживания и теоретических кривых распределения Вейбулла

 

В заключение отметим, что  имеется возможность анализировать  в качестве исходных  табулированные данные. Для этого нужно выбрать закладку Таблица времен жизни (Table of Survival Times) в диалоговом окне Таблицы и распределения времен жизни. В этом случае файл с табулированными данными должен содержать три переменные со следующей информацией:

а) нижняя граница временных интервалов;

б) количество цензурированных наблюдений;

в) число отказов (умерших) в каждом временном интервале.

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