SPSS script to Perform the Benjamini and Hochberg FDR procedure

[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