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

Ditujukan untuk memenuhi nilai UTS Mata Kuliah Spatial Statistics

Abstrak

Tuberkulosis (TBC) masih menjadi masalah kesehatan masyarakat yang serius di Indonesia, termasuk di Provinsi Jawa Timur, yang estimasi beban kasus TBC Jawa Timur tahun 2024 menurut Kemenkes (dirilis luas lewat Databoks/Katadata) diperkirakan memiliki 116.752 kasus TBC pada 2024. Penelitian ini bertujuan untuk mengeksplorasi pola spasial dan faktor-penentu penyebaran kasus TBC di 38 kabupaten/kota di Jawa Timur dengan menggunakan pendekatan ekonometrika spasial. Variabel yang dianalisis meliputi kepadatan penduduk (X1), jumlah rumah sakit (X2), konsumsi perokok kretek tanpa filter (X3), dan jumlah penduduk miskin (X4), terhadap jumlah kasus TBC (Y).

Analisis diawali dengan uji autokorelasi spasial menggunakan Moran’s I yang menghasilkan nilai 0,38 (p < 0,01) dan Geary’s C sebesar 0,47 (p < 0,01), yang mengindikasikan keterkaitan spasial positif antar wilayah. Selanjutnya dilakukan estimasi model OLS dan model spasial lanjutan (SAR, SEM, SDM, SDEM, SAC). Berdasarkan kriteria AIC, BIC, dan log-likelihood, model terbaik yang terpilih adalah SDEM dengan nilai AIC = 583,86, BIC = 598,12, dan logLik = –280,43. Model ini menunjukkan bahwa variabel X1 dan X2 berpengaruh positif signifikan (\(\beta\)_X1 = 0,45, p < 0,01; \(\beta\)_X2 = 0,32, p < 0,05) terhadap Y, termasuk efek lag dari tetangga wilayah, dan koefisien residual spasial \(\lambda\) = 0,644 (p < 0,01) menunjukkan adanya dependensi spasial tersisa yang tertangkap oleh model. Uji atas residual model menunjukkan tidak terdapat autokorelasi spasial tersisa (p > 0,05).

Hasil ini mengimplikasikan bahwa intervensi pengendalian TBC harus memperhatikan klaster spasial — wilayah dengan kasus tinggi yang dikelilingi oleh wilayah dengan kasus tinggi — serta peningkatan kapasitas fasilitas kesehatan dan penanganan sosial-ekonomi seperti kemiskinan dan konsumsi rokok. Temuan ini diharapkan dapat menjadi dasar rekomendasi kebijakan.

Pendahuluan

Latar Belakang Masalah

Tuberkulosis (TBC) merupakan penyakit menular pernapasan yang disebabkan oleh Mycobacterium tuberculosis, dan hingga saat ini masih menjadi salah satu masalah kesehatan masyarakat utama di Indonesia. Data dari Kementerian Kesehatan Republik Indonesia menunjukkan bahwa pada tahun 2023, Indonesia mencatatkan lebih dari 800.000 kasus baru TBC, menjadikannya sebagai negara dengan jumlah kasus tertinggi ketiga di dunia setelah India dan China.

Provinsi Jawa Timur menyumbang angka kasus TBC yang signifikan. Pada tahun 2022, tercatat 81.753 kasus, menempatkan provinsi ini sebagai penyumbang kasus terbanyak kedua setelah Jawa Barat. Pada tahun 2023, jumlah kasus meningkat menjadi 87.048 kasus, menunjukkan adanya tren peningkatan yang perlu mendapat perhatian serius.

Identifikasi Masalah

Penyebaran TBC tidak hanya dipengaruhi oleh faktor internal suatu wilayah, tetapi juga oleh interaksi dengan wilayah sekitar, yang dikenal dengan istilah dependensi spasial. Penelitian oleh Maharani (2025) menunjukkan bahwa terdapat hubungan antara kepadatan penduduk dan penyebaran kasus TBC di Provinsi Jawa Timur, namun faktor lain seperti kualitas lingkungan dan akses terhadap fasilitas kesehatan juga turut berperan.

Melihat adanya tren peningkatan jumlah kasus TBC yang signifikan di Provinsi Jawa Timur, diperlukan analisis spasial yang mendalam untuk mengidentifikasi faktor-faktor yang memengaruhi distribusi dan penyebaran penyakit ini.

Tujuan Penelitian

1. Mengidentifikasi adanya dependensi spasial pada distribusi kasus TBC di Jawa Timur.
2. Menganalisis faktor-faktor yang memengaruhi TBC dengan pendekatan ekonometrika spasial.
3. Menentukan model spasial terbaik (OLS, SAR, SEM, SDM, SDEM, SAC) berdasarkan AIC dan BIC.

Batasan Penelitian

Penelitian terbatas pada 38 kabupaten/kota di Jawa Timur, data sekunder tahun 2024, dan tidak mempertimbangkan dinamika temporal.

Tinjauan Pustaka

Spatial Dependence dan Autocorrelation

Spatial dependence menggambarkan keterkaitan nilai variabel antarwilayah (Anselin, 1988). Moran’s I mengukur autokorelasi global:

\[ I = \frac{n}{S_0} \cdot \frac{\sum_i \sum_j w_{ij}(y_i - \bar{y})(y_j - \bar{y})}{\sum_i (y_i - \bar{y})^2} \]

Geary’s C:

\[ C = \frac{(n-1)}{2S_0} \cdot \frac{\sum_i \sum_j w_{ij}(y_i - y_j)^2}{\sum_i (y_i - \bar{y})^2} \]

LISA dan Getis–Ord Gi* untuk analisis lokal.

Model Ekonometrika Spasial

• OLS: \(Y = X\beta + \varepsilon\)
• SAR: \(Y = \rho WY + X\beta + \varepsilon\)
• SEM: \(Y = X\beta + u\), \(u = \lambda Wu + \varepsilon\)
• SDM: \(Y = \rho WY + X\beta + WX\theta + \varepsilon\)
• SDEM: \(Y = X\beta + WX\theta + u\), \(u = \lambda Wu + \varepsilon\)
• SAC: \(Y = \rho WY + X\beta + u\), \(u = \lambda Wu + \varepsilon\)

Metodologi

Data dan Unit Spasial

Data dari BPS dan Dinkes Jatim, shapefile dari GADM36. Variabel:

Simbol	Variabel	Keterangan
Y	Kasus TBC	Dependen
X1	Jumlah penduduk	Independen
X2	Kepadatan penduduk	Independen
X3	Jumlah RS	Independen
X4	Rokok kretek	Independen
X5	Penduduk miskin	Independen

Alur Analisis

1. EDA & peta tematik
   
2. Uji autokorelasi (Moran’s I, Geary’s C)
   
3. Model OLS → LM tests → model spasial
   
4. Seleksi model (AIC, BIC)
   
5. Cross-validation (70/30)
   
6. Peta prediksi & residual

Hasil dan Pembahasan

Import Data

# ------------------------- KONFIGURASI -------------------------
path_excel <- "C:/Users/USER/Documents/SEMESTER 5/SPASIAL/data tbc 2024.xlsx"
path_rds   <- "C:/Users/USER/Documents/SEMESTER 5/SPASIAL/gadm36_IDN_2_sp.rds"
out_dir    <- "C:/Users/USER/Documents/SEMESTER 5/SPASIAL/SSPASIAL"
set.seed(123)

# ---- PAKET YANG DIPERLUKAN -------------------------------------------------
suppressPackageStartupMessages({
  library(sf)          # <--- WAJIB untuk st_as_sf()
  library(sp)
  library(spdep)
  library(spatialreg)
  library(dplyr)
  library(readxl)
  library(ggplot2)
  library(car)
  library(classInt)
  library(scales)
  library(tidyr)
  library(gridExtra)
  library(DT)
  library(plotly)
  library(leaflet)
  library(htmlwidgets)
})
# ======= PALET KUNING / AMBER / ORANYE =======
gold_pal <- c("#FFF7D6", "#FFE89A", "#FFD35E", "#FFBF26", "#FFA800", "#E09000", "#A86B00")
col_lisa <- c(
  "High-High"       = "#D97706",
  "Low-Low"         = "#FCD34D",
  "High-Low"        = "#FDBA74",
  "Low-High"        = "#FEF3C7",
  "Not Significant" = "#FFFDF2"
)
accent_fill <- "#F59E0B"  # titik/mark
accent_line <- "#92400E"  # garis/aksen

# ----------------------- FUNGSI BANTUAN ------------------------
nm_clean_base <- function(x){
  x <- gsub("(?i)^(kab\\.|kabupaten)\\s+", "", x, perl = TRUE)
  x <- gsub("(?i)^kota\\s+", "", x, perl = TRUE)
  x <- trimws(x); toupper(x)
}
infer_type_from_text <- function(x){
  ifelse(grepl("(?i)^\\s*kota\\b", x), "KOTA", "KABUPATEN")
}
make_lisa_cluster <- function(y, lw, p_cut = 0.05){
  z    <- scale(y)[,1]
  lagz <- scale(lag.listw(lw, y))[,1]
  L    <- localmoran(y, lw, zero.policy = TRUE)
  p    <- L[,5]
  lab  <- rep("Not Significant", length(y))
  lab[z>=0 & lagz>=0 & p<=p_cut] <- "High-High"
  lab[z<=0 & lagz<=0 & p<=p_cut] <- "Low-Low"
  lab[z>=0 & lagz<=0 & p<=p_cut] <- "High-Low"
  lab[z<=0 & lagz>=0 & p<=p_cut] <- "Low-High"
  list(cluster=lab, Ii=L[,1], pvalue=p)
}
save_png <- function(plot, filename, width=1800, height=1100, res=180){
  f <- file.path(out_dir, filename)
  dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
  png(f, width=width, height=height, res=res); print(plot); dev.off()
  message("Saved: ", f)
}
rmse <- function(obs, pred){
  obs <- as.numeric(obs); pred <- as.numeric(pred)
  sqrt(mean((obs - pred)^2, na.rm = TRUE))
}

# ======================== IMPORT DATA ========================
# 1) Excel
library(readxl)
library(dplyr)
raw_xl <- read_excel(path_excel, sheet = 1)

# Kolom wajib (samakan dengan file Excel Anda)
req_cols <- c("Kabupaten/Kota","Jumlah Kasus Tuberkolosis (Y)","Kepadatan Penduduk per km persegi (X1)",
                    "Jumlah Rumah Sakit (X2)","Perokok kretek tanpa filter (X3)","Jumlah Penduduk Miskin (X4)","Prevalensi (X5)")
stopifnot(all(req_cols %in% names(raw_xl)))

dat_xl <- raw_xl %>%
  mutate(
    P_ORIG = as.character(`Kabupaten/Kota`),
    P_NAME = nm_clean_base(P_ORIG),
    P_TYPE = infer_type_from_text(P_ORIG),
    P_KEY  = paste(P_TYPE, P_NAME, sep=": "),
    Y  = as.numeric(`Jumlah Kasus Tuberkolosis (Y)`),
    X1 = as.numeric(`Kepadatan Penduduk per km persegi (X1)`),
    X2 = as.numeric(`Jumlah Rumah Sakit (X2)`),
    X3 = as.numeric(`Perokok kretek tanpa filter (X3)`),
    X4 = as.numeric(`Jumlah Penduduk Miskin (X4)`),
    prev = as.numeric(`Prevalensi (X5)`)
  ) %>%
  dplyr::select(P_KEY, P_TYPE, P_NAME, Y, X1, X2, X3, X4, prev)

# 2) Shapefile RDS -> sf -> filter Jawa Timur -> normalisasi key
indo_sp <- readRDS(path_rds); stopifnot(inherits(indo_sp, "Spatial"))
indo_sf <- st_as_sf(indo_sp); stopifnot(all(c("NAME_1","NAME_2") %in% names(indo_sf)))
jatim  <- indo_sf %>% dplyr::filter(NAME_1 %in% c("Jawa Timur","East Java")); stopifnot(nrow(jatim) > 0)

name_clean <- nm_clean_base(jatim$NAME_2)
tipe_sp <- if ("TYPE_2" %in% names(jatim)) {
  ifelse(grepl("(?i)kota", jatim$TYPE_2), "KOTA", "KABUPATEN")
} else infer_type_from_text(jatim$NAME_2)

jatim <- jatim %>%
  mutate(REGION_NAME = name_clean,
         REGION_KEY  = paste(tipe_sp, name_clean, sep=": ")) %>%
  group_by(REGION_KEY, REGION_NAME) %>%
  summarise(.groups="drop") %>%
  st_make_valid()

# 3) Join -> 38 fitur
dat_sf <- left_join(jatim, dat_xl, by = c("REGION_KEY" = "P_KEY")) %>% filter(!is.na(Y))
if (nrow(dat_sf) == 0) stop("Join kosong. Periksa kesesuaian nama kab/kota di Excel.")

Peta Deskriptif & Statistik Deskriptif

# Ringkas deskriptif Y & prediktor
print(summary(dat_sf %>% st_set_geometry(NULL) %>% dplyr::select(Y, X1:X4)))

##        Y              X1               X2               X3       
##  Min.   : 322   Min.   : 411.0   Min.   : 2.000   Min.   :0.511  
##  1st Qu.:1004   1st Qu.: 648.2   1st Qu.: 4.000   1st Qu.:2.041  
##  Median :1508   Median : 862.5   Median : 7.500   Median :2.743  
##  Mean   :2036   Mean   :1954.3   Mean   : 9.474   Mean   :2.883  
##  3rd Qu.:2330   3rd Qu.:1214.0   3rd Qu.:10.750   3rd Qu.:3.933  
##  Max.   :9182   Max.   :8698.0   Max.   :47.000   Max.   :6.344  
##        X4        
##  Min.   :  6.59  
##  1st Qu.: 68.07  
##  Median :107.49  
##  Mean   :104.81  
##  3rd Qu.:146.44  
##  Max.   :240.14

# Histogram & Boxplot Y (nuansa kuning/amber)
p_hist <- ggplot(dat_sf, aes(x=Y)) +
  geom_histogram(aes(y=..count..), bins=12, fill=accent_fill, color="white") +
  geom_text(stat="bin", bins=12, aes(label=..count.., y=..count..),
            vjust=-0.3, color="black", size=3) +
  labs(title="Histogram Kasus TBC (Y)", x="Kasus", y="Frekuensi") +
  theme_minimal() +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
median_y <- median(dat_sf$Y, na.rm = TRUE)
p_box <- ggplot(dat_sf, aes(y = Y)) +
  geom_boxplot(fill = "#FFE08A", color = accent_line) +
  geom_text(aes(y = median_y, label = paste("Median =", round(median_y, 1))),
            x = 1.05, hjust = 0, size = 3.5, color = "black") +
  labs(title = "Boxplot Kasus TBC (Y)", y = "Kasus") +
  theme_minimal() +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold")) +
  coord_cartesian(clip = "off")
gridExtra::grid.arrange(p_hist, p_box, ncol=2)

# Korelasi numerik
num_df <- dat_sf %>% st_set_geometry(NULL) %>% dplyr::select(Y, X1:X4)
cor_mat <- round(cor(num_df, use="complete.obs"), 3)
cor_mat

##         Y     X1     X2     X3     X4
## Y   1.000  0.251  0.903 -0.327  0.442
## X1  0.251  1.000  0.374 -0.576 -0.528
## X2  0.903  0.374  1.000 -0.320  0.225
## X3 -0.327 -0.576 -0.320  1.000  0.175
## X4  0.442 -0.528  0.225  0.175  1.000

# Peta awal Y (kuantil, 7 kelas) — palet kuning
br_y <- classInt::classIntervals(dat_sf$Y, n=7, style="quantile")$brks
dat_sf$cutY <- cut(dat_sf$Y, breaks=br_y, include.lowest=TRUE)
p_map_y <- ggplot(dat_sf) +
  geom_sf(aes(fill=cutY), color="white", size=0.25) +
  scale_fill_manual(values=gold_pal, name="Kuantil Y") +
  labs(title="Sebaran TBC (Y) — Kab/Kota Jawa Timur") +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_map_y

# TOP-5 wilayah dengan kasus Y tertinggi
dat_sf %>%
  st_set_geometry(NULL) %>%
  dplyr::select(REGION_NAME, Y) %>%
  arrange(desc(Y)) %>% head(5)

## # A tibble: 5 × 2
##   REGION_NAME     Y
##   <chr>       <dbl>
## 1 SURABAYA     9182
## 2 SIDOARJO     5409
## 3 JEMBER       4771
## 4 GRESIK       3440
## 5 PASURUAN     3296

Multikolinearitas (VIF)

df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% dplyr::select(Y, X1:X4)
ols_vif  <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
car::vif(ols_vif)

##       X1       X2       X3       X4 
## 2.774826 1.633365 1.553310 1.978360

Uji Autokorelasi Spasial Global (Y)

# Bobot queen
nb <- spdep::poly2nb(as(dat_sf, "Spatial"), queen=TRUE)
lw <- spdep::nb2listw(nb, style="W", zero.policy=TRUE)

# Diagnostik konektivitas
deg <- sapply(nb, length)
c(Mean=mean(deg), Min=min(deg), Max=max(deg))

## Mean  Min  Max 
##  2.5  1.0  5.0

# Moran’s I & Geary’s C
moran.test(dat_sf$Y, lw, zero.policy=TRUE)

## 
##  Moran I test under randomisation
## 
## data:  dat_sf$Y  
## weights: lw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 3.118, p-value = 0.0009103
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.38238167       -0.03125000        0.01759825

geary.test(dat_sf$Y, lw, zero.policy=TRUE)

## 
##  Geary C test under randomisation
## 
## data:  dat_sf$Y 
## weights: lw  
## n reduced by no-neighbour observations 
## 
## Geary C statistic standard deviate = 2.7086, p-value = 0.003379
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.47210729        1.00000000        0.03798472

# Moran scatterplot (nuansa amber)
Y_std  <- scale(dat_sf$Y)[,1]
WY_std <- scale(lag.listw(lw, dat_sf$Y))[,1]
p_moran <- ggplot(data.frame(Y_std, WY_std), aes(x=Y_std, y=WY_std)) +
  geom_point(alpha=.9, size=2, color=accent_fill) +
  geom_smooth(method="lm", se=FALSE, linetype="dashed", color=accent_line) +
  geom_vline(xintercept=0, linetype="dotted") +
  geom_hline(yintercept=0, linetype="dotted") +
  labs(title="Moran Scatterplot (Y distandardisasi)", x="Y (std)", y="W * Y (std)") +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_moran

Analisis Lokal: LISA & Getis–Ord Gi* (Prevalensi)

lisa <- make_lisa_cluster(dat_sf$prev, lw, p_cut=0.05)
dat_sf$Ii      <- lisa$Ii
dat_sf$Pvalue  <- lisa$pvalue
dat_sf$Cluster <- factor(lisa$cluster,
                         levels=c("High-High","Low-Low","High-Low","Low-High","Not Significant"))

# Peta LISA (kuning–amber)
p_lisa <- ggplot(dat_sf) +
  geom_sf(aes(fill=Cluster), color="white", size=0.25) +
  scale_fill_manual(values=col_lisa) +
  labs(title="Peta Klaster LISA — Prevalensi TBC", fill="Klaster") +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_lisa

# Hotspot Getis–Ord Gi*
Gi <- spdep::localG(dat_sf$prev, lw)  # z-score
dat_sf$Gi_star <- as.numeric(Gi)
p_gi <- ggplot(dat_sf) +
  geom_sf(aes(fill = Gi_star), color="white", size=0.25) +
  scale_fill_gradient2(low="#8B5E00", mid="#FFF5CC", high="#FF8A00",
                       midpoint=0, name="Gi* z-score") +
  labs(title="Getis–Ord Gi* — Hotspot Prevalensi TBC") +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_gi

# Ringkas komposisi klaster
table(dat_sf$Cluster)

## 
##       High-High         Low-Low        High-Low        Low-High Not Significant 
##               0               2               0               0              36

Estimasi Model Spasial

# Data non-spasial untuk OLS
df_nogeo <- dat_sf %>% st_set_geometry(NULL) %>% dplyr::select(Y, X1:X4)

# OLS
ols <- lm(Y ~ X1 + X2 + X3 + X4, data=df_nogeo)
summary(ols)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1950.20  -347.31    27.31   265.29  1151.99 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -109.27858  421.77595  -0.259 0.797174    
## X1             0.06495    0.06959   0.933 0.357433    
## X2           149.31095   14.28364  10.453 5.27e-12 ***
## X3           -97.25111   78.66918  -1.236 0.225110    
## X4             8.43853    2.05134   4.114 0.000243 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 579.6 on 33 degrees of freedom
## Multiple R-squared:  0.8912, Adjusted R-squared:  0.878 
## F-statistic: 67.59 on 4 and 33 DF,  p-value: 1.991e-15

# Uji residual OLS terkait spasial
lm.morantest(ols, lw, zero.policy=TRUE)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## weights: lw
## 
## Moran I statistic standard deviate = 1.8304, p-value = 0.0336
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.20709784      -0.06953038       0.02284091

lm.LMtests(ols, lw, test=c("LMlag","LMerr","RLMlag","RLMerr"))

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
## 
## RSlag = 1.8334, df = 1, p-value = 0.1757
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
## 
## RSerr = 2.1706, df = 1, p-value = 0.1407
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
## 
## adjRSlag = 0.70664, df = 1, p-value = 0.4006
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = df_nogeo)
## test weights: listw
## 
## adjRSerr = 1.0438, df = 1, p-value = 0.3069

# SAR (lag Y)
sar <- spatialreg::lagsarlm(Y ~ X1 + X2 + X3 + X4,
                            data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
summary(sar)

## 
## Call:spatialreg::lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, 
##     listw = lw, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##         Min          1Q      Median          3Q         Max 
## -1829.71967  -307.87062    -0.26423   275.30382  1213.13159 
## 
## Type: lag 
## Regions with no neighbours included:
##  12 31 32 35 37 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -592.962223  415.925943 -1.4256    0.1540
## X1             0.115688    0.063689  1.8165    0.0693
## X2           135.774180   13.956756  9.7282 < 2.2e-16
## X3           -28.099166   76.353048 -0.3680    0.7129
## X4             8.730778    1.878674  4.6473 3.363e-06
## 
## Rho: 0.13435, LR test value: 2.5464, p-value: 0.11054
## Asymptotic standard error: 0.07243
##     z-value: 1.8549, p-value: 0.063609
## Wald statistic: 3.4407, p-value: 0.063609
## 
## Log likelihood: -291.735 for lag model
## ML residual variance (sigma squared): 271240, (sigma: 520.81)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 597.47, (AIC for lm: 598.02)
## LM test for residual autocorrelation
## test value: 1.7797, p-value: 0.18218

# SEM (error spasial)
sem <- spatialreg::errorsarlm(Y ~ X1 + X2 + X3 + X4,
                              data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
summary(sem)

## 
## Call:spatialreg::errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, 
##     listw = lw, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1969.664  -317.666    33.163   340.973   899.496 
## 
## Type: error 
## Regions with no neighbours included:
##  12 31 32 35 37 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -381.772182  384.300309 -0.9934   0.32050
## X1             0.126928    0.062753  2.0227   0.04311
## X2           136.189753   13.829124  9.8480 < 2.2e-16
## X3           -78.651500   76.908299 -1.0227   0.30647
## X4            10.607993    1.933421  5.4866 4.096e-08
## 
## Lambda: 0.22545, LR test value: 1.8464, p-value: 0.1742
## Asymptotic standard error: 0.17444
##     z-value: 1.2924, p-value: 0.19622
## Wald statistic: 1.6703, p-value: 0.19622
## 
## Log likelihood: -292.085 for error model
## ML residual variance (sigma squared): 273290, (sigma: 522.77)
## Number of observations: 38 
## Number of parameters estimated: 7 
## AIC: 598.17, (AIC for lm: 598.02)

# SDM (lag Y + lag X)
sdm <- spatialreg::lagsarlm(Y ~ X1 + X2 + X3 + X4,
                            data=dat_sf, listw=lw, type="mixed",
                            method="eigen", zero.policy=TRUE)
summary(sdm)

## 
## Call:spatialreg::lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, 
##     listw = lw, type = "mixed", method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1619.456  -260.774   -18.994   250.872   938.312 
## 
## Type: mixed 
## Regions with no neighbours included:
##  12 31 32 35 37 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -251.694006  458.362615 -0.5491   0.58293
## X1             0.119694    0.063953  1.8716   0.06126
## X2           124.554297   14.710171  8.4672 < 2.2e-16
## X3           -96.409475   94.494533 -1.0203   0.30760
## X4            11.466432    1.961856  5.8447 5.075e-09
## lag.X1        -0.214451    0.138199 -1.5518   0.12072
## lag.X2        12.228546   31.549488  0.3876   0.69831
## lag.X3        46.012344   89.662984  0.5132   0.60783
## lag.X4        -6.269699    2.781101 -2.2544   0.02417
## 
## Rho: 0.31115, LR test value: 4.2059, p-value: 0.040284
## Asymptotic standard error: 0.16
##     z-value: 1.9447, p-value: 0.051816
## Wald statistic: 3.7817, p-value: 0.051816
## 
## Log likelihood: -288.8023 for mixed model
## ML residual variance (sigma squared): 226290, (sigma: 475.7)
## Number of observations: 38 
## Number of parameters estimated: 11 
## AIC: 599.6, (AIC for lm: 601.81)
## LM test for residual autocorrelation
## test value: 0.54449, p-value: 0.46058

# SDEM (error spasial + lag X)
sdem <- spatialreg::errorsarlm(Y ~ X1 + X2 + X3 + X4,
                               data=dat_sf, listw=lw, etype="emixed",
                               method="eigen", zero.policy=TRUE)
summary(sdem)

## 
## Call:spatialreg::errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, 
##     listw = lw, etype = "emixed", method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1476.178  -260.940   -30.626   258.912  1021.531 
## 
## Type: error 
## Regions with no neighbours included:
##  12 31 32 35 37 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -524.873424  490.509301 -1.0701  0.284593
## X1             0.130491    0.066324  1.9675  0.049126
## X2           129.790950   13.716357  9.4625 < 2.2e-16
## X3           -49.676508   92.228435 -0.5386  0.590146
## X4            11.621600    1.783626  6.5157 7.234e-11
## lag.X1        -0.258342    0.147530 -1.7511  0.079925
## lag.X2        66.958801   20.522087  3.2628  0.001103
## lag.X3        -5.254513   94.046763 -0.0559  0.955444
## lag.X4        -3.327177    2.404541 -1.3837  0.166449
## 
## Lambda: 0.42338, LR test value: 6.3805, p-value: 0.011538
## Asymptotic standard error: 0.15079
##     z-value: 2.8077, p-value: 0.004989
## Wald statistic: 7.8834, p-value: 0.004989
## 
## Log likelihood: -287.715 for error model
## ML residual variance (sigma squared): 207290, (sigma: 455.29)
## Number of observations: 38 
## Number of parameters estimated: 11 
## AIC: 597.43, (AIC for lm: 601.81)

# SAC (lag + error)
sac <- spatialreg::sacsarlm(Y ~ X1 + X2 + X3 + X4,
                            data=dat_sf, listw=lw, method="eigen", zero.policy=TRUE)
summary(sac)

## 
## Call:spatialreg::sacsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = dat_sf, 
##     listw = lw, method = "eigen", zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1911.334  -311.503    21.299   255.595  1004.171 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -836.425337  439.421418 -1.9035  0.056979
## X1             0.165630    0.062683  2.6424  0.008233
## X2           128.622146   13.756721  9.3498 < 2.2e-16
## X3           -10.523155   83.616686 -0.1258  0.899851
## X4            10.670614    1.897647  5.6231 1.876e-08
## 
## Rho: 0.11911
## Asymptotic standard error: 0.073634
##     z-value: 1.6176, p-value: 0.10575
## Lambda: 0.20023
## Asymptotic standard error: 0.19044
##     z-value: 1.0514, p-value: 0.29308
## 
## LR test value: 3.9511, p-value: 0.13869
## 
## Log likelihood: -291.0327 for sac model
## ML residual variance (sigma squared): 258320, (sigma: 508.25)
## Number of observations: 38 
## Number of parameters estimated: 8 
## AIC: 598.07, (AIC for lm: 598.02)

Pemilihan Model Terbaik (AIC/BIC)

models <- list(OLS=ols, SAR=sar, SEM=sem, SDM=sdm, SDEM=sdem, SAC=sac)

get_wald_pvalue <- function(model) {
  summ <- summary(model)
  pval <- tryCatch({
    if (!is.null(summ$Wald$p.value)) summ$Wald$p.value
    else if (!is.null(summ$Wald_p.value)) summ$Wald_p.value
    else {
      txt <- capture.output(summary(model))
      line <- grep("Wald statistic", txt, value = TRUE)
      if (length(line) > 0) as.numeric(sub(".*p-value:\\s*([0-9.]+).*", "\\1", line)) else NA
    }
  }, error = function(e) NA)
  if (length(pval) == 0) pval <- NA
  pval
}

cmp <- do.call(rbind, lapply(names(models), function(nm) {
  m <- models[[nm]]
  data.frame(Model=nm, AIC=AIC(m), BIC=BIC(m), LogLik=as.numeric(logLik(m)),
             Wald_p=get_wald_pvalue(m), stringsAsFactors = FALSE)
}))

# Fokus model yang signifikan (jika tersedia), lalu urut AIC naik
cmp_sig <- cmp %>% dplyr::filter(is.na(Wald_p) | Wald_p < 0.05) %>% dplyr::arrange(AIC, BIC)
cmp

##                 Model      AIC      BIC    LogLik      Wald_p
## 1                 OLS 598.0165 607.8420 -293.0083          NA
## Wald statistic    SAR 597.4701 608.9332 -291.7350 0.063609137
## Wald statistic1   SEM 598.1701 609.6332 -292.0850 0.196215825
## Wald statistic2   SDM 599.6046 617.6180 -288.8023 0.051815594
## Wald statistic3  SDEM 597.4300 615.4434 -287.7150 0.004989006
## Wald statistic4   SAC 598.0654 611.1661 -291.0327 0.575395080

cmp_sig

##                 Model      AIC      BIC    LogLik      Wald_p
## Wald statistic3  SDEM 597.4300 615.4434 -287.7150 0.004989006
## 1                 OLS 598.0165 607.8420 -293.0083          NA

# Pilih yang AIC terendah dari cmp_sig bila ada; jika tidak, dari cmp
if (nrow(cmp_sig) > 0) {
  best_name <- cmp_sig$Model[1]
} else {
  best_name <- (cmp %>% dplyr::arrange(AIC, BIC))$Model[1]
}
best <- models[[best_name]]
cat("MODEL TERBAIK:", best_name, "\n")

## MODEL TERBAIK: SDEM

Prediksi, Residual, & Uji Residual (Model Terbaik)

# Vektor prediksi
if (best_name == "OLS") {
  dat_sf$Y_pred <- as.numeric(predict(best, newdata=df_nogeo))
} else {
  dat_sf$Y_pred <- as.numeric(fitted.values(best))
}
dat_sf$Residual <- dat_sf$Y - dat_sf$Y_pred

# Uji Moran pada residual model terbaik
spdep::moran.test(dat_sf$Residual, lw, zero.policy=TRUE)

## 
##  Moran I test under randomisation
## 
## data:  dat_sf$Residual  
## weights: lw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = -0.11792, p-value = 0.5469
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.04891648       -0.03125000        0.02244706

# Peta Prediksi (kuantil, palet kuning)
br_p <- classInt::classIntervals(dat_sf$Y_pred, n=7, style="quantile")$brks
dat_sf$cutPred <- cut(dat_sf$Y_pred, breaks = br_p, include.lowest = TRUE)
p_pred <- ggplot(dat_sf) +
  geom_sf(aes(fill = cutPred), color = "white", size = 0.25) +
  scale_fill_manual(values = gold_pal, name = "Kuantil Prediksi") +
  labs(title = paste0("Peta Prediksi Kasus TBC — Model ", best_name)) +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_pred

# Peta Residual (gradient amber)
p_res <- ggplot(dat_sf) +
  geom_sf(aes(fill = Residual), color = "white", size = 0.25) +
  scale_fill_gradient2(low = "#8B5E00", mid = "#FFF5CC", high = "#FF8A00",
                       midpoint = 0, name = "Residual") +
  labs(title = paste0("Peta Residual — Model ", best_name)) +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))
p_res

# Scatter Prediksi vs Aktual
df_sc <- data.frame(Aktual = dat_sf$Y, Prediksi = dat_sf$Y_pred)
rsq <- cor(df_sc$Aktual, df_sc$Prediksi)^2
rmse_all <- sqrt(mean((df_sc$Aktual - df_sc$Prediksi)^2))
ggplot(df_sc, aes(x = Prediksi, y = Aktual)) +
  geom_point(color = accent_fill, size = 2, alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, color = accent_line, linewidth = 0.8) +
  geom_abline(slope = 1, intercept = 0, color = "black", linetype = "dashed") +
  labs(
    title = sprintf("Prediksi vs Aktual — %s | R²=%.4f, RMSE=%.2f", best_name, rsq, rmse_all),
    x = "Prediksi Kasus TBC", y = "Aktual Kasus TBC"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(color = "#7C2D12", face = "bold"))

Kesimpulan dan Saran

Kesimpulan

Berangkat dari keseluruhan hasil, sebaran kasus TBC di Jawa Timur jelas tidak acak dan cenderung right-skewed: sebagian kecil wilayah menanggung beban jauh di atas rata-rata, dengan median 1.508,5 kasus, rerata 2.036 kasus, dan puncaknya 9.182 kasus di Kota Surabaya. Uji global menunjukkan adanya ketergantungan spasial yang nyata—Moran’s I sebesar 0,382 (p=0,0009) dan Geary’s C 0,472 (p=0,0034) yang menegaskan perlunya pemodelan spasial. Secara lokal, LISA dan Getis–Ord Gi* tidak menemukan banyak klaster kuat (hanya dua wilayah Low–Low dan tanpa High–High yang signifikan), sehingga pola mikro tampak heterogen meski ketergantungan globalnya ada.

OLS pada dasarnya sudah fit dengan R² 0,891 (adj-R² 0,878; F=67,59; p<0,001), tetapi residunya masih autokorelatif (Moran’s I 0,207; p=0,0336), sehingga model spasial pantas dipertimbangkan. Pada akhirnya, Spatial Durbin Error Model menjadi yang paling meyakinkan dengan AIC 597,43, BIC 615,44, log-likelihood –287,72, dan parameter error spasial λ=0,423 (p=0,005). Koefisien yang bermakna memperlihatkan pengaruh positif jumlah penduduk (X1=0,130; p=0,049), kepadatan (X2=129,79; p<0,001), kemiskinan (X4=11,62; p<0,001), serta efek tetangga untuk kepadatan (lag.X2=66,96; p=0,001), sementara jumlah rumah sakit tidak signifikan. Validasi mengindikasikan bahwa residual SDEM sudah bersih dari autokorelasi (Moran’s I –0,049; p=0,547) dan performa prediksi kuat dengan R² 0,923 dan RMSE 455,29. Secara substantif, temuan ini menegaskan peran kepadatan, kemiskinan, dan spillover antarwilayah bertetangga; koridor Surabaya–Sidoarjo–Gresik–Pasuruan muncul sebagai prioritas intervensi.

Saran

Berdasarkan itu, arah tindak lanjutnya menjadi cukup jelas. Intervensi perlu dipusatkan pada kawasan padat dan miskin dengan memperluas skrining aktif, mempercepat inisiasi pengobatan, serta menguatkan edukasi berhenti merokok—khususnya di koridor Surabaya–Sidoarjo–Gresik–Pasuruan tempat risiko terakumulasi. Mengingat adanya spillover, koordinasi lintas kabupaten/kota sebaiknya diperkokoh melalui integrasi data, pelacakan kontak lintas batas, dan penetapan paket kebijakan regional agar rantai penularan tidak lompat-lompat yurisdiksi. Karena jumlah rumah sakit tidak otomatis menekan kasus, fokusnya bukan pada kuantitas fasilitas, melainkan mutu layanan TBC: kesiapan tenaga, ketersediaan obat yang berkesinambungan, serta pencegahan loss to follow-up. Di saat yang sama, kebijakan lintas sektor—perumahan, ventilasi, sanitasi, dan pengendalian konsumsi rokok—perlu menyasar kantong urban-peri-urban yang berisiko tinggi. Untuk pengembangan ke depan, analisis panel multi-tahun dan penambahan variabel seperti kepadatan hunian, mobilitas harian, serta kepatuhan terapi akan membantu membangun sistem peringatan dini yang lebih tajam dan operasional bagi pengendalian TBC berbasis spasial.

Daftar Pustaka

• Anselin, L. (1988). Spatial Econometrics: Methods and Models. Springer.
• Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115.
• Dinas Kesehatan Provinsi Jawa Timur. (2025). Profil Kesehatan Provinsi Jawa Timur 2024.
• Kementerian Kesehatan RI. (2022). Program Penanggulangan Tuberkulosis.
• Rahmawati, N., Karno, F., & Hermanto, E. M. P. (2024). Analisis TBC Jawa Timur Menggunakan GWR. Indonesian Journal of Applied Statistics, 6(2), 116.
• PPTI. (2024). Rokok dan TBC: Dua Ancaman Serius bagi Gen-Z.
• WHO. (2017). Global Tuberculosis Report.

Lampiran

Dashboard terkait: https://fazilaazraanggina.shinyapps.io/DASHBOARD/

Laporan format PDF: https://drive.google.com/file/d/12udeMnHt2ovsJ-8VFvLrIFcFqTSJg5RO/view?usp=sharing

Syntax perhitungan: https://drive.google.com/file/d/1s74FW2VCM0Aospq6yYhtGoCFjkF-5mEL/view?usp=sharing

Tim Pengembang Dashboard:

Nadwah Khairunnisa Fazla Azra Anggina Gusti Ayu Anisa Eldina Stephani Julieta

```