Analisis Regresi Linier Berganda dengan R

Ringkasan Data

Sebelum masuk ke tahap pemodelan, sangat disarankan untuk melakukan eksplorasi data untuk lebih memahami kondisi data. Terdapat beberapa hal yang akan kita lakukan dan dimulai dengan mengecek ringkasan data. Hasil ringkasan data, memberikan sedikit gambaran bagaimana karakteristik masing-masing peubah. Misalkan untuk peubah respon kita yaitu ipm_2018, nilai terkecil pada data adalah 61,00 dan nilai terbesar adalah 86,11 denga nilai rataan adalah 71,23.

Kita juga perlu memastikan tidak ada data hilang. Jika ternyata terdapat data hilang, maka perlu dilakukan penanganan, entah dengan menghapus baris, kolom atau melakukan imputasi. Pada dataset yang kita gunakan kebetulan tidak terdapat data hilang sehingga dapat dilanjutkan dengan pengecekan lainnya.

# Melihat summary data
summary(data)

##     ipm_2018       PR_NO_LIS           PR_SAMPAH          PR_TINJA     
##  Min.   :61.00   Min.   :0.0000000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:67.91   1st Qu.:0.0000000   1st Qu.:0.00000   1st Qu.:0.1522  
##  Median :71.23   Median :0.0001314   Median :0.04982   Median :0.3169  
##  Mean   :71.98   Mean   :0.0015136   Mean   :0.07444   Mean   :0.3485  
##  3rd Qu.:74.92   3rd Qu.:0.0010234   3rd Qu.:0.10515   3rd Qu.:0.4942  
##  Max.   :86.11   Max.   :0.0246482   Max.   :0.35890   Max.   :0.9471  
##  PR_MKM_SUNGAI     PR_SUNGAI_LMBH      PR_KUMUH          PRA_1000     
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.00000   Min.   :0.1344  
##  1st Qu.:0.07685   1st Qu.:0.1428   1st Qu.:0.02484   1st Qu.:0.3792  
##  Median :0.12648   Median :0.2174   Median :0.07576   Median :0.4748  
##  Mean   :0.17818   Mean   :0.2682   Mean   :0.15779   Mean   :0.5107  
##  3rd Qu.:0.23558   3rd Qu.:0.3617   3rd Qu.:0.19812   3rd Qu.:0.6020  
##  Max.   :0.76471   Max.   :0.9649   Max.   :0.84091   Max.   :1.2772  
##     SD_1000          SM_1000       RS_PKS_PDK_1000  LIN_BID_POS_P_1000
##  Min.   :0.2093   Min.   :0.1627   Min.   :0.1487   Min.   :0.05515   
##  1st Qu.:0.4782   1st Qu.:0.2237   1st Qu.:0.2198   1st Qu.:0.26374   
##  Median :0.6318   Median :0.2760   Median :0.2641   Median :0.43536   
##  Mean   :0.6159   Mean   :0.3019   Mean   :0.3080   Mean   :0.43437   
##  3rd Qu.:0.7508   3rd Qu.:0.3433   3rd Qu.:0.3486   3rd Qu.:0.58097   
##  Max.   :1.1616   Max.   :0.7648   Max.   :0.8780   Max.   :0.92269   
##   APT_OBT_1000      DOK_DRG_1000        BID_1000       GZ_BURUK_1000    
##  Min.   :0.04144   Min.   :0.06402   Min.   :0.04001   Min.   :0.00000  
##  1st Qu.:0.13320   1st Qu.:0.13407   1st Qu.:0.31515   1st Qu.:0.03831  
##  Median :0.18011   Median :0.17146   Median :0.42429   Median :0.07618  
##  Mean   :0.21075   Mean   :0.23894   Mean   :0.42723   Mean   :0.12594  
##  3rd Qu.:0.24919   3rd Qu.:0.28077   3rd Qu.:0.52864   3rd Qu.:0.15679  
##  Max.   :1.09135   Max.   :1.07837   Max.   :1.90602   Max.   :0.93573

# mengecek missing value
# Melihat nilai NA
cek.na <- colSums(is.na(data))

# merubah menjadi dataframe agar lebih rapi
cek.na.df <- data.frame(Variabel=names(cek.na), "Jumlah NA"=cek.na, row.names = NULL)
print(cek.na.df)

##              Variabel Jumlah.NA
## 1            ipm_2018         0
## 2           PR_NO_LIS         0
## 3           PR_SAMPAH         0
## 4            PR_TINJA         0
## 5       PR_MKM_SUNGAI         0
## 6      PR_SUNGAI_LMBH         0
## 7            PR_KUMUH         0
## 8            PRA_1000         0
## 9             SD_1000         0
## 10            SM_1000         0
## 11    RS_PKS_PDK_1000         0
## 12 LIN_BID_POS_P_1000         0
## 13       APT_OBT_1000         0
## 14       DOK_DRG_1000         0
## 15           BID_1000         0
## 16      GZ_BURUK_1000         0

Boxplot Peubah Numerik

Semua variabel pada daset bertipe numerik, sehingga kita dapat membuat boxplot untuk masing-masing variabel. melalui boxplot kita dapat memperoleh informasi bagaimana data menyebar dan juga dapat digunakan unutk mengidentifikasi keberadaan outlier.

Dengan memanfaatkan fungsi lapply kita dapat membuat sekaligus semua boxplot dengan beberapa baris kode saja. Hasil visualisasi boxplot di bawah ini menunjukkan bahwa masing-masing variabel memiliki pola sebaran yang beragam. Sebagian besar terdapat titik-titik data yang diindikasikan sebagai outlier. Tentunya hal ini perlu dilihat dengan lebih seksama. Jika memang kondisi tersebut sesuatu yang wajar berarti tidak ada masalah. Namun jika memang dianggap ada yang tidak sesuai maka perlu dilakukan pengecekan dan penanganan kondisi tersebut.

Di sini, kita akan menganggap semua kondisi ini sesuai dengan fakta lapangan dan tidak akan melakukan penanganan apapun.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(cowplot)

plots <- lapply(names(data), function(var_x){
  p <- 
    ggplot(data) +
    aes(x=!!sym(var_x)) +
    geom_boxplot(lwd=0.2, color="black", fill = "darkgreen", alpha=0.8)+
    theme(axis.text = element_blank())
})

plot_grid(plotlist = plots)

Korelasi Antar Variabel

Melihat koerelasi antar variabel merupakan bagian penting. Hal ini dapat menjadi indikasi awal apakah terdapat kondisi multikolinieritas pada data yang kita miliki. Tentunya lebih jauh kita dapat melihat nilai VIF untuk lebih meyakinkan.

Dari plot korelasi berikut, terlihat beberapa kolom memiliki nilai korelasi yang tinggi, sehingga mungkin akan terdapat kondisi multikolinieritas. Adapun untuk korelasi antara peubah penjelas dengan peubah respon, terlihat beberapa peubah memiliki korelasi yang cukup tinggi seperti DOK_DRG_1000, SD_1000, PR_TINJA dan sebagainya.

Selain itu, korelasi ini dapat menjadi informasi awal untuk melihat bagaimana hubungan antara setiap variabe penjelas dengan variabel respon. Sebagai contoh, cukup mengejutkan melihat korelasi beberapa data yang tampak agak “aneh”. Misal, korelasi SD_1000 atau SM_1000 dengan nilai IPM memiliki korelasi yang negatif. Semakin banyak SD_1000 maka nilai IPM di Kabupaten/Kota tersebut semakin rendah. Sebaliknya Korelasi PR_KUMUH serta PR_SUNGAI_LMBH terhadap nilai IPM bernilai positif. Dimana artinya, semakin tinggi proporsi pemukiman kumuh, atau proporsi sungai yang tercemar limbah maka nilai IPM di kabupaten tersebut cenderung lebih tinggi juga.

Tentunya banyak hal yang mungkin menyebabkan kondisi tersebut terjadi, dan tentunya perlu dikaji secara lebih komprehensif. Namun untuk tulisan singkkat ini, kita tidak akan melihat jauh di luar teknis analisis regresi.

library(corrplot)

## Warning: package 'corrplot' was built under R version 4.2.3

## corrplot 0.92 loaded

# menghitung korelasi antar kolom
corr_matrix <- round(cor(data), 2)

# membuat plot korelasi
corrplot(corr_matrix, 
         type="lower",
         method = "color", 
         tl.cex = 0.6, 
         tl.col = "black",
         addCoef.col = "#2F2F2F",
         addCoefasPercent = FALSE,
         number.cex = 0.6,
         diag = FALSE)

Koefisien Regresi

Nilai ini menunjukkan nilai koefisien regresi ($\hat{\beta_i}$) untuk masing-masing variabel. Misalkan, untuk variabel PR_NO_LIS memiliki nilai koefisien -54,167 dan untuk SD_1000 bernilai 8,48. Tanda positif menunjukkan hubungan yang positif antara variabel penjelas dan responnya, sebaliknya tanda negatif menunjukkan hubungan negatif antara keduanya.

Statistik uji t

Komponen ini menunjukkan variabel penjelas mana saja yang signifikan. Jika nilai p-value lebih kecil dari nilai $\alpha$ tertentu, misal 0,05 maka peubah tersebut siginifikan pada tingkat keyakinan $100- $\alpha$. Informasi ini juga dapat dilihat dari kode signifikansi yang ditambahkan pada kolom terakhir. Tanda ‘***’ menyatakan variabel sangat signifikan bahkan untuk nilai $\alpha$ yang sangat kecil. Adapun untuk nilai yang standar dipakai yaitu $\alpha=0,05$ ditandai dengan simbol ‘.’ serta $\alpha=0,01$ ditandai dengan simbol ‘*’. Hasil di atas menunjukkan variabel mana saja yang signifikan atau tidak, dan terlihat bahwa banyak variabel yang tidak signifikan. Hasil ini tentunya perlu diuji lagi validitasnya dengan pengecekan asumsi serta kondisi yang tidak boleh dilanggar pada model regresi linier.

Statistik uji F

Jika statistik uji $t$ menunjukkan secara parsial signifikansi masing-masing variabel, maka statistik uji $F$ mengukur signifikansi model secara keseluruhan. Pada contoh ini, p-value $≈$ 0 ($p-value: < 2.2e-16$) menunjukkan model signifikan (bahkan dengan tingkat kesalahan yang sangat kecil, misal $alpha=0.00001$ atau lebih kecil lagi) . Dengan kata lain, terdapat setidaknya 1 variabel penjelas yang signifikan (atau setidaknya terdapat satu $i$ untuk $\text{i=1,2…,15}$ di mana nilai $\beta_i\neq0$. Hal ini tentunya juga sudah sesuai hasil pada uji parsial.

Keofisien Determinasi ($R^2$)

Nilai $R^2$ menunjukkan besarnya keragaman pada variabel respon yang dapat dijelaskan oleh variabel penjelas. Pada contoh ini nilai $R^2=0.8474$ artinya variabel-variabel penjelas pada model mampu menjelaskan sekitar 84,74 persen keragaman pada nilai IPM.

Adjusted $R^2$

Nilai adjusted $R^2$ adalah nilai $R^2$ yang diberikan penalti berupa banyaknya variabel penjelas yang digunakan. Ukuran ini biasanya digunakan saat membandingkan kinerja antar model dengan tetap mempertimbangkan kompleksitas model. Model yang lebih sederhana akan di penalti lebih kecil dibandingkan model yang kompleks. Sebagai contoh, Model A dengan 15 variabel memiliki nilai $R^2=0.84$ dan model B dengan hanya 5 variabel memiliki niai $R^2=0.82$. Jika dihitung nilai adjusted $R^2$, bisa jadi Model B akan memberikan nilai yang lebih tinggi dibandingkan Model A. Hal ini dapat dimaknai, dengan pengurangan kompleksitas model yang besar (dari 15 variabel menjadi hanya 5 variabel) model tetap memiliki kinerja yang baik dan hanya sedikit saja mengurangi $R^2$ (dari 0,84 menjadi 0,82). Jika terjadi situasi seperti ini mungkin kita akan lebih cenderung memilih model B dibandingkan A.

Linieritas

Premis utama dari analisis regresi linier berganda adalah adanya hubungan linear antara variabel respon dan variabel penjelas. Linearitas ini dapat diperiksa secara visual menggunakan scatterplot, yang seharusnya menunjukkan hubungan dengan kecenderungan garis lurus.

library(ggplot2)
library(cowplot)

var.respon <- "ipm_2018"
var.penjelas <- setdiff(names(data), var.respon)

plots <- lapply( var.penjelas, function(var_x){
  p <- 
    ggplot(data, aes_string(x = var_x, y = var.respon)) +
    aes(x=!!sym(var_x)) +
    geom_point(color = "blue") +
    theme_minimal() +
    theme(
      axis.text = element_blank(),
      axis.title = element_text(size = 8)
    )
})

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot_grid(plotlist = plots, nrow = 5, ncol = 3)

Berdasarkan scatter plot dapat dilihat pola hubungan antara variabel respon dengan variabel-variabel penjelasnya. Pada beberapa variabel penjelas, terlihat terdapat hubungan linier yang cukup kuat. Namun terdapat beberapa variabel yang sepertinya tidak berhubungan secara linier. Variabel-variabel ini yang kemungkinan tidak memiliki pengaruh signifikan.

Keputusan yang diambil tergantung dengan kepentingan peneliti, variabel yang tidak berhubungan linier seperti ini dapat dikeluarkan dari model ataupun dilakukan penanganan tertentu. Jika variabel-variabel tersebut misalkan tetap ingin dimasukkan, salah satu teknik yang dapat dicoba adalah transformasi data. Transformasi ini dapat berupa log natural, akar kuadrat, transformasi box-cox dan sebagainya.

Normalitas Residual

Analisis regresi linier berganda mengasumsikan bahwa residual (perbedaan antara nilai yang diamati dan nilai yang diprediksi) terdistribusi secara normal. Asumsi ini dapat dinilai dengan memeriksa histogram atau Q-Q plot dari residual, atau melalui uji statistik seperti uji Kolmogorov-Smirnov, Shapiro-Wilk dan sebagainya.

Pada hasil qq plot maupun uji formal menggunakan uji K-S (tidak tolak $H_0$), residual model memenuhi asumsi normalitas yaitu menyebar $\text{Normal}(0, \sigma^2)$.

# Mendapatkan residual dari model
model.res <- residuals(model_reg)

mean.res = round(mean(model.res), 4)
sd.res = round(sd(model.res), 4)

cat("Nilai tengah residual:", mean.res, "\n")

## Nilai tengah residual: 0

cat("Simpangan baku residual:", sd.res, "\n")

## Simpangan baku residual: 2.1323

# Membuat Q-Q plot dari residual
qqnorm(model.res)
qqline(model.res)

# Melakukan uji Kolmogorov-Smirnov pada residual
ks_result <- ks.test(model.res, "pnorm", mean = mean.res, sd = sd.res)
ks_result

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  model.res
## D = 0.049367, p-value = 0.9339
## alternative hypothesis: two-sided

Homoskedastisitas

Ragam dari residual harus konsisten di seluruh tingkat variabel penjelas. Scatterplot dari residual versus nilai prediksi tidak boleh menampilkan pola yang terlihat, seperti distribusi berbentuk kerucut ataupun corong, yang artinya menunjukkan heteroskedastisitas. Jika diperlukan, terdapat pula uji-uji statistik yang dgunakan untuk mengecek asumsi ini seperti uji Breusch-Pagan, uji White dan sebagainya. Jika terdapat kondisi heteroskedastisitas, maka dapat juga dilakukan penanganan melalui transformasi data.

Pada kasus ini, berdasarkan uji BP (tidak tolak $H_0$) dan plot sebaran nilai residual terlihat relatif menyebar acak untuk setiap nilai prediksi dan tidak ada kecenderungan membentuk pola kerucut. Dengan demikian dapat dikatakan tidak terdapat heteroskedastisitas.

library(lmtest)   # paket untuk uji model linier

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.2.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# Melakukan uji Breusch-Pagan
breusch_pagan_test <- bptest(model_reg)

breusch_pagan_test

## 
##  studentized Breusch-Pagan test
## 
## data:  model_reg
## BP = 15.26, df = 15, p-value = 0.4329

# Membuat plot residual vs nilai prediksi dengan scatter plot merah
plot(fitted(model_reg), model.res, 
     xlab = "Fitted Values", ylab = "Residuals", 
     main = "Residual vs Fitted Values", 
     col="red", pch = 20, cex=1.5)

# Menambahkan garis putus-putus merah di nilai 0
abline(h = 0, col = "blue", lty = 2)

Tidak Ada Multikolinieritas

Multikolinieritas (Ini sebenarnya bukanlah sebuah asumsi) merupakan masalah yang perlu ditangani. Analisis model dengan keberadaan multikolinieritas bisa menjadi tidak tepat. Ragam penduga menjadi lebih besar sehingga juga akan memperlebar selang kepercayaannya. Dampak dari selang kepercayaan yang lebih lebar yaitu uji dapat menjadi tidak signifikan meskipun seharusnya variabel tersebut memiliki pengaruh kuat. Hal ini juga tidak menutup kemungkinan sampai menyebabkan penduga keofisien menjadi berubah tanda (positif menjadi negatif atau sebaliknya). Jika terjadi demikian tentu hasil analisis menjadi sangat tidak reliable.

Keberadaan masalah multikolinieritas dapat dilihat dengan cara berikut:

• Matriks korelasi, di mana batas koefisien korelasi 0,80 dapat menjadi ‘warning’. Pada bagian sebelumnya, kita sudah menyajikan matriks korelasi dan menunjukkan beberapa variabel penjelas dengan korelasi cukup tinggi.

• Variance Inflation Factor (VIF), dengan nilai VIF di atas 10 menunjukkan adanya masalah multikolinieritas.

# Memuat paket yang diperlukan
library(car)

## Warning: package 'car' was built under R version 4.2.3

## Loading required package: carData

# Menghitung VIF untuk setiap variabel dalam model
vif_values <- vif(model_reg)

data.frame(vif_values)

##                    vif_values
## PR_NO_LIS            1.674004
## PR_SAMPAH            1.752889
## PR_TINJA             3.883222
## PR_MKM_SUNGAI        2.519156
## PR_SUNGAI_LMBH       2.248597
## PR_KUMUH             2.945363
## PRA_1000             1.840095
## SD_1000              4.638975
## SM_1000              2.323137
## RS_PKS_PDK_1000      5.283640
## LIN_BID_POS_P_1000   3.867567
## APT_OBT_1000         3.014481
## DOK_DRG_1000         4.545848
## BID_1000             2.377227
## GZ_BURUK_1000        1.272896

Jika terdapat masalah multikolinieritas tersebut kita mungkin perlu mengeluarkan salah satu atau lebih variabel-variabel dengan korelasi yang kuat. Berdasarkan hasil di atas, ternyata tidak terdapat indikasi masalah multikolinieritas antar variabel penjelas. Hal ini dapat meyakinkan kita bahwa analisis rergesi linier berganda yang kita lakukan sudah memenuhi asumsi-asumsi serta tidak terdapat masalah yang dapat mengganggu hasil dan interpretasinya.