version 3.5

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: rrBLUP
 Title: Ridge regression and other kernels for genomic selection
-Version: 3.4
+Version: 3.5
 Author: Jeffrey Endelman
 Maintainer: Jeffrey Endelman <j.endelman@gmail.com>
 Description: One application is to estimate marker effects by ridge
@@ -11,6 +11,6 @@ Description: One application is to estimate marker effects by ridge
 Suggests: multicore
 License: GPL-3
 LazyLoad: yes
-Packaged: 2012-02-06 18:10:17 UTC; jeffendelman
+Packaged: 2012-02-10 20:05:57 UTC; jeffendelman
 Repository: CRAN
-Date/Publication: 2012-02-06 21:06:52
+Date/Publication: 2012-02-14 07:20:59

MD5

@@ -1,10 +1,10 @@
-1a25d39223133ee17259798e2b9c61a8 *DESCRIPTION
+06639f9b86d7b7d138f356c360e8ed3a *DESCRIPTION
 9983baa53027ec7ecc5c44dc94318515 *NAMESPACE
-c049921a49fc22c3a9864722c8e3426c *NEWS
+c5119304691b321449950ff08d841a98 *NEWS
 3426a5afa2b38e987e53b354c267f488 *R/A.mat.R
-53311049c222a3c2bbe1c6d9f61c9a2d *R/GWA.R
+48dcf02533fd05c528f6c6ea2cd022fb *R/GWA.R
 f9b894abb6a3158ebfabe61486dccd94 *R/kinship.BLUP.R
-e6a37c2566e88f73c91bd18120dd5923 *R/mixed.solve.R
+b02c08fda6a873bb6a221650e7547a96 *R/mixed.solve.R
 b74e4e9e2c643e550a83c5ee4eb7e905 *inst/CITATION
 fbaf1e01452ed907bdb412ff35f3b067 *man/A.mat.Rd
 89ab7c69f6c5dbfca37b5220dd14525b *man/GWA.Rd

NEWS

@@ -1,36 +1,40 @@
+Changes in 3.5:
+
+* Labels now transfer from input to output with mixed.solve.
+
 Changes in 3.4:
 
 * Additional features for A.mat.
 
-Changes in 3.3:
-
-* Improved handling of missing genotype data in A.mat.
-* Removed impute function.  Imputation is now done through A.mat.
-
-Changes in 3.2:
-
-* Parallel computation enabled for GWA, A.mat, and kinship.BLUP.
-* Missing genotypic data permitted for GWA.
-
-Changes in 3.1:
-
-* Fixed error in regression coefficient for adjusted A.mat.
-
-Changes in 3.0:
-
-* Major changes to A.mat.  Read the help manual.
-
-Changes in 2.8:
-
-* Added NEWS file
-* Updated A.mat to calculate both the raw and adjusted UAR models.
-
-Changes in 2.7:
-
-* Updated mixed.solve to handle missing observations.
-* Updated A.mat to handle missing alleles.
-
-Changes in 2.5:
-
-* rrBLUP manuscript is published.  Included CITATION file.
-
+Changes in 3.3:
+
+* Improved handling of missing genotype data in A.mat.
+* Removed impute function.  Imputation is now done through A.mat.
+
+Changes in 3.2:
+
+* Parallel computation enabled for GWA, A.mat, and kinship.BLUP.
+* Missing genotypic data permitted for GWA.
+
+Changes in 3.1:
+
+* Fixed error in regression coefficient for adjusted A.mat.
+
+Changes in 3.0:
+
+* Major changes to A.mat.  Read the help manual.
+
+Changes in 2.8:
+
+* Added NEWS file
+* Updated A.mat to calculate both the raw and adjusted UAR models.
+
+Changes in 2.7:
+
+* Updated mixed.solve to handle missing observations.
+* Updated A.mat to handle missing alleles.
+
+Changes in 2.5:
+
+* rrBLUP manuscript is published.  Included CITATION file.
+

R/GWA.R

@@ -32,8 +32,8 @@ if (is.null(Z)) {Z <- diag(n)}
 stopifnot(nrow(Z)==n)
 stopifnot(ncol(Z)==n.line)
 
-Z <- Z[not.NA,]
-X <- X[not.NA,]
+Z <- as.matrix(Z[not.NA,])
+X <- as.matrix(X[not.NA,])
 n <- length(not.NA)
 y <- matrix(y[not.NA],n,1)
 Hinv <- mixed.solve(y,X=X,Z=Z,K=A.mat(G,n.core=n.core,min.MAF=min.MAF),return.Hinv=TRUE)$Hinv  
@@ -42,7 +42,8 @@ df <- p + 1
 if (length(which(is.na(G))) > 0) {missing=TRUE} else {missing=FALSE}
 
 score.calc <- function(G) {
-  scores <- rep(0,ncol(G))
+  scores <- array(0,ncol(G))
+  rownames(scores) <- colnames(G)
   for (i in 1:ncol(G)) {
     Gi <- G[,i]
   freq <- mean(Gi+1,na.rm=TRUE)/2
@@ -96,6 +97,7 @@ score.calc <- function(G) {
   } #ifelse Xsnp full rank
   } #if/else MAF < minMAF
   } #for i
+
   return(scores)
 } #end score.calc
 

R/mixed.solve.R

@@ -29,14 +29,14 @@ if (!is.null(K)) {
 	stopifnot(ncol(K) == m)
 }
 
-Z <- Z[not.NA,]
-X <- X[not.NA,]
+Z <- as.matrix(Z[not.NA,])
+X <- as.matrix(X[not.NA,])
 n <- length(not.NA)
 y <- matrix(y[not.NA],n,1)
 
 XtX <- crossprod(X, X)
 rank.X <- qr(XtX)$rank
-stopifnot(p == rank.X)
+if (rank.X < p) {stop("X not full rank")}
 XtXinv <- solve(XtX)
 S <- diag(n) - tcrossprod(X%*%XtXinv,X)
 if (n <= m + p) {
@@ -114,35 +114,47 @@ Vu.opt <- sum(omega.sq/(theta + lambda.opt))/df
 Ve.opt <- lambda.opt * Vu.opt
 Hinv <- U %*% (t(U)/(phi+lambda.opt))
 W <- crossprod(X,Hinv%*%X)
-beta <- solve(W,crossprod(X,Hinv%*%y))
+beta <- array(solve(W,crossprod(X,Hinv%*%y)))
+rownames(beta) <- colnames(X)
+
 if (is.null(K)) {
 	KZt <- t(Z)
 } else {
 	KZt <- tcrossprod(K,Z)
 }
 KZt.Hinv <- KZt %*% Hinv
-u <- KZt.Hinv %*% (y - X%*%beta)
+u <- array(KZt.Hinv %*% (y - X%*%beta))
+
+if (is.null(K)) {
+    rownames(u) <- colnames(Z)
+} else {
+    rownames(u) <- rownames(K)
+}
+
 LL = -0.5 * (soln$objective + df + df * log(2 * pi/df))
 if (!SE) {
   if (return.Hinv) {
-    list(Vu = Vu.opt, Ve = Ve.opt, beta = as.vector(beta), u = as.vector(u), LL = LL, Hinv = Hinv)
+    list(Vu = Vu.opt, Ve = Ve.opt, beta = beta, u = u, LL = LL, Hinv = Hinv)
   } else {
-    list(Vu = Vu.opt, Ve = Ve.opt, beta = as.vector(beta), u = as.vector(u), LL = LL)
+    list(Vu = Vu.opt, Ve = Ve.opt, beta = beta, u = u, LL = LL)
   }
 } else {
   Winv <- solve(W)
-  beta.SE <- sqrt(Vu.opt*diag(Winv))
+  beta.SE <- array(sqrt(Vu.opt*diag(Winv)))
+  rownames(beta.SE) <- rownames(beta)
   WW <- tcrossprod(KZt.Hinv,KZt)
   WWW <- KZt.Hinv%*%X
   if (is.null(K)) {
-	u.SE <- sqrt(Vu.opt * (rep(1,m) - diag(WW) + diag(tcrossprod(WWW%*%Winv,WWW))))	
+	u.SE <- array(sqrt(Vu.opt * (rep(1,m) - diag(WW) + diag(tcrossprod(WWW%*%Winv,WWW)))))	
   } else {
-	u.SE <- sqrt(Vu.opt * (diag(K) - diag(WW) + diag(tcrossprod(WWW%*%Winv,WWW))))
+	u.SE <- array(sqrt(Vu.opt * (diag(K) - diag(WW) + diag(tcrossprod(WWW%*%Winv,WWW)))))
   }
+  rownames(u.SE) <- rownames(u)
+
   if (return.Hinv) {
-    list(Vu = Vu.opt, Ve = Ve.opt, beta = as.vector(beta), beta.SE = as.vector(beta.SE), u = as.vector(u), u.SE = as.vector(u.SE), LL = LL, Hinv = Hinv)
+    list(Vu = Vu.opt, Ve = Ve.opt, beta = beta, beta.SE = beta.SE, u = u, u.SE = u.SE, LL = LL, Hinv = Hinv)
   } else {
-    list(Vu = Vu.opt, Ve = Ve.opt, beta = as.vector(beta), beta.SE = as.vector(beta.SE), u = as.vector(u), u.SE = as.vector(u.SE), LL = LL)
+    list(Vu = Vu.opt, Ve = Ve.opt, beta = beta, beta.SE = beta.SE, u = u, u.SE = u.SE, LL = LL)
   }
 }
 }