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