R help - Mar 2007 - Density estimation graphs

Dear all,

I'm struggling with a plot and would value any help!

I'm attempting to highlight a histogram and density plot to show a
proportion of cases above a threshold value. I wanted to cross-hatch the
area below the density curve. The breaks and bandwidth are deliberate
integer values because of the type of data I'm looking at.

I've managed to do this, but I don't think it is very good! It would be
difficult, for example, to do a cross-hatch using this technique.

allele.plot <- function(x, threshold=NULL, hatch.col='black',
hatch.border=hatch.col, lwd=par('lwd'),...) {
	h <- hist(x, breaks=max(x), plot=F)
	d <- density(x, bw=1)
	plot(d, lwd=lwd, ...)
	
	if (!is.null(threshold)) {
		d.t <- d$x>threshold
		d.x <- d$x[d.t]
		d.y <- d$y[d.t]
		d.l <- length(d.x)
		# draw all but first line of hatch
		for (i in 2:d.l) {
		lines(c(d.x[i],d.x[i]),c(0,d.y[i]),
			col=hatch.col,lwd=1)
		}
		# draw first line in hatch border colour
		lines(c(d.x[1],d.x[1]),c(0,d.y[1]),
			col=hatch.border,lwd=lwd)

		# and now re-draw density plot lines
		lines(d, lwd=lwd)
	}
}

# some pretend data
s8 = rnorm(100, 15, 5)
threshold = 19			# an arbitrary cut-off
allele.plot(s8, threshold,
hatch.col='grey',hatch.border='black')


Is there a better way? As always, I'm sure there's a one-liner rather
than my crude technique!

Best wishes,

Mark
-- 
Dr. Mark Wardle
Clinical research fellow and specialist registrar, Neurology
University Hospital Wales and Cardiff University, UK

Mark Wardle wrote:
> Dear all,
> 
> I'm struggling with a plot and would value any help!
> ...
> 
> Is there a better way? As always, I'm sure there's a one-liner
rather
> than my crude technique!
> 
As always, I've spent ages trying to sort this, and then the minute
after sending an email, I find the polygon() function.

Ignore previous message!

Best wishes,

Mark

On Mar 15, 2007, at 12:37 PM, Mark Wardle wrote:
> Dear all,
>
> I'm struggling with a plot and would value any help!
>
> I'm attempting to highlight a histogram and density plot to show a
> proportion of cases above a threshold value. I wanted to cross- 
> hatch the
> area below the density curve. The breaks and bandwidth are deliberate
> integer values because of the type of data I'm looking at.
>
> I've managed to do this, but I don't think it is very good! It  
> would be
> difficult, for example, to do a cross-hatch using this technique.
Don't know about a cross-hatch, but in general I use "polygon" for
highlighting areas like that:

allele.plot <- function(x, threshold=NULL, hatch.col='black',
hatch.border=hatch.col, lwd=par('lwd'),...) {
	h <- hist(x, breaks=max(x), plot=F)
	d <- density(x, bw=1)
	plot(d, lwd=lwd, ...)	
	if (!is.null(threshold)) {
		d.t <- d$x>threshold
		d.x <- d$x[d.t]
		d.y <- d$y[d.t]
		polygon(c(d.x[1],d.x,d.x[1]),c(0,d.y,0), col=hatch.col,lwd=1)
	}
}
# some pretend data
s8 = rnorm(100, 15, 5)
threshold = 19			# an arbitrary cut-off
allele.plot(s8, threshold,
hatch.col='grey',hatch.border='black')


Perhaps this can help a bit. Btw, what was d.l for?

Haris Skiadas
Department of Mathematics and Computer Science
Hanover College