%macro incid(data,group,relp,trm,time,out=); /* arguments are: dataset name, group variable(1,2,...), relapse indicator(1-yes, 0-no) trm indicator(1-yes, 0-no), dataset name for relapse incidence dataset name for trm incidence */ data lc_one; set &data; a=&time; b=&relp; c=&trm; d=&group; keep a b c d; proc sort data=lc_one; by descending a; /*proc sort data=lc_one; by d;*/ proc iml; use lc_one; read all into x; n=nrow(x); ngrp=max(x[ ,4]); gnum=J(1,ngrp,0); do k=1 to n; gnum[1,x[k,4]] = gnum[1,x[k,4]]+1; end; t=J(1,n,0); t[1,1]=x[1,1]; ntime=1; tnow=x[1,1]; do j=2 to n; if x[j,1] = t[1,k] then atrisk[k,ag]=atrisk[k,ag]+1; end; end; lfs=J(ntime,ngrp,-1); ci_rel=J(ntime,ngrp,-1); ci_trm=J(ntime,ngrp,-1); vci_rel=J(ntime,ngrp,-1); vci_trm=J(ntime,ngrp,-1); index=J(ntime,1,0); tt=t(t[1,1:ntime]); do ig=1 to ngrp; p=1; cr=0; cd=0; do j=1 to ntime ; index[j,1]=j; k=ntime-j+1; if atrisk[k,ig] >0 then do; cr=cr+relap[k,ig]*p/atrisk[k,ig]; cd=cd+trm[k,ig]*p/atrisk[k,ig]; p=p*(1-(trm[k,ig]+relap[k,ig])/atrisk[k,ig]); lfs[k,ig]=p; ci_rel[k,ig]=cr; ci_trm[k,ig]=cd; end; else do; lfs[k,ig]=.; ci_rel[k,ig]=.; ci_trm[k,ig]=.; end; end; do j=1 to ntime; know=ntime-j+1; vr=0; vd=0; if ci_rel[know,ig] = . then do; vci_rel[know,ig]=.; vci_trm[know,ig]=.;end; else do; do k=1 to ntime; jnow=ntime-k+1; if tt[jnow,1] <= tt[know,1] then do; wr=(trm[jnow,ig]+relap[jnow,ig])/atrisk[jnow,ig]**2; wr=wr*lfs[jnow,ig]*(ci_rel[know,ig]-ci_rel[jnow,ig])**2; q=(relap[jnow,ig]/atrisk[jnow,ig]**2)*lfs[jnow,ig]**2; q=q*(1-2*(ci_rel[know,ig]-ci_rel[jnow,ig])); vr=vr+wr+q; wd=(trm[jnow,ig]+relap[jnow,ig])/atrisk[jnow,ig]**2; wd=wd*lfs[jnow,ig]*(ci_trm[know,ig]-ci_trm[jnow,ig])**2; q=trm[jnow,ig]/atrisk[jnow,ig]**2*lfs[jnow,ig]**2; q=q*(1-2*(ci_trm[know,ig]-ci_trm[jnow,ig])); vd=vd+wd+q; end; end; end; vci_rel[know,ig]=sqrt(vr); vci_trm[know,ig]=sqrt(vd); end; end; nn=ngrp*ntime; yout=j(nn,6,0); k=0; do is=1 to ngrp; do it=1 to ntime; k=k+1; yout[k,1]=tt[it,1]; yout[k,2]=is; yout[k,3]=ci_rel[it,is]; yout[k,4]=vci_rel[it,is]; yout[k,5]=ci_trm[it,is]; yout[k,6]=vci_trm[it,is]; end; end; create dout from yout; append from yout; close dout; quit; data io; set dout; time=col1; group=col2; CI1=col3; SE_CI1=col4; CI2=col5; SE_CI2=col6; if CI1 =. then se_ci1=.; if CI2=. then se_ci2=.; drop col1-col6; proc sort data=io; by time; proc sort data=io; by group; proc print data=io; data &out; set io; %mend;