# Time-stamp: <2011-03-01 19:06:40 zaykind> (Dmitri Zaykin) # Chi.Com computes a pooled Z-test for the shared controls design # The combined p-value is converted back to two-sided # (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 (1 df chi-square p-values are assumed) # eff: vector of effect signs pooled.Z <- function(n, pvalues, eff) { nn <- length(n) nca <- nn-1 if(nca != length(pvalues)) { stop("pvalue vector length should be one less than the sample sizes vector") } if(nca != length(eff)) { stop("effect sign 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)) } cr <- matrix(nrow=nca, ncol=nca) diag(cr) <- 1 for(i in 2 : (nn-1)) { for(j in (i+1) : nn) { cr[i-1, j-1] <- cr[j-1, i-1] <- rho(n[1], n[i], n[j]) } } wgt1 <- (n[2:nn] / sum(n[2:nn]))^0.5 wgt3 <- sqrt(diag(solve(cr)))*wgt1 Rho3 <- array(0, nca*(nca-1)/2) ii <- 1 for(i in 2 : (nn-1)) { for(j in (i+1) : nn) { Rho3[ii] <- rho(n[1], n[i], n[j]) * wgt3[i-1] * wgt3[j-1] ii <- ii+1 } } varz3 <- sqrt(sum(wgt3^2) + 2*sum(Rho3)) zs <- array(0, nca) for(j in 1:nca) { zs.tmp <- qnorm(pvalues[j]/2, lower.tail=FALSE) zs[j] <- ifelse(eff[j] >= 0, zs.tmp, -zs.tmp) } nom <- sum(wgt3*zs) zstat <- nom / varz3 p1 <- pnorm(zstat, lower.tail=FALSE) p2 <- ifelse(p1 < 0.5, 2*p1, 2*(1-p1)) return (p2) }