final exam

Part 1: 不同行政區的速食店有不同的地理分布特徵嗎？(問答題)

(1) 圖 1是某城市各村里速食店密度的Moran Scatter Plot。簡述該圖表的 X軸與 Y軸所表示的意涵，以及其空間型態。（4分）

Moran Scatter Plot 用於分析空間自相關性。X 軸（Standardized Variable） 表示的是每個空間單元（在此為村里）的速食店密度原始觀測值（Variable X）。Y 軸（Standardized Lagged Variable） 則表示每個村里鄰近區域（空間鄰近單位）的標準化平均速食店密度，即空間上的鄰近值。根據圖 1 的散點分布和趨勢線，其空間型態顯示出顯著的正向自相關 (PositiveAutocorrelation)。這意味著速食店密度高的村里往往聚集在其他速食店密度高的村里附近 (High-High 聚集)，而速食店密度低的村里也傾向於聚集在其他速食店密度低的村里附近 (Low-Low 聚集)。

(2) 請簡述該圖趨勢線(Fit)的斜率值 (slope =0.982)，所代表的意涵。（3分）

在 Moran Scatter Plot 中，趨勢線（Fit line）的斜率實際上就是 Global Moran’s I 統計量 的值。該圖的趨勢線斜率為 0.982，非常接近 1，顯示速食店密度的空間分布具有正向自相關性。這表示：• 高速食店密度的村里，其周圍鄰近區域也傾向於擁有高密度；• 低密度區域則大多與其他低密度區域相鄰。因此，這種分布型態呈現出明顯的聚集（clustering）現象，即類似的值（不論高或低）會彼此聚集，而不是隨機或分散分布。

(3) 關於該城市各村里速食店密度的 Local Moran’s I統計量，請問其統計量可能會大於 0、等於 0或小於 0，並說明判斷的理由。（3分，理由正確才會給分）

根據圖 1 所示，Global Moran’s I 值為 0.982，表明整體上存在強烈的正向自相關或聚集。在 Moran Scatter Plot 中，這表示大多數散點落在 H-H (右上) 和 L-L (左下) 兩個象限。H-H 和 L-L 區域都代表了相似值的聚集。 因此，對於該城市各村里速食店密度的 Local Moran’s I 統計量，大多數村里的統計量可能會大於 0。這是因為 Global Moran’s I 的強烈正向值表示資料主要由 H-H 和 L-L 這兩類相似值聚集的區域所主導，而 Local Moran’s I 會因為這些聚集區域計算出較大的正值。圖中顯示明顯的同質性聚集現象（正向自相關），多數村里與鄰近地區的速食店密度趨勢一致。

Part 1:實作題（30分，每題10分）

1. 假設速食店的服務範圍是 1 公里，比較 A 區(文山+大安+中正) 與 B 區(信義+南港+松山) 這兩個地區的每一家速食店在服務可及範圍內，所涵蓋學校數量的平均值，是否有統計顯著差異。(需列出虛無假設與對立假設，統計檢定量，以及檢定的顯著水準等)。H0：μA 和 μB 具有顯著差異；H1：μA 和 μB 不具有顯著差異；α = 0.05

FF_joined <- st_join(food, vill["TOWN"])
TPV_A = dplyr::filter(FF_joined, TOWN %in% c("文山區", "大安區", "中正區"))
TPV_B = dplyr::filter(FF_joined, TOWN %in% c("信義區", "南港區", "松山區"))
buffer_A <- st_buffer(TPV_A, dist = 1000)
buffer_B <- st_buffer(TPV_B, dist = 1000)

count_schools <- function(buffers, schools) {
  sapply(1:nrow(buffers), function(i) {
    sum(st_intersects(buffers[i, ], schools, sparse = FALSE))
  })
}

schools_in_A <- count_schools(buffer_A, school)
schools_in_B <- count_schools(buffer_B, school)

buffer_A$school_count <- schools_in_A
buffer_B$school_count <- schools_in_B
t_test_result <- t.test(schools_in_A, schools_in_B, var.equal = FALSE)
t_test_result

    Welch Two Sample t-test

data:  schools_in_A and schools_in_B
t = 1.5977, df = 51, p-value = 0.1163
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1418044  1.2473762
sample estimates:
mean of x mean of y 
 3.870968  3.318182 

p-value = 0.1163 > 0.05，因此拒絕虛無假設(H0)，μA 和 μB 不具有顯著差異。

2. 以 500 公尺方格的空間單位建立網格，依照 contiguity 鄰近定義，計算並繪製速食店總數的 Moran’s correlogram，解讀其圖表資訊，並評估速食店的空間影響範圍。

taipei_boundary <- st_union(vill)

grid_500 <- st_make_grid(taipei_boundary, cellsize = 500, square = TRUE) %>% st_sf() %>%
  st_intersection(taipei_boundary)
grid_500$fastfood_count <- lengths(st_intersects(grid_500, food))

grid_sp <- as(grid_500, "Spatial")

nb <- poly2nb(grid_sp, queen = TRUE)

lw <- nb2listw(nb, style = "W")
grid_500$jittered_count <- jitter(grid_500$fastfood_count, amount = 0.1)
moran_corr <- sp.correlogram(nb, grid_500$jittered_count, order = 10, method = "I", style = "W", zero.policy = TRUE)

moran_df <- data.frame(
  order = 1:length(moran_corr$res),
  moran_I = sapply(moran_corr$res, function(x) x[1])
)
ggplot(moran_df, aes(x = order, y = moran_I)) +
  geom_line(color = "blue", size = 1) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(
    title = "Moran's I Correlogram for Fast Food Count (500m Grid)",
    x = "Order of Contiguity",
    y = "Moran's I"
  ) +
  scale_y_continuous(breaks = seq(-0.001, 0.2, by = 0.01)) + 
  theme_minimal()

在order ≦ 10，moran’s I皆大於0。因此，在大約5公里(500m*10)之前的範圍內，速食店呈現群聚分布，具空間自相關現象，而隨著距離越來越遠，moran’s I 趨近於0，速食店呈現隨機分布，無空間自相關現象。

3. 依照前一小題的網格為空間單位，依照 contiguity 鄰近定義，計算 Gi*統計量，分別繪製速食店家數與學校在 0.05 顯著水準的熱區分佈（需進行多重檢定的修正）。

grid_500m <- st_make_grid(taipei_boundary, cellsize = 500, square = TRUE) %>% st_sf() %>% 
  st_intersection(taipei_boundary)
grid_500m$grid_id <- 1:nrow(grid_500m)

fastfood_join <- st_join(food, grid_500m, left = FALSE)
fastfood_count <- fastfood_join %>%
  group_by(grid_id) %>%
  summarise(fastfood_n = n())

school_join <- st_join(school, grid_500m, left = FALSE)
school_count <- school_join %>%
  group_by(grid_id) %>%
  summarise(school_n = n())

grid_500m <- st_join(grid_500m, fastfood_count, by = "grid_id")
grid_500m <- st_join(grid_500m, school_count, by = "grid_id")
grid_500m$fastfood_n[is.na(grid_500m$fastfood_n)] <- 0
grid_500m$school_n[is.na(grid_500m$school_n)] <- 0
nb <- poly2nb(grid_500m, queen = TRUE)
lw <- nb2listw(nb, style = "W")
gi_fast <- localG(grid_500m$fastfood_n, lw)
gi_school <- localG(grid_500m$school_n, lw)

grid_500m$gi_fast <- as.numeric(gi_fast)
grid_500m$gi_school <- as.numeric(gi_school)

grid_500m$p_fast <- 2 * pnorm(-abs(grid_500m$gi_fast))
grid_500m$p_fast_adj <- p.adjust(grid_500m$p_fast, method = "fdr")
grid_500m$hotspot_fast <- ifelse(grid_500m$p_fast_adj < 0.05 & grid_500m$gi_fast > 0, "Hotspot",
                           ifelse(grid_500m$p_fast_adj < 0.05 & grid_500m$gi_fast < 0, "Coldspot", "Not Significant"))

grid_500m$p_school <- 2 * pnorm(-abs(grid_500m$gi_school))
grid_500m$p_school_adj <- p.adjust(grid_500m$p_school, method = "fdr")
grid_500m$hotspot_school <- ifelse(grid_500m$p_school_adj < 0.05 & grid_500m$gi_school > 0, "Hotspot",
                             ifelse(grid_500m$p_school_adj < 0.05 & grid_500m$gi_school < 0, "Coldspot", "Not Significant"))
tmap_mode("plot")

tm_shape(grid_500m) +
  tm_polygons(
    col = "hotspot_fast", 
    palette = c("red", "grey"),
    legend.show = TRUE,
    title = "Fast Food Hotspots (Gi*)"
  ) +
  tm_borders() +
  tm_layout(title = "速食店熱區 (Gi*) - FDR 修正", legend.outside = TRUE)

tm_shape(grid_500m) +
  tm_polygons(
    col = "hotspot_school", 
    palette = c("blue", "grey"),
    legend.show = TRUE,
    title = "School Hotspots (Gi*)"
  ) +
  tm_borders() +
  tm_layout(title = "學校熱區 (Gi*) - FDR 修正", legend.outside = TRUE)

Part 2：學校附近的速食店數量真的比較多嗎?

1. 請參考Transportation Research Part A, 45 (2011) 640–652（HW-09 的研讀教材）所使 用F(d)函數的作法，比較分別位於大安區、文山區與信義區的學校，到鄰近速食店的空間特徵。

target_vill <- vill %>% filter(TOWN %in% c("大安區", "文山區", "信義區"))
school_sel <- st_join(school, target_vill, join = st_within) %>% filter(!is.na(TOWN))
food_sel <- st_join(food, target_vill, join = st_within) %>% filter(!is.na(TOWN))

sf_to_ppp <- function(points_sf, window_polygon) {
  win <- as.owin(st_union(window_polygon))
  
  coords <- st_coordinates(points_sf)
  
  ppp_obj <- ppp(x = coords[,1], y = coords[,2], window = win)
  return(ppp_obj)
}

districts <- c("大安區", "文山區", "信義區")

fd_results <- list()

for (dist in districts) {
  dist_vill <- target_vill %>% filter(TOWN == dist)
  dist_school <- school_sel %>% filter(TOWN == dist)
  dist_food <- food_sel %>% filter(TOWN == dist)
  
  food_ppp <- sf_to_ppp(dist_food, dist_vill)
  
  school_coords <- st_coordinates(dist_school)
  marks <- rep(1, nrow(school_coords))
  school_ppp <- ppp(school_coords[,1], school_coords[,2], window = food_ppp$window, marks = marks)
  
  fd <- Fest(food_ppp, correction = "border")
  
  fd_results[[dist]] <- fd
}

plot(fd_results[["大安區"]], main = "F(d) function(大安,信義,文山)", col = "red")
plot(fd_results[["文山區"]], add = TRUE, col = "blue")

plot(fd_results[["信義區"]], add = TRUE, col = "green")
legend("bottomright", legend = districts, col = c("red", "blue", "green"), lty = 1)

par(mfrow = c(1, 3))  # 三區並排
for (dist in districts) {
  dist_vill <- target_vill %>% filter(TOWN == dist)
  dist_school <- school_sel %>% filter(TOWN == dist)
  dist_food <- food_sel %>% filter(TOWN == dist)

  win <- as.owin(st_union(dist_vill))
  school_coords <- st_coordinates(dist_school)
  food_coords <- st_coordinates(dist_food)

  school_ppp <- ppp(school_coords[,1], school_coords[,2], window = win)
  food_ppp   <- ppp(food_coords[,1], food_coords[,2], window = win)

  plot(school_ppp, main = paste(dist, "學校 vs 速食店"), cols = "blue", pch = 3)
  points(food_ppp, col = "red", pch = 16)
  legend("topright", legend = c("學校", "速食店"), col = c("blue", "red"), pch = c(3, 16))
}

districts <- c("大安區", "信義區", "文山區")

for (dist in districts) {
  
  dist_school <- school_sel %>% filter(TOWN == dist)
  dist_food   <- food_sel   %>% filter(TOWN == dist)

  distances <- st_distance(dist_school, dist_food) %>%
    apply(1, min) %>%
    as.numeric()

  d_values <- seq(0, 500, by = 10)
  F_empirical <- sapply(d_values, function(d) mean(distances <= d))

  n_random_points <- 500
  n_samples <- nrow(dist_school)
  n_sim <- 99

  bbox <- st_bbox(dist_school)
  random_points <- st_as_sf(
    data.frame(
      x = runif(n_random_points, bbox["xmin"], bbox["xmax"]),
      y = runif(n_random_points, bbox["ymin"], bbox["ymax"])
    ),
    coords = c("x", "y"),
    crs = st_crs(dist_school)
  )

  dist_food_proj <- st_transform(dist_food, 32618)
  random_points_proj <- st_transform(random_points, 32618)

  sim_results <- matrix(nrow = length(d_values), ncol = n_sim)

  for (i in 1:n_sim) {
    sampled_points <- random_points_proj %>%
      sample_n(n_samples)

    sim_distances <- st_distance(sampled_points, dist_food_proj) %>%
      apply(1, min) %>%
      as.numeric()

    sim_results[, i] <- sapply(d_values, function(d) mean(sim_distances <= d))
  }

  envelope_lower <- apply(sim_results, 1, min)
  envelope_upper <- apply(sim_results, 1, max)

  plot_data <- data.frame(
    d = d_values,
    Empirical = F_empirical,
    Lower = envelope_lower,
    Upper = envelope_upper
  )
  
  p <- ggplot(plot_data, aes(x = d)) +
    geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = "Simulation Envelope"),
                alpha = 0.7) +
    geom_line(aes(y = Empirical, color = "Empirical F(d)"), linewidth = 1) +
    scale_fill_manual(name = NULL, values = c("Simulation Envelope" = "grey80")) +
    scale_color_manual(name = NULL, values = c("Empirical F(d)" = "#D55E00")) +
    labs(
      x = "距離 (公尺)",
      y = "F(d)",
      title = paste0("F(d): ", dist, " 學校到速食店的可近性分析")
    ) +
    theme_minimal() +
    theme(
      panel.grid.minor = element_blank(),
      plot.title = element_text(hjust = 0.5, size = 14),
      legend.position = "bottom",
      legend.text = element_text(size = 10)
    )

  print(p)
} 

雖然從confidence envelope來看，三區皆呈現隨機分布，但從f-function的曲線可以看到大安區學校與速食店較偏向cluster分布,信義區學校與速食店較偏向random分布,文山區學校與速食店較偏向dispersion分布，顯現出大安區附近有較多速食店，反映三區在學校與速食店實際距離關係上的結構性差異，這些差異可能來自地理面積、建成環境密度等因素

2. 利用 Ripley’s K function, k(d)計算速食店在 distance (d) = 3000 公尺時的 K(3000) 以及 L(3000)，並使用 Monte Carlo significance test 計算 L(3000)的 95%信賴區間

food_proj <- st_transform(food, 3826)
vill_proj <- st_transform(vill, 3826)

taipei_win <- as.owin(st_union(vill_proj))
coords <- st_coordinates(food_proj)
food_ppp <- ppp(x = coords[,1], y = coords[,2], window = taipei_win)

K_result <- Kest(food_ppp, correction = "border")
K3000 <- K_result$border[which.min(abs(K_result$r - 3000))]
L3000 <- sqrt(K3000 / pi) - 3000

cat("K(3000) =", K3000, "\n")

K(3000) = 67640147 

cat("L(3000) =", L3000, "\n")

L(3000) = 1640.1 

taipei_win <- as.owin(st_union(vill))

set.seed(123)
env <- envelope(
  food_ppp,
  fun = Kest,
  nsim = 99,
  nrank = 5,
  correction = "border",
  savefuns = TRUE
)

Generating 99 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, 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.

env_L <- eval.fv(sqrt(env / pi)-3000)

L3000_lo <- env_L$lo[which.min(abs(env_L$r - 3000))]
L3000_hi <- env_L$hi[which.min(abs(env_L$r - 3000))]

cat("L(3000) 95% 信賴區間：", L3000_lo, "~", L3000_hi, "\n")

L(3000) 95% 信賴區間： -306.9529 ~ 282.8214 

plot(env_L, main = "L(d) Confidence Envelope")
abline(v = 3000, col = "blue", lty = 2)

從R-consule的計算結果可以看到L(3000)=1640.1遠大於信賴區間上限282.821，因此顯示距離3000公尺的尺度下，台北速食店分布呈現顯著「cluster」分布

從ci圖可以看到黑線(觀測值)整體高於envelope，印證了在3000公尺尺度下，台北速食店分布顯示顯著的cluster分布

3. 利用Bivariate F function分析方法，比較公立vs.私立學校，到鄰近速食店的空間特徵，並說明哪一種類型學校的附近，會有較多的速食店？（顯著水準=0.05）

pub = school[school$Type == "公立", ]
pri = school[school$Type == "私立", ]

coord_pub = st_coordinates(pub)
coord_pri = st_coordinates(pri)
coord_ff = st_coordinates(food)

windows = as.owin(st_union(vill))

pub_pp = as.ppp(coord_pub, windows)
pri_pp = as.ppp(coord_pri, windows)
ff_pp = as.ppp(coord_ff, windows)

NND_pub = nncross(pub_pp, ff_pp)
Fd_pub = ecdf(NND_pub$dist)

NND_pri = nncross(pri_pp, ff_pp)
Fd_pri = ecdf(NND_pri$dist)

n_pub = pub_pp$n
n_pri = pri_pp$n

F_sim_pub = matrix(nrow=99, ncol=n_pub)
for(i in 1:99){
  RND = rpoint(n_pub, win = windows)
  F_sim_pub[i, ] = nncross(RND, ff_pp)$dist
}
F_sim_s_pub = t(apply(F_sim_pub,1,sort))
F_sim_sort_pub = apply(F_sim_s_pub, 2, sort)

F_sim_pri = matrix(nrow=99, ncol=n_pri)
for(i in 1:99){
  RND = rpoint(n_pri, win = windows)
  F_sim_pri[i, ] = nncross(RND, ff_pp)$dist
}
F_sim_s_pri = t(apply(F_sim_pri,1,sort))
F_sim_sort_pri = apply(F_sim_s_pri, 2, sort)

plot(Fd_pub, col="blue", cex=0, main="F function: 公立學校 → 速食店", xlab="距離", ylab="F(d)", xlim=c(0,4000))
for(i in 1:99){ 
  lines(ecdf(F_sim_pub[i,]), col="grey", verticals=T, cex=0)
}

lines(ecdf(F_sim_sort_pub[5,]), col="black", lty=2, lwd=1.5)  
lines(ecdf(F_sim_sort_pub[95,]), col="black", lty=2, lwd=1.5) 

lines(Fd_pub, col="blue", lwd=3)

plot(Fd_pri, col="blue", cex=0, main="F function: 私立學校 → 速食店", xlab="距離", ylab="F(d)", xlim=c(0,4000))
for(i in 1:99){ 
  lines(ecdf(F_sim_pri[i,]), col="grey", verticals=T, cex=0)
}

lines(ecdf(F_sim_sort_pri[5,]), col="black", lty=2, lwd=1.5)  
lines(ecdf(F_sim_sort_pri[95,]), col="black", lty=2, lwd=1.5) 

lines(Fd_pri, col="blue", lwd=3)

plot(Fd_pub, col="blue", main="公立 vs. 私立學校 → 速食店", xlab="距離", ylab="F(d)", xlim=c(0,4000))
lines(Fd_pub, col="blue", lwd=3)

lines(Fd_pri, col="red", lwd=3)
legend("bottomright", legend = c("公立學校", "私立學校"), col = c("blue", "red"), lwd = 3)

從前兩張圖可以看到，在顯著水準0.05的情況下，台北市公私立學校附近的速食店皆呈現cluster分布，也就是兩種類型學校附近皆圍繞著速食店，而從第三張圖可以看到公立學校f-function的上升速度較快(較cluster)，換句話說，也就是公立學校附近的速食店較多