Model Basics with modelr

Memuat Package yang digunakan

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

Categorical Variables

Sebelumnya variabelya bersifat kontinu, sekarang jika variabelnya bersifat ktegorik.

Data yang digunakan adalah ini

sim2
## # A tibble: 40 x 2
##    x          y
##    <chr>  <dbl>
##  1 a      1.94 
##  2 a      1.18 
##  3 a      1.24 
##  4 a      2.62 
##  5 a      1.11 
##  6 a      0.866
##  7 a     -0.910
##  8 a      0.721
##  9 a      0.687
## 10 a      2.07 
## # ... with 30 more rows

prosesnya sama seperti yang sebelumnya(yang umum).

Berikut merupakan visualisasi dari data sim2.

ggplot(sim2) + 
  geom_point(aes(x, y))

Model dari sim2.

mod2 <- lm(y ~ x, data = sim2)
mod2
## 
## Call:
## lm(formula = y ~ x, data = sim2)
## 
## Coefficients:
## (Intercept)           xb           xc           xd  
##      1.1522       6.9639       4.9750       0.7588

Hasil prediksi berdasarkan model yang didapat.

grid <- sim2 %>%  data_grid(x) %>%  add_predictions(mod2) 
grid 
## # A tibble: 4 x 2
##   x      pred
##   <chr> <dbl>
## 1 a      1.15
## 2 b      8.12
## 3 c      6.13
## 4 d      1.91

Visualisasi data sim2 dengan hasil prediksi berdasarkan model yang didapat.

ggplot(sim2, aes(x)) +  
  geom_point(aes(y = y)) + 
  geom_point(data = grid,    
             aes(y = pred), 
             color = "red",  
             size = 4  )

Interaksi kontinu dan kategorik

Data yang digunakan.

sim3
## # A tibble: 120 x 5
##       x1 x2      rep      y    sd
##    <int> <fct> <int>  <dbl> <dbl>
##  1     1 a         1 -0.571     2
##  2     1 a         2  1.18      2
##  3     1 a         3  2.24      2
##  4     1 b         1  7.44      2
##  5     1 b         2  8.52      2
##  6     1 b         3  7.72      2
##  7     1 c         1  6.51      2
##  8     1 c         2  5.79      2
##  9     1 c         3  6.07      2
## 10     1 d         1  2.11      2
## # ... with 110 more rows

Berikut merupakan visualisasi dari data sim3.

ggplot(sim3, aes(x1, y)) +
  geom_point(aes(color = x2))

Model dari sim3 karena sebuah interaksi variabel kontinu dan kategorik, jadi ada dua kemungkinan yang mungkin terjadi, yaitu:

mod1 <- lm(y ~ x1 + x2, data = sim3) 
mod1
## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim3)
## 
## Coefficients:
## (Intercept)           x1          x2b          x2c          x2d  
##      1.8717      -0.1967       2.8878       4.8057       2.3596
mod2 <- lm(y ~ x1 * x2, data = sim3)
mod2
## 
## Call:
## lm(formula = y ~ x1 * x2, data = sim3)
## 
## Coefficients:
## (Intercept)           x1          x2b          x2c          x2d       x1:x2b  
##     1.30124     -0.09302      7.06938      4.43090      0.83455     -0.76029  
##      x1:x2c       x1:x2d  
##     0.06815      0.27728

Hasil prediksi berdasarkan kedua model yang didapat.

grid <- sim3 %>%  
  data_grid(x1, x2) %>%  
  gather_predictions(mod1, mod2)
grid 
## # A tibble: 80 x 4
##    model    x1 x2     pred
##    <chr> <int> <fct> <dbl>
##  1 mod1      1 a      1.67
##  2 mod1      1 b      4.56
##  3 mod1      1 c      6.48
##  4 mod1      1 d      4.03
##  5 mod1      2 a      1.48
##  6 mod1      2 b      4.37
##  7 mod1      2 c      6.28
##  8 mod1      2 d      3.84
##  9 mod1      3 a      1.28
## 10 mod1      3 b      4.17
## # ... with 70 more rows

Visualisasi data sim3 dengan hasil prediksi berdasarkan kedua model yang didapat.

ggplot(sim3, aes(x1, y, color = x2)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) + 
  facet_wrap(~ model)

Residual dari kedua model. Cari model dengan residual yang cenderung berkumpul di nol atau pilih model dengan residual paling kecil.

sim3 <- sim3 %>%  gather_residuals(mod1, mod2)
ggplot(sim3, aes(x1, resid, color = x2)) +
  geom_point() +  
  facet_grid(model ~ x2)

Dari dua model diatas didapat model paling bagus adalah model2 karena residualnya paling kecil.

Interaksi dua variabel kontinu

Data yang digunakan.

sim4
## # A tibble: 300 x 4
##       x1     x2   rep       y
##    <dbl>  <dbl> <int>   <dbl>
##  1    -1 -1         1  4.25  
##  2    -1 -1         2  1.21  
##  3    -1 -1         3  0.353 
##  4    -1 -0.778     1 -0.0467
##  5    -1 -0.778     2  4.64  
##  6    -1 -0.778     3  1.38  
##  7    -1 -0.556     1  0.975 
##  8    -1 -0.556     2  2.50  
##  9    -1 -0.556     3  2.70  
## 10    -1 -0.333     1  0.558 
## # ... with 290 more rows

Berikut merupakan visualisasi dari data sim4.

ggplot(sim4, aes(x1, x2)) +
  geom_tile(aes(fill = y))

Model dari sim4 karena sebuah interaksi dua variabel kontinu, jadi ada dua kemungkinan yang mungkin terjadi, yaitu:

mod1 <- lm(y ~ x1 + x2, data = sim4)
mod1
## 
## Call:
## lm(formula = y ~ x1 + x2, data = sim4)
## 
## Coefficients:
## (Intercept)           x1           x2  
##     0.03546      1.82167     -2.78252
mod2 <- lm(y ~ x1 * x2, data = sim4)
mod2
## 
## Call:
## lm(formula = y ~ x1 * x2, data = sim4)
## 
## Coefficients:
## (Intercept)           x1           x2        x1:x2  
##     0.03546      1.82167     -2.78252      0.95228

Hasil prediksi berdasarkan kedua model yang didapat.

grid <- sim4 %>%  
  data_grid(x1 = seq_range(x1, 5), 
            x2 = seq_range(x2, 5)  
            ) %>%  gather_predictions(mod1, mod2) 
grid 
## # A tibble: 50 x 4
##    model    x1    x2   pred
##    <chr> <dbl> <dbl>  <dbl>
##  1 mod1   -1    -1    0.996
##  2 mod1   -1    -0.5 -0.395
##  3 mod1   -1     0   -1.79 
##  4 mod1   -1     0.5 -3.18 
##  5 mod1   -1     1   -4.57 
##  6 mod1   -0.5  -1    1.91 
##  7 mod1   -0.5  -0.5  0.516
##  8 mod1   -0.5   0   -0.875
##  9 mod1   -0.5   0.5 -2.27 
## 10 mod1   -0.5   1   -3.66 
## # ... with 40 more rows

Pada data grid menggunakan seq_range(). Karena dari pada menggunakan setiap value pada variabel, akan lebih berguna pada kasus umum. Disini menggunakan 5 angka diantara maksimal dan minimal variabel x1 dan x2.

Macam-macam seq_range yang lain:

x <- c(0, 1) 
seq_range(x, n = 5)
## [1] 0.00 0.25 0.50 0.75 1.00
seq_range(x, n = 5, trim = 0.10) 
## [1] 0.050 0.275 0.500 0.725 0.950
seq_range(x, n = 5, trim = 0.25) 
## [1] 0.1250 0.3125 0.5000 0.6875 0.8750
seq_range(x, n = 5, trim = 0.50) 
## [1] 0.250 0.375 0.500 0.625 0.750
seq_range(x, n = 5)
## [1] 0.00 0.25 0.50 0.75 1.00
seq_range(x, n = 5, expand = 0.10)
## [1] -0.050  0.225  0.500  0.775  1.050
seq_range(x, n = 5, expand = 0.25) 
## [1] -0.1250  0.1875  0.5000  0.8125  1.1250
seq_range(x, n = 5, expand = 0.50) 
## [1] -0.250  0.125  0.500  0.875  1.250

Visualisasi data sim4 dengan hasil prediksi berdasarkan kedua model yang didapat pada seq_range.

ggplot(grid, aes(x1, pred, color = x2, group = x2)) +
  geom_line() + 
  facet_wrap(~ model) 

ggplot(grid, aes(x2, pred, color = x1, group = x1)) + 
  geom_line() + 
  facet_wrap(~ model)

Residual dari kedua model. Cari model dengan residual yang cenderung berkumpul di nol atau pilih model dengan residual paling kecil.

sim4 <- sim4 %>%  gather_residuals(mod1, mod2)

ggplot(sim4, aes(x1, resid, color = x2, group=x2)) +
   geom_point() + 
  facet_wrap(~ model)

ggplot(sim4, aes(x1, resid, group=model, color= model)) +
   geom_point() + 
  facet_wrap(~ x2)

a<-sim4%>%
  group_by(model)%>%
  summarise(mean=sum(resid))
a
## # A tibble: 2 x 2
##   model     mean
##   <chr>    <dbl>
## 1 mod1  4.23e-14
## 2 mod2  3.00e-14
a[1,2]>a[2,2]
##      mean
## [1,] TRUE

Dari dua model diatas didapat model paling bagus adalah model2 karena residualnya paling kecil (kayanya), menurut gambar mah gak jelas tapi menurut jumlah dari residual model 2 memiliki residual lebih kecil dari pada model 1.

TRANSFORMATIONS

Untuk data yang tidak linear dapat menggunakan I(). Untuk melihat persamaan yang dibentuk dari fitting lm dapat menggunakan model_matrix().

contoh:

df <- tribble(  ~y, ~x,   1,  1,   2,  2,   3,  3 ) 
df
## # A tibble: 3 x 2
##       y     x
##   <dbl> <dbl>
## 1     1     1
## 2     2     2
## 3     3     3
model_matrix(df, y ~  x) 
## # A tibble: 3 x 2
##   `(Intercept)`     x
##           <dbl> <dbl>
## 1             1     1
## 2             1     2
## 3             1     3
model_matrix(df, y ~ I(x^2) + x)
## # A tibble: 3 x 3
##   `(Intercept)` `I(x^2)`     x
##           <dbl>    <dbl> <dbl>
## 1             1        1     1
## 2             1        4     2
## 3             1        9     3
model_matrix(df, y ~ I(x^3)+I(x^2) + x)
## # A tibble: 3 x 4
##   `(Intercept)` `I(x^3)` `I(x^2)`     x
##           <dbl>    <dbl>    <dbl> <dbl>
## 1             1        1        1     1
## 2             1        8        4     2
## 3             1       27        9     3

pada fungsi nonlinear, dapat menggunakan pendekatan dengan polynomials atau dengan natural spline.

model_matrix(df, y ~ poly(x, 2))
## # A tibble: 3 x 3
##   `(Intercept)` `poly(x, 2)1` `poly(x, 2)2`
##           <dbl>         <dbl>         <dbl>
## 1             1     -7.07e- 1         0.408
## 2             1     -7.85e-17        -0.816
## 3             1      7.07e- 1         0.408
model_matrix(df, y ~ ns(x, 2)) 
## # A tibble: 3 x 3
##   `(Intercept)` `ns(x, 2)1` `ns(x, 2)2`
##           <dbl>       <dbl>       <dbl>
## 1             1       0           0    
## 2             1       0.566      -0.211
## 3             1       0.344       0.771

Berikut data yang digunakan, dapat dari bukunya.

sim5 <- tibble(x = seq(0, 3.5 * pi, length = 50),
               y = 4 * sin(x) + 
                 rnorm(length(x)) )
sim5
## # A tibble: 50 x 2
##        x       y
##    <dbl>   <dbl>
##  1 0     -0.325 
##  2 0.224  0.315 
##  3 0.449  0.0727
##  4 0.673  0.907 
##  5 0.898  6.01  
##  6 1.12   3.76  
##  7 1.35   3.56  
##  8 1.57   3.42  
##  9 1.80   4.52  
## 10 2.02   2.55  
## # ... with 40 more rows

Visualisasi data

ggplot(sim5, aes(x, y)) + 
  geom_point()

pendekatan model menggunakan natural spline.

mod1 <- lm(y ~ ns(x, 1), data = sim5) 
mod1
## 
## Call:
## lm(formula = y ~ ns(x, 1), data = sim5)
## 
## Coefficients:
## (Intercept)     ns(x, 1)  
##       1.789       -3.812
mod2 <- lm(y ~ ns(x, 2), data = sim5) 
mod2
## 
## Call:
## lm(formula = y ~ ns(x, 2), data = sim5)
## 
## Coefficients:
## (Intercept)    ns(x, 2)1    ns(x, 2)2  
##       1.906       -3.869       -2.238
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod3
## 
## Call:
## lm(formula = y ~ ns(x, 3), data = sim5)
## 
## Coefficients:
## (Intercept)    ns(x, 3)1    ns(x, 3)2    ns(x, 3)3  
##       4.793        2.911      -13.833       -3.720
mod4 <- lm(y ~ ns(x, 4), data = sim5) 
mod4
## 
## Call:
## lm(formula = y ~ ns(x, 4), data = sim5)
## 
## Coefficients:
## (Intercept)    ns(x, 4)1    ns(x, 4)2    ns(x, 4)3    ns(x, 4)4  
##       1.547       -9.291        7.440       -1.250      -10.435
mod5 <- lm(y ~ ns(x, 5), data = sim5)
mod5
## 
## Call:
## lm(formula = y ~ ns(x, 5), data = sim5)
## 
## Coefficients:
## (Intercept)    ns(x, 5)1    ns(x, 5)2    ns(x, 5)3    ns(x, 5)4    ns(x, 5)5  
##     -0.4997      -8.7192       4.2789       1.9943       4.9401     -10.7106

Visualisai model.

grid <- sim5 %>%  
  data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>% 
  gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")
ggplot(sim5, aes(x, y)) + 
  geom_point() + 
  geom_line(data = grid, color = "red") +
  facet_wrap(~ model)

MISSING VALUE

contoh:

df <- tribble(  ~x, ~y,  1, 2.2,  2, NA,  3, 3.5,  4, 8.3,  NA, 10 )

df
## # A tibble: 5 x 2
##       x     y
##   <dbl> <dbl>
## 1     1   2.2
## 2     2  NA  
## 3     3   3.5
## 4     4   8.3
## 5    NA  10
mod <- lm(y ~ x, data = df, na.action = na.exclude) 
mod
## 
## Call:
## lm(formula = y ~ x, data = df, na.action = na.exclude)
## 
## Coefficients:
## (Intercept)            x  
##     -0.2286       1.8357

untuk mengetahui berapa observasi yang digunakan, gunakan nobs.

nobs(mod)
## [1] 3

Bersambung