[Genabel-commits] r2056 - in pkg/MultiABEL: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 29 04:31:03 CEST 2016
Author: shenxia
Date: 2016-05-29 04:31:00 +0200 (Sun, 29 May 2016)
New Revision: 2056
Modified:
pkg/MultiABEL/R/MultiLoad.R
pkg/MultiABEL/R/MultiSummary.R
pkg/MultiABEL/R/Multivariate.R
pkg/MultiABEL/src/MultiSummaryLoop.f90
pkg/MultiABEL/src/MultiSummaryLoopInbred.f90
pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90
Log:
MultiSummary() output gets richer
Modified: pkg/MultiABEL/R/MultiLoad.R
===================================================================
--- pkg/MultiABEL/R/MultiLoad.R 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/R/MultiLoad.R 2016-05-29 02:31:00 UTC (rev 2056)
@@ -162,8 +162,9 @@
cat('\n')
}
R <- cor(Y)
+ vy <- colSums((t(t(Y) - colMeans(Y)))^2)/(nrow(Y) - 1)
rownames(bs) <- snp
- res <- list(gwa = bs, cor.pheno = R, cvg = cvg, snp = snp, nid = nrow(Y), nsnp = nsnp)
+ res <- list(gwa = bs, cor.pheno = R, var.pheno = vy, cvg = cvg, snp = snp, nid = nrow(Y), nsnp = nsnp)
class(res) <- 'multi.loaded'
return(res)
}
Modified: pkg/MultiABEL/R/MultiSummary.R
===================================================================
--- pkg/MultiABEL/R/MultiSummary.R 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/R/MultiSummary.R 2016-05-29 02:31:00 UTC (rev 2056)
@@ -95,6 +95,14 @@
cns[seq(3, m*3, 3)] <- pnames
res3 <- matrix(NA, k, m*3)
colnames(res3) <- cns
+ bnames <- paste0('est.coef.', rownames(x$cor.pheno))
+ snames <- paste0('se.coef.', rownames(x$cor.pheno))
+ pnames <- paste0('p.coef.', rownames(x$cor.pheno))
+ cns[seq(1, m*3, 3)] <- bnames
+ cns[seq(2, m*3, 3)] <- snames
+ cns[seq(3, m*3, 3)] <- pnames
+ resc <- matrix(NA, k, m*3)
+ colnames(resc) <- cns
}
if (type == 'outbred') {
scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
@@ -103,7 +111,7 @@
scan2 <- .Fortran('MultiSummaryLoop', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
f = as.numeric(f), betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno), D = matrix(0, m, m),
sdY = diag(sqrt(x$var.pheno)), invsdY = diag(1/sqrt(x$var.pheno)), sY = sqrt(x$var.pheno),
- b = betamat, s = betamat, pil = numeric(k), PACKAGE = "MultiABEL")
+ b = betamat, s = betamat, pil = numeric(k), coef = betamat, ss = betamat, PACKAGE = "MultiABEL")
res2 <- data.frame(p = pf(scan1$fstat, m, n - m - 1, lower.tail = FALSE),
beta.score = scan2$pil)
res2$se.score <- res2$beta.score/sqrt(qchisq(res2$p, 1, lower.tail = FALSE))
@@ -111,6 +119,9 @@
res3[,seq(2, m*3, 3)] <- scan2$s
res3[,seq(3, m*3, 3)] <- pchisq(scan2$b**2/scan2$s**2, 1, lower.tail = FALSE)
res <- cbind(res, res2, res3)
+ resc[,seq(1, m*3, 3)] <- scan2$coef
+ resc[,seq(2, m*3, 3)] <- sqrt(scan2$ss)
+ resc[,seq(3, m*3, 3)] <- pchisq(scan2$coef**2/scan2$ss, 1, lower.tail = FALSE)
} else if (type == 'inbred') {
scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
tmat = tmat, invR = solve(x$cor.pheno),
@@ -118,7 +129,7 @@
scan2 <- .Fortran('MultiSummaryLoopInbred', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
f = as.numeric(f), betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno), D = matrix(0, m, m),
sdY = diag(sqrt(x$var.pheno)), invsdY = diag(1/sqrt(x$var.pheno)), sY = sqrt(x$var.pheno),
- b = betamat, s = betamat, pil = numeric(k), PACKAGE = "MultiABEL")
+ b = betamat, s = betamat, pil = numeric(k), coef = betamat, ss = betamat, PACKAGE = "MultiABEL")
res2 <- data.frame(p = pf(scan1$fstat, m, n - m - 1, lower.tail = FALSE),
beta.score = scan2$pil)
res2$se.score <- res2$beta.score/sqrt(qchisq(res2$p, 1, lower.tail = FALSE))
@@ -126,6 +137,9 @@
res3[,seq(2, m*3, 3)] <- scan2$s
res3[,seq(3, m*3, 3)] <- pchisq(scan2$b**2/scan2$s**2, 1, lower.tail = FALSE)
res <- cbind(res, res2, res3)
+ resc[,seq(1, m*3, 3)] <- scan2$coef
+ resc[,seq(2, m*3, 3)] <- sqrt(scan2$ss)
+ resc[,seq(3, m*3, 3)] <- pchisq(scan2$coef**2/scan2$ss, 1, lower.tail = FALSE)
} else if (type == 'precise' & !is.null(vars)) {
scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
tmat = tmat, invR = solve(x$cor.pheno),
@@ -133,7 +147,7 @@
scan2 <- .Fortran('MultiSummaryLoopPrecise', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
varg = vars, betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno), D = matrix(0, m, m),
sdY = diag(sqrt(x$var.pheno)), invsdY = diag(1/sqrt(x$var.pheno)), sY = sqrt(x$var.pheno),
- b = betamat, s = betamat, pil = numeric(k), PACKAGE = "MultiABEL")
+ b = betamat, s = betamat, pil = numeric(k), coef = betamat, ss = betamat, PACKAGE = "MultiABEL")
res2 <- data.frame(p = pf(scan1$fstat, m, n - m - 1, lower.tail = FALSE),
beta.score = scan2$pil)
res2$se.score <- res2$beta.score/sqrt(qchisq(res2$p, 1, lower.tail = FALSE))
@@ -141,6 +155,9 @@
res3[,seq(2, m*3, 3)] <- scan2$s
res3[,seq(3, m*3, 3)] <- pchisq(scan2$b**2/scan2$s**2, 1, lower.tail = FALSE)
res <- cbind(res, res2, res3)
+ resc[,seq(1, m*3, 3)] <- scan2$coef
+ resc[,seq(2, m*3, 3)] <- sqrt(scan2$ss)
+ resc[,seq(3, m*3, 3)] <- pchisq(scan2$coef**2/scan2$ss, 1, lower.tail = FALSE)
} else if (type == 'direct') {
scan <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
tmat = tmat, invR = solve(x$cor.pheno),
@@ -150,5 +167,11 @@
stop('Wrong type of analysis!')
}
cat('Done.\n')
- return(res)
+ if (type %in% c('outbred', 'inbred', 'precise')) {
+ colnames(scan2$coef) <- rownames(x$cor.pheno)
+ rownames(scan2$coef) <- rownames(res)
+ return(list(scan = res, coef = resc))
+ } else {
+ return(list(scan = res, coef = NULL))
+ }
}
Modified: pkg/MultiABEL/R/Multivariate.R
===================================================================
--- pkg/MultiABEL/R/Multivariate.R 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/R/Multivariate.R 2016-05-29 02:31:00 UTC (rev 2056)
@@ -22,7 +22,7 @@
#' @author Xia Shen
#'
#' @references
-#' Xia Shen, ..., Gordan Lauc, Jim Wilson, Yurii Aulchenko (2015).
+#' Xia Shen, ..., Jim Wilson, Gordan Lauc, Yurii Aulchenko (2015).
#' Multi-omic-variate analysis identified novel loci associated with
#' compound N-Glycosylation of human Immunoglobulin G. \emph{Submitted}.
#'
@@ -72,12 +72,13 @@
m <- nrow(x$cor.pheno)
n <- rep(as.integer(x$nid), k)
betamat <- x$gwa[,seq(1, 2*m, 2)]
- scan <- .Fortran('MultiSummaryLoopPrecise', k = as.integer(k), m = as.integer(m), nn = n,
- betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno),
- E0 = diag(m), H = diag(m), coef = betamat,
- pil = numeric(k), fstat = numeric(k), vg = x$cvg, PACKAGE = "MultiABEL")
+ scan <- .Fortran('MultiSummaryLoopPrecise', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
+ varg = x$cvg, betamat = betamat, R = x$cor.pheno, invR = solve(x$cor.pheno), D = matrix(0, m, m),
+ sdY = diag(sqrt(x$var.pheno)), invsdY = diag(1/sqrt(x$var.pheno)), sY = sqrt(x$var.pheno),
+ b = betamat, s = betamat, pil = numeric(k), coef = betamat, PACKAGE = "MultiABEL")
cat(' OK\n')
res$Beta.S <- scan$pil
+ scan$fstat <- res$Beta.S/m/(1 - res$Beta.S)*(n - m - 1)
res$P <- pf(scan$fstat, m, n - m - 1, lower.tail = FALSE)
res$SE <- res$Beta.S/sqrt(qchisq(res$P, 1, lower.tail = FALSE))
coef <- scan$coef
Modified: pkg/MultiABEL/src/MultiSummaryLoop.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoop.f90 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/src/MultiSummaryLoop.f90 2016-05-29 02:31:00 UTC (rev 2056)
@@ -3,20 +3,20 @@
!
! Created by Xia Shen on 06/04/16.
-subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
+subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil, coef, ss)
implicit none
integer :: ii, i, j, k, m
- double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, D, A, invA
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, invE0, D, A, invA
double precision, dimension(m - 1, m - 1) :: zz, sdY1, R1
- double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY
+ double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY, tmp2, bb
double precision, dimension(m - 1, 1) :: beta1, alpha, sY1, cor1
double precision, dimension(1, m - 1) :: tbeta1
double precision, dimension(k) :: nn, f, pil
- double precision, dimension(k, m) :: betamat, b, s
+ double precision, dimension(k, m) :: betamat, b, s, coef, ss
double precision, dimension(m) :: lambda
- double precision, dimension(1, 1) :: abDalphabeta0, Vb
- double precision :: tmp, vg, abDalphabeta
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb, bDbeta0
+ double precision :: tmp, vg, abDalphabeta, bDbeta, gg
do j = 1, k
vg = 2*f(j)*(1 - f(j))
do i = 1, m
@@ -76,6 +76,21 @@
b(j,i) = ab(m,1)
s(j,i) = sqrt(Vb(1,1))
end do
+
+ tmp2 = matmul(invR, beta)*vg
+ do i = 1, m
+ coef(j,i) = tmp2(i,1)
+ bb(i,1) = tmp2(i,1)
+ D(i,i) = sY(i,1)**2
+ end do
+ gg = vg*nn(j)
+ bDbeta0 = matmul(matmul(transpose(bb), D), beta)
+ bDbeta = bDbeta0(1,1)*vg
+ invE0 = matmul(invsdY, matmul(invR, invsdY))
+ do i = 1, m
+ ss(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)/nn(j)
+ end do
+
H = matmul(beta, transpose(beta))*nn(j)*vg
HinvE0 = matmul(matmul(H, invsdY), matmul(invR, invsdY))/nn(j)
do i = 1, m
Modified: pkg/MultiABEL/src/MultiSummaryLoopInbred.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoopInbred.f90 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/src/MultiSummaryLoopInbred.f90 2016-05-29 02:31:00 UTC (rev 2056)
@@ -3,20 +3,20 @@
!
! Created by Xia Shen on 07/04/16.
-subroutine MultiSummaryLoopInbred(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
+subroutine MultiSummaryLoopInbred(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil, coef, ss)
implicit none
integer :: ii, i, j, k, m
- double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, D, A, invA
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, invE0, D, A, invA
double precision, dimension(m - 1, m - 1) :: zz, sdY1, R1
- double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY
+ double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY, tmp2, bb
double precision, dimension(m - 1, 1) :: beta1, alpha, sY1, cor1
double precision, dimension(1, m - 1) :: tbeta1
double precision, dimension(k) :: nn, f, pil
- double precision, dimension(k, m) :: betamat, b, s
+ double precision, dimension(k, m) :: betamat, b, s, coef, ss
double precision, dimension(m) :: lambda
- double precision, dimension(1, 1) :: abDalphabeta0, Vb
- double precision :: tmp, vg, abDalphabeta
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb, bDbeta0
+ double precision :: tmp, vg, abDalphabeta, bDbeta, gg
do j = 1, k
vg = f(j)*(1 - f(j))
do i = 1, m
@@ -76,6 +76,21 @@
b(j,i) = ab(m,1)
s(j,i) = sqrt(Vb(1,1))
end do
+
+ tmp2 = matmul(invR, beta)*vg
+ do i = 1, m
+ coef(j,i) = tmp2(i,1)
+ bb(i,1) = tmp2(i,1)
+ D(i,i) = sY(i,1)**2
+ end do
+ gg = vg*nn(j)
+ bDbeta0 = matmul(matmul(transpose(bb), D), beta)
+ bDbeta = bDbeta0(1,1)*vg
+ invE0 = matmul(invsdY, matmul(invR, invsdY))
+ do i = 1, m
+ ss(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)/nn(j)
+ end do
+
H = matmul(beta, transpose(beta))*nn(j)*vg
HinvE0 = matmul(matmul(H, invsdY), matmul(invR, invsdY))/nn(j)
do i = 1, m
Modified: pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90 2016-05-07 13:46:26 UTC (rev 2055)
+++ pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90 2016-05-29 02:31:00 UTC (rev 2056)
@@ -3,20 +3,20 @@
!
! Created by Xia Shen on 07/04/16.
-subroutine MultiSummaryLoopPrecise(k, m, nn, varg, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
+subroutine MultiSummaryLoopPrecise(k, m, nn, varg, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil, coef, ss)
implicit none
integer :: ii, i, j, k, m
- double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, D, A, invA
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, invE0, D, A, invA
double precision, dimension(m - 1, m - 1) :: zz, sdY1, R1
- double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY
+ double precision, dimension(m, 1) :: zg, beta, ab, alphabeta, sY, tmp2, bb
double precision, dimension(m - 1, 1) :: beta1, alpha, sY1, cor1
double precision, dimension(1, m - 1) :: tbeta1
double precision, dimension(k) :: nn, varg, pil
- double precision, dimension(k, m) :: betamat, b, s
+ double precision, dimension(k, m) :: betamat, b, s, coef, ss
double precision, dimension(m) :: lambda
- double precision, dimension(1, 1) :: abDalphabeta0, Vb
- double precision :: tmp, vg, abDalphabeta
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb, bDbeta0
+ double precision :: tmp, vg, abDalphabeta, bDbeta, gg
do j = 1, k
vg = varg(j)
do i = 1, m
@@ -76,6 +76,21 @@
b(j,i) = ab(m,1)
s(j,i) = sqrt(Vb(1,1))
end do
+
+ tmp2 = matmul(invR, beta)*vg
+ do i = 1, m
+ coef(j,i) = tmp2(i,1)
+ bb(i,1) = tmp2(i,1)
+ D(i,i) = sY(i,1)**2
+ end do
+ gg = vg*nn(j)
+ bDbeta0 = matmul(matmul(transpose(bb), D), beta)
+ bDbeta = bDbeta0(1,1)*vg
+ invE0 = matmul(invsdY, matmul(invR, invsdY))
+ do i = 1, m
+ ss(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)/nn(j)
+ end do
+
H = matmul(beta, transpose(beta))*nn(j)*vg
HinvE0 = matmul(matmul(H, invsdY), matmul(invR, invsdY))/nn(j)
do i = 1, m
More information about the Genabel-commits
mailing list