Dari bagian sebelumnya, mendapat ilham tentang model dan model building. Prinsip dasar dari membangun sebuah model melibat serangkaian proses tidak berurutan dan kontinuu terkait : * Eksplorasi — > Menemukan pattern * Pemilihan family model (berdasarkan kebutuhan dan kondisi riil data) * Model fitting * Repeat Satu hal yang perlu kita pahami adalah, carilah model paling optimal; yaitu model yang menghasilkan residual paling minimal. Membangun model hampir mirip dengan melukis. Kita harus tahu kapan untuk berhenti mencari kesempurnaan dari karya kita, dan berani menguji cobanya ke lapangan.
Untuk pembahasan ini, kita akan mulai dengan load seluruh data dan package yang kita butuhkan.
Jika kita ingat kembali pada dataset ‘diamonds’, ada satu fakta menarik di mana berlian berkualitas rendah (poor cuts, bad color, dan inferior clarity) memiliki harga lebih mahal daripada berlian berkualitas tinggi.
diamonds %>% ggplot(aes(cut, price)) +geom_boxplot()
diamonds %>% ggplot(aes(color,price)) + geom_boxplot()
diamonds %>% ggplot(aes(clarity, price)) + geom_boxplot()
Kenapa bisa seperti itu?
Jawabannya ada pada eksplorasi data carat vs price.
ggplot(diamonds, aes(carat, price)) + geom_hex(bins = 50)
# Jika bingung, silahkan lihat pada grafik garis berikut
ggplot(diamonds, aes(carat, price)) + geom_line()
Yup, grafik di atas menunjukkan hubungan linier positif antara carat dengan price diamond; dan hubungan itu sangat kuat. Kita mungkin bisa membuat hipotesis bahwa karat merupakan faktor siginigikan dalam penentuan harga diamond.
Namun untuk memperkuat analisa. Mari kita analisa bagaimana kondisi karat dalam berbagai kategori kualitas diamond
ggplot(diamonds, aes(cut, carat)) + geom_boxplot()
ggplot(diamonds, aes(clarity, carat)) + geom_boxplot()
ggplot(diamonds, aes(color, carat)) + geom_boxplot()
Pada boxplot ditunjukkan bahwa diamond dengan carat tinggi terdistribusi lebih banyak pada diamond dengan kualitas renda. Ini memperkuat hipotesis kita bahwa carat merupakan penentu utama dari harga diamond.
Setelah memiliki hipotesis secara visual, mari kita coba buktikan dugaan kita melalui modelling. Kita tidak bisa terlalu cepat mengambil kesimpulan carat sebagai satu-satunya variabel yang berpengaruh terhadap harga bukan?
Nah, sebelum lebih jauh, satu hal yg perlu kita pahami saat berhadapan dengan big data seperti diamonds, ada baiknya kita memperhatikan distribusi masing-masing data. (sebenarnya, lebih sempurna jika kita melakukan ini pada setiap proses analisa data)
Untuk melakukannya, kita bisa menggunakan fungsi summary().
summary(diamonds)
## 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 Very Good:12082 F: 9542 SI2 : 9194 Median :61.80
## Mean :0.7979 Premium :13791 G:11292 VS1 : 8171 Mean :61.75
## 3rd Qu.:1.0400 Ideal :21551 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
##
Dari summary ini kita bisa mengetahu gambaran umum tentang sebaran data. Namun ada satu parameter yang kurang di sini yaitu modus (mode). Untuk itu, mari kita buat terlebih dahulu fungsi modus kita yang kita beri nama modes().
# mode fucntion
modes <- function(x){
ux <- unique(x)
tab <- tabulate(match(x,ux))
ux[tab == max(tab)]
}
Setelah itu, kita mulai lakukan analisa modus.
modes(diamonds$carat)
## [1] 0.3
modes(diamonds$depth)
## [1] 62
modes(diamonds$table)
## [1] 56
modes(diamonds$price)
## [1] 605
modes(diamonds$x)
## [1] 4.37
modes(diamonds$y)
## [1] 4.34
modes(diamonds$z)
## [1] 2.7
Di sini kita bisa mengetahui pattern terkuat berada dari variabel ‘carat’ terhadap ‘price’.
ggplot(diamonds,aes(carat, price)) +
geom_hex(bins = 50)
Jika kita ingin mengevaluasi dampak attribut diamond lain terhadap harga, kita bisa membangun model untuk memisahkan efek dari ‘carat’ Caranya : * Pertama, kita transformasikan carat dan price menjadi log.
diamonds2 <- diamonds %>%
mutate(lprice = log2(price), lcarat = log2(carat))
diamonds2
## # A tibble: 53,940 x 12
## carat cut color clarity depth table price x y z lprice lcarat
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43 8.35 -2.12
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31 8.35 -2.25
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31 8.35 -2.12
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63 8.38 -1.79
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75 8.39 -1.69
## 6 0.24 Very G~ J VVS2 62.8 57 336 3.94 3.96 2.48 8.39 -2.06
## 7 0.24 Very G~ I VVS1 62.3 57 336 3.95 3.98 2.47 8.39 -2.06
## 8 0.26 Very G~ H SI1 61.9 55 337 4.07 4.11 2.53 8.40 -1.94
## 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 2.49 8.40 -2.18
## 10 0.23 Very G~ H VS1 59.4 61 338 4 4.05 2.39 8.40 -2.12
## # ... with 53,930 more rows
Apa gunanya transformasi log? Hal ini untuk mempermudah kita melihat hubungan antara carat dan price. Transformasi log akan membuat pattern linier dan memudahkan kita untuk bekerja.
ggplot(diamonds2, aes(lcarat, lprice)) + geom_hex(bins = 50)
Lihat? hasil transformasi log membuat hubungan linier kedua variabel di atas semakin jelas. Langkah selanjutnya dalam mengeluarkan pola linier dari data adalah dengan membuat hubungan carat vs price eksplisit melalui model fitting.
mod_diamonds <- lm(lprice~lcarat, data = diamonds2)
mod_diamonds
##
## Call:
## lm(formula = lprice ~ lcarat, data = diamonds2)
##
## Coefficients:
## (Intercept) lcarat
## 12.189 1.676
summary(mod_diamonds)
##
## Call:
## lm(formula = lprice ~ lcarat, data = diamonds2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.17606 -0.24456 -0.00853 0.24002 1.93023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.188841 0.001969 6190.9 <2e-16 ***
## lcarat 1.675817 0.001934 866.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3789 on 53938 degrees of freedom
## Multiple R-squared: 0.933, Adjusted R-squared: 0.933
## F-statistic: 7.51e+05 on 1 and 53938 DF, p-value: < 2.2e-16
Langkah selanjutnya adalah dengan mencoba memvisualisasikan model dengan sedikit perubahan.
# Membuat objek baru untuk visualisasi model menggunakan data_grid()
grid <- diamonds2 %>%
data_grid(carat) %>%
mutate(lcarat= log2(carat)) %>%
add_predictions(mod_diamonds, "lprice") %>%
mutate(price = 2 ^ lprice)
# lihat hasilnya
grid
## # A tibble: 273 x 4
## carat lcarat lprice price
## <dbl> <dbl> <dbl> <dbl>
## 1 0.2 -2.32 8.30 315.
## 2 0.21 -2.25 8.42 341.
## 3 0.22 -2.18 8.53 369.
## 4 0.23 -2.12 8.64 398.
## 5 0.24 -2.06 8.74 427.
## 6 0.25 -2 8.84 457.
## 7 0.26 -1.94 8.93 488.
## 8 0.27 -1.89 9.02 520.
## 9 0.28 -1.84 9.11 553.
## 10 0.29 -1.79 9.20 587.
## # ... with 263 more rows
Sekarang, kita visualisasikan model di atas
ggplot(diamonds2, aes(carat, price)) +
geom_hex(bins = 50) +
geom_line(data = grid, color ="red", size = 1)
Perhatikan bahwa dalam model di atas, kita tidak lagi menggunakan prediksi linier pada variabel ‘price’ (karena setelah mendapat nilai prediksi lprice, kita mengubahnya lagi menjadi ‘price’)
Apakah transformasi ini berdampak ? Kita bisa cek pada residualnya.
# Cek residual
diamonds2 <- diamonds2 %>% add_residuals(mod_diamonds, "lresid")
diamonds2 <- diamonds2 %>% mutate(resid = 2^lresid)
diamonds2
## # A tibble: 53,940 x 14
## carat cut color clarity depth table price x y z lprice lcarat
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43 8.35 -2.12
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31 8.35 -2.25
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31 8.35 -2.12
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63 8.38 -1.79
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75 8.39 -1.69
## 6 0.24 Very G~ J VVS2 62.8 57 336 3.94 3.96 2.48 8.39 -2.06
## 7 0.24 Very G~ I VVS1 62.3 57 336 3.95 3.98 2.47 8.39 -2.06
## 8 0.26 Very G~ H SI1 61.9 55 337 4.07 4.11 2.53 8.40 -1.94
## 9 0.22 Fair E VS2 65.1 61 337 3.87 3.78 2.49 8.40 -2.18
## 10 0.23 Very G~ H VS1 59.4 61 338 4 4.05 2.39 8.40 -2.12
## # ... with 53,930 more rows, and 2 more variables: lresid <dbl>, resid <dbl>
Untuk lebih mudah, kita lihat visualnya
ggplot(diamonds2, aes(lcarat, resid)) + geom_hex(bins = 50)
ggplot(diamonds2, aes(lcarat, lresid)) + geom_hex(bins = 50)
Kita bisa menggunakan ‘lresid’ pada attribute diamonds lain
ggplot(diamonds2, aes(cut, lresid)) + geom_boxplot()