CC BY
77
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2015 Управление, вычислительная техника и информатика № 2 (31)
УДК 519.2
DOI 10.17223/19988605/31/5
М.А. Семёнова, Е.В. Чимитова
КРИТЕРИИ ПРОВЕРКИ ГИПОТЕЗ О ПАРАМЕТРАХ ОБОБЩЕННЫХ МОДЕЛЕЙ ПРОПОРЦИОНАЛЬНЫХ ИНТЕНСИВНОСТЕЙ ПРИ НЕИЗВЕСТНОМ РАСПРЕДЕЛЕНИИ ВРЕМЕН ЖИЗНИ
Исследование выполнено при финансовой поддержке Министерства образования и науки Российской Федерации в рамках проектной части государственного задания в сфере научной деятельности № 2.541.2014K от 17.07.2014.
Рассматриваются вопросы построения вероятностных моделей пропорциональных интенсивностей Кокса и их обобщений - модели Ксая и SCE-модели в случае неизвестного базового распределения времен жизни. Разработан алгоритм оценивания регрессионных параметров и базовой функции риска с использованием функции частичного правдоподобия, описаны критерии проверки гипотез о параметрах моделей и критерий проверки выполнения предположения о пропорциональности рисков, предложенного М.С. Никулиным. Проведено исследование распределений статистик и мощности критериев отношения правдоподобия, Вальда и Никулина. Ключевые слова: модель пропорциональных интенсивностей; модель Ксая; SCE-модель; оценка максимального правдоподобия; предположение о пропорциональности рисков; критерий отношения правдоподобия; критерий Вальда; проверка адекватности.
В большинстве работ, посвященных исследованиям продолжительности жизни, при построении моделей выживаемости учитывается зависимость вероятности наступления системного события от значений ковариат [1-3]. При этом системное событие может представлять собой изменение определенных биохимических показателей, смерть тяжелобольного пациента, наступление ремиссии или рецидива заболевания при условии получения некоторого вида лечения или другие события. В качестве ковариат, в свою очередь, могут выступать как внутренние свойства объектов исследования (возраст, пол или наличие хронических заболеваний), так и условия проведения эксперимента (вид терапии или наличие вспомогательных видов лечения), которые могут оказывать влияние на время наступления исследуемого события.
Одной из первых моделей зависимости вероятности наступления системного события от ковариат является модель пропорциональных интенсивностей Кокса [4]. Данная модель получила широкую популярность благодаря двум неоспоримым преимуществам. Во-первых, модель пропорциональных ин-тенсивностей учитывает цензурированные наблюдения, наличие которых является типичным для задач анализа выживаемости. Действительно, не для всех наблюдаемых в исследовании объектов можно точно определить время наступления системного события, можно лишь утверждать, что системного события не произошло до некоторого момента времени, такие неполные наблюдения называют цензуриро-ванными. Во-вторых, для данной модели существует простая процедура непараметрического оценивания неизвестной базовой функции риска и регрессионных параметров [10].
Несмотря на все преимущества и насчитывающую десятилетия историю использования модели пропорциональных интенсивностей Кокса, вопросы о проверке предположения пропорциональности рисков и о корректности использования данной модели в случае непостоянного во времени отношения рисков наступления системного события при разных значениях ковариат остаются открытыми [5-7]. В [8] предложена модель Ксая, которая является обобщением модели пропорциональных интенсивно-стей Кокса и позволяет описывать пересекающиеся при разных значениях ковариат функции выживаемости, т.е. непропорциональные риски наступления системного события. Кроме этого, в [1, 9] приведена модель с пересечением функций выживаемости (SCE - simple cross-effect model), позволяющая описать не только пересекающиеся функции выживаемости, но и приближающиеся или отдаляющиеся друг от друга функции при разных значениях ковариат. Основной сложностью построения обобщенных мо-
делей является необходимость одновременного оценивания регрессионных параметров, в том числе обобщающих, и неизвестного базового распределения. Более того, при построении вероятностных моделей выживаемости требуется определение степени влияния ковариат на функцию выживаемости, для чего проверяется гипотеза о незначимости регрессионных параметров с использованием критерия отношения правдоподобия или критерия Вальда.
Таким образом, целью данной работы являются разработка алгоритма оценивания регрессионных параметров и базовой функции риска при построении полупараметрической модели Ксая и модели с пересечением функций выживаемости, исследование распределений статистик и мощности критериев проверки гипотезы о параметрах моделей и критерия проверки предположения о пропорциональности рисков.
1. Описание моделей
Пусть Тх - неотрицательная случайная величина, определяющая время наступления системного события, которое зависит от вектора ковариат х = [,х2,...,хт]т . Функция выживаемости определяется
соотношением Бх () = Р (Тх > t) = 1 - Ех (), а кумулятивная функция риска - выражением
í
Л х ^ ) = |х х (и )с1и = - 1п ( (t)) . Результаты эксперимента могу быть представлены в следующем виде:
о
(а, х1,51), ^2, х2,52),..., (^, х", 5„),
где п - объем выборки, х' - значение вектора ковариат для 7-го объекта, ti — ш1п(Т, С7} - время наступления системного события Т7 или момент цензурирования С7, 57 - индикатор цензурирования, 5' = 1{Т ^ С}, 7 — 1,...,п .
Модель пропорциональных интенсивностей определяется следующим соотношением [4]:
Лх (Р) = ехр(ртх)-Ло (0 , (1)
где р - т -мерный вектор параметров регрессии, Л0 ^) - базовая кумулятивная функция риска.
В соответствии с моделью пропорциональных интенсивностей отношение функций интенсивно-
X х-а (t) еХР (рТд )
сти при разных значениях ковариаты х не зависит от времени: — =--.
Хх=ь (t) ехр (рть)
Модель Ксая получена путем возведения базовой функции риска в степень ехр(уТх), позволяет описать непостоянное отношение функций интенсивности и имеет следующий вид:
Л х р, у) — ехр (рт х ){Л о « )}ехр(уТ х). (2)
Регрессионные параметры р и обобщающие регрессионные параметры у являются т -мерными.
При у — 0 данная модель является моделью пропорциональных интенсивностей, тогда как при у ф 0 функции интенсивности при разных значениях ковариат пересекаются [8]. Отношение интенсивностей
X х-а (0 = ехр (уТа ){Л о « ^Ч"
Х х—ь «"ехр (уТЬ ){Л 0 (t )}ехрИ-1
БСЕ-модель, предложенная в [9], позволяет получить не только пересекающиеся функции выживаемости, но и приближающиеся и отдаляющиеся друг от друга функции при разных значениях ковари-ат и может быть записана следующим образом:
Лх (;р, у) = (1 + ехр((Р + у)Т х)л0 (t))ехр(-у х) -1. (3)
монотонно.
Отношение функций интенсивности, соответствующих данной модели,
^x=fl (t) _ eXP + eXP ((Р + У)Т a)Л°
(t) aj
X x_* (t) /„ ТЛЛ чтл. / л\ехр(-УТь)
exp
монотонно. Обозначим отношение функций интенсив-
()( + ехр((р + у)Т ь)ло ^))
ности Xх=а (t)/Xх=ь (t) в точке t = 0 через с0. Тогда в зависимости от значения с0 функции выживаемости расходятся или приближаются друг к другу, если у < 0. Если же у > 0, то отношение функций интенсивности убывает с с0 до 0, т.е. функции выживаемости в этом случае при разных значениях ковариаты пересекаются на интервале (0, да).
2. Оценивание параметров
Для нахождения оценок неизвестных регрессионных параметров рассматриваемых моделей применяется метод максимального правдоподобия [10]. Логарифм функции правдоподобия для цензуриро-ванных справа данных в общем виде можно записать следующим образом:
ln L (0)_Ё8; ln
i _1
g (, Л о (t.), 0)
(4)
E g(,л° (t.),0)
=i, tj _
где 0 - вектор регрессионных параметров модели, функция g (x, ß)_ exp(ßrx) для модели пропорциональных интенсивностей, g (x, Л° (t), ß, у)_ exp ((ß + y)T x)(Л° (t))exp(y x) 1 для модели Ксая и
g(x,Л° (t),ß,у) _exp(ßTx)(l + exp((ß + y)T х)л° (t)) ( ) для SCE-модели.
Поскольку для модели пропорциональных интенсивностей функция g (x, ß) не зависит от базовой функции риска, то для нахождения оценок регрессионных параметров модели необходимо просто максимизировать логарифм функции правдоподобия: ß _ arg max lnL (ß). В [1°] предложена непараметриче-
ß
ская оценка базовой функции риска модели пропорциональных интенсивностей, которая может быть записана следующим образом:
Л° (t)_ E
i _1, t>t,
Sj E exp (ß T x)
/ j_1. tj щ
В [1] предложена итеративная процедура оценивания регрессионных параметров и базовой функции риска для семейства обобщенных моделей пропорциональных интенсивностей. Для модели Ксая и БСЕ-модели с зависящими от базовой функции риска функциями g (х, Л0 ^), р, у) данный алгоритм можно сформулировать следующим образом.
Чтобы оценить вектор регрессионных параметров 6 = |^рт, ут ^ обобщенных моделей пропорциональных интенсивностей, необходимо:
1. Положить к = 0, задать начальное приближение 60 = ^, ут ^ .
2. Оценить базовую функцию риска Л0 (/;66к), для этого:
а) упорядочить г различных значений полных наблюдений (5г- = 1) по возрастанию: w1 < <... < wг, задать равным числу наблюдений со значением {;
б) вычислить оценку функции риска
для
каждого
щ:
Л,
(о; 0к ) = о,
Л,
(0 к):
й1
£ ё (х, Ло (о; 0 к), 0 к)
1=1,
лЛо (щ+1;01) = лЛо (щ;0 к )+-
£ ё(х,Ло (щ;0к),0к) =1, ti >í,■ v 1
1 = 2,...,г -1.
3. Получить оценку 0 к+1 = а^шах1п Ь (0к).
4. Если
■'к+1 '
>е, задать к = к +1 и перейти на шаг 2, иначе - считать 0 к+1 найденной оценкой
регрессионных параметров.
В результате проведенного методами статистического моделирования исследования оценок регрессионных параметров модели пропорциональных интенсивностей, модели Ксая и 8СЕ-модели, полученных с использованием описанного алгоритма, показано, что с увеличением объема выборок уменьшается смещение и дисперсия получаемых оценок, тогда как при повышении степени цензурирования смещение и дисперсия оценок увеличиваются.
3. Критерии проверки гипотез о параметрах
В общем виде гипотеза о параметрах модели может быть записана как Но: 0 = 0о и проверена с помощью критерия отношения правдоподобия или критерия Вальда. Статистика критерия отношения правдоподобия имеет вид
ЬЯ = 2 (1п Ь (0)- 1п Ь (0о)) (5)
и асимптотически распределена по закону %2 с 5 степенями свободы, где 5 - количество оцениваемых параметров модели [13].
Статистика критерия Вальда для проверки гипотезы Но: 0 = 0о может быть записана следующим образом:
Ж = (0-0°)Т I(0)(0-0о), (6)
где I (0 ) = -
д21п Ь
(0)
д0 д0 1 1 У
- оценка информационной матрицы Фишера, 1, 1 = 1,...,5 . Статистика (6) также
асимптотически распределена по закону %2 с количеством степеней свободы, равным количеству оцениваемых параметров. Кроме этого, критерий Вальда позволяет проверять гипотезу о каждом параметре Но : 01 = 0о, ■ = 1,...,, в этом случае используется статистика
Ж =
(0- -0-)2 [ I-1(0)
(7)
где
1 (0)] ■ ■ _ диагональный элемент матрицы, обратной для матрицы I (0). Статистика (7) асимптотически распределена по закону %2 с одной степенью свободы [13].
Для получения оценки информационной матрицы Фишера вычислим производные по параметрам логарифма функции правдоподобия (3):
д 1п Ь
(0)
д0к
= £8г
к )
ё )к2 )
к )Ш - ё )дк Ь)'
, ё (^) = ё (х-, Л о (^), 0), к ) = £ ё (х1, Л о (^), 0),
д21п Ь
(0)
д0, д0
= £8г
1=1
1 д2 ё (^) 1 д2 к ) 1 дё )дк (^) + 1 дк )дё )
ё (^ )д0гд0к к (^ )д0гд0к ё2 (^) д0к д0г к2 ) д0к д0г
к,I = 1,...,5 .
й
В случае модели пропорциональных интенсивностей (1) выражения для вычисления первой и второй производных вспомогательной функции g (х, р) = ехр(ртх) имеют вид д(х Р)
ЗРк
• = "к ехр
(ртх) и
д 2 g (х, Р)
5Рг 5Рк
= х1хк ехр(ртх), соответственно, к,I = 1,...,т .
Далее обозначим для краткости через Л0 функцию Л0 ^). Тогда для модели Ксая (2) производные функции g (х,Л0,р,у) = ехр((р + у)т х)(Л0)ехр(у х=-1 по параметру р :
дg ( х, Л 0, Р, у)
ЗРк
д2 g (х, Л 0, Р, у)
дРг дРк
= хк ехр((р + у)т х)(Л0)ехр(утх)1, = ХЛ ехр ((р + у)т х )(Л 0 )ехр(ут х)-1,
по параметру у:
дg ( х, Л 0, Р, у)
дук
= хк ехр ((р + у)т х )(Л 0 )ехр(ут х)1 (1 + ехр (ут х =п Л0),
д2 g (х, Л 0, р, у)
дГ/дТк
смешанная производная по параметрам:
д2 g (х, Л 0, р, у)
= хх ехр ((р + у)т х)(Л0 )ехр(у х)-1 (1 + 3ехр (ут х)1п Л0 + ехр2 (ут х)1п2 Л0),
дуг дрк
= хх ехр ((р + у)т х)(Л0 )ехр(у х) 1 (1 + ехр (ут х)1пЛ0), к, I
= 1,...,т .
Для 8СЕ-модели (3) производные функции g(х,Л0,р,у) = ехр(ртх=1 + ехр((р + у)т х)л^ параметрам р и у приведены ниже:
ехр(-ут х)-1
по
дg (х, Л 0, р, у) дрк
= хк ехр
(рт х) (1 + ехр (рт х) Л0) (1 + ехр ((р + у)т х) Л0)
ехр(-ут х)-2
д g (x, Ло, р, у= = хл ехр ((Тх)( + ехр ((р + у=х)л0) у ^ 3 (1 + ехр(ртх)л0=1 + ехр ((р + у)т х— )-
дрг дрк
+ехр
(рт х) Л0 (2 + ехр (рт х) Л0- 3ехр (ут х — 2 ехр ((р + у)т х) Л0)
ехр (рт х)(1 - ехр (ут х= Л0 (ехр (-ут х) + ехр (рт х) Л0) 1п (1 + ехр ((р + у)т х) Л0)
= „ ехр(х=1+ехр((р + у)т х)л0)
/пт \Л . . ..//п . чт \ехр(-утх)-3
дТгдук
- = хх ехр (рт х) (^ + ехр ((р + у )т х) Л0) (ехр (ут х) -1-1(1 - ехр (у т х)) ехр (рт х) Л,
+ (ехр (-ут х) + ехр (рт х) Л0) 1п (1 + ехр ((р + у)т х) Л0 )| + + (1 + ехр ((р + у)т х) Л0) (1 + ехр ((р + у )Т х) Л0) - ехр ((р + у)Т х) Л0 -11
д g ду^Лр0^Р'^ = хЛ ехр(ртх)-1 + ехр(ртх-Л0-(1 + ехр((р + у)т х-Л0)СКр( ' х) 3 [ехр(ртх)
Л0 -
-2 ехр ((р + у= х=Л0 + (ехр (-ут х) + ехр (рт х =Л0 ) 1п (1 + ехр ((р + у= х=Л(
к, I = 1,..., т.
Как правило, при построении моделей зависимости функции выживаемости от ковариат возникает необходимость проверки гипотезы о равенстве нулю каждого параметра модели. Если гипотеза от-
вергается, то соответствующая данному параметру ковариата считается значимой. Важно отметить, что при построении обобщенных моделей значимость ковариаты может быть установлена на основании результатов проверки гипотезы о равенстве нулю как регрессионного параметра р, так и обобщающего параметра у, соответствующих данной ковариате.
В результате исследования методами статистического моделирования распределений О (| Но) статистик критериев отношения правдоподобия и критерия Вальда проверки гипотезы о параметрах модели пропорциональных интенсивностей и модели Ксая показано, что с увеличением объема выборок расстоя-
2
ние между эмпирическими распределениями статистик и соответствующим предельным % -распределением сокращается, распределения статистик не зависят от значений ковариат и регрессионных параметров. В случае проверки гипотезы о незначимости параметров р и у модели с пересечением функций выживаемости распределения статистики критерия Вальда оказываются далекими от предельного %2 -распределения даже при больших объемах выборок.
Не менее важной задачей, решаемой с использованием критериев проверки гипотез о параметрах обобщенных моделей, является проверка гипотезы о согласии с моделью пропорциональных интенсив-ностей Кокса. Для этого формулируется гипотеза о равенстве нулю только обобщающего параметра у .
Для исследования распределений статистик критериев проверки гипотезы Но: у = о о равенстве нулю обобщающего параметра 8СЕ-модели согласно модели пропорциональных интенсивностей Кокса с регрессионным параметром р = о,5 и базовым экспоненциальным распределением моделировались
цензурированные выборки объемом п = 1оо со скалярной ковариатой х = {о,1} , объем моделирования N = 1оооо . На рис. 1 и 2 представлены распределения статистик критерия отношения правдоподобия и критерия Вальда для разных степеней цензурирования, а также соответствующее предельное %2 -распределение.
С{1Я\Н0)
1.00 ■■
о.оо
¿Я
0.00
1.20
2.40
3.60
4.80
6.00
Рис 1. Эмпирические функции распределения статистики критерия отношения правдоподобия при проверке гипотезы о незначимости обобщающего параметра модели с пересечением функций выживаемости
Рис 2. Эмпирические функции распределения статистики критерия Вальда при проверке гипотезы о незначимости обобщающего параметра модели с пересечением функций выживаемости
На рис. 1 видно, что в случае полных выборок (без цензурирования) распределение статистики критерия отношения правдоподобия О (£ | Н0) близко к предельному %2 -распределению уже при
п = 100. Однако с ростом степени цензурирования расстояние от эмпирических распределений статистики данного критерия до предельного закона увеличивается. В свою очередь распределения статистики критерия Вальда значительно отклоняются от предельного %2 -распределения даже в случае проверки гипотезы по полным выборкам. Поэтому при проверке гипотезы о параметрах 8СЕ-модели по критерию Вальда использование предельного распределения для вычисления достигнутого уровня значимости может привести к неверному выводу. Данный факт является существенным недостатком критерия Вальда в сравнении с критерием отношения правдоподобия.
4. Критерий проверки предположения о пропорциональности рисков
Существует ряд методов проверки предположения о пропорциональности рисков при построении модели пропорциональных интенсивностей: графические методы, критерии, основанные на остатках, критерии Вальда и отношения правдоподобия для проверки незначимости добавленных в модель зависимых от времени ковариат и критерии против определенных конкурирующих гипотез. Одним из графических методов является сравнение графиков оценки Каплана-Мейера распределения остатков Кок-са-Снелла [14] с функцией стандартного экспоненциального распределения. Другой графический метод основан на сравнении наблюдаемых и ожидаемых кривых выживаемости: если они достаточно близки, то предположение пропорциональности рисков выполняется. Однако все графические методы основаны на субъективной оценке и могут быть использованы лишь для предварительной оценки адекватности модели.
Применение критерия проверки гипотезы о согласии при проверке предположения пропорциональности рисков является более объективным методом, чем графические методы, описанные ранее, так как позволяет получить значение статистик и достигнутые уровни значимости р. В существующей ли-
тературе рассматривается несколько критериев проверки предположения пропорциональности рисков. Модифицированный критерий, основанный на остатках Шонфельда, предложен в [7]. Для цензуриро-ванных данных в [5, 6] разработаны статистики проверки гипотезы о согласии с моделью пропорциональных интенсивностей. Однако применение данных критериев согласия для проверки предположения о пропорциональности рисков связано с необходимостью моделирования неизвестных распределений статистик [11] и / или с идентификацией базового распределения [12], что требует значительных вычислительных и временных ресурсов.
В настоящей работе рассматривается критерий, предложенный М.С. Никулиным в [9], для проверки гипотезы
Н0: Лх (/;р)= ехр(ртх)ло (/)
о модели пропорциональных интенсивностей против конкурирующей гипотезы о модели с пересечением функций выживаемости
Н1 : Л х (;р, у) = (1 + ехр ((р + у)т х) Л о (/^^ Х) -1.
Статистика критерия проверки гипотезы Н0 против гипотезы Н1 может быть записана в векторной форме следующим образом:
Т " (8)
Т = п-1и т Ъ-1и,
где и = [и„..,ит]Т , ик = £8.
- хк1п
(1 + ехР (р Т *))-Ц
к = 1,...,т , т - размерность вектора парамет-
ров модели пропорциональных интенсивностей; Ъ = Е** - Е*Е0лТ - ковариационная матрица вектора и,
Ео = - £8г п. =1
Ш - ЕЕТ
Е* = - £8г
п г=1
М) - ЕЕ т
¿о )
Е** = -£ §г'
п г=1
т - ееЕ т ¿о )
Ек =
¿1 к)
Ек =
¿1 к)
¿о )' к ¿о )'
¿о )= £ ехр (рт х]), ¿1 )= £ хк ехр (рт х*), ¿2 )= £ -X (х*)т ехр (рт х*),
*=1, щ
*=1, щ
*=1, о >и т
¿1 (и )= £ - хк ехр (рт х* )1п (1 + ехр (¡3т х* )Ло), ¿2 )= £ - х* (х*) ехр (рт х* )1п (1 + ехр (¡3т х* ),
¿2 (и )= £ -х1 (х*) ехр (рт х* )1п2 (1 + ехр (рт х* )Ло).
При справедливости проверяемой гипотезы статистика (8) подчиняется в пределе при п %2 -распределению с числом степеней свободы т .
На рис. 3 представлены распределения статистики (8) при объемах выборок п = 5о, 1оо, 2оо, зоо в
случае проверки гипотезы Но: Лх ^;р) = ехр(ртх)ло (/), компоненты вектора ковариат х1 ={о,1} и х2 = {о, 1,2,3}, значения параметров р1 = о,3 и р2 = о, 6.
Как видно из рис. 3, с увеличением объема выборок расстояние между эмпирическим и предельным распределением статистики (8) сокращается. Также было показано, что размерность вектора кова-риат и, следовательно, количество оцениваемых параметров влияют не только на число степеней свободы соответствующего предельного распределения, но и на близость эмпирического распределения статистики к теоретическому - чем меньше ковариат, тем ближе смоделированное распределение статистики к соответствующему предельному при том же объеме выборки.
Сравним мощность критерия проверки предположения о пропорциональности рисков, предложенного М.С. Никулиным, с мощностью критериев проверки гипотез о параметрах. Действительно, гипотеза о незначимости обобщающих параметров Но: у = о модели с пересечением функций выживаемости является гипотезой о выполнении предположения пропорциональности рисков, так как в случае
=1
у = 0 модель (3) является моделью пропорциональных интенсивностей. Рассмотрим следующие конкурирующие гипотезы:
1) Н1 : Лх (;р, у) = (1 + ехр ((р + у)т х) Ло ()
2) Н2 : Лх (;р,у) = (1 + ехр ((р + у)т х) Ло ())'
,С(Т\Н0)
-1, у = 0,5;
ехр-у х
-1, у = -0,5.
1.00--
0.80 --
0.60 --
0.40 --
0.20 --
0.00
0.00
2.00
4.00
6.00
8.00
10.00
Рис 3. Эмпирические функции распределения статистики Т критерия проверки предположения о пропорциональности рисков
Функции выживаемости, соответствующие гипотезам Н0, Н\ и Н2 при р = 0,5 для разных значений ковариаты х = {0,1}, представлены на рис. 4-6 соответственно.
Рис. 4. Функции выживаемости согласно Н 0
1.0
0.8
0.4
0.2
О.О
0.0
100
200
300
400
500
Рис 5. Функции выживаемости согласно Н
Рис 6. Функции выживаемости согласно Н
На рис. 5 видно, что при разных значениях ковариаты функции выживаемости, соответствующие модели с обобщающим параметром у = 0,5, пересекаются. Расстояние между функциями выживаемости
на рис. 6 при справедливости гипотезы Н12 больше аналогичного расстояния на рис. 4 при справедливости нулевой гипотезы о модели пропорциональных интенсивностей. В таблице представлены полученные оценки мощности рассматриваемых критериев для объемов выборок п = 100, 200 при заданном уровне значимости а = 0,1, объем моделирования N = 10000 .
Оценки мощности критериев
Критерий Н1 Н2
п = 100 п = 200 п = 100 п = 200
ЬК 0,40 0,57 0,14 0,34
W 0,01 0,20 0,27 0,45
Т 0,36 0,57 0,21 0,38
В случае конкурирующей гипотезы Н, соответствующей 8СЕ-модели с пересекающимися функциями выживаемости, наиболее мощным из рассмотренных критериев оказался критерий отношения правдоподобия. При объеме моделируемых выборок п = 100 мощность критерия Вальда оказалась меньше заданного уровня значимости, при увеличенном объеме выборок п = 200 критерий Вальда показал наименьшую мощность в сравнении с критерием отношения правдоподобия и критерием Никулина. В случае конкурирующей гипотезы Н2 с расходящимися функциями выживаемости критерии со статистиками ЬК и Т уступают по мощности критерию со статистикой W на всех рассмотренных объемах выборок, при этом мощность критерия Никулина больше мощности критерия отношения правдоподобия. Однако, учитывая смещенность критерия Вальда на паре конкурирующих гипотез Н0 и Н/, а также существенное отличие эмпирических распределений статистик О (£ | Н0) от соответствующего предель-2
ного х -распределения, для проверки предположения пропорциональности рисков рекомендуется применять критерий отношения правдоподобия со статистикой (5) и критерий Никулина со статистикой (8).
Заключение
В настоящей работе сформулирован итеративный алгоритм оценивания регрессионных параметров и базовой функции риска полупараметрической модели Ксая и 8СЕ-модели. В результате проведенного методами компьютерного моделирования исследования статистических свойств оценок параметров моделей показано, что с увеличением объема выборок уменьшаются смещение и дисперсия по-
лучаемых оценок, тогда как при повышении степени цензурирования смещение и дисперсия оценок увеличиваются.
На основании результатов исследования распределений статистик и мощности критериев отношения правдоподобия, Вальда и Никулина целесообразно рекомендовать использование критерия отношения правдоподобия и критерия Никулина для проверки гипотезы о выполнении предположения пропорциональности интенсивностей при построении полупараметрической модели пропорциональных интенсивностей Кокса. Применение критерия Вальда для проверки гипотезы как о значимости влияния ковариат, так и о согласии с моделью Кокса, сопряжено с возможными ошибками в определении достигнутого уровня значимости ввиду существенных отличий реальных распределений статистики при конечных объемах выборок от соответствующих предельных законов.
ЛИТЕРАТУРА
1. Bagdonavicius V., Nikulin M. Accelerated life models: modeling and statistical analysis. Boca Raton, Florida : Chapman &
Hall/CRC, 2002. 334 p.
2. Lee E., Wang J. Statistical methods for survival data analysis. 3rd. New Jersey : John Wiley & Sons, Inc., Hoboken, 2003. 534 p.
3. Semenova M., Bitukov A. Parametric models in the analysis of patients with multiple myeloma // Proceedings of the International
Workshop "Applied methods of statistical analysis. Applications in survival analysis, reliability and quality control". Novosibirsk : NSTU publisher, 2013. P. 250-256.
4. Cox D.R., Roy J. Regression models and life tables (with Discussion) // Journal of the Royal Statistical Society. 1972. Series B. V. 34.
P. 187-220.
5. Lin D.Y. Goodness-of-fit analysis for the Cox regression model based on a class of parameter estimators // JASA. 1991. V. 86.
P. 725-728.
6. Grambsch P., Therneau T.M. Proportional hazards tests and diagnostics based on weighted residuals // Biometrika. 1994. V. 81.
P. 515-526.
7. Harrell F.E. Regression modeling strategies with applications to linear models, logistic regression, and survival analysis. N.Y. :
Springer, 2002. 572 p.
8. Hsieh F. On heteroscedastic hazards regression models: theory and application // Journal of the Royal Statistical Society. 2001. Series
B. V. 63. P. 63-79.
9. Bagdonavicus V., Levuliene R., Nikulin M. Modeling and testing of presence of hazard rates crossing under censoring // Comm. in
Stat. Sim. and Comp. 2012. V. 41. P. 980-991.
10. Breslow N.E. Analysis of survival data under the proportional hazards model // International Statistical Review. 1975. V. 43. P. 4557.
11. Чимитова Е.В., Ведерникова М.А. Проверка адекватности модели пропорциональных интенсивностей Кокса по случайно цензурированным выборкам // Сборник научных трудов НГТУ. 2010. № 4(62). С. 103-108.
12. Balakrishnan N., Chimitova E., Galanova N., VedernikovaM. Testing goodness-of-fit of parametric AFT and PH models with residuals // Comm. in Stat. Sim. and Comp. 2013. V. 42. P. 1352-1367.
13. Айвазян С.А., Енюков И.С., Мешалкин Л.Д. Прикладная статистика: Основы моделирования и первичная обработка данных // Финансы и статистика. 1983. 471 с.
14. Kalbfleisch J.D., Prentice R.L. The statistical analysis of failure time data. N.Y. : John Wiley & Sons, Inc., 1980.
Семёнова Мария Александровна. E-mail: chimitova@corp.nstu.ru
Чимитова Екатерина Владимировна, канд. техн. наук, доцент. E-mail: vedernikova.m.a@gmail.com Новосибирский государственный технический университет
Поступила в редакцию 2 декабря 2014 г.
Semenova Maria A., Chimitova Ekaterina V. (Novosibirsk State Technical University, Russian Federation). Testing hypothesis of parameters of generalized proportional hazards models under unknown lifetime distribution. Keywords: proportional hazards model; Hsieh model; SCE-model; maximum likelihood estimation; proportional hazards assumption; likelihood ratio test; Wald test; goodness-of-fit testing
DOI 10.17223/19988605/31/5
The paper deals with the construction of the Cox proportional hazards model and its generalizations. The considered generalizations of the proportional hazards model are the Hsieh model and the simple cross-effect model (SCE-model), which allow decreasing, increasing or non-monotonic behavior of the ratio of hazard rate functions. The algorithm of the estimation of regression parameters and unknown baseline distribution for the generalized models is developed by using the partial likelihood function. The research on statistical properties of estimates carried out with computer simulations, has shown that the bias and the variance of obtained estimates decrease with the sample size growth. However, the bias and the variance of obtained estimates increase with the censoring degree growth.
The likelihood ratio test and the Wald test are used for testing hypothesis about parameters of considered models. In this paper, the expressions of elements of the matrix of the second partial derivatives by regression parameters of the Hsieh model and SCE-model have been obtained. In the case of testing hypothesis about parameters of proportional hazards model and the Hsieh model, the distributions G (S | H0) of likelihood ratio test statistic and the Wald statistic are independent of covariates values or regression parameters
values. The difference between the simulated test statistic distributions and the corresponding % -distributions decreases with the sample size growth. The dimension of the covariate vector and the number of estimated parameters affect not only the number of degrees of freedom of the limiting distribution, but also the closeness of simulated test statistic distributions to the theoretical distributions: the smaller the dimension of the covariate vector the smaller the difference between empirical and limiting distributions for the same sample size. In the case of testing hypothesis of insignificance of parameters P and y of SCE-model, the statistic distributions G (S
IH0) of
2
the Wald test are not close to the corresponding limiting % -distributions even for the large sample sizes.
Basing on the obtained results of research on statistic distributions and the power of considered tests, it is advisable to use the test proposed by M.S. Nikulin and likelihood ratio test for checking proportional hazard assumption against the competing hypothesis corresponding to the SCE-model. The application of the Wald test can result in inaccurate computation of p-value because empirical statistic distributions significantly differ from corresponding limiting distributions.
REFERENCES
1. Bagdonavicius, V. & Nikulin M. (2002) Accelerated life models: modeling and statistical analysis. Boca Raton, Florida: Chapman &
Hall/CRC.
2. Lee, E. & Wang, J. (2003) Statistical methods for survival data analysis. Hoboken, New Jersey: John Wiley & Sons, Inc.
3. Semenova, M. & Bitukov, A. (2013) Parametric models in the analysis of patients with multiple myeloma. Proc. of the International
Workshop "Applied methods of statistical analysis. Applications in survival analysis, reliability and quality control'. Novosibirsk: NSTU. pp. 250-256.
4. Cox, D.R. & Roy, J. (1972) Regression models and life tables (with Discussion). Journal of the Royal Statistical Society. Series B. 34.
pp. 187-220.
5. Lin, D.Y. (1991) Goodness-of-fit analysis for the Cox regression model based on a class of parameter estimators. JASA. 86. pp. 725-
728. DOI: 10.1080/01621459.1991.10475101
6. Grambsch, P. & Therneau, T.M. (1994) Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 81.
pp. 515-526. DOI: 10.1093/biomet/81.3.515
7. Harrell, F.E. (2002) Regression modeling strategies with applications to linear models, logistic regression, and survival analysis.
New York: Springer.
8. Hsieh, F. (2001) On heteroscedastic hazards regression models: theory and application. Journal of the Royal Statistical Society. Series
B. 63. pp. 63-79. DOI: http://dx.doi.org/10.1111/1467-9868.00276
9. Bagdonavicus, V., Levuliene, R. & Nikulin, M. (2012) Modeling and testing of presence of hazard rates crossing under censoring.
Communication in Statistics - Simulation and Computation. 41. pp. 980-991. DOI: 10.1080/03610918.2012.625758
10. Breslow, N.E. (1975) Analysis of survival data under the proportional hazards model. International Statistical Review. 43. pp. 4557.
11. Chimitova, E.V. & Vedernikova, M.A. (2010) Testing goodness-of-fit hypothesis with the proportional hazard Cox model by independent censored samples. Sbornik nauchnykh trudov NGTU. 4(62). 103-108. (In Russian).
12. Balakrishnan, N., Chimitova, E., Galanova, N. & Vedernikova, M. (2013) Testing goodness-of-fit of parametric AFT and PH models with residuals. Communication in Statistics - Simulation and Computation. 42. pp. 1352-1367. DOI: 10.1080/03610918.2012.659824
13. Ayvazyan, S.A., Enyukov, I.S. & Meshalkin, L.D. (1983) Prikladnaya statistika: Osnovy modelirovaniya ipervichnaya obrabotka dannykh [Applied statistic: Basis of modeling and data mining]. Moscow: Finansy i statistika.
14. Kalbfleisch, J.D. & Prentice, R.L. (1980) The statistical analysis of failure time data. New York: John Wiley & Sons, Inc.
