********************************************************************************; * Program: logrank_macro.sas; * Purpose: To perform a clustered logrank test; * input: /usr/path/simdata.sas7bdat; * definitions: datain = input datset; * group = variable for assigned treatment group; * cluster = cluster id; * censor = indicator for event=1 or censor=0; * time = observation time until event or censor; * * Copyright (C) 2009 Margaret Stedman; * * This program is free software: you can redistribute it and/or modify; * it under the terms of the GNU General Public License as published by; * the Free Software Foundation, either version 3 of the License, or; * any later version.; * * This program is distributed in the hope that it will be useful,; * but WITHOUT ANY WARRANTY; *without even the implied warranty of; * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the; * GNU General Public License for more details:; * ; ******************************************************************************; ******************************************************************************; ******Optional: Enter directory for input dataset here:*********************************; *libname ll "usr\path"; options nofmterr; %macro clusterlr (datain, time, cluster, group, censor); data one; set &datain; keep &cluster &group &time &censor; run; \**************** step 1 *******************************\; proc sql; select &group, count(distinct &cluster) as groupc label="Clusters", count(*) as groupn label="Observations" into :groupv1 - :groupv2, :groupc1 - :groupc2, :groupn1 - :groupn2 from one group by &group; quit; %let groups2= %eval(&groupn1 + 1); %let totaln = %eval(&groupn1 + &groupn2); %let totalc= %eval(&groupc1 + &groupc2); \************* step 2 *********************************\ proc sort data=one; by decending &time; run; \************ step 3 *********************************\ data cumm_y; set one; by descending &time ; if _n_=1 then do; yy=0; y_g1=0; y_g2=0; end; if first.&time then do; dd=0; d_g1=0; d_g2=0; end; retain yy dd y_g1 y_g2 d_g1 d_g2; yy=yy+1; if &group=&groupv1 then y_g1=y_g1+1; if &group=&groupv2 then y_g2=y_g2+1; if &censor = 1 then do; dd=dd+1; if &group=&groupv1 then d_g1=d_g1+1; if &group=&groupv2 then d_g2=d_g2+1; end; wt=y_g1*y_g2/yy; time= &time; if last.&time then output; keep time yy dd y_g1 y_g2 d_g1 d_g2 wt; run; proc sql; create table teststat as select sum((dd * wt**2) / y_g1 / yy) as vi_1, sum((dd * wt**2) / y_g2 / yy) as vi_2, sum(d_g1 * wt / y_g1 ) as test_1, sum(d_g2 * wt / y_g2 ) as test_2 from cumm_y; \************* step 4 ******************************************\ create table newoned as select a.time, a.wt, a.y_g1, a.y_g2, a.yy, b.&group as group, b.&cluster as cluster, b.&censor as censor, 1/a.y_g1 as y_g1_inv, 1/a.y_g2 as y_g2_inv, 1/a.yy as yy_inv from cumm_y a, one b where a.time=b.&time; quit; proc sort data=newoned; by group cluster; run; \*********** step 5 *********************************************\ proc iml; use newoned var{time censor yy_inv y_g1 y_g2 y_g1_inv y_g2_inv wt}; read all into newone; \*********** step 5.1 ******************************************\ A=newone[1:&groupn1 , {2 8 6}]; B=newone[&groups2 : &totaln , {2 8 7}]; C= A//B; D=C[,#]; use newoned var{cluster}; read all into csize; %do r=1 %to &totalc; LRgk_pt1p=D[loc(csize=&r)][+]; LRgk_pt1 = LRgk_pt1//LRgk_pt1p; %end; free A B C D; \************** step 5.2****************\ \************** step 5.2.1**************\; A=newone[,{2 4}]; B=newone[,{2 5}]; C=A[,#]; D=B[,#]; E=repeat(C,1,&groupc1 ); F=repeat(D,1,&groupc2 ); G= E||F; events=choose(G=0,0,1); free A B C D E F G; \************ step 5.2.2 ***************\ A=newone[,1]; B=repeat(A,1,&totaln); C=repeat(A`, &totaln ,1); D=(C>=B); %do s=1 %to &totalc; riskp=D[1:&totaln,loc(csize=&s)][,+]; risk= risk||riskp; %end; free A B C D; \*********** step 5.2.3 ******************\ A=newone[ ,{8 6 3}]; B=newone[ ,{8 7 3}]; C=A[,#]; D=B[,#]; E=repeat(C,1,&groupc1); F=repeat(D,1,&groupc2); weight=E||F; free A B C D E F; \********* step 5.2.4***********************\ combp2=events # risk # weight; sumpt2=combp2[+,]; LRgk_pt2= sumpt2`; \********** step 5.3 *********************\ LRgk=LRgk_pt1 - LRgk_pt2; group1=j(&groupc1 ,1,1); group2=j(&groupc2 ,1,2); group=group1//group2; cluster1=1: &groupc1; cluster2=1: &groupc2; cluster=cluster1`//cluster2`; fin=group || cluster || LRgk_pt1 || LRgk_pt2 || LRgk; create clustv_ijk from fin; append from fin; quit; data fin; set clustv_ijk; rename col1=group col2=cluster col3=LRgk_pt1 col4=LRgk_pt2 col5=LRgk; run; \*********** step 5.4 **************************\ ods listing; proc sql; create table ek as select group, sum(LRgk**2) as vk from fin group by group; create table eachtreat as select sum(vk) as clustered_v from ek; \************ step 5.5 ********************************\ title "final result"; select (a.vi_1 + a.vi_2) as ind_var label="independent variance", b.clustered_v label="clustered variance", (a.test_2-a.test_1) as diff label="logrank statistic",(calculated diff/sqrt(clustered_v)) as test label="clustered logrank test", (1-probnorm(abs(calculated test)))*2 as p label="clustered logrank p-value" from teststat a, eachtreat b; quit; %mend; *%clusterlr (datain=ll.simdata, time=time, cluster=drid, group=treat, censor=censor); quit;