Pemodelan Geographically Weigthed Regression (GWR Model) menggunakan R

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!