[Genabel-commits] r2053 - in pkg/MultiABEL: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat May 7 15:24:56 CEST 2016
Author: shenxia
Date: 2016-05-07 15:24:56 +0200 (Sat, 07 May 2016)
New Revision: 2053
Modified:
pkg/MultiABEL/ChangeLog
pkg/MultiABEL/R/MultiSummary.R
pkg/MultiABEL/R/load.summary.R
pkg/MultiABEL/src/MultiSummaryLoop.f90
pkg/MultiABEL/src/MultiSummaryLoopInbred.f90
pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90
Log:
Update to 1.1-5. See ChangeLog for details.
Modified: pkg/MultiABEL/ChangeLog
===================================================================
--- pkg/MultiABEL/ChangeLog 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/ChangeLog 2016-05-07 13:24:56 UTC (rev 2053)
@@ -1,5 +1,22 @@
-2015-02-25 15:52 xia
+2016-05-07 14:16 xia
+ * Updated to version 1.1-5
+ * MultiSummary() updated:
+ option type = 'direct' gives only p-values
+ with very fast calculation;
+ 'outbred', 'inbred', or 'precise' will
+ estimate the conditional genetic effects
+ for each of the traits
+ * load.summary() fixedN bug fixed
+
+2016-04-16 21:41 xia
+
+ * load.summary() updated by Yurii to allow
+ flexible input file column names and fixable
+ sample size by the user
+
+2016-02-25 15:52 xia
+
* Updated to version 1.1-4
* MultiSummary() implemented using T-stats
as option type = 'direct'
@@ -14,7 +31,7 @@
* MultiSummary() coefficients standard errors
provided for outbred HWE populations;
Traits can be partially selected for multi-
- trait scans.
+ trait scans
2015-07-31 15:54 xia
@@ -24,7 +41,7 @@
where the former loads individual-level
data and reports detailed single-trait
summary statistics, and the latter calls
- MultiSummary(type = 'precise') procedure.
+ MultiSummary(type = 'precise') procedure
2015-06-05 02:01 xia
Modified: pkg/MultiABEL/R/MultiSummary.R
===================================================================
--- pkg/MultiABEL/R/MultiSummary.R 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/R/MultiSummary.R 2016-05-07 13:24:56 UTC (rev 2053)
@@ -23,14 +23,12 @@
#' @author Xia Shen
#'
#' @references
-#' Xia Shen, Zheng Ning, Yakov Tsepilov, Masoud Shirali,
-#' Generation Scotland, Blair H. Smith, Lynne J. Hocking, Sandosh Padmanabhan, Caroline Hayward,
-#' David J. Porteous, Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2015).
-#' Simple multi-trait analysis identifies novel loci
-#' associated with growth and obesity measures. \emph{Submitted}.
+#' Xia Shen, Zheng Ning, Yakov Tsepilov, Peter K. Joshi,
+#' James F. Wilson, Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko (2016).
+#' Fast pleiotropic meta-analysis for genetic studies. \emph{Submitted}.
#'
#' @seealso
-#' \code{load.summary}, \code{Multivariate}
+#' \code{load.summary}
#'
#' @examples
#' \dontrun{
@@ -58,7 +56,7 @@
#' @aliases MultiSummary, multi.summary
#' @keywords multivariate, meta-analysis
#'
-`MultiSummary` <- function(x, index = NULL, type = 'outbred', vars = NULL) {
+`MultiSummary` <- function(x, index = NULL, type = 'direct', vars = NULL) {
if (class(x) != 'multi.summary') {
stop('Wrong data type!')
}
@@ -76,42 +74,81 @@
betamat <- x$gwa[,seq(1, 2*m, 2)]
semat <- x$gwa[,seq(2, 2*m, 2)]
tmat <- betamat/semat
+ bad <- which(rowSums(is.na(tmat)) > 0)
+ if (length(bad) > 0) {
+ betamat <- betamat[-bad,]
+ semat <- semat[-bad,]
+ tmat <- tmat[-bad,]
+ x$gwa <- x$gwa[-bad,]
+ k <- nrow(x$gwa)
+ }
dimnames(betamat) <- list(NULL, NULL)
- res <- data.frame(Marker = rownames(x$gwa), Freq = f, N = n, Beta.S = NA, SE = NA, P = NA)
+ res <- data.frame(marker = rownames(x$gwa), freq = f, n = n)
rownames(res) <- rownames(x$gwa)
+ if (type != 'direct') {
+ bnames <- paste0('beta.cond.', rownames(x$cor.pheno))
+ snames <- paste0('se.cond.', rownames(x$cor.pheno))
+ pnames <- paste0('p.cond.', rownames(x$cor.pheno))
+ cns <- rep(NA, m*3)
+ cns[seq(1, m*3, 3)] <- bnames
+ cns[seq(2, m*3, 3)] <- snames
+ cns[seq(3, m*3, 3)] <- pnames
+ res3 <- matrix(NA, k, m*3)
+ colnames(res3) <- cns
+ }
if (type == 'outbred') {
- scan <- .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),
- E0 = diag(m), H = diag(m), invE0 = diag(m), D = diag(m), s = betamat, coef = betamat,
- pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
+ scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
+ tmat = tmat, invR = solve(x$cor.pheno),
+ pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
+ 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")
+ 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))
+ res3[,seq(1, m*3, 3)] <- scan2$b
+ 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)
} else if (type == 'inbred') {
- scan <- .Fortran('MultiSummaryLoopInbred', k = as.integer(k), m = as.integer(m), nn = n,
- f = f, 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), PACKAGE = "MultiABEL")
+ scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
+ tmat = tmat, invR = solve(x$cor.pheno),
+ pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
+ 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")
+ 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))
+ res3[,seq(1, m*3, 3)] <- scan2$b
+ 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)
} else if (type == 'precise' & !is.null(vars)) {
- 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 = vars, PACKAGE = "MultiABEL")
+ scan1 <- .Fortran('MultiSummaryLoopDirect', k = as.integer(k), m = as.integer(m), nn = as.numeric(n),
+ tmat = tmat, invR = solve(x$cor.pheno),
+ pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
+ 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")
+ 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))
+ res3[,seq(1, m*3, 3)] <- scan2$b
+ 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)
} 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),
pil = numeric(k), fstat = numeric(k), PACKAGE = "MultiABEL")
+ res$p <- pf(scan$fstat, m, n - m - 1, lower.tail = FALSE)
}else {
stop('Wrong type of analysis!')
}
cat('Done.\n')
- res$Beta.S <- scan$pil
- 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))
- if (type != 'direct') {
- coef <- matrix(NA, k, m*2)
- coef[,seq(1, 2*m, 2)] <- scan$coef
- coef[,seq(2, 2*m, 2)] <- sqrt(scan$s)
- colnames(coef) <- colnames(x$gwa)[1:(2*m)]
- return(data.frame(res, coef))
- } else {
- return(res)
- }
+ return(res)
}
Modified: pkg/MultiABEL/R/load.summary.R
===================================================================
--- pkg/MultiABEL/R/load.summary.R 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/R/load.summary.R 2016-05-07 13:24:56 UTC (rev 2053)
@@ -20,22 +20,17 @@
#' @param est.var A logical value. If \code{FALSE}, each phenotypic variance is assumed to be known as 1.
#' If \code{TRUE}, each phenotypic variance will be estimated to adjust the summary statistics, so that
#' the corresponding phenoypic variance is 1.
-#' @param type A string gives the type of analysis. Default is \code{"outbred"}, referring to
-#' general outbred populations, following Hardy-Weinberg equilibrium. \code{"inbred"} refers to
-#' inbred populations, where no heterzygotes exists, namely, allele frequency = genotype frequency.
-#' \code{"precise"} refers to precise genotypic variance, especially when the individual-level data
-#' are available, for which the argument \code{vars} has to be given.
-#' @param vars A numeric vector gives the variance of the genotypes at each SNP, e.g. coded as 0, 1 and 2.
-#' Only used when \code{type = "precise"}.
#' @param columnNames A vector with names of columns containing necessary information in the input file;
-#' default values are c('snp','a1','freq','beta','se','n'). The values are case-insensitive.
+#' default values are c('snp','a1','freq','beta','se','n'). The values are case-insensitive. Note: check
+#' your allele definitions for different traits are based on the same strand!
#' @param fixedN sample size to assume across all analyses, when provided, this number will be used
#' (instead of the ones specified in the input files)
#'
#' @return The function returns a list of class \code{multi.summary}, containing two elements: \code{gwa}
-#' (the cleaned data to be processed in multi-trait GWAS) and \code{cor.pheno} (user input or estimated).
+#' (the cleaned data to be processed in multi-trait GWAS), \code{cor.pheno} (user input or estimated), and
+#' \code{var.pheno} (default or estimated).
#'
-#' @author Xia Shen, Yurii Aulchenko
+#' @author Xia Shen, Yurii S. Aulchenko
#'
#' @references
#' Xia Shen, Zheng Ning, Yakov Tsepilov, Masoud Shirali,
@@ -71,8 +66,8 @@
#' @aliases load.summary
#' @keywords multivariate, meta-analysis
#'
-`load.summary` <- function(files, cor.pheno = NULL, indep.snps = NULL, est.var = FALSE, type = 'outbred', vars = NULL,
- columnNames = c ('snp','a1','freq','beta','se','n'), fixedN = NULL ) {
+`load.summary` <- function(files, cor.pheno = NULL, indep.snps = NULL, est.var = FALSE,
+ columnNames = c('snp', 'a1', 'freq', 'beta', 'se', 'n'), fixedN = NULL) {
if (!all(is.character(files))) {
stop('files should be given as strings!')
}
@@ -88,30 +83,30 @@
stop('wrong dimensions of cor.pheno!')
}
}
- columnNames <- tolower( columnNames)
- if (!is.null( fixedN )) if (fixedN <= 0) {
- stop('fixedN should be a positive number')
- }
- if (is.null(fixedN)) { colNamLen = 6 } else { colNamLen = 5 }
- if (!is.character(columnNames)) {
- stop('columnNames should be character')
- }
- if (length( columnNames ) != colNamLen) {
- cat('columnNames should be a vector with',colNamLen,'elements')
- stop('... exiting')
- }
- if ( length(unique(columnNames)) != colNamLen ) {
- stop('elements of columnNames must be unique')
- }
- # column Names Translation
- cNT = list(
- 'snp' = columnNames[1],
- 'a1' = columnNames[2],
- 'freq'= columnNames[3],
- 'beta'= columnNames[4],
- 'se' = columnNames[5],
- 'n' = columnNames[6]
- )
+ columnNames <- tolower(columnNames)
+ if (!is.null(fixedN)) if (fixedN <= 0) {
+ stop('fixedN should be a positive number')
+ }
+ if (is.null(fixedN)) { colNamLen = 6 } else { colNamLen = 5 }
+ if (!is.character(columnNames)) {
+ stop('columnNames should be character')
+ }
+ if (length(columnNames) != colNamLen) {
+ cat('columnNames should be a vector with',colNamLen,'elements')
+ stop('... exiting')
+ }
+ if (length(unique(columnNames)) != colNamLen) {
+ stop('elements of columnNames must be unique')
+ }
+ # column Names Translation
+ cNT = list(
+ 'snp' = columnNames[1],
+ 'a1' = columnNames[2],
+ 'freq'= columnNames[3],
+ 'beta'= columnNames[4],
+ 'se' = columnNames[5],
+ 'n' = columnNames[6]
+ )
cat('loading data ...\n')
data <- c()
fn <- files # rev(files)
@@ -119,27 +114,26 @@
for (i in m:1) {
dd <- read.table(fn[i], header = TRUE, stringsAsFactors = FALSE)
colnames(dd) <- tolower(colnames(dd))
- currentColNames <- colnames(dd)
- if ( any( !( columnNames %in% currentColNames ) ) ) {
- cat('file column names do not match columnNames in ',fn[i],'... ')
- stop('exiting')
- }
- idx <- which(duplicated(dd[, cNT[['snp']] ]))
+ currentColNames <- colnames(dd)
+ if (any(!(columnNames %in% currentColNames))) {
+ stop('file column names do not match columnNames in ', fn[i], '... exiting!')
+ }
+ idx <- which(duplicated(dd[, cNT$snp]))
if (length(idx) > 0) {
data[[i]] <- dd[-idx,]
- rownames(data[[i]]) <- dd[ -idx , cNT[['snp']] ]
+ rownames(data[[i]]) <- dd[ -idx , cNT$snp]
} else {
data[[i]] <- dd
- rownames(data[[i]]) <- dd[, cNT[['snp']] ]
+ rownames(data[[i]]) <- dd[, cNT$snp]
}
if (est.var) {
- if (!is.null(fixedN)) {
- D <- dd[, cNT[['n']] ]*2*dd[, cNT[['freq']] ]*(1 - dd[, cNT[['freq']] ])
- vy <- D*dd[ , cNT[['se']] ]**2 + D*dd[, cNT[['beta']] ]**2/(dd[, cNT[['n']] ] - 1)
- } else {
- D <- fixedN*2*dd[, cNT[['freq']] ]*(1 - dd[, cNT[['freq']] ])
- vy <- D*dd[ , cNT[['se']] ]**2 + D*dd[, cNT[['beta']] ]**2/( fixedN - 1)
- }
+ if (is.null(fixedN)) {
+ D <- dd[,cNT$n]*2*dd[, cNT$freq]*(1 - dd[, cNT$freq])
+ vy <- D*dd[ , cNT$se]**2 + D*dd[, cNT$beta]**2/(dd[, cNT$n] - 1)
+ } else {
+ D <- fixedN*2*dd[, cNT$freq]*(1 - dd[, cNT$freq])
+ vy <- D*dd[ , cNT$se]**2 + D*dd[, cNT$beta]**2/(fixedN - 1)
+ }
dvy <- density(na.omit(vy))
vys[i] <- dvy$x[which.max(dvy$y)] #median(vy, na.rm = TRUE)
}
@@ -148,9 +142,9 @@
cat('\n')
if (est.var) cat('phenotypic variances are:', vys, '\n')
cat('checking markers ...\n')
- snps <- data[[1]][, cNT[['snp']] ]
+ snps <- data[[1]][, cNT$snp]
for (i in 2:m) {
- snps <- data[[i]][ data[[i]][, cNT[['snp']] ] %in% snps, cNT[['snp']] ]
+ snps <- data[[i]][ data[[i]][, cNT$snp] %in% snps, cNT$snp]
progress(i/m*100)
}
snps <- unique(snps)
@@ -163,25 +157,25 @@
cat('\n')
cat('correcting parameters ...\n')
for (i in 2:m) {
- if (any( data[[i]][, cNT[['a1']] ] != data[[1]][, cNT[['a1']] ] )) {
- adj <- 2*as.numeric( data[[i]][, cNT[['a1']] ] == data[[1]][, cNT[['a1']] ] ) - 1
- data[[i]][, cNT$beta ] <- data[[i]][, cNT$beta ]*adj
- data[[i]][, cNT$freq ] <- (adj == 1)*data[[i]][, cNT$freq] + (adj == -1)*(1 - data[[i]][,cNT$freq])
+ if (any(data[[i]][, cNT[['a1']]] != data[[1]][, cNT[['a1']]])) {
+ adj <- 2*as.numeric(data[[i]][, cNT[['a1']]] == data[[1]][, cNT[['a1']]]) - 1
+ data[[i]][, cNT$beta] <- data[[i]][, cNT$beta]*adj
+ data[[i]][, cNT$freq] <- (adj == 1)*data[[i]][, cNT$freq] + (adj == -1)*(1 - data[[i]][,cNT$freq])
}
progress(i/m*100)
}
cat('\n')
cat('adjusting sample size ... ')
n0 <- matrix(NA, nrow(data[[1]]), m)
- if (is.null(fixedN)) {
- for (i in 1:m) {
- n0[,i] <- data[[i]][,cNT$n]
- }
- } else {
- for (i in 1:m) {
- n0[,i] <- fixedN
- }
- }
+ if (is.null(fixedN)) {
+ for (i in 1:m) {
+ n0[,i] <- data[[i]][,cNT$n]
+ }
+ } else {
+ for (i in 1:m) {
+ n0[,i] <- fixedN
+ }
+ }
n <- apply(n0, 1, "min")
cat('done.\n')
cat('finalizing summary statistics ...\n')
@@ -219,7 +213,7 @@
dimnames(cor.pheno) <- list(files, files)
gwanames <- c(paste(rep(files, each = 2), rep(c('.beta', '.se'), m), sep = ''), 'f', 'n')
colnames(gwa0) <- gwanames
- dd <- list(gwa = gwa0, cor.pheno = cor.pheno)
+ dd <- list(gwa = gwa0, cor.pheno = cor.pheno, var.pheno = vys)
class(dd) <- 'multi.summary'
return(dd)
}
Modified: pkg/MultiABEL/src/MultiSummaryLoop.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoop.f90 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/src/MultiSummaryLoop.f90 2016-05-07 13:24:56 UTC (rev 2053)
@@ -1,49 +1,173 @@
!
! MultiSummaryLoop.f90
!
-! Created by Xia Shen on 5/29/15.
-!
+! Created by Xia Shen on 06/04/16.
-subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, E0, H, invE0, D, s, coef, pil, fstat)
+subroutine MultiSummaryLoop(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
implicit none
- integer :: i, j, k, m
- double precision, dimension(m, m) :: R, invR, E0, H, HinvE0, D, invE0
- double precision, dimension(1, m) :: beta, b
- double precision, dimension(m, 1) :: tmp
- double precision, dimension(k) :: nn, f, pil, fstat
- double precision, dimension(k, m) :: betamat, coef, s
+ integer :: ii, i, j, k, m
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, 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, 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(m) :: lambda
- double precision, dimension(1, 1) :: bDbeta0
- double precision :: vg, gg, bDbeta
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb
+ double precision :: tmp, vg, abDalphabeta
do j = 1, k
vg = 2*f(j)*(1 - f(j))
- E0 = nn(j)*R
- invE0 = invR/nn(j)
do i = 1, m
- beta(1,i) = betamat(j,i)
- D(i,i) = E0(i,i)
+ beta(i,1) = betamat(j,i)
end do
- tmp = matmul(invR, transpose(beta))*vg
do i = 1, m
- coef(j,i) = tmp(i,1)
- b(1,i) = tmp(i,1)
+ if (i == 1) then
+ R1 = R(2:m,2:m)
+ sdY1 = sdY(2:m,2:m)
+ beta1(1:(m - 1),1) = beta(2:m,1)
+ sY1(1:(m - 1),1) = sY(2:m,1)
+ cor1(1:(m - 1),1) = R(2:m,1)
+ else if (i == m) then
+ R1 = R(1:(m - 1),1:(m - 1))
+ sdY1 = sdY(1:(m - 1),1:(m - 1))
+ beta1(1:(m - 1),1) = beta(1:(m - 1),1)
+ sY1(1:(m - 1),1) = sY(1:(m - 1),1)
+ cor1(1:(m - 1),1) = R(1:(m - 1),m)
+ else
+ R1(1:(i - 1),1:(i - 1)) = R(1:(i - 1),1:(i - 1))
+ R1(i:(m - 1),i:(m - 1)) = R((i + 1):m,(i + 1):m)
+ R1(1:(i - 1),i:(m - 1)) = R(1:(i - 1),(i + 1):m)
+ R1(i:(m - 1),1:(i - 1)) = R((i + 1):m,1:(i - 1))
+ sdY1(1:(i - 1),1:(i - 1)) = sdY(1:(i - 1),1:(i - 1))
+ sdY1(i:(m - 1),i:(m - 1)) = sdY((i + 1):m,(i + 1):m)
+ sdY1(1:(i - 1),i:(m - 1)) = sdY(1:(i - 1),(i + 1):m)
+ sdY1(i:(m - 1),1:(i - 1)) = sdY((i + 1):m,1:(i - 1))
+ beta1(1:(i - 1),1) = beta(1:(i - 1),1)
+ beta1(i:(m - 1),1) = beta((i + 1):m,1)
+ sY1(1:(i - 1),1) = sY(1:(i - 1),1)
+ sY1(i:(m - 1),1) = sY((i + 1):m,1)
+ cor1(1:(i - 1),1) = R(1:(i - 1),i)
+ cor1(i:(m - 1),1) = R((i + 1):m,i)
+ end if
+ zg = beta1*vg
+ zz = matmul(matmul(sdY1, R1), sdY1)
+ tmp = sdY(i,i)
+ do ii = 1, m - 1
+ alpha(ii,1) = cor1(ii,1)/sY1(ii,1)*tmp
+ end do
+
+ alphabeta(1:(m - 1),1) = alpha(1:(m - 1),1)
+ alphabeta(m,1) = beta(i,1)
+ A(1:(m - 1),1:(m - 1)) = zz(1:(m - 1), 1:(m - 1))
+ A(1:(m - 1),m) = beta1(1:(m - 1), 1)*vg
+ tbeta1 = transpose(beta1)
+ A(m,1:(m - 1)) = tbeta1(1,1:(m - 1))*vg
+ A(m,m) = vg
+ do ii = 1, m
+ D(ii,ii) = A(ii,ii)
+ end do
+ invA = inverse(A, m)
+ ab = matmul(matmul(invA, D), alphabeta)
+ abDalphabeta0 = matmul(matmul(transpose(ab), D), alphabeta)
+ abDalphabeta = abDalphabeta0(1,1)
+ Vb = (sdY(i,i)**2 - abDalphabeta)/(nn(j) - 3)*invA(m, m)
+ b(j,i) = ab(m,1)
+ s(j,i) = sqrt(Vb(1,1))
end do
- gg = vg*nn(j)
- bDbeta0 = matmul(matmul(b, D), transpose(beta))
- bDbeta = bDbeta0(1,1)*vg
+ H = matmul(beta, transpose(beta))*nn(j)*vg
+ HinvE0 = matmul(matmul(H, invsdY), matmul(invR, invsdY))/nn(j)
do i = 1, m
- s(j,i) = (gg - bDbeta)/(nn(j) - m)*invE0(i,i)
- end do
- H = matmul(transpose(beta), beta)*nn(j)*vg
- HinvE0 = matmul(H, invR)/nn(j)
- do i = 1, m
lambda(i) = HinvE0(i,i)
end do
pil(j) = sum(lambda)
- fstat(j) = pil(j)/m/(1 - pil(j))*(nn(j) - m - 1)
end do
+contains
+
+function inverse(a, n)
+!============================================================
+! Inverse matrix
+! Method: Based on Doolittle LU factorization for Ax=b
+! Alex G. December 2009
+!-----------------------------------------------------------
+! input ...
+! a(n,n) - array of coefficients for matrix A
+! n - dimension
+! output ...
+! c(n,n) - inverse matrix of A
+! comments ...
+! the original matrix a(n,n) will be destroyed
+! during the calculation
+!===========================================================
+implicit none
+integer :: n
+double precision :: a(n,n)
+double precision :: c(n,n), inverse(n, n)
+double precision :: L(n,n), U(n,n), b(n), d(n), x(n)
+double precision :: coeff
+integer :: i, j, k
+
+! step 0: initialization for matrices L and U and b
+! Fortran 90/95 aloows such operations on matrices
+L=0.0
+U=0.0
+b=0.0
+
+! step 1: forward elimination
+do k=1, n-1
+ do i=k+1,n
+ coeff=a(i,k)/a(k,k)
+ L(i,k) = coeff
+ do j=k+1,n
+ a(i,j) = a(i,j)-coeff*a(k,j)
+ end do
+ end do
+end do
+
+! Step 2: prepare L and U matrices
+! L matrix is a matrix of the elimination coefficient
+! + the diagonal elements are 1.0
+do i=1,n
+ L(i,i) = 1.0
+end do
+! U matrix is the upper triangular part of A
+do j=1,n
+ do i=1,j
+ U(i,j) = a(i,j)
+ end do
+end do
+
+! Step 3: compute columns of the inverse matrix C
+do k=1,n
+ b(k)=1.0
+ d(1) = b(1)
+! Step 3a: Solve Ld=b using the forward substitution
+ do i=2,n
+ d(i)=b(i)
+ do j=1,i-1
+ d(i) = d(i) - L(i,j)*d(j)
+ end do
+ end do
+! Step 3b: Solve Ux=d using the back substitution
+ x(n)=d(n)/U(n,n)
+ do i = n-1,1,-1
+ x(i) = d(i)
+ do j=n,i+1,-1
+ x(i)=x(i)-U(i,j)*x(j)
+ end do
+ x(i) = x(i)/u(i,i)
+ end do
+! Step 3c: fill the solutions x(n) into column k of C
+ do i=1,n
+ c(i,k) = x(i)
+ end do
+ b(k)=0.0
+end do
+inverse = c
+end function inverse
+
end subroutine MultiSummaryLoop
Modified: pkg/MultiABEL/src/MultiSummaryLoopInbred.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoopInbred.f90 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/src/MultiSummaryLoopInbred.f90 2016-05-07 13:24:56 UTC (rev 2053)
@@ -1,38 +1,172 @@
!
-! MultiSummaryLoop.f90
+! MultiSummaryLoopInbred.f90
!
-! Created by Xia Shen on 5/29/15.
-!
+! Created by Xia Shen on 07/04/16.
-subroutine MultiSummaryLoopInbred(k, m, nn, f, betamat, R, invR, E0, H, coef, pil, fstat)
+subroutine MultiSummaryLoopInbred(k, m, nn, f, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
implicit none
- integer :: i, j, k, m
- double precision, dimension(m, m) :: R, invR, E0, H, HinvE0
- double precision, dimension(1, m) :: beta
- double precision, dimension(m, 1) :: tmp
- double precision, dimension(k) :: nn, f, pil, fstat
- double precision, dimension(k, m) :: betamat, coef
+ integer :: ii, i, j, k, m
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, 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, 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(m) :: lambda
- double precision :: vg
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb
+ double precision :: tmp, vg, abDalphabeta
do j = 1, k
vg = f(j)*(1 - f(j))
- E0 = nn(j)*R
do i = 1, m
- beta(1,i) = betamat(j,i)
+ beta(i,1) = betamat(j,i)
end do
- tmp = matmul(invR, transpose(beta))*vg
do i = 1, m
- coef(j,i) = tmp(i,1)
+ if (i == 1) then
+ R1 = R(2:m,2:m)
+ sdY1 = sdY(2:m,2:m)
+ beta1(1:(m - 1),1) = beta(2:m,1)
+ sY1(1:(m - 1),1) = sY(2:m,1)
+ cor1(1:(m - 1),1) = R(2:m,1)
+ else if (i == m) then
+ R1 = R(1:(m - 1),1:(m - 1))
+ sdY1 = sdY(1:(m - 1),1:(m - 1))
+ beta1(1:(m - 1),1) = beta(1:(m - 1),1)
+ sY1(1:(m - 1),1) = sY(1:(m - 1),1)
+ cor1(1:(m - 1),1) = R(1:(m - 1),m)
+ else
+ R1(1:(i - 1),1:(i - 1)) = R(1:(i - 1),1:(i - 1))
+ R1(i:(m - 1),i:(m - 1)) = R((i + 1):m,(i + 1):m)
+ R1(1:(i - 1),i:(m - 1)) = R(1:(i - 1),(i + 1):m)
+ R1(i:(m - 1),1:(i - 1)) = R((i + 1):m,1:(i - 1))
+ sdY1(1:(i - 1),1:(i - 1)) = sdY(1:(i - 1),1:(i - 1))
+ sdY1(i:(m - 1),i:(m - 1)) = sdY((i + 1):m,(i + 1):m)
+ sdY1(1:(i - 1),i:(m - 1)) = sdY(1:(i - 1),(i + 1):m)
+ sdY1(i:(m - 1),1:(i - 1)) = sdY((i + 1):m,1:(i - 1))
+ beta1(1:(i - 1),1) = beta(1:(i - 1),1)
+ beta1(i:(m - 1),1) = beta((i + 1):m,1)
+ sY1(1:(i - 1),1) = sY(1:(i - 1),1)
+ sY1(i:(m - 1),1) = sY((i + 1):m,1)
+ cor1(1:(i - 1),1) = R(1:(i - 1),i)
+ cor1(i:(m - 1),1) = R((i + 1):m,i)
+ end if
+ zg = beta1*vg
+ zz = matmul(matmul(sdY1, R1), sdY1)
+ tmp = sdY(i,i)
+ do ii = 1, m - 1
+ alpha(ii,1) = cor1(ii,1)/sY1(ii,1)*tmp
+ end do
+
+ alphabeta(1:(m - 1),1) = alpha(1:(m - 1),1)
+ alphabeta(m,1) = beta(i,1)
+ A(1:(m - 1),1:(m - 1)) = zz(1:(m - 1), 1:(m - 1))
+ A(1:(m - 1),m) = beta1(1:(m - 1), 1)*vg
+ tbeta1 = transpose(beta1)
+ A(m,1:(m - 1)) = tbeta1(1,1:(m - 1))*vg
+ A(m,m) = vg
+ do ii = 1, m
+ D(ii,ii) = A(ii,ii)
+ end do
+ invA = inverse(A, m)
+ ab = matmul(matmul(invA, D), alphabeta)
+ abDalphabeta0 = matmul(matmul(transpose(ab), D), alphabeta)
+ abDalphabeta = abDalphabeta0(1,1)
+ Vb = (sdY(i,i)**2 - abDalphabeta)/(nn(j) - 3)*invA(m, m)
+ b(j,i) = ab(m,1)
+ s(j,i) = sqrt(Vb(1,1))
end do
- H = matmul(transpose(beta), beta)*nn(j)*vg
- HinvE0 = matmul(H, invR)/nn(j)
+ H = matmul(beta, transpose(beta))*nn(j)*vg
+ HinvE0 = matmul(matmul(H, invsdY), matmul(invR, invsdY))/nn(j)
do i = 1, m
lambda(i) = HinvE0(i,i)
end do
pil(j) = sum(lambda)
- fstat(j) = pil(j)/m/(1 - pil(j))*(nn(j) - m - 1)
end do
+contains
+
+function inverse(a, n)
+!============================================================
+! Inverse matrix
+! Method: Based on Doolittle LU factorization for Ax=b
+! Alex G. December 2009
+!-----------------------------------------------------------
+! input ...
+! a(n,n) - array of coefficients for matrix A
+! n - dimension
+! output ...
+! c(n,n) - inverse matrix of A
+! comments ...
+! the original matrix a(n,n) will be destroyed
+! during the calculation
+!===========================================================
+implicit none
+integer :: n
+double precision :: a(n,n)
+double precision :: c(n,n), inverse(n, n)
+double precision :: L(n,n), U(n,n), b(n), d(n), x(n)
+double precision :: coeff
+integer :: i, j, k
+
+! step 0: initialization for matrices L and U and b
+! Fortran 90/95 aloows such operations on matrices
+L=0.0
+U=0.0
+b=0.0
+
+! step 1: forward elimination
+do k=1, n-1
+ do i=k+1,n
+ coeff=a(i,k)/a(k,k)
+ L(i,k) = coeff
+ do j=k+1,n
+ a(i,j) = a(i,j)-coeff*a(k,j)
+ end do
+ end do
+end do
+
+! Step 2: prepare L and U matrices
+! L matrix is a matrix of the elimination coefficient
+! + the diagonal elements are 1.0
+do i=1,n
+ L(i,i) = 1.0
+end do
+! U matrix is the upper triangular part of A
+do j=1,n
+ do i=1,j
+ U(i,j) = a(i,j)
+ end do
+end do
+
+! Step 3: compute columns of the inverse matrix C
+do k=1,n
+ b(k)=1.0
+ d(1) = b(1)
+! Step 3a: Solve Ld=b using the forward substitution
+ do i=2,n
+ d(i)=b(i)
+ do j=1,i-1
+ d(i) = d(i) - L(i,j)*d(j)
+ end do
+ end do
+! Step 3b: Solve Ux=d using the back substitution
+ x(n)=d(n)/U(n,n)
+ do i = n-1,1,-1
+ x(i) = d(i)
+ do j=n,i+1,-1
+ x(i)=x(i)-U(i,j)*x(j)
+ end do
+ x(i) = x(i)/u(i,i)
+ end do
+! Step 3c: fill the solutions x(n) into column k of C
+ do i=1,n
+ c(i,k) = x(i)
+ end do
+ b(k)=0.0
+end do
+inverse = c
+end function inverse
+
end subroutine MultiSummaryLoopInbred
Modified: pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90
===================================================================
--- pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90 2016-04-16 20:39:42 UTC (rev 2052)
+++ pkg/MultiABEL/src/MultiSummaryLoopPrecise.f90 2016-05-07 13:24:56 UTC (rev 2053)
@@ -1,40 +1,174 @@
!
-! MultiSummaryLoop.f90
+! MultiSummaryLoopPrecise.f90
!
-! Created by Xia Shen on 5/29/15.
-!
+! Created by Xia Shen on 07/04/16.
-subroutine MultiSummaryLoopPrecise(k, m, nn, betamat, R, invR, E0, H, coef, pil, fstat, vg)
+subroutine MultiSummaryLoopPrecise(k, m, nn, varg, betamat, R, invR, D, sdY, invsdY, sY, b, s, pil)
implicit none
- integer :: i, j, k, m
- integer, dimension(k) :: nn
- double precision, dimension(m, m) :: R, invR, E0, H, HinvE0
- double precision, dimension(1, m) :: beta
- double precision, dimension(m, 1) :: tmp
- double precision, dimension(k) :: pil, fstat, vg
- double precision, dimension(k, m) :: betamat, coef
+ integer :: ii, i, j, k, m
+ double precision, dimension(m, m) :: sdY, invsdY, R, invR, H, HinvE0, 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, 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(m) :: lambda
+ double precision, dimension(1, 1) :: abDalphabeta0, Vb
+ double precision :: tmp, vg, abDalphabeta
do j = 1, k
- !if (j < 2) then
- ! write(*,*) nn(j)
- !end if
- E0 = nn(j)*R
+ vg = varg(j)
do i = 1, m
- beta(1,i) = betamat(j,i)
+ beta(i,1) = betamat(j,i)
end do
- tmp = matmul(invR, transpose(beta))*vg(j)
do i = 1, m
- coef(j,i) = tmp(i,1)
+ if (i == 1) then
+ R1 = R(2:m,2:m)
+ sdY1 = sdY(2:m,2:m)
+ beta1(1:(m - 1),1) = beta(2:m,1)
+ sY1(1:(m - 1),1) = sY(2:m,1)
+ cor1(1:(m - 1),1) = R(2:m,1)
+ else if (i == m) then
+ R1 = R(1:(m - 1),1:(m - 1))
+ sdY1 = sdY(1:(m - 1),1:(m - 1))
+ beta1(1:(m - 1),1) = beta(1:(m - 1),1)
+ sY1(1:(m - 1),1) = sY(1:(m - 1),1)
+ cor1(1:(m - 1),1) = R(1:(m - 1),m)
+ else
+ R1(1:(i - 1),1:(i - 1)) = R(1:(i - 1),1:(i - 1))
+ R1(i:(m - 1),i:(m - 1)) = R((i + 1):m,(i + 1):m)
+ R1(1:(i - 1),i:(m - 1)) = R(1:(i - 1),(i + 1):m)
+ R1(i:(m - 1),1:(i - 1)) = R((i + 1):m,1:(i - 1))
+ sdY1(1:(i - 1),1:(i - 1)) = sdY(1:(i - 1),1:(i - 1))
+ sdY1(i:(m - 1),i:(m - 1)) = sdY((i + 1):m,(i + 1):m)
+ sdY1(1:(i - 1),i:(m - 1)) = sdY(1:(i - 1),(i + 1):m)
+ sdY1(i:(m - 1),1:(i - 1)) = sdY((i + 1):m,1:(i - 1))
+ beta1(1:(i - 1),1) = beta(1:(i - 1),1)
+ beta1(i:(m - 1),1) = beta((i + 1):m,1)
+ sY1(1:(i - 1),1) = sY(1:(i - 1),1)
+ sY1(i:(m - 1),1) = sY((i + 1):m,1)
+ cor1(1:(i - 1),1) = R(1:(i - 1),i)
+ cor1(i:(m - 1),1) = R((i + 1):m,i)
+ end if
+ zg = beta1*vg
+ zz = matmul(matmul(sdY1, R1), sdY1)
+ tmp = sdY(i,i)
+ do ii = 1, m - 1
+ alpha(ii,1) = cor1(ii,1)/sY1(ii,1)*tmp
+ end do
+
+ alphabeta(1:(m - 1),1) = alpha(1:(m - 1),1)
+ alphabeta(m,1) = beta(i,1)
+ A(1:(m - 1),1:(m - 1)) = zz(1:(m - 1), 1:(m - 1))
+ A(1:(m - 1),m) = beta1(1:(m - 1), 1)*vg
+ tbeta1 = transpose(beta1)
+ A(m,1:(m - 1)) = tbeta1(1,1:(m - 1))*vg
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 2053
More information about the Genabel-commits
mailing list