docs: def forward_lmm(snps,pheno,K=None,covs=None,qvalues=False,threshold = 5e-8, maxiter = 2,test='lrt',**kw_args):
      """
      univariate fixed effects test with forward selection
  
      Args:
          snps:   [N x S] SP.array of S SNPs for N individuals (test SNPs)
          pheno:  [N x 1] SP.array of 1 phenotype for N individuals
          K:      [N x N] SP.array of LMM-covariance/kinship koefficients (optional)
                          If not provided, then linear regression analysis is performed
          covs:   [N x D] SP.array of D covariates for N individuals
          threshold:      (float) P-value thrashold for inclusion in forward selection (default 5e-8)
          maxiter:        (int) maximum number of interaction scans. First scan is
                          without inclusion, so maxiter-1 inclusions can be performed. (default 2)
          test:           'lrt' for likelihood ratio test (default) or 'f' for F-test
  
      Returns:
          lm:             limix LMM object
          iadded:         array of indices of SNPs included in order of inclusion
          pvadded:        array of Pvalues obtained by the included SNPs in iteration
                          before inclusion
          pvall:   [maxiter x S] SP.array of Pvalues for all iterations
      """
  
      if K is None:
          K=SP.eye(snps.shape[0])
      if covs is None:
          covs = SP.ones((snps.shape[0],1))
  
      lm = simple_lmm(snps,pheno,K=K,covs=covs,test=test,**kw_args)
      pvall = SP.zeros((maxiter,snps.shape[1]))
      pv = lm.getPv()
      pvall[0:1,:]=pv
      imin= pv.argmin()
      niter = 1
      #start stuff
      iadded = []
      pvadded = []
      qvadded = []
      if qvalues:
          assert pv.shape[0]==1, "This is untested with the fdr package. pv.shape[0]==1 failed"
          qvall = SP.zeros((maxiter,snps.shape[1]))
          qv  = FDR.qvalues(pv)
          qvall[0:1,:] = qv
          score=qv.min()
      else:
          score=pv.min()
      while (score<threshold) and niter<maxiter:
          t0=time.time()
          iadded.append(imin)
          pvadded.append(pv[0,imin])
          if qvalues:
              qvadded.append(qv[0,imin])
          covs=SP.concatenate((covs,snps[:,imin:(imin+1)]),1)
          lm.setCovs(covs)
          lm.process()
          pv = lm.getPv()
          pvall[niter:niter+1,:]=pv
          imin= pv.argmin()
          if qvalues:
              qv = FDR.qvalues(pv)
              qvall[niter:niter+1,:] = qv
              score = qv.min()
          else:
              score = pv.min()
          t1=time.time()
          print(("finished GWAS testing in %.2f seconds" %(t1-t0)))
          niter=niter+1
      RV = {}
      RV['iadded']  = iadded
      RV['pvadded'] = pvadded
      RV['pvall']   = pvall
      if qvalues:
          RV['qvall'] = qvall
          RV['qvadded'] = qvadded
      return lm,RV