*********************************************
*        FRIEDMAN 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&3rd MACROS '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 25 (David J Sheskin, 2000; Chapman & Hall) *.

****************************************************************
* Friedman 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 MULTIFR (!POSITIONAL !CMDEND).
OMS /SELECT TABLES
 /IF SUBTYPE='Notes'
 /DESTINATION VIEWER=NO.
MATRIX.
PRINT /TITLE='FRIEDMAN TEST & MULTIPLE COMPARISON METHODS'.
GET rawdata /VAR=!1 /NAME=gnames /MISSING=OMIT.
COMPUTE b=NROW(rawdata).
COMPUTE k=NCOL(rawdata).
COMPUTE rdata=MAKE(b,k,0).
PRINT {k-1,b}
 /FORMAT='F8.0'
 /CLABELS='DF','N'
 /TITLE='Degrees of Freedom & Sample size'.
LOOP i=1 TO b.
- COMPUTE rdata(i,:)=RNKORDER(rawdata(i,:)).
END LOOP.
RELEASE rawdata.
COMPUTE ranks=CSUM(rdata).
COMPUTE VarRi = (MSSQ(rdata)-(MSUM(rdata))&**2/(b*k))/(b*(k-1)).
COMPUTE totalchi=(MSSQ(ranks)/b-(MSUM(rdata))&**2/(b*k))/VarRi.
COMPUTE tchisig=1-CHICDF(totalchi,k-1).
PRINT {totalchi,tchisig}
 /FORMAT='F8.4'
 /CLABELS='Chi^2','Sig.'
 /TITLE='Friedman test (asymptotic)'.
COMPUTE se_rankd=SQRT(VarRi*2*b).
COMPUTE drankd=MAKE(k-1,1,0).
COMPUTE compard=MAKE(k-1,3,' ').
LOOP i=1 TO k-1. 
- COMPUTE drankd(i)=ABS(ranks(1)-ranks(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 sranks=ranks.
COMPUTE sranks(GRADE(ranks))=ranks.
COMPUTE snames=gnames.
COMPUTE snames(GRADE(ranks))=gnames.
PRINT {sranks/b}
 /FORMAT='F8.2'
 /CNAME=snames
 /TITLE='Sorted Mean Ranks (in ascending order)'.
RELEASE rdata,gnames.
COMPUTE se_rank=SQRT(VarRi*2*b).
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(sranks(i)-sranks(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'.
END MATRIX.
OMSEND.
!ENDDEFINE.


* Macro 2 (Pairwise sign tests - with exact p-value) *.

DEFINE MULTISG (!POSITIONAL !CMDEND).
DATASET NAME Data.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='Sign Frequencies'
 /DESTINATION VIEWER=NO.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='Sign Test Statistics'
 /DESTINATION FORMAT=SAV
 OUTFILE='C:\Temp\PairwiseSignTestPvalues.sav'.
NPAR TEST
 /SIGN = !1.
OMSEND.
GET FILE='C:\Temp\PairwiseSignTestPvalues.sav' /DROP=Command_ TO Label_.
DATASET NAME PValues.
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 & Adjusted p-values'
  /MISSING = VARIABLE
  /CELLS = NONE.
OMSEND.
DATASET ACTIVATE Data.
DATASET CLOSE PValues.
!ENDDEFINE.

* Macro 3 (Pairwise Wilcoxon tests - with exact p-value: EXACT TESTS MODULE NEEDED)*.

DEFINE MULTIWX (!POSITIONAL !CMDEND).
PRESERVE. /* The output language MUST be English - see the SELECT IF below *.
SET OLANG=ENGLISH.
DATASET NAME Data.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='Wilcoxon Ranks'
 /DESTINATION VIEWER=NO.
OMS /SELECT TABLES
 /IF COMMANDS='NPar Tests' SUBTYPES='Wilcoxon Test Statistics'
 /DESTINATION FORMAT=SAV
  OUTFILE='C:\Temp\PairwiseWilcoxonTestPvalues.sav'.
NPAR TEST
  /WILCOXON = !1
  /METHOD = EXACT TIMER(0).
OMSEND.
OMS /SELECT TABLES
 /IF COMMANDS='Summarize' SUBTYPES='Case Processing Summary'
 /DESTINATION VIEWER=NO.
GET FILE='C:\Temp\PairwiseWilcoxonTestPvalues.sav' /DROP=Command_ TO Label_.
DATASET NAME PValues.
* This step is language dependent (English is assumed) *.
SELECT IF Var1='Exact Sig. (2-tailed)'.
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).
 SUMMARIZE
  /TABLES = pvalue bonferr TO simes
  /FORMAT = LIST NOCASENUM TOTAL
  /TITLE = 'Original & Adjusted p-values'
  /MISSING = VARIABLE
  /CELLS = NONE.
OMSEND.
DATASET ACTIVATE Data.
DATASET CLOSE Pvalues.
RESTORE.
!ENDDEFINE.


* Sample dataset: "Effectiveness of terbutaline and theophylline alone and
  in combination in exercise-induced bronchospasm" GG Shapiro, JJ McPhillips,
  K Smith, CT Furukawa, WE Pierson and CW Bierman (1981) Pediatrics 67: 508-513 *.
DATA LIST LIST/placebo terbut teofil terteo(4 F8.2).
BEGIN DATA
3.18 3.27 3.34 3.87
0.80 2.73 2.18 2.28
2.13 2.12 2.13 2.26
1.35 1.92 2.28 2.72
2.10 2.06 2.53 2.73
1.30 2.45 2.35 4.41
3.30 3.10 3.07 3.86
2.93 2.49 2.93 3.25
3.50 3.03 3.32 2.78
0.55 2.26 1.96 2.08
0.33 0.37 1.12 1.05
1.55 2.48 2.23 2.38
2.10 2.12 2.31 2.46
3.49 3.63 3.70 3.63
0.70 0.83 1.57 1.27
2.20 3.06 2.92 2.95
1.35 2.32 2.12 2.65
1.73 2.74 2.68 2.92
4.68 4.54 4.63 4.81
0.83 1.16 2.04 2.70
0.43 0.93 2.70 3.16
END DATA.
VAR LABEL placebo 'FEV(lt) after Placebo'.
VAR LABEL terbut  'FEV(lt) after Terbutaline'.
VAR LABEL teofil  'FEV(lt) after Theophylline'.
VAR LABEL terteo  'FEV(lt) after Terbutaline+Theophylline'.

* MACRO calls *.
MULTIFR placebo TO terteo.
MULTISG placebo TO terteo.
MULTIWX placebo TO terteo.