spatstat / spatstat.linnet Goto Github PK

Sub-package of spatstat containing code for linear networks

R 84.36% Makefile 0.24% C 15.40%

network-analysis spatial-data-analysis point-processes statistical-models statistical-inference statistical-analysis kernel-density-estimation density-estimation heat-equation

spatstat.linnet's Introduction

spatstat.linnet

Spatial analysis on a linear network, for the spatstat family

The original spatstat package has been split into several sub-packages (See spatstat/spatstat).

This package spatstat.linnet is one of the sub-packages. It contains the subset of the functionality of spatstat that deals with data on linear networks.

There is also an extension package spatstat.Knet which contains additional algorithms for linear networks.

Where to find data

Examples of datasets on linear networks are the point patterns chicago, dendrite and spiders provided in the spatstat.data package (available when spatstat.linnet is loaded) and the point pattern wacrashes provided in the extension package spatstat.Knet (which must be loaded separately).

Overview of `spatstat.linnet`

spatstat.linnet supports

Network geometry

creation of linear networks from coordinate data
extraction of networks from tessellations
modification of networks
interactive editing of networks
geometrical operations and measurement on networks
construction of the disc in the shortest-path metric
trees, tree branch labels, tree pruning

Point patterns on a network

creation of point patterns on a network from coordinate data
extraction of sub-patterns
shortest-path distance measurement

Covariates on a network

create pixel images and functions on a network
arithmetic operators for pixel images on a network
plot pixel images on a network (colour/thickness/perspective)
tessellation on a network

Simulation

completely random (uniform Poisson) point patterns on a network
nonuniform random (Poisson) point patterns on a network
Switzer-type point process
log-Gaussian Cox process

Exploratory analysis of point patterns on a network

kernel density estimation on a network
bandwidth selection
kernel smoothing on a network
estimation of intensity as a function of a covariate
ROC curves
Berman-Waller-Lawson test
CDF test
variable selection by Sufficient Dimension Reduction
K function on a network (shortest path or Euclidean distance)
pair correlation function on a network (shortest path or Euclidean distance)
inhomogeneous K function and pair correlation function
inhomogeneous F, G and J functions
simulation envelopes of summary functions

Parametric modelling and inference on a network

fit point process model on a network
fitted/predicted intensity
analysis of deviance for point process model
simulate fitted model

spatstat.linnet's People

Contributors

Stargazers

Watchers

Forkers

agila5 fairmiles

spatstat.linnet's Issues

Possible bug in `density.lpp` for networks with very short segments

This is an issue transferred from spatstat (spatstat/spatstat#133) because of the package splitup.

Several users have reported erroneous results from density.lpp/densityHeat in which the highest values of density occur at locations far away from the nearest data point. I do not yet have a minimal working example. The bug seems to arise when the network contains very short segments (e.g. shorter than 1 metre in a road network that is 10-100 km wide). This is being investigated.

Possible typo in the docs of `rhohat.lpp`

Dear developers, the documentation of rhohat.lpp says the following:

spatstat.linnet/man/rhohat.lpp.Rd

Lines 211 to 216 in c04b9d6

 This technique assumes that \code{covariate} has continuous values. 

 It is not applicable to covariates with categorical (factor) values 

 or discrete values such as small integers. 

 For a categorical covariate, use 

 \code{\link{intensity.quadratcount}} applied to the result of 

 \code{\link{quadratcount}(X, tess=covariate)}.

Maybe I don't understand something, but I think that the doc should be slightly updated since intensity.quadratcount or quadratcount cannot be run with lpp objects:

library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 2.2-2
#> Loading required package: spatstat.core
#> Loading required package: nlme
#> Loading required package: rpart
#> spatstat.core 2.3-0
#> Loading required package: spatstat.linnet
#> spatstat.linnet 2.3-0
#> 
#> spatstat 2.2-0       (nickname: 'That's not important right now') 
#> For an introduction to spatstat, type 'beginner'
Y <- runiflpp(30, simplenet)
quadratcount(Y)
#> Error in UseMethod("quadratcount"): no applicable method for 'quadratcount' applied to an object of class "c('lpp', 'ppx')"

^{Created on 2021-08-07 by the reprex package (v2.0.0)}

Bug in bw.relrisklpp (maybe) with duplicated points

Dear spatstat group, I'm sorry to bother you again. I'm writing this issue to report a problematic behaviour in bw.relrisklpp that occur with duplicated points. Reproducible example:

# packages
library(spatstat)
#> Carico il pacchetto richiesto: spatstat.data
#> Carico il pacchetto richiesto: spatstat.geom
#> spatstat.geom 2.0-0
#> Carico il pacchetto richiesto: spatstat.core
#> Carico il pacchetto richiesto: nlme
#> Carico il pacchetto richiesto: rpart
#> spatstat.core 2.0-0
#> Carico il pacchetto richiesto: spatstat.linnet
#> spatstat.linnet 2.0-0
#> 
#> spatstat 2.0-1       (nickname: 'Caution: contains small parts') 
#> For an introduction to spatstat, type 'beginner'

# simulate data
set.seed(1)
my_ppp1 <- ppp(
  x = c(0.2, 0.2, runif(5)), 
  y = c(0.2, 0.2, runif(5))
)
#> Warning: data contain duplicated points
my_lpp1 <- lpp(
  X = my_ppp1, 
  L = simplenet
)

# Estimate intensity at "points" following Rakshit et al. (2019)
density(my_lpp1, sigma = 0.2, distance = "euclidean", at = "points")
#> [1] 1.8888483 1.2261747 1.6113216 1.4250548 0.8453113 4.0378911
#> attr(,"sigma")
#> [1] 0.2

# It returns 6 estimates, i.e. one for each unique pair of coordinates. The
# problem is that I think that implies the following code fails:
my_ppp2 <- ppp(
  x = runif(5), 
  y = runif(5)
)
my_lpp2 <- lpp(
  X = my_ppp2, 
  L = simplenet
)

my_lpp <- superimpose(A = my_lpp1, B = my_lpp2)

# Fails
bw.relrisklpp(my_lpp, distance = "euclidean")
#> Error in `[<-.data.frame`(`*tmp*`, , j, value = structure(c(13.9870593985448, : replacement has 6 rows, data has 7

^{Created on 2021-03-18 by the reprex package (v1.0.0)}

I run the previous code in Rstudio with the debugger, and I think that the problem occurs here:

https://github.com/baddstats/spatstat.linnet/blob/2bc676e44d5db6ccc378a3183387ad58eb7457f3/R/relrisk.lpp.R#L62-L79

since the lapply() function returns a list with fewer observations than points (since there are duplicated coordinates in one of the spatial inputs). Unfortunately, I don't understand the internals of spatstat, so I have no idea how to fix that. Moreover, I'm also not extremely familiar with estimating relative risk on a road network, so I'm not sure if I should change something in case of duplicated points.

EDIT

I simply fixed the error on the real dataset by jittering the points by a small margin, but maybe you still want to add a warning message.

EDITv2

Unfortunately, that's not true since, even if the points are not identical after jittering, they are projected to the same point of the network (one of the edges near the boundary of the owin).

Bug in density.lpp with marks, duplicated points and weights

Dear spatstat developers, thank you very much for your work. I recently updated to spatstat v2.0.1, and I think I found a bug in density.lpp. Reprex:

# packages
library(spatstat)
#> Carico il pacchetto richiesto: spatstat.data
#> Carico il pacchetto richiesto: spatstat.geom
#> Registered S3 method overwritten by 'spatstat.geom':
#>   method     from
#>   print.boxx cli
#> spatstat.geom 1.65-5
#> Carico il pacchetto richiesto: spatstat.core
#> Carico il pacchetto richiesto: nlme
#> Carico il pacchetto richiesto: rpart
#> spatstat.core 1.65-5
#> Carico il pacchetto richiesto: spatstat.linnet
#> 
#> spatstat 2.0-1       (nickname: 'Caution: contains small parts') 
#> For an introduction to spatstat, type 'beginner'

# simulate ppp data in [0, 1] x [0, 1] 
my_ppp <- ppp(
  x = c(0.2, 0.2, seq(0.2, 0.9, length.out = 10)),
  y = c(0.2, 0.2, seq(0.9, 0.2, length.out = 10)), 
  marks = letters[1:12]
)

# simulate lpp data
my_lpp <- lpp(
  X = my_ppp, 
  L = simplenet
)

# plot
plot(my_lpp, use.marks = FALSE)

# density estimation
density(unmark(my_lpp), sigma = bw.scott.iso, distance = "euclidean", weights = runif(12))
#> Warning in m * weights: longer object length is not a multiple of shorter object
#> length
#> Error in pixellate.ppp(x, ..., padzero = TRUE): length(weights) == npoints(x) || length(weights) == 1 is not TRUE

^{Created on 2021-03-15 by the reprex package (v1.0.0)}

Considering the warning message, I think that the error is related to the following code:

https://github.com/baddstats/spatstat.linnet/blob/99f77051a76c787f7fc5ce189a42fc4db0bb2b22/R/density.lpp.R#L23-L27

In that code, length(m) is equal to the number of unique points in X (i.e. 11), while length(weights) is equal to the number of points in X (i.e. 12).

Session info

sessioninfo::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value                       
#>  version  R version 4.0.4 (2021-02-15)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  Italian_Italy.1252          
#>  ctype    Italian_Italy.1252          
#>  tz       Europe/Berlin               
#>  date     2021-03-15                  
#> 
#> - Packages -------------------------------------------------------------------
#>  package         * version date       lib source        
#>  abind             1.4-5   2016-07-21 [1] CRAN (R 4.0.3)
#>  assertthat        0.2.1   2019-03-21 [1] CRAN (R 4.0.4)
#>  backports         1.2.1   2020-12-09 [1] CRAN (R 4.0.3)
#>  cli               2.3.1   2021-02-23 [1] CRAN (R 4.0.4)
#>  crayon            1.4.1   2021-02-08 [1] CRAN (R 4.0.4)
#>  curl              4.3     2019-12-02 [1] CRAN (R 4.0.4)
#>  deldir            0.2-10  2021-02-16 [1] CRAN (R 4.0.4)
#>  digest            0.6.27  2020-10-24 [1] CRAN (R 4.0.4)
#>  ellipsis          0.3.1   2020-05-15 [1] CRAN (R 4.0.4)
#>  evaluate          0.14    2019-05-28 [1] CRAN (R 4.0.4)
#>  fansi             0.4.2   2021-01-15 [1] CRAN (R 4.0.4)
#>  fs                1.5.0   2020-07-31 [1] CRAN (R 4.0.4)
#>  glue              1.4.2   2020-08-27 [1] CRAN (R 4.0.4)
#>  goftest           1.2-2   2019-12-02 [1] CRAN (R 4.0.3)
#>  highr             0.8     2019-03-20 [1] CRAN (R 4.0.4)
#>  htmltools         0.5.1.1 2021-01-22 [1] CRAN (R 4.0.4)
#>  httr              1.4.2   2020-07-20 [1] CRAN (R 4.0.4)
#>  knitr             1.31    2021-01-27 [1] CRAN (R 4.0.4)
#>  lattice           0.20-41 2020-04-02 [2] CRAN (R 4.0.4)
#>  lifecycle         1.0.0   2021-02-15 [1] CRAN (R 4.0.4)
#>  magrittr          2.0.1   2020-11-17 [1] CRAN (R 4.0.4)
#>  Matrix            1.3-2   2021-01-06 [2] CRAN (R 4.0.4)
#>  mgcv              1.8-33  2020-08-27 [2] CRAN (R 4.0.4)
#>  mime              0.10    2021-02-13 [1] CRAN (R 4.0.4)
#>  nlme            * 3.1-152 2021-02-04 [2] CRAN (R 4.0.4)
#>  pillar            1.5.1   2021-03-05 [1] CRAN (R 4.0.4)
#>  pkgconfig         2.0.3   2019-09-22 [1] CRAN (R 4.0.4)
#>  polyclip          1.10-0  2019-03-14 [1] CRAN (R 4.0.3)
#>  purrr             0.3.4   2020-04-17 [1] CRAN (R 4.0.4)
#>  R6                2.5.0   2020-10-28 [1] CRAN (R 4.0.4)
#>  reprex            1.0.0   2021-01-27 [1] CRAN (R 4.0.4)
#>  rlang             0.4.10  2020-12-30 [1] CRAN (R 4.0.4)
#>  rmarkdown         2.7     2021-02-19 [1] CRAN (R 4.0.4)
#>  rpart           * 4.1-15  2019-04-12 [2] CRAN (R 4.0.4)
#>  sessioninfo       1.1.1   2018-11-05 [1] CRAN (R 4.0.4)
#>  spatstat        * 2.0-1   2021-03-13 [1] CRAN (R 4.0.4)
#>  spatstat.core   * 1.65-5  2021-02-01 [1] CRAN (R 4.0.4)
#>  spatstat.data   * 2.0-0   2021-02-10 [1] CRAN (R 4.0.4)
#>  spatstat.geom   * 1.65-5  2021-01-15 [1] CRAN (R 4.0.4)
#>  spatstat.linnet * 1.65-3  2021-02-05 [1] CRAN (R 4.0.4)
#>  spatstat.sparse   1.2-1   2021-02-06 [1] CRAN (R 4.0.4)
#>  spatstat.utils    2.0-0   2021-01-24 [1] CRAN (R 4.0.4)
#>  stringi           1.5.3   2020-09-09 [1] CRAN (R 4.0.3)
#>  stringr           1.4.0   2019-02-10 [1] CRAN (R 4.0.4)
#>  styler            1.3.2   2020-02-23 [1] CRAN (R 4.0.4)
#>  tensor            1.5     2012-05-05 [1] CRAN (R 4.0.3)
#>  tibble            3.1.0   2021-02-25 [1] CRAN (R 4.0.4)
#>  utf8              1.2.1   2021-03-12 [1] CRAN (R 4.0.4)
#>  vctrs             0.3.6   2020-12-17 [1] CRAN (R 4.0.4)
#>  withr             2.4.1   2021-01-26 [1] CRAN (R 4.0.4)
#>  xfun              0.22    2021-03-11 [1] CRAN (R 4.0.4)
#>  xml2              1.3.2   2020-04-23 [1] CRAN (R 4.0.4)
#>  yaml              2.2.1   2020-02-01 [1] CRAN (R 4.0.3)
#> 
#> [1] C:/Users/Utente/Documents/R/win-library/4.0
#> [2] C:/Program Files/R/R-4.0.4/library

Bug in arithmetic for images on a network

If Z is a pixel image on a linear network (class "linim") then Z/3 is a pixel image in two dimensions (class "im") instead of a pixel image on the network, as intended. This is a bug in Math.linim.

The bug has been fixed in the development version 2.2-0.002.

In the meantime, use eval.linim(Z/3) to get the correct behaviour.

Problem with the inclusion of lintess objects in lppm

Dear all, I think I've found a possible problem with the inclusion of lintess objects in a model estimated by lppm. For example:

library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 3.2-2
#> Loading required package: spatstat.random
#> spatstat.random 3.1-5
#> Loading required package: spatstat.explore
#> Loading required package: nlme
#> spatstat.explore 3.2-1
#> Loading required package: spatstat.model
#> Loading required package: rpart
#> spatstat.model 3.2-4
#> Loading required package: spatstat.linnet
#> spatstat.linnet 3.1-1
#> 
#> spatstat 3.0-6 
#> For an introduction to spatstat, type 'beginner'

toy_points <- rpoislpp(20, simplenet)
toy_tess <- lintess(
  L = simplenet, 
  df = data.frame(
    seg = 1:10, 
    t0 = 0, 
    t1 = 1, 
    tile = c("A", "B")
  )
)

lppm(toy_points ~ toy_tess)
#> Error in mpl.get.covariates(covariates, P, Pname, covfunargs): Each entry in the list 'covariates' should be an image, a function, a window, a tessellation or a single number

A possible workaround I've found:

toy_linfun <- as.linfun(toy_tess)
lppm(toy_points ~ toy_linfun)
#> Point process model on linear network
#>  Fitted to point pattern dataset 'toy_points'
#> 
#> Nonstationary Poisson process
#> 
#> Log intensity:  ~toy_linfun
#> 
#> Fitted trend coefficients:
#> (Intercept) toy_linfunB 
#>    3.311431    0.123307 
#> 
#>             Estimate      S.E.    CI95.lo  CI95.hi Ztest       Zval
#> (Intercept) 3.311431 0.1601282  2.9975856 3.625276   *** 20.6798800
#> toy_linfunB 0.123307 0.2176698 -0.3033181 0.549932        0.5664863
#> 
#> Domain: Linear network with 10 vertices and 10 lines
#> Enclosing window: rectangle = [0, 1] x [0, 1] units

^{Created on 2023-07-12 with reprex v2.0.2}

Am I doing something wrong in the first call to lppm? Is the workaround ok? In any case, thank you very much for maintaining this amazing set of packages.

How to change resolution of Zimage in rhohat.lpp

Dear developers, I'm not sure if the following reasoning does make sense, but I was wondering if it's possible to modify the code behind evalCovar.lppm to adjust the resolution of the object Zimage when creating a linim object starting from a linfun object. In particular, I noticed that, at the moment, the argument eps is ignored in evalCovar.lppm, and I was wondering if it can be used when creating the linim object named Zimage.

I tried to implement that idea in a separate fork of this repository and this is the result:

library(spatstat.linnet)
set.seed(1)
Y <- runiflpp(50, domain(chicago))
rhoY <- rhohat(Y, "y", eps = 5)
plot(predict(rhoY))

rhoY <- rhohat(Y, "y", eps = 20)
plot(predict(rhoY))

^{Created on 2021-08-24 by the reprex package (v2.0.0)}

(the current version returns a 128 x 128 image in both cases).

I'm not sure the resolution of Zimage can be modified like that, and I'm not sure I can "reuse" the argument eps (since it defines the grid used to evaluate the covariate). In any case, if you agree, I could complete the implementation and create a pull request.

Sorry for all these issues.

	This technique assumes that \code{covariate} has continuous values.
	It is not applicable to covariates with categorical (factor) values
	or discrete values such as small integers.
	For a categorical covariate, use
	\code{\link{intensity.quadratcount}} applied to the result of
	\code{\link{quadratcount}(X, tess=covariate)}.