Spatial analysis
Bogdan Jasiński
2025-06-04
# Load required packages
library(spatstat)## Loading required package: spatstat.data## Loading required package: spatstat.univar## spatstat.univar 3.1-3## Loading required package: spatstat.geom## spatstat.geom 3.4-1## Loading required package: spatstat.random## spatstat.random 3.4-1## Loading required package: spatstat.explore## Loading required package: nlme## spatstat.explore 3.4-3## Loading required package: spatstat.model## Loading required package: rpart## spatstat.model 3.3-6## Loading required package: spatstat.linnet## spatstat.linnet 3.2-6##
## spatstat 3.3-3
## For an introduction to spatstat, type 'beginner'library(sp)
library(sf)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(maps)
# Set up environment
options(scipen = 999)
set.seed(123) Data preparation
dunkin <- read.csv("dunkin_stores.csv", stringsAsFactors = FALSE)
starbucks <- read.csv("directory.csv", stringsAsFactors = FALSE)
dunkin_ma <- dunkin[dunkin$state == "MA", ]
starbucks_ma <- starbucks[starbucks$State.Province == "MA" &
starbucks$Brand == "Starbucks", ]
dunkin_clean <- data.frame(
brand = "Dunkin",
lat = as.numeric(dunkin_ma$loc_lat),
lon = as.numeric(dunkin_ma$loc_long)
)## Warning in data.frame(brand = "Dunkin", lat = as.numeric(dunkin_ma$loc_lat), :
## NAs introduced by coercion
## Warning in data.frame(brand = "Dunkin", lat = as.numeric(dunkin_ma$loc_lat), :
## NAs introduced by coercionstarbucks_clean <- data.frame(
brand = "Starbucks",
lat = starbucks_ma$Latitude,
lon = starbucks_ma$Longitude
)
coffee_data <- rbind(dunkin_clean, starbucks_clean)
coffee_data <- coffee_data[!is.na(coffee_data$lat) & !is.na(coffee_data$lon), ]
coffee_data <- unique(coffee_data)
cat("\nData summary:\n")##
## Data summary:table(coffee_data$brand)##
## Dunkin Starbucks
## 776 219Create spatial objects
ma_map <- map("state", regions = "massachusetts", plot = FALSE, fill = TRUE)
ma_coords <- cbind(ma_map$x, ma_map$y)
ma_coords <- ma_coords[!is.na(ma_coords[,1]), ]
ma_poly <- Polygon(ma_coords)
ma_sp <- SpatialPolygons(list(Polygons(list(ma_poly), "MA")))
proj4string(ma_sp) <- CRS("+proj=longlat +datum=WGS84")
ma_sf <- st_as_sf(ma_sp)
ma_planar <- st_transform(ma_sf, crs = 32619)
coffee_sf <- st_as_sf(coffee_data,
coords = c("lon", "lat"),
crs = 4326) # WGS84
coffee_planar <- st_transform(coffee_sf, crs = 32619)
coffee_clipped <- st_intersection(coffee_planar, ma_planar)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesCreate point pattern object
W <- as.owin(ma_planar)
coords <- st_coordinates(coffee_clipped)
coffee_ppp <- ppp(x = coords[,1],
y = coords[,2],
window = W,
marks = as.factor(coffee_clipped$brand))
coffee_ppp <- rjitter(coffee_ppp, retry = TRUE)
coffee_ppp <- rescale(coffee_ppp, 1000, "km")
print(summary(coffee_ppp))## Marked planar point pattern: 903 points
## Average intensity 0.04291856 points per square km
##
## Coordinates are given to 13 decimal places
##
## Multitype:
## frequency proportion intensity
## Dunkin 712 0.7884828 0.033840550
## Starbucks 191 0.2115172 0.009078013
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 254 20826.50000 0.990000
## polygon 2 16 35.83020 0.001700
## polygon 3 4 2.63848 0.000125
## polygon 4 8 59.68730 0.002840
## polygon 5 30 115.20200 0.005480
## enclosing rectangle: [127.4444, 421.4184] x [4566.069, 4748.639] km
## (294 x 182.6 km)
## Window area = 21039.8 square km
## Unit of length: 1 km
## Fraction of frame area: 0.392Basic visualization
plot(coffee_ppp,
main = "Coffee Shops in Massachusetts",
cols = c("orange", "green"),
pch = 16, cex = 0.8)
First-order analysis
coffee_split <- split(coffee_ppp)
overall_intensity <- intensity(coffee_ppp)
cat("Overall intensity:", round(overall_intensity, 4), "stores per km²\n")## Overall intensity: 0.0338 0.0091 stores per km²brand_intensities <- intensity(coffee_split)
cat("Brand intensities:\n")## Brand intensities:for(i in 1:length(brand_intensities)) {
cat(names(brand_intensities)[i], ":", round(brand_intensities[i], 4), "stores per km²\n")
}## Dunkin : 0.0338 stores per km²
## Starbucks : 0.0091 stores per km²density_overall <- density(coffee_ppp)
plot(density_overall, main = "Overall Store Density")
par(mfrow = c(1, 2))
plot(density(coffee_split$Dunkin), main = "Dunkin' Density")
plot(density(coffee_split$Starbucks), main = "Starbucks Density")
par(mfrow = c(1, 1))
quad_test <- quadrat.test(coffee_ppp, nx = 4, ny = 4)## Warning: Some expected counts are small; chi^2 approximation may be inaccuratecat("Quadrat test p-value:", round(quad_test$p.value, 4), "\n")## Quadrat test p-value: 0if(quad_test$p.value < 0.05) {
cat("Pattern is NOT spatially random (p < 0.05)\n")
} else {
cat("Pattern appears spatially random (p >= 0.05)\n")
}## Pattern is NOT spatially random (p < 0.05)Second-order analysis
K_func <- Kest(coffee_ppp)
plot(K_func, main = "K-function (Overall Pattern)")
L_func <- Lest(coffee_ppp)
plot(L_func, main = "L-function (Overall Pattern)")
K_cross <- Kcross(coffee_ppp, i = "Dunkin", j = "Starbucks")
plot(K_cross, main = "Cross K-function: Dunkin ↔ Starbucks")
L_cross <- Lcross(coffee_ppp, i = "Dunkin", j = "Starbucks")
plot(L_cross, main = "Cross L-function: Dunkin ↔ Starbucks")
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2)
The overall spatial patterns show that the trans function (red dashed line) exhibits the strongest clustering, with values significantly higher than the Poisson random benchmark (blue line), while the other patterns show moderate clustering that increases with distance.
The cross-correlation analysis between Dunkin’ and Starbucks locations revealed a significant spatial relationship between the two coffee chains, with the iso function (black line) showing the strongest cross-clustering.
Statistical testing
# Test for spatial independence between brands
# H0: brands are spatially independent
# H1: brands show spatial dependence
independence_test <- envelope(coffee_ppp, Lcross,
nsim = 39,
i = "Dunkin", j = "Starbucks",
simulate = expression(rlabel(coffee_ppp)),
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(independence_test,
main = "Test for Spatial Independence\n(Random Labeling)")
if(any(independence_test$obs > independence_test$hi |
independence_test$obs < independence_test$lo)) {
cat("REJECT H0: Brands show spatial dependence\n")
} else {
cat("FAIL TO REJECT H0: No evidence of spatial dependence\n")
}## REJECT H0: Brands show spatial dependencecsr_test <- envelope(coffee_ppp, Lest, nsim = 39, global = TRUE)## Generating 39 simulations of CSR ...
## 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.
##
## Done.plot(csr_test, main = "Test against Complete Spatial Randomness")
Dunkin’ and Starbucks show spatial independence in their location patterns. Coffee shop locations are significantly clustered, not randomly distributed.
Simple point process modeling
model_simple <- ppm(coffee_ppp ~ marks)
summary(model_simple)## Point process model
## Fitted to data: coffee_ppp
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.formula(Q = coffee_ppp ~ marks)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Marked planar point pattern: 903 points
## Average intensity 0.0429 points per square km
## Multitype:
## frequency proportion intensity
## Dunkin 712 0.788 0.03380
## Starbucks 191 0.212 0.00908
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 254 20826.50000 0.990000
## polygon 2 16 35.83020 0.001700
## polygon 3 4 2.63848 0.000125
## polygon 4 8 59.68730 0.002840
## polygon 5 30 115.20200 0.005480
## enclosing rectangle: [127.4444, 421.4184] x [4566.069, 4748.639] km
## (294 x 182.6 km)
## Window area = 21039.8 square km
## Unit of length: 1 km
## Fraction of frame area: 0.392
##
## Dummy quadrature points:
## 70 x 70 grid of dummy points, plus 4 corner points
## dummy spacing: 4.199629 x 2.608138 km
##
## Original dummy parameters: =
## Marked planar point pattern: 5041 points
## Average intensity 0.24 points per square km
## Multitype:
## frequency proportion intensity
## Dunkin 2260 0.448 0.107
## Starbucks 2780 0.552 0.132
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 254 20826.50000 0.990000
## polygon 2 16 35.83020 0.001700
## polygon 3 4 2.63848 0.000125
## polygon 4 8 59.68730 0.002840
## polygon 5 30 115.20200 0.005480
## enclosing rectangle: [127.4444, 421.4184] x [4566.069, 4748.639] km
## (294 x 182.6 km)
## Window area = 21039.8 square km
## Unit of length: 1 km
## Fraction of frame area: 0.392
## Quadrature weights:
## (counting weights based on 70 x 70 array of rectangular tiles)
## All weights:
## range: [0.261, 11] total: 41900
## Weights on data points:
## range: [0.261, 5.48] total: 2950
## Weights on dummy points:
## range: [0.261, 11] total: 38900
## --------------------------------------------------------------------------------
## FITTED :
##
## Stationary multitype Poisson process
## Possible marks:
## Dunkin Starbucks
## ---- Intensity: ----
##
## Log intensity: ~marks
##
## Intensities:
## beta_Dunkin beta_Starbucks
## 0.033994787 0.009119389
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -3.381548 0.03747658 -3.455001 -3.308095 *** -90.23096
## marksStarbucks -1.315804 0.08148678 -1.475516 -1.156093 *** -16.14746
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) marksStarbucks
## -3.381548 -1.315804
##
## Fitted exp(theta):
## (Intercept) marksStarbucks
## 0.03399479 0.26825845model_trend <- ppm(coffee_ppp ~ marks + x + y)
summary(model_trend)## Point process model
## Fitted to data: coffee_ppp
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.formula(Q = coffee_ppp ~ marks + x + y)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Marked planar point pattern: 903 points
## Average intensity 0.0429 points per square km
## Multitype:
## frequency proportion intensity
## Dunkin 712 0.788 0.03380
## Starbucks 191 0.212 0.00908
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 254 20826.50000 0.990000
## polygon 2 16 35.83020 0.001700
## polygon 3 4 2.63848 0.000125
## polygon 4 8 59.68730 0.002840
## polygon 5 30 115.20200 0.005480
## enclosing rectangle: [127.4444, 421.4184] x [4566.069, 4748.639] km
## (294 x 182.6 km)
## Window area = 21039.8 square km
## Unit of length: 1 km
## Fraction of frame area: 0.392
##
## Dummy quadrature points:
## 70 x 70 grid of dummy points, plus 4 corner points
## dummy spacing: 4.199629 x 2.608138 km
##
## Original dummy parameters: =
## Marked planar point pattern: 5041 points
## Average intensity 0.24 points per square km
## Multitype:
## frequency proportion intensity
## Dunkin 2260 0.448 0.107
## Starbucks 2780 0.552 0.132
##
## Window: polygonal boundary
## 5 separate polygons (no holes)
## vertices area relative.area
## polygon 1 254 20826.50000 0.990000
## polygon 2 16 35.83020 0.001700
## polygon 3 4 2.63848 0.000125
## polygon 4 8 59.68730 0.002840
## polygon 5 30 115.20200 0.005480
## enclosing rectangle: [127.4444, 421.4184] x [4566.069, 4748.639] km
## (294 x 182.6 km)
## Window area = 21039.8 square km
## Unit of length: 1 km
## Fraction of frame area: 0.392
## Quadrature weights:
## (counting weights based on 70 x 70 array of rectangular tiles)
## All weights:
## range: [0.261, 11] total: 41900
## Weights on data points:
## range: [0.261, 5.48] total: 2950
## Weights on dummy points:
## range: [0.261, 11] total: 38900
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary multitype Poisson process
## Possible marks:
## Dunkin Starbucks
## ---- Intensity: ----
##
## Log intensity: ~marks + x + y
##
## Fitted trend coefficients:
## (Intercept) marksStarbucks x y
## -50.606351170 -1.315804483 0.011267725 0.009393028
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -50.606351170 4.3717695991 -59.17486213 -42.03784021 ***
## marksStarbucks -1.315804483 0.0814867871 -1.47551565 -1.15609332 ***
## x 0.011267725 0.0006140288 0.01006425 0.01247120 ***
## y 0.009393028 0.0009141743 0.00760128 0.01118478 ***
## Zval
## (Intercept) -11.57571
## marksStarbucks -16.14746
## x 18.35048
## y 10.27488
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) marksStarbucks x y
## -50.606351170 -1.315804483 0.011267725 0.009393028
##
## Fitted exp(theta):
## (Intercept) marksStarbucks
## 0.0000000000000000000001051819 0.2682584269663113585302482988
## x y
## 1.0113314449927264249140534957 1.0094372813789493914526929075anova(model_simple, model_trend, test = "LR")## Analysis of Deviance Table
##
## Model 1: ~marks Poisson
## Model 2: ~marks + x + y Poisson
## Npar Df Deviance Pr(>Chi)
## 1 2
## 2 4 2 383.3 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals_plot <- plot(residuals(model_trend, type = "pearson"),
main = "Model Residuals")
The model systematically over-predicts locations in central Massachusetts and under-predicts eastern areas for both coffee chains, suggesting no covariates or spatial effects in the modeling approach.
Summary
cat("\n=== ANALYSIS SUMMARY ===\n")##
## === ANALYSIS SUMMARY ===cat("Total stores analyzed:", npoints(coffee_ppp), "\n")## Total stores analyzed: 903cat("Study area:", round(area(coffee_ppp), 1), "km²\n")## Study area: 21039.8 km²cat("Overall store density:", round(overall_intensity, 3), "stores/km²\n")## Overall store density: 0.034 0.009 stores/km²cat("\nBrand comparison:\n")##
## Brand comparison:cat("- Dunkin':", npoints(coffee_split$Dunkin), "stores\n")## - Dunkin': 712 storescat("- Starbucks:", npoints(coffee_split$Starbucks), "stores\n")## - Starbucks: 191 stores