<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/introducing-the-kernelheaping-package/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.inwt-statistics.com/read-blog/introducing-the-kernelheaping-package-512.html"> INWT-Blog-RBloggers</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p>In this blog article I’d like to introduce the univariate kernel density estimation for heaped (i.e. rounded or interval censored) data with the Kernelheaping package.</p>
<p>It is not unusual to have interval censored data such as in income surveys due to anonymisation or simplification issues. However, a simple task like plotting a density may fail completely for rounded or interval censored data. When using class centers as input, the density estimate might get bumpy or spiky around the center points. This problem gets even worse for larger sample sizes with decreasing bandwidth. One might try to manually increase the bandwidth until one obtains a sufficiently smooth estimate. However, this will result in oversmoothing as we will see in our example. The <a rel="nofollow" href="https://cran.r-project.org/package=Kernelheaping" title="Kernelheaping package on CRAN" target="_blank">Kernelheaping package, available on CRAN</a>, delivers an algorithm to obtain nonparametric density estimates for interval censored data.</p>
</div>
<div>
<h2>Example: Kernel density estimation for interval data without correction</h2>
<p>In the following chunk I’ll simulate 5000 observations from a log-normal distribution. Then I plot the histogram together with a conventional kernel density estimate.</p>
</div>
<pre><pre><code>library(&quot;ggplot2&quot;) library(&quot;dplyr&quot;) </code></pre></pre>
<pre><pre><code>set.seed(0) dat &lt;- data.frame(x = rlnorm(5000, meanlog = 7.5, sdlog = 0.5)) ggplot(dat, aes(x, y = ..density..)) +   geom_histogram(position = &quot;identity&quot;, color = &quot;black&quot;, fill = &quot;white&quot;) +   geom_density(fill = &quot;blue&quot;, alpha = 0.7, linetype = 0) +   ylim(c(0, 0.00055)) +   xlim(c(0, 10000)) </code></pre></pre>
<div>
<p><img alt="" width="450" src="https://www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_plot_1.svg"></p>
</div>
<div>
<p>Now let’s split our data into 14 classes imitating a typical income survey.</p>
</div>
<pre><pre><code>classes &lt;- c(0, 500, 1000, 1500, 2000, 2500, 3000, 4000, 5000, 6000, 8000,              10000, 15000, Inf) dat$xclass &lt;- cut(dat$x, breaks = classes)</code></pre></pre>
<div>
<p>If I now use the class centers as input, the result is a very unrealistic, spiky density.</p>
</div>
<pre><pre><code>classCenters &lt;- (classes[-1] - classes[-length(classes)]) /   2 + classes[-length(classes)] classCenters[length(classCenters)] &lt;- 2 * classCenters[length(classCenters) - 1] levels(dat$xclass) &lt;- classCenters dat$xclassCenter &lt;- as.numeric(as.character(dat$xclass)) ggplot(dat, aes(xclassCenter, y = ..density..)) +   geom_histogram(breaks = classes, color = &quot;black&quot;, fill = &quot;white&quot;) +   geom_density(fill = &quot;blue&quot;, alpha = 0.7, linetype = 0) </code></pre></pre>
<div>
<p><img alt="" width="450" src="https://www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_plot_2.svg"></p>
</div>
<div>
<p>I’ve got the option to adjust the density by manipulating the adjust parameter. But getting a nice and smooth density shape requires manipulation to such an extent that we end up with serious oversmoothing.</p>
</div>
<pre><pre><code>ggplot(dat, aes(xclassCenter, y = ..density..)) +   geom_histogram(breaks = classes, color = &quot;black&quot;, fill = &quot;white&quot;) +   geom_density(adjust = 2.5, fill = &quot;blue&quot;, alpha = 0.7, linetype = 0) +   ylim(c(0, 0.00055)) +   xlim(c(0, 10000))</code></pre></pre>
<div>
<p><img alt="" width="450" src="https://www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_plot_3.svg"></p>
</div>
<div>
<h2>Kernel density estimation for interval data with the Kernelheaping package</h2>
<p>Now, we install and load the Kernelheaping package. Its <em>dclass </em>function implements an iterative SEM (Stochastic Expectation Maximization) algorithm.</p>
</div>
<pre><pre><code>install.packages(&quot;Kernelheaping&quot;) library(&quot;Kernelheaping&quot;)</code></pre></pre>
<div>
<p>It can be called with any bandwidth selector available to the <em>density </em>function in the stats package and also provides optional boundary correction for positive only data.</p>
</div>
<pre><pre><code>densityEst &lt;- dclass(xclass = dat$xclass, classes = classes, burnin = 50,                      samples = 100, evalpoints = 1000)</code></pre></pre>
<pre><pre><code>classes &lt;- data.frame(xclass = densityEst$xclass) estimates &lt;- data.frame(Mestimates = densityEst$Mestimates, gridx = densityEst$gridx) ggplot() +   geom_histogram(data = classes, aes(xclass, y = ..density..),                   breaks = densityEst$classes, color = &quot;black&quot;, fill = &quot;white&quot;) +   geom_area(data = estimates, aes(gridx, Mestimates),              fill = &quot;blue&quot;, alpha = 0.7, size = 1) +   ylim(c(0, 0.00055)) +   xlim(c(0, 10000)) </code></pre></pre>
<div>
<p><img alt="" width="450" src="https://www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_plot_4.svg"></p>
</div>
<div>
<p>The resulting density estimate is smooth and very close to the estimate on the original unclassed data. The package is also able to handle the case of survey participants rounding values by their own, each by a different level of coarseness. Please check out the <em>dheaping </em>function with given examples.</p>
<p>Soon there will be a second part to this article, where I’ll show how to compute bivariate densities given area-level data with the <a rel="nofollow" href="https://cran.r-project.org/web/packages/Kernelheaping/index.html" title="Kernelheaping Package on CRAN" target="_blank">Kernelheaping package</a>.</p>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://www.inwt-statistics.com/read-blog/introducing-the-kernelheaping-package-512.html"> INWT-Blog-RBloggers</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/kqmrp_cq1Yk" height="1"></a></div>