From 59d3603db4fb6c040d599b2915320cb7c6f3676a Mon Sep 17 00:00:00 2001 From: Jeffrey Endelman Date: Fri, 10 Feb 2012 00:00:00 +0000 Subject: [PATCH] version 3.5 --- DESCRIPTION | 6 ++--- MD5 | 8 +++--- NEWS | 68 ++++++++++++++++++++++++++----------------------- R/GWA.R | 8 +++--- R/mixed.solve.R | 36 +++++++++++++++++--------- 5 files changed, 72 insertions(+), 54 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index cc1113f..7dcb9af 100755 --- a/DESCRIPTION +++ b/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 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 diff --git a/MD5 b/MD5 index ebdaf22..102e087 100644 --- a/MD5 +++ b/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 diff --git a/NEWS b/NEWS index ceebee8..4ed1c08 100755 --- a/NEWS +++ b/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. + diff --git a/R/GWA.R b/R/GWA.R index 901bae1..2090e77 100755 --- a/R/GWA.R +++ b/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 diff --git a/R/mixed.solve.R b/R/mixed.solve.R index cdaa0f6..26df3dd 100755 --- a/R/mixed.solve.R +++ b/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) } } }