*********************************************
*         COCHRAN POST-HOC METHODS          *
*********************************************
* (c) Marta Garcia-Granero (07/2008)        *
* mgarciagranero@gmail.com                  *
* Downloaded from: http://gjyp.nl/marta/    *
* Feel free to use or modify this code, but *
* acknowledge the author and the web site   *
*********************************************

* Important: For the 2nd MACRO 'C:\Temp' folder must exist *.
PRESERVE.
SET RESULTS=NONE MESSAGES=NONE ERRORS=NONE.
HOST COMMAND=['IF not exist C:\Temp MD C:\Temp'].
RESTORE.

* Macro 1: asymptotic methods *.

* See HANDBOOK OF PARAMETRIC AND NONPARAMETRIC STATISTICAL PROCEDURES
  Test 26 (David J Sheskin, 2000; Chapman & Hall) *.

****************************************************************
* Cochran test & Multiple Comparisons for large samples        *
* Tukey-Kramer & SNK methods available for up to 20 groups,    *
* Dunnett test for 12 groups (1 control & 11 treatment groups) *
* Critical values for Dunnett test obtained from: Dunnett, CW  *
* (1964) "New tables for multiple comparisons with a control"  *
* Biometrics 20, 482-91                                        *
* Critical Values of the Studentized Range Statistic obtained  *
* from http://cse.niaes.affrc.go.jp/miwa/probcalc/             *
****************************************************************

DEFINE MULTICOC (!POSITIONAL !CMDEND).
PRESERVE.
SET MXLOOPS=1000. /* Increase it if sample size is GT 1000 *.
MATRIX.
PRINT /TITLE='COCHRAN TEST & MULTIPLE COMPARISON METHODS'.
GET data /VAR=!1 /NAME=gnames /MISSING=OMIT.
DO IF (MMIN(data) EQ 0) AND (MMAX(data) EQ 1).
COMPUTE b=NROW(data).
COMPUTE k=NCOL(data).
COMPUTE ci=CSUM(data).
COMPUTE fi=RSUM(data).
COMPUTE total=MSUM(data).
COMPUTE valuelab={'Negative','Positive'}.
PRINT {T(b-ci),T(ci)}
 /FORMAT='F8.0'
 /CNAMES=valuelab
 /RNAMES=gnames
 /TITLE='          Frequencies (*)'.
PRINT /TITLE='(*) 1 is treated as Positive'.
COMPUTE q=(k-1)*(k*MSUM(ci&**2)-total&**2)/(k*total-MSUM(fi&**2))).
COMPUTE qsig=1-CHICDF(q,(k-1)).
PRINT {k-1,b}
 /FORMAT='F8.0'
 /CLABEL='      DF','       N'
 /TITLE='Degrees of Freedom & Sample Size'.
PRINT {q,qsig}
 /FORMAT='F8.3'
 /CLABEL='       Q','    Sig.'
 /TITLE='Cochran statistic & asymptotic significance (*)'.
COMPUTE discord=0.
LOOP i=1 TO b.
- DO IF (fi(i) NE 0) AND (fi(i) NE k).
-  COMPUTE discord=discord+1.
- END IF.
END LOOP.
DO IF discord*k LT 24.
- PRINT /TITLE='(*) Not enough discordant patterns for asymptotic testing (use McNemar for MC)'.
ELSE.
- PRINT /TITLE='(*) Enough discordant patterns for asymptotic testing'.
END IF.
COMPUTE se_rankd=SQRT(2*(k*total-MSUM(fi&**2))/(k*(k-1))).
COMPUTE drankd=MAKE(k-1,1,0).
COMPUTE compard=MAKE(k-1,3,' ').
LOOP i=1 TO k-1. 
- COMPUTE drankd(i)=ABS(ci(1)-ci(i+1)).
- COMPUTE compard(i,:)={gnames(1),'  vs    ',gnames(i+1)}.
END LOOP.
COMPUTE dunnett={1.96,2.21,2.35,2.41,2.51,2.57,2.65,2.72,2.77,2.83,2.91;  
                 2.58,2.79,2.92,3.00,3.06,3.11,3.19,3.25,3.29,3.35,3.42}.
COMPUTE dunnett5=dunnett(1,k-1)*se_rankd.
COMPUTE dunnett1=dunnett(2,k-1)*se_rankd.
COMPUTE sigd=MAKE(k-1,1,' ').
LOOP i=1 TO k-1.
- DO IF drankd(i) GE dunnett1.
-  COMPUTE sigd(i)='   <0.01'.
- ELSE IF drankd(i) GE dunnett5.
-  COMPUTE sigd(i)='   <0.05'.
- ELSE.
-  COMPUTE sigd(i)='      NS'.
- END IF.
END LOOP.
PRINT {compard,sigd}
 /FORMAT='A8'
 /TITLE='Comparisons (adjusted): Dunnett method (1st group is reference)'.
RELEASE se_rankd,drankd,compard,dunnett,dunnett5,dunnett1,sigd.
COMPUTE sci=ci.
COMPUTE sci(GRADE(ci))=ci.
COMPUTE snames=gnames.
COMPUTE snames(GRADE(ci))=gnames.
RELEASE data,gnames,ci.
COMPUTE num=sci.
COMPUTE den=MAKE(1,k,b).
COMPUTE p=num/den .
COMPUTE z = 1.95996398.
COMPUTE zsq = z*z.
COMPUTE x1 = 2*num+zsq .
COMPUTE x2 = z*(zsq+4*num&*(1-p))&**0.5 .
COMPUTE x3 = 2*(den+zsq) .
COMPUTE x4 = (x1 - x2) / x3 .
COMPUTE x5 = (x1 + x2) / x3 .
PRINT {T(p),T(x4),T(x5)}
 /FORMAT='F8.3'
 /RNAME=snames
 /CLABELS='Point','Lower95%','Upper95%'
 /TITLE='Sorted Proportions (in ascending order) & 95% CI'.
RELEASE num,den,p,z,zsq,x1,x2,x3,x4,x5.
COMPUTE se_rank=SQRT(2*(k*total-MSUM(fi&**2))/(k*(k-1))).
COMPUTE drank=MAKE(k*(k-1)/2,1,0).
COMPUTE compar=MAKE(k*(k-1)/2,3,' ').
COMPUTE sig=MAKE(k*(k-1)/2,1,' ').
COMPUTE count=0.
LOOP i=1 TO k-1.
- LOOP j=i+1 TO k.
-  COMPUTE count=count+1.
-  COMPUTE drank(count)=ABS(sci(i)-sci(j)).
-  COMPUTE compar(count,:)={snames(i),'  vs    ',snames(j)}.
- END LOOP.
END LOOP.
COMPUTE pvalue=2*(1-CDFNORM(drank/se_rank)).
LOOP i=1 TO count.
- DO IF pvalue(i) LE 0.01 .
-  COMPUTE sig(i)='   <0.01'.
- ELSE IF pvalue(i) LE 0.05.
-  COMPUTE sig(i)='   <0.05'.
- ELSE.
-  COMPUTE sig(i)='      NS'.
- END IF.
END LOOP.
PRINT {compar,sig}
 /FORMAT='A8'
 /TITLE='Comparisons (unadjusted): LSD test'.
COMPUTE apvalue=1-(1-pvalue)&**count.
LOOP i=1 TO count.
- DO IF apvalue(i) LE 0.01.
-  COMPUTE sig(i)='   <0.01'.
- ELSE IF apvalue(i) LE 0.05.
-  COMPUTE sig(i)='   <0.05'.
- ELSE.
-  COMPUTE sig(i)='      NS'.
- END IF.
END LOOP.
PRINT {compar,sig}
 /FORMAT='A8'
 /TITLE='Comparisons (adjusted): Dunn-Sidak test'.
COMPUTE studrang={2.772,3.314,3.633,3.858,4.030,4.170,4.286,4.387,4.474,
                  4.552,4.622,4.685,4.743,4.796,4.845,4.891,4.934,4.974,5.012;  
                  3.643,4.120,4.403,4.603,4.757,4.882,4.987,5.078,5.157,
                  5.227,5.290,5.348,5.400,5.448,5.493,5.535,5.574,5.611,5.645}.
COMPUTE ctukey5=studrang(1,k-1)*se_rank/SQRT(2).
COMPUTE ctukey1=studrang(2,k-1)*se_rank/SQRT(2).
LOOP i=1 TO count.
- DO IF drank(i) GE ctukey1.
-  COMPUTE sig(i)='   <0.01'.
- ELSE IF drank(i) GE ctukey5.
-  COMPUTE sig(i)='   <0.05'.
- ELSE.
-  COMPUTE sig(i)='      NS'.
- END IF.
END LOOP.
PRINT {compar,sig}
 /FORMAT='A8'
 /TITLE='Comparisons (adjusted): Tukey-Kramer method'.
COMPUTE snk5=MAKE(k-1,1,0).
COMPUTE snk1=MAKE(k-1,1,0).
LOOP i=1 TO k-1.
- COMPUTE snk5(i)=studrang(1,k-i)*se_rank/SQRT(2).
- COMPUTE snk1(i)=studrang(2,k-i)*se_rank/SQRT(2).
END LOOP.
COMPUTE step=MAKE(k*(k-1)/2,1,0).
COMPUTE count=0.
LOOP i=1 TO k.
- LOOP j=i+1 TO k.
-  COMPUTE count=count+1.
-  COMPUTE step(count)=k+i-j.
- END LOOP.
END LOOP.
LOOP i=1 TO count.
- DO IF drank(i) GE snk1(step(i)).
-  COMPUTE sig(i)='   <0.01'.
- ELSE IF drank(i) GE snk5(step(i)).
-  COMPUTE sig(i)='   <0.05'.
- ELSE.
-  COMPUTE sig(i)='      NS'.
- END IF.
END LOOP.
PRINT {compar,sig}
 /FORMAT='A8'
 /TITLE='Comparisons (adjusted): Step-down SNK method'.
ELSE.
- PRINT /TITLE='Warning: Data are not 0/1 coded'.
END IF.
END MATRIX.
RESTORE.
!ENDDEFINE.


* Macro 2: pairwise McNemar tests with exact p-values *.

DEFINE MULTIMCN (!POSITIONAL !CMDEND).
DATASET NAME Data.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='McNemar Crosstabs'
 /DESTINATION VIEWER=NO.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='McNemar Test Statistics'
 /DESTINATION FORMAT=SAV
  OUTFILE='C:\Temp\PairwiseMcNemarTestPvalues.sav'.
NPAR TEST
   /MCNEMAR = !1.
OMSEND.
GET FILE='C:\Temp\PairwiseMcNemarTestPvalues.sav' /DROP=Command_ TO Label_.
DATASET NAME PValues.
SELECT IF Var1 NE 'N'.
PRESERVE.
SET ERRORS=NONE RESULTS=NONE.
FLIP VAR=ALL/NEWNAME=Var1  .
DELETE VARIABLES CASE_LBL.
RENAME VARIABLE (ALL=pvalue).
RANK pvalue /n into N.
RESTORE.
COMPUTE id = $CASENUM.
FORMAT id (F2.0).
SORT CASES BY  pvalue (A) .
COMPUTE pos = $CASENUM.
FORMAT pos (F2.0).
COMPUTE bonferr=MIN(pvalue*n,1).
COMPUTE sidak=1-(1-pvalue)**n.
COMPUTE holm = MIN(1,(n-pos+1)*pvalue).
IF (holm LT LAG(holm)) holm = LAG(holm).
COMPUTE downsidk = 1-(1-pvalue)**(n-pos+1).
IF (downsidk LT LAG(downsidk)) downsidk = LAG(downsidk).
COMPUTE finner = 1-(1-pvalue)**(n/pos).
IF (finner LT LAG(finner)) finner = LAG(finner).
COMPUTE cn = cn+1/pos.
LEAVE cn.
SORT CASES BY pos(D).
IF cn LT LAG(cn) cn = LAG(cn).
COMPUTE hommel = MIN(1,cn*n*pvalue/pos).
IF (hommel GT LAG(hommel)) hommel = LAG(hommel).
COMPUTE hochberg = (n-pos+1)*pvalue.
IF (hochberg GT LAG(hochberg)) hochberg = LAG(hochberg).
COMPUTE simes = n*pvalue/pos.
IF (simes GT LAG(simes)) simes = LAG(simes).
EXECUTE.
DELETE VARIABLES pos,n,cn.
FORMAT pvalue bonferr to simes (F9.4).
VARIABLE LABELS id 'Nr.' /pvalue 'Original p-value'
        /bonferr 'One-step Bonferroni' /sidak 'One-step Sidak'
        /holm 'Step-down Holm' /downsidk 'Step-down Dunn-Sidak'
        /finner 'Step-down Finner' /hommel 'Step-up Hommel'
        /hochberg 'Step-up Hochberg' /simes 'Step-up Simes'.
SORT CASES BY id (A).
OMS /SELECT TABLES
 /IF COMMANDS='Summarize' SUBTYPES='Case Processing Summary'
 /DESTINATION VIEWER=NO.
 SUMMARIZE
  /TABLES = pvalue bonferr TO simes
  /FORMAT = LIST NOCASENUM TOTAL
  /TITLE = 'Original and adjusted p-values'
  /MISSING = VARIABLE
  /CELLS = NONE.
OMSEND.
DATASET ACTIVATE Data.
DATASET CLOSE PValues.
!ENDDEFINE.

* Sample dataset (replace by your own): data MUST be 0/1 coded *.
DATA LIST FREE /march june august (3 F8.0).
BEGIN DATA
0 1 0 0 0 0 0 1 0 1 0 0 1 1 0 1 0 0
0 0 0 1 1 0 1 1 0 1 0 0 1 0 0 1 1 0
0 1 1 0 1 1 0 0 0 0 1 0 0 1 0 1 0 0
1 0 0 1 0 0 1 0 1 1 1 1 1 0 0 1 1 0
END DATA.
VALUE LABEL march TO august 0'Negative' 1'Positive'.

DOCUMENT'Janzen, D. H. 1967. Interaction of the bull's-horn acacia (Acacia
         cornigera L.) with an ant inhabitant (Pseudomyrmex ferruginea F. Smith)
         in eastern Mexico. University of Kansas Science Bulletin 47:315-558.'.

* MACRO calls *.
MULTICOC march TO august.
MULTIMCN march TO august.