[R] Selectively shading areas under two density curves

[Jim Lemon]

Ian Watson wrote:
> Dear Listers,
> 
> I am drawing a plot of two density curves, for male and female incomes. I would 
> like to shade/hatch/color (whatever) the areas under the curves which are 
> distinctive for each gender. This is the code I have tried so far:
> 
> 
> m <- density(topmal.d$y, bw = "sj")
> f <- density(topfem.d$y, bw = "sj")
> par(mfrow = c(1,1))
> plot(x = c(0,400), y = c(0,0.02), type = "n",
>      bty = "l", xlab = "Annual earnings (in $thousands)", ylab = "Density")
> polygon(m, col = "blue", border = "blue", density = 20, angle = -45)
> polygon(f, col = "red", border = "red", density = 20, angle = 45)
> 
> 
> Without the data I realise you may not be able to see my goal. But basically the 
>   m density and the f density overlap for most of the plot, but areas "bulge" 
> out on the left and right of the overlapped area. What I'm wanting to do is 
> shade/hatch/color etc these areas which are unique to each density, that is, 
> which are outside the overlapped area.
> 
> My code is successful at hatching the bulging areas, but leaves me with a double 
> hatched area for the overlap, which is distracting. If I could turn this double 
> hatched area white, that would achieve my goal, though ideally I would like to 
> be able to specify something like "shade only areas under the m density curve 
> which are not also under the f density curve (and shade only under the f density 
> curve those areas not under the m density curve".
> 
Hi Ian,
Once you get used to the vagaries of "density" output, this isn't too 
hard, as you only have one crossing. I think this gets you what you want:

minc<-rnorm(1000,mean=65000,sd=15000)
finc<-rnorm(1000,mean=55000,sd=15000)
m <- density(minc, bw = "sj")
f <- density(finc, bw = "sj")
plot(m$x,m$y,main="Male vs female income distribution",
  xlab="",ylab="Density",type="l",xaxt="n")
library(plotrix)
axis.mult(1,at=c(0,20000,40000,60000,80000,100000,120000),mult=1000,
  labels=c("0","20","40","60","80","100","120"),mult.label="Annual 
income/1000")
lines(f$x,f$y)
fontop.end<-max(which(f$y>m$y))
xend<-length(f$x)
polygon(c(f$x[1:fontop.end],rev(m$x[1:fontop.end])),
  c(f$y[1:fontop.end],rev(m$y[1:fontop.end])),col="red")
fontop.end<-fontop.end+1
polygon(c(f$x[fontop.end:xend],rev(m$x[fontop.end:xend])),
  c(f$y[fontop.end:xend],rev(m$y[fontop.end:xend])),col="blue")

Jim