# Time-stamp: <2011-03-01 19:04:21 zaykind> (Dmitri Zaykin)
# See comments above "conditional.Z" which evaluates
# \Pr(P_j < p_j | other p-values)
# for the shared control design
# (Zaykin, Kozbur 2010; http://www.ncbi.nlm.nih.gov/pubmed/20976797)

Params <- function (cv, p) {
  s11 <- cv[1,1]
  s12 <- cv[1,][-1]
  s21 <- cv[,1][-1]
  s22 <- cv[-1,-1]
  s22i <- solve(cv[-1,-1])
  m <- s12 %*% s22i %*% p
  v <- s11 - s12 %*% s22i %*% s21
  return ( c(m, v) )
}

rho <- function(n0, ni, nj) {
  sqrt(1/(1+n0/ni) * 1/(1+n0/nj))
}

# 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
# idx: denotes index of the p-value for which the function will do the
# adjustment (first p-value is assumed if idx is not specified)
conditional.Z <- function(n, pvalues, eff, idx=1) {
  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")
  }

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

  iRh <- solve(cr)
  fct2 <- 1/diag(iRh)
  fct <- 1/sqrt(fct2)
  wV <- diag(fct) %*% cr %*% t(diag(fct))
  
  zs <- array(0, nca)
  for(j in 1:nca) {
    p2z <- ifelse(eff[j] < 0, 1-pvalues[j]/2, pvalues[j]/2)
    zs[j] <- qnorm(1 - p2z)
  }

  zw <- zs*fct
  x <- wV
  
  if(idx>1) {
    x[, c(1,idx)] <- x[, c(idx,1)]
    x[c(1,idx), ] <- x[c(idx,1), ]
    zw[c(1,idx)] <- zw[c(idx,1)]
  }
  mv <- Params(x, zw[-1])
  p2 <- pchisq(zw[1]^2/(mv[2]), ncp=mv[1]^2, df=1, lower.tail=FALSE)
  return (p2)
}