Nothing
# Version: 07-07-2013
# Changes:
# Added the options "rejLine" and "alpha"
# Version: 25-05-2013, Daniel Fischer
# Changes:
# Added the rej.lty parameter in order to change the line type of the rejection line
# Added the rej.col parameter in order to change the color of the rejection line
# Changed the default of crit to NULL, so that by default the lower part of the plot isn't plotted!
# rej.para accepts now also vectors as entry and plots then several rejection lines.
# Fixed a bug in the matrix input (couldn't handle data.frame input)
# 07-07-2013: Changed the Bonferroni slope (it depends now also on the alpha value
# 21-09-2013: Added the option for a legend
# Added automatical colors, if not enough are given.
rejectionPlot <- function(X, lCol="red", xlim=NULL, crit=NULL, rejLine=NULL, alpha=0.01, rej.lty=c("dotted"), rej.col="black", incLegend=NULL){
# Input check for the option 'crit'
if(!is.null(crit)){
crit <- match.arg(crit,c("distance","ratio"))
} else {
crit <- "NULL"
}
# Input check for the option 'rejLine'
if(!is.null(rejLine)) rejLine <- match.arg(rejLine,c("bh","bonferroni","simes"))
# X-axis dimension setting
if(is.null(xlim)) xlim <- c(0,1)
# If the input is a data.frame bring it into matrix format
if(is.data.frame(X)) X <- as.matrix(X)
# Logical, was the input 'X' originally a matrix?
wasMatrix <- is.matrix(X) && dim(X)[1]>1
# In case of vector input, transform it into a matrix with one row
if(!is.matrix(X)) X <- t(as.matrix(X))
# check if we have for each row in the input 'X' a different color
if(length(lCol)!=dim(X)[1])
{
warning("Too less colors given, colors have automatically be assigned!")
lCol <- 1:dim(X)[1]
}
# Get the required p-value grid (in order to calcualte the obs/exp plot - we pool the p-values over all given vectors, in order to calcualte the most dense one)
temp <- unique(c(0,sort(as.vector(X)),1))
sigTests <- matrix(NA,ncol=length(temp),nrow=dim(X)[1]+1)
sigTests[1,] <- temp
for(i in 1:length(temp))
{
temp <- (X <= sigTests[1,i])
for(j in 2:dim(sigTests)[1])
{
sigTests[j,i] <- sum(temp[j-1,])
}
}
for(i in 2:dim(sigTests)[1])
{
sigTests[i,] <- sigTests[i,]/dim(X)[2]
}
ratios <- matrix(NA,ncol=ncol(sigTests),nrow=nrow(sigTests)-1)
distances <- matrix(NA,ncol=ncol(sigTests),nrow=nrow(sigTests)-1)
for(i in 1:dim(ratios)[1])
{
ratios[i,] <- sigTests[i+1,]/sigTests[1,]
distances[i,] <- sigTests[i+1,] - sigTests[1,]
}
ratios[1] <- 0
ylim <- c(0,max(sigTests[,sum(sigTests[1,]<=xlim[2])]))
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
if(crit=="NULL"){
par(oma=c(2,0,1,0),
mar=c(4,4,1,1))
} else {
nf <- layout(matrix(c(1,2),ncol=1),c(4,4), c(3,1), TRUE)
par(oma=c(2,0,1,0),
mar=c(4,4,1,0))
}
if (crit=="NULL"){
plot(c(0,sigTests[1,]),c(0,sigTests[2,]),col=lCol[1],type="l",xlab="Expected Ratio",ylab="Observed Ratio",xlim=xlim,ylim=ylim)
} else {
plot(c(0,sigTests[1,]),c(0,sigTests[2,]),col=lCol[1],type="l",ylab="Observed Ratio",xlab=" ",xlim=xlim,ylim=ylim)
}
lines(c(-10,10),c(-10,10),type="l")
if(wasMatrix)
{
for(i in 3:dim(sigTests)[1])
{
lines(c(0,sigTests[1,]),c(0,sigTests[i,]),type="l",col=lCol[i-1])
}
}
# Add the different rejection lines, if requested:
if(!is.null(rejLine))
{
if(rejLine=="bh"|rejLine=="simes")
{
for(i in 1:length(alpha)){
lines(c(0,1),c(0,(1+alpha[i])/alpha[i]),lty=rej.lty, col=rej.col)
}
} else if(rejLine=="bonferroni") {
for(i in 1:length(alpha)){
lines(c(0,1),c(0,dim(X)[2]/alpha[i]),lty=rej.lty, col=rej.col)
}
}
}
# Add a possible Legend
if(!is.null(incLegend)){
# legend(incLegend,legend=paste("Row",lCol), col=lCol, pch=20, lty=1,bg="white")#, trace=TRUE)
legend(incLegend,legend=rownames(X), col=lCol, pch=20, lty=1,bg="white")#, trace=TRUE)
}
if(crit=="ratio")
{
plot(c(-1,2),c(1,1),type="l",xlim=xlim,ylim=c(min(ratios),max(ratios)),xlab="Expected Ratio",ylab="Ratio")
for(i in 1:dim(ratios)[1])
{
lines(sigTests[1,],ratios[i,],col=lCol[i])
}
} else if(crit=="distance") {
plot(c(-1,2),c(1,1),type="l",xlim=xlim,ylim=c(min(distances),max(distances)),xlab="Expected Ratio",ylab="Distance")
for(i in 1:dim(ratios)[1])
{
lines(sigTests[1,],distances[i,],col=lCol[i])
}
lines(c(-10,10),c(0,0),lty="dotted")
}
}
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.