* PHREG4.SAS
* xref:
* input: SPRINGS.DAT
* output:
* does: - Illustrate the computational algorithm for partial likelihood
          in SAS IML.
        - time: no ties.
        - one covariate only.
        - no stratification.
        - Compare to proc phreg.
***************************************************;
filename INF 'SPRINGS.DAT';
options linesize=176;
***************************************************;
data A;
  infile INF;
  input time delta group stress;
  label
    time  = 'Cycles to failure (x1000)    ..'
    delta = 'Event: 0=censored 1=observed ..'
    group = 'Group: 1,...,6 stress levels ..'
    stress= 'Stress (N/mm^2)              ..'
  ;
  if (group >= 2) and (group <= 4);
  run;
***************************************************;
proc sort data=A nodupkey; by time; run;
* just to make sure no tied times exist;
***************************************************;
proc phreg data=A;
  model time * delta(0) = group;
  run;
***************************************************;


* sort the data properly before starting IML;
proc sort data=A; by descending time delta; run;
***************************************************;
proc iml worksize=1024;
reset nocenter   name pagesize=65 linesize=178;
***************************************************;

start failonly (t, x, y, d);

  * subset (t, x, y) only d = 1;
  i = loc(d);
  t = t[i,];
  x = x[i,];
  y = y[i,];

finish;

***************************************************;

start phreg (beta, score , info, time, delta, x);

  * Proportional Hazards Model Fitting;
  * time, delta, x:  n by 1;
  * sorts time, score, info by descending time;
  * time, score, info: a subset (delta=1) is returned;
  * delta, x: permuted;

  run sort3d (time, delta, x);
  beta = 0;
  oldbeta = 0;
  niter = 0;

  do until (converge | (niter=10)); * main loop;
    psi=exp(x#beta);
    s0 = cusum(psi);
    s1 = cusum(psi#x);
    s2 = cusum(psi#x#x);
    score = (delta#(x - s1/s0));
    info  = (delta#(s2/s0 - (s1/s0)##2));
    beta = beta + sum(score)  / sum(info);
    converge = (abs(beta-oldbeta) < 1.0e-4);
    oldbeta = beta;
    niter = niter + 1;
  end;

  if (^ converge) then print "No convergence";
  else run failonly (time, score, info, delta);      * subset delta = 1;

finish;

***************************************************;

* main program ->;

  use A;
  read all var {time delta group stress};


  run phreg (beta, score , info, time, delta, group);

  v = 1 / sum(info);
  se = sqrt(v);

  print beta se,  time score info;