Anggota Kelompok 1
- G1401201004 Alfidhia Rahman Nasa Juhanda
- G1401201050 Hanung Safrizal
- G1401201064 Muhammad Irsyad Robbani
- G1401201081 Zakiyah Azzahra
- G1401201092 Muhammad Zhillan Zakiyyan
lapply(c("kableExtra","tidyverse","plm","tinytex","tseries","lmtest", "broom", "stargazer", "RColorBrewer", "googlesheets4", "corrplot", "car"), library, character.only = T)[[1]]## [1] "kableExtra" "stats" "graphics" "grDevices" "utils"
## [6] "datasets" "methods" "base"Data
Data yang digunakan adalah data Jumlah Tindak Pidana per 100 ribu penduduk dengan 7 Peubah bebas. Data ini berjenis data panel dengan individunya adalah 34 Provinsi di Indonesia selama 6 tahun (Tahun 2015-2020) dengan peubah-peubah sebagai berikut:
PDN : Jumlah tindak pidana per 100.000 Penduduk (Peubah Respons)
KMIS : Tingkat crimekinan (%)
RLS : Rata-rata Lama Sekolah (Tahun)
IPM : Indeks Pembangunan Manusia (Poin)
TPT : Tingkat Pengangguran Terbuka (%)
PDRB : Produk Domestik Regional Bruto Per Kapita (Ribu Rupiah)
AHH : Angka Harapan Hidup (Tahun)
PKD : Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun)
Selain itu, digunakan pula peubah Jumlah Penduduk dan Luas Wilayah untuk manipulasi data. Untuk lebih jelasnya, dapat dilihat dari output berikut:
gs4_deauth()
crime <- read_sheet("https://docs.google.com/spreadsheets/d/1nJUZTW7l52XB77-zQpurDjVw0tYyTWCF141CCRZ27yg/edit#gid=1513708653")## ✔ Reading from "Kelompok 1 ADP".## ✔ Range 'table'.colnames(crime) <- c("PROV", "TAHUN", "PENDUDUK", "LUAS", "PDN", "KMIS", "RLS", "IPM", "TPT", "PDRB", "AHH", "PKD")
crime$PROV <- as.factor(crime$PROV)
crime$TAHUN <- as.factor(crime$TAHUN)
str(crime)## tibble [204 × 12] (S3: tbl_df/tbl/data.frame)
## $ PROV : Factor w/ 34 levels "ACEH","BALI",..: 1 34 32 26 8 33 4 19 17 18 ...
## $ TAHUN : Factor w/ 6 levels "2015","2016",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ PENDUDUK: num [1:204] 5388 14798 5546 6951 3604 ...
## $ LUAS : num [1:204] 57956 72981 42013 87024 50058 ...
## $ PDN : num [1:204] 149 231 150 123 134 147 172 92 133 127 ...
## $ KMIS : num [1:204] 15.21 8.95 6.42 6.93 7.78 ...
## $ RLS : num [1:204] 9.33 9.54 8.99 9.14 8.55 ...
## $ IPM : num [1:204] 72 71.8 72.4 72.7 71.3 ...
## $ TPT : num [1:204] 6 5.81 6.07 5.62 4.7 4.71 3.58 4.47 4.3 8.16 ...
## $ PDRB : num [1:204] 31633 54979 43826 114167 57958 ...
## $ AHH : num [1:204] 70 69.2 69.5 71.7 71.2 ...
## $ PKD : num [1:204] 9492 10420 10733 10675 10392 ...Dapat terlihat bahwa terdapat 204 total observasi, yaitu 34 individu dengan periode 6 tahun. Semua peubah berjenis numerik kontinu.
Pre-processing Data
Penanganan Missing Value
colSums(is.na(crime))## PROV TAHUN PENDUDUK LUAS PDN KMIS RLS IPM
## 0 0 0 0 5 0 0 0
## TPT PDRB AHH PKD
## 0 0 0 0Dari output di atas, dapat terlihat bahwa terdapat 5 missing value di peubah respons, yaitu PDN. kita akan menggunakan imputasi Group Mean Imputation untuk mengatasi masalah ini dengan menjadikan setiap provinsi adalah grup. Sehingga, kita perlu tahu provinsi apa saja yang memiliki missing value.
crime %>%
select(PROV, TAHUN, PDN) %>%
filter(is.na(PDN))num <- crime %>% mutate(no = row_number()) %>% filter(is.na(PDN)) %>% pull(no)Dapat terlihat bahwa provinsi Kalimantan Utara dan Sulawesi Barat merupakan provinsi yang memiliki nilai missing value, sehingga kita perlu menghitung rataan nilai PDN dari kedua provinsi tersebut.
crime %>%
na.omit() %>%
select(PROV, PDN) %>%
group_by(PROV) %>%
summarise("Rata-rata"=mean(PDN)) %>%
filter(PROV %in% c("KALIMANTAN UTARA", "SULAWESI BARAT"))Kemudian, input rata-rata yang didapatkan di atas terhadap nilai yang hilang sebelumnya.
crime <- crime %>%
mutate(PDN = ifelse(is.na(PDN),
ifelse(PROV == "KALIMANTAN UTARA",
110, ifelse(PROV == "SULAWESI BARAT",
140.75,PDN)), PDN))Sehingga didapat:
crime %>%
select(PROV, TAHUN, PDN) %>%
slice(as.numeric(num))Masalah missing value terselesaikan.
Feature Engineering
Data yang diperoleh pada sumber adalah Jumlah tindak pidana yang telah terboboti oleh Jumlah penduduk setiap provinsi. Analisis pada peubah ini telah dilakukan sebelumnya oleh Febrianti (2022). Pada penelitian ini, kami mengusung pembobotan tambahan pada peubah yang sama, yaitu dengan pembobotan pada luas wilayah setiap provinsi. Syntax yang digunakan adalah sebagai berikut:
crime <- crime %>%
mutate(PDN_AREA=PDN/LUAS*100) %>%
select(PROV,TAHUN,PDN,PDN_AREA,everything(),-PENDUDUK,-LUAS)
crimeMetodologi
Metodologi yang digunakan adalah sebagai berikut:
Melakukan pre-pocessing data, yaitu menanganani missing value dan melakukan feature engineering dengan pembobotan pada luas wilayah di setiap provinsi pada data jumlah tindak pidana.
Melakukan eksplorasi data.
Melakukan pemodelan regresi panel CEM, FEM, dan REM.
Melakukan uji Chow, Uji Hausman, dan Uji Lagrange Multiplier untuk memilih model panel terbaik di antara model CEM, FEM, dan REM.
Memeriksa heterogenitas spasial dengan uji Breusch Pagan pada model terpilih. Pemodelan GWPR mengasumsikan adanya heterogenitas spasial.
Menghitung bandwidth optimum dan matriks pembobot spasial.
Mengestimasi nilai parameter pemodelan GWPR.
Melakukan pengujian serentak pemodelan GWPR dengan statistik uji F dan pengujian parsial dengan statistik uji t
Melakukan perbandingan model antara regresi panel dan GWPR dengan kriteria R^2.
Interpretasi model terbaik yang diperoleh.
Eksplorasi Data
Keragaman Antar Waktu
ggplot(data = crime, aes(x = TAHUN, y = PDN_AREA)) +
geom_line() +
labs(x = "Tahun", y = "Jumlah Pidana") +
theme(legend.position = "none") +
theme_bw()
Keragaman antar Individu
gplots::plotmeans(PDN ~ `PROV`, main="Keragaman Antar Jumlah Pidana / 100rb Penduduk", crime, n.label = F, xlab = "a")
gplots::plotmeans(PDN_AREA ~ `PROV`, main="Keragaman Antar Jumlah Pidana / 100rb Penduduk / 100km2", crime, n.label = F, xlab = "a")
Pada plot keragaman antarwaktu untuk peubah jumlah tindak pidana per
100.000 penduduk per 100km persegi (PDN_AREA) mengalami penurunan jumlah
seiring waktunya. Kemudian, pada plot keragaman antarprovinsi untuk
peubah jumlah tindak pidana per 100.000 penduduk (PDN) terlihat cukup
beragam nilainya. Sementara untuk peubah PDN_AREA persegi terlihat
perbedaan nilai yang signifikan pada wilayah DKI Jakarta dan DI
Yogyakarta. Sebaran nilai pada plot kedua tersebut tidak terlihat
perbedaan yang jauh pada provinsi lainnya.
Correlation Plot
Sebelumnya disebutkan bahwa kami menggunakan pembobotan tambahan pada peubah respons (PDN). Alasan kami memilih pembobotan tambahan akan dijelaskan lebih lanjut pada plot korelasi di bawah.
corrplot(cor(crime[-c(1:2)]), method="color", type = "upper", tl.pos = 'tp')
corrplot(cor(crime[-c(1:2)]), method="number", diag=F, add=T, type = "lower",
tl.pos = 'n', cl.pos = 'n')
Dari plot korelasi di atas, dapat terlihat bahwa Jumlah Pidana yang
tidak terboboti luas wilayah memiliki korelasi pada peubah penjelas yang
relatif rendah (|r|<= 0.3). Sedangkan, pada Jumlah Pidana yang
terboboti oleh Luas Wilayah memiliki korelasi yang relatif lebih tinggi
jika dibanding Jumlah pidana yang tidak terboboti Luas Wilayah. Oleh
karena itu, akan lebih baik jika digunakan peubah respons yang terboboti
PDN_AREA jika dibandingkan PDN. Sehingga:
Time Series Plot
Perbandingan PDN dan PDN_AREA
crime %>%
mutate(TAHUN=as.numeric(as.character(TAHUN))) %>%
pivot_longer(cols = c(PDN, PDN_AREA), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = TAHUN, y = value, color = PROV)) +
geom_line() +
geom_point(size=2, shape=19, fill="red") +
facet_wrap(~variable, scales = "free_y") +
labs(title = "Jumlah Tindak Pidana / 100rb Penduduk vs\nJumlah Tindak Pidana / 100rb Penduduk / 100km2", x = "Year", y = "Value", color = "Province") +
theme(legend.position = "bottom", plot.title = element_text(face = "bold", hjust = 0.5))
Plot PDN
ggplot(data=crime, aes(x=TAHUN, y=PDN, group = PROV, colour = PROV))+ theme_bw()+
geom_line(size=1.2) +
geom_point(size=3, shape=19, fill="red") +
labs(colour="PROV", title = "Jumlah tindak pidana per 100.000 Penduduk", subtitle = "Tahun 2015-2020") +
theme(plot.title = element_text(face = "bold"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Plot PDN_AREA
ggplot(crime, aes(x=TAHUN, y=PDN_AREA, group = PROV, colour = PROV))+ theme_bw()+
geom_line(size=1.2) +
geom_point( size=3, shape=19, fill="red") +
labs(colour="PROV", title = "Jumlah Tindak Pidana / 100rb Penduduk / 100km2", subtitle = "Tahun 2015-2020") +
theme(plot.title = element_text(face = "bold"))
Plot PDN AREA tanpa DKI Jakarta dan DI Yogyakarta
crime %>%
filter(PROV!="DKI JAKARTA" & PROV!="D I YOGYAKARTA") %>%
ggplot(aes(x=TAHUN, y=PDN_AREA, group = PROV, colour = PROV))+ theme_bw()+
geom_line(size=1.2) +
geom_point( size=3, shape=19, fill="red") +
labs(colour="PROV", title = "Jumlah Tindak Pidana / 100rb Penduduk / 100km2", subtitle = "Tahun 2015-2020") +
theme(plot.title = element_text(face = "bold"))
Pada gambar di atas, terlihat bahwa perkembangan PDN_AREA setiap provinsi Indonesia tahun 2015-2020 cukup beragam. Ada provinsi tertentu yang meningkat, menurun, atau stagnan nilai peubah PDN_AREA provinsi Indonesia tahun 2015-2020 hanya terlihat pada provinsi DKI Jakarta dan DI Yogyakarta saja yang mana perkembanganya menurun. Untuk lebih jelasnya, dibuatlah gambar 7 yang mana gambar 7 merupakan plot yang sama dengan gambar 6. Hanya saja gambar 7 tidak memasukkan kategori provinsi DKI Jakarta dan DI Yogyakarta agar provinsi dengan PDN_AREA yang rendah dapat terlihat bagaimana perkembangannya. Terlihat bahwa nilai PDN_AREA masing-masing provinsi mayoritas mengalami penurunan pada rentang waktu 2015-2020.
Map Chart
Map di bawah ini dibuat pada Aplikasi Tableau Public, merupakan sebaran rata-rata peubah PDN dan PDN_AREA di Indonesia pada selama tahun 2015 hingga 2020.
Map chart untuk rataan PDN 2015-2020 untuk masing-masing provinsi menunjukkan bahwa di beberapa wilayah pelosok Sumatera, Sulawesi dan Papua masih cukup tinggi banyaknya tindak pidana yang terjadi. Kemudian, setelah peubah PDN diboboti dengan wilayah per 100km persegi, sebagian besar wilayah pelosok banyak tindak pidana terboboti wilayahnya cukup kecil (lebih kecil dari 1). Meskipun begitu, masih terlihat beberapa wilayah yang cukup tinggi seperti DKI Jakarta, DI Yogyakarta, Kepulauan Riau, Gorontalo, Bali, hingga Bengkulu.
Seleksi Peubah Multikolinearitas
# Menghapus peubah PDN karena sudah tidak akan digunakan
crime <- crime %>% select(-PDN)
check_model <- lm(PDN_AREA ~ KMIS + RLS + IPM + TPT + PDRB + AHH + PKD, crime)
data.frame(t(vif(check_model)))Dari hasil di atas, dapat terlihat bahwa peubah IPM memiliki VIF yang sangat tinggi, yaitu 38.78 yang mana lebih besar dari nol. Hal ini melanggar asumsi Multikolinearitas sehingga peubah ini akan dihilangkan.
crime <- crime %>% select(-IPM)
check_model <- lm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, crime)
data.frame(t(vif(check_model)))Setelah menghilangkan peubah IPM, didapat semua nilai VIF peubah menjadi << 10, artinya asumsi multikolinearitas terpenuhi.
Common-Effect Model
Common Effect Model (CEM) atau disebut juga pooled regression model merupakan salah satu model dalam regresi data panel yang mengkombinasikan data cross-section dan time series menjadi satu kesatuan (pooldata) yang kemudian diestimasi modelnya dengan teknik OLS. Akibatnya, perbedaan baik antar individu maupun antar waktu tidak dapat diamati. Dengan kata lain, dalam pendekatan ini dimensi individu maupun waktu tidak diperhatikan. Diasumsikan bahwa perilaku data antar Individu sama dalam berbagai kurun waktu. Persamaan pool model dapat dituliskan sebagai berikut:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+\varepsilon_{it} \\ \end{aligned} \]
dengan \(y_{ij}\) merupakan variabel dependen individu ke-i dan waktu ke-j. \(\alpha\) adalah intersep, \(x_{p,ij}\) adalah variabel independen ke-p individu ke-i dan waktu ke-j, \(\beta_k\) koefisien variabel bebas, dan \(\varepsilon_{ij}\) merupakan galat individu ke-i dan waktu ke-j.
Selain asumsi kesamaan pengaruh individu dan waktu, asumsi lain yang diperlukan pada model ini adalah asumsi normalitas, homogenitas, dan kebebasan galat.
Model CEM
Program R menyediakan package plm yang dapat memodelkan
data menjadi model CEM.
Dalam pemodelan ini, dapat dibuat hipotesis:
\[ \begin{aligned} H_0 &: \text{CEM adalah model yang tidak layak}\\ H_1 &: \text{CEM adalah model yang layak} \end{aligned} \]
library(plm)
cem <- plm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, data = crime, model = "pooling")
summary(cem)## Pooling Model
##
## Call:
## plm(formula = PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
## data = crime, model = "pooling")
##
## Balanced Panel: n = 34, T = 6, N = 204
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -5.72244 -1.14021 -0.18869 0.79534 15.12686
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -7.4014e+00 5.9327e+00 -1.2476 0.213673
## KMIS 2.1825e-01 3.9638e-02 5.5060 1.136e-07 ***
## RLS 6.9626e-01 2.5559e-01 2.7242 0.007026 **
## TPT -1.7381e-01 1.0775e-01 -1.6131 0.108330
## PDRB 2.6495e-05 5.0511e-06 5.2453 4.010e-07 ***
## AHH -1.7906e-01 8.3613e-02 -2.1415 0.033463 *
## PKD 1.1957e-03 1.2684e-04 9.4265 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2934.7
## Residual Sum of Squares: 1096.4
## R-Squared: 0.62639
## Adj. R-Squared: 0.61501
## F-statistic: 55.0488 on 6 and 197 DF, p-value: < 2.22e-16Pada dasarnya, pendugaan model yang dilakukan oleh metode CEM menghasilkan model yang sama persis dengan metode Regresi Linear Berganda dengan memandang setiap baris pada data adalah individu yang berbeda-beda. Model regresi yang dihasilkan tertulis di atas, dengan interpretasinya adalah:
Setiap kenaikan satu satuan
KMISmenyebabkanPDN_AREAmeningkat sebesar0.21825dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
RLSmenyebabkanPDN_AREAmeningkat sebesar0.6926dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
TPTmenyebabkanPDN_AREAmenurun sebesar0.17381dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PDRBmenyebabkanPDN_AREAmeningkat sebesar2.6495e-05dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
AHHmenyebabkanPDN_AREAmenurun sebesar0.17906dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PKDmenyebabkanPDN_AREAmeningkat sebesar1.1957e-03dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk peubah
KMIS,RLS,PDRB,AHH, danPKD< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.12645 atau 12.645% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
Meskipun sebagian besar parameter model signifikan, tetapi nilai R2 model kurang baik dan fakta bahwa model CEM mengabaikan pengaruh individu dan waktu tidak bisa diabaikan. Kedua faktor tersebut bisa berpengaruh terhadap keheterogenan sisaan pada model. Oleh karena itu, perlu dilakukan perbandingan model CEM dengan FEM dan REM untuk memastikan model terbaik.
Fixed-Effect Model
Fixed-Effect model merupakan salah satu jenis model dalam regresi data panel yang mempertimbangkan adanya pengaruh individu dan waktu. Terdapat asumsi bahwa intersep antar individu dan waktu berbeda tetapi koefisien regresi konstan untuk semua individu dan waktu. Model yang memiliki pengaruh dari salah satu dari individu atau waktu disebut model sisaan satu arah dan model yang memiliki pengaruh dari keduanya disebut model sisaan dua arah. Pendugaan pada model ini biasa dilakukan dengan model Within atau Least Square Dummy Variable (LSDV). Pada penelitian ini, hanya digunakan model within.
Persamaan model Fixed-Effect :
- Pengaruh Individu
\[ \begin{aligned} y_{it} = \alpha + f_i +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
- Pengaruh Waktu
\[ \begin{aligned} y_{it} = \alpha + \lambda_t +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
- Model Dua arah (pengaruh individu dan waktu)
\[ \begin{aligned} y_{it} = \alpha + f_i + \lambda_t +\beta_1x_{1,it} + \beta_2x_{2,it} + ... + \beta_kx_{k,it} + \varepsilon_{it} \end{aligned} \]
FEM Individual
Efek komponen sisaaan satu arah pada individu.
fem_ind <- plm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, index = c("PROV", "TAHUN"), model = "within", effect= "individual", data = crime)
summary(fem_ind)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
## data = crime, effect = "individual", model = "within", index = c("PROV",
## "TAHUN"))
##
## Balanced Panel: n = 34, T = 6, N = 204
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.453977 -0.176466 -0.013324 0.171397 2.908483
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## KMIS -2.0378e-02 1.0946e-01 -0.1862 0.852544
## RLS 1.5999e+00 5.5946e-01 2.8597 0.004792 **
## TPT -2.2009e-01 9.0523e-02 -2.4313 0.016122 *
## PDRB -7.7429e-05 8.5420e-06 -9.0645 3.662e-16 ***
## AHH -3.0887e-01 3.4891e-01 -0.8852 0.377324
## PKD -3.0543e-04 2.5732e-04 -1.1870 0.236953
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 106.82
## Residual Sum of Squares: 57.409
## R-Squared: 0.46255
## Adj. R-Squared: 0.33474
## F-statistic: 23.5243 on 6 and 164 DF, p-value: < 2.22e-16Interpretasi :
Setiap kenaikan satu satuan
KMISmenyebabkanPDN_AREAmenurun sebesar0.020378dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
RLSmenyebabkanPDN_AREAmeningkat sebesar1.5999dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
TPTmenyebabkanPDN_AREAmenurun sebesar0.22009dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PDRBmenyebabkanPDN_AREAmenurun sebesar7.7429e-05dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
AHHmenyebabkanPDN_AREAmenurun sebesar0.30887dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PKDmenyebabkanPDN_AREAmenurun sebesar3.0543e-04dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk peubah
RLS,TPT, danPDRB< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.46225 atau 46.225% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
summary(fixef(fem_ind, effect="individual"))## Estimate Std. Error t-value Pr(>|t|)
## ACEH 14.233 21.226 0.6705 0.50345
## BALI 18.482 21.927 0.8429 0.40050
## BANTEN 17.737 21.130 0.8394 0.40245
## BENGKULU 15.299 21.187 0.7221 0.47126
## D I YOGYAKARTA 22.604 23.001 0.9828 0.32717
## DKI JAKARTA 51.759 21.403 2.4183 0.01669 *
## GORONTALO 17.789 21.180 0.8399 0.40218
## JAMBI 17.485 21.737 0.8044 0.42234
## JAWA BARAT 17.756 22.361 0.7941 0.42831
## JAWA TENGAH 18.469 23.360 0.7906 0.43029
## JAWA TIMUR 18.801 22.226 0.8459 0.39886
## KALIMANTAN BARAT 17.009 21.897 0.7768 0.43841
## KALIMANTAN SELATAN 16.292 20.854 0.7812 0.43579
## KALIMANTAN TENGAH 16.312 21.232 0.7683 0.44341
## KALIMANTAN TIMUR 25.722 22.294 1.1538 0.25028
## KALIMANTAN UTARA 21.671 22.018 0.9842 0.32646
## KEP. BANGKA BELITUNG 18.657 21.683 0.8604 0.39082
## KEPULAUAN RIAU 22.799 20.714 1.1006 0.27267
## LAMPUNG 16.563 21.829 0.7588 0.44908
## MALUKU 11.845 19.709 0.6010 0.54866
## MALUKU UTARA 12.996 20.406 0.6369 0.52510
## NUSA TENGGARA BARAT 16.044 20.716 0.7745 0.43976
## NUSA TENGGARA TIMUR 13.928 20.946 0.6649 0.50702
## PAPUA 17.731 21.124 0.8393 0.40250
## PAPUA BARAT 19.225 20.657 0.9307 0.35338
## RIAU 21.192 21.615 0.9804 0.32831
## SULAWESI BARAT 14.795 19.981 0.7405 0.46007
## SULAWESI SELATAN 17.720 21.620 0.8196 0.41363
## SULAWESI TENGAH 15.746 20.782 0.7577 0.44971
## SULAWESI TENGGARA 15.486 21.680 0.7143 0.47608
## SULAWESI UTARA 18.161 21.502 0.8446 0.39957
## SUMATERA BARAT 15.664 20.873 0.7504 0.45406
## SUMATERA SELATAN 16.988 21.504 0.7900 0.43065
## SUMATERA UTARA 15.076 20.601 0.7318 0.46532
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nilai diatas adalah nilai pengaruh konstan dari masing-masing individu yang pada model dapat dituliskan sebagai \(f_i\).
FEM TIME
Efek komponen sisaaan satu arah pada waktu.
fem_time <- plm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, index = c("PROV", "TAHUN"), model = "within", effect= "time", data = crime)
summary(fem_time)## Oneway (time) effect Within Model
##
## Call:
## plm(formula = PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
## data = crime, effect = "time", model = "within", index = c("PROV",
## "TAHUN"))
##
## Balanced Panel: n = 34, T = 6, N = 204
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -5.0372730 -0.9663536 0.0034536 0.7862123 13.1063295
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## KMIS 2.2939e-01 3.7023e-02 6.1959 3.458e-09 ***
## RLS 1.0310e+00 2.4524e-01 4.2043 4.013e-05 ***
## TPT -3.2468e-01 1.0437e-01 -3.1109 0.002149 **
## PDRB 2.6992e-05 4.7139e-06 5.7261 3.910e-08 ***
## AHH -1.8149e-01 7.8086e-02 -2.3242 0.021161 *
## PKD 1.2167e-03 1.1833e-04 10.2818 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2925.2
## Residual Sum of Squares: 928.18
## R-Squared: 0.68269
## Adj. R-Squared: 0.66451
## F-statistic: 68.8482 on 6 and 192 DF, p-value: < 2.22e-16Interpretasi:
Setiap kenaikan satu satuan
KMISmenyebabkanPDN_AREAmeningkat sebesar0.22939dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
RLSmenyebabkanPDN_AREAmeningkat sebesar1.031dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
TPTmenyebabkanPDN_AREAmenurun sebesar0.32648dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PDRBmenyebabkanPDN_AREAmeningkat sebesar2.6992e-05dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
AHHmenyebabkanPDN_AREAmenurun sebesar0.18149dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PKDmenyebabkanPDN_AREAmeningkat sebesar1.2167e-03dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk seluruh koefisien peubah < 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.68269 atau 68.269% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
summary(fixef(fem_time, effect="time"))## Estimate Std. Error t-value Pr(>|t|)
## 2015 -8.1078 5.5630 -1.4575 0.14662
## 2016 -8.7157 5.5607 -1.5674 0.11867
## 2017 -9.5108 5.5578 -1.7112 0.08865 .
## 2018 -10.1735 5.5643 -1.8284 0.06905 .
## 2019 -10.7981 5.5820 -1.9344 0.05453 .
## 2020 -10.5003 5.6058 -1.8731 0.06257 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nilai diatas adalah nilai pengaruh konstan dari masing-masing waktu yang pada model dapat dituliskan sebagai \(\lambda_t\).
FEM TWO WAYS
Efek komponen sisaaan satu arah pada individu.
fem_twoways <- plm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, index = c("PROV", "TAHUN"), model = "within", effect= "twoways", data = crime)
summary(fem_twoways)## Twoways effects Within Model
##
## Call:
## plm(formula = PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
## data = crime, effect = "twoways", model = "within", index = c("PROV",
## "TAHUN"))
##
## Balanced Panel: n = 34, T = 6, N = 204
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.34461974 -0.21115029 0.00011108 0.19612635 2.64700693
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## KMIS -2.3596e-03 1.1019e-01 -0.0214 0.98294
## RLS 1.8633e+00 8.5443e-01 2.1807 0.03067 *
## TPT -1.2729e-01 1.0033e-01 -1.2688 0.20638
## PDRB -8.0578e-05 8.8120e-06 -9.1441 2.736e-16 ***
## AHH -5.3199e-01 4.6077e-01 -1.1546 0.25000
## PKD -9.2309e-04 4.2894e-04 -2.1520 0.03291 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 97.319
## Residual Sum of Squares: 54.356
## R-Squared: 0.44147
## Adj. R-Squared: 0.28691
## F-statistic: 20.9458 on 6 and 159 DF, p-value: < 2.22e-16Interpretasi :
Setiap kenaikan satu satuan
KMISmenyebabkanPDN_AREAmenurun sebesar2.3596e-03dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
RLSmenyebabkanPDN_AREAmeningkat sebesar1.8633dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
TPTmenyebabkanPDN_AREAmenurun sebesar0.12729dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PDRBmenyebabkanPDN_AREAmenurun sebesar8.0578e-05dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
AHHmenyebabkanPDN_AREAmenurun sebesar0.53199dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PKDmenyebabkanPDN_AREAmenurun sebesar9.2309e-04dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk peubah
RLS,PDRB, danPKD< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.44147 atau 44.147% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
data.frame(summary(fixef(fem_twoways, effect="twoways")))Nilai diatas adalah nilai pengaruh konstan dari gabungan masing-masing individu dan waktu yang pada model dapat dituliskan sebagai \(f_i+\lambda_t\).
Perbandingan Kebaikan FEM Individu, Waktu, dan Twoways
Hipotesis pada uji ini yaitu :
- Untuk menguji pengaruh individu
\[ \begin{aligned} H_{0} &: \text{Tidak ada pengaruh spesifik individu} \\ H_{1} &: \text{Ada pengaruh spesifik individu} \end{aligned} \]
- Untuk menguji pengaruh waktu
\[ \begin{aligned} H_{0} &: \text{Tidak ada pengaruh spesifik waktu} \\ H_{1} &: \text{Ada pengaruh spesifik waktu} \end{aligned} \]
Statsitik Uji yang digunakan :
- Pengaruh individu
\[ \begin{aligned} F_{Hitung} = \frac{JKG_{MG} - JKG_{MPTI}}{JKG_{MPTI}} . \frac{NT-N-K}{N-1} \end{aligned} \]
- Pengaruh Waktu
\[ \begin{aligned} F_{Hitung} = \frac{JKG_{MG} - JKG_{MPTW}}{JKG_{MPTW}} . \frac{N-T-K}{T-1} \end{aligned} \]
Uji Signifikasi Pengaruh Individu / Waktu / Two ways
Uji pengaruh individu
plmtest(fem_twoways, type = "bp", effect = "individual")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 237.67, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsUji pengaruh waktu
plmtest(fem_twoways, type = "bp", effect = "time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 40.608, df = 1, p-value = 1.861e-10
## alternative hypothesis: significant effectsUji pengaruh twoways
plmtest(fem_twoways, type = "bp", effect = "twoways")##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 278.28, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectsDari ketiga uji diatas didapat nilai p-value yang lebih kecil dari 0.05 yang berarti keputusannya adalah Tak Tolak \(H_0\) atau terdapat pengaruh yang signifikan baik dari individu, waktu, atau keduanya secara bersamaan. Sehingga, kita perlu melihat metrik lain dalam menentukan metode terpilih.
Nilai Kebaikan Model
# Sum Squared Error
dsse <- data.frame(Individu=sum(fem_ind$residuals^2),Time=sum(fem_time$residuals^2),Twoways=sum(fem_twoways$residuals^2))
# AIC
lsdv_ind <- lm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD + PROV, crime)
lsdv_time <- lm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD + TAHUN, crime)
lsdv_twoways <- lm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD + PROV + TAHUN, crime)
daic <- data.frame(Individu=AIC(lsdv_ind),Time=AIC(lsdv_time),Twoways=AIC(lsdv_twoways))
# MAPE
mape <- function(actual, forecast) {
mean(abs((actual - forecast) / actual)) * 100
}
dmape <- data.frame(Individu=mape(crime$PDN_AREA,predict(fem_ind)),
Time=mape(crime$PDN_AREA,predict(fem_time)),
Twoways=mape(crime$PDN_AREA,predict(fem_twoways)))
# BIC
dbic <- data.frame(Individu=BIC(lsdv_ind),Time=BIC(lsdv_time),Twoways=BIC(lsdv_twoways))
# R-squared
ind <- summary(fem_ind)
time <- summary(fem_time)
tways <- summary(fem_twoways)
drsq <- data.frame(Individu=ind$r.squared,Time=time$r.squared,Twoways=tways$r.squared)
# Perbandingan
compare <- t(rbind(dsse,daic,dbic,dmape,drsq))
colnames(compare) <- c("SSE","AIC","BIC", "MAPE","R-Squared","Adj R-Squared")
compare## SSE AIC BIC MAPE R-Squared Adj R-Squared
## Individu 57.40917 402.2721 538.3150 618.7911 0.4625523 0.3347446
## Time 928.18393 914.0093 957.1449 873.7714 0.6826914 0.6645123
## Twoways 54.35551 401.1218 553.7554 619.0807 0.4414685 0.2869063Dari hasil di atas, dapat terlihat bahwa mesipun model pengaruh waktu merupakan model yang paling signifikan dengan nilai R-Squared yang besar jika dibandingkan dengan yang lain, tetapi error yang didapat oleh model waktu sangatlah tinggi. Oleh karena itu, akan dipilih model Twoways yang memiliki nilai error relatif lebih rendah daripada model Individu.
FEM vs CEM
Untuk membandingkan kebaikan model FEM atau CEM, uji Chow akan digunakan. Uji Chow adalah uji formal untuk menguji ada atau tidaknya pengaruh spesifik individu dan pengaruh spesifik waktu pada data panel.
\[ \begin{aligned} H_0 &: \text{Model Common Effect} \\ H_1 &: \text{Model Fixed Effect} \end{aligned} \]
\(H_0\) akan ditolak apabila
p-value lebih kecil dari alpha yaitu 0.05
pooltest(cem, fem_twoways)##
## F statistic
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## F = 80.216, df1 = 38, df2 = 159, p-value < 2.2e-16
## alternative hypothesis: unstabilityNilai p-value sebesar \(2.2e-16\) berarti keputusan Tolak
H0 atau model fixed effect lebih layak digunakan pada data ini.
Random-Effect Model
Model pengaruh acak mengasumsikan tidak ada korelasi antara pengaruh spesifik individu dan pengaruh spesifik waktu dengan peubah bebas. Asumsi ini membuat komponen sisaan dari pengaruh spesifik individu dan pengaruh spesifik waktu dimasukan kedalam sisaan. Model sisaan satu arah yaitu:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ f_{i} +\varepsilon_{it}\ \text{(pengaruh Individu)} \\ \end{aligned} \]
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ \lambda_{t} +\varepsilon_{it} \\ \end{aligned} \]
Model sisaan dua arah yaitu:
\[ \begin{aligned} y_{it} &= \alpha + \beta_{1}x_{1,it}+ \beta_{2}x_{2,it}+...+\beta_{k}x_{k,it}+ f_{i}+ \lambda_{t} +\varepsilon_{it} \\ \end{aligned} \] Dengan \(\alpha\) merupakan intersep model dan \(f_i\) dan \(\lambda_t\) merupakan pengaruh acak masing-masing individu dan waktu. Pendugaan parameter model pengaruh acak menggunakan Generalized Least Square (Baltagi 2008).
Generalized Least Square
Analisis data panel tidak sepenuhnya mengikuti prosedur regresi biasa yang tidak dapat menjelaskan keragaman data. Perlunya mempertimbangkan keragaman data tersebut untuk dapat mengetahui dalam penggalian lebih lanjut terhadap informasi yang dikandung data. Salah satu cara untuk melakukan hal tersebut yaitu dengan metode Generalized Least Square (GLS). GLS dipakai sebagai estimator analisis data panel dengan persamaan sistem karena informasi heterogenitas antara individu dan waktu digunakan sebagai informasi yang penting untuk menghasilkan parameter \(\hat{\beta}\) Persamaan dalam pendugaan parameter GLS adalah sebagai berikut.
\[ \begin{aligned} \hat {\beta} &= (X'(\sum)^{-1} X )^{-1} X'(\sum) ^{-1} y \\ \end{aligned} \]
Metode GLS memasukkan struktur matriks ragam peragam residual pada parameter \(\hat{\beta}\) (Ekananda 2016).
Pendugaan Parameter
Jika pada model FEM mengasumsikan bahwa nilai dari pengaruh spesifik individu bersifat fixed pada setiap unit, maka pada model random effect diasumsikan bahwa pengaruh tersebut bersifat acak. Dalam hal ini misalnya model one way individu pada data Grundfeld, maka ke-10 firm memiliki nilai rataan umum intersep yang sama, sedangkan perbedaan intersep antar masing-masing unit firm direfleksikan pada error termnya.
rem_gls <- plm(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, data = crime,
index = c("PROV", "TAHUN"),
effect = "twoways", model = "random", random.method = "nerlove")
summary(rem_gls)## Twoways effects Random Effect Model
## (Nerlove's transformation)
##
## Call:
## plm(formula = PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
## data = crime, effect = "twoways", model = "random", random.method = "nerlove",
## index = c("PROV", "TAHUN"))
##
## Balanced Panel: n = 34, T = 6, N = 204
##
## Effects:
## var std.dev share
## idiosyncratic 0.2664 0.5162 0.004
## individual 61.9916 7.8735 0.994
## time 0.1274 0.3569 0.002
## theta: 0.9732 (id) 0.7592 (time) 0.7591 (total)
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.667620 -0.230679 -0.031907 0.137530 3.787819
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 3.0779e+00 2.2719e+01 0.1355 0.8922353
## KMIS -1.7684e-02 9.9666e-02 -0.1774 0.8591664
## RLS 2.0957e+00 6.0135e-01 3.4849 0.0004922 ***
## TPT -1.2347e-01 9.3177e-02 -1.3251 0.1851276
## PDRB -7.3490e-05 7.8922e-06 -9.3118 < 2.2e-16 ***
## AHH -1.4412e-01 3.5020e-01 -0.4115 0.6806785
## PKD -3.9659e-04 2.9043e-04 -1.3655 0.1720944
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 99.894
## Residual Sum of Squares: 60.882
## R-Squared: 0.39053
## Adj. R-Squared: 0.37196
## Chisq: 126.23 on 6 DF, p-value: < 2.22e-16Interpretasi :
Setiap kenaikan satu satuan
KMISmenyebabkanPDN_AREAmenurun sebesar0.017684dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
RLSmenyebabkanPDN_AREAmeningkat sebesar2.0957dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
TPTmenyebabkanPDN_AREAmenurun sebesar0.12347dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PDRBmenyebabkanPDN_AREAmenurun sebesar7.3490e-05dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
AHHmenyebabkanPDN_AREAmenurun sebesar0.14412dengan menganggap peubah lain konstan.Setiap kenaikan satu satuan
PKDmenyebabkanPDN_AREAmenurun sebesar3.9659e-04dengan menganggap peubah lain konstan.
Untuk semua p-value di bawah 5%, maka terdapat cukup bukti untuk menolak \(H_0\) (kurang dari 5% kemungkinan bahwa \(H_0\) benar).
p-value uji-F < 2.22e-16 << 5%. Artinya, dapat dikatakan bahwa model layak untuk digunakan.
p-value uji-T untuk peubah
RLSdanPDRB< 5%. Artinya, dapat dikatakan bahwa peubah tersebut berpengaruh signifikan terhadap model.
\(R^2\) merupakan nilai yang menunjukan kekuatan hubungan antara peubah bebas dan peubah respons. Semakin tinggi nilai \(R^2\), maka semakin baik model tersebut.
- \(R^2\) pada model sebesar 0.39053 atau 39.053% yang mengindikasikan bahwa peubah yang digunakan kurang baik dalam menggambarkan model.
Perbandingan Kebaikan REM
Uji Pengaruh Individu
#efek individu
plmtest(rem_gls,type = "bp", effect="individu")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 237.67, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsKarena nilai p-value < 0.05, maka tolak \(H_0\) atau pengaruh individu signifikan terhadap model.
Uji pengaruh waktu
#efek waktu
plmtest(rem_gls,type = "bp", effect="time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 40.608, df = 1, p-value = 1.861e-10
## alternative hypothesis: significant effectsKarena nilai p-value < 0.05, maka tolak \(H_0\) atau pengaruh waktu signifikan terhadap model. Artinya, model twoways tepat digunakan untuk metode REM.
Uji pengaruh twoways
#efek waktu
plmtest(rem_gls,type = "bp", effect="twoways")##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 278.28, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectstidy_ranef <- tidy(ranef(rem_gls, effect="individual"))
colnames(tidy_ranef) <- c("Provinsi", "Pengaruh Acak Individu")
tidy_ranefOutput diatas merupakan pengaruh acak dari setiap unit individu, nilai tersebut tersebut menunjukkan seberapa besar perbedaan nilai komponen error acak masing-masing unit indvidu terhadap nilai intersep umum.
tidy_ranef <- tidy(ranef(rem_gls, effect="time"))
colnames(tidy_ranef) <- c("Tahun", "Pengaruh Acak Waktu")
tidy_ranefOutput diatas merupakan pengaruh acak dari setiap unit waktu, nilai tersebut tersebut menunjukkan seberapa besar perbedaan nilai komponen error acak masing-masing unit waktu terhadap nilai intersep umum.
FEM vs REM
Uji Haussman
Pada bagian sebelumnya, cukup bukti yang menyatakan bahwa model FEM
lebih baik daripada model CEM. Oleh karena itu, tahap selanjutnya akan
dibandingkan kebaikan model FEM jika dibanding model REM. Uji yang
digunakan adalah uji Hausman dengan menggunakan fungsi
phtest(FEM,REM).
Pengujian ini dilakukan untuk memilih model antara model pengaruh tetap atau model pengaruh acak yang sesuai untuk menggambarkan suatu data panel. Uji Hausman didasarkan pada perbedaan penduga model pengaruh tetap \(\hat{\beta}_{MPT}\) dengan penduga model pengaruh acak \(\hat{\beta}_{MPA}\). Kedua penduga konsisten dalam kondisi \(H_0\), tetapi \(\hat{\beta}_{MPA}\) akan bias dan tidak konsisten pada \(H_1\). Pengujian hipotesisnya adalah:
\(H_0\): Model pengaruh acak (REM) adalah model yang tepat
\(H_1\): Model pengaruh tetap (FEM) adalah model yang tepat
#fem dengan rem gls
phtest(fem_twoways, rem_gls)##
## Hausman Test
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 6.0755, df = 6, p-value = 0.4148
## alternative hypothesis: one model is inconsistentKarena nilai p-value yang diperoleh lebih dari dari alpha 5%, maka Terima H0 yang mengindikasikan bahwa model yang sesuai adalah model REM.
UJi Diagnostik Sisaan
Normalitas
ks.test(rem_gls$residuals, "pnorm",
mean=mean(rem_gls$residuals),
sd=sd(rem_gls$residuals))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: rem_gls$residuals
## D = 0.18177, p-value = 2.797e-06
## alternative hypothesis: two-sidedHistogram
# Histogram
ggplot(as.data.frame(rem_gls$residuals), aes(x = rem_gls$residuals)) +
geom_histogram(aes(y = after_stat(density)), color = "white", fill = "steelblue") +
geom_density(color = "red", linewidth = 1) +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Uji Autokorelasi
pbgtest(rem_gls)##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
##
## data: PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD
## chisq = 62.514, df = 6, p-value = 1.387e-11
## alternative hypothesis: serial correlation in idiosyncratic errorsUji Kehomogenan
bptest(rem_gls)##
## studentized Breusch-Pagan test
##
## data: rem_gls
## BP = 55.643, df = 6, p-value = 3.437e-10Kesimpulan uji diagnostik :
- Sisaan tidak menyebar normal (Asumsi kenormalan tidak terpenuhi)
- Sisaan memiliki autokorelasi (Asumsi kebebasan belum terpenuhi)
- Sisaan tidak homogen (Asumsi kehomogenan belum terpenuhi)
Terlihat bahwa asumsi kehomogenan tidak terpenuhi. Oleh karena itu, diperlukan penanganan yaitu dengan menggunakan Analisis Geographically Weighted Panel Regression.
Geographically Weighted Panel Regression
Ide utama analisis Geographically Weighted Panel Regression (GWPR) sama halnya dengan analisis GWR cross-sectional. Analisis GWPR mengasumsikan runtutan waktu dari observasi pada sebuah lokasi geografis merupakan realisasi dari sebuah proses smooth spatiotemporal. Proses ini mengikuti sebuah distribusi observasi terdekat lebih berhubungan dari pada observasi yang jauh. Analisis GWPR bertujuan menggabungkan secara keseluruhan lokasi dan observasi (Yu 2010). Pendugaan parameter model GWPR menggunakan pendekatan Weighted Least Square (WLS) seperti pada model GWR.
Input Latitude Longitude
library(jsonlite)
wilayah <- fromJSON("https://raw.githubusercontent.com/yusufsyaifudin/wilayah-indonesia/master/data/list_of_area/provinces.json")
wilayahdata.sp.GWPR <- crime %>%
mutate(LAT=rep(wilayah$latitude,6), LONG=rep(wilayah$longitude,6)) %>%
select(PROV, TAHUN, LAT, LONG, everything())
data.sp.GWPRModeling GWPR
library(GWmodel)
coordinates(data.sp.GWPR)=3:4
class(data.sp.GWPR)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"head(data.sp.GWPR,34)## coordinates PROV TAHUN PDN_AREA KMIS RLS TPT
## 1 (4.36855, 97.0253) ACEH 2020 0.25709159 15.21 9.33 6.00
## 2 (2.19235, 99.38122) SUMATERA UTARA 2020 0.31651974 8.95 9.54 5.81
## 3 (-1.34225, 100.0761) SUMATERA BARAT 2020 0.35703328 6.42 8.99 6.07
## 4 (0.50041, 101.5476) RIAU 2020 0.14134087 6.93 9.14 5.62
## 5 (-1.61157, 102.7797) JAMBI 2020 0.26768862 7.78 8.55 4.70
## 6 (-3.12668, 104.0931) SUMATERA SELATAN 2020 0.16049361 12.82 8.24 4.71
## 7 (-3.51868, 102.536) BENGKULU 2020 0.86348286 15.17 8.84 3.58
## 8 (-4.8555, 105.0273) LAMPUNG 2020 0.26571318 12.55 8.05 4.47
## 9 (-2.75775, 107.5839) KEP. BANGKA BELITUNG 2020 0.80978759 4.71 8.06 4.30
## 10 (-0.15478, 104.5804) KEPULAUAN RIAU 2020 1.54845569 6.03 10.12 8.16
## 11 (6.1745, 106.8227) DKI JAKARTA 2020 15.81301486 4.61 11.13 8.05
## 12 (-6.88917, 107.6405) JAWA BARAT 2020 0.08197240 8.16 8.55 9.09
## 13 (-7.30324, 110.0044) JAWA TENGAH 2020 0.09451021 11.63 7.69 5.34
## 14 (7.7956, 110.3695) D I YOGYAKARTA 2020 6.38335222 12.54 9.55 3.98
## 15 (-6.96851, 113.98) JAWA TIMUR 2020 0.09204349 11.28 7.78 4.72
## 16 (-6.44538, 106.1376) BANTEN 2020 0.65197684 6.28 8.89 9.32
## 17 (-8.23566, 115.1224) BALI 2020 1.03805151 4.12 8.95 3.44
## 18 (-8.12179, 117.637) NUSA TENGGARA BARAT 2020 0.89918761 14.10 7.31 3.63
## 19 (-8.56568, 120.6979) NUSA TENGGARA TIMUR 2020 0.18063102 21.06 7.63 3.46
## 20 (-0.13224, 111.0969) KALIMANTAN BARAT 2020 0.05159293 7.21 7.37 5.14
## 21 (-1.49958, 113.2903) KALIMANTAN TENGAH 2020 0.04037391 5.04 8.59 3.96
## 22 (-2.94348, 115.3756) KALIMANTAN SELATAN 2020 0.50588178 4.61 8.29 4.21
## 23 (0.78844, 116.242) KALIMANTAN TIMUR 2020 0.07747935 6.37 9.77 6.80
## 24 (2.72594, 116.911) KALIMANTAN UTARA 2020 0.19346025 7.11 9.00 5.34
## 25 (0.65557, 124.0902) SULAWESI UTARA 2020 1.81393230 7.70 9.49 6.36
## 26 (-1.69378, 120.8089) SULAWESI TENGAH 2020 0.28945062 12.99 8.83 3.35
## 27 (-3.64467, 119.9472) SULAWESI SELATAN 2020 0.31037633 8.86 8.38 6.01
## 28 (-3.54912, 121.728) SULAWESI TENGGARA 2020 0.21277881 11.35 9.04 3.84
## 29 (0.71862, 122.4556) GORONTALO 2020 1.90102753 15.41 7.82 3.79
## 30 (-2.49745, 119.3919) SULAWESI BARAT 2020 0.74461583 11.19 7.89 2.86
## 31 (-3.11884, 129.4208) MALUKU 2020 0.64586223 17.72 9.93 7.14
## 32 (0.63012, 127.972) MALUKU UTARA 2020 0.21574298 6.88 9.04 4.62
## 33 (-1.38424, 132.9025) PAPUA BARAT 2020 0.31858533 21.54 7.60 6.79
## 34 (-3.98857, 138.3485) PAPUA 2020 0.06519639 26.72 6.69 3.85
## PDRB AHH PKD
## 1 31633.38 69.975 9492
## 2 54979.04 69.150 10420
## 3 43825.66 69.520 10733
## 4 114166.90 71.650 10675
## 5 57957.73 71.170 10392
## 6 53842.74 69.930 10652
## 7 36552.50 69.370 10380
## 8 39290.33 70.695 9982
## 9 52023.40 70.680 12794
## 10 123464.79 69.990 14209
## 11 262615.17 72.950 18227
## 12 43236.51 73.150 10845
## 13 36964.78 74.405 10930
## 14 37693.64 75.025 14015
## 15 56640.82 71.345 11601
## 16 52729.40 70.005 11964
## 17 52015.45 72.155 13929
## 18 25183.56 66.510 10351
## 19 20056.99 67.055 7598
## 20 39622.24 70.735 8930
## 21 57145.08 69.750 11154
## 22 44100.79 68.690 12032
## 23 161798.85 74.375 11728
## 24 143533.29 72.535 8756
## 25 50521.13 71.745 10791
## 26 66306.27 68.745 9335
## 27 55675.03 70.620 11079
## 28 49718.15 71.340 9331
## 29 35693.26 68.115 10020
## 30 32836.75 65.110 9168
## 31 25094.36 66.025 8732
## 32 33069.32 68.375 8032
## 33 73932.60 66.050 8086
## 34 46416.36 65.835 69543.4.1. Menentukan Fungsi Pembobot Spasial Terbaik
Untuk melakukan pemodelan GWPR, tahap awal yang harus dilakukan adalah mengukur jarak euclidean antar provinsi. Data koordinat geografis, yaitu latitude dan longitude, digunakan untuk mengukur jarak euclidean. Data koordinat geografis yang digunakan dalam penelitian ini adalah centroid dari masing-masing provinsi dengan menggunakan software R. Setelah mendapatkan jarak euclidean, tahap selanjutnya adalah memilih fungsi pembobot spasial yang paling sesuai untuk mendapatkan matriks pembobot spasial yang akan digunakan dalam pemodelan GWPR.
# menggunakan 3 fungsi kernel: gaussian, eksponensial, bisquare )adaptive dan fixed)
#adaptive bisquare
bwd.GWPR.bisquare.ad <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, data = data.sp.GWPR, approach = "CV", kernel = "bisquare", adaptive=T)## Adaptive bandwidth: 133 CV score: 693.3964
## Adaptive bandwidth: 90 CV score: 556.5396
## Adaptive bandwidth: 62 CV score: 430.2603
## Adaptive bandwidth: 46 CV score: 369.3776
## Adaptive bandwidth: 35 CV score: 324.3215
## Adaptive bandwidth: 29 CV score: 187.1796
## Adaptive bandwidth: 24 CV score: 88.422
## Adaptive bandwidth: 22 CV score: 88.422hasil.GWPR.bisquare.ad <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, data = data.sp.GWPR, bw = bwd.GWPR.bisquare.ad, kernel = "bisquare", adaptive=T)
#adaptive gaussian
bwd.GWPR.gaussian.ad <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD, data = data.sp.GWPR, approach = "CV", kernel = "gaussian", adaptive=T)## Adaptive bandwidth: 133 CV score: 1073.73
## Adaptive bandwidth: 90 CV score: 932.2855
## Adaptive bandwidth: 62 CV score: 861.8887
## Adaptive bandwidth: 46 CV score: 740.2562
## Adaptive bandwidth: 35 CV score: 641.6354
## Adaptive bandwidth: 29 CV score: 572.2988
## Adaptive bandwidth: 24 CV score: 514.914
## Adaptive bandwidth: 22 CV score: 514.914hasil.GWPR.gaussian.ad <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, bw = bwd.GWPR.gaussian.ad, kernel = "gaussian", adaptive=T)
#adaptive exponential
bwd.GWPR.exponential.ad <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, approach = "CV", kernel = "exponential", adaptive=T)## Adaptive bandwidth: 133 CV score: 895.3918
## Adaptive bandwidth: 90 CV score: 783.4486
## Adaptive bandwidth: 62 CV score: 730.2009
## Adaptive bandwidth: 46 CV score: 655.3079
## Adaptive bandwidth: 35 CV score: 600.9486
## Adaptive bandwidth: 29 CV score: 558.1274
## Adaptive bandwidth: 24 CV score: 518.1572
## Adaptive bandwidth: 22 CV score: 518.1572hasil.GWPR.exponential.ad <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, bw = bwd.GWPR.exponential.ad, kernel = "exponential", adaptive=T)
#fixed bisquare
bwd.GWPR.bisquare.fix <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, approach = "CV", kernel = "bisquare", adaptive=F)## Fixed bandwidth: 26.05943 CV score: 973.6662
## Fixed bandwidth: 16.10883 CV score: 809.8809
## Fixed bandwidth: 9.959026 CV score: 521.4118
## Fixed bandwidth: 6.158237 CV score: 119.7549
## Fixed bandwidth: 3.809221 CV score: 169846.3
## Fixed bandwidth: 7.610009 CV score: 117.1879
## Fixed bandwidth: 8.507254 CV score: 353.8743
## Fixed bandwidth: 7.055482 CV score: 117.4174
## Fixed bandwidth: 7.952726 CV score: 308.7523
## Fixed bandwidth: 7.398199 CV score: 113.4041
## Fixed bandwidth: 7.267292 CV score: 113.4004
## Fixed bandwidth: 7.186388 CV score: 114.3983
## Fixed bandwidth: 7.317294 CV score: 113.1673
## Fixed bandwidth: 7.348197 CV score: 113.1696
## Fixed bandwidth: 7.298195 CV score: 113.2215
## Fixed bandwidth: 7.329098 CV score: 113.1551
## Fixed bandwidth: 7.336393 CV score: 113.1556
## Fixed bandwidth: 7.324589 CV score: 113.1578
## Fixed bandwidth: 7.331884 CV score: 113.1545
## Fixed bandwidth: 7.333607 CV score: 113.1547
## Fixed bandwidth: 7.33082 CV score: 113.1546hasil.GWPR.bisquare.fix <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, bw = bwd.GWPR.bisquare.fix, kernel = "bisquare", adaptive=F)
#fixed gaussian
bwd.GWPR.gaussian.fix <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, approach = "CV", kernel = "gaussian", adaptive=F)## Fixed bandwidth: 26.05943 CV score: 1170.719
## Fixed bandwidth: 16.10883 CV score: 1079.048
## Fixed bandwidth: 9.959026 CV score: 917.7705
## Fixed bandwidth: 6.158237 CV score: 734.8031
## Fixed bandwidth: 3.809221 CV score: 569.5653
## Fixed bandwidth: 2.357449 CV score: 203.5027
## Fixed bandwidth: 1.460204 CV score: 460.5427
## Fixed bandwidth: 2.911976 CV score: 391.4167
## Fixed bandwidth: 2.014732 CV score: 102.411
## Fixed bandwidth: 1.802921 CV score: 70.63039
## Fixed bandwidth: 1.672015 CV score: 186.1903
## Fixed bandwidth: 1.883825 CV score: 63.42011
## Fixed bandwidth: 1.933827 CV score: 74.38713
## Fixed bandwidth: 1.852923 CV score: 61.60032
## Fixed bandwidth: 1.833824 CV score: 63.08311
## Fixed bandwidth: 1.864726 CV score: 61.73331
## Fixed bandwidth: 1.845628 CV score: 61.90642
## Fixed bandwidth: 1.857431 CV score: 61.56244
## Fixed bandwidth: 1.860218 CV score: 61.59454
## Fixed bandwidth: 1.855709 CV score: 61.56364
## Fixed bandwidth: 1.858496 CV score: 61.56978
## Fixed bandwidth: 1.856774 CV score: 61.56099
## Fixed bandwidth: 1.856367 CV score: 61.56127
## Fixed bandwidth: 1.857025 CV score: 61.56126
## Fixed bandwidth: 1.856618 CV score: 61.56099hasil.GWPR.gaussian.fix <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, bw = bwd.GWPR.gaussian.fix, kernel = "gaussian", adaptive=F)
#fixed exponential
bwd.GWPR.exponential.fix <- bw.gwr(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, approach = "CV", kernel = "exponential", adaptive=F)## Fixed bandwidth: 26.05943 CV score: 1013.952
## Fixed bandwidth: 16.10883 CV score: 909.6848
## Fixed bandwidth: 9.959026 CV score: 783.5012
## Fixed bandwidth: 6.158237 CV score: 649.328
## Fixed bandwidth: 3.809221 CV score: 523.2297
## Fixed bandwidth: 2.357449 CV score: 404.5852
## Fixed bandwidth: 1.460204 CV score: 165.9064
## Fixed bandwidth: 0.9056765 CV score: 84.65157
## Fixed bandwidth: 0.5629596 CV score: 529.4583
## Fixed bandwidth: 1.117487 CV score: 74.65
## Fixed bandwidth: 1.248393 CV score: 103.3616
## Fixed bandwidth: 1.036583 CV score: 64.81533
## Fixed bandwidth: 0.9865809 CV score: 65.47438
## Fixed bandwidth: 1.067485 CV score: 67.36285
## Fixed bandwidth: 1.017484 CV score: 64.269
## Fixed bandwidth: 1.00568 CV score: 64.40019
## Fixed bandwidth: 1.024779 CV score: 64.37273
## Fixed bandwidth: 1.012975 CV score: 64.27412
## Fixed bandwidth: 1.02027 CV score: 64.2926
## Fixed bandwidth: 1.015762 CV score: 64.26456
## Fixed bandwidth: 1.014697 CV score: 64.26575
## Fixed bandwidth: 1.016419 CV score: 64.26533
## Fixed bandwidth: 1.015355 CV score: 64.26466hasil.GWPR.exponential.fix <- gwr.basic(PDN_AREA ~ KMIS + RLS + TPT + PDRB + AHH + PKD,
data = data.sp.GWPR, bw = bwd.GWPR.exponential.fix, kernel = "exponential", adaptive=F)
AIC <- c(hasil.GWPR.bisquare.ad$GW.diagnostic$AIC,
hasil.GWPR.gaussian.ad$GW.diagnostic$AIC,
hasil.GWPR.exponential.ad$GW.diagnostic$AIC,
hasil.GWPR.bisquare.fix$GW.diagnostic$AIC,
hasil.GWPR.gaussian.fix$GW.diagnostic$AIC,
hasil.GWPR.exponential.fix$GW.diagnostic$AIC)
R2 <- c(hasil.GWPR.bisquare.ad$GW.diagnostic$gw.R2,
hasil.GWPR.gaussian.ad$GW.diagnostic$gw.R2,
hasil.GWPR.exponential.ad$GW.diagnostic$gw.R2,
hasil.GWPR.bisquare.fix$GW.diagnostic$gw.R2,
hasil.GWPR.gaussian.fix$GW.diagnostic$gw.R2,
hasil.GWPR.exponential.fix$GW.diagnostic$gw.R2)
AICc <- c(hasil.GWPR.bisquare.ad$GW.diagnostic$AICc,
hasil.GWPR.gaussian.ad$GW.diagnostic$AICc,
hasil.GWPR.exponential.ad$GW.diagnostic$AICc,
hasil.GWPR.bisquare.fix$GW.diagnostic$AICc,
hasil.GWPR.gaussian.fix$GW.diagnostic$AICc,
hasil.GWPR.exponential.fix$GW.diagnostic$AICc)
R2adj <- c(hasil.GWPR.bisquare.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.gaussian.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.exponential.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.bisquare.fix$GW.diagnostic$gwR2.adj,
hasil.GWPR.gaussian.fix$GW.diagnostic$gwR2.adj,
hasil.GWPR.exponential.fix$GW.diagnostic$gwR2.adj)
RSS <- c(hasil.GWPR.bisquare.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.gaussian.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.exponential.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.bisquare.fix$GW.diagnostic$RSS.gw,
hasil.GWPR.gaussian.fix$GW.diagnostic$RSS.gw,
hasil.GWPR.exponential.fix$GW.diagnostic$RSS.gw)Perbandingan
tabel <- cbind(RSS,R2,R2adj,AIC,AICc)
rownames(tabel) <- c("Adaptive Bisquare", "Adaptive Gaussian",
"Adaptive Exponential", "Fixed Bisquare",
"Fixed Gaussian", "Fixed Exponential")
data.frame(tabel)bwd.GWPR.exponential.fix## [1] 1.015762Output di atas menunjukkan bahwa fungsi pembobot Fixed kernel Exponential memiliki nilai error dan AIC paling rendah dengan nilai R2 dan R2 adj paling tinggi, sehingga fungsi pembobot yang terpilih adalah Fixed Kernel Exponential. Fungsi pembobot Fixed Exponential memiliki nilai CV paling minimum yaitu 64,2646 dengan nilai banwidth yaitu 1,01535.
Dengan demikian, matriks pembobot spasial dapat ditentukan dengan menggunakan fungsi pembobot Fixed Exponential. Matriks pembobot spasial ini akan digunakan untuk melakukan pemodelan GWPR untuk menduga nilai parameter model. Nilai dugaan parameter dalam pemodelan GWPR akan berbeda untuk tiap provinsi tetapi sama (konstan) untuk tiap tahun.
Pendugaan Parameter
#Penduga parameter
parameter.GWPR <- as.data.frame(hasil.GWPR.exponential.fix$SDF[1:34,1:7])[-c(8,9)]
parameter.GWPR %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything())Model GWPR dengan Jumlah Pidana per Jumlah Penduduk per Luas Wilayah Tertinggi
parameter_R2 <- as.data.frame(hasil.GWPR.exponential.fix$SDF[1:34,c(1:7,27)])[-9:-10]
parameter_R2 %>%
mutate(PROV=crime$PROV[1:34]) %>%
filter(PROV %in% c("DKI JAKARTA", "D I YOGYAKARTA", "SULAWESI UTARA")) %>%
select(PROV, everything())R2 setiap provinsi
parameter_R2 %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, Local_R2) %>%
arrange(-Local_R2)Signifikasi parameter (p-value)
p_value <- gwr.t.adjust(hasil.GWPR.bisquare.ad)$results$p
p_value <- as.data.frame(ifelse(p_value<=0.05, "Signifikan", "Tidak"))[1:34,]
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything())Filter Provinsi yang mempunyai signifikansi pada peubah Kemiskinan
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$KMIS_p)Filter Provinsi yang mempunyai signifikansi peubah Rata-rata Lama Sekolah
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$RLS_p)Filter Provinsi yang mempunyai signifikansi peubah Tingkat Pengangguran Terbuka
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$TPT_p)Peubah TPT tidak signifikan terhadap Jumlah Pidana per Penduduk per Luas Wilayah
Filter Provinsi yang mempunyai signifikansi peubah PDRB
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$PDRB_p)Filter Provinsi yang mempunyai signifikansi Angka Harapan Hidup
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$AHH_p)Filter Provinsi yang mempunyai signifikansi Pengeluaran Per Kapita
p_value %>%
mutate(PROV=crime$PROV[1:34]) %>%
select(PROV, everything()) %>%
filter("Signifikan" == p_value$PKD_p)Uji Kesesuaian Model
Uji kesesuaian model dilakukan untuk mengetahui apakah terdapat perbedaan yang signifikan antara model regresi global dan model GWPR, dengan hipotesisnya sebagai berikut:
\(H_0:β_k (u_i,v_i)=β_k\) (Tidak ada perbedaan antara model regresi panel dengan GWPR)
\(H_1:β_k (u_i,v_i )≠β_k\) (Ada perbedaan antara model regresi panel dengan GWPR)
# Uji Kesesuaian Model
# Uji F
rss_rem <- sum(rem_gls$residuals^2)
rss_gwpr <- hasil.GWPR.exponential.fix$GW.diagnostic$RSS.gw
df_rem <- rem_gls$df.residual
df_gwpr <- hasil.GWPR.exponential.fix$GW.diagnostic$edf
Fhit=(rss_gwpr/df_gwpr)/(rss_rem/df_rem)
Fhit## [1] 0.2912719pf(Fhit,df_gwpr,df_rem)## [1] 2.512821e-11Berdasarkan hasil uji, didapatkan nilai F-hitung sebesar 0.29 << F-tabel (1.318) dengan p-value sebesar 2.51e-11. Oleh karena itu, H0 ditolak sehingga dapat disimpulkan bahwa cukup bukti untuk menyatakan ada perbedaan antara model regresi panel dengan GWPR.
comparison <- data.frame(R2=c(hasil.GWPR.bisquare.fix$GW.diagnostic$gw.R2,summary(rem_gls)$r.squared[1]),
R2adj=c(hasil.GWPR.bisquare.fix$GW.diagnostic$gwR2.adj,summary(rem_gls)$r.squared[2]),
RSS=c(hasil.GWPR.bisquare.fix$GW.diagnostic$RSS.gw,sum(summary(rem_gls)$residuals^2)))
rownames(comparison) <- c("GWPR", "REM Two Ways")
comparisonDari hasil di atas, model GWPR memiliki nilai R^2 dan R^2 adj terbesar dan RSS terkecil artinya model GWPR memiliki goodness of fit yang lebih baik daripada model REM sehingga model GWPR terpilih menjadi model terbaik. Selanjutnya, akan dilihat model GWPR beserta signifikansi peubah di 3 Provinsi dengan Jumlah Pidana per Jumlah Penduduk per Luas Wilayah Tertinggi.
Kesimpulan dan Saran
Dari hasil analisis yang telah dilakukan, dapat disimpulkan bahwa peubah-peubah yang berpengaruh terhadap Jumlah Pidana per penduduk per luas wilayah adalah Tingkat Kemiskinan, Rata-rata lama sekolah, PDRB per Kapita, Angka Harapan Hidup, Pengeluaran per Kapita Disesuaikan. Nilai R2 terendah adalah sebesar 0.88 pada Provinsi Kepulauan Riau yang menunjukan bahwa model telah cukup baik dalam memprediksi Jumlah pidana per penduduk per luas wilayah. Hampir semua model provinsi tidak memiliki peubah yang signifikan terhadap respons menunjukan bahwa model buruk dalam interpretasi pendugaan parameter setiap peubah.
Pada analisis disebutkan bahwa faktor yang berpengaruh adalah faktor individu dan faktor waktu. Oleh karena itu, disarankan untuk penelitian selanjutnya menggunakan metode yang digunakan adalah Geographically and Temporally Weighted Regression dengan mempertimbangkan pengaruh waktu. Pengaruh peubah-peubah yang digunakan berpengaruh secara signifikan umumnya pada Provinsi DKI Jakarta dan DI Yogyakarta. Oleh karena itu, disarankan kepada pemerintah untuk mengupayakan langkah-langkah efektif dalam menangani tingkat kriminalitas di kedua provinsi tersebut. Cara yang dapat dilakukan adalah dengan mengentas Kemiskinan, meningkatkan mutu pendidikan, meningkatkan Pendapatan dan Pengeluaran per Kapita, hingga memperbaiki kualitas hidup masyarakat.