Wstęp

W niniejszej pracy poddano analizie wzorzec punktowy w postaci lokalizacji dyskontów Lidl i Biedronka na terenie miasta stołecznego Warszawy.

Przegląd literatury

Teoretyczną podstawą rozważań może być model lokalizacji liniowej (Hotelling, 1929) z mikroekonomii. Zgodnie z tym modelem, sklepy różnych sieci będą lokalizować się w podobnej lokalizacji, bo takie działanie zapewni im maksymalizację zysku. Gdy sprzedawcy znajdują się blisko siebie, możemy mówić o równowadze Nasha z teorii gier.

Temat jest również rozpoznany w literaturze dotyczącej analizy wzorców punktowych i liniowych. W pracy pt. “Grocery Retail Location Patterns in Brno: Clustering, Inequality and Street Network Centrality” (Kunc et al, 2025) autorzy przeanalizowali wzorce lokowania się sklepów w czeskim mieście Brno. Badanie wykazało, że lokalizacje sklepów mają tendencje do klastrowania, jednak znaleziono różnicę między lokowaniem się sieci z kapitałem krajowym (centra miast) i zagranicznym (obrzeża). Podobna praca z Korei - “Exploring spatial distribution of retail stores in South Korea” (Lee, 2024) - również wskazała na klastrowanie sklepów. W jej przypadku wykorzystano podobne metody do zastosowanych przeze mnie (funkcja L pochodna do K Ripleya), natomiast w czeskiej pracy zdecydowano sie na estymację rozkładu z wykorzystaniem funckji jądra.

Mimo, że obserwujemy wzrost zainteresowania niniejszym tematem w ostatnim czasie, to podobne prace pojawiają się od dłuższego czasu. Często cytowany jest artykuł “Spatial point pattern analysis of urban retail stores: the case of twelve large- and medium-sized Greek cities” (Tsoutsos et Photis, 2026), którego wnioski są zbieżne z powyższymi. Z kolei jedna z pierwszych tego typu prac “Point pattern analysis: an application to the loyalty networks of chain-stores” pojawiła się już prawie 25 lat temu (Giovani et Lorenzo, 2002). Można znaleźć również próby zastosowania tych metod na polskim rynku, jednak dla innych rodzajów punktów.

Ten projekt będzie przeniesienim powyższych doświadczeń na polski rynek dyskontów spożywczych.

Analiza eksploracyjna danych

Lokalizacje sklepów pobrałem z serwisu OpenStreetMaps z wykorzystaniem pakietu osmdata. Należy zwrócić uwagę na wybranie działajacego serwera API overpass, inaczej kwerendy nie będą się wykonywać.

get_overpass_url()
## [1] "https://overpass.kumi.systems/api/interpreter"
new_url <- "https://overpass-api.de/api/interpreter"
set_overpass_url(new_url)
get_overpass_url()
## [1] "https://overpass-api.de/api/interpreter"

Z kolei kontur Warszawy wykorzystywany do wizualizacji (shapefile) to przetworzone dane z Geoportalu Krajowego geoportal.gov.pl.

warsaw_bb <- getbb("Warsaw")

poviat <- st_read("poviat/A02_Granice_powiatow.shp")
## Reading layer `A02_Granice_powiatow' from data source 
##   `E:\RStudio\awpl testy\poviat\A02_Granice_powiatow.shp' using driver `ESRI Shapefile'
## Simple feature collection with 380 features and 34 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 14.12288 ymin: 49.00204 xmax: 24.14578 ymax: 54.83642
## Geodetic CRS:  ETRS89
poviat_wgs84 <- st_transform(poviat, crs = 4326) # WGS84 spherical projections (as Google)
poviat_mercator <- st_transform(poviat_wgs84, crs = 3857) # mercator

warsaw_map <- poviat_mercator[poviat_mercator$JPT_NAZWA_=="powiat Warszawa",]

Na potrzeby dalszej analizy wykonujemy konwersję CRS na 3857, czyli na mapę typu Mercator używaną między innymi przez Google Maps. Analiza wzroców punktowych opiera się na rzeczywistych odległościach w metrach (a nie koordynatach), dlatego ten sposób projekcji będzie najlepszy.

sklepy <- getbb("Warszawa") %>%
  opq() %>%
  add_osm_feature(key = "shop", value = "supermarket") %>%
  osmdata_sf()

punkty <- st_centroid(sklepy$osm_polygons)
punkty <- st_transform(punkty, crs = 3857)
punkty <- st_intersection(punkty, warsaw_map)
lidle <- punkty[punkty$name == "Lidl",]
biedry <- punkty[punkty$name == "Biedronka",]

Drobna uwaga techniczna - pobieranie danych może chwile zająć lub momentami nie działać poprawnie, ze względu na przeciążenie serwera. Należy wtedy zmienić serwer lub spróbować ponownie później.

ggplot() +
  geom_sf(data = warsaw_map) + 
  geom_sf(data = lidle, aes(color = "Lidl")) +
  geom_sf(data = biedry, aes(color = "Biedronka")) +
  scale_color_manual(name = "Supermarkety", values = c("Lidl" = "blue", "Biedronka" = "red"))

Analiza ilościowa

Plan działania

Zostanie przeprowadzona analiza bez i z podziałem na etykietki, w wyniku której potwierdzimy charakter obserwowanego wzorca punktowego.

Wzorzec bez podziału na etykietki

Analiza pierwsego rzędu

summary(sklepy.ppp.um)
## Planar point pattern:  160 points
## Average intensity 0.3093342 points per square km
## 
## Coordinates are given to 12 decimal places
## 
## Window: polygonal boundary
## single connected closed polygon with 4259 vertices
## enclosing rectangle: [1422.9789, 1451.6043] x [4179.58, 4209.699] km
##                      (28.63 x 30.12 km)
## Window area = 517.24 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6
intensity(sklepy.ppp.um)
## [1] 0.3093342
npoints(sklepy.ppp.um)/area(sklepy.ppp.um)
## [1] 0.3093342
## quadrats
qdr <- quadratcount(sklepy.ppp.um) # by default 5*5
plot(qdr)

qdr <- quadratcount(sklepy.ppp.um, nx=3, ny=3)
intensity(qdr)
##                      x
## y                     [1.42e+03,1.43e+03] (1.43e+03,1.44e+03]
##   (4.2e+03,4.21e+03]            0.3229686           0.2982102
##   (4.19e+03,4.2e+03]            0.5549433           0.2415089
##   [4.18e+03,4.19e+03]           0.1699644           0.3044925
##                      x
## y                     (1.44e+03,1.45e+03]
##   (4.2e+03,4.21e+03]            0.0000000
##   (4.19e+03,4.2e+03]            0.2673112
##   [4.18e+03,4.19e+03]           0.1488145
intensity(qdr, image=T, dimyx = 200) # pixel image, dimyx - to set resolution
## real-valued pixel image
## 200 x 200 pixel array (ny, nx)
## enclosing rectangle: [1423, 1451.6] x [4179.6, 4209.7] km
plot(qdr)

qdr.test <- quadrat.test(sklepy.ppp.um, nx=3, ny=3)
qdr.test # p-value < 0.05, pattern is inhomogenous
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  sklepy.ppp.um
## X2 = 22.973, df = 8, p-value = 0.006798
## alternative hypothesis: two.sided
## 
## Quadrats: 9 tiles (irregular windows)
qdr.test <- quadrat.test(sklepy.ppp.um, nx=3, ny=3,
                         alternative = "clustered")
qdr.test
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  sklepy.ppp.um
## X2 = 22.973, df = 8, p-value = 0.003399
## alternative hypothesis: clustered
## 
## Quadrats: 9 tiles (irregular windows)
qdr.test <- quadrat.test(sklepy.ppp.um, nx=3, ny=3,
                         alternative = "regular")
qdr.test
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  sklepy.ppp.um
## X2 = 22.973, df = 8, p-value = 0.9966
## alternative hypothesis: regular
## 
## Quadrats: 9 tiles (irregular windows)
plot(sklepy.ppp.um)
plot(qdr.test, add=T, col="blue")

## kernel smoothing
d <- density(sklepy.ppp.um)
summary(d)
## real-valued pixel image
## 128 x 128 pixel array (ny, nx)
## enclosing rectangle: [1422.979, 1451.604] x [4179.58, 4209.699] km
## dimensions of each pixel: 0.224 x 0.2353044 km
## Image is defined on a subset of the rectangular grid
## Subset area = 516.173497102695 square km
## Subset area fraction = 0.599
## Pixel values (inside window):
##  range = [0.1491663, 0.5698746]
##  integral = 162.8499
##  mean = 0.3154945
par(mfrow=c(2,2))
plot(d, main="Density") # graph 1 - density
contour(d, main="Contour plot") # graph 2 - contour plot
plot(d, main="Density + \n contour plot") 
contour(d, add=T) # graph 3 - combination of 1 and 2
persp(d, theta = 30, phi = 30, col = "lightblue", main="Perspective plot") # graph 4 - perspective plot

par(mfrow=c(1,1))

cdf.test(sklepy.ppp.um, "x") # covariate given in ""
## 
##  Spatial Kolmogorov-Smirnov test of CSR in two dimensions
## 
## data:  covariate 'x' evaluated at points of 'sklepy.ppp.um' 
##      and transformed to uniform distribution under CSR
## D = 0.12539, p-value = 0.01306
## alternative hypothesis: two-sided
berman.test(sklepy.ppp.um, "x") # Z1 default
## 
##  Berman Z1 test of CSR in two dimensions
## 
## data:  covariate 'x' evaluated at points of 'sklepy.ppp.um'
## Z1 = 0.013979, p-value = 0.9888
## alternative hypothesis: two-sided
berman.test(sklepy.ppp.um, "x", "Z2")
## 
##  Berman Z2 test of CSR in two dimensions
## 
## data:  covariate 'x' evaluated at points of 'sklepy.ppp.um' 
##   and transformed to uniform distribution under CSR
## Z2 = -2.6004, p-value = 0.009312
## alternative hypothesis: two-sided

Wszystkie testy wskazują na występowanie niehomogenicznego wzorca klastrowanego.

Analiza drugiego rzędu

plot(sklepy.ppp.um)

plot(sklepy.ppp.um)
plot(Kinhom(sklepy.ppp.um))

plot(pcfinhom(sklepy.ppp.um))

plot(Linhom(sklepy.ppp.um))

plot(Finhom(sklepy.ppp.um))

plot(Ginhom(sklepy.ppp.um))

plot(Jinhom(sklepy.ppp.um))

## inhom

D <- density(sklepy.ppp.um)
E.inhom <- envelope(sklepy.ppp.um, 
                    Kinhom, 
                    simulate=expression(rpoispp(D)))
## Generating 99 simulations by evaluating expression  ...
## 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.
## 
## Done.
par(mfrow=c(1,2))
plot(D)
plot(E.inhom)

par(mfrow=c(1,1))


E.ptwise <- envelope(sklepy.ppp.um, 
                     Linhom, 
                     simulate=expression(rpoispp(D)),
                     nsim = 99, 
                     fix.n = T)
## Generating 99 simulations by evaluating expression  ...
## 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.
## 
## Done.
plot(E.ptwise)

# Global
E.g <- envelope(sklepy.ppp.um, 
                Linhom, 
                nsim=39, # p-v = 1/(1+nsim)
                simulate=expression(rpoispp(D)),
                rank=1, # min and max to use
                global=TRUE)
## Generating 78 simulations by evaluating expression (39 to estimate the mean and 
## 39 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.
## 
## Done.
plot(E.g)

mad.test(sklepy.ppp.um, Linhom, nsim=19, verbose=T, use.theo=T) # around 1 min
## Generating 19 simulated realisations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.
## 
##  Maximum absolute deviation test of CSR
##  Monte Carlo test based on 19 simulations
##  Summary function: L[inhom](r)
##  Reference function: theoretical
##  Alternative: two.sided
##  Interval of distance values: [0, 7.15633830376646] km
##  Test statistic: Maximum absolute deviation
##  Deviation = observed minus theoretical
## 
## data:  sklepy.ppp.um
## mad = 0.37225, rank = 1, p-value = 0.05
dclf.test(sklepy.ppp.um, Linhom, nsim=19, use.theo=T) # as above
## Generating 19 simulated realisations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.
## 
##  Diggle-Cressie-Loosmore-Ford test of CSR
##  Monte Carlo test based on 19 simulations
##  Summary function: L[inhom](r)
##  Reference function: theoretical
##  Alternative: two.sided
##  Interval of distance values: [0, 7.15633830376646] km
##  Test statistic: Integral of squared absolute deviation
##  Deviation = observed minus theoretical
## 
## data:  sklepy.ppp.um
## u = 0.20983, rank = 1, p-value = 0.05
clarkevans.test(sklepy.ppp.um, alternative = "clustered")
## 
##  Clark-Evans test
##  CDF correction
##  Z-test
## 
## data:  sklepy.ppp.um
## R = 0.79637, p-value = 4.166e-07
## alternative hypothesis: clustered (R < 1)
#hopskel.test(sklepy.ppp.um)

Przeprowadzane analizy wskauzują, że wzorzec jest klastrowany. Przykładowe interpretacje wykresów:

  • K Ripleya powyżej f. Poissonowskiej wskazuje na klastrowanie
  • podobna interpretacja dla PCF i L
  • G powyzej Poissona - mam bliżej do sąsiada niż w przypadku procesu losowego
  • F poniżej Poissona - mam dalej do pustej przestzeni niż w przypadku procesu losowego
  • iloraz J < 1 wskazuje na klastrowanie

Również testy oparte na kopertach (envelopes) oraz te tradycyjne (p-value <= 0.05) wskazują na klastrowanie.

Wzorzec z podziałem na etykiety sieciowe

Analiza pierwszego rzędu

summary(sklepy.ppp.marked)
## Marked planar point pattern:  97 points
## Average intensity 0.1875338 points per square km
## 
## Coordinates are given to 12 decimal places
## 
## Multitype:
##           frequency proportion  intensity
## Biedronka        61   0.628866 0.11793360
## Lidl             36   0.371134 0.06960019
## 
## Window: polygonal boundary
## single connected closed polygon with 4259 vertices
## enclosing rectangle: [1422.9789, 1451.6043] x [4179.58, 4209.699] km
##                      (28.63 x 30.12 km)
## Window area = 517.24 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6
is.multitype(sklepy.ppp.marked)
## [1] TRUE
best_bw <- bw.diggle(sklepy.ppp.marked)
plot(density(split(sklepy.ppp.marked), sigma = best_bw))

plot(relrisk(sklepy.ppp.marked, sigma = best_bw, relative = FALSE)) 

plot(relrisk(sklepy.ppp.marked, sigma = best_bw, relative = TRUE, control = "Biedronka"))

segregation.test.ppp(sklepy.ppp.marked, nsim = 99, verbose = FALSE)
## 
##  Monte Carlo test of spatial segregation of types
## 
## data:  sklepy.ppp.marked
## T = 0.043759, p-value = 0.88

Przy p-value rzędu 0.88 nie mamy podstaw do odrzucenia hipotezy zerowej o braku segregacji między sieciami sklepów.

Analiza drugiego rzędu

plot(alltypes(sklepy.ppp.marked, Kcross))

plot(alltypes(sklepy.ppp.marked, pcfcross))

plot(alltypes(sklepy.ppp.marked, markconnect))

densities_split <- solapply(split(sklepy.ppp.marked), density, sigma = bw.diggle)

K_adjusted <- Kcross.inhom(
  sklepy.ppp.marked, 
  i = "Biedronka", 
  j = "Lidl", 
  lambdaI = densities_split$Biedronka, 
  lambdaJ = densities_split$Lidl
)
plot(K_adjusted)

plot(alltypes(sklepy.ppp.marked, Kdot))

plot(alltypes(sklepy.ppp.marked, Jdot))

# CSRI
aE <- alltypes(sklepy.ppp.marked, Lcross.inhom, nsim = 99, envelope = TRUE, verbose = FALSE, global = TRUE)
plot(aE, . - r ~ r)

# Separate Spatial Profiles
E_shift <- envelope(sklepy.ppp.marked, Lcross.inhom, nsim = 99, i = "Biedronka", j = "Lidl", verbose = FALSE,
                    simulate = expression(rshift(sklepy.ppp.marked, radius = 1.5, edge="erode")))
plot(E_shift, . - r ~ r)

# Complete Random Labeling
E_label <- envelope(sklepy.ppp.marked, Lcross.inhom, nsim = 99, i = "Biedronka", j = "Lidl", verbose = FALSE,
                    simulate = expression(rlabel(sklepy.ppp.marked)))
plot(E_label, . - r ~ r)

Na podstawie wykresu dla K Ripleya i pochodnych, możemy potwierdzić, że Biedronki i Lidle zajmują często podobne lokalizacje.

Wnioski

Z przeprowadzanej analizy wynika, że sklepy jako wzorzec punktowy podlegają klastrowaniu, natomiast sama etykietka (sieć) sklepu jest już niezależną zmienną losową. Powyższe wynioski zgodne są zarówno z teoretycznym modelem Hotellinga, jak i wynikami z cytowanych prac naukowych.

Potencjalne możliwości dalszej rozbudowy niniejszego projektu to: 1) estymacja modelu procesu punktowego (PPM, point process model) z doborem odpowiednich zmiennych zależnych, 2) wykorzystanie procesu Gibbsa, 3) rozbudowa zbioru danych na więcej etykietek lub większy/inny obszar.

Bibliografia

  1. Hotelling, H. (1929). Stability in Competition. The Economic Journal, 39(153), 41–57. https://doi.org/10.2307/2224214
  2. Kunc, J., Križan, F., Novotná, M., Výbošťok, J., Charvátová, K. and Bilková, K. (2025), Grocery Retail Location Patterns in Brno: Clustering, Inequality and Street Network Centrality. Tijds. voor econ. en Soc. Geog.. https://doi.org/10.1111/tesg.70027
  3. Lee, C. (2024). Exploring spatial distribution of retail stores in South Korea. Bulletin of Geography. Socio-economic Series, 66(66):47-55. DOI: http://doi.org/10.12775/bgss-2024-003
  4. TSOUTSOS, Michail-Christos, and Yorgos N. PHOTIS. 2020. “Spatial Point Pattern Analysis of Urban Retail Stores: The Case of Twelve Large- and Medium-Sized Greek Cities”. European Journal of Geography 11 (4). https://doi.org/10.48088/ejg.m.mar.11.4.36.63.
  5. Rabino, Giovanni; Mastrangelo, Lorenzo (2002) : Point pattern analysis: an application to the loyalty networks of chain-stores, 42nd Congress of the European Regional Science Association: “From Industry to Advanced Services - Perspectives of European Metropolitan Regions”, August 27th - 31st, 2002, Dortmund, Germany, European Regional Science Association (ERSA), Louvain-la-Neuve at: https://hdl.handle.net/10419/115827