Description Usage Arguments Details Value Author(s) Examples

For a set of phenotypes, given summary results for association tests for each single phenotype, over a large fixed set of marker genotypes (e.g. GWAS results), this function calculates a multi-phenotype association test for each marker that is equivalent to using the subject-specific data to perform a multivariate analysis of variance for each marker.

1 |

`z` |
a matrix of association Z statistics, with rows corresponding to markers and columns corresponding to phenotypes. |

The function is named for the close correspondence with Hotelling's T2 test.

It is the user's responsibility to ensure that the columns of
“`z`

” are marginally well calibrated, i.e. over a set of
null markers, each column of “`z`

” should be standard
normal marginal to the other columns of “`z`

”.

Rank deficient situations are currently not handled.

A list with two elements (both of length `nrow(z)`

).
“`T2`

” contains the test statistics, which are chi-squared
with `ncol(z)`

d.f. under the null, and “`pval`

”
contains the P-values.

Toby Johnson Toby.x.Johnson@gsk.com

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | ```
## Not run:
nsubj <- 400 # number of subjects
nsnp <- 4000 # intended number of SNPs in GWAS
snps <- replicate(nsnp, rbinom(nsubj, 2, rbeta(1, 1.2, 1.2)))
## simulate with random allele frequencies
snps <- snps[ , apply(snps, 2, var) > 0]
nsnp <- ncol(snps) # number after filtering monomorphic
beta <- matrix(rnorm(30, 0, 0.1), ncol = 3)
## matrix of effects of 10 snps on 3 phenotypes
y1 <- rnorm(nsubj)
y2 <- .1*y1 + rnorm(nsubj)
y3 <- .1*y1 + .3*y2 + rnorm(nsubj) # simulate correlated errors
y <- cbind(y1, y2, y3) + snps[ , 1:10] %*% beta
## wlog the truely associated snps are 1:10
rm(y1, y2, y3)
astats <- function(ii) {
with(list(snp = snps[ , ii]),
c(coef(summary(lm(y[ , 1] ~ snp)))["snp", 3],
coef(summary(lm(y[ , 2] ~ snp)))["snp", 3],
coef(summary(lm(y[ , 3] ~ snp)))["snp", 3],
summary(manova(y ~ snp))$stats["snp", 6]))
}
system.time(gwas <- t(sapply(1:nsnp, astats)))
## columns 1-3 contain single phenotype Z statistics
## column 4 contains manova P-value
pv <- multipheno.T2(gwas[ , 1:3])$pval
plot(-log10(gwas[ , 4]), -log10(pv), type = "n",
xlab = "MANOVA -log10(P)", ylab = "Summary statistic -log10(P)", las = 1)
points(-log10(gwas[-(1:10), 4]), -log10(pv[-(1:10)]))
points(-log10(gwas[1:10, 4]), -log10(pv[1:10]), cex = 2, pch = 21, bg = "red")
legend("topleft", pch = c(1, 21), pt.cex = c(1, 2), pt.bg = c("white", "red"),
legend = c("null SNPs", "associated SNPs"))
## nb approximation gets better as nsnp becomes large, even with small nsubj
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.