R/sim_pants.R

In jdreyf/PANTS: Pathway analysis via network smoothing

#' Simulate \code{pants} to estimate size & power
#' 
#' Simulate \code{pants} to estimate size (effect=0) & power (effect>0) with \code{score_fcn=identity} vs. 
#' H1: "greater". Does not support \code{type=="mediation"}.
#' 
#' @param effect.v Numeric vector of log fold-changes or percent of phentotypes to add.
#' @inheritParams ezlimma::sim_barfield
#' @inheritParams pants
#' @seealso \code{\link{sim_pants_mediation}}

# contr
# test mult pwys for efficiency
sim_pants <- function(Gmat, phenotype, type=c("contrasts", "correlation"), 
                      effect.v=c(0, 0.2), alpha=0.05, nsim=10**3,
                      nperm=10**3, ncores=1, seed=1, verbose=TRUE, ker=NULL){
  type <- match.arg(type)
  if (is.null(names(phenotype))) names(phenotype) <- paste0("s", 1:length(phenotype))
  
  prop.sig.mat <- matrix(NA, nrow=nsim, ncol=length(effect.v), 
                         dimnames=list(paste0("sim", 1:nsim), paste0("eff_", effect.v)))
  
  if (type=="contrasts"){
    contrast.v <- c(vs=paste(unique(phenotype)[2], unique(phenotype)[1], sep="-"))
  } else {
    contrast.v <- NULL
  }
  
  set.seed(seed)
  for (sim in 1:nsim){
    for (ev in effect.v){
      obj.test <- matrix(stats::rnorm(n=nrow(Gmat)*length(phenotype)), ncol=length(phenotype), nrow=nrow(Gmat),
                         dimnames=list(rownames(Gmat), names(phenotype)))
      if (ev > 0){
        if (type=="contrasts"){
          obj.test[, phenotype == unique(phenotype)[2]] <- obj.test[, phenotype == unique(phenotype)[2]] + ev
        } else {
          obj.test <- obj.test + ev*matrix(phenotype, nrow=nrow(obj.test), ncol=ncol(obj.test), byrow = TRUE)
        }
      }
      res <- pants(object=obj.test, Gmat=Gmat, phenotype=phenotype, type=type, score_fcn=abs,
                   contrast.v=contrast.v, nperm=nperm, seed=sample(10**5, size = 1), ker=ker, ncores=ncores)
      prop.sig.mat[sim, paste0("eff_", ev)] <- mean(res$pwy.stats[, "p"] < alpha)
    }
    if (sim %% 100 == 0) cat("sim: ", sim, "\n")
  }
  (prop.sig.mat <- rbind(avg=colMeans(prop.sig.mat), prop.sig.mat))
}