Latar Belakang

Analisis regresi adalah satu metode machine learning yang bertujuan untuk memprediksi target variabel berdasarkan variabel-variabel pendukung sesuai dengan kajian-kajian yang mendasarinya. Analisis regresi cukup populer dan banyak diaplikasikan pada banyak disiplin ilmu mulai dari kajian dalam ruang lingkup ekonomi dan bisnis hingga fenomena sehari-hari.

Salah satu kajian dalam analisis regresi adalah dengan mempertimbangkan data spasial. Model regresi spasial dapat melakukan eksplorasi suatu kasus yang berhubungan dengan aspek kewilayahan. Contohnya kasus kelahiran balita di suatu wilayah, kasus tindak pidana di provinsi A, kasus penyebaran penyakit di wilayah yang berdekatan dll. Berbeda dengan pemodelan analisis regresi pada umumnya, pemodelan dengan menggunakan data spasial mempertimbangkan aspek kewilayahan yang memberikan pengaruh lingkungan dan karakteristik yang muncul dari suatu wilayah.

Data spasial adalah data cross-section yang memiliki struktur geografis berupa longitude dan latitude dan dapat divisualisasikan dalam bentuk peta. Berdasarkan Hukum I Tobler, setiap hal memiliki keterkaitan dengan hal lainnya namun yang lebih berdekatan memiliki keterkaitan yang lebih dari lainnya ¹. Hukum ini mendasari pemodelan data spasial yang mempertimbangkan aspek kewilayahan yang berdekatan.

Regresi Terboboti Geografis

Regresi Terboboti Geografis (RTG) adalah metode yang dapat digunakan untuk menganalisis aspek heterogenitas spasial antar variabel. Model regresi terboboti geografis berangkat dari model regresi linear klasik dengan melibatkan hubungan antar variabel berdasarkan aspek ruang. Model yang dihasilkan menggunakan regresi linear klasik merupakan model global yang berarti pendugaan parameter yang dihasilkan berlaku untuk seluruh observasi. Regresi linear klasik dengan Ordinary Least Square (OLS) tidak mempertimbangkan aspek ruang atau efek spasial yang timbul. Hal ini dapat mengakibatkan interpretasi yang keliru terhadap model yang dihasilkan. Simpson’s Paradox adalah salah satu fenomena yang muncul pada kasus statistik yang dapat mengakibatkan kekeliruan pengambilan kesimpulan². Fenomena ini terjadi ketika suatu group data yang dianalisis antara secara terpisah dan secara penggabungan memiliki hasil yang kontradiktif.

Dalam pemodelan kasus regresi linear klasik yang melibatkan ruang, model yang dihasilkan dapat berhadapan dengan Simpson’s Paradox karena bersinggungan dengan data yang melalui agregasi. Setiap observasi dalam pemodelan regresi yang melibatkan ruang memiliki kemungkinan menghasilkan pendugaan parameter yang berbeda ketika pemodelan yang dilakukan bersifat global. Regresi terboboti geografis dapat menangani permasalahan yang muncul akibat misinterpretasi saat melakukan pemodelan regresi linear klasik karena adanya keragaman spasial dalam data³. Output yang muncul pada pemodelan regresi terboboti spasial adalah setiap observasi yang berbasis ruang memiliki model regresi masing-masing.

Model regresi global dapat dinyatakan sebagai berikut:

\(y_{i} = \beta_{0} + \sum_{k=1}^p \beta_{k}x_{ik} + \varepsilon_{i}\)

Melalui regresi terboboti geografis, model regresi tersebut dapat digunakan untuk menduga parameter lokal dengan model umum sebagai berikut:

\(y_{i} = \beta_{0}(u_{i},v_{i}) + \sum_{k=1}^p \beta_{k}(u_{i},v_{i})x_{ik} + \varepsilon_{i}\)

dengan:

\(y_{i}\) : Variabel target untuk setiap lokasi ke-i

\((u_{i},v_{i})\) : Titik koordinat lokasi ke-\(i\) dalam bentuk longitude dan latitude

\(\beta_{k}(u_{i},v_{i})\) : Nilai koefisien regresi untuk variabel prediktor ke-\(k\) dan lokasi ke-\(i\)

\(\varepsilon_{i}\) : Nilai residual regresi antara variabel prediktor dengan variabel target di lokasi ke-\(i\)

Model regresi yang dihasilkan melalui RTG berguna ketika hendak dilakukan eksplorasi suatu variabel beserta pengambilan keputusannya pada daerah tersebut berdasarkan variabel-variabel yang diteliti. Karena setiap daerah memiliki model masing-masing, maka pengambilan kebijakan pada satu daerah akan berbeda dengan daerah lainnya.

Untuk mendeteksi aspek keragaman spasial yang muncul dapat dilakukan menggunakan uji Breusch-Pagan dengan statistik uji sebagai berikut:

\(BP = \frac{1}{2}h^{‘}y(y^{‘}y)^{-1}y^{‘}h\)

dengan

\(h = (\frac{\epsilon_{i}^{2}}{\sigma^{2}} – 1)\)

Vektor observasi variabel target dinyatakan sebagai \(y\) sedangkan \(\epsilon_{i}^{2}\) adalah kuadrat error pengamatan ke-\(i\) dan \(\sigma^{2}\) adalah varians dari \(\epsilon_{i}\)

Setelah dilakukan uji Breush-Pagan dan hasil yang ditunjukan adalah keragaman spasial yang signifikan maka pemodelan dengan regresi terboboti geografis dapat dilakukan. Pada model regresi terboboti geografis setiap lokasi memiliki parameter persamaan yang berbeda sehingga banyaknya model regresi yang diduga adalah sebagai obervasi yang digunakan dalam data.

Pembobot Spasial

Dalam melakukan analisis spasial, dibutuhkan suatu pembobot spasial berdasarkan aspek ketetanggaan atau aspek jarak. Pembentukan pembobot spasial ini bersifat penting karena akan mempengaruhi hasil analisis regresi spasial yang dibentuk. Pembobotan ini akan disimpan ke dalam suatu matriks W yang kemudian bernama matriks pembobot spasial. Adapun beberapa konsep pembobotan spasial yang biasa digunakan diantaranya:

Critical cut of Neighbourhood
Dua lokasi yang memiliki koordinat easting-northing \(s_{i}\) dan \(s_{j}\) dianggap bertetangga jika \(0\le d_{ij} \le d^{*}\) dengan \(d_{ij}\) adalah jarak antarlokasi dan \(d^{*}\) mewakili jarak critical cut off.
\(k^{th}\) Nearest Neighboar (Tetangga terdekat)
Melibatkan lokasi terdekat dengan membertimbangkan nilai \(k\) yang dipilih. Untuk lokasi \(i\) yang terdekat dengan lokasi \(j\) akan dipilih jika lokasi \(i\) masih dalam urutan \(k\) jarak terdekat dalam unit pengamatan.
Contiguity Based Neighbourhood (Ketetanggaan berbasis persinggungan)
Dua lokasi berkoordinat \(s_{i}\) dan \(s_{j}\) disebut bertetangga jika saling berbagi batas wilayah.

Dalam melakukan pembobotan untuk regresi terboboti geografis, akan digunakan suatu fungsi kernal yaitu fungsi pembobot kernel Gaussian dengan persamaan sebagai berikut:

\(w_{ij} = exp \ [\frac{-1}{2}(\frac{d_{ij}}{b})^2]\)

\(d_{ij}\) adalah jarak antara dua lokasi yaitu \(i\) dan \(j\) dan \(b\) adalah lebar jendela (bandwidth). Bandwidth dapat dianalogikan sebagai suatu titik dalam lingkaran yang memiliki efek mempengaruhi area sekitar yang masih berada dalam area tersebut. Hal ini mempengaruhi pendugaan parameternya di lokasi ke-\(i\). Bandwith terbagi menjadi dua yaitu memiliki efek sama pada setiap lokasi pengamatan (fixed) dan berbeda pada setiap lokasi pengamatan (adaptif).

Metode yang dapat digunakan untuk menentukan nilai bandwidth optimum adalah dengan menggunakan validasi silang (Cross Validation). Adapun persamaan untuk mendapatkan CV adalah sebagai berikut:

\(CV = \Sigma_{i=1}^{n}[y_{i} – \hat{y}_{\ne i}(b)]^2\)

\(\hat{y}_{\ne i}(b)\) adalah nilai dugaan \(y_{i}\) pada observasi di lokasi ke-\(i\) dihilangkan dalam pendugaan. Bandwidth optimum akan dicari melalui serangkaian iterasi sehingga didapatkan bandwidth bernilai minimum.

Persiapan Data

Pada tulisan ini akan digunakan metode regresi terboboti geografis pada data tingkat kriminalitas yang di pulau Jawa. Data ini merupakan kumpulan data yang berasal dari publikasi Badan Pusat Statistik (BPS) untuk tahun 2018 dari 6 provinsi yaitu Banten⁴, Jawa Barat⁵, DKI Jakarta⁶ ⁷, Jawa Tengah⁸, DI Yogyakarta⁹ dan Jawa Timur¹⁰. Pemilihan variabel-variabel prediktor yang digunakan berdasarkan penelitian yang dilakukan oleh Putra (2019)¹¹.

library(geojsonio)
library(sp)
library(rgdal)
library(tidyverse)
library(leaflet)
library(RColorBrewer)
library(spatialreg)
library(spdep)
library(htmlwidgets)

Data_Tingkat_Kriminalitas <- read.csv("data_input/rtg/Data Risiko Kejatahan Jawa.csv")

glimpse(Data_Tingkat_Kriminalitas)

#> Rows: 112
#> Columns: 9
#> $ Provinsi           <chr> "Banten", "Banten", "Banten", "Banten", "Banten", "~
#> $ Tipe.Administratif <chr> "Kota", "Kota", "Kabupaten", "Kabupaten", "Kabupate~
#> $ Nama.Daerah        <chr> "Cilegon", "Kota Serang", "Lebak", "Pandeglang", "S~
#> $ TK                 <dbl> 1418.02, 15.34, 18.98, 20.26, 13.65, 9.15, 1110.00,~
#> $ KP                 <dbl> 2458, 2541, 378, 440, 866, 3649, 19757, 19212, 1590~
#> $ PDRB               <dbl> 224, 43, 21, 22, 48, 35, 169, 692, 262, 156, 272, 3~
#> $ PPM                <dbl> 3.25, 5.36, 8.41, 9.61, 4.30, 5.18, 3.39, 3.59, 2.8~
#> $ TPT                <dbl> 9.33, 8.16, 7.69, 8.33, 12.77, 9.70, 5.00, 6.64, 6.~
#> $ IPM                <dbl> 72.65, 71.68, 63.37, 64.34, 65.93, 71.59, 80.88, 81~

Penjelasan untuk setiap variabel dari data yang akan digunakan adalah:

Provinsi : Nama Provinsi
Tipe Administratif: Tipe Kabupaten/Kota
Nama.Daerah : Nama Kabupaten/Kota
TK `: Tingkat Kriminalitas/Risiko Penduduk Terjadi Tindak Pidana per 100.000 Penduduk
KP : Kepadatan Penduduk
PDRB : Produk Domestik Regional Bruto (juta rupiah)
PPM : Persentase Penduduk Miskin
TPT : Tingkat Pengangguran Terbuka
IPM : Indeks Pembangunan Manusia

Membaca Data SHP

Hal pertama yang harus disiapkan adalah menyiapkan file dalam bentuk peta daerah observasi yang hendak dikaji. Adapun peta tersebut dapat didapatkan di GADM. Setelah mendapatkan file tersebut dalam bentuk SHP, dapat dilakukan pemilihan objek observasi sesuai yang dibutuhkan. Dalam melakukan ini kita bisa gunakan Mapshaper ¹².

Pulau_Jawa <- readOGR(dsn = "data_input/Pulau_Jawa", layer = "Pulau_Jawa")

#> OGR data source with driver: ESRI Shapefile
#> Source: "C:\Users\User\Documents\rtg\algoritmablog\content\blog\data_input\Pulau_Jawa", layer: "Pulau_Jawa"
#> with 120 features
#> It has 13 fields

dsn berfungsi sebagai directory tempat penyimpanan file SHP dan layer adalah nama file SHP yang akan digunakan.

Selanjutnya file tersebut dikonversikan ke dalam bentuk sf yang akan digunakan sebagai dasar visualisasi menggunakan leaflet.

jawa_sf <- sf::st_as_sf(Pulau_Jawa)

Menggabungkan Data SHP

Data yang sudah diimport baik dalam bentuk SHP dan CSV perlu digabungkan untuk keperluan data visualisasi yang akan dilakukan pada tahap selanjutnya. Variabel NAME_2 akan digunakan sebagai dasar penggabungan berdasarkan nama kabupaten/kota. Selain itu perlu juga dilakukan penamaan ulang pada nama provinsi sehingga menampilkan informasi yang lebih baik.

data_pulau_jawa <- jawa_sf %>%
  left_join(Data_Tingkat_Kriminalitas, by = c("NAME_2"="Nama.Daerah")) %>%
  mutate(Provinsi = str_replace(Provinsi, "Jakarta Raya", "DKI Jakarta")) %>%
  mutate(Provinsi = str_replace(Provinsi, "Yogyakarta", "DI Yogyakarta"))

Karena peta yang dihasilkan pada setiap visualisasinya bersifat interaktif, ukuran file yang dihasilkan akan besar sehingga untuk efisiensi memori bisa dilakukan pengecilan ukuran file sebagai berikut:

data_pulau_jawa <- data_pulau_jawa %>%
  st_simplify(preserveTopology = T, dTolerance = 0.0001)

Eksplorasi Data

Sebaran Tingkat Kriminalitas Pulau Jawa

Selanjutnya, akan digunakan library leaflet untuk membuat visualisai sebaran tingkat kriminalitas di tiap kabupaten kota di pulau Jawa. Sebelum itu akan diseting elemen-elemen pendukung visualisasi tersebut melingkupi pemilahan warna dan popup teks.

mybin <-  c(min(data_pulau_jawa$TK), (fivenum(data_pulau_jawa$TK)))

mypalette <- colorBin(palette = "Spectral", domain = data_pulau_jawa, na.color = "transparent", bins = mybin)


mytext <- paste(

  data_pulau_jawa$TYPE_2," ", data_pulau_jawa$NAME_2, "<br/>",
  "Provinsi ", data_pulau_jawa$Provinsi, "<br/>",
  "Rasio Kejahatan ", data_pulau_jawa$TK, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_kriminalitas <-
  leaflet(data_pulau_jawa) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette(data_pulau_jawa$TK),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette,
            values = ~data_pulau_jawa$TK,
            title = "Tingkat Kriminalitas",
            labFormat = labelFormat(),
            opacity = 1)

peta_kriminalitas

summary(Data_Tingkat_Kriminalitas)

#>    Provinsi         Tipe.Administratif Nama.Daerah              TK
#>  Length:112         Length:112         Length:112         Min.   :   0.400
#>  Class :character   Class :character   Class :character   1st Qu.:   0.975
#>  Mode  :character   Mode  :character   Mode  :character   Median :  29.500
#>                                                           Mean   : 139.929
#>                                                           3rd Qu.: 119.500
#>                                                           Max.   :1490.000
#>        KP               PDRB             PPM              TPT
#>  Min.   :  278.4   Min.   : 18.21   Min.   : 2.830   Min.   : 1.430
#>  1st Qu.:  778.6   1st Qu.: 24.18   1st Qu.: 6.770   1st Qu.: 3.353
#>  Median : 1098.0   Median : 30.51   Median : 8.550   Median : 4.490
#>  Mean   : 3001.9   Mean   : 60.62   Mean   : 9.585   Mean   : 5.228
#>  3rd Qu.: 2954.2   3rd Qu.: 57.17   3rd Qu.:11.990   3rd Qu.: 6.942
#>  Max.   :19757.0   Max.   :692.00   Max.   :21.210   Max.   :12.770
#>       IPM
#>  Min.   :61.00
#>  1st Qu.:67.93
#>  Median :71.04
#>  Mean   :71.76
#>  3rd Qu.:74.63
#>  Max.   :86.11

Median tingkat kriminalitas di pulau Jawa adalah 29.5 yang bermakna terdapat risiko 29.5 per 100,000 penduduk terjadi tindak pidana dengan visualisasi setiap provinsi sebagai berikut:

library(ggridges)
library(viridis)
library(hrbrthemes)

ridge_chart <- data_pulau_jawa %>%
  filter(TK >= 0) %>%
  mutate(Provinsi =fct_reorder(Provinsi,TK))

ridge_chart %>%
  ggplot(aes(x = TK, y = Provinsi, fill = ..x..)) +
  geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01) +
  scale_fill_viridis(name = "Risiko", option = "C") +
  labs(title = 'Tingkat Kriminalitas Pulau Jawa') +
  theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 3)
    )

Data_Tingkat_Kriminalitas %>%
  filter(!is.na(TK)) %>%
  arrange(TK) %>%
  tail(20) %>%
  mutate(Nama.Daerah=factor(Nama.Daerah, Nama.Daerah)) %>%
  ggplot( aes(x=Nama.Daerah, y=TK) ) +
    geom_segment( aes(x=Nama.Daerah,xend=Nama.Daerah, y=0, yend=TK), color="#3B1877") +
    geom_point(size=3, color="#E69A8D") +
    coord_flip() +
    theme_ipsum() +
    theme(
      panel.grid.minor.y = element_blank(),
      panel.grid.major.y = element_blank(),
      legend.position="none"
    ) +
    xlab("") +
    ylab("Tingkat Kriminalitas 2018 per 100,000 penduduk")

Interpretasi untuk tingkat kriminalitas di pulau Jawa adalah Provinsi DKI Jakarta memiliki tingkat kriminalitas tertinggi dibanding provinsi lain. Pada tingkat kabupaten/kota non DKI Jakarta, kota Cilegon memiliki tingkat kriminalitas 2018 cukup ekstrim dibanding kabupaten/kota non DKI Jakarta. Di Provinsi Jawa Timur, Kota Blitar dan Kota Madiun memiliki tingkat kriminalitas di atas 500 dan merupakan 2 daerah dengan tingkat kriminalitas tertinggi di Jawa Timur.

Hubungan antar variabel

library(GGally)

ggcorr(Data_Tingkat_Kriminalitas, label = TRUE, label_size = 2.9, hjust = 1, layout.exp = 2)

Pada plot korelasi hubungan antar variabel terlihat variabel TK (Tingkat Kriminalitas) memiliki hubungan positif dengan PDRB (Produk Domestik Regional Bruto). Hal ini mengindikasikan semakin tinggi tingkat kriminalitas, PDRB di lokasi tersebut juga tinggi.

Pemodelan Regresi Global

Teknik Multiple Linear Regression

Sebelum dilakukan pemodelan dengan regresi terboboti geografis, akan dilakukan pemodelan dengan teknik multiple linear regression dengan Ordinary Least Square (OLS) untuk mendapatkan model global dari tingkat kriminalitas di pulau Jawa.

library(stats)
library(car)
library(lmtest)

model_global <- lm(TK ~ KP + PDRB + PPM + TPT + IPM,data = Data_Tingkat_Kriminalitas)

summary(model_global)

#>
#> Call:
#> lm(formula = TK ~ KP + PDRB + PPM + TPT + IPM, data = Data_Tingkat_Kriminalitas)
#>
#> Residuals:
#>     Min      1Q  Median      3Q     Max
#> -683.21  -63.76  -23.23   33.30  941.15
#>
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 562.973003 477.430520   1.179   0.2410
#> KP            0.018202   0.007423   2.452   0.0158 *
#> PDRB          2.211645   0.260209   8.499 1.31e-13 ***
#> PPM          -7.092031   6.509865  -1.089   0.2784
#> TPT         -13.044204   9.992780  -1.305   0.1946
#> IPM          -6.627663   5.960117  -1.112   0.2687
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 198 on 106 degrees of freedom
#> Multiple R-squared:  0.5999, Adjusted R-squared:  0.581
#> F-statistic: 31.79 on 5 and 106 DF,  p-value: < 2.2e-16

AIC(model_global)

#> [1] 1510.229

Pendugaan parameter regresi menghasilkan dua variabel yang signifikan terhadap variabel tingkat kriminalitas yaitu Produk Domestik Regional Bruto (PDRB) dan Kepadatan Penduduk (KP). Variabel KP dan PDRB memiliki nilai koefisien positif sehingga tingkat kriminalitas akan meningkat seiring pertambahan kepadatan penduduk dan produk domestik regional bruto.

Adjusted R-Squared yang dihasilkan adalah 58% yang berarti variabel-variabel yang dilibatkan mampu menjelaskan 58% varians tingkat kriminalitas di Pulau Jawa sedangkan 42% sisanya berasal dari faktor lain.

Uji Multikolinearitas

vif(model_global)

#>       KP     PDRB      PPM      TPT      IPM
#> 2.790513 1.509029 2.067562 1.655423 2.897423

Variabel-variabel yang digunakan dalam pemodelan tidak menunjukan multikolinearitas ditandai dengan nilai VIF yang kurang dari 5.

Uji Keragaman Spasial

Untuk mengetahui signifikansi keragaman spasial dalam observasi akan dilakukan uji keragaman spasial menggunakan uji Breusch-Pagan error dari model global yang sudah terbentuk. Uji ini sama dengan uji homoskedastisitas dalam pengujian dengan teknik multiple linear regression.

Hipotesis

\(H_{0}\): Ragam dari error bersifat homogen

\(H_{1}\): Ragam dari error tidak homogen

bptest(model_global,studentize = FALSE)

#>
#>  Breusch-Pagan test
#>
#> data:  model_global
#> BP = 202.11, df = 5, p-value < 2.2e-16

Dengan melakukan uji Breusch-Pagan didapatkan hasil tolak \(H_{0}\) (p-value < 5%) yang berarti ragam dari error yang didapatkan dalam model tidak homogen. Hal ini disebabkan karena terdapat efek spasial dalam data. Dengan munculnya faktor keragaman spasial mengindikasikan setiap daerah memiliki karakteristiknya masing-masing sehingga tidak lagi digunakan model regressi global dengan teknik multiple linear regression. Model yang akan digunakan untuk mengidentifikasi karakteristik masing-masing daerah adalah regresi terboboti geografis.

Pemodelan Regresi Terboboti Geografis

Mencari Titik Koordinat Pusat Setiap Wilayah

Langkah awal pembentukan model regresi terboboti geografis adalah dengan menentukan letak pusat wilayah setiap lokasi observasi. Fungsi tersebut bernama coordinates, hasil penentuan titik pusat tersebut disimpan dalam objek titik_centroid. Sebelumnya, data_pulau_jawa akan dikonversi menjadi format sp serta membuang beberapa lokasi yang tidak memiliki data tingkat kriminalitas.

data_pulau_jawa_sp <- sf:::as_Spatial(data_pulau_jawa %>% drop_na(TK))

titik_centroid <- coordinates(data_pulau_jawa_sp)

Pemodelan Lokal dengan Regresi Terboboti Geografis

Setelah didapatkan titik centroid masing-masing daerah, dilakukan pemodelan regresi terboboti geografis. Library yang digunakan untuk analisis spasial RTG adalah GWmodel

library(GWmodel)

bwd<-bw.gwr(TK ~ KP + PDRB + PPM + TPT + IPM,data = data_pulau_jawa_sp,kernel = "gaussian", adaptive = TRUE)

#> Adaptive bandwidth: 76 CV score: 5714020
#> Adaptive bandwidth: 55 CV score: 5603460
#> Adaptive bandwidth: 40 CV score: 6726455
#> Adaptive bandwidth: 62 CV score: 5521878
#> Adaptive bandwidth: 68 CV score: 5595912
#> Adaptive bandwidth: 59 CV score: 5564769
#> Adaptive bandwidth: 64 CV score: 5502422
#> Adaptive bandwidth: 65 CV score: 5483381
#> Adaptive bandwidth: 66 CV score: 5532648
#> Adaptive bandwidth: 64 CV score: 5502422
#> Adaptive bandwidth: 65 CV score: 5483381

Pembobotan pada fungsi kernel Gaussian tersebut menggunakan bobot kernel adaptif yang berarti terdapat bandwidth yang berbeda pada setiap lokasi pengamatan.

hasil<-gwr.basic(TK ~ KP + PDRB + PPM + TPT + IPM,data = data_pulau_jawa_sp,bw=bwd, kernel = "gaussian", adaptive = TRUE)

hasil

#>    ***********************************************************************
#>    *                       Package   GWmodel                             *
#>    ***********************************************************************
#>    Program starts at: 2021-10-19 14:58:41
#>    Call:
#>    gwr.basic(formula = TK ~ KP + PDRB + PPM + TPT + IPM, data = data_pulau_jawa_sp,
#>     bw = bwd, kernel = "gaussian", adaptive = TRUE)
#>
#>    Dependent (y) variable:  TK
#>    Independent variables:  KP PDRB PPM TPT IPM
#>    Number of data points: 112
#>    ***********************************************************************
#>    *                    Results of Global Regression                     *
#>    ***********************************************************************
#>
#>    Call:
#>     lm(formula = formula, data = data)
#>
#>    Residuals:
#>     Min      1Q  Median      3Q     Max
#> -683.21  -63.76  -23.23   33.30  941.15
#>
#>    Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)
#>    (Intercept) 562.973003 477.430520   1.179   0.2410
#>    KP            0.018202   0.007423   2.452   0.0158 *
#>    PDRB          2.211645   0.260209   8.499 1.31e-13 ***
#>    PPM          -7.092031   6.509865  -1.089   0.2784
#>    TPT         -13.044204   9.992780  -1.305   0.1946
#>    IPM          -6.627663   5.960117  -1.112   0.2687
#>
#>    ---Significance stars
#>    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>    Residual standard error: 198 on 106 degrees of freedom
#>    Multiple R-squared: 0.5999
#>    Adjusted R-squared: 0.581
#>    F-statistic: 31.79 on 5 and 106 DF,  p-value: < 2.2e-16
#>    ***Extra Diagnostic information
#>    Residual sum of squares: 4154935
#>    Sigma(hat): 194.3505
#>    AIC:  1510.229
#>    AICc:  1511.306
#>    BIC:  1450.288
#>    ***********************************************************************
#>    *          Results of Geographically Weighted Regression              *
#>    ***********************************************************************
#>
#>    *********************Model calibration information*********************
#>    Kernel function: gaussian
#>    Adaptive bandwidth: 65 (number of nearest neighbours)
#>    Regression points: the same locations as observations are used.
#>    Distance metric: Euclidean distance metric is used.
#>
#>    ****************Summary of GWR coefficient estimates:******************
#>                    Min.    1st Qu.     Median    3rd Qu.     Max.
#>    Intercept 180.884287 283.972572 357.225181 387.831878 520.5507
#>    KP          0.012381   0.013115   0.015854   0.017046   0.0220
#>    PDRB        0.604867   0.889585   1.395199   2.516883   2.5735
#>    PPM        -8.015244  -6.518996  -5.768897  -4.174312  -3.6043
#>    TPT       -15.921081 -14.944228 -13.313235 -11.725363 -11.4185
#>    IPM        -5.614398  -4.305718  -3.952752  -2.521792  -0.6963
#>    ************************Diagnostic information*************************
#>    Number of data points: 112
#>    Effective number of parameters (2trace(S) - trace(S'S)): 10.52591
#>    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 101.4741
#>    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 1489.56
#>    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 1476.542
#>    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 1397.594
#>    Residual sum of squares: 3219284
#>    R-square value:  0.6900096
#>    Adjusted R-square value:  0.6575343
#>
#>    ***********************************************************************
#>    Program stops at: 2021-10-19 14:58:41

Output hasil tersebut menampilkan perbandingan model antara dua jenis pemodelan yakni pemodelan menggunakan regresi global dan regresi terboboti geografis. Nilai koefiesien yang dihasilkan menggunakan regresi terboboti geografis menunjukan koefisien positif untuk variabel kepadatan penduduk (KP) dan produk domestik regional bruto (PDRB).

Dengan menggunakan dua jenis pemodelan tersebut dapat dilihat perbedaan ukuran kebaikan antar model. Ukuran kebaikan ini dapat diketahui melalui nilai \(R^2\) dan AIC sebagai berikut:

r2_global <- 0.581
aic_global <- AIC(model_global)
r2_gwr <- 0.657
aic_gwr <- 1476.542


komparasi_model <- data.frame(Jenis_Model = c("MKT", "RTG"), Adjusted_R2 = c(r2_global*100, r2_gwr*100), AIC = c(aic_global,aic_gwr))

komparasi_model

#>   Jenis_Model Adjusted_R2      AIC
#> 1         MKT        58.1 1510.229
#> 2         RTG        65.7 1476.542

Model yang dihasilkan dengan regresi terboboti geografis (RTG) menghasilkan model yang lebih baik dibanding model regresi dengan Ordinary Least Square (OLS). Hal ini dapat terlihat dari nilai Adjusted \(R^2\) yang lebih besar dan nilai AIC yang lebih kecil pada model RTG dibandingkan model OLS.

Penduga Parameter, P-Value Observasi dan R Square Lokal

Pemodelan menggunakan RTG menghasilkan sejumlah model sesuai jumlah observasi yang digunakan. Dalam hal ini model RTG menghasilkan 112 model yang merupakan model masing-masing daerah observasi. Setiap observasi atau wilayah memiliki koefisien parameternya masing-masing. Adapun setiap koefisien paramater dari tiap variabel dapat diakses menggunakan fungsi berikut:

parameter <- as.data.frame(hasil$SDF[1:6])

glimpse(parameter)

#> Rows: 112
#> Columns: 6
#> $ Intercept <dbl> 412.9709, 407.4182, 400.1325, 411.2708, 407.2777, 399.4753, ~
#> $ KP        <dbl> 0.01330009, 0.01317198, 0.01300549, 0.01325055, 0.01316647, ~
#> $ PDRB      <dbl> 2.547144, 2.552279, 2.554020, 2.543360, 2.552339, 2.559178, ~
#> $ PPM       <dbl> -7.996424, -8.002756, -8.006775, -8.015244, -8.006302, -7.97~
#> $ TPT       <dbl> -11.73880, -11.72544, -11.75332, -11.79172, -11.72725, -11.6~
#> $ IPM       <dbl> -4.543578, -4.467307, -4.361653, -4.508615, -4.464381, -4.36~

Setiap observasi memiliki nilai signifikansi masing-masing untuk setiap variabel. Nilai-nilai ini dapat diakses melalui fungsi berikut:

p.value.GWR <- gwr.t.adjust(hasil)$results$p

glimpse(p.value.GWR)

#>  num [1:112, 1:6] 0.383 0.391 0.401 0.385 0.391 0.402 0.408 0.411 0.414 0.416 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:112] "1" "2" "3" "4" ...
#>   ..$ : chr [1:6] "Intercept_p" "KP_p" "PDRB_p" "PPM_p" ...

Sebagaimana model regresi pada umumnya, setiap model regresi pada tiap observasi memiliki \(R^2\) masing-masing. Nilai-nilai ini dapat pula diakses melalui fungsi berikut:

local.R2<-hasil$SDF$Local_R2

glimpse(local.R2)

#>  num [1:112] 0.74 0.74 0.739 0.74 0.74 ...

Eksplorasi Koefisien Regresi Terboboti Geografis

Setelah didapatkan koefisien regresi untuk setiap variabel pada tiap observasi, dapat dilakukan visualisasi untuk melihat pembobotan koefisien

nama.daerah <- data_pulau_jawa_sp$NAME_2

hasil.GWR <- as.data.frame(cbind(nama.daerah,parameter,p.value.GWR,local.R2))

glimpse(hasil.GWR)

#> Rows: 112
#> Columns: 14
#> $ nama.daerah <chr> "Cilegon", "Kota Serang", "Lebak", "Pandeglang", "Serang",~
#> $ Intercept   <dbl> 412.9709, 407.4182, 400.1325, 411.2708, 407.2777, 399.4753~
#> $ KP          <dbl> 0.01330009, 0.01317198, 0.01300549, 0.01325055, 0.01316647~
#> $ PDRB        <dbl> 2.547144, 2.552279, 2.554020, 2.543360, 2.552339, 2.559178~
#> $ PPM         <dbl> -7.996424, -8.002756, -8.006775, -8.015244, -8.006302, -7.~
#> $ TPT         <dbl> -11.73880, -11.72544, -11.75332, -11.79172, -11.72725, -11~
#> $ IPM         <dbl> -4.543578, -4.467307, -4.361653, -4.508615, -4.464381, -4.~
#> $ Intercept_p <dbl> 0.383, 0.391, 0.401, 0.385, 0.391, 0.402, 0.408, 0.411, 0.~
#> $ KP_p        <dbl> 0.083, 0.086, 0.090, 0.084, 0.086, 0.090, 0.092, 0.093, 0.~
#> $ PDRB_p      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
#> $ PPM_p       <dbl> 0.223, 0.224, 0.226, 0.222, 0.224, 0.228, 0.232, 0.234, 0.~
#> $ TPT_p       <dbl> 0.240, 0.242, 0.242, 0.238, 0.242, 0.246, 0.248, 0.249, 0.~
#> $ IPM_p       <dbl> 0.439, 0.448, 0.460, 0.442, 0.448, 0.460, 0.465, 0.468, 0.~
#> $ local.R2    <dbl> 0.7403801, 0.7398869, 0.7391092, 0.7402497, 0.7399046, 0.7~

Data berisi paramater-parameter tersebut harus digabung dengan data bertipe sf untuk keperluan visualisasi.

gwr_gabung <- left_join(jawa_sf,hasil.GWR,by=c("NAME_2"="nama.daerah")) %>%
  mutate(Provinsi = str_replace(NAME_1, "Jakarta Raya", "DKI Jakarta")) %>%
  mutate(Provinsi = str_replace(NAME_1, "Yogyakarta", "DI Yogyakarta"))

colSums(is.na(gwr_gabung))

#>       GID_0      NAME_0       GID_1      NAME_1   NL_NAME_1       GID_2
#>           0           0           0           0         120           0
#>      NAME_2   VARNAME_2   NL_NAME_2      TYPE_2   ENGTYPE_2        CC_2
#>           0         120         120           0           0           0
#>      HASC_2   Intercept          KP        PDRB         PPM         TPT
#>           2           8           8           8           8           8
#>         IPM Intercept_p        KP_p      PDRB_p       PPM_p       TPT_p
#>           8           8           8           8           8           8
#>       IPM_p    local.R2    geometry    Provinsi
#>           8           8           0           0

Pada data frame tersebut terdapat beberapa data NA yang harus dihilangkan. Selain itu hanya akan dipilih variabel-variabel untuk mengurangi redudansi informasi.

gwr_clean <- gwr_gabung %>%
  select(c(Provinsi, NAME_2,Intercept, KP, PDRB, PPM, TPT, IPM, Intercept_p, KP_p, PDRB_p, PPM_p, TPT_p, IPM_p, local.R2, geometry)) %>%
  drop_na(Intercept)


gwr_clean <- gwr_clean %>%
  st_simplify(preserveTopology = T, dTolerance = 0.001)


colSums(is.na(gwr_clean))

#>    Provinsi      NAME_2   Intercept          KP        PDRB         PPM
#>           0           0           0           0           0           0
#>         TPT         IPM Intercept_p        KP_p      PDRB_p       PPM_p
#>           0           0           0           0           0           0
#>       TPT_p       IPM_p    local.R2    geometry
#>           0           0           0           0

Pemetaan Koefisien Kepadatan Penduduk (KP)

Koefisien kepadatan penduduk merupakan variabel yang signifikan menggunakan regresi global. Dengan model regresi terboboti geografis dapat dihasilkan koefisien masing-masing wilayah dengan visualisasi sebagai berikut:

mybin_kp <-  (fivenum(gwr_clean$KP))

mypalette_kp <- colorBin(palette = "YlOrBr", domain = gwr_clean, na.color = "transparent", bins = mybin_kp)


mytext_kp <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Koefisien Kepadatan Penduduk ", gwr_clean$KP, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_koef_kp <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_kp(gwr_clean$KP),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_kp,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette_kp,
            values = ~ gwr_clean$KP,
            title = "Koefisien Kepadatan Penduduk",
            labFormat = labelFormat(),
            opacity = 1)


peta_koef_kp

Melalui visualisasi tersebut dapat terlihat koefisien kepadatan penduduk di provinsi Jawa Timur memiliki bobot yang lebih tinggi dibanding provinsi lainnya.

Pemetaan Koefisien Produk Domestik Regional Bruto (PDRB)

mybin_pdrb <-  (fivenum(gwr_clean$PDRB))

mypalette_pdrb <- colorBin(palette = "PuRd", domain = gwr_clean, na.color = "transparent", bins = mybin_pdrb)


mytext_pdrb <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Koefisien PDRB ", gwr_clean$PDRB, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_koef_pdrb <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_pdrb(gwr_clean$PDRB),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_pdrb,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette_pdrb,
            values = ~ gwr_clean$PDRB,
            title = "Koefisien Kepadatan Penduduk",
            labFormat = labelFormat(),
            opacity = 1)

peta_koef_pdrb

Pada visualisasi menggunakan koefisien Produk Domestik Regional Bruto (PDRB) didapat hasil yang berbeda dengan koefisien kepadatan penduduk. Pada visualisasi ini, bobot koefisien PDRB yang dihasilkan lebih besar pada pulau Jawa bagian barat meliputi DKI Jakarta, Banten, dan Jawa Barat diikuti Jawa Tengah, DI Yogyakarta dan sebagian wilayah Jawa Timur.

Eksplorasi Nilai Signifikansi Variabel

Pemetaan Signifikansi Variabel Kepadatan Penduduk (KP)

Setelah dilakukan pemetaan koefisien kepadatan penduduk untuk seluruh wilayah, akan dilakukan pemetaan di daerah mana saja variabel kepadatan penduduk bernilai signifikan (p-value < 5%).

gwr_clean <- gwr_clean %>% mutate(sig_kp_5 = ifelse(KP_p <= 0.05,"Signifikan", "Tidak Signifikan")) %>%
  mutate(sig_kp_10 = ifelse(KP_p <= 0.10, "Signifikan", "Tidak Signifikan"))

mybin_kp_sig <- c(min(gwr_clean$KP_p), 0.05, max(gwr_clean$KP_p))

mypalette_kp_sig <- colorBin(palette = "GnBu", domain = gwr_clean, na.color = "transparent", bins = mybin_kp_sig)


mytext_kp_sig <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Signifikansi Variabel 5% : "," ", gwr_clean$sig_kp_5, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_koef_kp_sig <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_kp_sig(gwr_clean$KP_p),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_kp_sig,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette_kp_sig,
            values = ~ gwr_clean$KP_p,
            title = "Signifikansi Variabel",
            labFormat = labelFormat(),
            opacity = 1)


peta_koef_kp_sig

Terlihat tidak semua daerah memiliki p-value yang signifikan untuk variabel Kepadatan Penduduk. Daerah-daerah yang memiliki nilai signifikan berada pada provinsi Jawa Tengah yaitu Kabupaten Pemalang, Kabupaten Banjarnegara dan Kabupaten Kendal. Sementara di Jawa Timur, nilai P-value yang signifikan beberapa di beberapa wilayah diantaranya Kabupaten Jember, Kabupaten Banyuwangi, dan Kabupaten Pasuruan.

Pemetaan Signifikansi Variabel Produk Domestik Regional Bruto

Selain variabel Kepadatan Penduduk, pada regresi global didapatkan hasil variabel Produk Domestik Regional Bruto juga bernilai signifikan. Berikut adalah visualisasi untuk variabel PDRB.

gwr_clean <- gwr_clean %>% mutate(sig_pdrb = ifelse(PDRB_p <= 0.05,"Signifikan", "Tidak Signifikan"))


mybin_pdrb_sig <- c(min(gwr_clean$PDRB_p), 0.05, max(gwr_clean$PDRB_p))

mypalette_pdrb_sig <- colorBin(palette = "YlGn", domain = gwr_clean, na.color = "transparent", bins = mybin_pdrb_sig)


mytext_pdrb_sig <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Signifikansi Variabel 5% : "," ", gwr_clean$sig_pdrb, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_pdrb_sig <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_pdrb_sig(gwr_clean$PDRB_p),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_pdrb_sig,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette_pdrb_sig,
            values = ~ gwr_clean$PDRB_p,
            title = "Signifikansi Variabel",
            labFormat = labelFormat(),
            opacity = 1)

peta_pdrb_sig

Hampir semua observasi pada variabel PDRB bernilai signifikan pada taraf nyata 5%. Namun variabel ini tidak bernilai signifikan untuk sebagian wilayah Jawa Timur seperti Kabupaten Madiun, Kabupaten Bojonegoro dan Kabupaten Tuban.

Eksplorasi R2 Lokal

Setiap observasi dalam model RTG memiliki R2 masing-masing, hasil ini dapat pula dipetakan sebagai berikut:

gwr_clean <- gwr_clean %>% mutate(local.R2 = local.R2 * 100)

mybin_r2 <- (fivenum(gwr_clean$local.R2))

mypalette_r2 <- colorBin(palette = "YlOrBr", domain = gwr_clean, na.color = "transparent", bins = mybin_r2)


mytext_r2 <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Nilai R2 : "," ", round(gwr_clean$local.R2,2),"%", "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_r2_lokal <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_r2(gwr_clean$local.R2),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_r2,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            pal = mypalette_r2,
            values = ~ gwr_clean$local.R2,
            title = "Nilai R2",
            labFormat = labelFormat(),
            opacity = 1)

peta_r2_lokal

Nilai \(R^2\) yang lebih tinggi berada di pulau Jawa bagian barat meliputi DKI Jakarta, Banten dan Jawa Barat diikuti provinsi Jawa Tengah, DI Yogyakarta, Jawa Timur bagian barat dan terakhir Jawa Timur bagian Timur.

Sebagai perbandingan, nilai \(R^2\) di Kabupaten Bogor sebesar 73.71% yang bermakna variabel-variabel prediktor mampu menjelaskan keragaman tingkat kriminalitas sebesar 73.71% sedangkan sisanya dijelaskan diluar variabel prediktor. Hal ini berbeda dengan Kabupaten Banyuwangi yang memiliki \(R^2\) lebih kecil yaitu 63.01%.

Pendekatan Analisis Klaster

Untuk meringkas pendugaan parameter masing-masing yang dihasilkan oleh RTG, akan dilakukan pendekatan analisis klaster. Analisis klaster akan menghasilkan rangkuman klaster-klaster berdasarkan karakteristik setiap observasi yang mirip. Pengklasteran akan dilakukan pada setiap penduga parameter pada setiap variabel prediktor.

Kemungkinan Pengklasteran

library(ggplot2)

ggplot(gwr_clean, aes(KP, PDRB, color = Provinsi)) +
    geom_point(alpha = 0.8) + theme_minimal()

Visualiasi tersebut memperlihatkan sejumlah observasi memiliki kemungkinan membentuk suatu klaster berdasarkan koefisien parameter Kepadatan Penduduk dan Produk Domestik Regional Bruto.

Penentuan Klaster Optimum Berdasarkan Pendugaan Parameter

Untuk menentukan jumlah klaster optimum, akan dilakukan dua metode yaitu “wss” dan “silhouette” menggunakan metode elbow. Nilai optimum yang muncul akan digunakan sebagai “k” pada analisis klaster. Sebelumnya variabel-variabel prediktor akan dilakukan subsetting.

gwr_clean_sp <- as.data.frame(gwr_clean)

gwr_clean_koef <- gwr_clean_sp %>%
   select(KP, PDRB, PPM, TPT, IPM)

library(factoextra)

fviz_nbclust(gwr_clean_koef, kmeans, method = "wss")

fviz_nbclust(gwr_clean_koef, kmeans, method = "silhouette")

Berdasarkan metode elbow diketahui klaster optimum yang muncul adalah 4 karena jumlah kuadarat setelah klaster ke 4 tidak lagi mengalami penurunan drastis.

set.seed(11)
gwr_klaster <- kmeans(gwr_clean_koef, 4)

gwr_klaster

#> K-means clustering with 4 clusters of sizes 34, 21, 29, 28
#>
#> Cluster means:
#>           KP      PDRB       PPM       TPT       IPM
#> 1 0.01280933 2.5446312 -7.379492 -11.64877 -4.178523
#> 2 0.02028626 0.9893856 -5.930103 -15.72078 -4.009320
#> 3 0.01614251 1.7541582 -4.110653 -12.74133 -3.832290
#> 4 0.01587982 0.8119999 -4.903384 -14.45671 -1.501139
#>
#> Clustering vector:
#>   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
#>   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
#>  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40
#>   1   1   1   1   1   1   1   1   1   1   1   1   1   1   3   3   3   4   3   3
#>  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60
#>   3   3   4   3   4   3   3   4   3   3   3   3   3   3   4   3   3   3   3   4
#>  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80
#>   3   3   4   4   4   3   3   4   3   2   2   2   4   4   2   2   2   4   4   4
#>  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100
#>   4   4   2   2   2   2   4   2   4   4   2   2   4   4   4   2   2   4   2   2
#> 101 102 103 104 105 106 107 108 109 110 111 112
#>   2   2   2   2   4   4   4   3   4   3   3   3
#>
#> Within cluster sum of squares by cluster:
#> [1] 20.19429 24.72763 36.68086 31.34728
#>  (between_SS / total_SS =  85.0 %)
#>
#> Available components:
#>
#> [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#> [6] "betweenss"    "size"         "iter"         "ifault"

gwr_klaster$cluster <- as.factor(gwr_klaster$cluster)

fviz_cluster(object = gwr_klaster,
             data = gwr_clean_koef)

gwr_clean_koef$cluster <- gwr_klaster$cluster

gwr_clean_koef %>%
   group_by(cluster) %>%
   summarise_all(.funs = "mean")

#> # A tibble: 4 x 6
#>   cluster     KP  PDRB   PPM   TPT   IPM
#>   <fct>    <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1       0.0128 2.54  -7.38 -11.6 -4.18
#> 2 2       0.0203 0.989 -5.93 -15.7 -4.01
#> 3 3       0.0161 1.75  -4.11 -12.7 -3.83
#> 4 4       0.0159 0.812 -4.90 -14.5 -1.50

Terdapat 4 klaster yang terbentuk dengan karakteristik masing-masing. Klaster pertama menunjukan koefisien PDRB terbesar dibanding klaster lain sementara pada klaster kedua memiliki karakteristik koefisien KP yang lebih besar dibanding klaster lainnya.

gwr_clean$clusterp <- gwr_klaster$cluster

gwr_clean <- gwr_clean %>% mutate(klaster_p = as.numeric(clusterp))

Pemetaan Analisis Klaster Berdasarkan Koefisien Variabel Prediktor

mybin_k <- c(1,2,3,4,5)

mypalette_k <- colorBin(palette = c("Yellow", "Green", "Turquoise", "Navajowhite"), domain = gwr_clean, na.color = "transparent", bins = mybin_k)


mytext_k <- paste(

  gwr_clean$TYPE_2," ", gwr_clean$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Anggota Klaster : "," ", gwr_clean$klaster_p, "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_klaster_p <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_k(gwr_clean$klaster_p),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.7,
    label = mytext_k,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            colors = c("Yellow", "Green", "Turquoise", "Navajowhite"),
            labels= c("Klaster 1", "Klaster 2", "Klaster 3", "Klaster 4"),
            values = ~ gwr_clean$klaster_p,
            title = "Anggota Klaster",
            labFormat = labelFormat(),
            opacity = 1)


peta_klaster_p

Peta tersebut memperlihatkan pemetaan berdasarkan anggota klaster masing-masing. Terlihat klaster pertama berasal dari provinsi Jawa Barat, DKI Jakarta, dan Banten. Klaster kedua terdiri dari sebagian besar kabupaten/kota di provinsi Jawa Tengah dan DI Yogyakarta. Sementara di Provinsi Jawa Timur terdapat dua klaster dengan karakteristik berbeda.

Penentuan Klaster Optimum Berdasarkan Signifikansi Variabel

Selain menggunakan koefisien setiap variabel prediktor, analisis dapat juga dilakukan menggunakan signifikansi variabel prediktor. Seperti halnya pada pengklasteran menggunakan koefisien antar variabel, akan juga dilakukan penentuan klaster optimum berdasarkan signifikansi antar variabel prediktor.

gwr_clean_p <- as.data.frame(gwr_clean)

gwr_clean_p <- gwr_clean_p %>%
   select(KP_p, PDRB_p, PPM_p, TPT_p, IPM_p)

fviz_nbclust(gwr_clean_p, kmeans, method = "wss")

fviz_nbclust(gwr_clean_p, kmeans, method = "silhouette")

Berbeda dengan jenis pengklasteran sebelumnya, pada pengklasteran menggunakan signifikansi variabel prediktor terdapat 6 klaster yang optimum.

set.seed(11)
gwr_klaster_p <- kmeans(gwr_clean_p, 6)

gwr_klaster_p

#> K-means clustering with 6 clusters of sizes 29, 18, 16, 11, 16, 22
#>
#> Cluster means:
#>         KP_p      PDRB_p     PPM_p     TPT_p     IPM_p
#> 1 0.09806897 0.000000000 0.2563448 0.2523103 0.4846552
#> 2 0.05183333 0.000000000 0.5122778 0.2348333 0.4705000
#> 3 0.07256250 0.031125000 0.3865000 0.2110625 0.5890000
#> 4 0.02345455 0.005454545 0.3606364 0.1513636 0.4185455
#> 5 0.06668750 0.003812500 0.5491875 0.2196250 0.6509375
#> 6 0.12445455 0.076590909 0.4500909 0.2216818 0.8364545
#>
#> Clustering vector:
#>   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
#>   1   1   1   1   1   1   1   1   1   1   1   1   1   3   1   3   1   1   3   1
#>  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40
#>   1   1   1   1   3   1   1   3   1   1   1   1   1   1   2   2   2   6   5   2
#>  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60
#>   2   5   5   2   6   2   2   5   5   2   2   2   5   5   5   2   2   2   2   5
#>  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80
#>   5   5   6   6   6   2   2   6   2   4   4   3   6   6   4   3   4   6   6   6
#>  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100
#>   6   6   3   3   3   4   3   4   6   6   3   3   6   6   6   4   3   6   4   4
#> 101 102 103 104 105 106 107 108 109 110 111 112
#>   3   4   4   3   6   6   6   5   5   5   5   5
#>
#> Within cluster sum of squares by cluster:
#> [1] 0.06387917 0.05586911 0.12999663 0.03678127 0.05745300 0.14756682
#>  (between_SS / total_SS =  88.6 %)
#>
#> Available components:
#>
#> [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#> [6] "betweenss"    "size"         "iter"         "ifault"

gwr_klaster_p$klaster_s <- as.numeric(gwr_klaster_p$cluster)

fviz_cluster(object = gwr_klaster_p,
             data = gwr_clean_p)

gwr_clean_p$klaster_s <- gwr_klaster_p$cluster

gwr_clean_p %>%
   group_by(klaster_s) %>%
   summarise_all(.funs = "mean")

#> # A tibble: 6 x 6
#>   klaster_s   KP_p  PDRB_p PPM_p TPT_p IPM_p
#>       <int>  <dbl>   <dbl> <dbl> <dbl> <dbl>
#> 1         1 0.0981 0       0.256 0.252 0.485
#> 2         2 0.0518 0       0.512 0.235 0.471
#> 3         3 0.0726 0.0311  0.386 0.211 0.589
#> 4         4 0.0235 0.00545 0.361 0.151 0.419
#> 5         5 0.0667 0.00381 0.549 0.220 0.651
#> 6         6 0.124  0.0766  0.450 0.222 0.836

Hampir semua klaster yang muncul memperlihatkan nilai signifikan untuk p-value <= 0.05 pada variabel PDRB. Hal ini tidak terjadi pada klaster keenam. Selain itu, wilayah-wilayah yang masuk di klaster keenam juga tidak memperlihatkan nilai p-value yang signifikan untuk variabel KP. Klaster keempat adalah klaster dengan p-value signifikan baik pada variabel KP maupun PDRB.

Pemetaan Analisis Klaster Berdasarkan Signifikansi Variabel prediktor

gwr_clean$cluster_s <- gwr_klaster_p$cluster

mybin_p <- c(1,2,3,4,5,6,7)

mypalette_p <- colorBin(palette = c("yellow", "green", "Turquoise", "blue", "red","brown"), domain = gwr_clean_p, na.color = "transparent", bins = mybin_p)


mytext_p <- paste(

  gwr_clean$NAME_2," ", gwr_clean_p$NAME_2, "<br/>",
  "Provinsi ", gwr_clean$Provinsi, "<br/>",
  "Anggota Klaster : "," ", gwr_clean_p$klaster_s
  , "<br/>",
  sep = "") %>%
  lapply(htmltools::HTML)

peta_klaster_s <- leaflet(gwr_clean) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor = ~ mypalette_p(gwr_clean$cluster_s),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "2",
    fillOpacity = 0.5,
    label = mytext_p,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend("bottomleft",
            colors =c("Yellow", "Green", "Turquoise", "Blue", "Red","Brown"),
            values = ~ gwr_clean$cluster_s,
            labels = c("Klaster 1", "Klaster 2", "Klaster 3", "Klaster 4","Klaster 5", "Klaster 6"),
            title = "Anggota Klaster",
            labFormat = labelFormat(),
            opacity = 1)

peta_klaster_s

Pada pemetaan tersebut muncul enam klaster dengan karakteristik masing-masing. Pada klaster pertama didominasi oleh kabupaten/kota dari provinsi Jawa Barat, DKI Jakarta dan Banten. Klaster ketiga merupakan kabupaten/kota yang terletak antara Jawa Barat dan Jawa Tengah. Klaster keempat yang memiliki nilai p-value signifikan pada variabel KP dan PDRB berasal dari provinsi Jawa Timur bagian Timur.

Kesimpulan

Model regresi yang dihasilkan dengan regresi terboboti geografis dapat dijadikan sebagai bahan eksplorasi pada suatu daerah yang memiliki karakteristik berbeda satu sama lain. Berbeda dengan metode multiple linear regression yang bersifat global dan menghasilkan satu persamaan model, regresi terboboti geografis mampu menjelaskan karakteristik suatu daerah yang beragam.

Regresi terboboti geografis baik digunakan ketika ingin melakukan eksplorasi serangkaian variabel terhadap suatu daerah dan bagaimana antar daerah memiliki variabel prediktor signifikan yang tidak sama. Hasil eksplorasi tersebut dapat dijadikan pengambilan kebijakan di lokasi mana suatu variabel prediktor harus diperhatikan. Pada pemodelan regresi terboboti geografis yang telah dilakukan terdapat variabel Kepadatan Penduduk (KP) dan Produk Domestik Regional Bruto (PDRB) memiliki koefisien positif dan memiliki p-value signifikan pada taraf nyata 5%. Variabel Kepadatan Penduduk dan Produk Domestik Regional Bruto yang signifikan tidak berlaku di semua daerah. Terdapat beberapa daerah yang memiliki variabbel PDRB yang signifikan namun tidak signifikan di variabel Kepadatan Penduduk. Nilai R2 yang dihasilkan di pulau Jawa bagian barat lebih tinggi dibandingkan Jawa bagian timur. Hal ini menunjukan variabel-variabel prediktor yang digunakan mampu menjelaskan variabel tingkat kriminalitas lebih baik pada daerah di barat pulau Jawa.

Hasil regresi terboboti geografis dapat didekati dengan analisis klaster untuk merangkum karakteristik daerah yang mirip berdasarkan koefisien parameter yang dihasilkan. Selain itu dengan menggunakan pendekatan analisis klaster dapat dilakukan eksplorasi daerah observasi yang memiliki karakteristik yang mirip pada tingkat signifikansi variabel. Berdasarkan parameter model regresi yang dihasilkan pulau Jawa terbagi menjadi 4 klaster. Pada tingkat signifkansi variabel prediktor, terdapat 6 klaster yang membagi pulau Jawa berdasarkan antar daerahnya. Setiap daerah yang berdekatan cenderung berada dalam satu klaster. Daerah-daerah yang berada di perlintasan antar provinsi cenderung membentuk suatu klaster tertentu.

Penelitian Lebih Lanjut

Terdapat beberapa macam fungsi kernel lain yang bisa digunakan untuk melakukan analisis RTG. Beberapa diantaranya adalah kernel Tricube, kernel kuadrat, dan kernel eksponensial. Penggunaan pembobotan kernel yang berbeda dapat dikaji lebih lanjut untuk melihat kemampuan pembobot dalam menghasilkan pemodelan yang lebih baik. Selain itu terdapat beberapa macam model lain dalam pemodelan RTG yaitu RTG Logistik, RTG Lasso dan RTG Ridge. Setiap teknik pemodelan RTG ini memiliki karakter masing-masing dan dapat dilakukan sebagai rujukan perbandingan teknik terbaik yang dapat diterapkan pada suatu data.

Sumber Rujukan

https://osgis.org/2020/06/what-is-toblers-first-law-of-geography/↩︎
https://anatomisebiostats.com/biostatistics-blog/simpsons-paradox-hidden-bias/↩︎
Fotheringham AS, Brunsdon C, Charlton M. 2002. Geographically Weighted Regression : the Analysis of Spatially Varying Relationships. West Sussex (UK): John Wiley and Sons LTD.↩︎
[BPS]. Badan Pusat Statistik.2020.Provinsi Banten Dalam Angka 2020. Serang(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2021.Provinsi Jawa Barat Dalam Angka 2020. Bandung(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2020.Provinsi DKI Jakarta Dalam Angka 2020. Jakarta(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2018.Statistik Kesejahteraan Rakyat Provinsi DKI Jakarta 2018. Jakarta(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2021.Provinsi Jawa Timur Dalam Angka 2021. Surabaya(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2021.Provinsi DI Yogyakarta Dalam Angka 2021. Yogyakarta(ID) : BPS ↩︎
[BPS]. Badan Pusat Statistik.2021.Provinsi Jawa Tengah Dalam Angka 2020. Semarang(ID) : BPS ↩︎
Zikalta P. 2019. Model Regresi Terboboti Geografis Temporal Kekar untuk Tingkat Kriminalitas di Provinsi Jawa Tengah dan Provinsi Jawa Timur. [tesis]. Bogor (ID): Institut Pertanian Bogor.↩︎
https://algotech.netlify.app/blog/creating-choropleth-with-mapshaper-and-r/↩︎