[docs]defmultipletests(pvals,alpha=0.05,method='hs',is_sorted=False,returnsorted=False):""" Test results and p-value correction for multiple tests Parameters ---------- pvals : array_like, 1-d uncorrected p-values. Must be 1-dimensional. alpha : float FWER, family-wise error rate, e.g. 0.1 method : str Method used for testing and adjustment of pvalues. Can be either the full name or initial letters. Available methods are: - `bonferroni` : one-step correction - `sidak` : one-step correction - `holm-sidak` : step down method using Sidak adjustments - `holm` : step-down method using Bonferroni adjustments - `simes-hochberg` : step-up method (independent) - `hommel` : closed method based on Simes tests (non-negative) - `fdr_bh` : Benjamini/Hochberg (non-negative) - `fdr_by` : Benjamini/Yekutieli (negative) - `fdr_tsbh` : two stage fdr correction (non-negative) - `fdr_tsbky` : two stage fdr correction (non-negative) is_sorted : bool If False (default), the p_values will be sorted, but the corrected pvalues are in the original order. If True, then it assumed that the pvalues are already sorted in ascending order. returnsorted : bool not tested, return sorted p-values instead of original sequence Returns ------- reject : ndarray, boolean true for hypothesis that can be rejected for given alpha pvals_corrected : ndarray p-values corrected for multiple tests alphacSidak: float corrected alpha for Sidak method alphacBonf: float corrected alpha for Bonferroni method Notes ----- There may be API changes for this function in the future. Except for 'fdr_twostage', the p-value correction is independent of the alpha specified as argument. In these cases the corrected p-values can also be compared with a different alpha. In the case of 'fdr_twostage', the corrected p-values are specific to the given alpha, see ``fdrcorrection_twostage``. The 'fdr_gbs' procedure is not verified against another package, p-values are derived from scratch and are not derived in the reference. In Monte Carlo experiments the method worked correctly and maintained the false discovery rate. All procedures that are included, control FWER or FDR in the independent case, and most are robust in the positively correlated case. `fdr_gbs`: high power, fdr control for independent case and only small violation in positively correlated case **Timing**: Most of the time with large arrays is spent in `argsort`. When we want to calculate the p-value for several methods, then it is more efficient to presort the pvalues, and put the results back into the original order outside of the function. Method='hommel' is very slow for large arrays, since it requires the evaluation of n partitions, where n is the number of p-values. """importgcpvals=np.asarray(pvals)alphaf=alpha# Notation ?ifnotis_sorted:sortind=np.argsort(pvals)pvals=np.take(pvals,sortind)ntests=len(pvals)alphacSidak=1-np.power((1.-alphaf),1./ntests)alphacBonf=alphaf/float(ntests)ifmethod.lower()in['b','bonf','bonferroni']:reject=pvals<=alphacBonfpvals_corrected=pvals*float(ntests)elifmethod.lower()in['s','sidak']:reject=pvals<=alphacSidakpvals_corrected=1-np.power((1.-pvals),ntests)elifmethod.lower()in['hs','holm-sidak']:alphacSidak_all=1-np.power((1.-alphaf),1./np.arange(ntests,0,-1))notreject=pvals>alphacSidak_alldelalphacSidak_allnr_index=np.nonzero(notreject)[0]ifnr_index.size==0:# nonreject is empty, all rejectednotrejectmin=len(pvals)else:notrejectmin=np.min(nr_index)notreject[notrejectmin:]=Truereject=~notrejectdelnotrejectpvals_corrected_raw=1-np.power((1.-pvals),np.arange(ntests,0,-1))pvals_corrected=np.maximum.accumulate(pvals_corrected_raw)delpvals_corrected_rawelifmethod.lower()in['h','holm']:notreject=pvals>alphaf/np.arange(ntests,0,-1)nr_index=np.nonzero(notreject)[0]ifnr_index.size==0:# nonreject is empty, all rejectednotrejectmin=len(pvals)else:notrejectmin=np.min(nr_index)notreject[notrejectmin:]=Truereject=~notrejectpvals_corrected_raw=pvals*np.arange(ntests,0,-1)pvals_corrected=np.maximum.accumulate(pvals_corrected_raw)delpvals_corrected_rawgc.collect()elifmethod.lower()in['sh','simes-hochberg']:alphash=alphaf/np.arange(ntests,0,-1)reject=pvals<=alphashrejind=np.nonzero(reject)ifrejind[0].size>0:rejectmax=np.max(np.nonzero(reject))reject[:rejectmax]=Truepvals_corrected_raw=np.arange(ntests,0,-1)*pvalspvals_corrected=np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]delpvals_corrected_rawelifmethod.lower()in['ho','hommel']:# we need a copy because we overwrite it in a loopa=pvals.copy()forminrange(ntests,1,-1):cim=np.min(m*pvals[-m:]/np.arange(1,m+1.))a[-m:]=np.maximum(a[-m:],cim)a[:-m]=np.maximum(a[:-m],np.minimum(m*pvals[:-m],cim))pvals_corrected=areject=a<=alphafelifmethod.lower()in['fdr_bh','fdr_i','fdr_p','fdri','fdrp']:# delegate, call with sorted pvalsreject,pvals_corrected=fdrcorrection(pvals,alpha=alpha,method='indep',is_sorted=True)elifmethod.lower()in['fdr_by','fdr_n','fdr_c','fdrn','fdrcorr']:# delegate, call with sorted pvalsreject,pvals_corrected=fdrcorrection(pvals,alpha=alpha,method='n',is_sorted=True)elifmethod.lower()in['fdr_tsbky','fdr_2sbky','fdr_twostage']:# delegate, call with sorted pvalsreject,pvals_corrected=fdrcorrection_twostage(pvals,alpha=alpha,method='bky',is_sorted=True)[:2]elifmethod.lower()in['fdr_tsbh','fdr_2sbh']:# delegate, call with sorted pvalsreject,pvals_corrected=fdrcorrection_twostage(pvals,alpha=alpha,method='bh',is_sorted=True)[:2]elifmethod.lower()in['fdr_gbs']:#adaptive stepdown in Gavrilov, Benjamini, Sarkar, Annals of Statistics 2009## notreject = pvals > alphaf / np.arange(ntests, 0, -1) #alphacSidak## notrejectmin = np.min(np.nonzero(notreject))## notreject[notrejectmin:] = True## reject = ~notrejectii=np.arange(1,ntests+1)q=(ntests+1.-ii)/ii*pvals/(1.-pvals)pvals_corrected_raw=np.maximum.accumulate(q)#up requirementdpvals_corrected=np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]delpvals_corrected_rawreject=pvals_corrected<=alphaelse:raiseValueError('method not recognized')ifpvals_correctedisnotNone:#not necessary anymorepvals_corrected[pvals_corrected>1]=1ifis_sortedorreturnsorted:returnreject,pvals_corrected,alphacSidak,alphacBonfelse:pvals_corrected_=np.empty_like(pvals_corrected)pvals_corrected_[sortind]=pvals_correcteddelpvals_correctedreject_=np.empty_like(reject)reject_[sortind]=rejectreturnreject_,pvals_corrected_,alphacSidak,alphacBonf

[docs]deffdrcorrection(pvals,alpha=0.05,method='indep',is_sorted=False):'''pvalue correction for false discovery rate This covers Benjamini/Hochberg for independent or positively correlated and Benjamini/Yekutieli for general or negatively correlated tests. Both are available in the function multipletests, as method=`fdr_bh`, resp. `fdr_by`. Parameters ---------- pvals : array_like set of p-values of the individual tests. alpha : float error rate method : {'indep', 'negcorr'} is_sorted : bool If False (default), the p_values will be sorted, but the corrected pvalues are in the original order. If True, then it assumed that the pvalues are already sorted in ascending order. Returns ------- rejected : ndarray, bool True if a hypothesis is rejected, False if not pvalue-corrected : ndarray pvalues adjusted for multiple hypothesis testing to limit FDR Notes ----- If there is prior information on the fraction of true hypothesis, then alpha should be set to alpha * m/m_0 where m is the number of tests, given by the p-values, and m_0 is an estimate of the true hypothesis. (see Benjamini, Krieger and Yekuteli) The two-step method of Benjamini, Krieger and Yekutiel that estimates the number of false hypotheses will be available (soon). Method names can be abbreviated to first letter, 'i' or 'p' for fdr_bh and 'n' for fdr_by. '''pvals=np.asarray(pvals)ifnotis_sorted:pvals_sortind=np.argsort(pvals)pvals_sorted=np.take(pvals,pvals_sortind)else:pvals_sorted=pvals# aliasifmethodin['i','indep','p','poscorr']:ecdffactor=_ecdf(pvals_sorted)elifmethodin['n','negcorr']:cm=np.sum(1./np.arange(1,len(pvals_sorted)+1))#corrected thisecdffactor=_ecdf(pvals_sorted)/cm## elif method in ['n', 'negcorr']:## cm = np.sum(np.arange(len(pvals)))## ecdffactor = ecdf(pvals_sorted)/cmelse:raiseValueError('only indep and negcorr implemented')reject=pvals_sorted<=ecdffactor*alphaifreject.any():rejectmax=max(np.nonzero(reject)[0])reject[:rejectmax]=Truepvals_corrected_raw=pvals_sorted/ecdffactorpvals_corrected=np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]delpvals_corrected_rawpvals_corrected[pvals_corrected>1]=1ifnotis_sorted:pvals_corrected_=np.empty_like(pvals_corrected)pvals_corrected_[pvals_sortind]=pvals_correcteddelpvals_correctedreject_=np.empty_like(reject)reject_[pvals_sortind]=rejectreturnreject_,pvals_corrected_else:returnreject,pvals_corrected

[docs]deffdrcorrection_twostage(pvals,alpha=0.05,method='bky',iter=False,is_sorted=False):'''(iterated) two stage linear step-up procedure with estimation of number of true hypotheses Benjamini, Krieger and Yekuteli, procedure in Definition 6 Parameters ---------- pvals : array_like set of p-values of the individual tests. alpha : float error rate method : {'bky', 'bh') see Notes for details * 'bky' - implements the procedure in Definition 6 of Benjamini, Krieger and Yekuteli 2006 * 'bh' - the two stage method of Benjamini and Hochberg iter : bool Returns ------- rejected : ndarray, bool True if a hypothesis is rejected, False if not pvalue-corrected : ndarray pvalues adjusted for multiple hypotheses testing to limit FDR m0 : int ntest - rej, estimated number of true hypotheses alpha_stages : list of floats A list of alphas that have been used at each stage Notes ----- The returned corrected p-values are specific to the given alpha, they cannot be used for a different alpha. The returned corrected p-values are from the last stage of the fdr_bh linear step-up procedure (fdrcorrection0 with method='indep') corrected for the estimated fraction of true hypotheses. This means that the rejection decision can be obtained with ``pval_corrected <= alpha``, where ``alpha`` is the original significance level. (Note: This has changed from earlier versions (<0.5.0) of statsmodels.) BKY described several other multi-stage methods, which would be easy to implement. However, in their simulation the simple two-stage method (with iter=False) was the most robust to the presence of positive correlation TODO: What should be returned? '''pvals=np.asarray(pvals)ifnotis_sorted:pvals_sortind=np.argsort(pvals)pvals=np.take(pvals,pvals_sortind)ntests=len(pvals)ifmethod=='bky':fact=(1.+alpha)alpha_prime=alpha/factelifmethod=='bh':fact=1.alpha_prime=alphaelse:raiseValueError("only 'bky' and 'bh' are available as method")alpha_stages=[alpha_prime]rej,pvalscorr=fdrcorrection(pvals,alpha=alpha_prime,method='indep',is_sorted=True)r1=rej.sum()if(r1==0)or(r1==ntests):returnrej,pvalscorr*fact,ntests-r1,alpha_stagesri_old=r1whileTrue:ntests0=1.0*ntests-ri_oldalpha_star=alpha_prime*ntests/ntests0alpha_stages.append(alpha_star)#print ntests0, alpha_starrej,pvalscorr=fdrcorrection(pvals,alpha=alpha_star,method='indep',is_sorted=True)ri=rej.sum()if(notiter)orri==ri_old:breakelifri<ri_old:# prevent cycles and endless loopsraiseRuntimeError(" oops - should not be here")ri_old=ri# make adjustment to pvalscorr to reflect estimated number of Non-Null cases# decision is then pvalscorr < alpha (or <=)pvalscorr*=ntests0*1.0/ntestsifmethod=='bky':pvalscorr*=(1.+alpha)ifnotis_sorted:pvalscorr_=np.empty_like(pvalscorr)pvalscorr_[pvals_sortind]=pvalscorrdelpvalscorrreject=np.empty_like(rej)reject[pvals_sortind]=rejreturnreject,pvalscorr_,ntests-ri,alpha_stageselse:returnrej,pvalscorr,ntests-ri,alpha_stages

[docs]classNullDistribution(object):""" Estimate a Gaussian distribution for the null Z-scores. The observed Z-scores consist of both null and non-null values. The fitted distribution of null Z-scores is Gaussian, but may have non-zero mean and/or non-unit scale. Parameters ---------- zscores : array_like The observed Z-scores. null_lb : float Z-scores between `null_lb` and `null_ub` are all considered to be true null hypotheses. null_ub : float See `null_lb`. estimate_mean : bool If True, estimate the mean of the distribution. If False, the mean is fixed at zero. estimate_scale : bool If True, estimate the scale of the distribution. If False, the scale parameter is fixed at 1. estimate_null_proportion : bool If True, estimate the proportion of true null hypotheses (i.e. the proportion of z-scores with expected value zero). If False, this parameter is fixed at 1. Attributes ---------- mean : float The estimated mean of the empirical null distribution sd : float The estimated standard deviation of the empirical null distribution null_proportion : float The estimated proportion of true null hypotheses among all hypotheses References ---------- B Efron (2008). Microarrays, Empirical Bayes, and the Two-Groups Model. Statistical Science 23:1, 1-22. Notes ----- See also: http://nipy.org/nipy/labs/enn.html#nipy.algorithms.statistics.empirical_pvalue.NormalEmpiricalNull.fdr """def__init__(self,zscores,null_lb=-1,null_ub=1,estimate_mean=True,estimate_scale=True,estimate_null_proportion=False):# Extract the null z-scoresii=np.flatnonzero((zscores>=null_lb)&(zscores<=null_ub))iflen(ii)==0:raiseRuntimeError("No Z-scores fall between null_lb and null_ub")zscores0=zscores[ii]# Number of Z-scores, and null Z-scoresn_zs,n_zs0=len(zscores),len(zscores0)# Unpack and transform the parameters to the natural scale, hold# parameters fixed as specified.defxform(params):mean=0.sd=1.prob=1.ii=0ifestimate_mean:mean=params[ii]ii+=1ifestimate_scale:sd=np.exp(params[ii])ii+=1ifestimate_null_proportion:prob=1/(1+np.exp(-params[ii]))returnmean,sd,probfromscipy.stats.distributionsimportnormdeffun(params):""" Negative log-likelihood of z-scores. The function has three arguments, packed into a vector: mean : location parameter logscale : log of the scale parameter logitprop : logit of the proportion of true nulls The implementation follows section 4 from Efron 2008. """d,s,p=xform(params)# Mass within the central regioncentral_mass=(norm.cdf((null_ub-d)/s)-norm.cdf((null_lb-d)/s))# Probability that a Z-score is null and is in the central regioncp=p*central_mass# Binomial termrval=n_zs0*np.log(cp)+(n_zs-n_zs0)*np.log(1-cp)# Truncated Gaussian term for null Z-scoreszv=(zscores0-d)/srval+=np.sum(-zv**2/2)-n_zs0*np.log(s)rval-=n_zs0*np.log(central_mass)return-rval# Estimate the parametersfromscipy.optimizeimportminimize# starting values are mean = 0, scale = 1, p0 ~ 1mz=minimize(fun,np.r_[0.,0,3],method="Nelder-Mead")mean,sd,prob=xform(mz['x'])self.mean=meanself.sd=sdself.null_proportion=prob# The fitted null density function

[docs]defpdf(self,zscores):""" Evaluates the fitted empirical null Z-score density. Parameters ---------- zscores : scalar or array_like The point or points at which the density is to be evaluated. Returns ------- The empirical null Z-score density evaluated at the given points. """zval=(zscores-self.mean)/self.sdreturnnp.exp(-0.5*zval**2-np.log(self.sd)-0.5*np.log(2*np.pi))