Testing hypothesis of parameters of generalized proportional hazards models under unknown lifetime distribution.pdf В большинстве работ, посвященных исследованиям продолжительности жизни, при построении моделей выживаемости учитывается зависимость вероятности наступления системного события от значений ковариат [1-3]. При этом системное событие может представлять собой изменение определенных биохимических показателей, смерть тяжелобольного пациента, наступление ремиссии или рецидива заболевания при условии получения некоторого вида лечения или другие события. В качестве ковариат, в свою очередь, могут выступать как внутренние свойства объектов исследования (возраст, пол или наличие хронических заболеваний), так и условия проведения эксперимента (вид терапии или наличие вспомогательных видов лечения), которые могут оказывать влияние на время наступления исследуемого события. Одной из первых моделей зависимости вероятности наступления системного события от ковариат является модель пропорциональных интенсивностей Кокса [4]. Данная модель получила широкую популярность благодаря двум неоспоримым преимуществам. Во-первых, модель пропорциональных ин-тенсивностей учитывает цензурированные наблюдения, наличие которых является типичным для задач анализа выживаемости. Действительно, не для всех наблюдаемых в исследовании объектов можно точно определить время наступления системного события, можно лишь утверждать, что системного события не произошло до некоторого момента времени, такие неполные наблюдения называют цензуриро-ванными. Во-вторых, для данной модели существует простая процедура непараметрического оценивания неизвестной базовой функции риска и регрессионных параметров [10]. Несмотря на все преимущества и насчитывающую десятилетия историю использования модели пропорциональных интенсивностей Кокса, вопросы о проверке предположения пропорциональности рисков и о корректности использования данной модели в случае непостоянного во времени отношения рисков наступления системного события при разных значениях ковариат остаются открытыми [5-7]. В [8] предложена модель Ксая, которая является обобщением модели пропорциональных интенсивно-стей Кокса и позволяет описывать пересекающиеся при разных значениях ковариат функции выживаемости, т.е. непропорциональные риски наступления системного события. Кроме этого, в [1, 9] приведена модель с пересечением функций выживаемости (SCE - simple cross-effect model), позволяющая описать не только пересекающиеся функции выживаемости, но и приближающиеся или отдаляющиеся друг от друга функции при разных значениях ковариат. Основной сложностью построения обобщенных моделей является необходимость одновременного оценивания регрессионных параметров, в том числе обобщающих, и неизвестного базового распределения. Более того, при построении вероятностных моделей выживаемости требуется определение степени влияния ковариат на функцию выживаемости, для чего проверяется гипотеза о незначимости регрессионных параметров с использованием критерия отношения правдоподобия или критерия Вальда. Таким образом, целью данной работы являются разработка алгоритма оценивания регрессионных параметров и базовой функции риска при построении полупараметрической модели Ксая и модели с пересечением функций выживаемости, исследование распределений статистик и мощности критериев проверки гипотезы о параметрах моделей и критерия проверки предположения о пропорциональности рисков. 1. Описание моделей Пусть TX - неотрицательная случайная величина, определяющая время наступления системного события, которое зависит от вектора ковариат х = \,x2,...,xm]т . Функция выживаемости определяется соотношением SX (t) = P (TX > t) = 1 - FX ((), а кумулятивная функция риска - выражением t Лх (t) = JXх (u)du = -In(Sx (t)) . Результаты эксперимента могу быть представлены в следующем виде: о (tl, X1, Si ), (t2, X2,52),..., (t„, Xn, S„ ) , где n - объем выборки, X - значение вектора ковариат для /-го объекта, t/ = min{T, C/} - время наступления системного события T или момент цензурирования C/, 5/ - индикатор цензурирования, 5/ = 1{T ^ C}, i = 1,...,n . Модель пропорциональных интенсивностей определяется следующим соотношением [4]: Лх (t;Р) = exp(ртх)-Ло (t), (1) где р - m -мерный вектор параметров регрессии, Л0 (t) - базовая кумулятивная функция риска. В соответствии с моделью пропорциональных интенсивностей отношение функций интенсивно- Xх=а (t) exp (pTfl) сти при разных значениях ковариаты х не зависит от времени: - =-т-f. Xx=fe (t) exp (ртй) Модель Ксая получена путем возведения базовой функции риска в степень exp(yTx), позволяет описать непостоянное отношение функций интенсивности и имеет следующий вид: Л х (t; Р, y) = exp (рт х ){Л 0 (t )}exp(yT х). (2) Регрессионные параметры р и обобщающие регрессионные параметры у являются m -мерными. При y = 0 данная модель является моделью пропорциональных интенсивностей, тогда как при y ф 0 функции интенсивности при разных значениях ковариат пересекаются [8]. Отношение интенсивностей X х=a (t) = exp (yTa ){Л 0 (t )}exp(yTfl)-1 Xx=*(t^exp(yTb){лЛ0 (t)}exp(yTb)-1 SCE-модель, предложенная в [9], позволяет получить не только пересекающиеся функции выживаемости, но и приближающиеся и отдаляющиеся друг от друга функции при разных значениях ковари-ат и может быть записана следующим образом: Лх (t;P,y) = (1 + exp((P + y)T х)л0 (t))exp(-y X)-1. (3) монотонно. Отношение функций интенсивности, соответствующих данной модели, Xx_a(t) _ exp(pT°+ exp((P + y)T a)ло (t ))) aj X x_b (t) /„ т,\Л ч'Л. / лН^) exp монотонно. Обозначим отношение функций интенсив- ()( + exp((P + y)t б)лq (t))e ности Xx_a (t)/Xx_b (t) в точке t _ Q через cQ. Тогда в зависимости от значения cQ функции выживаемости расходятся или приближаются друг к другу, если y < 0. Если же y > 0, то отношение функций интенсивности убывает с с0 до 0, т.е. функции выживаемости в этом случае при разных значениях ковариаты пересекаются на интервале (0, да). 2. Оценивание параметров Для нахождения оценок неизвестных регрессионных параметров рассматриваемых моделей применяется метод максимального правдоподобия [10]. Логарифм функции правдоподобия для цензуриро-ванных справа данных в общем виде можно записать следующим образом: (4) g (, л о (ti), е) ln L (9)_Ё8г ln i _1 Е g(,ло (t.),е) =1, tj ^ _ где е - вектор регрессионных параметров модели, функция g (x, p)_ exp(prx) для модели пропорциональных интенсивностей, g (x, лQ (t), p, y)_ exp ((P + y)t x)(лQ (t))exp(Y x) 1 для модели Ксая и g(x,лQ (t),p,y) _exp(pTx)(i + exp((P + y)t x)лQ (t)) ( ) для SCE-модели. Поскольку для модели пропорциональных интенсивностей функция g (x, p) не зависит от базовой функции риска, то для нахождения оценок регрессионных параметров модели необходимо просто максимизировать логарифм функции правдоподобия: p _ argmaxlnL (p). В [10] предложена непараметричер ская оценка базовой функции риска модели пропорциональных интенсивностей, которая может быть записана следующим образом: 87 ЕЕ exp (( т x) / j_1, 'j щ лло (t)_ Е i _1, t>t, В [1] предложена итеративная процедура оценивания регрессионных параметров и базовой функции риска для семейства обобщенных моделей пропорциональных интенсивностей. Для модели Ксая и SCE-модели с зависящими от базовой функции риска функциями g (x, л0 (t), p, y) данный алгоритм можно сформулировать следующим образом. Чтобы оценить вектор регрессионных параметров е _ |^pT, yt J обобщенных моделей пропорциональных интенсивностей, необходимо: 1. Положить k _ 0, задать начальное приближение Э0 _ ^p', Y' J . 2. Оценить базовую функцию риска j\q (t;Эk), для этого: а) упорядочить r различныx значений полных наблюдений (8i _ 1) по возрастанию: w1 < w2 t, v 1 £ g (x, лЛо (w/; ek), ek) =1,tj>ti v 1 3. Получить оценку 0 k+1 = argmaxln L (ek). >e, задать k = k +1 и перейти на шаг 2, иначе - считать e k+1 найденной оценкой 4. Если Jk+1 ' регрессионных параметров. В результате проведенного методами статистического моделирования исследования оценок регрессионных параметров модели пропорциональных интенсивностей, модели Ксая и SCE-модели, полученных с использованием описанного алгоритма, показано, что с увеличением объема выборок уменьшается смещение и дисперсия получаемых оценок, тогда как при повышении степени цензурирования смещение и дисперсия оценок увеличиваются. 3. Критерии проверки гипотез о параметрах В общем виде гипотеза о параметрах модели может быть записана как Но: e = 0о и проверена с помощью критерия отношения правдоподобия или критерия Вальда. Статистика критерия отношения правдоподобия имеет вид LR = 2 (lnL (e)- lnL (e° )) (5) и асимптотически распределена по закону %1 с 5 степенями свободы, где 5 - количество оцениваемых параметров модели [13]. Статистика критерия Вальда для проверки гипотезы Но: e = 0о может быть записана следующим образом: w=(e-e°)T i (e)(e-e°), (6) - оценка информационной матрицы Фишера, i, j = 1,...,s . Статистика (6) также (e) 52 ln L где i (e )=- 5ei 5e, 1 J асимптотически распределена по закону %2 с количеством степеней свободы, равным количеству оцениваемых параметров. Кроме этого, критерий Вальда позволяет проверять гипотезу о каждом параметре Но : ei = 0о, i = 1,...,sV~, в этом случае используется статистика (e i -ei )2 [ i-1(e) (7) Wi = где 1 (e)]■■ _ диагональный элемент матрицы, обратной для матрицы I (e). Статистика (7) асимптотически распределена по закону %2 с одной степенью свободы [13]. Для получения оценки информационной матрицы Фишера вычислим производные по параметрам логарифма функции правдоподобия (3): , g (ti) = g (x, Л о (ti), e), h (ti) = £ g (xj, Л о (ti), e), (e) h (ti )Ш - g (ti )5h (ti) 5 ln L h (ti ) = Zsi g (ti )h2 (ti) 5ek В случае модели пропорциональных интенсивностей (1) выражения для вычисления первой и второй производных вспомогательной функции g (x, p)_ exp(pTx) имеют вид dg(x P) (pTx) и d 2 g (x, P) dPi 5Pk • _ xk exp dPk _ xixk exp()тx), соответственно, k,i _ 1,...,m . Далее обозначим для краткости через л0 функцию л0 (t). Тогда для модели Ксая (2) производные функции g (x,лQ,p, y) _ exp()) + y)T x)(лQ)exp{Y по параметру p : _ xk exp((P + y)t x)(лq )exp(y'x)"1, _ Xl4 exp((P + y)t x)(лq)exp(y'x)-1, dg ( x, л q, p, y) 3Pk d2 g (x, л q, p, y) dPi dPk по параметру у: _ xk exp ((P + y)t x )(л q )exp(Y'x-_ (1 + exp ( x )ln ло), dg ( x, л q, p, y) dYk _ xx exp ((p +y)t x)(ло )exp(Y x)_1 (1 + 3exp (T x)ln лq + exp2 (T x)ln2 л0), dYidYk смешанная производная по параметрам: d2 g (x, л q, p, y) d2 g (x, л о, p, y) _ 1,...,m . _ xixk exp((P + y)t x)^)exp(y x) 1 (1 + exp(tx)lnло), k,l dYi dPk Для SCE-модели (3) производные функции g(x,л0,р,y)_exp()тx)(i + exp((P + Y- x-Q параметрам p и у приведены ниже: exp(-YT x)-1 по (T x) (1 + exp (T x) лQ) (1 + exp (() + y)T x) лq ) exp(-YT x)-2 dg (x, л q, p, y) dPk _ xk exp d g (x, л^ p, y) _ xx exp (p^x)( + exp )(p + Y_x)ло) Y x) 3 (1 + exp(Tx)лq)( + exp ((p + Y_ x))dPi dPk (pr x) лQ (2 + exp (pT x) лQ - 3exp (yt x) - 2 exp ((p + y)t x) лQ) exp (pT x)( - exp (yt x)- (exp (-yt x) + exp (pT x) лQ) ln (1 + exp )(P + Y_ x) лQ) +exp M)^ _ „ exp(prx)+exp()_ x)лq)^ d2g(x,ло,p,y) ../„t \/, . . Jtn . чт \exp(-ytx)-3 • _ xixk exp(ртx) ( + exp((p + y)T x)лq ) (exp(yTx) -1) {(1 - exp(ytx))exp(ртx)л, + (exp (-y' x) + exp (рт x) л0) ln (1 + exp )(P + Y_ x) л0 )| + + (1 + exp ((P + y)' x) ло) {ln ( + exp ((p + y)T x) лq ) - exp ((p + y)T x) ло -1J dYidYk Aq d g (лр0^Р'^ _ xx exp (рт x))1 + exp (рТ x) -(1 + exp ({p + Y- x)Aq ' ^ 3 [exp (рт x) k,i _ 1,...,m. -2 exp)(P + y_ x) A0 + (exp ({x) + exp (ртx_A0) ln (1 + exp )(P +y)T x_A( Как правило, при построении моделей зависимости функции выживаемости от ковариат возникает необходимость проверки гипотезы о равенстве нулю каждого параметра модели. Если гипотеза отвергается, то соответствующая данному параметру ковариата считается значимой. Важно отметить, что при построении обобщенных моделей значимость ковариаты может быть установлена на основании результатов проверки гипотезы о равенстве нулю как регрессионного параметра р, так и обобщающего параметра у, соответствующих данной ковариате. В результате исследования методами статистического моделирования распределений G (S | Но) статистик критериев отношения правдоподобия и критерия Вальда проверки гипотезы о параметрах модели пропорциональных интенсивностей и модели Ксая показано, что с увеличением объема выборок расстоя- 2 ние между эмпирическими распределениями статистик и соответствующим предельным % -распределением сокращается, распределения статистик не зависят от значений ковариат и регрессионных параметров. В случае проверки гипотезы о незначимости параметров р и у модели с пересечением функций выживаемости распределения статистики критерия Вальда оказываются далекими от предельного %2 -распределения даже при больших объемах выборок. Не менее важной задачей, решаемой с использованием критериев проверки гипотез о параметрах обобщенных моделей, является проверка гипотезы о согласии с моделью пропорциональных интенсив-ностей Кокса. Для этого формулируется гипотеза о равенстве нулю только обобщающего параметра у . Для исследования распределений статистик критериев проверки гипотезы Но: у = о о равенстве нулю обобщающего параметра SCE-модели согласно модели пропорциональных интенсивностей Кокса с регрессионным параметром р = о,5 и базовым экспоненциальным распределением моделировались цензурированные выборки объемом n = 1оо со скалярной ковариатой x = {о,1} , объем моделирования N = 1оооо . На рис. 1 и 2 представлены распределения статистик критерия отношения правдоподобия и критерия Вальда для разных степеней цензурирования, а также соответствующее предельное %2 -распределение. LR G{LR\H0) 1.00 ■■ 0.00 4.80 0.00 1.20 2.40 3.60 6.00 Рис 1. Эмпирические функции распределения статистики критерия отношения правдоподобия при проверке гипотезы о незначимости обобщающего параметра модели с пересечением функций выживаемости Рис 2. Эмпирические функции распределения статистики критерия Вальда при проверке гипотезы о незначимости обобщающего параметра модели с пересечением функций выживаемости На рис. 1 видно, что в случае полных выборок (без цензурирования) распределение статистики критерия отношения правдоподобия G (S | H0) близко к предельному %2 -распределению уже при n = 100. Однако с ростом степени цензурирования расстояние от эмпирических распределений статистики данного критерия до предельного закона увеличивается. В свою очередь распределения статистики критерия Вальда значительно отклоняются от предельного %2 -распределения даже в случае проверки гипотезы по полным выборкам. Поэтому при проверке гипотезы о параметрах SCE-модели по критерию Вальда использование предельного распределения для вычисления достигнутого уровня значимости может привести к неверному выводу. Данный факт является существенным недостатком критерия Вальда в сравнении с критерием отношения правдоподобия. 4. Критерий проверки предположения о пропорциональности рисков Существует ряд методов проверки предположения о пропорциональности рисков при построении модели пропорциональных интенсивностей: графические методы, критерии, основанные на остатках, критерии Вальда и отношения правдоподобия для проверки незначимости добавленных в модель зависимых от времени ковариат и критерии против определенных конкурирующих гипотез. Одним из графических методов является сравнение графиков оценки Каплана-Мейера распределения остатков Кок-са-Снелла [14] с функцией стандартного экспоненциального распределения. Другой графический метод основан на сравнении наблюдаемых и ожидаемых кривых выживаемости: если они достаточно близки, то предположение пропорциональности рисков выполняется. Однако все графические методы основаны на субъективной оценке и могут быть использованы лишь для предварительной оценки адекватности модели. Применение критерия проверки гипотезы о согласии при проверке предположения пропорциональности рисков является более объективным методом, чем графические методы, описанные ранее, так как позволяет получить значение статистик и достигнутые уровни значимости p. В существующей литературе рассматривается несколько критериев проверки предположения пропорциональности рисков. Модифицированный критерий, основанный на остатках Шонфельда, предложен в [7]. Для цензуриро-ванных данных в [5, 6] разработаны статистики проверки гипотезы о согласии с моделью пропорциональных интенсивностей. Однако применение данных критериев согласия для проверки предположения о пропорциональности рисков связано с необходимостью моделирования неизвестных распределений статистик [11] и / или с идентификацией базового распределения [12], что требует значительных вычислительных и временных ресурсов. В настоящей работе рассматривается критерий, предложенный М.С. Никулиным в [9], для проверки гипотезы Н о: Л х (t ;р )= exp (р т х )Ло (t) о модели пропорциональных интенсивностей против конкурирующей гипотезы о модели с пересечением функций выживаемости Hi : Л х (t ;р, у) = (l + exp ((р + у)т х) Л о (t ^^ Х) -1. Статистика критерия проверки гипотезы Но против гипотезы H1 может быть записана в векторной форме следующим образом: Т " (8) T = n-1U т D-1U, k = 1,...,m , m - размерность вектора парамет (1+ exP (р Т ^))-Ц где U = [U1,...,U„]Т , Uk =£5i - 4ln =1 ров модели пропорциональных интенсивностей; D = Е** - Е*Ео1Ет - ковариационная матрица вектора U, 'Щ - EET Sо (ti) Е* = - Z8i ni=1 M) - eee т S (ti) Ео = - £8i ni=1 Е** = -£ §i n i=1 SkM - eeE т Sо (ti) S1 (ti) Ek = Sо (ti)' k Sо (ti)' Sо(ti)= £ exp(ртxJ), S1 (ti)= £ xJ exp(ртxJ), S2 )= £ -xJ (xJ)т exp(|TX), J=1, tj Щ J=1, tj Щ J=1, tJ >ti т S1 (ti )= £ - xJ exp (рт xJ )ln (1 + exp (рт xJ )), S2 (tt )= £ - xJ (xJ) exp (рт xJ )ln (1 + exp (рт xJ ), j=1,tj щ v ' j=1, tj щ v ' S2 (tt )= £ -xJ (xJ )т exp ((3т xJ )ln2 (1 + exp (рт xJ )Ло). j=1, tj щ V У При справедливости проверяемой гипотезы статистика (8) подчиняется в пределе при n %2 -распределению с числом степеней свободы m . На рис. 3 представлены распределения статистики (8) при объемах выборок n = 5о, 1оо, 2оо, зоо в случае проверки гипотезы Но: Лх (t;р) = exp(ртх)ло (t), компоненты вектора ковариат х1 ={о,1} и х2 = {о, 1,2,3}, значения параметров р1 = о,3 и р2 = о, 6. Как видно из рис. 3, с увеличением объема выборок расстояние между эмпирическим и предельным распределением статистики (8) сокращается. Также было показано, что размерность вектора кова-риат и, следовательно, количество оцениваемых параметров влияют не только на число степеней свободы соответствующего предельного распределения, но и на близость эмпирического распределения статистики к теоретическому - чем меньше ковариат, тем ближе смоделированное распределение статистики к соответствующему предельному при том же объеме выборки. Сравним мощность критерия проверки предположения о пропорциональности рисков, предложенного М.С. Никулиным, с мощностью критериев проверки гипотез о параметрах. Действительно, гипотеза о незначимости обобщающих параметров Но: у = о модели с пересечением функций выживаемости является гипотезой о выполнении предположения пропорциональности рисков, так как в случае S1 (ti) Ek = у = 0 модель (3) является моделью пропорциональных интенсивностей. Рассмотрим следующие конкурирующие гипотезы: 1) H1 : Лх (;р, у) = (l + exp ((р + у)т х) Ло () 2) H2 : Лх (;р,у) = (1 + exp ((р + у)т х) Ло ())' ,G(T\H0) -1, у = 0,5; exp-у х -1, у = -0,5. 1.00-- 0.80 -- 0.60 -- 0.40 -- 0.20 -- 2.00 4.00 6.00 8.00 10.00 Рис 3. Эмпирические функции распределения статистики T критерия проверки предположения о пропорциональности рисков 0.00 0.00 Функции выживаемости, соответствующие гипотезам H0, H и H12 при р = 0,5 для разных значений ковариаты х = {0,1}, представлены на рис. 4-6 соответственно. Рис. 4. Функции выживаемости согласно H 0 500 Рис 5. Функции выживаемости согласно H 1.0 0.8 0.4 0.2 100 200 300 400 Рис 6. Функции выживаемости согласно H1 0.Q 0.0 На рис. 5 видно, что при разных значениях ковариаты функции выживаемости, соответствующие модели с обобщающим параметром у = 0,5, пересекаются. Расстояние между функциями выживаемости на рис. 6 при справедливости гипотезы H12 больше аналогичного расстояния на рис. 4 при справедливости нулевой гипотезы о модели пропорциональных интенсивностей. В таблице представлены полученные оценки мощности рассматриваемых критериев для объемов выборок n = 100, 200 при заданном уровне значимости а = 0,1, объем моделирования N = 10000 . Оценки мощности критериев Критерий H H2 n = 100 n = 200 n = 100 n = 200 LR 0,40 0,57 0,14 0,34 W 0,01 0,20 0,27 0,45 T 0,36 0,57 0,21 0,38 В случае конкурирующей гипотезы H, соответствующей SCE-модели с пересекающимися функциями выживаемости, наиболее мощным из рассмотренных критериев оказался критерий отношения правдоподобия. При объеме моделируемых выборок n = 100 мощность критерия Вальда оказалась меньше заданного уровня значимости, при увеличенном объеме выборок n = 200 критерий Вальда показал наименьшую мощность в сравнении с критерием отношения правдоподобия и критерием Никулина. В случае конкурирующей гипотезы H12 с расходящимися функциями выживаемости критерии со статистиками LR и T уступают по мощности критерию со статистикой W на всех рассмотренных объемах выборок, при этом мощность критерия Никулина больше мощности критерия отношения правдоподобия. Однако, учитывая смещенность критерия Вальда на паре конкурирующих гипотез H0 и н/, а также существенное отличие эмпирических распределений статистик G (S | H0) от соответствующего предель-2 ного х -распределения, для проверки предположения пропорциональности рисков рекомендуется применять критерий отношения правдоподобия со статистикой (5) и критерий Никулина со статистикой (8). Заключение В настоящей работе сформулирован итеративный алгоритм оценивания регрессионных параметров и базовой функции риска полупараметрической модели Ксая и SCE-модели. В результате проведенного методами компьютерного моделирования исследования статистических свойств оценок параметров моделей показано, что с увеличением объема выборок уменьшаются смещение и дисперсия получаемых оценок, тогда как при повышении степени цензурирования смещение и дисперсия оценок увеличиваются. На основании результатов исследования распределений статистик и мощности критериев отношения правдоподобия, Вальда и Никулина целесообразно рекомендовать использование критерия отношения правдоподобия и критерия Никулина для проверки гипотезы о выполнении предположения пропорциональности интенсивностей при построении полупараметрической модели пропорциональных интенсивностей Кокса. Применение критерия Вальда для проверки гипотезы как о значимости влияния ковариат, так и о согласии с моделью Кокса, сопряжено с возможными ошибками в определении достигнутого уровня значимости ввиду существенных отличий реальных распределений статистики при конечных объемах выборок от соответствующих предельных законов.

Скачать электронную версию публикации