Analisis Kasus Tuberkulosis di Kabupaten/Kota
Provinsi Jawa Timur Tahun 2024

Ditujukan Guna Memenuhi UAS Mata Kuliah Spasial



Disusun oleh:
Khalisha Syahla Putri Marindra (140610230008)


Dosen Pengampu:
I Gede Nyoman Mindra Jaya, Ph.D


PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

Abstrak

Tuberkulosis (TBC) masih menjadi masalah kesehatan masyarakat yang signifikan dengan distribusi kasus yang tidak merata secara geografis. Penelitian ini bertujuan untuk menganalisis pola spasial dan faktor-faktor yang memengaruhi insidensi TBC di tingkat kabupaten/kota menggunakan pendekatan analisis spasial dan pemodelan regresi berbasis lokasi. Data yang digunakan meliputi jumlah kasus TBC, insidensi rate (IR) per 100.000 penduduk, serta beberapa variabel penjelas pada 38 wilayah observasi.

Analisis spasial menunjukkan bahwa distribusi IR TBC tidak bersifat acak dan membentuk pola pengelompokan spasial. Hasil Standardized Morbidity Ratio (SMR) dan Local Indicator of Spatial Autocorrelation (LISA) mengidentifikasi adanya wilayah berisiko tinggi (hotspot), terutama di Mojokerto dan sekitarnya. Pemodelan Geographically Weighted Regression (GWR) menunjukkan bahwa pengaruh variabel penjelas terhadap IR TBC bervariasi antar wilayah, dengan tingkat kecocokan model yang berbeda secara spasial. Selain itu, hasil interpolasi menggunakan Inverse Distance Weighting (IDW) dan Ordinary Kriging menggambarkan permukaan risiko TBC yang kontinu dan konsisten dengan temuan klaster spasial.

Secara keseluruhan, hasil penelitian ini menegaskan adanya heterogenitas dan ketergantungan spasial dalam distribusi TBC, sehingga analisis berbasis lokasi menjadi penting untuk memahami variasi risiko antar wilayah.

Kata kunci: Tuberkulosis, Spasial, Jawa Timur

BAB I

LATAR BELAKANG

Tuberkulosis (TBC) merupakan salah satu penyakit menular yang masih menjadi permasalahan kesehatan masyarakat di Indonesia hingga saat ini. Penyakit ini disebabkan oleh bakteri Mycobacterium tuberculosis yang menular melalui droplet udara dan terutama menyerang paru-paru [1]. Indonesia masih termasuk dalam negara dengan beban TBC tertinggi di dunia, dengan jumlah kasus yang besar dan penularan yang masih aktif di berbagai wilayah [2]. Kondisi ini menuntut upaya pengendalian yang tidak hanya bersifat klinis, tetapi juga berbasis wilayah melalui pendekatan analisis spasial untuk memahami pola penyebaran penyakit secara geografis [3].

Provinsi Jawa Timur merupakan salah satu provinsi dengan jumlah kasus Tuberkulosis yang tinggi di Indonesia. Berdasarkan data Dinas Kesehatan Provinsi Jawa Timur, hingga tahun 2024 tercatat lebih dari 70 ribu kasus TBC yang tersebar di berbagai kabupaten/kota [4]. Namun, jumlah kasus tersebut tidak terdistribusi secara merata antarwilayah. Beberapa daerah, khususnya wilayah perkotaan dan kawasan dengan kepadatan penduduk tinggi seperti Kota Surabaya dan kabupaten penyangga sekitarnya, menunjukkan jumlah kasus yang jauh lebih tinggi dibandingkan daerah lain [5]. Variasi ini mengindikasikan adanya perbedaan risiko penularan TBC yang berkaitan erat dengan karakteristik spasial suatu wilayah.

Analisis spasial menjadi penting karena mampu menggambarkan distribusi geografis kasus TBC, mengidentifikasi wilayah dengan risiko tinggi (hotspot), serta mendeteksi adanya autokorelasi spasial antarwilayah [6]. Melalui pendekatan ini, perbedaan jumlah kasus antar kabupaten/kota tidak hanya dipahami sebagai variasi angka semata, tetapi juga sebagai hasil dari pengaruh lingkungan dan keterkaitan spasial antarwilayah. Pendekatan spasial memungkinkan pemanfaatan peta tematik dan metode statistik spasial untuk mengungkap pola persebaran penyakit yang tidak dapat dijelaskan secara optimal oleh analisis non-spasial [3,6].

Penelitian ini berfokus pada analisis spasial kasus Tuberkulosis di tingkat kabupaten/kota di Provinsi Jawa Timur tahun 2024. Penelitian ini bertujuan untuk menggambarkan pola persebaran geografis kasus TBC serta menganalisis keterkaitan spasial antarwilayah administratif. Hasil analisis diharapkan dapat memberikan pemahaman yang lebih mendalam mengenai pola distribusi TBC di Jawa Timur dan menjadi dasar dalam perencanaan intervensi kesehatan masyarakat yang lebih efektif, terarah, dan berbasis wilayah.

BAB II

TINJAUAN PUSTAKA

a. Spatial Dependence

Spatial Dependence atau ketergantungan spasial merupakan kondisi di mana nilai suatu variabel pada suatu lokasi dipengaruhi oleh nilai variabel yang sama di lokasi lain yang berdekatan [7]. Artinya, data yang berdekatan secara geografis cenderung memiliki karakteristik yang mirip, sehingga observasi di satu wilayah tidak bersifat independen [8]. Dalam analisis spasial, spatial dependence menunjukkan adanya pola keterkaitan antarwilayah yang perlu diperhitungkan agar hasil analisis tidak bias [9]. Jika wilayah dengan jumlah kasus tinggi berdekatan dengan wilayah lain yang juga memiliki jumlah kasus tinggi, menunjukkan adanya spatial dependence positif [10]. Sebaliknya, wilayah dengan nilai tinggi dikelilingi oleh wilayah dengan nilai rendah, maka menunjukkan adanya spatial dependence negatif [7]. Ketergantungan spasial dapat diukur menggunakan metode seperti Moran’s I, Geary’s C, dan LISA (Local Indicators of Spatial Association) untuk mendeteksi pola autokorelasi spasial baik secara global maupun lokal. Spatial dependence sangat penting dalam analisis statistik spasial [8], karena dapat membantu peneliti mengidentifikasi pola penyebaran dan klasterisasi fenomena yang terjadi di ruang geografis [9].

b. Autokorelasi Spasial

Autokorelasi spasial adalah konsep yang menggambarkan hubungan atau keterkaitan antara nilai suatu variabel di suatu lokasi dengan nilai variabel yang sama di lokasi lain yang berdekatan secara geografis [10]. Dalam analisis spasial, autokorelasi spasial digunakan untuk mengidentifikasi apakah suatu fenomena memiliki pola yang teratur, acak, atau menyebar [11]. Jika wilayah dengan nilai tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki nilai tinggi, maka terjadi autokorelasi spasial positif [12]. Sebaliknya, jika wilayah dengan nilai tinggi dikelilingi oleh wilayah dengan nilai rendah, maka terjadi autokorelasi spasial negatif [10]. Untuk mengukurnya autokorelasi spasial digunakan Moran’s I, Geary’s C, dan Getis-Ord Gi*, yang dapat digunakan untuk mendeteksi adanya pola global maupun lokal pada data spasial [13].

c. Model Spatial Econometrics

Model Spatial Econometrics merupakan pendekatan analisis yang digunakan ketika data menunjukkan adanya keterkaitan antarwilayah atau efek spasial yang signifikan [14]. Tiga model yang umum digunakan dalam analisis ini adalah Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), dan Spatial Durbin Model (SDM) [15]. Model SAR digunakan ketika terdapat pengaruh langsung antarwilayah, yaitu nilai variabel dependen pada suatu lokasi dipengaruhi oleh nilai variabel dependen pada lokasi lain yang berdekatan [14]. Dengan kata lain, model SAR memperhitungkan adanya efek spasial dalam variabel terikat melalui matriks pembobot spasial (spatial weight matrix) yang merepresentasikan hubungan kedekatan antarwilayah [15].

Sementara itu, SEM digunakan ketika ketergantungan spasial tidak muncul secara langsung pada variabel dependen, melainkan pada komponen galat atau error model [16]. Model ini mengasumsikan bahwa kesalahan pada suatu wilayah dapat berkorelasi dengan kesalahan pada wilayah lain yang berdekatan akibat adanya faktor-faktor tak teramati yang bersifat spasial. Dengan demikian, SEM berperan dalam mengatasi masalah bias dan inefisiensi estimasi yang sering muncul pada model regresi klasik Ordinary Least Squares (OLS) ketika asumsi independensi error tidak terpenuhi [15].

Adapun Spatial Durbin Model (SDM) merupakan model yang lebih umum dan fleksibel karena mengakomodasi ketergantungan spasial baik pada variabel dependen maupun pada variabel independen [16]. Dalam model ini, nilai variabel dependen suatu wilayah tidak hanya dipengaruhi oleh nilai variabel dependen wilayah lain yang berdekatan, tetapi juga oleh nilai variabel independen dari wilayah tetangga melalui lag spasial variabel penjelas. Dengan karakteristik tersebut, SDM mampu menangkap efek limpahan (spillover effects) antarwilayah secara lebih komprehensif, baik efek langsung maupun tidak langsung. Selain itu, SDM sering dianggap sebagai model yang mencakup SAR dan SEM sebagai kasus khusus, sehingga dapat digunakan sebagai pendekatan awal untuk mengidentifikasi bentuk ketergantungan spasial yang paling sesuai dalam suatu analisis.

Ketiga model tersebut memiliki peran penting dalam analisis data spasial karena mampu menangkap hubungan geografis yang tidak dapat dijelaskan oleh model regresi klasik. Penggunaan model Spatial Econometrics memungkinkan pemahaman yang lebih akurat terhadap fenomena yang tersebar secara spasial, terutama dalam konteks analisis wilayah seperti distribusi penyakit, ekonomi regional, dan perencanaan kebijakan berbasis spasial [16].

BAB III

METODOLOGI PENELITIAN

1. Jenis dan Sumber Data

Tuberkulosis (TBC) adalah penyakit menular yang disebabkan oleh bakteri Mycobacterium tuberculosis yang umumnya menyerang paru-paru.Penelitian ini menggunakan data sekunder yang bersifat cross-section pada tahun 2024. Data dianalisis pada tingkat kabupaten/kota di Provinsi Jawa Timur.. Data diperoleh dari publikasi resmi Badan Pusat Statistik (BPS). Unit analisis dalam penelitian ini adalah 38 kabupaten/kota di Provinsi Jawa Timur, di mana setiap unit direpresentasikan sebagai satu entitas spasial dengan batas wilayah administratif yang jelas.

2. Variabel Penelitian

Selain variabel dependen berupa jumlah kasus TBC tahun 2024 (Y), penelitian ini juga melibatkan enam variabel independen, yaitu:

X1: Kepadatan Penduduk,

X2: Jumlah Penduduk Miskin,

X3: Konsumsi Rokok Tanpa Filter,

Pemilihan Provinsi Jawa Timur sebagai wilayah penelitian didasarkan pada tingginya variasi jumlah kasus TBC antarwilayah, sehingga memungkinkan untuk mengidentifikasi adanya pola ketergantungan spasial yang signifikan antar kabupaten/kota. Penelitian ini menggunakan pendekatan spatial econometrics untuk mengidentifikasi adanya efek spasial pada persebaran TBC. Tahapan analisis meliputi:

  1. Eksplorasi Data Spasial, yaitu menyajikan peta persebaran kasus TBC untuk mengenali pola awal.

3. Uji Autokorelasi Spasial

3.1 Moran’s I (Global)

Uji Moran’s I digunakan untuk mendeteksi autokorelasi spasial global:

\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \frac{\sum_i \sum_j w_{ij}(x_i - \bar{x})(x_j - \bar{x})} {\sum_i (x_i - \bar{x})^2} \]

3.2 LISA (Local Indicators of Spatial Association)

LISA digunakan untuk mengidentifikasi klaster spasial lokal (High-High, Low-Low, High-Low, Low-High).

4. Model Regresi Spasial

4.1 Ordinary Least Squares (OLS)

\[ y = X\beta + \varepsilon \]

4.2 Spatial Autoregressive Model (SAR)

\[ y = \rho Wy + X\beta + \varepsilon \]

4.3 Spatial Error Model (SEM)

\[ y = X\beta + u,\quad u = \lambda Wu + \varepsilon \]

4.4 Spatial Durbin Model (SDM)

\[ y = \rho Wy + X\beta + WX\theta + \varepsilon \]

5. Geographically Weighted Regression (GWR)

Model GWR digunakan untuk menangkap variasi hubungan lokal antarwilayah:

\[ y_i = \beta_0(u_i,v_i) + \sum_k \beta_k(u_i,v_i)x_{ik} + \varepsilon_i \]

Parameter diestimasi secara lokal dengan bandwidth adaptif.

6. Interpolasi Spasial

6.1 Inverse Distance Weighting (IDW)

\[ \hat{Z}(s_0) = \sum_{i=1}^n w_i Z(s_i) \]

\[ w_i = \frac{1}{d_{i0}^p} \]

9.2 Ordinary Kriging

\[ Z(s) = \mu + \varepsilon(s) \]

Ketidakpastian prediksi diukur melalui varians kriging, dan validasi dilakukan menggunakan RMSE cross-validation.

7. Perangkat Lunak

Seluruh analisis dilakukan menggunakan RStudio

8. Alur Penelitian

Secara umum, alur penelitian dimulai dari pengumpulan dan pengolahan data, pembuatan peta deskriptif, pengujian autokorelasi spasial, estimasi model SEM, SAR, OLS, pemilihan model terbaik, Analaisis, GWR, dan Analasis Interpolasi hingga interpretasi hasil dalam bentuk peta tematik. Alur ini disusun agar analisis spasial dapat menggambarkan hubungan antarwilayah secara menyeluruh.

BAB IV

HASIL DAN PEMBAHASAN

4.1 Peta Deskriptif Persebaran Kasus Tuberkulosis di Jawa Timur Tahun 2024

library(sf)
library(readxl)
library(dplyr)
library(ggplot2)

clean_name <- function(x) tolower(gsub("\\s+|kabupaten|kota|,", "", x))

to_num_safe <- function(v) {
  if (is.numeric(v)) return(v)
  v <- as.character(v)
  v <- gsub("\\.", "", v)   # buang ribuan titik
  v <- gsub(",", ".", v)    # desimal koma -> titik
  suppressWarnings(as.numeric(v))
}

path_zip   <- "GADM PETA JATIM.zip"               # contoh: "data/GADM PETA JATIM.zip"
path_excel <- "Kasus Penyakit Menurut Kabupaten.xlsx"  # contoh: "data/Kasus....xlsx"
tahun_pilih <- 2024

tmpdir <- file.path(tempdir(), paste0("shp_", as.integer(Sys.time())))
dir.create(tmpdir, showWarnings = FALSE, recursive = TRUE)
utils::unzip(path_zip, exdir = tmpdir)

shp_files <- list.files(tmpdir, pattern = "\\.shp$", full.names = TRUE, recursive = TRUE, ignore.case = TRUE)
stopifnot(length(shp_files) > 0)

shp <- sf::st_read(shp_files[1], quiet = TRUE) |> sf::st_make_valid()

prov_col <- names(shp)[grepl("NAME_1|prov", names(shp), ignore.case = TRUE)]
if (length(prov_col) > 0) {
  shp <- shp |> dplyr::filter(grepl("Jawa Timur", .data[[prov_col[1]]], ignore.case = TRUE))
}

if (is.na(sf::st_crs(shp))) sf::st_crs(shp) <- 4326
stopifnot(nrow(shp) > 0)
stopifnot("NAME_2" %in% names(shp))

df <- readxl::read_excel("C:/Users/khalisha marindra/Downloads/Kasus Penyakit Menurut Kabupaten_Kota dan Jenis Penyakit di Provinsi Jawa Timur, 2024.xlsx")
names(df) <- trimws(gsub("\\s+", "_", names(df)))

x_idx <- grep("^X", names(df))
if (length(x_idx) > 0) {
  for (nm in names(df)[x_idx]) df[[nm]] <- to_num_safe(df[[nm]])
}
if ("Kasus_TBC" %in% names(df)) df$Kasus_TBC <- to_num_safe(df$Kasus_TBC)

stopifnot("Tahun" %in% names(df))
df <- df |> dplyr::filter(Tahun == tahun_pilih)
stopifnot(nrow(df) > 0)

shp$.key <- clean_name(shp$NAME_2)

if ("Kabupaten_Kota" %in% names(df)) {
  df$.key <- clean_name(df$Kabupaten_Kota)
} else if ("Kabupaten" %in% names(df)) {
  df$.key <- clean_name(df$Kabupaten)
} else if ("Kota" %in% names(df)) {
  df$.key <- clean_name(df$Kota)
} else {
  stop("Kolom nama wilayah tidak ditemukan di Excel (Kabupaten_Kota/Kabupaten/Kota).")
}

joined <- dplyr::left_join(shp, df, by = ".key") |> sf::st_as_sf()
stopifnot(nrow(joined) > 0)
stopifnot("Kasus_TBC" %in% names(joined))
centroids <- sf::st_point_on_surface(joined)

ggplot(joined) +
  geom_sf(aes(fill = Kasus_TBC), color = "grey60", linewidth = 0.2) +
  geom_sf_text(data = centroids, aes(label = NAME_2), size = 2.8, color = "black") +
  scale_fill_viridis_c(option = "C", direction = -1, na.value = "grey95") +
  theme_minimal() +
  labs(
    title = paste0("Peta Kasus TBC Jawa Timur (", tahun_pilih, ")"),
    fill  = "Kasus TBC"
  )

Peta kasus Tuberkulosis (TBC) di Provinsi Jawa Timur tahun 2024 menunjukkan bahwa jumlah kasus TBC tidak terdistribusi secara merata antar kabupaten/kota. Wilayah dengan jumlah kasus tertinggi cenderung terkonsentrasi di kawasan perkotaan, terutama Kota Surabaya, yang ditandai dengan warna paling gelap pada peta. Selain Surabaya, beberapa kabupaten dengan jumlah kasus relatif tinggi juga terlihat di wilayah Jember, Sidoarjo, dan Malang, yang merupakan daerah dengan jumlah penduduk besar dan aktivitas sosial ekonomi yang tinggi.

Sebaliknya, wilayah di bagian barat dan selatan Jawa Timur, seperti Pacitan, Trenggalek, dan beberapa kabupaten di wilayah Mataraman, menunjukkan jumlah kasus TBC yang relatif lebih rendah. Perbedaan ini mengindikasikan bahwa kepadatan penduduk, tingkat urbanisasi, serta intensitas interaksi sosial berperan penting dalam variasi jumlah kasus TBC antarwilayah.

Namun demikian, peta ini masih merepresentasikan jumlah kasus, sehingga wilayah dengan populasi besar cenderung memiliki jumlah kasus yang lebih tinggi.

4.2Uji Autokorelasi Spasial Global dan Lokal (Moran’s I dan LISA)

library(sf)
library(dplyr)
library(spdep)
library(ggplot2)

shp_path   <- "data_shp/JATIM_KABKOTA.shp"  # atau hasil unzip shapefile
excel_path <- "data_excel/Kasus_TBC.xlsx"  # file excel

# Baca shapefile & excel
shp <- sf::st_read("C:/Users/khalisha marindra/Downloads/GADM PETA JATIM/GADM_JATIM_KABKOTA.shp", quiet = TRUE) %>%
  sf::st_make_valid()

df <- readxl::read_excel("C:/Users/khalisha marindra/Downloads/Kasus Penyakit Menurut Kabupaten_Kota dan Jenis Penyakit di Provinsi Jawa Timur, 2024.xlsx")
names(df) <- trimws(gsub("\\s+", "_", names(df)))

clean_name <- function(x) tolower(gsub("\\s+|kabupaten|kota|,", "", x))

tahun_pakai <- 2024
df <- df %>% dplyr::filter(Tahun == tahun_pakai)

# Buat key join
shp$.key <- clean_name(shp$NAME_2)

if ("Kabupaten_Kota" %in% names(df)) {
  df$.key <- clean_name(df$Kabupaten_Kota)
} else if ("Kabupaten" %in% names(df)) {
  df$.key <- clean_name(df$Kabupaten)
} else if ("Kota" %in% names(df)) {
  df$.key <- clean_name(df$Kota)
} else {
  stop("Kolom nama wilayah tidak ditemukan di excel (Kabupaten_Kota/Kabupaten/Kota).")
}

# Join
sf_join <- dplyr::left_join(shp, df, by = ".key") %>% sf::st_as_sf()

stopifnot("Kasus_TBC" %in% names(sf_join))

ok <- !is.na(sf_join$Kasus_TBC)
sf_m <- sf_join[ok, ]

# Neighbors & listw
nb_full <- spdep::poly2nb(sf_join, queen = TRUE)
nb_sub  <- subset(nb_full, ok)
lw_sub  <- spdep::nb2listw(nb_sub, style = "W", zero.policy = TRUE)

# Variabel analisis
y <- sf_m$Kasus_TBC

moran_res <- spdep::moran.test(y, lw_sub, zero.policy = TRUE)

moran_I  <- unname(moran_res$estimate[["Moran I statistic"]])
moran_EI <- unname(moran_res$estimate[["Expectation"]])
moran_p  <- moran_res$p.value

cat("GLOBAL Moran's I\n")
## GLOBAL Moran's I
cat("Moran I   :", round(moran_I, 4), "\n")
## Moran I   : 0.3738
cat("E(I)      :", round(moran_EI, 4), "\n")
## E(I)      : -0.027
cat("p-value   :", signif(moran_p, 4), "\n\n")
## p-value   : 0.0002389
# =========================
# 4) LISA (Local Moran)
# =========================
alpha_lisa <- 0.05  # ubah kalau mau

local <- spdep::localmoran(y, lw_sub, zero.policy = TRUE)

# Z-score & spatial lag (untuk klasifikasi HH/LL/HL/LH)
z     <- as.numeric(scale(y))
lag_z <- spdep::lag.listw(lw_sub, z, zero.policy = TRUE)

sig <- local[, 5] < alpha_lisa

cluster <- rep("Not Sig", length(z))
cluster[z > 0 & lag_z > 0 & sig] <- "High-High"
cluster[z < 0 & lag_z < 0 & sig] <- "Low-Low"
cluster[z > 0 & lag_z < 0 & sig] <- "High-Low"
cluster[z < 0 & lag_z > 0 & sig] <- "Low-High"

sf_m$LISA <- factor(cluster, levels = c("High-High","Low-Low","High-Low","Low-High","Not Sig"))

# =========================
# 5) PETA LISA
# =========================
ggplot(sf_m) +
  geom_sf(aes(fill = LISA), color = "grey60", linewidth = 0.2) +
  theme_minimal() +
  labs(
    title = paste0(
      "Peta LISA Kasus TBC (", tahun_pakai, ") | Global Moran's I = ",
      round(moran_I, 3), " (p = ", signif(moran_p, 3), ")"
    ),
    fill = "LISA Cluster"
  )

Hasil analisis autokorelasi spasial menggunakan Global Moran’s I menunjukkan nilai Moran’s I sebesar 0,3738 dengan p-value = 0,000239. Nilai Moran’s I yang positif mengindikasikan adanya autokorelasi spasial positif yang kuat pada distribusi kasus Tuberkulosis antar kabupaten/kota di Provinsi Jawa Timur. Hal ini menunjukkan bahwa wilayah dengan jumlah kasus TBC tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki jumlah kasus tinggi, sehingga distribusi kasus tidak bersifat acak secara geografis.

Berdasarkan hasil Local Indicators of Spatial Association (LISA), teridentifikasi klaster High–High yang signifikan di wilayah Kota Surabaya dan kabupaten/kota sekitarnya. Klaster High–High menunjukkan bahwa wilayah dengan jumlah kasus TBC tinggi dikelilingi oleh wilayah dengan jumlah kasus tinggi pula, sehingga dapat dikategorikan sebagai hotspot penularan TBC di Jawa Timur. Keberadaan klaster ini mengindikasikan adanya proses penularan yang dipengaruhi oleh kedekatan geografis, kepadatan penduduk, serta interaksi sosial antarwilayah.

Sebagian besar wilayah lainnya berada dalam kategori Not Significant, yang menunjukkan tidak adanya autokorelasi spasial lokal yang signifikan pada wilayah tersebut. Kondisi ini mengindikasikan bahwa variasi jumlah kasus TBC di wilayah-wilayah tersebut relatif tidak dipengaruhi oleh wilayah tetangganya secara langsung.

Secara keseluruhan, hasil Moran’s I dan LISA menunjukkan bahwa penyebaran kasus TBC di Jawa Timur memiliki pola spasial yang terklaster, khususnya di kawasan perkotaan dan sekitarnya.

4.3 Perbandingan Kinerja Model Spasial

library(spatialreg)
library(broom)
library(DT)

data_model0 <- sf_m %>%
  mutate(
    Kasus_TBC = as.numeric(Kasus_TBC),
    X1 = as.numeric(X1),
    X2 = as.numeric(X2),
    X3 = as.numeric(X3),
    X4 = as.numeric(X4)
  )

needed <- c("Kasus_TBC","X1","X2","X3","X4")
ok <- complete.cases(st_drop_geometry(data_model0)[, needed, drop = FALSE])
sf_m <- data_model0[ok, ]
stopifnot(nrow(sf_m) >= 10)

df <- st_drop_geometry(sf_m)

nb_full <- spdep::poly2nb(sf_m, queen = TRUE)
lw_full <- spdep::nb2listw(nb_full, style = "W", zero.policy = TRUE)

f <- Kasus_TBC ~ X1 + X2 + X3 + X4

ols <- lm(f, data = df)

sar <- spatialreg::lagsarlm(
  f, data = df, listw = lw_full,
  zero.policy = TRUE
)

sem <- spatialreg::errorsarlm(
  f, data = df, listw = lw_full,
  zero.policy = TRUE
)

sdm <- spatialreg::lagsarlm(
  f, data = df, listw = lw_full,
  Durbin = TRUE, zero.policy = TRUE
)

safe_mse <- function(obs, pred) {
  ok2 <- !is.na(obs) & !is.na(pred)
  if (sum(ok2) == 0) return(NA_real_)
  mean((obs[ok2] - pred[ok2])^2)
}

safe_rmse <- function(obs, pred) sqrt(safe_mse(obs, pred))

safe_mape <- function(obs, pred) {
  ok2 <- !is.na(obs) & !is.na(pred) & obs != 0
  if (sum(ok2) == 0) return(NA_real_)
  mean(abs((obs[ok2] - pred[ok2]) / obs[ok2])) * 100
}

obs <- df$Kasus_TBC

pred_ols <- fitted(ols)
pred_sar <- fitted(sar)
pred_sem <- fitted(sem)
pred_sdm <- fitted(sdm)

hasil_compare <- tibble::tibble(
  n_model = nrow(df),
  Model = c("OLS", "SAR", "SEM", "SDM"),
  AIC  = c(AIC(ols), AIC(sar), AIC(sem), AIC(sdm)),
  MSE  = c(safe_mse(obs, pred_ols),
           safe_mse(obs, pred_sar),
           safe_mse(obs, pred_sem),
           safe_mse(obs, pred_sdm)),
  RMSE = c(safe_rmse(obs, pred_ols),
           safe_rmse(obs, pred_sar),
           safe_rmse(obs, pred_sem),
           safe_rmse(obs, pred_sdm)),
  MAPE = c(safe_mape(obs, pred_ols),
           safe_mape(obs, pred_sar),
           safe_mape(obs, pred_sem),
           safe_mape(obs, pred_sdm))
) %>% 
  arrange(AIC)

# -------------------------
# Tampilkan tabel perbandingan
# -------------------------
DT::datatable(
  hasil_compare,
  rownames = FALSE,
  options = list(pageLength = 10, scrollX = TRUE)
)
# -------------------------
# (Opsional) ringkasan model
# -------------------------
cat("\n\n=== Ringkasan SAR ===\n")
## 
## 
## === Ringkasan SAR ===
print(summary(sar))
## 
## Call:spatialreg::lagsarlm(formula = f, data = df, listw = lw_full, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1762.249  -217.218    24.923   255.478  1305.031 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -1.2042e+03  4.3872e+02 -2.7448  0.006056
## X1           1.9370e+00  2.0232e-01  9.5743 < 2.2e-16
## X2           2.2426e-01  5.6131e-02  3.9954  6.46e-05
## X3           7.2335e-01  2.3823e+00  0.3036  0.761410
## X4          -3.5224e+01  7.3548e+01 -0.4789  0.631991
## 
## Rho: 0.31777, LR test value: 7.5103, p-value: 0.0061347
## Asymptotic standard error: 0.087052
##     z-value: 3.6504, p-value: 0.00026184
## Wald statistic: 13.325, p-value: 0.00026184
## 
## Log likelihood: -291.2926 for lag model
## ML residual variance (sigma squared): 258650, (sigma: 508.58)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 596.59, (AIC for lm: 602.1)
## LM test for residual autocorrelation
## test value: 3.4941, p-value: 0.061588
cat("\n\n=== Ringkasan SEM ===\n")
## 
## 
## === Ringkasan SEM ===
print(summary(sem))
## 
## Call:spatialreg::errorsarlm(formula = f, data = df, listw = lw_full, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1706.207  -219.775   -24.713   338.803  1219.229 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -318.715253  374.684584 -0.8506    0.3950
## X1             2.055856    0.220939  9.3051 < 2.2e-16
## X2             0.237174    0.055068  4.3069 1.655e-05
## X3            -0.453731    2.676897 -0.1695    0.8654
## X4           -97.217593   78.876618 -1.2325    0.2178
## 
## Lambda: 0.39502, LR test value: 4.7108, p-value: 0.029973
## Asymptotic standard error: 0.16074
##     z-value: 2.4575, p-value: 0.013991
## Wald statistic: 6.0393, p-value: 0.013991
## 
## Log likelihood: -292.6923 for error model
## ML residual variance (sigma squared): 273530, (sigma: 523)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 599.38, (AIC for lm: 602.1)
cat("\n\n=== Ringkasan SDM ===\n")
## 
## 
## === Ringkasan SDM ===
print(summary(sdm))
## 
## Call:spatialreg::lagsarlm(formula = f, data = df, listw = lw_full, 
##     Durbin = TRUE, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1770.621  -196.030    88.407   253.932  1111.562 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -227.39004  723.89444 -0.3141 0.7534296
## X1             1.96745    0.21755  9.0437 < 2.2e-16
## X2             0.21329    0.05519  3.8647 0.0001112
## X3            -0.55735    2.74124 -0.2033 0.8388834
## X4            -8.55029   83.04470 -0.1030 0.9179947
## lag.X1        -0.23717    0.49606 -0.4781 0.6325797
## lag.X2        -0.16284    0.18429 -0.8836 0.3769018
## lag.X3        -1.78501    3.30518 -0.5401 0.5891525
## lag.X4      -146.50508  121.69605 -1.2039 0.2286435
## 
## Rho: 0.43623, LR test value: 6.3947, p-value: 0.011446
## Asymptotic standard error: 0.15017
##     z-value: 2.905, p-value: 0.0036723
## Wald statistic: 8.4391, p-value: 0.0036723
## 
## Log likelihood: -289.3943 for mixed model
## ML residual variance (sigma squared): 227310, (sigma: 476.77)
## Number of observations: 38 
## Number of parameters estimated: 11 
## AIC: 600.79, (AIC for lm: 605.18)
## LM test for residual autocorrelation
## test value: 1.3574, p-value: 0.24398

Berdasarkan hasil perbandingan kinerja model regresi, diperoleh bahwa model Spatial Autoregressive (SAR) memiliki nilai AIC paling kecil (596,59) dibandingkan model SEM, SDM, dan OLS. Nilai AIC yang lebih kecil menunjukkan bahwa model SAR memberikan keseimbangan terbaik antara tingkat kecocokan model dan kompleksitas parameter, sehingga secara keseluruhan model SAR dapat dianggap sebagai model paling cocok untuk menjelaskan variasi jumlah kasus Tuberkulosis di Jawa Timur.

Secara keseluruhan, hasil ini menunjukkan bahwa pendekatan spasial memberikan peningkatan kinerja model yang signifikan dibandingkan model lain. Model SAR direkomendasikan sebagai model terbaik berdasarkan kriteria AIC, sedangkan SDM unggul dari sisi ketepatan prediksi. Pemilihan model akhir dapat disesuaikan dengan tujuan analisis, yaitu apakah lebih menekankan pada interpretasi hubungan spasial (SAR) atau akurasi prediksi (SDM).

4.4 Analisis Geographically Weighted Regression

library(spgwr)
library(lmtest)
library(car)

pick_first_sf <- function(...) {
  cand <- list(...)
  for (obj in cand) if (inherits(obj, "sf")) return(obj)
  stop("Tidak ditemukan objek bertipe sf. Pastikan sudah ada sf_join/joined/sf_m/Jatim_sf di environment.")
}

sf_data <- pick_first_sf(
  if (exists("sf_m")) sf_m else NULL,
  if (exists("sf_join")) sf_join else NULL,
  if (exists("joined")) joined else NULL,
  if (exists("Jatim_sf")) Jatim_sf else NULL,
  if (exists("Jatim_Kabkota")) Jatim_Kabkota else NULL
)
find_col <- function(df, patterns) {
  nm <- names(df)
  idx <- which(sapply(patterns, \(p) any(grepl(p, nm, ignore.case = TRUE))))
  if (length(idx) == 0) return(NA_character_)
  # ambil match pertama dari pattern pertama yang ketemu
  for (p in patterns) {
    hit <- nm[grepl(p, nm, ignore.case = TRUE)]
    if (length(hit) > 0) return(hit[1])
  }
  NA_character_
}


col_kasus <- if ("Kasus_TBC" %in% names(sf_data)) "Kasus_TBC" else
  find_col(sf_data, c("^kasus", "tbc"))

col_x1 <- if ("X1" %in% names(sf_data)) "X1" else
  find_col(sf_data, c("^x1$", "penduduk", "populasi", "jumlah_penduduk"))

col_x2 <- if ("X2" %in% names(sf_data)) "X2" else
  find_col(sf_data, c("^x2$", "kepadatan", "density"))

col_x3 <- if ("X3" %in% names(sf_data)) "X3" else
  find_col(sf_data, c("^x3$", "miskin", "kemiskinan", "poverty"))

col_x4 <- if ("X4" %in% names(sf_data)) "X4" else
  find_col(sf_data, c("^x4$", "rokok", "cig", "smok"))

# validasi minimal untuk bikin IR_100k + GWR
if (is.na(col_x1)) stop("Kolom X1 (penduduk) tidak ditemukan. Sesuaikan mapping col_x1.")
if (is.na(col_x2)) stop("Kolom X2 tidak ditemukan. Sesuaikan mapping col_x2.")
if (is.na(col_x3)) stop("Kolom X3 tidak ditemukan. Sesuaikan mapping col_x3.")
if (is.na(col_x4)) stop("Kolom X4 tidak ditemukan. Sesuaikan mapping col_x4.")
if (is.na(col_kasus)) stop("Kolom Kasus_TBC tidak ditemukan. Sesuaikan mapping col_kasus.")

# rename agar seragam
sf_data <- sf_data %>%
  dplyr::rename(
    Kasus_TBC = all_of(col_kasus),
    X1 = all_of(col_x1),
    X2 = all_of(col_x2),
    X3 = all_of(col_x3),
    X4 = all_of(col_x4)
  )

# pastikan numeric + hitung IR_100k kalau belum ada
sf_data <- sf_data %>%
  dplyr::mutate(
    Kasus_TBC = suppressWarnings(as.numeric(Kasus_TBC)),
    X1 = suppressWarnings(as.numeric(X1)),
    X2 = suppressWarnings(as.numeric(X2)),
    X3 = suppressWarnings(as.numeric(X3)),
    X4 = suppressWarnings(as.numeric(X4)),
    Populasi = X1 * 1000,
    IR_100k = ifelse(!is.na(Populasi) & Populasi > 0 & !is.na(Kasus_TBC),
                     (Kasus_TBC / Populasi) * 100000,
                     NA_real_)
  )

needed <- c("IR_100k","X1","X2","X3","X4")
if (!all(needed %in% names(sf_data))) stop("Kolom target belum lengkap setelah standardisasi.")
# complete-case
sf_m2 <- sf_data %>%
  dplyr::filter(
    complete.cases(
      sf::st_drop_geometry(.) %>%
        dplyr::select(all_of(needed))
    )
  )


if (nrow(sf_m2) < 10) stop("Data valid < 10 baris (complete-case). Cek NA/variabel.")

pts_sf <- sf::st_point_on_surface(sf_m2)
pts_utm <- tryCatch(sf::st_transform(pts_sf, 32749), error=function(e) pts_sf)
coords  <- sf::st_coordinates(pts_utm)

dat <- pts_utm %>%
  sf::st_drop_geometry() %>%
  dplyr::transmute(
    y  = IR_100k,
    x1 = X1, x2 = X2, x3 = X3, x4 = X4,
    easting  = coords[,1],
    northing = coords[,2]
  )

sp::coordinates(dat) <- ~ easting + northing
sp::proj4string(dat) <- sp::CRS(sp::proj4string(as(pts_utm, "Spatial")))


ols <- lm(y ~ x1 + x2 + x3 + x4, data = dat)

sh_ols <- shapiro.test(residuals(ols))     # normalitas residual
bp     <- lmtest::bptest(ols)             # heteroskedastisitas
vif_v  <- car::vif(ols)                   # multikolinearitas

# Moran residual OLS (kNN k=4)
xy <- cbind(dat$easting, dat$northing)
nb_knn <- spdep::knn2nb(spdep::knearneigh(xy, k = 4))
lw_knn <- spdep::nb2listw(nb_knn, style = "W", zero.policy = TRUE)
mor_ols <- spdep::moran.test(residuals(ols), lw_knn, zero.policy = TRUE)

bw_adapt <- spgwr::gwr.sel(y ~ x1 + x2 + x3 + x4, data = dat, adapt = TRUE)
## Adaptive q: 0.381966 CV score: 185579.9 
## Adaptive q: 0.618034 CV score: 189889.1 
## Adaptive q: 0.236068 CV score: 185919.2 
## Adaptive q: 0.330599 CV score: 184527.8 
## Adaptive q: 0.3138362 CV score: 184111.4 
## Adaptive q: 0.2841314 CV score: 183462.3 
## Adaptive q: 0.2657728 CV score: 185218.8 
## Adaptive q: 0.2954776 CV score: 183393.4 
## Adaptive q: 0.2918024 CV score: 183229.6 
## Adaptive q: 0.2902648 CV score: 183158.7 
## Adaptive q: 0.2879221 CV score: 183214.4 
## Adaptive q: 0.2897535 CV score: 183134.9 
## Adaptive q: 0.2894023 CV score: 183125.9 
## Adaptive q: 0.2888369 CV score: 183159.2 
## Adaptive q: 0.2890561 CV score: 183146.2 
## Adaptive q: 0.2894721 CV score: 183121.8 
## Adaptive q: 0.289535 CV score: 183124.6 
## Adaptive q: 0.2894721 CV score: 183121.8
gwr_fit <- spgwr::gwr(
  y ~ x1 + x2 + x3 + x4,
  data      = dat,
  adapt     = bw_adapt,
  hatmatrix = TRUE,
  se.fit    = TRUE
)

sdf <- as.data.frame(gwr_fit$SDF)

# tempel balik hasil lokal ke sf untuk peta kalau perlu
sf_m2$gwr_x1      <- sdf$x1
sf_m2$gwr_x2      <- sdf$x2
sf_m2$gwr_x3      <- sdf$x3
sf_m2$gwr_x4      <- sdf$x4
sf_m2$gwr_localR2 <- sdf$localR2
sf_m2$gwr_pred    <- sdf$pred
sf_m2$gwr_resid   <- sdf$gwr.e

sh_gwr  <- shapiro.test(sf_m2$gwr_resid)
mor_gwr <- spdep::moran.test(sf_m2$gwr_resid, lw_knn, zero.policy = TRUE)

safe1 <- function(x) {
  # pastikan output 1 nilai saja (kalau NULL / panjang 0 -> NA)
  if (is.null(x) || length(x) == 0) return(NA)
  x <- x[1]
  x
}

safe_num1 <- function(expr) {
  # jalankan ekspresi; kalau error -> NA; kalau hasil vektor -> ambil 1
  out <- tryCatch(eval.parent(substitute(expr)), error = function(e) NA)
  safe1(out)
}

ringkas <- tibble::tibble(
  Item = c(
    "n (complete-case)",
    "Bandwidth adapt (GWR)",
    "AICc (GWR)",
    "R2 (GWR results)",
    "Shapiro p (OLS residual)",
    "Breusch-Pagan p (OLS)",
    "Max VIF (OLS)",
    "Moran p (OLS residual)",
    "Shapiro p (GWR residual)",
    "Moran p (GWR residual)"
  ),
  Value = c(
    safe_num1(nrow(sf_m2)),
    safe_num1(bw_adapt),
    safe_num1(gwr_fit$results$AICc),
    safe_num1(gwr_fit$results$r.squared),
    safe_num1(sh_ols$p.value),
    safe_num1(bp$p.value),
    safe_num1(max(vif_v, na.rm = TRUE)),
    safe_num1(mor_ols$p.value),
    safe_num1(sh_gwr$p.value),
    safe_num1(mor_gwr$p.value)
  )
)

ringkas
# ============================================================
# GWR LENGKAP (1 CHUNK): Data -> OLS + Uji Asumsi -> GWR -> Ringkasan -> Peta
# ============================================================

# ---- Packages ----
library(sf)
library(dplyr)
library(sp)
library(spgwr)
library(lmtest)
library(car)
library(spdep)
library(ggplot2)
library(tibble)
library(patchwork)

pick_first_sf <- function(...) {
  cand <- list(...)
  for (obj in cand) if (!is.null(obj) && inherits(obj, "sf")) return(obj)
  stop("Tidak ditemukan objek bertipe sf. Pastikan ada: sf_m / sf_join / joined / Jatim_sf / Jatim_Kabkota.")
}

sf_data <- pick_first_sf(
  if (exists("sf_m")) sf_m else NULL,
  if (exists("sf_join")) sf_join else NULL,
  if (exists("joined")) joined else NULL,
  if (exists("Jatim_sf")) Jatim_sf else NULL,
  if (exists("Jatim_Kabkota")) Jatim_Kabkota else NULL
)

find_col <- function(df, patterns) {
  nm <- names(df)
  for (p in patterns) {
    hit <- nm[grepl(p, nm, ignore.case = TRUE)]
    if (length(hit) > 0) return(hit[1])
  }
  NA_character_
}

# mapping kolom (otomatis; ubah manual kalau perlu)
col_kasus <- if ("Kasus_TBC" %in% names(sf_data)) "Kasus_TBC" else find_col(sf_data, c("^kasus", "tbc"))
col_x1    <- if ("X1" %in% names(sf_data)) "X1" else find_col(sf_data, c("^x1$", "penduduk", "populasi", "jumlah_penduduk"))
col_x2    <- if ("X2" %in% names(sf_data)) "X2" else find_col(sf_data, c("^x2$", "kepadatan", "density"))
col_x3    <- if ("X3" %in% names(sf_data)) "X3" else find_col(sf_data, c("^x3$", "miskin", "kemiskinan", "poverty"))
col_x4    <- if ("X4" %in% names(sf_data)) "X4" else find_col(sf_data, c("^x4$", "rokok", "cig", "smok"))

if (is.na(col_kasus)) stop("Kolom Kasus_TBC tidak ditemukan. Cek nama kolom kasus di data kamu.")
if (is.na(col_x1))    stop("Kolom X1 tidak ditemukan. Cek kolom penduduk (ribu) di data kamu.")
if (is.na(col_x2))    stop("Kolom X2 tidak ditemukan. Cek kolom kepadatan di data kamu.")
if (is.na(col_x3))    stop("Kolom X3 tidak ditemukan. Cek kolom kemiskinan di data kamu.")
if (is.na(col_x4))    stop("Kolom X4 tidak ditemukan. Cek kolom rokok di data kamu.")

sf_data <- sf_data %>%
  dplyr::rename(
    Kasus_TBC = all_of(col_kasus),
    X1 = all_of(col_x1),
    X2 = all_of(col_x2),
    X3 = all_of(col_x3),
    X4 = all_of(col_x4)
  ) %>%
  dplyr::mutate(
    Kasus_TBC = suppressWarnings(as.numeric(Kasus_TBC)),
    X1 = suppressWarnings(as.numeric(X1)),
    X2 = suppressWarnings(as.numeric(X2)),
    X3 = suppressWarnings(as.numeric(X3)),
    X4 = suppressWarnings(as.numeric(X4)),
    Populasi = X1 * 1000,
    IR_100k  = ifelse(!is.na(Populasi) & Populasi > 0 & !is.na(Kasus_TBC),
                      (Kasus_TBC / Populasi) * 100000,
                      NA_real_)
  )

needed <- c("IR_100k","X1","X2","X3","X4")

sf_m2 <- sf_data %>%
  dplyr::filter(
    complete.cases(sf::st_drop_geometry(.)[, needed, drop = FALSE])
  )

if (nrow(sf_m2) < 10) stop("Data valid (complete-case) < 10. Cek NA / variabel.")


pts_sf  <- sf::st_point_on_surface(sf_m2)
pts_utm <- tryCatch(sf::st_transform(pts_sf, 32749), error=function(e) pts_sf)
coords  <- sf::st_coordinates(pts_utm)

dat <- pts_utm %>%
  sf::st_drop_geometry() %>%
  dplyr::transmute(
    y  = IR_100k,
    x1 = X1, x2 = X2, x3 = X3, x4 = X4,
    easting  = coords[,1],
    northing = coords[,2]
  )

sp::coordinates(dat) <- ~ easting + northing
sp::proj4string(dat) <- sp::CRS(sp::proj4string(as(pts_utm, "Spatial")))


ols   <- lm(y ~ x1 + x2 + x3 + x4, data = dat)
sh_ols <- shapiro.test(residuals(ols))
bp     <- lmtest::bptest(ols)
vif_v  <- car::vif(ols)

# Moran residual OLS (kNN k=4)
xy <- cbind(dat$easting, dat$northing)
nb_knn <- spdep::knn2nb(spdep::knearneigh(xy, k = 4))
lw_knn <- spdep::nb2listw(nb_knn, style = "W", zero.policy = TRUE)
mor_ols <- spdep::moran.test(residuals(ols), lw_knn, zero.policy = TRUE)


bw_adapt <- spgwr::gwr.sel(y ~ x1 + x2 + x3 + x4, data = dat, adapt = TRUE)
## Adaptive q: 0.381966 CV score: 185579.9 
## Adaptive q: 0.618034 CV score: 189889.1 
## Adaptive q: 0.236068 CV score: 185919.2 
## Adaptive q: 0.330599 CV score: 184527.8 
## Adaptive q: 0.3138362 CV score: 184111.4 
## Adaptive q: 0.2841314 CV score: 183462.3 
## Adaptive q: 0.2657728 CV score: 185218.8 
## Adaptive q: 0.2954776 CV score: 183393.4 
## Adaptive q: 0.2918024 CV score: 183229.6 
## Adaptive q: 0.2902648 CV score: 183158.7 
## Adaptive q: 0.2879221 CV score: 183214.4 
## Adaptive q: 0.2897535 CV score: 183134.9 
## Adaptive q: 0.2894023 CV score: 183125.9 
## Adaptive q: 0.2888369 CV score: 183159.2 
## Adaptive q: 0.2890561 CV score: 183146.2 
## Adaptive q: 0.2894721 CV score: 183121.8 
## Adaptive q: 0.289535 CV score: 183124.6 
## Adaptive q: 0.2894721 CV score: 183121.8
gwr_fit <- spgwr::gwr(
  y ~ x1 + x2 + x3 + x4,
  data      = dat,
  adapt     = bw_adapt,
  hatmatrix = TRUE,
  se.fit    = TRUE
)

sdf <- as.data.frame(gwr_fit$SDF)

# tempel balik hasil lokal ke sf (untuk peta)
sf_m2$gwr_x1      <- sdf$x1
sf_m2$gwr_x2      <- sdf$x2
sf_m2$gwr_x3      <- sdf$x3
sf_m2$gwr_x4      <- sdf$x4
sf_m2$gwr_localR2 <- sdf$localR2
sf_m2$gwr_pred    <- sdf$pred
sf_m2$gwr_resid   <- sdf$gwr.e

# Uji residual GWR
sh_gwr  <- shapiro.test(sf_m2$gwr_resid)
mor_gwr <- spdep::moran.test(sf_m2$gwr_resid, lw_knn, zero.policy = TRUE)

# ============================================================
# 7) Tabel ringkas (ANTI ERROR ukuran beda)
# ============================================================
safe1 <- function(x) {
  if (is.null(x) || length(x) == 0) return(NA)
  x[1]
}

ringkas <- tibble::tibble(
  Item = c(
    "n (complete-case)",
    "Bandwidth adapt (GWR)",
    "AICc (GWR)",
    "R2 (GWR results)",
    "Shapiro p (OLS residual)",
    "Breusch-Pagan p (OLS)",
    "Max VIF (OLS)",
    "Moran p (OLS residual)",
    "Shapiro p (GWR residual)",
    "Moran p (GWR residual)"
  ),
  Value = c(
    safe1(nrow(sf_m2)),
    safe1(bw_adapt),
    safe1(gwr_fit$results$AICc),
    safe1(gwr_fit$results$r.squared),
    safe1(sh_ols$p.value),
    safe1(bp$p.value),
    safe1(tryCatch(max(vif_v, na.rm = TRUE), error=function(e) NA)),
    safe1(mor_ols$p.value),
    safe1(sh_gwr$p.value),
    safe1(mor_gwr$p.value)
  )
)

cat("\n================ OLS SUMMARY ================\n")
## 
## ================ OLS SUMMARY ================
print(summary(ols))
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -105.329  -45.328   -4.799   26.201  169.222 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 198.928880  47.370817   4.199  0.00019 ***
## x1           -0.045064   0.024295  -1.855  0.07256 .  
## x2            0.032064   0.007191   4.459 8.99e-05 ***
## x3            0.299565   0.299899   0.999  0.32512    
## x4          -14.487678   8.848735  -1.637  0.11108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.26 on 33 degrees of freedom
## Multiple R-squared:  0.6522, Adjusted R-squared:   0.61 
## F-statistic: 15.47 on 4 and 33 DF,  p-value: 3.181e-07
cat("\n================ RINGKASAN (UNTUK LAPORAN) ================\n")
## 
## ================ RINGKASAN (UNTUK LAPORAN) ================
print(ringkas)
## # A tibble: 10 × 2
##    Item                         Value
##    <chr>                        <dbl>
##  1 n (complete-case)         38      
##  2 Bandwidth adapt (GWR)      0.289  
##  3 AICc (GWR)               437.     
##  4 R2 (GWR results)          NA      
##  5 Shapiro p (OLS residual)   0.0773 
##  6 Breusch-Pagan p (OLS)      0.00262
##  7 Max VIF (OLS)              3.34   
##  8 Moran p (OLS residual)     0.386  
##  9 Shapiro p (GWR residual)   0.168  
## 10 Moran p (GWR residual)     0.468
p_r2 <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_localR2), color="grey60", linewidth=0.25) +
  scale_fill_viridis_c(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Local R²", fill="Local R²")

p_x1 <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_x1), color="grey60", linewidth=0.25) +
  scale_fill_viridis_c(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Koef X1", fill="Koef X1")

p_x2 <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_x2), color="grey60", linewidth=0.25) +
  scale_fill_viridis_c(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Koef X2", fill="Koef X2")

p_x3 <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_x3), color="grey60", linewidth=0.25) +
  scale_fill_viridis_c(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Koef X3", fill="Koef X3")

p_x4 <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_x4), color="grey60", linewidth=0.25) +
  scale_fill_viridis_c(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Koef X4", fill="Koef X4")

p_res <- ggplot(sf_m2) +
  geom_sf(aes(fill = gwr_resid), color="grey60", linewidth=0.25) +
  scale_fill_gradient2(na.value="grey95") +
  theme_minimal() + labs(title="GWR: Residual", fill="Residual")

cat("\n================ PETA GWR (GRID DASHBOARD) ================\n")
## 
## ================ PETA GWR (GRID DASHBOARD) ================
print((p_r2 | p_x1 | p_x2) / (p_x3 | p_x4 | p_res))

Model Geographically Weighted Regression (GWR) digunakan untuk menganalisis variasi spasial hubungan antara faktor demografi dan sosial ekonomi terhadap angka insidensi TBC per 100.000 penduduk (IR_100k) di Provinsi Jawa Timur tahun 2024. Hasil analisis menunjukkan bahwa hubungan antara variabel independen dan IR TBC tidak bersifat homogen, melainkan berbeda antar kabupaten/kota. Hal ini menegaskan bahwa model global kurang mampu menangkap kompleksitas spasial persebaran TBC, sehingga pendekatan GWR lebih sesuai untuk menggambarkan dinamika lokal.

Berdasarkan peta Local R², terlihat bahwa kemampuan model dalam menjelaskan variasi IR TBC bervariasi secara spasial. Wilayah barat dan barat daya Jawa Timur memiliki nilai Local R² yang relatif lebih tinggi, yang mengindikasikan bahwa variabel-variabel dalam model mampu menjelaskan variasi IR TBC dengan baik di wilayah tersebut. Sebaliknya, wilayah timur dan kepulauan menunjukkan nilai Local R² yang lebih rendah, yang mengindikasikan adanya faktor lain di luar model yang turut memengaruhi tingginya atau rendahnya IR TBC di wilayah tersebut.

Koefisien lokal jumlah penduduk (X1) menunjukkan nilai negatif di seluruh wilayah Jawa Timur, dengan intensitas yang lebih besar di wilayah barat. Hal ini mengindikasikan bahwa peningkatan jumlah penduduk cenderung berasosiasi dengan penurunan IR TBC, terutama di wilayah dengan konsentrasi penduduk yang tinggi. Kondisi ini dapat mencerminkan peran wilayah perkotaan dengan populasi besar yang umumnya memiliki fasilitas kesehatan, sistem deteksi, dan pelaporan kasus yang lebih baik dibandingkan wilayah dengan jumlah penduduk lebih kecil.

Variabel kepadatan penduduk (X2) menunjukkan hubungan positif terhadap IR TBC di seluruh wilayah, meskipun dengan kekuatan yang berbeda. Pengaruh kepadatan penduduk terhadap peningkatan IR TBC lebih kuat di wilayah barat dan tengah Jawa Timur dibandingkan wilayah timur. Hal ini sejalan dengan teori epidemiologi penyakit menular, di mana kepadatan penduduk yang tinggi meningkatkan intensitas kontak antarindividu dan memperbesar peluang penularan TBC.

Koefisien lokal jumlah penduduk miskin (X3) juga menunjukkan hubungan positif dengan IR TBC, terutama di wilayah barat dan selatan Jawa Timur. Temuan ini mengindikasikan bahwa kemiskinan merupakan determinan penting dalam persebaran TBC, karena berkaitan erat dengan kondisi hunian yang kurang layak, status gizi yang rendah, serta keterbatasan akses terhadap layanan kesehatan. Sebaliknya, pengaruh kemiskinan terhadap IR TBC relatif lebih lemah di wilayah timur dan kepulauan.

Sementara itu, koefisien lokal konsumsi rokok filter (X4) menunjukkan hubungan negatif terhadap IR TBC di seluruh wilayah Jawa Timur, dengan variasi spasial yang cukup signifikan. Pengaruh negatif yang lebih kuat ditemukan di wilayah barat dan tengah Jawa Timur, sedangkan wilayah timur menunjukkan pengaruh yang lebih lemah. Hasil ini perlu ditafsirkan secara hati-hati, karena kemungkinan konsumsi rokok filter berperan sebagai proksi karakteristik sosial ekonomi atau gaya hidup wilayah tertentu, dan tidak serta-merta mencerminkan hubungan kausal langsung terhadap penurunan risiko TBC.

Secara keseluruhan, hasil GWR menunjukkan bahwa faktor-faktor yang memengaruhi IR TBC di Provinsi Jawa Timur tahun 2024 bersifat spasial dan kontekstual. Variasi pengaruh antar wilayah menegaskan pentingnya pendekatan kebijakan pengendalian TBC yang berbasis wilayah, di mana intervensi kesehatan masyarakat perlu disesuaikan dengan karakteristik lokal masing-masing kabupaten/kota untuk mencapai efektivitas yang optimal.

4.5 Analisis Interpolasi (IDW & Krigging

library(gstat)
library(patchwork)

pick_first_sf <- function(...) {
  cand <- list(...)
  for (obj in cand) if (!is.null(obj) && inherits(obj, "sf")) return(obj)
  stop("Tidak ditemukan objek bertipe sf. Pastikan ada: sf_m2/sf_m/sf_join/joined/Jatim_sf/Jatim_Kabkota.")
}

sf_data <- pick_first_sf(
  if (exists("sf_m2")) sf_m2 else NULL,
  if (exists("sf_m")) sf_m else NULL,
  if (exists("sf_join")) sf_join else NULL,
  if (exists("joined")) joined else NULL,
  if (exists("Jatim_sf")) Jatim_sf else NULL,
  if (exists("Jatim_Kabkota")) Jatim_Kabkota else NULL
)

find_col <- function(df, patterns) {
  nm <- names(df)
  for (p in patterns) {
    hit <- nm[grepl(p, nm, ignore.case = TRUE)]
    if (length(hit) > 0) return(hit[1])
  }
  NA_character_
}

col_kasus <- if ("Kasus_TBC" %in% names(sf_data)) "Kasus_TBC" else find_col(sf_data, c("^kasus", "tbc"))
col_x1    <- if ("X1" %in% names(sf_data)) "X1" else find_col(sf_data, c("^x1$", "penduduk", "populasi", "jumlah_penduduk"))

if (!is.na(col_kasus) && col_kasus != "Kasus_TBC") sf_data <- sf_data %>% dplyr::rename(Kasus_TBC = all_of(col_kasus))
if (!is.na(col_x1)    && col_x1    != "X1")      sf_data <- sf_data %>% dplyr::rename(X1 = all_of(col_x1))

if (!("IR_100k" %in% names(sf_data))) {
  if (!("Kasus_TBC" %in% names(sf_data)) || !("X1" %in% names(sf_data))) {
    stop("IR_100k belum ada dan tidak bisa dihitung karena Kasus_TBC/X1 tidak ditemukan.")
  }
  sf_data <- sf_data %>%
    dplyr::mutate(
      Kasus_TBC = suppressWarnings(as.numeric(Kasus_TBC)),
      X1 = suppressWarnings(as.numeric(X1)),
      Populasi = X1 * 1000,
      IR_100k  = ifelse(!is.na(Populasi) & Populasi > 0 & !is.na(Kasus_TBC),
                        (Kasus_TBC / Populasi) * 100000,
                        NA_real_)
    )
}

sf_m <- sf_data %>% dplyr::filter(!is.na(IR_100k))
if (nrow(sf_m) < 5) stop("Data valid untuk interpolasi < 5. Cek IR_100k / NA.")
grid_n    <- 120   # mirip slider dashboard (50-200). makin besar makin berat
idp_power <- 2     # IDW power

pts_sf <- sf::st_point_on_surface(sf_m)
pts_m  <- tryCatch(sf::st_transform(pts_sf, 32749), error=function(e) pts_sf)
sf_m_m <- tryCatch(sf::st_transform(sf_m, 32749), error=function(e) sf_m)

pts_sp <- as(pts_m, "Spatial")
pts_sp$IR_100k <- sf_m_m$IR_100k

# grid (centers)
bb <- sf::st_bbox(sf_m_m)
grd_sf <- sf::st_make_grid(sf::st_as_sfc(bb), n = c(grid_n, grid_n), what = "centers") %>%
  sf::st_as_sf()
sf::st_crs(grd_sf) <- sf::st_crs(sf_m_m)
grd_sp <- as(grd_sf, "Spatial")


idw_res <- gstat::idw(IR_100k ~ 1, locations = pts_sp, newdata = grd_sp, idp = idp_power)
## [inverse distance weighted interpolation]
idw_df  <- as.data.frame(idw_res)
names(idw_df)[names(idw_df) == "var1.pred"] <- "pred"

tab_idw <- tibble::tibble(
  Method = "IDW",
  Parameter = c("Power (idp)", "Grid N (N x N)", "n titik input"),
  Value = c(idp_power, paste0(grid_n, " x ", grid_n), length(pts_sp))
)

p_idw <- ggplot(idw_df, aes(x = coords.x1, y = coords.x2, fill = pred)) +
  geom_raster(interpolate = TRUE) +
  geom_sf(data = sf_m_m, fill = NA, color = "grey60", linewidth = 0.25, inherit.aes = FALSE) +
  theme_minimal() +
  labs(title = "Interpolasi IDW (p=2) — Prediksi IR_100k", fill = "IR_100k")

# 1) Variogram empiris
vgm_emp <- gstat::variogram(IR_100k ~ 1, data = pts_sp)

# 2) Fit variogram (model spherical) - sama spirit dashboard kamu
vgm_fit <- tryCatch(
  gstat::fit.variogram(
    vgm_emp,
    model = gstat::vgm(
      psill  = stats::var(pts_sp$IR_100k, na.rm = TRUE),
      model  = "Sph",
      range  = max(vgm_emp$dist, na.rm = TRUE) / 2,
      nugget = 0
    ),
    fit.method = 2
  ),
  error = function(e) NULL
)
if (is.null(vgm_fit)) stop("Fitting variogram gagal. Cek CRS/data/NA.")

# 3) Kriging
kr <- gstat::krige(IR_100k ~ 1, locations = pts_sp, newdata = grd_sp, model = vgm_fit)
## [using ordinary kriging]
kr_df <- as.data.frame(kr)
names(kr_df)[names(kr_df) == "var1.pred"] <- "pred"
names(kr_df)[names(kr_df) == "var1.var"]  <- "var"

# 4) LOOCV RMSE (opsional tapi sesuai ringkasan dashboard)
cv <- tryCatch(
  gstat::krige.cv(IR_100k ~ 1, locations = pts_sp, model = vgm_fit, nfold = length(pts_sp)),
  error = function(e) NULL
)
rmse <- if (!is.null(cv)) sqrt(mean((cv$observed - cv$var1.pred)^2, na.rm = TRUE)) else NA_real_

tab_krig <- tibble::tibble(
  Method = "Ordinary Kriging",
  Parameter = c("Model", "Nugget", "Partial Sill (psill)", "Range", "Grid N (N x N)", "n titik input", "CV RMSE (LOOCV)"),
  Value = c(as.character(vgm_fit$model[2]),
            vgm_fit$psill[1],
            vgm_fit$psill[2],
            vgm_fit$range[2],
            paste0(grid_n, " x ", grid_n),
            length(pts_sp),
            rmse)
)

p_kr_pred <- ggplot(kr_df, aes(x = coords.x1, y = coords.x2, fill = pred)) +
  geom_raster(interpolate = TRUE) +
  geom_sf(data = sf_m_m, fill = NA, color = "grey60", linewidth = 0.25, inherit.aes = FALSE) +
  theme_minimal() +
  labs(title = "Ordinary Kriging — Prediksi IR_100k", fill = "IR_100k")

p_kr_var <- ggplot(kr_df, aes(x = coords.x1, y = coords.x2, fill = var)) +
  geom_raster(interpolate = TRUE) +
  geom_sf(data = sf_m_m, fill = NA, color = "grey60", linewidth = 0.25, inherit.aes = FALSE) +
  theme_minimal() +
  labs(title = "Ordinary Kriging — Varian (Uncertainty)", fill = "Varian")

# Variogram plot (panel diagnostik)
vgm_df <- as.data.frame(vgm_emp)
p_vgm <- ggplot(vgm_df, aes(x = dist, y = gamma)) +
  geom_point() +
  theme_minimal() +
  labs(title = "Empirical Variogram", x = "Distance", y = "Semivariance")


cat("\n================ RINGKASAN IDW ================\n")
## 
## ================ RINGKASAN IDW ================
print(tab_idw)
## # A tibble: 3 × 3
##   Method Parameter      Value    
##   <chr>  <chr>          <chr>    
## 1 IDW    Power (idp)    2        
## 2 IDW    Grid N (N x N) 120 x 120
## 3 IDW    n titik input  38
cat("\n================ RINGKASAN KRIGING ================\n")
## 
## ================ RINGKASAN KRIGING ================
print(tab_krig)
## # A tibble: 7 × 3
##   Method           Parameter            Value           
##   <chr>            <chr>                <chr>           
## 1 Ordinary Kriging Model                Sph             
## 2 Ordinary Kriging Nugget               10994.1702583597
## 3 Ordinary Kriging Partial Sill (psill) 79.1936137455226
## 4 Ordinary Kriging Range                200041.191804332
## 5 Ordinary Kriging Grid N (N x N)       120 x 120       
## 6 Ordinary Kriging n titik input        38              
## 7 Ordinary Kriging CV RMSE (LOOCV)      105.875378770452
cat("\n================ PETA INTERPOLASI (DASHBOARD GRID) ================\n")
## 
## ================ PETA INTERPOLASI (DASHBOARD GRID) ================
# Baris 1: IDW prediksi | Kriging prediksi
# Baris 2: Kriging varian | Variogram
print((p_idw | p_kr_pred) / (p_kr_var | p_vgm))

Hasil interpolasi Inverse Distance Weighting (IDW) dengan parameter power (idp) = 2 dan resolusi grid 120 × 120 menunjukkan bahwa sebaran nilai angka insidensi TBC per 100.000 penduduk (IR_100k) di Provinsi Jawa Timur membentuk pola konsentrasi lokal di sekitar wilayah dengan nilai pengamatan tinggi. Peta interpolasi IDW memperlihatkan adanya beberapa hotspot IR TBC yang dominan di wilayah tengah dan barat Jawa Timur, sementara wilayah lain menunjukkan nilai yang relatif lebih rendah. Karena metode IDW sangat bergantung pada jarak titik pengamatan, pola yang dihasilkan cenderung menampilkan perubahan nilai yang tajam di sekitar titik data dan lebih halus di area yang jauh dari titik observasi. Dengan demikian, IDW lebih merepresentasikan pola lokal berbasis kedekatan spasial tanpa mempertimbangkan struktur autokorelasi spasial secara eksplisit.

Sementara itu, hasil Ordinary Kriging memberikan gambaran sebaran IR TBC yang lebih halus dan terstruktur dibandingkan IDW. Peta prediksi Kriging menunjukkan adanya gradien spasial yang relatif konsisten, di mana wilayah dengan IR TBC tinggi tidak hanya terfokus pada titik tertentu, tetapi menyebar mengikuti pola spasial yang dipengaruhi oleh struktur variogram. Hal ini mengindikasikan bahwa metode Kriging mampu menangkap ketergantungan spasial antar wilayah secara lebih baik dibandingkan IDW, sehingga menghasilkan estimasi yang lebih stabil pada area yang tidak memiliki data pengamatan langsung.

Berdasarkan hasil pemodelan variogram pada Ordinary Kriging, digunakan model spherical (Sph) dengan nilai nugget yang relatif besar, yang mengindikasikan adanya variasi skala kecil atau kesalahan pengukuran yang cukup dominan pada data IR TBC. Nilai partial sill yang lebih kecil dibandingkan nugget menunjukkan bahwa komponen variasi spasial terstruktur masih ada, namun tidak mendominasi keseluruhan variasi data. Nilai range sekitar 200 km mengindikasikan bahwa autokorelasi spasial IR TBC masih bertahan hingga jarak tersebut, sehingga wilayah-wilayah dalam radius tersebut cenderung memiliki nilai IR TBC yang saling berkaitan.

Hasil validasi silang (Leave-One-Out Cross Validation) pada Ordinary Kriging menghasilkan nilai RMSE sekitar 105,88, yang menunjukkan tingkat kesalahan prediksi yang masih cukup besar. Hal ini mengindikasikan bahwa meskipun Kriging mampu menangkap pola spasial secara lebih baik dibandingkan IDW, akurasi prediksi masih dipengaruhi oleh keterbatasan jumlah titik pengamatan dan tingginya variasi lokal kasus TBC antar wilayah. Oleh karena itu, hasil interpolasi sebaiknya digunakan sebagai gambaran umum pola spasial dan bukan sebagai estimasi nilai absolut yang presisi.

Secara keseluruhan, interpolasi IDW dan Ordinary Kriging sama-sama menunjukkan adanya pola spasial persebaran IR TBC di Provinsi Jawa Timur tahun 2024, namun dengan karakteristik yang berbeda. Metode IDW lebih menonjolkan pengaruh lokal dari titik pengamatan terdekat, sedangkan Ordinary Kriging memberikan estimasi yang lebih halus dan mempertimbangkan struktur spasial data. Dengan demikian, hasil interpolasi ini dapat dimanfaatkan sebagai alat eksplorasi spasial untuk mengidentifikasi wilayah berisiko tinggi TBC, namun tetap perlu dikombinasikan dengan analisis regresi spasial untuk penarikan kesimpulan kebijakan yang lebih kuat.

BAB V

KESIMPULAN DAN SARAN

Berdasarkan hasil analisis spasial kasus Tuberkulosis (TBC) di Provinsi Jawa Timur tahun 2024, dapat disimpulkan bahwa distribusi kasus TBC menunjukkan pola yang tidak merata antar kabupaten/kota. Wilayah dengan jumlah penduduk besar dan tingkat kepadatan tinggi, khususnya kawasan perkotaan dan daerah penyangga metropolitan, cenderung memiliki jumlah kasus absolut yang lebih tinggi dibandingkan wilayah lain. Hasil uji autokorelasi spasial mengindikasikan adanya keterkaitan spasial antarwilayah, di mana daerah dengan nilai kasus tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki risiko tinggi. Hal ini menunjukkan bahwa penularan TBC tidak terjadi secara acak, melainkan membentuk pola spasial tertentu yang dipengaruhi oleh faktor lingkungan dan karakteristik wilayah. Selain itu, penerapan model spasial (SAR, SEM, dan SDM) menunjukkan bahwa sar menjadi model terbaik Secara keseluruhan, hasil penelitian diharapkan dapat menjadi dasar bagi pemerintah daerah dan pemangku kebijakan dalam merancang strategi pengendalian TBC yang lebih efektif, terarah, dan berbasis risiko wilayah.

DAFTAR PUSTAKA

[1] World Health Organization. Global Tuberculosis Report 2023. Geneva: WHO; 2023. [2]Kementerian Kesehatan RI. Laporan Program Penanggulangan Tuberkulosis (TBC) di Indonesia. Jakarta: Direktorat Jenderal Pencegahan dan Pengendalian Penyakit; 2024. [3]Pfeiffer DU, Robinson TP, Stevenson M, et al. Spatial Analysis in Epidemiology. Oxford: Oxford University Press; 2008. [4]Dinas Kesehatan Provinsi Jawa Timur. Profil Kesehatan Provinsi Jawa Timur Tahun 2023 (Data Update 2024). Surabaya: Dinkes Jatim; 2024. [5]Badan Pusat Statistik Provinsi Jawa Timur. Provinsi Jawa Timur Dalam Angka 2024. Surabaya: BPS Jatim; 2024. [6]Anselin L. Interactive techniques and exploratory spatial data analysis. In: Longley PA, Goodchild MF, Maguire DJ, Rhind DW, editors. Geographical Information Systems: Principles, Techniques, Management and Applications. 2nd ed. Hoboken: John Wiley & Sons; 2010. Tobler WR. A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography. 1970;46:234-240. [7]Goodchild MF. Spatial Autocorrelation. Norwich: Geo Books; 1986. [8]Haining RP. Spatial Data Analysis: Theory and Practice. Cambridge: Cambridge University Press; 2003. [9]Cliff AD, Ord JK. Spatial Processes: Models & Applications. London: Pion; 1981. [10]Odland J. Spatial Autocorrelation. Beverly Hills: Sage Publications; 1988. [11]Moran PAP. Notes on Continuous Stochastic Phenomena. Biometrika. 1950;37(1/2):17-23. [12]Getis A, Ord JK. The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis. 1992;24(3):189-206. [13]Anselin L. Spatial Econometrics: Methods and Models. Dordrecht: Kluwer Academic Publishers; 1988. [14]LeSage J, Pace RK. Introduction to Spatial Econometrics. Boca Raton: CRC Press Taylor & Francis Group; 2009. [15]Elhorst JP. Spatial Econometrics: From Cross-Sectional Data to Spatial [16]Panels. Berlin: Springer; 2014.