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

Создайте псевдослучайные данные в течение времени жизни, сопоставленного с этими девятью моделями. Создайте 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])

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

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

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

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

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

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

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

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_{ij}{b}_j)$$.

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

$$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

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

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

Данные совпадают с теоретической линией для вероятностей много больше 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