SPSS script to Perform the Benjamini and Hochberg FDR procedure

This procedure may also be computed using the [:FAQ/pvs:spreadsheet here]. Keselman HJ, Cribbie R and Holland B (1999) compared several pairwise comparison tests and found the FDR approach to be the most powerful when more than 6 pairwise tests were performed e.g. when comparing all possible pairs of means for five groups.

The FDR is a stepdown procedure using a prescribed alpha usually 0.05 for testing in advance of the p-value testing and so the above spreadsheet should be used for hypothesis testing of multiple p-values. See [:FAQ/SpssBonferroni: here] for a warning about using exact p-values instead of statements of form p<alpha.

To compute, however, the exact p-value which would produce specific p-values of form p= rather than p<alpha (for assessing sensitivity to see how close a p-value's unadjusted version is to statistical significance) this [attachment:fdrexact.xls spreadsheet] can be used.

Reference

Keselman HJ, Cribbie R and Holland B (1999) The pairwise multiple comparison multiplicity problem: an alternative approach to familywise and comparisonwise type I error control Psychological methods 4(1) 58-69.

[COPY AND PASTE THE BELOW BOX SYNTAX INTO A SPSS SYNTAX WINDOW; SELECT ALL AND RUN. EDIT THE INPUT DATA AS REQUIRED]

DATA LIST FREE / pvals .
BEGIN DATA
0.021
0.001
0.017
0.041
0.005
0.036
0.042
0.023
0.07
0.1
END DATA .

COMPUTE alpha = 0.05 .

COMPUTE index = $CASENUM .
EXECUTE .

SORT CASES BY pvals (A) .

COMPUTE rank = index .
COMPUTE ref = alpha * rank / 10 .
COMPUTE diff = pvals - ref .
COMPUTE comp = diff LT 0 .
CREATE ccomp = CSUM(comp) .
EXECUTE .

MATRIX .
GET rnk / VARIABLES = rank .
GET pvals / VARIABLES = pvals .
GET ccomp / VARIABLES = ccomp .
GET alpha / VARIABLES = alpha .
GET diff / VARIABLES = diff .
COMPUTE ccmpmx = pvals LE CMAX(((CMAX(ccomp) &* (diff LE 0)) EQ ccomp) &* pvals) .
SAVE {rnk, pvals, alpha, ccmpmx} / OUTFILE = * / VARIABLES = index,pvals,alpha,fdrsig .
END MATRIX .

SORT CASES BY index (A) .

EXPORT OUTFILE = out .
IMPORT FILE = out / KEEP = pvals alpha fdrsig .

FORMAT pvals (F5.3) .
FORMAT alpha (F5.3) .
FORMAT fdrsig (F1.0) .

VALUE LABELS fdrsig 1 'Sig' 0 'Nonsig' .

SET RESULTS = LISTING .

DO IF $CASENUM EQ 1 .
PRINT EJECT /'P Value' 1 'FDR Criterion' 9 'FDR Significance' 24 .
END IF .
PRINT / pvals (T3 F5.3) alpha (T17, F5.3) fdrsig (T39, F1.0) .
EXECUTE .

This produces the following output:

P Value FDR Criterion  FDR Significance
   .021          .050                 1
   .001          .050                 1
   .017          .050                 1
   .041          .050                 0
   .005          .050                 1
   .036          .050                 0
   .042          .050                 0
   .023          .050                 1
   .070          .050                 0
   .100          .050                 0

R Code to Perform the Benjamini and Hochberg FDR procedure

[Insert P values in any order]

pvalue <- c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1) 
sorted.pvalue<-sort(pvalue) 
sorted.pvalue 
j.alpha <- (1:10)*(.05/10) 
diff <- sorted.pvalue-j.alpha 
neg.diff <- diff[diff<0] 
pos.diff <- neg.diff[length(neg.diff)] 
index <- diff==pos.diff 
p.cutoff <-sorted.pvalue[index] 
p.cutoff 
p.sig <- pvalue[pvalue <= p.cutoff] 
p.sig

R Output

> pvalue<-c(0.021,0.001,0.017,0.041,0.005,0.036,0.042,0.023,0.07,0.1) 
> sorted.pvalue <- sort(pvalue) 
> sorted.pvalue

[sorted P values in ascending order]
[1] 0.001 0.005 0.017 0.021 0.023 0.036 0.041 0.042 0.070 0.100 

> j.alpha <- (1:10) * (0.05/10) 
> diff <- sorted.pvalue - j.alpha 
> neg.diff <- diff[diff < 0] 
> pos.diff <- neg.diff[length(neg.diff)] 
> index <- diff == pos.diff 
> p.cutoff <- sorted.pvalue[index] 
> p.cutoff 
[1] 0.023 

> p.sig <- pvalue[pvalue <= p.cutoff] > p.sig

[significant p values by B-H method]: 
[1] 0.021 0.001 0.017 0.005 0.023