Modul 3: Analisis Peramalan Kausal (Regresi Berganda)

Objectives Praktikum

  1. Praktikan mampu memahami konsep-konsep dan prinsip dasar pada analisis peramalan kausal (regresi berganda)
  2. Praktikan mampu menjalankan prosedur analisis peramalan kausal dengan R dan RStudio.
  3. Praktikan mampu melakukan interpretasi hasil analisis peramalan kausal (regresi berganda) dalam konteks PWK.

Review Teori

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.

    • Setiap nilai variabel prediktor X ada suatu distribusi probabilitas dari nilai-nilai kriteria Y. Untuk setiap distribusi Y, satu atau lebih diambil sampel secara random
    • Variansi distribusi Y semuanya sama satu dengan yang lain (homosedasticity).
    • Rata-rata distribusi Y jatuh pada garis regresi \(\mu_Y = \alpha+\beta X\), di mana \(\mu_Y\) adalah rata-rata distribusi Y untuk nilai X tertentu, \(\beta\) adalah kelelerangan garis regresi, \(\alpha\) perpotongan garis regresi dengan sumbu Y.
  • 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.

  • Membuat daftar panjang variabel prediktor. Variabel prediktor atau independen dapat diperoleh dari teori dan logika hasil studi literatur, serta arahan pakar.
  • Membuat daftar pendek variabel prediktor. Persamaan regresi hendaknya mempunyai variabel prediktor sesedikit mungkin sehingga variabel prediktor yang tidak terlalu berperan harus dihapus.

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.

Multikolinearity

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.

  • Bila multicollinearity terjadi dalam model regresi berganda, solusi least square tidak mungkin dihasilkan. Hal ini disebabkan karena solusi least square dipengaruhi oleh error yang berlebihan.
  • Kestabilan model regresi berganda dipengaruhi oleh multicollinearity. Standard error koefisien dipengaruhi oleh korelasi variabel bebas, semakin besar korelasi, dan koefisien semakin tidak stabil.
  • Dapat mempengaruhi hasil peramalan, yaitu hasil peramalan menjadi over estimate.

Multikolinearitas dapat dihindari dengan beberapa alternatif cara yang diantaranya adalah sebagai berikut.

  • Menyisihkan atau membuang data ekstrim dalam suatu variabel.
  • Memperbesar ukuran sampel.
  • Memasukkan persamaan tambahan atau data baru ke dalam model.
  • Menghubungkan data cross section dan data time series.
  • Membuang salah satu atau beberapa variabel yang berkorelasi tinggi.
  • Mentransformasikan variabel.
  • Melakukan teknik pemodelan yang lebih kompleks terlebih dahulu seperti analisis faktor sehingga variabel-variabel yang menyebabkan multikolinearitas akan difiltrasi terlebih dahulu.

Dummy Variable

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.

Handson Praktikum

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:

  • Nama Kabupaten/Kota
  • Jumlah penduduk
  • Jumlah produksi padi
  • Jumlah ternak yang keluar
  • Jumlah tenaga kerja industri kecil menengah
  • Jumlah tenaga kerja industri besar
  • Jumlah unit usaha industri kecil menengah
  • Jumlah unit usaha industri besar
  • Jumlah hotel dan akomodasi
  • PDRB

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

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

Rumusan Masalah

Masalah yang akan diteliti adalah sebagai berikut.

  • Model regresi mana yang tepat untuk menjelaskan PDRB di Jawa Barat? Jelaskan alasannya!
  • Variabel-variabel apa saja yang paling baik untuk menjelaskan PDRB di Jawa Barat?
  • Bagaimana interpretasi dari model regresi yang terbentuk? (pola dan arah hubungan antara PDRB dengan variabel-variabel yang memengaruhinya)
  • Bagaimana validitas dan signifikansi dari model yang terbentuk?
  • Perkirakan PDRB suatu kota di Jawa Barat dengan jumlah ternak yang keluar dari Provinsi Jawa Barat sebanyak 100 ribu, jumlah penduduknya 10 juta jiwa jumlah penduduk dan jumlah unit usaha besar 100 (Tingkat kepercayaan 95% dan Standar Error = 1346308)

Descriptive Statistics

# 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.

Missing values

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

Correlation

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
Matriks Korelasi Variabel Numerik
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)

Regresi Seluruh Variabel

1. Regresi 8 Variabel

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:

    • Estimasi koefisien untuk setiap variabel independen dan intercept diberikan di kolom 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.
    • Variabel 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)\)

Post-Estimation Model Regresi (_All)

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.

2. Regresi 7 variabel

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

Regresi dengan Stepwise

1. Regresi dengan Stepwise Forward atau Variables Entered

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.

  • AIC (Akaike Information Criteria), SBC (Schwarz Bayesian Criterion), dan SBIC adalah metrik untuk membandingkan model; semakin kecil nilainya, semakin baik modelnya.
  • R2 adalah koefisien determinasi yang menunjukkan seberapa baik data sesuai dengan model regresi. Adj. R2 adalah versi yang disesuaikan dari R2 yang mempertimbangkan jumlah prediktor dalam model.

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.

  • R adalah akar kuadrat dari R-squared.
  • R-Squared adalah proporsi varians dalam variabel dependen yang dapat dijelaskan oleh variabel independen dalam model. Nilainya 0.831 menunjukkan bahwa model menjelaskan 83.1% varians.
  • Adj. R-Squared menyesuaikan R-squared berdasarkan jumlah prediktor dan jumlah sampel, nilainya 0.662.
  • RMSE (Root Mean Square Error) menunjukkan ukuran rata-rata kesalahan yang dihasilkan oleh model, yaitu 59721467219518.180.
  • MAE (Mean Absolute Error) adalah rata-rata absolut dari kesalahan dalam model, yaitu 5220740.445.

ANOVA

Tabel ini menunjukkan Analisis Varians yang membandingkan model dengan model tanpa prediktor. Hal ini dilakukan untuk menguji apakah prediktor secara keseluruhan signifikan.

  • DF adalah derajat kebebasan.
  • Sum of Squares adalah variasi yang dijelaskan oleh model dan variasi residual (yang tidak dijelaskan oleh model).
  • Mean Square adalah Sum of Squares dibagi dengan derajat kebebasan.
  • F adalah statistik F, yang digunakan untuk menentukan apakah secara keseluruhan variabel independen berdampak signifikan pada variabel dependen.
  • Sig. adalah nilai p dari F-statistik, di sini menunjukkan bahwa model kita memiliki prediktor yang secara signifikan mempengaruhi variabel respons (dengan nilai p sangat kecil, hampir mendekati 0).

Parameter Estimates

Ini adalah tabel yang memberikan estimasi untuk parameter (koefisien) dari model regresi.

  • model menunjukkan variabel independen.
  • Beta adalah koefisien estimasi.
  • Std. Error menunjukkan standar kesalahan dari koefisien.
  • Std. Beta adalah koefisien beta yang telah distandarisasi.
  • t adalah nilai t, digunakan untuk menentukan apakah masing-masing koefisien berdampak signifikan pada variabel dependen.
  • Sig. adalah nilai p dari tes t, di mana nilai kurang dari 0.05 umumnya menunjukkan bahwa koefisien secara statistik signifikan. Di sini, ‘jml_trnkk_klr’ dan ‘_jml_pndkk’ keduanya signifikan.

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

  • Grafik ini menunjukkan peningkatan dalam R-squared dari model dasar tanpa prediktor ke model akhir dengan dua prediktor. Nilai R-squared meningkat dari 0 ke sekitar 0.69, menunjukkan bahwa model final dengan kedua prediktor menjelaskan sekitar 69% dari variabilitas dalam variabel dependen, yang merupakan peningkatan yang signifikan dari model dasar.

Adjusted R-Square Plot

  • Plot ini mirip dengan plot R-squared tetapi menunjukkan Adjusted R-squared, yang juga meningkat dari 0 menjadi sekitar 0.662. Ini menyesuaikan untuk jumlah prediktor dalam model dan mencegah overfitting. Peningkatan ini menunjukkan bahwa setelah menyesuaikan untuk jumlah variabel, model masih menjelaskan sebagian besar variabilitas dalam data.

Akaike Information Criteria (AIC) Plot

  • Plot AIC menggambarkan penurunan nilai AIC saat lebih banyak variabel ditambahkan ke model. Dimulai dengan nilai tinggi pada model dasar dan menurun secara signifikan setelah menambahkan kedua prediktor, yang menunjukkan bahwa model final lebih disukai dibandingkan model hanya dengan satu prediktor atau tanpa prediktor.

Root Mean Squared Error (RMSE) Plot

  • RMSE menunjukkan ukuran kesalahan model. Plot RMSE ini menunjukkan penurunan yang sangat signifikan dalam kesalahan ketika prediktor pertama ditambahkan, dan penurunan lebih lanjut namun lebih kecil ketika prediktor kedua ditambahkan. Ini menunjukkan bahwa model menjadi lebih akurat dalam memprediksi variabel dependen dengan penambahan setiap prediktor.

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.

Post-Estimation Regresi Stepwise Forward

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.

2. Regresi dengan Stepwise Backward

# 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.

Regresi Robust dengan Stepwise

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.

1. Robust dengan Stepwise Forward

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.
  • Iterasi 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.
  • Jika 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.
  • Jika tidak ada variabel yang memiliki p-value di bawah ambang batas, loop 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
Post-Estimation Regresi Robust Stepwise Forward

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.

2. Robust dengan Stepwise Backward

# 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.

Pemilihan Model

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).

Interpretasi Model

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:

  • Min: Nilai residual terkecil adalah -21763868, ini menunjukkan bahwa nilai prediksi paling rendah meleset dari nilai sebenarnya sebanyak itu.
  • 1Q (Kuartil Pertama): 25% dari residuals adalah lebih kecil dari atau sama dengan -2212674.
  • Median: Nilai tengah dari residuals adalah -375302, ini adalah perbedaan median antara nilai prediksi dan nilai observasi.
  • 3Q (Kuartil Ketiga): 75% dari residuals adalah lebih kecil dari atau sama dengan 1935965.
  • Max: Nilai residual terbesar adalah 19637308, ini menunjukkan bahwa nilai prediksi paling tinggi meleset dari nilai sebenarnya sebanyak itu.

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