R/simulation.R

In osbosb/DEClust: DE analysis based on Clustering

generateSimData <- function() {
  Ngene <- 10000
  alpha.true <- rnorm(Ngene, 4, 1)
  beta.true <- matrix(c(0, 0, 0, 1, -.5, -.5, -.5, 1, -.5, -.5, -.5, 1), 3, 4)
  Ngroup <- nrow(beta.true)
  K.true <- ncol(beta.true)
  norm.true <- rnorm(9, 0, 0)
  normalizer <- matrix(rep(norm.true, each = Ngene), Ngene, 9)
  phi.true <- rep(0.1, Ngene)
  treatment <- gl(3, 3)
  depatterns <- c(rep(1, 9500), rep(2, 200), rep(3, 200), rep(4, 100))
  pi <- table(depatterns) / length(depatterns)
  lambda <- rep(alpha.true, times = 9) + as.vector(t(beta.true[treatment, depatterns])) + rep(norm.true, times = Ngene)
  lambda <- exp(lambda)
  n <- rnbinom(Ngene * 9, size = rep(1 / phi.true, times = 9), mu = lambda)
  dim(n) <- c(Ngene, 9)

  logL <- vapply(1:K.true, function(i) {
    lambda <- rep(alpha.true, times = 9) + as.vector(t(beta.true[treatment, rep(i, Ngene)])) + rep(norm.true, times = Ngene)
    lambda <- exp(lambda)
    logL <- dnbinom(n, size = rep(1 / phi.true, times = 9), mu = lambda, log = TRUE)
    dim(logL) <- c(Ngene, 9)
    apply(logL, 1, sum)
  }, vector("numeric", Ngene))
  logL <- sweep(logL, 2, log(pi), "+")
  logL <- sum(apply(logL, 1, function(x) {
    log(sum(exp(x - max(x)))) + max(x)
  }))
  return(list(count = n, treatment = treatment, pi = pi, alpha = alpha.true,
              beta = beta.true, phi = phi.true, normalizer = normalizer, logL = logL,
              cluster = depatterns, normalizer = norm.true, k = K.true))
}

#sim <- generateSimData()
#K <- 2:6
#decls <- DEClust(sim$count, sim$treatment, K, normalizer = sim$normalizer, phi = sim$phi)
##mymb <- mymbcluster(sim$count, sim$treatment, 4, normalizer = sim$normalizer, phi = sim$phi)
#aic <- sapply(decls, "[[", "AIC")
#bic <- sapply(decls, "[[", "BIC")
#logL <- sapply(decls, "[[", "logL")
#
#plot(0, 0, type = "n", xlim = range(K), ylim = range(c(logL, sim$logL)), xlab = "k", ylab = "logL")
#lines(K, logL)
#lines(K, rep(sim$logL, length(K)), col = "red")