Overlapping density plots

density plots.r

### r code to support: http://econometricsense.blogspot.com/2011/01/flexibility-of-r-graphics.html

library(colorspace) # package for rainbow_hcl function
 
 
ds <- rbind(data.frame(dat=KyCropsAndSubsidies[,][,"LogAcres"], grp="All"),
            data.frame(dat=KyCropsAndSubsidies[,][KyCropsAndSubsidies$subsidy_in_millions > 2.76,"LogAcres"], grp=">median"),
            data.frame(dat=KyCropsAndSubsidies[,][KyCropsAndSubsidies$subsidy_in_millions <= 2.76,"LogAcres"], grp="<=median"))
 
 
# histogram and density for all ears
hs <- hist(ds[ds$grp=="All",1], main="", xlab="LogAcres", col="grey90", ylim=c(0, 25), breaks="fd", border=TRUE)
 
dens <- density(ds[ds$grp=="All",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[1])
 
# density for above median subsidies
dens <- density(ds[ds$grp==">median",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[2])
 
# density for below median subsidies
dens <- density(ds[ds$grp=="<=median",1], na.rm=TRUE)
rs <- max(hs$counts)/max(dens$y)
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(3)[3])
 
# Add a rug to illustrate density.
 
rug(ds[ds$grp==">median", 1], col=rainbow_hcl(3)[2])
rug(ds[ds$grp=="<=median", 1], col=rainbow_hcl(3)[3])
 
# Add a legend to the plot.
 
legend("topright", c("All", ">median", "<=media"), bty="n", fill=rainbow_hcl(3))
 
# Add a title to the plot. 
 
title(main="Distribution of Acres Planted by Subsidies Recieved Above or Below Median", sub=paste("Created Using R Statistical Package"))

# histogram and density estimates for all data 
hs <- hist(trade_by_yr$logTrade,main="", xlab="trade", col="grey90", ylim=c(0, 95), breaks="fd", border=TRUE)  # histogram 
 
dens <- density(trade_by_yr$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density 
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[1]) # plot densiy 
 
# density estimates for year 2000 trade data
 
y2000 <- trade_by_yr[trade_by_yr$year==2000,] # subset data for year
dens <- density(y2000$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density  
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[2]) # plot densiy 
 
# density estimates for year 2004 trade data
 
y2004 <- trade_by_yr[trade_by_yr$year==2004,]  # subset data for year 
dens <- density(y2004$logTrade)  # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density 
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[3]) # plot densiy 
 
# densty estimates for year 2008 trade data
 
y2008 <- trade_by_yr[trade_by_yr$year==2008,] # subset data for year
dens <- density(y2008$logTrade) # density 
rs <- max(hs$counts)/max(dens$y)  # rescale/mormalize density  
lines(dens$x, dens$y*rs, type="l", col=rainbow_hcl(4)[4]) # plot densiy 
 
# Add a legend to the plot.
 
legend("topright", c("All", "2000", "2004", "2008"), bty="n", fill=rainbow_hcl(4))
 
# Add a title to the plot. 
 
title(main="Distribution of Total World Trade Volume by Country by Year", sub=paste("Created Using R Statistical Package"))