41023224L李姿瑾_FinalReport

school_sf <- st_read("./DATA/SCHOOL/SCHOOL.shp",options = "encoding=Big5" )
Vill_sf <- st_read("./DATA/Taipei_Vill/Taipei_Vill.shp",options = "encoding=Big5" )
FF_sf <- st_read("./DATA/Tpe_Fastfood/Tpe_Fastfood.shp",options = "encoding=Big5" )

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

問答題

1. 此圖中的 X 軸表示各村里速食店密度的標準化值，Y軸則表示鄰近村里的速食店密度之平均值，即空間滯後值，也經過標準化。由於圖中大多數點分布於第一象限（高密度鄰近高密度）與第三象限（低密度鄰近低密度），顯示該城市速食店密度具有明顯的空間群聚特徵，屬於正空間自相關的型態。

2. 該圖中的趨勢線斜率值為 0.982，代表整體速食店密度與鄰近區域速食店密度之間具有高度正向關聯，此斜率即為 Global Moran’s I 值。此數值越接近 1，表示該城市的速食店密度分布具有越強的正空間自相關，代表高密度區附近傾向出現其他高密度區，低密度區則鄰近其他低密度區，呈現出明顯的空間群聚現象。

3. 基於圖 1，我認為其 Local Moran’s I 統計量是大於 0，因為圖中多數點分布於第一象限（高值鄰近高值）與第三象限（低值鄰近低值），顯示該城市的速食店密度呈現正向空間自相關，Local Moran’s I 值會反映其群聚的程度，因此在這種情形下其統計量會大於 0。

實作題

1.

p >= 0.05，無法拒絕虛無假設，表示無顯著差異。

FF_vill_sf <- st_join(FF_sf, Vill_sf, join = st_within, left = FALSE)

a_ff <- FF_vill_sf %>% 
  filter(TOWN %in% c("文山區", "大安區", "中正區"))

# B區：信義、南港、松山
b_ff <- FF_vill_sf %>% 
  filter(TOWN %in% c("信義區", "南港區", "松山區"))

a_buffer <- st_buffer(a_ff, dist = 1000)
b_buffer <- st_buffer(b_ff, dist = 1000)

a_school_counts <- sapply(1:nrow(a_buffer), function(i) {
  sum(st_intersects(a_buffer[i,], school_sf, sparse = FALSE))
})

b_school_counts <- sapply(1:nrow(b_buffer), function(i) {
  sum(st_intersects(b_buffer[i,], school_sf, sparse = FALSE))
})

t_result <- t.test(a_school_counts, b_school_counts, var.equal = FALSE)

cat("統計檢定摘要\n")

統計檢定摘要

cat("虛無假設 H₀：A區與B區速食店涵蓋學校數平均值無差異\n")

虛無假設 H₀：A區與B區速食店涵蓋學校數平均值無差異

cat("對立假設 H₁：A區與B區速食店涵蓋學校數平均值有差異\n")

對立假設 H₁：A區與B區速食店涵蓋學校數平均值有差異

cat("A區平均學校數：", mean(a_school_counts), "\n")

A區平均學校數： 3.870968 

cat("B區平均學校數：", mean(b_school_counts), "\n")

B區平均學校數： 3.318182 

cat("t 值：", round(t_result$statistic, 3), "\n")

t 值： 1.598 

cat("p 值：", round(t_result$p.value, 4), "\n")

p 值： 0.1163 

if(t_result$p.value < 0.05){
  cat("結果：p < 0.05，拒絕虛無假設，表示 A區與B區涵蓋學校數量有顯著差異。\n")
} else {
  cat("結果：p >= 0.05，無法拒絕虛無假設，表示無顯著差異。\n")
}

結果：p >= 0.05，無法拒絕虛無假設，表示無顯著差異。

2.

圖中為以速食店數量為變數所繪製的 Moran’s I Correlogram，X 軸代表鄰近距離的標準化單位（distance classes），Y 軸為對應距離下的 Moran’s I 自相關指標值。從圖中可以觀察到，前三個距離階層的 Moran’s I 值為正，且達統計顯著水準（紅點標記），顯示速食店於鄰近空間內具有明顯的群聚現象。而隨著距離增加，Moran’s I 呈現持續上升但仍偏低，表示速食店的空間群聚效應較為局部，整體空間結構趨向微弱正相關。基於此圖，可推估速食店的空間影響範圍約在前三個距離階層內（約 1000–1500 公尺），此範圍內的空間單元較可能因鄰近速食店影響而呈現分布集聚現象。

grid <- st_make_grid(Vill_sf, cellsize = 500, square = TRUE)
grid_sf <- st_sf(grid_id = 1:length(grid), geometry = grid)

ff_count <- st_join(grid_sf, FF_sf, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(count = n()) %>%
  st_as_sf()
View(school_sf)
View(school_sf)

nb <- poly2nb(ff_count, queen = TRUE)
lw <- nb2listw(nb, style = "W")

correlog <- sp.correlogram(nb, ff_count$count, order = 10,
                           method = "I", style = "W", zero.policy = TRUE)

correlog_mat <- matrix(unlist(correlog$res), ncol = 5, byrow = TRUE)
correlog_df <- as.data.frame(correlog_mat)
colnames(correlog_df) <- c("lag", "I", "expected", "sd", "p.value")

plot(correlog_df$lag, correlog_df$I, type = "b",
     xlab = "distance classes", ylab = "Moran I statistic",
     ylim = c(min(correlog_df$I) - 0.05, max(correlog_df$I) + 0.05))

abline(h = 0, lty = 2)

sig_lags <- correlog_df$p.value < 0.05
points(correlog_df$lag[sig_lags], correlog_df$I[sig_lags],
       col = "red", pch = 19)

3.

grid_sf$count_ff <- ff_count$count
grid_sf$count_school <- st_join(grid_sf, school_sf, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(count = n()) %>%
  pull(count)

grid_sf$count_school[is.na(grid_sf$count_school)] <- 0  # 替代 NA 為 0

nb <- poly2nb(grid_sf, queen = TRUE)
lw <- nb2listw(nb, style = "W")

gi_ff <- localG(grid_sf$count_ff, lw)
gi_school <- localG(grid_sf$count_school, lw)

p_ff <- 2 * pnorm(-abs(gi_ff))
p_school <- 2 * pnorm(-abs(gi_school))

p_ff_fdr <- p.adjust(p_ff, method = "fdr")
p_school_fdr <- p.adjust(p_school, method = "fdr")

grid_sf$ff_hotspot <- ifelse(p_ff_fdr < 0.05, gi_ff, 0)
grid_sf$school_hotspot <- ifelse(p_school_fdr < 0.05, gi_school, 0)

tm_shape(grid_sf) +
  tm_fill("ff_hotspot", palette = "-RdBu", title = "Fast food Gi*",
          style = "cont", midpoint = NA) +
  tm_borders() +
  tm_layout(main.title = "速食店熱區分析", legend.outside = TRUE)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_fill()`: instead of `style = "cont"`, use fill.scale = `tm_scale_continuous()`.
ℹ Migrate the argument(s) 'midpoint', 'palette' (rename to 'values') to 'tm_scale_continuous(<HERE>)'
For small multiples, specify a 'tm_scale_' for each multiple, and put them in a list: 'fill'.scale =
list(<scale1>, <scale2>, ...)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the visual variable `fill` namely
'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run `cols4all::c4a_gui()` to explore
them. The old palette name "-RdBu" is named "rd_bu" (in long format "brewer.rd_bu")

tm_shape(grid_sf) +
  tm_fill("school_hotspot", palette = "-RdBu", title = "School Gi*",
          style = "cont", midpoint = NA) +
  tm_borders() +
  tm_layout(main.title = "學校熱區分析", legend.outside = TRUE)

── tmap v3 code detected ─────────────────────────────────────────────────────────────────────────────────────
[v3->v4] `tm_fill()`: instead of `style = "cont"`, use fill.scale = `tm_scale_continuous()`.
ℹ Migrate the argument(s) 'midpoint', 'palette' (rename to 'values') to 'tm_scale_continuous(<HERE>)'
For small multiples, specify a 'tm_scale_' for each multiple, and put them in a list: 'fill'.scale =
list(<scale1>, <scale2>, ...)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the visual variable `fill` namely
'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run `cols4all::c4a_gui()` to explore
them. The old palette name "-RdBu" is named "rd_bu" (in long format "brewer.rd_bu")

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

實作題

1.

此圖比較大安區、文山區與信義區學校至最近速食店的距離累積分布情形，顯示信義區的學校可及性最高，約八成學校在 500 公尺內就能抵達速食店；大安區次之，約六成學校在同樣距離內；而文山區則可及性最低，多數學校距離速食店超過 500 公尺。整體而言，信義區學校暴露於速食環境的程度最高，文山區則相對較低。

SH_vill_sf <- st_join(school_sf, Vill_sf, join = st_within, left = FALSE)

school_da <- SH_vill_sf %>% filter(TOWN == "大安區")
school_ws <- SH_vill_sf %>% filter(TOWN == "文山區")
school_xy <- SH_vill_sf %>% filter(TOWN == "信義區")

min_dist <- function(school_pts, fastfood_pts) {
  dists <- st_distance(school_pts, fastfood_pts)
  apply(dists, 1, min) %>% as.numeric()
}

dist_da <- min_dist(school_da, FF_sf)
dist_ws <- min_dist(school_ws, FF_sf)
dist_xy <- min_dist(school_xy, FF_sf)

d_range <- seq(0, 2000, by = 10)

get_Fd <- function(dist_vec, d_range) {
  sapply(d_range, function(d) mean(dist_vec <= d))
}

Fd_da <- get_Fd(dist_da, d_range)
Fd_ws <- get_Fd(dist_ws, d_range)
Fd_xy <- get_Fd(dist_xy, d_range)

plot(d_range, Fd_da, type = "l", col = "blue", lwd = 3,
     xlab = "Distance (m)", ylab = "Cumulative Proportion",
     main = "F(d): 學校到速食店距離累積分布", ylim = c(0, 1))

lines(d_range, Fd_ws, col = "green", lwd = 3)

lines(d_range, Fd_xy, col = "red", lwd = 3)

text(2000, Fd_da[length(Fd_da)], "大安區", col = "blue", pos = 4)

text(2000, Fd_ws[length(Fd_ws)], "文山區", col = "green", pos = 4)
text(2000, Fd_xy[length(Fd_xy)], "信義區", col = "red", pos = 4)

legend("bottomright", legend = c("大安區", "文山區", "信義區"),
       col = c("blue", "green", "red"), lty = 1, lwd = 3)

2.

根據 Ripley’s L 函數分析結果，速食店在 3000 公尺距離下的L值為4374.25，明顯高於 Monte Carlo 模擬所得之 95% 信賴區間（上限為 3204.45），顯示其空間分布存在顯著的群聚現象。並顯示速食店在台北市區內傾向集中開立於特定地區，而非隨機分布。

ff_sp <- st_transform(FF_sf, 3826)
win <- as.owin(st_bbox(ff_sp))  
coords <- st_coordinates(ff_sp)
pp <- ppp(x = coords[,1], y = coords[,2], window = win)

警告： data contain duplicated points

Kest_result <- Kest(pp, r = seq(0, 3000, by = 100), correction = "Ripley")

Lest_result <- Lest(pp, r = seq(0, 3000, by = 100), correction = "Ripley")

envelope_result <- envelope(pp, fun = Lest, nsim = 99, rank = 1,
                            correction = "Ripley", r = seq(0, 3000, by = 100))

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.

plot(envelope_result, main = "L(r) with 95% Confidence Interval",
     legend = FALSE)
abline(a = 0, b = 1, lty = 2, col = "gray") 

K3000 <- Kest_result$iso[which(Kest_result$r == 3000)]
L3000 <- Lest_result$iso[which(Lest_result$r == 3000)]
L3000_lo <- envelope_result$lo[which(envelope_result$r == 3000)]
L3000_hi <- envelope_result$hi[which(envelope_result$r == 3000)]

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

K(3000)： 60111493 

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

L(3000)： 4374.252 

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

L(3000) 95% 信賴區間： 2810.895 ~ 3202.906 

3.

Bivariate F function 分析公立與私立學校到速食店的空間可及性，並進一步結合 CSR（完全隨機分布）模擬，以顯著水準=0.05 檢定兩者是否有顯著差異。從第一張圖顯示，私立學校的 F(d) 累積曲線明顯快於公立學校，表示其在較短距離內即有較多速食店可達，可及性較高。為檢定此差異是否具統計意義，分別對私立與公立學校進行 CSR 模擬（第二與第三張圖）。

結果顯示，私立學校的 F(d) 曲線明顯高於 CSR 模擬的 95% 信賴區間上限（圖二），表示速食店在空間上有顯著集中於私立學校周邊的現象（p < 0.05）。相對地，公立學校的 F(d) 曲線多數落在 CSR 模擬區間內（圖三），表示其鄰近速食店的分布趨近於隨機，無顯著群聚特性。綜上所述，私立學校的周邊擁有顯著較多的速食店，可及性高於公立學校，且達統計顯著水準的0.05。

school_public <- SH_vill_sf %>% filter(Type == "公立")
school_private <- SH_vill_sf %>% filter(Type == "私立")

win <- as.owin(st_bbox(school_sf))
coords_public <- st_coordinates(st_transform(school_public, 3826))
coords_private <- st_coordinates(st_transform(school_private, 3826))
coords_fastfood <- st_coordinates(st_transform(FF_sf, 3826))

pp_public <- ppp(coords_public[,1], coords_public[,2], window = win)
pp_private <- ppp(coords_private[,1], coords_private[,2], window = win)
pp_fastfood <- ppp(coords_fastfood[,1], coords_fastfood[,2], window = win)

警告： data contain duplicated points

F_public <- Fest(pp_public, pp_fastfood, correction = "rs", r = seq(0, 2000, by = 25))
F_private <- Fest(pp_private, pp_fastfood, correction = "rs", r = seq(0, 2000, by = 25))

plot(F_public$r, F_public$rs, type = "l", col = "blue", lwd = 2,
     xlab = "距離 d (公尺)", ylab = "F(d)",
     main = "公立 vs 私立學校到速食店 F(d)", ylim = c(0, 1))  # 加這行！

lines(F_private$r, F_private$rs, col = "red", lwd = 2)

legend("bottomright", legend = c("公立學校", "私立學校"),
       col = c("blue", "red"), lty = 1, lwd = 2)

F_env_private <- envelope(pp_private, fun = Fest,
                          simulate = expression(runifpoint(pp_fastfood$n, win = win)),
                          nsim = 99, rank = 1,
                          correction = "rs", r = seq(0, 2000, by = 25))

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.

F_env_public <- envelope(pp_public, fun = Fest,
                         simulate = expression(runifpoint(pp_fastfood$n, win = win)),
                         nsim = 99, rank = 1,
                         correction = "rs", r = seq(0, 2000, by = 25))

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(F_env_private, main = "私立學校到速食店 F(d) 與 CSR 模擬 95% 信賴區間",
     xlab = "距離 d (公尺)", ylab = "F(d)", legend = FALSE)

lines(F_private$r, F_private$rs, col = "red", lwd = 2)

legend("bottomright", legend = c("私立學校", "CSR模擬區間"),
       col = c("red", "gray"), lty = c(1,1,NA), lwd = 2)

plot(F_env_public, main = "公立學校到速食店 F(d) 與 CSR 模擬 95% 信賴區間",
     xlab = "距離 d (公尺)", ylab = "F(d)", legend = FALSE)

lines(F_public$r, F_public$rs, col = "blue", lwd = 2)

legend("bottomright", legend = c("公立學校", "CSR模擬區間"),
       col = c("blue", "gray"), lty = c(1,1,NA), lwd = 2)

Part3

將 Vill_sf 合併為「臺北市整體邊界」

Taipei_city <- Vill_sf %>%
  summarise(geometry = st_union(geometry)) %>%
  st_make_valid() 
Taipei_city$COUNTY <- "臺北市"

school_sf增加COUNTY、TOWN、VILLAGE欄位

school_data_joined <- st_join(school_data, Vill_sf, join = st_within)
school_sf2 <- school_data_joined %>%
  select(COUNTY,TOWN,VILLAGE,Type,Name,SID)

將Vill_sf從村里層級，變成行政區層級

TOWN_sf <- Vill_sf %>%
  group_by(TOWN) %>%
  summarise(geometry = st_union(geometry)) %>%
  ungroup() %>%
  st_make_valid()  
TOWN_sf$COUNTY <- "臺北市"
TOWN_line <- TOWN_sf %>%
  st_cast("MULTILINESTRING")

STEP 1：計算每所學校的周邊速食店數量

school_buff <- st_buffer(school_sf2, dist = 500)

# 計算每個學校 buffer 內有幾間速食店
school_ff <- st_join(school_buff, FF_sf, join = st_intersects) %>%
  group_by(Name) %>%
  summarise(ff_count = n()) %>%
  st_drop_geometry()

# 合併筆數進原始學校點
school_data <- school_sf2 %>%
  left_join(school_ff, by = "Name") %>%
  mutate(ff_count = replace_na(ff_count, 0))

STEP 2：建立鄰接矩陣（以學校為點位k-nearest neighbor）

coords <- st_coordinates(st_centroid(school_data))

nb <- knn2nb(knearneigh(coords, k = 5))
lw <- nb2listw(nb, style = "W")

STEP 3：計算 Local G* 值（z-score）

gi_star <- localG(school_data$ff_count, lw)

school_data$GiZ <- as.numeric(gi_star)

# 確保 GiZ_class 分類完成
school_data <- school_data %>%
  mutate(GiZ_class = case_when(
    GiZ >= 2.58  ~ "極顯著熱點",  # p < 0.01
    GiZ >= 1.96  ~ "顯著熱點",    # p < 0.05
    GiZ <= -2.58 ~ "極顯著冷點",  # p < 0.01
    GiZ <= -1.96 ~ "顯著冷點",    # p < 0.05
    TRUE         ~ "非顯著"
  ))

# 畫出底圖 + 熱點分類地圖
ggplot() +
  geom_sf(data = Vill_sf, fill = "grey95", color = "grey80") +    # 行政區底圖
  

  geom_sf(data = school_data, aes(fill = GiZ_class), shape = 21, color = "black") +
  geom_sf(data = TOWN_line, color = "black", size = 0.5) +   # 加上邊界線
  theme_minimal()+
  
  geom_sf(data = school_data, aes(fill = GiZ_class), color = "black", size = 2, shape = 21) +
  scale_fill_manual(
    values = c(
      "極顯著熱點" = "#800026",
      "顯著熱點"   = "#FC4E2A",
      "非顯著"     = "white",
      "顯著冷點"   = "#2C7FB8",
      "極顯著冷點" = "#084081"
    ),
    name = "Local G* 分類"
  ) +
  labs(
    title = "臺北市國小周邊速食店 Local G* 熱點分析",
    caption = "註：分類依照 GiZ z-score 門檻（±1.96、±2.58）"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 14),
    legend.position = "right"
  )

NA
NA

ggplot(data.frame(GiZ = GiZ_vector), aes(x = GiZ)) +
  geom_histogram(binwidth = 0.3, fill = "lightblue", color = "black") +
  geom_vline(xintercept = c(-2.58, -1.96, 1.96, 2.58), 
             color = "darkgray", linetype = "dotted") +
  scale_x_continuous(
    breaks = c(-2.58, -1.96, 0, 1.96, 2.58),
    labels = c("-2.58", "-1.96", "0", "1.96", "2.58")
  ) +
  labs(title = "臺北市國小 Local G* Z-score 分布圖",
       x = "GiZ 值", y = "學校數量") +
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5)  
  )