An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110920/c1534f70/attachment.pl>
kernel density estimation - scaling factor problem?
2 messages · Myles Falconer, Rolf Turner
(1) You make an awfully complicated issue of the calculations.
Why not just:
require(spatstat)
W <- owin(c(0,10),c(0,10))
X <- ppp(x=c(2,7,4,5,8,8,1,3,3,9),
y=c(6,5,1,1,4,5,7,1,9,8),
marks=c(10,7,5,9,12,14,8,11,7,10),
window=W)
ddd <- density(X,sigma=1.4,weights=marks(X),edge=FALSE)
Doing range(ddd$v) gives
[1] 0.008501562 2.454624379
so the highest density is about 2.45.
(2) Why do you say that 2.4 is too small? You have a total of 93
apples over 100 square metres, giving a mean density of 0.93.
So 2.4 as a maximum density seems commensurate with this.
(3) The integral of ddd is a bit on the small side; 77.01386. This
seems to be due to your peculiar insistence on no edge correction.
(4) Your value of sigma may be a tad large; if we let density() choose
sigma for itself we get 1.25.
(5) Doing
eee <- density(X,weights=marks(X),diggle=TRUE)
yields an integral of exactly 93 and a maximal density of about 3.128.
Does that make you happier?
cheers,
Rolf Turner
On 21/09/11 09:43, Myles Falconer wrote:
Hi all, I'm a little new to R, so any help is appreciated. I'm having a problem interpreting the output from my kernel density plot. The density numbers seem smaller than I would expect based on my observed data. Here's an example dataset. Let's pretend I counted apples on apple trees in an orchard. Let's say the units are in metres. library(maptools) library(gstat) #make observation window x1<-seq(0,10, by=1) y1<-seq(0,10, by=1) df<-data.frame(x=rep(x1, each=length(y1)),y=rep(y1, length(x1)), z=1) gridded(df)<-~x+y df<-as(df, "SpatialGridDataFrame") dfowin<-as(df,"owin") #make a dataset to fit in the owin. rx<-c(2,7,4,5,8,8,1,3,3,9) ry<-c(6,5,1,1,4,5,7,1,9,8) apples<-c(10,7,5,9,12,14,8,11,7,10) rdf<-data.frame(x=rx,y=ry,apples=count) coordinates(rdf)<-~x+y #as.ppp rdf1<-as(rdf["count"],"ppp") rdf2<-ppp(x=rdf1$x,y=rdf1$y,marks=rdf1$marks,window=dfowin) plot(rdf2) plot(density(rdf2, sigma=1.4, weights=rdf2$marks, edge=FALSE)) So the highest density in the orchard is only 2.4 apples per metre... that doesn't seem right... Does it?
___________________________________________ Myles Falconer, Ontario Region Project Biologist, Bird Studies Canada, 115 Front St., P.O. Box 160, Port Rowan, ON N0E 1M0 Canada Email: mfalconer at birdscanada.org Ph: 1-888-448-2473 ext. 165 Local Ph: 519-586-3531 ext. 165 BSC website: http://www.bsc-eoc.org/