#' Simulate \code{pants} with mediation to estimate size & power
#'
#' Simulate \code{pants} with mediation to estimate size (effect=0) & power (effect>0).
#'
#' @inheritParams ezlimma::sim_barfield
#' @inheritParams pants
#' @inheritParams sim_pants
#' @seealso \code{\link{sim_pants}}
# test mult pwys for efficiency
sim_pants_mediation <- function(Gmat, exposure, effect.v=c(0, 0.2), alpha=0.05, nsim=10**3,
nperm=10**3, ncores=1, seed=1, verbose=TRUE, ker=NULL){
prop.sig.mat <- matrix(NA, nrow=nsim, ncol=length(effect.v),
dimnames=list(paste0("sim", 1:nsim), paste0("eff_", effect.v)))
set.seed(seed)
err <- stats::rnorm(n=length(exposure), sd=stats::sd(exposure)/2)
phenotype <- exposure + err
for (sim in 1:nsim){
for (ev in effect.v){
obj.test <- matrix(stats::rnorm(n=nrow(Gmat)*length(phenotype), sd=1), ncol=length(phenotype), nrow=nrow(Gmat),
dimnames=list(rownames(Gmat), names(phenotype)))
if (ev > 0){
obj.test <- t(t(obj.test) + phenotype + ev*err)
}
# need to randomize seed to get variation in prop sig
res <- pants(object=obj.test, exposure=exposure, phenotype=phenotype, Gmat=Gmat, nperm=nperm, type="mediation",
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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.