{"id":407,"date":"2014-01-01T22:10:41","date_gmt":"2014-01-02T03:10:41","guid":{"rendered":"http:\/\/homepages.uc.edu\/~yaozo\/wordpress\/?p=407"},"modified":"2014-01-01T22:10:41","modified_gmt":"2014-01-02T03:10:41","slug":"shading-a-kernel-density-plot-between-two-points","status":"publish","type":"post","link":"https:\/\/zhuoyao.net\/index.php\/2014\/01\/01\/shading-a-kernel-density-plot-between-two-points\/","title":{"rendered":"Shading a kernel density plot between two points."},"content":{"rendered":"
\n\n\n\n\n
\n
\n

I frequently use kernel density plots to illustrate distributions. These are easy and fast to create in R like so:<\/p>\n

set.seed(1)\ndraws <- rnorm(100)^2\ndens <- density(draws)\nplot(dens)\n#or in one line like this: plot(density(rnorm(100)^2))<\/code><\/pre>\n
Which gives me this nice little PDF:<\/p>\n
\"It's<\/p>\n
I’d like to shade the area under the PDF from the 75th to 95th percentiles. It’s easy to calculate the points using the\u00a0quantile<\/code>\u00a0function:<\/p>\n
q75 <- quantile(draws, .75)\nq95 <- quantile(draws, .95)<\/code><\/pre>\n
But how do I shade the the area between\u00a0q75<\/code>\u00a0and\u00a0q95<\/code>?<\/p>\n<\/div>\n
<\/td>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n
<\/a><\/p>\n
\n\n\n\n\n
<\/td>\n\n
\n

With the\u00a0polygon()<\/code>\u00a0function, see its help page and I believe we had similar questions here too.<\/p>\n

You need to find the index of the quantile values to get the actual\u00a0(x,y)<\/code>\u00a0pairs.<\/p>\n

Edit:<\/em>\u00a0Here you go:<\/p>\n

x1 <- min(which(dens$x >= q75))  \nx2 <- max(which(dens$x <  q95))\nwith(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=\"gray\"))<\/code><\/pre>\n
Output (added by JDL)<\/p>\n
\"alt<\/p>\n<\/div>\n<\/td>\n<\/tr>\n
<\/td>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n
<\/div>\n

<\/a><\/p>\n

\n\n\n\n\n
<\/td>\n\n
\n

Another solution:<\/p>\n

dd <- with(dens,data.frame(x,y))\nlibrary(ggplot2)\nqplot(x,y,data=dd,geom=\"line\")+\n  geom_ribbon(data=subset(dd,x>q75 & x<q95),aes(ymax=y),ymin=0,\n              fill=\"red\",colour=NA,alpha=0.5)<\/code><\/pre>\n
Result:\u00a0\"alt<\/p>\n<\/div>\n<\/td>\n<\/tr>\n
<\/td>\n\n
\n\n\n\n\n
<\/td>\n\n
\n

An expanded solution:<\/p>\n

If you wanted to shade both tails (copy & paste of Dirk’s code) and use known x values:<\/p>\n

set.seed(1)\ndraws <- rnorm(100)^2\ndens <- density(draws)\nplot(dens)\n\nq2     <- 2\nq65    <- 6.5\nqn08   <- -0.8\nqn02   <- -0.2\n\nx1 <- min(which(dens$x >= q2))  \nx2 <- max(which(dens$x <  q65))\nx3 <- min(which(dens$x >= qn08))  \nx4 <- max(which(dens$x <  qn02))\n\nwith(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=\"gray\"))\nwith(dens, polygon(x=c(x[c(x3,x3:x4,x4)]), y= c(0, y[x3:x4], 0), col=\"gray\"))<\/code><\/pre>\n
Result:<\/p>\n
\"2-tailed<\/p>\n<\/div>\n<\/td>\n<\/tr>\n
<\/td>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n

<\/a><\/p>\n

\n\n\n\n
<\/td>\n\n
\n

This question needs a\u00a0lattice<\/code>\u00a0answer. Here’s a very basic one, simply adapting the method employed by Dirk and others:<\/p>\n

#Set up the data\nset.seed(1)\ndraws <- rnorm(100)^2\ndens <- density(draws)\n\n#Put in a simple data frame   \nd <- data.frame(x = dens$x, y = dens$y)\n\n#Define a custom panel function;\n# Options like color don't need to be hard coded    \nshadePanel <- function(x,y,shadeLims){\n    panel.lines(x,y)\n    m1 <- min(which(x >= shadeLims[1]))\n    m2 <- max(which(x <= shadeLims[2]))\n    tmp <- data.frame(x1 = x[c(m1,m1:m2,m2)], y1 = c(0,y[m1:m2],0))\n    panel.polygon(tmp$x1,tmp$y1,col = \"blue\")\n}\n\n#Plot\nxyplot(y~x,data = d, panel = shadePanel, shadeLims = c(1,3))<\/code><\/pre>\n
\"enter<\/p>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
I frequently use kernel density plots to illustrate distributions. These are easy and fast to create in R like so: set.seed(1) draws <- rnorm(100)^2 dens… <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[20],"tags":[],"class_list":["post-407","post","type-post","status-publish","format-standard","hentry","category-r"],"_links":{"self":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/407","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/comments?post=407"}],"version-history":[{"count":0,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/407\/revisions"}],"wp:attachment":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/media?parent=407"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/categories?post=407"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/tags?post=407"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}