Metode peramalan kausal (regresi berganda) adalah metode analisis untuk meramalkan nilai suatu variabel berdasarkan pola hubungan variabel tersebut dengan variabel lainnya yang mempengaruhi. Oleh karena metode peramalan ini didasarkan pada pola hubungan antarvariabel, maka bentuk hubungan variabel yang diramalkan dengan variabel yang diduga mempengaruhi merupakan hal krusial dalam metode analisis ini.
Dalam peramalan kausal, terdapat sejumlah asumsi yang digunakan. Asumsi-asumsi yang digunakan adalah sebagai berikut ini.
Variabel yang akan diramal menunjukkan suatu hubungan sebab akibat dengan satu atau lebih variabel yang lain.
Bersifat kuantitatif dan terdapat hubungan sebab akibat. Variabel yang digunakan terdiri dari variabel bebas dan terikat, bentuk hubungan fungsional dapat berupa time series maupun cross section.
Memenuhi asumsi model regresi. Secara khusus, dalam model regresi terdapat sejumlah asumsi seperti berikut ini.
Memenuhi kriteria estimasi hubungan, yaitu kriteria least square.
Memenuhi keakuratan model, yang terdiri dari kesalahan baku dan proporsi variansi yang dijelaskan model.
Memenuhi keberartian persamaan regresi.
Secara umum, analisis regresi terbagi menjadi analisis regresi sederhana dan analisis regresi berganda. Persamaan regresi sederhana dapat dinyatakan dalam persamaan \(Y = a + bX\). Sedangkan regresi berganda merupakan perluasan dari regresi sederhana. Regresi berganda menyatakan hubungan antara satu variabel dependen dengan satu atau lebih variabel independen. Persamaan regresi berganda dapat dinyatakan seperti berikut:
\(Y' = a + b_1X_1 + b_2X_2 + \ldots + b_kX_k\),
dengan \(Y'\) adalah nilai prediksi dari variabel; \(a\) merupakan koefisien yang ditentukan berdasarkan data sampel; dan \(X_1, X_2, \ldots, X_k\) merupakan variabel prediktor.
Tahapan penting dalam analisis regresi adalah seleksi variabel prediktor. Prosedur dalam seleksi variabel prediktor adalah sebagai berikut.
Seleksi variabel prediktor dapat dilakukan dalam berbagai cara, antara lain adalah plot variabel bebas terhadap variabel terikat, hitung inter-korelasi, analisis faktor melalui prosedur lain diantaranya “stepwise”. Prosedur stepwise terdiri dari prosedur forward dan prosedur backward. a. Backward Elimination Metode ini dimulai dengan memasukkan semua variabel terlebih dahulu, kemudian dilakukan analisis. Variabel yang tidak layak masuk dalam regresi akan dikeluarkan satu per satu. b. Forward Elimination Metode ini dilakukan dengan memasukkan variabel independen satu per satu.
Dalam regresi terdapat permasalahan multicollinierarity. Multikolinearity adalah suatu masalahperhitungan yang muncul jika kondisi-kondisi seperti ini terjadi karena: (1) dua variabel bebas berkorelasi sempurna; (2) dua variabel bebas hampir berkorelasi sempurna; (3) kombinasi linear dari beberapa variabel bebas berkorelasi sempurna dengan variabel bebas lainnya; (4) kombinasi linear dari suatu sub-himpunan variabel bebas berkorelasi sempurna atau mendekati sempurna dengan suatu kombinasi linear dari sub-himpunan variabel bebas lainnya.
Multikolinearitas harus dihindari dengan alasan berikut.
Multikolinearitas dapat dihindari dengan beberapa alternatif cara yang diantaranya adalah sebagai berikut.
Dalam membuat persamaan regresi, variabel independen yang bersifat kualitatif (data nominal atau ordinal) tetap dapat digunakan. Hanya saja data yang bersifat kualitatif tidak dapat digunakan secara langsung. Variabel tersebut harus diubah menjadi variabel dummy. Variabel dummy dinyatakan dalam nilai 1 dan 0, nilai 1 untuk nilai yang diinginkan, dan nilai 0 untuk nilai yang tidak diinginkan. Suatu variabel dummy equivalen dengan variabel regresor baru. Bila digunakan variabel kualitatif dengan n kriteria, maka gunakan n variabel dummy untuk menunjukkan n kriteria yang berbeda.
Pada pratikum kali ini akan digunakan data Provinsi Jawa Barat pada tahun 2006. Provinsi Jawa Barat memiliki 25 kabupaten dan kota. Variabel-variabel yang digunakan adalah:
Dari 10 variabel diatas, jumlah PDRB tahun 2006 merupakan variabel dependen sedangkan variabel lainnya digunakan sebagai variabel prediktor. Untuk mengetahui apa saja variabel yang mempengaruhi jumlah PDRB pada tahun 2006 di Provinsi Jawa Barat, maka digunakan analisis regresi berganda. Metode yang digunakan dalam analisis regresi berganda kali ini adalah metode Stepwise. Secara umum, cara kerja dari metode tersebut adalah dengan memasukkan variabel bebas satu per satu dan variabel yang tidak memiliki korelasi dengan variabel dependen dapat dikeluarkan.
Adapun tujuan dari pengolahan data tersebut adalah untuk membuat model regresi yang mengestimasikan besaran PDRB, berdasarkan variabel-variabel lainnya yang memiliki hubungan dengan variabel PDRB.
Data bisa didapatkan pada link berikut DATA. Berikut adalah data yang digunakan.
data <- read.csv("data.csv")
data
## kab_kota jml_pnddk jml_prod_padi jml_trnk_klr
## 1 kabupaten bogor 4100934 411511 17304
## 2 sukabumi 2224993 677328 82773
## 3 cianjur 2098644 686619 8699
## 4 kabupaten bandung 4263934 604540 13164
## 5 garut 2321070 617922 3864
## 6 kabupaten tasikmalaya 1693479 645065 0
## 7 ciamis 1542661 552347 3334
## 8 kuningan 1096848 331802 11863
## 9 kabupaten cirebon 2107918 430794 531
## 10 majalengka 1191490 531572 105
## 11 sumedang 1067361 383070 469
## 12 indramayu 1760286 1137958 35604
## 13 subang 1421973 962898 20743
## 14 purwakarta 770660 176960 1441
## 15 karawang 1985574 919843 0
## 16 kabupaten bekasi 1953380 529093 157852
## 17 kota bogor 844778 9219 0
## 18 sukabumi 287760 22844 9380
## 19 kota bandung 2315895 24529 8999
## 20 kota cirebon 281089 3496 0
## 21 kota bekasi 1994850 10406 12222
## 22 depok 1373860 7668 469
## 23 cimahi 493693 2930 0
## 24 kota tasikmalaya 594158 73773 939
## 25 banjar 173576 33030 0
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_kclmngh
## 1 1026 41 3346
## 2 1233 85 600
## 3 423 45 0
## 4 531 80 291
## 5 1055 81 129
## 6 13 2 0
## 7 10 2 0
## 8 10 3 0
## 9 573 23 170
## 10 69 8 0
## 11 1692 3 150
## 12 0 0 0
## 13 5 1 0
## 14 230 11 244
## 15 95 5 662
## 16 242 4 988
## 17 36 3 687
## 18 249 19 0
## 19 1324 99 129
## 20 57 7 0
## 21 232 15 139
## 22 56 7 102
## 23 0 0 0
## 24 32 6 0
## 25 205 25 0
## jml_unt_ush_nd_bsr jml_htl_akomodasi PDRB
## 1 26 117 38182120.0
## 2 1 95 11324257.0
## 3 0 103 10776519.0
## 4 5 47 25494164.0
## 5 8 93 13697884.0
## 6 0 8 7253242.0
## 7 0 218 9247507.0
## 8 0 42 4573373.0
## 9 3 13 9760061.0
## 10 0 9 5129025.5
## 11 1 19 7048211.0
## 12 0 22 23591254.0
## 13 0 76 9063475.0
## 14 3 19 8531292.0
## 15 5 22 25804278.0
## 16 6 8 57017160.0
## 17 4 49 5391919.0
## 18 0 32 2402017.5
## 19 3 253 34792184.0
## 20 0 44 6840256.0
## 21 3 16 19226332.0
## 22 1 5 7541666.0
## 23 0 3 7227777.5
## 24 0 32 4617522.0
## 25 0 8 973962.8
Berdasarkan studi kasus tersebut, maka penyelesaian yang dapat dilakukan terdiri dari beberapa tahapan yang diantaranya adalah: (1) merumuskan masalah; (2) mengolah data menggunakan R; (3) menganalisis dan menginterpretasikan hasil analisis
Masalah yang akan diteliti adalah sebagai berikut.
# summary function
summary(data)
## kab_kota jml_pnddk jml_prod_padi jml_trnk_klr
## Length:25 Min. : 173576 Min. : 2930 Min. : 0
## Class :character 1st Qu.: 844778 1st Qu.: 24529 1st Qu.: 105
## Mode :character Median :1542661 Median : 411511 Median : 3334
## Mean :1598435 Mean : 391489 Mean : 15590
## 3rd Qu.:2098644 3rd Qu.: 617922 3rd Qu.: 12222
## Max. :4263934 Max. :1137958 Max. :157852
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_kclmngh
## Min. : 0.0 Min. : 0 Min. : 0.0
## 1st Qu.: 32.0 1st Qu.: 3 1st Qu.: 0.0
## Median : 205.0 Median : 7 Median : 102.0
## Mean : 375.9 Mean :23 Mean : 305.5
## 3rd Qu.: 531.0 3rd Qu.:25 3rd Qu.: 244.0
## Max. :1692.0 Max. :99 Max. :3346.0
## jml_unt_ush_nd_bsr jml_htl_akomodasi PDRB
## Min. : 0.00 Min. : 3.00 Min. : 973963
## 1st Qu.: 0.00 1st Qu.: 13.00 1st Qu.: 6840256
## Median : 1.00 Median : 32.00 Median : 9063475
## Mean : 2.76 Mean : 54.12 Mean :14220298
## 3rd Qu.: 3.00 3rd Qu.: 76.00 3rd Qu.:19226332
## Max. :26.00 Max. :253.00 Max. :57017160
Dari descriptive statistics tersebut, perlu dilihat apakah terdapat missing value. Tabeldescriptive statistics ini diperuntukkan untuk melihat gambaran karakteristik data yang dimiliki, dalam hal ini rata-rata dan standar deviasi dari masing-masing variabelnya.
Cek terlebih dahulu missing values setiap kolom.
# check missing values in data
missing_values <- colSums(is.na(data))
# menampilkan missing values
print("Missing values:")
## [1] "Missing values:"
print(missing_values)
## kab_kota jml_pnddk jml_prod_padi
## 0 0 0
## jml_trnk_klr jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr
## 0 0 0
## jml_unt_ush_nd_kclmngh jml_unt_ush_nd_bsr jml_htl_akomodasi
## 0 0 0
## PDRB
## 0
Tidak ada missing values
Kali ini kita akan melihat korelasi seluruh variabel.
# Menghitung matriks korelasi
numeric_vars <- sapply(data, is.numeric) # Menentukan variabel numerik
data_numeric <- data[, numeric_vars] # Subset data untuk variabel numerik saja
cor_matrix <- cor(data_numeric)
Notes: Output bisa dilihat di halaman terpisah untuk lebih mudah. Namun, untuk melihat dalam bentuk HTML bisa dilihat di bawah ini.
library(knitr)
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.2
# Menggunakan kable untuk membuat tabel dasar
kable_cor <- kable(cor_matrix, format = "html", caption = "Matriks Korelasi Variabel Numerik") %>%
# Menambahkan styling dengan kableExtra
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
# Menampilkan tabel
kable_cor
jml_pnddk | jml_prod_padi | jml_trnk_klr | jml_tngkrj_inc_kclmngh | jml_tngkrj_inc_bsr | jml_unt_ush_nd_kclmngh | jml_unt_ush_nd_bsr | jml_htl_akomodasi | PDRB | |
---|---|---|---|---|---|---|---|---|---|
jml_pnddk | 1.0000000 | 0.4438861 | 0.2185172 | 0.4335774 | 0.5885246 | 0.5657576 | 0.6589290 | 0.3434810 | 0.6354552 |
jml_prod_padi | 0.4438861 | 1.0000000 | 0.2742949 | 0.0249897 | 0.0448039 | 0.0643892 | 0.0554494 | 0.0622070 | 0.2643122 |
jml_trnk_klr | 0.2185172 | 0.2742949 | 1.0000000 | 0.1105341 | 0.0995877 | 0.2677504 | 0.1289454 | -0.0308428 | 0.6607597 |
jml_tngkrj_inc_kclmngh | 0.4335774 | 0.0249897 | 0.1105341 | 1.0000000 | 0.6766148 | 0.3108878 | 0.3769838 | 0.3955271 | 0.2596851 |
jml_tngkrj_inc_bsr | 0.5885246 | 0.0448039 | 0.0995877 | 0.6766148 | 1.0000000 | 0.1542704 | 0.2852532 | 0.5477687 | 0.2785888 |
jml_unt_ush_nd_kclmngh | 0.5657576 | 0.0643892 | 0.2677504 | 0.3108878 | 0.1542704 | 1.0000000 | 0.9400460 | 0.1528709 | 0.5695537 |
jml_unt_ush_nd_bsr | 0.6589290 | 0.0554494 | 0.1289454 | 0.3769838 | 0.2852532 | 0.9400460 | 1.0000000 | 0.2012609 | 0.5753913 |
jml_htl_akomodasi | 0.3434810 | 0.0622070 | -0.0308428 | 0.3955271 | 0.5477687 | 0.1528709 | 0.2012609 | 1.0000000 | 0.2355098 |
PDRB | 0.6354552 | 0.2643122 | 0.6607597 | 0.2596851 | 0.2785888 | 0.5695537 | 0.5753913 | 0.2355098 | 1.0000000 |
Tabel matriks korelasi menjelaskan korelasi masing-masing variabel dengan variabel lainnya. Pada baris diagonal terdapat angka 1,000, hal ini berarti korelasi antarvariabel dengan dirinya sendiri sangat kuat. Jika nilai korelasi > 0,5 maka hubungannya cukup kuat, sedangkan bila < 0,5 hubungannya lemah. Apabila terdapat dua variabel independen dengan korelasi sangat kuat, maka perlu berhati-hati karena dapat memunculkan fenomena multikolinearitas. Akan lebih baik apabila hanya salah satu variabel yang digunakan.
Tabel tersebut juga menunjukkan korelasi antar variabel, baik itu antara variabel dependen dengan variabel independen dan juga antar variabel independennya. Variabel dependen (PDRB) memiliki korelasi paling besar dengan variabel jumlah penduduk yaitu sebesar 0,635 dan dengan jumlah ternak yang keluar yaitu 0,661. Berdasarkan korelasi ini ada kemungkinan kedua variabel independen dapat menjelaskan variabel dependen (PDRB). Akan tetapi hal tersebut harus diperkuat melalui pengujian selanjutnya.
Berdasarkan tabel tersebut dapat dilihat bahwa adanya korelasi yang kuat antar variabel independen yaitu variabel jumlah unit usaha menengah kecil dan jumlah unit usaha besar dengan koefisien korelasi 0,94. Tingginya korelasi tersebut dapat mengindikasikan adanya multikolineariti antara kedua variabel tersebut. Konsekuensinya adalah hanya salah satu variabel saja yang boleh diambil dalam model.
Hal ini bisa dilihat lebih jelas jika menggunakan heatmap correlation.
# Install and load the ggcorrplot package
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.3.2
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.2
# Compute and show the result
ggcorrplot(cor_matrix, hc.order = TRUE, type = "lower",
lab = TRUE)
Prosedur default adalah prosedur forward tanpa persyaratan nilai siginifikansi. Sehingga seluruh variabel dimasukkan ke dalam model regresi.
# Menghapus kolom pertama (Kab/Kota)
data_reg <- data[,-1]
# agar output tidak menggunakan notes saintifik seperti 'e+'
options(scipen = 999)
# Melakukan multiple regression dengan seluruh variabel independen yang ada
model_all <- lm(formula = PDRB ~ . ,data = data_reg) # . artinya adalah seluruh variabel independen
summary(model_all)
##
## Call:
## lm(formula = PDRB ~ ., data = data_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9234970 -2455815 -635103 1854901 15179635
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1308508.286 3126097.230 0.419 0.681091
## jml_pnddk 6.127 3.022 2.027 0.059645 .
## jml_prod_padi -5.472 5.751 -0.951 0.355514
## jml_trnk_klr 269.620 53.129 5.075 0.000113 ***
## jml_tngkrj_inc_kclmngh -112.258 4321.453 -0.026 0.979597
## jml_tngkrj_inc_bsr -135635.250 94151.668 -1.441 0.168979
## jml_unt_ush_nd_kclmngh -12413.056 7956.540 -1.560 0.138294
## jml_unt_ush_nd_bsr 2064914.621 1075149.056 1.921 0.072789 .
## jml_htl_akomodasi 42709.476 28416.468 1.503 0.152323
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7278000 on 16 degrees of freedom
## Multiple R-squared: 0.7999, Adjusted R-squared: 0.6999
## F-statistic: 7.996 on 8 and 16 DF, p-value: 0.000236
Call: Ini menunjukkan fungsi yang dipanggil
(lm
) dan formula yang digunakan dalam model, yaitu
PDRB ~ .
, yang berarti PDRB
adalah variabel
dependen dan .
mewakili semua variabel lain dalam
data
yang digunakan sebagai variabel independen.
Residuals: Memberikan ringkasan statistik dari residu, yaitu selisih antara nilai yang diamati dari variabel dependen dan nilai yang diprediksi oleh model. Residu memiliki rentang yang sangat lebar dari minimum -9234970 hingga maksimum 15179635, yang mungkin menunjukkan variabilitas yang tinggi dalam model atau outlier.
Coefficients:
Estimate
. Intercept (konstanta) adalah
1.309e+06, yang merupakan nilai prediksi untuk PDRB
ketika
semua variabel independen adalah 0.Std. Error
menunjukkan kesalahan standar dari estimasi
koefisien, yang mengukur variabilitas estimasi koefisien.t value
adalah rasio dari estimasi koefisien terhadap
kesalahan standarnya. Nilai t yang tinggi (positif atau negatif)
menunjukkan bahwa koefisien tersebut secara statistik signifikan berbeda
dari 0.Pr(>|t|)
menunjukkan p-value dari uji t untuk
hipotesis bahwa koefisien tersebut adalah 0. Nilai yang rendah
menunjukkan bukti yang kuat melawan hipotesis nol dan menunjukkan bahwa
variabel tersebut penting dalam model.jml_trnk_klr
memiliki tiga bintang
(***
), menunjukkan bahwa variabel ini sangat signifikan
secara statistik dalam model (p < 0.001).Signif. codes: Legenda yang menunjukkan tingkat signifikansi dari koefisien berdasarkan p-value mereka.
Residual standard error: Ini adalah perkiraan standar deviasi dari residu dan memberikan ukuran rata-rata seberapa jauh data yang diamati jatuh dari garis regresi yang diestimasi.
Multiple R-squared: Koefisien determinasi, yang
mengukur seberapa baik model regresi menjelaskan variabilitas dari
variabel dependen. Nilai ini adalah 0.7999, yang berarti sekitar 79.99%
dari variabilitas PDRB
dapat dijelaskan oleh
variabel-variabel independen dalam model.
Adjusted R-squared: Versi penyesuaian dari R-squared yang mempertimbangkan jumlah variabel dalam model. Nilai ini lebih kecil daripada R-squared dan merupakan ukuran yang lebih baik untuk membandingkan model dengan jumlah prediktor yang berbeda.
F-statistic: Uji statistik untuk mengevaluasi apakah setidaknya satu koefisien regresi dalam model secara signifikan berbeda dari 0. Dengan F-statistik sebesar 7.996 dan p-value yang sangat rendah (0.000236), kita memiliki bukti yang kuat untuk menyimpulkan bahwa model secara keseluruhan signifikan secara statistik.
Persamaan regresi: \(PDRB = 1308508.286 + (6.127 * Jumlah penduduk) - (5.472 * Jumlah produksi padi) +(269.620 * Jumlah ternak yang keluar) - (112.258 * Jumlah tenaga kerja industri kecil menengah) -(135635.250 * Jumlah tenaga kerja industri besar) - (12413.056 * Jumlah unit usaha industri kecil menengah) +(2064941.621 * Jumlah unit usaha industri besar) + (42709.476 * Jumlah hotel dan akomodasi)\)
Test Heteroskedasticity: Untuk menguji
heteroskedastisitas dalam model regresi, salah satu metode yang paling
umum digunakan adalah uji Breusch-Pagan. Uji ini memeriksa apakah
varians residu (kesalahan prediksi) berhubungan dengan varians dari
variabel independen. Di R, kita dapat melakukan uji Breusch-Pagan
menggunakan fungsi bptest()
dari paket
lmtest
.
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# melakukan uji Breusch-Pagan
bptest_all <- bptest(model_all)
# Menampilkan hasil uji Breusch-Pagan
bptest_all
##
## studentized Breusch-Pagan test
##
## data: model_all
## BP = 8.9261, df = 8, p-value = 0.3486
Asumsi model regresi, variansi distribusi variabel dependen semuanya sama satu dengan yang lain (Homosedasticity). Hipotesis nol, menyatakan bahwa variansi konstan (homo). Hal tersebut dapat dilihat dari nilai-p, apabila kurang dari 0,05 artinya hipotesa nol ditolak, artinya terdapat heteroskedasticity. Pada kasus ini, p bernilai 0.3486 artinya tidak terdapat heteroskedasticity.
Tes VIF: Uji Variance Inflation Factor (VIF) digunakan untuk mendeteksi adanya multikolinearitas di antara variabel independen dalam model regresi. Multikolinearitas terjadi ketika variabel independen berkorelasi erat satu sama lain, yang dapat menyebabkan masalah dalam estimasi koefisien regresi. VIF mengukur seberapa banyak varians dari suatu estimasi koefisien regresi meningkat karena adanya kolinearitas.
library(car)
## Warning: package 'car' was built under R version 4.3.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.2
# VIF test pada model
vif_model_all <- vif(model_all)
# Menampilkan hasil VIF
print(vif_model_all)
## jml_pnddk jml_prod_padi jml_trnk_klr
## 4.401629 1.778473 1.514696
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_kclmngh
## 2.062389 3.782425 13.584754
## jml_unt_ush_nd_bsr jml_htl_akomodasi
## 15.068421 1.506562
Hasil di atas menunjukkan nilai VIF yang dihasilkan oleh masing-masing variabel. Pada kasus ini, terdapat multikolinearitas yakni variabel jumlah unit usaha besar dan variabel jumlah unit usaha kecil menengah. Karena kedua variabel tersebut memiliki nilai VIF yang berbeda signifikan dengan nilai VIF variabel lainnya.
Distribusi Residual
# library untuk density plot
library(ggplot2)
# Menghitung residu dari model regresi
model_all_residuals = model_all$residuals
# Membuat density plot dengan histogram dari residu
ggplot(data.frame(model_all_residuals), aes(x=model_all_residuals)) +
geom_histogram(aes(y=..density..), binwidth = diff(range(model_all_residuals))/30, fill="skyblue", alpha=0.5) +
geom_density(color="blue", fill="blue", alpha=0.3) +
labs(title="Density and Histogram Plot of Residuals", x="Residuals", y="Density")
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Garis kurva density yang menyerupai bentuk lonceng menandakan bahwa residu cenderung terdistribusi secara normal, yang merupakan asumsi penting dalam regresi linear. Namun, ada beberapa penyimpangan dari kurva normal sempurna, seperti terlihat dari ekor pada kedua sisi kurva.
Untuk menguji normalitas secara langsung, mari kita Q-Q (Quantile) plot
Q-Q Plot
# Plot the residuals
qqnorm(model_all_residuals)
# Plot the Q-Q line
qqline(model_all_residuals)
Grafik menunjukkan persebaran data residu observasi terhadap inverse normal. Apabila jarak plot dengan garis semakin kecil, artinya model regresi semakin baik.
Di sebagian besar jangkauan data, poin-poin tampak mengikuti garis referensi yang mengindikasikan bahwa residu memiliki distribusi yang dekat dengan distribusi normal. Namun, tetap ada indikasi outliers melihat bahwa masih ada poin-poin yang jauh dari garis.
Melihat scatter plot hasil prediksi
# melakukan predict pada train data (seluruh data di awal tadi)
predicted_values <- predict(model_all)
# Buat dataframe dari nilai asli dan prediksi
comparison_df <- data.frame(Actual = data_reg$PDRB, Predicted = predicted_values)
# Buat scatter plot
ggplot(comparison_df, aes(x = Actual, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") + # Garis y=x untuk referensi
labs(title = "Scatter Plot of Actual vs. Predicted", x = "Actual Values", y = "Predicted Values") +
theme_minimal()
Grafik menunjukkan persebaran data hasil regesi dibandingkan dengan data hasil observasi. Semakin baik apabila hasil plot membentuk garis diagonal (45 derajat).
Kesimpulan
Dari hasil post estimasi di atas, bisa disimpulkan bahwa model kurang baik karena adanya multikolinearitas (sesuai hasil VIF). Maka dari itu, kita akan mencoba menghapus salah satu variabel yang menyebabkan hal tersebut.
Berdasarkan prosedur regresi seluruh variabel yang dilakukan sebelumnya, diindikasikan terdapat multikorelasi antara variabel jumlah usaha besar dan usaha kecil menengah. Maka dari itu perlu dihapuskan salah satu variabelnya. Pada kasus ini dihapus variael jumlah usaha kecil menengah:
# Melakukan multiple regression dengan seluruh variabel kecuali jumlah usaha kecil
model_minus_one <- lm(formula = PDRB ~ jml_pnddk + jml_prod_padi + jml_trnk_klr + jml_tngkrj_inc_kclmngh + jml_tngkrj_inc_bsr + jml_unt_ush_nd_bsr + jml_htl_akomodasi,
data = data_reg)
summary(model_minus_one)
##
## Call:
## lm(formula = PDRB ~ jml_pnddk + jml_prod_padi + jml_trnk_klr +
## jml_tngkrj_inc_kclmngh + jml_tngkrj_inc_bsr + jml_unt_ush_nd_bsr +
## jml_htl_akomodasi, data = data_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15192412 -2709949 -765396 1738084 15871239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1359116.425 3255094.985 0.418 0.681513
## jml_pnddk 6.382 3.143 2.031 0.058226 .
## jml_prod_padi -5.008 5.980 -0.837 0.414012
## jml_trnk_klr 227.333 47.583 4.778 0.000175 ***
## jml_tngkrj_inc_kclmngh -798.012 4476.680 -0.178 0.860626
## jml_tngkrj_inc_bsr -86476.279 92388.444 -0.936 0.362377
## jml_unt_ush_nd_bsr 530636.371 452445.761 1.173 0.257041
## jml_htl_akomodasi 35116.523 29153.423 1.205 0.244883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7579000 on 17 degrees of freedom
## Multiple R-squared: 0.7695, Adjusted R-squared: 0.6746
## F-statistic: 8.106 on 7 and 17 DF, p-value: 0.0002167
Coefficients (Koefisien):
Intercept: Nilai PDRB
ketika semua
variabel independen adalah 0 adalah sekitar 1,359,116.425.
jml_pnddk
: Koefisien positif (sekitar 6.382)
menunjukkan bahwa dengan setiap peningkatan satu unit pada ‘jumlah
penduduk’, PDRB
meningkat sekitar 6.382, meskipun ini tidak
signifikan secara statistik (p-value > 0.05).
jml_prod_padi
: Koefisien negatif (sekitar -5.008)
mengindikasikan bahwa peningkatan ‘jumlah produksi padi’ cenderung
dikaitkan dengan penurunan kecil pada PDRB
, tapi ini juga
tidak signifikan secara statistik.
jml_trnk_klr
: Setiap peningkatan satu unit pada
‘jumlah ternak yang keluar’ dikaitkan dengan peningkatan
PDRB
sekitar 227.333, dan ini signifikan secara statistik
(p-value < 0.01).
…dan seterusnya untuk variabel-variabel lainnya. Penting untuk mencatat bahwa koefisien yang signifikan secara statistik memberikan indikasi hubungan yang lebih meyakinkan antara variabel independen dengan variabel dependen.
Signif. codes: Menunjukkan level signifikansi
dari koefisien, dengan ***
berarti p < 0.001,
**
untuk p < 0.01, dan *
untuk p <
0.05.
Residual standard error: Sebuah perkiraan dari standar deviasi residu, menunjukkan seberapa jauh data yang diamati menyimpang dari model yang diestimasi.
Multiple R-squared: Menunjukkan bahwa model
menjelaskan sekitar 76.95% dari variabilitas dalam
PDRB
.
Adjusted R-squared: Penyesuaian R-squared untuk jumlah variabel dalam model menunjukkan nilai sekitar 67.46%, yang merupakan indikator lebih baik dari kecocokan model karena mempertimbangkan jumlah prediktor.
F-statistic: Mengukur keseluruhan signifikansi model, dengan nilai F sekitar 8.106 dan p-value sangat rendah (0.0002167), menunjukkan bahwa model secara keseluruhan signifikan secara statistik.
Dengan menggunakan koefisien estimasi dari output, persamaan regresi linier dapat ditulis sebagai:
\[ PDRB = 1359116.425 + 6.382 * (jumlah penduduk) - 5.008 * (jumlah produksi padi) +227.333 * (jumlah ternak yang keluar) - 798.012 * (jumlah tenaga kerja industri kecil menengah) -86476.279 * (jumlah tenaga kerja industri besar) + 360336.319 * (jumlah unit usaha industri kecil menengah) +35116.523 * (jumlah hotel dan akomodasi) \]
Memeriksa kembali multikolinearitas
# VIF test pada model
vif_model_minus_one <- vif(model_minus_one)
# Menampilkan hasil VIF
print(vif_model_minus_one)
## jml_pnddk jml_prod_padi jml_trnk_klr
## 4.388750 1.773715 1.120449
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_bsr
## 2.041051 3.358771 2.460902
## jml_htl_akomodasi
## 1.462369
Prosedur forward digunakan untuk seleksi variabel dengan cara menambahkan variabel satu per satu (dengan mensyaratkan level signifikansi variabel saat dimasukkan ke dalam analisis regresi).
# library untuk stepwise
library(olsrr)
## Warning: package 'olsrr' was built under R version 4.3.3
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
##
## rivers
# Menerapkan stepwise forward regression dengan menggunakan model_all sebegai model yang diiterasi
stepwise_forward <- ols_step_forward_p(model_all, data = data_reg, p_val = 0.05)
# Menampilkan hasil model setelah regresi stepwise forward
stepwise_forward
##
##
## Stepwise Summary
## ---------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## ---------------------------------------------------------------------------
## 0 Base Model 894.036 896.474 820.769 0.00000 0.00000
## 1 jml_trnk_klr 881.692 885.348 808.568 0.43660 0.41211
## 2 jml_pnddk 868.769 873.644 797.870 0.68984 0.66165
## ---------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------------------------
## R 0.831 RMSE 7249475.227
## R-Squared 0.690 MSE 59721467129518.180
## Adj. R-Squared 0.662 Coef. Var 54.345
## Pred R-Squared 0.151 AIC 868.769
## MAE 5220740.445 SBC 873.644
## --------------------------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -----------------------------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -----------------------------------------------------------------------------------------
## Regression 2922282249482486.500 2 1461141124741243.250 24.466 0.0000
## Residual 1313872276849400.000 22 59721467129518.180
## Total 4236154526331886.500 24
## -----------------------------------------------------------------------------------------
##
## Parameter Estimates
## --------------------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## --------------------------------------------------------------------------------------------------------
## (Intercept) 302669.552 2898460.874 0.104 0.918 -5708370.393 6313709.496
## jml_trnk_klr 211.571 46.970 0.548 4.504 0.000 114.161 308.981
## jml_pnddk 6.644 1.568 0.516 4.238 0.000 3.393 9.894
## --------------------------------------------------------------------------------------------------------
Penjelasan Output
Stepwise Summary
Tabel pertama adalah ringkasan stepwise yang menunjukkan model pada setiap langkah (step) dan statistik pemilihan modelnya. Anda memiliki 3 model (base_model, jml_trnkk_klr, dan _jml_pndkk) dengan berbagai metrik.
Dari tabel ini, kita dapat melihat bahwa penambahan variabel ‘jml_trnkk_klr’ dan ‘_jml_pndkk’ meningkatkan nilai R2 dan Adj. R2 dan menurunkan nilai AIC, SBC, dan SBIC, yang menunjukkan perbaikan dalam model.
Final Model Output
Bagian ini memberikan ringkasan statistik dari model regresi yang final.
ANOVA
Tabel ini menunjukkan Analisis Varians yang membandingkan model dengan model tanpa prediktor. Hal ini dilakukan untuk menguji apakah prediktor secara keseluruhan signifikan.
Parameter Estimates
Ini adalah tabel yang memberikan estimasi untuk parameter (koefisien) dari model regresi.
Interpretasi Umum
Model regresi yang dihasilkan dari stepwise forward regression tampaknya memperbaiki fitur statistik dengan menambahkan dua prediktor, ‘jml_trnkk_klr’ dan ‘_jml_pndkk’, karena keduanya signifikan dan meningkatkan R-squared model. RMSE dan MAE sangat besar, yang mungkin menunjukkan kesalahan model yang besar, tetapi ini bisa juga karena skala data.
Kita juga bisa melihat secara details proses stepwise untuk setiap
penambahan variabelnya melalui parameter metrics
stepwise_forward$metrics
## step variable r2 adj_r2 aic sbc sbic mallows_cp
## 1 1 jml_trnk_klr 0.4366034 0.4121079 881.6916 885.3482 808.5681 24.052093
## 2 2 jml_pnddk 0.6898432 0.6616471 868.7689 873.6444 797.8698 5.801738
## rmse
## 1 9770640
## 2 7249475
Hasil tersebut bisa divisualisasikan, terutama dalam indikator R-Square, Adjusted R-Square, AIC (Akaike Information Criteria), dan RMSE (Root Mean Squared Error).
plot(stepwise_forward)
R-Square Plot
Adjusted R-Square Plot
Akaike Information Criteria (AIC) Plot
Root Mean Squared Error (RMSE) Plot
Secara keseluruhan, plot tersebut menunjukkan bahwa dengan menambahkan prediktor ‘jml_trnkk_klr’ dan ‘_jml_pndkk’ ke dalam model, kualitas fit meningkat, sebagaimana dicerminkan oleh peningkatan R-squared dan Adjusted R-squared serta penurunan dalam AIC dan RMSE. Ini menunjukkan bahwa kedua variabel tersebut memberikan informasi yang signifikan yang meningkatkan kemampuan prediktif model regresi.
Untuk menyimpan hasil modelnya kita lakukan perintah di bawah ini.
model_stepwise_forward <- stepwise_forward$model
summary(model_stepwise_forward)
##
## Call:
## lm(formula = paste(response, "~", paste(preds, collapse = " + ")),
## data = l)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21272484 -4300057 -523034 3645258 17199940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 302669.552 2898460.874 0.104 0.917779
## jml_trnk_klr 211.571 46.970 4.504 0.000176 ***
## jml_pnddk 6.644 1.568 4.238 0.000337 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7728000 on 22 degrees of freedom
## Multiple R-squared: 0.6898, Adjusted R-squared: 0.6616
## F-statistic: 24.47 on 2 and 22 DF, p-value: 0.000002555
Berdasarkan output di atas, persamaan regresi dari model adalah:
\[ Y = 302669.552 + 211.571 \times \text{jumlah ternak yang keluar} + 6.644 \times \text{jumlah penduduk} \]
Sekarang mari kita interpretasi koefisien tersebut:
Intercept ( \(\beta_0\) ): Nilainya adalah 302,669.552. Ini merupakan nilai prediksi untuk \(Y\) ketika semua variabel independen (\(X\)) adalah 0. Namun, interpretasi intercept bisa tidak bermakna jika 0 bukan merupakan nilai yang mungkin atau masuk akal untuk variabel independen dalam konteks masalah Anda.
Koefisien untuk jml_trnkk_klr
(
\(\beta_1\) ): Nilainya adalah 211.571.
Ini berarti bahwa dengan setiap peningkatan satu unit dalam
jml_trnkk_klr
, nilai prediksi \(Y\) meningkat sebanyak 211.571 unit, asumsi
variabel lain tetap konstan.
Koefisien untuk jml_pndkk
( \(\beta_2\) ): Nilainya adalah 6.644. Ini
berarti bahwa dengan setiap peningkatan satu unit dalam
jml_pndkk
, nilai prediksi \(Y\) meningkat sebanyak 6.644 unit, asumsi
variabel lain tetap konstan.
Koefisien-koefisien ini adalah signifikan pada level 0.05, seperti yang ditunjukkan oleh nilai p (\(p < 0.05\)) yang ditandai dengan bintang-bintang di output.
Terakhir, nilai R-squared dari 0.6898 menunjukkan bahwa sekitar 68.98% variabilitas dalam variabel dependen dijelaskan oleh model. Sedangkan nilai Adjusted R-squared adalah 0.6616, menunjukkan sekitar 66.16% dari variabilitas dijelaskan ketika menyesuaikan untuk jumlah variabel dalam model. Nilai F-statistik yang besar dan nilai p yang sangat kecil untuk F-statistik menunjukkan bahwa model secara keseluruhan memiliki signifikansi statistik yang kuat.
Tes Heteroskedasticity
# melakukan uji Breusch-Pagan
bptest_stepwise_forward <- bptest(model_stepwise_forward)
# Menampilkan hasil uji Breusch-Pagan
bptest_stepwise_forward
##
## studentized Breusch-Pagan test
##
## data: model_stepwise_forward
## BP = 6.2264, df = 2, p-value = 0.04446
Hipotesis nol, menyatakan bahwa variansi konstan (homo). Hal tersebut dapat dilihat dari nilai-p, apabila kurang dari 0,05 artinya hipotesa nol ditolak, artinya terdapat heteroskedasticity. Pada kasus ini, p bernilai 0.04446 artinya terdapat heteroskedasticity.
Tes VIF
# VIF test pada model
vif_model_stepwise_forward <- vif(model_stepwise_forward)
# Menampilkan hasil VIF
print(vif_model_stepwise_forward)
## jml_trnk_klr jml_pnddk
## 1.050144 1.050144
Hasilnya tidak terdapat nilai VIF yang terlalu besar menandakan bahwa tidak ada multikolinearitas
Distribusi Residual
# Menghitung residu dari model regresi
model_stepwise_forward_residuals = model_stepwise_forward$residuals
# Membuat density plot dengan histogram dari residu
ggplot(data.frame(model_stepwise_forward_residuals), aes(x=model_stepwise_forward_residuals)) +
geom_histogram(aes(y=..density..), binwidth = diff(range(model_stepwise_forward_residuals))/30, fill="skyblue", alpha=0.5) +
geom_density(color="blue", fill="blue", alpha=0.3) +
labs(title="Density and Histogram Plot of Residuals (Stepwise Forward)", x="Residuals", y="Density")
Grafik ini menunjukkan distribusi residu dari model. Idealnya, residu harus mendekati distribusi normal (Gaussian), yang akan terlihat seperti kurva lonceng simetris pada plot density. Tampak bahwa residu tidak sepenuhnya simetris dan terdapat beberapa lonjakan yang tinggi. Ini bisa menunjukkan adanya outlier atau bahwa model tidak sepenuhnya menangkap pola dalam data. Kurva kepadatan (garis biru) menunjukkan beberapa penyimpangan dari bentuk lonceng normal, dengan beberapa ekor pada kedua sisi dari puncak.
Secara keseluruhan, plot tersebut menunjukkan bahwa residu model tidak sepenuhnya normal, mengindikasikan bahwa model mungkin tidak mencakup semua pola dalam data atau bahwa ada variabel yang tidak termasuk dalam model yang dapat mempengaruhi variabel respons.
Q-Q Plot
# Plot the residuals
qqnorm(model_stepwise_forward_residuals)
# Plot the Q-Q line
qqline(model_stepwise_forward_residuals)
Grafik menunjukkan persebaran data residu observasi terhadap inverse normal. Apabila jarak plot dengan garis semakin kecil, artinya model regresi semakin baik. Dari hasil di atas, terlihat bahwa hanya sebagaian poin-poin yang bersinggungan langsung dengan garis. Hal ini bisa mengindikasikan bahwa distribusi mendekati distribusi normal, namun tetap terdapat outlier.
Melihat scatter plot hasil prediksi
# melakukan predict pada train data (seluruh data di awal tadi)
predicted_values <- predict(model_stepwise_forward)
# Buat dataframe dari nilai asli dan prediksi
comparison_df <- data.frame(Actual = data_reg$PDRB, Predicted = predicted_values)
# Buat scatter plot
ggplot(comparison_df, aes(x = Actual, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") + # Garis y=x untuk referensi
labs(title = "Scatter Plot of Actual vs. Predicted (Stepwise Forward)", x = "Actual Values", y = "Predicted Values") +
theme_minimal()
Grafik menunjukkan persebaran data hasil regesi dibandingkan dengan data hasil observasi. Semakin baik apabila hasil plot membentuk garis diagonal (45 derajat).
Kesimpulan
Dari berbagai hasil post-estimasi di atas, bisa disimpulkan bahwa model masih kurang baik karena terdapat heteroskedasticity yang melanggar asumsi model regresi, meskipun performa model cukup baik.
# Menerapkan stepwise backward regression dengan menggunakan model_all sebegai model yang diiterasi
stepwise_backward <- ols_step_backward_p(model_all, data = data_reg, p_val = 0.05)
# Menampilkan hasil model setelah regresi stepwise backward
stepwise_backward
##
##
## Stepwise Summary
## -------------------------------------------------------------------------------------
## Step Variable AIC SBC SBIC R2 Adj. R2
## -------------------------------------------------------------------------------------
## 0 Full Model 869.811 882.000 808.356 0.79991 0.69987
## 1 jml_tngkrj_inc_kclmngh 867.812 878.782 805.231 0.79990 0.71751
## 2 jml_prod_padi 867.203 876.954 802.490 0.78845 0.71794
## 3 jml_htl_akomodasi 868.024 876.556 800.714 0.76319 0.70087
## 4 jml_tngkrj_inc_bsr 867.135 874.448 798.678 0.75243 0.70291
## 5 jml_unt_ush_nd_kclmngh 866.557 872.652 797.079 0.73793 0.70049
## 6 jml_unt_ush_nd_bsr 868.769 873.644 797.870 0.68984 0.66165
## -------------------------------------------------------------------------------------
##
## Final Model Output
## ------------------
##
## Model Summary
## --------------------------------------------------------------------------------
## R 0.831 RMSE 7249475.227
## R-Squared 0.690 MSE 59721467129518.227
## Adj. R-Squared 0.662 Coef. Var 54.345
## Pred R-Squared 0.151 AIC 868.769
## MAE 5220740.445 SBC 873.644
## --------------------------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
## AIC: Akaike Information Criteria
## SBC: Schwarz Bayesian Criteria
##
## ANOVA
## -----------------------------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## -----------------------------------------------------------------------------------------
## Regression 2922282249482485.000 2 1461141124741242.500 24.466 0.0000
## Residual 1313872276849401.000 22 59721467129518.227
## Total 4236154526331886.000 24
## -----------------------------------------------------------------------------------------
##
## Parameter Estimates
## --------------------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## --------------------------------------------------------------------------------------------------------
## (Intercept) 302669.552 2898460.874 0.104 0.918 -5708370.393 6313709.496
## jml_pnddk 6.644 1.568 0.516 4.238 0.000 3.393 9.894
## jml_trnk_klr 211.571 46.970 0.548 4.504 0.000 114.161 308.981
## --------------------------------------------------------------------------------------------------------
Hasil yang didapatkan sama namun dengan step yang berbeda. Berikut adalah model yang dihasilkan.
model_stepwise_backward <- stepwise_backward$model
summary(model_stepwise_backward)
##
## Call:
## lm(formula = paste(response, "~", paste(c(include, cterms), collapse = " + ")),
## data = l)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21272484 -4300057 -523034 3645258 17199940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 302669.552 2898460.874 0.104 0.917779
## jml_pnddk 6.644 1.568 4.238 0.000337 ***
## jml_trnk_klr 211.571 46.970 4.504 0.000176 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7728000 on 22 degrees of freedom
## Multiple R-squared: 0.6898, Adjusted R-squared: 0.6616
## F-statistic: 24.47 on 2 and 22 DF, p-value: 0.000002555
Interpretasi dan post-estimation akan sama dengan hasil sebelumnya.
Untuk menghilangkan kasus heteroskedasticity yang terdapat pada
regresi stepwise, dapat dilakukan analisis regresi menggunakan
robust yaitu dengan terjelaskannya standar error dan
membuang outlier data yang menyebabkan kasus heteroskedasticity terjadi.
Dalam R, regression robust ini dilakukan menggunakan fungsi
rlm()
dalam paket MASS
.
Berikut adalah inisialisasi variabel dan threshold yang digunakan
library(MASS) # untuk rlm()
##
## Attaching package: 'MASS'
## The following object is masked from 'package:olsrr':
##
## cement
# Inisialisasi
current_vars <- character(0)
remaining_vars <- setdiff(names(data_reg), "PDRB")
model_robust_forward <- NULL
p_value_threshold <- 0.05 # Atau threshold lain yang Anda tentukan
# melakukan iterasi penambahan variabel secara forward selection
repeat {
p_values <- numeric(length(remaining_vars))
names(p_values) <- remaining_vars
# melakukan iterasi running rlm() setiap variabel independen pada data_reg ditambahkan
for (var in remaining_vars) {
formula <- as.formula(paste("PDRB ~", paste(c(current_vars, var), collapse = "+")))
model <- rlm(formula, data = data_reg)
model_summary <- summary(model)
coefficients <- model_summary$coefficients[, "Value"]
std_error <- model_summary$coefficients[, "Std. Error"]
t_stats <- coefficients / std_error
p_values[var] <- 2 * pt(-abs(t_stats[names(t_stats) == var]), df = model_summary$df[2])
}
# Cari variabel dengan p-value terendah
min_p_value <- min(p_values)
min_p_var <- names(which.min(p_values))
# Jika p-value di bawah threshold, tambahkan ke model dan ulangi, jika tidak, berhenti
if (min_p_value < p_value_threshold) {
current_vars <- c(current_vars, min_p_var)
model_robust_forward <- rlm(as.formula(paste("PDRB ~", paste(current_vars, collapse = "+"))), data = data_reg)
remaining_vars <- setdiff(remaining_vars, min_p_var)
} else {
break
}
}
## Warning in rlm.default(x, y, weights, method = method, wt.method = wt.method, :
## 'rlm' failed to converge in 20 steps
## Warning in rlm.default(x, y, weights, method = method, wt.method = wt.method, :
## 'rlm' failed to converge in 20 steps
## Warning in rlm.default(x, y, weights, method = method, wt.method = wt.method, :
## 'rlm' failed to converge in 20 steps
summary(model_robust_forward)
##
## Call: rlm(formula = as.formula(paste("PDRB ~", paste(current_vars,
## collapse = "+"))), data = data_reg)
## Residuals:
## Min 1Q Median 3Q Max
## -21763868 -2212674 -375302 1935965 19637308
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 1746500.0257 1724927.9424 1.0125
## jml_trnk_klr 265.0803 26.2601 10.0944
## jml_pnddk 3.9406 1.1553 3.4108
## jml_unt_ush_nd_bsr 632286.1861 218576.3098 2.8927
##
## Residual standard error: 3281000 on 21 degrees of freedom
Penjelasan setiap baris:
Inisialisasi
Pada bagian ini kita menyiapkan lingkungan untuk proses seleksi:
current_vars
: Sebuah vektor kosong untuk menyimpan
variabel yang telah dipilih.remaining_vars
: Semua variabel lain di
data_reg
kecuali PDRB
, yang merupakan variabel
dependen. Variabel-variabel ini adalah kandidat untuk ditambahkan ke
model.model_robust_forward
: Variabel ini akan menyimpan model
akhir setelah proses seleksi selesai.p_value_threshold
: Ini adalah nilai ambang batas yang
Anda tetapkan untuk p-value. Jika p-value suatu variabel di bawah ambang
ini, variabel tersebut dianggap signifikan.Iterasi Penambahan Variabel
Proses seleksi variabel dilakukan dengan loop
repeat
:
p_values
: Vektor untuk menyimpan p-value untuk setiap
variabel yang tersisa untuk ditambahkan ke model. Ini diinisialisasi
dengan jumlah elemen yang sama dengan jumlah
remaining_vars
.for
menguji setiap variabel yang tersisa dengan
menambahkannya ke model satu per satu.formula
: Ini adalah rumus regresi yang diperbarui
setelah menambahkan variabel baru dari remaining_vars
.model
: Model rlm
yang di-fit dengan
formula
yang diperbarui.model_summary
: Ringkasan dari model yang baru
di-fit.coefficients
dan std_error
: Ini adalah
koefisien estimasi dan kesalahan standar dari model. Diperlukan untuk
menghitung t-statistik.t_stats
: Ini dihitung untuk setiap koefisien untuk
menentukan nilai t-statistik mereka.p_values
: Di sini, p-value dihitung dari t-statistik
dan derajat kebebasan (df) model. P-values ini disimpan kembali ke
vektor p_values
dengan nama variabel yang sesuai.Setelah semua variabel diuji, proses seleksi berlanjut:
min_p_value
dan min_p_var
: Ini adalah
p-value terendah yang ditemukan dan nama variabel yang terkait dengan
p-value tersebut.min_p_value
kurang dari
p_value_threshold
, variabel min_p_var
ditambahkan ke current_vars
, yang berarti variabel tersebut
dianggap signifikan dan termasuk dalam model.model_robust_forward
: Model ini diperbarui dengan
formula yang termasuk semua variabel yang telah dipilih sejauh ini.remaining_vars
: Variabel yang baru saja ditambahkan ke
model dihapus dari daftar variabel yang tersisa.repeat
dihentikan, menandakan bahwa tidak ada lagi
variabel signifikan yang bisa ditambahkan.Model Akhir
Setelah loop selesai, model_robust_forward
berisi model
regresi robust yang termasuk semua variabel yang dipilih melalui proses
forward selection.
summary(model_robust_forward)
: Ini menampilkan
ringkasan statistik dari model akhir, termasuk koefisien estimasi,
kesalahan standar, dan statistik lainnya. Namun, perlu diingat bahwa
p-values harus dihitung dan ditambahkan secara manual, sebagaimana
ditunjukkan dalam contoh kode sebelumnya.Pada akhirnya, Anda mendapatkan model dengan variabel yang dianggap
signifikan berdasarkan p-value yang dihitung secara manual, karena
regresi robust (rlm
) tidak secara otomatis menyediakan
p-values dalam outputnya.
Hasil summary di atas, tidak terdapat nilai p-value untuk memperlihatkan signifikansi setiap variabelnya meskipun sudah diseleksi berdasarkan p-value tersebut. Mari kita tambahkan.
# Misalkan `final_model` adalah model robust regression Anda yang telah di-fit
model_robust_forward_p <- summary(model_robust_forward)
# Mendapatkan koefisien dan standar error dari model final
final_coefficients <- model_robust_forward_p$coefficients[, "Value"]
final_std_error <- model_robust_forward_p$coefficients[, "Std. Error"]
# Menghitung t-statistics
final_t_stats <- final_coefficients / final_std_error
# Menghitung p-values
final_p_values <- 2 * pt(-abs(final_t_stats), df = model_robust_forward_p$df[2])
# Menambahkan p-values ke summary
model_robust_forward_p$coefficients <- cbind(model_robust_forward_p$coefficients, "Pr(>|t|)" = final_p_values)
# Menampilkan ringkasan model yang diperbarui
print(model_robust_forward_p)
##
## Call: rlm(formula = as.formula(paste("PDRB ~", paste(current_vars,
## collapse = "+"))), data = data_reg)
## Residuals:
## Min 1Q Median 3Q Max
## -21763868 -2212674 -375302 1935965 19637308
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1746500.0257 1724927.9424 1.0125 0.3228
## jml_trnk_klr 265.0803 26.2601 10.0944 0.0000
## jml_pnddk 3.9406 1.1553 3.4108 0.0026
## jml_unt_ush_nd_bsr 632286.1861 218576.3098 2.8927 0.0087
##
## Residual standard error: 3281000 on 21 degrees of freedom
# Menghitung Total Sum of Squares (TSS)
y_obs <- data_reg$PDRB # Mengganti dengan nama variabel dependen yang sebenarnya
y_mean <- mean(y_obs)
TSS <- sum((y_obs - y_mean)^2)
# Menghitung Residual Sum of Squares (RSS) dari model final
RSS <- sum(residuals(model_robust_forward)^2)
# Menghitung pseudo R-squared
pseudo_R_squared <- 1 - RSS / TSS
pseudo_R_squared
## [1] 0.7192884
Pada regresi menggunakan robust, uji post estimasi hettest tidak perlu digunakan karena dengan robust suda menghilangkan heteroskedasticity tersebut.
Tes VIF
# Menghitung VIF
vif_robust_forward <- vif(model_robust_forward)
# Melihat hasil VIF
print(vif_robust_forward)
## jml_trnk_klr jml_pnddk jml_unt_ush_nd_bsr
## 1.050585 1.825900 1.768112
Hasilnya tidak terdapat nilai VIF yang terlalu besar menandakan bahwa tidak ada multikolinearitas
Distribusi Residual
# Menghitung residu dari model regresi
model_robust_forward_residuals = model_robust_forward$residuals
# Membuat density plot dengan histogram dari residu
ggplot(data.frame(model_robust_forward_residuals), aes(x=model_robust_forward_residuals)) +
geom_histogram(aes(y=..density..), binwidth = diff(range(model_robust_forward_residuals))/30, fill="skyblue", alpha=0.5) +
geom_density(color="blue", fill="blue", alpha=0.3) +
labs(title="Density and Histogram Plot of Residuals (Robust Stepwise Forward)", x="Residuals", y="Density")
Hasilnya lebih mendekati distribusi normal
Q-Q Plot
# Plot the residuals
qqnorm(model_robust_forward_residuals)
# Plot the Q-Q line
qqline(model_robust_forward_residuals)
Lebih banyak poin-poin yang mendekati garis.
Melihat scatter plot hasil prediksi
# melakukan predict pada train data (seluruh data di awal tadi)
predicted_values <- predict(model_robust_forward)
# Buat dataframe dari nilai asli dan prediksi
comparison_df <- data.frame(Actual = data_reg$PDRB, Predicted = predicted_values)
# Buat scatter plot
ggplot(comparison_df, aes(x = Actual, y = Predicted)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, color = "red") + # Garis y=x untuk referensi
labs(title = "Scatter Plot of Actual vs. Predicted (Robust Stepwise Forward)", x = "Actual Values", y = "Predicted Values") +
theme_minimal()
Kesimpulan
Dari berbagai hasil post-estimasi di atas, bisa disimpulkan bahwa model menjadi lebih baik karena sudah menghilangkan heteroskedasticity dengan performa yang lebih baik juga.
# Inisialisasi
current_vars <- setdiff(names(data_reg), "PDRB")
model_robust_backward <- rlm(as.formula(paste("PDRB ~", paste(current_vars, collapse = "+"))), data = data_reg)
repeat {
p_values <- setNames(numeric(length(current_vars)), current_vars)
for (var in current_vars) {
formula <- as.formula(paste("PDRB ~", paste(current_vars, collapse = "+"))) # Gunakan model dengan semua variabel saat ini
model <- rlm(formula, data = data_reg)
model_summary <- summary(model)
coefficients <- model_summary$coefficients[, "Value"]
std_error <- model_summary$coefficients[, "Std. Error"]
t_stats <- coefficients / std_error
# Hitung p-value untuk setiap variabel tanpa menghapusnya dari model
p_values[var] <- 2 * pt(-abs(t_stats[names(t_stats) == var]), df = model_summary$df[2])
}
# Cari variabel dengan p-value tertinggi
max_p_value <- max(p_values)
max_p_var <- names(which.max(p_values))
# Jika p-value di atas threshold, hapus dari model dan ulangi, jika tidak, berhenti
if (max_p_value > p_value_threshold && length(current_vars) > 1) { # Pastikan setidaknya satu variabel tetap
current_vars <- setdiff(current_vars, max_p_var)
model_robust_backward <- rlm(as.formula(paste("PDRB ~", paste(current_vars, collapse = "+"))), data = data_reg)
} else {
break
}
}
summary(model_robust_backward)
##
## Call: rlm(formula = as.formula(paste("PDRB ~", paste(current_vars,
## collapse = "+"))), data = data_reg)
## Residuals:
## Min 1Q Median 3Q Max
## -21763868 -2212674 -375302 1935965 19637308
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 1746500.0257 1724927.9424 1.0125
## jml_pnddk 3.9406 1.1553 3.4108
## jml_trnk_klr 265.0803 26.2601 10.0944
## jml_unt_ush_nd_bsr 632286.1861 218576.3098 2.8927
##
## Residual standard error: 3281000 on 21 degrees of freedom
Hasilnya sama dengan forward selection di awal. Tidak perlu melanjutkan proses post-estimation.
Kita review ulang seluruh hasil post-estimation dari seluruh model yang kita gunakan
All Variables
summary_model_all <- summary(model_all)
# Mengambil nilai R-squared
r_squared <- summary_model_all$r.squared
# Mengambil nilai adjusted R-squared
adj_r_squared <- summary_model_all$adj.r.squared
# Menampilkan nilai R-squared dan adjusted R-squared
print(paste("R-squared:", r_squared))
## [1] "R-squared: 0.79991284031587"
print(paste("Adjusted R-squared:", adj_r_squared))
## [1] "Adjusted R-squared: 0.699869260473805"
# fungsi mengambil p-value
overall_p <- function(my_model) {
f <- summary(my_model)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
# hasil p-value model
print(paste("Prob > F:", overall_p(model_all)))
## [1] "Prob > F: 0.000235967824418562"
bptest_all
##
## studentized Breusch-Pagan test
##
## data: model_all
## BP = 8.9261, df = 8, p-value = 0.3486
Tidak terdapat heteroskedasticity.
vif_model_all
## jml_pnddk jml_prod_padi jml_trnk_klr
## 4.401629 1.778473 1.514696
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_kclmngh
## 2.062389 3.782425 13.584754
## jml_unt_ush_nd_bsr jml_htl_akomodasi
## 15.068421 1.506562
Terdapat multikolinearitas
Pengurangan 1 variabel
summary_model_minus_one <- summary(model_minus_one)
# Mengambil nilai R-squared
r_squared <- summary_model_minus_one$r.squared
# Mengambil nilai adjusted R-squared
adj_r_squared <- summary_model_minus_one$adj.r.squared
# Menampilkan nilai R-squared dan adjusted R-squared
print(paste("R-squared:", r_squared))
## [1] "R-squared: 0.769475391511165"
print(paste("Adjusted R-squared:", adj_r_squared))
## [1] "Adjusted R-squared: 0.674553493898115"
# hasil p-value model
print(paste("Prob > F:", overall_p(model_minus_one)))
## [1] "Prob > F: 0.000216710731102793"
bptest_minus_one <- bptest(model_minus_one)
bptest_minus_one
##
## studentized Breusch-Pagan test
##
## data: model_minus_one
## BP = 11.504, df = 7, p-value = 0.1181
Tidak terdapat heteroskedasticity.
vif_model_minus_one
## jml_pnddk jml_prod_padi jml_trnk_klr
## 4.388750 1.773715 1.120449
## jml_tngkrj_inc_kclmngh jml_tngkrj_inc_bsr jml_unt_ush_nd_bsr
## 2.041051 3.358771 2.460902
## jml_htl_akomodasi
## 1.462369
Tidak terdapat multikolinearitas
Stepwise Forward = Stepwise Backward
# Mengambil nilai R-squared
r_squared <- summary(model_stepwise_forward)$r.squared
# Mengambil nilai adjusted R-squared
adj_r_squared <- summary(model_stepwise_forward)$adj.r.squared
# Menampilkan nilai R-squared dan adjusted R-squared
print(paste("R-squared:", r_squared))
## [1] "R-squared: 0.689843165851863"
print(paste("Adjusted R-squared:", adj_r_squared))
## [1] "Adjusted R-squared: 0.661647090020214"
# hasil p-value model
print(paste("Prob > F:", overall_p(model_stepwise_forward)))
## [1] "Prob > F: 0.00000255502353981603"
bptest_stepwise_forward
##
## studentized Breusch-Pagan test
##
## data: model_stepwise_forward
## BP = 6.2264, df = 2, p-value = 0.04446
Terdapat heteroskedasticity
vif_model_stepwise_forward
## jml_trnk_klr jml_pnddk
## 1.050144 1.050144
Tidak terdapat multikolinearitas
Robust Forward = Robust Backward
pseudo_R_squared
## [1] 0.7192884
vif_robust_forward
## jml_trnk_klr jml_pnddk jml_unt_ush_nd_bsr
## 1.050585 1.825900 1.768112
Dari hasil, kita bisa memilih model hasil robust (baik forward maupun backward).
summary(model_robust_forward)
##
## Call: rlm(formula = as.formula(paste("PDRB ~", paste(current_vars,
## collapse = "+"))), data = data_reg)
## Residuals:
## Min 1Q Median 3Q Max
## -21763868 -2212674 -375302 1935965 19637308
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 1746500.0257 1724927.9424 1.0125
## jml_trnk_klr 265.0803 26.2601 10.0944
## jml_pnddk 3.9406 1.1553 3.4108
## jml_unt_ush_nd_bsr 632286.1861 218576.3098 2.8927
##
## Residual standard error: 3281000 on 21 degrees of freedom
Residuals:
Residuals merupakan perbedaan antara nilai observasi sebenarnya dan nilai yang diprediksi oleh model. Pada output Anda, kita memiliki:
Distribusi residual yang lebar seperti ini bisa mengindikasikan model yang kurang fit dengan data atau adanya outlier yang signifikan.
Coefficients:
Ini adalah estimasi koefisien untuk variabel independen dalam model regresi.
Intercept: Estimasi untuk intercept (konstanta) adalah 1746500.0257 dengan standard error 1742497.9424 dan t-value 1.0125. Nilai t ini menunjukkan bahwa intercept tidak secara statistik signifikan (asumsi tingkat signifikansi standar 0.05), karena biasanya kita mencari nilai t yang lebih besar dari 2 atau kurang dari -2 untuk signifikansi.
jml_trnk_klr: Koefisien untuk variabel
jml_trnk_klr
adalah 265.0803 dengan standard error 26.2601
dan t-value 10.0944. Ini menunjukkan bahwa variabel ini signifikan
secara statistik dalam model.
jml_pnddk: Koefisien untuk
jml_pnddk
adalah 3.9406 dengan standard error 1.1553 dan
t-value 3.4108, yang juga menandakan signifikansi statistik.
jml_unt_ush_nd_bsr: Koefisien untuk
jml_unt_ush_nd_bsr
adalah 632286.1861 dengan standard error
218576.3098 dan t-value 2.8927, menunjukkan signifikansi
statistik.
Residual standard error:
Ini adalah perkiraan deviasi standar dari residuals, yang adalah 3281000 pada 21 derajat kebebasan (jumlah observasi dikurangi jumlah estimasi koefisien).
Formula Regresi:
Dengan menggunakan markdown, formula regresi bisa ditulis sebagai:
\[ Y = 1746500.0257 + 265.0803 * jml_trnk_klr + 3.9406 * jml_pnddk + 632286.1861 * jml_unt_ush_nd_bsr \]
Dimana: - Y
adalah variabel dependen yang ingin
diprediksi (PDRB). - jml_trnk_klr
, jml_pnddk
,
dan jml_unt_ush_nd_bsr
adalah variabel independen dalam
model.
Di sini, Y
merupakan variabel dependen, yang artinya
nilai Y
diprediksi berdasarkan kombinasi dari
variabel-variabel independen dalam model. Setiap koefisien memberikan
informasi tentang seberapa banyak kita mengharapkan Y
untuk
berubah, dengan asumsi bahwa variabel independen lainnya tetap
konstan.
Intercept (1746500.0257):
Intercept, atau konstanta, mewakili nilai rata-rata yang diharapkan
untuk Y
ketika semua variabel independen
(jml_trnk_klr
, jml_pnddk
,
jml_unt_ush_nd_bsr
) bernilai nol. Meskipun dalam banyak
konteks praktis, interpretasi ini tidak selalu masuk akal jika nilai nol
dari variabel independen tidak mungkin atau tidak masuk akal.
Koefisien jml_trnk_klr
(265.0803):
Koefisien ini mengindikasikan bahwa dengan setiap peningkatan satu
unit pada jumlah ternak keluar
, nilai PDRB diperkirakan
akan meningkat sebesar 265.0803, asumsi variabel lainnya tetap.
Koefisien jml_pnddk
(3.9406):
Koefisien ini menunjukkan bahwa dengan setiap peningkatan satu unit
pada jumlah penduduk
, nilai PDRB diperkirakan akan
meningkat sebesar 3.9406, asumsi variabel lainnya konstan.
Koefisien jml_unt_ush_nd_bsr
(632286.1861):
Koefisien ini menyatakan bahwa untuk setiap peningkatan satu unit
dalam jumlah unit usaha besar
, nilai PDRB diperkirakan akan
meningkat sebesar 632286.1861, dengan asumsi bahwa variabel independen
lainnya tidak berubah.
Pentingnya Koefisien:
Tingkat signifikansi statistik dari koefisien diukur oleh nilai t dan standar error mereka: - Koefisien dengan nilai t tinggi dan standar error yang rendah dianggap sebagai prediktor yang signifikan dan kuat dalam model. - Nilai t di atas 2 atau di bawah -2 (secara tradisional pada tingkat kepercayaan 95%) sering digunakan sebagai aturan praktis untuk menentukan signifikansi statistik. Selain itu, hasil stepwise menunjukkan bahwa seluruh variabel signifikan pada tingkat kepercayaan 95% juga/
Dalam konteks ini, semua variabel independen kecuali intercept
tampaknya memiliki kontribusi signifikan terhadap nilai Y
berdasarkan nilai t mereka yang relatif tinggi. Namun, karena standard
error yang besar pada intercept, kita mungkin memiliki beberapa keraguan
tentang kestabilan estimasi intercept tersebut.
Jika jumlah ternak yang keluar dari Provinsi Jawa Barat sebanyak 100.000, jumlah penduduknya 10.000.000 jiwa jumlah penduduk dan jumlah unit usaha besar 100
# Membuat data frame baru dengan nilai variabel independen yang ingin diprediksi
new_data <- data.frame(jml_trnk_klr = 100000, jml_pnddk = 10000000, jml_unt_ush_nd_bsr = 100)
# Menggunakan fungsi predict() untuk mendapatkan estimasi dari model regresi robust
predicted_values <- predict(model_robust_forward, newdata = new_data)
# Menampilkan nilai prediksi
print(predicted_values)
## 1
## 130889337
Berikut adalah hasilnya:
\[ PDRB = 1746500.0257 + 265.0803 * 100000 + 3.9406 * 10000000 + 632286.1861 * 100 = 130889337 \]
PDRB yang diprediksi adalah senilai Rp 130.889.337