 |
| Geographically Weigthed Regression (GWR) dengan R |
Halo teman-teman, kemarin kita telah sejenak break dari pemodelan dengan membahas fungsi sprintf() yang ada di R. Kali ini kita akan melanjutkan kembali pembahasan model-model statistik dengan R. Yang akan kita bahas kali ini adalah bagaimana menerapkan model Geographically Weigthed Regression atau biasa disingkat model GWR.
GWR merupakan sebuah model alternatif bila di dalam pemodelan regresi linier (dalam parameter) kita terganggu oleh asumsi heteroskedastisitas spasial atau asumsi homoskedastisitas residual model tidak terpenuhi. Mengapa disebut heteroskedastisitas spasial? Ya karena di dalam data kita terdiri atas amatan-amatan lokasi yang masing-masing memiliki variabel pembobot spasial dan di dalam praktiknya masing-masing lokasi ini ditunjukkan oleh posisi astronomis longitute dan latitude.
Kalau di model regresi Robust yang sebelumnya kita bahas, ia merupakan model alternatif bila residual modelnya terganggu asumsi normalitas dan atau heteroskedastisitas, tapi setiap amatannya bukan terdiri atas lokasi-lokasi amatan. Bila ternyata di dalam data kita memiliki lokasi-lokasi amatan yang ditunjukkan oleh longitude dan latitude, maka model GWR bisa juga menjadi alternatif model untuk analisis.
Baik, kita kembali ke model GWR. Bila di dalam model regresi dengan metode estimasi parameter Ordinary Least Square (OLS) hanya melihat pengaruh variabel independen terhadap variabel dependennya baik secara parsial maupun simultan, maka di model GWR kita dapat melihat seberapa besar pengaruh variabel independen terhadap variabel dependen di sebuah wilayah atau lokasi amatan dan seberapa besar pengaruh variabel independen seluruh lokasi terhadap variabel dependennya. Kalau diistilahkan, di dalam model GWR ini kita akan mengenal model lokal dan model global. Artinya, pada masing-masing lokasi atau wilayah akan memiliki model regresi GWR yang berbeda dan secara global juga memiliki sebuah model GWR.
Secara mendasar, karakteristik dari model GWR ini ditunjukkan oleh matriks pembobot dengan menggunakan beberapa teknik dengan menggunakan fungsi Kernel, yaitu fixed kernel Gaussian, adaptive kernel Gaussian, fixed kernel Bisquare, adaptive kernel Bisquare atau lainnya. Menurut Fotheringham et.al (2002), fungsi matriks pembobot ini memiliki bandwith (h) yang berbeda untuk setiap lokasi. Bandwith ini sebenarnya merupakan parameter penghalus non-negatif berbentuk lingkaran yang memiliki radius sebesar b dari titik pusat lokasi sebagai dasar penentuan bobot model regresi di setiap lokasi amatan. Dengan demikian, untuk pengamatan yang dekat dengan lokasi atau wilayah X misalkan, maka ia akan lebih berpengaruh terhadap pembentukan parameter lokasi atau wilayah X tersebut.
Dalam kasus model GWR, nilai bandwith yang kecil menyebabkan varians membesar sehingga diperlukan metode-metode bagaimana memilih bandwith sedemikian rupa sehingga mampu memperkecil varians. Sebab varians yang besar (tidak homogen) menyebabkan penduga parameter yang dihasilkan tidak akurat. Adapun metode pemilihan bandwith yang optimum ini ada beberapa, yaitu dengan berdasarkan nilai Cross Validation (CV), Akaike Information Criteria (AIC), General Cross Validation (GCV), atau Bayessian Information Criteria (BIC).
Secara matematis, bentuk model GWR secara umum dapat dituliskan sebagai berikut:
 |
| Bentuk umum model Geographically Weigthed Regression (GWR) |
Baik, itu sekilas teori mengenai GWR teman-teman. Selanjutnya kita akan mempraktikkan bagaimana memodelkan GWR dengan R. Data yang akan kita gunakan dalam praktik kali ini adalah data mengenai jumlah kematian bayi, indeks ketersediaan pangan, dan jumlah persalinan yang dibantu tenaga kesehatan. Data tersebut bersumber dari Profil Kesehatan Jawa Timur dan kemeterian Pertanian tahun 2020 yang dapat teman-teman unduh pada tautan berikut. Setelah datanya telah siap, pemodelan GWR dapat mengikuti beberapa code berikut:
Code:
#Mengimport Data
library(readxl)
dataku <- read_excel("C:/Users/Joko Ade/Downloads/datagwr.xlsx")
#Sesuaikan dengan letak penyimpanan file
#melihat sekilas struktur data
str(dataku)
Hasil:
tibble [38 x 6] (S3: tbl_df/tbl/data.frame)
$ Wilayah : chr [1:38] "Pacitan" "Ponorogo" "Trenggalek" "Tulungagung" ...
$ long : num [1:38] 111 111 112 112 112 ...
$ lat : num [1:38] -8.13 -7.87 -8.18 -8.91 -8.1 ...
$ JKB : num [1:38] 62 139 38 158 115 163 86 174 338 117 ...
$ IP : num [1:38] 74.1 76.2 71.6 76.5 76.1 ...
$ Salinakes: num [1:38] 6.55 10.62 8.99 14.3 15.52 ...
Code:
#Melakukan Seleksi data dan melihat ringkasan data
dataku <- dataku[,c(-1)]
summary(dataku)
#Mengattach data
attach(dataku)
Hasil:
long lat JKB IP Salinakes
Min. :111.1 Min. :-8.912 Min. : 6.0 Min. :58.40 Min. : 2.103
1st Qu.:112.0 1st Qu.:-7.971 1st Qu.: 44.5 1st Qu.:67.11 1st Qu.: 8.729
Median :112.5 Median :-7.712 Median : 96.0 Median :75.11 Median :14.629
Mean :112.6 Mean :-7.691 Mean :101.8 Mean :73.58 Mean :15.062
3rd Qu.:113.1 3rd Qu.:-7.463 3rd Qu.:148.8 3rd Qu.:79.41 3rd Qu.:17.820
Max. :114.4 Max. :-6.895 Max. :338.0 Max. :90.11 Max. :43.534
Code:
#Pemodelan regresi OLS
ols <- lm(JKB~IP+Salinakes)
#Ringkasan model OLS
summary(ols)
#Melihat AIC
AIC(ols)
Hasil:
Call:
lm(formula = JKB ~ IP + Salinakes)
Residuals:
Min 1Q Median 3Q Max
-123.989 -32.292 1.592 29.144 105.588
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 210.3187 76.3647 2.754 0.00927 **
IP -2.5151 1.0633 -2.365 0.02368 *
Salinakes 5.0789 0.8212 6.184 4.42e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 49.12 on 35 degrees of freedom
Multiple R-squared: 0.5261, Adjusted R-squared: 0.499
F-statistic: 19.42 on 2 and 35 DF, p-value: 2.114e-06
AIC: [1] 408.6804
Code:
#Uji Asumsi Normalitas Residual
ks.test(ols$residuals, ecdf(ols$residuals))
shapiro.test(residuals(ols))
#Uji Non Autokorelasi Residual
dwtest(ols)
#Uji Non Multikolinieritas Variabel Bebas Model
ols_vif_tol(ols)
#Uji Homoskedastisitas Residual Model
bptest(ols)
Hasil:
One-sample Kolmogorov-Smirnov test
data: ols$residuals
D = 0.026316, p-value = 1
alternative hypothesis: two-sided
Shapiro-Wilk normality test
data: residuals(ols)
W = 0.99049, p-value = 0.9836
Durbin-Watson test
data: ols
DW = 1.6189, p-value = 0.09199
alternative hypothesis: true autocorrelation is greater than 0
Variables Tolerance VIF
1 IP 0.9317424 1.073258
2 Salinakes 0.9317424 1.073258
studentized Breusch-Pagan test
data: ols
BP = 11.072, df = 2, p-value = 0.003943
Interpretasi: terlihat bahwa model regresi memenuhi semua uji asumsi klasik kecuali uji homoskedastisitas karena nilai p-value Breush-Pagan 0,0003943 atau < 0,05 (alpha penelitian), maka model GWR dapat menjadi alternatif bila setiap amatan kita adalah lokasi dengan memiliki data longitude dan latitude
Code:
#Karena Terganggu Heteroskedastisitas Opsi Lain Bisa Menggunakan GWR
#Install dan aktivasi beberapa package untuk pemodelan GWR
install.packages("spgwr")
install.packages("sp")
install.packages("spData")
library(spgwr)
library(sp)
library(spData)
#Pemodelan GWR
#Mencari Bandwith, misalkan Bisquare Adaptive
band_bidap <- gwr.sel(formula = JKB~IP+Salinakes, data = dataku,
coords = cbind(dataku$long, dataku$lat), adapt = TRUE, gweight = gwr.bisquare)
Hasil:
Adaptive q: 0.381966 CV score: 71986.42
Adaptive q: 0.618034 CV score: 85144.32
Adaptive q: 0.236068 CV score: 79166.29
Adaptive q: 0.3985704 CV score: 72813.99
Adaptive q: 0.3493846 CV score: 71763.77
Adaptive q: 0.3061015 CV score: 78009.31
Adaptive q: 0.3617677 CV score: 71361.09
Adaptive q: 0.3639216 CV score: 71345.02
Adaptive q: 0.3650074 CV score: 71341.32
Adaptive q: 0.3658301 CV score: 71340.32
Adaptive q: 0.3659504 CV score: 71340.29
Adaptive q: 0.3659911 CV score: 71340.29
Adaptive q: 0.3660318 CV score: 71340.3
Adaptive q: 0.3659911 CV score: 71340.29 ----> nilai bandwith yang digunakan 0,3659911
Code:
#pemodelan GWR
modelfix <- gwr(formula = JKB~IP+Salinakes, data = dataku,
coords = cbind(dataku$long, dataku$lat), bandwidth = band_bidap,
hatmatrix = TRUE)
#Pengecekan koefisien model
modelfix$lm$coefficients
#Mengecek signifikansi variabel dalam model
LMZ.F3GWR.test(modelfix)
Hasil:
(Intercept) IP Salinakes
210.318714 -2.515098 5.078869
Leung et al. (2000) F(3) test
F statistic Numerator d.f. Denominator d.f. Pr(>)
(Intercept) 5.9564 13.6357 24.734 6.602e-05 ***
IP 5.7204 12.0966 24.734 0.0001238 ***
Salinakes 3.4095 13.6284 24.734 0.0039280 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Interpretasi: terlihat seluruh variabel independen berpengaruh signifikan statistik terhadap variabel dependen
Code:
#Melihat AIC
modelfix$results[8]
Hasil:
$AICh
[1] 355.3489
Code:
#cek Anova modelnya perbandingan OLS dan GWR
#bila signifikan berarti model GWR lebih baik dari OLS
anova(modelfix)
#mencari nilai F tabel
qf(p = 0.05, df1 = anova(modelfix)[2,1], df2 = anova(modelfix)[3,1],
lower.tail = FALSE)
#Jika nilai |F hitung| > F tabel tolak H0 artinya GWR lebih baik dari OLS
# Atau Mencari p-value nilai F Anova, jika < 0.05 tolak H0 artinya GWR lebih baik dari OLS
pf(q = anova(modelfix)[3,4], df1 = anova(modelfix)[2,1],
df2 = anova(modelfix)[3,1], lower.tail = FALSE)
Hasil:
Analysis of Variance Table
Df Sum Sq Mean Sq F value
OLS Residuals 3.000 84451
GWR Improvement 16.169 67176 4154.7
GWR Residuals 18.831 17275 917.4 4.5289
Nilai F-Tabel: [1] 2.217443
Nilai p-value: [1] 0.00118195
Code:
#Selain itu uji model Global dapat dilakukan dengan uji F LMZ
LMZ.F1GWR.test(modelfix)
Hasil:
Leung et al. (2000) F(1) test
data: modelfix
F = 0.38019, df1 = 24.734, df2 = 35.000, p-value = 0.007164
alternative hypothesis: less
sample estimates:
SS OLS residuals SS GWR residuals
84451.27 17274.95
Code:
#Melihat Intersep model lokal
modelfix$SDF$`(Intercept)`
#melihat Koefisien IP lokal
modelfix$SDF$IP
Hasil:
[1] -165.493877 -109.437424 -158.089894 -1046.972846 -67.279399 20.705466 212.146159 483.523836
[9] 555.629387 836.886416 552.258880 659.191891 487.971285 158.777521 137.298122 120.523552
[17] 111.351302 6.488024 -57.322304 -62.434970 -23.977988 68.010960 395.432418 244.031644
[25] 201.522895 150.840227 124.501058 259.844691 353.003696 -26.771901 -67.279399 86.833596
[33] 343.234161 156.035701 146.394134 -66.493632 162.328836 58.899471
[1] 2.00291333 1.42004407 2.23550425 16.07065781 1.33894888 -0.04054907 -2.36473580 -6.71981226
[9] -7.97684097 -11.37571776 -7.54229252 -8.63917298 -7.19002370 -1.92665900 -1.74410450 -1.46782969
[17] -1.33347504 0.06484369 0.80541727 0.76526999 0.29533918 -0.72940634 -4.83500655 -3.08973403
[25] -2.62090249 -2.14313147 -2.06405256 -4.44891112 -5.27246568 0.56831402 1.33894888 -0.77049139
[33] -4.85706069 -1.97571435 -1.83141695 0.88599258 -2.16050614 -0.45529833
Oke, demikian sedikit pembahasan mengenai penerapan model Geographically Weigthed Regression (GWR) dengan R. Jangan lupa share dan bisa komentar bila ada pertanyaan di kolom komentar. Dan pastinya, simak dan ikuti terus unggahan-unggahan berikutnya dalam blog ini. Selamat mempraktikkan!