Point pattern classification with spatstat: how to choose the right bandwidth?

Question

I'm still trying to find the best way to classify bivariate point patterns:

Point pattern classification with spatstat: what am I doing wrong?

I now analysed 110 samples of my dataset using @Adrian's suggestion with sigma=bw.diggle (as I wanted an automatic bandwidth selection). f is a "resource selection function" (RSF) which describes the relationship between the intensity of the Cancer point process and the covariate (here kernel density of Immune):

Cancer <- split(cells)[["tumor"]]
Immune <- split(cells)[["bcell"]]
Dimmune <- density(Immune,sigma=bw.diggle)
f <- rhohat(Cancer, Dimmune)

I am in doubt about some results I've got. A dozen of rho-functions looked weird (disrupted, single peak). After changing to default sigma=NULL or sigma=bw.scott (which are smoother) the functions became "better" (see examples below). I also experimented with the following manipulations:

  cells # bivariate point pattern with marks "tumor" and "bcell"
  o.marks<-cells$marks  # original marks
    #A) randomly re-assign original marks
  a.marks <- sample(cells$marks)
    #B) replace marks randomly with a 50/50 proportion
  b.marks<-as.factor(sample(c("tumor","bcell"), replace=TRUE, size=length(o.marks))) 
    #C) random (homogenious?) pattern with the original number of points 
  randt<-runifpoint(npoints(subset(cells,marks=="tumor")),win=cells$window) 
  randb<-runifpoint(npoints(subset(cells,marks=="bcell")),win=cells$window)
  cells<-superimpose(tumor=randt,bcell=randb)
    #D) tumor points are associated with bcell points (is "clustered" a right term?)
  Cancer<-rpoint(npoints(subset(cells,marks=="tumor")),Dimmune,win=cells$window)
    #E) tumor points are segregated from bcell points
  reversedD<-Dimmune
  density.scale.v<-sort(unique((as.vector(Dimmune$v)[!is.na(as.vector(Dimmune$v))]))) # density scale
  density.scale.v.rev<-rev(density.scale.v)# reversed density scale
  new.image.v<-Dimmune$v
  # Loop over matrix
  for(row in 1:nrow(Dimmune$v)) {
    for(col in 1:ncol(Dimmune$v)) {
      if (is.na(Dimmune$v[row, col])==TRUE){next}
      number<-which(density.scale.v==Dimmune$v[row, col])
      new.image.v[row, col]<-density.scale.v.rev[number]}
  }
  reversedD$v<-new.image.v # reversed density
  Cancer<-rpoint(npoints(subset(cells,marks=="tumor")),reversedD,win=cells$window)

A better way to generate inverse density heatmaps is given by @Adrian in his post below.

I could not generate rpoint patterns for the bw.diggle density as it produced negative numbers.Thus I replaced the negatives Dimmune$v[which(Dimmune$v<0)]<-0 and could run rpoint then. As @Adrian explained in the post below, this is normal and can be solved easier by using a density.ppp option positive=TRUE.

I first used bw.diggle, because hopskel.test indicarted "clustering" for all my patterns. Now I'm going to use bw.scott for my analysis but can this decision be somehow justified? Is there a better method besides "RSF-function is looking weird"?

some examples:

sample10: sample20: sample110:

score 1 · Accepted Answer · answered Jul 16 '21 at 09:04

That is a lot of questions!

Please try to ask only one question per post.

But here are some answers to your technical questions about spatstat.

Negative values: The help for density.ppp explains that small negative values can occur because of numerical effects. To force the density values to be non-negative, use the argument positive=TRUE in the call to density.ppp. For example density(Immune, bw.diggle, positive=TRUE).

Reversed image: to reverse the ordering of values in an image Z you can use the following code:

V <- Z
A <- order(Z[])
V[][A] <- Z[][rev(A)]

Then V is the order-reversed image.

Thank you @Adrian for the technical advice! Sorry, I obviously did not read the help carefully enough. The code for reversing is much neater as my bulky ones. — Sergej S, Jul 16 '21 at 09:56

score 1 · Answer 2 · answered Jul 17 '21 at 01:11

Tips for your code:

to generate a random point pattern with the same number of points and in the same window as an existing point pattern X, use Y <- runifpoint(ex=X).
To extract the marks of a point pattern X, use a <- marks(X). To assign new marks to a point pattern X, use marks(X) <- b.
to randomly permute the marks attached to the points in a point pattern X, use Y <- rlabel(X).
to assign new marks to a point pattern X where the new marks are drawn randomly-with-replacement from a given vector of values m, use Y <- rlabel(X, m, permute=FALSE).

Point pattern classification with spatstat: how to choose the right bandwidth?

2 Answers2