### R functions ### rm(list=ls()) dyn.load("Concordance.dll") sim.multinom <- function(p){ return( c(rmultinom(1,1,p)) ) } sim.all <- function(k, p, pi){ m <- dim(p)[1] z <- apply(p, 1, sim.multinom) w <- rowSums(z)>m*pi return( c(w, (w[2] || w[3])) ) } sim.es.all <- function(p, pi, B=1000){ return( rowMeans( sapply(1:B, sim.all, p, pi) ) ) } #modified from package mclust unmap <- function(classification){ n <- length(classification) u <- sort(unique(classification)) labs <- as.character(u) k <- length(u) z <- matrix(0, n, k) for (j in 1:k) z[classification == u[j], j] <- 1 dimnames(z) <- list(NULL, labs) return(z) } em.normal.univariate <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ n <- as.integer(length(data)) g <- as.integer(nlevels(as.factor(class))) z <- unmap(class) pi <- double(g) mu <- double(g) sigma <- double(g) loglik <- double(1) convergence <- integer(1) if(restriction>0){ if(length(constrain)==g){ results <- .C("em_normal_univariate", as.double(data), as.double(z), n, pi, mu, sigma, as.integer(g), loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) }else{ print("Error with constrain!") return(0) } }else{ results <- .C("em_normal_univariate", as.double(data), as.double(z), n, pi, mu, sigma, as.integer(g), loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) } print(paste("convergence within", results[[12]], "run?", results[[13]], sep=" ")) return(list(model="univariate", convergence=results[[13]], pi=results[[4]], mu=results[[5]], sigma=results[[6]], loglik=results[[8]], class=apply(array(results[[2]], dim=c(n,g)),1,order,decreasing=T)[1,], z=array(results[[2]], dim=c(n,g)))) } em.normal.complete.discordant <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ r1 <- em.normal.univariate(data[,1], class[,1], tol, restriction, constrain, iteration) r2 <- em.normal.univariate(data[,2], class[,2], tol, restriction, constrain, iteration) return(list(model="CD", convergence=c(r1\$convergence, r2\$convergence), pi=rbind(r1\$pi, r2\$pi), mu_sigma=rbind(r1\$mu, r1\$sigma), nu_tau=rbind(r2\$mu, r2\$sigma), loglik=r1\$loglik+r2\$loglik, class=cbind(r1\$class, r2\$class), z1=r1\$z, z2=r2\$z)) } em.normal.complete.concordant <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ n <- as.integer(dim(data)[1]) g <- as.integer(nlevels(as.factor(class))) z <- unmap(class) pi <- double(g) mu <- double(g) sigma <- double(g) nu <- double(g) tau <- double(g) loglik <- double(1) convergence <- integer(1) if(restriction>0){ if(length(constrain)==g){ results <- .C("em_normal_complete_concordant", as.double(data[,1]), as.double(data[,2]), as.double(z), n, pi, mu, sigma, nu, tau, as.integer(g), loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) }else{ print("Error with constrain!") return(0) } }else{ results <- .C("em_normal_complete_concordant", as.double(data[,1]), as.double(data[,2]), as.double(z), n, pi, mu, sigma, nu, tau, as.integer(g), loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) } print(paste("convergence within", results[[15]], "run?", results[[16]], sep=" ")) return(list(model="CC", convergence=results[[16]], pi=results[[5]], mu_sigma=rbind(results[[6]], results[[7]]), nu_tau=rbind(results[[8]], results[[9]]), loglik=results[[11]], class=apply(array(results[[3]], dim=c(n,g)),1,order,decreasing=T)[1,], z=array(results[[3]], dim=c(n,g)))) } em.normal.partial.concordant <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ n <- as.integer(dim(data)[1]) g <- as.integer(nlevels(as.factor(class))) yl.outer <- function(k, zx, zy){ return( c(zx[k,] %o% zy[k,]) ) } yl.diag <- function(k, z){ return( c(diag(z[k,])) ) } zx <- unmap(class[,1]) zy <- unmap(class[,2]) zxy <- sapply(1:dim(zx)[1], yl.outer, zx, zy) pi <- double(g*g) mu <- double(g) sigma <- double(g) nu <- double(g) tau <- double(g) loglik <- double(1) convergence <- integer(1) if(restriction>0){ if(length(constrain)==g){ results <- .C("em_normal_partial_concordant", as.double(data[,1]), as.double(data[,2]), as.double(t(zxy)), n, pi, mu, sigma, nu, tau, g, loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) }else{ print("Error with constrain!") return(0) } }else{ results <- .C("em_normal_partial_concordant", as.double(data[,1]), as.double(data[,2]), as.double(t(zxy)), n, pi, mu, sigma, nu, tau, g, loglik, as.double(tol), as.integer(restriction), as.integer(constrain), as.integer(iteration), convergence) } print(paste("convergence within", results[[15]], "run?", results[[16]], sep=" ")) return(list(model="PCD", convergence=results[[16]], pi=t(array(results[[5]],dim=c(g,g))), mu_sigma=rbind(results[[6]], results[[7]]), nu_tau=rbind(results[[8]], results[[9]]), loglik=results[[11]], class=apply(array(results[[3]], dim=c(n,g*g)),1,order,decreasing=T)[1,], z=array(results[[3]], dim=c(n,g*g)))) } concordance.test.discordant.partial <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ r0 <- em.normal.complete.discordant(data, class, tol, restriction, constrain, iteration) r1 <- em.normal.partial.concordant(data, class, tol, restriction, constrain, iteration) test <- NA if(min(c(r0\$convergence, r1\$convergence))>0){ test <- r1\$loglik - r0\$loglik } return(test) } concordance.test.concordant.partial <- function(data, class, tol=0.000001, restriction=0, constrain=0, iteration=1000){ r0 <- em.normal.complete.concordant(data, apply(class,1,min), tol, restriction, constrain, iteration) r1 <- em.normal.partial.concordant(data, class, tol, restriction, constrain, iteration) test <- NA if(min(c(r0\$convergence, r1\$convergence))>0){ test <- r1\$loglik - r0\$loglik } return(test) } sim.normal.complete.discordant <- function(n, g, pi.x, mu.x, sigma.x, pi.y, mu.y, sigma.y){ data.x <- sample(1:g, n, replace=T, prob=pi.x) data.y <- sample(1:g, n, replace=T, prob=pi.y) for(i in 1:g){ l <- length(data.x[data.x==i]) if(l>0){ data.x[data.x==i] <- rnorm(l, mean=mu.x[i], sd=sqrt(sigma.x[i])) } l <- length(data.y[data.y==i]) if(l>0){ data.y[data.y==i] <- rnorm(l, mean=mu.y[i], sd=sqrt(sigma.y[i])) } } return(cbind(data.x, data.y)) } sim.normal.complete.concordant <- function(n, g, pi, mu.x, sigma.x, mu.y, sigma.y){ data <- sample(1:g, n, replace=T, prob=pi) data.x <- data data.y <- data for(i in 1:g){ l <- length(data[data==i]) if(l>0){ data.x[data==i] <- rnorm(l, mean=mu.x[i], sd=sqrt(sigma.x[i])) data.y[data==i] <- rnorm(l, mean=mu.y[i], sd=sqrt(sigma.y[i])) } } return(cbind(data.x, data.y)) } ### Running samples ### #simulate paired z-scores n <- 3000 g <- 3 set.seed(123) pdata <- sim.normal.complete.concordant(n, g, pi=c(0.6,0.2,0.2), c(0,-1.5,2), c(1,1.5,2), c(0,-2,1.5), c(1,2,1.5)) class <- cbind(rep(0, dim(pdata)[1]), rep(0, dim(pdata)[1])) class[pdata[,1]>0+1,1] <- 2 class[pdata[,1]<0-1,1] <- 1 class[pdata[,2]>0+1,2] <- 2 class[pdata[,2]<0-1,2] <- 1 #estimate models for parametric bootstrap cd <- em.normal.complete.discordant(pdata, class, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) cc <- em.normal.complete.concordant(pdata, apply(class,1,min), tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) pd <- em.normal.partial.concordant(pdata, class, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) #test of CC/CD t.dp <- concordance.test.discordant.partial(pdata, class, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) t.cp <- concordance.test.concordant.partial(pdata, class, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) #significance for CD: reject CD since t.dp > all t0.dp n <- dim(cd\$class)[1] g <- dim(cd\$mu_sigma)[2] t0.dp <- double(1000) i <- 1 while(i<=length(t0.dp)){ sdata <- sim.normal.complete.discordant(n, g, cd\$pi[1,], cd\$mu_sigma[1,], cd\$mu_sigma[2,], cd\$pi[2,], cd\$nu_tau[1,], cd\$nu_tau[2,]) sclass <- cbind(rep(0, dim(sdata)[1]), rep(0, dim(sdata)[1])) sclass[sdata[,1]>0+1,1] <- 2 sclass[sdata[,1]<0-1,1] <- 1 sclass[sdata[,2]>0+1,2] <- 2 sclass[sdata[,2]<0-1,2] <- 1 t0.dp[i] <- concordance.test.discordant.partial(sdata, sclass, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) print(c(t.dp, i, t0.dp[i])) i <- i+1 } length(t0.dp) sum(!is.na(t0.dp)) sum(!is.na(t0.dp) & (t0.dp>=t.dp)) sum(!is.na(t0.dp) & (t0.dp>=t.dp)) / sum(!is.na(t0.dp)) #significance for CC: not reject CC since t.cp ~ t0.cp n <- length(cc\$class) g <- dim(cc\$mu_sigma)[2] t0.cp <- double(1000) i <- 1 while(i<=length(t0.cp)){ sdata <- sim.normal.complete.concordant(n, g, cc\$pi, cc\$mu_sigma[1,], cc\$mu_sigma[2,], cc\$nu_tau[1,], cc\$nu_tau[2,]) sclass <- cbind(rep(0, dim(sdata)[1]), rep(0, dim(sdata)[1])) sclass[sdata[,1]>0+1,1] <- 2 sclass[sdata[,1]<0-1,1] <- 1 sclass[sdata[,2]>0+1,2] <- 2 sclass[sdata[,2]<0-1,2] <- 1 t0.cp[i] <- concordance.test.concordant.partial(sdata, sclass, tol=0.001, restriction=1, constrain=c(0,-1,1), iteration=1000) print(c(t.cp, i, t0.cp[i])) i <- i+1 } length(t0.cp) sum(!is.na(t0.cp)) sum(!is.na(t0.cp) & (t0.cp>=t.cp)) sum(!is.na(t0.cp) & (t0.cp>=t.cp)) / sum(!is.na(t0.cp)) #confirmed model cc <- em.normal.complete.concordant(pdata, apply(class,1,min), tol=0.000001, restriction=1, constrain=c(0,-1,1), iteration=10000) #random gene set analysis plot(pdata[,1], pdata[,2], main="A random gene set", xlab="z-score (data set 1)", ylab="z-score (data set 2)", pch=19, col="gray75", cex.lab=1.5, cex.axis=1.5, cex=1.25) abline(h=0, col="gray50", lwd=1, lty=3) abline(v=0, col="gray50", lwd=1, lty=3) set.seed(123) h <- sample(1:dim(pdata)[1], 100) lines(pdata[h,1], pdata[h,2], type="p", pch=19, col="gray0", cex=1.5) print( sim.es.all(cc\$z[h,], cc\$pi) )