Загрузите образцы данных.

load('lightbulb.mat');

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

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

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

b = 4.7262

Оценка коэффициента опасности eb $${}^{}$$= 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')

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