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