/**************************************************************************************** ** PROGRAM: rankcorr_mmer.sas ** ** PURPOSE: This program calculates observed and measurement-error corrected ** Spearman rank correlation coefficients and 95% confidence ** intervals. ** ** REQUIREMENT: A SAS data with the following variables: ** 1. id ** 2. Variable with surrogate exposure, such as ffq. ** 3. Goldstandard exposure such as diet records. Goldstandard variables ** can have repeated measure. They must be named with ordering ** number attached, for example, vitc1, vitc2, vitc3, vitc4. Also ** they must have the same number of repeated measure for each ** subject or the program will stop. ** ** MACRO VARIABLES: ** 1. Name of the SAS data set. ** 2. Name of the surrogate exposure. ** vitc1, vitc2, vitc3 and vitc4, enter vitc. ** 3. Number of repeated measure for goldstandard. For example, ** for vitc1, vitc2, vitc3, and vitc4, enter 4. ** 4. The name of the data set and the names of the variables entered ** in as macro variables (data=, ffq=, dr=) must be consistent with ** the names stored in the SAS data set. ** ** USAGE: Enter "%rankcorr(data=drdat, ffq=vc, dr=vitc, repeated=4 ) ** ********************************************************************************************/ %macro rankcorr_mmer(data=, ffq=, dr=, repeated=); options ps=72 nosource; libname in "."; %let pi=3.1415926; data check; set in.&data end=eof; /*check if all dr have same number of repeated measure*/ %do i=1 %to &repeated; retain ndr&i 0; if &dr&i ne . then ndr&i=ndr&i+1; %end; retain nffq 0; if &ffq ne . then nffq=nffq+1; if eof; data check; set check; array kndr (*) nffq ndr2-ndr&repeated; kk=0; do i=1 to dim(kndr); if ndr1-kndr(i)^=0 then kk=1; end; call symput('flag',kk); run; %if &flag eq 1 %then %do; data _null_; put '************************************************************************'; put '*** PROGRAM STOPPED: ***'; put '*** THE NUMBER OF THE SURROGATE EXPOSURE AND THE NUMBER OF RETEADTED ***'; put '*** GOLD STANDARD MEASURE MUST BE THE SAME ***'; put '************************************************************************'; %goto out; %end; %if &flag eq 0 %then %do; data ppp; set in.&data; proc means noprint; var &dr.1-&dr&repeated; output out=results mean=mean&dr.1-mean&dr&repeated std =std&dr.1-std&dr&repeated min =min&dr.1-min&dr&repeated max =max&dr.1-max&dr&repeated; /*calculate correrations between gold standard and each repeated measure*/ data allout; data _null_; call symput('num','1'); %macro process(num); data ppp numm(keep=nn); set ppp; array intake(nn) &dr.1-&dr&repeated; nn=# h=intake; hstar=&ffq; proc rank data=ppp out=rank; var h hstar; ranks rh rhstar; data getn; set rank end=eof; retain n 0; n=n+1; if eof then call symput("n",trim(left(n))); data posa; set rank; hhat=probit(rh/(&n+1)); hhstar=probit(rhstar/(&n+1)); one=1; proc corr pearson outp=corrs1 noprint; var hhat hhstar; proc corr pearson data=posa outp=corrs2 noprint; var rh rhstar; data corrs1; set corrs1; one=1; if _TYPE_='CORR' and _NAME_='hhat'; r=hhstar; keep one r; data corrs2; set corrs2; one=1; if _TYPE_='CORR' and _NAME_='rh' ; rs=rhstar; keep one rs; data final; merge corrs1 corrs2 end=eof; by one; zhat=0.5*log((1+r)/(1-r)); z1=zhat-1.96/sqrt(&n-3); z2=zhat+1.96/sqrt(&n-3); r1=(exp(2*z1)-1)/(exp(2*z1)+1); r2=(exp(2*z2)-1)/(exp(2*z2)+1); rs1=6/&pi*arsin(r1/2); rs2=6/&pi*arsin(r2/2); ps=6/&pi*arsin(r/2); data allout; set allout final; data setnum; set numm; nn=nn+1; call symput('num',nn); stop: %mend process; %macro macgen(repeat); %do i=1 %to &repeat; %str(%process); %end; %mend macgen; %macgen(&repeated) data outprt; set allout end=eof; if r=. then delete; keep one ps rs1 rs2; /*calculate correlation between gold standard and mean of repeated measure*/ data ppp; set ppp; proc rank data=ppp out=rank1; var &dr.1-&dr&repeated &ffq; ranks rh1-rh&repeated rhstar; data getn; set rank1 end=eof; retain n 0; n=n+1; if eof then call symput("n",trim(left(n))); data posa; set rank1 end=eof; %do i=1 %to &repeated; hhat&i=probit(rh&i/(&n+1)); %end; array hhat(*) hhat1-hhat&repeated; hhstar=probit(rhstar/(&n+1)); hhatmean=mean(of hhat(*)); proc corr pearson outp=corrs1 noprint; var hhatmean hhstar; data corrs1; set corrs1; one=1; if _TYPE_='CORR' and _NAME_='hhatmean'; mr=hhstar; data mfinal; set corrs1; mzhat=0.5*log((1+mr)/(1-mr)); mz1=mzhat-1.96/sqrt(&n-3); mz2=mzhat+1.96/sqrt(&n-3); mr1=(exp(2*mz1)-1)/(exp(2*mz1)+1); mr2=(exp(2*mz2)-1)/(exp(2*mz2)+1); mrs1=6/&pi*arsin(mr1/2); mrs2=6/&pi*arsin(mr2/2); mrs=6/&pi*arsin(mr/2); keep one mrs mrs1 mrs2; /*taking measurement error into account*/ data tfinal; set posa; array hh hhat1 hhat2 hhat3 hhat4; do over hh; hhat=hh; output; end; proc sort; by id; proc glm outstat=outst noprint; model hhat=id; absorb id; data stat; set outst end=eof; retain modelms errorms modeldf modelss errordf errorss; if _SOURCE_='ERROR' and _TYPE_='ERROR' then do; errordf=df; errorss=ss; errorms=errorss/errordf; end; if _n_=2 and _SOURCE_='id' and _TYPE_='SS1' then do; modeldf=df; modelss=ss; modelms=modelss/modeldf; end; one=1; if eof then output; data all; merge corrs1 stat; by one; within=errorms; betwen=(modelms-errorms)/&repeated; lambda=within/betwen; ri=betwen/(betwen+within); ro=mr; rt=ro*sqrt(1+lambda/&repeated); aa=(1-ro**2)**2/(&n*ro**2); bb=(1-ri)**2*(1+(&repeated-1)*ri)**2; cc=2*ri**4*(&repeated+lambda)**2*&repeated*(&repeated-1)*(&n-1); se=rt/(1-rt**2)*sqrt(aa-bb/cc); zt=0.5*log((1+rt)/(1-rt)); z1=zt-1.96*se; z2=zt+1.96*se; p1=(exp(2*z1)-1)/(exp(2*z1)+1); p2=(exp(2*z2)-1)/(exp(2*z2)+1); pst=6/&pi*arsin(rt/2); pst1=6/&pi*arsin(p1/2); pst2=6/&pi*arsin(p2/2); keep one errorms betwen lambda pst pst1 pst2 ri; /*calculations end*/ /*following part is to make a table*/ data transresults; set results; one=1; length num $4; %do i=1 %to &repeated; num=&i; rmean&dr = round(mean&dr.&i); rstd&dr = round(std&dr.&i); rmin&dr = round(min&dr.&i); rmax&dr = round(max&dr.&i); keep one num rmean&dr rstd&dr rmin&dr rmax&dr; output; %end; data left(drop=one); merge transresults outprt; by one; array psrs ps rs1 rs2; do over psrs; psrs=round(psrs,.001); end; data right(drop=one ln rn hf); merge mfinal all; by one; array mrspst mrs1 mrs2 mrs errorms betwen pst pst1 pst2 lambda ri; do over mrspst; mrspst=round(mrspst,.001); end; /*rename so they align with ps,rs1 and rs2*/ rename mrs=ps mrs1=rs1 mrs2=rs2; length ln $1; length rn $2; length hf $1; length num $4; ln=1; rn=&repeated; hf='-'; num=compress(ln||hf||rn); options ls=256 nocenter; proc sql; select num label='No' format=$char4., rmean&dr label='Mean' format=6., rstd&dr label='SD' format=4., rmin&dr label='Lower Range' format=6., rmax&dr label='Upper Range' format=6., ps label='Observed Spearman rank correlation' format=5.3, rs1 format=5.3 label='95% CI (lower)', rs2 format=5.3 label='95% CI (upper)' from left outer union corr select num label='No' format=$char4., ps label='Observed Spearman rank correlation' format=5.3, rs1 format=5.3 label='95% CI (lower)', rs2 format=5.3 label='95% CI (upper)', pst format=5.3 label='Corrected Spearman rank correlation', pst1 format=5.3 label='95% CI (lower)', pst2 format=5.3 label='95% CI (upper)', errorms format=5.3 label='sig^2', betwen format=5.3 label='siga^2', lambda format=5.3, ri format=5.3 label='ICC' from right; run; %end; %out: %mend rankcorr_mmer; %*rankcorr_mmer(data=example, ffq=repvitcn, dr=vitcno, repeated=4);