# 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)
}