docs: def fdrcorrection(pvals, alpha=0.05, method='indep'):
      '''
      NOTE: This function was copied from
      statsmodels.sandbox.stats.multicomp.fdrcorrection0, from statsmodels
      version 0.5.0.
  
      This is to avoid requiring all of statsmodels to be a dependency for
      metaseq, just for this function.
  
  
  
  
      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')
  
      Returns
      -------
      rejected : array, bool
          True if a hypothesis is rejected, False if not
      pvalue-corrected : array
          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)
  
      pvals_sortind = np.argsort(pvals)
      pvals_sorted = pvals[pvals_sortind]
      sortrevind = pvals_sortind.argsort()
  
      if method in ['i', 'indep', 'p', 'poscorr']:
          ecdffactor = _ecdf(pvals_sorted)
      elif method in ['n', 'negcorr']:
          cm = np.sum(1./np.arange(1, len(pvals_sorted)+1))  # corrected this
          ecdffactor = _ecdf(pvals_sorted) / cm
  #    elif method in ['n', 'negcorr']:
  #        cm = np.sum(np.arange(len(pvals)))
  #        ecdffactor = ecdf(pvals_sorted)/cm
      else:
          raise ValueError('only indep and necorr implemented')
      reject = pvals_sorted <= ecdffactor*alpha
      if reject.any():
          rejectmax = max(np.nonzero(reject)[0])
          reject[:rejectmax] = True
  
      pvals_corrected_raw = pvals_sorted / ecdffactor
      pvals_corrected = np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]
      pvals_corrected[pvals_corrected > 1] = 1
      return reject[sortrevind], pvals_corrected[sortrevind]