<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-kernel-heaping-package-iii/"></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-kernel-heaping-package-iii.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 the second part of this blog series, I showed how to compute spatial kernel density estimates based on area-level data. The Kernelheaping package also supports boundary-corrected kernel density estimation, which allows us to exclude certain areas, where we know that the density must be zero. One example is estimating the population density where we like to exclude uninhabited areas such as lakes, forests, parks etc. The Kernelheaping package employs a boundary correction method, where each single kernel is restricted to the area of interest. We continue with our example of elderly people in Berlin from part two:</p>
</div>
<pre><code><pre>library(maptools) library(dplyr) library(fields) library(ggplot2) library(RColorBrewer) library(Kernelheaping) library(rgeos) library(rgdal)</pre></code></pre>
<div>
<p>Again, we load a shapefile with the administrative districts, available from: <a rel="nofollow" href="https://www.statistik-berlin-brandenburg.de/opendata/RBS_OD_LOR_2015_12.zip" target="_blank">https://www.statistik-berlin-brandenburg.de/opendata/RBS_OD_LOR_2015_12.zip</a></p>
</div>
<pre><code><pre>data &lt;- read.csv2(&quot;EWR201512E_Matrix.csv&quot;) berlin &lt;- readOGR(&quot;RBS_OD_LOR_2015_12/RBS_OD_LOR_2015_12.shp&quot;) berlin &lt;- spTransform(berlin, CRS(&quot;+proj=longlat +datum=WGS84&quot;))</pre></code></pre>
<div>
<p>We load an OpenStreetMap file including shapes or polygons with information on uninhabited areas such as lakes, rivers, forests and parks: <a rel="nofollow" href="https://daten.berlin.de/datensaetze/openstreetmap-daten-f%C3%BCr-berlin" target="_blank">https://daten.berlin.de/datensaetze/openstreetmap-daten-für-berlin</a></p>
</div>
<pre><code><pre>berlinN &lt;- readOGR(&quot;berlin-latest-free.shp/gis_osm_landuse_a_free_1.shp&quot;) # land  berlinWater &lt;- readOGR(&quot;berlin-latest-free.shp/gis_osm_water_a_free_1.shp&quot;) # water</pre></code></pre>
<div>
<p>We specifically exlude residential areas and split the shapefile into the two remaining categories (“Nature” and “Other”):</p>
</div>
<pre><code><pre>table(berlinN@data$fclass) berlinN &lt;- berlinN[!(berlinN@data$fclass == &quot;residential&quot;), ] berlinGreen &lt;- berlinN[(berlinN@data$fclass %in%  	c(&quot;forest&quot;, &quot;grass&quot;, &quot;nature_reserve&quot;, &quot;park&quot;,    	  &quot;cemetery&quot;, &quot;allotments&quot;, &quot;farm&quot;, &quot;meadow&quot;,    	  &quot;orchard&quot;, &quot;vineyard&quot;, &quot;heath&quot;)), ] berlinOther &lt;- berlinN[!(berlinN@data$fclass %in%  	c(&quot;forest&quot;, &quot;grass&quot;, &quot;nature_reserve&quot;, &quot;park&quot;,  	  &quot;cemetery&quot;, &quot;allotments&quot;, &quot;farm&quot;, &quot;meadow&quot;,  	  &quot;orchard&quot;, &quot;vineyard&quot;, &quot;heath&quot;)), ]</pre></code></pre>
<div>
<p>These shapes are very complicated with many polygons. Thus we simplify them with the <span>gSimplify()</span> function from the <span>rgeos</span> package:</p>
</div>
<pre><code><pre>berlinGreen &lt;- spTransform(gSimplify(berlinGreen, tol = 0.0005, topologyPreserve = FALSE), 	                   CRS(&quot;+proj=longlat +datum=WGS84&quot;)) 	 berlinOther &lt;- spTransform(gSimplify(berlinOther, tol = 0.0005, topologyPreserve = FALSE),                            CRS(&quot;+proj=longlat +datum=WGS84&quot;)) 	 berlinWater &lt;- spTransform(gSimplify(berlinWater, tol = 0.0005, topologyPreserve = FALSE),                            CRS(&quot;+proj=longlat +datum=WGS84&quot;))</pre></code></pre>
<div>
<p>For the <span>dshapebivr()</span> and <span>dshapebivrProp()</span> functions we need a single shapefile; therefore we have to unite the water, nature and other shapefiles:</p>
</div>
<pre><code><pre>berlinUnInhabitated &lt;- gUnion(gSimplify(gUnion(berlinGreen, berlinOther), tol = 0.0005), 	                      berlinWater)</pre></code></pre>
<div>
<p>Now we perform the same data preparation steps as in the previous part and estimate the boundary-corrected density of people between 65 and 80 in Berlin. The shapefile of uninhabited areas now goes into the <span>deleteshapes</span> argument:</p>
</div>
<pre><code><pre> dataIn &lt;- cbind(do.call(rbind, lapply(berlin@polygons, function(x) x@labpt)),                    data$E_E65U80) est &lt;- dshapebivr(data = dataIn,                    burnin = 5,                    samples = 15,                    adaptive = FALSE,                    shapefile = berlin,                    deleteShapes = berlinUnInhabitated,                    gridsize = 325,                    boundary = TRUE)</pre></code></pre>
<div>
<p>To plot the map in <span>ggplot2</span>, we need to perform some additional data preparation steps:</p>
</div>
<pre><code><pre>berlin@data$id &lt;- as.character(berlin@data$PLR) berlin@data$E_E65U80 &lt;- data$E_E65U80 berlinPoints &lt;- fortify(berlin, region = &quot;id&quot;) berlin@data$E_E65U80density &lt;- berlin@data$E_E65U80 / (gArea(berlin, byid = TRUE) / 1000000) berlinDf &lt;- left_join(berlinPoints, berlin@data, by = &quot;id&quot;) kData &lt;- data.frame(expand.grid(long = est$Mestimates$eval.points[[1]],                                 lat = est$Mestimates$eval.points[[2]]),                     Density = est$Mestimates$estimate %>% as.vector) %>%    filter(Density > 0)</pre></code></pre>
<div>
<p>Now, we are able to plot the density together with the administrative districts and uninhabited areas of different types:</p>
</div>
<pre><code><pre>ggplot(kData) +   geom_raster(aes(long, lat, fill = Density)) +    ggtitle(&quot;Bivariate density of Inhabitants between 65 and 80 years&quot;) +   scale_fill_gradientn(colours = c(&quot;#FFFFFF&quot;, &quot;coral1&quot;))+   geom_polygon(fill = &quot;grey20&quot;, data = fortify(gIntersection(berlin, berlinOther)),                aes(long, lat, group = group), alpha = 0.25) +   geom_polygon(fill = &quot;darkolivegreen3&quot;, data = fortify(gIntersection(berlin, berlinGreen)),                aes(long, lat, group = group), alpha = 0.25) +   geom_polygon(fill = &quot;deepskyblue3&quot;, data = fortify(gIntersection(berlin, berlinWater)),                aes(long, lat, group = group), alpha = 0.25) +   geom_path(color = &quot;#000000&quot;, data = berlinDf, size = 0.1,             aes(long, lat, group = group)) +   coord_quickmap()</pre></code></pre>
<div>
<p><img alt="" src="https://i0.wp.com/www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_3_BivariateDensity.png?w=450&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<h2>Smooth Estimates of Proportion</h2>
<p>One may not only estimate the density, but also the proportion of a certain group relative to the overall population. The Kernelheaping package provides the <span>dshapebivrProp()</span> function which smoothly estimates the spatial proportion using a Nadaraya-Watson-type estimator. Naturally, it includes boundary correction as well. We use another open data example for Berlin on inhabitants with migration background from <a rel="nofollow" href="https://daten.berlin.de/datensaetze/einwohnerinnen-und-einwohner-mit-migrationshintergrund-berlin-lor-planungsr%C3%A4umen" target="_blank">https://daten.berlin.de/datensaetze/einwohnerinnen-und-einwohner-mit-migrationshintergrund-berlin-lor-planungsräumen</a></p>
<p>First, we load the dataset and merge the area ids such that they fit with the shapefile of Berlin:</p>
</div>
<pre><code><pre>berlinMigration &lt;- read.csv2(&quot;EWRMIGRA201512H_Matrix.csv&quot;) berlinMigration$RAUMID &lt;- as.character(berlinMigration$RAUMID) berlinMigration$RAUMID[nchar(berlinMigration$RAUMID) == 7] &lt;-   paste0(&quot;0&quot;, berlinMigration$RAUMID[nchar(berlinMigration$RAUMID) == 7]) berlinMigration &lt;- berlinMigration[order(berlinMigration$RAUMID), ] </pre></code></pre>
<div>
<p>We model the spatial proportion of inhabitants with Turkish migration background. For the proportion, a fourth column with the total number of people in that area is necessary:</p>
</div>
<pre><code><pre> dataTurk &lt;- cbind(do.call(rbind, lapply(berlin@polygons, function(x) x@labpt)),                    berlinMigration$HK_Turk,                    berlinMigration$MH_E)</pre></code></pre>
<div>
<p>We estimate the proportion with the <span>dshapebivrProp()</span> function now:</p>
</div>
<pre><code><pre>estTurk &lt;- dshapebivrProp(data = dataTurk,                            burnin = 5,                            samples = 10,                            adaptive = FALSE,                            deleteShapes = berlinUnInhabitated,                                            shapefile = berlin,                            gridsize = 325,                            boundary = TRUE,                            numChains = 4,                            numThreads = 4)</pre></code></pre>
<div>
<p>Now we can plot these proportions:</p>
</div>
<pre><code><pre>gridBerlin &lt;- expand.grid(long = estTurk$Mestimates$eval.points[[1]],                           lat = estTurk$Mestimates$eval.points[[2]]) kDataTurk &lt;- data.frame(gridBerlin,                          Proportion = estTurk$proportions %>% as.vector) %>%    filter(Proportion > 0) ggplot(kDataTurk) +    geom_raster(aes(long, lat, fill = Proportion)) +    ggtitle(&quot;Proportion of inhabitants with turkish migration background &quot;) +    scale_fill_gradientn(colours = c(&quot;#FFFFFF&quot;, &quot;coral1&quot;)) +    geom_polygon(fill = &quot;grey20&quot;, data = fortify(gIntersection(berlin, berlinOther)),                aes(long, lat, group = group), alpha = 0.25) +    geom_polygon(fill = &quot;darkolivegreen3&quot;, data = fortify(gIntersection(berlin, berlinGreen)),                aes(long, lat, group = group), alpha = 0.25) +    geom_polygon(fill = &quot;deepskyblue3&quot;, data = fortify(gIntersection(berlin, berlinWater)),                aes(long, lat, group = group), alpha = 0.25) +    geom_path(color = &quot;#000000&quot;, data = berlinDf, size = 0.1,                aes(long, lat, group = group)) +     coord_quickmap()</pre></code></pre>
<div>
<p><img alt="" src="https://i2.wp.com/www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_3_PropTurkInhab.png?w=450&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<h2>Hotspot Estimation</h2>
<p>Spatial kernel density estimates are a great tool to identify subpopulation hotspots. Three different countries / regions of origin are compared: Arabian countries, countries of the former Soviet Union and Poland. We perform the usual data preparation and estimation steps:</p>
</div>
<pre><code><pre> dataArab &lt;- cbind(do.call(rbind, lapply(berlin@polygons, function(x) x@labpt)),                    berlinMigration$HK_Arab) dataSU &lt;- cbind(do.call(rbind, lapply(berlin@polygons, function(x) x@labpt)),                  berlinMigration$HK_EheSU) dataPol &lt;- cbind(do.call(rbind, lapply(berlin@polygons, function(x) x@labpt)),                   berlinMigration$HK_Polen) estArab &lt;- dshapebivr(data = dataArab, burnin = 5, samples = 10, adaptive = FALSE,                       shapefile = berlin, gridsize = 325, boundary = TRUE) estSU &lt;- dshapebivr(data = dataSU, burnin = 5, samples = 10, adaptive = FALSE,                     shapefile = berlin, gridsize = 325, boundary = TRUE) estPol &lt;- dshapebivr(data = dataPol, burnin = 5, samples = 10, adaptive = FALSE,                      shapefile = berlin, gridsize = 325, boundary = TRUE) gridBerlin &lt;- expand.grid(long = estArab$Mestimates$eval.points[[1]],                           lat = estArab$Mestimates$eval.points[[2]])</pre></code></pre>
<div>
<p>Now we use the 97.5% quantile of the inhabited area to define hotspots:</p>
</div>
<pre><code><pre>kDataArab &lt;- data.frame(gridBerlin,                          Density = estArab$Mestimates$estimate %>% as.vector) %>%     filter(Density > 0) %>%   filter(Density > quantile(Density, 0.975)) %>%    mutate(Density = &quot;Arabian countries&quot;) kDataSU &lt;- data.frame(gridBerlin,                        Density = estSU$Mestimates$estimate %>% as.vector) %>%    filter(Density > 0) %>%   filter(Density > quantile(Density, 0.975)) %>%    mutate(Density = &quot;Former Soviet Union&quot;) kDataPol &lt;- data.frame(gridBerlin,                         Density = estPol$Mestimates$estimate %>% as.vector) %>%     filter(Density > 0) %>%   filter(Density > quantile(Density, 0.975)) %>%    mutate(Density = &quot;Poland&quot;) </pre></code></pre>
<div>
<p>Now, we display the hotspots of all three population subgroups in a single plot:</p>
</div>
<pre><code><pre>ggplot() +   geom_raster(aes(long, lat), fill = &quot;#FFFFFF&quot;, data = kData, alpha = 0.6) +    geom_raster(aes(long, lat, fill = Density), data = kDataArab, alpha = 0.6) +    geom_raster(aes(long, lat, fill = Density), data = kDataSU, alpha = 0.6) +    geom_raster(aes(long, lat, fill = Density), data = kDataPol, alpha = 0.6) +   scale_fill_manual(guide_legend(title = &quot;&quot;), values = c(&quot;#f8eb4a&quot;, &quot;#DD9123&quot;, &quot;#8A3B89&quot;)) +   ggtitle(&quot;Hotspots of Inhabitants With Different Migration Background&quot;) +   geom_polygon(fill = &quot;grey20&quot;, data = fortify(gIntersection(berlin, berlinOther)),                aes(long, lat, group = group), alpha = 0.25) +   geom_polygon(fill = &quot;darkolivegreen3&quot;, data = fortify(gIntersection(berlin, berlinGreen)),                aes(long, lat, group = group), alpha = 0.25) +   geom_polygon(fill = &quot;deepskyblue3&quot;, data = fortify(gIntersection(berlin, berlinWater)),                aes(long, lat, group = group), alpha = 0.25) +   geom_path(color = &quot;#000000&quot;, data = berlinDf, size = 0.1,                aes(long, lat, group = group)) +   coord_quickmap() +   theme(legend.position = &quot;top&quot;)</pre></code></pre>
<div>
<p><img alt="" src="https://i0.wp.com/www.inwt-statistics.com/files/INWT/images_blog/Kernelheaping/Kernelheaping_3_subPopHotspots.png?w=450&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<p>Further parts of the article series <strong>Introducing the Kernelheaping Package</strong>:</p>
<ul>
<li><a rel="nofollow" href="https://www.inwt-statistics.com/read-blog/introducing-the-kernelheaping-package-512.html" target="_blank">Part 1: Univariate Kernel Density Estimation for Heaped Data</a></li>
<li><a rel="nofollow" href="https://www.inwt-statistics.com/read-blog/introducing-the-kernelheaping-package-ii.html" target="_blank">Part 2: Bivariate Kernel Density Estimation for Heaped Data</a></li>
<li><strong>Part 3: Boundary Corrected Kernel Density Estimation</strong></li>
</ul></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-kernel-heaping-package-iii.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/7gphpF_JBQ8" height="1"></a></div>