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   ;
/* 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 ;
* 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);
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=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) ;
gm=j(nrg,nrt,0) ;
vh=j(nrg,nrt,0) ;
li=j(nrg,1,0) ;
 h =j(nz ,nrt,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) ;
chi=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;
  do j = 1 to nrow(x) ;
e1=exp(-lh[g,k]*x[j,2] ) ;
fh[g,k]=fh[g,k] + e1;
gm[g,k]=gm[g,k] +x[j,2]*exp(-lh[g,k]*x[j,2]) ;
    do l=1 to nz; chi[l,k]=chi[l,k] +z[j,l]*x[j,2]*e1/ng[+] ;end;
  end;
   fh[g,k]=fh[g,k]/ ng[+]  ;
   gm[g,k]=gm[g,k]/ng[+] ;
 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;
  do j = 1 to nrow(x) ;
e1=exp(-lh[g,k]*x[j,2] ) ;
fh[g,k]=fh[g,k] + e1;
gm[g,k]=gm[g,k] +x[j,2]*e1 ;
    do l=1 to nz; chi[l,k]=chi[l,k] +z[j,l]*x[j,2]*e1/ng[+] ;end;
end;
   fh[g,k]=fh[g,k]/ ng[+]  ;
   gm[g,k]=gm[g,k]/ng[+] ;
                             end;
k=k+1;
     end;
 p1=p2+1;  j1=j2+1;
if g=1 then do;sa1=s1; ha=h; chia=chi; end;
if g=2 then do; sb1=s1;  hb=h; chib=chi ;end;
 end;
   ah0=j(nrg,1,0) ;
omga=j(nrg,1,0);
p1=2 ; p2=0;
  do g=1 to nrg;
if g=1 then do ;h=ha;  chi=chia; end;
if g=2 then do; h=hb;  chi=chib; end;
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]*gm[g,k-1]**2*(tpp-tp[p-1])**2;
           do l=1 to nz;
phi[l]= phi[l] +(gm[g,k-1]*h[l,k-1]-lh[g,k-1]*chi[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]*gm[g,k-1]*gm[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 ,gm ,chia,chib ;
*print vh;
*print  ha,hb ;
print ah0;
print phia '      ' phib;
print omga '   '   v '   ' czi ;
quit;