/*  SAS code for using lifereg to fit a
    Weibull model to assess the association
    between HPA staining and survial times of 
    breast cancer patients.  This code is 
    posted as  weibull2.sas  */

data set1;
  input d y x;
  datalines;
  1   5  1
  1   8  1
  1  10  1
  1  13  1
  1  18  1
  1  24  1
  1  26  1
  1  26  1
  1  31  1
  1  35  1
  1  40  1
  1  41  1
  1  48  1
  1  50  1
  1  59  1
  1  61  1
  1  68  1
  1  71  1
  0  76  1
  0 105  1
  0 107  1
  0 109  1
  1 113  1
  0 116  1
  1 118  1
  1 143  1
  0 154  1
  0 162  1
  0 188  1
  0 212  1
  0 217  1
  0 225  1
  1  23  0
  1  47  0
  1  69  0
  0  70  0
  0  71  0
  0 100  0
  0 101  0
  1 148  0
  1 181  0
  0 198  0
  0 208  0
  0 212  0
  0 224  0
  run;

  proc lifereg covout outest=set2;
    model y*d(0)= x / dist=weibull covb itprint
                    maxiter=30 converge=.0001;
    output out=set3 censored=c cdf=cdf predicted=p
           quantiles = .1 to .9 by .1 std=std;
   run;

proc print data=set2; run;

proc print data=set3; run; 


/* Keep only the columns containing the
   estimates and the covariance matrix */

data set2; set set2;
 drop _MODEL_ _NAME_ _TYPE_ _DIST_ _STATUS_ _LNLIKE_ y;
 run;


proc iml;
start btrans;

use set2;       /* Enter the data  */
read all into zz;

n=ncol(zz);  /* Number of parameters */
nr=nrow(zz);

b=zz[1, ];   /* Parameter estimates */

v = zz[2:nr, ];  /* Covariance matrix */

/* transform to parameters for weibull
   model in the notes */

a = - b / b[1,n];
a[1,n] = 1/b[1,n]; 
a[1,1] = exp(-b[1,1]/b[1,n]); 


/* Apply delta method to get a large 
   sample covraiance matrix for a */
G = diag(j(1,n,-a[1,n])) ;
G[1,1] =  -a[1,1]*a[1,n];
G[ ,n] =  -t(a)*a[1,n];
G[1,n] = -a[1,1]*a[1,n]*log(a[1,1]);
G[n,n] =  -(a[1,n])##(2);

parms = t(a);
covparms = G*V*t(G);
stdparms = sqrt(vecdiag(covparms));
print parms stdparms, covparms;

finish;

run btrans;