* PHREG6.SAS
* xref:
* input: SPRINGS.DAT
* output:
* does: - Illustrate the computational algorithm for partial likelihood
          in SAS IML.
        - time: no ties.
        - 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);
  g3 = (group = 3);
  g4 = (group = 4);
  run;
***************************************************;
* sort the data properly before starting IML;
proc sort data=A nodupkey; by descending time delta; run;
* make sure no tied times exist;
***************************************************;
proc phreg data=A;
  model time * delta(0) = g3 g4;
  run;
***************************************************;
proc iml worksize=1024;
reset nocenter   name pagesize=65 linesize=178;
***************************************************;

start failonly (t, x, d);

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

finish;

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

start phreg (beta, score, info, avar, time, delta, z);

  * Proportional Hazards Model Fitting
  * Expects data sorted by descending time and ascending delta.
  * time, delta:  n by 1
  * z:  n by p
  * beta <- p x 1 estimates
  * info  <- p x p info matrix
  * avar  <- p x p asymptotic var. est.
  * score <- d x p residuals     (delta=1)
  * time <- d x 1 subset (delta=1)  in desc. order
  ;

  n = nrow(z);
  p = ncol(z);
  beta  = j(p, 1, 0);
  score = j(n, p, 0);   * score process residuals;
  oldbeta = beta;
  niter = 0;

  do until (converge | (niter=10)); * main loop;
    niter = niter + 1;
    eta = z * beta;     * linear predictor;
    psi = exp( eta );   * weights;
    s0 = cusum( psi );  * total weight;
    s1 = j(p, 1, 0);    * 1st moment;
    s2 = j(p, p, 0);    * 2nd moment;
    info = s2;
    n2ll = -2 * sum(delta#(eta - log(s0)));     * -2 log Lik.;
    if (niter=1) then print '-2 LOG L without covariates: ' n2ll;
    do i = 1 to n;
      zi = t(z[i,]);
      psiizi = psi[i] * zi;
      s1 = s1 + psiizi;
      s2 = s2 + psiizi * t(zi);
      if (delta[i]) then
      do;
        zbar = s1 / (s0[i]);         * weighted average;
        score[i,] = t(zi - zbar);    * residual;
        info      = info + s2/(s0[i]) - zbar* t(zbar); * weighted var.;
      end;
    end;
    u = t( score[+,] );         * the total score vector;
    avar = inv(info);
    beta = beta + avar * u;
    converge = (abs(beta-oldbeta) < 1.0e-4);
    print niter n2ll oldbeta beta u;
    oldbeta = beta;
  end;

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

finish;

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

* main program ->;

  use A;
  read all var {time delta};
  read all var {g3 g4} into z;

  run phreg (beta, score, info, avar, time, delta, z);

  se = sqrt(vecdiag(avar));

  print beta se,  time score;