Analiza lokalizacji automatów paczkowych (APM) w Warszawie przy wykorzystaniu Analizy wzorca punktowego oraz liniowego
Artur Stankiewicz
04.2026 - 06.2026
Część teoretyczna
Punktem wyjścia analizy stanowi opisywany podczas średnio-zaawansowanego kursu mikroekonomii model Hottelinga opisany przez Harolda Hottelinga w 1929 roku. Zgodnie z opisanym przez niego modelem, pomimo tego, że optymalna społecznie oraz “gospodarczo” lokalizacja biznesu wskazuje, że firmy powinny lokować swoje punkty w pewnej odległości między sobą, tak aby maksymalizować liczbę klientów, a tym samym, uzyskiwany zysk. Hotteling w 1929 doszedł jednak do wniosku, że optymalnym rozwiązaniem do firm będzie lokować się w jednym miejscu. Pojawić się może jednak pytanie, to co w przypadku, gdy rozpatrując model Hottelinga pojawi się jeszcze trzeci, czwarty… wniosek będzie analogiczny, pomimo tego, że optimum społeczne przy 3 punktach wskazuje, że optymalne odległości powinny wynosić 1/6, 1/2 oraz 5/6 (Hotteling, 1929).
Podobny wzorzec można ujrzeć na rynku automatów paczkowych (zwanych skrótowo APM). Słynny mem “zebranie paczkomatów w Bydgoszczy” nie bez powodu zrobił furorę w sieci, klienci sami zauważają, że punkty APM są często zebrane w jednym miejscu, a nie “porozrzucane” w przestrzenii. Takie zachowanie jest zgodne z opisywanym modelem Hottelinga, firmy alokując punkty, pośrednio, grają w grę stworzoną przez Hottelinga, choć NIE MUSI to wcale oznaczać, że zachowują się w ten sposób dlatego, że zmusza ich optymalny wybór stworzony przez Hottelinga.
Literatura empiryczna dostarcza ciekawych wniosków dotyczących alokacji punktów APM, wskazując, że zmienne takie jak dochód, dostęp do internetu, rozbudowana infrastruktura drogowa i miejska oraz zmienne demograficzne wpływają na rozmieszczenie punktów APM przez firmy (Schaefer i Figliozzi, 2021). Jak wskazują Schaefer i Figliozzi (2021) w przeprowadzonym badaniu, wykorzystując do tego metodę Warda (ciekawe, metoda z zakresu ML, która ma na celu wykazanie klastrów. Coś w stylu LISA ze Spatial) punkty APM zazwyczaj są rozmieszczone bliżej dróg, często przy sklepach. Co ciekawe, przeprowadzili oni badanie dla Portland, miasta w Stanie Zjednoczonym dla punktów Amazona, który ma całkowicie inne podejście co do alokacji punktów niż jest to w Warszawie, lokując swoje punkty częściej w często odwiedzanych miejscach. Autorzy wskazują, że rynek APM jest przykładem market failure, wykluczając obszary o niższych dochodach i populacji. Osobiście, uważam, że podejście Amazona w tej kwestii jest uzasadnione. Amazon oferuje subskrypcje Amazon Prime, która oferuje darmowe dostawy do domów od minimalnej wartości zamówienia, może patrząc na rynek zrozumieli, że nie ma większego sensu lokować APM w miejscach, w których preferowaną usługą jest dostawa do domu.
Inne podejście, prezentuje Xinwei Kang (2024), który używając metody ML: DBSCAN oraz OLS/GWR wyznaczył determinanty alokacji punktów APM. Nie zastosował on jednak PPA. Wnioski: gęstość sklepów, restauracji oraz przystanków autobusowych pozytywnie wpływa na intensywność punktów APM w danej jednostce terytorialnej. Co ciekawe, im więcej punktów APM, tym mniejsza populacja w danym regionie.
Cel pracy i dlaczego ten temat jest ważny
Z punktu widzenia firm alokujących swoje punkty APM, wykorzystanie analizy wzorca punktowego może prowadzić do kluczowych dla biznesu wniosku, między innymi: w których miejscach mamy mniejszą dominację rynkową, jakie kowariaty są istotne statystycznie, która firma najczęściej lokuje swoje punkty koło naszych punktów oraz jaka jest optymalna trasa dostawcy, minimalizująca odległość.
Celem pracy jest odpowiedzieć na niektóre z powyższych pytań: estymacja modelu przy użyciu kowariat jest jako further research :). Pośrednio, praca stanowi również wprowadzenie do analizy wzorca punktowego dla niewtajemniczonych
Wstęp
Warszawa dysponuje jedną z najgęstszych sieci punktów odbioru paczek w Polsce. Na mapie Warszawy znajduje się 3148 automatów paczkowych (APM) oraz 3223 punktów odbioru w punktach usługowych (PUDO), firm takich jak: InPost, DPD, Orlen, Poczta Polska oraz DHL. Łącznie stanowi to 6371 punktów odbioru paczek (OOH). Wszechobecna sieć automatów paczkowych rodzi kluczowe pytania:
- czy automaty lokalizowane są czysto losowe?
- czy istnieje przestrzenna zależność między automatami?
- czy firmy preferują określone lokalizacje? W których miejscach firmy mają rynkową dominacje nad innymi firmami?
- od czego zależy alokacja automatów?
Na te, oraz wiele innych pytań, odpowiem podczas przygotowanej analizy. Zapraszam do śledzenia RPubs, na którym można zobaczyć od przysłowiowej “kuchni” jak wygląda proces analizy wzorca punktowego na przykładzie automatów paczkowych (APM).
Niniejsza praca służy w celach czysto edukacyjnych w zakresie użycia analizy wzorca punktowego i liniowego (Point Pattern Analysis, PPA) przy użycia języka R. Do analizy wykorzystano publicznie dostępne dane na stronie dostawców (stan na kwiecień 2026).
Wzorce punktowe możemy podzielić na 3 typy:
#Hm... zazwyczaj tytuł pojawia się bezpośrednio nad wykresem... do dalszego dopracowania :). Raportuję użycie AI, nie znałem tej funkcji z pakietu spatstat
library(spatstat)
set.seed(2115)
W<-owin(c(0,2), c(0,2))
par(mfrow=c(1,3), mar=c(1, 1, 1, 1))
plot(rpoispp(lambda=100, win=W), main="Losowy (CSR)", pch=16, cex=1, legend=FALSE)
plot(rSSI(r=0.08, n=100, win=W), main="Regularny", pch=16, cex=1, legend=FALSE)
plot(rMatClust(kappa=10, r=0.1, mu=10, win=W), main="Skupiony", pch=16, cex=1, legend=FALSE)
Kod
Kod przygotowany na podstawie kursu Analiza wzorca punktowego i liniowego w R, prowadzonego przez dr Katerynę Zabarinę, na Wydziale Nauk Ekonomicznych Uniwersytetu Warszawskiego.
Marked oraz unmarked PPP
requiredPackages<-c("sf", "spatstat", "rgugik", "automap",
"magrittr", "sp", "rnaturalearth",
"ggplot2", "tidyverse", "viridis", "osrm")
for(i in requiredPackages)
{if(!require(i, character.only=TRUE)) install.packages(i)}
for(i in requiredPackages)
{library(i, character.only=TRUE)}
Sys.setenv(LANG="en")
windowsFonts(TNR=windowsFont("Times New Roman"))Zacznijmy od wczytania mapy oraz zbioru danych z punktami APM
#Pobieramy mapę Warszawy
warszawa_granica<-suppressWarnings(borders_get(TERYT="1465")) #Pobieramy mapę Warszawy z GUGiK
warszawa_granica<-st_transform(warszawa_granica, crs=2180) #Ustawiamy CRS (układ planetarny) na 2180, czyli odpowiedni dla Warszawy
st_crs(warszawa_granica) #Sprawdzamy czy jest CRS 2180 po zmianie
W<-as.owin(warszawa_granica)
class(W) #Sprawdzamy jaki to obiekt
area_km2<-as.numeric(st_area(warszawa_granica))/1e6 #Powierzchnia Warszawy w km²
rsc<-sqrt(area.owin(W)/area_km2) #Zmieniamy z metrów na kilometry
#Wczytujemy zbiór danych z punktami APM
APM<-read.csv("Lista_APM.csv", fileEncoding="UTF-8-BOM") #Wczytujemy dane z pliku CSV
APM_sf<-st_as_sf(APM, coords=c("lon", "lat"), crs=4326) #Kodujemy koordynaty punktów APM jako długość i szerokość geograficzną
APM_sf<-st_transform(APM_sf, crs=2180) #Przechodzimy z CRS 4326 na 2180
st_crs(APM_sf) #Sprawdzamy czy jest CRS 2180
APM_waw<-st_intersection(APM_sf, warszawa_granica) #Przycinamy dataset do Warszawy## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesSprawdźmy ile mamy automatów paczkowych w zbiorze danych z podziałem na marki
##
## DHL DPD INPOST ORLEN POCZTA POLSKA
## 417 848 1413 419 51p_mapa<-ggplot() +
geom_sf(data=warszawa_granica, fill="grey95", color="black", linewidth=0.6) +
geom_sf(data=APM_waw, aes(color=marka), size=1.5, alpha=0.7) +
scale_color_viridis_d(option="turbo", name="Marka") +
labs(title="Punkty APM w Warszawie", subtitle=paste("Łącznie", nrow(APM_waw), "punktów")
) +
theme_minimal() +
theme(plot.title=element_text(size=14, face="bold"))
print(p_mapa)
Wyświetlmy podsumowanie wzorca unmarked, czyli takiego, który nie różnicuje punktów ze względu na markę
#Unmarked, wszystkie punkty APM bez rozróżnienia marki
APM.ppp.um<-ppp(
x=st_coordinates(APM_waw)[, 1], #Współrzędne X
y=st_coordinates(APM_waw)[, 2], #Współrzędne Y
window=W) #Punkty w Warszawie
APM.ppp.um<-rjitter(APM.ppp.um, 0.05, retry=TRUE) #Przesuwamy lekko punkty żeby nie było duplikatów (o 0.05 metra).
APM.ppp.um<-rescale.ppp(APM.ppp.um, rsc, "km") #Zmieniamy z metrów na kilometry
summary(APM.ppp.um) #Wyświetlmy podsumowanie wzorca unmarked## Planar point pattern: 3148 points
## Average intensity 6.092221 points per square km
##
## Coordinates are given to 13 decimal places
##
## Window: polygonal boundary
## single connected closed polygon with 3389 vertices
## enclosing rectangle: [626.506, 655.2603] x [472.2295, 502.1724] km
## (28.75 x 29.94 km)
## Window area = 516.724 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6Wyświetlmy też podsumowanie wzorca marked, czyli takiego, który różnicuje punkty ze względu na markę
#Marka jako cecha punktów APM
APM.ppp.m<-ppp(
x=st_coordinates(APM_waw)[, 1],
y=st_coordinates(APM_waw)[, 2],
window=W,
marks=as.factor(APM_waw$marka))
APM.ppp.m<-rjitter(APM.ppp.m, 0.5)
APM.ppp.m<-rescale.ppp(APM.ppp.m, rsc, "km")
summary(APM.ppp.m) #Wyświetlamy podsumowanie wzorca## Marked planar point pattern: 3148 points
## Average intensity 6.092221 points per square km
##
## Coordinates are given to 13 decimal places
##
## Multitype:
## frequency proportion intensity
## DHL 417 0.13246510 0.80700650
## DPD 848 0.26937740 1.64110700
## INPOST 1413 0.44885640 2.73453300
## ORLEN 419 0.13310040 0.81087700
## POCZTA POLSKA 51 0.01620076 0.09869863
##
## Window: polygonal boundary
## single connected closed polygon with 3389 vertices
## enclosing rectangle: [626.506, 655.2603] x [472.2295, 502.1724] km
## (28.75 x 29.94 km)
## Window area = 516.724 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6Narysujmy mapę z podziałem na markę
split_ppp_m<-split(APM.ppp.m) #Dzielimy wzorzec marked na osobne ppp dla każdej marki
par(mfrow=c(2,3), family="TNR", cex.main=1.4, cex.sub=1.1)
plot(split_ppp_m[[1]], main="DHL", pch=16, cex=0.5, cols="red")
mtext(paste(npoints(split_ppp_m[[1]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[2]], main="DPD", pch=16, cex=0.5, cols="green")
mtext(paste(npoints(split_ppp_m[[2]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[3]], main="InPost", pch=16, cex=0.5, cols="blue")
mtext(paste(npoints(split_ppp_m[[3]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[4]], main="Orlen", pch=16, cex=0.5, cols="orange")
mtext(paste(npoints(split_ppp_m[[4]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[5]], main="Poczta Polska", pch=16, cex=0.5, cols="black")
mtext(paste(npoints(split_ppp_m[[5]]), "APM"), side=1, line=2.5, cex=1.0)
plot.new() #Puste pole na legendę
legend("center", legend=c("DHL","DPD","InPost","Orlen","Poczta Polska"),
col=c("red","green","blue","orange","black"),
pch=16, cex=1.3, title="Marka", bty="n", pt.cex=1.5)
Funkcje K, L, F, G oraz J. Koperty oraz testy MAD/DCLF/Clark-Evans
Zanim przejdziemy do funkcji… przypomnijmy sobie czym się różni funkcja homogeniczna od niehomogenicznej. Najprościej rzecz ujmując:
- funkcja homogeniczna zakłada, że intensywność punktów w przestrzeni jest taka sama.
- Odwrotnie jest w przypadku funkcji niehomogenicznej, która uwzględnia również występującą intensywność w przestrzeni pomiędzy punktami.
A jak należy interpretować te funkcje? Charakterystyczne są:
- pokrycia z funkcją Poissona, co może sugerować, że mamy do czynienia z CSR
- funkcje nad lub pod funkcją Poissona, co sugeruje, że mamy do czynienia ze wzorcem (skupiony, a może regularny?)
- peaków, przecięć, pokryć z funkcją Poissona od określonej długości promienia
Zwróćcie proszę uwagę na różne osie x oraz y na wykresach
Funkcja K Ripleya mierzy ile punktów wpada do koła o promieniu r wokół losowo wybranego punktu zwiększając stopniowo promień r. Następnie porównujemy z rozkładem losowym (CSR, Poisson)
- K(r) powyżej linii Poisson sugeruje skupienie
- K(r) poniżej linii Poisson sugeruje regularność
Wniosek: funkcja K Ripleya sugeruje, że mamy do czynienia ze skupieniem.
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
K<-Kest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(K, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Kin<-Kinhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Kin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Funkcja L to transformacja funkcji K: L(r)=pierwiastek(K(r)/pi)
- L(r) powyżej linii Poisson sugeruje skupienie
- L(r) poniżej linii Poisson sugeruje regularność
Wniosek: analogiczny do funkcji K Ripleya :)
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
L<-Lest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(L, main="Homogeniczna",
xlab="r (km)", ylab="L(r)", legendpos="bottomright")
Lin<-Linhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Lin, main="Niehomogeniczna",
xlab="r (km)", ylab="L(r)", legendpos="bottomright")
Funkcja F (empty space) to funkcja, która mierzy odległość od pustego w przestrzeni punktu do najbliższego punktu APM.
- F(r) powyżej linii Poisson sugeruje regularność
- F(r) poniżej linii Poisson sugeruje skupienie
Wniosek: F(r) poniżej linii Poissona sugeruje skupienie. Występuje więcej pustej przestrzeni między punktami APM, co jest charakterystyczne dla wzorca skupionego.
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
Ffun<-Fest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Ffun, main="Homogeniczna",
xlab="r (km)", ylab="F(r)", legendpos="bottomright")
Fin<-Finhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Fin, main="Niehomogeniczna",
xlab="r (km)", ylab="F(r)", legendpos="bottomright")
Funkcja G (nearest neighbour) to funkcja, która mierzy odległość od punktu APM w przestrzeni do najbliższego punktu APM.
- G(r) powyżej linii Poisson sugeruje skupienie
- G(r) poniżej linii Poisson sugeruje regularność
Wniosek: sugeruje, że punkty stoją bliżej siebie niż gdyby to było w przypadku rozkładu Poissona
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
G<-Gest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(G, main="Homogeniczna",
xlab="r (km)", ylab="G(r)", legendpos="bottomright")
Gin<-Ginhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Gin, main="Niehomogeniczna",
xlab="r (km)", ylab="G(r)", legendpos="bottomright")
Funkcja J służy do porównania czy mamy bardziej do czynienia z pustą przestrzenią (empty space, funkcja F) czy najbliższego sąsiada (nearest neighbour, funkcja G)
- J(r) powyżej 1 wskazuje, że silniejszy jest efekt funkcji F (regularność)
- J(r) poniżej 1 wskazuje, że silniejszy jest efekt funkcji G (skupienie)
Wniosek: zdecydowanie silniejszy efekt funkcji G względem funkcji F
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
J<-Jest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(J, main="Homogeniczna",
xlab="r (km)", ylab="J(r)", legendpos="bottomright")
Jin<-Jinhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Jin, main="Niehomogeniczna",
xlab="r (km)", ylab="J(r)", legendpos="bottomright")
Narysowaliśmy funkcje i wstępnie postawiliśmy wnioski. Dla potwierdzenia wniosków posłużmy się innymi testami. Zacznijmy od kopert.
W kopertach interesuje nas czy empiryczna funkcja, w naszym przypadku funkcja G, wychodzi poza szary obszar koperty. Szary zakres wynika z otrzymanych wartości minimum oraz maksimum funkcji CSR w drodze przeprowadzonej symulacji Monte Carlo. Jeśli funkcja empiryczna:
- wychodzi poza szary obszar koperty: istnieją podstawy do odrzucenia hipotezy zerowej o CSR
- nie wychodzi poza szary obszar koperty: brak podstaw do odrzucenia hipotezy zerowej o CSR
Funkcja kopert pozwala nam na wybranie następujących atrybutów:
- nsim liczba symulacji Monte Carlo. Wybór nsim będzie potem definiował p-value. Dla kopert punktowych: 2/(nsim+1), dla kopert globalnych: 1/(nsim+1). Zazwyczaj wybieramy nsim spośród 19, 39, 99 (w zależności jaki poziom istotności nas interesuje. Im wyższe nsim tym lepiej)
- fix.n liczba punktów ma być zgodna z liczbą punktów w naszym zbiorze danych
- correction dla ponad 3000 punktów automatycznie ustawione na “border”
- r możemy sami wyznaczyć maksymalny promień oraz co ile ma się on zmieniać
- global koperta globalna (TRUE) lub punktowa (FALSE)
- typ funkcji możemy wybrać czy funkcja ma być homogeniczna/niehomogeniczna oraz typ funkcji
Wyróżniamy dwa typy kopert:
Koperta punktowa, która oblicza kopertę dla funkcji G pozwala określić w jakich odległościach dokładnie punkty APM się “grupują” (lub wręcz odwrotnie)
set.seed(2115)
n_pacz<-npoints(APM.ppp.um)
E.ptwise<-envelope(APM.ppp.um, Gest,
nsim=39, #p-value = 5%
fix.n=TRUE,
correction="border",
rmax=3,
r=seq(0, 3, by=0.05),
verbose=FALSE)
par(family="TNR", bty="l")
plot(E.ptwise, main=("Koperta punktowa"),
xlab="r (km)", ylab="G(r)", legendpos="bottomright")
Koperta globalna, która oblicza kopertę dla funkcji G wybiera jeden maksymalną oraz minimalną wartość z symulacji Monte Carlo i stosuje ją dla całej długości promienia r. Jest po prostu bardziej restrykcyjna względem koperty punktowej
set.seed(2115)
E.global<-envelope(APM.ppp.um, Gest,
nsim=19, #p-value = 5%
rank=1,
global=TRUE,
correction="border",
fix.n=TRUE,
rmax=3,
r=seq(0, 3, by=0.05),
verbose=FALSE)
par(family="TNR", bty="l")
plot(E.global, main=paste0("Koperta globalna"),
xlab="r (km)", ylab="G(r)", legendpos="bottomright")
Wniosek z obu kopert jest jasny: punkty APM rozmieszczone są w przestrzeni zgodnie ze skupionym wzorcem.
Możemy też przeprowadzić test MAD i DCLF, które potwierdzają wnioski z kopert. Niestety, przy tak dużej liczbie punktów test się nie wykona
mad.test(APM.ppp.um, Linhom, nsim=19, verbose=FALSE)
dclf.test(APM.ppp.um, Linhom, nsim=19, verbose=FALSE)Teraz testy Clark-Evans. Postanowiłem sprawdzić dwie hipotezy alternatywne:
- Hipoteza zerowa H0: punkty są zgodne z czysto losowym wzorcem (CSR)
- H1(a): punkty są zgodne ze skupionym wzorcem
- H1(b): punkty są zgodne z regularnym wzorcem
Na podstawie poniższego testu istnieją podstawy do odrzucenia hipotezy zerowej o CSR na rzecz hipotezy alternatywnej H1(a) o skupieniu (p-value=2,2e-16)
##
## Clark-Evans test
## CDF correction
## Z-test
##
## data: APM.ppp.um
## R = 0.48321, p-value < 2.2e-16
## alternative hypothesis: clustered (R < 1)Na podstawie poniższego testu nie ma podstaw do odrzucenia hipotezy zerowej o CSR na rzecz hipotezy alternatywnej H1(b) o regularności (p-value=1)
##
## Clark-Evans test
## CDF correction
## Z-test
##
## data: APM.ppp.um
## R = 0.48321, p-value = 1
## alternative hypothesis: regular (R > 1)Wniosek: testy wskazują jednoznacznie: punkty APM nie są rozmieszczone losowo, a zgodnie ze skupionym wzorcem.
Intensywność oraz relatywne ryzyko
Intensywność to najprościej rzecz ujmując liczba punktów przypadająca na daną jednostkę przestrzenną (na przykład kilometr kwadratowy). Intensywność możemy przedstawić na dwa sposoby:
1. podzielić mapę Warszawy na tak zwany “grid”, czyli siatkę i przypisać każdej siatce punkty APM.
#Dzielimy mapę Warszawy na GRID, czyli siatkę 9x9
qdr<-quadratcount(APM.ppp.um, nx=9, ny=9)
par(family="TNR", bty="l")
plot(qdr, main="Liczba punktów APM w kwadratach")
a następnie stworzyć mapę ciepła (punkty na km^2).
#Intensywność w każdym polu (punkty/km²)
par(family="TNR", bty="l")
plot(intensity(qdr, image=TRUE), main="Intensywność punktów APM (punkty/km^2)")
2. narysować mapę ciepła używając do tego KDE (Kernel Density Estimation). Zacznijmy od wyboru sigmy. W praktyce, wartość sigmy oznacza rozmycie punktu w przestrzeni. Można wybrać dowolną wartość. Narysujmy 6 przykładowych sigm i zobaczmy jakie ma to przełożenie na mapę intensywności. Przy:
- Małej sigmie rozmycie jest małe, widać lokalne punkty
- Dużej sigmie rozmycie jest duże, widać bardziej globalny trend
Co wybrać? To zależy na czym chcemy się dokładnie skupić: globalne trendy czy lokalna sytuacja?
Zwróćcie proszę uwagę na różne przedziały legendy
sigmy<-c(0.1, 0.5, 0.7, 1, 1.5, 2)
par(mfrow=c(2,3), mar=c(1,1,2,1), family="TNR")
for(s in sigmy)
{plot(density(APM.ppp.um, sigma=s),
main=paste0("Sigma=", s,"km"),
ribbon=TRUE)}
Postanowiłem wybrać sigmę=1. Oczywiście, istnieje również funkcja bw.diggle, która wybierze optymalną sigmę automatycznie, problemem jest jedynie fakt, że mapy pomiędzy firmami stają się wtedy nieporównywalne. Sprawdźmy jaki poziom sigm proponuje bw.diggle:
## DHL DPD INPOST ORLEN POCZTA POLSKA
## 0.7104178 0.2532182 0.4149965 0.7104178 1.7162568Dlatego dla ujednolicenia map, będę posługiwał się jedną wspólną sigmą.
Narysujmy teraz mapy intensywności dla punktów APM łącznie oraz z podziałem na marki.
Zwróćcie proszę uwagę na różne przedziały legendy
sigma<-1
par(mfrow=c(2,3), mar=c(1,1,2,1), family="TNR")
plot(density(APM.ppp.um, sigma=sigma), main="Łącznie", ribbon=TRUE)
for(i in 1:5) {
plot(density(split_ppp_m[[i]], sigma=sigma),
main=names(split_ppp_m)[i],
ribbon=TRUE)
}
I teraz najciekawsze: relatywne ryzyko. Przejdziemy przez każdy sposób wyrażenia relatywnego ryzyka żeby lepiej zrozumieć na czym tak właściwie to polega oraz jak należy interpretować otrzymany wynik. Przejdziemy przez 3 warianty, zaczynając od najprostszego idąc do tego najogólniejszego. Wszystko zależy od tego jakie parametry funkcji relrisk z pakietu Spatstat wybierzemy, możemy wybrać:
- relative. Wartość TRUE obliczy nam iloraz
intensywności jednej firmy względem drugiej. Wartości w tym przypadku
mogą przyjmować wartości od 0 do nieskończoności.
Wartość FALSE obliczy nam prawdopodobieństwo dzieląc intensywność jednej firmy względem sumy intensywności obu firm. Przyjmuje wartości od 0 do 1. - case/control. Określa, którą firmę wybierzemy jako tą znajdującą się w liczniku (case) lub mianowniku (control).
- przykładowo: dla firm InPost (case) oraz DPD (control) będzie to liczone tak:
FALSE Intensywność InPost / (Intensywność InPost +
Intensywność DPD)
TRUE Intensywność InPost / Intensywność DPD
#Musimy na początek wyznaczyć oddzielny subset, w którym będą znajdowały się tylko dwie marki
APM_sub2<-APM.ppp.m[marks(APM.ppp.m) %in% c("INPOST", "DPD")]
marks(APM_sub2)<-droplevels(marks(APM_sub2))Wariant 1
Ile razy częściej w danym miejscu występuje DPD względem InPost? W tym przypadku atrybuty funkcji to relative=TRUE oraz control=“INPOST”. Otrzymujemy następującą mapę:
rr_1<-relrisk(APM_sub2, sigma=sigma, relative=TRUE, control="INPOST")
par(family="TNR")
plot(rr_1, main="Int_DPD / Int_InPost")
Na podstawie mapy widzimy, że na północnym zachodzie Warszawy są obszary, gdzie intensywność DPD jest 1,5 razy większa niż InPost.
Wariant 2
Prawdopodobieństwo, że napotkany automat w danym miejscu należy do DPD (spośród automatów InPost oraz DPD). W tym przypadku mamy atrybut relative=FALSE, czyli w mianowniku otrzymamy sumę intensywności InPost + DPD. Zwróćcie uwagę na parametr case=“DPD”, to on sprawia, że DPD trafia do licznika.
rr_2<-relrisk(APM_sub2, sigma=sigma, relative=FALSE, case="DPD")
par(family="TNR")
plot(rr_2, main="Int_DPD / (Int_DPD + Int_InPost)")
Wniosek? Dokładnie ten sam jak poprzednio, zmieniła się tylko skala. Zależało mi na przedstawieniu czym różni się funkcja relative ustawiona na TRUE względem FALSE, przy dwóch firmach zmienia się jedynie skala. W przypadku Wariantu 2 prawdopodobieństwo napotkania w tym miejscu automatu DPD wynosi około 60-70%.
Wariant 3
W wariancie trzecim mamy więcej marek niż dwie. Funkcja relrisk narysuje nam tym razem 5 map (bo tyle mamy marek w zbiorze). Możemy interpretować wynik jako dominację rynkową.
- Przykładowo: czy firma posiada wyższą dominację w danym regionie niż w innych, w których posiada swoje APM?
- Odpowiedź: DHL ma wyższą dominację rynkową na obrzeżach miasta niż w centrum
Zwróćcie proszę uwagę na różne przedziały legendy
rr_3<-relrisk(APM.ppp.m, sigma=sigma)
par(family="TNR")
plot(rr_3, main="Int_Marka / Suma intensywności wszystkich APM")
Analiza odległości, segregacja oraz dystans z punktu widzenia konsumenta
Test segregacji przestrzennej
Test segregacji sprawdza czy marki APM rozmieszczone są losowo względem siebie (H0: losowe etykietowanie), czy preferują określone lokalizacje w przestrzeni.
##
## Monte Carlo test of spatial segregation of types
##
## data: APM.ppp.m
## T = 10.828, p-value = 0.05Wniosek:
Odległości do najbliższego sąsiada
Macierz zwraca odległość do najbliższego sąsiada z podziałem na markę
## DHL DPD INPOST ORLEN POCZTA POLSKA
## [1,] 0.2113222 0.26528205 0.2165693 0.25166230 0.7374301
## [2,] 0.3322191 0.13976803 0.3199745 0.19682665 0.7018047
## [3,] 0.6510884 0.64815109 0.6182601 0.64329991 1.1938274
## [4,] 0.3643505 0.38707012 0.3601018 0.69089334 1.7101166
## [5,] 1.1173159 0.35250568 0.1585622 0.13729013 2.1040918
## [6,] 0.0730599 0.01460516 0.1732513 0.01121624 0.5936821## DHL DPD INPOST ORLEN POCZTA POLSKA
## [1,] 2446 2208 751 2812 3135
## [2,] 2530 2189 530 2992 3125
## [3,] 2495 1893 665 3044 3144
## [4,] 2548 2068 931 2784 3117
## [5,] 2545 1886 87 2747 3126
## [6,] 2558 1570 82 2993 3106nncorr
- unnormalised prawdopodobieństwo, że APM i jego sąsiad należą do tej samej marki
- normalised unnormalised podzielone przez wartość oczekiwaną. (1=brak korelacji), poniżej 1 oznacza, że firmy wolą nie lokować blisko siebie swoich punktów APM
## unnormalised normalised
## 0.2291200 0.7401439Tabela sąsiedztwa, kto z kim najczęściej sąsiaduje
## neighbour
## point DHL DPD INPOST ORLEN POCZTA POLSKA
## DHL 4 145 111 140 17
## DPD 168 133 366 167 14
## INPOST 192 429 580 193 19
## ORLEN 141 143 123 8 4
## POCZTA POLSKA 21 12 12 6 0## DHL DPD INPOST ORLEN POCZTA POLSKA
## 417 848 1413 419 51Dystans z punktu widzenia dostawcy (OSRM)
library(leaflet)
#Sortownia InPost (Warszawa)
sortownia_sf<-st_as_sf(
data.frame(name='Sortownia InPost',
lat=52.3110908105985,
lon=21.01836853299125),
coords=c('lon','lat'), crs=4326)
#6 automatów paczkowych InPost w okolicach Woli
paczkomaty_df<-data.frame(
name=paste0('APM_', 1:6),
lat=c(52.310627066358265, 52.3069925344919, 52.30505075049821,
52.304027343587585, 52.314837524287334, 52.31743472285266),
lon=c(20.969839717946318, 20.978540768513582, 20.984677662732807,
20.97845493798414, 20.9656232499619, 20.96785484782913))
paczkomaty_sf<-st_as_sf(paczkomaty_df, coords=c('lon','lat'), crs=4326)
#Wszystkie przystanki razem (sortownia + 6 APM)
all_stops<-rbind(sortownia_sf[, 'name'], paczkomaty_sf[, 'name'])
#Optymalna trasa kuriera (OSRM Trip = TSP)
trips<-osrmTrip(loc=all_stops, overview='full')
mytrip<-trips[[1]]$trip
#Mapa interaktywna (leaflet)
leaflet() %>%
addTiles() %>%
addPolylines(data=mytrip, color='blue', weight=3,
opacity=0.8, label='Trasa kuriera') %>%
addMarkers(data=sortownia_sf,
popup='Sortownia InPost',
label='Sortownia') %>%
addCircleMarkers(data=paczkomaty_sf,
color='red', radius=8, fillOpacity=0.9,
popup=paczkomaty_df$name,
label=paczkomaty_df$name)Wnioski
Na podstawie przeprowadzonej analizy, można wyciągnąć następujące wnioski * firmy alokują punkty bliżej siebie niż wynikałoby to z CSR (patrz, funkcja G i koperta dla G). Istnieją również testy statystyczne potwierdzające nasze wnioski * chcąc przedstawić mapę intensywności lub relatywnego ryzyka musimy posłużyć się sigmą, która w praktyce oznacza stopień rozmycia punktu * biznesowo: DHL ma wyższą dominację rynkową na obrzeżach miasta niż w centrum * segregation test wykazał: p-value = 0.05. Przydałoby się zwiększyć liczbę symulacji. Przyjmując poziom istotności = 10%, odrzucając H0 i wskazując, że firmy wolą pewne obszary. * możemy narysować interaktywną mapę, która wskaże optymalną trasę przez wyznaczone punkty