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

load('lightbulb.mat');

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

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

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 не указывает, что элемент (лампа) подвергается цензуре.

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

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 отображает отчет относительно итоговой итерации. Изменение максимального количества итераций не влияло на содействующую оценку.

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

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

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

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

[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]

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

Постройте оценку Cox базовой функции оставшегося в живых вместе с известной функцией Weibull.

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')

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