/* Interactive SAS program for calculating critical values of multiple comparison procedures. For details see the order of events below and the short description of the cooresponding QBATCH program. Author : Frank Bretz and Paul Somerville Contact: bretz@ifgb.uni-hannover.de Date : 11.05.99 - 1.0 - start version - 1.1 - skiped 02.06.00 - 1.2 - inclusion of several new swtypes (7-15) References: Somerville (1997, 1998, 1999), Somerville and Bretz (2000, 2001) */ /* Order of events: info1 swtype info2 conf info3 k info4 ndenom info5 irep if(irep.lt.0)then info6 givense else info7 mocar if(mocar.lt.0)then info8 irep info9 mocar info10 swcov if(swcov.eq.0)then info11 vc(i,j) if((swcov.eq.0).and.((swtype.eq.4).or.(swtype.eq.5)))then info12 icontrol if(swtype.eq.0)then info13 mm info14 contr if(swtype.eq.10)then info16 iflo ighi if (swtype=11) or (swtype=12) or (swtype=13) or (swtype=14) or (swtype=15) then info17 znn info15 seed */ %let givense=0; %let swcov=1; %let irep=0; %let mocar=-1; %let mm=0; %let icontrol=.; %let iflo=.; %let ighi=.; %window info1 #10 @10 "This program calculates the constants necessary for" #11 @10 "obtaining simultaneous confidence intervals.For two-" #12 @10 "sided confidence intervals the form is" #13 @10 "(x(i)-x(j)) plus or minus sd(x(i)-x(j))*q/sqrt(2)" #15 @10 "The program has the capability of calculating q" #16 @10 "(or several q s for Hsu s MCB) for the following cases" #17 @10 " 1 Tukey" #18 @10 " 2 Bofinger" #19 @10 " 3 Hsu s MCB" #20 @10 " 4 Dunnett s comparisons with a control (1 sided)" #21 @10 " 5 Dunnett s comparisons with a control (2 sided)" #22 @10 " 6 OSRT" #23 @10 " 7 Successive ordered treatments (one-sided)" #24 @10 " 8 Successive ordered treatments (two-sided)" #25 @10 " 9 umbrella contrasts" #26 @10 " 10 generalized umbrella contrasts" #27 @10 " 11 Williams type contrast" #28 @10 " 12 Marcus type contrast" #29 @10 " 13 Hirotsu test" #30 @10 " 14 Test of McDermott and Mudholkar" #31 @10 " 15 Isootonic contrast" #32 @10 " 0 Arbitrary set of contrasts (1 and 2 sided)" #37 @10 "Enter your choice and press ENTER" +2 swtype 2 attr=underline; %window info2 #10 @10 "Enter desired confidence as a decimal (i.e .95) and" #11 @10 "press ENTER." +2 conf 5 attr=underline; %window info3 #10 @10 "Enter no. of populations or treatments, press enter." +2 k 2 attr=underline; %window info4 #10 @10 "Give df for error estimate (if the variance is assumed" #11 @10 "known type -1) and press enter" +2 ndenom 4 attr=underline; %window info5 #10 @10 "The calculation of q involves Monte Carlo methods." #11 @10 "You may select the number of random directions to use." #12 @10 "10000 is frequently sufficient, or you may specify the" #13 @10 "standard error of the calculated estimate for q" #15 @10 "To specify the standard error type -1 and press enter" #16 @10 "For no. of random directions, type 0, press enter." #18 @10 "Enter your choice and press ENTER" +2 irep 3 attr=underline; %window info6 #10 @10 "The smaller the standard error, the longer the running" #11 @10 "time of the program will be. Some choices are .01," #12 @10 ".005, or .001. Type choice, press ENTER" +2 givense 7 attr=underline; %window info7 #10 @10 "Type number of random directions you wish to use and" #11 @10 "press enter. (The standard error of your estimate will" #12 @10 "be calculated using an empirical formula.)" #14 @10 "An alternate procedure is to make n independent esti-" #15 @10 "mates of q, each with nn random directions. The stand-" #16 @10 "ard error of estimated q is calculated using the n" #17 @10 "estimates. The alternate procedure is especially use-" #18 @10 "ful when the estimates of treatment means are" #19 @10 "correlated and an accurate formula for the standard" #20 @10 "error of the estimate of q is not available. " #22 @10 "If the alternate procedure is desired, type -1 and" #23 @10 "press enter" +2 mocar 7 attr=underline; %window info8 #10 @10 "Give n and press enter." +2 irep 7 attr=underline; %window info9 #10 @10 "Give nn and press enter." +2 mocar 7 attr=underline; %window info10 #10 @10 "The program will calculate the constant q for an" #11 @10 "arbitrary variance-covariance matrix (correlated case)" #12 @10 "or the uncorrelated case." #14 @10 " 0 correlated or unequal var case" #15 @10 " 1 uncorrelated equal variance case" #16 @10 "Enter 0 or 1 and press enter" +2 swcov 2 attr=underline; %window info11 #10 @10 "Input the lower triangular portion of the variance covariance matrix." #12 @10 "Press enter after each value" +2 vc_i 7 attr=underline; %macro mvar; %display info10; %if (&swcov=0) %then %do; %do i=1 %to &k; %do j=1 %to &i; %display info11; vc[&i,&j]=&vc_i; %end; %end; %end; %else %put 'hallo'; %mend; %window info12 #10 @10 "Type number of the control population and press enter." +2 icontrol 2 attr=underline; %window info13 #10 @10 "A confidence interval of the form (g,+infinity) for" #11 @10 "-x2+x4 is input as 0 -1 0 +1 (one line)" #12 @10 "A confidence interval of the form (-infinity,g) for" #14 @10 "-x2+x4 is input as 0 +1 0 -1 (one line)" #14 @10 "A confidence interval of the form (-g,g) for -x2+x4" #15 @10 "requires both (two lines)" #16 @10 "Type no. of lines required for contrasts," #27 @10 "press enter" +2 mm 2 attr=underline; %window info14 #10 @10 "Type the contrasts, pressing enter after each value" +2 contr_i 7 attr=underline; %window info16a #10 @10 "Give iflo and press enter." +2 iflo 2 attr=underline; %window info16b #10 @10 "Give ighi and press enter." +2 ighi 2 attr=underline; %window info17 #10 @10 "Type the sample sizes, pressing enter after each value" +2 znn_i 3 attr=underline; %macro mcontr; %if (&swtype=0) %then %do; %do i=1 %to &mm; %do j=1 %to &k; %display info14; contr[&i,&j]=&contr_i; %end; %end; %end; %else %if (&swtype=9) %then %do; %display info16a; %display info16b; %end; %else %if (&swtype=11) or (&swtype=12) or (&swtype=13) or (&swtype=14) or (&swtype=15) %then %do; %do j=1 %to &k; %display info17; znn[&j]=&znn_i; %end; %end; %mend; %window info15 #10 @10 "Type seed for random no generator, press enter" +2 seed 10 attr=underline; %macro main; %display info1; %display info2; %display info3; %display info4; %display info5; %if (&irep<0) %then %display info6; %else %do; %display info7; %if (&mocar<0) %then %do; %display info8; %display info9; %end; %end; %if (&swtype=0) %then %display info13; %display info15; %mend; %macro icontr; %if (&swcov=0) %then %do; %if ((&swtype=4) | (&swtype=5)) %then %display info12; %end; %if (&swtype=9) %then %display info12; %mend; %main; proc iml symsize=250000; swtype = &swtype; conf = &conf; k = &k; ndenom = &ndenom; irep = &irep; mocar = &mocar; givense = &givense; vc=j(k,k,.); %mvar; %icontr; swcov = &swcov; icontrol = &icontrol; do i=1 to k-1; do j=i+1 to k; vc[i,j]=vc[j,i]; end; end; if (swtype=0) then do; mm = &mm; contr=j(mm,k,.); end; znn=j(1,k,.); %mcontr; ighi = &ighi; iflo = &iflo; seed = &seed; /* swcov = 0 ; swtype = 3 ; conf = .9 ; givense= 999 ; k = 4 ; seed = 8763 ; ndenom = 25 ; mocar = 1000 ; irep = 100 ; mm = 99 ; */ nq = 25 ; /* Optional Parameters vc={76.807 36.148 -20.366 7.2384, 36.148 27.373 15.176 21.198, -20.366 15.176 64.583 40.605, 7.2384 21.198 40.605 31.239}; contr=.; icontrol=.; */ /*********************************************************** *********************************************************** ***********************************************************/ /* This function gives the logarithm of the gamma function when the argument is a positive integer or an integer divided by 2 */ start gln(a); i1=int(a); c2=0; c1=1; gln=0; if (a-i1)=0 then i1=i1-2; else do; c1=-.5; gln=log(sqrt(arcos(-1))); end; do i=1 to i1; z=c1+i; x=log(z); gln=gln+x; end; return(gln); finish; /* This function gives the ratio of gamma(a+b) to gamma(a)*gamma(b) */ start beta(a,b); beta=gln(a+b)-gln(a)-gln(b); beta= exp(beta); return(beta); finish; /* Calculates k-1 rows of Helmert matrix */ start helmert(k,k1); hm=j(k1,k,.); do i=1 to k1; zi=i; do j=k-i to k; hm[i,j]=1/sqrt(zi*(zi+1)); end; do j=1 to k1-i; hm[i,j]=0; end; hm[i,k-i]=-zi*hm[i,k-i]; end; return(hm); finish; /* SUBROUTINE TO ENUMERATE ALL OF THE PAIRS OF FACES */ start facelist(mm) global(k, ifaceb,ifacee); ifacee=j(1,mm,.); ifaceb=j(1,mm,.); i=0; ie=2; do until(ie=k+1); ib=1; krit=0; do until(krit=1); i=i+1; ifaceb[i]=ib; ifacee[i]=ie; if (ib^=ie-1) then ib=ib+1; else krit=1; end; ie=ie+1; end; finish; /* Generation of Contrast Matrices */ start contrast(swtype,mm,icontrol) global(znn,ighi,iflo,k,k1,ifaceb,ifacee); contr=j(mm,k,0); if (swtype=1 | swtype=2 | swtype=6) then run facelist(mm); if swtype=1 then do; imm=mm/2; do i=1 to imm; j=i+imm; contr[i,ifaceb[i]]=1; contr[i,ifacee[i]]=-1; contr[j,ifaceb[i]]=-1; contr[j,ifacee[i]]=1; end; end; if (swtype=2 | swtype=6) then do; do i=1 to mm; contr[i,ifaceb[i]]=1; contr[i,ifacee[i]]=-1; end; end; if (swtype=3 | swtype=4) then do; do i=1 to mm; contr[i,icontrol]=-1; end; do i=1 to icontrol-1; contr[i,i]=1; end; do i=icontrol+1 to mm; contr[i,i+1]=1; end; if icontrol<=mm then contr[icontrol,icontrol+1]=1; end; if swtype=5 then do; do i=1 to k1; j=i+k1; contr[i,icontrol]=-1; contr[j,icontrol]=1; end; do i=1 to icontrol-1; j=i+k1; contr[i,i]=1; contr[j,i]=-1; end; do i=icontrol+1 to k1; j=i+k1; contr[i,i+1]=1; contr[j,i+1]=-1; end; if icontrol<=k1 then do; contr[icontrol,icontrol+1]=1; contr[icontrol+k1,icontrol+1]=-1; end; end; if (swtype=7 | swtype=8) then do; do i=1 to k-1; contr[i,i]=1; contr[i,i+1]=-1; end; end; if (swtype=8) then do; do i=k to mm; contr[i,i-k+1]=-1; contr[i,i-k+2]=1; end; end; if (swtype=9) then do; mm1=icontrol*(icontrol-1)/2; mm2=(k+1-icontrol)*(k-icontrol)/2; if (mm1 ^= 0) then do; ifacee=j(1,mm,0); ifaceb=j(1,mm,0); i=0; ie=2; do until(ie=icontrol+1); ib=1; do until(ib=ie); i=i+1; ifaceb[i]=ib; ifacee[i]=ie; ib=ib+1; end; ie=ie+1; end; do jj=1 to mm1; contr[jj,ifaceb[jj]]=-1; contr[jj,ifacee[jj]]=1; end; end; if (mm2 ^= 0) then do; ifacee=j(1,mm,0); ifaceb=j(1,mm,0); i=0; ie=2; do until(ie=k-icontrol+2); ib=1; do until(ib=ie); i=i+1; ifaceb[i]=ib; ifacee[i]=ie; ib=ib+1; end; ie=ie+1; end; do jj=1 to mm2; contr[jj+mm1,ifaceb[jj]+icontrol-1]=1; contr[jj+mm1,ifacee[jj]+icontrol-1]=-1; end; end; end; if (swtype=10) then do; in=0; do i=1 to ighi-1; do j=i+1 to ighi; in=in+1; contr[in,i]=1; contr[in,j]=-1; end; end; do j=iflo to k-1; do i=j+1 to k; in=in+1; contr[in,j]=-1; contr[in,i]=1; end; end; end; if swtype=11 then contr=cm_will(znn); if swtype=12 then contr=cm_marc(znn); if swtype=13 then contr=cm_hirot(znn); if swtype=14 then contr=cm_mcder(znn); if swtype=15 then contr=cm_iso(znn); print contr; return(contr); finish; start cm_will(n); k1=ncol(n)-1; c=j(1,k1,.); cm=j(k1,k1,0); cm=j(k1,1,-1)||cm; do i=1 to k1; x=sum(n[k1-i+2:k1+1]); do j=1 to i; cm[i,k1-j+2]=n[k1-j+2]/x; end; end; return(cm); finish; start cm_marc(n); k=ncol(n); cm1=j(k-1,k,0); cm2=cm1; do i=1 to k-1; cm1[i,i+1:k]=t(n[i+1:k]/sum(n[i+1:k])); end; do i=1 to k-1; cm2[i,1:i]=t(n[1:i]/sum(n[1:i])); end; row=k*(k-1)/2; cm=j(row,k,0); index=1; do i=1 to k-1; do j=1 to i; cm[index,]=cm1[i,]-cm2[j,]; index=index+1; end; end; return(cm); finish; start cm_hirot(n); k=ncol(n); cm=j(k-1,k,0); do i=1 to k-1; do j=1 to k; if j>i then cm[i,j]=n[j]/sum(n[i+1:k]); else cm[i,j]=-n[j]/sum(n[1:i]); end; end; return(cm); finish; start cm_mcder(n); k=ncol(n); cm=j(k-1,k,0); do i=1 to k-1; do j=1 to k; if j <= i then cm[i,j]=-n[j]/sum(n[1:i]); else if j=i+1 then cm[i,j]=1; end; end; return(cm); finish; /******************************************************************/ /* The module SCHNITT cuts the place STELLE out of the vector VEC */ /* and returns a vector which dimension dimension is reduced by 1 */ /******************************************************************/ start schnitt(vec,stelle); col=ncol(vec); if stelle=col then neu=vec[,1:col-1]; else do; neu=vec[,1:col-1]; neu[stelle:col-1]=vec[,stelle+1:col]; end; return(neu); finish; /******************************************************************/ /* The main module to calculate the contrast coefficients */ /******************************************************************/ start cm_iso(n); k=ncol(n)-1; anz=2##k; cm=j(2##k-1,k+1,.); mat=j(anz,k,0); mat[2,k]=1; /*****************************************/ /* Generation of all 2^k-1 possibilities */ /* 0: two subsequent groups are pooled */ /* 1: subsequent groups are not pooled */ /*****************************************/ do i=2 to k; mat[2##(i-1)+1:2##i,k-i+1]=j(2##(i-1),1,1); mat[2##(i-1)+1:2##i,k-i+2:k]=mat[1:2##(i-1),k-i+2:k]; end; mat=mat[2:2##k,]; row=nrow(mat); n1=j(1,k,n[1]); n2=j(1,k,sum(n[2:k+1])); do j=2 to k; n1[j]=sum(n[1:j]); n2[j]=sum(n[j+1:k+1]); end; /*************************************************/ /* Beginning of main loop: each row = one of the */ /* 2^k-1 possibilities is evaluated on its own */ /*************************************************/ do i=1 to row; y=mat[i,]; /* Amalgamated contrast vector */ x=j(1,k,1); /* Weight vector according to the number of amalgamations */ s=mat[i,+]; ssd=j(1,s+1,0); nn=-y#sqrt(n1#n2/n[+]); nn_help=nn; do j=1 to ncol(nn); nn[j]=nn_help[ncol(nn)-j+1]; end; krit=0; j=1; do until(krit=1); if nn[j]=0 then nn=schnitt(nn,j); else j=j+1; if j>ncol(nn) then krit=1; end; ssd[2:s+1]=nn; ssd_help=ssd; do j=1 to ncol(ssd); ssd[j]=ssd_help[ncol(ssd)-j+1]; end; krit=0; j=1; y_help=y; n_gew=n; do until(krit=1); if y_help[j]=0 then do; n_gew[j]=n_gew[j]+n_gew[j+1]; n_gew=schnitt(n_gew,j+1); y_help=schnitt(y_help,j); end; else j=j+1; if all(y_help)=1 then krit=1; end; krit=0; j=1; do until(krit=1); /* Actual amalgamation */ if y[j]=0 then do; if y[j]=y[j+1] then do; y=schnitt(y,j+1); /* Pooling, i.e. reduction of contrast vector */ x=schnitt(x,j+1); /* Reduction des weight vector */ x[j]=x[j]+1; /* Increase of weight by 1, if pooling happend */ end; else j=j+1; end; else j=j+1; if j>=ncol(y) then krit=1; end; do j=1 to ncol(y); /* Correction of weights */ if y[j]=0 then x[j]=x[j]+1; end; j=1; do while(j<(ncol(y)-1)); /* Special case. if a non-pooling group lies between */ if y[j]=0 then do; /* two groups to be pooled, e.g. 010 */ l=j+1; do until(y[l]=0 | l>=ncol(y)); l=l+1; end; if y[l]=0 then do; y=schnitt(y,j+1); x=schnitt(x,j+1); end; end; j=j+1; end; if i=row then do; /* Special case, if no pooling at all */ y=y||{1}; /* (always the last of the 2^k-1 cases) */ x=x||{1}; end; c=j(1,s+1,ssd[1]/n_gew[1]); do j=2 to s+1; c[j]=(ssd[j]-sum(c[1:j-1]#n_gew[1:j-1]))/n_gew[j]; end; contrast=c; cm[i,1:x[1]]=contrast[1]; /* Solution of amalgamation */ do j=2 to ncol(y); cm[i,sum(x[1:j-1])+1:sum(x[1:j])]=contrast[j]; end; end; cm=cm#n; return(cm); finish; /*********************************************************** *********************************************************** ***********************************************************/ /* This routine does the binning of the reciprocal distances */ start fillbin2(recdist,nr) global(pl, ibin); i=nr; krit=0; do until(krit=1 | i=0); if recdist>= pl[i] then do; j=nr+1-i; ibin[j]=ibin[j]+1; krit=1; end; else i=i-1; end; finish; /* This routine calculates proportions in the bins */ start binprop(points,nq) global(ibin, prop); do i=2 to nq; ibin[i]=ibin[i-1]+ibin[i]; end; do i=1 to nq; prop[i]=1-ibin[i]/points; end; finish; /* This routine calculates the probability of correct statements outside of region of radius q. */ start fringe(nr,prop) global(wl,pl,zm,zn,fconst); x=0; do i=1 to nr; j=nr-i+1; cc=zm*log(1+zn*pl[i]*pl[i])/2; cd=zn*log(1+1/zn/pl[i]/pl[i])/2; xx=exp(-log(pl[i])-cc-cd); x=x+xx*wl[i]*prop[j]; end; fringe=x*fconst; *print wl pl fconst; return(fringe); finish; /* Function obtains confidence corresponding to values of q, k1 and ndenom */ start prob(q,nq,prop) global(zm,zn); qq=sqrt(2)/q; run gauleg(0,qq,nq); if(zn>0) then do; x=fringe(nq,prop); qqdk1=q*q/zm/2; poff=probf(qqdk1,zm,zn); end; else do; x=fringinf(nq,zm,prop); qqdk1=q*q/2; poff=probchi(qqdk1,zm); end; prob2=x+poff; return(prob2); finish; /* Calculation of g(v) for Multivariate Normal Distribution */ start fringinf(nr,q,prop) global(wl,pl); x=0; do i=1 to nr; j=nr-i+1; ce=-log(gamma(q/2))-(q/2-1)*log(2)-(q+1)*log(pl[i])-1/pl[i]/pl[i]/2; xx=exp(ce); x=x+xx*wl[i]*prop[j]; end; fringinf=x; return(fringinf); finish; /* Determination of Nodes and Weigths for Gauss-Legendre */ start gauleg(x1,x2,n) global(pl,wl); x=j(1,n,.); w=j(1,n,.); eps2=1E-09; m=(n+1)/2; xm=(x2+x1)/2; xl=(x2-x1)/2; pi=2*arsin(1); do i=1 to m; z=cos(pi*(i-0.25)/(n+0.5)); do until (abs(z-z1)0 then do; tt=coc; mocar=int(max(add,256)); coc=qvalue(nq,mocar); coc=(10000*tt+mocar*coc)/(10000+mocar); mocar=mocar+10000; end; zmean = coc; se=givense; if add<0 then se=exp(beta[i+1,1]+beta[i+1,2]*coc+beta[i+1,2]/zk)*10; end; if j=0 then do; coc=qvalue(nq,mocar); zmean=sx/zrep + coc; zk=k; se=exp(beta[i+1,1]+beta[i+1,2]*coc+beta[i+1,2]/zk)*1000/sqrt(zmocar); end; if j=1 then do; mocar=zmocar; sx=0; ssx=0; do njj=1 to irep; qqq=qvalue(nq,mocar); *print qqq; qans[njj]=qqq; coc=qans[1]; cm=qans[njj]-coc; sx=sx+cm; ssx=ssx+cm*cm; end; sd=(ssx-sx*sx/zrep)/(zrep-1); se=sqrt(sd/zrep); sd=sqrt(sd); zmean=sx/zrep+coc; mocar=zrep*mocar; end; * print zmean; finish; /******************************************************************** ******************************************************************** ******************************************************************** ********************************************************************/ zmocar=mocar; if (swtype^=4 & swtype^=5 & swtype^=9) then icontrol=1; if (swtype=1) then mm=k*(k-1); if (swtype=2 | swtype=6) then mm=k*(k-1)/2; if (swtype=3 | swtype=4) then mm=k-1; if (swtype=5) then mm=2*(k-1); if (swtype=7) then mm=k-1; if (swtype=8) then mm=2*(k-1); if (swtype=9) then do; mm1=icontrol*(icontrol-1)/2; mm2=(k+1-icontrol)*(k-icontrol)/2; mm=mm1+mm2; end; if (swtype=10) then mm=ighi*(ighi-1)/2+(k-iflo+1)*(k-iflo)/2; if (swtype=11) then mm=k-1; if (swtype=12) then mm=k*(k-1)/2; if (swtype=13) then mm=k-1; if (swtype if (swtype=15) then mm=2**(k-1)-1; if (swtype^=0) then contr=j(mm,k,.); if (swtype^=10) then do; iflo=.; ighi=.; end; brakl = .0001; braku = 5; xacc = .0001; zero=0; zone=1; two=2; iseed=abs(seed); jseed=iseed; /* Random Number Generation */ iy=0; iran=j(1,97,.); mr=714025; ia=1366; ic=150889; rm=1/mr; do j=1 to 97; iseed=mod(ia*iseed+ic,mr); iran[j]=iseed; end; irr=iseed; k1=k-1; zm=k1; zn=ndenom; zmd2=zm/2; zmnd2=(zm+zn)/2; zmdzn=zm/zn; znd2=zn/two; if ndenom>0 then cnlogn=zn*log(zn)/2; sx=0; ssx=0; zrep=irep; dc=j(1,k1,.); iii=10*((irep-1)/10)+10; qans=j(1,iii,0); lines=(irep-1)/10+1; if ndenom>0 then fconst=2*beta(zmd2,znd2); if swcov^=0 then vc=I(k); hm=helmert(k,k1); vcc=hm*vc*t(hm); c=t(root(vcc))+1E-12; p=t(hm)*c; pp=p*t(p); wl=j(1,nq,.); pl=j(1,nq,.); prop=j(1,nq,.); ibin=j(1,nq,.); if (icontrol^=k & swtype=3 & swcov=0) then do; do until(icontrol>k); if ((swtype=4 | swtype=5) & swcov=1) then icontrol=1; if swtype^=0 then contr=contrast(swtype,mm,icontrol); gg=contr*p; dt=gg*t(gg); do i=1 to mm; dt[i,i]=sqrt(dt[i,i]); end; do i=1 to mm; do j=1 to k1; gg[i,j]=gg[i,j]/dt[i,i]; end; end; reset noname; if (swtype=3 & swcov=0 & icontrol^=1) then do; run estimate(swtype,nq,jjj); print 'q-value corresponding to population' icontrol ' is ' zmean; print mocar ' random directions used. ' 'swtype=' swtype; print 'Standard error of estimate is ' se; end; else do; if zrep<1 then jjj=-1; if zrep=1 then jjj=0; if zrep>1 then jjj=1; run estimate(swtype,nq,jjj); if (zrep<1 & swtype=zero) then do; print 'The estimated standard error may be very conservative'; print 'unless the sample means are highly correlated.'; end; print 'The q-value is ' zmean '. The seed is' jseed; print k ' populations,' conf ' confidence.'; if ndenom>0 then print 'df for variance estimate is ' ndenom; else print 'variance is assumed known'; print mocar ' random directions used. ' 'swtype=' swtype; print 'Standard error of estimate is ' se; end; icontrol=icontrol + 1; end; end; else do; if ((swtype=4 | swtype=5) & swcov=1) then icontrol=1; if swtype^=0 then contr=contrast(swtype,mm,icontrol); gg=contr*p; dt=gg*t(gg); do i=1 to mm; dt[i,i]=sqrt(dt[i,i]); end; do i=1 to mm; do j=1 to k1; gg[i,j]=gg[i,j]/dt[i,i]; end; end; reset noname; if (swtype=3 & swcov=0 & icontrol^=1) then do; run estimate(swtype,nq,jjj); print 'q-value corresponding to population' icontrol ' is ' zmean; print mocar ' random directions used. ' 'swtype=' swtype; print 'Standard error of estimate is ' se; end; else do; if zrep<1 then jjj=-1; if zrep=1 then jjj=0; if zrep>1 then jjj=1; run estimate(swtype,nq,jjj); if (zrep<1 & swtype=zero) then do; print 'The estimated standard error may be very conservative'; print 'unless the sample means are highly correlated.'; end; print 'The q-value is ' zmean '. The seed is' jseed; print k ' populations,' conf ' confidence.'; if ndenom>0 then print 'df for variance estimate is ' ndenom; else print 'variance is assumed known'; print mocar ' random directions used. ' 'swtype=' swtype; print 'Standard error of estimate is ' se; end; end; quit;