[R-sig-Geo] kernel density estimation - scaling factor problem?

[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/