[R] kernel density estimation for univariate data using splancs

Kate Zinszer kate.zinszer at mail.mcgill.ca
Tue Aug 18 23:37:18 CEST 2009


Hi,

I previously received help in extract data from a shapefile and now my question is about kernel density estimation.  My objective is to have 3 kernel density plots; 2 for the each set of cases and the 3rd is the difference in kernel densities between the 2 sets of cases.  Previously, I used the density function from the stats package, which worked but I wanted finer control of the bandwidth.  Therefore, I was advised to switch to the package splancs.   I've been trying to copy the example provided in the help file but cannot produce the same results for my data.  I've pasted my code below and I am not sure how to define the grd1 (cellsize and cells.dim).  Here is my code:  

coords <- S at polygons[[1]]@Polygons[[1]]@coords 
coord <- as.points(coords) 
cases1.xy <- as.points(mergedis$jit.x.x, mergedis$pre.fsa.shp.Y2)
grd1 <- GridTopology(cellcentre.offset=c(50281627.2, 580082), cellsize=c(0.2, 0.2), cells.dim=c(75,100))
k1000 <- spkernel2d(cases1.xy, coord, h0=1000, grd1)
k5000 <- spkernel2d(cases1.xy, coord, h0=5000, grd1)
if (.sp_lt_0.9()) {
  df <- AttributeList(list(k1000=k1000, k5000=k5000))
} else {
  df <- data.frame(k1000=k1000, k5000=k5000)
}
kernels <- SpatialGridDataFrame(grd1, data=df)
spplot(kernels, checkEmptyRC=FALSE, col.regions=terrain.colors(16), cuts=15)

Using this code, I get NAs for all values of k1000 and k5000.  If this did work, can I then subtract the spplots from one another to have the difference between the two sets of cases?

Any advice would be greatly appreciated! 
Kate



. You need to check the @hole
slot of the ... at Polygons[[1]] object if you care whether it's a hole
or an island!

Barry



More information about the R-help mailing list