options nocenter  ps=500  ls=80;
title ' Carcinoma of the oropharynx' ;
data org; infile 'karri.dat';
input case inst sex group grade age cond site t n ed status time ;
if sex=1 then sx= 0; if sex=2 then sx=1 ;
if grade=9 then grade=2;
if cond=0 or cond=9 then cond=1;
s1=0; s2=0;
if site=2 then s1=1;  if site=4 then s2=1 ;
 rename status=d   sx=z1 grade=z2 age=z3 cond=z4 s1=z5 s2=z6 t=z7
         n=z8 time=t   ;
data org; set org;
Z1 = z1 - 0.2358974;
Z2 = z2 - 1.9282051;
Z3 = z3 - 60.4358974;
Z4 = z4 - 1.2974359;
Z5 = z5 - 0.3282051;
Z6 = z6 - 0.3384615;
Z7 = z7 - 3.1179487;
Z8 = z8 - 1.9230769;
/* proc ttest;class group; var age;
proc freq; tables GROup*(sx grade cond site t n status)/nopercent
 nocol norow;
endsas;
data xb  ; infile 'zoker.dat';
input group t d z1 z2;
ebz=z1+z2;
*/
* ------------------  ;
*  varible t = response (failure time) .
   varible d = censor  .
   variables z1 z2 ...   are covariates .   ;
*  -------------------------------------------  ;

  proc phreg data=org covout outest=cov ;
  model t*d(0)=z1 -z8       ;
    strata group;
 Id group ;
output out=xb xbeta=xb      ;
*
     THE names XB and COV are nessary .
   @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@   ;
data xb; set xb;    * xb=xb -0.074909*group;
ebz=exp(xb);

  proc sort data=xb; by group t ;
* proc print data=cov;
data tp; set xb ;
  if d=1;
keep group t;
/*proc sort data=tp ;by group t; */
*proc print;
proc  means data=tp noprint ; by group;
output out=np n=n ;
data NP;set NP; keep n;
*proc print ;
proc means data=xb  noprint; by group;
output out=ng n=n ;
data NG;set NG; keep n;
*proc print ;
DATA x; set xb ;    keep t ebz  ;
DATA z;set xb ; keep z1 -z8 ;
  data cov;set cov ; if _type_='COV  ' ; keep z1 -z8;
    data cov;set cov(obs=8);
/*
data cov;input z1 z2;cards;
4 0
0 4
;
  */
proc IML;
use tp ; read all var { t } into tp ;
use NP ; read all into np ;
use NG ; read all into ng ;
use  X ; read all into x ;
use  Z ; read all into z ;
use cov; read all into sig;
 print      np ng    ;
print sig;
* &&&&&&&&&&&&&&&&&&  ;
*tu=4.99 ; tu=1277.5 ;               *    ++++++++++++++++++++++++++++;
* &&&&&&&&&&&&&&&&&&  ;
nz=ncol(z) ;
nrg=nrow(np) ;      nrt=max(np) ;
s0=j(nrg,nrt,0) ;
fh=j(nrg,nrt,0) ;
lh=j(nrg,nrt,0) ;
vh=j(nrg,nrt,0) ;
 h =j(nz ,nrt,0) ;
li=j(nrg,1,0) ;
p1=1    ; p2=0;
j1=1      ;j2=0;
  do   g=1 to nrow(ng) ;
*     ----------------------- ;
p2=p2+np[g] ;
j2=j2 +ng[g] ;
  k=1;
s1=j(nz,nrt,0) ;
      do   p=p1 to p2    ;
If (tp[p] <= tu) then li[g]= li[g] +1;
*                                        ;
  do    j=j1 to j2 ;
   s0[g,k]= s0[g,k] + x[j,2]*(   x[j,1] >= tp[p] ) ;
         do  l=1 to nz;
   s1[l,k]= s1[l,k] + x[j,2]*z[j,l]*(   x[j,1] >= tp[p] ) ;
                   end;
          end;
   s0[g,k]=s0[g,k]/ng[g];
do l=1 to nz;   s1[l,k]=s1[l,k]/ng[g]; end;
*  ***********************   ;
IF k=1 then do;
lh[g,k]= (1.0/s0[g,k])/ng[g] ;
vh[g,k]=(1.0/s0[g,k]**2 )/ng[g];
do l=1 to nz; h[l,k]= (s1[l,k]/s0[g,k]**2 )/ng[g];  end;
fh[g,k]=exp(-lh[g,k] ) ;
        end;
else;  do;
lh[g,k]=  lh[g,k-1] +(1.0/s0[g,k])/ng[g] ;
vh[g,k]= vh[g,k-1] +(1.0/s0[g,k]**2 )/ng[g];
do l=1 to nz; h[l,k]=  h[l,k-1] +(s1[l,k]/s0[g,k]**2 )/ng[g]; end;
fh[g,k]=exp(-lh[g,k] ) ;       end;
k=k+1;
     end;
 p1=p2+1;  j1=j2+1;
if g=1 then do;sa1=s1; ha=h; end;
if g=2 then do; sb1=s1;  hb=h; end;
 end;
                     *print fh;
   ah0=j(nrg,1,0) ;
omga=j(nrg,1,0);
p1=2 ; p2=0;
  do g=1 to nrg;
if g=1 then h=ha;
if g=2 then h=hb;
phi=j(nz,1,0) ;
k=2;
       ah0[g]=tp[p1-1] ;
     do  p=p1 to p2+li[g] +1 ;
 tpp= tu;
if k < li[g]+1 then tpp=tp[p] ;
ah0[g]= ah0[g] +fh[g,k-1]*(tpp   -tp[p-1] ) ;
 omga[g]= omga[g]+ vh[g,k-1]*fh[g,k-1]**2*(tpp-tp[p-1])**2;
           do l=1 to nz;
  phi[l]= phi[l] +fh[g,k-1]*h[l,k-1]*(tpp - tp[p-1]) ;  end;
  k=k+1;
  end;
p2=p2 +np[g] ;
p1=p2+2;
if g=1 then phia=phi ;
if g=2 then phib=phi ;
end;
              *      print omga;
p1=2; p2=0;
     do g=1 to nrg ;
 k=2;
       do p=p1 to (p2+li[g] )  ;
 kq=k+1;
         do q=(p +1) to (p2+li[g]+1) ;
 tpq= tu;
if kq< li[g]+1 then tpq=tp[q] ;
  omga[g]=omga[g] + vh[g,k -1]*fh[g,k-1]*fh[g,kq-1]*
                        ( tp[p] -tp[p-1])*(tpq -tp[q-1])*2.0 ;
 kq=kq+1;
    end ;    k =k +1;     end;

p2=p2 +np[g] ;
p1=p2+2;
       end;
 v=(phia - phib)`*sig*(phia- phib)+omga[1]/ng[1] +omga[2]/ng[2] ;
   czi=(ah0[1] -ah0[2] )/sqrt(v) ;
*  ------------- ;
print tu li;
*print s0;
*print sa1, sb1;
*print lh;
*print fh;
*print vh;
*print  ha,hb ;
print ah0;
print phia,phib;
print omga   v czi ;
quit;