Features:
Fast Kernel Density calculation using hexagonal grid and plotting of result.
Edge correction including Jones-Diggllle algorithm as described in Jones, M.C. (1993) Simple boundary corrections for kernel density estimation. Statistics and Computing 3, 135--146.
Default bandwidth is the diagonal normal scale bandwidth
Using bei spatial dataset from spatstat.explore
library(spatstat.data)
data=bei
Calculating kernel density using hexagonal grid
#specify the x, y vectors
density = hexDensity(x=data$x, y=data$y)
#or just let hexDensity figure it out.
density = hexDensity(data)
Comparing to density.ppp by spatstat which use square-grid. (eps is to ensure square instead of rectangular grid)
library(ks)
bandwidth = sqrt(diag(Hns.diag(cbind(data$x,data$y)))) #Set bandwidth to be the same default plug-in bandwidth in hexDensity.
library(spatstat.explore)
#eps variable is used to turn the grid square instead of rectangle
density = density.ppp(data,sigma=bandwidth, eps=diff(range(data$x))/128)
plot.im(density,col=colorRampPalette(viridis::viridis(11)))
Comparison to SpatialKDE package, which can also do hexagonal kernel density but really slow to compute and plot. Selected "bandwidth" and "cell size" values are chosen to best fit with the above examples but may not match perfectly. Note that SpatialKDE does not have option for different bandwidth values in different directions and does not have edge correction.
library(SpatialKDE)
library(dplyr)
library(sp)
library(sf)
library(tmap)
#Prepare data
data <- data.frame(bei) %>%
st_as_sf(coords = c("x", "y"), dim = "XY") %>%
st_set_crs(28992) %>%
select()
cell_size <- 12
band_width <- 160
#Create grid
grid <- data %>%
create_grid_hexagonal(cell_size = cell_size, side_offset = band_width)
#Calculate KDE
kde <- data %>%
kde(band_width = band_width, kernel = "quartic", grid = grid)
#Plot
tm_shape(kde) +
tm_polygons(col = "kde_value",style="cont", palette = "viridis", title = "KDE Estimate",legend.show=FALSE)
Using spatial transcriptomic dataset. This example use MERFISH mouse hypothalamic preoptic region dataset from MerfishData package on bioconductor. Finding density for cells of 'inhibitory' cell class.
library(MerfishData)
spe = MouseHypothalamusMoffitt2018()
#filter for just Inhibitory cells on z-layer -0.14.
cdat = data.frame(colData(spe),spatialCoords(spe))
cdat = subset(cdat, cell_class!= "Ambiguous",select = -c(cell_id,sample_id,sex,behavior,neuron_cluster_id))
cdat = subset(cdat,z == -0.14)
cdat.inhibitory = subset(cdat, cell_class == ("Inhibitory"))
#hexDensity
density.inhibitory = hexDensity(cdat.inhibitory$x,cdat.inhibitory$y)
plotHexDensity(density.inhibitory)
#comparison with density.ppp by spatstat
#need to convert into ppp object
cdat.inhibitory = ppp(cdat.inhibitory$x,cdat.inhibitory$y,window = owin(range(cdat.inhibitory$x),range(cdat.inhibitory$y)))
bandwidth = sqrt(diag(Hns.diag(cbind(cdat.inhibitory$x,cdat.inhibitory$y))))
density.inhibitory = density.ppp(cdat.inhibitory,sigma=bandwidth, eps=diff(range(cdat.inhibitory$x))/128)
plot.im(density.inhibitory, col=colorRampPalette(viridis::viridis(11)))
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4
