# Time-stamp: <2011-03-01 19:02:48 zaykind> (Dmitri Zaykin) # depends on "CompQuadForm"; try install.packages("CompQuadForm") library(CompQuadForm) # Chi.Com computes a pooled chi-square test for the shared controls design # (Zaykin, Kozbur 2010; http://www.ncbi.nlm.nih.gov/pubmed/20976797) # # n: sample sizes vector (n0, n1, ...). n0 is the shared group sample size # pvalues: p-values vector (chi-square p-values) # dfs: degrees of freedom (1 if not specified) pooled.Chi <- function(n, pvalues, dfs=1) { xs <- qchisq(pvalues, df=dfs, lower.tail = FALSE) nn <- length(n) nca <- nn-1 if(nca != length(xs)) { stop("pvalue vector length should be one less than the sample sizes vector") } rho <- function(n0, ni, nj) { sqrt(1/(1+n0/ni) * 1/(1+n0/nj)) } Rh <- matrix(nrow=nca, ncol=nca) diag(Rh) <- 1 ii <- 1 for(i in 2 : (nn-1)) { for(j in (i+1) : nn) { Rh[i-1, j-1] <- Rh[j-1, i-1] <- rho(n[1], n[i], n[j]) ii <- ii+1 } } dmp <- 0 wgt1 <- ( (n[2:nn]+dmp) / (sum(n[2:nn]) + dmp) )^0.5 wgt3 <- sqrt(diag(solve(Rh)))*wgt1 wRh3 <- diag(wgt3^0.5) %*% Rh %*% t(diag(wgt3^0.5)) eX3 <- eigen(wRh3)\$values err <- 0; res <- 1 cc <- sum(wgt3*xs) pv <- farebrother(cc, lambda=eX3, h=rep(dfs,nca))\$res return(pv) }