R/simulateData.R

In ChongWu-Biostat/MiSPU: Microbiome Based Sum of Powered Score (MiSPU) Tests

###############################################################
# Chong Wu
# email: wuxx0845@umn.edu
###############################################################
# Intent:
#   This function suppresses the following notes generated by "R CMD check":
#   - "Note: no visible binding for global variable '.->ConfigString'"
#   - "Note: no visible binding for '<<-' assignment to 'ConfigString'"
# Usage:
#   Add the following right in the beginning of the .r file (before the Reference
#   class is defined in the sourced .r file):
#   suppressBindingNotes(c(".->ConfigString","ConfigString"))
#suppressBindingNotes <- function(variablesMentionedInNotes) {
#    for(variable in variablesMentionedInNotes) {
#        assign(variable,NULL, envir = .GlobalEnv)
#    }
#}

#suppressBindingNotes(c(".->ConfigString","ConfigString"))

simulateData <- function(nSam=100, s=12,ncluster = 20,mu = 1000, size = 25) {
    
    data(throat.tree,envir = environment())
    data(dd,envir = environment())
    
    tree <- get("throat.tree", envir  = environment())
    dd1 = get("dd",envir  = environment())
    tree.dist <- cophenetic(tree)
    obj <- pam(as.dist(tree.dist), ncluster,diss =TRUE)
    clustering <- obj$clustering
    otu.ids <- tree$tip.label
    
    p.est = dd1$pi
    names(p.est) <- names(dd1$pi)
    theta <- dd1$theta
    gplus <- (1 - theta) / theta
    p.est <- p.est[otu.ids]
    g.est <- p.est * gplus
    p.clus <- sort(tapply(p.est, clustering, sum), decreasing=T)
    scale2 = function(x)as.numeric(scale(x))
    
    
    comm <- matrix(0, nSam, length(g.est))
    rownames(comm) <- 1:nrow(comm)
    colnames(comm) <- names(g.est)
    # comm.p hold the underlying proportions
    comm.p <- comm
    nSeq <- rnbinom(nSam, mu = mu, size = size)
    for (i in 1:nSam) {
        comm.p[i, ] <- rdirichlet(1, g.est)[1, ]
        comm[i, ] <- rmultinom(1, nSeq[i], prob=comm.p[i, ])[, 1]
    }
    
    otu.ids <- names(which(clustering == s))
    
    # No additional covariates in this case.
    OTU = comm[, otu.ids]
    return(list(informative.OTU = OTU,whole.OTU = OTU))
}