Import Library
library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.4.3
library(spdep)
## Warning: package 'spdep' was built under R version 4.4.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.4.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.4.3
## Loading required package: Matrix
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption
library(GWmodel)
## Warning: package 'GWmodel' was built under R version 4.4.3
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 4.4.3
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(viridis)
## Warning: package 'viridis' was built under R version 4.4.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.4.3
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
library(mapview)
## Warning: package 'mapview' was built under R version 4.4.3
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.4.3
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(geodata)
## Warning: package 'geodata' was built under R version 4.4.3
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.4.3
## terra 1.8.60
##
## Attaching package: 'terra'
## The following object is masked from 'package:tidyr':
##
## extract
library(stringr)
## Warning: package 'stringr' was built under R version 4.4.3
Persiapkan Data
# Ambil peta administrasi Indonesia level 2 (kab/kota)
indo_shp <- geodata::gadm(country = "IDN", level = 2, path = tempdir())
indo_shp
## class : SpatVector
## geometry : polygons
## dimensions : 502, 13 (geometries, attributes)
## extent : 95.00971, 141.0194, -11.00761, 6.076941 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : GID_2 GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1
## type : <chr> <chr> <chr> <chr> <chr> <chr>
## values : IDN.1.2_1 IDN Indonesia IDN.1_1 Aceh NA
## IDN.1.1_1 IDN Indonesia IDN.1_1 Aceh NA
## IDN.1.3_1 IDN Indonesia IDN.1_1 Aceh NA
## NAME_2 VARNAME_2 NL_NAME_2 TYPE_2 (and 3 more)
## <chr> <chr> <chr> <chr>
## Aceh Barat NA NA Kabupaten
## Aceh Barat Daya NA NA Kabupaten
## Aceh Besar NA NA Kabupaten
# Konversi SpatVector -> sf agar kompatibel dengan dplyr dan tmap
indo_sf <- st_as_sf(indo_shp)
# Filter hanya Provinsi Jawa Barat
jabar <- indo_sf %>%
filter(NAME_1 %in% c("West Java", "Jawa Barat")) %>%
mutate(KABKOT = NAME_2)
# Normalisasi huruf besar untuk join data
jabar$KABKOT <- str_to_upper(jabar$KABKOT)
# Import data TBC (gunakan file CSV Anda)
# Format CSV:
# KABKOT;Y;X1;X2;X3;X4
# Y = Jumlah Kasus TBC
# X1 = Jumlah Puskesmas
# X2 = % RT dengan Air Layak
# X3 = % RT dengan Sanitasi Layak
# X4 = Indeks Kualitas Udara
data_tbc = read.csv(file.choose(), sep = ";", dec = ",")
str(data_tbc)
## 'data.frame': 27 obs. of 6 variables:
## $ Kabupaten.Kota: chr "Bogor" "Sukabumi" "Cianjur" "Bandung" ...
## $ Y : int 28527 12251 11924 13678 9404 4795 3299 3995 9414 4451 ...
## $ X1 : int 101 58 47 62 67 40 37 37 60 32 ...
## $ X2 : num 93.9 89.3 91.6 96.3 85.4 ...
## $ X3 : num 88.8 100 70 96.8 91.3 ...
## $ X4 : num 75.7 91.3 89.8 82.6 85.6 ...
data_tbc$KABKOT <- str_to_upper(data_tbc$Kabupaten.Kota)
# Gabungkan data CSV dan shapefile
data_spatial <- jabar %>%
left_join(data_tbc, by = "KABKOT")
# Cek struktur data hasil gabungan
glimpse(data_spatial)
## Rows: 27
## Columns: 21
## $ GID_2 <chr> "IDN.9.2_1", "IDN.9.1_1", "IDN.9.3_1", "IDN.9.4_1", "ID…
## $ GID_0 <chr> "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN",…
## $ COUNTRY <chr> "Indonesia", "Indonesia", "Indonesia", "Indonesia", "In…
## $ GID_1 <chr> "IDN.9_1", "IDN.9_1", "IDN.9_1", "IDN.9_1", "IDN.9_1", …
## $ NAME_1 <chr> "Jawa Barat", "Jawa Barat", "Jawa Barat", "Jawa Barat",…
## $ NL_NAME_1 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ NAME_2 <chr> "Bandung", "Bandung Barat", "Banjar", "Bekasi", "Bogor"…
## $ VARNAME_2 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ NL_NAME_2 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TYPE_2 <chr> "Kabupaten", "Kabupaten", "Kota", "Kabupaten", "Kabupat…
## $ ENGTYPE_2 <chr> "Regency", "Regency", "City", "Regency", "Regency", "Re…
## $ CC_2 <chr> "3204", "3217", "3279", "3216", "3201", "3207", "3203",…
## $ HASC_2 <chr> "ID.JR.BU", "ID.JR.BB", "ID.JR.BA", "ID.JR.BS", "ID.JR.…
## $ KABKOT <chr> "BANDUNG", "BANDUNG BARAT", "BANJAR", "BEKASI", "BOGOR"…
## $ Kabupaten.Kota <chr> "Bandung", "Bandung Barat", "Banjar", "Bekasi", "Bogor"…
## $ Y <int> 13678, 4474, 1435, 14788, 28527, 3299, 11924, 4434, 941…
## $ X1 <int> 62, 32, 10, 51, 101, 37, 47, 13, 60, 38, 67, 49, 50, 80…
## $ X2 <dbl> 96.28, 89.25, 96.97, 96.76, 93.93, 96.42, 91.62, 99.60,…
## $ X3 <dbl> 96.79, 94.91, 99.11, 97.00, 88.77, 95.16, 69.99, 78.58,…
## $ X4 <dbl> 82.63, 73.60, 89.58, 71.25, 75.65, 91.39, 89.83, 75.07,…
## $ geometry <GEOMETRY [°]> POLYGON ((107.6257 -7.29998..., POLYGON ((107.…
#Bersihkan data yang kosong
data_spatial <- data_spatial %>%
filter(!st_is_empty(geometry)) %>%
filter(!is.na(Y) & !is.na(X1) & !is.na(X2) & !is.na(X3) & !is.na(X4))
# Pilih hanya kolom numerik yang ingin dikorelasikan
data_numerik <- data_spatial %>%
st_drop_geometry() %>%
select(Y, X1, X2, X3, X4)
# Hitung matriks korelasi Pearson
cor_matrix <- cor(data_numerik, use = "complete.obs", method = "pearson")
# Tampilkan hasil korelasi
cor_matrix
## Y X1 X2 X3 X4
## Y 1.00000000 0.85316544 -0.02143465 -0.10935991 -0.40969516
## X1 0.85316544 1.00000000 -0.22668926 0.04247369 -0.18870956
## X2 -0.02143465 -0.22668926 1.00000000 0.13228690 -0.16130848
## X3 -0.10935991 0.04247369 0.13228690 1.00000000 -0.01313847
## X4 -0.40969516 -0.18870956 -0.16130848 -0.01313847 1.00000000
# Visualisasi korelasi dengan heatmap
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.4.3
## corrplot 0.95 loaded
corrplot(cor_matrix, method = "color", type = "upper",
tl.col = "black", tl.srt = 45, addCoef.col = "black",
number.cex = 0.7, diag = FALSE)
Sebaran Kasus TBC di Kabupaten/Kota Jawa Barat di Tahun 2024
tm_shape(data_spatial) +
tm_polygons(
fill = "Y",
fill.scale = tm_scale(values = "Reds"),
fill.legend = tm_legend(title = "Jumlah Kasus TBC")
) +
tm_borders() +
tm_title("Sebaran Kasus TBC di Kabupaten/Kota Jawa Barat (2024)") +
tm_layout(legend.outside = TRUE)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"
## Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.
##
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
# Buat matriks bobot spasial Queen contiguity
nb_queen <- poly2nb(data_spatial, queen = TRUE)
lw_queen <- nb2listw(nb_queen, style = "W")
# Cek apakah ada nilai hilang di variabel Y dan X
colSums(is.na(data_spatial[, c("Y", "X1", "X2", "X3", "X4")]))
## Y X1 X2 X3 X4 geometry
## 0 0 0 0 0 0
# Hapus baris yang memiliki NA di variabel Y atau X
data_spatial <- data_spatial %>%
filter(!is.na(Y) & !is.na(X1) & !is.na(X2) & !is.na(X3) & !is.na(X4))
# Pastikan sudah bersih
summary(data_spatial$Y)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1435 4438 5977 8610 12169 28527
sum(is.na(data_spatial$Y))
## [1] 0
# --- Sinkronisasi data dan matriks bobot ---
# 1️⃣ Pastikan tidak ada geometry kosong
data_spatial <- data_spatial %>% filter(!st_is_empty(geometry))
# 2️⃣ Pastikan tidak ada missing value di Y dan variabel lainnya
data_spatial <- data_spatial %>%
filter(!is.na(Y) & !is.na(X1) & !is.na(X2) & !is.na(X3) & !is.na(X4))
# 3️⃣ Buat ulang matriks bobot QUEEN berdasarkan data terbaru
nb_queen <- poly2nb(data_spatial, queen = TRUE)
lw_queen <- nb2listw(nb_queen, style = "W")
# 4️⃣ Cek konsistensi jumlah observasi
length(data_spatial$Y)
## [1] 26
length(nb_queen)
## [1] 26
Uji Autokorelasi Spasial
# Uji Global Moran's I
moran_global <- moran.test(data_spatial$Y, lw_queen)
moran_global
##
## Moran I test under randomisation
##
## data: data_spatial$Y
## weights: lw_queen
##
## Moran I statistic standard deviate = 2.9246, p-value = 0.001725
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.35013827 -0.04000000 0.01779552
# Local Moran (LISA)
local_moran <- localmoran(data_spatial$Y, lw_queen)
data_spatial$LISA <- local_moran[,1]
local_moran
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 -0.032537196 -0.0290959179 0.078694605 -0.01226725 0.9902123958
## 2 -0.150089784 -0.0193758710 0.065182088 -0.51198553 0.6086611399
## 3 1.079067980 -0.0583119640 1.427703651 0.95188949 0.3411530312
## 4 1.707315198 -0.0432364393 0.328638227 3.05362629 0.0022609353
## 5 1.464864472 -0.4493523724 0.569615570 2.53629758 0.0112031489
## 6 0.714778850 -0.0319491953 0.134023258 2.03972856 0.0413773688
## 7 0.333033865 -0.0124416452 0.042150724 1.68273120 0.0924271385
## 8 -0.425103023 -0.0197524724 0.153822814 -1.03352358 0.3013589395
## 9 -0.092339181 -0.0007325053 0.004163072 -1.41977705 0.1556725987
## 10 -0.119458371 -0.0001315980 0.001639276 -2.94721529 0.0032064987
## 11 -0.001643886 -0.0007144005 0.004060250 -0.01458701 0.9883616665
## 12 0.216670086 -0.0094597994 0.053293647 0.97953564 0.3273153892
## 13 0.631661520 -0.0231747577 0.128751852 1.82497109 0.0680053927
## 14 -0.301648711 -0.1099499036 0.777450660 -0.21741174 0.8278874824
## 15 1.175476655 -0.0264821400 0.204814421 2.65588529 0.0079100517
## 16 1.477837089 -0.0077765339 0.200617544 3.31681630 0.0009104946
## 17 -0.103748227 -0.0235109719 0.596913359 -0.10385333 0.9172857461
## 18 -0.530601311 -0.0299948075 0.756473094 -0.57557211 0.5649044344
## 19 0.516834629 -0.0181225396 0.221684993 1.13618974 0.2558771623
## 20 0.377488755 -0.0241238677 0.187027370 0.92865579 0.3530674994
## 21 0.406043190 -0.0195919739 0.065894552 1.65810774 0.0972957219
## 22 -0.480512788 -0.0180682107 0.076880919 -1.66782473 0.0953505251
## 23 0.070511086 -0.0063939482 0.021794545 0.52093180 0.6024142821
## 24 0.621612676 -0.0150179461 0.117517459 1.85710444 0.0632962962
## 25 0.169210798 -0.0267533102 0.089323332 0.65568323 0.5120279599
## 26 0.378870683 -0.0164849084 0.070257010 1.49156798 0.1358124364
## attr(,"call")
## localmoran(x = data_spatial$Y, listw = lw_queen)
## attr(,"class")
## [1] "localmoran" "matrix" "array"
## attr(,"quadr")
## mean median pysal
## 1 High-Low High-Low High-Low
## 2 Low-Low Low-High Low-High
## 3 Low-Low Low-Low Low-Low
## 4 High-High High-High High-High
## 5 High-High High-High High-High
## 6 Low-Low Low-Low Low-Low
## 7 High-High High-High High-High
## 8 Low-High Low-High Low-High
## 9 High-Low High-Low High-Low
## 10 Low-High High-High Low-High
## 11 High-Low High-Low High-Low
## 12 Low-Low Low-Low Low-Low
## 13 High-High High-High High-High
## 14 High-Low High-Low High-Low
## 15 High-High High-High High-High
## 16 High-High High-High High-High
## 17 Low-Low Low-High Low-High
## 18 Low-High Low-High Low-High
## 19 Low-Low Low-Low Low-Low
## 20 Low-Low Low-Low Low-Low
## 21 Low-Low Low-Low Low-Low
## 22 Low-High Low-High Low-High
## 23 Low-Low High-Low Low-Low
## 24 High-High High-High High-High
## 25 Low-Low Low-Low Low-Low
## 26 Low-Low Low-Low Low-Low
#Peta LISA
tm_shape(data_spatial) +
tm_polygons(
fill = "LISA",
fill.scale = tm_scale(values = "RdYlBu"),
fill.legend = tm_legend(title = "Local Moran's I (LISA)")
) +
tm_borders() +
tm_title("Peta LISA Kasus TBC di Jawa Barat") +
tm_layout(legend.outside = TRUE)
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "RdYlBu" is named
## "brewer.rd_yl_bu"
## Multiple palettes called "rd_yl_bu" found: "brewer.rd_yl_bu", "matplotlib.rd_yl_bu". The first one, "brewer.rd_yl_bu", is returned.
##
## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)
Uji Lagrange Multiplier Test (SAR vs SEM)
ols_model <- lm(Y ~ X1 + X2 + X3 + X4, data = data_spatial)
lm_tests <- lm.LMtests(ols_model, lw_queen, test = "all")
## Please update scripts to use lm.RStests in place of lm.LMtests
lm_tests
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial)
## test weights: listw
##
## RSerr = 9.3529, df = 1, p-value = 0.002226
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial)
## test weights: listw
##
## RSlag = 18.325, df = 1, p-value = 1.862e-05
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial)
## test weights: listw
##
## adjRSerr = 1.7062, df = 1, p-value = 0.1915
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial)
## test weights: listw
##
## adjRSlag = 10.679, df = 1, p-value = 0.001084
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial)
## test weights: listw
##
## SARMA = 20.031, df = 2, p-value = 4.469e-05
Estimasi Model SAR
sar_model <- lagsarlm(
Y ~ X1 + X2 + X3 + X4,
data = data_spatial,
listw = lw_queen,
method = "eigen"
)
summary(sar_model)
##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial,
## listw = lw_queen, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2964.834 -1112.827 -77.856 981.781 4548.365
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7111.828 10610.883 -0.6702 0.502705
## X1 247.941 17.651 14.0465 < 2.2e-16
## X2 167.510 93.263 1.7961 0.072478
## X3 -68.125 23.739 -2.8697 0.004109
## X4 -92.440 51.785 -1.7851 0.074253
##
## Rho: 0.29934, LR test value: 15.582, p-value: 7.8985e-05
## Asymptotic standard error: 0.079661
## z-value: 3.7577, p-value: 0.00017147
## Wald statistic: 14.12, p-value: 0.00017147
##
## Log likelihood: -231.5509 for lag model
## ML residual variance (sigma squared): 3108700, (sigma: 1763.2)
## Number of observations: 26
## Number of parameters estimated: 7
## AIC: 477.1, (AIC for lm: 490.68)
## LM test for residual autocorrelation
## test value: 0.43646, p-value: 0.50883
Estimasi Model SEM
sem_model <- errorsarlm(
Y ~ X1 + X2 + X3 + X4,
data = data_spatial,
listw = lw_queen,
method = "eigen"
)
summary(sem_model)
##
## Call:errorsarlm(formula = Y ~ X1 + X2 + X3 + X4, data = data_spatial,
## listw = lw_queen, method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2571.25 -1255.50 -422.94 1247.70 5013.68
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7066.247 11567.157 -0.6109 0.54127
## X1 226.324 16.037 14.1125 < 2e-16
## X2 184.085 101.778 1.8087 0.07050
## X3 -66.077 25.868 -2.5544 0.01064
## X4 -68.433 58.184 -1.1761 0.23954
##
## Lambda: 0.59413, LR test value: 10.648, p-value: 0.001102
## Asymptotic standard error: 0.16372
## z-value: 3.6288, p-value: 0.00028472
## Wald statistic: 13.168, p-value: 0.00028472
##
## Log likelihood: -234.0182 for error model
## ML residual variance (sigma squared): 3454500, (sigma: 1858.6)
## Number of observations: 26
## Number of parameters estimated: 7
## AIC: 482.04, (AIC for lm: 490.68)