Загрузите выборочные данные.

load('lightbulb.mat');

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

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

b = coxphfit(lightbulb(:,2),lightbulb(:,1), ...
'Censoring',lightbulb(:,3))

b = 4.7262

Оценка коэффициента опасности $$e^b$$ = 112.8646. Это означает, что опасность для лампочек накаливания в 112,86 раза превышает опасность для лампочек флуоресценции.

Загрузите выборочные данные.

load('lightbulb.mat');

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

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

b = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3))

b = 4.7262

Отобразите параметры управления по умолчанию для алгоритма coxphfit используется для оценки коэффициентов.

statset('coxphfit')

ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 200
          MaxIter: 100
           TolBnd: 1.0000e-06
           TolFun: 1.0000e-08
       TolTypeFun: []
             TolX: 1.0000e-08
         TolTypeX: []
          GradObj: []
         Jacobian: []
        DerivStep: []
      FunValCheck: []
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

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

coxphopt = statset('coxphfit');
coxphopt.Display = 'final';
coxphopt.MaxIter = 50;

Выполняйте coxphfit с новыми параметрами алгоритма.

b = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3),'Options',coxphopt)

Successful convergence: Norm of gradient less than OPTIONS.TolFun

b = 4.7262

coxphfit отображает отчет об итоговой итерации. Изменение максимального количества итераций не повлияло на оценку коэффициента.

Сгенерируйте данные Вейбула в зависимости от X предиктора.

rng('default') % for reproducibility
X = 4*rand(100,1);
A = 50*exp(-0.5*X); 
B = 2;
y = wblrnd(A,B);

Значения отклика генерируются из распределения Вейбула с параметром формы в зависимости от переменной предиктора X и параметр шкалы 2.

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

[b,logL,H,stats] = coxphfit(X,y);
[b logL]

ans = 1×2

    0.9409 -331.1479

Оценка коэффициента составляет 0,9409, и журнал значение правдоподобия составляет -331,1479.

Запросите статистику модели.

stats

stats = struct with fields:
                    covb: 0.0158
                    beta: 0.9409
                      se: 0.1256
                       z: 7.4889
                       p: 6.9462e-14
                   csres: [100x1 double]
                  devres: [100x1 double]
                 martres: [100x1 double]
                  schres: [100x1 double]
                 sschres: [100x1 double]
                  scores: [100x1 double]
                 sscores: [100x1 double]
    LikelihoodRatioTestP: 6.6613e-16

Ковариационная матрица оценок коэффициентов, covb, содержит только одно значение, которое равно отклонению оценки коэффициента в этом примере. Оценка коэффициента, beta, то же, что и b и равен 0,9409. Стандартная ошибка оценки коэффициента, se, равен 0,1256, что является квадратным корнем отклонения 0,0158. $$z$$-статистический, z, есть beta/se = 0.9409/0.1256 = 7.4880. Значение p, p, указывает, что эффект X является значительным.

Постройте график оценки Кокса базовой функции выживания вместе с известной функцией Вейбула.

stairs(H(:,1),exp(-H(:,2)),'LineWidth',2)
xx = linspace(0,100);
line(xx,1-wblcdf(xx,50*exp(-0.5*mean(X)),B),'color','r','LineWidth',2)
xlim([0,50])
legend('Estimated Survivor Function','Weibull Survivor Function')

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