GWR SAR
teguh
2026-10-11
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
# ---- 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)1. DATA
data <- data.frame(read_excel(file_data, sheet = sheet_data))
y <- data[, 5] # respon: TFR
x <- data[, c(6, 7, 9, 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 0.334 0.299 0.316 0.334 0.314 ...
## $ 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 ...2. OLS
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.055565 -0.021499 -0.001962 0.016526 0.093182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2364334 0.0465002 5.085 2.57e-06 ***
## RLS 0.0024485 0.0044653 0.548 0.585066
## TPAKP 0.0012863 0.0005249 2.450 0.016567 *
## TPT -0.0048266 0.0032404 -1.489 0.140498
## MISKIN.... 0.0049704 0.0013841 3.591 0.000581 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03276 on 76 degrees of freedom
## Multiple R-squared: 0.4186, Adjusted R-squared: 0.388
## F-statistic: 13.68 on 4 and 76 DF, p-value: 1.894e-083. BANDWIDTH GWR (PEMILIHAN KERNEL)
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: 0.08278404
## Fixed bandwidth: 564865.3 CV score: 0.07721869
## Fixed bandwidth: 349218.9 CV score: 0.06918309
## Fixed bandwidth: 215942.1 CV score: 0.05921385
## Fixed bandwidth: 133572.5 CV score: 0.05017147
## Fixed bandwidth: 82665.27 CV score: 0.04482401
## Fixed bandwidth: 51202.88 CV score: 0.04320759
## Fixed bandwidth: 31758.06 CV score: 0.04866014
## Fixed bandwidth: 63220.45 CV score: 0.04341538
## Fixed bandwidth: 43775.62 CV score: 0.04369163
## Fixed bandwidth: 55793.18 CV score: 0.04318265
## Fixed bandwidth: 913788.5 CV score: 0.08883346
## Fixed bandwidth: 564865.3 CV score: 0.08346086
## Fixed bandwidth: 349218.9 CV score: 0.0721317
## Fixed bandwidth: 215942.1 CV score: 0.05790676
## Fixed bandwidth: 133572.5 CV score: 0.04626845
## Fixed bandwidth: 82665.27 CV score: 0.05214777
## Fixed bandwidth: 165034.9 CV score: 0.05114479
## Fixed bandwidth: 114127.7 CV score: 0.04366026
## Fixed bandwidth: 102110.1 CV score: 0.04555428
## Fixed bandwidth: 121554.9 CV score: 0.04456849
## Fixed bandwidth: 109537.4 CV score: 0.04493185
## Fixed bandwidth: 116964.6 CV score: 0.04399161
## Fixed bandwidth: 112374.3 CV score: 0.04346589
## Fixed bandwidth: 111290.7 CV score: 0.0433499
## Fixed bandwidth: 110621 CV score: 0.04490323
## Fixed bandwidth: 111704.6 CV score: 0.04339382
## Fixed bandwidth: 913788.5 CV score: 0.07517516
## Fixed bandwidth: 564865.3 CV score: 0.06106176
## Fixed bandwidth: 349218.9 CV score: 0.05044693
## Fixed bandwidth: 215942.1 CV score: 0.05595139
## Fixed bandwidth: 431588.5 CV score: 0.05306428
## Fixed bandwidth: 298311.7 CV score: 0.04744146
## Fixed bandwidth: 266849.3 CV score: 0.05838481
## Fixed bandwidth: 317756.5 CV score: 0.04830623
## Fixed bandwidth: 286294.1 CV score: 0.04677532
## Fixed bandwidth: 278866.9 CV score: 0.05560715
## Fixed bandwidth: 290884.4 CV score: 0.04708051
## Fixed bandwidth: 283457.2 CV score: 0.05493043
## Fixed bandwidth: 288047.5 CV score: 0.04689615
## Fixed bandwidth: 285210.5 CV score: 0.05472604
## Fixed bandwidth: 286963.8 CV score: 0.04682168
## Fixed bandwidth: 913788.5 CV score: 0.07731267
## Fixed bandwidth: 564865.3 CV score: 0.06287393
## Fixed bandwidth: 349218.9 CV score: 0.05055393
## Fixed bandwidth: 215942.1 CV score: 0.06188474
## Fixed bandwidth: 431588.5 CV score: 0.0537769
## Fixed bandwidth: 298311.7 CV score: 0.04713675
## Fixed bandwidth: 266849.3 CV score: 0.0592232
## Fixed bandwidth: 317756.5 CV score: 0.04834538
## Fixed bandwidth: 286294.1 CV score: 0.05593707
## Fixed bandwidth: 305738.9 CV score: 0.04762516
## Fixed bandwidth: 293721.4 CV score: 0.05522601
## Fixed bandwidth: 301148.6 CV score: 0.04733994
## Fixed bandwidth: 296558.3 CV score: 0.0550267
## Fixed bandwidth: 299395.3 CV score: 0.04721736bw_kernel## exponential gaussian bisquare tricube
## 55793.18 111290.70 286294.12 298311.68bandwith.GWR <- unname(bw_kernel[1]) # bandwidth kernel terpilih
h_gwr <- bandwith.GWR4. PENGUJIAN ASUMSI
# ---- 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.085734, p-value = 0.5616
## alternative hypothesis: two-sidedprint(shapiro.test(sisaan_ols)) # lebih tepat untuk n kecil dan parameter diduga##
## Shapiro-Wilk normality test
##
## data: sisaan_ols
## W = 0.96141, p-value = 0.01546# ---- 4.2 Heterogenitas spasial (Breusch-Pagan) ----
BPTest_Spasial <- bptest(OLS, studentize = FALSE)
BPTest_Spasial##
## Breusch-Pagan test
##
## data: OLS
## BP = 10.014, df = 4, p-value = 0.04018# ---- 4.3 Multikolinearitas (VIF) ----
Multi <- vif(OLS)
Multi## RLS TPAKP TPT MISKIN....
## 1.833854 1.354483 2.206433 1.605330# ---- 4.4 Korelasi antar variabel ----
Korelasi <- corrplot(cor(data.frame(TFR = y, x)), method = "number")
5. GWR
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
## Min. :-0.04657 Min. :-0.022236 Min. :-2.057e-03
## 1st Qu.: 0.19580 1st Qu.:-0.006839 1st Qu.: 1.688e-06
## Median : 0.28684 Median :-0.001852 Median : 6.967e-04
## Mean : 0.28640 Mean : 0.002150 Mean : 5.500e-04
## 3rd Qu.: 0.37988 3rd Qu.: 0.011345 3rd Qu.: 1.157e-03
## Max. : 0.55817 Max. : 0.035220 Max. : 2.709e-03
## X4 X5
## Min. :-0.0193255 Min. :-0.004200
## 1st Qu.:-0.0103342 1st Qu.: 0.001343
## Median :-0.0062097 Median : 0.003205
## Mean :-0.0040757 Mean : 0.004084
## 3rd Qu.: 0.0008578 3rd Qu.: 0.006443
## Max. : 0.0220731 Max. : 0.012530# ---- 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)
# ---- 5.3 Kebaikan model ----
R.Square.Adj.GWR <- gwr$R_Square_Adj.GWR
R.Square.JKR.GWR <- gwr$R_Square.GWR
AIC.GWR <- gwr$AIC6. MATRIKS PEMBOBOT SPASIAL W, MORAN’S I, DAN UJI LM
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 TFR ----
lw <- mat2listw(W_lag, style = "W")
moran_tfr <- moran.test(data$TFR, lw, randomisation = TRUE, alternative = "two.sided")
print(moran_tfr)##
## Moran I test under randomisation
##
## data: data$TFR
## weights: lw
##
## Moran I statistic standard deviate = 11.358, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.530026989 -0.012500000 0.002281404# ---- 6.3 Uji LM spasial pada OLS global ----
hasil_lm <- lm.spatial.tests(OLS, W_lag)
keputusan <- pilih.model.lm(hasil_lm, alpha = alpha)
print(hasil_lm, digits = 5)## Uji Statistik df p.value
## 1 LMerr 49.1964 1 2.3157e-12
## 2 LMlag 61.0042 1 5.6952e-15
## 3 RLMerr 2.7745 1 9.5775e-02
## 4 RLMlag 14.5823 1 1.3417e-04cat("=> Keputusan LM:", keputusan, "\n")## => Keputusan LM: SAR7. GWR-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 = 68328.26, min ESS = 1.33, q = 9# (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 = 68328.26 (CV = 0.04082)# ---- 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.5498 Min. :-0.09014
## 1st Qu.:-176467 1st Qu.:-270339 1st Qu.: 0.1269 1st Qu.: 0.01529
## Median : 8909 Median :-104925 Median : 0.6269 Median : 0.08106
## Mean : 25406 Mean : -26631 Mean : 0.5018 Mean : 0.13925
## 3rd Qu.: 172975 3rd Qu.: 281036 3rd Qu.: 0.7343 3rd Qu.: 0.33472
## Max. : 588654 Max. : 687412 Max. : 1.2580 Max. : 0.45962
## RLS TPAKP TPT
## Min. :-0.0189923 Min. :-2.270e-03 Min. :-0.0170834
## 1st Qu.:-0.0086700 1st Qu.: 6.494e-05 1st Qu.:-0.0086037
## Median :-0.0036049 Median : 5.840e-04 Median :-0.0061715
## Mean :-0.0007472 Mean : 5.039e-04 Mean :-0.0045891
## 3rd Qu.: 0.0075244 3rd Qu.: 9.933e-04 3rd Qu.:-0.0002186
## Max. : 0.0268045 Max. : 2.260e-03 Max. : 0.0159560
## MISKIN....
## Min. :-0.003606
## 1st Qu.: 0.001441
## Median : 0.002627
## Mean : 0.003141
## 3rd Qu.: 0.004714
## Max. : 0.010531Rho_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.074cat(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 TPT MISKIN....
## 0.741 0.148 0.049 0.148 0.185 0.185# ---- 7.5 Uji kebaikan model (F-test) ----
print(round(fit$F_test, 5))## F df1 df2 p
## 0.80167 39.39985 75.00000 0.22596# ---- 7.6 Uji signifikansi parameter lokal (Uji F2) ----
# Pastikan fit sudah di-generate dengan fungsi yang baru dimodifikasi
hasil_f2 <- uji_f2_spasial(fit = fit, alpha = alpha)
print(hasil_f2)## Variabel F2_Statistik df1 df2 p_value Status
## 1 Rho 5.53 6.71 39.4 < 0.001 Signifikan
## 2 Intercept 7.41 7.00 39.4 < 0.001 Signifikan
## 3 RLS 5.63 7.92 39.4 < 0.001 Signifikan
## 4 TPAKP 2.79 7.80 39.4 0.01569 Signifikan
## 5 TPT 5.66 9.08 39.4 < 0.001 Signifikan
## 6 MISKIN.... 5.10 7.32 39.4 < 0.001 Signifikan8. PERBANDINGAN MODEL: OLS vs GWR vs GWR-SAR
gwr0 <- prog.gwr.sar(y, x, coord, W_lag, h_gwr, kernel, lag = FALSE) # GWR tanpa lag (AICc sebanding)
k_ols <- ncol(x) + 1
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),
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 AICc
## 1 OLS 0.41863 0.38803 -310.797
## 2 GWR 0.94375 0.73446 -62.279
## 3 GWR-SAR 0.92806 0.72290 -125.588