%macro homotest(data,strata,group,time,censor,factors,outdata=); data lc_dat1; set &data; %if &strata=0 %then %str(lc_str=1;);%else %str(lc_str=&strata;); %if &group=0 %then %str(lc_grp=_n_;);%else %str(lc_grp=&group;); lc_cen=&censor; lc_t=&time; keep lc_str lc_grp lc_cen lc_t &factors; proc phreg covout outest=covmat data=lc_dat1 noprint; model lc_t*lc_cen(0)=&factors; output out=lc_suv logsurv=lc_cmh resmart=lc_rm xbeta=lc_l /method=ch order=data; strata lc_str; data lc_dat2; merge lc_dat1 lc_suv; proc sort data=lc_dat2; by descending lc_t lc_cen; data lc_dat3; set lc_dat2; lc_expx=exp(lc_l); lc_baseh=-lc_cmh/exp(lc_l); keep lc_str lc_grp lc_cen lc_t lc_expx lc_baseh lc_rm &factors; proc iml; use lc_dat3; read all into x; do i=1 to ncol(x); x=x[loc(x[,i]^=.),];end; nc=ncol(x)-7; rowx=nrow(x); strindex=J(nrow(x),1,1);strn=1; do i=2 to nrow(x); true=0; do j=1 to i-1 until(true=1); if x[i,nc+1]=x[j,nc+1] then do;true=1;strindex[i]=strindex[j]; end;end; if true=0 then do;strn=strn+1;strindex[i]=strn;end; end; x=x||strindex;free strindex; iv=J(rowx,1,0); jv=J(rowx,1,0); iv[1]=1; jv[1]=1; gnum=1; true=0; do i=2 to rowx; true=0; gidx=0; do j=1 to i-1; if x[i,nc+2]=x[j,nc+2] then do; true = 1; iv[i]=iv[j]; gidx=gidx+1; end; end; if true=0 then do; gnum=gnum+1; iv[i]=gnum; end; jv[i]=gidx+1; end; gindex=J(gnum,1,0); do m=1 to rowx; gindex[iv[m]]=gindex[iv[m]]+1; end; xrow=J(gnum,max(gindex),0); do i=1 to gnum; do j=1 to gindex[i]; do m=1 to rowx; if iv[m]=i then if jv[m]=j then xrow[i,j]=m; end;end;end; free gidx iv jv; death=sum(x[,nc+3]); d=0; dindex=J(death,1,0); dnum=J(death,strn,0); do m=rowx to 1 by -1; if x[m,nc+3]=1 then if d=0 then do; d=d+1; dindex[d]=m;dnum[d,x[m,nc+8]]=dnum[d,x[m,nc+8]]+1; end; else if x[m,nc+4]>x[dindex[d],nc+4] then do; d=d+1; dindex[d]=m;dnum[d,x[m,nc+8]]=dnum[d,x[nc+8]]+1; end; else if x[m,nc+4]=x[dindex[d],nc+4] then dnum[d,x[m,nc+8]]=dnum[d,x[m,nc+8]]+1; end; start y(i,j,t) global(xrow,x,nc); if x[xrow[i,j],nc+4] >= t then return(1); else return(0); finish; start p(i,j,n,s0) global(xrow,x,nc); return(y(i,j,x[n,nc+4])*x[xrow[i,j],nc+6]/s0[n]); finish; start stra(i,j,h) global(xrow,x,nc); if x[xrow[i,j],nc+8]=h then return(1); else return(0); finish; m_=J(gnum,d,0); do i=1 to gnum; do k=1 to d; do j=1 to gindex[i];true=0; if x[xrow[i,j],nc+4] <= x[dindex[k],nc+4] then m_[i,k]=m_[i,k]+x[xrow[i,j],nc+3]-x[xrow[i,j],nc+7]*x[xrow[i,j],nc+6]; else if x[dindex[k],nc+8]=x[xrow[i,j],nc+8] then m_[i,k]=m_[i,k]-x[dindex[k],nc+7]*x[xrow[i,j],nc+6]; else do m=dindex[k] to rowx until(true=1); if x[m,nc+8]=x[xrow[i,j],nc+8] then do; true=1; m_[i,k]=m_[i,k]-x[m,nc+7]*x[xrow[i,j],nc+6]; end;end; end;end;end; cr=J(strn,1,0); hat_i=J(strn,1,0);theta=J(nc,strn,0); do h=1 to strn; p_=J(gnum,d,0);q=J(gnum,d,0); s_0=J(rowx,1,0); do n=1 to rowx; do i=1 to gnum; do j=1 to gindex[i]; s_0[n]=s_0[n]+stra(i,j,h)*y(i,j,x[n,nc+4])*x[xrow[i,j],nc+6]; end; end;end; do i=1 to gnum; do k=1 to d; do j=1 to gindex[i]; p_[i,k]=p_[i,k]+stra(i,j,h)*p(i,j,dindex[k],s_0); end;end;end; do k=1 to d; do i=1 to gnum; cr[h]=cr[h]+dnum[k,h]*p_[i,k]**2; end; end; do i=1 to gnum; do g=1 to gnum; q[i,1]=q[i,1]+p_[g,1]**2;end; q[i,1]=2*(q[i,1]-p_[i,1]); do k=2 to d; do g=1 to gnum; q[i,k]=q[i,k]+p_[g,k]*(p_[g,k]-m_[g,k-1]); end; q[i,k]=2*(q[i,k]+m_[i,k-1]-p_[i,k]); end; end; do i=1 to gnum; do k=1 to d; hat_i[h]=hat_i[h]+q[i,k]**2*p_[i,k]*dnum[k,h]; end;end; do v=1 to nc;do i=1 to gnum;do k=1 to d;do j=1 to gindex[i]; theta[v,h]=theta[v,h]+ q[i,k]*dnum[k,h]*stra(i,j,h)*p(i,j,dindex[k],s_0)*x[xrow[i,j],v]; end;end;end;end; end; free m_ p_ q; tmp=J(gnum,1,0); do i=1 to gnum; do j=1 to gindex[i]; tmp[i]=tmp[i]+x[xrow[i,j],nc+5]; end; end; score=sum(tmp#tmp)-death+sum(cr); free cr tmp; use covmat; read all into cm; var=sum(hat_i)+(theta[,+])`*cm[2:nrow(cm),1:ncol(cm)-1]*theta[,+]; out=J(1,2,0);out[1]=score;out[2]=var; create lc_outds from out [colname={'score_t' 'variance'}]; append from out; close lc_outds; quit; data lc_outda; set lc_outds; test=score_t/sqrt(variance); p_value=2*(1-CDF('NORMAL',abs(test))); proc print data=lc_outda; %if %length(&outdata) ^=0 %then %str(data=&strata;set lc_outda); %mend;