model_2
Galih Dwika Putra R
1/22/2022
Model Building
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.
Diamonds
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.3modes(diamonds$depth)## [1] 62modes(diamonds$table)## [1] 56modes(diamonds$price)## [1] 605modes(diamonds$x)## [1] 4.37modes(diamonds$y)## [1] 4.34modes(diamonds$z)## [1] 2.7Di 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 rowsApa 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.676summary(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-16Langkah 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 rowsSekarang, 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()
Many Models
Pada bagian ini, kita akan mempelajari tiga ide yang membantu untuk mengolah model big data. * Menggunakan banyak model sederhana untuk memahami dataset yang kompleks * Menggunakan list-columns untuk menyimpan struktur data abritrari dalam data frame. Sebagai contoh, list-column akan memungkinkan kita memiliki sebuah kolom yang memiliki model linier * Menggunakan package ‘broom’ untuk mengubah model menjadi data yang tidy. Ini adalah teknik yang sangat berguna untuk bekerja dengan banyak model.
Kita akan belajar bagaimana untuk : * Kenapa menaruh list dalam data frame bisa dianggap valid * 3 cara utama untuk membuat list-column * Cara mengembalikan list-column menjadi atomic vector * dan menata data menggunkan broom
Untuk mempelajari bagian ini, package awal yang dibutuhkan adalah
library(tidyverse)
library(modelr)Sebagai awalan, kita akan mulai dengan menggunakan dataset bernama ‘gapminder’.
library(gapminder)## Warning: package 'gapminder' was built under R version 4.1.2gapminder## # A tibble: 1,704 x 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Afghanistan Asia 1952 28.8 8425333 779.
## 2 Afghanistan Asia 1957 30.3 9240934 821.
## 3 Afghanistan Asia 1962 32.0 10267083 853.
## 4 Afghanistan Asia 1967 34.0 11537966 836.
## 5 Afghanistan Asia 1972 36.1 13079460 740.
## 6 Afghanistan Asia 1977 38.4 14880372 786.
## 7 Afghanistan Asia 1982 39.9 12881816 978.
## 8 Afghanistan Asia 1987 40.8 13867957 852.
## 9 Afghanistan Asia 1992 41.7 16317921 649.
## 10 Afghanistan Asia 1997 41.8 22227415 635.
## # ... with 1,694 more rowsMari kita buat pertanyaan riset terlebih dahulu : * Bagaimana ekspektasi usia (lifeExp) berubah sepanjang waktu (year) untuk masing-masing negara(country)?
Seperti biasa, praktik termudah adalah dengan membuat plot untuk menjawab pertanyaan tersebut
gapminder %>%
ggplot(aes(year, lifeExp, group = country)) +
geom_line(alpha = 1/3)
Pada data di atas,overall, lifeExp mengalami kenaikan sepanjang waktu. Namun dapat kita lihat juga beberapa negara yang memiliki trend berbeda dari kebanyakan. Pertanyaan : Bagaimana cara agar kita bisa melihat pertumbuhan negara-negara seperti itu?
Jika kita hanya membahas satu negara, tentu penyusunannya lebih masuk akal.
nz <- filter(gapminder, country == "New Zealand")
nz %>% ggplot(aes(year, lifeExp)) +
geom_line() +
ggtitle("Full data = ")
nz_mod <- lm(lifeExp~year, data = nz)
nz %>% add_predictions(nz_mod) %>%
ggplot(aes(year, pred)) +
geom_line() +
ggtitle("Linier trend +")
nz %>%
add_residuals(nz_mod) %>%
ggplot(aes(year, resid)) +
geom_line() +
geom_hline(yintercept = 0, color = "white", size = 3) +
ggtitle("Remaining pattern")
Lalu, bagaimana cara kita menarik data negara-negara dengan trend berbeda secara keseluruhan? R memberikan tool berguna yang memungkinkan kita untuk menata data frame khusus di dalam data frame lain, dan itu menggunakan fungsi nest().
by_country <- gapminder %>% group_by(country, continent) %>%
nest()
by_country## # A tibble: 142 x 3
## # Groups: country, continent [142]
## country continent data
## <fct> <fct> <list>
## 1 Afghanistan Asia <tibble [12 x 4]>
## 2 Albania Europe <tibble [12 x 4]>
## 3 Algeria Africa <tibble [12 x 4]>
## 4 Angola Africa <tibble [12 x 4]>
## 5 Argentina Americas <tibble [12 x 4]>
## 6 Australia Oceania <tibble [12 x 4]>
## 7 Austria Europe <tibble [12 x 4]>
## 8 Bahrain Asia <tibble [12 x 4]>
## 9 Bangladesh Asia <tibble [12 x 4]>
## 10 Belgium Europe <tibble [12 x 4]>
## # ... with 132 more rowsDengan menggunakan nest(), akan muncul satu kolom baru bernama’data’ yang cukup unik karena berisikan list dari data-data frame lain.
Contoh, jika kita take out salah satu row, maka akan muncul dataframe berisi variabel lengkap untuk negara masing-masing.
by_country$data[1]## [[1]]
## # A tibble: 12 x 4
## year lifeExp pop gdpPercap
## <int> <dbl> <int> <dbl>
## 1 1952 28.8 8425333 779.
## 2 1957 30.3 9240934 821.
## 3 1962 32.0 10267083 853.
## 4 1967 34.0 11537966 836.
## 5 1972 36.1 13079460 740.
## 6 1977 38.4 14880372 786.
## 7 1982 39.9 12881816 978.
## 8 1987 40.8 13867957 852.
## 9 1992 41.7 16317921 649.
## 10 1997 41.8 22227415 635.
## 11 2002 42.1 25268405 727.
## 12 2007 43.8 31889923 975.Lihat?
Setiap elemen dalam kolom ‘data’, adalah dataframe lain yang berisikan seluruh data lengkap masing-masing negara. Apa dampaknya? Dengan dataframe masing-masing negara kini sudah tersusun dalam bentuk list di dalam variabel ‘data’, maka kita bisa menerapkan hal yang sama pada list pada umumnya.
Yep… fungsi map; fitur yang menerapkan satu fungsi pada seluruh elemen list secara simultan. Inilah yang disebut dengan ‘nested dataframe’.
Sekarang, saatnya kita melakukan fitting model. Buat sebuah fungsi yang berisikan model fitting yang kita inginkan.
country_model <- function(df){
lm(lifeExp~year, data = df)
}Lalu, kita terapkan fungsi map. Kita bisa membuat objek baru untuk menyimpan hasil model, atau menjadikannya kolom baru dalam dataframe by_country.
by_country <- by_country %>%
mutate(model = map(data, country_model))
by_country## # A tibble: 142 x 4
## # Groups: country, continent [142]
## country continent data model
## <fct> <fct> <list> <list>
## 1 Afghanistan Asia <tibble [12 x 4]> <lm>
## 2 Albania Europe <tibble [12 x 4]> <lm>
## 3 Algeria Africa <tibble [12 x 4]> <lm>
## 4 Angola Africa <tibble [12 x 4]> <lm>
## 5 Argentina Americas <tibble [12 x 4]> <lm>
## 6 Australia Oceania <tibble [12 x 4]> <lm>
## 7 Austria Europe <tibble [12 x 4]> <lm>
## 8 Bahrain Asia <tibble [12 x 4]> <lm>
## 9 Bangladesh Asia <tibble [12 x 4]> <lm>
## 10 Belgium Europe <tibble [12 x 4]> <lm>
## # ... with 132 more rowsCoba kita cek,
by_country$model[1]## [[1]]
##
## Call:
## lm(formula = lifeExp ~ year, data = df)
##
## Coefficients:
## (Intercept) year
## -507.5343 0.2753Seperti yang kita lihat, saat kita call salah satu elemen model, akan muncul hasil dari model fitting tersebut.
Menjadikan hasil estimasi model sebagai kolom dari dataframe by_country merupakan praktik yang dianjukan, karena untuk menjaga keakuratan data. Karena saat terpisah, dan kita melakukan manipulasi pada salah satu vector, hasil yang ada akan tertukar dan kita memperolah prediksi yang salah.
Sekarang, setelah selesai dengan model yg berisi model fitting masing-masing negara dalam kolom ‘data’, bagaimana kita menghitung residualnya.
Seperti yang kita tahu, ‘add_residuals(data, function)’ Nah dalam dataframe ini, kita akan menerapkan mapping menggunakan dua vector, yaitu ‘data’ yang berisi data frame masing-masing negara, dan ‘model’ yang menjadi vector fungsi.
by_country <- by_country %>%
mutate(resids = map2(data,model,add_residuals))
by_country## # A tibble: 142 x 5
## # Groups: country, continent [142]
## country continent data model resids
## <fct> <fct> <list> <list> <list>
## 1 Afghanistan Asia <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 2 Albania Europe <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 3 Algeria Africa <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 4 Angola Africa <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 5 Argentina Americas <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 6 Australia Oceania <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 7 Austria Europe <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 8 Bahrain Asia <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 9 Bangladesh Asia <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## 10 Belgium Europe <tibble [12 x 4]> <lm> <tibble [12 x 5]>
## # ... with 132 more rows# kita cek
by_country$resids[1]## [[1]]
## # A tibble: 12 x 5
## year lifeExp pop gdpPercap resid
## <int> <dbl> <int> <dbl> <dbl>
## 1 1952 28.8 8425333 779. -1.11
## 2 1957 30.3 9240934 821. -0.952
## 3 1962 32.0 10267083 853. -0.664
## 4 1967 34.0 11537966 836. -0.0172
## 5 1972 36.1 13079460 740. 0.674
## 6 1977 38.4 14880372 786. 1.65
## 7 1982 39.9 12881816 978. 1.69
## 8 1987 40.8 13867957 852. 1.28
## 9 1992 41.7 16317921 649. 0.754
## 10 1997 41.8 22227415 635. -0.534
## 11 2002 42.1 25268405 727. -1.54
## 12 2007 43.8 31889923 975. -1.22Oke… lalu sekarang bagaimana cara memvisualkan prediksi dari masing-masing negara tersebut? Akan sulit jika menggunakan ggplot pada dataframe by_country begitu saja yang terdiri dari data struktur list-coloumn.
Solusi paling efisien adalah dengan membuka nest data frame khusus data ‘resids’; caranya dengan menggunakan function unnest()
resids <- unnest(by_country,resids)
resids## # A tibble: 1,704 x 9
## # Groups: country, continent [142]
## country continent data model year lifeExp pop gdpPercap resid
## <fct> <fct> <list> <list> <int> <dbl> <int> <dbl> <dbl>
## 1 Afghanistan Asia <tibble ~ <lm> 1952 28.8 8.43e6 779. -1.11
## 2 Afghanistan Asia <tibble ~ <lm> 1957 30.3 9.24e6 821. -0.952
## 3 Afghanistan Asia <tibble ~ <lm> 1962 32.0 1.03e7 853. -0.664
## 4 Afghanistan Asia <tibble ~ <lm> 1967 34.0 1.15e7 836. -0.0172
## 5 Afghanistan Asia <tibble ~ <lm> 1972 36.1 1.31e7 740. 0.674
## 6 Afghanistan Asia <tibble ~ <lm> 1977 38.4 1.49e7 786. 1.65
## 7 Afghanistan Asia <tibble ~ <lm> 1982 39.9 1.29e7 978. 1.69
## 8 Afghanistan Asia <tibble ~ <lm> 1987 40.8 1.39e7 852. 1.28
## 9 Afghanistan Asia <tibble ~ <lm> 1992 41.7 1.63e7 649. 0.754
## 10 Afghanistan Asia <tibble ~ <lm> 1997 41.8 2.22e7 635. -0.534
## # ... with 1,694 more rowsSekarang, kita kembali memiliki dataframe yang benar… Kita bisa melakukan visual
resids %>%
ggplot(aes(year, resid)) +
geom_line(aes(group = country), alpha = 1/3) +
facet_wrap(~continent)
Pada plot residual per continent di atas, kita bisa melihat adanya yang tidak beres di Afrika. Untuk mengetahui yang mana yang salah, kita bisa menggunakan ukuran kualitas model (R square, adj.R Square, p.value, dll)
Untuk mencari ukuran kualitas model dengan memanfaatkan tool dari broom,kita bisa menggunakan broom::glance
broom::glance(nz_mod)## # A tibble: 1 x 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.954 0.949 0.804 205. 0.0000000541 1 -13.3 32.6 34.1
## # ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Nah, skg kita terapkan ‘glance’ pada dataframe by_country untuk mencari kualitas model per negara.
glance <- by_country %>%
mutate(glance = map(model,broom::glance)) %>%
unnest(glance, .drop = TRUE)## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.glance## # A tibble: 142 x 17
## # Groups: country, continent [142]
## country continent data model resids r.squared adj.r.squared sigma statistic
## <fct> <fct> <lis> <lis> <list> <dbl> <dbl> <dbl> <dbl>
## 1 Afghani~ Asia <tib~ <lm> <tibb~ 0.948 0.942 1.22 181.
## 2 Albania Europe <tib~ <lm> <tibb~ 0.911 0.902 1.98 102.
## 3 Algeria Africa <tib~ <lm> <tibb~ 0.985 0.984 1.32 662.
## 4 Angola Africa <tib~ <lm> <tibb~ 0.888 0.877 1.41 79.1
## 5 Argenti~ Americas <tib~ <lm> <tibb~ 0.996 0.995 0.292 2246.
## 6 Austral~ Oceania <tib~ <lm> <tibb~ 0.980 0.978 0.621 481.
## 7 Austria Europe <tib~ <lm> <tibb~ 0.992 0.991 0.407 1261.
## 8 Bahrain Asia <tib~ <lm> <tibb~ 0.967 0.963 1.64 291.
## 9 Banglad~ Asia <tib~ <lm> <tibb~ 0.989 0.988 0.977 930.
## 10 Belgium Europe <tib~ <lm> <tibb~ 0.995 0.994 0.293 1822.
## # ... with 132 more rows, and 8 more variables: p.value <dbl>, df <dbl>,
## # logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>Dengan menggunakan broom::glance sekarang kita punya parameter kualitas model untuk masing-masing negara yang tersusun menjadi dataframe. Struktur ini memungkinkan kita untuk melihat dan memfilter mana negara yang tidak sesuai dengan model.
glance %>% arrange(r.squared)## # A tibble: 142 x 17
## # Groups: country, continent [142]
## country continent data model resids r.squared adj.r.squared sigma statistic
## <fct> <fct> <lis> <lis> <list> <dbl> <dbl> <dbl> <dbl>
## 1 Rwanda Africa <tib~ <lm> <tibb~ 0.0172 -0.0811 6.56 0.175
## 2 Botswana Africa <tib~ <lm> <tibb~ 0.0340 -0.0626 6.11 0.352
## 3 Zimbabwe Africa <tib~ <lm> <tibb~ 0.0562 -0.0381 7.21 0.596
## 4 Zambia Africa <tib~ <lm> <tibb~ 0.0598 -0.0342 4.53 0.636
## 5 Swazila~ Africa <tib~ <lm> <tibb~ 0.0682 -0.0250 6.64 0.732
## 6 Lesotho Africa <tib~ <lm> <tibb~ 0.0849 -0.00666 5.93 0.927
## 7 Cote d'~ Africa <tib~ <lm> <tibb~ 0.283 0.212 3.93 3.95
## 8 South A~ Africa <tib~ <lm> <tibb~ 0.312 0.244 4.74 4.54
## 9 Uganda Africa <tib~ <lm> <tibb~ 0.342 0.276 3.19 5.20
## 10 Congo, ~ Africa <tib~ <lm> <tibb~ 0.348 0.283 2.43 5.34
## # ... with 132 more rows, and 8 more variables: p.value <dbl>, df <dbl>,
## # logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>, df.residual <int>,
## # nobs <int>Terlihat pada dataframe di atas, menunjukkan bahwa seluruh negara di benua Africa tidak sesuai dengan model yang kita buat.
Tiga alat dalam broom
- broom::glance(model) —> menghasilkan baris untuk setiap model, memberikan model summary.
- broom::tidy(model) —> Menghasilkan satu baris untuk setiap koefisien model. Setiap kolom memberi infomrasi soal estimasi dan variabilitinya.
- broom::augment(model, data) –> Mengerluarkan satu baris untuk tiap baris dalam data. menambah value extra seperti residual, dan influence statistic.