libname t 'c:\tas'; options ps=85 ls=100 mprint; proc printto log='c:\pval' new; run; ods listing close; ods results off; ****************************************************************************** The SIMULATIONS program below will produce 50 observed and adjusted p-values in about 8 minutes on an on an H-P workstation with dual processors (3.20 GHz) and 1 Gig Ram. 10 observed and 10 adjusted p-values take about 2 minutes, 100 in about 16 minutes, and 1000, as in paper, about 2 1/2 hours. Change 2nd line from top and 7th line from bottom to suite your preference. ******************************************************************************; /* Author: Eugene Komaroff */ ********* BEGIN SIMULATIONS **************************************************; %MACRO DOVER; %DO MMM=1 %TO 50; /* change to 1000 for closer approximation to results in paper */ %let seed=%eval(&mmm-1)+123456*&mmm; data c1 c2; do sub=1 to 100; true=rannor(&seed); cor=0.70; v1=sqrt(cor)*true+sqrt(1-cor)*rannor(&seed); gp=1; output c1; v2=sqrt(cor)*true+sqrt(1-cor)*rannor(&seed); gp=2; output c2; end; run; data c3; merge c1 c2; by sub; if v1 > 0; dif = v2-v1; run; proc ttest data=c3; var dif; ods output ttests=t; run; data t2; set t; rename probt = pobs df=dfobs; do i = 1 to 100; output; end; run; proc means data = c3 mean std var; var v2 v1; output out=m (where = (_STAT_ = 'MEAN' or _STAT_ = 'STD')); run; %macro DOIT (seed1=, seed2=, num=); data n; set m; retain ave; if _STAT_ = 'MEAN' then do; ave = 1/2 * (v1 + v2); end; if _STAT_ = 'STD' then do; var = 1/2 * (v1**2 + v2**2); end; std = sqrt(var); keep ave std; if _n_ = 2; run; data n2; set n; do sub=1 to 100; s1=ave+std*rannor(&seed1); s2=ave+std*rannor(&seed2); output; end; keep s1 s2; run; data n3; set n2; if s1 > 0; dif=s2-s1; run; proc ttest data=n3; var dif; ods output TTests = p# ods output Statistics = m# run; data final; set p&num %if &num > 1 %then %do; final %end;; keep df tValue probt; run; %MEND DOIT; %macro runit(r=); %do j = 1 %to &r; %let m=%eval(&j-1)+1; %doit (seed1 = &m, seed2=&m+91, num=&j); %end; %mend; %runit(r=100) data padj_&MMM; merge t2 final; if probt <= pobs then flag=1; else flag=0; dataset=&MMM; run; %end; %MEND DOVER; %DOVER; ****Appending simulated data sets**************; %MACRO THEEND; data padj; set %do m=1 %to 50; /* change to 1000 for closer approximation to results in paper */ padj_&m %end; ; run; %MEND; %THEEND; *******END SIMULATIONS**************************************************************************; *******BEGIN ANALYSIS***************************************************************************; proc printto log=log; run; proc freq data=padj; table dataset * flag / out=flag sparse; /**Cumulate: probt (p_i) <= pobs then flag=1; else flag=0**/ run; data flag2; set flag; if flag=1 then padj = percent; run; data flag3; set flag2; by dataset; if last.dataset; keep dataset padj; run; data pobs; set padj; by dataset; if first.dataset; keep dataset pobs; run; proc sort data=flag3; by dataset; run; proc sort data=pobs; by dataset; run; data p; merge flag3 pobs; by dataset; run; ods listing; options ps=50; proc chart data=p; vbar pobs; run; proc chart data=p; vbar padj; run; proc plot data=p; plot padj * pobs; run;