Name:
Location: United States

Friday, June 30, 2006

Babies 1: qqplot with marginal histograms

# A function for drawing QQ-plot of 2 distributions along with marginal histograms

qq.hist = function(x=rnorm(1000), y=rnorm(1000), xlab="X", ylab="Y")
{

# calculate quantiles
xq = quantile(x, probs = seq(0,1,0.01))
yq = quantile(y, probs = seq(0,1,0.01))

# calculate histogram values
minval = floor(min(c(x,y)))-1
maxval = ceiling(max(c(x,y)))+1
xhist = hist(x, breaks=seq(minval,maxval,0.5), plot=FALSE)
yhist = hist(y, breaks=seq(minval,maxval,0.5), plot=FALSE)

top <- max(c(xhist$counts, yhist$counts))
xrange <- c(minval,maxval)
yrange <- c(minval,maxval)

# layout
nf <- layout(matrix(c(2,4,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
layout.show(nf)

# draw the qq-plot
par(mar=c(4,4,1,1))
qqplot(x, y, xlim=xrange, ylim=yrange, pch=".", cex=2, main="", xlab = xlab, ylab=ylab, cex.lab=1.5)
abline(h=0,v=0,col="blue")

# histogram for the x-axis
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes=F, ylim=c(0, top), space=0)


# histogram for the y-axis
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes=F, xlim=c(0, top), space=0, horiz=TRUE)


}