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; 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 cn; set xb; keep d; 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 cn; read all into d ; 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) ; n=ng[+]; tpd=j(n,1,0); td=0; do j=1 to n; if j=ng[1]+1 then td=0; if d[j]=1 then do; tpd[j]=tp[j] ; td=tp[j] ; end; else tpd[j]=td ; end; * ----------------- ; 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) ; qb=j(nz,1,0) ; bh=j(nz,nz,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; * *********************** ; k=k+1; end; p1=p2+1; j1=j2+1; if g=1 then sa1=s1; if g=2 then sb1=s1; end; * ------------------------ ; j1=1 ;j2=0; rb=j(n,nrt,0); qbj=j(nrt,nz,0); do g=1 to nrg ; j2=j2 +ng[g] ; k=1; if g=1 then s1=sa1; if g=2 then s1=sb1; ch=j(nz,nrt,0) ; do j=j1 to j2 ; do l=1 to nz; ab=0; kj=1; lj=j1; do while ( x[lj,1] <= x[j,1] ) ; ab =ab+ d[lj]*(x[j ,2]/s0[g,kj])*(z[ j,l] -s1[l,kj]/s0[g,kj] ); kj=kj+1; lj=lj +1; if lj > j2 then go to skip;end; skip: qb[l]=d[j]*(z[j,l] -s1[l,k]/s0[g,k]) -ab/ng[g] ; end; qbj[k,]=qb` ; * <<<<<<<<<<<<<<< >>>>>' ; do m=j1 to j2; lm=j1; ar=0; km=1; do while ( x[lm,1] <= x[m,1] & x[lm,1]<= x[j,1]) ; ar= ar +d[lm]*x[j,2]/s0[g,km]**2 ; lm=lm+1 ; km=km+1 ; if lm >j2 then go to sof; end; sof: rb[m,k]= d[j]*(x[j,1] <= x[m,1])/s0[g,k] -ar/ng[g] ; end; * ----------------------------------------------------; IF k=1 then do; lh[g,k]= d[j]*(1.0/s0[g,k])/ng[g] ; do l=1 to nz; h[l,k]= d[j]*(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] +d[j]*(1.0/s0[g,k])/ng[g] ; do l=1 to nz; h[l,k]= h[l,k-1] +d[j]*(s1[l,k]/s0[g,k]**2 )/ng[g]; end; fh[g,k]=exp(-lh[g,k] ) ; end; * ---------------; k=k+1 ; bh= bh +qb*qb` ; end; ch1=j(nz,1,0) ; jk=1 ; do j=j1 to j2 ; do k=1 to ng[g] ; do l=1 to nz; ch[l,jk]= ch[l,jk] + qbj[k,l]*rb[j, k]/ng[g] ; end; end; if j=j1 then do l=1 to nz; ch1[l]=ch[l,1]; end ; if d[j]=1 then do l=1 to nz; ch1[l]=ch[l,jk]; end; if d[j]=0 then do l=1 to nz; ch[l,jk]=ch1[l] ; end; * ++++++++++++++++; jk=jk+1; end; j1=j2+1; if g=1 then do;cha=ch; ha=h; end; if g=2 then do; chb=ch; hb=h; end; end; * ; do k=1 to ng[1]; rb1=rb[1,k]; do j=1 to ng[1]; if d[j]=1 then rb1=rb[j,k]; if d[j]=0 then rb[j,k]=rb1 ; end; end; * ; do k=1 to ng[2]; rb1=rb[ng[1]+1,k]; do j=ng[1]+1 to n ; if d[j]=1 then rb1=rb[j,k]; if d[j]=0 then rb[j,k]=rb1 ; end; end; bh=bh/ng[+] ; ggh= n*sig*bh*sig ; print ' >> B << ' bh; print ggh; * -------------- ; ah0=j(nrg,1,0) ; p1=2 ; p2=0; do g=1 to nrg; if g=1 then do; ch=cha; h=ha; end; if g=2 then do; ch=chb; h=hb; end; pi=j(nz,1,0); phi=j(nz,1,0) ; k=2; nd=d[p1-1] ; ah0[g]=d[p1-1]*tpd[p1-1] ; do p=p1 to p2+li[g] +1 ; tpp= tu; dp=1; nd= nd +d[p] ; if k < li[g]+1 then do; tpp=tpd[p] ;dp=d[p]; end; ah0[g]= ah0[g] +dp*fh[g,k-1]*(tpp -tpd[p-1] ) ; if nd >= 2 then do l=1 to nz; phi[l]= phi[l] +dp*fh[g,k-1]*h[l,k-1]*(tpp - tpd[p-1]) ; pi[l]= pi[l] +dp*fh[g,k-1]*ch[l,k-1]*(tpp - tpd[p-1]) ; end; k=k+1; end; p2=p2 +np[g] ; p1=p2+2; if g=1 then do; phia=phi ; pia=pi; end; if g=2 then do; phib=phi ; pib=pi; end; end; * print omga; omga=j(nrg,1,0); p1=2; p2=0; do g=1 to nrg ; do j=1 to ng[g] ; k=2; nd=d[p1-1]; omj=0; do p=p1 to (p2+li[g]+1 ) ; nd=nd +d[p] ; tpp= tu; dp=1; if k < li[g]+1 then do; tpp=tpd[p] ; dp=d[p];end; if nd >=2 then omj=omj + dp*rb[p-1,j ]*fh[g,k-1]*(tpp - tpd[p-1]) ; k=k+1; end ; omga[g] = omga[g] +omj**2 ; end ; omga[g] = omga[g]/ng[g] ; p2=p2 +np[g] ; p1=p2+2; end; vma=omga[1]/ng[1] +phia`*ggh*phia - 2*phia`*sig*pia ; vmb=omga[2]/ng[2] +phib`*ggh*phib - 2*phib`*sig*pib ; v=(phia - phib)`*sig*(pia- pib) ; v= omga[1]/ng[1] +omga[2]/ng[2] +(phia-phib)`*ggh*(phia- phib) - 2*v ; 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 pia ' ' pib; print vma vmb ' ' omga ' ' v ' ' czi ; quit;