PEMODELAN LOCALLY COMPENSATED RIDGE IMPROVED GEOGRAPHICALLY TEMPORALLY WEIGHTED REGRESSION PADA DATA TINGKAT FERTILITAS (TFR) SULAWESI DAFTAR ISI 0. Pengaturan dan paket 1. Data 2. OLS 3. Bandwidth GWR (pemilihan kernel) 4. Pengujian asumsi (normalitas, heterogenitas spasial, multikolinearitas) 5. GWR 6. Matriks pembobot spasial W, Moran’s I, dan uji LM 7. GWR-SAR: bandwidth, estimator, dan estimasi 8. Perbandingan model
rm(list = ls())
# ==============================================================================
# 0. PENGATURAN DAN PAKET
# ==============================================================================
# ---- 0.1 Lokasi file (sesuaikan) ----
dir_kerja <- "C:/Users/ASUS/Downloads"
file_func <- file.path(dir_kerja, "FUNCTION_LCR-IGTWRne.R")
file_data <- file.path(dir_kerja, "Data_TFR_Sulawesi_2025.xlsx")
sheet_data <- "Data Gabungan"
# ---- 0.2 Pengaturan model ----
alpha <- 0.05 # taraf signifikansi
# ---- 0.3 Paket ----
# install.packages(c("readxl", "openxlsx", "lmtest", "car", "corrplot",
# "spdep", "MASS", "sp", "spgwr", "GWmodel")) # cukup sekali
paket <- c("readxl", # membaca data Excel
"openxlsx", # menyimpan hasil ke Excel
"lmtest", # bptest()
"car", # vif()
"corrplot", # plot korelasi
"spdep", # mat2listw(), moran.test()
"MASS", # ginv()
"sp", "spgwr", "GWmodel") # gold() untuk bandwidth CV
invisible(lapply(paket, library, character.only = TRUE))
## Warning: package 'readxl' was built under R version 4.5.3
## Warning: package 'openxlsx' was built under R version 4.5.3
## Warning: package 'lmtest' was built under R version 4.5.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Warning: package 'car' was built under R version 4.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.3
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
## Warning: package 'spdep' was built under R version 4.5.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.5.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
## Warning: package 'MASS' was built under R version 4.5.3
## Warning: package 'sp' was built under R version 4.5.3
## Warning: package 'spgwr' was built under R version 4.5.3
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
## Warning: package 'GWmodel' was built under R version 4.5.3
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 4.5.3
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.5.3
## Welcome to GWmodel version 2.4-2.
# ---- 0.4 Fungsi buatan sendiri ----
source(file_func)
data <- data.frame(read_excel(file_data, sheet = sheet_data))
y <- data[, 5] # respon: TFR
x <- data[, c(6:10)] # prediktor: RLS, TPAKP, TPT, MISKIN
coord <- data[, 3:4] # koordinat: Longitude, Latitude
n <- nrow(data)
D <- jarak.euclid(coord) # matriks jarak (dipakai di bagian 6 dan 7)
str(data)
## 'data.frame': 81 obs. of 10 variables:
## $ prov : chr "Sulawesi Utara" "Sulawesi Utara" "Sulawesi Utara" "Sulawesi Utara" ...
## $ kabkot : chr "Bolaang Mongondow" "Minahasa" "Kepulauan Sangihe" "Kepulauan Talaud" ...
## $ Longitude : num 282361 373638 450241 588654 335866 ...
## $ Latitude : num 299607 359465 613607 687412 341272 ...
## $ TFR : num 2.16 1.99 2.07 2.16 2.06 2.02 2.19 2.06 2.24 2.27 ...
## $ RLS : num 8.49 10.25 8.7 9.97 9.53 ...
## $ TPAKP : num 54.2 48.9 54.5 50.8 50 ...
## $ TPAKL : num 87.6 80.1 81.2 81 84.4 ...
## $ TPT : num 4.97 7.74 2.64 3.4 4.77 6.8 4.07 2.41 3.38 2.91 ...
## $ MISKIN....: num 6.98 5.88 10.91 8.34 8.26 ...
OLS <- lm(y ~ ., data = data.frame(y = y, x))
summary(OLS)
##
## Call:
## lm(formula = y ~ ., data = data.frame(y = y, x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34952 -0.13151 -0.00616 0.07845 0.50553
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.784685 0.601391 4.630 1.5e-05 ***
## RLS 0.003450 0.024733 0.140 0.88943
## TPAKP 0.008101 0.002895 2.798 0.00653 **
## TPAKL -0.012853 0.006009 -2.139 0.03570 *
## TPT -0.034674 0.018505 -1.874 0.06486 .
## MISKIN.... 0.025759 0.007585 3.396 0.00110 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1793 on 75 degrees of freedom
## Multiple R-squared: 0.4327, Adjusted R-squared: 0.3949
## F-statistic: 11.44 on 5 and 75 DF, p-value: 3.182e-08
kernel_uji <- c("exponential", "gaussian", "bisquare", "tricube")
bw_kernel <- sapply(kernel_uji, function(k)
bwd.gwr.twr(y = y, x = x, long.lat = coord, spasial = TRUE, kernel = k))
## Fixed bandwidth: 913788.5 CV score: 2.5205
## Fixed bandwidth: 564865.3 CV score: 2.340359
## Fixed bandwidth: 349218.9 CV score: 2.103019
## Fixed bandwidth: 215942.1 CV score: 1.845772
## Fixed bandwidth: 133572.5 CV score: 1.64356
## Fixed bandwidth: 82665.27 CV score: 1.518837
## Fixed bandwidth: 51202.88 CV score: 1.463511
## Fixed bandwidth: 31758.06 CV score: 1.688484
## Fixed bandwidth: 63220.45 CV score: 1.468245
## Fixed bandwidth: 43775.62 CV score: 1.509792
## Fixed bandwidth: 55793.18 CV score: 1.460464
## Fixed bandwidth: 58630.15 CV score: 1.461897
## Fixed bandwidth: 54039.85 CV score: 1.46072
## Fixed bandwidth: 56876.81 CV score: 1.460768
## Fixed bandwidth: 55123.47 CV score: 1.460446
## Fixed bandwidth: 913788.5 CV score: 2.725857
## Fixed bandwidth: 564865.3 CV score: 2.526895
## Fixed bandwidth: 349218.9 CV score: 2.166979
## Fixed bandwidth: 215942.1 CV score: 1.850646
## Fixed bandwidth: 133572.5 CV score: 1.637793
## Fixed bandwidth: 82665.27 CV score: 1.560676
## Fixed bandwidth: 51202.88 CV score: 2.571048
## Fixed bandwidth: 102110.1 CV score: 1.552873
## Fixed bandwidth: 114127.7 CV score: 1.56719
## Fixed bandwidth: 94682.83 CV score: 1.549737
## Fixed bandwidth: 90092.53 CV score: 1.55978
## Fixed bandwidth: 97519.8 CV score: 1.548459
## Fixed bandwidth: 99273.13 CV score: 1.549094
## Fixed bandwidth: 96436.17 CV score: 1.548582
## Fixed bandwidth: 98189.51 CV score: 1.548586
## Fixed bandwidth: 97105.89 CV score: 1.548456
## Fixed bandwidth: 913788.5 CV score: 2.229115
## Fixed bandwidth: 564865.3 CV score: 1.914661
## Fixed bandwidth: 349218.9 CV score: 1.702733
## Fixed bandwidth: 215942.1 CV score: 2.246674
## Fixed bandwidth: 431588.5 CV score: 1.850383
## Fixed bandwidth: 298311.7 CV score: 1.690326
## Fixed bandwidth: 266849.3 CV score: 2.568912
## Fixed bandwidth: 317756.5 CV score: 1.682905
## Fixed bandwidth: 329774.1 CV score: 1.676433
## Fixed bandwidth: 337201.3 CV score: 1.686054
## Fixed bandwidth: 325183.8 CV score: 1.692446
## Fixed bandwidth: 332611 CV score: 1.680005
## Fixed bandwidth: 328020.7 CV score: 1.67429
## Fixed bandwidth: 326937.1 CV score: 1.672986
## Fixed bandwidth: 326267.4 CV score: 1.694112
## Fixed bandwidth: 327351 CV score: 1.673481
## Fixed bandwidth: 326681.3 CV score: 1.672685
## Fixed bandwidth: 326523.2 CV score: 1.694514
## Fixed bandwidth: 326779 CV score: 1.6728
## Fixed bandwidth: 326620.9 CV score: 1.672614
## Fixed bandwidth: 913788.5 CV score: 2.294018
## Fixed bandwidth: 564865.3 CV score: 1.961853
## Fixed bandwidth: 349218.9 CV score: 1.732014
## Fixed bandwidth: 215942.1 CV score: 2.569846
## Fixed bandwidth: 431588.5 CV score: 1.938748
## Fixed bandwidth: 298311.7 CV score: 1.719732
## Fixed bandwidth: 266849.3 CV score: 3.160788
## Fixed bandwidth: 317756.5 CV score: 1.718661
## Fixed bandwidth: 329774.1 CV score: 1.730365
## Fixed bandwidth: 310329.2 CV score: 1.720281
## Fixed bandwidth: 322346.8 CV score: 1.721571
## Fixed bandwidth: 314919.5 CV score: 1.718216
## Fixed bandwidth: 313166.2 CV score: 1.718566
## Fixed bandwidth: 316003.2 CV score: 1.718246
bw_kernel
## exponential gaussian bisquare tricube
## 55123.47 97105.89 326620.91 314919.54
bandwith.GWR <- unname(bw_kernel[1]) # bandwidth kernel terpilih
h_gwr <- bandwith.GWR
# ---- 4.1 Normalitas sisaan ----
sisaan_ols <- residuals(OLS)
print(ks.test(sisaan_ols, "pnorm", mean = mean(sisaan_ols), sd = sd(sisaan_ols)))
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan_ols
## D = 0.10947, p-value = 0.2663
## alternative hypothesis: two-sided
print(shapiro.test(sisaan_ols)) # lebih tepat untuk n kecil dan parameter diduga
##
## Shapiro-Wilk normality test
##
## data: sisaan_ols
## W = 0.96985, p-value = 0.05296
# ---- 4.2 Heterogenitas spasial (Breusch-Pagan) ----
BPTest_Spasial <- bptest(OLS, studentize = FALSE)
BPTest_Spasial
##
## Breusch-Pagan test
##
## data: OLS
## BP = 16.376, df = 5, p-value = 0.005848
# ---- 4.3 Multikolinearitas (VIF) ----
Multi <- vif(OLS)
Multi
## RLS TPAKP TPAKL TPT MISKIN....
## 1.879287 1.376575 1.395072 2.403514 1.610396
# ---- 4.4 Korelasi antar variabel ----
Korelasi <- corrplot(cor(data.frame(TFR = y, x)), method = "number")
gwr <- prog.gwr(y = y, x = x, long.lat = coord, bwd = bandwith.GWR, kernel = "exponential")
# ---- 5.1 Estimasi parameter ----
Estimasi_GWR <- data.frame(gwr$beta)
summary(Estimasi_GWR)
## X1 X2 X3 X4
## Min. :-1.360 Min. :-0.15172 Min. :-0.009103 Min. :-0.039108
## 1st Qu.: 1.218 1st Qu.:-0.02596 1st Qu.:-0.001397 1st Qu.:-0.014609
## Median : 1.780 Median : 0.01010 Median : 0.004026 Median :-0.004396
## Mean : 2.074 Mean : 0.01907 Mean : 0.003078 Mean :-0.002695
## 3rd Qu.: 3.197 3rd Qu.: 0.07061 3rd Qu.: 0.006717 3rd Qu.: 0.010542
## Max. : 5.955 Max. : 0.15846 Max. : 0.014776 Max. : 0.031342
## X5 X6
## Min. :-0.12970 Min. :-0.03059
## 1st Qu.:-0.05374 1st Qu.: 0.01044
## Median :-0.03359 Median : 0.01877
## Mean :-0.03424 Mean : 0.02288
## 3rd Qu.:-0.01170 3rd Qu.: 0.03670
## Max. : 0.10758 Max. : 0.07435
# ---- 5.2 Matriks jarak dan pembobot ----
Matriks_Jarak <- data.frame(gwr$jarak)
Matriks_Pembobot_W <- data.frame(gwr$bobot)
# write.csv(Matriks_Pembobot_W, file = "Matriks_Pembobot_W.csv", row.names = FALSE)
# ---- Uji-t lokal GWR biasa ----
inf_gwr <- inferensi.gwr(gwr, y, x)
sig_gwr <- uji.signifikansi(inf_gwr, id = data[, c("prov", "kabkot")], alpha = alpha)
print(sig_gwr$ringkasan, digits = 4)
## Parameter Min Median Mean Max Mean_t N_Signif Prop_Signif
## 1 Intercept -1.359595 1.779968 2.074339 5.95455 1.9803 36 0.4444
## 2 RLS -0.151719 0.010100 0.019073 0.15846 0.2934 9 0.1111
## 3 TPAKP -0.009103 0.004026 0.003078 0.01478 0.8446 13 0.1605
## 4 TPAKL -0.039108 -0.004396 -0.002695 0.03134 -0.2131 17 0.2099
## 5 TPT -0.129701 -0.033587 -0.034245 0.10758 -0.9911 15 0.1852
## 6 MISKIN.... -0.030593 0.018774 0.022879 0.07435 1.4218 33 0.4074
## N_Positif N_Negatif Status
## 1 36 0 Signifikan di sebagian lokasi
## 2 5 4 Signifikan di sebagian lokasi
## 3 13 0 Signifikan di sebagian lokasi
## 4 7 10 Signifikan di sebagian lokasi
## 5 0 15 Signifikan di sebagian lokasi
## 6 30 3 Signifikan di sebagian lokasi
# Buat workbook baru
wb <- createWorkbook()
# Tambahkan sheet untuk Ringkasan
addWorksheet(wb, sheetName = "Ringkasan")
writeData(wb, sheet = "Ringkasan", x = sig_gwr$ringkasan)
# Tambahkan sheet untuk Detail Lokal
addWorksheet(wb, sheetName = "Lokal")
writeData(wb, sheet = "Lokal", x = sig_gwr$lokal)
# Simpan file Excel
saveWorkbook(wb, file = "Uji_Sig_Parameter0.xlsx", overwrite = TRUE)
cat("File 'Uji_Sig_Parameter0.xlsx' berhasil disimpan.\n")
## File 'Uji_Sig_Parameter0.xlsx' berhasil disimpan.
kernel<- "exponential"
# ---- 6.1 W lag: kernel, diagonal 0, row-normalized ----
W_lag <- buat.W.lag(D, h_gwr, kernel)
# ---- 6.2 Moran's I pada sisaan GWR ----
sisaan_gwr <- drop((diag(n) - gwr$hat_matrix) %*% as.numeric(y))
lw <- mat2listw(W_lag, style = "W")
moran_gwr <- moran.test(sisaan_gwr, lw, randomisation = TRUE, alternative = "two.sided")
print(moran_gwr)
##
## Moran I test under randomisation
##
## data: sisaan_gwr
## weights: lw
##
## Moran I statistic standard deviate = -1.9899, p-value = 0.0466
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## -0.105070194 -0.012500000 0.002164029
# ---- 6.3 Uji LM spasial pada sisaan GWR ----
hasil_lm <- lm.spatial.tests.gwr(gwr, y, W_lag)
keputusan <- pilih.model.lm(hasil_lm, alpha = alpha)
print(hasil_lm, digits = 5)
## Uji Statistik df p.value
## 1 LMerr 4.16291 1 4.1318e-02
## 2 LMlag 43.65755 1 3.9117e-11
## 3 RLMerr 0.28697 1 5.9217e-01
## 4 RLMlag 39.78161 1 2.8401e-10
cat("=> Keputusan LM:", keputusan, "\n")
## => Keputusan LM: SAR
q <- 2 * ncol(x) + 1 # jumlah kolom instrumen
# ---- 7.1 Bandwidth dan estimator ----
# (a) Awal: WLS
bw.sar <- bwd.gwr.sar(y, x, coord, W_lag, kernel, "wls")
estimator <- "wls"
ESS_wls <- hitung.ESS(D, bw.sar$bw, kernel)
cat(sprintf("WLS: h = %.2f, min ESS = %.2f, q = %d\n", bw.sar$bw, min(ESS_wls), q))
## WLS: h = 69670.95, min ESS = 1.36, q = 11
# (b) 2SLS dipakai hanya jika ESS minimum pada bandwidth BARU masih > q
if (min(ESS_wls) > q) {
bw2 <- bwd.gwr.sar(y, x, coord, W_lag, kernel, "2sls")
ESS2 <- hitung.ESS(D, bw2$bw, kernel)
cat(sprintf("2SLS: h = %.2f, min ESS = %.2f\n", bw2$bw, min(ESS2)))
if (min(ESS2) > q) {
bw.sar <- bw2
estimator <- "2sls"
} else {
message("ESS 2SLS <= q pada bandwidth baru -> tetap WLS")
}
}
# (c) Cek batas grid
if (bw.sar$idx_grid %in% c(1, nrow(bw.sar$grid)))
warning("Optimum jatuh di batas grid: perlebar rentang bandwidth.")
# ---- 7.2 Estimasi akhir ----
fit <- prog.gwr.sar(y, x, coord, W_lag, bw.sar$bw, kernel, estimator)
cat(sprintf("=> Estimator = %s, h optimum = %.2f (CV = %.5f)\n",
fit$estimator, bw.sar$bw, bw.sar$cv))
## => Estimator = wls, h optimum = 69670.95 (CV = 1.38291)
# ---- 7.3 Parameter lokal ----
Estimasi_GWRSAR <- data.frame(coord, fit$beta)
print(summary(Estimasi_GWRSAR))
## Longitude Latitude Rho Intercept
## Min. :-286460 Min. :-504107 Min. :-0.6308 Min. :-0.62847
## 1st Qu.:-176467 1st Qu.:-270339 1st Qu.: 0.4194 1st Qu.: 0.08699
## Median : 8909 Median :-104925 Median : 0.6130 Median : 0.47627
## Mean : 25406 Mean : -26631 Mean : 0.5800 Mean : 0.79280
## 3rd Qu.: 172975 3rd Qu.: 281036 3rd Qu.: 0.7739 3rd Qu.: 1.49696
## Max. : 588654 Max. : 687412 Max. : 1.2689 Max. : 3.15215
## RLS TPAKP TPAKL
## Min. :-0.116678 Min. :-0.009176 Min. :-0.038878
## 1st Qu.:-0.040113 1st Qu.:-0.001596 1st Qu.:-0.011760
## Median : 0.017348 Median : 0.002685 Median :-0.002065
## Mean : 0.005107 Mean : 0.002432 Mean :-0.001311
## 3rd Qu.: 0.041101 3rd Qu.: 0.006074 3rd Qu.: 0.013714
## Max. : 0.171802 Max. : 0.011913 Max. : 0.029306
## TPT MISKIN....
## Min. :-0.133296 Min. :-0.01798
## 1st Qu.:-0.051126 1st Qu.: 0.01001
## Median :-0.026650 Median : 0.01482
## Mean :-0.030466 Mean : 0.01890
## 3rd Qu.:-0.005126 3rd Qu.: 0.02578
## Max. : 0.077838 Max. : 0.05975
Rho_lokal <- fit$beta[, "Rho"]
cat(sprintf("Proporsi |rho_i| >= 1 (melanggar stasioneritas): %.3f\n",
mean(abs(Rho_lokal) >= 1)))
## Proporsi |rho_i| >= 1 (melanggar stasioneritas): 0.099
cat(sprintf("Proporsi lokasi signifikan (alpha = %g):\n", alpha))
## Proporsi lokasi signifikan (alpha = 0.05):
print(round(colMeans(fit$pval < alpha), 3))
## Rho Intercept RLS TPAKP TPAKL TPT MISKIN....
## 0.556 0.099 0.062 0.210 0.173 0.198 0.259
# Jika Anda tetap ingin melihat Uji t Lokal per wilayah:
uji_sig <- uji.signifikansi(fit, id = data[, c("prov", "kabkot")], alpha = alpha)
print(uji_sig$ringkasan, digits = 4)
## Parameter Min Median Mean Max Mean_t N_Signif
## 1 Rho -0.630806 0.613029 0.579957 1.26887 1.94081 45
## 2 Intercept -0.628469 0.476271 0.792801 3.15215 0.97684 8
## 3 RLS -0.116678 0.017348 0.005107 0.17180 -0.04305 5
## 4 TPAKP -0.009176 0.002685 0.002432 0.01191 0.81033 17
## 5 TPAKL -0.038878 -0.002065 -0.001311 0.02931 -0.17833 14
## 6 TPT -0.133296 -0.026650 -0.030466 0.07784 -0.97919 16
## 7 MISKIN.... -0.017978 0.014823 0.018901 0.05975 1.22429 21
## Prop_Signif N_Positif N_Negatif Status
## 1 0.55556 45 0 Signifikan di sebagian lokasi
## 2 0.09877 8 0 Signifikan di sebagian lokasi
## 3 0.06173 0 5 Signifikan di sebagian lokasi
## 4 0.20988 17 0 Signifikan di sebagian lokasi
## 5 0.17284 3 11 Signifikan di sebagian lokasi
## 6 0.19753 0 16 Signifikan di sebagian lokasi
## 7 0.25926 21 0 Signifikan di sebagian lokasi
# Buat workbook baru
wb <- createWorkbook()
# Tambahkan sheet untuk Ringkasan
addWorksheet(wb, sheetName = "Ringkasan")
writeData(wb, sheet = "Ringkasan", x = uji_sig$ringkasan)
# Tambahkan sheet untuk Detail Lokal
addWorksheet(wb, sheetName = "Lokal")
writeData(wb, sheet = "Lokal", x = uji_sig$lokal)
# Simpan file Excel
saveWorkbook(wb, file = "Uji_Sig_Parameter.xlsx", overwrite = TRUE)
cat("File 'Uji_Sig_Parameter.xlsx' berhasil disimpan.\n")
## File 'Uji_Sig_Parameter.xlsx' berhasil disimpan.
gwr0 <- prog.gwr.sar(y, x, coord, W_lag, h_gwr, kernel, lag = FALSE) # GWR tanpa lag (AICc sebanding)
k_ols <- ncol(x) + 1
mse_ols <- mean(residuals(OLS)^2)
mse_gwr <- mean(gwr0$residuals^2)
mse_sar <- mean(fit$residuals^2)
perbandingan <- data.frame(
Model = c("OLS", "GWR", "GWR-SAR"),
R2 = c(summary(OLS)$r.squared, gwr0$R2, fit$R2),
R2adj = c(summary(OLS)$adj.r.squared, gwr0$R2adj, fit$R2adj),
MSE = c(mse_ols, mse_gwr, mse_sar),
RMSE = sqrt(c(mse_ols, mse_gwr, mse_sar)),
AICc = c(aicc.fun(sum(residuals(OLS)^2), n - k_ols, k_ols, n), gwr0$AICc, fit$AICc))
print(perbandingan, digits = 5)
## Model R2 R2adj MSE RMSE AICc
## 1 OLS 0.43271 0.39489 0.0297525 0.172489 -33.066
## 2 GWR 0.95732 0.74143 0.0022382 0.047310 299.583
## 3 GWR-SAR 0.94098 0.73740 0.0030956 0.055638 189.377