Fast food spatial analysis

WSTĘP

W literaturze występuje wieloletni trend badania wpływu bliskości punktów fast food do szkół na otyłość dzieci do nich uczęszczających, a także ich skłonnośc do korzystania z tych punktów. Większość badań wiążących lokację fast foodów z bliskością szkół wychodzi z tego punktu.

Jednym z głównych badań było badanie Simona i in. (2008) analizowało bliskość lokali fast food względem szkół publicznych w hrabstwie Los Angeles. Wykazano, że 23,3% szkół miało co najmniej jeden taki lokal w promieniu 400 metrów, a blisko 65% w promieniu 800 metrów. Restauracje te znajdowały się znacznie bliżej szkół średnich niż podstawowych czy gimnazjów. W obszarach komercyjnych szkoły w uboższych dzielnicach były znacznie bardziej narażone na sąsiedztwo fast foodów niż placówki w bogatszych rejonach. Autorzy sugerują, że otoczenie szkół może niweczyć wysiłki na rzecz poprawy żywienia uczniów, zwłaszcza w społecznościach o niskich dochodach.

Drugim badaniem analizującym bliskość szkół i fast foodów było badnie Kwate i Loh (2010) Skupiono się w nim na przestrzennym skupieniu fast foodów wokół szkół publicznych i prywatnych w Nowym Jorku. Ustalono, że publiczne szkoły średnie są najbardziej narażone na bliskość takich lokali, często przekraczającą dwukrotnie wartości wynikające z przypadku. Zaobserwowano też, że szkoły z większym odsetkiem czarnoskórych uczniów miały gęstszą sieć fast foodów w swojej okolicy. Wyniki te wskazują na nierówny dostęp do zdrowego otoczenia żywieniowego wokół nowojorskich placówek edukacyjnych.

W nieniejszej pracy skupimy się nie na szczegółowym rozpatrywnaniu różnych typów szkół, ale generalnym zjawisku. Sprawdzimy też czy istnieje podobny trend w bliskości fast foodów do szkół w Warszawie.

Biblioteki

library(osmdata)

  Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(tidyverse)

  ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
  ✔ dplyr     1.2.0     ✔ readr     2.2.0
  ✔ forcats   1.0.1     ✔ stringr   1.6.0
  ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
  ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
  ✔ purrr     1.2.1

  ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
  ✖ dplyr::filter() masks stats::filter()
  ✖ dplyr::lag()    masks stats::lag()
  ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggmap)

  ℹ Google's Terms of Service: <https://mapsplatform.google.com>
    Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service>
    OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles>
  ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

library(maptiles)
library(tidyterra)

  
  Attaching package: 'tidyterra'
  
  The following object is masked from 'package:stats':
  
      filter

library(sf)

  Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(spatstat)

  Loading required package: spatstat.data
  Loading required package: spatstat.univar
  spatstat.univar 3.1-7
  Loading required package: spatstat.geom
  spatstat.geom 3.7-2
  Loading required package: spatstat.random
  spatstat.random 3.4-5
  Loading required package: spatstat.explore
  Loading required package: nlme
  
  Attaching package: 'nlme'
  
  The following object is masked from 'package:dplyr':
  
      collapse
  
  spatstat.explore 3.8-0
  Loading required package: spatstat.model
  Loading required package: rpart
  spatstat.model 3.6-1
  Loading required package: spatstat.linnet
  spatstat.linnet 3.4-1
  
  spatstat 3.5-1 
  For an introduction to spatstat, type 'beginner'

library(naniar)
library(rgugik)
library(dixon)

  Loading required package: splancs
  Loading required package: sp
  
  Spatial Point Pattern Analysis Code in S-Plus
   
   Version 2 - Spatial and Space-Time analysis
  
  
  Attaching package: 'splancs'
  
  The following object is masked from 'package:dplyr':
  
      tribble
  
  The following object is masked from 'package:tidyr':
  
      tribble
  
  The following object is masked from 'package:tibble':
  
      tribble

Dane

Dane o punktach fast food i szkołach zostały pobrane z OpenStreetMap poprzez Overpass API. Obiekty oznaczone innaczej niż jako punkty zmieniono sprowadzono do ich centroidów i analizowano jako punkty. Obiekty, którym nie przypisano name, a w przypadku fast foodów — ani name, ani cuisine, zostały usunięte z analizy po wizualnej inspecji, ponieważ nie wnosiły wartościowych infromacji. Po ograniczeniu do granic administracyjnych Warszawy zbiór danych zawiera 1445 punkty fast food i 854 szkół. Wzorzec punktowy reprojektowano na Mercatora (EPSG:3857) i przeskalowano na kilometry używając powierzchni Warszawy równej 517.24 km².

Pobieramy fast foody.

# fast food
warsaw_bb <- getbb("Warsaw, Poland")
warsaw_ff <- warsaw_bb %>% 
  opq(timeout = 300) %>%
  add_osm_feature(key = "amenity", value = "fast_food") %>%
  osmdata_sf()

Obiekty zapisane w innym formacie niż punkty sprowadzamy do centroidów.

ff_poly <- warsaw_ff$osm_polygons %>% select(osm_id, name, cuisine, geometry)
ff_point <- warsaw_ff$osm_points %>% select(osm_id, name, cuisine, geometry)

ff_marks <- rbind(ff_poly, ff_point) %>%
  st_centroid() %>% 
  mutate(type = "Fast Food")
  Warning: st_centroid assumes attributes are constant over geometries
nrow(ff_marks)
  [1] 3094

Usuwamy obiekty, które nie mają ani nazwy, ani kuchni.

plot(ff_marks)

ff_marks <- ff_marks %>% filter(!is.na(name) | !is.na(cuisine))
nrow(ff_marks)

  [1] 1543

Pobieramy szkoły

# schools
warsaw_schools <- opq(getbb("Warsaw, Poland")) %>%
  add_osm_feature(key = "amenity", value = "school") %>%
  osmdata_sf()

Sprowadzamy niepunkty do centroidów.

school_point <- warsaw_schools$osm_points %>% select(osm_id, name, geometry)
school_lines <- warsaw_schools$osm_lines %>% select(osm_id, name, geometry)
school_poly <- warsaw_schools$osm_polygons %>% select(osm_id, name, geometry)
school_multipoly <- warsaw_schools$osm_multipolygons %>% select(osm_id, name, geometry)

school_marks <- bind_rows(school_point, school_lines, school_poly, school_multipoly) %>%
  st_centroid() %>%                          # lines / polygons / multipolygons → point
  mutate(type = "School")
  Warning: st_centroid assumes attributes are constant over geometries
nrow(school_marks)
  [1] 12315

Liczba szkół (12315) zdecydowanie przewyższa realną wartość (ok. 1 tys.). Prawdopobonie te same obiekty mapowano kilkakrotnie na różne sposoby, co potwierdza poniższy wykres - w wielu przypadkach wiele punktów nakłada się na siebie, a po przefiltrowaniu są reprezentowane pojedynczo:

plot(school_marks)

Liczba szkół jest też zbliżona do realnej (https://um.warszawa.pl/edukacja-2025).

school_marks <- school_marks %>% filter(!is.na(name))
nrow(school_marks)

  [1] 940

Łączyny dane o fast foodach i szkołach.

marks <- bind_rows(school_marks, ff_marks) %>%
  st_set_crs(4326)                           # make sure CRS is explicit

Ograniczamu do granic administracyjnych Warszawy.

warsaw_borders <- opq("Warsaw, Poland") %>%
  add_osm_feature(key = "boundary", value = "administrative") %>%
  add_osm_feature(key = "admin_level", value = "8") %>%
  osmdata_sf()
warsaw_borders <- warsaw_borders$osm_multipolygons[1,]

marks_waw <- st_filter(marks, warsaw_borders)

marks_waw %>% st_drop_geometry() %>% count(type)

         type    n
  1 Fast Food 1421
  2    School  854

Przekształcamy dane pod analizę wzorca punktowego.

# reproject to planar (Mercator
proj_marks_waw <- st_transform(marks_waw, crs = 3857)   # ← marks (schools + FF)

# Warsaw administrative boundary from OSM
warsaw_poly <- st_transform(warsaw_borders, crs = 3857)

# create owin
W <- as.owin(warsaw_poly)

# extract coordinates
coords <- st_coordinates(proj_marks_waw)

# MARKED ppp
warsaw_ppp <- ppp(
  x     = coords[, 1],
  y     = coords[, 2],
  window = W,
  marks = as.factor(proj_marks_waw$type)   # ← 'type' column: "School" / "Fast Food"
)
warsaw_ppp <- rjitter(warsaw_ppp, 0.2) #rozrzut

# rescale to km
rsc        <- sqrt(area.owin(W) / 517.24) # 517.24 km² = Warsaw area
warsaw_ppp <- rescale.ppp(warsaw_ppp, rsc, "km")

# separate point patterns for marks
warsaw_ppp_ff <- warsaw_ppp[warsaw_ppp$marks == "Fast Food"]
warsaw_ppp_ff <- warsaw_ppp_ff %>% unmark()
warsaw_ppp_sc <- warsaw_ppp[warsaw_ppp$marks == "School"]
warsaw_ppp_sc <- warsaw_ppp_sc %>% unmark()

WZORZEC PUNKTOWY

Spojrzenie na mapę ujawnia silną koncentrację fast foodów w centrum Warszawy z mniejszym zagęszczeniem w ościennych dzielnicach. Szkoły wydają się bardziej równomiernie rozłożone, co odzwierciedla ich rolę jako infrastuktry publicznej, która powinna odwzorowywać raczej zagęszczenie ludności niż komercyjną logikę. Skupimy się najpierw na rozkładzie fast foodów.

ggplot() +
  geom_sf(data = warsaw_borders, fill = "white") +
  geom_sf(data = marks_waw, aes(color = type), size = 0.4) +
  scale_color_manual(name = "Type", values = c("Fast Food" = "steelblue", "School" = "tomato")) +
  theme_minimal()

Unmarked

Analiza siatki kwadratów ujawnia, że intensywność fast foodów nie jest identyczna na całym analizowanym obszarze, a skupia się głównie w centrum. Potwierdza to formalny test równoliczności (X² = 2481.9, df = 21, p < 0.001). Wynik ten jest zbierzny z przypuszczeniami z obserwacji punktów.

par(
  mfrow    = c(1, 2),
  mar      = c(.5, .5, .5, .5), # marginesy
  pty      = "s",           # to robi, że wykresy są tej samej wielkości
  cex.main = 0.9
)

qdr <- quadratcount(warsaw_ppp_ff, nx = 5, ny = 5)
img_intensity <- intensity(qdr, image = TRUE)

plot(img_intensity, values0 = FALSE, ribbon = FALSE, main = "Intensity — count", col = NA)
plot(qdr, add = TRUE, col = "blue")
plot(img_intensity, ribbon = FALSE, main = "Intensity — heatmap")

par(mfrow = c(1, 1), pty = "m")

quadrat.test(warsaw_ppp_ff, nx = 5, ny = 5, alternative = "clustered")

  Warning: Some expected counts are small; chi^2 approximation may be inaccurate

  
    Chi-squared test of CSR using quadrat counts
  
  data:  warsaw_ppp_ff
  X2 = 2481.9, df = 21, p-value < 2.2e-16
  alternative hypothesis: clustered
  
  Quadrats: 22 tiles (irregular windows)

Sprawdzamy jak wygląda intensywność dla innych gęstości siatki.

par(
  mfrow = c(2,2),
  mar = c(.5, .5, .5, .5),   # marginesy
  cex.main = 0.9
)
plot(intensity(quadratcount(warsaw_ppp_ff, nx = 5, ny = 5), 
               image=TRUE), main = "Intensity 5x5")
plot(intensity(quadratcount(warsaw_ppp_ff, nx = 10, ny = 10), 
               image=TRUE), main = "Intensity 10x10")
plot(intensity(quadratcount(warsaw_ppp_ff, nx = 15, ny = 15), 
               image=TRUE), main = "Intensity 15x15")
plot(intensity(quadratcount(warsaw_ppp_ff, nx = 20, ny = 20), 
               image=TRUE), main = "Intensity 20x20")

par(mfrow = c(1, 1))

Estymator jądrowy gęstości (Gaussian, Scott’s bandwidth = 1.06 km) potwierdza centralne skupienie widoczene w siatce kwadratów.

par(
  mfrow = c(1,1),
  mar = c(.5, .5, .5, .5)   # marginesy
)
# bandwidth
bw_scott <- bw.scott(warsaw_ppp_ff)
# Kernel density
d_scott <- density(warsaw_ppp_ff, sigma = bw_scott[1], leaveoneout = FALSE, positive = TRUE)
plot(d_scott, main = "Fast food KDE", kernel = "gaussian")

## distances between fast foods
dist_fun <- distfun(warsaw_ppp_ff)
dist_map <- distmap(warsaw_ppp_ff)

par(mfrow = c(1, 2))
plot(dist_fun, main = "distfun – fast food")
plot(dist_map, main = "distmap – fast food")

par(mfrow = c(1, 1))

Wszystkie statystyki drugiego rzędu (funkcje K, G, F, J) wskazują na klastrowanie fast foofów. Obserwowane funkcje K i G leżą nad funkcją, którą generowałby proces Poissona (CSR) na wszystkich odległościach. Podobnie funkcja J jest przyjmuje wartości znacznie poniżej 1 dla wszystkich odległości, co wzkazuje na klastrowanie. Co ciekawe funkcja F w większości leży poniżej liini wskazującej na CSR, co wskazywałoby na klastrowanie, jednak dla przedziału od 0 do ok. 0.2km leży ona powyżej tej linii, co sugeruje, że w takiej skali występuje więcej pustej przestrzeni, niż wynikałoby to z losowego rozłożenia punktów.

all.ff <- allstats(warsaw_ppp_ff)
plot(all.ff, main = " ")

# Envelope
# Formal tests
mad <- mad.test(warsaw_ppp_ff,  Lest,   nsim = 49, use.theo = TRUE)

  Generating 49 simulated realisations of CSR  ...
  1, 2,  [7:52 remaining] 3,
   [7:21 remaining] 4,  [7:19 remaining] 5,  [7:04 remaining] 6,
   [6:57 remaining] 7,  [6:42 remaining] 8,  [6:42 remaining] 9,
   [6:31 remaining] 10,  [6:20 remaining] 11,  [6:08 remaining] 12,
   [5:57 remaining] 13,  [5:47 remaining] 14,  [5:38 remaining] 15,
   [5:28 remaining] 16,  [5:16 remaining] 17,  [5:06 remaining] 18,
   [4:57 remaining] 19,  [4:49 remaining] 20,  [4:38 remaining] 21,
   [4:31 remaining] 22,  [4:22 remaining] 23,  [4:12 remaining] 24,
   [4:03 remaining] 25,  [3:53 remaining] 26,  [3:44 remaining] 27,
   [3:34 remaining] 28,  [3:25 remaining] 29,  [3:15 remaining] 30,
   [3:05 remaining] 31,  [2:55 remaining] 32,  [2:46 remaining] 33,
   [2:36 remaining] 34,  [2:26 remaining] 35,  [2:16 remaining] 36,
   [2:06 remaining] 37,  [1:56 remaining] 38,  [1:47 remaining] 39,
   [1:37 remaining] 40,  [1:27 remaining] 41,  [1:17 remaining] 42,
   [1:08 remaining] 43,  [58 sec remaining] 44,  [48 sec remaining] 45,
   [39 sec remaining] 46,  [29 sec remaining] 47,  [19 sec remaining] 48,
   [10 sec remaining] 
  49.
  
  Done.

dclf <- dclf.test(warsaw_ppp_ff, Lest,   nsim = 49, use.theo = TRUE)

  Generating 49 simulated realisations of CSR  ...
  1, 2,  [7:23 remaining] 3,
   [7:28 remaining] 4,  [7:10 remaining] 5,  [7:01 remaining] 6,
   [6:53 remaining] 7,  [6:42 remaining] 8,  [6:32 remaining] 9,
   [6:22 remaining] 10,  [6:13 remaining] 11,  [6:01 remaining] 12,
   [5:51 remaining] 13,  [5:45 remaining] 14,  [5:36 remaining] 15,
   [5:26 remaining] 16,  [5:17 remaining] 17,  [5:07 remaining] 18,
   [4:57 remaining] 19,  [4:48 remaining] 20,  [4:37 remaining] 21,
   [4:27 remaining] 22,  [4:18 remaining] 23,  [4:07 remaining] 24,
   [3:58 remaining] 25,  [3:48 remaining] 26,  [3:39 remaining] 27,
   [3:29 remaining] 28,  [3:21 remaining] 29,  [3:11 remaining] 30,
   [3:02 remaining] 31,  [2:52 remaining] 32,  [2:43 remaining] 33,
   [2:34 remaining] 34,  [2:24 remaining] 35,  [2:15 remaining] 36,
   [2:05 remaining] 37,  [1:55 remaining] 38,  [1:45 remaining] 39,
   [1:36 remaining] 40,  [1:26 remaining] 41,  [1:16 remaining] 42,
   [1:07 remaining] 43,  [58 sec remaining] 44,  [48 sec remaining] 45,
   [38 sec remaining] 46,  [29 sec remaining] 47,  [19 sec remaining] 48,
   [10 sec remaining] 
  49.
  
  Done.

mad

  
    Maximum absolute deviation test of CSR
    Monte Carlo test based on 49 simulations
    Summary function: L(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 7.15641113464284] km
    Test statistic: Maximum absolute deviation
    Deviation = observed minus theoretical
  
  data:  warsaw_ppp_ff
  mad = 2.8776, rank = 1, p-value = 0.02

dclf

  
    Diggle-Cressie-Loosmore-Ford test of CSR
    Monte Carlo test based on 49 simulations
    Summary function: L(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 7.15641113464284] km
    Test statistic: Integral of squared absolute deviation
    Deviation = observed minus theoretical
  
  data:  warsaw_ppp_ff
  u = 40.437, rank = 1, p-value = 0.02

Marked

Porównanie jądrowych estymatorów gęstości punktów fast food i szkół wskazuje, że w znaczenij mierze pokrywają się, oba typy koncentrują się w centrum, choć koncentracja fast foodów jest intensywniejsza.

# Density by type
d_split <- density(split(warsaw_ppp), sigma = bw.scott(unmark(warsaw_ppp)))
par(mfrow = c(1, 2))
plot(d_split$`Fast Food`, main = "Fast food density")
plot(d_split$School, main = "School density")

par(mfrow = c(1, 1))

Analiza względnego ryzyka potwierdza, że punkty fast food są nadreprezentowane względem szkół w centrum.

# Relative risk
par(mar = c(.5, .5, .5, .5))
probs_rr   <- relrisk.ppp(warsaw_ppp, sigma = bw.scott(warsaw_ppp), relative = FALSE)
plot(probs_rr, main = "P(type | location)")

risk_ff    <- relrisk.ppp(warsaw_ppp, sigma = bw.scott(warsaw_ppp),
                          casecontrol = TRUE, relative = TRUE,
                          control = "School")
plot(risk_ff, main = "Relative risk: Fast Food vs School")

Funcja K wskazuje, że punkty fast food ze szkołami, dla dowolnie przyjętego promienia w obszar wpada większa liczba szkół niż by to wynikało z losowego rozmieszczenia.

#Cross K
plot(Kcross(warsaw_ppp, i = "Fast Food", j = "School"),
     main = "Kcross  FF → School")

Odległośc do najbliższego sąsiada innego typu pokazuje, że medianowy punkt fast food leży 0.24 km od najbliższej szkoły, podczas gdy medianowa szkoła leży 0.26 km od najbliższego lokalu fast food. Rozkłady obu odległości są prawoskośne, co wskazuje, że zdecydowana większość typów znajduje się w pobliiżu przeciwnego typu.

# NN distances between types
# cross-distances: FF to nearest school
nn_ff_sc <- nncross(warsaw_ppp_ff, warsaw_ppp_sc)
summary(nn_ff_sc$dist)

      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003705 0.155570 0.239266 0.301609 0.380119 2.000167

hist(nn_ff_sc$dist, breaks = 30, col = "tomato", border = "white",
     main = "FF – dist to nearest school (km)", xlab = "km")

nn_sc_ff <- nncross(warsaw_ppp_sc, warsaw_ppp_ff)
summary(nn_sc_ff$dist)

      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  0.003705 0.145989 0.261456 0.377914 0.424768 4.245132

hist(nn_sc_ff$dist, breaks = 30, col = "steelblue", border = "white",
     main = "School – dist to nearest FF (km)", xlab = "km")

Zarówno test Dixona, jak i test segregacji wskazują na istotny związek przestrzenny między punktami sprzedaży fast foodów a szkołami (p < 0.05), co pozwala odrzucić hipotezę o niezależności przestrzennej tych dwóch typów.

# Marktable & Dixon test
m.table <- marktable(warsaw_ppp, N = 1, collapse = TRUE)
m.table

             neighbour
  point       Fast Food School
    Fast Food      1220    201
    School          263    591

rowSums(m.table)

  Fast Food    School 
       1421       854

# Dixon test for spatial segregation via NN contingency table
dixon.res <- dixon(as.data.frame(warsaw_ppp))

  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
  21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
  41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
  61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
  81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
  99.

dixon.res$tablaC

                               df  Chi-sq P.asymp   P.rand
  Overall segregation            2 443.07       0        0
  From  Fast Food                1 339.65       0        0
  From  School                   1 302.89       0        0

# Segregation test
segregation.test.ppp(warsaw_ppp, nsim = 49)

  Computing observed value... Done.
  Computing 49 simulated values...

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  1,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  2,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  3,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  4,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  5,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  6,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  7,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  8,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  9,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  10,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  11,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  12,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  13,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  14,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  15,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  16,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  17,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  18,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  19,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  20,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  21,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  22,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  23,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  24,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  25,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  26,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  27,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  28,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  29,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  30,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  31,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  32,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  33,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  34,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  35,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  36,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  37,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  38,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  39,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  40,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  41,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  42,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  43,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  44,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  45,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  46,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  47,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  48,

  Warning: Negative Likelihood Cross-Validation criterion was minimised at
  right-hand end of interval [0.0181, 2.33]; use arguments 'hmin', 'hmax' to
  specify a wider interval for bandwidth 'sigma'

  
  49.
  Done.

  
    Monte Carlo test of spatial segregation of types
  
  data:  warsaw_ppp
  T = 135.84, p-value = 0.02

W przeciwieństwie do wyników powyższych testów wyniki analizy funkcji Lcross przy losowym przesunięciu sugeruje, że typy są niezależne. Funkcja leży w granicach koperty na wszystkihc odległościach. Oznacza to, że lokale fast food nie lokują się konkretnie blisko szkół ponad to, co byłoby oczekiwane z niezależnego koncentrowania się typów w centrum. Odzwierciedla to zapewne gętość zaludnienia.

# Envelopes – cross functions
# Random shift test (independence of two types)
# rectangular window needed for rshift
Window(warsaw_ppp) <- owin(range(warsaw_ppp$x) + c(-0.1, 0.1),
                           range(warsaw_ppp$y) + c(-0.1, 0.1))


E_cross_global <- envelope(warsaw_ppp, Lcross,
                           nsim     = 49,
                           i        = "Fast Food", j = "School",
                           simulate = expression(rshift(warsaw_ppp, radius = 1)),
                           global   = TRUE)
  Generating 98 simulations by evaluating expression (49 to estimate the mean and 
  49 to calculate envelopes) ...
  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
  21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
  41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
  61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
  81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 
  98.
  
  Done.
plot(E_cross_global, main = "Lcross global envelope – rshift")

# Restoring proper window
Window(warsaw_ppp) <- rescale(
  as.owin(st_transform(warsaw_poly, crs = 3857)), rsc, "km")

MODELE

Aby formalnie zbadać, czy bliskość szkoły wpływa na intensywność występowania lokali fast foodów, dopasowano szereg modeli procesu punktowego Poissona, w których zmienną onjaśnianą stanowiła liczba lokali fast foodowych. Zdefiniowano dwie przestrzenne współzmienne: odległość od najbliższej szkoły (dist_school_im) oraz estymator jądrowy gęstości szkół (school_dens_im), przedstawione jako pixel image w oknie obejmującym Warszawę. Wszystkie modele zawierają wycentrowany układ współrzędnych w celu zapewnienia stabilności numerycznej polynomials.

# Rescalling
ctr <- centroid.owin(Window(warsaw_ppp))
warsaw_ppp_ff <- shift(warsaw_ppp_ff, c(-ctr$x, -ctr$y))
warsaw_ppp_sc <- shift(warsaw_ppp_sc, c(-ctr$x, -ctr$y))
warsaw_ppp <- shift(warsaw_ppp, c(-ctr$x, -ctr$y))

# Covariates as pixel images
dist_school_im <- as.im(distfun(warsaw_ppp_sc), W = Window(warsaw_ppp_ff))
school_dens_im <- density(warsaw_ppp_sc,
                          sigma = bw.scott(warsaw_ppp_sc),
                          leaveoneout = FALSE, positive = TRUE)

plot(dist_school_im,  main = "Distance to nearest school (km)")

plot(school_dens_im,  main = "School kernel density")

Przed dopasowaniem modelu parametrycznego sprawdzamy nieparametryczne związki między intensywnością lokali fast fooda a współzmiennymi uzywając oszacowania rhohat. W przypadku dystansu do szkół widać związek kwadratowy. Intensywność jest najwyższa ok. 100-200 metrów od szkół i spada zarówno z wyższą odległością, jak i niższą. Wskazuje to na potrzebę użycia kwadratu zmiennej w specyfikacji modelu. W przypadku gęstości szkół widzimy monotoniczny dodatni związek z intensywnością fast foodów.

# rhohat: empirical intensity vs covariate
rho_dist  <- rhohat(warsaw_ppp_ff, dist_school_im)
rho_sdens <- rhohat(warsaw_ppp_ff, school_dens_im)

par(mfrow = c(2, 1))
plot(rho_dist,  main = "rhohat – dist to school")
plot(rho_sdens, main = "rhohat – school density")

par(mfrow = c(1, 1))

Na podstawie AIC porównano 10 specyfickaji modelu, od modelu zerowego do pełnej specyfikacji zawieracjącej polynomial 3. rzędu i obie współzmienne. Test ilorazu wiarogodności potwierdza, że poly(3) jest lepiej dopasowany niż poly(2) (p < .001), co uzasadnia dodatkowe parametry.

# A. Poisson models – fast food as DV
 
# Null
m0 <- ppm(warsaw_ppp_ff ~ 1)

# Linear trend
m1 <- ppm(warsaw_ppp_ff ~ x + y)

# Polynomial trends
m2 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 2))
m3 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3))

# Distance to school
m4 <- ppm(warsaw_ppp_ff ~ dist_school_im)
m5 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im)

# School density
m6 <- ppm(warsaw_ppp_ff ~ school_dens_im)
m7 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3) + school_dens_im)

# Both covariates
m8 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im + school_dens_im)

m9 <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im + I(dist_school_im^2) + school_dens_im)

Uwzględnienie związku parabolicznego okazało się zbędne ponieważ jego współczynnik jest nieistotny, a kryterium informacyjne nie poprawiło się. Na podstawie kryterium informacyjnego AIC najlepszy jest więc model M8: polynom(x, y, 3) + dist_school_im + school_dens_im.

# Model comparison
as.data.frame(
  cbind(
    sapply(list(m0, m1, m2, m3, m4, m5, m6, m7, m8, m9), AIC),
    c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8", "M9"))) %>% 
  rename(AIC = V1, Model = V2) %>% mutate(AIC = as.numeric(AIC)) %>% 
  arrange(AIC)

             AIC Model
  1  -2629.57084    M8
  2  -2627.69749    M9
  3  -2414.29486    M7
  4  -2403.41608    M5
  5  -2207.92053    M6
  6  -2041.39898    M3
  7  -1968.79637    M2
  8  -1481.83678    M4
  9   -319.05875    M1
  10   -28.15119    M0

anova(m2, m3, test = "LR")

  Analysis of Deviance Table
  
  Model 1: ~x + y + I(x^2) + I(x * y) + I(y^2)   Poisson
  Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3)   Poisson
    Npar Df Deviance  Pr(>Chi)    
  1    6                          
  2   10  4   80.603 < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m5, m8, test = "LR")

  Analysis of Deviance Table
  
  Model 1: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3) + dist_school_im      Poisson
  Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3) + dist_school_im + school_dens_im     Poisson
    Npar Df Deviance  Pr(>Chi)    
  1   11                          
  2   12  1   228.16 < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m8, m9, test = "LR")

  Analysis of Deviance Table
  
  Model 1: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3) + dist_school_im + school_dens_im     Poisson
  Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3) + dist_school_im + I(dist_school_im^2) + school_dens_im   Poisson
    Npar Df Deviance Pr(>Chi)
  1   12                     
  2   13  1  0.12666   0.7219

par(mar = c(0,0,0,0))
plot(simulate(m8, nsim = 4), main = "Poisson simulation")

  Generating 4 simulated patterns ...1, 2, 3, 
  4.

par(mar = c(5,4,4,2))

# GOF for best Poisson model
best_pois <- m8   

diagnose.ppm(best_pois)

  Model diagnostics (raw residuals)
  Diagnostics available:
    four-panel plot
    mark plot 
    smoothed residual field
    x cumulative residuals
    y cumulative residuals
    sum of all residuals
  sum of raw residuals in entire window = -2.353e-11
  area of entire window = 517.2
  quadrature area = 516.5
  range of smoothed field =  [-0.5211, 0.2791]

qqplot.ppm(best_pois)

  Extracting model information...Evaluating trend...done.
  Simulating 100 realisations... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
  21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
  41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
  61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
  81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 
  100.
  
  Diagnostic info:
   simulated patterns contained an average of 1403.99 points.
  Calculating quantiles...averaging.....Done.

envelope(best_pois, Linhom, nsim = 49, global = TRUE) |> plot()

  Generating 98 simulated patterns (49 to estimate the mean and 49 to calculate 
  envelopes) ...
  1, 2, 
  [23:04 remaining, estimate finish 2026-05-21 04:27:41]
  3, 
  [21:45 remaining, estimate finish 2026-05-21 04:26:35]
  4, 
  [21:30 remaining, estimate finish 2026-05-21 04:26:35]
  5, 
  [21:19 remaining, estimate finish 2026-05-21 04:26:37]
  6, 
  [20:55 remaining, estimate finish 2026-05-21 04:26:26]
  7, 
  [20:21 remaining, estimate finish 2026-05-21 04:26:05]
  8, 
  [20:05 remaining, estimate finish 2026-05-21 04:26:02]
  9, 
  [19:58 remaining, estimate finish 2026-05-21 04:26:09]
  10, 
  [19:51 remaining, estimate finish 2026-05-21 04:26:16]
  11, 
  [19:56 remaining, estimate finish 2026-05-21 04:26:37]
  12, 
  [19:40 remaining, estimate finish 2026-05-21 04:26:34]
  13, 
  [19:27 remaining, estimate finish 2026-05-21 04:26:35]
  14, 
  [19:17 remaining, estimate finish 2026-05-21 04:26:39]
  15, 
  [19:02 remaining, estimate finish 2026-05-21 04:26:38]
  16, 
  [18:52 remaining, estimate finish 2026-05-21 04:26:42]
  17, 
  [18:30 remaining, estimate finish 2026-05-21 04:26:32]
  18, 
  [18:22 remaining, estimate finish 2026-05-21 04:26:39]
  19, 
  [18:08 remaining, estimate finish 2026-05-21 04:26:38]
  20, 
  [17:53 remaining, estimate finish 2026-05-21 04:26:38]
  21, 
  [17:38 remaining, estimate finish 2026-05-21 04:26:36]
  22, 
  [17:23 remaining, estimate finish 2026-05-21 04:26:34]
  23, 
  [17:09 remaining, estimate finish 2026-05-21 04:26:33]
  24, 
  [16:51 remaining, estimate finish 2026-05-21 04:26:28]
  25, 
  [16:35 remaining, estimate finish 2026-05-21 04:26:25]
  26, 
  [16:20 remaining, estimate finish 2026-05-21 04:26:23]
  27, 
  [16:05 remaining, estimate finish 2026-05-21 04:26:22]
  28, 
  [15:53 remaining, estimate finish 2026-05-21 04:26:23]
  29, 
  [15:42 remaining, estimate finish 2026-05-21 04:26:27]
  30, 
  [15:28 remaining, estimate finish 2026-05-21 04:26:26]
  31, 
  [15:16 remaining, estimate finish 2026-05-21 04:26:29]
  32, 
  [15:01 remaining, estimate finish 2026-05-21 04:26:27]
  33, 
  [14:48 remaining, estimate finish 2026-05-21 04:26:28]
  34, 
  [14:33 remaining, estimate finish 2026-05-21 04:26:26]
  35, 
  [14:21 remaining, estimate finish 2026-05-21 04:26:28]
  36, 
  [14:06 remaining, estimate finish 2026-05-21 04:26:27]
  37, 
  [13:51 remaining, estimate finish 2026-05-21 04:26:25]
  38, 
  [13:40 remaining, estimate finish 2026-05-21 04:26:28]
  39, 
  [13:25 remaining, estimate finish 2026-05-21 04:26:26]
  40, 
  [13:12 remaining, estimate finish 2026-05-21 04:26:28]
  41, 
  [13:00 remaining, estimate finish 2026-05-21 04:26:31]
  42, 
  [12:48 remaining, estimate finish 2026-05-21 04:26:32]
  43, 
  [12:35 remaining, estimate finish 2026-05-21 04:26:34]
  44, 
  [12:21 remaining, estimate finish 2026-05-21 04:26:33]
  45, 
  [12:07 remaining, estimate finish 2026-05-21 04:26:33]
  46, 
  [11:55 remaining, estimate finish 2026-05-21 04:26:36]
  47, 
  [11:41 remaining, estimate finish 2026-05-21 04:26:36]
  48, 
  [11:28 remaining, estimate finish 2026-05-21 04:26:38]
  49, 
  [11:14 remaining, estimate finish 2026-05-21 04:26:38]
  50, 
  [11:01 remaining, estimate finish 2026-05-21 04:26:40]
  51, 
  [10:49 remaining, estimate finish 2026-05-21 04:26:43]
  52, 
  [10:35 remaining, estimate finish 2026-05-21 04:26:41]
  53, 
  [10:20 remaining, estimate finish 2026-05-21 04:26:40]
  54, 
  [10:06 remaining, estimate finish 2026-05-21 04:26:40]
  55,  [9:52 remaining] 56,  [9:38 remaining] 57,  [9:24 remaining] 58,  [9:11 remaining] 59,  [8:57 remaining] 60,  [8:44 remaining] 61,  [8:30 remaining] 62,  [8:16 remaining] 63,  [8:02 remaining] 64,  [7:48 remaining] 65,  [7:35 remaining] 66,  [7:21 remaining] 67,  [7:08 remaining] 68,  [6:55 remaining] 69,  [6:41 remaining] 70,  [6:27 remaining] 71,  [6:13 remaining] 72,  [5:59 remaining] 73,  [5:45 remaining] 74,
   [5:32 remaining] 75,  [5:18 remaining] 76,  [5:04 remaining] 77,  [4:51 remaining] 78,  [4:37 remaining] 79,  [4:23 remaining] 80,  [4:09 remaining] 81,  [3:55 remaining] 82,  [3:41 remaining] 83,  [3:27 remaining] 84,  [3:14 remaining] 85,  [3:00 remaining] 86,  [2:46 remaining] 87,  [2:32 remaining] 88,  [2:18 remaining] 89,  [2:05 remaining] 90,  [1:51 remaining] 91,  [1:37 remaining] 92,  [1:23 remaining] 93,  [1:09 remaining] 94,
   [55 sec remaining] 95,  [41 sec remaining] 96,  [28 sec remaining] 97,  [14 sec remaining] 
  98.
  
  Done.

mad.test(best_pois,  Linhom, nsim = 49)

  Generating 49 simulated patterns  ...
  1, 2, 
  [10:43 remaining, estimate finish 2026-05-21 04:38:10]
  3, 
  [10:44 remaining, estimate finish 2026-05-21 04:38:26]
  4, 
  [10:43 remaining, estimate finish 2026-05-21 04:38:39]
  5, 
  [10:42 remaining, estimate finish 2026-05-21 04:38:54]
  6, 
  [10:14 remaining, estimate finish 2026-05-21 04:38:39]
  7, 
  [10:04 remaining, estimate finish 2026-05-21 04:38:44]
  8,  [9:45 remaining] 9,  [9:32 remaining] 10,  [9:18 remaining] 11,  [9:01 remaining] 12,  [8:48 remaining] 13,  [8:34 remaining] 14,  [8:20 remaining] 15,  [8:02 remaining] 16,  [7:47 remaining] 17,  [7:31 remaining] 18,  [7:14 remaining] 19,  [7:00 remaining] 20,  [6:47 remaining] 21,  [6:29 remaining] 22,  [6:16 remaining] 23,  [6:01 remaining] 24,  [5:46 remaining] 25,  [5:33 remaining] 26,  [5:19 remaining] 27,
   [5:05 remaining] 28,  [4:52 remaining] 29,  [4:38 remaining] 30,  [4:24 remaining] 31,  [4:10 remaining] 32,  [3:56 remaining] 33,  [3:42 remaining] 34,  [3:29 remaining] 35,  [3:15 remaining] 36,  [3:01 remaining] 37,  [2:47 remaining] 38,  [2:33 remaining] 39,  [2:19 remaining] 40,  [2:05 remaining] 41,  [1:52 remaining] 42,  [1:38 remaining] 43,  [1:24 remaining] 44,  [1:10 remaining] 45,  [56 sec remaining] 46,  [42 sec remaining] 47,
   [28 sec remaining] 48,  [14 sec remaining] 
  49.
  
  Done.

  
    Maximum absolute deviation test of fitted Poisson model
    Monte Carlo test based on 49 simulations
    Summary function: L[inhom](r)
    Reference function: sample mean
    Alternative: two.sided
    Interval of distance values: [0, 7.15641113464284] km
    Test statistic: Maximum absolute deviation
    Deviation = leave-one-out
  
  data:  best_pois
  mad = 0.42362, rank = 1, p-value = 0.02

Test dopasowania modelu ujawnia, iż po uwzględnieniu współzmiennych i trendów wciąż pozostaje klastrowanie na niskich odległościahc. Z tego względu wyestymujemy model Geyera, który dopuszcza interakcję między punktami na ustalonej odległości. Promień interakcji dobrano na podstawie maksymalizacji pseudowiarogodności i wyniósł on 0.161km.

# Gibbs models – fast food
# Profile likelihood – Geyer r
prof_geyer <- profilepl(
  data.frame(r = seq(0.001, 1.0, by = 0.01), sat = 2),
  Geyer,
  warsaw_ppp_ff ~ polynom(x,y,3) 
)

  (computing rbord)

  comparing 100 models...

  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
  21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
  41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
  61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
  81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 
  100.

  fitting optimal model...

  done.

plot(prof_geyer, main = "Profile PL – Geyer r")

r_gey <- prof_geyer$param[which.max(prof_geyer$prof), 1]
r_gey

  [1] 0.161

m_geyer <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3), Geyer(r = r_gey, sat = 2))
# gamma > 1 → clustering tendency

m_geyer_dist <- ppm(warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im + school_dens_im, Geyer(r = r_gey, sat = 2))

# Compare all
sapply(list(m_geyer, m_geyer_dist), AIC)

  [1] -3701.000 -3938.586

anova(m_geyer, m_geyer_dist, test = "LR")

  Warning: anova.ppm now computes the *adjusted* deviances when the models are
  not Poisson processes.

  Analysis of Deviance Table
  
  Model 1: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3)   Geyer
  Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * y^2) + I(y^3) + dist_school_im + school_dens_im     Geyer
    Npar Df AdjDeviance  Pr(>Chi)    
  1   11                             
  2   13  2      200.98 < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Simulations from Gibbs models (visual GOF)
par(mar = c(0,0,0,0))
plot(simulate(m_geyer_dist, nsim = 4), main = "Geyer simulation")

  Generating 4 simulated patterns ...1, 2, 3, 
  4.

par(mar = c(5,4,4,2))

# Envelopes and tests against Gibbs
envelope(m_geyer_dist, Linhom, nsim = 49, global = TRUE) |> plot()

  Generating 98 simulated patterns (49 to estimate the mean and 49 to calculate 
  envelopes) ...
  1, 2, 
  [23:45 remaining, estimate finish 2026-05-21 10:51:59]
  3, 
  [25:28 remaining, estimate finish 2026-05-21 10:54:00]
  4, 
  [26:16 remaining, estimate finish 2026-05-21 10:55:05]
  5, 
  [25:38 remaining, estimate finish 2026-05-21 10:54:44]
  6, 
  [25:58 remaining, estimate finish 2026-05-21 10:55:22]
  7, 
  [25:49 remaining, estimate finish 2026-05-21 10:55:31]
  8, 
  [25:25 remaining, estimate finish 2026-05-21 10:55:23]
  9, 
  [24:51 remaining, estimate finish 2026-05-21 10:55:04]
  10, 
  [24:26 remaining, estimate finish 2026-05-21 10:54:55]
  11, 
  [24:15 remaining, estimate finish 2026-05-21 10:55:01]
  12, 
  [23:31 remaining, estimate finish 2026-05-21 10:54:31]
  13, 
  [23:07 remaining, estimate finish 2026-05-21 10:54:22]
  14, 
  [22:48 remaining, estimate finish 2026-05-21 10:54:19]
  15, 
  [22:21 remaining, estimate finish 2026-05-21 10:54:06]
  16, 
  [22:06 remaining, estimate finish 2026-05-21 10:54:08]
  17, 
  [21:45 remaining, estimate finish 2026-05-21 10:54:02]
  18, 
  [21:24 remaining, estimate finish 2026-05-21 10:53:56]
  19, 
  [21:03 remaining, estimate finish 2026-05-21 10:53:49]
  20, 
  [20:48 remaining, estimate finish 2026-05-21 10:53:51]
  21, 
  [20:30 remaining, estimate finish 2026-05-21 10:53:49]
  22, 
  [20:13 remaining, estimate finish 2026-05-21 10:53:47]
  23, 
  [20:00 remaining, estimate finish 2026-05-21 10:53:51]
  24, 
  [19:43 remaining, estimate finish 2026-05-21 10:53:50]
  25, 
  [19:23 remaining, estimate finish 2026-05-21 10:53:44]
  26, 
  [19:07 remaining, estimate finish 2026-05-21 10:53:44]
  27, 
  [18:50 remaining, estimate finish 2026-05-21 10:53:42]
  28, 
  [18:32 remaining, estimate finish 2026-05-21 10:53:40]
  29, 
  [18:20 remaining, estimate finish 2026-05-21 10:53:45]
  30, 
  [18:05 remaining, estimate finish 2026-05-21 10:53:47]
  31, 
  [17:49 remaining, estimate finish 2026-05-21 10:53:47]
  32, 
  [17:33 remaining, estimate finish 2026-05-21 10:53:47]
  33, 
  [17:15 remaining, estimate finish 2026-05-21 10:53:44]
  34, 
  [17:02 remaining, estimate finish 2026-05-21 10:53:48]
  35, 
  [16:48 remaining, estimate finish 2026-05-21 10:53:51]
  36, 
  [16:34 remaining, estimate finish 2026-05-21 10:53:54]
  37, 
  [16:18 remaining, estimate finish 2026-05-21 10:53:54]
  38, 
  [16:01 remaining, estimate finish 2026-05-21 10:53:53]
  39, 
  [15:47 remaining, estimate finish 2026-05-21 10:53:56]
  40, 
  [15:30 remaining, estimate finish 2026-05-21 10:53:54]
  41, 
  [15:16 remaining, estimate finish 2026-05-21 10:53:57]
  42, 
  [14:59 remaining, estimate finish 2026-05-21 10:53:57]
  43, 
  [14:43 remaining, estimate finish 2026-05-21 10:53:57]
  44, 
  [14:27 remaining, estimate finish 2026-05-21 10:53:56]
  45, 
  [14:12 remaining, estimate finish 2026-05-21 10:53:59]
  46, 
  [13:57 remaining, estimate finish 2026-05-21 10:54:00]
  47, 
  [13:42 remaining, estimate finish 2026-05-21 10:54:02]
  48, 
  [13:25 remaining, estimate finish 2026-05-21 10:54:01]
  49, 
  [13:09 remaining, estimate finish 2026-05-21 10:54:02]
  50, 
  [12:53 remaining, estimate finish 2026-05-21 10:54:02]
  51, 
  [12:36 remaining, estimate finish 2026-05-21 10:54:00]
  52, 
  [12:19 remaining, estimate finish 2026-05-21 10:53:58]
  53, 
  [12:04 remaining, estimate finish 2026-05-21 10:54:00]
  54, 
  [11:48 remaining, estimate finish 2026-05-21 10:54:01]
  55, 
  [11:33 remaining, estimate finish 2026-05-21 10:54:02]
  56, 
  [11:18 remaining, estimate finish 2026-05-21 10:54:05]
  57, 
  [11:03 remaining, estimate finish 2026-05-21 10:54:07]
  58, 
  [10:47 remaining, estimate finish 2026-05-21 10:54:07]
  59, 
  [10:30 remaining, estimate finish 2026-05-21 10:54:06]
  60, 
  [10:13 remaining, estimate finish 2026-05-21 10:54:05]
  61,  [9:56 remaining] 62,  [9:40 remaining] 63,  [9:24 remaining] 64,  [9:08 remaining] 65,  [8:52 remaining] 66,  [8:36 remaining] 67,  [8:20 remaining] 68,  [8:03 remaining] 69,  [7:47 remaining] 70,  [7:30 remaining] 71,  [7:15 remaining] 72,  [6:58 remaining] 73,  [6:42 remaining] 74,  [6:26 remaining] 75,  [6:09 remaining] 76,  [5:53 remaining] 77,  [5:37 remaining] 78,  [5:22 remaining] 79,  [5:06 remaining] 80,
   [4:50 remaining] 81,  [4:34 remaining] 82,  [4:18 remaining] 83,  [4:01 remaining] 84,  [3:45 remaining] 85,  [3:29 remaining] 86,  [3:13 remaining] 87,  [2:57 remaining] 88,  [2:41 remaining] 89,  [2:25 remaining] 90,  [2:08 remaining] 91,  [1:52 remaining] 92,  [1:36 remaining] 93,  [1:20 remaining] 94,  [1:04 remaining] 95,  [48 sec remaining] 96,  [32 sec remaining] 97,  [16 sec remaining] 
  98.
  
  Done.

mad.test(m_geyer_dist,  Linhom, nsim = 49)

  Generating 49 simulated patterns  ...
  1, 2, 
  [11:54 remaining, estimate finish 2026-05-21 11:06:37]
  3, 
  [11:34 remaining, estimate finish 2026-05-21 11:06:32]
  4, 
  [11:31 remaining, estimate finish 2026-05-21 11:06:46]
  5, 
  [11:40 remaining, estimate finish 2026-05-21 11:07:13]
  6, 
  [11:32 remaining, estimate finish 2026-05-21 11:07:21]
  7, 
  [11:25 remaining, estimate finish 2026-05-21 11:07:32]
  8, 
  [11:06 remaining, estimate finish 2026-05-21 11:07:28]
  9, 
  [10:46 remaining, estimate finish 2026-05-21 11:07:24]
  10, 
  [10:32 remaining, estimate finish 2026-05-21 11:07:27]
  11, 
  [10:15 remaining, estimate finish 2026-05-21 11:07:25]
  12,  [9:56 remaining] 13,  [9:40 remaining] 14,  [9:29 remaining] 15,  [9:16 remaining] 16,  [8:56 remaining] 17,  [8:38 remaining] 18,  [8:19 remaining] 19,  [8:02 remaining] 20,  [7:46 remaining] 21,  [7:31 remaining] 22,  [7:15 remaining] 23,  [7:02 remaining] 24,  [6:46 remaining] 25,  [6:30 remaining] 26,  [6:13 remaining] 27,  [5:55 remaining] 28,  [5:39 remaining] 29,  [5:23 remaining] 30,  [5:08 remaining] 31,
   [4:51 remaining] 32,  [4:35 remaining] 33,  [4:19 remaining] 34,  [4:02 remaining] 35,  [3:46 remaining] 36,  [3:30 remaining] 37,  [3:13 remaining] 38,  [2:57 remaining] 39,  [2:41 remaining] 40,  [2:25 remaining] 41,  [2:09 remaining] 42,  [1:53 remaining] 43,  [1:37 remaining] 44,  [1:21 remaining] 45,  [1:04 remaining] 46,  [48 sec remaining] 47,  [32 sec remaining] 48,  [16 sec remaining] 
  49.
  
  Done.

  
    Maximum absolute deviation test of fitted Gibbs model
    Monte Carlo test based on 49 simulations
    Summary function: L[inhom](r)
    Reference function: sample mean
    Alternative: two.sided
    Interval of distance values: [0, 7.15641113464284] km
    Test statistic: Maximum absolute deviation
    Deviation = leave-one-out
  
  data:  m_geyer_dist
  mad = 0.39581, rank = 1, p-value = 0.02

Mimo uwzględnienia interakcji między punktami na którkich odległościach pozostaje nadmierne względem CSR klastrowanie.

# Full summary table
all_models <- list(
  "Null" = m0,
  "Linear" = m1,
  "Poly2" = m2,
  "Poly3" = m3,
  "dist_school" = m4,
  "Poly3+dist" = m5,
  "school_dens" = m6,
  "Poly3+sdens" = m7,
  "Poly3+dist+sdens" = m8,
  "Poly3+dist+dist2+sdens" = m9,
  "Geyer" = m_geyer,
  "Geyer_full" = m_geyer_dist
)
aic_tab <- data.frame(
  model = names(all_models),
  AIC   = sapply(all_models, AIC),
  npar  = sapply(all_models, function(m) length(coef(m)))
)
print(aic_tab[order(aic_tab$AIC), ])

                                          model         AIC npar
  Geyer_full                         Geyer_full -3938.58578   13
  Geyer                                   Geyer -3700.99956   11
  Poly3+dist+sdens             Poly3+dist+sdens -2629.57084   12
  Poly3+dist+dist2+sdens Poly3+dist+dist2+sdens -2627.69749   13
  Poly3+sdens                       Poly3+sdens -2414.29486   11
  Poly3+dist                         Poly3+dist -2403.41608   11
  school_dens                       school_dens -2207.92053    2
  Poly3                                   Poly3 -2041.39898   10
  Poly2                                   Poly2 -1968.79637    6
  dist_school                       dist_school -1481.83678    2
  Linear                                 Linear  -319.05875    3
  Null                                     Null   -28.15119    1

summary(m8)

  Point process model
  Fitted to data: warsaw_ppp_ff
  Fitting method: maximum likelihood (Berman-Turner approximation)
  Model was fitted using glm()
  Algorithm converged
  Call:
  ppm.formula(Q = warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im + 
      school_dens_im)
  Edge correction: "border"
    [border correction distance r = 0 ]
  --------------------------------------------------------------------------------
  Quadrature scheme (Berman-Turner) = data + dummy + weights
  
  Data pattern:
  Planar point pattern:  1421 points
  Average intensity 2.75 points per square km
  Window: polygonal boundary
  single connected closed polygon with 1745 vertices
  enclosing rectangle: [-13.375867, 15.249778] x [-14.642966, 15.476199] km
                       (28.63 x 30.12 km)
  Window area = 517.24 square km
  Unit of length: 1 km
  Fraction of frame area: 0.6
  
  Dummy quadrature points:
       80 x 80 grid of dummy points, plus 4 corner points
       dummy spacing: 0.3578206 x 0.3764896 km
  
  Original dummy parameters: =
  Planar point pattern:  3953 points
  Average intensity 7.64 points per square km
  Window: polygonal boundary
  single connected closed polygon with 1745 vertices
  enclosing rectangle: [-13.375867, 15.249778] x [-14.642966, 15.476199] km
                       (28.63 x 30.12 km)
  Window area = 517.24 square km
  Unit of length: 1 km
  Fraction of frame area: 0.6
  Quadrature weights:
       (counting weights based on 80 x 80 array of rectangular tiles)
  All weights:
    range: [0.00481, 0.135] total: 516
  Weights on data points:
    range: [0.00481, 0.0674]    total: 45.2
  Weights on dummy points:
    range: [0.00481, 0.135] total: 471
  --------------------------------------------------------------------------------
  FITTED :
  
  Nonstationary Poisson process
  
  ---- Intensity: ----
  
  Log intensity: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x 
  * y^2) + I(y^3) + dist_school_im + school_dens_im
  Model depends on external covariates 'dist_school_im' and 'school_dens_im'
  Covariates provided:
    dist_school_im: im
    school_dens_im: im
  
  Fitted trend coefficients:
     (Intercept)              x              y         I(x^2)       I(x * y) 
    0.7120427274  -0.0602434370  -0.0894984662  -0.0054653562   0.0043080633 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
   -0.0026374866   0.0001862507   0.0009568807  -0.0006473063   0.0002505739 
  dist_school_im school_dens_im 
   -1.6875117961   0.4159937846 
  
                      Estimate         S.E.       CI95.lo       CI95.hi Ztest
  (Intercept)     0.7120427274 0.1501087700  4.178349e-01  1.006251e+00   ***
  x              -0.0602434370 0.0169760012 -9.351579e-02 -2.697109e-02   ***
  y              -0.0894984662 0.0129208556 -1.148229e-01 -6.417405e-02   ***
  I(x^2)         -0.0054653562 0.0015604532 -8.523788e-03 -2.406924e-03   ***
  I(x * y)        0.0043080633 0.0024554504 -5.045311e-04  9.120658e-03      
  I(y^2)         -0.0026374866 0.0014974807 -5.572495e-03  2.975216e-04      
  I(x^3)          0.0001862507 0.0001736040 -1.540068e-04  5.265083e-04      
  I(x^2 * y)      0.0009568807 0.0003480210  2.747721e-04  1.638989e-03    **
  I(x * y^2)     -0.0006473063 0.0003399462 -1.313589e-03  1.897597e-05      
  I(y^3)          0.0002505739 0.0001236771  8.171161e-06  4.929766e-04     *
  dist_school_im -1.6875117961 0.1270184972 -1.936463e+00 -1.438560e+00   ***
  school_dens_im  0.4159937846 0.0269657469  3.631419e-01  4.688457e-01   ***
                       Zval
  (Intercept)      4.743512
  x               -3.548741
  y               -6.926667
  I(x^2)          -3.502416
  I(x * y)         1.754490
  I(y^2)          -1.761283
  I(x^3)           1.072848
  I(x^2 * y)       2.749491
  I(x * y^2)      -1.904143
  I(y^3)           2.026032
  dist_school_im -13.285559
  school_dens_im  15.426748
  
  ----------- gory details -----
  
  Fitted regular parameters (theta):
     (Intercept)              x              y         I(x^2)       I(x * y) 
    0.7120427274  -0.0602434370  -0.0894984662  -0.0054653562   0.0043080633 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
   -0.0026374866   0.0001862507   0.0009568807  -0.0006473063   0.0002505739 
  dist_school_im school_dens_im 
   -1.6875117961   0.4159937846 
  
  Fitted exp(theta):
     (Intercept)              x              y         I(x^2)       I(x * y) 
       2.0381504      0.9415353      0.9143897      0.9945496      1.0043174 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
       0.9973660      1.0001863      1.0009573      0.9993529      1.0002506 
  dist_school_im school_dens_im 
       0.1849792      1.5158764

summary(m_geyer_dist)

  Point process model
  Fitted to data: warsaw_ppp_ff
  Fitting method: maximum pseudolikelihood (Berman-Turner approximation)
  Model was fitted using glm()
  Algorithm converged
  Call:
  ppm.formula(Q = warsaw_ppp_ff ~ polynom(x, y, 3) + dist_school_im + 
      school_dens_im, interaction = Geyer(r = r_gey, sat = 2))
  Edge correction: "border"
    [border correction distance r = 0.322 ]
  --------------------------------------------------------------------------------
  Quadrature scheme (Berman-Turner) = data + dummy + weights
  
  Data pattern:
  Planar point pattern:  1421 points
  Average intensity 2.75 points per square km
  Window: polygonal boundary
  single connected closed polygon with 1745 vertices
  enclosing rectangle: [-13.375867, 15.249778] x [-14.642966, 15.476199] km
                       (28.63 x 30.12 km)
  Window area = 517.24 square km
  Unit of length: 1 km
  Fraction of frame area: 0.6
  
  Dummy quadrature points:
       80 x 80 grid of dummy points, plus 4 corner points
       dummy spacing: 0.3578206 x 0.3764896 km
  
  Original dummy parameters: =
  Planar point pattern:  3953 points
  Average intensity 7.64 points per square km
  Window: polygonal boundary
  single connected closed polygon with 1745 vertices
  enclosing rectangle: [-13.375867, 15.249778] x [-14.642966, 15.476199] km
                       (28.63 x 30.12 km)
  Window area = 517.24 square km
  Unit of length: 1 km
  Fraction of frame area: 0.6
  Quadrature weights:
       (counting weights based on 80 x 80 array of rectangular tiles)
  All weights:
    range: [0.00481, 0.135] total: 516
  Weights on data points:
    range: [0.00481, 0.0674]    total: 45.2
  Weights on dummy points:
    range: [0.00481, 0.135] total: 471
  --------------------------------------------------------------------------------
  FITTED :
  
  Nonstationary Geyer saturation process
  
  ---- Trend: ----
  
  Log trend: ~x + y + I(x^2) + I(x * y) + I(y^2) + I(x^3) + I(x^2 * y) + I(x * 
  y^2) + I(y^3) + dist_school_im + school_dens_im
  Model depends on external covariates 'dist_school_im' and 'school_dens_im'
  Covariates provided:
    dist_school_im: im
    school_dens_im: im
  
  Fitted trend coefficients:
     (Intercept)              x              y         I(x^2)       I(x * y) 
    3.042414e-01  -4.290323e-02  -5.978295e-02  -3.944874e-03   1.993108e-03 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
   -2.944440e-03  -2.807688e-05   6.841272e-04  -6.208957e-04   5.712484e-05 
  dist_school_im school_dens_im 
   -1.116069e+00   2.751963e-01 
  
                      Estimate         S.E.       CI95.lo       CI95.hi Ztest
  (Intercept)     3.042414e-01 0.1805325003 -4.959584e-02  0.6580785567      
  x              -4.290323e-02 0.0169120147 -7.605017e-02 -0.0097562866     *
  y              -5.978295e-02 0.0162791966 -9.168959e-02 -0.0278763119   ***
  I(x^2)         -3.944874e-03 0.0018327345 -7.536967e-03 -0.0003527802     *
  I(x * y)        1.993108e-03 0.0025749043 -3.053611e-03  0.0070398281      
  I(y^2)         -2.944440e-03 0.0017148389 -6.305463e-03  0.0004165824      
  I(x^3)         -2.807688e-05 0.0001945530 -4.093937e-04  0.0003532399      
  I(x^2 * y)      6.841272e-04 0.0003498385 -1.543612e-06  0.0013697979      
  I(x * y^2)     -6.208957e-04 0.0004124937 -1.429369e-03  0.0001875771      
  I(y^3)          5.712484e-05 0.0002130338 -3.604138e-04  0.0004746635      
  dist_school_im -1.116069e+00 0.1503999672 -1.410848e+00 -0.8212907406   ***
  school_dens_im  2.751963e-01 0.0290975024  2.181662e-01  0.3322263516   ***
  Interaction     6.549735e-01 0.0422615707  5.721423e-01  0.7378046313   ***
                       Zval
  (Intercept)     1.6852442
  x              -2.5368489
  y              -3.6723527
  I(x^2)         -2.1524524
  I(x * y)        0.7740514
  I(y^2)         -1.7170360
  I(x^3)         -0.1443148
  I(x^2 * y)      1.9555516
  I(x * y^2)     -1.5052247
  I(y^3)          0.2681492
  dist_school_im -7.4206749
  school_dens_im  9.4577291
  Interaction    15.4980864
  
   ---- Interaction: -----
  
  Interaction: Geyer saturation process
  Interaction distance: 0.161
  Saturation parameter: 2
  Fitted interaction parameter gamma:   1.9250915
  
  Relevant coefficients:
  Interaction 
    0.6549735 
  
  ----------- gory details -----
  
  Fitted regular parameters (theta):
     (Intercept)              x              y         I(x^2)       I(x * y) 
    3.042414e-01  -4.290323e-02  -5.978295e-02  -3.944874e-03   1.993108e-03 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
   -2.944440e-03  -2.807688e-05   6.841272e-04  -6.208957e-04   5.712484e-05 
  dist_school_im school_dens_im    Interaction 
   -1.116069e+00   2.751963e-01   6.549735e-01 
  
  Fitted exp(theta):
     (Intercept)              x              y         I(x^2)       I(x * y) 
       1.3555962      0.9580041      0.9419690      0.9960629      1.0019951 
          I(y^2)         I(x^3)     I(x^2 * y)     I(x * y^2)         I(y^3) 
       0.9970599      0.9999719      1.0006844      0.9993793      1.0000571 
  dist_school_im school_dens_im    Interaction 
       0.3275648      1.3167891      1.9250915

Mimo niedoskonałości dopasowania modeli współczynniki wzspózmiennych związanych ze szkołami były istotne w obu modelach, co potwierdza, że fast foody lokują się w pobliżu szkół.

PODSUMOWANIE

Wyniki analizy są zbieżne z wcześniejszymi pracami nt. związków przestrzennych lokacji fast foodów i szkół. W Warszawie występuje klastrowanie między fast foodami oraz istotny związaek przestrzenny między fast foodami i szkołami. Każde 100m odległości od odszło wiąże się ze spadkiem intensywności fast foodów o 10.6%. Gęstość szkół ma komplementrany pozywny efekt. Nie można jednak swierdzić, że lokalizacja fast foodów zależy od lokalizacji szkół, raczej obie wynikają z podobnego procesu generującego dane, np. z gęstości zaludnienia, która powinna determinować zarówno lokalizacje szkół jak i fast foodów.

BIBLIOGRAFIA

1. Simon, P. A., Kwan, D., Angelescu, A., Shih, M., & Fielding, J. E. (2008). Proximity of fast food restaurants to schools: Do neighborhood income and type of school matter? Preventive Medicine, 47(3), 284–288,.
2. Kwate, N. O. A., & Loh, J. M. (2010). Separate and unequal: The influence of neighborhood and school characteristics on spatial proximity between fast food and schools. Preventive Medicine, 51, 153–156.