Rmodel
Arwan Zhagi
2/12/2021
Model Basics with rmodel
Memuat Package yang digunakan
library(tidyverse)
library(modelr)
options(na.action = na.warn)
library(nycflights13)
library(lubridate)
library(splines)
library(MASS)Model Building
Dataset
flights## # A tibble: 336,776 x 19
## year month day dep_time sched_dep_time dep_delay arr_time sched_arr_time
## <int> <int> <int> <int> <int> <dbl> <int> <int>
## 1 2013 1 1 517 515 2 830 819
## 2 2013 1 1 533 529 4 850 830
## 3 2013 1 1 542 540 2 923 850
## 4 2013 1 1 544 545 -1 1004 1022
## 5 2013 1 1 554 600 -6 812 837
## 6 2013 1 1 554 558 -4 740 728
## 7 2013 1 1 555 600 -5 913 854
## 8 2013 1 1 557 600 -3 709 723
## 9 2013 1 1 557 600 -3 838 846
## 10 2013 1 1 558 600 -2 753 745
## # ... with 336,766 more rows, and 11 more variables: arr_delay <dbl>,
## # carrier <chr>, flight <int>, tailnum <chr>, origin <chr>, dest <chr>,
## # air_time <dbl>, distance <dbl>, hour <dbl>, minute <dbl>, time_hour <dttm>Banyaknya penerbangan berdasarkan waktu (per hari).
daily <- flights %>%
mutate(date = make_date(year, month, day)) %>%
group_by(date) %>%
summarize(n = n())
daily## # A tibble: 365 x 2
## date n
## <date> <int>
## 1 2013-01-01 842
## 2 2013-01-02 943
## 3 2013-01-03 914
## 4 2013-01-04 915
## 5 2013-01-05 720
## 6 2013-01-06 832
## 7 2013-01-07 933
## 8 2013-01-08 899
## 9 2013-01-09 902
## 10 2013-01-10 932
## # ... with 355 more rowsVisualisasi data banyaknya penerbangan berdasarkan waktu.
ggplot(daily, aes(date, n)) +
geom_line()
Pada grafik diatas terlihat terdapat efek sebuah hari dalam seminggu yang sangat kuat dan mendominasi sehingga grafik berbentuk pola yang curam.
Distribusi penerbangan berdasarkan hari dalam seminggu.
daily <- daily %>%
mutate(wday = wday(date, label = TRUE)) %>%
mutate(month=month(date, label = T))
ggplot(daily, aes(wday, n)) +
geom_boxplot()
Pada grafik diatas terlihat bahwa pada weekend terdapat lebih sedikit penerbangan, hal ini karena banyak perjalan adalah untuk bisnis. Efek tersebut terlihat jelas pada hari sabtu.
Untuk menghilangkan efek kuat tersebut adalah menggunkan model tersebut.
Model
mod <- lm(n ~ wday, data = daily)
mod##
## Call:
## lm(formula = n ~ wday, data = daily)
##
## Coefficients:
## (Intercept) wday.L wday.Q wday.C wday^4 wday^5
## 922.595 -83.322 -155.111 -62.834 -80.126 -4.967
## wday^6
## -16.934Prediksi berdasarkan model.
grid <- daily %>% data_grid(wday) %>%
add_predictions(mod, "n")
grid## # A tibble: 7 x 2
## wday n
## <ord> <dbl>
## 1 Sun 891.
## 2 Mon 975.
## 3 Tue 951.
## 4 Wed 963.
## 5 Thu 966.
## 6 Fri 967.
## 7 Sat 745.Distribusi penerbangan berdasarkan hari daam seminggu dengan prediksinya.
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid, color = "red", size = 4)
Residual
daily <- daily %>%
add_residuals(mod)
daily## # A tibble: 365 x 5
## date n wday month resid
## <date> <int> <ord> <ord> <dbl>
## 1 2013-01-01 842 Tue Jan -109.
## 2 2013-01-02 943 Wed Jan -19.7
## 3 2013-01-03 914 Thu Jan -51.8
## 4 2013-01-04 915 Fri Jan -52.5
## 5 2013-01-05 720 Sat Jan -24.6
## 6 2013-01-06 832 Sun Jan -59.5
## 7 2013-01-07 933 Mon Jan -41.8
## 8 2013-01-08 899 Tue Jan -52.4
## 9 2013-01-09 902 Wed Jan -60.7
## 10 2013-01-10 932 Thu Jan -33.8
## # ... with 355 more rowsVisualisasi dari residual
daily %>%
ggplot(aes(date, resid)) +
geom_ref_line(h = 0) +
geom_line()
Pada visualisasi diatas terlihat penyimpangan pada perkiraan banyaknya penerbangan.
Pada visualisasi tersebut juga berguna karena telah menghapus banyak efek dari hari dalam seminggu sebelumnya. Namun disini masih tersisa beberapa pola curam.
Pada model tersebut terlihat gagal pada juni, karena terlihat pola curam.
Agar mudah dilihat, graph dibuat untuk tiap-tiap hari dalam seminggu.
ggplot(daily, aes(date, resid, color = wday)) +
geom_ref_line(h = 0) +
geom_line()
Data dengan residual lebih atau kurang dari 100.
daily %>% filter(resid < -100 )## # A tibble: 11 x 5
## date n wday month resid
## <date> <int> <ord> <ord> <dbl>
## 1 2013-01-01 842 Tue Jan -109.
## 2 2013-01-20 786 Sun Jan -105.
## 3 2013-05-26 729 Sun May -162.
## 4 2013-07-04 737 Thu Jul -229.
## 5 2013-07-05 822 Fri Jul -145.
## 6 2013-09-01 718 Sun Sep -173.
## 7 2013-11-28 634 Thu Nov -332.
## 8 2013-11-29 661 Fri Nov -306.
## 9 2013-12-24 761 Tue Dec -190.
## 10 2013-12-25 719 Wed Dec -244.
## 11 2013-12-31 776 Tue Dec -175.daily %>% filter(resid > 100 )## # A tibble: 1 x 5
## date n wday month resid
## <date> <int> <ord> <ord> <dbl>
## 1 2013-11-30 857 Sat Nov 112.Dilihat dari data diatas penyimpangan yang terjadi pada visualisasi sebelumnya adalah karena hari-hari tersebut merupakan hari-hari besar atau hari libur.
Berikut trend menurut pola.
daily %>% ggplot(aes(date, resid)) +
geom_ref_line(h = 0) +
geom_line(color = "grey50") +
geom_smooth(se = FALSE, span = 0.20) ## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Pada grafik sebelumnya terlihat bahwa hari sabtu (warna kuning), terlihat pola musiman. Pada bulan juni - september jumlah penerbangan lebih besar dari prediksi (expected) dan pada bulan september - Januari jumlah penerbangan lebih kecil dari predeksi (expected).
Berikut visualisasi hari sabtu setiap bulannya
daily %>%
filter(wday == "Sat") %>%
ggplot(aes(date, n)) +
geom_point() +
geom_line() +
scale_x_date( NULL,
date_breaks = "1 month",
date_labels = "%b" )
Visualisasi hari sabtu berdasarkan three school term.
term <- function(date) {
cut(date,
breaks = ymd(20130101, 20130605, 20130825, 20140101),
labels = c("spring", "summer", "fall")
)
}
daily <- daily %>%
mutate(term = term(date))
daily %>%
filter(wday == "Sat") %>%
ggplot(aes(date, n, color = term)) +
geom_point(alpha = 1/3) +
geom_line() +
scale_x_date( NULL,
date_breaks = "1 month",
date_labels = "%b" )
Distribusi penerbangan berdasarkan hari daam seminggu berdasarkan three school term sebelumnya.
daily %>%
ggplot(aes(wday, n, color = term)) +
geom_boxplot()
Pada garfik diatas terlihat terdapat variasi(perbedaan) yang signifikan antara masing-masing term. Jadi disini dipisahkan pada hari dalam seminggu berdasarkan term.
Berikut perbandingan model sebelumnya dengan model yang diperbaiki dan juga visualisasinya.
mod1 <- lm(n ~ wday, data = daily)
mod2 <- lm(n ~ wday * term, data = daily)
mod1##
## Call:
## lm(formula = n ~ wday, data = daily)
##
## Coefficients:
## (Intercept) wday.L wday.Q wday.C wday^4 wday^5
## 922.595 -83.322 -155.111 -62.834 -80.126 -4.967
## wday^6
## -16.934mod2##
## Call:
## lm(formula = n ~ wday * term, data = daily)
##
## Coefficients:
## (Intercept) wday.L wday.Q wday.C
## 912.8684 -71.4674 -161.0112 -65.6215
## wday^4 wday^5 wday^6 termsummer
## -82.0997 -1.3425 -12.4501 37.6922
## termfall wday.L:termsummer wday.Q:termsummer wday.C:termsummer
## 3.7598 -6.2215 26.4121 26.7451
## wday^4:termsummer wday^5:termsummer wday^6:termsummer wday.L:termfall
## 23.5933 -13.1343 -5.0026 -31.9048
## wday.Q:termfall wday.C:termfall wday^4:termfall wday^5:termfall
## -0.6959 -8.8349 -9.8570 -3.2107
## wday^6:termfall
## -8.2061daily %>%
gather_residuals(without_term = mod1, with_term = mod2) %>%
ggplot(aes(date, resid, color = model)) +
geom_line(alpha = 0.75)
Distribusi penerbangan berdasarkan hari daam seminggu dengan prediksinya pada model yang telah diperbaiki.
grid <- daily %>%
data_grid(wday, term) %>%
add_predictions(mod2, "n")
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid, color = "red") +
facet_wrap(~ term)
Seharusnya hasil prediksi mendekati rata-rata. Pada visualisasi diatas terlihat prediksi jauh dari rata-rata. Hal ini dikarenakan oleh adanya outlier yang cukup banyak. Masalah outlier dapat diatasi dengan menggunakan M estimator atau menggunakan perintah rlm pada package MASS.
mod3 <- MASS::rlm(n ~ wday * term, data = daily)
mod3## Call:
## rlm(formula = n ~ wday * term, data = daily)
## Converged in 17 iterations
##
## Coefficients:
## (Intercept) wday.L wday.Q wday.C
## 922.4857006 -79.3059001 -153.7810509 -67.7749553
## wday^4 wday^5 wday^6 termsummer
## -74.9294649 -6.4803220 -9.7257891 32.9796852
## termfall wday.L:termsummer wday.Q:termsummer wday.C:termsummer
## 1.9847024 9.5345903 12.3444851 17.0750064
## wday^4:termsummer wday^5:termsummer wday^6:termsummer wday.L:termfall
## 11.0812708 -3.5357217 0.3378038 -31.7324552
## wday.Q:termfall wday.C:termfall wday^4:termfall wday^5:termfall
## -33.1644078 -23.8270526 -22.5287983 -4.9156541
## wday^6:termfall
## -3.3798444
##
## Degrees of freedom: 365 total; 344 residual
## Scale estimate: 17.6summary(mod3)##
## Call: rlm(formula = n ~ wday * term, data = daily)
## Residuals:
## Min 1Q Median 3Q Max
## -347.537 -12.899 3.915 11.500 160.101
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 922.4857 1.5914 579.6839
## wday.L -79.3059 4.2202 -18.7920
## wday.Q -153.7811 4.2136 -36.4961
## wday.C -67.7750 4.2081 -16.1056
## wday^4 -74.9295 4.2229 -17.7437
## wday^5 -6.4803 4.1961 -1.5444
## wday^6 -9.7258 4.2011 -2.3151
## termsummer 32.9797 2.7178 12.1349
## termfall 1.9847 2.3615 0.8404
## wday.L:termsummer 9.5346 7.2110 1.3222
## wday.Q:termsummer 12.3445 7.1875 1.7175
## wday.C:termsummer 17.0750 7.2040 2.3702
## wday^4:termsummer 11.0813 7.1885 1.5415
## wday^5:termsummer -3.5357 7.1969 -0.4913
## wday^6:termsummer 0.3378 7.1550 0.0472
## wday.L:termfall -31.7325 6.2480 -5.0788
## wday.Q:termfall -33.1644 6.2524 -5.3043
## wday.C:termfall -23.8271 6.2399 -3.8185
## wday^4:termfall -22.5288 6.2606 -3.5985
## wday^5:termfall -4.9157 6.2318 -0.7888
## wday^6:termfall -3.3798 6.2550 -0.5403
##
## Residual standard error: 17.58 on 344 degrees of freedomgrid2 <- daily %>%
data_grid(wday, term) %>%
add_predictions(mod3, "n")
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid2, color = "red") +
facet_wrap(~ term)+
labs(title = "Model Sesudah Diperbaiki")
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid, color = "red") +
facet_wrap(~ term)+
labs(title = "Model Sebelum Diperbaiki")
Visualisasi residual pada model yang telah diperbaiki.
daily <- daily %>%
add_residuals(mod3, "resid")
daily %>%
ggplot(aes(date, resid)) +
geom_hline(yintercept = 0, size = 2, color = "blue") +
geom_line()
Pada grafik diatas terlihat bahwa model berhasil menghilangkan pola hari dalam seminggu (Residual mendaki nol). Berikut perbandingan antara model menggunakan lm dan rlm.
daily %>%
gather_residuals(menggunakan_lm = mod2, menggunakan_rlm = mod3) %>%
ggplot(aes(date, resid, color = model)) +
geom_hline(yintercept = 0, size = 2, color = "blue") +
geom_line()
APABILA DIBUAT MENJADI SUATU FUNCTION
Untuk mencari model. Untuk function wday gunakan wday2 atau yang lain selain wday, karena function wday sudah ada di library lubridate/tidyverse.
daily <- flights %>%
mutate(date = make_date(year, month, day)) %>%
group_by(date) %>%
summarize(n = n()) ## `summarise()` ungrouping output (override with `.groups` argument)term <- function(date) {
cut(date,
breaks = ymd(20130101, 20130605, 20130825, 20140101),
labels = c("spring", "summer", "fall")
)
}
wday2 <- function(x) {
wday(x, label = TRUE)
}
mod4 <- rlm(n ~ wday2(date) * term(date), data = daily)
mod4## Call:
## rlm(formula = n ~ wday2(date) * term(date), data = daily)
## Converged in 17 iterations
##
## Coefficients:
## (Intercept) wday2(date).L
## 922.4857006 -79.3059001
## wday2(date).Q wday2(date).C
## -153.7810509 -67.7749553
## wday2(date)^4 wday2(date)^5
## -74.9294649 -6.4803220
## wday2(date)^6 term(date)summer
## -9.7257891 32.9796852
## term(date)fall wday2(date).L:term(date)summer
## 1.9847024 9.5345903
## wday2(date).Q:term(date)summer wday2(date).C:term(date)summer
## 12.3444851 17.0750064
## wday2(date)^4:term(date)summer wday2(date)^5:term(date)summer
## 11.0812708 -3.5357217
## wday2(date)^6:term(date)summer wday2(date).L:term(date)fall
## 0.3378038 -31.7324552
## wday2(date).Q:term(date)fall wday2(date).C:term(date)fall
## -33.1644078 -23.8270526
## wday2(date)^4:term(date)fall wday2(date)^5:term(date)fall
## -22.5287983 -4.9156541
## wday2(date)^6:term(date)fall
## -3.3798444
##
## Degrees of freedom: 365 total; 344 residual
## Scale estimate: 17.6Untuk memasukan variabel pada dataset gunakan function:
compute_vars <- function(data) {
data %>%
mutate( term = term(date),
wday = wday(date, label = TRUE)
)
}
daily<-compute_vars(daily)
daily## # A tibble: 365 x 4
## date n term wday
## <date> <int> <fct> <ord>
## 1 2013-01-01 842 spring Tue
## 2 2013-01-02 943 spring Wed
## 3 2013-01-03 914 spring Thu
## 4 2013-01-04 915 spring Fri
## 5 2013-01-05 720 spring Sat
## 6 2013-01-06 832 spring Sun
## 7 2013-01-07 933 spring Mon
## 8 2013-01-08 899 spring Tue
## 9 2013-01-09 902 spring Wed
## 10 2013-01-10 932 spring Thu
## # ... with 355 more rowsSebelumnya model diperbaiki dengan pengetahuan tentang three term school. Sebagai alternatif dapat digunakan natural spline agar pengetahuan kita dapat masuk dalam model secara explisit.
mod <- MASS::rlm(n ~ wday * ns(date, 5), data = daily)
mod## Call:
## rlm(formula = n ~ wday * ns(date, 5), data = daily)
## Converged in 14 iterations
##
## Coefficients:
## (Intercept) wday.L wday.Q wday.C
## 839.1368183 -58.0966417 -177.5614443 -76.1745048
## wday^4 wday^5 wday^6 ns(date, 5)1
## -106.5456685 14.7851980 -18.4661221 90.3704238
## ns(date, 5)2 ns(date, 5)3 ns(date, 5)4 ns(date, 5)5
## 133.7193537 41.3816719 180.6645268 28.9177394
## wday.L:ns(date, 5)1 wday.Q:ns(date, 5)1 wday.C:ns(date, 5)1 wday^4:ns(date, 5)1
## -35.0804683 11.6485897 -6.8750934 23.7107143
## wday^5:ns(date, 5)1 wday^6:ns(date, 5)1 wday.L:ns(date, 5)2 wday.Q:ns(date, 5)2
## -21.0684733 6.6102716 19.8439518 66.9609288
## wday.C:ns(date, 5)2 wday^4:ns(date, 5)2 wday^5:ns(date, 5)2 wday^6:ns(date, 5)2
## 54.5911251 59.4655330 -21.2689663 8.7438786
## wday.L:ns(date, 5)3 wday.Q:ns(date, 5)3 wday.C:ns(date, 5)3 wday^4:ns(date, 5)3
## -104.5352071 -83.6926248 -77.0186820 -36.3447345
## wday^5:ns(date, 5)3 wday^6:ns(date, 5)3 wday.L:ns(date, 5)4 wday.Q:ns(date, 5)4
## -25.0391121 0.8459285 -29.5596025 58.6435656
## wday.C:ns(date, 5)4 wday^4:ns(date, 5)4 wday^5:ns(date, 5)4 wday^6:ns(date, 5)4
## 32.5922420 61.0040478 -41.3647943 13.1681565
## wday.L:ns(date, 5)5 wday.Q:ns(date, 5)5 wday.C:ns(date, 5)5 wday^4:ns(date, 5)5
## 39.3654330 59.7652461 39.0232787 22.0796539
## wday^5:ns(date, 5)5 wday^6:ns(date, 5)5
## 2.1498306 -4.9029910
##
## Degrees of freedom: 365 total; 323 residual
## Scale estimate: 10.7daily %>%
data_grid(wday, date = seq_range(date, n = 13)) %>%
add_predictions(mod)%>%
ggplot(aes(date, pred, color = wday)) +
geom_line() +
geom_point()