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