modelr

Model Basics with modelr

Memuat Package yang digunakan

library(tidyverse)
library(modelr) 
options(na.action = na.warn) 

Basic of How models work.

A Simple Model

Data yang digunakan adalah data Sim1, dengan data sebaga berikut:

sim1

## # A tibble: 30 x 2
##        x     y
##    <int> <dbl>
##  1     1  4.20
##  2     1  7.51
##  3     1  2.13
##  4     2  8.99
##  5     2 10.2 
##  6     2 11.3 
##  7     3  7.36
##  8     3 10.5 
##  9     3 10.5 
## 10     4 12.4 
## # ... with 20 more rows

Berikut merupakan visualisasi dari sim1 (scatter plot).

ggplot(sim1, aes(x, y)) + 
  geom_point()+
  scale_x_continuous(breaks = rep(1:10))

Berikut merupakan 250 pola/model yang mungkin dari data tersebut. Walaupun kalau dilihat pasti ketebak linear sih.

models <- tibble( 
  a1 = runif(250, -20, 40), 
  a2 = runif(250, -5, 5)
  )

ggplot(sim1, aes(x, y)) +  
  geom_abline(    
    aes(intercept = a1, slope = a2),  
    data = models, alpha = 1/4  ) +
  geom_point()

Kemudian diambil 10 pola/model terbaik (yang mempunyai jarak antara pola dengan data paling kecil). Warnanya semakin biru muda maka jaraknya semakin dekat.

model1 <- function(a, data) { 
  a[1] + data$x * a[2] 
} 


measure_distance <- function(mod, data){
  diff <- data$y - model1(mod, data)
  sqrt(mean(diff ^ 2))
} 


sim1_dist <- function(a1, a2) { 
  measure_distance(c(a1, a2), sim1) 
  }

models <- models %>%
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist)) 

ggplot(sim1, aes(x, y)) +  
  geom_point(size = 2, color = "grey30") + 
  geom_abline(aes(intercept = a1, slope = a2, color = -dist),
              data = filter(models, rank(dist) <= 10)  
              )

Berikut merupakan grafik versus (scatter plot) antara a1 (konstanta) dan a2 (koefisien x). Bulatan merah adalah menunjukkan a1 dan a2 10 pola/ model sebelumnya.

ggplot(models, aes(a1, a2)) + 
  geom_point(data = filter(models, rank(dist) <= 10),   
             size = 4, color = "red"  ) +  
  geom_point(aes(colour = -dist))

Kemudian dari Plot sebelumnya dibuat beberapa grid, dan didapat 10 plot terbaik (yang mempunyai jarak terdekat).

grid <- expand.grid(a1 = seq(-5, 20, length = 25), 
                    a2 = seq(1, 3, length = 25)  ) %>%  
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))


grid %>%  ggplot(aes(a1, a2)) +
  geom_point(data = filter(grid, rank(dist) <= 10),
             size = 4, colour = "red"  ) + 
  geom_point(aes(color = -dist))

Berikut plot 10 pola/model terbaik berdasarkan grafik sebelumnya.

ggplot(sim1, aes(x, y)) +
  geom_point(size = 2, color = "grey30") + 
  geom_abline(aes(intercept = a1, slope = a2, color = -dist),
              data = filter(grid, rank(dist) <= 10)  
              )

Dan berikut merupakan pola/model terbaik.

best <- optim(c(0, 0), measure_distance, data = sim1)

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, color = "grey30") + 
  geom_abline(intercept = best$par[1], slope = best$par[2]
              )

Dan berikut merupakan nilai a1 dan a2

best$par

## [1] 4.222248 2.051204

---

Padahal mah, karena terlihat linear jadi sebenernya proses diatas bisa satu baris.

sim1_mod<-lm(y~x, data = sim1)
coef(lm(y ~ x, data = sim1)
     )

## (Intercept)           x 
##    4.220822    2.051533

---

Keseluruhan proses sebelumnya dapat dilakukan dengan cara pendekatan, seperti berikut:

grid <- sim1 %>%  
  data_grid(x) 


grid <- grid %>% 
  add_predictions(sim1_mod) 


ggplot(sim1, aes(x)) +  geom_point(aes(y = y)) +
  geom_line(aes(y = pred),  
            data = grid,
            colour = "red",    
            size = 1  )

Proses diatas lebih umum karena dapat digunakan pada kasus yang lebih kompleks.

---

RESIDUAL

sim1 <- sim1 %>%  
  add_residuals(sim1_mod) 
sim1 

## # A tibble: 30 x 3
##        x     y    resid
##    <int> <dbl>    <dbl>
##  1     1  4.20 -2.07   
##  2     1  7.51  1.24   
##  3     1  2.13 -4.15   
##  4     2  8.99  0.665  
##  5     2 10.2   1.92   
##  6     2 11.3   2.97   
##  7     3  7.36 -3.02   
##  8     3 10.5   0.130  
##  9     3 10.5   0.136  
## 10     4 12.4   0.00763
## # ... with 20 more rows

ggplot(sim1, aes(resid)) +
  geom_freqpoly(binwidth = 0.5)

ggplot(sim1, aes(x, resid)) +
  geom_ref_line(h = 0) +  geom_point()

---

BERSAMBUNG …

---