Mercurial > repos > dereeper > mlmm

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/MLMM.pl	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,155 @@
+#!/usr/bin/perl
+
+use strict;
+use Switch;
+use Getopt::Long;
+use Bio::SeqIO;
+
+
+my $usage = qq~Usage:$0 <args> [<opts>]
+where <args> are:
+    -g, --geno         <Genotype input>
+    -i, --info         <SNP information. Genome position.>
+    -p, --pheno        <Phenotype input>
+    -o, --out          <output name>
+    -d, --directory    <directory for MLMM R libraries>
+    -s, --step_number  <number of steps. Maximum: 20. Default: 10>
+    -m, --method       <Method: mbonf or extBIC. Default: mbonf>
+~;
+$usage .= "\n";
+
+my ($geno,$map,$pheno,$out,$dir,$steps,$method);
+
+
+GetOptions(
+	"geno=s"      => \$geno,
+	"info=s"      => \$map,
+	"pheno=s"     => \$pheno,
+	"out=s"       => \$out,
+	"dir=s"       => \$dir,
+	"steps=s"     => \$steps,
+	"method=s"    => \$method
+);
+
+
+die $usage
+  if ( !$geno || !$map || !$pheno || !$out || !$dir || !$steps || !$method);
+
+my $max_steps = 10;
+my $plot_opt = "mbonf";
+if ($method && $method ne 'mbonf' && $method ne 'extBIC')
+{
+	print "Aborted: Method must be mbonf or extBIC.\n";
+	exit;
+}
+else
+{
+	$plot_opt = $method;
+}
+if ($steps && $steps !~/\d+/ && $steps > 20 && $steps < 2)
+{
+	print "Aborted: Number of steps must be greater than 2 and lower than 20.\n";
+	exit;
+}
+else
+{
+	$max_steps = $steps;
+}
+
+
+my $chunk = 2;
+
+my $RSCRIPT_EXE = "Rscript";
+my $R_DIR = $dir;
+
+
+my $head_trait = `head -1 $pheno`;
+my @headers_traits = split(/\t/,$head_trait);
+my $trait_name = $headers_traits[1];
+
+
+open( my $RCMD, ">rscript" ) or throw Error::Simple($!);
+
+
+
+
+print $RCMD "Y_file <- \"" . $pheno . "\"\n";
+print $RCMD "X_file <- \"" . $geno . "\"\n";
+if($map)
+{
+	print $RCMD "map_file <- \"$map\"\n";
+	print $RCMD "map <- read.table(map_file, sep = \"\\t\", header = T)\n";
+}
+print $RCMD "mlmm_data = list()\n";
+print $RCMD "mlmm_data\$chunk <- $chunk\n";
+print $RCMD "mlmm_data\$maxsteps <- $max_steps\n";
+print $RCMD "genot <- read.table(X_file, sep = \"\\t\", header = T)\n";
+print $RCMD "genot_mat <- as.matrix(genot[, 2:ncol(genot)])\n";
+print $RCMD "rownames(genot_mat) <- genot\$Ind_id\n";
+
+print $RCMD "phenot <- read.table(Y_file, sep = \"\\t\", header = T)\n";
+
+
+
+# missing data imputation
+print $RCMD "genot_imp <- genot_mat\n";
+print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n";
+print $RCMD "for (i in 1:ncol(genot_imp)){genot_imp[is.na(genot_imp[,i]), i] <- average[i]}\n";
+
+# kinship matrix computation
+print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n";
+print $RCMD "stdev <- apply(genot_imp, 2, sd)\n";
+print $RCMD "genot_stand <- sweep(sweep(genot_imp, 2, average, \"-\"), 2, stdev, \"/\")\n";
+print $RCMD "K_mat <- (genot_stand %*% t(genot_stand)) / ncol(genot_stand)\n";
+print $RCMD "write.table(K_mat, '$out.kinship', sep='\\t', dec='.', quote=F, col.names=T, row.names=T)\n";
+
+print $RCMD "source(\"" . $R_DIR. "/mlmm.r\")\n";
+print $RCMD "source(\"" . $R_DIR. "/emma.r\")\n";
+
+# mlmm
+print $RCMD "mygwas <- mlmm(Y = phenot\$$trait_name, X = genot_imp, K = K_mat, nbchunks=mlmm_data\$chunk, maxsteps=mlmm_data\$maxsteps)\n";
+
+# plots
+print $RCMD "pdf('$out.pdf')\n";
+print $RCMD "plot_step_table(mygwas, \"h2\")\n";
+print $RCMD "plot_step_table(mygwas, \"extBIC\")\n";
+print $RCMD "plot_step_table(mygwas, \"maxpval\")\n";
+print $RCMD "plot_step_RSS(mygwas)\n";
+# for (my $j = 1; $j <= ($max_steps - 1); $j++)
+# {
+	# print $RCMD "plot_fwd_GWAS(mygwas, step = $j, snp_info = map, pval_filt = 0.1)\n";
+# }
+print $RCMD "plot_opt_GWAS(mygwas, opt = \"extBIC\", snp_info = map, pval_filt = 0.1)\n";
+print $RCMD "plot_opt_GWAS(mygwas, opt = \"mbonf\", snp_info = map, pval_filt = 0.1)\n";
+#print $RCMD "qqplot_fwd_GWAS(mygwas, nsteps = mlmm_data\$maxsteps)\n";
+print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"extBIC\")\n";
+print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"mbonf\")\n";
+
+# outputs
+print $RCMD "write.table(mygwas\$RSSout, '$out.rss', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n";
+print $RCMD "write.table(mygwas\$step_table, '$out.steptable', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n";
+
+$plot_opt = "\$opt_" . $plot_opt;
+print $RCMD "pval = mygwas" . $plot_opt . "\$out\n";
+print $RCMD "colnames(pval) = c(\"Marker_name\", \"Pvalue\")\n";
+print $RCMD "info_tmp = map\n";
+print $RCMD "colnames(info_tmp) = c(\"Marker_name\", \"Chr\", \"Pos\")\n";
+print $RCMD "res_asso = pval\n";
+print $RCMD qq~
+if(exists("info_tmp")){
+	res_asso = merge(info_tmp, res_asso, by="Marker_name")
+	if( !is.element("Trait", colnames(info_tmp)) ){
+		m = matrix(data="traitname", ncol=1, nrow=nrow(res_asso), dimnames=list(c(), c("Trait")))
+		res_asso = cbind(m, res_asso)
+	}
+}
+~;
+print $RCMD "res_asso = res_asso[order(res_asso[, \"Trait\"], res_asso[, \"Chr\"], res_asso[, \"Pos\"]), ]\n";
+print $RCMD "write.table(res_asso, '$out.res_asso', sep='\t', dec='.', quote=F, col.names=T, row.names=F)\n";
+close($RCMD);
+
+system("$RSCRIPT_EXE --vanilla rscript");
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/MLMM.sh	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,28 @@
+#!/bin/bash
+geno=$1
+map=$2
+pheno=$3
+steps=$4
+method=$5
+output=$6
+pdf=$7
+kinship=$8
+rss=$9
+step_table=${10}
+log=${11}
+
+directory=`dirname $0`
+mkdir tmpdir$$
+cp -rf $geno tmpdir$$/geno
+cp -rf $map tmpdir$$/map
+cp -rf $pheno tmpdir$$/pheno
+
+perl $directory/MLMM.pl -g tmpdir$$/geno -i tmpdir$$/map -p tmpdir$$/pheno -s $steps -m $method -o tmpdir$$/output -d $directory/source_library >>$log 2>&1
+
+mv tmpdir$$/output.pdf $pdf
+mv tmpdir$$/output.kinship $kinship
+mv tmpdir$$/output.res_asso $output
+mv tmpdir$$/output.rss $rss
+mv tmpdir$$/output.steptable $step_table
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/MLMM.xml	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,61 @@
+<tool id="mlmm" name="MLMM" version="1.0">
+	<description>GWAS using Multi-Locus Mixed-Model (MLMM)</description>
+	<requirements>
+		<requirement type="binary">Rscript</requirement>
+	</requirements>
+	<command interpreter="bash">./MLMM.sh $geno $map $pheno $steps $method $output $pdf $kinship $rss $step_table $log
+    </command>
+	<inputs>
+		<param format="txt" name="geno" type="data" label="Genotype matrix" help="NxM, N = individuals in line, M = Markers in columns, Genotype coded in 0,1,2"/>
+		<param type="data" format="txt" name="map" label="SNP Information file" help="3 columns: SNP, Chrom, Pos"/>
+		<param format="txt" name="pheno" type="data" label="Phenotype matrix" help="NxT, N = individuals in line, T = Trait in columns (Phenot1, Phenot2...)"/>
+		<param type="text" name="steps" label="Maximum number of steps for the forward approach" value="10"/>
+		<param name="method" type="select">
+                        <option value="extBIC">EBIC</option>
+                        <option value="mbonf" selected="True">MBonf</option>
+                </param>
+	</inputs>
+	<outputs>
+		<data format="txt" name="output" label="Association results"/>
+		<data format="txt" name="kinship" label="Kinship matrix"/>
+		<data format="pdf" name="pdf" label="PDF Graphical outputs"/>
+		<data format="txt" name="rss" label="RSS"/>
+		<data format="txt" name="step_table" label="Step Table"/>
+		<data format="txt" name="log" label="Log file"/>
+	</outputs>
+	<help>
+
+
+.. class:: infomark
+
+**Program encapsulated in Galaxy by Southgreen**
+
+.. class:: infomark
+
+**MLMM version 1.0**
+
+-----
+
+==============
+ Please cite:
+==============
+
+"An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations.", **Segura V, Vilhjalmsson BJ, Platt A, Korte A, Seren U, Long Q, Nordborg M.**, Nature Genetics, 44: 825-830, 2012.
+
+-----
+
+===========
+ Overview:
+===========
+
+MLMM is an efficient multi-locus mixed-model approach for genome-wide association studies in structured populations.
+
+-----
+
+For further informations, please visite the MLMM_ website.
+
+
+.. _MLMM: https://sites.google.com/site/vincentosegura/mlmm
+	</help>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/source_library/emma.mlmm.r	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,1274 @@
+emma.kinship <- function(snps, method="additive", use="all") {
+  n0 <- sum(snps==0,na.rm=TRUE)
+  nh <- sum(snps==0.5,na.rm=TRUE)
+  n1 <- sum(snps==1,na.rm=TRUE)
+  nNA <- sum(is.na(snps))
+
+  stopifnot(n0+nh+n1+nNA == length(snps))
+
+  if ( method == "dominant" ) {
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))
+    snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]
+  }
+  else if ( method == "recessive" ) {
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))
+    snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]
+  }
+  else if ( ( method == "additive" ) && ( nh > 0 ) ) {
+    dsnps <- snps
+    rsnps <- snps
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))
+    dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))
+    rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]
+    snps <- rbind(dsnps,rsnps)
+  }
+
+  if ( use == "all" ) {
+    mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))
+    snps[is.na(snps)] <- mafs[is.na(snps)]
+  }
+  else if ( use == "complete.obs" ) {
+    snps <- snps[rowSums(is.na(snps))==0,]
+  }
+
+  n <- ncol(snps)
+  K <- matrix(nrow=n,ncol=n)
+  diag(K) <- 1
+
+  for(i in 2:n) {
+    for(j in 1:(i-1)) {
+      x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])
+      K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))
+      K[j,i] <- K[i,j]
+    }
+  }
+  return(K)
+}
+
+emma.eigen.L <- function(Z,K,complete=TRUE) {
+  if ( is.null(Z) ) {
+    return(emma.eigen.L.wo.Z(K))
+  }
+  else {
+    return(emma.eigen.L.w.Z(Z,K,complete))
+  }
+}
+
+emma.eigen.L.wo.Z <- function(K) {
+  eig <- eigen(K,symmetric=TRUE)
+  return(list(values=eig$values,vectors=eig$vectors))
+}
+
+emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {
+  if ( complete == FALSE ) {
+    vids <- colSums(Z)>0
+    Z <- Z[,vids]
+    K <- K[vids,vids]
+  }
+  eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)
+  return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))
+}
+
+emma.eigen.R <- function(Z,K,X,complete=TRUE) {
+  if ( ncol(X) == 0 ) {
+    return(emma.eigen.L(Z,K))
+  }
+  else if ( is.null(Z) ) {
+    return(emma.eigen.R.wo.Z(K,X))
+  }
+  else {
+    return(emma.eigen.R.w.Z(Z,K,X,complete))
+  }
+}
+
+emma.eigen.R.wo.Z <- function(K, X) {
+  n <- nrow(X)
+  q <- ncol(X)
+  S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)
+  eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)
+  stopifnot(!is.complex(eig$values))
+  return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))
+}
+
+emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {
+  if ( complete == FALSE ) {
+    vids <-  colSums(Z) > 0
+    Z <- Z[,vids]
+    K <- K[vids,vids]
+  }
+  n <- nrow(Z)
+  t <- ncol(Z)
+  q <- ncol(X)
+
+  SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)
+  eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)
+  if ( is.complex(eig$values) ) {
+    eig$values <- Re(eig$values)
+    eig$vectors <- Re(eig$vectors)
+  }
+  qr.X <- qr.Q(qr(X))
+  return(list(values=eig$values[1:(t-q)],
+              vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),
+                complete=TRUE)[,c(1:(t-q),(t+1):n)]))
+}
+
+emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {
+  n <- length(xi)
+  delta <- exp(logdelta)
+  return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) )
+}
+
+emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
+  t <- length(xi.1)
+  delta <- exp(logdelta)
+#  stopifnot(length(lambda) == length(etas.1))
+  return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )
+}
+
+emma.delta.ML.dLL.wo.Z <- function(logdelta, lambda, etas, xi) {
+  n <- length(xi)
+  delta <- exp(logdelta)
+  etasq <- etas*etas
+  ldelta <- lambda+delta
+  return( 0.5*(n*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/(xi+delta))) )
+}
+
+emma.delta.ML.dLL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
+  t <- length(xi.1)
+  delta <- exp(logdelta)
+  etasq <- etas.1*etas.1
+  ldelta <- lambda+delta
+  return( 0.5*(n*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/(xi.1+delta))+(n-t)/delta) ) )
+}
+
+emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) {
+  nq <- length(etas)
+  delta <-  exp(logdelta)
+  return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas*etas/(lambda+delta))))-sum(log(lambda+delta))) )
+}
+
+emma.delta.REML.LL.w.Z <- function(logdelta, lambda, etas.1, n, t, etas.2.sq ) {
+  tq <- length(etas.1)
+  nq <- n - t + tq
+  delta <-  exp(logdelta)
+  return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(lambda+delta))+(n-t)*logdelta)) )
+}
+
+emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) {
+  nq <- length(etas)
+  delta <- exp(logdelta)
+  etasq <- etas*etas
+  ldelta <- lambda+delta
+  return( 0.5*(nq*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) )
+}
+
+emma.delta.REML.dLL.w.Z <- function(logdelta, lambda, etas.1, n, t1, etas.2.sq ) {
+  t <- t1
+  tq <- length(etas.1)
+  nq <- n - t + tq
+  delta <- exp(logdelta)
+  etasq <- etas.1*etas.1
+  ldelta <- lambda+delta
+  return( 0.5*(nq*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/ldelta)+(n-t)/delta)) )
+}
+
+emma.MLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL, eig.R = NULL)
+{
+  n <- length(y)
+  t <- nrow(K)
+  q <- ncol(X)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X) == n)
+
+  if ( det(crossprod(X,X)) == 0 ) {
+    warning("X is singular")
+    return (list(ML=0,delta=0,ve=0,vg=0))
+  }
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.wo.Z(K)
+    }
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.wo.Z(K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE)
+    Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n-q,m)
+    LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Xis)))
+    dLL <- 0.5*delta*(n*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Xis))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.R$values,etas,eig.L$values))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.R$values,etas,eig.L$values))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas, xi=eig.L$values)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.R$values, etas, eig.L$values))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.w.Z(Z,K)
+    }
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.w.Z(Z,K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+    etas.1 <- etas[1:(t-q)]
+    etas.2 <- etas[(t-q+1):(n-q)]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE)
+    Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t-q,m)
+    #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas))
+    dLL <- 0.5*delta*(n*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.R$values, etas.1, eig.L$values, n, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.R$values+maxdelta))/n
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/n
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.MLE.noX <- function(y, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL)
+{
+  n <- length(y)
+  t <- nrow(K)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.wo.Z(K)
+    }
+    etas <- crossprod(eig.L$vectors,y)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n,m)
+    LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Xis)))-colSums(log(Xis)))
+    dLL <- 0.5*delta*(n*colSums(Etasq/(Xis*Xis))/colSums(Etasq/Xis)-colSums(1/Xis))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    #print(dLL)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.L$values,etas,eig.L$values))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.L$values,etas,eig.L$values))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        #if ( ( dLL[i]*dLL[i+1] < 0 ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas=etas, xi=eig.L$values)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.L$values, etas, eig.L$values))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.w.Z(Z,K)
+    }
+    etas <- crossprod(eig.L$vectors,y)
+    etas.1 <- etas[1:t]
+    etas.2 <- etas[(t+1):n]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t,m)
+    #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas))
+    dLL <- 0.5*delta*(n*(colSums(Etasq/(Xis*Xis))+etas.2.sq/(delta*delta))/(colSums(Etasq/Xis)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.L$values, etas.1, eig.L$values, n, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.L$values+maxdelta))/n
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.L$values+maxdelta))+etas.2.sq/maxdelta)/n
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.REMLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL, eig.R = NULL) {
+  n <- length(y)
+  t <- nrow(K)
+  q <- ncol(X)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X) == n)
+
+  if ( det(crossprod(X,X)) == 0 ) {
+    warning("X is singular")
+    return (list(REML=0,delta=0,ve=0,vg=0))
+  }
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.wo.Z(K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n-q,m)
+    LL <- 0.5*((n-q)*(log((n-q)/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas)))
+    dLL <- 0.5*delta*((n-q)*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Lambdas))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.REML.LL.wo.Z(llim,eig.R$values,etas))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.REML.LL.wo.Z(ulim,eig.R$values,etas))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root,eig.R$values, etas))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.w.Z(Z,K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+    etas.1 <- etas[1:(t-q)]
+    etas.2 <- etas[(t-q+1):(n-q)]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t-q,m)
+    dLL <- 0.5*delta*((n-q)*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Lambdas)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.REML.LL.w.Z(llim,eig.R$values,etas.1,n,t,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.REML.LL.w.Z(ulim,eig.R$values,etas.1,n,t,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.REML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, n=n, t1=t, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.REML.LL.w.Z(r$root,eig.R$values, etas.1, n, t, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.R$values+maxdelta))/(n-q)
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/(n-q)
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(REML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.ML.LRT <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) {
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  g <- nrow(ys)
+  n <- ncol(ys)
+  m <- nrow(xs)
+  t <- ncol(xs)
+  q0 <- ncol(X0)
+  q1 <- q0 + 1
+
+  if ( !ponly ) {
+    ML1s <- matrix(nrow=m,ncol=g)
+    ML0s <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+  ML0 <- vector(length=g)
+
+  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X0) == n)
+
+  if ( sum(is.na(ys)) == 0 ) {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    for(i in 1:g) {
+      ML0[i] <- emma.MLE(ys[i,],X0,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R0)$ML
+    }
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          stats[i,] <- rep(NA,g)
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          ML1s[i,] <- rep(NA,g)
+          ML0s[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          ML1s[i,] <- ML1s[i-1,]
+          ML0s[i,] <- ML0s[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0[vids,,drop=FALSE],xs[i,vids])
+          eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X)
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          nr <- sum(vrows)
+          X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids]%*%t(xs[i,vids,drop=FALSE]))
+          eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X)
+        }
+
+        for(j in 1:g) {
+          if ( nv == t ) {
+            MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+#            MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+            if (!ponly) {
+              ML1s[i,j] <- MLE$ML
+              vgs[i,j] <- MLE$vg
+              ves[i,j] <- MLE$ve
+            }
+            stats[i,j] <- 2*(MLE$ML-ML0[j])
+
+          }
+          else {
+            if ( is.null(Z) ) {
+              eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids])
+              MLE0 <- emma.MLE(ys[j,vids],X0[vids,,drop=FALSE],K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vids],X,K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0)
+            }
+            else {
+              if ( nr == n ) {
+                MLE1 <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids],K[vids,vids])
+                MLE0 <- emma.MLE(ys[j,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0)
+                MLE1 <- emma.MLE(ys[j,vrows],X,K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0)
+              }
+            }
+            if (!ponly) {
+              ML1s[i,j] <- MLE1$ML
+              ML0s[i,j] <- MLE0$ML
+              vgs[i,j] <- MLE1$vg
+              ves[i,j] <- MLE1$ve
+            }
+            stats[i,j] <- 2*(MLE1$ML-MLE0$ML)
+          }
+        }
+        if ( ( nv == t ) && ( !ponly ) ) {
+          ML0s[i,] <- ML0
+        }
+        ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE)
+      }
+    }
+  }
+  else {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    for(i in 1:g) {
+      vrows <- !is.na(ys[i,])
+      if ( is.null(Z) ) {
+        ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)$ML
+      }
+      else {
+        vids <- colSums(Z[vrows,]>0)
+
+        ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp)$ML
+      }
+    }
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          stats[i,] <- rep(NA,g)
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          ML1s[i,] <- rep(NA,g)
+          ML0s[,i] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          ML1s[i,] <- ML1s[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+#          print(j)
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                MLE <- emma.MLE(ys[j,],X,K,NULL,ngrids,llim,ulim,esp,eig.L,eig.R1)
+              }
+              else {
+                MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)
+              }
+            }
+            else {
+              if ( nr == n ) {
+                MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vtids,vtids],Z[vrows,vtids],ngrids,llim,ulim,esp)
+              }
+            }
+
+            if (!ponly) {
+              ML1s[i,j] <- MLE$ML
+              vgs[i,j] <- MLE$vg
+              ves[i,j] <- MLE$ve
+            }
+            stats[i,j] <- 2*(MLE$ML-ML0[j])
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L(NULL,K[vtids,vtids])
+              MLE0 <- emma.MLE(ys[j,vtids],X0[vtids,,drop=FALSE],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vtids],X[vtids,],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0)
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids]))
+              eig.L0 <- emma.eigen.L(Z[vtrows,vtids],K[vtids,vtids])
+              MLE0 <- emma.MLE(ys[j,vtrows],X0[vtrows,,drop=FALSE],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vtrows],X[vtrows,],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0)
+            }
+            if (!ponly) {
+              ML1s[i,j] <- MLE1$ML
+              vgs[i,j] <- MLE1$vg
+              ves[i,j] <- MLE1$ve
+              ML0s[i,j] <- MLE0$ML
+            }
+            stats[i,j] <- 2*(MLE1$ML-MLE0$ML)
+          }
+        }
+        if ( ( nv == t ) && ( !ponly ) ) {
+          ML0s[i,] <- ML0
+        }
+        ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,ML1s=ML1s,ML0s=ML0s,stats=stats,vgs=vgs,ves=ves))
+  }
+}
+
+emma.test <- function(ys, xs, K, Z=NULL, x0s = NULL, X0 = NULL, dfxs = 1, dfx0s = 1, use.MLE = FALSE, use.LRT = FALSE, ngrids = 100, llim = -10, ulim = 10, esp=1e-10, ponly = FALSE)
+{
+  stopifnot (dfxs > 0)
+
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  nx <- nrow(xs)/dfxs
+
+  if ( is.null(x0s) ) {
+    dfx0s = 0
+    x0s <- matrix(NA,0,ncol(xs))
+  }
+  # X0 automatically contains intercept. If no intercept is to be used,
+  #    X0 should be matrix(nrow=ncol(ys),ncol=0)
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  stopifnot(Z == NULL) # The case where Z is not null is not implemented
+
+  ny <- nrow(ys)
+  iy <- ncol(ys)
+  ix <- ncol(xs)
+
+  stopifnot(nrow(K) == ix)
+  stopifnot(ncol(K) == ix)
+  stopifnot(nrow(X0) == iy)
+
+  if ( !ponly ) {
+    LLs <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  dfs <- matrix(nrow=m,ncol=g)
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+
+  # The case with no missing phenotypes
+  if ( sum(is.na(ys)) == 0 ) {
+    if ( ( use.MLE ) || ( !use.LRT ) ) {
+      eig.L0 <- emma.eigen.L(Z,K)
+    }
+    if ( dfx0s == 0 ) {
+      eig.R0 <- emma.eigen.R(Z,K,X0)
+    }
+    x.prev <- NULL
+
+    for(i in 1:ix) {
+      x1 <- t(xs[(dfxs*(i-1)+1):(dfxs*i),,drop=FALSE])
+      if ( dfxs0 == 0 ) {
+        x0 <- X0
+      }
+      else {
+        x0 <- cbind(t(x0s[(dfx0s*(i-1)+1):(dfx0s*i),,drop=FALSE]),X0)
+      }
+      x <-  cbind(x1,x0)
+      xvids <- rowSums(is.na(x) == 0)
+      nxv <- sum(xvids)
+      xv <- x[xvids,,drop=FALSE]
+      Kv <- K[xvids,xvids,drop=FALSE]
+      yv <- ys[j,xvids]
+
+      if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          dfs[i,] <- dfs[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          stats[i,] <- stats[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        eig.R1 = emma.eigen.R.wo.Z(Kv,xv)
+
+        for(j in 1:iy) {
+          if ( ( use.MLE ) || ( !use.LRT ) ) {
+            if ( nxv < t ) {
+              # NOTE: this complexity can be improved by avoiding eigen computation for identical missing patterns
+              eig.L0v <- emma.eigen.L.wo.Z(Kv)
+            }
+            else {
+              eig.L0v <- eig.L0
+            }
+          }
+
+          if ( use.MLE ) {
+            MLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1)
+            stop("Not implemented yet")
+          }
+          else {
+            REMLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1)
+            if ( use.LRT ) {
+              stop("Not implemented yet")
+            }
+            else {
+              U <- eig.L0v$vectors * matrix(sqrt(1/(eig.L0v$values+REMLE$delta)),t,t,byrow=TRUE)
+              dfs[i,j] <- length(eig.R1$values)
+              yt <- crossprod(U,yv)
+              xt <- crossprod(U,xv)
+              ixx <- solve(crossprod(xt,xt))
+              beta <- ixx%*%crossprod(xt,yt)
+              if ( dfxs == 1 ) {
+                stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+              }
+              else {
+                model.m <- c(rep(1,dfxs),rep(0,ncol(xv)-dfxs))
+                stats[i,j] <-
+                  crossprod(crossprod(solve(crossprod(crossprod(iXX,model.m),
+                                                      model.m)),
+                                      model.m*beta),model.m*beta)
+
+              }
+              if ( !ponly ) {
+                vgs[i,j] <- REMLE$vg
+                ves[i,j] <- REMLE$ve
+                REMLs[i,j] <- REMLE$REML
+              }
+            }
+          }
+        }
+        if ( dfxs == 1 ) {
+          ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+        }
+        else {
+          ps[i,] <- pf(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+        }
+      }
+    }
+  }
+  # The case with missing genotypes - not implemented yet
+  else {
+    stop("Not implemented yet")
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          dfs[i,] <- dfs[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids,drop=FALSE]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            yv <- ys[j,vrows]
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X[vrows,,drop=FALSE])
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtids]
+              nr <- sum(vtids)
+              REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtids,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))
+              eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtrows]
+              nr <- sum(vtrows)
+              REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtrows,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            yt <- crossprod(U,yv)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))
+  }
+}
+
+emma.REML.t <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) {
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  g <- nrow(ys)
+  n <- ncol(ys)
+  m <- nrow(xs)
+  t <- ncol(xs)
+  q0 <- ncol(X0)
+  q1 <- q0 + 1
+
+  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X0) == n)
+
+  if ( !ponly ) {
+    REMLs <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  dfs <- matrix(nrow=m,ncol=g)
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+
+  if ( sum(is.na(ys)) == 0 ) {
+    eig.L <- emma.eigen.L(Z,K)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if ( !ponly ) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          stats[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          dfs[i,] <- dfs[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          stats[i,] <- stats[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0[vids,,drop=FALSE],xs[i,vids])
+          eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X)
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X)
+        }
+
+        for(j in 1:g) {
+          if ( nv == t ) {
+            REMLE <- emma.REMLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+            if ( is.null(Z) ) {
+              U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),t,t,byrow=TRUE)
+              dfs[i,j] <- nv - q1
+            }
+            else {
+              U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              dfs[i,j] <- n - q1
+            }
+            yt <- crossprod(U,ys[j,])
+            Xt <- crossprod(U,X)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids])
+              nr <- sum(vids)
+              yv <- ys[j,vids]
+              REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],NULL,ngrids,llim,ulim,esp,eig.R1)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              dfs[i,j] <- nr - q1
+            }
+            else {
+              eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids,drop=FALSE],K[vids,vids])
+              yv <- ys[j,vrows]
+              nr <- sum(vrows)
+              tv <- sum(vids)
+              REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],Z[vrows,vids,drop=FALSE],ngrids,llim,ulim,esp,eig.R1)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-tv)),nr,nr,byrow=TRUE)
+              dfs[i,j] <- nr - q1
+            }
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if (!ponly) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  else {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          dfs[i,] <- dfs[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids,drop=FALSE]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            yv <- ys[j,vrows]
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X[vrows,,drop=FALSE])
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtids]
+              nr <- sum(vtids)
+              REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtids,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))
+              eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtrows]
+              nr <- sum(vtrows)
+              REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtrows,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            yt <- crossprod(U,yv)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/source_library/emma.r	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,1274 @@
+emma.kinship <- function(snps, method="additive", use="all") {
+  n0 <- sum(snps==0,na.rm=TRUE)
+  nh <- sum(snps==0.5,na.rm=TRUE)
+  n1 <- sum(snps==1,na.rm=TRUE)
+  nNA <- sum(is.na(snps))
+
+  stopifnot(n0+nh+n1+nNA == length(snps))
+
+  if ( method == "dominant" ) {
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))
+    snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]
+  }
+  else if ( method == "recessive" ) {
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))
+    snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]
+  }
+  else if ( ( method == "additive" ) && ( nh > 0 ) ) {
+    dsnps <- snps
+    rsnps <- snps
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))
+    dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]
+    flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))
+    rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]
+    snps <- rbind(dsnps,rsnps)
+  }
+
+  if ( use == "all" ) {
+    mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))
+    snps[is.na(snps)] <- mafs[is.na(snps)]
+  }
+  else if ( use == "complete.obs" ) {
+    snps <- snps[rowSums(is.na(snps))==0,]
+  }
+
+  n <- ncol(snps)
+  K <- matrix(nrow=n,ncol=n)
+  diag(K) <- 1
+
+  for(i in 2:n) {
+    for(j in 1:(i-1)) {
+      x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])
+      K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))
+      K[j,i] <- K[i,j]
+    }
+  }
+  return(K)
+}
+
+emma.eigen.L <- function(Z,K,complete=TRUE) {
+  if ( is.null(Z) ) {
+    return(emma.eigen.L.wo.Z(K))
+  }
+  else {
+    return(emma.eigen.L.w.Z(Z,K,complete))
+  }
+}
+
+emma.eigen.L.wo.Z <- function(K) {
+  eig <- eigen(K,symmetric=TRUE)
+  return(list(values=eig$values,vectors=eig$vectors))
+}
+
+emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {
+  if ( complete == FALSE ) {
+    vids <- colSums(Z)>0
+    Z <- Z[,vids]
+    K <- K[vids,vids]
+  }
+  eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)
+  return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))
+}
+
+emma.eigen.R <- function(Z,K,X,complete=TRUE) {
+  if ( ncol(X) == 0 ) {
+    return(emma.eigen.L(Z,K))
+  }
+  else if ( is.null(Z) ) {
+    return(emma.eigen.R.wo.Z(K,X))
+  }
+  else {
+    return(emma.eigen.R.w.Z(Z,K,X,complete))
+  }
+}
+
+emma.eigen.R.wo.Z <- function(K, X) {
+  n <- nrow(X)
+  q <- ncol(X)
+  S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)
+  eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)
+  stopifnot(!is.complex(eig$values))
+  return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))
+}
+
+emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {
+  if ( complete == FALSE ) {
+    vids <-  colSums(Z) > 0
+    Z <- Z[,vids]
+    K <- K[vids,vids]
+  }
+  n <- nrow(Z)
+  t <- ncol(Z)
+  q <- ncol(X)
+
+  SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)
+  eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)
+  if ( is.complex(eig$values) ) {
+    eig$values <- Re(eig$values)
+    eig$vectors <- Re(eig$vectors)
+  }
+  qr.X <- qr.Q(qr(X))
+  return(list(values=eig$values[1:(t-q)],
+              vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),
+                complete=TRUE)[,c(1:(t-q),(t+1):n)]))
+}
+
+emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {
+  n <- length(xi)
+  delta <- exp(logdelta)
+  return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) )
+}
+
+emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
+  t <- length(xi.1)
+  delta <- exp(logdelta)
+#  stopifnot(length(lambda) == length(etas.1))
+  return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )
+}
+
+emma.delta.ML.dLL.wo.Z <- function(logdelta, lambda, etas, xi) {
+  n <- length(xi)
+  delta <- exp(logdelta)
+  etasq <- etas*etas
+  ldelta <- lambda+delta
+  return( 0.5*(n*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/(xi+delta))) )
+}
+
+emma.delta.ML.dLL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
+  t <- length(xi.1)
+  delta <- exp(logdelta)
+  etasq <- etas.1*etas.1
+  ldelta <- lambda+delta
+  return( 0.5*(n*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/(xi.1+delta))+(n-t)/delta) ) )
+}
+
+emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) {
+  nq <- length(etas)
+  delta <-  exp(logdelta)
+  return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas*etas/(lambda+delta))))-sum(log(lambda+delta))) )
+}
+
+emma.delta.REML.LL.w.Z <- function(logdelta, lambda, etas.1, n, t, etas.2.sq ) {
+  tq <- length(etas.1)
+  nq <- n - t + tq
+  delta <-  exp(logdelta)
+  return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(lambda+delta))+(n-t)*logdelta)) )
+}
+
+emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) {
+  nq <- length(etas)
+  delta <- exp(logdelta)
+  etasq <- etas*etas
+  ldelta <- lambda+delta
+  return( 0.5*(nq*sum(etasq/(ldelta*ldelta))/sum(etasq/ldelta)-sum(1/ldelta)) )
+}
+
+emma.delta.REML.dLL.w.Z <- function(logdelta, lambda, etas.1, n, t1, etas.2.sq ) {
+  t <- t1
+  tq <- length(etas.1)
+  nq <- n - t + tq
+  delta <- exp(logdelta)
+  etasq <- etas.1*etas.1
+  ldelta <- lambda+delta
+  return( 0.5*(nq*(sum(etasq/(ldelta*ldelta))+etas.2.sq/(delta*delta))/(sum(etasq/ldelta)+etas.2.sq/delta)-(sum(1/ldelta)+(n-t)/delta)) )
+}
+
+emma.MLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL, eig.R = NULL)
+{
+  n <- length(y)
+  t <- nrow(K)
+  q <- ncol(X)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X) == n)
+
+  if ( det(crossprod(X,X)) == 0 ) {
+    warning("X is singular")
+    return (list(ML=0,delta=0,ve=0,vg=0))
+  }
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.wo.Z(K)
+    }
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.wo.Z(K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE)
+    Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n-q,m)
+    LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Xis)))
+    dLL <- 0.5*delta*(n*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Xis))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.R$values,etas,eig.L$values))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.R$values,etas,eig.L$values))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas, xi=eig.L$values)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.R$values, etas, eig.L$values))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.w.Z(Z,K)
+    }
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.w.Z(Z,K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+    etas.1 <- etas[1:(t-q)]
+    etas.2 <- etas[(t-q+1):(n-q)]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE)
+    Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t-q,m)
+    #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas))
+    dLL <- 0.5*delta*(n*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.R$values, etas.1, eig.L$values, n, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.R$values+maxdelta))/n
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/n
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.MLE.noX <- function(y, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL)
+{
+  n <- length(y)
+  t <- nrow(K)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.wo.Z(K)
+    }
+    etas <- crossprod(eig.L$vectors,y)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Xis <- matrix(eig.L$values,n,m) + matrix(delta,n,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n,m)
+    LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Xis)))-colSums(log(Xis)))
+    dLL <- 0.5*delta*(n*colSums(Etasq/(Xis*Xis))/colSums(Etasq/Xis)-colSums(1/Xis))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    #print(dLL)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.L$values,etas,eig.L$values))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.L$values,etas,eig.L$values))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        #if ( ( dLL[i]*dLL[i+1] < 0 ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas=etas, xi=eig.L$values)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.L$values, etas, eig.L$values))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.L) ) {
+      eig.L <- emma.eigen.L.w.Z(Z,K)
+    }
+    etas <- crossprod(eig.L$vectors,y)
+    etas.1 <- etas[1:t]
+    etas.2 <- etas[(t+1):n]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Xis <- matrix(eig.L$values,t,m) + matrix(delta,t,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t,m)
+    #LL <- 0.5*(n*(log(n/(2*pi))-1-log(colSums(Etasq/Lambdas)+etas.2.sq/delta))-colSums(log(Xis))+(n-t)*log(deltas))
+    dLL <- 0.5*delta*(n*(colSums(Etasq/(Xis*Xis))+etas.2.sq/(delta*delta))/(colSums(Etasq/Xis)+etas.2.sq/delta)-(colSums(1/Xis)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.L$values,etas.1,eig.L$values,n,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.L$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.L$values, etas.1, eig.L$values, n, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.L$values+maxdelta))/n
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.L$values+maxdelta))+etas.2.sq/maxdelta)/n
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.REMLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
+  esp=1e-10, eig.L = NULL, eig.R = NULL) {
+  n <- length(y)
+  t <- nrow(K)
+  q <- ncol(X)
+
+#  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X) == n)
+
+  if ( det(crossprod(X,X)) == 0 ) {
+    warning("X is singular")
+    return (list(REML=0,delta=0,ve=0,vg=0))
+  }
+
+  if ( is.null(Z) ) {
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.wo.Z(K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,n-q,m) + matrix(delta,n-q,m,byrow=TRUE)
+    Etasq <- matrix(etas*etas,n-q,m)
+    LL <- 0.5*((n-q)*(log((n-q)/(2*pi))-1-log(colSums(Etasq/Lambdas)))-colSums(log(Lambdas)))
+    dLL <- 0.5*delta*((n-q)*colSums(Etasq/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(1/Lambdas))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.REML.LL.wo.Z(llim,eig.R$values,etas))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.REML.LL.wo.Z(ulim,eig.R$values,etas))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas)
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root,eig.R$values, etas))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+  else {
+    if ( is.null(eig.R) ) {
+      eig.R <- emma.eigen.R.w.Z(Z,K,X)
+    }
+    etas <- crossprod(eig.R$vectors,y)
+    etas.1 <- etas[1:(t-q)]
+    etas.2 <- etas[(t-q+1):(n-q)]
+    etas.2.sq <- sum(etas.2*etas.2)
+
+    logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
+    m <- length(logdelta)
+    delta <- exp(logdelta)
+    Lambdas <- matrix(eig.R$values,t-q,m) + matrix(delta,t-q,m,byrow=TRUE)
+    Etasq <- matrix(etas.1*etas.1,t-q,m)
+    dLL <- 0.5*delta*((n-q)*(colSums(Etasq/(Lambdas*Lambdas))+etas.2.sq/(delta*delta))/(colSums(Etasq/Lambdas)+etas.2.sq/delta)-(colSums(1/Lambdas)+(n-t)/delta))
+
+    optlogdelta <- vector(length=0)
+    optLL <- vector(length=0)
+    if ( dLL[1] < esp ) {
+      optlogdelta <- append(optlogdelta, llim)
+      optLL <- append(optLL, emma.delta.REML.LL.w.Z(llim,eig.R$values,etas.1,n,t,etas.2.sq))
+    }
+    if ( dLL[m-1] > 0-esp ) {
+      optlogdelta <- append(optlogdelta, ulim)
+      optLL <- append(optLL, emma.delta.REML.LL.w.Z(ulim,eig.R$values,etas.1,n,t,etas.2.sq))
+    }
+
+    for( i in 1:(m-1) )
+      {
+        if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) )
+        {
+          r <- uniroot(emma.delta.REML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, n=n, t1=t, etas.2.sq = etas.2.sq )
+          optlogdelta <- append(optlogdelta, r$root)
+          optLL <- append(optLL, emma.delta.REML.LL.w.Z(r$root,eig.R$values, etas.1, n, t, etas.2.sq ))
+        }
+      }
+#    optdelta <- exp(optlogdelta)
+  }
+
+  maxdelta <- exp(optlogdelta[which.max(optLL)])
+  maxLL <- max(optLL)
+  if ( is.null(Z) ) {
+    maxva <- sum(etas*etas/(eig.R$values+maxdelta))/(n-q)
+  }
+  else {
+    maxva <- (sum(etas.1*etas.1/(eig.R$values+maxdelta))+etas.2.sq/maxdelta)/(n-q)
+  }
+  maxve <- maxva*maxdelta
+
+  return (list(REML=maxLL,delta=maxdelta,ve=maxve,vg=maxva))
+}
+
+emma.ML.LRT <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) {
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  g <- nrow(ys)
+  n <- ncol(ys)
+  m <- nrow(xs)
+  t <- ncol(xs)
+  q0 <- ncol(X0)
+  q1 <- q0 + 1
+
+  if ( !ponly ) {
+    ML1s <- matrix(nrow=m,ncol=g)
+    ML0s <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+  ML0 <- vector(length=g)
+
+  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X0) == n)
+
+  if ( sum(is.na(ys)) == 0 ) {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    for(i in 1:g) {
+      ML0[i] <- emma.MLE(ys[i,],X0,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R0)$ML
+    }
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          stats[i,] <- rep(NA,g)
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          ML1s[i,] <- rep(NA,g)
+          ML0s[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          ML1s[i,] <- ML1s[i-1,]
+          ML0s[i,] <- ML0s[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0[vids,,drop=FALSE],xs[i,vids])
+          eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X)
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          nr <- sum(vrows)
+          X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids]%*%t(xs[i,vids,drop=FALSE]))
+          eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X)
+        }
+
+        for(j in 1:g) {
+          if ( nv == t ) {
+            MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+#            MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+            if (!ponly) {
+              ML1s[i,j] <- MLE$ML
+              vgs[i,j] <- MLE$vg
+              ves[i,j] <- MLE$ve
+            }
+            stats[i,j] <- 2*(MLE$ML-ML0[j])
+
+          }
+          else {
+            if ( is.null(Z) ) {
+              eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids])
+              MLE0 <- emma.MLE(ys[j,vids],X0[vids,,drop=FALSE],K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vids],X,K[vids,vids],NULL,ngrids,llim,ulim,esp,eig.L0)
+            }
+            else {
+              if ( nr == n ) {
+                MLE1 <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids],K[vids,vids])
+                MLE0 <- emma.MLE(ys[j,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0)
+                MLE1 <- emma.MLE(ys[j,vrows],X,K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp,eig.L0)
+              }
+            }
+            if (!ponly) {
+              ML1s[i,j] <- MLE1$ML
+              ML0s[i,j] <- MLE0$ML
+              vgs[i,j] <- MLE1$vg
+              ves[i,j] <- MLE1$ve
+            }
+            stats[i,j] <- 2*(MLE1$ML-MLE0$ML)
+          }
+        }
+        if ( ( nv == t ) && ( !ponly ) ) {
+          ML0s[i,] <- ML0
+        }
+        ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE)
+      }
+    }
+  }
+  else {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    for(i in 1:g) {
+      vrows <- !is.na(ys[i,])
+      if ( is.null(Z) ) {
+        ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)$ML
+      }
+      else {
+        vids <- colSums(Z[vrows,]>0)
+
+        ML0[i] <- emma.MLE(ys[i,vrows],X0[vrows,,drop=FALSE],K[vids,vids],Z[vrows,vids],ngrids,llim,ulim,esp)$ML
+      }
+    }
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          stats[i,] <- rep(NA,g)
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          ML1s[i,] <- rep(NA,g)
+          ML0s[,i] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          ML1s[i,] <- ML1s[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+#          print(j)
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                MLE <- emma.MLE(ys[j,],X,K,NULL,ngrids,llim,ulim,esp,eig.L,eig.R1)
+              }
+              else {
+                MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vrows,vrows],NULL,ngrids,llim,ulim,esp)
+              }
+            }
+            else {
+              if ( nr == n ) {
+                MLE <- emma.MLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.L,eig.R1)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                MLE <- emma.MLE(ys[j,vrows],X[vrows,],K[vtids,vtids],Z[vrows,vtids],ngrids,llim,ulim,esp)
+              }
+            }
+
+            if (!ponly) {
+              ML1s[i,j] <- MLE$ML
+              vgs[i,j] <- MLE$vg
+              ves[i,j] <- MLE$ve
+            }
+            stats[i,j] <- 2*(MLE$ML-ML0[j])
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L(NULL,K[vtids,vtids])
+              MLE0 <- emma.MLE(ys[j,vtids],X0[vtids,,drop=FALSE],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vtids],X[vtids,],K[vtids,vtids],NULL,ngrids,llim,ulim,esp,eig.L0)
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids]))
+              eig.L0 <- emma.eigen.L(Z[vtrows,vtids],K[vtids,vtids])
+              MLE0 <- emma.MLE(ys[j,vtrows],X0[vtrows,,drop=FALSE],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0)
+              MLE1 <- emma.MLE(ys[j,vtrows],X[vtrows,],K[vtids,vtids],Z[vtrows,vtids],ngrids,llim,ulim,esp,eig.L0)
+            }
+            if (!ponly) {
+              ML1s[i,j] <- MLE1$ML
+              vgs[i,j] <- MLE1$vg
+              ves[i,j] <- MLE1$ve
+              ML0s[i,j] <- MLE0$ML
+            }
+            stats[i,j] <- 2*(MLE1$ML-MLE0$ML)
+          }
+        }
+        if ( ( nv == t ) && ( !ponly ) ) {
+          ML0s[i,] <- ML0
+        }
+        ps[i,] <- pchisq(stats[i,],1,lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,ML1s=ML1s,ML0s=ML0s,stats=stats,vgs=vgs,ves=ves))
+  }
+}
+
+emma.test <- function(ys, xs, K, Z=NULL, x0s = NULL, X0 = NULL, dfxs = 1, dfx0s = 1, use.MLE = FALSE, use.LRT = FALSE, ngrids = 100, llim = -10, ulim = 10, esp=1e-10, ponly = FALSE)
+{
+  stopifnot (dfxs > 0)
+
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  nx <- nrow(xs)/dfxs
+
+  if ( is.null(x0s) ) {
+    dfx0s = 0
+    x0s <- matrix(NA,0,ncol(xs))
+  }
+  # X0 automatically contains intercept. If no intercept is to be used,
+  #    X0 should be matrix(nrow=ncol(ys),ncol=0)
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  stopifnot(Z == NULL) # The case where Z is not null is not implemented
+
+  ny <- nrow(ys)
+  iy <- ncol(ys)
+  ix <- ncol(xs)
+
+  stopifnot(nrow(K) == ix)
+  stopifnot(ncol(K) == ix)
+  stopifnot(nrow(X0) == iy)
+
+  if ( !ponly ) {
+    LLs <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  dfs <- matrix(nrow=m,ncol=g)
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+
+  # The case with no missing phenotypes
+  if ( sum(is.na(ys)) == 0 ) {
+    if ( ( use.MLE ) || ( !use.LRT ) ) {
+      eig.L0 <- emma.eigen.L(Z,K)
+    }
+    if ( dfx0s == 0 ) {
+      eig.R0 <- emma.eigen.R(Z,K,X0)
+    }
+    x.prev <- NULL
+
+    for(i in 1:ix) {
+      x1 <- t(xs[(dfxs*(i-1)+1):(dfxs*i),,drop=FALSE])
+      if ( dfxs0 == 0 ) {
+        x0 <- X0
+      }
+      else {
+        x0 <- cbind(t(x0s[(dfx0s*(i-1)+1):(dfx0s*i),,drop=FALSE]),X0)
+      }
+      x <-  cbind(x1,x0)
+      xvids <- rowSums(is.na(x) == 0)
+      nxv <- sum(xvids)
+      xv <- x[xvids,,drop=FALSE]
+      Kv <- K[xvids,xvids,drop=FALSE]
+      yv <- ys[j,xvids]
+
+      if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          dfs[i,] <- dfs[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          stats[i,] <- stats[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        eig.R1 = emma.eigen.R.wo.Z(Kv,xv)
+
+        for(j in 1:iy) {
+          if ( ( use.MLE ) || ( !use.LRT ) ) {
+            if ( nxv < t ) {
+              # NOTE: this complexity can be improved by avoiding eigen computation for identical missing patterns
+              eig.L0v <- emma.eigen.L.wo.Z(Kv)
+            }
+            else {
+              eig.L0v <- eig.L0
+            }
+          }
+
+          if ( use.MLE ) {
+            MLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1)
+            stop("Not implemented yet")
+          }
+          else {
+            REMLE <- emma.REMLE(yv,xv,Kv,NULL,ngrids,llim,ulim,esp,eig.R1)
+            if ( use.LRT ) {
+              stop("Not implemented yet")
+            }
+            else {
+              U <- eig.L0v$vectors * matrix(sqrt(1/(eig.L0v$values+REMLE$delta)),t,t,byrow=TRUE)
+              dfs[i,j] <- length(eig.R1$values)
+              yt <- crossprod(U,yv)
+              xt <- crossprod(U,xv)
+              ixx <- solve(crossprod(xt,xt))
+              beta <- ixx%*%crossprod(xt,yt)
+              if ( dfxs == 1 ) {
+                stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+              }
+              else {
+                model.m <- c(rep(1,dfxs),rep(0,ncol(xv)-dfxs))
+                stats[i,j] <-
+                  crossprod(crossprod(solve(crossprod(crossprod(iXX,model.m),
+                                                      model.m)),
+                                      model.m*beta),model.m*beta)
+
+              }
+              if ( !ponly ) {
+                vgs[i,j] <- REMLE$vg
+                ves[i,j] <- REMLE$ve
+                REMLs[i,j] <- REMLE$REML
+              }
+            }
+          }
+        }
+        if ( dfxs == 1 ) {
+          ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+        }
+        else {
+          ps[i,] <- pf(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+        }
+      }
+    }
+  }
+  # The case with missing genotypes - not implemented yet
+  else {
+    stop("Not implemented yet")
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          dfs[i,] <- dfs[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids,drop=FALSE]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            yv <- ys[j,vrows]
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X[vrows,,drop=FALSE])
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtids]
+              nr <- sum(vtids)
+              REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtids,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))
+              eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtrows]
+              nr <- sum(vtrows)
+              REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtrows,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            yt <- crossprod(U,yv)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))
+  }
+}
+
+emma.REML.t <- function(ys, xs, K, Z=NULL, X0 = NULL, ngrids=100, llim=-10, ulim=10, esp=1e-10, ponly = FALSE) {
+  if ( is.null(dim(ys)) || ncol(ys) == 1 ) {
+    ys <- matrix(ys,1,length(ys))
+  }
+  if ( is.null(dim(xs)) || ncol(xs) == 1 ) {
+    xs <- matrix(xs,1,length(xs))
+  }
+  if ( is.null(X0) ) {
+    X0 <- matrix(1,ncol(ys),1)
+  }
+
+  g <- nrow(ys)
+  n <- ncol(ys)
+  m <- nrow(xs)
+  t <- ncol(xs)
+  q0 <- ncol(X0)
+  q1 <- q0 + 1
+
+  stopifnot(nrow(K) == t)
+  stopifnot(ncol(K) == t)
+  stopifnot(nrow(X0) == n)
+
+  if ( !ponly ) {
+    REMLs <- matrix(nrow=m,ncol=g)
+    vgs <- matrix(nrow=m,ncol=g)
+    ves <- matrix(nrow=m,ncol=g)
+  }
+  dfs <- matrix(nrow=m,ncol=g)
+  stats <- matrix(nrow=m,ncol=g)
+  ps <- matrix(nrow=m,ncol=g)
+
+  if ( sum(is.na(ys)) == 0 ) {
+    eig.L <- emma.eigen.L(Z,K)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if ( !ponly ) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          stats[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          dfs[i,] <- dfs[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          stats[i,] <- stats[i-1,]
+        }
+        ps[i,] <- ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0[vids,,drop=FALSE],xs[i,vids])
+          eig.R1 = emma.eigen.R.wo.Z(K[vids,vids],X)
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids]))
+          X <- cbind(X0[vrows,,drop=FALSE],Z[vrows,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          eig.R1 = emma.eigen.R.w.Z(Z[vrows,vids],K[vids,vids],X)
+        }
+
+        for(j in 1:g) {
+          if ( nv == t ) {
+            REMLE <- emma.REMLE(ys[j,],X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+            if ( is.null(Z) ) {
+              U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),t,t,byrow=TRUE)
+              dfs[i,j] <- nv - q1
+            }
+            else {
+              U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              dfs[i,j] <- n - q1
+            }
+            yt <- crossprod(U,ys[j,])
+            Xt <- crossprod(U,X)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              eig.L0 <- emma.eigen.L.wo.Z(K[vids,vids])
+              nr <- sum(vids)
+              yv <- ys[j,vids]
+              REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],NULL,ngrids,llim,ulim,esp,eig.R1)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              dfs[i,j] <- nr - q1
+            }
+            else {
+              eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vids,drop=FALSE],K[vids,vids])
+              yv <- ys[j,vrows]
+              nr <- sum(vrows)
+              tv <- sum(vids)
+              REMLE <- emma.REMLE(yv,X,K[vids,vids,drop=FALSE],Z[vrows,vids,drop=FALSE],ngrids,llim,ulim,esp,eig.R1)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-tv)),nr,nr,byrow=TRUE)
+              dfs[i,j] <- nr - q1
+            }
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if (!ponly) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  else {
+    eig.L <- emma.eigen.L(Z,K)
+    eig.R0 <- emma.eigen.R(Z,K,X0)
+
+    x.prev <- vector(length=0)
+
+    for(i in 1:m) {
+      vids <- !is.na(xs[i,])
+      nv <- sum(vids)
+      xv <- xs[i,vids]
+
+      if ( ( mean(xv) <= 0 ) || ( mean(xv) >= 1 ) ) {
+        if (!ponly) {
+          vgs[i,] <- rep(NA,g)
+          ves[i,] <- rep(NA,g)
+          REMLs[i,] <- rep(NA,g)
+          dfs[i,] <- rep(NA,g)
+        }
+        ps[i,] = rep(1,g)
+      }
+      else if ( identical(x.prev, xv) ) {
+        if ( !ponly ) {
+          stats[i,] <- stats[i-1,]
+          vgs[i,] <- vgs[i-1,]
+          ves[i,] <- ves[i-1,]
+          REMLs[i,] <- REMLs[i-1,]
+          dfs[i,] <- dfs[i-1,]
+        }
+        ps[i,] = ps[i-1,]
+      }
+      else {
+        if ( is.null(Z) ) {
+          X <- cbind(X0,xs[i,])
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.wo.Z(K,X)
+          }
+        }
+        else {
+          vrows <- as.logical(rowSums(Z[,vids,drop=FALSE]))
+          X <- cbind(X0,Z[,vids,drop=FALSE]%*%t(xs[i,vids,drop=FALSE]))
+          if ( nv == t ) {
+            eig.R1 = emma.eigen.R.w.Z(Z,K,X)
+          }
+        }
+
+        for(j in 1:g) {
+          vrows <- !is.na(ys[j,])
+          if ( nv == t ) {
+            yv <- ys[j,vrows]
+            nr <- sum(vrows)
+            if ( is.null(Z) ) {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE)
+              }
+              else {
+                eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              if ( nr == n ) {
+                REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)
+                U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE)
+              }
+              else {
+                vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))
+                eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+                REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+                U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              }
+              dfs[i,j] <- nr-q1
+            }
+
+            yt <- crossprod(U,yv)
+            Xt <- crossprod(U,X[vrows,,drop=FALSE])
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+          }
+          else {
+            if ( is.null(Z) ) {
+              vtids <- vrows & vids
+              eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtids]
+              nr <- sum(vtids)
+              REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtids,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            else {
+              vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids
+              vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))
+              eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])
+              yv <- ys[j,vtrows]
+              nr <- sum(vtrows)
+              REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)
+              U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)
+              Xt <- crossprod(U,X[vtrows,,drop=FALSE])
+              dfs[i,j] <- nr-q1
+            }
+            yt <- crossprod(U,yv)
+            iXX <- solve(crossprod(Xt,Xt))
+            beta <- iXX%*%crossprod(Xt,yt)
+            if ( !ponly ) {
+              vgs[i,j] <- REMLE$vg
+              ves[i,j] <- REMLE$ve
+              REMLs[i,j] <- REMLE$REML
+            }
+            stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)
+
+          }
+        }
+        ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE)
+      }
+    }
+  }
+  if ( ponly ) {
+    return (ps)
+  }
+  else {
+    return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/source_library/mlmm.r	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,477 @@
+##############################################################################################################################################
+###MLMM - Multi-Locus Mixed Model
+###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH
+#######
+#
+##note: require EMMA
+#library(emma)
+#source('emma.r')
+#
+##REQUIRED DATA & FORMAT
+#
+#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names
+#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=number of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names
+#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names
+#each of these data being sorted in the same way, according to the individual name
+#
+##FOR PLOTING THE GWAS RESULTS
+#SNP INFORMATION - snp_info: a data frame having at least 3 columns:
+# - 1 named 'SNP', with SNP names (same as colnames(X)),
+# - 1 named 'Chr', with the chromosome number to which belong each SNP
+# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.
+#######
+#
+##FUNCTIONS USE
+#save this file somewhere on your computer and source it!
+#source('path/mlmm.r')
+#
+###FORWARD + BACKWARD ANALYSES
+#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)
+#X,Y,K as described above
+#nbchunks: an integer defining the number of chunks of X to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory
+#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,
+#			however to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.
+#			It's value must be specified as an integer >= 3
+#
+###RESULTS
+#
+##STEPWISE TABLE
+#mygwas$step_table
+#
+##PLOTS
+#
+##PLOTS FORM THE FORWARD TABLE
+#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))
+#
+##RSS PLOT
+#plot_step_RSS(mygwas)
+#
+##GWAS MANHATTAN PLOTS
+#
+#FORWARD STEPS
+#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)
+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#
+#OPTIMAL MODELS
+#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria
+#
+#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#
+##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST
+#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)
+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#chrom is an integer specifying the chromosome on which the region of interest is
+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
+#
+#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#chrom is an integer specifying the chromosome on which the region of interest is
+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
+#
+##QQPLOTS of pvalues
+#qqplot_fwd_GWAS(mygwas,nsteps)
+#nsteps=maximum number of forward steps to be displayed
+#
+#qqplot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'))
+#
+##############################################################################################################################################
+
+mlmm<-function(Y,X,K,nbchunks,maxsteps) {
+
+n<-length(Y)
+m<-ncol(X)
+
+stopifnot(ncol(K) == n)
+stopifnot(nrow(K) == n)
+stopifnot(nrow(X) == n)
+stopifnot(nbchunks >= 2)
+stopifnot(maxsteps >= 3)
+
+#INTERCEPT
+
+Xo<-rep(1,n)
+
+#K MATRIX NORMALISATION
+
+K_norm<-(n-1)/sum((diag(n)-matrix(1,n,n)/n)*K)*K
+rm(K)
+
+#step 0 : NULL MODEL
+cof_fwd<-list()
+cof_fwd[[1]]<-as.matrix(Xo)
+colnames(cof_fwd[[1]])<-'Xo'
+
+mod_fwd<-list()
+mod_fwd[[1]]<-emma.REMLE(Y,cof_fwd[[1]],K_norm)
+
+herit_fwd<-list()
+herit_fwd[[1]]<-mod_fwd[[1]]$vg/(mod_fwd[[1]]$vg+mod_fwd[[1]]$ve)
+
+RSSf<-list()
+RSSf[[1]]<-'NA'
+
+RSS_H0<-list()
+RSS_H0[[1]]<-'NA'
+
+df1<-1
+df2<-list()
+df2[[1]]<-'NA'
+
+Ftest<-list()
+Ftest[[1]]<-'NA'
+
+pval<-list()
+pval[[1]]<-'NA'
+
+fwd_lm<-list()
+
+cat('null model done! pseudo-h=',round(herit_fwd[[1]],3),'\n')
+
+#step 1 : EMMAX
+
+M<-solve(chol(mod_fwd[[1]]$vg*K_norm+mod_fwd[[1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[1]])
+fwd_lm[[1]]<-summary(lm(Y_t~0+cof_fwd_t))
+Res_H0<-fwd_lm[[1]]$residuals
+Q_<-qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[1]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[2]]<-unlist(RSS)
+RSS_H0[[2]]<-sum(Res_H0^2)
+df2[[2]]<-n-df1-ncol(cof_fwd[[1]])
+Ftest[[2]]<-(rep(RSS_H0[[2]],length(RSSf[[2]]))/RSSf[[2]]-1)*df2[[2]]/df1
+pval[[2]]<-pf(Ftest[[2]],df1,df2[[2]],lower.tail=FALSE)
+
+cof_fwd[[2]]<-cbind(cof_fwd[[1]],X[,colnames(X) %in% names(which(RSSf[[2]]==min(RSSf[[2]]))[1])])
+colnames(cof_fwd[[2]])<-c(colnames(cof_fwd[[1]]),names(which(RSSf[[2]]==min(RSSf[[2]]))[1]))
+mod_fwd[[2]]<-emma.REMLE(Y,cof_fwd[[2]],K_norm)
+herit_fwd[[2]]<-mod_fwd[[2]]$vg/(mod_fwd[[2]]$vg+mod_fwd[[2]]$ve)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
+
+cat('step 1 done! pseudo-h=',round(herit_fwd[[2]],3),'\n')
+
+#FORWARD
+
+for (i in 3:(maxsteps)) {
+if (herit_fwd[[i-2]] < 0.01) break else {
+
+M<-solve(chol(mod_fwd[[i-1]]$vg*K_norm+mod_fwd[[i-1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[i-1]])
+fwd_lm[[i-1]]<-summary(lm(Y_t~0+cof_fwd_t))
+Res_H0<-fwd_lm[[i-1]]$residuals
+Q_ <- qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[i-1]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[i]]<-unlist(RSS)
+RSS_H0[[i]]<-sum(Res_H0^2)
+df2[[i]]<-n-df1-ncol(cof_fwd[[i-1]])
+Ftest[[i]]<-(rep(RSS_H0[[i]],length(RSSf[[i]]))/RSSf[[i]]-1)*df2[[i]]/df1
+pval[[i]]<-pf(Ftest[[i]],df1,df2[[i]],lower.tail=FALSE)
+
+cof_fwd[[i]]<-cbind(cof_fwd[[i-1]],X[,colnames(X) %in% names(which(RSSf[[i]]==min(RSSf[[i]]))[1])])
+colnames(cof_fwd[[i]])<-c(colnames(cof_fwd[[i-1]]),names(which(RSSf[[i]]==min(RSSf[[i]]))[1]))
+mod_fwd[[i]]<-emma.REMLE(Y,cof_fwd[[i]],K_norm)
+herit_fwd[[i]]<-mod_fwd[[i]]$vg/(mod_fwd[[i]]$vg+mod_fwd[[i]]$ve)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)}
+cat('step ',i-1,' done! pseudo-h=',round(herit_fwd[[i]],3),'\n')}
+rm(i)
+
+##gls at last forward step
+M<-solve(chol(mod_fwd[[length(mod_fwd)]]$vg*K_norm+mod_fwd[[length(mod_fwd)]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[length(mod_fwd)]])
+fwd_lm[[length(mod_fwd)]]<-summary(lm(Y_t~0+cof_fwd_t))
+
+Res_H0<-fwd_lm[[length(mod_fwd)]]$residuals
+Q_ <- qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[length(mod_fwd)]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[length(mod_fwd)+1]]<-unlist(RSS)
+RSS_H0[[length(mod_fwd)+1]]<-sum(Res_H0^2)
+df2[[length(mod_fwd)+1]]<-n-df1-ncol(cof_fwd[[length(mod_fwd)]])
+Ftest[[length(mod_fwd)+1]]<-(rep(RSS_H0[[length(mod_fwd)+1]],length(RSSf[[length(mod_fwd)+1]]))/RSSf[[length(mod_fwd)+1]]-1)*df2[[length(mod_fwd)+1]]/df1
+pval[[length(mod_fwd)+1]]<-pf(Ftest[[length(mod_fwd)+1]],df1,df2[[length(mod_fwd)+1]],lower.tail=FALSE)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
+
+##get max pval at each forward step
+max_pval_fwd<-vector(mode="numeric",length=length(fwd_lm))
+max_pval_fwd[1]<-0
+for (i in 2:length(fwd_lm)) {max_pval_fwd[i]<-max(fwd_lm[[i]]$coef[2:i,4])}
+rm(i)
+
+##get the number of parameters & Loglikelihood from ML at each step
+mod_fwd_LL<-list()
+mod_fwd_LL[[1]]<-list(nfixed=ncol(cof_fwd[[1]]),LL=emma.MLE(Y,cof_fwd[[1]],K_norm)$ML)
+for (i in 2:length(cof_fwd)) {mod_fwd_LL[[i]]<-list(nfixed=ncol(cof_fwd[[i]]),LL=emma.MLE(Y,cof_fwd[[i]],K_norm)$ML)}
+rm(i)
+
+cat('backward analysis','\n')
+
+##BACKWARD (1st step == last fwd step)
+
+dropcof_bwd<-list()
+cof_bwd<-list()
+mod_bwd <- list()
+bwd_lm<-list()
+herit_bwd<-list()
+
+dropcof_bwd[[1]]<-'NA'
+cof_bwd[[1]]<-as.matrix(cof_fwd[[length(mod_fwd)]][,!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]])
+colnames(cof_bwd[[1]])<-colnames(cof_fwd[[length(mod_fwd)]])[!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]]
+mod_bwd[[1]]<-emma.REMLE(Y,cof_bwd[[1]],K_norm)
+herit_bwd[[1]]<-mod_bwd[[1]]$vg/(mod_bwd[[1]]$vg+mod_bwd[[1]]$ve)
+M<-solve(chol(mod_bwd[[1]]$vg*K_norm+mod_bwd[[1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_bwd_t<-crossprod(M,cof_bwd[[1]])
+bwd_lm[[1]]<-summary(lm(Y_t~0+cof_bwd_t))
+
+rm(M,Y_t,cof_bwd_t)
+
+for (i in 2:length(mod_fwd)) {
+dropcof_bwd[[i]]<-(colnames(cof_bwd[[i-1]])[2:ncol(cof_bwd[[i-1]])])[which(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])==min(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])))]
+cof_bwd[[i]]<-as.matrix(cof_bwd[[i-1]][,!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]])
+colnames(cof_bwd[[i]])<-colnames(cof_bwd[[i-1]])[!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]]
+mod_bwd[[i]]<-emma.REMLE(Y,cof_bwd[[i]],K_norm)
+herit_bwd[[i]]<-mod_bwd[[i]]$vg/(mod_bwd[[i]]$vg+mod_bwd[[i]]$ve)
+M<-solve(chol(mod_bwd[[i]]$vg*K_norm+mod_bwd[[i]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_bwd_t<-crossprod(M,cof_bwd[[i]])
+bwd_lm[[i]]<-summary(lm(Y_t~0+cof_bwd_t))
+rm(M,Y_t,cof_bwd_t)}
+
+rm(i)
+
+##get max pval at each backward step
+max_pval_bwd<-vector(mode="numeric",length=length(bwd_lm))
+for (i in 1:(length(bwd_lm)-1)) {max_pval_bwd[i]<-max(bwd_lm[[i]]$coef[2:(length(bwd_lm)+1-i),4])}
+max_pval_bwd[length(bwd_lm)]<-0
+
+##get the number of parameters & Loglikelihood from ML at each step
+mod_bwd_LL<-list()
+mod_bwd_LL[[1]]<-list(nfixed=ncol(cof_bwd[[1]]),LL=emma.MLE(Y,cof_bwd[[1]],K_norm)$ML)
+for (i in 2:length(cof_bwd)) {mod_bwd_LL[[i]]<-list(nfixed=ncol(cof_bwd[[i]]),LL=emma.MLE(Y,cof_bwd[[i]],K_norm)$ML)}
+rm(i)
+
+cat('creating output','\n')
+
+##Forward Table: Fwd + Bwd Tables
+#Compute parameters for model criteria
+BIC<-function(x){-2*x$LL+(x$nfixed+1)*log(n)}
+extBIC<-function(x){BIC(x)+2*lchoose(m,x$nfixed-1)}
+
+fwd_table<-data.frame(step=ncol(cof_fwd[[1]])-1,step_=paste('fwd',ncol(cof_fwd[[1]])-1,sep=''),cof='NA',ncof=ncol(cof_fwd[[1]])-1,h2=herit_fwd[[1]]
+	,maxpval=max_pval_fwd[1],BIC=BIC(mod_fwd_LL[[1]]),extBIC=extBIC(mod_fwd_LL[[1]]))
+for (i in 2:(length(mod_fwd))) {fwd_table<-rbind(fwd_table,
+	data.frame(step=ncol(cof_fwd[[i]])-1,step_=paste('fwd',ncol(cof_fwd[[i]])-1,sep=''),cof=paste('+',colnames(cof_fwd[[i]])[i],sep=''),ncof=ncol(cof_fwd[[i]])-1,h2=herit_fwd[[i]]
+	,maxpval=max_pval_fwd[i],BIC=BIC(mod_fwd_LL[[i]]),extBIC=extBIC(mod_fwd_LL[[i]])))}
+
+rm(i)
+
+bwd_table<-data.frame(step=length(mod_fwd),step_=paste('bwd',0,sep=''),cof=paste('-',dropcof_bwd[[1]],sep=''),ncof=ncol(cof_bwd[[1]])-1,h2=herit_bwd[[1]]
+	,maxpval=max_pval_bwd[1],BIC=BIC(mod_bwd_LL[[1]]),extBIC=extBIC(mod_bwd_LL[[1]]))
+for (i in 2:(length(mod_bwd))) {bwd_table<-rbind(bwd_table,
+	data.frame(step=length(mod_fwd)+i-1,step_=paste('bwd',i-1,sep=''),cof=paste('-',dropcof_bwd[[i]],sep=''),ncof=ncol(cof_bwd[[i]])-1,h2=herit_bwd[[i]]
+	,maxpval=max_pval_bwd[i],BIC=BIC(mod_bwd_LL[[i]]),extBIC=extBIC(mod_bwd_LL[[i]])))}
+
+rm(i,BIC,extBIC,max_pval_fwd,max_pval_bwd,dropcof_bwd)
+
+fwdbwd_table<-rbind(fwd_table,bwd_table)
+
+#RSS for plot
+mod_fwd_RSS<-vector()
+mod_fwd_RSS[1]<-sum((Y-cof_fwd[[1]]%*%fwd_lm[[1]]$coef[,1])^2)
+for (i in 2:length(mod_fwd)) {mod_fwd_RSS[i]<-sum((Y-cof_fwd[[i]]%*%fwd_lm[[i]]$coef[,1])^2)}
+mod_bwd_RSS<-vector()
+mod_bwd_RSS[1]<-sum((Y-cof_bwd[[1]]%*%bwd_lm[[1]]$coef[,1])^2)
+for (i in 2:length(mod_bwd)) {mod_bwd_RSS[i]<-sum((Y-cof_bwd[[i]]%*%bwd_lm[[i]]$coef[,1])^2)}
+expl_RSS<-c(1-sapply(mod_fwd_RSS,function(x){x/mod_fwd_RSS[1]}),1-sapply(mod_bwd_RSS,function(x){x/mod_bwd_RSS[length(mod_bwd_RSS)]}))
+h2_RSS<-c(unlist(herit_fwd),unlist(herit_bwd))*(1-expl_RSS)
+unexpl_RSS<-1-expl_RSS-h2_RSS
+plot_RSS<-t(apply(cbind(expl_RSS,h2_RSS,unexpl_RSS),1,cumsum))
+
+#GLS pvals at each step
+pval_step<-list()
+pval_step[[1]]<-list(out=data.frame('SNP'=colnames(X),'pval'=pval[[2]]),cof='')
+for (i in 2:(length(mod_fwd))) {pval_step[[i]]<-list(out=rbind(data.frame(SNP=colnames(cof_fwd[[i]])[-1],'pval'=fwd_lm[[i]]$coef[2:i,4]),
+	data.frame(SNP=colnames(X)[-which(colnames(X) %in% colnames(cof_fwd[[i]]))],'pval'=pval[[i+1]])),cof=colnames(cof_fwd[[i]])[-1])}
+
+#GLS pvals for best models according to extBIC and mbonf
+
+opt_extBIC<-fwdbwd_table[which(fwdbwd_table$extBIC==min(fwdbwd_table$extBIC))[1],]
+opt_mbonf<-(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),])[which(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof==max(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof))[1],]
+bestmodel_pvals<-function(model) {if(substr(model$step_,start=0,stop=3)=='fwd') {
+		pval_step[[as.integer(substring(model$step_,first=4))+1]]} else if (substr(model$step_,start=0,stop=3)=='bwd') {
+		cof<-cof_bwd[[as.integer(substring(model$step_,first=4))+1]]
+		mixedmod<-emma.REMLE(Y,cof,K_norm)
+		M<-solve(chol(mixedmod$vg*K_norm+mixedmod$ve*diag(n)))
+		Y_t<-crossprod(M,Y)
+		cof_t<-crossprod(M,cof)
+		GLS_lm<-summary(lm(Y_t~0+cof_t))
+		Res_H0<-GLS_lm$residuals
+		Q_ <- qr.Q(qr(cof_t))
+		RSS<-list()
+		for (j in 1:(nbchunks-1)) {
+		X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+		RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+		rm(X_t)}
+		X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof)-1))])
+		RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+		rm(X_t,j)
+		RSSf<-unlist(RSS)
+		RSS_H0<-sum(Res_H0^2)
+		df2<-n-df1-ncol(cof)
+		Ftest<-(rep(RSS_H0,length(RSSf))/RSSf-1)*df2/df1
+		pval<-pf(Ftest,df1,df2,lower.tail=FALSE)
+		list(out=rbind(data.frame(SNP=colnames(cof)[-1],'pval'=GLS_lm$coef[2:(ncol(cof)),4]),
+		data.frame('SNP'=colnames(X)[-which(colnames(X) %in% colnames(cof))],'pval'=pval)),cof=colnames(cof)[-1])} else {cat('error \n')}}
+opt_extBIC_out<-bestmodel_pvals(opt_extBIC)
+opt_mbonf_out<-bestmodel_pvals(opt_mbonf)
+
+list(step_table=fwdbwd_table,pval_step=pval_step,RSSout=plot_RSS,bonf_thresh=-log10(0.05/m),opt_extBIC=opt_extBIC_out,opt_mbonf=opt_mbonf_out)}
+
+plot_step_table<-function(x,type){
+	if (type=='h2') {plot(x$step_table$step,x$step_table$h2,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='h2')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='maxpval'){plot(x$step_table$step,-log10(x$step_table$maxpval),type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='-log10(max_Pval)')
+	abline(h=x$bonf_thresh,lty=2)
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='BIC'){plot(x$step_table$step,x$step_table$BIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='BIC')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='extBIC'){plot(x$step_table$step,x$step_table$extBIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='EBIC')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+	else {cat('error! \n argument type must be one of h2, maxpval, BIC, extBIC')}}
+
+plot_step_RSS<-function(x){
+	op<-par(mar=c(5, 5, 2, 2))
+	plot(0,0,xlim=c(0,nrow(x$RSSout)-1),ylim=c(0,1),xlab='step',ylab='%var',col=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,3],0,0), col='brown1', border=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,2],0,0), col='forestgreen', border=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,1],0,0), col='dodgerblue4', border=0)
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)
+	par(op)}
+
+plot_GWAS<-function(x) {
+	output_<-x$out[order(x$out$Pos),]
+	output_ok<-output_[order(output_$Chr),]
+
+	maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200)))
+	plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2))
+#	plot_col<-c('blue','darkgreen','red','cyan','purple')
+	size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x
+	a<-rep(maxpos[1],size[1])
+	b<-rep(plot_col[1],size[1])
+		for (i in 2:max(unique(output_ok$Chr))){
+	a<-c(a,rep(maxpos[i],size[i]))
+	b<-c(b,rep(plot_col[i],size[i]))}
+
+	output_ok$xpos<-output_ok$Pos+a
+	output_ok$col<-b
+	output_ok$col[output_ok$SNP %in% x$cof]<-'red'
+
+	d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2
+
+	plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome')
+	axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))
+	abline(h=x$bonf_thresh,lty=3,col='black')}
+
+plot_region<-function(x,chrom,pos1,pos2){
+	region<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)
+	region$col<- if (chrom %% 2 == 0) {'gray60'} else {'gray10'}
+	region$col[which(region$SNP %in% x$cof)]<-'red'
+	plot(region$Pos,-log10(region$pval),type='p',pch=20,main=paste('chromosome',chrom,sep=''),xlab='position (bp)',ylab='-log10(pval)',col=region$col,xlim=c(pos1,pos2))
+	abline(h=x$bonf_thresh,lty=3,col='black')}
+
+
+plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {
+	stopifnot(step<=length(x$pval_step))
+	output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
+	plot_GWAS(output)}
+
+plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {
+	stopifnot(step<=length(x$pval_step))
+	output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
+	plot_region(output,chrom,pos1,pos2)}
+
+
+plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {
+	if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
+		plot_GWAS(output)}
+	else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
+		plot_GWAS(output)}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {
+	if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
+		plot_region(output,chrom,pos1,pos2)}
+	else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
+		plot_region(output,chrom,pos1,pos2)}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+
+qqplot_fwd_GWAS<-function(x,nsteps){
+	stopifnot(nsteps<=length(x$pval_step))
+	e<--log10(ppoints(nrow(x$pval_step[[1]]$out)))
+	ostep<-list()
+	ostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))
+	for (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}
+
+	maxp<-ceiling(max(unlist(ostep)))
+
+	plot(e,ostep[[1]],type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))
+	abline(0,1,col="dark grey")
+
+	for (i in 2:nsteps) {
+	par(new=T)
+	plot(e,ostep[[i]],type='b',pch=20,cex=0.8,col=i,axes='F',xlab='',ylab='',xlim=c(0,max(e)+1),ylim=c(0,maxp))}
+	legend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),'cof',sep=' '))
+}
+
+qqplot_opt_GWAS<-function(x,opt){
+	if (opt=='extBIC') {
+		e<--log10(ppoints(nrow(x$opt_extBIC$out)))
+		o<--log10(sort(x$opt_extBIC$out$pval))
+		maxp<-ceiling(max(o))
+		plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to extBIC'))
+		abline(0,1,col="dark grey")}
+	else if (opt=='mbonf') {
+		e<--log10(ppoints(nrow(x$opt_mbonf$out)))
+		o<--log10(sort(x$opt_mbonf$out$pval))
+		maxp<-ceiling(max(o))
+		plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to mbonf'))
+		abline(0,1,col="dark grey")}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/source_library/mlmm1.r	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,486 @@
+##############################################################################################################################################
+###MLMM - Multi-Locus Mixed Model
+###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH
+#######
+#
+##note: require EMMA
+#library(emma)
+#source('emma.r')
+#
+##REQUIRED DATA & FORMAT
+#
+#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names
+#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=n umber of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names
+#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names
+#each of these data being sorted in the same way, according to the individual name
+#
+##FOR PLOTING THE GWAS RESULTS
+#SNP INFORMATION - snp_info: a data frame having at least 3 columns:
+# - 1 named 'SNP', with SNP names (same as colnames(X)),
+# - 1 named 'Chr', with the chromosome number to which belong each SNP
+# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.
+#######
+#
+##FUNCTIONS USE
+#save this file somewhere on your computer and source it!
+#source('path/mlmm.r')
+#
+###FORWARD + BACKWARD ANALYSES
+#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)
+#X,Y,K as described above
+#nbchunks: an integer defining the number of chunks of to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory
+#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,
+#			however to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.
+#			It's value must be specified as an integer >= 3
+#
+###RESULTS
+#
+##STEPWISE TABLE
+#mygwas$step_table
+#
+##PLOTS
+#
+##PLOTS FORM THE FORWARD TABLE
+#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))
+#
+##RSS PLOT
+#plot_step_RSS(mygwas)
+#
+##GWAS MANHATTAN PLOTS
+#
+#FORWARD STEPS
+#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)
+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#
+#OPTIMAL MODELS
+#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria
+#
+#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#
+##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST
+#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)
+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#chrom is an integer specifying the chromosome on which the region of interest is
+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
+#
+#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)
+#snp_info as described above
+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot
+#chrom is an integer specifying the chromosome on which the region of interest is
+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info
+#
+##QQPLOTS of pvalues
+#qqplot_fwd_GWAS(mygwas,nsteps)
+#nsteps=maximum number of forward steps to be displayed
+#
+#qqplot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'))
+#
+##############################################################################################################################################
+
+mlmm<-function(Y,X,K,nbchunks,maxsteps) {
+
+n<-length(Y)
+m<-ncol(X)
+
+stopifnot(ncol(K) == n)
+stopifnot(nrow(K) == n)
+stopifnot(nrow(X) == n)
+stopifnot(nbchunks >= 2)
+stopifnot(maxsteps >= 3)
+
+#INTERCEPT
+
+Xo<-rep(1,n)
+
+#K MATRIX NORMALISATION
+
+K_norm<-(n-1)/sum((diag(n)-matrix(1,n,n)/n)*K)*K
+rm(K)
+
+#step 0 : NULL MODEL
+cof_fwd<-list()
+cof_fwd[[1]]<-as.matrix(Xo)
+colnames(cof_fwd[[1]])<-'Xo'
+
+mod_fwd<-list()
+mod_fwd[[1]]<-emma.REMLE(Y,cof_fwd[[1]],K_norm)
+
+herit_fwd<-list()
+herit_fwd[[1]]<-mod_fwd[[1]]$vg/(mod_fwd[[1]]$vg+mod_fwd[[1]]$ve)
+
+RSSf<-list()
+RSSf[[1]]<-'NA'
+
+RSS_H0<-list()
+RSS_H0[[1]]<-'NA'
+
+df1<-1
+df2<-list()
+df2[[1]]<-'NA'
+
+Ftest<-list()
+Ftest[[1]]<-'NA'
+
+pval<-list()
+pval[[1]]<-'NA'
+
+fwd_lm<-list()
+
+cat('null model done! pseudo-h=',round(herit_fwd[[1]],3),'\n')
+
+#step 1 : EMMAX
+
+M<-solve(chol(mod_fwd[[1]]$vg*K_norm+mod_fwd[[1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[1]])
+fwd_lm[[1]]<-summary(lm(Y_t~0+cof_fwd_t))
+Res_H0<-fwd_lm[[1]]$residuals
+Q_<-qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[1]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[2]]<-unlist(RSS)
+RSS_H0[[2]]<-sum(Res_H0^2)
+df2[[2]]<-n-df1-ncol(cof_fwd[[1]])
+Ftest[[2]]<-(rep(RSS_H0[[2]],length(RSSf[[2]]))/RSSf[[2]]-1)*df2[[2]]/df1
+pval[[2]]<-pf(Ftest[[2]],df1,df2[[2]],lower.tail=FALSE)
+
+cof_fwd[[2]]<-cbind(cof_fwd[[1]],X[,colnames(X) %in% names(which(RSSf[[2]]==min(RSSf[[2]]))[1])])
+colnames(cof_fwd[[2]])<-c(colnames(cof_fwd[[1]]),names(which(RSSf[[2]]==min(RSSf[[2]]))[1]))
+mod_fwd[[2]]<-emma.REMLE(Y,cof_fwd[[2]],K_norm)
+herit_fwd[[2]]<-mod_fwd[[2]]$vg/(mod_fwd[[2]]$vg+mod_fwd[[2]]$ve)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
+
+cat('step 1 done! pseudo-h=',round(herit_fwd[[2]],3),'\n')
+
+#FORWARD
+
+for (i in 3:(maxsteps)) {
+if (herit_fwd[[i-2]] < 0.01) break else {
+
+M<-solve(chol(mod_fwd[[i-1]]$vg*K_norm+mod_fwd[[i-1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[i-1]])
+fwd_lm[[i-1]]<-summary(lm(Y_t~0+cof_fwd_t))
+Res_H0<-fwd_lm[[i-1]]$residuals
+Q_ <- qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[i-1]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[i-1]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[i]]<-unlist(RSS)
+RSS_H0[[i]]<-sum(Res_H0^2)
+df2[[i]]<-n-df1-ncol(cof_fwd[[i-1]])
+Ftest[[i]]<-(rep(RSS_H0[[i]],length(RSSf[[i]]))/RSSf[[i]]-1)*df2[[i]]/df1
+pval[[i]]<-pf(Ftest[[i]],df1,df2[[i]],lower.tail=FALSE)
+
+cof_fwd[[i]]<-cbind(cof_fwd[[i-1]],X[,colnames(X) %in% names(which(RSSf[[i]]==min(RSSf[[i]]))[1])])
+colnames(cof_fwd[[i]])<-c(colnames(cof_fwd[[i-1]]),names(which(RSSf[[i]]==min(RSSf[[i]]))[1]))
+mod_fwd[[i]]<-emma.REMLE(Y,cof_fwd[[i]],K_norm)
+herit_fwd[[i]]<-mod_fwd[[i]]$vg/(mod_fwd[[i]]$vg+mod_fwd[[i]]$ve)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)}
+cat('step ',i-1,' done! pseudo-h=',round(herit_fwd[[i]],3),'\n')}
+rm(i)
+
+##gls at last forward step
+M<-solve(chol(mod_fwd[[length(mod_fwd)]]$vg*K_norm+mod_fwd[[length(mod_fwd)]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_fwd_t<-crossprod(M,cof_fwd[[length(mod_fwd)]])
+fwd_lm[[length(mod_fwd)]]<-summary(lm(Y_t~0+cof_fwd_t))
+
+Res_H0<-fwd_lm[[length(mod_fwd)]]$residuals
+Q_ <- qr.Q(qr(cof_fwd_t))
+
+RSS<-list()
+for (j in 1:(nbchunks-1)) {
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t)}
+X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof_fwd[[length(mod_fwd)]])])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof_fwd[[length(mod_fwd)]])-1))])
+RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+rm(X_t,j)
+
+RSSf[[length(mod_fwd)+1]]<-unlist(RSS)
+RSS_H0[[length(mod_fwd)+1]]<-sum(Res_H0^2)
+df2[[length(mod_fwd)+1]]<-n-df1-ncol(cof_fwd[[length(mod_fwd)]])
+Ftest[[length(mod_fwd)+1]]<-(rep(RSS_H0[[length(mod_fwd)+1]],length(RSSf[[length(mod_fwd)+1]]))/RSSf[[length(mod_fwd)+1]]-1)*df2[[length(mod_fwd)+1]]/df1
+pval[[length(mod_fwd)+1]]<-pf(Ftest[[length(mod_fwd)+1]],df1,df2[[length(mod_fwd)+1]],lower.tail=FALSE)
+rm(M,Y_t,cof_fwd_t,Res_H0,Q_,RSS)
+
+##get max pval at each forward step
+max_pval_fwd<-vector(mode="numeric",length=length(fwd_lm))
+max_pval_fwd[1]<-0
+for (i in 2:length(fwd_lm)) {max_pval_fwd[i]<-max(fwd_lm[[i]]$coef[2:i,4])}
+rm(i)
+
+##get the number of parameters & Loglikelihood from ML at each step
+mod_fwd_LL<-list()
+mod_fwd_LL[[1]]<-list(nfixed=ncol(cof_fwd[[1]]),LL=emma.MLE(Y,cof_fwd[[1]],K_norm)$ML)
+for (i in 2:length(cof_fwd)) {mod_fwd_LL[[i]]<-list(nfixed=ncol(cof_fwd[[i]]),LL=emma.MLE(Y,cof_fwd[[i]],K_norm)$ML)}
+rm(i)
+
+cat('backward analysis','\n')
+
+##BACKWARD (1st step == last fwd step)
+
+dropcof_bwd<-list()
+cof_bwd<-list()
+mod_bwd <- list()
+bwd_lm<-list()
+herit_bwd<-list()
+
+dropcof_bwd[[1]]<-'NA'
+cof_bwd[[1]]<-as.matrix(cof_fwd[[length(mod_fwd)]][,!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]])
+colnames(cof_bwd[[1]])<-colnames(cof_fwd[[length(mod_fwd)]])[!colnames(cof_fwd[[length(mod_fwd)]]) %in% dropcof_bwd[[1]]]
+mod_bwd[[1]]<-emma.REMLE(Y,cof_bwd[[1]],K_norm)
+herit_bwd[[1]]<-mod_bwd[[1]]$vg/(mod_bwd[[1]]$vg+mod_bwd[[1]]$ve)
+M<-solve(chol(mod_bwd[[1]]$vg*K_norm+mod_bwd[[1]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_bwd_t<-crossprod(M,cof_bwd[[1]])
+bwd_lm[[1]]<-summary(lm(Y_t~0+cof_bwd_t))
+
+rm(M,Y_t,cof_bwd_t)
+
+for (i in 2:length(mod_fwd)) {
+dropcof_bwd[[i]]<-(colnames(cof_bwd[[i-1]])[2:ncol(cof_bwd[[i-1]])])[which(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])==min(abs(bwd_lm[[i-1]]$coef[2:nrow(bwd_lm[[i-1]]$coef),3])))]
+cof_bwd[[i]]<-as.matrix(cof_bwd[[i-1]][,!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]])
+colnames(cof_bwd[[i]])<-colnames(cof_bwd[[i-1]])[!colnames(cof_bwd[[i-1]]) %in% dropcof_bwd[[i]]]
+mod_bwd[[i]]<-emma.REMLE(Y,cof_bwd[[i]],K_norm)
+herit_bwd[[i]]<-mod_bwd[[i]]$vg/(mod_bwd[[i]]$vg+mod_bwd[[i]]$ve)
+M<-solve(chol(mod_bwd[[i]]$vg*K_norm+mod_bwd[[i]]$ve*diag(n)))
+Y_t<-crossprod(M,Y)
+cof_bwd_t<-crossprod(M,cof_bwd[[i]])
+bwd_lm[[i]]<-summary(lm(Y_t~0+cof_bwd_t))
+rm(M,Y_t,cof_bwd_t)}
+
+rm(i)
+
+##get max pval at each backward step
+max_pval_bwd<-vector(mode="numeric",length=length(bwd_lm))
+for (i in 1:(length(bwd_lm)-1)) {max_pval_bwd[i]<-max(bwd_lm[[i]]$coef[2:(length(bwd_lm)+1-i),4])}
+max_pval_bwd[length(bwd_lm)]<-0
+
+##get the number of parameters & Loglikelihood from ML at each step
+mod_bwd_LL<-list()
+mod_bwd_LL[[1]]<-list(nfixed=ncol(cof_bwd[[1]]),LL=emma.MLE(Y,cof_bwd[[1]],K_norm)$ML)
+for (i in 2:length(cof_bwd)) {mod_bwd_LL[[i]]<-list(nfixed=ncol(cof_bwd[[i]]),LL=emma.MLE(Y,cof_bwd[[i]],K_norm)$ML)}
+rm(i)
+
+cat('creating output','\n')
+
+##Forward Table: Fwd + Bwd Tables
+#Compute parameters for model criteria
+BIC<-function(x){-2*x$LL+(x$nfixed+1)*log(n)}
+extBIC<-function(x){BIC(x)+2*lchoose(m,x$nfixed-1)}
+
+fwd_table<-data.frame(step=ncol(cof_fwd[[1]])-1,step_=paste('fwd',ncol(cof_fwd[[1]])-1,sep=''),cof='NA',ncof=ncol(cof_fwd[[1]])-1,h2=herit_fwd[[1]]
+	,maxpval=max_pval_fwd[1],BIC=BIC(mod_fwd_LL[[1]]),extBIC=extBIC(mod_fwd_LL[[1]]))
+for (i in 2:(length(mod_fwd))) {fwd_table<-rbind(fwd_table,
+	data.frame(step=ncol(cof_fwd[[i]])-1,step_=paste('fwd',ncol(cof_fwd[[i]])-1,sep=''),cof=paste('+',colnames(cof_fwd[[i]])[i],sep=''),ncof=ncol(cof_fwd[[i]])-1,h2=herit_fwd[[i]]
+	,maxpval=max_pval_fwd[i],BIC=BIC(mod_fwd_LL[[i]]),extBIC=extBIC(mod_fwd_LL[[i]])))}
+
+rm(i)
+
+bwd_table<-data.frame(step=length(mod_fwd),step_=paste('bwd',0,sep=''),cof=paste('-',dropcof_bwd[[1]],sep=''),ncof=ncol(cof_bwd[[1]])-1,h2=herit_bwd[[1]]
+	,maxpval=max_pval_bwd[1],BIC=BIC(mod_bwd_LL[[1]]),extBIC=extBIC(mod_bwd_LL[[1]]))
+for (i in 2:(length(mod_bwd))) {bwd_table<-rbind(bwd_table,
+	data.frame(step=length(mod_fwd)+i-1,step_=paste('bwd',i-1,sep=''),cof=paste('-',dropcof_bwd[[i]],sep=''),ncof=ncol(cof_bwd[[i]])-1,h2=herit_bwd[[i]]
+	,maxpval=max_pval_bwd[i],BIC=BIC(mod_bwd_LL[[i]]),extBIC=extBIC(mod_bwd_LL[[i]])))}
+
+rm(i,BIC,extBIC,max_pval_fwd,max_pval_bwd,dropcof_bwd)
+
+fwdbwd_table<-rbind(fwd_table,bwd_table)
+
+#RSS for plot
+mod_fwd_RSS<-vector()
+mod_fwd_RSS[1]<-sum((Y-cof_fwd[[1]]%*%fwd_lm[[1]]$coef[,1])^2)
+for (i in 2:length(mod_fwd)) {mod_fwd_RSS[i]<-sum((Y-cof_fwd[[i]]%*%fwd_lm[[i]]$coef[,1])^2)}
+mod_bwd_RSS<-vector()
+mod_bwd_RSS[1]<-sum((Y-cof_bwd[[1]]%*%bwd_lm[[1]]$coef[,1])^2)
+for (i in 2:length(mod_bwd)) {mod_bwd_RSS[i]<-sum((Y-cof_bwd[[i]]%*%bwd_lm[[i]]$coef[,1])^2)}
+expl_RSS<-c(1-sapply(mod_fwd_RSS,function(x){x/mod_fwd_RSS[1]}),1-sapply(mod_bwd_RSS,function(x){x/mod_bwd_RSS[length(mod_bwd_RSS)]}))
+h2_RSS<-c(unlist(herit_fwd),unlist(herit_bwd))*(1-expl_RSS)
+unexpl_RSS<-1-expl_RSS-h2_RSS
+plot_RSS<-t(apply(cbind(expl_RSS,h2_RSS,unexpl_RSS),1,cumsum))
+
+#GLS pvals at each step
+pval_step<-list()
+pval_step[[1]]<-list(out=data.frame('SNP'=colnames(X),'pval'=pval[[2]]),cof='')
+for (i in 2:(length(mod_fwd))) {pval_step[[i]]<-list(out=rbind(data.frame(SNP=colnames(cof_fwd[[i]])[-1],'pval'=fwd_lm[[i]]$coef[2:i,4]),
+	data.frame(SNP=colnames(X)[-which(colnames(X) %in% colnames(cof_fwd[[i]]))],'pval'=pval[[i+1]])),cof=colnames(cof_fwd[[i]])[-1])}
+
+#GLS pvals for best models according to extBIC and mbonf
+
+opt_extBIC<-fwdbwd_table[which(fwdbwd_table$extBIC==min(fwdbwd_table$extBIC))[1],]
+opt_mbonf<-(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),])[which(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof==max(fwdbwd_table[which(fwdbwd_table$maxpval<=0.05/m),]$ncof))[1],]
+bestmodel_pvals<-function(model) {if(substr(model$step_,start=0,stop=3)=='fwd') {
+		pval_step[[as.integer(substring(model$step_,first=4))+1]]} else if (substr(model$step_,start=0,stop=3)=='bwd') {
+		cof<-cof_bwd[[as.integer(substring(model$step_,first=4))+1]]
+		mixedmod<-emma.REMLE(Y,cof,K_norm)
+		M<-solve(chol(mixedmod$vg*K_norm+mixedmod$ve*diag(n)))
+		Y_t<-crossprod(M,Y)
+		cof_t<-crossprod(M,cof)
+		GLS_lm<-summary(lm(Y_t~0+cof_t))
+		Res_H0<-GLS_lm$residuals
+		Q_ <- qr.Q(qr(cof_t))
+		RSS<-list()
+		for (j in 1:(nbchunks-1)) {
+		X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j-1)*round(m/nbchunks)+1):(j*round(m/nbchunks))])
+		RSS[[j]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+		rm(X_t)}
+		X_t<-crossprod(M %*% (diag(n)-tcrossprod(Q_,Q_)),(X[,!colnames(X) %in% colnames(cof)])[,((j)*round(m/nbchunks)+1):(m-(ncol(cof)-1))])
+		RSS[[nbchunks]]<-apply(X_t,2,function(x){sum(lsfit(x,Res_H0,intercept = FALSE)$residuals^2)})
+		rm(X_t,j)
+		RSSf<-unlist(RSS)
+		RSS_H0<-sum(Res_H0^2)
+		df2<-n-df1-ncol(cof)
+		Ftest<-(rep(RSS_H0,length(RSSf))/RSSf-1)*df2/df1
+		pval<-pf(Ftest,df1,df2,lower.tail=FALSE)
+		list(out=rbind(data.frame(SNP=colnames(cof)[-1],'pval'=GLS_lm$coef[2:(ncol(cof)),4]),
+		data.frame('SNP'=colnames(X)[-which(colnames(X) %in% colnames(cof))],'pval'=pval)),cof=colnames(cof)[-1])} else {cat('error \n')}}
+opt_extBIC_out<-bestmodel_pvals(opt_extBIC)
+opt_mbonf_out<-bestmodel_pvals(opt_mbonf)
+
+list(step_table=fwdbwd_table,pval_step=pval_step,RSSout=plot_RSS,bonf_thresh=-log10(0.05/m),opt_extBIC=opt_extBIC_out,opt_mbonf=opt_mbonf_out)}
+
+plot_step_table<-function(x,type){
+	if (type=='h2') {plot(x$step_table$step,x$step_table$h2,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='h2')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='maxpval'){plot(x$step_table$step,-log10(x$step_table$maxpval),type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='-log10(max_Pval)')
+	abline(h=x$bonf_thresh,lty=2)
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='BIC'){plot(x$step_table$step,x$step_table$BIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='BIC')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+    else if (type=='extBIC'){plot(x$step_table$step,x$step_table$extBIC,type='b',lty=2,pch=20,col='darkblue',xlab='step',ylab='EBIC')
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)}
+	else {cat('error! \n argument type must be one of h2, maxpval, BIC, extBIC')}}
+
+plot_step_RSS<-function(x){
+	op<-par(mar=c(5, 5, 2, 2))
+	plot(0,0,xlim=c(0,nrow(x$RSSout)-1),ylim=c(0,1),xlab='step',ylab='%var',col=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,3],0,0), col='brown1', border=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,2],0,0), col='forestgreen', border=0)
+	polygon(c(0:(nrow(x$RSSout)-1),(nrow(x$RSSout)-1),0), c(x$RSSout[,1],0,0), col='dodgerblue4', border=0)
+	abline(v=(nrow(x$step_table)/2-0.5),lty=2)
+	par(op)}
+
+plot_GWAS<-function(x) {
+	output_<-x$out[order(x$out$Pos),]
+	output_ok<-output_[order(output_$Chr),]
+
+	maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200)))
+	plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2))
+#	plot_col<-c('blue','darkgreen','red','cyan','purple')
+	size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x
+	a<-rep(maxpos[1],size[1])
+	b<-rep(plot_col[1],size[1])
+		for (i in 2:max(unique(output_ok$Chr))){
+	a<-c(a,rep(maxpos[i],size[i]))
+	b<-c(b,rep(plot_col[i],size[i]))}
+
+	output_ok$xpos<-output_ok$Pos+a
+	output_ok$col<-b
+	output_ok$col[output_ok$SNP %in% x$cof]<-'red'
+
+	d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2
+
+	plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome')
+	axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))
+	abline(h=x$bonf_thresh,lty=3,col='black')
+
+
+    if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) {
+      text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7)
+      legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
+      } else {
+      legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
+      }
+}
+
+plot_region<-function(x,chrom,pos1,pos2){
+	region<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)
+	region$col<- if (chrom %% 2 == 0) {'gray60'} else {'gray10'}
+	region$col[which(region$SNP %in% x$cof)]<-'red'
+	plot(region$Pos,-log10(region$pval),type='p',pch=20,main=paste('chromosome',chrom,sep=''),xlab='position (bp)',ylab='-log10(pval)',col=region$col,xlim=c(pos1,pos2))
+	abline(h=x$bonf_thresh,lty=3,col='black')}
+
+
+plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {
+	stopifnot(step<=length(x$pval_step))
+	output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
+	plot_GWAS(output)}
+
+plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {
+	stopifnot(step<=length(x$pval_step))
+	output<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by='SNP'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)
+	plot_region(output,chrom,pos1,pos2)}
+
+
+plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {
+	if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
+		plot_GWAS(output)}
+	else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
+		plot_GWAS(output)}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {
+	if (opt=='extBIC') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by='SNP'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)
+		plot_region(output,chrom,pos1,pos2)}
+	else if (opt=='mbonf') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by='SNP'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)
+		plot_region(output,chrom,pos1,pos2)}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+
+qqplot_fwd_GWAS<-function(x,nsteps){
+	stopifnot(nsteps<=length(x$pval_step))
+	e<--log10(ppoints(nrow(x$pval_step[[1]]$out)))
+	ostep<-list()
+	ostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))
+	for (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}
+
+	maxp<-ceiling(max(unlist(ostep)))
+
+	plot(e,ostep[[1]],type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))
+	abline(0,1,col="dark grey")
+
+	for (i in 2:nsteps) {
+	par(new=T)
+	plot(e,ostep[[i]],type='b',pch=20,cex=0.8,col=i,axes='F',xlab='',ylab='',xlim=c(0,max(e)+1),ylim=c(0,maxp))}
+	legend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),'cof',sep=' '))
+}
+
+qqplot_opt_GWAS<-function(x,opt){
+	if (opt=='extBIC') {
+		e<--log10(ppoints(nrow(x$opt_extBIC$out)))
+		o<--log10(sort(x$opt_extBIC$out$pval))
+		maxp<-ceiling(max(o))
+		plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to extBIC'))
+		abline(0,1,col="dark grey")}
+	else if (opt=='mbonf') {
+		e<--log10(ppoints(nrow(x$opt_mbonf$out)))
+		o<--log10(sort(x$opt_mbonf$out$pval))
+		maxp<-ceiling(max(o))
+		plot(e,o,type='b',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste('optimal model according to mbonf'))
+		abline(0,1,col="dark grey")}
+	else {cat('error! \n opt must be extBIC or mbonf')}}
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mlmm/source_library/plot_MLMM.Ago.r	Thu Jul 02 05:42:38 2015 -0400
@@ -0,0 +1,36 @@
+## fonction ###########################################################################
+plot_MLMM.Ago<-function(x) {
+	output1 <-x$out[order(x$out$Pos),]
+	output_ok<-output1[order(output1$Chr),]
+
+	maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200)))
+	plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2))
+#	plot_col<-c('blue','darkgreen','red','cyan','purple')
+	size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x
+	a<-rep(maxpos[1],size[1])
+	b<-rep(plot_col[1],size[1])
+		for (i in 2:max(unique(output_ok$Chr))){
+	a<-c(a,rep(maxpos[i],size[i]))
+	b<-c(b,rep(plot_col[i],size[i]))}
+
+	output_ok$xpos<-output_ok$Pos+a
+	output_ok$col<-b
+	output_ok$col[output_ok$mk=='qtl']<-'cyan'
+	output_ok$col[output_ok$SNP %in% x$cof]<-'red'
+
+	d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2
+
+	plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome')
+	axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))
+  xline(output_ok$xpos[output_ok$mk=='qtl'], col='cyan', lwd=0.1)
+	abline(h=x$bonf_thresh,lty=3,col='black')
+
+
+    if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) {
+      text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7)
+      legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
+      } else {
+      legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
+      }
+
+}