Unexpected result with weights in spatstat::density.psp()

Question

I am trying to use kernel density smoothing to map the intensity of possible pest escape from vehicle traffic. Each route has been broken down into straight lines with each line having an integer attribute for the number of times the segment was travelled upon. However, when I use this attribute as the weight in kernel density smoothing, the weights don't seem to be used.

I've created a simplified reprex below with two abutting straight lines. Can anyone explain to me how I can make density.psp() account for the fact that one segment has an attribute 2x the magnitude of the other?

Many thanks for your help,

Josh

    # Load packages 
    library(spatstat)
    
    # Create the data frame. Note that coordinates are projected to the
    # North America Lambert Conformal Conic system. https://epsg.io/102009
    my_lines <- data.frame(
      fx = c(1252365.22479882, 1233600.0015391),
      fy = c(510853.438463705, 626675.859171899),
      tx = c(1233600.0015391, 1218256.03484937),
      ty = c(626675.859171899, 721347.256108354),
      attrib = c(100, 50)
    ) 
    
    # Creat the observation window
    my_owin <- owin(
      xrange = c(1200000, 1270000),
      yrange = c(500000,  730000)
    ) 
    
    # Create the psp object with 'attrib' as the mark. 
    my_psp <- psp(
      x0 = my_lines$fx, y0 = my_lines$fy,
      x1 = my_lines$tx, y1 = my_lines$ty,
      window = my_owin, marks = my_lines$attrib
    )
    
    # Create the KDS image with 'attrib' as the weights
    kdi_image_weighted <-
      density.psp(
        my_psp,
        kernel = "gaussian",
        sigma = 4000,
        weights = my_psp$attrib,
        edge = FALSE
      )
    
    # Create the KDS image without weights
    kdi_image_unweighted <-
      density.psp(
        my_psp,
        kernel = "gaussian",
        sigma = 4000,
        edge = FALSE
      )
    
    # Plot the weighted and unweighted KDS images. Note that they are the same 
    # despite one being weighted. 
    plot(kdi_image_weighted) 
    plot(kdi_image_unweighted)

Hi Josh Persi, it seems that `weights = my_psp$num_trips` in the `kdi_image_weighted <- density.psp(...)` code does not refer to anything. Maybe, try with `weights = my_psp$marks`. Cheers — lovalery, Feb 07 '22 at 18:29
Hi @lovalery, thanks for catching that! My code was supposed to read `weights = my_psp$attrib` after I renamed the `num_trips` column to make this reprex more general. However, switching `weights = my_psp$marks` works and provides the expected output! Thank you so much! — Josh Persi, Feb 07 '22 at 18:39

score 3 · Accepted Answer · answered Feb 08 '22 at 02:57

Short answer:

Use marks() to extract the mark values from an object in the spatstat package. Example:

 Z <- density(my_psp, weights=marks(my_psp))

Long answer:

This problem is not related to the function density.psp.

Firstly there is a misunderstanding about names of objects. In your example, an object my_psp is created as psp(.... marks=my_lines$attrib). Then the subsequent code expects my_psp to contain an element my_psp$attrib. This is not how evaluation works in R.

When you type the command psp(...., marks=my_lines$attrib), the expression my_lines$attrib will be evaluated (by extracting the element called attrib from the object my_lines). After evaluating these data, the system forgets where it got the data from; the name attrib is forgotten. The resulting data (that was stored in my_lines$attrib) will be passed to the function psp as the argument marks.

Secondly you are using the $ operator to extract data from an object (in this case an object of class "psp" in the spatstat package). This is OK for simple objects such as lists, but for objects created by packages it is unsafe, because the internal structure could change when the package is updated.

Your example code expects the object my_psp to contain an element with a particular name, but there's no reason why it should. The function psp() creates an object of class "psp" which has some weird internal structure, described in help(psp.object).

To extract the mark values from an object in the spatstat package, it is currently possible to use $marks, but this could change, and it is strongly recommended to use the function marks().

Unexpected result with weights in spatstat::density.psp()

1 Answers1