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

Скрипт Open Live Script

В этом примере показано, как соответствовать и анализировать Cox пропорциональная модель опасностей использование CoxModel объект. Cox пропорциональные модели опасностей относится ко времени жизни или данным времени отказа. Для фона смотрите то, Что Анализ Выживания? и Cox Пропорциональная Модель Опасностей.

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

Ванна, интенсивность отказов которой высока вначале, уменьшается до низкого уровня, затем поднимается к постоянному уровню
Логарифмически увеличение: $\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 object. The axes object contains 3 objects of type line. These objects represent Hazard 1, Hazard 2, Hazard 3.

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

Создайте псевдослучайные данные в течение времени жизни, сопоставленного с этими девятью моделями. Создайте 9 000 выборок, выбранных случайным образом (приблизительно 1 000 из каждого типа) путем инвертирования кумулятивных распределений. (Для получения дополнительной информации смотрите, что Инверсия преобразовывает выборку).

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"]);

Подбирайте модель Cox

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

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 object. The axes object 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 object. The axes object 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 object. The axes object 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 предикторов уровня, вспомните определение модели Cox:

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

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

$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 object. The axes object 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 предиктор значения.

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

supports = coxMdl.ProportionalHazardsPValueGlobal

supports = 0.9730

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

Исследуйте отношения опасности

Вычислите отношение опасности для значений предиктора 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 object. The axes object 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

Смотрите также

fitcox | CoxModel

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