Introduction

Dokumen ini mengeksplorasi kumpulan data yang berisi harga dan atribut untuk sekitar 54.000 berlian bulat. Ada 53.940 berlian dalam dataset dengan 10 fitur (karat, potongan, warna, kejelasan, kedalaman, tabel, harga, x, y, dan z). Sebagian besar variabel bersifat numerik, tetapi variabel potongan(cut), warna(color), dan kejernihan(clarity) merupakan variabel faktor terurut dengan tingkatan di dalamnya. Tentang mata uang untuk kolom harga(price) di dataset ini menggunakan satuan dollar($). Dan Tentang kolom x,y, dan z adalah pengukuran berlian sebagai (( x: panjang dalam mm, y: lebar dalam mm,z: kedalaman dalam mm ))

Dalam pembahasan kali ini, kita akan menggunakan dataset yang telah disediakan oleh Kaggle yang bernama “Diamond dataset” .

Objektif:

Menganalisis faktor-faktor apa saja yang memengaruhi harga berlian.
Membandingkan beberapa model regresi linear dengan metode berbeda.
Mencari model regresi linear paling akurat dengan memperhatikan nilai R squared, AIC dan nilai error.

Data Preperation

Read data / importing data

Sebelum menganalisa dataset, kita akan memanggil dataset “Diamond” dari diamond.csv dan disimpan kedalam object bernama dia, kemudian akan dilakukan Exploratory Data Analyst terhadap dataset tersebut.

dia <- read.csv("diamonds.csv")

Import Library

Pertama kita akan memanggil library yang dibutuhkan, karena kita akan melakukan sedikit visualisasi data dan membuat model regresi linear, maka kita akan menggunkan library sebagai berikut.

library(dplyr)
library(ggplot2)
library(GGally)
library(MLmetrics)
library(performance)
library(lmtest)
library(leaps)
library(car)
library(corrplot)
library(ggpubr)

Inspeksi Data

Selanjutnya kita akan melihat dataset yang sudah kita panggil, gunakan head untuk melihat 6 baris pertama.

head(dia)

Kita juga dapat melihat isi dataset dengan menggunakan str

str(dia)

#> 'data.frame':    53940 obs. of  11 variables:
#>  $ X      : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ carat  : num  0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
#>  $ cut    : chr  "Ideal" "Premium" "Good" "Premium" ...
#>  $ color  : chr  "E" "E" "E" "I" ...
#>  $ clarity: chr  "SI2" "SI1" "VS1" "VS2" ...
#>  $ depth  : num  61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
#>  $ table  : num  55 61 65 58 58 57 57 55 61 61 ...
#>  $ price  : int  326 326 327 334 335 336 336 337 337 338 ...
#>  $ x      : num  3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
#>  $ y      : num  3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
#>  $ z      : num  2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

Dari hasil tinjauan dataset di atas, ternyata dataset terdiri dari 53.940 baris dan 11 kolom. Beberapa kolom ada yang bertipe integer, numerik, dan character.

Deskripsi Data

Dari inspeksi data yang dilakukan terdapat beberapa kolom yang memiliki deskripsi sebagai berikut:

price harga dalam US dollars ($326–$18,823)
carat berat dari berlian (0.2–5.01)
cut kualitas dari hasil potongan (Fair, Good, Very Good, Premium, Ideal), potongan mengacu pada ringkasan proporsi yang memengaruhi kualitas seperti kecemerlangan, kilauan, dll.
color warna berlian, dari J (buruk, agak kekuningan) to D (terbaik, tanpa warna)
clarity ukuran seberapa jernih sebuah berlian dari unsur lain yang terkandung (I1 (terburuk), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (terbaik))
x panjang dalam mm (0–10.74)
y lebar dalam mm (0–58.9)
z tebal dalam mm (0–31.8)
depth total ketebalam dalam persentase = z / mean(x, y) = 2 * z / (x + y) (43–79)
table lebar bagian atas berlian terhadap lebar keseluruhan (43–95)

Data Wrangling

Cek Duplikat

Mengecek apakah terdapat data yang memiliki kesamaan nilai atau duplikat.

#cek duplicate
sum(duplicated(dia))

#> [1] 0

Ternyata tidak ada data yang berisi duplikat value, selanjutnya cek apakah terdapat missing values.

Mengecek kolom Missing Values

colSums(is.na(dia))

#>       X   carat     cut   color clarity   depth   table   price       x       y 
#>       0       0       0       0       0       0       0       0       0       0 
#>       z 
#>       0

dia[]

Ternyata tidak terdapat kolom yang mengandung missing values.

Cek Unique

Mengecek nilai unik di setiap kolom untuk membantu dalam penentuan tipe data terutama tipe data kategorik/factor. Karena ingin mengecek semua kolomnya sekaligus maka digunakan perintah lapply().

lengths(lapply(dia, unique))

#>       X   carat     cut   color clarity   depth   table   price       x       y 
#>   53940     273       5       7       8     184     127   11602     554     552 
#>       z 
#>     375

Cek Tipe Data

Melakukan inspeksi tipe data, untuk memastikan tipe data dari setiap kolomnya sudah sesuai.

glimpse(dia)

#> Rows: 53,940
#> Columns: 11
#> $ X       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
#> $ carat   <dbl> 0.23, 0.21, 0.23, 0.29, 0.31, 0.24, 0.24, 0.26, 0.22, 0.23, 0.…
#> $ cut     <chr> "Ideal", "Premium", "Good", "Premium", "Good", "Very Good", "V…
#> $ color   <chr> "E", "E", "E", "I", "J", "J", "I", "H", "E", "H", "J", "J", "F…
#> $ clarity <chr> "SI2", "SI1", "VS1", "VS2", "SI2", "VVS2", "VVS1", "SI1", "VS2…
#> $ depth   <dbl> 61.5, 59.8, 56.9, 62.4, 63.3, 62.8, 62.3, 61.9, 65.1, 59.4, 64…
#> $ table   <dbl> 55, 61, 65, 58, 58, 57, 57, 55, 61, 61, 55, 56, 61, 54, 62, 58…
#> $ price   <int> 326, 326, 327, 334, 335, 336, 336, 337, 337, 338, 339, 340, 34…
#> $ x       <dbl> 3.95, 3.89, 4.05, 4.20, 4.34, 3.94, 3.95, 4.07, 3.87, 4.00, 4.…
#> $ y       <dbl> 3.98, 3.84, 4.07, 4.23, 4.35, 3.96, 3.98, 4.11, 3.78, 4.05, 4.…
#> $ z       <dbl> 2.43, 2.31, 2.31, 2.63, 2.75, 2.48, 2.47, 2.53, 2.49, 2.39, 2.…

Kolom yang tipe datanya akan diubah: - cut, color, clarity -> Factor

Kolom yang dibuang: - X

Mengubah sekaligus beberapa tipe data dengan menggunakan perintah mutate_at

dia <- dia %>% 
  select(-c(X)) %>% 
  mutate_at(vars(cut, color, clarity), as.factor)

Mengecek kembali tipe data:

# Cek kembali struktur data
glimpse(dia)

#> Rows: 53,940
#> Columns: 10
#> $ carat   <dbl> 0.23, 0.21, 0.23, 0.29, 0.31, 0.24, 0.24, 0.26, 0.22, 0.23, 0.…
#> $ cut     <fct> Ideal, Premium, Good, Premium, Good, Very Good, Very Good, Ver…
#> $ color   <fct> E, E, E, I, J, J, I, H, E, H, J, J, F, J, E, E, I, J, J, J, I,…
#> $ clarity <fct> SI2, SI1, VS1, VS2, SI2, VVS2, VVS1, SI1, VS2, VS1, SI1, VS1, …
#> $ depth   <dbl> 61.5, 59.8, 56.9, 62.4, 63.3, 62.8, 62.3, 61.9, 65.1, 59.4, 64…
#> $ table   <dbl> 55, 61, 65, 58, 58, 57, 57, 55, 61, 61, 55, 56, 61, 54, 62, 58…
#> $ price   <int> 326, 326, 327, 334, 335, 336, 336, 337, 337, 338, 339, 340, 34…
#> $ x       <dbl> 3.95, 3.89, 4.05, 4.20, 4.34, 3.94, 3.95, 4.07, 3.87, 4.00, 4.…
#> $ y       <dbl> 3.98, 3.84, 4.07, 4.23, 4.35, 3.96, 3.98, 4.11, 3.78, 4.05, 4.…
#> $ z       <dbl> 2.43, 2.31, 2.31, 2.63, 2.75, 2.48, 2.47, 2.53, 2.49, 2.39, 2.…

Exploratory Data Analysis

Dalam proses EDA kali ini, kita akan meninjau masing masing kolomnya, dan menemukan apakah ada keganjilan dari masing-masing kolom. Tidak hanya itu kita akan melakukan visualisasi data untuk melihat hubungan dari kolom prediktor dan kolom target.

summary(dia)

#>      carat               cut        color        clarity          depth      
#>  Min.   :0.2000   Fair     : 1610   D: 6775   SI1    :13065   Min.   :43.00  
#>  1st Qu.:0.4000   Good     : 4906   E: 9797   VS2    :12258   1st Qu.:61.00  
#>  Median :0.7000   Ideal    :21551   F: 9542   SI2    : 9194   Median :61.80  
#>  Mean   :0.7979   Premium  :13791   G:11292   VS1    : 8171   Mean   :61.75  
#>  3rd Qu.:1.0400   Very Good:12082   H: 8304   VVS2   : 5066   3rd Qu.:62.50  
#>  Max.   :5.0100                     I: 5422   VVS1   : 3655   Max.   :79.00  
#>                                     J: 2808   (Other): 2531                  
#>      table           price             x                y         
#>  Min.   :43.00   Min.   :  326   Min.   : 0.000   Min.   : 0.000  
#>  1st Qu.:56.00   1st Qu.:  950   1st Qu.: 4.710   1st Qu.: 4.720  
#>  Median :57.00   Median : 2401   Median : 5.700   Median : 5.710  
#>  Mean   :57.46   Mean   : 3933   Mean   : 5.731   Mean   : 5.735  
#>  3rd Qu.:59.00   3rd Qu.: 5324   3rd Qu.: 6.540   3rd Qu.: 6.540  
#>  Max.   :95.00   Max.   :18823   Max.   :10.740   Max.   :58.900  
#>                                                                   
#>        z         
#>  Min.   : 0.000  
#>  1st Qu.: 2.910  
#>  Median : 3.530  
#>  Mean   : 3.539  
#>  3rd Qu.: 4.040  
#>  Max.   :31.800  
#>

Di kolom x, y dan z yang mana ini merupakan kolom ukuran berlian, dari hasil summary ditemukan kalau nilai minimum dari ketiga kolom ini adalah 0, dan nilai 0 ini tidak logis dalam sebuah ukuran sebuah dimensi. Untuk itu data dengan nilai 0 di kolom tersebut harus dihilangkan.

Pertama kita cek ada berapa baris data yang ukuran berliannya bernilai 0

dia %>% 
  filter(x == 0 | y == 0 | z == 0)

Ternyata ada 20 baris, kemudian bisa dihilangkan dengan cara berikut

dia_clean<- dia %>% 
  filter(!(x == 0 | y == 0 | z == 0))

Untuk kolom cut dan clarity, masing-masing kolom ini memiliki urutan dari tingkatan kategorinya yang salah. Dan ini harus diperbaiki karena tingkatan dalam masing-masing kolom tersebut memiliki andil dalam menaikkan harga/price suatu berlian. Maka dari itu kita akan mengurutkan kategorinya berdasarkan urutan yang benar.

# kolom cut
dia_clean$cut <- factor(dia_clean$cut, levels = c("Fair", "Good", "Very Good", "Premium", "Ideal"), ordered = F)
levels(dia_clean$cut)

#> [1] "Fair"      "Good"      "Very Good" "Premium"   "Ideal"

# kolom clarity
dia_clean$clarity <- factor(dia_clean$clarity, c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF"), ordered = F)
levels(dia_clean$clarity)

#> [1] "I1"   "SI2"  "SI1"  "VS2"  "VS1"  "VVS2" "VVS1" "IF"

Balik lagi ke hasil summary tadi, ternyata tampak ada ukuran berlian yang sangat besar terlihat di kolom lebar dan tebal, hal ini bisa menjadi outlier pada dataset ini.

Hubungan Carat dengan Price

ggplot(dia_clean, aes(x = carat, y = price)) + 
  geom_point()

Dari hasil plot di atas, dapat langsung di lihat bahwa hubungan antara carat dan price ada kemungkinan tidak linier. Hal ini harus diperhatikan ketika kita mulai membangun model terutama model regresi.

ggplot(dia_clean, aes(x = carat, y = price))+
  geom_point()+
  geom_smooth()

Bentuk kurva di atas memberikan bukti tambahan bahwa hubungan antara carat dan price bersifat eksponensial, tetapi ketika ditinjau pada carat berlian yang lebih besar dari 2 hubungannya justru melandai. Hal ini mungkin tidak bagus untuk model yang akan kita buat kedepannya.

Hal ini dapat terjadi karena kemungkinan adanya batas atas pada harga berlian dengan carat tertentu di dataset ini, ini dapat dilihat semacam ada batasan harga berlian tidak lebih dari $19000. Kemungkinan kumpulan dataset ini hanya menyertakan berlian yang berada di bawah harga tertentu. Jadi, untuk berlian yang ukurannya sangat besar dalam hal ini carat yang besar, ada bias harga yang terjadi antara ukuran berlian dan harganya.

Hal ini juga mungkin terjadi karena data untuk berlian yang memiliki carat tinggi dalam dataset ini hanya sekumpulan berlian dengan kualitas rendah, sehingga untuk menemukan penyebabnya kita akan memecah scatter plot di atas berdasarkan kategori cut, color dan clarity untuk melihat kategori dari variabel/kolom mana yang mempengaruhi kurva korelasi antara price dan carat sehingga bisa melandai.

Hubungan Carat dengan Price berdasarkan Cut

Untuk plot pertama kita tinjau berdasarkan cut,

ggplot(data = dia_clean) +
  geom_point(mapping = aes(x = carat, y = price, color = cut)) +
  facet_wrap(~cut) +
  labs(title = "Hubungan Price dan Carat Pada Tiap Jenis Cut")

Dari plot di atas tampak kalau semua tipe cut memberikan sedikit outlier pada scatter plot. Kita bisa lihat peningkatan di tipe cut lain selain tipe cut fair bisa dibilang bergerak secara exponensial, sedangkan pada tipe cut fair ada beberapa outlier dimana isinya adalah berlian dengan harga carat tinggi tapi harganya justru tidak wajar, hal ini yang menyebabkan garis di scatter plot utama tadi melandai.

Hubungan Carat dengan Price berdasarkan Color

ggplot(data = dia_clean)+
  geom_point(mapping = aes(x = carat, y = price, color = color))+
  facet_wrap(~color) +
  labs(title = "Hubungan Price dan Carat Pada Tiap Jenis Color")

Sedangkan dilihat dari color, ternyata semua jenis color memiliki outlier. Dan apabila ditinjau dari tipe colornya, tipe J merupakan yang terburuk di dataset ini, dan ternyata color J ini yang memberi pengaruh outlier lebih banyak dibanding yang tipe color yang lain.

Hubungan Carat dengan Price berdasarkan Clarity

ggplot(data = dia_clean)+
  geom_point(mapping = aes(x = carat, y = price, color = clarity))+
  facet_wrap(~clarity) +
  labs(title = "Hubungan Price dan Carat Pada Tiap Jenis Clarity")

Dan terakhir untuk clarity, hanya di tipe I1 yang memberikan sumbangsih outlier terbanyak pada scatter plot utama dan VS2 terlihat hanya 1 outlier. Terlihat di clarity tipe I1 memiliki berlian dengan carat yang besar tetapi peningkatan harganya tidak sebanding dengan berlian dengan clarity tipe lain.

Dari penjabaran scatter plot berdasarkan tipe di setiap kolom/variabel cut, color dan clarity, terlihat jelas bahwa penyebab utama kenapa garis pada scatter plot utama tadi melandai karena banyaknya outlier atau banyaknya berlian dengan carat yang besar tetapi dari segi kualitas justru buruk dalam dataset kali ini.

Karena outlier-outlier tersebut akan memengaruhi model kita dalam memprediksi data baru nantinya, maka kita bisa menghapus beberapa outlier tersebut agar kedepannya model kita tidak salah memprediksi untuk berlian dengan carat yang besar.

Dengan menggunakan bantuan code boxplot.stats(nama_data)$out kita bisa megambil semua outlier yang ada di data.

dia_no_out <- dia_clean %>% 
  filter(!carat %in% boxplot.stats(dia_clean$carat)$out)

Kemudian kita cek kembali scatter plot antra carat dengan price.

ggplot(dia_no_out, aes(x = carat, y = price)) + 
  geom_point() + 
  geom_smooth()

Sekarang tampak garis melandai yang disebabkan oleh outlier sudah hilang.

dia_no_out

Hubungan Price dengan variabel lainnya

plot1 <- ggplot(dia_no_out, aes(x = depth, y = price)) + 
  geom_point()

plot2 <- ggplot(dia_no_out, aes(x = table, y = price)) + 
  geom_point()

plot3 <- ggplot(dia_no_out, aes(x = x, y = price)) + 
  geom_point()

plot4 <- ggplot(dia_no_out, aes(x = y, y = price)) + 
  geom_point()

plot5 <- ggplot(dia_no_out, aes(x = z, y = price)) + 
  geom_point()

ggarrange(plot1, plot2, plot3, plot4, plot5, ncol = 3, nrow = 2)

Tampak ada beberapa outlier terutama pada masing-masing variabel terhadap price, akan tetapi pada table dan depth nilainya itu dalam bentuk persen atau hasil perhitungann jadi kita tidak akan menghapus outliernya, tetapi pada y dan z akan di hapus karena hasil pengukuran, hal ini karena bisa jadi akan mempengaruhi model.

length(boxplot.stats(dia_no_out$x)$out)

#> [1] 0

length(boxplot.stats(dia_no_out$y)$out)

#> [1] 2

length(boxplot.stats(dia_no_out$z)$out)

#> [1] 3

Tampak hanya ada 5 outlier dari variabel x, y, dan z.

dia_no_out1 <- dia_no_out %>% 
  filter(!x %in% boxplot.stats(x)$out) %>% 
  filter(!y %in% boxplot.stats(y)$out) %>% 
  filter(!z %in% boxplot.stats(z)$out)

Mendefinisikan Business Problem

Dalam kasus kali ini, kita akan memprediksi nilai pada kolom price berdasarkan semua atau beberapa variabel prediktor yang ada. Kita akan menentukan:

variable target (y): price.
variable prediktor (x): semua atau beberapa variabel selain kolom price.

Note untuk kasus ini kita akan mencari variabel prediktor mana yang akan menghasilkan model regresi linear terbaik.

Membuat Model Multiple Linear Regression

Linear regression dengan lebih dari satu prediktor bisa meningkatkan performa model karena lebih banyak informasi yang dapat menjelaskan target.

Kita akan membuat beberapa model linear regresi dan membandingkan nilai R-Squarednya.

Simple Linear Regression

Akan dibuat model regresi linear sederhana dimana variable targetnya price dan variabel prediktornya berdasarkan korelasi terkuat terhadap variabel targetnya. Untuk itu kita cek korelasi untuk tiap prediktor terhadap variabel targetnya.

ggcorr(dia_no_out1, label = TRUE, label_round = 3)

Karena carat merupakan variabel prediktor dengan nilai korelasi paling tinggi terhadap variabel target dibanding variabel prediktor lainnya, maka kita bisa mengambilnya sebagai variabel prediktor di model regresi linear sederhana.

model_simple <- lm(formula = price ~ carat, data = dia_no_out1)
summary(model_simple)

#> 
#> Call:
#> lm(formula = price ~ carat, data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -8101.2  -775.1   -22.0   519.5 12769.8 
#> 
#> Coefficients:
#>             Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept) -2222.80      13.45  -165.3 <0.0000000000000002 ***
#> carat        7687.52      15.82   485.8 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1451 on 52031 degrees of freedom
#> Multiple R-squared:  0.8194, Adjusted R-squared:  0.8194 
#> F-statistic: 2.36e+05 on 1 and 52031 DF,  p-value: < 0.00000000000000022

ggplot(dia_no_out, aes(x = carat, y = price)) + 
  geom_point() + 
  geom_smooth(method = "lm")

Insight:

Karena p-value 0.00000000000000022 < 0.05, sehingga disimpulkan bahwa carat memiliki pengaruh yang signifikan terhadap peningkatan harga/price sebuah berlian.
Ketika carat = 0, maka nilai Price kita sebesar -$2222.80. Hal ini logis karena carat tidak mungkin bernilai 0.
Ketika carat naik satu satuan, maka price naik (positif) sebesar $7687.52
Dari hasil model diperoleh Adj R-Squared = 0.8194, Artinya sebesar 81.94% variabel Price bisa dijelaskan oleh variabel Carat. Sisanya 18.06% dijelaskan oleh variabel lain yang tidak diketahui.
Nilai adjR2 yang dihasilkan sudah baik, akan tetapi kita akan menguji model lain dengan penambahan faktor variabel lainnya.

Multiple Linear Regression

Model dengan semua prediktor

model_all <- lm(formula = price ~ ., data = dia_no_out1)
summary(model_all)

#> 
#> Call:
#> lm(formula = price ~ ., data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7586.2  -534.2  -155.1   348.6 10866.6 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)   2205.234    646.404   3.412             0.000646 ***
#> carat        13954.180     66.162 210.909 < 0.0000000000000002 ***
#> cutGood        376.133     31.121  12.086 < 0.0000000000000002 ***
#> cutVery Good   500.165     30.133  16.599 < 0.0000000000000002 ***
#> cutPremium     592.929     29.460  20.127 < 0.0000000000000002 ***
#> cutIdeal       629.758     30.810  20.440 < 0.0000000000000002 ***
#> colorE        -218.988     15.877 -13.792 < 0.0000000000000002 ***
#> colorF        -280.400     16.075 -17.443 < 0.0000000000000002 ***
#> colorG        -493.410     15.765 -31.298 < 0.0000000000000002 ***
#> colorH       -1019.767     16.867 -60.458 < 0.0000000000000002 ***
#> colorI       -1535.520     19.177 -80.072 < 0.0000000000000002 ***
#> colorJ       -2281.876     24.118 -94.614 < 0.0000000000000002 ***
#> claritySI2    1747.080     41.561  42.036 < 0.0000000000000002 ***
#> claritySI1    2673.923     41.218  64.872 < 0.0000000000000002 ***
#> clarityVS2    3261.670     41.383  78.816 < 0.0000000000000002 ***
#> clarityVS1    3574.692     41.951  85.210 < 0.0000000000000002 ***
#> clarityVVS2   3928.479     43.055  91.243 < 0.0000000000000002 ***
#> clarityVVS1   3974.729     44.155  90.018 < 0.0000000000000002 ***
#> clarityIF     4300.185     47.444  90.636 < 0.0000000000000002 ***
#> depth           25.262      9.459   2.671             0.007572 ** 
#> table          -34.828      2.641 -13.186 < 0.0000000000000002 ***
#> x            -1805.687     96.588 -18.695 < 0.0000000000000002 ***
#> y             1119.108     98.946  11.310 < 0.0000000000000002 ***
#> z            -2239.327    150.500 -14.879 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 999.1 on 52009 degrees of freedom
#> Multiple R-squared:  0.9144, Adjusted R-squared:  0.9144 
#> F-statistic: 2.416e+04 on 23 and 52009 DF,  p-value: < 0.00000000000000022

💡 Hasil interpretasi dari summary(model_all):

Significant Predictor:

Variabel yang signifikan mempengaruhi price adalah semua variabel prediktor yakni carat, cut, color, clarity, depth, table, x, y, dan z.
Tidak ada variabel yang tidak signifikan

Slope for Numeric

Interpretasi coefficient untuk prediktor numerik:
- Ketika variabel carat naik 1 satuan, maka variabel price naik sebesar $13954.180 dengan variabel prediktor lainnya tetap.
- Ketika variabel depth naik 1 satuan, maka variabel price naik sebesar $25.262 dengan variabel prediktor lainnya tetap.
- Ketika variabel table naik 1 satuan, maka variabel price turun sebesar $34.828 dengan variabel prediktor lainnya tetap.
- Koefisien cutIdeal = 629.758, artinya berlian dengan potongan/cut tingkat Ideal (best) meningkatkan price sebesar $629.758 apabila dibandingan dengan potongan/cut tingkat Fair (worst).
- Koefisien colorJ = -2281.876, artinya berlian dengan warna/color tingkat J (worst) menurunkan price sebesar $2281.876 apabila dibandingan dengan warna/color berlian dengan tingkat D (best).
- Koefisien clarityIF = 4300.185, artinya berlian dengan kejernihan/calrity tingkat IF (best) menaikkan price sebesar $4300.185 apabila dibandingan dengan kejernihan/calrity tingkat I1 (worst).

Goodness of Fit

Adjusted R-squared: 0.9144 artinya model kita berhasil menangkap informasi/varians sebesar 91.44% untuk memprediksi variabel target.

Model dengan prediktor berdasarkan korelasi

Sebelum membuat model Multiple Linear Regression, kita perlu melihat korelasi antara variabel prediktornya

ggcorr(dia_no_out1, label = TRUE, label_round = 3)

Hasil di atas hanya menampilkan korelasi pada variabel dengan tipe data numerik.

Tampak prediktor x, y dan z memiliki korelasi yang sangat kuat terhadap prediktor carat, hal ini tidak kita inginkan dalam pembuatan model regresi linear.

Dalam kasus regresi linear kita tidak mau adanya korelasi yang kuat antar variabel prediktor hal ini akan mengganggu hasil dari model kita karena terdapatnya multikolinearitas.

Karena kita tidak ingin nilai x, y, dan z termuat di dalam model, kita kana membuat model multiple linear regression tanap ketiga variabel prediktor tersebut.

Selanjutnya kita bisa membuat model regresi yang baru tanpa memasukkan variabel tersebut

model_no_xyz <- lm(formula = price ~ .-x -y -z, data = dia_no_out1)
summary(model_no_xyz)

#> 
#> Call:
#> lm(formula = price ~ . - x - y - z, data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5559.0  -645.6  -187.9   425.7 10390.6 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)  -3734.671    351.078 -10.638 < 0.0000000000000002 ***
#> carat         8783.929     12.833 684.457 < 0.0000000000000002 ***
#> cutGood        521.436     32.418  16.085 < 0.0000000000000002 ***
#> cutVery Good   676.149     31.132  21.719 < 0.0000000000000002 ***
#> cutPremium     715.817     31.172  22.964 < 0.0000000000000002 ***
#> cutIdeal       765.449     32.271  23.719 < 0.0000000000000002 ***
#> colorE        -219.039     16.846 -13.002 < 0.0000000000000002 ***
#> colorF        -316.183     17.049 -18.545 < 0.0000000000000002 ***
#> colorG        -513.313     16.725 -30.692 < 0.0000000000000002 ***
#> colorH        -997.025     17.891 -55.728 < 0.0000000000000002 ***
#> colorI       -1491.312     20.338 -73.326 < 0.0000000000000002 ***
#> colorJ       -2232.200     25.581 -87.261 < 0.0000000000000002 ***
#> claritySI2    1706.376     44.026  38.758 < 0.0000000000000002 ***
#> claritySI1    2660.524     43.644  60.959 < 0.0000000000000002 ***
#> clarityVS2    3315.610     43.826  75.654 < 0.0000000000000002 ***
#> clarityVS1    3640.555     44.409  81.978 < 0.0000000000000002 ***
#> clarityVVS2   4074.712     45.549  89.458 < 0.0000000000000002 ***
#> clarityVVS1   4171.351     46.695  89.332 < 0.0000000000000002 ***
#> clarityIF     4512.524     50.154  89.973 < 0.0000000000000002 ***
#> depth          -13.943      3.835  -3.635             0.000278 ***
#> table          -28.271      2.797 -10.108 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1060 on 52012 degrees of freedom
#> Multiple R-squared:  0.9037, Adjusted R-squared:  0.9036 
#> F-statistic: 2.439e+04 on 20 and 52012 DF,  p-value: < 0.00000000000000022

Tampak hasil dari R-Squarednya lebih rendah 1 % dari model dengan semua variabel predictor.

Model dengan bantuan step_wise

Backward

model_all_backward <- step(object = model_all,
                       direction = "backward",
                       trace = FALSE) #tidak mengeluarkan output
summary(model_all_backward)

#> 
#> Call:
#> lm(formula = price ~ carat + cut + color + clarity + depth + 
#>     table + x + y + z, data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7586.2  -534.2  -155.1   348.6 10866.6 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)   2205.234    646.404   3.412             0.000646 ***
#> carat        13954.180     66.162 210.909 < 0.0000000000000002 ***
#> cutGood        376.133     31.121  12.086 < 0.0000000000000002 ***
#> cutVery Good   500.165     30.133  16.599 < 0.0000000000000002 ***
#> cutPremium     592.929     29.460  20.127 < 0.0000000000000002 ***
#> cutIdeal       629.758     30.810  20.440 < 0.0000000000000002 ***
#> colorE        -218.988     15.877 -13.792 < 0.0000000000000002 ***
#> colorF        -280.400     16.075 -17.443 < 0.0000000000000002 ***
#> colorG        -493.410     15.765 -31.298 < 0.0000000000000002 ***
#> colorH       -1019.767     16.867 -60.458 < 0.0000000000000002 ***
#> colorI       -1535.520     19.177 -80.072 < 0.0000000000000002 ***
#> colorJ       -2281.876     24.118 -94.614 < 0.0000000000000002 ***
#> claritySI2    1747.080     41.561  42.036 < 0.0000000000000002 ***
#> claritySI1    2673.923     41.218  64.872 < 0.0000000000000002 ***
#> clarityVS2    3261.670     41.383  78.816 < 0.0000000000000002 ***
#> clarityVS1    3574.692     41.951  85.210 < 0.0000000000000002 ***
#> clarityVVS2   3928.479     43.055  91.243 < 0.0000000000000002 ***
#> clarityVVS1   3974.729     44.155  90.018 < 0.0000000000000002 ***
#> clarityIF     4300.185     47.444  90.636 < 0.0000000000000002 ***
#> depth           25.262      9.459   2.671             0.007572 ** 
#> table          -34.828      2.641 -13.186 < 0.0000000000000002 ***
#> x            -1805.687     96.588 -18.695 < 0.0000000000000002 ***
#> y             1119.108     98.946  11.310 < 0.0000000000000002 ***
#> z            -2239.327    150.500 -14.879 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 999.1 on 52009 degrees of freedom
#> Multiple R-squared:  0.9144, Adjusted R-squared:  0.9144 
#> F-statistic: 2.416e+04 on 23 and 52009 DF,  p-value: < 0.00000000000000022

Tampak hasil r-square dari step-wise tipe backward sebesar 0.9144.

Forward

modeldia_none <- lm(price ~ 1, dia_no_out1)
model_all_forward <- step(object = modeldia_none,
                      direction = "forward", 
                      scope = list(lower = modeldia_none, upper = model_all),
                      trace = F) #tidak mengeluarkan output
summary(model_all_forward)

#> 
#> Call:
#> lm(formula = price ~ carat + clarity + color + z + table + x + 
#>     cut + y + depth, data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7586.2  -534.2  -155.1   348.6 10866.6 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)   2205.234    646.404   3.412             0.000646 ***
#> carat        13954.180     66.162 210.909 < 0.0000000000000002 ***
#> claritySI2    1747.080     41.561  42.036 < 0.0000000000000002 ***
#> claritySI1    2673.923     41.218  64.872 < 0.0000000000000002 ***
#> clarityVS2    3261.670     41.383  78.816 < 0.0000000000000002 ***
#> clarityVS1    3574.692     41.951  85.210 < 0.0000000000000002 ***
#> clarityVVS2   3928.479     43.055  91.243 < 0.0000000000000002 ***
#> clarityVVS1   3974.729     44.155  90.018 < 0.0000000000000002 ***
#> clarityIF     4300.185     47.444  90.636 < 0.0000000000000002 ***
#> colorE        -218.988     15.877 -13.792 < 0.0000000000000002 ***
#> colorF        -280.400     16.075 -17.443 < 0.0000000000000002 ***
#> colorG        -493.410     15.765 -31.298 < 0.0000000000000002 ***
#> colorH       -1019.767     16.867 -60.458 < 0.0000000000000002 ***
#> colorI       -1535.520     19.177 -80.072 < 0.0000000000000002 ***
#> colorJ       -2281.876     24.118 -94.614 < 0.0000000000000002 ***
#> z            -2239.327    150.500 -14.879 < 0.0000000000000002 ***
#> table          -34.828      2.641 -13.186 < 0.0000000000000002 ***
#> x            -1805.687     96.588 -18.695 < 0.0000000000000002 ***
#> cutGood        376.133     31.121  12.086 < 0.0000000000000002 ***
#> cutVery Good   500.165     30.133  16.599 < 0.0000000000000002 ***
#> cutPremium     592.929     29.460  20.127 < 0.0000000000000002 ***
#> cutIdeal       629.758     30.810  20.440 < 0.0000000000000002 ***
#> y             1119.108     98.946  11.310 < 0.0000000000000002 ***
#> depth           25.262      9.459   2.671             0.007572 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 999.1 on 52009 degrees of freedom
#> Multiple R-squared:  0.9144, Adjusted R-squared:  0.9144 
#> F-statistic: 2.416e+04 on 23 and 52009 DF,  p-value: < 0.00000000000000022

Both

model_all_both <- step(object = modeldia_none,
                   direction = "both",
                   scope = list(upper = model_all),
                   trace = F)
summary(model_all_both)

#> 
#> Call:
#> lm(formula = price ~ carat + clarity + color + z + table + x + 
#>     cut + y + depth, data = dia_no_out1)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -7586.2  -534.2  -155.1   348.6 10866.6 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)   2205.234    646.404   3.412             0.000646 ***
#> carat        13954.180     66.162 210.909 < 0.0000000000000002 ***
#> claritySI2    1747.080     41.561  42.036 < 0.0000000000000002 ***
#> claritySI1    2673.923     41.218  64.872 < 0.0000000000000002 ***
#> clarityVS2    3261.670     41.383  78.816 < 0.0000000000000002 ***
#> clarityVS1    3574.692     41.951  85.210 < 0.0000000000000002 ***
#> clarityVVS2   3928.479     43.055  91.243 < 0.0000000000000002 ***
#> clarityVVS1   3974.729     44.155  90.018 < 0.0000000000000002 ***
#> clarityIF     4300.185     47.444  90.636 < 0.0000000000000002 ***
#> colorE        -218.988     15.877 -13.792 < 0.0000000000000002 ***
#> colorF        -280.400     16.075 -17.443 < 0.0000000000000002 ***
#> colorG        -493.410     15.765 -31.298 < 0.0000000000000002 ***
#> colorH       -1019.767     16.867 -60.458 < 0.0000000000000002 ***
#> colorI       -1535.520     19.177 -80.072 < 0.0000000000000002 ***
#> colorJ       -2281.876     24.118 -94.614 < 0.0000000000000002 ***
#> z            -2239.327    150.500 -14.879 < 0.0000000000000002 ***
#> table          -34.828      2.641 -13.186 < 0.0000000000000002 ***
#> x            -1805.687     96.588 -18.695 < 0.0000000000000002 ***
#> cutGood        376.133     31.121  12.086 < 0.0000000000000002 ***
#> cutVery Good   500.165     30.133  16.599 < 0.0000000000000002 ***
#> cutPremium     592.929     29.460  20.127 < 0.0000000000000002 ***
#> cutIdeal       629.758     30.810  20.440 < 0.0000000000000002 ***
#> y             1119.108     98.946  11.310 < 0.0000000000000002 ***
#> depth           25.262      9.459   2.671             0.007572 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 999.1 on 52009 degrees of freedom
#> Multiple R-squared:  0.9144, Adjusted R-squared:  0.9144 
#> F-statistic: 2.416e+04 on 23 and 52009 DF,  p-value: < 0.00000000000000022

summary(model_all_backward)$adj.r.squared

#> [1] 0.9143889

summary(model_all_forward)$adj.r.squared

#> [1] 0.9143889

summary(model_all_both)$adj.r.squared

#> [1] 0.9143889

Ternyata dari ketiga jenis step-wise, ketiga model menghasilkan nilai adjusted R-squared yang sama yakni 0.9143889. Hasil ini sama dengan hasil adj. R-squared di model_all.

Model dengan bantuan regsubsets

library(mltools)
library(data.table)

regs <- regsubsets(price ~ ., dia_clean, nbest=10)
plot(regs, scale="adjr", main="All possible regression: ranked by Adjusted R-squared")

Dilihat dari hasil adj r-squared dari kombinasi prediktor diatas, tampak hasil paling tinggi hanya di 0.89, nilai ini lebih rendah dari hasil adj r-squared model sebelumnya.

Evaluasi Model

Performa Model

Salah satu cara untuk melihat performa model yaitu dengan melihat nilai RMSE dari model kita.

RMSE adalah bentuk akar kuadrat dari MSE. Karena sudah diakarkan, maka interpretasinya kurang lebih sama dengan MAE. RMSE dapat digunakan jika kita lebih concern dengan error yang sangat besar.

\[ RMSE = \sqrt{\frac{1}{n} \sum (\hat y - y)^2} \] Sebelum menghitung nilai RMSE dengan menggunakan perintah RMSE, kita harus menyimpan hasil prediksi dari model kita dalam dataframe.

Kita akan fokus ke satu model untuk menghitung RMSE yaitu model_no_xyz.

dia_pred <- dia_no_out1
dia_pred$pred_price <- predict(model_no_xyz, dia_no_out1)

# Menghitung RMSE dari library MLmetrics
RMSE(y_pred = dia_pred$pred_pric,
     y_true = dia_no_out1$price)

#> [1] 1059.843

range(dia_no_out1$price)

#> [1]   326 18818

Interpretasi RMSE: Nilai RMSE dari model_no_xyz masih berada di dalam rentang data aktualnya ini artinya model masih memberikan kesalahan yang cukup tinggi untuk prediksi apabila dibandingkan dengan nilai aktual target variabelnya. Oleh karena itu masih perlu perbaikan model.

Kita akan membandingkan performa keenam model regression sebelumnya menggunakan fungsi compare_performance dari packages performance

library(performance)
#install.packages("performance")
comparison <- compare_performance(model_simple, model_all, model_no_xyz, model_all_backward, model_all_forward, model_all_both)

as.data.frame(comparison)

Kesimpulan: Model terbaik berdasarkan adjusted r-squared terbesar, AIC terkecil dan RMSE terkecil adalah - adjusted r-squared terbesar -> model all. ini untuk model terbaik - AIC terkecil -> model all - RMSE terkecil -> model_all. ini untuk melihat prediksi terbaik

Apabila ingin melakukan prediksi gunakan model_all

Akan tetapi karena kita tidak ingin adanya multikolinearitas dalam model kita, maka kita akan menggunakan model_no_xyz. Dan model ini juga bagus karena adj. R-squared nya hanya memiliki selisih 1% dengan model_all

Uji Asumsi

Asumsi Linear Regression

Sebagai salah satu model statistik, linear regression adalah model yang ketat asumsi. Berikut beberapa asumsi yang harus dicek untuk memastikan apakah model yang kita buat dianggap sebagai Best Linear Unbiased Estimator (BLUE) model, yaitu model yang dapat memprediksi data baru secara konsisten.

Asumsi model linear regression:

Linearity
Normality of Residuals
Homoscedasticity of Residuals
No Multicollinearity

Linearity

Linearity artinya target variabel dengan prediktornya memiliki hubungan yang linear atau hubungannya bersifat garis lurus.

Untuk menentukan asumsi linearitas dari sebuah multiple linear regression dapat dilakukan dengan membuat plot residual vs fitted values. Plot ini merupakan scatter plot dengan sumbu x adalah nilai fitted values (hasil prediksi variabel respon / nilai prediksi) dan sumbu y adalah nilai residual / nilai error

plot(model_no_xyz, which = 1)
abline(h = 1200, col = "green")
abline(h = -1200, col = "green")

Kesimpulan: karena garis merah sudah keluar dari cakupan toleransi kita, sehingga model_no_xyz adalah model yang tidak linear.

Normality of Residuals

Model linear regression diharapkan menghasilkan error yang berdistribusi normal. Dengan begitu, error lebih banyak berkumpul di sekitar angka nol.

Visualisasi histogram residual menggunakan fungsi hist()

# histogram residual
hist(model_no_xyz$residuals)

Uji statistik dengan ad.test()

Shapiro-Wilk hypothesis test:

H0: error berdistribusi normal / p-value>0.05
H1: error TIDAK berdistribusi normal

Kondisi yang diharapkan: H0

# shapiro test dari residual
library(nortest)
ad.test(model_no_xyz$residuals)

#> 
#>  Anderson-Darling normality test
#> 
#> data:  model_no_xyz$residuals
#> A = 1269.6, p-value < 0.00000000000000022

Kesimpulan: Karena p-values < 0.05, sehingga tolak H0 atau error TIDAK berdistribusi normal.

Homoscedasticity of Residuals

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola. Kondisi ini disebut juga sebagai homoscedasticity

Cara melakukan pengecekan asumsi homoscedasticity:

Uji statistik dengan bptest() dari package lmtest

Breusch-Pagan hypothesis test:

H0: error menyebar konstan atau homoscedasticity
H1: error menyebar TIDAK konstan atau heteroscedasticity

Kondisi yang diharapkan: H0

# bptest dari model
library(lmtest)
bptest(model_no_xyz)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_no_xyz
#> BP = 6330.3, df = 20, p-value < 0.00000000000000022

Kesimpulan: karena nilai p-value dari bptest < 0.05, sehingga tolak H0 atau error menyebar TIDAK konstan atau heteroscedasticity. Sehingga asumsi Homoscedasticity TIDAK terpenuhi.

No Multicollinearity

Multicollinearity adalah kondisi adanya korelasi antar prediktor yang kuat. Hal ini tidak diinginkan karena menandakan prediktor redundan pada model, yang seharusnya dapat dipilih salah satu saja dari variable yang hubungannya amat kuat tersebut. Harapannya tidak terjadi multicollinearity

Uji VIF (Variance Inflation Factor) dengan fungsi vif() dari package car: * nilai VIF > 10: terjadi multicollinearity pada model * nilai VIF < 10: tidak terjadi multicollinearity pada model

Kondisi yang diharapkan: VIF < 10

# vif dari model
library(car)
vif(model_no_xyz)

#>             GVIF Df GVIF^(1/(2*Df))
#> carat   1.232804  1        1.110317
#> cut     1.911708  4        1.084370
#> color   1.126551  6        1.009980
#> clarity 1.256039  7        1.016416
#> depth   1.369454  1        1.170237
#> table   1.789993  1        1.337906

Kesimpulan: Dari uji VIF, prediktor di model_no_xyz lolos uji asumsi multicolinearity (karena tidak ada nilai GVIF >10).

Karena dari nilai RMSE dan uji asumsi keduanya belum terpenuhi, maka kita akan melakukan tuning pada model

Perbaikan Model

Dari uji asumsi yang dilakukan, model_no_xyz menunjukan ketidaklinearan, Tidak berdistribusi normal, dan heteroscedasticity. Dengan demikian akan dilakukan transformasi data untuk variabel prediktor dan respon. Salah satu jenis transformasi yang dapat dilakukan adalah transformasi logaritmik.

Tuning Model

Karena variabel predikotr depth dari hasil summary model_no_xyz memberikan p-value yang terbesar dibanding variabel prediktor lainnya maka akan dipertimbangkan untuk melakukan transformasi logaritmik pada variabel price dan variabel depth. Dan juga dilihat dari scatter plor antara price dan carat kita akan menggunakan logaritmik pada Price untuk membuat hubungan lebih linier.

ggplot(dia_no_out1, aes(x = carat, y = log(price))) + 
  geom_point() + 
  geom_smooth()

Kemudian kita akan membuat model baru berdasarkan asumsi tadi.

model_tuning <- lm(formula = log(price) ~ log(carat) + cut + color + clarity + log(depth) + table, data = dia_no_out1)
summary(model_tuning)

#> 
#> Call:
#> lm(formula = log(price) ~ log(carat) + cut + color + clarity + 
#>     log(depth) + table, data = dia_no_out1)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.76627 -0.08589 -0.00154  0.08106  1.90155 
#> 
#> Coefficients:
#>                Estimate Std. Error  t value            Pr(>|t|)    
#> (Intercept)   8.0411303  0.1313519   61.218 <0.0000000000000002 ***
#> log(carat)    1.8895819  0.0011758 1606.998 <0.0000000000000002 ***
#> cutGood       0.0752561  0.0040270   18.688 <0.0000000000000002 ***
#> cutVery Good  0.1119749  0.0038572   29.030 <0.0000000000000002 ***
#> cutPremium    0.1367829  0.0038582   35.453 <0.0000000000000002 ***
#> cutIdeal      0.1561895  0.0039919   39.127 <0.0000000000000002 ***
#> colorE       -0.0548194  0.0020979  -26.130 <0.0000000000000002 ***
#> colorF       -0.0961137  0.0021239  -45.253 <0.0000000000000002 ***
#> colorG       -0.1612427  0.0020827  -77.419 <0.0000000000000002 ***
#> colorH       -0.2512101  0.0022244 -112.932 <0.0000000000000002 ***
#> colorI       -0.3764010  0.0025271 -148.945 <0.0000000000000002 ***
#> colorJ       -0.5105059  0.0031785 -160.611 <0.0000000000000002 ***
#> claritySI2    0.3811098  0.0054817   69.524 <0.0000000000000002 ***
#> claritySI1    0.5467700  0.0054328  100.643 <0.0000000000000002 ***
#> clarityVS2    0.6979403  0.0054569  127.901 <0.0000000000000002 ***
#> clarityVS1    0.7678174  0.0055295  138.858 <0.0000000000000002 ***
#> clarityVVS2   0.9031873  0.0056752  159.145 <0.0000000000000002 ***
#> clarityVVS1   0.9757102  0.0058211  167.617 <0.0000000000000002 ***
#> clarityIF     1.0709476  0.0062522  171.292 <0.0000000000000002 ***
#> log(depth)   -0.0257261  0.0291976   -0.881               0.378    
#> table        -0.0004313  0.0003483   -1.238               0.216    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.132 on 52012 degrees of freedom
#> Multiple R-squared:  0.9815, Adjusted R-squared:  0.9815 
#> F-statistic: 1.381e+05 on 20 and 52012 DF,  p-value: < 0.00000000000000022

Tampak hasil R-Squared jauh lebih tinggi yakni 0.9815

Dan terdapat variabel yang tidak signifikan yakni log(depth) dan table, sehingga variabel ini bisa dihapus dari model, karena tidak berpengaruh banyak ke nilai target.

Sehingga bisa dibuat model baru tanpa variabel log(depth) dan table.

model_tuning_new <- lm(formula = log(price) ~ log(carat) + cut + color + clarity , data = dia_no_out1)
summary(model_tuning_new)

#> 
#> Call:
#> lm(formula = log(price) ~ log(carat) + cut + color + clarity, 
#>     data = dia_no_out1)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.76911 -0.08589 -0.00159  0.08104  1.90376 
#> 
#> Coefficients:
#>               Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)   7.908547   0.006070 1302.89 <0.0000000000000002 ***
#> log(carat)    1.889447   0.001170 1614.75 <0.0000000000000002 ***
#> cutGood       0.076010   0.003958   19.20 <0.0000000000000002 ***
#> cutVery Good  0.113246   0.003688   30.71 <0.0000000000000002 ***
#> cutPremium    0.137951   0.003655   37.75 <0.0000000000000002 ***
#> cutIdeal      0.158335   0.003617   43.77 <0.0000000000000002 ***
#> colorE       -0.054833   0.002098  -26.14 <0.0000000000000002 ***
#> colorF       -0.096092   0.002124  -45.24 <0.0000000000000002 ***
#> colorG       -0.161234   0.002082  -77.44 <0.0000000000000002 ***
#> colorH       -0.251234   0.002223 -113.00 <0.0000000000000002 ***
#> colorI       -0.376447   0.002526 -149.04 <0.0000000000000002 ***
#> colorJ       -0.510567   0.003177 -160.69 <0.0000000000000002 ***
#> claritySI2    0.381276   0.005479   69.59 <0.0000000000000002 ***
#> claritySI1    0.546918   0.005431  100.70 <0.0000000000000002 ***
#> clarityVS2    0.698117   0.005454  128.00 <0.0000000000000002 ***
#> clarityVS1    0.768032   0.005525  139.01 <0.0000000000000002 ***
#> clarityVVS2   0.903400   0.005671  159.30 <0.0000000000000002 ***
#> clarityVVS1   0.975926   0.005817  167.78 <0.0000000000000002 ***
#> clarityIF     1.071260   0.006245  171.55 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.132 on 52014 degrees of freedom
#> Multiple R-squared:  0.9815, Adjusted R-squared:  0.9815 
#> F-statistic: 1.535e+05 on 18 and 52014 DF,  p-value: < 0.00000000000000022

Ternyata adj R-squared dari model ini tetap sama dengan model_tuning yaitu sebesar 0.9815

Performa Model

Kita akan membandingkan performa kedua model regression yang telah dituning menggunakan fungsi compare_performance dari packages performance

library(performance)
#install.packages("performance")
comparison1 <- compare_performance(model_tuning, model_tuning_new)

as.data.frame(comparison1)

Dapat dilihat bahwa dari kedua model, ternyata model_tuning_new memiliki kriteria performa yang lebih unggul dari model_tuning yaitu Adj R-squared terbesar dan AIC lebih kecil, walaupun RMSE-nya lebih besar tapi hanya beda 0.000002.

Sehingga disimpulkan bahwa model regresi terbaik untuk memodelkan data pada kasus ini adalah model_tuniing_new.

Uji Asumsi

Linearity

Linearity artinya target variabel dengan prediktornya memiliki hubungan yang linear atau hubungannya bersifat garis lurus.

plot(model_tuning_new, which = 1)
abline(h = 0.2, col = "green")
abline(h = -0.2, col = "green")

Kesimpulan: karena garis merah berada dalam cakupan toleransi kita, sehingga model_no_xyz adalah model yang linear.

Normality of Residuals

Model linear regression diharapkan menghasilkan error yang berdistribusi normal. Dengan begitu, error lebih banyak berkumpul di sekitar angka nol.

Visualisasi histogram residual menggunakan fungsi hist()

# histogram residual
  hist(model_tuning_new$residuals)

Uji statistik dengan ad.test()

Shapiro-Wilk hypothesis test:

H0: error berdistribusi normal / p-value>0.05
H1: error TIDAK berdistribusi normal

Kondisi yang diharapkan: H0

# shapiro test dari residual
library(nortest)
ad.test(model_tuning_new$residuals)

#> 
#>  Anderson-Darling normality test
#> 
#> data:  model_tuning_new$residuals
#> A = 29.407, p-value < 0.00000000000000022

Kesimpulan: Karena p-values < 0.05, sehingga tolak H0 atau error TIDAK berdistribusi normal.

Homoscedasticity of Residuals

Diharapkan error yang dihasilkan oleh model menyebar secara acak atau dengan variasi konstan. Apabila divisualisasikan maka error tidak berpola. Kondisi ini disebut juga sebagai homoscedasticity

Cara melakukan pengecekan asumsi homoscedasticity:

Uji statistik dengan bptest() dari package lmtest

Breusch-Pagan hypothesis test:

H0: error menyebar konstan atau homoscedasticity
H1: error menyebar TIDAK konstan atau heteroscedasticity

Kondisi yang diharapkan: H0

# bptest dari model
library(lmtest)
bptest(model_tuning_new)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_tuning_new
#> BP = 1121.6, df = 18, p-value < 0.00000000000000022

Kesimpulan: karena nilai p-value dari bptest < 0.05, sehingga tolak H0 atau error menyebar TIDAK konstan atau heteroscedasticity. Sehingga asumsi Homoscedasticity TIDAK terpenuhi.

No Multicollinearity

Uji VIF (Variance Inflation Factor) dengan fungsi vif() dari package car: * nilai VIF > 10: terjadi multicollinearity pada model * nilai VIF < 10: tidak terjadi multicollinearity pada model

Kondisi yang diharapkan: VIF < 10

# vif dari model
library(car)
vif(model_tuning_new)

#>                GVIF Df GVIF^(1/(2*Df))
#> log(carat) 1.241616  1        1.114278
#> cut        1.096451  4        1.011576
#> color      1.110500  6        1.008772
#> clarity    1.280887  7        1.017840

Kesimpulan: Dari uji VIF, prediktor di model_no_xyz lolos uji asumsi multicolinearity (karena tidak ada nilai GVIF >10).

Hasil dari uji asumsi menunjukkan bahwa 2 dari 4 asumsi masih tidak terpenuhi sehingga bisa dikatakan data berlian tidak terlalu cocok untuk menggunakan model regresi ini.

model_tuning_new <- lm(formula = log(price) ~ log(carat) + cut + color + clarity , data = dia_no_out1)
summary(model_tuning_new)

#> 
#> Call:
#> lm(formula = log(price) ~ log(carat) + cut + color + clarity, 
#>     data = dia_no_out1)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.76911 -0.08589 -0.00159  0.08104  1.90376 
#> 
#> Coefficients:
#>               Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)   7.908547   0.006070 1302.89 <0.0000000000000002 ***
#> log(carat)    1.889447   0.001170 1614.75 <0.0000000000000002 ***
#> cutGood       0.076010   0.003958   19.20 <0.0000000000000002 ***
#> cutVery Good  0.113246   0.003688   30.71 <0.0000000000000002 ***
#> cutPremium    0.137951   0.003655   37.75 <0.0000000000000002 ***
#> cutIdeal      0.158335   0.003617   43.77 <0.0000000000000002 ***
#> colorE       -0.054833   0.002098  -26.14 <0.0000000000000002 ***
#> colorF       -0.096092   0.002124  -45.24 <0.0000000000000002 ***
#> colorG       -0.161234   0.002082  -77.44 <0.0000000000000002 ***
#> colorH       -0.251234   0.002223 -113.00 <0.0000000000000002 ***
#> colorI       -0.376447   0.002526 -149.04 <0.0000000000000002 ***
#> colorJ       -0.510567   0.003177 -160.69 <0.0000000000000002 ***
#> claritySI2    0.381276   0.005479   69.59 <0.0000000000000002 ***
#> claritySI1    0.546918   0.005431  100.70 <0.0000000000000002 ***
#> clarityVS2    0.698117   0.005454  128.00 <0.0000000000000002 ***
#> clarityVS1    0.768032   0.005525  139.01 <0.0000000000000002 ***
#> clarityVVS2   0.903400   0.005671  159.30 <0.0000000000000002 ***
#> clarityVVS1   0.975926   0.005817  167.78 <0.0000000000000002 ***
#> clarityIF     1.071260   0.006245  171.55 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.132 on 52014 degrees of freedom
#> Multiple R-squared:  0.9815, Adjusted R-squared:  0.9815 
#> F-statistic: 1.535e+05 on 18 and 52014 DF,  p-value: < 0.00000000000000022

Interpretasi Model

💡 Hasil interpretasi dari summary(model_tuning_all):

Significant Predictor:

Variabel yang signifikan mempengaruhi price dalam model ini adalah variabel prediktor log(carat), cut, color, dan clarity.
Tidak ada variabel yang tidak signifikan dalam model ini

Slope for Numeric

Interpretasi coefficient untuk prediktor numerik:
- Ketika variabel carat naik 1 satuan, maka variabel price naik sebesar 1.889447 kali dengan variabel prediktor lainnya tetap.
- Koefisien cutIdeal = e^(0.158335) = 1.171559, artinya berlian dengan potongan/cut tingkat Ideal (best) meningkatkan price sebesar 1.171559 kali apabila dibandingan dengan potongan/cut tingkat Fair (worst).
- Koefisien colorJ = e^(-0.510567) = 0.6001552, artinya berlian dengan warna/color tingkat J (worst) menurunkan price sebesar 0.6001552 kali apabila dibandingan dengan warna/color berlian dengan tingkat D (best).
- Koefisien clarityIF = e^(1.071260), artinya berlian dengan kejernihan/calrity tingkat IF (best) menaikkan price sebesar 2.919055 kali apabila dibandingan dengan kejernihan/calrity tingkat I1 (worst).

Goodness of Fit

Adjusted R-squared: 0.9815 artinya model kita berhasil menangkap informasi/varians sebesar 98.15% untuk memprediksi variabel target.

Kesimpulan

Dari hasil analisis dataset diamond di atas, terdapat empat faktor yang mempengaruhi harga berlian. Faktor utama tersebut adalah karat berlian, warna, potongan dan kejernihannya. Namun, karat tampaknya menjadi salah satu faktor yang paling berpengaruh terhadap harga berlian.
Adanya batas atas pada harga berlian dengan carat tertentu di dataset ini, membuat semacam batasan harga berlian yang tidak lebih dari $19000. Kemungkinan kumpulan dataset ini hanya menyertakan berlian yang berada di bawah harga tertentu. Jadi, untuk berlian yang ukurannya sangat besar dalam hal ini carat yang besar, ada bias harga yang terjadi antara ukuran berlian dan harganya. Dan setelah melakukan analisis, ternyata data berlian memiliki banyak outlier yaitu berlian dengan carat yang besar tetapi dari segi kualitas justru buruk dalam dataset kali ini.
Dari hasil model regresi linear bisa dilihat bahwa semakin baik jenis potongan, warna, dan kejernihan maka semakin naik pula harga berlian tersebut. Dan carat merupakan variabel prediktor yang sangat signifikan dalam meningkatkan harga suatu berlian
Ukuran berlian dalam model kali ini tidak begitu signifikan dalam mengubah harga berlian, hal ini dikarenakan variabel-variabel tersebut sudah terwakili di variabel carat.
Untuk hasil model multiple linear regression yang dibuat, tampak model_tuning_new merupakan model terbaik dibanding model-model yang lain. Akan tetapi setelah melakuakan uji asumsi, 2 dari 4 asumsi regresi linear tidak terpenuhi, sehingga model tersebut tidak cocok digunakan untuk dataset kali ini apabila model akan digunakan dalam kebutuhan penelitian atau research di ranah pendidikan. Kita bisa mengabaikan ketidak-terpenuhan dari uji asusmi jika dataset ini kita analisis untuk keperluan bisnis.

Mencari Model Linear Regression terbaik untuk dataset Diamond

LBB: Regression Model

Tri Saputra Ungko

March 18, 2023

Introduction

Data Preperation

Read data / importing data

Import Library

Inspeksi Data

Deskripsi Data

Data Wrangling

Cek Duplikat

Mengecek kolom Missing Values

Cek Unique

Cek Tipe Data

Exploratory Data Analysis

Hubungan Carat dengan Price

Hubungan Carat dengan Price berdasarkan Cut

Hubungan Carat dengan Price berdasarkan Color

Hubungan Carat dengan Price berdasarkan Clarity

Hubungan Price dengan variabel lainnya

Mendefinisikan Business Problem

Membuat Model Multiple Linear Regression

Simple Linear Regression

Multiple Linear Regression

Model dengan semua prediktor

Significant Predictor:

Slope for Numeric

Goodness of Fit

Model dengan prediktor berdasarkan korelasi

Model dengan bantuan step_wise

Backward

Forward

Both

Model dengan bantuan regsubsets

Evaluasi Model

Performa Model

Uji Asumsi

Asumsi Linear Regression

Linearity

Normality of Residuals

Homoscedasticity of Residuals

No Multicollinearity

Perbaikan Model

Tuning Model

Performa Model

Uji Asumsi

Linearity

Normality of Residuals

Homoscedasticity of Residuals

No Multicollinearity

Interpretasi Model

Significant Predictor:

Slope for Numeric

Goodness of Fit

Kesimpulan