/*************** Cox Proportional Hazards  ******************/


/* coxfit

Coxfit estimates the parameters from a
Cox proportional hazards model.


Written by Dan Powers (12/17/98),
Department of Sociology
University of Texas at Austin
E-mail: dpowers@mail.la.utexas.edu

This code is written and submitted for public, non-commercial use.  This
code will reproduce results from standard packages in a timely manner but
there are no performance guarantees.  If you use this code extensively in
a research project, please cite the author.  This code can be used alone,
however, there are many canned packages that can estimate Cox proportional
hazard  models. The gauss coxfit procedure was designed to be used as part
of the M-Step in an EM-algorithm to estimate the variance in a gamma
distributed random effect as described by Lancater (1990).  A companion
program (mfEM_cox.arc) provides the gauss code to do this.


Proc COXFIT:  Gauss procedure to estimate the Cox proportional
                 hazards model using Breslow-Peto adjustment for
                 ties.

Inputs:    b  -  starting coeficient vector (kx1)
           t  -  waiting times              (nx1)
           d  -  event=1/censoring=0        (nx1)
           x  -  covariates                 (nxk)
        offs  -  offset                     (nx1)
        tol   -  convergence tolerance      (1x1)
      maxit   -  maximum iterations         (1x1)

Output
       Mod    -  column[1]=integrated hazards (Aalen's Estimate)
                 for use in computing gamma frailties.  (nx1)
              -  column[2]=original id number (after a few sort
                 procedures we need to check this).     (nx1)
         b    -  estimates                              (kx1)
         H    -  information matrix                     (kxk)

*/

proc(4)=coxfit(b,t,d,x,offs,tol,maxit);

local dta, ndta, n_obs, n_fail, n_uni, id, i,
      ii, j, jj, k, tf, tf_u, df, xf, sx, dd, iter, db,
      new_b, U, H, lam, sum_d, sum_x, mc, inc, idx, R,
      sum_id,ntf,rsdta,tidx, tlam, haz, Ihaz, risk, LogL,
      modout;

/* first concatenate the input data */

dta = t~d~x~offs;
  k = cols(x);
 id = seqa(1,1,rows(dta));

/* reverse sort by failure-time for computational ease */

ndta = rev(sortc(dta,1));
   t = ndta[.,1];
   d = ndta[.,2];
   x = ndta[.,3:3+(k-1)];
offs = ndta[.,3+k];
n_obs  = rows(ndta);  /* number of observations    */
n_fail = sumc(d);     /* number of failures/events */


/*
Some calculations work only over failure times.
*/

tf     = selif(t,d.==1);   /* failure times                            */
df     = selif(d,d.==1);   /* ones                                     */
xf     = selif(x,d.==1);   /* covariate vector                         */
tf_u   = rev(unique(tf,1));/* unique failure times                     */
n_uni  = rows(tf_u);


/*    We need some quantities based on the
      unique failure times, such as:
dd  - number of ties at each unique failure time
sx  - the sum of covariate vector for tied obs.

*/


ii=0;
 j=0;
sx =zeros(n_uni,k);
dd =zeros(n_uni,1);
jj =zeros(n_uni,1);
sum_x  = zeros(1,k);
sum_d  = 0;

i=1;
  do until i > n_fail;
        if i eq 1;
  sum_x  = xf[1,.];
  sum_d  = df[1];
           ii=1;
            j=1;
        else;
          if tf[i] ne tf[i-1];
sum_x  = xf[i,.];
sum_d  = df[i];
           ii = ii+1;
            j=1;
        elseif tf[i] eq tf[i-1];
sum_x  = sum_x+xf[i,.];
sum_d  = sum_d+df[i];
            j=j+1;
           ii=ii;
        endif;
       endif;
       jj[ii]=j;
sx[ii,.] = sum_x;
dd[ii]   = sum_d;
i=i+1;
endo;

/* Get Risk Set Indicies */

idx=zeros(n_uni,1);
 i=1;
  do until i > n_uni;
     idx[i]=sumc(ones(n_obs,1).*(t .>= tf_u[i]));
  i=i+1;
endo;


mc=maxc(jj);
"NUMBER OF OBSERVATIONS                        : ";;n_obs;?;
"NUMBER OF FAILURES                            : ";;n_fail;?;
"PERCENT CENSORED                              ; ";;1-n_fail/n_obs;?;
"NUMBER OF UNIQUE FAILURE TIMES                : ";;n_uni;?;
"MAXIMUM NUMBER OF TIES PER UNIQUE FAILURE TIME: ";;mc;?;
"TOTAL NUMBER OF TIED FAILURE TIMES            : ";;
sumc(dd.*(dd.>1));?;
"PROPORTION OF TIES TO FAILURES                : ";;
sumc(dd.*(dd.>1))./n_fail;?;


/* if many ties call FIXUP before estimating this model */

/* NEWTON-STEP */

 db  = 1;
iter = 0;
do until abs(db) < tol or iter > maxit;

{lam,U,H,R} = CALC(b,x,sx,offs,k,n_uni,dd,idx);

    new_b = b  + solpd(U,H);
       db = new_b-b;
        b = new_b;
iter=iter+1;
endo;

"COX MODEL CONVERGED WITH TOLERANCE ";; format /rd (5,5); meanc(abs(db));;
"IN ";;format /rd (5,0);iter;;"ITERATIONS";
format /rd (12,8);

/*
   These are quantites to be output along with results.
   We need to get the times into ascending order to make
   this come out right.

   Sort by time and bring id and x along.  Next, sort by
  unique time, bring along estimated hazard.
*/

tidx  = sortc(t~d~id~x,1);
tlam  = sortc(tf_u~lam,1);

/*
   Get quantities needed to compute Aalen's estimate of the
   Integrated Hazard. We will also need this to compute the
   likelihood function for this model.
*/

risk  = exp(tidx[.,4:4+k-1]*b);

/* Loop over ALL observation to get Hazard and Cumulative Hazard */

Ihaz=zeros(n_obs,1);
 haz=zeros(n_obs,1);

i=1;
   do until i > n_obs;
          haz[i] = risk[i].*maxc(tlam[.,2].*(tlam[.,1] .== tidx[i,1]));
         Ihaz[i] = risk[i].*sumc(tlam[.,2].*(tlam[.,1] .<= tidx[i,1]));
   i=i+1;
endo;

/*
   Bookkeeping - rearrange by original id number so that Ihaz[i]
                 is matched up with the right observation in the
                 original data.
*/

LogL = sumc(sx*b - dd.*ln(R));

modout=sortc(Ihaz~haz~id,3);

"(Partial) Log-Likelihood ";;LogL;

retp(modout[.,1],modout[.,2],H,b);
endp;



/* PROC CALC:  risk-set loop to do main work */


proc(4)=CALC(b,x,sx,offs,npar,n_uni,dd,idx);

  local i,ii,j,xi, A, rf, rfx, rfxpx, xpx,
        rr,rx, rxx, U ,H ,lambda;

       xi = exp(x*b + ln(offs) );
    rf    = cumsumc(xi);
    rfx   = cumsumc(x.*xi);
    rfxpx = cumsumc(x*~x.*xi);


/* score and hessian */

rx=zeros(n_uni,npar);
rxx=zeros(n_uni,npar*npar);
rr=zeros(n_uni,1);

/* Score and Information  */

i = 1;
  do until i > n_uni;
ii=idx[i];
rx[i,.] = rfx[ii,.];
rxx[i,.]= rfxpx[ii,.];
rr[i]   = rf[ii];
i=i+1;
endo;

U = sumc(sx - dd.*rx./rr);

xpx=zeros(n_uni,npar*npar);

A=(rx./rr)*~(rx./rr);

j=1;
  do until j > npar*npar;
    xpx[.,j]=dd.*(rxx[.,j]./rr - A[.,j]);
   j=j+1;
  endo;

H=reshape(sumc(xpx),npar,npar);

/*
  Lambda is the estimated baseline hazard at each failure time -
  to be used to calculated the Integrated hazard and Log L.
*/

lambda = dd./rr;

retp(lambda,U,H,rr);
endp;

/*
   PROC fixup - replaces tied failure times with time + a small
   increment.  For use in situations with extreme ties.
   Another model would be preferred.
*/


Proc(1)=fixup(t);
  local n, nt, i, inc, ties;
n=rows(t);
nt=t;
ties=0;
   i=1;
 inc = .001;
     do until i > n;
           if t[i] ne t[i-1]; nt[i]=t[i];
       elseif t[i] eq t[i-1]; nt[i]=nt[i]+inc;
   ties=ties+1;
 endif;
   inc= inc + .001;
i=i+1;
endo;
"Ties " ;ties;
retp(nt);
endp;






