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
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
