Объект модели пропорциональных опасностей Кокса

Открыть Live Script

В этом примере показано, как подогнать и проанализировать модель пропорциональных рисков Кокса с помощью CoxModel объект. Модели пропорциональных рисков Кокса относятся к данным о времени жизни или времени отказа. Для справки смотрите Что такое анализ выживания? и модель пропорциональных рисков Кокса.

Сгенерируйте данные для трех моделей в течение жизни со следующими типами коэффициентов опасности. Эти модели соответствуют трем уровням расслоения; см. Расширение модели пропорциональных опасностей Кокса.

Ванна, чей уровень отказа высок в начале, уменьшается до низкого уровня, затем поднимается к постоянному уровню
Логарифмическое увеличение: $\log (x) / 10$
Постоянный $1 / 4$

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

1
1/20
1/100

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

Постройте график трех скоростей опасности, когда значение предиктора 1.

t = linspace(1,40);
plot(t,hazard(t),t,log(t)/10,t,1/4*ones(size(t)))
legend('Hazard 1','Hazard 2','Hazard 3','Location','northeast')
ylim([0 1])

Figure contains an axes. The axes contains 3 objects of type line. These objects represent Hazard 1, Hazard 2, Hazard 3.

Создайте данные для подбора кривой

Создайте псевдослучайные данные для времени жизни, связанного с девятью моделями. Создайте 9000 выборок, выбранных случайным образом (около 1000 каждого типа), инвертировав совокупные распределения. (Для получения дополнительной информации см. Раздел «Выборка обратного преобразования»).

N = 9000;
rng default
mults = [1;1/20;1/100]; % Three predictors
strata = randi(3,N,1); % Three strata
t1 = zeros(N,1);
a0 = randi(3,N,1); % Predictor
a1 = mults(a0);
v1 = rand(N,1);
for i = 1:N
    switch strata(i)
        case 1 % Bathtub
            t1(i) = zeropt(a1(i),v1(i));
        case 2 % Logarithmic
            t1(i) = zeroptold(a1(i),v1(i));
        case 3 % Constant
            t1(i) = 1 + exprnd(4/a1(i));
    end
end

Поместите данные в таблицу.

a3 = categorical(a1);
tbldata = table(t1,a3,strata,'VariableNames',["Lifetime" "Predictors" "Strata"]);

Подгонка модели Кокса

Подбор модели пропорциональных рисков Кокса к данным таблицы.

coxMdl = fitcox(tbldata,'Lifetime ~ Predictors',"Stratification",'Strata');

Построение графика выживания

Постройте график вероятности выживания как функции времени для значения предиктора 1 и три уровня расслоения.

oo = categorical(1);
plotSurvival(coxMdl,[oo;oo;oo],[1;2;3])
xlim([1,30])

Figure contains an axes. The axes with title Estimated Survival Probability contains 3 objects of type stair. These objects represent X=1, Stratification=1, X=1, Stratification=2, X=1, Stratification=3.

Постройте график вероятности выживания как функции времени для значения предиктора 1/20 и три уровня расслоения.

tt = categorical(1/20);
figure
plotSurvival(coxMdl,[tt;tt;tt],[1;2;3])
xlim([1,200])

Figure contains an axes. The axes with title Estimated Survival Probability contains 3 objects of type stair. These objects represent X=0.05, Stratification=1, X=0.05, Stratification=2, X=0.05, Stratification=3.

Постройте график вероятности выживания как функции времени для значения предиктора 1/100 и три уровня расслоения.

uu = categorical(1/100);
figure
plotSurvival(coxMdl,[uu;uu;uu],[1;2;3])
xlim([1,1000])

Figure contains an axes. The axes with title Estimated Survival Probability contains 3 objects of type stair. These objects represent X=0.01, Stratification=1, X=0.01, Stratification=2, X=0.01, Stratification=3.

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

Анализируйте подгонку

Исследуйте коэффициенты подобранной модели.

disp(coxMdl.Coefficients)

                        Beta        SE       zStat     pValue
                       ______    ________    ______    ______

    Predictors_0.05    1.5301    0.031783    48.143      0   
    Predictors_1       4.5593    0.052149    87.427      0

Заметьте, что pValue значения 0, что означает, что указанный предиктор Beta значения не равны нулю.

Просмотрите доверительные интервалы для коэффициентов модели на уровне значимости 0,01.

coefci(coxMdl,0.01)

ans = 2×2

    1.4483    1.6120
    4.4249    4.6936

Чтобы вывести опасность для предиктора уровня 0.01, вспомним определение модели Кокса:

$h (X_{i}, t) = h_{0} (t) \exp (\sum x_{i j} b_{j})$ .

The fitcox функция использует фиктивные переменные со ссылочной группой для обработки категориальных данных. В этом случае функция обрабатывает предиктор уровня 0.01 как ссылочную группу и кодирует предиктор как все 0s при подборе модели. Если вы вводите все 0 в функцию опасности, вы получаете

$h ([0, 0], t) = h_{0} (t) \exp (0 * b 1 + 0 * b 2) = h_{0} (t) \exp (0) = h_{0} (t) .$

$h_{0} (t)$ - базовая опасность, которая хранится в coxMdl.Hazard. Поэтому, чтобы получить опасность для предиктора уровня 0.01, можно изучить coxMdl.Hazard. Постройте график кумулятивной опасности базового уровня для трех уровней расслоения.

figure
hold on
for i = 1:3
    pred3 = find(coxMdl.Hazard(:,3) == i); % Find indices of stratification i
    plot(coxMdl.Hazard(pred3,1),coxMdl.Hazard(pred3,2))
end
xlabel('Time')
ylabel('Cumulative Hazard')
xlim([1 300])
legend('X = 0.01, Stratification 1',...
    'X = 0.01, Stratification 2',...
    'X = 0.01, Stratification 3','Location','northwest')
hold off

Figure contains an axes. The axes contains 3 objects of type line. These objects represent X = 0.01, Stratification 1, X = 0.01, Stratification 2, X = 0.01, Stratification 3.

Совокупная опасность для других значений предиктора $\exp (B e t a)$ умножение на базовую совокупную опасность, где $B e t a$ - выходной коэффициент.

disp(exp(coxMdl.Coefficients.Beta))

    4.6188
   95.5127

Эти относительные значения опасности близки к теоретическим значениям для данных, которые были сгенерированы умножителями 1, 1/20 и 1/100. Базовое значение соответствует умножителю 1/100, поэтому теоретические умножители равны 5 и 100.

Просмотрите linhyptest таблица.

linhyptest(coxMdl)

ans=2×2 table
         Predictor         pValue
    ___________________    ______

    {'Empty Model'    }      0   
    {'Predictors_0.05'}      0

Снова, модель требует 1/20 предиктор значения и предиктор 1 значения.

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

supports = coxMdl.ProportionalHazardsPValueGlobal

supports = 0.9730

Нулевая гипотеза для этого теста состоит в том, что данные взяты из модели пропорциональных рисков Кокса. Чтобы отвергнуть эту гипотезу, supports должен быть меньше 0,05 или какой-либо другой малый уровень значимости. Статистика указывает на поддержку гипотезы о том, что данные согласуются с моделью Кокса.

Изучение коэффициентов опасности

Вычислите коэффициент опасности для значений предиктора 1, 1/20, и 1/100 по сравнению с базовым уровнем 0 для трех уровней расслоения.

hazards = hazardratio(coxMdl,...
    categorical([1;1;1;1/20;1/20;1/20;1/100;1/100;1/100]),...
    [1;2;3;1;2;3;1;2;3],'Baseline',0)

hazards = 9×1

   95.5127
   95.5127
   95.5127
    4.6188
    4.6188
    4.6188
    1.0000
    1.0000
    1.0000

Коэффициенты опасности близки к их теоретическим значениям 100, 5 и 1, как объяснено в предыдущем разделе. Коэффициенты опасности не зависят от уровня расслоения.

Насколько хорошо постоянный уровень расслоения опасности соответствует теории?

Уровень расслоения 3 имеет постоянную скорость опасности 1/4. Теоретически, постоянная скорость опасности означает экспоненциальную функцию выживания, логарифм которой является прямой линией. Постройте график выживания расслоения уровня 3 с помощью значения предиктора 1/20, что приводит к экспоненциальной скорости 1/80.

tt = categorical(1/20);
h = figure;
axes1 = axes('Parent',h);
plotSurvival(coxMdl,axes1,tt,3);
xlim([1 400]);
axes1.YScale = 'log';
hold on
tms = linspace(1,400);
semilogy(tms,exp(-tms/80),'r--')
hold off

Figure contains an axes. The axes with title Estimated Survival Probability contains 2 objects of type stair, line. This object represents X=0.05, Stratification=3.

Данные совпадают с теоретической линией для вероятностей значительно выше 1/100.

Вспомогательные функции

Эта функция создает коэффициент опасности для ванны.

function h = hazard(x)
h = exp(-2*(x - 1)) + (1 + tanh(x/10 - 3));
end

Эта функция создает интеграл уровня опасности ванны от 1 на x.

function eh = exphazard(x)
eh = 1/2 - exp(-2*(x-1))/2;
eh2 = (10*log(cosh(x/10 - 3)) - 10*log(cosh(1/10 - 3)) + x - 1);
eh = eh + eh2;
end

Эта функция решает для корня совокупной скорости опасности с умножителем a на уровень v.

function zz = zeropt(a,v)
zz = fzero(@(x)(1 - exp(-a*exphazard(x)) - v),[1 100*max(1,1/a)]);
end

Эта функция создает интеграл логарифмической скорости опасности с умножителем 1/10 от 1 на x.

function cr = cumrisk(x)
cr = 1/10*(x.*(log(x) - 1) + 1);
end

Эта функция решает для корня совокупной скорости опасности с умножителем a на уровень v.

function zz = zeroptold(a,v)
zz = fzero(@(x)(1 - exp(-a*cumrisk(x)) - v),[1 50*max(1,1/a)]);
end

См. также

CoxModel | fitcox

Документация