Logo Search packages:      
Sourcecode: survival version File versions

coxscho.c

/*  SCCS @(#)coxscho.c  5.2 10/27/98
** Return the Schoenfeld residuals.
**
**  the input parameters are
**
**       nused        :number of people
**       nvar         :number of covariates
**       y(3,n)       :start, stop, and event for each subject
**       covar(nv,n)  :covariates for person i.
**                        Note that S sends this in column major order.
**       strata(n)    :marks the strata.  Will be 1 if this person is the
**                       last one in a strata.  If there are no strata, the
**                       vector can be identically zero, since the nth person's
**                       value is always assumed to be = to 1.
**       score(n)     :the risk score for the subject
**       method       : =1 if the Efron method was used
**
**  returned parameters
**       covar        :for each death, the row is replaced with the residuals
**
**  work arrays
**       a(nvar)
**       a2(nvar)
**       mean(nvar)
**
**  the 3 vectors a, a2 and mean are passed as a single
**    vector of storage, and then broken out.
**
**  the data must be sorted by ascending time within strata, deaths before
**          living within tied times.
*/
#include <math.h>
#include "survS.h"
#include "survproto.h"

void coxscho(int   *nusedx,    int   *nvarx,    double *y, 
           double *covar2,    double *score,    int   *strata,  
           int   *method2,   double *work)
{
    int i,k,person;
    int     nused, nvar;
    double **covar;
    double *a;
    double *a2;
    double *mean;
    double  denom, weight;
    double  time;
    double  temp;
    double     method;
    double   deaths;
    double efron_wt;
    double  *start,
          *stop,
          *event;

    nused = *nusedx;
    nvar  = *nvarx;
    method= *method2;
    /*
    **  Set up the ragged arrays
    */
    covar= dmatrix(covar2, nused, nvar);
    a = work;
    a2= a+nvar;
    mean = a2+nvar;
    start =y;
    stop  =y + nused;
    event =y + nused +nused;

    /*
    ** Now walk through the data
    */
    for (person=0; person<nused;) {
      if (event[person]==0) person++;
      else {
          /*
          ** compute the mean over the risk set and over the deaths (a & a2)
          */
          denom =0;
          efron_wt =0;
          for (i=0; i<nvar; i++) {
            a[i] =0;
            a2[i]=0;
            }
          time = stop[person];
          deaths=0;
          for (k=person; k<nused; k++) {
            if (start[k] < time) {
                weight = score[k];
                denom += weight;
                for (i=0; i<nvar; i++) {
                  a[i] += weight*covar[i][k];
                  }
                if (stop[k]==time && event[k]==1) {
                  deaths += 1;
                  efron_wt += weight*event[k];
                  for (i=0; i<nvar; i++) a2[i]+= weight*covar[i][k];
                  }
                 }
            if (strata[k]==1) break;
            }

          /*
          ** Compute the mean at this time point
          */
          for (i=0; i<nvar; i++) mean[i] =0;
          for (k=0; k<deaths; k++) {
            temp = method *k/deaths;
            for (i=0; i<nvar; i++)
                mean[i] += (a[i] - temp*a2[i])/(deaths*(denom -temp*efron_wt));
            }
          /*
          ** Compute the residual(s) for this time point
          */
          for (k=person; k<nused && stop[k]==time; k++) {
            if (event[k]==1) {
                for (i=0; i<nvar; i++) {
                  covar[i][k] -= mean[i];
                  }
                }
            person++;
            if (strata[k]==1) break;
            }
          }
      }
    }

Generated by  Doxygen 1.6.0   Back to index