Сравнение оценок параметров специальной модели кривой дожития для выборки с интервальным цензурированием Текст научной статьи по специальности «Математика»

Сер. 10. 2009. Вып. 2

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 519.233.2

А. И. Коробейников

СРАВНЕНИЕ ОЦЕНОК ПАРАМЕТРОВ СПЕЦИАЛЬНОЙ МОДЕЛИ КРИВОЙ ДОЖИТИЯ ДЛЯ ВЫБОРКИ С ИНТЕРВАЛЬНЫМ ЦЕНЗУРИРОВАНИЕМ

1. Введение. Пусть X - случайная величина, имеющая смысл времени до наступления какого-либо события, обычно называемого отказом. Функцией дожития называется зависимость вероятности отсутствия отказа к некоторому моменту времени t:

S (t) = P (X >t).

Мы будем рассматривать специальную модель кривой дожития в виде произведения экспоненты на косинус:

S (ж; г), т) = 1 — F (ж; г/, т) = exp (—r/x) cos > V > 0, 0 < х < г, (1)

здесь F (x) - функция распределения случайной величины X.

Эта модель была впервые предложена в работе [1]. Формальным основанием для ее использования в биометрии служит разработанный в [2] лагранжево-гамильтоновый механизм для биологической систематики, согласно которому компоненты в виде произведения экспоненты на косинус являются решением системы, описывающей взаимодействие факторов органа и организма. Модель применялась для описания динамики хронического гломерулонефрита [1], раневых процессов [3], хронического генерализованного пародонтита [4].

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

Как правило, наличие цензурирования ведет к достаточно большой потере информации о выборке. Особенно сильно это проявляется при работе с параметрическими моделями: заметная смещенность оценок параметров даже при оценивании по выборкам достаточно большого размера, влияние механизма цензурирования на статистические свойства оценок, отсутствие асимптотической нормальности и т. п. В таких ситуациях значительный интерес представляют робастные оценки параметров, устойчивые к существенной потере информации о выборке, возможно, за счет некоторой потери эффективности. Например, в работе [5] исследуются оценки, минимизирующие расстояние Хеллинджера между теоретическим и эмпирическим распределениями,

Коробейников Антон Иванович — аспирант кафедры статистического моделирования математикомеханического факультета Санкт-Петербургского государственного университета. Научный руководитель: докт. физ.-мат. наук, проф. С. М. Ермаков. Количество опубликованных работ: 6. Научные направления: теория вероятностей, математическая статистика. E-mail: [email protected].

© А. И. Коробейников, 2009

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

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

2. Модель интервального цензурирования. Как было отмечено раньше, наблюдать непосредственно время отказа X, как правило, невозможно по тем или иным объективным причинам. Часто вместо этого можно наблюдать пару случайных величин (Ь, Н), такую, что

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

Так, если обозначить через X момент рецидива, У - момент последней диагностики, то величина (Ь,Н) может быть задана как

Для независимых X и У эта модель известна как интервальное цензурирование первого типа [6] и пара (Ь, Н) обычно представлена в виде информации о текущем состоянии (У, I), где I - индикатор интересуемого события.

Естественным обобщением данной модели является интервальное цензурирование к-го типа [7]. В таком случае предполагается наличие к отдельных моментов диагностики У\ < ■■■ < УЬ, не зависящих от X, и доступна информация, какому из интервалов (—^о, У1], ..., (Уь, +то) принадлежит реализация X.

На практике, однако, нечасто можно наблюдать возможность применения модели интервального цензурирования к-го вида при к ^ 2 в чистом виде: маловероятно, что число наблюдений состояния индивида будет всегда одинаковым. Отдельные индивиды могут выбывать из процесса наблюдения до окончания эксперимента, по тем или иным причинам пропускать диагностику и т. д.

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

Дадим формальное описание модели интервального цензурирования. Для этого достаточно определить зависимость пары наблюдаемых случайных величин (Ь, Н) от ненаблюдаемой X. Нам понадобится набор вспомогательных случайных величин.

Пусть К - неотрицательная целочисленная случайная величина. Для к = 1, 2,... положим Ук = {Ук,}^=1 - набор случайных величин и Укд < ••• < Ук,к. Кроме того, положим формально Ук,0 = —ж, Ук,к+1 = +ж. Всюду далее будем предполагать, что набор (К, У1, У2,...) от X не зависит.

Рассмотрим событие {К = к}. Для набора Ук введем вектор индикаторов Г (Ук) = (Г1,...,Гк+1) с Г,- = I {X ^ Ук —1}. Положим J = ^ к=1 Г,. Наконец, определим (Ь,Е) = (Ук, з-1 ,Ук , з).

Фактически набор Ук задает случайное разбиение числовой оси на полуинтервалы: М = (—ж, Ук, 1] и (Ук, 1, Ук,2] и ••• и (Ук,к, +ж). Случайная величина J - это номер полуинтервала, которому принадлежит X: X € (Ук,з-1,Ук, з], а (Ь,Е) - его концы.

Итак, при фиксированном к имеем

’ (—ж, Ук,\), X < Ук,1,

(ЬК)={(Ук’1 ,Ук,2) , Ук,1 <X ^ Уk,2, (2)

, (Ук,к , +ж) , Ук,к < X.

Величина X имеет смысл случайного времени до наступления некоторого события (в частности, времени элиминации объекта, времени наработки на отказ). Случайная величина К - это количество наблюдений за состоянием индивида, а случайный набор Ук - набор временных точек, в которых производится наблюдение за состоянием индивида. Например, периодический контроль состояния индивида на протяжении всего периода проведения эксперимента может быть представлен в виде следующей модели:

з-1 ( з-1 Л

Ук,3 =^2^, К = шах < з ^ 1:^Аг <Т>, (3)

г=1 I г=1 )

где Т - общая продолжительность наблюдений, а Аг - промежутки времени между отдельными моментами контроля состояния индивида.

3. Оценки параметров модели. Пусть X1,...,Xn - п независимых одинаково распределенных случайных величин - набор времен до наступления интересуемого события. Функция распределения Г для XI задается формулой (1). Пусть (Ь1, Д1),..., (Ьп, Яп) - наблюдаемые цензурированные случайные величины, определенные в соответствии с формулой (2).

Запишем эмпирическую функцию правдоподобия для выборки с интервальным цензурированием (формально полагаем Г (—ж) = 0, Г (+ж) = 1):

п п

Ьп (Г) = П [Г(Щ) — Г(Ц)] = П Рр (Iз). (4)

з=1 з=1

Ее логарифм равен

п

1п (Г) = ^ё Рр I). (5)

з=1

Здесь символом Рр (13) обозначена вероятность события {X, € I3 } для интервала I,.

Всюду далее будем обозначать символом в вектор параметров модели: в = (ц, т), соответственно функцию распределения - Гд. Пространство параметров 0 представляет собой декартово произведение полупрямой на вещественную прямую: 0 = М+ х М.

Рассмотрим разные способы получения оценок двумерного параметра в:

1) Оценивание методом максимального правдоподобия:

П

вп = arg max ln (Fe) = arg ma^ log Pe (Ij). (6)

9ев 9ев

2) Оценивание при помощи следующей двухшаговой процедуры:

а) вычисление оценки максимального правдоподобия для функции распределения

Fe

Fn = arg max ln(G), (7)

GeG

где G - класс всех функций распределения с носителем [0, ж).

б) получение оценки вп посредством минимизации расстояния между Fn и параметрическим семейством F = {Fe : в e 0}

вп = argmind (Fn,F) = argmind (Fn,Fe) (8)

feF v ' eee v '

для некоторой «функции расстояния» d (■, ■).

Статистические свойства (например, состоятельность) оценок en неизвестны, что вынуждает искать альтернативные методы построения оценок параметров модели. Отметим, что в некоторых так называемых нерегулярных случаях оценки максимального правдоподобия могут не обладать свойством состоятельности [8].

Способ построения оценок en был выбран нами не случайно. Пусть B - борелевская сигма-алгебра на R. Введем на B меру ц:

оо к

М (B) = Y, P (K = k)^2 P (Yk,j e B | K = k), B eB. (9)

k=l j = 1

Пусть p - метрика пространства L1 (м), т. е. для F1, F2 e G:

p (Fi,F2)= IFi(t) - F2 (t)| dM(t),

J 0

тогда при E(K) < ж имеет место сходимость p(Fn, F) ^ 0 (п.н.) [9].

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

Сами по себе задачи оптимизации (6), (8) не являются чем-то нестандартным: мы имеем дело с конечномерной задачей оптимизации в евклидовом пространстве. Задача (7) существенно отличается от этих двух задач (например, непосредственно из условия не ясно, при каких условиях существует хотя бы приближенное решение задачи), остановимся на ней подробнее.

4. Построение непараметрической оценки функции распределения. Для нахождения оценки Fn непосредственно из уравнения (7) необходимо максимизировать логарифм функции правдоподобия ln(F) по функциональному пространству G, тем самым, задача оптимизации (7) является бесконечномерной.

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

4.1. Редукция.

Определение 1. Интервал H назовем внутренним подынтервалом (в англоязычной литературе иногда используется термин «максимальное пересечение» — maximal intersection) для набора интервалов {Ij}™=1, если выполняются следующие условия:

1) H представим в виде пересечения интервалов f]jEa Ij для некоторого набора индексов а С{1,...,п};

2) для всякого номера j либо H П Ij = H, либо H П Ij = 0.

Заметим, что если H, H' - два различных внутренних интервала, то H П H' = 0.

В работе [11] рассматривается теоретико-графовый подход к построению набора всех внутренних подинтервалов для данного набора интервалов {Ij}j=1. Показано, что внутренние подинтервалы соответствуют максимальным кликам графа пересечений для данного набора интервалов.

Поведение оценки максимального правдоподобия Fn описывается леммами из [11].

Лемма 1. Пусть H = {Hl,..., Hm} - набор всех внутренних подинтервалов для {Ij}™=1- Тогда любая функция распределения из G, возрастающая вне H, не может быть оценкой максимального правдоподобия для функции распределения F.

Лемма 2. При фиксированных значениях Pf (Hi) функция правдоподобия (5) не зависит от поведения функции распределения F внутри каждого интервала Hi.

Эти леммы имеют очень важный смысл, с точки зрения процедуры вычисления оценки Fn. Лемма 1 позволяет свести сложную бесконечномерную оптимизационную задачу (7) по всему пространству G к аналогичной, но уже по некоторому подпространству G. Лемма 2, с одной стороны, фактически означает, что оценкой максимального правдоподобия для F является целый класс эквивалентности функций распределения и нет никакой дополнительной информации для того, чтобы выбрать конкретного представителя из этого класса, с другой - позволяет при построении оценки Fn ограничиться поиском исключительно в классе скачкообразных на H функций и тем самым получить конечномерную оптимизационную задачу.

Действительно, имеется ровно m интервалов, где функция Fn может производить скачки, потому для ее определения достаточно найти величину этих скачков.

Для i = 1,. ..,m положим pi = Pf (Hi) и p = (pl,... ,pm). Тогда для j = 1,...,n имеем

После этого можно переписать выражение для функции правдоподобия (4) и ее логарифма (5) в терминах величин рі и а^:

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

m

m

n

n

m

j=1

j=i Li=i

n

n

m

1) Редукция. Вычисление набора Н всех внутренних подинтервалов. Это можно сделать за О (п 1с^ п) операций и О (п) памяти [12].

2) Оптимизация. Решение конечномерной оптимизационной задачи

р = а^шах 1п (рі,... ,рт) = а^шах 1п (р), (10)

V р^Т

здесь

т ^

р е Мт : рі > 0, і = 1,...,т;^2 Рі = 1 ґ .

і=1 )

Таким образом, получаем т-мерную выпуклую задачу оптимизации.

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

Теорема 1. Величина р, определенная уравнением (10), существует.

Доказательство. Полагая к^(0) = —ж, продолжим логарифм по непрерывности на расширенную числовую прямую К = К и {—ж}. Тогда логарифм функции правдоподобия 1п (р) является непрерывной функцией, заданной на непустом компакте V С Кт и принимающей значения в К. Тем самым существует максимум Іп (р) на V, и величина р, являющаяся точкой максимума Іп (р), в действительности существует. □

4.2. Оптимизация. ЕЫ-алгоритм. Для решения задачи оптимизации (10) будем применять ЕМ-алгоритм [13]. В силу лемм 1, 2, «полные данные» представляют собой знание, какому интервалу Ні принадлежит наблюдение с номером і. Это, естественно, не эквивалентно знанию точного времени отказа X^: в некотором смысле механизм цензурирования диктует выбор информативной модели для полных данных.

Обозначим [іі^ = I {Xj е Ні} и пі = ^п=і Ріі. Тогда естественной моделью для

полных данных является мультиномиальная модель на Н с вектором вероятностей (р1,... ,рт). Запишем логарифм (ненаблюдаемой) функции правдоподобия для полных данных іп (рі,..., Рт):

т

іп(Р) = іп (Рі,...,Рт) = ^2 пі 1сё Pi, р =(Рі,...,Рт) ер. (11)

і=і

Общая схема ЕЫ-алгоритма. ЕМ-алгоритм решает задачу максимимизации логарифма функции правдоподобия для неполных данных (10) неявно, применяя для этого логарифм функции правдоподобия для полных данных (11). Так как величины пі не наблюдаются, следовательно, не наблюдается и іп(р), то используется условное математическое ожидание іп (р) по отношению к наблюдаемым случайным величинам aij. Приведем (к + 1)-тую итерацию ЕМ-алгоритма:

Е-шаг. Вычислить функцию Q (р; р(к)), где

Я(р; р(к)) = Ер(к) (іп(р) \Н).

М-шаг. Выбрать в качестве следующего приближения р(к+1) точку максимума функции Q (р; р(к)):

Q(Уk+l); р(к)) > Q(p; р(к)) , Ур еР.

Всюду далее символом Ер(к) (•) будет обозначаться математическое ожидание с использованием вектора параметров р(к).

Е- и М-шаги алгоритма последовательно выполняются до тех пор, пока последовательность р(к) не достигнет стационарной точки.

Теорема 2. Применение процедуры ЕМ-алгоритма приводит к увеличению значения функции правдоподобия, т. е.

1и(р{к+1)) > ¡и(р(к)), к = 0,1, 2,... .

Доказательство. Следует из общей теории ЕМ-алгоритма [13, 14]. □ ЕМ-алгоритм для оценивания Ри. Для того чтобы применить ЕМ-алгоритм к максимизации функции ¡П (р), необходимо получить выражение для функции

Я (р; р(к)).

Заметим, что ¡С(р) линейна по щ, которые, в свою очередь, являются простой суммой [3^. Из определения [3^ непосредственно имеем

Ер(к) (^ \ Н) = Ер(к) (I {Xj е Ні }\Н)= Р X е Ні \ Xj е ^ ) =

(к)

^2 р.

і=1

(к)

Отсюда, суммируя по і = 1,...,п, получаем

П П

Ні(р(к)) = Ер(к) (щ \Н) =^2 Ер№) (^ \Н)=]Г

j=l j=l

(к)

Ет (к)

¿=1 ач Рі .

После этих приготовлений несложно вывести явное выражение для функции Q (р; р(к))

т

Q (р; р(к)) = Ер« (іп (р) \Н) = £ Ні (р(к) )1с%рі.

=і

М-шаг ЕМ-алгоритма заключается в максимизации функции Q (р; р(к)) по пространству V. Этот максимум несложно получить в явной форме. Рассмотрим функцию

І (рі,...,рп,^) = ^ні (р(к) )^рі — хЕрі.

=і

=і

Приравняем частные производные /(рі,...,рп,Х) по рі,... ,рт к нулю:

д/ Ціір{к)) л „ ч Иі(р{к))

—— =---------------Л = 0 =----------.

дрі рі X

Суммируя рі по і = 1,...,т, находим

Ні (р(кї)

X

=>

X = ^2 Ні (р(к)) = п.

=і =і =і Окончательно для пересчета р имеем явную формулу

(й+1) = М»(р(й)) = 1 * п п

j=l

рі

(к)

Ет

¿=1а ij р

(к)

(12)

Замечание 1. EM-алгоритм, вообще говоря, не всегда сходится к Fn .В работе [15] приведены условия типа Куна-Таккера, позволяющие проверить, действительно ли данный вектор (pi,... ,pm) является решением задачи (10).

EM-алгоритм был выбран из-за простоты описания и реализации, он хорошо работал в интересующих нас случаях и на практике всегда сходился к оценке максимального правдоподобия Fn.

Существуют другие методы решения оптимизационной задачи (10). Среди наиболее часто используемых отметим метод ICM (Iterative Convex Minorant) [17]. Детальный обзор методов вычисления непараметрической оценки максимального правдоподобия функции распределения можно найти в [11].

5. Экспериментальное сравнение свойств оценок. Интересующие нас свойства оценок (состоятельность и скорость сходимости к истинным значениям параметров) изучались на модельных выборках.

Для моделирования значения параметров были выбраны совпадающими с оценками, полученными в работе [18] для реальных данных из кардиологии: ц = 0.125, т = 16. Цензурирование производилось по схеме, представленной в формуле (3): промежутки времени Ai между отдельными моментами контроля состояния индивида имели экспоненциальное распределение со средним 1. Общая продолжительность наблюдений T равнялась 8, таким образом, случайная величина K — 1 имела пуассоновское распределение с параметром А = 8.

На рис. 1 изображены модельная кривая дожития и ее непараметрическая оценка (размер выборки n = 1000). Отметим, что скачок в точке x = 8 связан с отсутствием какой-либо информации о времени отказа для тех индивидов, у которых оно больше 8.

Рис. 1. Модельная кривая дожития (точечная линия) и ее оценка (сплошная)

Объем выборки n варьировался от 500 до 5000 с шагом 50 и от 5000 до 10 000 с шагом 100 (в связи с высокой трудоемкостью вычисления оценок при больших объемах выборки). Моделирование производилось 100 раз, соответственно для каждого значения n получалась выборка из m = 100 оценок.

5.1. Вычислительная процедура. Все вычисления производились в пакете R версии 2.7.0 [19] на четырехпроцессорном компьютере под управлением ОС Linux. Ввиду значительной вычислительной сложности задачи (время, необходимое для получения одной реализации оценок параметров, достигало нескольких минут для достаточно больших объемов выборки n) моделирование производилось параллельно на всех процессорах компьютера. Для организации процедуры распределения задач по отдельным процессорам использовались пакеты snow и Rmpi [20].

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

Оценки f]n, тп. Несложно видеть, что целевая функция в задаче оптимизации (6) не является вогнутой по паре переменных (п,т). Следовательно, возможно существование нескольких локальных максимумов логарифма функции правдоподобия. Однако функция (5) оказывается вогнутой по п.

Лемма 3. При фиксированном значении аргумента т логарифм функции правдоподобия (5) является вогнутой функцией аргумента п.

Доказательство. Действительно,

nn

ln (п,Т) = Elog ЬП) (п,Т) = Е П (п,т), j=i j=i

Ln] (V, т) = F(Rj) ~ F(Lj) = exP (~vLj) cos - exp (-rjRj) cos (J^Rj) ■

( j ) ( j )

Функция In = log Ln вогнута по п, так как

д2 ,(i) = ехр (-?7(Дд + Lj)) cos (^:Дд) cos (Дд- - ¿д-)2 q

дг12 " (exp (-r/Rj) cos (fpRj) - exp (-r/Lj) cos (^Lj))2

И ln (п, t) является вогнутой функцией аргумента п как сумма вогнутых. □

Для нахождения глобального максимума функции (5) применялась следующая процедура мультистарта:

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

а) перебор значения т по узлам сетки на отрезке [T, T] с некоторым шагом

h: = T + kh, 0 ^ k ^ (T — T)/h;

б) поиск точки максимума функции ди(п) = ln (п,ти). Обозначим щ = arg max ди(п).

2. Результатом шага 1 является набор точек (пи, ти), к = 0,1,... .В силу леммы 3, точка пи является (возможно, не единственной) точкой глобального максимума функции ди (п). Точки (пи ,ти) использовались в качестве начального приближения для поиска максимума логарифма функции правдоподобия (5) по двум переменным при помощи встроенного в пакет R алгоритма L-BFGS-B [21]. Результатом работы алгоритма является точка (локального) максимума (п*и,тк).

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

Оценки fjn, тп. Как было отмечено выше (лемма 2), оценка максимального правдоподобия Fn, вообще говоря, не единственная: решением оптимизационной задачи (10) является целый класс эквивалентности функций распределения. Однако всегда можно указать такую функцию F+, принадлежащую данному классу, что для любой оценки максимального правдоподобия Fn выполняется Fn ^ F+ всюду: для этого достаточно взять кусочно-постоянную функции со скачками на левых концах внутренних подин-тервалов соответственно.

Использовалась функция расстояния d (■, ■) в виде

( т 2\ 1/2

d (f, к) = ^2 F a) - т+ a)) ) , (13)

где ai - левый конец г-го внутреннего подинтервала.

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

Процедура вычисления оценок естественным образом включает в себя такие шаги:

1. Вычисление оценки максимального правдоподобия Fn при помощи EM-алгоритма. Итерации (12) продолжаются до тех пор, пока результат текущего шага перестанет отличаться от предыдущего. Таким образом, критерием прекращения

7 ( (k) (k-lA 2 .

итераций на шаге к является условие > ,i=i [pi — pi I < e для некоторого малого

значения e. Для вычисления оценок использовали e = 10-10.

После этого для полученного вектора (pi,... ,pm) проверялись условия типа Куна-Таккера из работы [15]. Предполагается, что при их невыполнении EM-алгоритм запус-

( (0) (0)° /

кается заново с другим значением вектора начальных данных I pi ,... ,pm I (отметим,

что в рассматриваемом случае этого не наблюдалось никогда).

2. Решение задачи оптимизации (8). Поскольку функция (13) не является выпуклой, то применялась процедура мультистарта с перебором начальной точки по узлам сетки на прямоугольнике [T, Т] х [0, fj\. Для нахождения точки минимума при фиксированном начальном приближении использовался встроенный в пакет R алгоритм L-BFGS-B.

Оценка сверху fj параметра п находилась следующим образом. Заметим, что для q = P (X > T) имеет место

(пТ\ 1 1

q=l-Fe(T) = exp(-r]T)cos l^-j =>т) = -- log—^ ~f 1оё«-

Величину q можно оценить при помощи F+. Отметим, что точка Т является регулярной (в смысле [9]), т. е. для меры л из формулы (9) имеет место

Л {(Т — e,T\} > 0, Уе > 0,

и F+ (Т) ——a.s. Fq (Т) [9] (вообще говоря, для выбранной модели цензурирования все точки полуинтервала (0,Т\ являются регулярными). Таким образом, 1 -F+ (Т) служит

состоятельной оценкой для q и можно положить fj = — ^ log ^1 — F^ (Т)^ .

5.2. Скорость сходимости и дисперсия оценок. Для исследования состоятельности и скорости сходимости оценок строились графики зависимости оценки стандартного отклонения оценок fjn,тп и fjn,fn от объема выборки п. Предполагали, что оценки параметров имеют типичную скорость сходимости порядка Op (1 /л/п), поэтому для построения зависимостей оценки были центрированы и нормированы на л/п (для сравнения на рис. 2, а изображена зависимость стандартного отклонения оценок параметров без дополнительной нормировки).

Из изображения зависимостей стандартного отклонения оценок от объема выборки на рис. 2 очевидно, что наличие интервального цензурирования не изменило скорость сходимости оценок максимального правдоподобия fjn. Оценки fn демонстрируют такие же скорость сходимости и сравнимую дисперсию.

Для оценок тп, тп естественно ожидать достаточно большую дисперсию: общая продолжительность наблюдений Т < т, а, значит, достаточно заметная часть выборки представляла собой полуинтервалы вида (8, +го) (при данных значениях параметров и механизме цензурирования такие наблюдения составляли порядка 36% от общего объема выборки).

2000 4000 6000 8000

Рис. 2. Зависимость стандартного отклонения оценок fjn и г)п от объема выборки

Заметим, кроме того, что носитель распределения зависит от параметра т, поэтому в случае оценивания методом максимального правдоподобия в ситуации без цензурирования оценка тп демонстрировала бы сверхсходимость (т. е. стандартное отклонение оценки имело бы порядок Ор (1 /п) против обычного Ор (1 /л/п)). Наличие же цензурирования привело к пропаданию сверхсходимости.

На рис. 3 изображена зависимость л/п std (fn) и л/п std (fn) от объема выборки.

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

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

5.3. Чувствительность оценок к выделяющимся наблюдениям. Кроме наличия цензурирования, для данных типа времени жизни характерна неоднородность:

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

1501-

100

50

2000 4000 6000 8000 10 000

п

Рис. 3. Зависимость y/n std(fn) и л/rïstd (тп) от объема выборки

. std(v77(xn-T0)) □ std(vn(xn -т0))

________________I__

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

Типичным случаем является неоднородность по моменту начала наблюдений: часть выборки получает дополнительный сдвиг. Для сравнения оценок в таком случае моделировалась смесь распределений: с параметров п = 0.125, т = 16 и п = 0.125, т = 17 с дополнительным сдвигом влево на 1. Таким образом, правый конец носителя распределений компонентов смеси совпадал. Были использованы выборки из m = 500 оценок, каждая из которых была получена по выборке из n = 1000 наблюдений, при этом доля наблюдений, взятых из второй компоненты смеси, составила 0.05.

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

Для сравнения оценок использовалась относительная эффективность [22, 23]

' D(«2)/E(«|)

По результатам моделирования были получены следующие величины относительной эффективности оценок:

RE (Пп, Пп) « 1.2, RE (Т, т„) « 1.9.

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

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

Автор благодарит рецензента за ценные замечания, существенно улучшившие статью.

Литература

1. Барт А. Г., Бондаренко Б. Б., Бойко В. И. Математический анализ течения ХГН // Гломерулонефрит / под ред. С. И. Рябова. М.: Наука, 1980. С. 213-225.

2. Калинин О. М. О единых математических трактовках в биологической систематике и динамике популяций и о связи диффузии с нелинейными уравнениями // Проблемы кибернетики. 1972. Вып. 25. С. 107-117.

3. Барт А. Г. Анализ медико-биологических систем. Метод частично-обратных функций. СПб.: Изд-во С.-Петерб. ун-та, 2003. 280 с.

4. Мадай Д. Ю., Барт А. Г., Рыжак Г. А. и др. Репаративная эффективность вилона при лечении больных пожилого и старческого возраста с хроническим генерализованным пародонтитом и сопутствующими возрастными соматическими заболеваниями. Великий Новгород: Изд-во Новгород. гос. ун-та. им. Ярослава Мудрого, 2006. 26 с.

5. Beran R. J. Minimum Hellinger distance estimates for parametric models // Ann. Stat. 1977. Vol. 5. P. 445—463.

6. Groeneboom P., Wellner J. A. Information Bounds and Nonparametric Maximum Likelihood Estimation. Basel: Birkhauser, 1992. 126 p.

7. Wellner J. A. Interval censoring case 2: alternative hypotheses // Analysis of Censored Data / H. Koul, J. V. Deshpande (eds.). Proc. of the Workshop on Analysis of Censored Data. IMS Lecture Notes, Monographs Ser. 1995. Vol. 27. P. 271—291.

8. Cheng R. C. H., Traylor L. Non-Regular maximum likelihood problems // J. R. Stat. Soc. Ser. B. 1995. N 57. P. 3-44.

9. Schick A., Yu Q. Consistency of the GMLE with mixed case interval censored data // Scand. J. Stat. 2000. Vol. 27. P. 45-55.

10. Turnbull B. W. The empirical distribution function with arbitrarily grouped, censored and truncated data // J. R. Stat. Soc. Ser. B. 1976. Vol. 38. P. 290—295.

11. Gentleman R., Vandal A. C. Computational algorithms for censored data problems using intersection graphs // J. Comput. & Graph. Stat. 2001. Vol. 10. P. 403-421.

12. Maathuis M. H. Reduction algorithm for the NPMLE for the distribution of bivariate interval-censored data // J. Comput. Graph. Stat. 2005. Vol. 14. P. 352—362.

13. Dempster A. P., Laird N. M., Rubin D. B. Maximum likelihood data from incomplete data via the EM algorithm // J. R. Stat. Soc. Ser. B. 1977. N 39. P. 1-38.

14. McLachlan G. J., Krishnan T. The EM Algorithm and Extensions. New York: Wiley, 1997. 274 p.

15. Gentleman R., Geyer C. J. Maximum likelihood for interval censored data: Consistency and computation // Biometrika. 1994. Vol. 81. P. 618-623.

16. Nettleton D. Convergence properties of the EM algorithm in constrained parameter spaces // Canad. J. of Stat. 1999. N 27. P. 639-648.

17. Jongbloed G. The iterative convex minorant algorithm for nonparametric estimation // J. Comput. Graph. Stat. 1998. Vol. 7. P. 301—321.

18. Коробейников А. И., Бондаренко Б. Б., Алексеева Н. П. Оценка критических периодов в кривых дожития // Тез. Всерос. науч.-практ. конференции «Высокотехнологичные методы диагностики и лечения заболеваний сердца и крови». Артериальная гипертензия. 2008. Т. 14, № 1. С. 26.

19. R Development Core Team. R: A language and environment for statistical computing // R Found. for Stat. Comp. Vienna, Austria, 2008. URL: http://www.R-project.org.

20. Hao Yu Rmpi: Interface (Wrapper) to MPI (Message-Passing Interface) // R package version 0.5-5. 2007. URL: http://www.stats.uwo.ca/faculty/yu/Rmpi.

21. Zhu C., Byrd R. H., Nocedal J. Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization // ACM Transactions on Mathematical Software. 1997. Vol. 23, N 4. P. 550-560.

22. Huber P. J. Robust statistical procedures // SIAM Regional Conference Ser. in Appl. Math. 1977. N 27. 56 p.

23. Хьюберт П. Дж. Робастность в статистике / пер. с англ.; под ред. И. Г. Журбенко. М.: Мир, 1984. 303 с.

Статья рекомендована к печати проф. Л. А. Петросяном.

Статья принята к печати 25 декабря 2008 г.