Nothing
hapPower <- function(hapRun.result, alpha = 0.05){
##
## Power simulations
##
##
#
## Extract information from hapRun object
.hap.run <- hapRun.result
.hapfunc <- attributes(hapRun.result)$hapfunc
#
## Extract p-values from .hap.run
if(.hapfunc == "haplin"){
.n.sim <- nrow(unique(.hap.run["sim.no"]))
if(any(is.na(.hap.run["haplos"]))) .hap.run <- .hap.run[-which(is.na(.hap.run["haplos"])),]
.no.files <- nrow(unique(.hap.run["sim.no"]))
.hap.run <- cbind(.hap.run["haplos"],.hap.run["reference"],.hap.run["pv.overall"],.hap.run[grep("p.value",names(.hap.run))])
.reference <- unique(.hap.run[which(.hap.run["reference"]=="ref"),which(colnames(.hap.run)=="haplos")])
if(length(.reference)> 1) stop("More than one haplotype have been chosen as reference", call. = F)
.hap.run <- .hap.run[,-which(colnames(.hap.run)=="reference")]
.pvalues <- lapply(unique(.hap.run$haplos), function(x){
.pvalues <- .hap.run[.hap.run$haplos == x,]
.pvalues$haplos <- NULL
.pvalues
})
.power <- cbind(haplos=unique(.hap.run$haplos),as.data.frame(t(sapply(.pvalues, function(x) colMeans(x <= alpha, na.rm = TRUE)))))
.power[.power$haplos==.reference,which(is.na(.power[.power$haplos==.reference,]))] <- "ref"
if(any(is.na(.hap.run[-which(.hap.run$haplos==.reference),]))) warning("Some p-values were estimated to NA. These were removed from the power calculations", call. = F)
}else if(.hapfunc == "haplinSlide"){
.n.sim <- length(.hap.run)
.pvalues <- sapply(.hap.run, function(x){
if(is.null(attr(x, "Error"))) x$pval.alt.corr
else NA
})
.no.files <- sum(sapply(.hap.run, function(x){
if(is.null(attr(x, "Error"))) TRUE
else FALSE
}))
if(length(which(!is.na(.pvalues)))!= .no.files) warning("Some p-values were estimated to NA. These were removed from the power calculations", call. = F)
.power <- mean(.pvalues <= alpha, na.rm = TRUE)
}else if(.hapfunc == "haplinStrat"){
.n.sim <- length(.hap.run)
.error <- sapply(.hap.run, function(x){
if(is.null(attr(x, "Error"))) TRUE
else FALSE
})
.hap.run[which(!.error)] <- NULL
.pvalues <- lapply(.hap.run, function(x){
.pvalues <- x$gxe.test$pval
names(.pvalues) <- x$gxe.test$gxe.test
.pvalues <- .pvalues[-which(names(.pvalues)=="haplo.freq")]
})
.pvalues <- do.call("cbind",.pvalues)
.no.files <- length(.hap.run)
if(any(is.na(.pvalues))) warning("Some p-values were estimated to NA. These were removed from the power calculations", call. = F)
.power <- rowMeans(.pvalues <= alpha, na.rm = TRUE)
}
#
## Output
.names <- names(.power)
.names <- sub("p.value","power",.names)
.names <- sub("pv.overall","overall.power",.names)
names(.power) <- .names
if(.hapfunc == "haplinSlide"){
names(.power) <- "Power (corrected for multiple testing)"
warning("\"haplinSlide\" is only partially implemented, and the results should be interpreted with caution. Please see the \"hapRun\" help file for more information", call. = F)
}
#
cat(paste("\nThe power was calculated using ", .no.files, " of ", .n.sim, " files\n", sep = ""))
#
return(.power)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.