Import Pakages
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.1
library(psych)
## Warning: package 'psych' was built under R version 4.5.1
library(sf)
## Warning: package 'sf' was built under R version 4.5.1
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.1
##
## 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(tmap)
## Warning: package 'tmap' was built under R version 4.5.1
library(spdep)
## Warning: package 'spdep' was built under R version 4.5.1
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.5.1
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(GWmodel)
## Warning: package 'GWmodel' was built under R version 4.5.1
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 4.5.1
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.5.1
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
Import Data
data <- read_excel("C:/PERKULIAHAN/SEMESTER 7/Kapita Selekta Sains Data/UTS/Data IPP.xlsx")
head(data)
## # A tibble: 6 × 7
## PROVINSI ipp miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aceh 60 14.5 0.296 1178. 67.8 299
## 2 Bali 62.7 4.25 0.362 270. 80.9 21
## 3 Banten 53.3 6.17 0.368 997. 89.1 342
## 4 Bengkulu 57.3 14.0 0.333 703. 77.6 119
## 5 Daerah Isti… 73.3 11.0 0.449 279. 79.0 5
## 6 DKI Jakarta 52 4.44 0.431 109. 87.0 337
Statistika Deskriptif Biasa (Non Spasial)
str(data)
## tibble [34 × 7] (S3: tbl_df/tbl/data.frame)
## $ PROVINSI : chr [1:34] "Aceh" "Bali" "Banten" "Bengkulu" ...
## $ ipp : num [1:34] 60 62.7 53.3 57.3 73.3 ...
## $ miskin : num [1:34] 14.5 4.25 6.17 14.04 11.04 ...
## $ gini_ratio : num [1:34] 0.296 0.362 0.368 0.333 0.449 0.431 0.417 0.343 0.425 0.369 ...
## $ dau_pendidikan : num [1:34] 1178 270 997 703 279 ...
## $ internet : num [1:34] 67.8 80.9 89.1 77.6 79 ...
## $ putus_sekolah_sma: num [1:34] 299 21 342 119 5 337 168 191 688 166 ...
describe(data[ , c("ipp", "miskin", "gini_ratio", "dau_pendidikan", "internet", "putus_sekolah_sma")])
## vars n mean sd median trimmed mad min max
## ipp 1 34 56.09 4.02 55.42 55.60 2.59 51.17 73.33
## miskin 2 34 10.09 5.18 8.43 9.41 4.26 4.25 26.03
## gini_ratio 3 34 0.34 0.05 0.34 0.34 0.05 0.24 0.45
## dau_pendidikan 4 34 1128.71 929.60 844.74 975.05 586.39 109.23 4078.64
## internet 5 34 76.19 5.93 76.63 76.41 4.87 60.78 89.10
## putus_sekolah_sma 6 34 259.21 255.98 179.50 218.07 166.79 1.26 924.00
## range skew kurtosis se
## ipp 22.16 2.30 7.48 0.69
## miskin 21.78 1.13 0.89 0.89
## gini_ratio 0.20 0.25 -0.45 0.01
## dau_pendidikan 3969.41 1.55 1.99 159.43
## internet 28.32 -0.38 0.26 1.02
## putus_sekolah_sma 922.74 1.37 0.94 43.90
Peta Tematik IPP
shapefile_data <- st_read("C:/PERKULIAHAN/SEMESTER 7/Kapita Selekta Sains Data/UTS/Batas Provinsi Desember 2019 Dukcapil/batas_provinsi_desember_2019_dukcapil.shp")
## Reading layer `batas_provinsi_desember_2019_dukcapil' from data source
## `C:\PERKULIAHAN\SEMESTER 7\Kapita Selekta Sains Data\UTS\Batas Provinsi Desember 2019 Dukcapil\batas_provinsi_desember_2019_dukcapil.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## z_range: zmin: 2.65e-05 zmax: 2.65e-05
## Geodetic CRS: WGS 84
print(shapefile_data)
## Simple feature collection with 34 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## z_range: zmin: 2.65e-05 zmax: 2.65e-05
## Geodetic CRS: WGS 84
## First 10 features:
## OBJECTID PROVINSI SHAPE_LENG SHAPE_AREA AKB
## 1 1 ACEH 27.455786 4.62543602 19.41
## 2 2 BALI 6.026646 0.45871702 13.26
## 3 3 BANTEN 9.282228 0.76491115 13.83
## 4 4 BENGKULU 11.706367 1.63012892 0.00
## 5 5 DAERAH ISTIMEWA YOGYAKARTA 3.342892 0.26012975 10.90
## 6 6 DKI JAKARTA 3.123689 0.05342567 10.38
## 7 7 GORONTALO 11.013237 0.97760000 29.47
## 8 8 JAMBI 11.835072 3.97771124 0.00
## 9 9 JAWA BARAT 11.614950 3.03278484 13.56
## 10 10 JAWA TENGAH 15.456349 2.81983779 12.77
## geometry
## 1 MULTIPOLYGON Z (((97.39178 ...
## 2 MULTIPOLYGON Z (((115.1251 ...
## 3 MULTIPOLYGON Z (((105.5498 ...
## 4 MULTIPOLYGON Z (((102.3862 ...
## 5 MULTIPOLYGON Z (((110.8198 ...
## 6 MULTIPOLYGON Z (((106.8768 ...
## 7 MULTIPOLYGON Z (((121.4254 ...
## 8 MULTIPOLYGON Z (((104.4071 ...
## 9 MULTIPOLYGON Z (((108.685 -...
## 10 MULTIPOLYGON Z (((108.8835 ...
unique(shapefile_data$PROVINSI)
## [1] "ACEH" "BALI"
## [3] "BANTEN" "BENGKULU"
## [5] "DAERAH ISTIMEWA YOGYAKARTA" "DKI JAKARTA"
## [7] "GORONTALO" "JAMBI"
## [9] "JAWA BARAT" "JAWA TENGAH"
## [11] "JAWA TIMUR" "KALIMANTAN BARAT"
## [13] "KALIMANTAN SELATAN" "KALIMANTAN TENGAH"
## [15] "KALIMANTAN TIMUR" "KALIMANTAN UTARA"
## [17] "KEPULAUAN BANGKA BELITUNG" "KEPULAUAN RIAU"
## [19] "LAMPUNG" "MALUKU"
## [21] "MALUKU UTARA" "NUSA TENGGARA BARAT"
## [23] "NUSA TENGGARA TIMUR" "PAPUA"
## [25] "PAPUA BARAT" "RIAU"
## [27] "SULAWESI BARAT" "SULAWESI SELATAN"
## [29] "SULAWESI TENGAH" "SULAWESI TENGGARA"
## [31] "SULAWESI UTARA" "SUMATERA BARAT"
## [33] "SUMATERA SELATAN" "SUMATERA UTARA"
unique(data$PROVINSI)
## [1] "Aceh" "Bali"
## [3] "Banten" "Bengkulu"
## [5] "Daerah Istimewa Yogyakarta" "DKI Jakarta"
## [7] "Gorontalo" "Jambi"
## [9] "Jawa Barat" "Jawa Tengah"
## [11] "Jawa Timur" "Kalimantan Barat"
## [13] "Kalimantan Selatan" "Kalimantan Tengah"
## [15] "Kalimantan Timur" "Kalimantan Utara"
## [17] "Kepulauan Bangka Belitung" "Kepulauan Riau"
## [19] "Lampung" "Maluku"
## [21] "Maluku Utara" "Nusa Tenggara Barat"
## [23] "Nusa Tenggara Timur" "Papua"
## [25] "Papua Barat" "Riau"
## [27] "Sulawesi Barat" "Sulawesi Selatan"
## [29] "Sulawesi Tengah" "Sulawesi Tenggara"
## [31] "Sulawesi Utara" "Sumatera Barat"
## [33] "Sumatera Selatan" "Sumatera Utara"
data$PROVINSI <- toupper(data$PROVINSI)
setdiff(shapefile_data$PROVINSI, data$PROVINSI)
## character(0)
setdiff(data$PROVINSI, shapefile_data$PROVINSI)
## character(0)
indo_data <- shapefile_data %>%
left_join(data, by = "PROVINSI")
tm_shape(indo_data) +
tm_polygons("ipp",
palette = "YlGnBu",
title = "Indeks Pembangunan Pemuda (IPP)") +
tm_layout(main.title = "Persebaran IPP di Indonesia",
legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGnBu" is named
## "brewer.yl_gn_bu"
## Multiple palettes called "yl_gn_bu" found: "brewer.yl_gn_bu", "matplotlib.yl_gn_bu". The first one, "brewer.yl_gn_bu", is returned.
Uji Global Moran’s I dan Breusch-Pagan Test
indo_data <- st_make_valid(indo_data)
# Membentuk matriks ketetanggaan (neighbors & weights)
coords <- st_centroid(st_geometry(indo_data)) %>% st_coordinates()
knn_nb <- knn2nb(knearneigh(coords, k = 4)) # k = 4 contohnya
listw_knn <- nb2listw(knn_nb, style = "W", zero.policy = TRUE)
# Uji Global Moran's I
moran_res <- moran.test(indo_data$ipp, listw_knn, zero.policy = TRUE)
print(moran_res)
##
## Moran I test under randomisation
##
## data: indo_data$ipp
## weights: listw_knn
##
## Moran I statistic standard deviate = 0.43468, p-value = 0.3319
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.008266999 -0.030303030 0.007873540
# Breusch-Pagan test
ols_model <- lm(ipp ~ miskin + gini_ratio + dau_pendidikan + internet + putus_sekolah_sma, data = indo_data)
bp_test <- bptest(ols_model)
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: ols_model
## BP = 11.743, df = 5, p-value = 0.03848
Best Bandwidth
Fungsi untuk mencari Best Bandwidth
# Mengonversi sf ke Spatial
indo_sp <- as(indo_data, "Spatial")
# Formula analisis
formula_ipp <- ipp ~ miskin + gini_ratio + dau_pendidikan + internet + putus_sekolah_sma
# Membuat fungsi pencarian bandwidth
custom_bw_search <- function(formula, data, kernel, adaptive = TRUE) {
best_bw <- NULL
tryCatch({
best_bw <- bw.gwr(formula,
data = data,
approach = "AICc",
kernel = kernel,
adaptive = adaptive)
}, error = function(e) {
cat("Error in bw.gwr:", e$message, "\n")
})
return(best_bw)
}
# Kernel yang akan diuji
kernel_types1 <- c("bisquare")
kernel_types2 <- c("gaussian")
best_bandwidths <- list()
best_bandwidths_fixed <- list()
a. Bandwidth Adaptive Bisquare
for (kernel_type in kernel_types1) {
best_bandwidth <- custom_bw_search(formula_ipp, indo_sp, kernel_type, adaptive = TRUE)
best_bandwidths[[kernel_type]] <- best_bandwidth
}
## Adaptive bandwidth (number of nearest neighbours): 28 AICc value: 207.814
## Adaptive bandwidth (number of nearest neighbours): 25 AICc value: 212.7475
## Adaptive bandwidth (number of nearest neighbours): 30 AICc value: 204.6235
## Adaptive bandwidth (number of nearest neighbours): 31 AICc value: 201.6446
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 199.9734
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 199.9734
for (kernel_type in kernel_types1) {
cat("Best adaptive bandwidth for", kernel_type, "kernel:", best_bandwidths[[kernel_type]], "\n")
}
## Best adaptive bandwidth for bisquare kernel: 32
b. Bandwidth Fixed Bisquare
for (kernel_type in kernel_types1) {
best_bandwidth <- custom_bw_search(formula_ipp, indo_sp, kernel_type, adaptive = FALSE)
best_bandwidths_fixed[[kernel_type]] <- best_bandwidth
}
## Fixed bandwidth: 26.42036 AICc value: 201.3115
## Fixed bandwidth: 16.33194 AICc value: 220.2038
## Fixed bandwidth: 32.65534 AICc value: 198.4297
## Fixed bandwidth: 36.50877 AICc value: 197.8225
## Fixed bandwidth: 38.89032 AICc value: 197.5784
## Fixed bandwidth: 40.3622 AICc value: 197.4622
## Fixed bandwidth: 41.27187 AICc value: 197.402
## Fixed bandwidth: 41.83408 AICc value: 197.3688
## Fixed bandwidth: 42.18154 AICc value: 197.3498
## Fixed bandwidth: 42.39629 AICc value: 197.3385
## Fixed bandwidth: 42.52901 AICc value: 197.3318
## Fixed bandwidth: 42.61103 AICc value: 197.3277
## Fixed bandwidth: 42.66172 AICc value: 197.3252
## Fixed bandwidth: 42.69306 AICc value: 197.3237
## Fixed bandwidth: 42.71242 AICc value: 197.3228
## Fixed bandwidth: 42.72439 AICc value: 197.3222
## Fixed bandwidth: 42.73178 AICc value: 197.3218
## Fixed bandwidth: 42.73635 AICc value: 197.3216
## Fixed bandwidth: 42.73918 AICc value: 197.3215
## Fixed bandwidth: 42.74092 AICc value: 197.3214
for (kernel_type in kernel_types1) {
cat("Best fixed bandwidth for", kernel_type, "kernel:", best_bandwidths_fixed[[kernel_type]], "\n")
}
## Best fixed bandwidth for bisquare kernel: 42.74092
c. Bandwidth Fixed Gaussian
for (kernel_type in kernel_types2) {
best_bandwidth <- custom_bw_search(formula_ipp, indo_sp, kernel_type, adaptive = FALSE)
best_bandwidths_fixed[[kernel_type]] <- best_bandwidth
}
## Fixed bandwidth: 26.42036 AICc value: 197.4153
## Fixed bandwidth: 16.33194 AICc value: 197.79
## Fixed bandwidth: 32.65534 AICc value: 197.5051
## Fixed bandwidth: 22.56692 AICc value: 197.3907
## Fixed bandwidth: 20.18537 AICc value: 197.4298
## Fixed bandwidth: 24.0388 AICc value: 197.3926
## Fixed bandwidth: 21.65725 AICc value: 197.3976
## Fixed bandwidth: 23.12913 AICc value: 197.3898
## Fixed bandwidth: 23.4766 AICc value: 197.3904
## Fixed bandwidth: 22.91439 AICc value: 197.3899
for (kernel_type in kernel_types2) {
cat("Best fixed bandwidth for", kernel_type, "kernel:", best_bandwidths_fixed[[kernel_type]], "\n")
}
## Best fixed bandwidth for gaussian kernel: 23.12913
d. Bandwidth Adaptive Gaussian
for (kernel_type in kernel_types2) {
best_bandwidth <- custom_bw_search(formula_ipp, indo_sp, kernel_type, adaptive = TRUE)
best_bandwidths[[kernel_type]] <- best_bandwidth
}
## Adaptive bandwidth (number of nearest neighbours): 28 AICc value: 196.8704
## Adaptive bandwidth (number of nearest neighbours): 25 AICc value: 196.0434
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 197.1092
## Adaptive bandwidth (number of nearest neighbours): 25 AICc value: 196.0434
for (kernel_type in kernel_types2) {
cat("Best adaptive bandwidth for", kernel_type, "kernel:", best_bandwidths[[kernel_type]], "\n")
}
## Best adaptive bandwidth for gaussian kernel: 25
Pembentukan Model
a. Model Fixed Gaussian
kernel_type1 <- "gaussian"
gwr_model_gaussfix <- gwr.basic(formula_ipp, data = indo_sp,
bw = 23.12913,
kernel = kernel_type1,
adaptive = FALSE)
print(gwr_model_gaussfix)
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-10-16 21:29:45.231591
## Call:
## gwr.basic(formula = formula_ipp, data = indo_sp, bw = 23.12913,
## kernel = kernel_type1, adaptive = FALSE)
##
## Dependent (y) variable: ipp
## Independent variables: miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## Number of data points: 34
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0825 -1.8675 -0.9381 2.1737 11.7657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.0949007 9.8170418 6.020 1.74e-06 ***
## miskin -0.0560155 0.1494157 -0.375 0.7106
## gini_ratio 40.7790915 16.2343653 2.512 0.0181 *
## dau_pendidikan -0.0005748 0.0007956 -0.723 0.4760
## internet -0.1902551 0.1311815 -1.450 0.1581
## putus_sekolah_sma -0.0051399 0.0031168 -1.649 0.1103
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 3.739 on 28 degrees of freedom
## Multiple R-squared: 0.265
## Adjusted R-squared: 0.1337
## F-statistic: 2.019 on 5 and 28 DF, p-value: 0.1066
## ***Extra Diagnostic information
## Residual sum of squares: 391.4779
## Sigma(hat): 3.497668
## AIC: 193.5691
## AICc: 197.8768
## BIC: 194.9382
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: gaussian
## Fixed bandwidth: 23.12913
## 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 54.82734420 55.35472802 56.30880076 57.19059883 59.5421
## miskin -0.14002192 -0.05171164 -0.00178854 0.04966324 0.1248
## gini_ratio 39.57440756 41.47671544 42.16094660 43.03458599 43.8183
## dau_pendidikan -0.00066333 -0.00065115 -0.00061847 -0.00058783 -0.0005
## internet -0.18817485 -0.17587570 -0.16799879 -0.15922663 -0.1456
## putus_sekolah_sma -0.00603783 -0.00591191 -0.00543154 -0.00505091 -0.0042
## ************************Diagnostic information*************************
## Number of data points: 34
## Effective number of parameters (2trace(S) - trace(S'S)): 7.54586
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 26.45414
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 197.3898
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 183.021
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 166.3208
## Residual sum of squares: 354.2663
## R-square value: 0.3348546
## Adjusted R-square value: 0.1376727
##
## ***********************************************************************
## Program stops at: 2025-10-16 21:29:45.315535
b. Model Adaptive Gaussian
kernel_type2 <- "gaussian"
gwr_model_gaussadap <- gwr.basic(formula_ipp, data = indo_sp,
bw = 25,
kernel = kernel_type2,
adaptive = TRUE)
print(gwr_model_gaussadap)
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-10-16 21:29:45.359055
## Call:
## gwr.basic(formula = formula_ipp, data = indo_sp, bw = 25, kernel = kernel_type2,
## adaptive = TRUE)
##
## Dependent (y) variable: ipp
## Independent variables: miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## Number of data points: 34
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0825 -1.8675 -0.9381 2.1737 11.7657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.0949007 9.8170418 6.020 1.74e-06 ***
## miskin -0.0560155 0.1494157 -0.375 0.7106
## gini_ratio 40.7790915 16.2343653 2.512 0.0181 *
## dau_pendidikan -0.0005748 0.0007956 -0.723 0.4760
## internet -0.1902551 0.1311815 -1.450 0.1581
## putus_sekolah_sma -0.0051399 0.0031168 -1.649 0.1103
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 3.739 on 28 degrees of freedom
## Multiple R-squared: 0.265
## Adjusted R-squared: 0.1337
## F-statistic: 2.019 on 5 and 28 DF, p-value: 0.1066
## ***Extra Diagnostic information
## Residual sum of squares: 391.4779
## Sigma(hat): 3.497668
## AIC: 193.5691
## AICc: 197.8768
## BIC: 194.9382
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: gaussian
## Adaptive bandwidth: 25 (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 48.50653224 50.68812548 53.21848914 55.41083127 59.1951
## miskin -0.10522091 -0.04512274 0.11138802 0.18556609 0.2111
## gini_ratio 39.82510733 41.59449455 43.51020958 45.97283028 49.9863
## dau_pendidikan -0.00079224 -0.00072809 -0.00067891 -0.00061635 -0.0005
## internet -0.19052251 -0.15945671 -0.13635727 -0.11933483 -0.1010
## putus_sekolah_sma -0.00724741 -0.00665606 -0.00593829 -0.00478868 -0.0045
## ************************Diagnostic information*************************
## Number of data points: 34
## Effective number of parameters (2trace(S) - trace(S'S)): 8.728306
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 25.27169
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 196.0434
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 179.7769
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 164.8852
## Residual sum of squares: 315.3179
## R-square value: 0.4079813
## Adjusted R-square value: 0.1950863
##
## ***********************************************************************
## Program stops at: 2025-10-16 21:29:45.511036
c. Model Fixed Bisquare
kernel_type3 <- "bisquare"
gwr_model_bisfix <- gwr.basic(formula_ipp, data = indo_sp,
bw = 42.74092,
kernel = kernel_type3,
adaptive = FALSE)
print(gwr_model_bisfix)
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-10-16 21:29:45.555346
## Call:
## gwr.basic(formula = formula_ipp, data = indo_sp, bw = 42.74092,
## kernel = kernel_type3, adaptive = FALSE)
##
## Dependent (y) variable: ipp
## Independent variables: miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## Number of data points: 34
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0825 -1.8675 -0.9381 2.1737 11.7657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.0949007 9.8170418 6.020 1.74e-06 ***
## miskin -0.0560155 0.1494157 -0.375 0.7106
## gini_ratio 40.7790915 16.2343653 2.512 0.0181 *
## dau_pendidikan -0.0005748 0.0007956 -0.723 0.4760
## internet -0.1902551 0.1311815 -1.450 0.1581
## putus_sekolah_sma -0.0051399 0.0031168 -1.649 0.1103
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 3.739 on 28 degrees of freedom
## Multiple R-squared: 0.265
## Adjusted R-squared: 0.1337
## F-statistic: 2.019 on 5 and 28 DF, p-value: 0.1066
## ***Extra Diagnostic information
## Residual sum of squares: 391.4779
## Sigma(hat): 3.497668
## AIC: 193.5691
## AICc: 197.8768
## BIC: 194.9382
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: bisquare
## Fixed bandwidth: 42.74092
## 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 49.92489823 53.14777306 55.23050258 56.58730157 59.0014
## miskin -0.17276597 -0.05078207 0.02325179 0.11895875 0.3103
## gini_ratio 38.73044777 41.84080638 42.58400458 43.55488575 44.7163
## dau_pendidikan -0.00074267 -0.00069316 -0.00063824 -0.00059484 -0.0004
## internet -0.17787502 -0.17157659 -0.15859699 -0.13803012 -0.0942
## putus_sekolah_sma -0.00689214 -0.00625379 -0.00555851 -0.00503985 -0.0037
## ************************Diagnostic information*************************
## Number of data points: 34
## Effective number of parameters (2trace(S) - trace(S'S)): 8.028703
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 25.9713
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 197.3214
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 182.1044
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 166.2303
## Residual sum of squares: 341.5436
## R-square value: 0.3587418
## Adjusted R-square value: 0.1525662
##
## ***********************************************************************
## Program stops at: 2025-10-16 21:29:45.6511
d. Model Adaptive Bisquare
kernel_type4 <- "bisquare"
gwr_model_bisadapt <- gwr.basic(formula_ipp, data = indo_sp,
bw = 32,
kernel = kernel_type4,
adaptive = TRUE)
print(gwr_model_bisadapt)
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-10-16 21:29:45.686956
## Call:
## gwr.basic(formula = formula_ipp, data = indo_sp, bw = 32, kernel = kernel_type4,
## adaptive = TRUE)
##
## Dependent (y) variable: ipp
## Independent variables: miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## Number of data points: 34
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0825 -1.8675 -0.9381 2.1737 11.7657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.0949007 9.8170418 6.020 1.74e-06 ***
## miskin -0.0560155 0.1494157 -0.375 0.7106
## gini_ratio 40.7790915 16.2343653 2.512 0.0181 *
## dau_pendidikan -0.0005748 0.0007956 -0.723 0.4760
## internet -0.1902551 0.1311815 -1.450 0.1581
## putus_sekolah_sma -0.0051399 0.0031168 -1.649 0.1103
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 3.739 on 28 degrees of freedom
## Multiple R-squared: 0.265
## Adjusted R-squared: 0.1337
## F-statistic: 2.019 on 5 and 28 DF, p-value: 0.1066
## ***Extra Diagnostic information
## Residual sum of squares: 391.4779
## Sigma(hat): 3.497668
## AIC: 193.5691
## AICc: 197.8768
## BIC: 194.9382
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: bisquare
## Adaptive bandwidth: 32 (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 40.62984016 43.73764260 46.22734302 48.98907960 59.3907
## miskin -0.18898724 0.00108198 0.27741644 0.40702707 0.4563
## gini_ratio 34.52113779 39.05060201 44.30051841 49.62716034 58.0914
## dau_pendidikan -0.00097247 -0.00085993 -0.00078190 -0.00071585 -0.0003
## internet -0.15663146 -0.09765541 -0.06230331 -0.05755196 -0.0459
## putus_sekolah_sma -0.00898504 -0.00837856 -0.00686074 -0.00383266 -0.0027
## ************************Diagnostic information*************************
## Number of data points: 34
## Effective number of parameters (2trace(S) - trace(S'S)): 11.97666
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 22.02334
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 199.9734
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 176.0864
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 167.2714
## Residual sum of squares: 263.5633
## R-square value: 0.5051521
## Adjusted R-square value: 0.2232453
##
## ***********************************************************************
## Program stops at: 2025-10-16 21:29:45.846245
Evaluasi Model
# Membuat data frame kosong untuk menyimpan hasil evaluasi
results_df <- data.frame(
model = character(),
kernel = character(),
adaptive = logical(),
aic = numeric(),
aicc = numeric(),
bic = numeric(),
rsquared = numeric(),
adj_rsquared = numeric(),
stringsAsFactors = FALSE
)
# Membuat fungsi untuk ekstrak hasil evaluasi dari model GWR
extract_results <- function(model_obj, model_name, kernel_type, adaptive_flag) {
return(data.frame(
model = model_name,
kernel = kernel_type,
adaptive = adaptive_flag,
aic = model_obj$GW.diagnostic$AIC,
aicc = model_obj$GW.diagnostic$AICc,
bic = model_obj$GW.diagnostic$BIC,
rsquared = model_obj$GW.diagnostic$gw.R2,
adj_rsquared = model_obj$GW.diagnostic$gwR2.adj,
stringsAsFactors = FALSE
))
}
# Menambahkan semua model ke results_df
results_df <- rbind(results_df,
extract_results(gwr_model_bisadapt, "GWR Adaptive Bisquare", "bisquare", TRUE),
extract_results(gwr_model_bisfix, "GWR Fixed Bisquare", "bisquare", FALSE),
extract_results(gwr_model_gaussadap, "GWR Adaptive Gaussian", "gaussian", TRUE),
extract_results(gwr_model_gaussfix, "GWR Fixed Gaussian", "gaussian", FALSE)
)
# Menampilkan hasil final
print(results_df)
## model kernel adaptive aic aicc bic rsquared
## 1 GWR Adaptive Bisquare bisquare TRUE 176.0864 199.9734 167.2714 0.5051521
## 2 GWR Fixed Bisquare bisquare FALSE 182.1044 197.3214 166.2303 0.3587418
## 3 GWR Adaptive Gaussian gaussian TRUE 179.7769 196.0434 164.8852 0.4079813
## 4 GWR Fixed Gaussian gaussian FALSE 183.0210 197.3898 166.3208 0.3348546
## adj_rsquared
## 1 0.2232453
## 2 0.1525662
## 3 0.1950863
## 4 0.1376727
# Copy dari results_df
df <- results_df
# Ranking AIC/AICc/BIC (lebih kecil lebih baik)
df$rank_AIC <- rank(df$aic, ties.method = "min")
df$rank_AICc <- rank(df$aicc, ties.method = "min")
df$rank_BIC <- rank(df$bic, ties.method = "min")
# Ranking R2 dan Adjusted R2 (lebih besar lebih baik)
df$rank_R2 <- rank(-df$rsquared, ties.method = "min")
df$rank_AdjR2 <- rank(-df$adj_rsquared, ties.method = "min")
# Menghitung total ranking (semakin kecil = semakin baik)
df$total_score <- df$rank_AIC + df$rank_AICc + df$rank_BIC + df$rank_R2 + df$rank_AdjR2
# Mengurutkan berdasarkan total_score
df_ranked <- df[order(df$total_score), ]
print(df_ranked)
## model kernel adaptive aic aicc bic rsquared
## 3 GWR Adaptive Gaussian gaussian TRUE 179.7769 196.0434 164.8852 0.4079813
## 1 GWR Adaptive Bisquare bisquare TRUE 176.0864 199.9734 167.2714 0.5051521
## 2 GWR Fixed Bisquare bisquare FALSE 182.1044 197.3214 166.2303 0.3587418
## 4 GWR Fixed Gaussian gaussian FALSE 183.0210 197.3898 166.3208 0.3348546
## adj_rsquared rank_AIC rank_AICc rank_BIC rank_R2 rank_AdjR2 total_score
## 3 0.1950863 2 1 1 2 2 8
## 1 0.2232453 1 4 4 1 1 11
## 2 0.1525662 3 2 2 3 3 13
## 4 0.1376727 4 3 3 4 4 18
Model GWR Adaptive Gaussian
# Model terbaik
gwr_model_best <- gwr_model_gaussadap
# Ekstrak Koefisien & T-Value
head(gwr_model_best$SDF@data[, c("miskin", "gini_ratio", "dau_pendidikan", "internet", "putus_sekolah_sma")])
## miskin gini_ratio dau_pendidikan internet putus_sekolah_sma
## 1 0.10354079 39.82511 -0.0006547596 -0.1508077 -0.005925672
## 2 0.04979091 48.29174 -0.0007060717 -0.1397292 -0.005992706
## 3 0.20669622 45.97156 -0.0007538010 -0.1165061 -0.007239812
## 4 0.16925371 42.55838 -0.0007073129 -0.1306901 -0.006680413
## 5 0.18926652 49.98629 -0.0007922402 -0.1108380 -0.007247414
## 6 0.19035671 45.97325 -0.0007464462 -0.1192071 -0.007081634
head(gwr_model_best$SDF@data[, c("miskin_TV", "gini_ratio_TV", "dau_pendidikan_TV", "internet_TV", "putus_sekolah_sma_TV")])
## miskin_TV gini_ratio_TV dau_pendidikan_TV internet_TV putus_sekolah_sma_TV
## 1 0.6589513 2.536258 -0.8629603 -1.1804298 -1.976267
## 2 0.3202900 3.087601 -0.9185193 -1.0941556 -1.993704
## 3 1.1383707 2.864438 -0.9887145 -0.8855715 -2.371208
## 4 0.9870493 2.675870 -0.9316402 -1.0037668 -2.206603
## 5 1.0472810 3.115630 -1.0277012 -0.8459355 -2.367054
## 6 1.0719730 2.879720 -0.9802955 -0.9131993 -2.330006
Uji Parsial
# Signifikansi Parsial
gwr_model_best$SDF@data$signif_miskin <- ifelse(abs(gwr_model_best$SDF@data$miskin_TV) > 1.96, "Signifikan", "Tidak")
gwr_model_best$SDF@data$signif_gini <- ifelse(abs(gwr_model_best$SDF@data$gini_ratio_TV) > 1.96, "Signifikan", "Tidak")
gwr_model_best$SDF@data$signif_dau_pendidikan <- ifelse(abs(gwr_model_best$SDF@data$dau_pendidikan_TV) > 1.96, "Signifikan", "Tidak")
gwr_model_best$SDF@data$signif_internet <- ifelse(abs(gwr_model_best$SDF@data$internet_TV) > 1.96, "Signifikan", "Tidak")
gwr_model_best$SDF@data$signif_putus_sma <- ifelse(abs(gwr_model_best$SDF@data$putus_sekolah_sma_TV) > 1.96, "Signifikan", "Tidak")
uji_parsial_df <- gwr_model_best$SDF@data[, c("miskin", "miskin_TV", "signif_miskin",
"gini_ratio", "gini_ratio_TV", "signif_gini",
"dau_pendidikan", "dau_pendidikan_TV", "signif_dau_pendidikan",
"internet", "internet_TV", "signif_internet",
"putus_sekolah_sma", "putus_sekolah_sma_TV", "signif_putus_sma")]
uji_parsial_df$PROVINSI <- indo_data$PROVINSI
uji_parsial_df <- uji_parsial_df[, c("PROVINSI",
"miskin", "miskin_TV", "signif_miskin",
"gini_ratio", "gini_ratio_TV", "signif_gini",
"dau_pendidikan", "dau_pendidikan_TV", "signif_dau_pendidikan",
"internet", "internet_TV", "signif_internet",
"putus_sekolah_sma", "putus_sekolah_sma_TV", "signif_putus_sma")]
head(uji_parsial_df)
## PROVINSI miskin miskin_TV signif_miskin gini_ratio
## 1 ACEH 0.10354079 0.6589513 Tidak 39.82511
## 2 BALI 0.04979091 0.3202900 Tidak 48.29174
## 3 BANTEN 0.20669622 1.1383707 Tidak 45.97156
## 4 BENGKULU 0.16925371 0.9870493 Tidak 42.55838
## 5 DAERAH ISTIMEWA YOGYAKARTA 0.18926652 1.0472810 Tidak 49.98629
## 6 DKI JAKARTA 0.19035671 1.0719730 Tidak 45.97325
## gini_ratio_TV signif_gini dau_pendidikan dau_pendidikan_TV
## 1 2.536258 Signifikan -0.0006547596 -0.8629603
## 2 3.087601 Signifikan -0.0007060717 -0.9185193
## 3 2.864438 Signifikan -0.0007538010 -0.9887145
## 4 2.675870 Signifikan -0.0007073129 -0.9316402
## 5 3.115630 Signifikan -0.0007922402 -1.0277012
## 6 2.879720 Signifikan -0.0007464462 -0.9802955
## signif_dau_pendidikan internet internet_TV signif_internet
## 1 Tidak -0.1508077 -1.1804298 Tidak
## 2 Tidak -0.1397292 -1.0941556 Tidak
## 3 Tidak -0.1165061 -0.8855715 Tidak
## 4 Tidak -0.1306901 -1.0037668 Tidak
## 5 Tidak -0.1108380 -0.8459355 Tidak
## 6 Tidak -0.1192071 -0.9131993 Tidak
## putus_sekolah_sma putus_sekolah_sma_TV signif_putus_sma
## 1 -0.005925672 -1.976267 Signifikan
## 2 -0.005992706 -1.993704 Signifikan
## 3 -0.007239812 -2.371208 Signifikan
## 4 -0.006680413 -2.206603 Signifikan
## 5 -0.007247414 -2.367054 Signifikan
## 6 -0.007081634 -2.330006 Signifikan
Peta Tematik Variabel Signifikan
# Menggabungkan hasil uji parsial dengan shapefile
indo_gwr <- indo_data %>%
left_join(uji_parsial_df, by = "PROVINSI")
# Kolom kombinasi signifikan
indo_gwr <- indo_gwr %>%
mutate(
Kombinasi_Signifikan = case_when(
signif_gini == "Signifikan" & signif_putus_sma == "Signifikan" ~ "1. Gini Ratio & SMA Dropout Rate",
signif_gini == "Signifikan" & signif_putus_sma == "Tidak" ~ "2. Gini Ratio Only",
TRUE ~ NA_character_
)
)
# Mendefinisikan level dan palet warna
custom_levels <- c(
"1. Gini Ratio & SMA Dropout Rate",
"2. Gini Ratio Only"
)
custom_palette_final <- c(
"1. Gini Ratio & SMA Dropout Rate" = "#F4A6A6",
"2. Gini Ratio Only" = "#8DB5E2"
)
indo_gwr$Kombinasi_Signifikan <- factor(
indo_gwr$Kombinasi_Signifikan,
levels = custom_levels
)
# Peta tematik
map_kluster_final <- tm_shape(indo_gwr) +
tm_polygons(
fill = "Kombinasi_Signifikan",
palette = custom_palette_final,
title = "Kombinasi_Signifikan",
colorNA = "#FFFFFF",
textNA = "Missing"
) +
tm_borders(col = "grey40", lwd = 0.5) +
tm_layout(
main.title = "Significant Variable Clusters Against IPP (GWR)",
main.title.size = 1.2,
legend.outside = TRUE,
legend.title.size = 1,
frame = FALSE,
bg.color = "white"
)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values'), 'colorNA'
## (rename to 'value.na'), 'textNA' (rename to 'label.na') to fill.scale =
## tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
## a list: 'fill.scale = list(<scale1>, <scale2>, ...)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
map_kluster_final
