BB-SGoF: A correction for SGoF when the tests are correlated

Jacobo de Uña-Álvarez –University of Vigo

February 24, 2012

SGoF multitesting procedure (Carvajal-Rodríguez et al, 2009; de Uña-Álvarez, 2011) makes a decision on the number of effects among those p-values smaller than a significance level α. It is a useful alternative to FDR-based methods (e.g. Benjamini-Hochberg) when the number of tests in large, the proportion of effects is small, and the effects are weak to moderate. Unfortunately, SGoF is very sensitive to correlation among the tests and, indeed, it may be very anticonservative in dependent scenarios.

BB-SGoF (from Beta-Binomial SGoF, de Uña-Álvarez 2012) is a correction of SGoF for correlated tests. It is assumed that there exist k independent blocks of correlated p-values, where k is unknown. As SGoF, BB-SGoF makes a decision on the number of effects with p-values smaller than α, but depending on the number of blocks k and the within block correlation.

Here it is an application of BB-SGoF to Hedenfalk’s 3170 two-sample tests for gene expression levels, where α=0.05. The first Figure shows the statistical evidence on gene correlation: most of the p-values of Tarone’s test fall below the significance threshold. It suggest that there exist a number of blocks of correlated genes. The second Figure is the decision plot: it reveals that there are 389 true effects in the most conservative situation, which corresponds to k=266 blocks. Ordinary SGoF for independent tests would declare 412 effects, so BB-SGoF correction for dependence discards 23 potential discoveries.

Interestingly, BB-SGoF decision is still much more powerful than that obtained from standard methods which control the FDR at 5%. Indeed, Benjamini-Hochberg FDR-based method at that level gives for this data set only 157 discoveries (de Uña-Álvarez, 2011), which are less than half the discoveries declared by BB-SGoF.

$Descripción: Descripción: Descripción: Descripción: Descripción: Descripción: G:\bbsgof\hedenfalk_tarone_effects.eps$

References:

Carvajal Rodríguez, A., de Uña Álvarez, J., Rolán Álvarez, E. (2009). A new multitest correction (SGoF) that increases its statistical power when increasing the number of tests, BMC Bioinformatics 10:209, pp. 1-14.

de Uña-Álvarez J (2011) On the statistical properties of SGoF multitesting method. Statistical Applications in Genetics and Molecular Biology, Vol. 10, Iss. 1, Article 18.

de Uña-Álvarez J (2012) The Beta-Binomial SGoF method for multiple dependent tests. Statistical Applications in Genetics and Molecular Biology, Vol. 11, Iss. 3, Article 14.

R code for the application of BB-SGoF (updated: April 27^th, 2012):

[Remark, May 20^th, 2013: sgof package for R available at CRAN: http://cran.r-project.org/web/packages/sgof/ . This

package performs original SGoF and BB-SGoF methods, and it allows for the computation of adjusted p-values too.]

####code below may be run to analyze Hedenfalk data by starting with:

## library(qvalue)

## data(hedenfalk)

## u=hedenfalk

#####

#the only input you need is the sequence of p-values (u)

g=0.05 #the gamma parameter (how many effects are there among those p-values below gamma?)

a=0.05 #the alpha parameter (FWER control, weak sense)

v=as.numeric(u<=g)

n=length(v)

p=mean(v)

kmin=2 #minimum number of blocks

jmin=kmin-1

kmax=min(n%/%6,1000) #maximum number of blocks

jmax=kmax-1

Zvalue=vector(length=jmax)

pvalue=vector(length=jmax)

prob=vector(length=jmax)

rho=vector(length=jmax)

low=vector(length=jmax)

beta1=vector(length=jmax) #logit transform of prob

beta2=vector(length=jmax) #logit transfor of rho

sebeta2=vector(length=jmax) #std.err beta2

sebeta1=vector(length=jmax) #std.err beta1

library(VGAM)

for (j in jmin:jmax) {

k=j+1

A=vector(length=k)

n1=n%/%k

for (i in 1:k){A[i]=sum(v[((i-1)*n1+1):(i*n1)])}

A[k]=sum(v[((k-1)*n1+1):n])

B=c(rep(n1,k-1),length(v[((k-1)*n1+1):n]))

fit=vglm(cbind(A,B-A)~1,betabinomial)

pmodel=exp(coef(fit)[1])/(1+exp(coef(fit)[1]))

#Tarone statistic:#############

S=sum((A-p*B)^2)/(p*(1-p))

Zvalue[j]=(S-sum(B))/sqrt(2*sum(B*(B-1)))

pvalue[j]=1-pnorm(abs(Zvalue[j]))

################################

prob[j]=pmodel

beta1[j]=coef(fit)[1]

beta2[j]=coef(fit)[2]

sebeta1[j]=sqrt(vcov(fit)[1,1])

sebeta2[j]=sqrt(vcov(fit)[2,2])

rho[j]=exp(coef(fit)[2])/(1+exp(coef(fit)[2]))

low[j]=coef(fit)[1]-sqrt(vcov(fit))[1,1]*qnorm(1-a)

}

tlow=exp(low)/(1+exp(low))-g

effects=n*tlow

###control for too large standard errors when fitting the beta-binomial:

se=(log(prob/(1-prob))-low)/qnorm(1-a)

####se is equal to sebeta1

sec=length(se[se>25*median(se)|prob==0|prob==1])

sec #number of cases with a poor fit

#########################################

Mini=min(tlow[se<=25*median(se)&prob!=0&prob!=1])

mineffects=n*Mini #effects declared by BB-SGoF based on kN=k^N

Aux2=min(which(tlow==Mini))

kN=Aux2+1 #number of blocks for BB-SGoF decision

###################################################################

SGoF=n*(mean(v)-g)-n*sqrt(mean(v)*(1-mean(v))/n)*qnorm(1-a)+1

#effects declared by original SGoF for independent tests

#################################

#####numerical output of interest:

SGoF #effects declared by original SGoF (for independent tests)

mineffects #effects declared by BB-SGoF (most conservative decision, dependent tests)

kN #number of blocks

prob[Aux2] #estimated mean of the (random) probability of being a significant p-value

rho[Aux2];pvalue[Aux2] #estimated within-block correlation and significance

beta1[Aux2]

beta2[Aux2]

sebeta1[Aux2]

sebeta2[Aux2]

####logit transformation of p and rho, and standard errors

########for plotting the beta density we need:

alfa=(1-rho[Aux2])*prob[Aux2]/rho[Aux2]

beta=(1-rho[Aux2])*(1-prob[Aux2])/rho[Aux2]

ss=seq(0.001,0.999,by=0.001)

##################################

ww=(jmin:jmax)[se<=25*median(se)&prob!=0&prob!=1]+1 #index vector for non-poor fit of the betabinomial model

par(mfrow=c(2,2)) #plotting your results:

plot(ww,pvalue[se<=25*median(se)&prob!=0&prob!=1],ylab="pvalue",xlab="blocks",main="Tarone's test");abline(h=0.05,col=2,lty=2)

plot(ww,rho[se<=25*median(se)&prob!=0&prob!=1],ylab="rho",xlab="blocks",main="Within-block correlation")

plot(effects[se<=25*median(se)&prob!=0&prob!=1],ylab="effects",xlab="blocks",main="Decision plot");abline(h=SGoF,col=2,lty=2)

plot(ss,dbeta(ss,alfa,beta),type="l",ylab="density",xlab="probability",main="Beta density"); abline(v=prob[Aux2],lty=2)

################################# this is the BB-SGoF 2x2 plot you get for Hedenfalk data:

Cuadro de texto:

####and these are the numerical outputs:

> #################################

> SGoF #effects declared by original SGoF (for independent tests)

[1] 412.0839

> minieffects #effects declared by BB-SGoF (most conservative decision, dependent tests)

[1] 389.1544

> kN #number of blocks

[1] 266

> prob[Aux2] #estimated mean of the (random) probability of being a significant p-value

[1] 0.1856618

> rho[Aux2];pvalue[Aux2] #estimated within-block correlation and significance

[1] 0.02711878

[1] 4.863077e-11

> beta1[Aux2]

[1] -1.478449

> beta2[Aux2]

[1] -3.580035

> sebeta1[Aux2]

[1] 0.05333666

> sebeta2[Aux2]

[1] 0.3868749

> ####logit transformation of p and rho, and standard errors