TutoRial R for GLM V2.0

1 Kenapa menggunakan R?

“Out of all the back half of the Top 20 languages, R has shown the most consistent upwards movement over time.”-RedMonk Programming Language Rankings

“R is a programming language and free software environment for statistical computing and graphics supported by the R Foundation for Statistical Computing”-wikipedia


“RStudio is an integrated development environment (IDE) for R”-RStudio.com

2 Seputar R

• Dibuat oleh statistikawan pada tahun 1992
  
• Open Sources
  
• High-Level Language
  
• R Interpreted Language
  
• R Case Sensitive
  
• Salah satu bahasa pemrograman yang paling sering digunakan dalam bidang Data Science
  

---

3 Fungsi-fungsi Primitif pada R

3.1 Assignment variable

Pada R kita dapat menyimpan sebuah objek kedalam variable menggunakan simbol <- atau simbol =

variable_1 <- 3
variable_1 = 3

3.2 Tipe data sederhana

Disini kita menyimpan angka 1 kedalam variable_1, untuk mengecek tipe data yang disimpan, kita menggunakan fungsi typeof

typeof(variable_1)

[1] "double"

Tipe" value yang dapat disimpan adalah :
1. numeric, double (bilangan real)
2. integer (bilangan bulat)
3. logical (TRUE atau FALSE)
4. chr (karakter atau string, biasa diperlakukann jugasebagai factor)
5. Factor (kategorikal seperti : warna , jenis kelamin, merek, dll)

Kita dapat mengganti tipe value dengan fungsi as.integer,as.factor / factor,as.numeric,dll

as.character(variable_1)

[1] "3"

3.3 Conditional Expression

3.3.1 If Else

Fungsi if dalam R :

if (10 + 5 < 10){
print("berikut contoh if yang berhasil")
}else{
print("berikut contoh else yang berhasil")
}

[1] "berikut contoh else yang berhasil"

3.3.2 For dan While loop

Iterasi di R dilakukan dengan 2 cara, yakni while dan for

x <- 0
while (x < 10){
x <- x + 1
}
print(x)

[1] 10

rm(x)

Kita dapat membuat vektor integer terurut dengan cara :

-1:20

 [1] -1  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

for (x in 1:10){
  print(x)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10

Selain iterasi dalam sebuah selang integer, kita juga dapat melakukan iterasi terhadap komponen vektor

variable_2 <- c(2,3,5,7,9,11,13,17,19,23,29)
for (m in variable_2){
  print(m)
}

[1] 2
[1] 3
[1] 5
[1] 7
[1] 9
[1] 11
[1] 13
[1] 17
[1] 19
[1] 23
[1] 29

---

3.4 Struktur data dan slicing

Struktur data di R yang umum adalah Vektor, Matrix dan Data frame.

3.4.1 Vector

Vektor dapat dianalogikan sebagai array, sebuah vektor harus memiliki tipe data yang sama untuk setiap komponennya.

#inisiasi vektor
vektor_a <- c("satu","dua","tiga","empat")
print(vektor_a)

[1] "satu"  "dua"   "tiga"  "empat"

#Kita dapat menghapus variable menggunakan fungsi rm(nama_Variable)
rm(vektor_a) 

3.4.2 Matrix

Matriks merupakan struktur data yang memiliki baris dan kolom, setiap elemen pada Matriks harus memiliki tipe yang sama.

matriks_a <- matrix(c(1,2,3,4),ncol = 2) #Pendefinisian matrix
print(matriks_a)

     [,1] [,2]
[1,]    1    3
[2,]    2    4

rm(matriks_a)

3.4.3 Data Frame

Data frame dapat dianalogikan sebagai tabel di excel, data frame memiliki nama kolom dan setiap kolom.

df_a <- data.frame(X1 = c(18,12,33,17,16,20),X2=c("M","F","M","M","M","F"),X3 = c(1,1,0,0,1,0))
print(df_a)

3.5 Slicing

Slicing adalah proses mengambil sebagian anggota dari sebuah data, slicing yang dasar pada R menggunakan operator [ ].
Berikut contoh slicing pada vektor

var_2 <- 1:10
print(var_2[5]) #Mengambil element ke 5 pada vektor var_2, counting di R dimulai dari indeks 1, BUKAN 0

[1] 5

print(var_2[1:3]) #Mengambil elemen pertama hingga ketiga dari vektor var_2

[1] 1 2 3

print(var_2[c(1,3,7)]) #Mengambil elemen ke 1, 3 dan 7 dari var_2

[1] 1 3 7

print(var_2[c(-9,-10,-2)]) # Membuang elemen ke 9, 10 dan 2

[1] 1 3 4 5 6 7 8

Slicing pada matriks menggunakan sistem double indeks, indeks kosong setara dengan memilih seluruh index.
Berikut contoh slicing pada matriks

var_3 <- round(matrix(runif(5*4), ncol=5)*100)
print(var_3)

     [,1] [,2] [,3] [,4] [,5]
[1,]   42   42   32   98   92
[2,]   78   84   41   55   18
[3,]   44    3   49   16   64
[4,]   55   59    8   32   74

print(var_3[3,])

[1] 44  3 49 16 64

print(var_3[1:3,-2])

     [,1] [,2] [,3] [,4]
[1,]   42   32   98   92
[2,]   78   41   55   18
[3,]   44   49   16   64

Slicing pada data frame mirip dengan slicing pada matriks.
Kita dapat mengambil sebuah kolom dari dataframe menggunakan nama kolom tersebut dengan operator $

df_a$X2

[1] M F M M M F
Levels: F M

Untuk menghapus sebuah kolom dari data frame, masukkan NULL kedalam kolom yang ingin kita hapus

df_a$X2 <- NULL

---

4 Packages

Package adalah kumpulan fungsi atau data yang terkompilasi dalam 1 file. Umumnya package digunakan untuk mempermudah kita dalam melakukan beberapa pekerjaan.
Beberapa package yang terkenal :
- caret(Classification and Regresion Training, untuk membuat model machine learning)
- GGplot (Untuk melakukan visualisasi data)
- tidyr (untuk mempermudah pengolahan dataframe)

install.packages("tidyverse") #Untuk install packages, cukup satu kali dijalankan
library(tidyverse) #Untuk menggunakan package caret

Jika kita ingin menggunakan satu fungsi saja dari suatu package, lebih baik memanggil fungsi tersebut secara spesifik dari package daripada memanggil seluruh fungsi dari package. Contoh :

tidyr::spread()

Selain melalui perintah install.packages, kita juga dapat melakukan instalasi packages melalui opsi install packages

---

5 Input file

R dapat menerima data dalam berbagai format (CSV,Excel,JSON, hingga terhubung ke SQL) Tipe data yang paling umum digunakan adalah csv.
“Comma Separated Values atau CSV adalah suatu format data dalam basis data di mana setiap record dipisahkan dengan tanda koma (,) atau titik koma (;)” - Wikipedia.

data <- read.csv("data1.csv",sep=",") #Untuk file csv, sep(separator) adalah tanda yang digunakan untuk memisahkan tiap cell dalam data
str(data)

'data.frame':   27 obs. of  2 variables:
 $ A: num  77 91 88 29 91.5 77 80.5 86.5 92 86.5 ...
 $ B: num  83 50 76.5 79 73 81.5 76.5 85.5 82.5 82.5 ...

6 Statistika Deskriptif

6.1 Sari Numerik

untuk 6 fungsi dibawah, input yang diterima adalah vector

sum(data$A)     #Mengeluarkan penjumlahan elemen-elemen vector

[1] 2267

unique(data$A)  #Menghasilkan elemen-elemen vektor yang unik

 [1] 77.0 91.0 88.0 29.0 91.5 80.5 86.5 92.0 30.0 90.0 82.0 93.5 86.0 87.5 90.5

length(data$A)  #Menampilkan panjang vector

[1] 27

mean(data$A)    #Menampilkan rata-rata dari elemen vector

[1] 83.96296

median(data$A)  #Menampilkan median dari elemen vector

[1] 90

var(data$A)     #Menampilan variansi dari elemen vector

[1] 266.4793

sd(data$A)      #Menampilkan standar deviasi dari elemen vector 

[1] 16.3242

summary(data)   #Dapat menerima input berupa vector atau data frame, menghasilan beberapa sari numerik

       A               B        
 Min.   :29.00   Min.   :50.00  
 1st Qu.:86.25   1st Qu.:77.00  
 Median :90.00   Median :80.50  
 Mean   :83.96   Mean   :79.41  
 3rd Qu.:91.50   3rd Qu.:82.75  
 Max.   :93.50   Max.   :89.50  

6.2 Visualisasi Sederhana

hist(data$B,breaks = 10) #Histogram-Input : vector, breaks menyatakan banyaknya bar yang diinginkan

boxplot(data$A,data$B,col = "blue")#Box Plot-Input : bisa beberapa vector, col menyatakan warna

plot(data$A,data$B,main = "Scatter Plot",xlab = "kolom A",ylab = "Kolom B")#Scatter Plot-Input:1 atau 2 vector,main menyatakan judul dari plot

stem(data$B,scale = 1)#Stem Leaf plot - Input : vector,scale menyatakan ukuran plot

  The decimal point is 1 digit(s) to the right of the |

  5 | 0
  5 | 
  6 | 
  6 | 
  7 | 33
  7 | 7777799
  8 | 0011122233333
  8 | 566
  9 | 0

---

7 Regresi Linear

7.1 Analisis DataSet mtcars(Motor Trend Car Road Tests)

Dataset yang akan kita gunakan adalah mtcars.
Ketika menginput dataset, 3 hal pertama yang umumnya dilakukan adalah melihat head() untuk melihat beberapa record paling atas
str() untuk melihat struktur dari dataset
summary() untuk melihat sari numerik dari setiap variabel di dataset

Description

The data was extracted from the 1974 Motor Trend US magazine, and comprises fuel consumption and 10 aspects of automobile design and performance for 32 automobiles (1973-74 models).

A data frame with 32 observations on 11 (numeric) variables.
[, 1] mpg Miles/(US) gallon
[, 2] cyl Number of cylinders
[, 3] disp Displacement (cu.in.)
[, 4] hp Gross horsepower
[, 5] drat Rear axle ratio
[, 6] wt Weight (1000 lbs)
[, 7] qsec 1/4 mile time
[, 8] vs Engine (0 = V-shaped, 1 = straight)
[, 9] am Transmission (0 = automatic, 1 = manual)
[,10] gear Number of forward gears
[,11] carb Number of carburetors

attach(mtcars)
head(mtcars)

str(mtcars)

'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

summary(mtcars)

      mpg             cyl             disp             hp             drat      
 Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0   Min.   :2.760  
 1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5   1st Qu.:3.080  
 Median :19.20   Median :6.000   Median :196.3   Median :123.0   Median :3.695  
 Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7   Mean   :3.597  
 3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0   3rd Qu.:3.920  
 Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0   Max.   :4.930  
       wt             qsec             vs               am        
 Min.   :1.513   Min.   :14.50   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000   1st Qu.:0.0000  
 Median :3.325   Median :17.71   Median :0.0000   Median :0.0000  
 Mean   :3.217   Mean   :17.85   Mean   :0.4375   Mean   :0.4062  
 3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000   3rd Qu.:1.0000  
 Max.   :5.424   Max.   :22.90   Max.   :1.0000   Max.   :1.0000  
      gear            carb      
 Min.   :3.000   Min.   :1.000  
 1st Qu.:3.000   1st Qu.:2.000  
 Median :4.000   Median :2.000  
 Mean   :3.688   Mean   :2.812  
 3rd Qu.:4.000   3rd Qu.:4.000  
 Max.   :5.000   Max.   :8.000  

Ketika kita membangun model regresi linear, coba periksa terlebih dahulu korelasi antar variabel.2 hal yang ingin kita lihat adalah :
1. Korelasi setiap peubah terjadap response
2. adanya multikolienaritas antar prediktor

cor(mtcars,method = "spearman") #menghasilakn matrix korelasi

            mpg        cyl       disp         hp        drat         wt
mpg   1.0000000 -0.9108013 -0.9088824 -0.8946646  0.65145546 -0.8864220
cyl  -0.9108013  1.0000000  0.9276516  0.9017909 -0.67888119  0.8577282
disp -0.9088824  0.9276516  1.0000000  0.8510426 -0.68359210  0.8977064
hp   -0.8946646  0.9017909  0.8510426  1.0000000 -0.52012499  0.7746767
drat  0.6514555 -0.6788812 -0.6835921 -0.5201250  1.00000000 -0.7503904
wt   -0.8864220  0.8577282  0.8977064  0.7746767 -0.75039041  1.0000000
qsec  0.4669358 -0.5723509 -0.4597818 -0.6666060  0.09186863 -0.2254012
vs    0.7065968 -0.8137890 -0.7236643 -0.7515934  0.44745745 -0.5870162
am    0.5620057 -0.5220712 -0.6240677 -0.3623276  0.68657079 -0.7377126
gear  0.5427816 -0.5643105 -0.5944703 -0.3314016  0.74481617 -0.6761284
carb -0.6574976  0.5800680  0.5397781  0.7333794 -0.12522294  0.4998120
            qsec         vs          am       gear        carb
mpg   0.46693575  0.7065968  0.56200569  0.5427816 -0.65749764
cyl  -0.57235095 -0.8137890 -0.52207118 -0.5643105  0.58006798
disp -0.45978176 -0.7236643 -0.62406767 -0.5944703  0.53977806
hp   -0.66660602 -0.7515934 -0.36232756 -0.3314016  0.73337937
drat  0.09186863  0.4474575  0.68657079  0.7448162 -0.12522294
wt   -0.22540120 -0.5870162 -0.73771259 -0.6761284  0.49981205
qsec  1.00000000  0.7915715 -0.20333211 -0.1481997 -0.65871814
vs    0.79157148  1.0000000  0.16834512  0.2826617 -0.63369482
am   -0.20333211  0.1683451  1.00000000  0.8076880 -0.06436525
gear -0.14819967  0.2826617  0.80768800  1.0000000  0.11488698
carb -0.65871814 -0.6336948 -0.06436525  0.1148870  1.00000000

pairs(mtcars) #menghasilkan matrix scatter plot

Untuk beberapa section kedepan, kita akan fokus melihat hubungan disp dan mpg karena keduanya numerik dan memiliki korelasi yang cukup kuat.

plot(mtcars$disp,mtcars$mpg)

Variabel Response disini adalah mpg.
Disini data akan kita modelkan dengan GLM (residu diasumsikan berdistribusi normal)

---

7.2 Model Regresi Linear Sederhana

Misalkan kita ingin menciptakan model regresi \[\hat{Y_i} = \hat{\beta_0}+ \hat{\beta_1}x\]

model <- glm(mpg ~ disp,mtcars,family = "gaussian")
model <- lm(mpg ~ disp,mtcars) 
summary(model)

Call:
lm(formula = mpg ~ disp, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8922 -2.2022 -0.9631  1.6272  7.2305 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
disp        -0.041215   0.004712  -8.747 9.38e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.251 on 30 degrees of freedom
Multiple R-squared:  0.7183,    Adjusted R-squared:  0.709 
F-statistic: 76.51 on 1 and 30 DF,  p-value: 9.38e-10

Coba bandingkan model yang dibuat oleh fungsi lm dan glm !

summary(model)

Call:
lm(formula = mpg ~ disp, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8922 -2.2022 -0.9631  1.6272  7.2305 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
disp        -0.041215   0.004712  -8.747 9.38e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.251 on 30 degrees of freedom
Multiple R-squared:  0.7183,    Adjusted R-squared:  0.709 
F-statistic: 76.51 on 1 and 30 DF,  p-value: 9.38e-10

{plot(mtcars$disp,mtcars$mpg)
abline(model)}

7.3 Model Regresi Polynomial

Misalkan kita ingin menciptakan model regresi \[\hat{Y_i} = \hat{\beta_0}+ \hat{\beta_1}x + \hat{\beta_2}x^2 + \hat{\beta_3}x^3 \]

model_poly <- lm(mpg ~ poly(disp,3),mtcars) 
summary(model_poly)

Call:
lm(formula = mpg ~ poly(disp, 3), data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0896 -1.5653 -0.3619  1.4368  4.7617 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     20.0906     0.3931  51.112  < 2e-16 ***
poly(disp, 3)1 -28.4410     2.2235 -12.791 3.26e-13 ***
poly(disp, 3)2   9.1524     2.2235   4.116 0.000308 ***
poly(disp, 3)3  -9.7446     2.2235  -4.382 0.000150 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.224 on 28 degrees of freedom
Multiple R-squared:  0.8771,    Adjusted R-squared:  0.8639 
F-statistic: 66.58 on 3 and 28 DF,  p-value: 7.347e-13

x = seq(0,500,0.1)
y= predict(model_poly,list(disp = x))
{plot(mtcars$disp,mtcars$mpg,col="red",pch = 16)
abline(model,col="blue")
lines(x,y,col="green")}

Perhatikan bahwa regresi polinom dengan pangkat yang lebih tinggi menghasilkan kurva regresi yang lebih cocok ke titik-titik data. Akan tetapi, polinom dengan pangkat yang terlalu tinggi dapat menyebabkan overfit.

7.4 Multiple Linear Regression/Model Regresi Berganda

Misalkan kita ingin menciptakan model regresi \[\hat{Y_i} = \hat{\beta_0}+ \hat{\beta_1}x_1 + \hat{\beta_2}x_2 +...+ \hat{\beta_p}x_p\]

model_multi <- lm(mpg ~.,mtcars) 
summary(model_multi)

Call:
lm(formula = mpg ~ ., data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4506 -1.6044 -0.1196  1.2193  4.6271 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 12.30337   18.71788   0.657   0.5181  
cyl         -0.11144    1.04502  -0.107   0.9161  
disp         0.01334    0.01786   0.747   0.4635  
hp          -0.02148    0.02177  -0.987   0.3350  
drat         0.78711    1.63537   0.481   0.6353  
wt          -3.71530    1.89441  -1.961   0.0633 .
qsec         0.82104    0.73084   1.123   0.2739  
vs           0.31776    2.10451   0.151   0.8814  
am           2.52023    2.05665   1.225   0.2340  
gear         0.65541    1.49326   0.439   0.6652  
carb        -0.19942    0.82875  -0.241   0.8122  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.65 on 21 degrees of freedom
Multiple R-squared:  0.869, Adjusted R-squared:  0.8066 
F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07

7.5 Stepwise Regresion

Dalam memilih variabel yang ingin ditambahkan/dibuang dalam stepwise regresion, sebenarnya ada banyak criterion/kriteria yang dapat dipakai.
Buku walpole menggunakan criterion \(R^2\) dalam menyeleksi variabel.

model_step <- stepAIC(lm(mpg~.,mtcars),direction = "both")

Start:  AIC=70.9
mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb

       Df Sum of Sq    RSS    AIC
- cyl   1    0.0799 147.57 68.915
- vs    1    0.1601 147.66 68.932
- carb  1    0.4067 147.90 68.986
- gear  1    1.3531 148.85 69.190
- drat  1    1.6270 149.12 69.249
- disp  1    3.9167 151.41 69.736
- hp    1    6.8399 154.33 70.348
- qsec  1    8.8641 156.36 70.765
<none>              147.49 70.898
- am    1   10.5467 158.04 71.108
- wt    1   27.0144 174.51 74.280

Step:  AIC=68.92
mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb

       Df Sum of Sq    RSS    AIC
- vs    1    0.2685 147.84 66.973
- carb  1    0.5201 148.09 67.028
- gear  1    1.8211 149.40 67.308
- drat  1    1.9826 149.56 67.342
- disp  1    3.9009 151.47 67.750
- hp    1    7.3632 154.94 68.473
<none>              147.57 68.915
- qsec  1   10.0933 157.67 69.032
- am    1   11.8359 159.41 69.384
+ cyl   1    0.0799 147.49 70.898
- wt    1   27.0280 174.60 72.297

Step:  AIC=66.97
mpg ~ disp + hp + drat + wt + qsec + am + gear + carb

       Df Sum of Sq    RSS    AIC
- carb  1    0.6855 148.53 65.121
- gear  1    2.1437 149.99 65.434
- drat  1    2.2139 150.06 65.449
- disp  1    3.6467 151.49 65.753
- hp    1    7.1060 154.95 66.475
<none>              147.84 66.973
- am    1   11.5694 159.41 67.384
- qsec  1   15.6830 163.53 68.200
+ vs    1    0.2685 147.57 68.915
+ cyl   1    0.1883 147.66 68.932
- wt    1   27.3799 175.22 70.410

Step:  AIC=65.12
mpg ~ disp + hp + drat + wt + qsec + am + gear

       Df Sum of Sq    RSS    AIC
- gear  1     1.565 150.09 63.457
- drat  1     1.932 150.46 63.535
<none>              148.53 65.121
- disp  1    10.110 158.64 65.229
- am    1    12.323 160.85 65.672
- hp    1    14.826 163.35 66.166
+ carb  1     0.685 147.84 66.973
+ vs    1     0.434 148.09 67.028
+ cyl   1     0.414 148.11 67.032
- qsec  1    26.408 174.94 68.358
- wt    1    69.127 217.66 75.350

Step:  AIC=63.46
mpg ~ disp + hp + drat + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- drat  1     3.345 153.44 62.162
- disp  1     8.545 158.64 63.229
<none>              150.09 63.457
- hp    1    13.285 163.38 64.171
+ gear  1     1.565 148.53 65.121
+ cyl   1     1.003 149.09 65.242
+ vs    1     0.645 149.45 65.319
+ carb  1     0.107 149.99 65.434
- am    1    20.036 170.13 65.466
- qsec  1    25.574 175.67 66.491
- wt    1    67.572 217.66 73.351

Step:  AIC=62.16
mpg ~ disp + hp + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- disp  1     6.629 160.07 61.515
<none>              153.44 62.162
- hp    1    12.572 166.01 62.682
+ drat  1     3.345 150.09 63.457
+ gear  1     2.977 150.46 63.535
+ cyl   1     2.447 150.99 63.648
+ vs    1     1.121 152.32 63.927
+ carb  1     0.011 153.43 64.160
- qsec  1    26.470 179.91 65.255
- am    1    32.198 185.63 66.258
- wt    1    69.043 222.48 72.051

Step:  AIC=61.52
mpg ~ hp + wt + qsec + am

       Df Sum of Sq    RSS    AIC
- hp    1     9.219 169.29 61.307
<none>              160.07 61.515
+ disp  1     6.629 153.44 62.162
+ carb  1     3.227 156.84 62.864
+ drat  1     1.428 158.64 63.229
- qsec  1    20.225 180.29 63.323
+ cyl   1     0.249 159.82 63.465
+ vs    1     0.249 159.82 63.466
+ gear  1     0.171 159.90 63.481
- am    1    25.993 186.06 64.331
- wt    1    78.494 238.56 72.284

Step:  AIC=61.31
mpg ~ wt + qsec + am

       Df Sum of Sq    RSS    AIC
<none>              169.29 61.307
+ hp    1     9.219 160.07 61.515
+ carb  1     8.036 161.25 61.751
+ disp  1     3.276 166.01 62.682
+ cyl   1     1.501 167.78 63.022
+ drat  1     1.400 167.89 63.042
+ gear  1     0.123 169.16 63.284
+ vs    1     0.000 169.29 63.307
- am    1    26.178 195.46 63.908
- qsec  1   109.034 278.32 75.217
- wt    1   183.347 352.63 82.790

*Untuk forward selection, buka dokumentasi ?MASS:stepAIC, object diisi dengan model regresi kosonglm(y~1), scope diisi dengan list(upper= model_Regresi_paling_komplit,lower = model_kosong)

7.6 Interaksi

Model interaksi dibangun ketika terdapat multikolienaritas(prediktor tidak saling bebas).
Model yang kita dapatkan dari stepwise regresion adalah model menggunakan 2 prediktor, yakni wt dan cyl.
Dari section 6.1, kita dapat melihat bahwa wt dan cyl memiliki koefisien korelasi sebesar 0.856.
Oleh karena itu, kita akan membangun model \[\hat{Y_i} = \hat{\beta_0}+ \hat{\beta_1}x_1 + \hat{\beta_2}x_2 + \hat{\beta_3}x_1 x_2 \] atau

\[\hat{Y_i} = \hat{\beta_0}+ \hat{\beta_1}x_1 + (\hat{\beta_2} + \hat{\beta_3}x_1) x_2 \] Membuat model dengan interaksi

model_interaksi <- lm(mpg ~ wt+ cyl + wt*cyl,mtcars)
summary(model_interaksi)

Call:
lm(formula = mpg ~ wt + cyl + wt * cyl, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2288 -1.3495 -0.5042  1.4647  5.2344 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  54.3068     6.1275   8.863 1.29e-09 ***
wt           -8.6556     2.3201  -3.731 0.000861 ***
cyl          -3.8032     1.0050  -3.784 0.000747 ***
wt:cyl        0.8084     0.3273   2.470 0.019882 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.368 on 28 degrees of freedom
Multiple R-squared:  0.8606,    Adjusted R-squared:  0.8457 
F-statistic: 57.62 on 3 and 28 DF,  p-value: 4.231e-12

7.7 Menguji asumsi pada model regresi linear

#Asumsi pertama, galat berdistribusi normal
shapiro.test(model_step$residuals)

    Shapiro-Wilk normality test

data:  model_step$residuals
W = 0.93745, p-value = 0.06341

#Asumsi kedua, variansi konstan
plot(model_step)

8 Tipe data kategorikal

Mari lihat kembali dataset mtcars

#attach(mtcars) harap di run sekali saja
str(mtcars)

'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

Melalui keterangan data, vs dan am adalah variabel kategorikal, akan tetapi R membaca kedua variabel tersebut sebagai variabel numerik.
Oleh karena itu kita akan mengubah kedua variabel tersebut menjadi variabel kategorikal/faktor di R

data <- mtcars
data$vs <- factor(data$vs)
data$am <- factor(data$am)
str(data)

'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : Factor w/ 2 levels "0","1": 1 1 2 2 1 2 1 2 2 2 ...
 $ am  : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

Sekarang kita coba buat model regresi multiple, dengan variabel kategori yang sudah kita koreksi

model <- lm(mpg~.,data)
summary(model)

Call:
lm(formula = mpg ~ ., data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4506 -1.6044 -0.1196  1.2193  4.6271 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 12.30337   18.71788   0.657   0.5181  
cyl         -0.11144    1.04502  -0.107   0.9161  
disp         0.01334    0.01786   0.747   0.4635  
hp          -0.02148    0.02177  -0.987   0.3350  
drat         0.78711    1.63537   0.481   0.6353  
wt          -3.71530    1.89441  -1.961   0.0633 .
qsec         0.82104    0.73084   1.123   0.2739  
vs1          0.31776    2.10451   0.151   0.8814  
am1          2.52023    2.05665   1.225   0.2340  
gear         0.65541    1.49326   0.439   0.6652  
carb        -0.19942    0.82875  -0.241   0.8122  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.65 on 21 degrees of freedom
Multiple R-squared:  0.869, Adjusted R-squared:  0.8066 
F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07

vs1 berarti jika nilai vs=1, maka nilai mpg bertambah sebesar 0.31776. jika vs=0, berarti nilai mpg tidak bertambah.

9 dplyr dan tidyr

dplyr dan tidyr adalah packages yang digunakan untuk melakukan data wrangling (merapikan data).
Untuk informasi lebih lanjut kunjungi https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf

library(tidyr)
library(dplyr)

Attaching package: <U+393C><U+3E31>dplyr<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:

    filter, lag

The following objects are masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    intersect, setdiff, setequal, union

dt <- read.csv("tidy.csv")
head(dt)

Pada package dplyr, dikenalkan pipe operation %>%.
f(g(a),b) \(\equiv\) a %>% g() %>% f(b)

dt %>% select(Tanggal,mobil)

dt %>% filter(mobil>70)

dt %>% gather("kendaraan","jumlah",2:4)

dt %>% separate(Tanggal,c("tanggal","bulan","tahun"),sep ="/")

10 Data cleansing sederhana

Dua jenis data kotor yang umum dijumpai adalah data kotor akibat salah input dan akibat nilai kosong(Missing Value).
Untuk memperbaiki kesalahan akibat salah input, tidak ada metode umum perbaikan(tergantung kasus).
Untuk masalah missing values, metode paling umum adalah melakukan subtitusi nilai yang hilang dengan median dari variabel tersebut

mis <- dt$bus #Mengambil kolom yang ingin kita perbaiki
mis[is.na(mis)] <- median(mis,na.rm = T) #Melakukan subtitusi elemen missing value dengan median dari variabel tersebut. na.rm = T berarti observasi yang NA (Not Avaiable)/missing kita abaikan dalam perhitungan median.
dt$bus <- mis
rm(mis)

Untuk pembahasan lebih detail mengenai data cleansing, ada di modul “Basic Data Cleansing”

11 Regresi Logistik

Regresi logistik adalah regresi yang response nya merupakan nilai kategorikal (antara 0 atau 1), merupakan salah satu GLM, dimana response mengikuti distribusi binomial.

\[Pr(y|x) = \frac{1}{1+e^{-(\beta x)}} \]

11.1 Analisis dataset Credit Scoring

credit <- read.csv("npl_train.csv",dec=".",sep = ",")
str(credit)

'data.frame':   92 obs. of  11 variables:
 $ jumlah_kartu                 : int  2 2 3 4 2 2 2 4 2 4 ...
 $ outstanding                  : int  36158 268691 6769149 3496732 9402085 6227439 3906290 9534837 4145065 1818606 ...
 $ limit_kredit                 : num  7.0e+06 1.0e+07 2.8e+07 2.1e+07 1.0e+07 8.0e+07 4.0e+06 2.0e+07 5.0e+06 7.0e+06 ...
 $ tagihan                      : int  23437 254564 4159779 111231 6099283 2081248 2043682 3692028 4021399 1765911 ...
 $ total_pemakaian_retail       : num  94 1012 0 2536660 2666558 ...
 $ sisa_tagihan_tidak_terbayar  : int  26323 0 0 581334 5951865 4613435 3314046 7881069 4122425 1627786 ...
 $ total_pemakaian              : num  94 1012 0 2536660 2666558 ...
 $ sisa_tagihan_per_jumlah_kartu: num  13162 0 0 145334 2975933 ...
 $ sisa_tagihan_per_limit       : num  0.00376 0 0 0.02768 0.59519 ...
 $ total_pemakaian_per_limit    : num  1.34e-05 1.01e-04 0.00 1.21e-01 2.67e-01 ...
 $ flag_kredit_macet            : int  0 0 0 0 0 0 0 0 0 0 ...

11.2 Penentuan base level

Kita dapat memiih levels mana dari suatu factor yang dijadikan sebagai base level

a <- factor(c("a","b","c"))
a <- relevel(a,3)#Kita menjadikan elemen ketiga menjadi base level

11.3 Pembuatan Regresi logistik

model_rl <- glm(flag_kredit_macet~.,credit,family = "binomial")
summary(model_rl)

Call:
glm(formula = flag_kredit_macet ~ ., family = "binomial", data = credit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6908  -0.4744  -0.2848  -0.1391   2.1794  

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)
(Intercept)                   -6.231e-01  1.815e+00  -0.343    0.731
jumlah_kartu                  -8.945e-01  7.706e-01  -1.161    0.246
outstanding                   -3.186e-07  2.533e-07  -1.258    0.208
limit_kredit                   1.825e-08  2.326e-08   0.785    0.433
tagihan                        2.732e-07  1.973e-07   1.385    0.166
total_pemakaian_retail        -1.616e-06  1.767e-06  -0.915    0.360
sisa_tagihan_tidak_terbayar    1.619e-07  2.131e-07   0.759    0.448
total_pemakaian                2.190e-06  1.994e-06   1.098    0.272
sisa_tagihan_per_jumlah_kartu -1.416e-07  2.040e-07  -0.694    0.488
sisa_tagihan_per_limit         1.096e+00  1.162e+00   0.943    0.346
total_pemakaian_per_limit     -1.106e+01  8.432e+00  -1.311    0.190

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 70.966  on 90  degrees of freedom
Residual deviance: 53.021  on 80  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 75.021

Number of Fisher Scoring iterations: 7

11.4 Confusion Matrix

Untuk melihat confusion matrix dari hasil prediksi model, gunakan fungsi confussionMatrix dari package

prediksi <- predict(model_rl,credit,type = "response") #Harus menggunakan type = "response" untuk prediksi probabilitas
prediksi_biner <- ifelse(prediksi > 0.5,1,0) #Menggunakan threshold 0.5 untuk menentukan true/false
library(caret)

Loading required package: lattice
Loading required package: ggplot2

Attaching package: <U+393C><U+3E31>ggplot2<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>mtcars<U+393C><U+3E32>:

    mpg

confusionMatrix(factor(prediksi_biner),factor(credit$flag_kredit_macet))

Confusion Matrix and Statistics

          Reference
Prediction  0  1
         0 78  9
         1  1  3
                                         
               Accuracy : 0.8901         
                 95% CI : (0.8072, 0.946)
    No Information Rate : 0.8681         
    P-Value [Acc > NIR] : 0.33223        
                                         
                  Kappa : 0.3309         
 Mcnemar's Test P-Value : 0.02686        
                                         
            Sensitivity : 0.9873         
            Specificity : 0.2500         
         Pos Pred Value : 0.8966         
         Neg Pred Value : 0.7500         
             Prevalence : 0.8681         
         Detection Rate : 0.8571         
   Detection Prevalence : 0.9560         
      Balanced Accuracy : 0.6187         
                                         
       'Positive' Class : 0              
                                         

11.5 ROC Curve

library(pROC)
roc_var <- roc(response = credit$flag_kredit_macet,predictor = prediksi,smooth = F,auc=T) 
#Untuk plot ROC yang halus, atur smooth menjadi T
{par(pty = "s") #Supaya hasil plot berbentuk persegi
plot.roc(roc_var,auc.polygon=T,print.auc = T,legacy.axes = T,add=F)}

12 Generalized Linear Model

Dalam pembuatan model GLM, ada 2 hal yang harus kita tentukan. Yakni distribusi dan linkfunction. Berikut distribusi beserta link function bawaannya.

-binomial(link = “logit”)
-gaussian(link = “identity”)
-Gamma(link = “inverse”)
-inverse.gaussian(link = “1/mu^2”)
-poisson(link = “log”)
-quasi(link = “identity”, variance = “constant”)
-quasibinomial(link = “logit” )
-quasipoisson(link = “log”)

model1 <- glm(y~x,data,family = "Gamma") #Pembuatan model GLM distribusi Gamma dengan linkfunction inverse
model2 <- glm(y~x,data,family= binomial(link = "probit")) #Pembuatan model GLM distribusi binomial/bernoulli dengan link function probit

*Untuk pendefinisian offset dan weights, silahkan baca dokumentasi ?glm