Этот пример показывает, как оценить и построить совокупные функции опасности и выживания для различных групп.
Загрузите выборочные данные.
load('readmissiontimes.mat')
Данные имеют время реадмиссии пациентов с информацией об их полу, возрасте, весе, статусе курения и цензуре. Это моделируемые данные.
Создайте матрицу времени реадмиссии и цензуры для каждого пола.
female = [ReadmissionTime(Sex==1),Censored(Sex==1)]; male = [ReadmissionTime(Sex==0),Censored(Sex==0)];
Постройте график оценки Каплана-Мейера совокупной функции распределения для пациентов женского и мужского пола.
figure() ecdf(gca,female(:,1),'Censoring',female(:,2)); hold on [f,x] = ecdf(male(:,1),'Censoring',male(:,2)); stairs(x,f,'--r') hold off legend('female','male','Location','SouthEast')
Сравните функции выживания для пациентов женского и мужского пола.
figure() ax1 = gca; ecdf(ax1,female(:,1),'Censoring',female(:,2),'function','survivor'); hold on [f,x] = ecdf(male(:,1),'Censoring',male(:,2),'function','survivor'); stairs(x,f,'--r') legend('female','male')
Этот рисунок показывает, что время реадмиссии для пациентов мужского пола меньше, чем для пациентов женского пола.
Подгонка распределений Weibull ко времени реадмиссии пациентов женского и мужского пола.
pd = fitdist(female(:,1),'wbl','Censoring',female(:,2))
pd =
WeibullDistribution
Weibull distribution
A = 12.5593 [10.749, 14.6745]
B = 1.99834 [1.56489, 2.55185]
pd2 = fitdist(male(:,1),'wbl','Censoring',male(:,2))
pd2 =
WeibullDistribution
Weibull distribution
A = 4.63991 [3.91039, 5.50551]
B = 1.94422 [1.48496, 2.54552]
pd2 = fitdist(male(:,1),'wbl','Censoring',male(:,2))
pd2 =
WeibullDistribution
Weibull distribution
A = 4.63991 [3.91039, 5.50551]
B = 1.94422 [1.48496, 2.54552]
Постройте график функций выживания Вейбула для пациентов женского и мужского пола при предполагаемых функциях выживания.
plot(0:1:25,1-cdf('wbl',0:1:25,12.5593,1.99834),'-.') plot(0:1:25,1-cdf('wbl',0:1:25,4.63991,1.94422),':r') hold off legend('Festimated','Mestimated','FWeibull','MWeibull')
Распределение Weibull обеспечивает хорошую подгонку для данных.
Рассчитать совокупную функцию опасности для полов и соответствовать совокупной функции опасности Вейбула.
figure() [f,x] = ecdf(female(:,1),'Censoring',female(:,2),... 'function','cumhazard'); plot(x,f) hold on plot(x,cumsum(pdf(pd,x)./(1-cdf(pd,x))),'-.') [f,x] = ecdf(male(:,1),'Censoring',male(:,2),... 'function','cumhazard'); plot(x,f,'--r') plot(x,cumsum(pdf(pd2,x)./(1-cdf(pd2,x))),':r') legend('Festimated','FWeibull','Mestimated','MWeibull',... 'Location','North')
