Relationship vs Causal Relationship
Tidak semua hubungan (relationship) berupa hubungan sebab-akibat
Penentuan suatu hubungan bersifat sebab-akibat memerlukan well-argued position dari bidang ilmu terkait.
Alat Analisis Keterkaitan
Ditentukan oleh :
Skala pengukuran data/variable
Jenis hubungan antar variabel
| Relationship | Numerik | Kategorik |
|---|---|---|
| Numerik | Korelasi Pearson, Spearman | Tabel Ringkasan Korelasi Biserial |
| Kategorik | Tabel Ringkasan Korelasi Biserial | Spearman (ordinal), Chi-Square Korelasi Tetrachoric |
Causal Relationship
Digunakan jika variabel response dan prediktor memiliki hubungan sebab-akibat.
| Y \ X | Numerik | Kategorik |
|---|---|---|
| Numerik | Regresi Linear | ANOVA |
| Kategorik |
|
|
Jika kedua variabel memiliki hubungan sebab akibat tetapi tidak linear, bukan berarti kedua variabel tersebut sepenuhnya tidak memiliki hubungan sebab akibat. Terdapat hubungan sebab akibat lain yang mungkin terjadi yaitu eksponensial, kubik, kuadratik, logaritma, dll. Tetapi pada hubungan sebab akibat jenis ini tidak digunakan istilah korelasi. Korelasi hanya digunakan untuk hubungan sebab akibat yang linear.
Koefisien Korelasi
tidak menggambarkan hubungan sebab akibat
nilainya berkisar antara -1 dan 1
tanda (+) / (-) menunjukkan arah hubungan yaitu (+) searah dan (-) berlawanan arah
Pearson’s Coef of Correlation menunjukkan linear relationship
Spearman’s Coef of Correlation (rank correlation) menunjukkan trend relationship
Linear relationship (Parametrik) -> Pearson Correlation
\[ r_{xy} = \frac {S_{xy}}{S_xS_y} \]
\[ S_{xy} = \frac {\sum (x_i - \bar{x})(y_i - \bar{y})}{n-1} \]
dengan \[ S_x = \sqrt {\frac {\sum (x_i - \bar{x})^2} {n-1}} \]
\[ S_y = \sqrt {\frac {\sum (y_i - \bar{y})^2} {n-1}} \]
Trend Relationship -> Rank Correlation -> Spearman Correlation (NonParametrik)
\[ \theta = \frac {\sum ((R_i - \bar{R})(S_i - \bar{S}))} {\sqrt {\sum \limits_i (R_i - \bar{R})^2 \sum (S_i - \bar{S})^2} } \]
R : peringkat dari X
S : peringkat dari Y
\[ \bar{R} : rataan \ peringkat \ X \]
\[ \bar{S} : rataan \ peringkat \ Y \]
Uji Korelasi
Angka korelasi yang diperoleh dari sejumlah pengamatan merupakan angka contoh; Angka korelasi yang sebenarnya sebagai suatu paramater tidak diketahui, kecuali apabila dilakukan sensus.
Uji Korelasi adalah uji atas keberartian angka korelasi, yaitu uji atas
\(H_0 : \rho = 0\)dan \(H_1 : \rho \neq 0\) tolak H0 apabila \(|t| > t_{\alpha /2, n-2}\) dengan \(t = \frac {r\sqrt {n-2}}{\sqrt {1-r^2}}\)
Mendapatkan model hubungan antar variabel
Menganalisis hubungan/pengaruh antara satu atau lebih variabel numerik terhadap sebuah variabel numerik lain.
Menduga nilai suatu variabel berdasarkan nilai variabel lainnya.
Model Umum: \(Y = b_0 + b_1 X_1 + b_2 X_2 + \ ...\ +b_kX_k\) dengan Y merupakan dependent variabel, X sebagai independect/predictor variables, dan b merupakan regression coeficients.
Metode Kuadrat Terkecil: pendugaan parameter pada regresi didapat dengan meminimumkan jumlah kuadrat galat.
(KEBAIKAN MODEL REGRESI) Dilihat dari nilai koefisien determinasi \(R^2\) merupakan ukuran seberapa besar keragaman dari peubah respon \(y\) dapat dijelaskan oleh model peubah penjelas \(x\). Nilainya antara \(0-100\%\) , semakin mendekati \(100\%\) maka nilai semakin bagus.
Cara meningkatkan \(R^2\) yaitu cek asumsi yang diperlukan dalam regresi melalui evaluasi residual (diagnostik terhadap plot-plot residual) dan tambah peubah penjelas -> analisis regresi linear berganda (multiple linear regression).
Asumsi - asumsi dalam Regresi
Nilai mean dari variabel-variabel Y dimodelkan secara akurat oleh fungsi linier dari variabel-variabel \(X\)
Istilah galat acak, \(\varepsilon\) , diasumsikan menyebar normal dengan nilai tengah nol dan memiliki ragam yang konstan, \(\sigma^2\)
Galat bersifat independen/saling bebas.
QQ Plot -> uji kenormalan sisaan
Homogenity Test -> uji kehomogenan ragam sisaan. Dengan H0 (sisaan homogen) dan H1 (sisaan tidak homogen).
Detecting Autocorrelation
The Durbin Watson Test \(d = \frac {\sum \limits_{t=2}^{t=n} {(\hat {u}_t - \hat{u}_{t-1})^2}} {\sum \limits_{t=1}^{t=n} {\hat {u_t}^2}}\)
It is simply the ration of the sum of squared differences in successive residuals to the RSS.
The number of observation is n-1 as one observation is lost in taking successive differences.
Model \(Y_i = \beta_0 + \beta_1 X_{1i} + \ ...\ \beta_k X_{ki} + \varepsilon_|\).
Banyaknya variabel dummy yang dibentuk dari \(k\) variabel kategorik adalah \(k-1\) dengan \(k\) adalah kategorik dalam variabel tersebut.
Pengujian Parameter
ANOVA digunakan untuk menguji secara simultan pengaruh seluruh X terhadap Y
Hipotesis Nol: Model regresi tidak sesuai untuk data jika dibandingkan dengan model baseline. Dan b1 = b2 = … = bk = 0.
Hipotesis Alternatif: Model regresi sesuai untuk data jika dibandingkan dengan model baseline. Dan Tidak semua b sama dengan nol
Jika nilai sig kecil, disimpulkan tolak H0. dengan kata lain, jika nilai sig kecil berarti ada X yang berpengaruh terhadap Y. Demikian pula sebaliknya. -> UJI Parsial -> uji-t
Multikolinearitas
Masalah multikolinieritas terjadi pada regresi berganda jika variabel-variabel X saling berkorelasi.
Hal ini akan mempengaruhi ragam dari dugaan koefisien regresi
Variabel X yang dianggap penting kemungkinan akan tidak signifikan
Pendugaan dari koefisien regresi menjadi tidak benar, misalnya koefisien memiliki tanda negatif padahal dalam hubungan x dan y sebenarnya adalah positif.
| Bagaimana cara mendeteksi multikolinearitas? | Bagaimana cara mengatasi mengatasi multikolimearitas? |
|---|---|
| Periksa korelasi antar variabel X | Jika kita ingin memilih variabel X dimana hanya X yang signifikan akan memasuki model: Gunakan prosedur penyeleksian variabel, seperti forward, backward, stepwise |
| Hitunglah nilai Variance Inflation Factor (VIF) dimana VIF = \((1-R_j^ 2)^{-1}\) . \(R_j^2\) adalah kuadrat dari regresi dimana \(Xj\) merupakan variabel respon dan variabel X lainnya menjadi predictor. | Jika kita ingin mempertahankan konfigurasi variabel X yang akan memasuki model: Gunakan metoda estimasi diluar metoda kuadrat terkecil, seperti Ridge Regression, Principal Component Regression, Partial Least Square. |
| Jika VIF lebih besar dari 10 biasanya ada masalah multikolinieritas. |
Metode Penyeleksian Variabel
Forward (Hanya masuk)
Backward (hanya keluar)
Stepwise (Masuk dan keluar)
library(readr)
library(corrplot)
## corrplot 0.92 loaded
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Loading required package: carData
library(MASS)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ purrr 1.0.2
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ dplyr::select() masks MASS::select()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
setwd("D:/Kuliah/Mat/TSA Kominfo/Praktikum")
data <- read.csv("4. limit_kredit.csv")
head(data, 15)
## NAMA EVENT STATUS JK USIA STATUS_NIKAH JML_TANGGUNGAN
## 1 A123145 36 NEW LAKI LAKI 44 MENIKAH 2
## 2 A118717 48 NEW LAKI LAKI 53 MENIKAH 2
## 3 A125680 24 NEW PEREMPUAN 28 MENIKAH 2
## 4 A14798 60 NEW LAKI LAKI 45 MENIKAH 3
## 5 A142852 24 NEW PEREMPUAN 44 MENIKAH 1
## 6 A96802 24 NEW PEREMPUAN 51 MENIKAH 3
## 7 A10617 60 NEW LAKI LAKI 28 TIDAK MENIKAH 0
## 8 A134638 36 NEW PEREMPUAN 32 MENIKAH 0
## 9 A147552 60 NEW LAKI LAKI 43 MENIKAH 1
## 10 A43766 18 NEW LAKI LAKI 56 MENIKAH 3
## 11 A105309 36 NEW LAKI LAKI 29 MENIKAH 1
## 12 A257422 48 NEW LAKI LAKI 44 MENIKAH 4
## 13 A49034 60 NEW LAKI LAKI 41 MENIKAH 2
## 14 A54294 54 NEW PEREMPUAN 53 MENIKAH 2
## 15 A313803 6 TOP LAKI LAKI 27 TIDAK MENIKAH 0
## STATUS_TT PENDIDIKAN WILAYAH PEKERJAAN LAMA_BEKERJA
## 1 Milik Sendiri SMA Balikpapan Pegawai swasta 14
## 2 Milik Sendiri SMA Balikpapan Pegawai BUMN 25
## 3 Milik Sendiri SMA Balikpapan Pegawai swasta 6
## 4 Milik Sendiri SMA Medan Pegawai swasta 15
## 5 Milik Orang Tua S1 Jakarta 5 Pegawai swasta 3
## 6 Milik Sendiri SMA Jakarta 4 PNS 10
## 7 Rumah Dinas DIPLOMA Medan Pegawai BUMN 5
## 8 Milik Orang Tua S1 Jakarta 5 Pegawai swasta 1
## 9 Milik Sendiri SMA Jakarta 5 Pegawai swasta 3
## 10 Milik Sendiri SMA Palembang PNS 10
## 11 Milik Orang Tua SMA Balikpapan Pegawai swasta 2
## 12 Milik Sendiri SMA Jakarta 5 Pegawai swasta 3
## 13 Sewa SMA Jakarta 3 Pegawai swasta 13
## 14 Milik Sendiri S1 Jakarta 3 Pegawai BUMN 25
## 15 Milik Sendiri DIPLOMA Semarang Pegawai swasta 4
## PENGHASILAN DBR WAKTU_KREDIT LIMIT_KREDIT STATUS_KREDIT PENDIDIKAN_KODE
## 1 2829819 0.45 36 23000000 LANCAR 3
## 2 23757496 0.37 48 300000000 LANCAR 3
## 3 2733000 0.45 24 16000000 LANCAR 3
## 4 3000000 0.42 60 50000000 LANCAR 3
## 5 4195200 0.35 24 28000000 LANCAR 5
## 6 8883124 0.49 24 90000000 LANCAR 3
## 7 4943520 0.51 60 100000000 LANCAR 4
## 8 3749812 0.23 36 18500000 LANCAR 5
## 9 5000000 0.20 60 40000000 LANCAR 3
## 10 10672327 0.15 18 25000000 LANCAR 3
## 11 5368339 0.28 36 35000000 LANCAR 3
## 12 4257454 0.35 48 50000000 MACET 3
## 13 4088257 0.47 60 63000000 LANCAR 3
## 14 13329450 0.55 54 300000000 LANCAR 5
## 15 3267000 0.38 60 36700000 LANCAR 4
## PENDIDIKAN2
## 1 3. SMA
## 2 3. SMA
## 3 3. SMA
## 4 3. SMA
## 5 5. S1
## 6 3. SMA
## 7 4. DIPLOMA
## 8 5. S1
## 9 3. SMA
## 10 3. SMA
## 11 3. SMA
## 12 3. SMA
## 13 3. SMA
## 14 5. S1
## 15 4. DIPLOMA
str(data)
## 'data.frame': 500 obs. of 19 variables:
## $ NAMA : chr "A123145" "A118717" "A125680" "A14798" ...
## $ EVENT : int 36 48 24 60 24 24 60 36 60 18 ...
## $ STATUS : chr "NEW" "NEW" "NEW" "NEW" ...
## $ JK : chr "LAKI LAKI" "LAKI LAKI" "PEREMPUAN" "LAKI LAKI" ...
## $ USIA : int 44 53 28 45 44 51 28 32 43 56 ...
## $ STATUS_NIKAH : chr "MENIKAH" "MENIKAH" "MENIKAH" "MENIKAH" ...
## $ JML_TANGGUNGAN : int 2 2 2 3 1 3 0 0 1 3 ...
## $ STATUS_TT : chr "Milik Sendiri" "Milik Sendiri" "Milik Sendiri" "Milik Sendiri" ...
## $ PENDIDIKAN : chr "SMA" "SMA" "SMA" "SMA" ...
## $ WILAYAH : chr "Balikpapan" "Balikpapan" "Balikpapan" "Medan" ...
## $ PEKERJAAN : chr "Pegawai swasta" "Pegawai BUMN" "Pegawai swasta" "Pegawai swasta" ...
## $ LAMA_BEKERJA : int 14 25 6 15 3 10 5 1 3 10 ...
## $ PENGHASILAN : int 2829819 23757496 2733000 3000000 4195200 8883124 4943520 3749812 5000000 10672327 ...
## $ DBR : num 0.45 0.37 0.45 0.42 0.35 0.49 0.51 0.23 0.2 0.15 ...
## $ WAKTU_KREDIT : int 36 48 24 60 24 24 60 36 60 18 ...
## $ LIMIT_KREDIT : int 23000000 300000000 16000000 50000000 28000000 90000000 100000000 18500000 40000000 25000000 ...
## $ STATUS_KREDIT : chr "LANCAR" "LANCAR" "LANCAR" "LANCAR" ...
## $ PENDIDIKAN_KODE: int 3 3 3 3 5 3 4 5 3 3 ...
## $ PENDIDIKAN2 : chr "3. SMA" "3. SMA" "3. SMA" "3. SMA" ...
# Karena terdapat peubah kategorik yang tersimpan sebagai numerik
# maka sebaiknya dikonversi menjadi faktor
data$PENDIDIKAN_KODE <- as.factor(data$PENDIDIKAN_KODE)
# Summary data
summary(data)
## NAMA EVENT STATUS JK
## Length:500 Min. : 6.00 Length:500 Length:500
## Class :character 1st Qu.: 24.00 Class :character Class :character
## Mode :character Median : 36.00 Mode :character Mode :character
## Mean : 45.34
## 3rd Qu.: 60.00
## Max. :120.00
## USIA STATUS_NIKAH JML_TANGGUNGAN STATUS_TT
## Min. :22.00 Length:500 Min. :0.000 Length:500
## 1st Qu.:30.00 Class :character 1st Qu.:0.000 Class :character
## Median :34.00 Mode :character Median :1.000 Mode :character
## Mean :36.54 Mean :1.472
## 3rd Qu.:43.00 3rd Qu.:2.000
## Max. :59.00 Max. :5.000
## PENDIDIKAN WILAYAH PEKERJAAN LAMA_BEKERJA
## Length:500 Length:500 Length:500 Min. : 1.000
## Class :character Class :character Class :character 1st Qu.: 3.000
## Mode :character Mode :character Mode :character Median : 6.000
## Mean : 7.892
## 3rd Qu.:11.000
## Max. :29.000
## PENGHASILAN DBR WAKTU_KREDIT LIMIT_KREDIT
## Min. : 2005020 Min. :0.0400 Min. : 12.00 Min. : 5000000
## 1st Qu.: 3437854 1st Qu.:0.2500 1st Qu.: 36.00 1st Qu.: 25000000
## Median : 4594916 Median :0.3400 Median : 48.00 Median : 40000000
## Mean : 6387706 Mean :0.3482 Mean : 55.54 Mean : 74359000
## 3rd Qu.: 7021339 3rd Qu.:0.4400 3rd Qu.: 60.00 3rd Qu.: 95500000
## Max. :65459300 Max. :0.6000 Max. :120.00 Max. :517000000
## STATUS_KREDIT PENDIDIKAN_KODE PENDIDIKAN2
## Length:500 1: 1 Length:500
## Class :character 2: 4 Class :character
## Mode :character 3:247 Mode :character
## 4: 95
## 5:144
## 6: 9
Untuk melihat korelasi antara 2 peubah dapat menggunakan
fungsi cor. Jika ingin menghitung korelasi Spearman
khususnya pada data rank, maka dapat mengatur argumen untuk
parameter method dengan nilai "spearman".
# Melihat korelasi antara 2 variabel
cat("*************** KORELASI ***************\n\n")
## *************** KORELASI ***************
cat("Penghasilan & Limit Kredit -->",
round(cor(data$PENGHASILAN, data$LIMIT_KREDIT), 3), "\n")
## Penghasilan & Limit Kredit --> 0.628
cat("Penghasilan & Lama Bekerja -->",
round(cor(data$PENGHASILAN, data$LAMA_BEKERJA), 3), "\n")
## Penghasilan & Lama Bekerja --> 0.185
cat("Waktu kredit & Limit Kredit -->",
round(cor(data$WAKTU_KREDIT, data$LIMIT_KREDIT), 3), "\n")
## Waktu kredit & Limit Kredit --> 0.593
cat("Pendidikan & Penghasilan -->",
round(cor(as.integer(data$PENDIDIKAN_KODE), data$PENGHASILAN, method ="spearman"), 3), "\n")
## Pendidikan & Penghasilan --> 0.222
cat("Pendidikan & Limit Kredit -->",
round(cor(as.integer(data$PENDIDIKAN_KODE), data$LIMIT_KREDIT, method ="spearman"), 3),"\n")
## Pendidikan & Limit Kredit --> 0.124
cat("Usia & Lama Bekerja -->",
round(cor(data$USIA, data$LAMA_BEKERJA), 3), "\n")
## Usia & Lama Bekerja --> 0.653
cat("\n****************************************\n")
##
## ****************************************
Fungsi cor juga dapat digunakan untuk mengukur korelasi
dari banyak variabel sekaligus (numerik) dimana output yang dihasilkan
berupa matriks korelasi.
# Membuat matriks Korelasi untuk setiap variabel numerik
library(tidyverse)
num.data <-data %>% select_if(is.numeric) # Mengambil kolom-kolom numerik saja
corr_matrix <- cor(num.data)
round(corr_matrix, 3)
## EVENT USIA JML_TANGGUNGAN LAMA_BEKERJA PENGHASILAN DBR
## EVENT 1.000 0.070 0.055 0.111 0.095 0.122
## USIA 0.070 1.000 0.618 0.653 0.301 -0.088
## JML_TANGGUNGAN 0.055 0.618 1.000 0.437 0.159 -0.116
## LAMA_BEKERJA 0.111 0.653 0.437 1.000 0.185 -0.021
## PENGHASILAN 0.095 0.301 0.159 0.185 1.000 -0.147
## DBR 0.122 -0.088 -0.116 -0.021 -0.147 1.000
## WAKTU_KREDIT 0.578 -0.047 0.016 0.096 0.080 0.372
## LIMIT_KREDIT 0.339 0.197 0.089 0.243 0.628 0.344
## WAKTU_KREDIT LIMIT_KREDIT
## EVENT 0.578 0.339
## USIA -0.047 0.197
## JML_TANGGUNGAN 0.016 0.089
## LAMA_BEKERJA 0.096 0.243
## PENGHASILAN 0.080 0.628
## DBR 0.372 0.344
## WAKTU_KREDIT 1.000 0.593
## LIMIT_KREDIT 0.593 1.000
library(corrplot)
# Menampilkan Nilai Korelasi dalam bentuk visual
corrplot(corr_matrix, method="circle", type="lower", diag=F)
# Menampilkan korelasi untuk beberapa variabel
subs <- c("PENGHASILAN", "LAMA_BEKERJA", "LIMIT_KREDIT", "USIA")
cor(data[subs])
## PENGHASILAN LAMA_BEKERJA LIMIT_KREDIT USIA
## PENGHASILAN 1.0000000 0.1845696 0.6275672 0.3013116
## LAMA_BEKERJA 0.1845696 1.0000000 0.2429656 0.6534017
## LIMIT_KREDIT 0.6275672 0.2429656 1.0000000 0.1972260
## USIA 0.3013116 0.6534017 0.1972260 1.0000000
cat("\n")
corrplot(cor(data[subs]), method="circle", type="lower", diag=F)
Model regresi linear sederhana yaitu model regresi dimana terdapat 1 peubah respon dengan hanya 1 peubah penjelas. Contohnya, model untuk menghitung LIMIT KREDIT (Sebagai peubah respon) berdasarkan informasi mengenai PENGHASILAN (Sebagai peubah penjelas : X)
Untuk menyederhanakan data yang digunakan pada model regresi ini,
maka kita dapat menyeleksi kolom-kolom yang akan digunakan saja dan
menyimpannya pada variabel baru data.reg.1.
(Proses ini opsional, kita dapat saja tetap menggunakan data secara keseluruhan walaupun model nantinya hanya menggunakan 1 peubah bebas)
# standar R
data.reg.1 <- data[c("LIMIT_KREDIT", "PENGHASILAN")]
# # Slicing menggunakan dplyr
# data.reg.1 <- data %>% select(LIMIT_KREDIT, PENGHASILAN)
head(data.reg.1)
## LIMIT_KREDIT PENGHASILAN
## 1 23000000 2829819
## 2 300000000 23757496
## 3 16000000 2733000
## 4 50000000 3000000
## 5 28000000 4195200
## 6 90000000 8883124
Model regresi linear sederhana yang pertama yaitu menggunakan data apa adanya, baik untuk peubah LIMIT_KREDIT maupun PENGHASILAN.
Membuat Model
Untuk membuat model regresi linear dapat menggunakan fungsi built-in
pada R yaitu lm. Argumen pada parameter pertama adalah
formula yang menunjukkan model regresi yang akan dibuat.
Misal, pada contoh berikut ini,
LIMIT_KREDIT ~ PENGHASILAN menunjukkan bahwa model yang
dibuat yaitu model linear dengan peubah respon adalah
LIMIT_KREDIT dan peubah penjelas adalah
PENGHASILAN. Sementara data yang digunakan adalah
data.reg.1
# Model Regresi Lienar Sederhana (Y : Limit Kredit, X : Penghasilan)
reg.1 = lm(LIMIT_KREDIT ~ PENGHASILAN, data = data.reg.1)
summary(reg.1)
##
## Call:
## lm(formula = LIMIT_KREDIT ~ PENGHASILAN, data = data.reg.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -241548318 -30563160 -19034577 15820323 253846561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.101e+07 4.059e+06 5.177 3.28e-07 ***
## PENGHASILAN 8.351e+00 4.643e-01 17.988 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 61980000 on 498 degrees of freedom
## Multiple R-squared: 0.3938, Adjusted R-squared: 0.3926
## F-statistic: 323.6 on 1 and 498 DF, p-value: < 2.2e-16
Berdasarkan summary model, dapat diketahui bahwa
peubah PENGHASILAN memiliki pengaruh yang signifikan di
dalam model (P-value sangat kecil ***). Namun, jika dilihat dari
nilai R2, sepertinya performa model belum cukup baik dimana
Peubah PENGHASILAN hanya mampu menjelaskan sekitar 39%
keragaman dari peubah respon.
Scatter Plot dan Garis Regresi
library(ggplot2)
# Membuat scatter plot
scatter_plot <- ggplot(data = data.reg.1, aes(x = PENGHASILAN, y = LIMIT_KREDIT)) +
# Menambahkan titik-titik scatter plot
geom_point(color="darkred", size=4, alpha=0.3) +
# Menambahkan garis regresi
geom_smooth(method = "lm",
formula = y ~ x,
se = T, # Menambahkan C.I
color = "blue")
# Menampilkan scatter plot
scatter_plot
Jika melihat scatter Plot, tidak tergambarkans ecara jelas hubungan
antara peubah PENGHASILAN dan LIMIT_KREDIT.
Sebagian besar data menumpuk pada nilai-nilai kecil, namun ada data-data
lainnya yang menyebar pada rentang nilai yang sangat besar. Hal ini
mengindikasikan juga bahwa terdapat kemenjuluran yang besar pada kedua
peubah serta kemungkinan besar asumsi-asumsi model linear tidak
terpenuhi.
Plot Normal Q-Q
Dari Normal Q-Q Plot maupun uji formal menggunakan Shapiro-Wilk maka dapat disimpulkan dengan sangat jelas bahwa Residual model tidak menyebar Normal.
Asumsi Normalitas pada model regresi linear : Residual menyebar Normal(0,σ2)
# Menghitung nilai residual dari model reg.1
res <- residuals(reg.1)
# Membuat Kurva Normal QQ Plot
qqnorm(res)
qqline(res, col="red")
# Mengecek apakah residual menyebar Normal
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.84589, p-value < 2.2e-16
Plot Sisaan (Residual Plot)
Asumsi lainnya di dalam model regresi linear adalah Homegenitas yaitu ragam dari residual harus konstan atau homogen.
Berdasarkan plot sisaan dapat kita lihat pula bahwa nilai-nilai residual tidak menyebar secara acak, dimana terdapat pola seperti corong yang mengindikasikan semakin besar nilai Y maka residual-nya juga semakin besar.
fits<-fitted(reg.1)
plot(fits, res, col="orange", cex=1.5, lwd=2)
abline(h = 0, col = "red", lty = 5)
Pengecekan Homogenitas
Terkait dengan pola sebaran residual, dengan menggunakan Uji Formal Breusch-Pagan, Goldfeld-Quant, serta Non-Constant Variance Score dapat disimpulkan bahwa terdapat masalah heteroskedastisitas atau ragam residual tidak homogen.
# Breusch-Pagan test (Tolak Ho : Heteroskedastisitas)
lmtest::bptest(reg.1)
##
## studentized Breusch-Pagan test
##
## data: reg.1
## BP = 160.94, df = 1, p-value < 2.2e-16
cat("\n")
# Goldfeld-Quandt test (Tolak Ho : Heteroskedastisitas)
lmtest::gqtest(reg.1)
##
## Goldfeld-Quandt test
##
## data: reg.1
## GQ = 1.4721, df1 = 248, df2 = 248, p-value = 0.001213
## alternative hypothesis: variance increases from segment 1 to 2
cat("\n")
# ncvTest
car::ncvTest(reg.1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 464.2134, Df = 1, p = < 2.22e-16
Plot-plot rsidual dapat juga ditampilkan hanya dengan menggunakan fungsi `plot` dengan memasukkan model yang dibuat sebagai nilai argumennya.
# dengan memasukkan model regresi sebagai argumen fungsi plot
# dapat menampilkan beberapa plot residual secara langsung
plot(reg.1)
Autokorelasi
Menggunakan uji Durbin-Watson dapat disimpulkan tidka terdapat autokorelasi pada data
lmtest::dwtest(reg.1)
##
## Durbin-Watson test
##
## data: reg.1
## DW = 1.9919, p-value = 0.4643
## alternative hypothesis: true autocorrelation is greater than 0
Mengingat sebaran data pada dua peubah dengan kemenjuluran ke kanan yang sangat besar, maka dapat dicobakan transformasi: khususnya untuk peubah Y (misal log natural).
Pada contoh ini akan dilakukan transformasi selain untuk peubah Y, juga peubah X dengan harapan dapat meningkatkan performa model.
Membuat Model
data.reg.2 <- log(data.reg.1)
reg.2 <- lm(LIMIT_KREDIT~PENGHASILAN, data=data.reg.2)
summary(reg.2)
##
## Call:
## lm(formula = LIMIT_KREDIT ~ PENGHASILAN, data = data.reg.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.22586 -0.44672 -0.05501 0.51649 1.86892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.34445 0.83985 2.792 0.00545 **
## PENGHASILAN 0.99311 0.05427 18.299 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6872 on 498 degrees of freedom
## Multiple R-squared: 0.4021, Adjusted R-squared: 0.4009
## F-statistic: 334.9 on 1 and 498 DF, p-value: < 2.2e-16
Hasil pemodelan regresi linear sederhana dengan peubah yang sudah ditrasnformasi secara umum tidak berbeda dengan sebelumnya, dimana nilai R2 hanya sidikit lebih besar dari sebelumnya yaitu sebesar 0,4021.
Scatter Plot dan Garis Regresi
Melihat scatterplot data hasil transformasi, pola hubungan linear dari kedua peubah menjadi lebih jelas terlihat walaupun memang masih cukup menyebar.
# Membuat scatter plot dengan ggplot2
scatter_plot <- ggplot(data = data.reg.2, aes(x = PENGHASILAN, y = LIMIT_KREDIT)) +
# Menambahkan titik-titik scatter plot
geom_point(color="orange", size=4, alpha=0.6) +
# Menambahkan garis regresi
geom_smooth(method = "lm",
formula = y ~ x,
se = T, # Menambahkan C.I
color = "blue")
# Menampilkan scatter plot
scatter_plot
Plot Normal Q-Q
Hasil uji Shapiro-Wilk pada model dengan data yang sudah ditransformasi masih menunjukkan residual tidak menyebar Normal. Namun jika melihat dari Normal Q-Q Plot, sebaran residual lebih mendekati Normal dibandingkan sebelumnya.
# Menghitung nilai residual dari model reg.2
res <- residuals(reg.2)
# Membuat Kurva Normal QQ Plot
qqnorm(res)
qqline(res, col="red")
# Mengecek apakah residual menyebar Normal
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98735, p-value = 0.0002497
Plot Sisaan (Residual Plot)
fits<-fitted(reg.2)
plot(fits, res, col="orange", cex=1.5, lwd=2)
abline(h = 0, col = "red", lty = 5)
Pengecekan Homogenitas
hasil uji Breusch-Pagan, Goldfeld-Quandt serta Non-constant Variance Score menunjukkan masih terdapat masalah Heteroskedastisitas atau ragam yang tidak homogen. Jika melihat dari plot residual, nilai-nilai residual sudah terlihat menyebar acak tanpa adanya pola-pola tertentu.
# Breusch-Pagan test (Tolak Ho : Heteroskedastisitas)
lmtest::bptest(reg.1)
##
## studentized Breusch-Pagan test
##
## data: reg.1
## BP = 160.94, df = 1, p-value < 2.2e-16
cat("\n")
# Goldfeld-Quandt test (Tolak Ho : Heteroskedastisitas)
lmtest::gqtest(reg.1)
##
## Goldfeld-Quandt test
##
## data: reg.1
## GQ = 1.4721, df1 = 248, df2 = 248, p-value = 0.001213
## alternative hypothesis: variance increases from segment 1 to 2
cat("\n")
# ncvTest
car::ncvTest(reg.1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 464.2134, Df = 1, p = < 2.22e-16
Plot-plot residual menggunakan fungsi plot:
plot(reg.2)
Autokorelasi
lmtest::dwtest(reg.2)
##
## Durbin-Watson test
##
## data: reg.2
## DW = 1.9595, p-value = 0.325
## alternative hypothesis: true autocorrelation is greater than 0
Selain inferensi, tujuan lain dari pembuatan model adalah untuk melakukan prediksi. Dalam konteks machine learning, sering kali performa model dalam memprediksi menjadi perhatian yang lebih utama dibandingkan inferensi.
Untuk melakukan predisi dapat menggunakan fungsi
predict, dimana argumen pertama adalah model yang digunakan
untuk memprediksi, serta argumen kedua, yaitu data-data peubah X yang
akan menjadi dasar untuk memprediksi nilai peubah Y.
# PREDICT
# Misal terdapat dua orang dengan PENGASILAN SEBESAR Rp 15.000.000 dan Rp 9.000.000,
# berdasarkan 3 model yang sudah dibuat, berapakah prediksi LIMIT KREDIT yang dapat diberikan
income <- data.frame("PENGHASILAN"=c(15000000, 9000000))
# Model Reg.1
limit.1 <- predict(reg.1, income)
print(t(round(limit.1)))
## 1 2
## [1,] 146281380 96174605
# Orang 1 : 146.281.380
# Orang 2 : 96.174.605
Sedikit berbeda dengan model 1, Nilai peubah Y pada model 2 adalah
hasil transformasi logaritma. Sehingga output dari fungsi
predict juga akan dalam bentuk nilai logaritmanya. Oleh
karena itu, untuk mendapatkan nilai aslinya (Misal: Rupiah) maka nilai Y
harus di-inverse. Pada contoh ini inverse dari log natural adalah
eksponensial (eypred).
# Transformasi data sesuai model
log.income <- log(income)
# Model Reg.2
limit.2 <- predict(reg.2, log.income)
# Karena peubah Y ditransformasi
# untuk memperoleh hasil dalam Rp harus di lakukan inverse : exp(Y)
limit.2.rp <- exp(limit.2)
print(t(round(limit.2.rp)))
## 1 2
## [1,] 139586974 84047379
# Orang 1 : 139.586.974
# Orang 2 : 84.047.379
Di dalam contoh sebelumnya, nilai \(R^2\) model tidak begitu tinggi, hal ini mengindikasikan perlunya tambahan peubah penjelas lain agar dapat menjelaskan keragaman peubah Y dengan lebih baik. Model regresi linear dengan lebih dari 1 peubah penjelas disebut model regresi linear berganda.
Pada contoh ini akan digunakan 8 peubah penjelas seperti yang tersaji pada sintaks berikut ini.
columns <- c("LIMIT_KREDIT", "PENGHASILAN", "PEKERJAAN", "LAMA_BEKERJA", "STATUS_TT",
"USIA", "STATUS_NIKAH", "JML_TANGGUNGAN", "WAKTU_KREDIT")
# slicing data yang diperlukan
data.mul.1 <- data[columns]
head(data.mul.1)
## LIMIT_KREDIT PENGHASILAN PEKERJAAN LAMA_BEKERJA STATUS_TT USIA
## 1 23000000 2829819 Pegawai swasta 14 Milik Sendiri 44
## 2 300000000 23757496 Pegawai BUMN 25 Milik Sendiri 53
## 3 16000000 2733000 Pegawai swasta 6 Milik Sendiri 28
## 4 50000000 3000000 Pegawai swasta 15 Milik Sendiri 45
## 5 28000000 4195200 Pegawai swasta 3 Milik Orang Tua 44
## 6 90000000 8883124 PNS 10 Milik Sendiri 51
## STATUS_NIKAH JML_TANGGUNGAN WAKTU_KREDIT
## 1 MENIKAH 2 36
## 2 MENIKAH 2 48
## 3 MENIKAH 2 24
## 4 MENIKAH 3 60
## 5 MENIKAH 1 24
## 6 MENIKAH 3 24
Formula LIMIT_KREDIT ~ . menunjukkan bahwa model regrsi
yang dibuat yaitu menggunakan peubah LIMIT_KREDIT sebagai
peubah respon, dan tanda “titik” menunjukkan bahwa seluruh peubah
lainnya pada data.mul.1 digunakan sebagai peubah
penjelas.
Multikolinearitas
Untuk melihat ada tidaknya kondisi multikolinearitas pada peubah-peubah penjelas dapat menggunakan ukuran VIF (Variance Inflation Factor). Nilai VIF di atas 5 (atau referensi lain menyebut di atas 10) mengindikasikan adanya kondisi multikolinearitas.
Untuk menghitung nilai VIF dalam bahasa R, dapat menggunakan fungsi
vif yang tersedia pada package car
mul.reg.1 <- lm(LIMIT_KREDIT ~ ., data=data.mul.1)
summary(mul.reg.1)
##
## Call:
## lm(formula = LIMIT_KREDIT ~ ., data = data.mul.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -205151888 -19357136 -1433490 15484055 185612046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.035e+08 3.624e+07 -2.856 0.00448 **
## PENGHASILAN 7.568e+00 3.621e-01 20.898 < 2e-16 ***
## PEKERJAANPegawai swasta -1.530e+07 7.305e+06 -2.094 0.03678 *
## PEKERJAANPNS -2.725e+07 9.223e+06 -2.954 0.00329 **
## LAMA_BEKERJA 1.177e+06 3.970e+05 2.963 0.00319 **
## STATUS_TTMilik Orang Tua 1.717e+07 3.252e+07 0.528 0.59771
## STATUS_TTMilik Sendiri 1.602e+07 3.261e+07 0.491 0.62344
## STATUS_TTRumah Dinas 3.044e+07 3.312e+07 0.919 0.35852
## STATUS_TTSewa 2.119e+07 3.277e+07 0.647 0.51820
## USIA 3.648e+05 3.638e+05 1.003 0.31656
## STATUS_NIKAHMENIKAH 2.546e+07 2.102e+07 1.211 0.22655
## STATUS_NIKAHTIDAK MENIKAH 2.913e+07 2.138e+07 1.362 0.17375
## JML_TANGGUNGAN -3.824e+06 2.251e+06 -1.699 0.09004 .
## WAKTU_KREDIT 1.500e+06 7.120e+04 21.064 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43310000 on 486 degrees of freedom
## Multiple R-squared: 0.7111, Adjusted R-squared: 0.7034
## F-statistic: 92.04 on 13 and 486 DF, p-value: < 2.2e-16
car::vif(mul.reg.1)
## GVIF Df GVIF^(1/(2*Df))
## PENGHASILAN 1.245991 1 1.116240
## PEKERJAAN 1.247807 2 1.056907
## LAMA_BEKERJA 1.867946 1 1.366728
## STATUS_TT 1.612642 4 1.061554
## USIA 2.600554 1 1.612623
## STATUS_NIKAH 2.057287 2 1.197633
## JML_TANGGUNGAN 2.157633 1 1.468888
## WAKTU_KREDIT 1.112715 1 1.054853
plot(mul.reg.1)
Stepwise regression adalah suatu teknik dalam analisis regresi yang digunakan untuk memilih peubah-peubah penjelas yang akan dimasukkan atau dihapus dari model regresi berdasarkan sejumlah kriteria statistik tertentu. Tujuan dari stepwise regression adalah untuk membangun model regresi yang paling relevan dan efisien dengan memilih subset terbaik dari variabel independen yang tersedia. Teknik ini umumnya digunakan untuk menghindari overfitting dan untuk meningkatkan interpretabilitas model.
stepwise.reg <- MASS::stepAIC(mul.reg.1, direction="both")
## Start: AIC=17597.77
## LIMIT_KREDIT ~ PENGHASILAN + PEKERJAAN + LAMA_BEKERJA + STATUS_TT +
## USIA + STATUS_NIKAH + JML_TANGGUNGAN + WAKTU_KREDIT
##
## Df Sum of Sq RSS AIC
## - STATUS_TT 4 8.9351e+15 9.2070e+17 17595
## - STATUS_NIKAH 2 3.6740e+15 9.1544e+17 17596
## - USIA 1 1.8857e+15 9.1365e+17 17597
## <none> 9.1176e+17 17598
## - JML_TANGGUNGAN 1 5.4127e+15 9.1717e+17 17599
## - PEKERJAAN 2 1.6486e+16 9.2825e+17 17603
## - LAMA_BEKERJA 1 1.6474e+16 9.2824e+17 17605
## - PENGHASILAN 1 8.1930e+17 1.7311e+18 17916
## - WAKTU_KREDIT 1 8.3243e+17 1.7442e+18 17920
##
## Step: AIC=17594.65
## LIMIT_KREDIT ~ PENGHASILAN + PEKERJAAN + LAMA_BEKERJA + USIA +
## STATUS_NIKAH + JML_TANGGUNGAN + WAKTU_KREDIT
##
## Df Sum of Sq RSS AIC
## - USIA 1 2.1840e+15 9.2288e+17 17594
## - STATUS_NIKAH 2 5.9488e+15 9.2665e+17 17594
## <none> 9.2070e+17 17595
## - JML_TANGGUNGAN 1 5.3533e+15 9.2605e+17 17596
## + STATUS_TT 4 8.9351e+15 9.1176e+17 17598
## - PEKERJAAN 2 1.7196e+16 9.3789e+17 17600
## - LAMA_BEKERJA 1 1.4886e+16 9.3558e+17 17601
## - WAKTU_KREDIT 1 8.3012e+17 1.7508e+18 17914
## - PENGHASILAN 1 8.5573e+17 1.7764e+18 17921
##
## Step: AIC=17593.83
## LIMIT_KREDIT ~ PENGHASILAN + PEKERJAAN + LAMA_BEKERJA + STATUS_NIKAH +
## JML_TANGGUNGAN + WAKTU_KREDIT
##
## Df Sum of Sq RSS AIC
## - STATUS_NIKAH 2 5.2503e+15 9.2813e+17 17593
## - JML_TANGGUNGAN 1 3.5920e+15 9.2647e+17 17594
## <none> 9.2288e+17 17594
## + USIA 1 2.1840e+15 9.2070e+17 17595
## + STATUS_TT 4 9.2334e+15 9.1365e+17 17597
## - PEKERJAAN 2 1.7272e+16 9.4015e+17 17599
## - LAMA_BEKERJA 1 2.9461e+16 9.5234e+17 17608
## - WAKTU_KREDIT 1 8.4366e+17 1.7665e+18 17917
## - PENGHASILAN 1 9.0368e+17 1.8266e+18 17933
##
## Step: AIC=17592.67
## LIMIT_KREDIT ~ PENGHASILAN + PEKERJAAN + LAMA_BEKERJA + JML_TANGGUNGAN +
## WAKTU_KREDIT
##
## Df Sum of Sq RSS AIC
## <none> 9.2813e+17 17593
## + STATUS_NIKAH 2 5.2503e+15 9.2288e+17 17594
## + USIA 1 1.4855e+15 9.2665e+17 17594
## + STATUS_TT 4 1.1330e+16 9.1680e+17 17595
## - JML_TANGGUNGAN 1 7.9330e+15 9.3606e+17 17595
## - PEKERJAAN 2 1.8892e+16 9.4702e+17 17599
## - LAMA_BEKERJA 1 2.7654e+16 9.5579e+17 17605
## - WAKTU_KREDIT 1 8.5298e+17 1.7811e+18 17917
## - PENGHASILAN 1 9.0542e+17 1.8336e+18 17931
# Menampilkan summary model hasil stepwise
summary(stepwise.reg)
##
## Call:
## lm(formula = LIMIT_KREDIT ~ PENGHASILAN + PEKERJAAN + LAMA_BEKERJA +
## JML_TANGGUNGAN + WAKTU_KREDIT, data = data.mul.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -211420202 -20222348 -1895159 15305945 184188958
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.328e+07 9.333e+06 -4.637 4.54e-06 ***
## PENGHASILAN 7.616e+00 3.473e-01 21.930 < 2e-16 ***
## PEKERJAANPegawai swasta -1.779e+07 7.239e+06 -2.458 0.014333 *
## PEKERJAANPNS -2.865e+07 9.201e+06 -3.113 0.001957 **
## LAMA_BEKERJA 1.278e+06 3.335e+05 3.833 0.000143 ***
## JML_TANGGUNGAN -3.527e+06 1.718e+06 -2.053 0.040624 *
## WAKTU_KREDIT 1.464e+06 6.879e+04 21.286 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43390000 on 493 degrees of freedom
## Multiple R-squared: 0.706, Adjusted R-squared: 0.7024
## F-statistic: 197.3 on 6 and 493 DF, p-value: < 2.2e-16
mul.reg.log <- lm(log(LIMIT_KREDIT) ~ ., data=data.mul.1)
summary(mul.reg.log)
##
## Call:
## lm(formula = log(LIMIT_KREDIT) ~ ., data = data.mul.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.68710 -0.28196 0.01652 0.30232 1.18451
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.601e+01 4.087e-01 39.168 < 2e-16 ***
## PENGHASILAN 6.281e-08 4.083e-09 15.383 < 2e-16 ***
## PEKERJAANPegawai swasta -2.833e-01 8.236e-02 -3.440 0.000632 ***
## PEKERJAANPNS -2.165e-01 1.040e-01 -2.082 0.037860 *
## LAMA_BEKERJA 7.325e-03 4.477e-03 1.636 0.102411
## STATUS_TTMilik Orang Tua -2.585e-01 3.666e-01 -0.705 0.481073
## STATUS_TTMilik Sendiri -2.330e-01 3.677e-01 -0.634 0.526626
## STATUS_TTRumah Dinas -6.657e-02 3.734e-01 -0.178 0.858584
## STATUS_TTSewa -1.607e-01 3.695e-01 -0.435 0.663831
## USIA 5.312e-03 4.102e-03 1.295 0.195963
## STATUS_NIKAHMENIKAH 4.337e-01 2.371e-01 1.830 0.067904 .
## STATUS_NIKAHTIDAK MENIKAH 4.612e-01 2.411e-01 1.913 0.056334 .
## JML_TANGGUNGAN -1.155e-02 2.539e-02 -0.455 0.649214
## WAKTU_KREDIT 1.966e-02 8.028e-04 24.491 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4884 on 486 degrees of freedom
## Multiple R-squared: 0.7053, Adjusted R-squared: 0.6974
## F-statistic: 89.46 on 13 and 486 DF, p-value: < 2.2e-16
Misalkan terdapat 2 pelanggan baru yang berencana mengambil kredit. Berdasarkan profil pelanggan tersebut hitunglah perkiraan limit kredit yang dapat diberikan menggunakan 3 model yang telah dibangun sebelumnya.
# informasi calon pelanggan baru
data.baru <- data.frame(
PENGHASILAN = c(9000000, 15000000),
PEKERJAAN = c("Pegawai swasta", "Pegawai BUMN"),
LAMA_BEKERJA = c(14, 6),
STATUS_TT = c("Milik Sendiri", "Milik Sendiri"),
USIA = c(36, 25),
STATUS_NIKAH = c("MENIKAH", "TIDAK MENIKAH"),
JML_TANGGUNGAN = c(3, 0),
WAKTU_KREDIT = c(48, 36) # dalam bulan
)
print(t(data.baru))
## [,1] [,2]
## PENGHASILAN "9.0e+06" "1.5e+07"
## PEKERJAAN "Pegawai swasta" "Pegawai BUMN"
## LAMA_BEKERJA "14" " 6"
## STATUS_TT "Milik Sendiri" "Milik Sendiri"
## USIA "36" "25"
## STATUS_NIKAH "MENIKAH" "TIDAK MENIKAH"
## JML_TANGGUNGAN "3" "0"
## WAKTU_KREDIT "48" "36"
# Prediksi dengan Model 1 (Tanpa Transformasi)
limit.kredit.1 <- predict(mul.reg.1, data.baru)
cat(">>> BESAR LIMIT KREDIT (model 1):\n")
## >>> BESAR LIMIT KREDIT (model 1):
print(t(limit.kredit.1))
## 1 2
## [1,] 80899306 125322607
# Prediksi dengan Model 2 (Stepwise)
limit.kredit.2 <- predict(stepwise.reg, data.baru)
cat("\n\n>>> BESAR LIMIT KREDIT (model 2):\n")
##
##
## >>> BESAR LIMIT KREDIT (model 2):
print(t(limit.kredit.2))
## 1 2
## [1,] 85071604 131344055
# Prediksi dengan Model 3 (tRANSFORM Y)
limit.kredit.3 <- predict(mul.reg.log, data.baru)
# Karena nilai Y dalam bentuk Log natural,
# maka untuk memperoleh nilai dalam satuan asli
# perlu dilakukan inverse log natural (e^y)
cat("\n\n>>> BESAR LIMIT KREDIT (model 3):\n")
##
##
## >>> BESAR LIMIT KREDIT (model 3):
print(t(exp(limit.kredit.3)))
## 1 2
## [1,] 48281805 69854080
Menurut hasil pada model 1, orang pertama dapat diberikan kredit hingga 80,9 juta rupiah dan orang kedua sebesar 125,3 juta rupiah.
Sementara menurut model 2, memberikan limit yang sedikit lebih tinggi yaitu Rp 85,1 juta rupiah dan 131,3 juta rupiah.
Model 3 memberikan perkiraan limit kredit yang cukup rendah dibandingkan dua model sebelumnya, dimana orang pertama dapat mengambil kredit hingga 48,3 juta rupiah dan orang kedua sekitar 69,8 juta rupiah.
Model regresi logistik menggunakan peubah penjelas, baik kategorik atau kontinu, untuk memprediksi peluang dari hasil yang spesifik.
Regresi logistik dirancang untuk menggambarkan peluang yang terkait dengan nilai-nilai peubah respon.
Model logit (log odds) \(logit[\pi (x)] = log (\frac {\pi (x)}{1-\pi (x)}) = \alpha +\beta x\) sehingga, \(\pi(x) = \frac {exp (\alpha + \beta x)} {1+exp(\alpha + \beta x)}\). Dengan \(X\) sebagai peubah penjelas kuantitatif, \(Y\) peubah respons biner, \(\pi (x)\) peluang sukses peubah X.
Odds akan meningkat secara sebesar \(e^\beta\) untuk setiap kenaikan 1 unit \(X\)
\(e^\beta\) merupakan odds ratio (OR) yaitu \(OR = \frac {odds (X = x+1}{odds (X=x)}\)
An odds ratio indicates how much are likely, with respects to odds, a certain event occurs in one group relative to its occurrence in another group.
Nominal scales assign numbers to categories as labels with no ordering implied by the numbers.
Ordinal scales use numbers to indicate rank ordering on a single attribute.
To predict nominal outcomes using logistics regression, we use Multinomial Logistics Regression.
Multiple Logistics Regression
Also known as: Multi-class logistic regression or Maximum entropy classifier
This is a form of linear regression analysis conducted when the dependent variable is nominal with more than two levels.
t is used to describe data and to explain the relationship between one dependent nominal variable and one or more continuous-level (interval or ratio scale) independent variables.
The multinomial logistic regression estimates a separate binary logistic regression model for each dummy variables.
The result is M-1 binary logistic regression models.
Each model conveys the effect of predictors on the probability of success in that category, in comparison to the reference category.
Multinomial Logit Model \(log (\frac {\pi_j} {\pi_J}) = \alpha_j + \beta_j x,\ \ \ j = 1, ..., J-1\) where \(x\) denotes the explanotry variables, \(j\) denotes the level of nominal response variable, \(\alpha\) and \(\beta\) denotes the model coefficients. This is also known as Baseline-Category Logits. When \(J=2\) , this model simplifies to a single equation for \(log(\pi_1 / \pi_2) = logit(\pi_1)\) , resultiong in ordinary logistic regression for binary responses.
Limitation of Multinomial Logistics Regression
The category to which an outcome belongs to, does not assume any order in it.
However, in reality, we come across problems where categories have a natural order.
The category to which an outcome belongs to, does not assume any order in it. However, in reality, we come across problems where categories have a natural order. Use Ordinal Regression.
Ordinal Logistics Regression Model (OLM)
This is the most common model for ordinal outcomes.
The logit version also called: The proportional odds model or the cumulative logit model.
OLM also known as: parallel regression model or grouped continuous model.
Cummulative Logit Models for Ordinal Responses
A cumulative probability for Y is the probability that Y falls at or belowa particular point.
For outcome category j , the cumulative probability is:
\(P(Y\leq j) =\pi_1 + ... + \pi_j\), \(j = 1, …, J\) .
The logits of the cumulative probabilities are:
\(logit [P(Y\leq j)] = log [\frac{P(Y \leq j)} {1-P(Y \leq j)}] = log [ \frac{\pi_1 + ... + \pi_j} {\pi_{j+1} + ...+ \pi_j}]\), \(j = 1, …, J-1\) .
Cummulative Logit Models with Proportional Odds Property
A model for cumulative logit j looks like a binary logistic regression model in which categories 1–j combine to form a single category and categories j + 1 to J form a second category.
For an explanatory variable x, the model:
\(logit [P(Y\leq j)] =\alpha_j + \beta x\), \(j = 1, …, J-1\) .
Penyiapan Data
library(caret)
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(dplyr)
library(ggplot2)
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following object is masked from 'package:caret':
##
## predictors
## The following object is masked from 'package:car':
##
## logit
## The following object is masked from 'package:lmtest':
##
## lrtest
library(corrplot)
setwd("D:/Kuliah/Mat/TSA Kominfo/Praktikum")
data <- read.csv("4. obesity_dataset.csv")
head(data, 15)
## Gender Age Height Weight family_history_with_overweight FAVC FCVC NCP
## 1 Female 21 1.62 64.0 yes no 2 3
## 2 Female 21 1.52 56.0 yes no 3 3
## 3 Male 23 1.80 77.0 yes no 2 3
## 4 Male 27 1.80 87.0 no no 3 3
## 5 Male 22 1.78 89.8 no no 2 1
## 6 Male 29 1.62 53.0 no yes 2 3
## 7 Female 23 1.50 55.0 yes yes 3 3
## 8 Male 22 1.64 53.0 no no 2 3
## 9 Male 24 1.78 64.0 yes yes 3 3
## 10 Male 22 1.72 68.0 yes yes 2 3
## 11 Male 26 1.85 105.0 yes yes 3 3
## 12 Female 21 1.72 80.0 yes yes 2 3
## 13 Male 22 1.65 56.0 no no 3 3
## 14 Male 41 1.80 99.0 no yes 2 3
## 15 Male 23 1.77 60.0 yes yes 3 1
## CAEC SMOKE CH2O SCC FAF TUE CALC MTRANS
## 1 Sometimes no 2 no 0 1 no Public_Transportation
## 2 Sometimes yes 3 yes 3 0 Sometimes Public_Transportation
## 3 Sometimes no 2 no 2 1 Frequently Public_Transportation
## 4 Sometimes no 2 no 2 0 Frequently Walking
## 5 Sometimes no 2 no 0 0 Sometimes Public_Transportation
## 6 Sometimes no 2 no 0 0 Sometimes Automobile
## 7 Sometimes no 2 no 1 0 Sometimes Motorbike
## 8 Sometimes no 2 no 3 0 Sometimes Public_Transportation
## 9 Sometimes no 2 no 1 1 Frequently Public_Transportation
## 10 Sometimes no 2 no 1 1 no Public_Transportation
## 11 Frequently no 3 no 2 2 Sometimes Public_Transportation
## 12 Frequently no 2 yes 2 1 Sometimes Public_Transportation
## 13 Sometimes no 3 no 2 0 Sometimes Public_Transportation
## 14 Sometimes no 2 no 2 1 Frequently Automobile
## 15 Sometimes no 1 no 1 1 Sometimes Public_Transportation
## NObeyesdad
## 1 Normal_Weight
## 2 Normal_Weight
## 3 Normal_Weight
## 4 Overweight_Level_I
## 5 Overweight_Level_II
## 6 Normal_Weight
## 7 Normal_Weight
## 8 Normal_Weight
## 9 Normal_Weight
## 10 Normal_Weight
## 11 Obesity_Type_I
## 12 Overweight_Level_II
## 13 Normal_Weight
## 14 Obesity_Type_I
## 15 Normal_Weight
str(data)
## 'data.frame': 2111 obs. of 17 variables:
## $ Gender : chr "Female" "Female" "Male" "Male" ...
## $ Age : num 21 21 23 27 22 29 23 22 24 22 ...
## $ Height : num 1.62 1.52 1.8 1.8 1.78 1.62 1.5 1.64 1.78 1.72 ...
## $ Weight : num 64 56 77 87 89.8 53 55 53 64 68 ...
## $ family_history_with_overweight: chr "yes" "yes" "yes" "no" ...
## $ FAVC : chr "no" "no" "no" "no" ...
## $ FCVC : num 2 3 2 3 2 2 3 2 3 2 ...
## $ NCP : num 3 3 3 3 1 3 3 3 3 3 ...
## $ CAEC : chr "Sometimes" "Sometimes" "Sometimes" "Sometimes" ...
## $ SMOKE : chr "no" "yes" "no" "no" ...
## $ CH2O : num 2 3 2 2 2 2 2 2 2 2 ...
## $ SCC : chr "no" "yes" "no" "no" ...
## $ FAF : num 0 3 2 2 0 0 1 3 1 1 ...
## $ TUE : num 1 0 1 0 0 0 0 0 1 1 ...
## $ CALC : chr "no" "Sometimes" "Frequently" "Frequently" ...
## $ MTRANS : chr "Public_Transportation" "Public_Transportation" "Public_Transportation" "Walking" ...
## $ NObeyesdad : chr "Normal_Weight" "Normal_Weight" "Normal_Weight" "Overweight_Level_I" ...
Terdapat 7 Kelas pada dataset ini seperti yang tertera pada output kode berikut.
unique(data$NObeyesdad)
## [1] "Normal_Weight" "Overweight_Level_I" "Overweight_Level_II"
## [4] "Obesity_Type_I" "Insufficient_Weight" "Obesity_Type_II"
## [7] "Obesity_Type_III"
ob.levels <- c('Insufficient_Weight', 'Normal_Weight', 'Overweight_Level_I', 'Overweight_Level_II',
'Obesity_Type_I', 'Obesity_Type_II', 'Obesity_Type_III')
ob.levels.s <- c("Insufficient", "Normal", "Over_Lv1", "Over_Lv2", "Obes_Lv1", "Obes_Lv2", "Obes_Lv3")
Untuk alasan merapikan dan mempersingkat nama pada masing-masing kategori maka akan dilakukan perubahan nama pada setiap kategorinya.
ubah_nama_level <- function(level) {
# Cari indeks teks dalam vektor
indeks <- which(ob.levels == level) #["NOrmal Weight",.....]
#indeks = 2
return (ob.levels.s[indeks])
}
# merubah nama setiap kategori berdasarkan on.leves.s
data$NObeyesdad <- sapply(data$NObeyesdad, ubah_nama_level)
# Merubah features kategorik menggunakan label encoding
data$FAVC <- as.integer(as.factor(data$FAVC))
data$CAEC <- as.integer(as.factor(data$CAEC))
data$CALC <- as.integer(as.factor(data$CALC))
data$MTRANS <- as.integer(as.factor(data$MTRANS))
data$Gender <- factor(data$Gender)
data$family_history_with_overweight <- factor(data$family_history_with_overweight)
data$SMOKE <- factor(data$SMOKE)
data$SCC <- factor(data$SCC)
# Mengubah peubah Y menjadi factor dengan urutan (level)
data$NObeyesdad <- factor(data$NObeyesdad, levels=ob.levels.s)
head(data)
## Gender Age Height Weight family_history_with_overweight FAVC FCVC NCP CAEC
## 1 Female 21 1.62 64.0 yes 1 2 3 4
## 2 Female 21 1.52 56.0 yes 1 3 3 4
## 3 Male 23 1.80 77.0 yes 1 2 3 4
## 4 Male 27 1.80 87.0 no 1 3 3 4
## 5 Male 22 1.78 89.8 no 1 2 1 4
## 6 Male 29 1.62 53.0 no 2 2 3 4
## SMOKE CH2O SCC FAF TUE CALC MTRANS NObeyesdad
## 1 no 2 no 0 1 3 4 Normal
## 2 yes 3 yes 3 0 4 4 Normal
## 3 no 2 no 2 1 2 4 Normal
## 4 no 2 no 2 0 2 5 Over_Lv1
## 5 no 2 no 0 0 4 4 Over_Lv2
## 6 no 2 no 0 0 4 1 Normal
Secara umum dataset yang dimiliki memiliki jumlah yang relatif balance untuk setiap kategorinya. Sehingga tidak perlu melakukan penanganan untuk menyeimbangkan data.
# Jumlah observasi untuk masing-masing kelas
data.frame(table(data$NObeyesdad))
## Var1 Freq
## 1 Insufficient 272
## 2 Normal 287
## 3 Over_Lv1 290
## 4 Over_Lv2 290
## 5 Obes_Lv1 351
## 6 Obes_Lv2 297
## 7 Obes_Lv3 324
Statistika Deskriptif
count.by.gender <- table(data$NObeyesdad, data$Gender)
data_frame(count.by.gender)
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## ℹ Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## # A tibble: 7 × 1
## count.by.gender[,"Female"] [,"Male"]
## <int> <int>
## 1 173 99
## 2 141 146
## 3 145 145
## 4 103 187
## 5 156 195
## 6 2 295
## 7 323 1
count.by.famhist <- table(data$NObeyesdad, data$family_history_with_overweight)
data.frame(count.by.famhist)
## Var1 Var2 Freq
## 1 Insufficient no 146
## 2 Normal no 132
## 3 Over_Lv1 no 81
## 4 Over_Lv2 no 18
## 5 Obes_Lv1 no 7
## 6 Obes_Lv2 no 1
## 7 Obes_Lv3 no 0
## 8 Insufficient yes 126
## 9 Normal yes 155
## 10 Over_Lv1 yes 209
## 11 Over_Lv2 yes 272
## 12 Obes_Lv1 yes 344
## 13 Obes_Lv2 yes 296
## 14 Obes_Lv3 yes 324
count.by.smoke <- table(data$NObeyesdad, data$SMOKE)
data.frame(count.by.smoke)
## Var1 Var2 Freq
## 1 Insufficient no 271
## 2 Normal no 274
## 3 Over_Lv1 no 287
## 4 Over_Lv2 no 285
## 5 Obes_Lv1 no 345
## 6 Obes_Lv2 no 282
## 7 Obes_Lv3 no 323
## 8 Insufficient yes 1
## 9 Normal yes 13
## 10 Over_Lv1 yes 3
## 11 Over_Lv2 yes 5
## 12 Obes_Lv1 yes 6
## 13 Obes_Lv2 yes 15
## 14 Obes_Lv3 yes 1
count.by.scc <- table(data$NObeyesdad, data$SCC)
data.frame(count.by.scc)
## Var1 Var2 Freq
## 1 Insufficient no 250
## 2 Normal no 257
## 3 Over_Lv1 no 253
## 4 Over_Lv2 no 286
## 5 Obes_Lv1 no 349
## 6 Obes_Lv2 no 296
## 7 Obes_Lv3 no 324
## 8 Insufficient yes 22
## 9 Normal yes 30
## 10 Over_Lv1 yes 37
## 11 Over_Lv2 yes 4
## 12 Obes_Lv1 yes 2
## 13 Obes_Lv2 yes 1
## 14 Obes_Lv3 yes 0
Pembagian Data
Tahapan selanjutnya adalah pembagian dataset menjadi 2 bagian yaitu data latih dan data uji. Data latih akan digunakan untuk membangun model. Sementara data uji akan digunakan untuk evaluasi dan pengukuran performa model pada data baru.
# Membagi data Latih dan Data Uji
library(caret)
set.seed(123) # Untuk hasil yang dapat direproduksi
# Fungsi createDataPartition dari paket caret untuk membagi data
# dengan menghasilkan indeks data latih dan data uji
trainIndex <- createDataPartition(data$NObeyesdad, p = 0.75,
list = FALSE,
times = 1)
# Buat data latih dan data uji berdasarkan indeks yang dihasilkan
data.train <- data[trainIndex, ]
data.test <- data[-trainIndex, ]
# melihat keterwakilan masing-masing kategori
# pada data latih maupun data uji
cbind("Train"=table(data.train$NObeyesdad), "Test"=table(data.test$NObeyesdad))
## Train Test
## Insufficient 204 68
## Normal 216 71
## Over_Lv1 218 72
## Over_Lv2 218 72
## Obes_Lv1 264 87
## Obes_Lv2 223 74
## Obes_Lv3 243 81
Pembuatan Model
# pengaturan agar informasi warning tidak tercetak
options(warn = -1)
library(VGAM)
model.1 <- vglm(NObeyesdad ~ ., data=data.train, family="propodds")
summary(model.1)
##
## Call:
## vglm(formula = NObeyesdad ~ ., family = "propodds", data = data.train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 34.241026 2.058222 16.636 < 2e-16 ***
## (Intercept):2 29.334225 1.977170 14.836 < 2e-16 ***
## (Intercept):3 26.112495 1.950249 13.389 < 2e-16 ***
## (Intercept):4 22.749482 1.899227 11.978 < 2e-16 ***
## (Intercept):5 17.638708 1.850752 9.531 < 2e-16 ***
## (Intercept):6 13.115666 1.836558 7.141 9.24e-13 ***
## GenderMale -0.866568 0.179717 -4.822 1.42e-06 ***
## Age 0.009587 0.013609 0.704 0.481143
## Height -33.732558 1.439910 -23.427 < 2e-16 ***
## Weight 0.377887 0.012520 30.183 < 2e-16 ***
## family_history_with_overweightyes 0.548120 0.194446 2.819 0.004819 **
## FAVC 0.065387 0.193491 0.338 0.735415
## FCVC 0.070794 0.125744 0.563 0.573437
## NCP 0.161940 0.078754 2.056 0.039755 *
## CAEC 0.332652 0.090936 3.658 0.000254 ***
## SMOKEyes -0.357705 0.407013 -0.879 0.379481
## CH2O -0.086401 0.103896 -0.832 0.405628
## SCCyes -0.099355 0.288968 -0.344 0.730975
## FAF -0.058516 0.079157 -0.739 0.459761
## TUE -0.119395 0.102124 -1.169 0.242353
## CALC -0.240213 0.112360 -2.138 0.032525 *
## MTRANS 0.166312 0.062519 2.660 0.007810 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 6
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]),
## logitlink(P[Y>=4]), logitlink(P[Y>=5]), logitlink(P[Y>=6]), logitlink(P[Y>=7])
##
## Residual deviance: 1191.529 on 9494 degrees of freedom
##
## Log-likelihood: NA on 9494 degrees of freedom
##
## Number of Fisher scoring iterations: 15
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## 'Height', 'Weight'
##
##
## Exponentiated coefficients:
## GenderMale Age
## 4.203921e-01 1.009633e+00
## Height Weight
## 2.239423e-15 1.459198e+00
## family_history_with_overweightyes FAVC
## 1.729998e+00 1.067572e+00
## FCVC NCP
## 1.073360e+00 1.175790e+00
## CAEC SMOKEyes
## 1.394661e+00 6.992794e-01
## CH2O SCCyes
## 9.172260e-01 9.054209e-01
## FAF TUE
## 9.431633e-01 8.874568e-01
## CALC MTRANS
## 7.864601e-01 1.180941e+00
Prediksi (Data Train)
Hasil prediksi pada model regresi logistik ordinal adalah besarnya peluang data tersebut masuk ke setiap kelas atau kategori. Penentuan data tersebut masuk ke kelas yang mana yaitu dengan memilih kelas dengan nilai peluang terbesar.
# Memprediksi data TRAIN
pred.train <- predict(model.1, newdata = data.train, type="response")
round(head(pred.train), 4)
## Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## 2 0.0031 0.2904 0.6189 0.0843 0.0033 0e+00 0
## 3 0.0131 0.6293 0.3358 0.0209 0.0008 0e+00 0
## 4 0.0004 0.0453 0.4997 0.4265 0.0279 2e-04 0
## 6 0.4613 0.5302 0.0082 0.0003 0.0000 0e+00 0
## 7 0.0013 0.1462 0.6652 0.1794 0.0079 0e+00 0
## 8 0.5814 0.4133 0.0051 0.0002 0.0000 0e+00 0
# merubah hasil prediksi (peluang) menjadi kategori
# berdasarkan kategori dengan peluang terbesar
pred.train.cat <- factor(ob.levels.s[max.col(pred.train)])
# prediksi kategori untuk 6 data pertama (train)
head(pred.train.cat)
## [1] Over_Lv1 Normal Over_Lv1 Normal Over_Lv1
## [6] Insufficient
## Levels: Insufficient Normal Obes_Lv1 Obes_Lv2 Obes_Lv3 Over_Lv1 Over_Lv2
# membuat confusion Matrix hasil data train
cm.train <- confusionMatrix(data = pred.train.cat,
reference = data.train$NObeyesdad)
# performa model secara keseluruhan (data train)
round(data.frame(cm.train$overall), 4)
## cm.train.overall
## Accuracy 0.9079
## Kappa 0.8925
## AccuracyLower 0.8926
## AccuracyUpper 0.9217
## AccuracyNull 0.1665
## AccuracyPValue 0.0000
## McnemarPValue NaN
# performa model untuk setiap kategori (data train)
# transpose : (agar tabel tidak terlalu lebar)
round(t(cm.train$byClass), 4)
## Class: Insufficient Class: Normal Class: Over_Lv1
## Sensitivity 0.9853 0.7500 0.9404
## Specificity 0.9797 0.9956 0.9642
## Pos Pred Value 0.8777 0.9643 0.8071
## Neg Pred Value 0.9978 0.9619 0.9902
## Precision 0.8777 0.9643 0.8071
## Recall 0.9853 0.7500 0.9404
## F1 0.9284 0.8437 0.8686
## Prevalence 0.1286 0.1362 0.1375
## Detection Rate 0.1267 0.1021 0.1293
## Detection Prevalence 0.1444 0.1059 0.1602
## Balanced Accuracy 0.9825 0.8728 0.9523
## Class: Over_Lv2 Class: Obes_Lv1 Class: Obes_Lv2
## Sensitivity 0.8853 0.9659 0.9193
## Specificity 0.9905 0.9962 0.9780
## Pos Pred Value 0.9369 0.9808 0.8723
## Neg Pred Value 0.9819 0.9932 0.9867
## Precision 0.9369 0.9808 0.8723
## Recall 0.8853 0.9659 0.9193
## F1 0.9104 0.9733 0.8952
## Prevalence 0.1375 0.1665 0.1406
## Detection Rate 0.1217 0.1608 0.1293
## Detection Prevalence 0.1299 0.1639 0.1482
## Balanced Accuracy 0.9379 0.9811 0.9486
## Class: Obes_Lv3
## Sensitivity 0.9012
## Specificity 0.9888
## Pos Pred Value 0.9359
## Neg Pred Value 0.9822
## Precision 0.9359
## Recall 0.9012
## F1 0.9182
## Prevalence 0.1532
## Detection Rate 0.1381
## Detection Prevalence 0.1475
## Balanced Accuracy 0.9450
PREDIKSI & EVALUASI (DATA TEST)
Hasil yang sudah ditampilkan sebelumnya merupakan informasi mengenai performa model untuk data latih. Atau dengan kata lain menunjukkan pefroma model pada data yang digunakan untuk membangun model tersebut.
Untuk mengetahui performa dengan lebih fair, dapat dilakukan menggunakan data lainnya yang belum pernah dilihat model selama pelatihan. Dalam hal ini adalah data uji.
# Memprediksi data TEST
pred.test <- predict(model.1, newdata = data.test, type="response")
round(head(pred.test), 4)
## Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## 1 0.0020 0.2126 0.6580 0.1224 0.0050 0e+00 0
## 5 0.0002 0.0217 0.3375 0.5825 0.0578 4e-04 0
## 17 0.0001 0.0172 0.2889 0.6210 0.0723 5e-04 0
## 26 0.6091 0.3862 0.0045 0.0002 0.0000 0e+00 0
## 27 0.0716 0.8409 0.0837 0.0037 0.0001 0e+00 0
## 29 0.0061 0.4484 0.4998 0.0440 0.0016 0e+00 0
# merubah hasil prediksi (peluang) menjadi kategori
# berdasarkan kategori dengan peluang terbesar
pred.test.cat <- factor(ob.levels.s[max.col(pred.test)])
# prediksi kategori untuk 6 data pertama (test)
head(pred.test.cat)
## [1] Over_Lv1 Over_Lv2 Over_Lv2 Insufficient Normal
## [6] Over_Lv1
## Levels: Insufficient Normal Obes_Lv1 Obes_Lv2 Obes_Lv3 Over_Lv1 Over_Lv2
# membuat confusion Matrix (data test)
cm.test <- confusionMatrix(data = pred.test.cat,
reference = data.test$NObeyesdad)
# performa model secara keseluruhan (data test)
round(data.frame(cm.test$overall), 4)
## cm.test.overall
## Accuracy 0.8990
## Kappa 0.8821
## AccuracyLower 0.8700
## AccuracyUpper 0.9235
## AccuracyNull 0.1657
## AccuracyPValue 0.0000
## McnemarPValue NaN
# Performa model pada setiap kategori (data test)
# transpose : (agar tabel tidak terlalu lebar)
round(t(cm.test$byClass), 4)
## Class: Insufficient Class: Normal Class: Over_Lv1
## Sensitivity 1.0000 0.7465 0.9028
## Specificity 0.9781 0.9956 0.9603
## Pos Pred Value 0.8718 0.9636 0.7831
## Neg Pred Value 1.0000 0.9617 0.9842
## Precision 0.8718 0.9636 0.7831
## Recall 1.0000 0.7465 0.9028
## F1 0.9315 0.8413 0.8387
## Prevalence 0.1295 0.1352 0.1371
## Detection Rate 0.1295 0.1010 0.1238
## Detection Prevalence 0.1486 0.1048 0.1581
## Balanced Accuracy 0.9891 0.8710 0.9315
## Class: Over_Lv2 Class: Obes_Lv1 Class: Obes_Lv2
## Sensitivity 0.8333 0.9885 0.8919
## Specificity 0.9890 0.9909 0.9823
## Pos Pred Value 0.9231 0.9556 0.8919
## Neg Pred Value 0.9739 0.9977 0.9823
## Precision 0.9231 0.9556 0.8919
## Recall 0.8333 0.9885 0.8919
## F1 0.8759 0.9718 0.8919
## Prevalence 0.1371 0.1657 0.1410
## Detection Rate 0.1143 0.1638 0.1257
## Detection Prevalence 0.1238 0.1714 0.1410
## Balanced Accuracy 0.9111 0.9897 0.9371
## Class: Obes_Lv3
## Sensitivity 0.9136
## Specificity 0.9865
## Pos Pred Value 0.9250
## Neg Pred Value 0.9843
## Precision 0.9250
## Recall 0.9136
## F1 0.9193
## Prevalence 0.1543
## Detection Rate 0.1410
## Detection Prevalence 0.1524
## Balanced Accuracy 0.9500
# confusion matrix (data test)
cm.test$table
## Reference
## Prediction Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## Insufficient 68 10 0 0 0 0 0
## Normal 0 53 2 0 0 0 0
## Over_Lv1 0 8 65 10 0 0 0
## Over_Lv2 0 0 5 60 0 0 0
## Obes_Lv1 0 0 0 2 86 2 0
## Obes_Lv2 0 0 0 0 1 66 7
## Obes_Lv3 0 0 0 0 0 6 74
Regresi Logistik Multinomial (Multinomial Logistic Regression) adalah salah satu metode analisis statistik yang digunakan untuk memodelkan hubungan antara variabel independen (prediktor) dan variabel dependen yang terdiri dari lebih dari dua kategori atau kelas. Ini adalah perluasan bentuk dari regresi logistik binomial, yang hanya digunakan untuk memodelkan variabel dependen dengan dua kategori (biner) saja. Dalam regresi logistik multinomial, variabel dependen memiliki tiga atau lebih kategori yang tidak berurutan (tidak memiliki tingkatan atau tidak memntingkan tingkatannya).
Output dari model regresi logistik multinomial adalah kumpulan regresi logistik biner yaitu sebanyak K-1 Model. (K: banyaknya kategori)
Note :
Untuk menyederhanakan contoh pada bagian ini, mari kita anggap saja data
sebelumnya yaitu data obesitas sebagai kategori multinomial. Dimana
tidak terdapat urutan atau tingkatan pada setiap kategori. Kita hanya
tertarik untuk mengetahui seseorang masuk ke dalam kategori apa
berdasarkan informasi-informasi yang tersedia.
Pembuatan Model
model.mul <- vglm(NObeyesdad ~ ., family=multinomial(refLevel="Normal"), data=data.train)
summary(model.mul)
##
## Call:
## vglm(formula = NObeyesdad ~ ., family = multinomial(refLevel = "Normal"),
## data = data.train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.026e+02 2.615e+01 -3.924 8.72e-05 ***
## (Intercept):2 1.087e+02 2.215e+01 4.909 9.15e-07 ***
## (Intercept):3 2.325e+02 2.960e+01 7.853 4.08e-15 ***
## (Intercept):4 3.480e+02 4.032e+01 8.630 < 2e-16 ***
## (Intercept):5 3.839e+02 5.750e+01 6.678 2.43e-11 ***
## (Intercept):6 3.489e+02 1.552e+02 NA NA
## GenderMale:1 -3.681e-01 1.884e+00 -0.195 0.84509
## GenderMale:2 -2.865e+00 1.336e+00 -2.145 0.03196 *
## GenderMale:3 -2.276e+00 1.589e+00 -1.432 0.15200
## GenderMale:4 -4.974e+00 2.591e+00 -1.920 0.05489 .
## GenderMale:5 4.158e+00 6.922e+00 0.601 0.54806
## GenderMale:6 -1.956e+01 1.604e+01 NA NA
## Age:1 5.338e-02 1.798e-01 0.297 0.76653
## Age:2 7.865e-02 7.906e-02 0.995 0.31981
## Age:3 2.185e-03 9.586e-02 0.023 0.98182
## Age:4 7.833e-02 1.686e-01 0.465 0.64226
## Age:5 8.286e-01 3.910e-01 2.119 0.03407 *
## Age:6 -4.231e-02 1.340e+00 -0.032 0.97481
## Height:1 1.238e+02 2.733e+01 4.529 5.93e-06 ***
## Height:2 -1.433e+02 2.493e+01 -5.748 9.03e-09 ***
## Height:3 -2.855e+02 3.218e+01 -8.872 < 2e-16 ***
## Height:4 -4.371e+02 4.343e+01 -10.063 < 2e-16 ***
## Height:5 -5.517e+02 5.403e+01 -10.211 < 2e-16 ***
## Height:6 -5.075e+02 9.769e+01 -5.195 2.05e-07 ***
## Weight:1 -2.029e+00 4.083e-01 -4.969 6.72e-07 ***
## Weight:2 1.807e+00 3.056e-01 5.913 3.35e-09 ***
## Weight:3 3.463e+00 3.864e-01 8.963 < 2e-16 ***
## Weight:4 5.076e+00 4.871e-01 10.420 < 2e-16 ***
## Weight:5 6.392e+00 5.622e-01 11.369 < 2e-16 ***
## Weight:6 6.072e+00 5.950e-01 10.205 < 2e-16 ***
## family_history_with_overweightyes:1 1.552e+00 1.405e+00 1.105 0.26931
## family_history_with_overweightyes:2 -3.016e-01 9.262e-01 -0.326 0.74467
## family_history_with_overweightyes:3 2.839e+00 1.312e+00 2.163 0.03056 *
## family_history_with_overweightyes:4 3.569e+00 2.315e+00 1.542 0.12318
## family_history_with_overweightyes:5 2.130e-01 1.278e+01 0.017 0.98671
## family_history_with_overweightyes:6 -1.038e+01 1.699e+01 -0.611 0.54128
## FAVC:1 -1.446e+00 1.994e+00 -0.725 0.46847
## FAVC:2 4.954e-01 1.474e+00 0.336 0.73673
## FAVC:3 -2.330e+00 1.724e+00 -1.352 0.17647
## FAVC:4 -1.578e+00 2.241e+00 -0.704 0.48135
## FAVC:5 7.797e-02 1.221e+01 0.006 0.99491
## FAVC:6 6.383e+00 2.366e+01 0.270 0.78736
## FCVC:1 1.301e+00 1.135e+00 1.146 0.25172
## FCVC:2 -1.101e+00 7.874e-01 -1.398 0.16197
## FCVC:3 -1.565e+00 1.108e+00 -1.412 0.15788
## FCVC:4 -1.645e+00 1.538e+00 -1.069 0.28504
## FCVC:5 -2.927e+00 2.642e+00 -1.108 0.26802
## FCVC:6 3.312e-01 1.081e+01 0.031 0.97556
## NCP:1 4.406e-01 7.342e-01 0.600 0.54840
## NCP:2 3.539e-01 5.305e-01 0.667 0.50475
## NCP:3 -4.491e-02 6.620e-01 -0.068 0.94592
## NCP:4 -1.498e-01 9.697e-01 -0.154 0.87726
## NCP:5 -1.059e+00 1.469e+00 -0.721 0.47105
## NCP:6 5.239e+00 7.336e+00 0.714 0.47509
## CAEC:1 -4.636e-01 5.434e-01 -0.853 0.39366
## CAEC:2 1.399e+00 5.298e-01 2.642 0.00825 **
## CAEC:3 1.626e+00 7.933e-01 2.050 0.04039 *
## CAEC:4 1.771e+00 1.063e+00 1.666 0.09563 .
## CAEC:5 3.127e+00 1.950e+00 1.604 0.10882
## CAEC:6 1.562e+00 7.355e+00 0.212 0.83185
## SMOKEyes:1 -5.386e+00 2.840e+01 NA NA
## SMOKEyes:2 -4.073e+00 2.818e+00 -1.445 0.14843
## SMOKEyes:3 -6.446e+00 3.737e+00 -1.725 0.08456 .
## SMOKEyes:4 -4.605e+00 5.479e+00 -0.841 0.40062
## SMOKEyes:5 -7.996e+00 2.213e+01 -0.361 0.71783
## SMOKEyes:6 -9.580e-01 3.623e+01 -0.026 0.97890
## CH2O:1 8.150e-01 1.003e+00 0.813 0.41628
## CH2O:2 -2.513e-01 7.833e-01 -0.321 0.74837
## CH2O:3 5.749e-02 9.471e-01 0.061 0.95160
## CH2O:4 -2.570e-01 1.281e+00 -0.201 0.84096
## CH2O:5 -1.728e+00 2.319e+00 -0.745 0.45618
## CH2O:6 4.560e-02 5.828e+00 0.008 0.99376
## SCCyes:1 -6.156e-01 2.219e+00 -0.277 0.78146
## SCCyes:2 4.981e+00 1.685e+00 2.956 0.00312 **
## SCCyes:3 5.965e+00 2.201e+00 2.710 0.00673 **
## SCCyes:4 8.138e+00 9.355e+00 0.870 0.38437
## SCCyes:5 7.126e+00 2.303e+01 0.309 0.75696
## SCCyes:6 8.502e+00 2.136e+01 0.398 0.69061
## FAF:1 -8.439e-01 7.373e-01 -1.145 0.25238
## FAF:2 -3.346e-01 4.737e-01 -0.706 0.47994
## FAF:3 -8.778e-01 6.212e-01 -1.413 0.15765
## FAF:4 -1.159e+00 1.018e+00 -1.138 0.25508
## FAF:5 -4.093e+00 1.868e+00 -2.191 0.02847 *
## FAF:6 -2.750e+00 6.432e+00 -0.428 0.66892
## TUE:1 6.670e-02 9.662e-01 0.069 0.94496
## TUE:2 2.731e-01 6.320e-01 0.432 0.66570
## TUE:3 1.236e+00 8.369e-01 1.477 0.13966
## TUE:4 8.708e-01 1.217e+00 0.716 0.47416
## TUE:5 -4.333e-01 1.791e+00 -0.242 0.80882
## TUE:6 -1.485e+00 8.178e+00 -0.182 0.85590
## CALC:1 9.498e-02 9.622e-01 0.099 0.92137
## CALC:2 3.622e-01 6.298e-01 0.575 0.56520
## CALC:3 -1.320e+00 7.770e-01 -1.699 0.08939 .
## CALC:4 -1.661e+00 1.171e+00 -1.418 0.15610
## CALC:5 -2.501e+00 3.005e+00 -0.832 0.40524
## CALC:6 -9.624e-01 1.201e+01 NA NA
## MTRANS:1 -1.712e-01 4.935e-01 -0.347 0.72864
## MTRANS:2 -1.822e-02 4.284e-01 -0.043 0.96608
## MTRANS:3 -2.735e-01 5.184e-01 -0.528 0.59774
## MTRANS:4 2.292e-01 9.978e-01 0.230 0.81835
## MTRANS:5 1.970e+00 1.405e+00 1.402 0.16085
## MTRANS:6 5.631e-01 6.361e+00 0.089 0.92946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 6
##
## Names of linear predictors: log(mu[,1]/mu[,2]), log(mu[,3]/mu[,2]),
## log(mu[,4]/mu[,2]), log(mu[,5]/mu[,2]), log(mu[,6]/mu[,2]), log(mu[,7]/mu[,2])
##
## Residual deviance: 132.6246 on 9414 degrees of freedom
##
## Log-likelihood: -66.3123 on 9414 degrees of freedom
##
## Number of Fisher scoring iterations: 22
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):2', '(Intercept):5', '(Intercept):6', 'GenderMale:6', 'Height:1', 'Height:3', 'Height:4', 'Height:5', 'Weight:2', 'SMOKEyes:1', 'CALC:6'
##
##
## Reference group is level 2 of the response
Prediksi (Data Latih)
# Memprediksi data TRAIN
pred.train.mul <- predict(model.mul, newdata = data.train, type="response")
round(head(pred.train.mul), 4)
## Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## 2 0e+00 0.8972 0.1028 0.00 0 0 0
## 3 0e+00 0.9988 0.0012 0.00 0 0 0
## 4 0e+00 0.0000 0.3200 0.68 0 0 0
## 6 6e-04 0.9994 0.0000 0.00 0 0 0
## 7 0e+00 0.6045 0.3955 0.00 0 0 0
## 8 1e-03 0.9990 0.0000 0.00 0 0 0
pred.train.mul.cat <- factor(ob.levels.s[max.col(pred.train.mul)])
head(pred.train.mul.cat)
## [1] Normal Normal Over_Lv2 Normal Normal Normal
## Levels: Insufficient Normal Obes_Lv1 Obes_Lv2 Obes_Lv3 Over_Lv1 Over_Lv2
# membuat confusion Matrix (data train)
cm.mul.train <- confusionMatrix(data = pred.train.mul.cat,
reference = data.train$NObeyesdad)
# performa model secara keseluruhan (data test)
round(data.frame(cm.mul.train$overall), 4)
## cm.mul.train.overall
## Accuracy 0.9887
## Kappa 0.9867
## AccuracyLower 0.9821
## AccuracyUpper 0.9933
## AccuracyNull 0.1665
## AccuracyPValue 0.0000
## McnemarPValue NaN
# Performa model pada setiap kategori (data train)
# transpose : (agar tabel tidak terlalu lebar)
round(t(cm.mul.train$byClass), 4)
## Class: Insufficient Class: Normal Class: Over_Lv1
## Sensitivity 1.0000 0.9861 0.9633
## Specificity 0.9986 0.9993 0.9942
## Pos Pred Value 0.9903 0.9953 0.9633
## Neg Pred Value 1.0000 0.9978 0.9942
## Precision 0.9903 0.9953 0.9633
## Recall 1.0000 0.9861 0.9633
## F1 0.9951 0.9907 0.9633
## Prevalence 0.1286 0.1362 0.1375
## Detection Rate 0.1286 0.1343 0.1324
## Detection Prevalence 0.1299 0.1349 0.1375
## Balanced Accuracy 0.9993 0.9927 0.9787
## Class: Over_Lv2 Class: Obes_Lv1 Class: Obes_Lv2
## Sensitivity 0.9679 1.0000 1.0000
## Specificity 0.9949 1.0000 1.0000
## Pos Pred Value 0.9679 1.0000 1.0000
## Neg Pred Value 0.9949 1.0000 1.0000
## Precision 0.9679 1.0000 1.0000
## Recall 0.9679 1.0000 1.0000
## F1 0.9679 1.0000 1.0000
## Prevalence 0.1375 0.1665 0.1406
## Detection Rate 0.1330 0.1665 0.1406
## Detection Prevalence 0.1375 0.1665 0.1406
## Balanced Accuracy 0.9814 1.0000 1.0000
## Class: Obes_Lv3
## Sensitivity 1.0000
## Specificity 1.0000
## Pos Pred Value 1.0000
## Neg Pred Value 1.0000
## Precision 1.0000
## Recall 1.0000
## F1 1.0000
## Prevalence 0.1532
## Detection Rate 0.1532
## Detection Prevalence 0.1532
## Balanced Accuracy 1.0000
# confusion matrix (data test)
cm.mul.train$table
## Reference
## Prediction Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## Insufficient 204 2 0 0 0 0 0
## Normal 0 213 1 0 0 0 0
## Over_Lv1 0 1 210 7 0 0 0
## Over_Lv2 0 0 7 211 0 0 0
## Obes_Lv1 0 0 0 0 264 0 0
## Obes_Lv2 0 0 0 0 0 223 0
## Obes_Lv3 0 0 0 0 0 0 243
Prediksi & Evaluasi (Data Test)
# Memprediksi data TEST
pred.test.mul <- predict(model.mul, newdata = data.test, type="response")
round(head(pred.test.mul), 4)
## Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## 1 0e+00 0.6638 0.3361 0.0002 0e+00 0 0
## 5 0e+00 0.0000 0.0004 0.9995 1e-04 0 0
## 17 0e+00 0.0000 0.0030 0.9970 0e+00 0 0
## 26 3e-04 0.9997 0.0000 0.0000 0e+00 0 0
## 27 0e+00 1.0000 0.0000 0.0000 0e+00 0 0
## 29 0e+00 0.0436 0.9563 0.0000 0e+00 0 0
pred.test.mul.cat <- factor(ob.levels.s[max.col(pred.test.mul)])
head(pred.test.mul.cat)
## [1] Normal Over_Lv2 Over_Lv2 Normal Normal Over_Lv1
## Levels: Insufficient Normal Obes_Lv1 Obes_Lv2 Obes_Lv3 Over_Lv1 Over_Lv2
# membuat confusion Matrix (data test)
cm.mul.test <- confusionMatrix(data = pred.test.mul.cat,
reference = data.test$NObeyesdad)
# performa model secara keseluruhan (data test)
round(data.frame(cm.mul.test$overall), 4)
## cm.mul.test.overall
## Accuracy 0.9581
## Kappa 0.9511
## AccuracyLower 0.9372
## AccuracyUpper 0.9736
## AccuracyNull 0.1657
## AccuracyPValue 0.0000
## McnemarPValue NaN
# Performa model pada setiap kategori (data test)
# transpose : (agar tabel tidak terlalu lebar)
round(t(cm.mul.test$byClass), 4)
## Class: Insufficient Class: Normal Class: Over_Lv1
## Sensitivity 0.9853 0.8873 0.9444
## Specificity 0.9869 0.9956 0.9868
## Pos Pred Value 0.9178 0.9692 0.9189
## Neg Pred Value 0.9978 0.9826 0.9911
## Precision 0.9178 0.9692 0.9189
## Recall 0.9853 0.8873 0.9444
## F1 0.9504 0.9265 0.9315
## Prevalence 0.1295 0.1352 0.1371
## Detection Rate 0.1276 0.1200 0.1295
## Detection Prevalence 0.1390 0.1238 0.1410
## Balanced Accuracy 0.9861 0.9415 0.9656
## Class: Over_Lv2 Class: Obes_Lv1 Class: Obes_Lv2
## Sensitivity 0.9167 0.9885 0.9865
## Specificity 0.9934 0.9954 0.9956
## Pos Pred Value 0.9565 0.9773 0.9733
## Neg Pred Value 0.9868 0.9977 0.9978
## Precision 0.9565 0.9773 0.9733
## Recall 0.9167 0.9885 0.9865
## F1 0.9362 0.9829 0.9799
## Prevalence 0.1371 0.1657 0.1410
## Detection Rate 0.1257 0.1638 0.1390
## Detection Prevalence 0.1314 0.1676 0.1429
## Balanced Accuracy 0.9550 0.9920 0.9910
## Class: Obes_Lv3
## Sensitivity 0.9877
## Specificity 0.9977
## Pos Pred Value 0.9877
## Neg Pred Value 0.9977
## Precision 0.9877
## Recall 0.9877
## F1 0.9877
## Prevalence 0.1543
## Detection Rate 0.1524
## Detection Prevalence 0.1543
## Balanced Accuracy 0.9927
# confusion matrix (data test)
cm.mul.test$table
## Reference
## Prediction Insufficient Normal Over_Lv1 Over_Lv2 Obes_Lv1 Obes_Lv2 Obes_Lv3
## Insufficient 67 6 0 0 0 0 0
## Normal 1 63 1 0 0 0 0
## Over_Lv1 0 2 68 4 0 0 0
## Over_Lv2 0 0 3 66 0 0 0
## Obes_Lv1 0 0 0 2 86 0 0
## Obes_Lv2 0 0 0 0 1 73 1
## Obes_Lv3 0 0 0 0 0 1 80