PENDAHULUAN

Tujuan analisis ini adalah membuat model untuk memprediksi kekuatan kompresi (strength) berdasarkan komposisi campurannya. Data yang digunakan adalah data Concrete Analysis.

Yang ingin dicapai pada pembuatan model:

  1. Metrics :

SET UP

library(tidyverse)
library(GGally)
library(caret)
library(car)
library(MLmetrics)

DATA PREPROCESSING

concrete_train <- read.csv("data-train.csv")
concrete_test <- read.csv("data-test.csv")

glimpse(concrete_train)
## Rows: 825
## Columns: 10
## $ id          <chr> "S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10…
## $ cement      <dbl> 540.0, 540.0, 332.5, 332.5, 198.6, 380.0, 380.0, 475.0, 19…
## $ slag        <dbl> 0.0, 0.0, 142.5, 142.5, 132.4, 95.0, 95.0, 0.0, 132.4, 132…
## $ flyash      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ water       <dbl> 162, 162, 228, 228, 192, 228, 228, 228, 192, 192, 228, 228…
## $ super_plast <dbl> 2.5, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0…
## $ coarse_agg  <dbl> 1040.0, 1055.0, 932.0, 932.0, 978.4, 932.0, 932.0, 932.0, …
## $ fine_agg    <dbl> 676.0, 676.0, 594.0, 594.0, 825.5, 594.0, 594.0, 594.0, 82…
## $ age         <int> 28, 28, 270, 365, 360, 365, 28, 28, 90, 28, 28, 90, 90, 36…
## $ strength    <dbl> 79.99, 61.89, 40.27, 41.05, 44.30, 43.70, 36.45, 39.29, 38…
length(unique(concrete_train$id))
## [1] 825

Keterangan variabel:

id: Id of each cement mixture, cement: The amount of cement (Kg) in a m3 mixture, slag: The amount of blast furnace slag (Kg) in a m3 mixture, flyash: The amount of fly ash (Kg) in a m3 mixture, water: The amount of water (Kg) in a m3 mixture, super_plast: The amount of Superplasticizer (Kg) in a m3 mixture, coarse_agg: The amount of Coarse Aggreagate (Kg) in a m3 mixture, fine_agg: The amount of Fine Aggreagate (Kg) in a m3 mixture, age: the number of resting days before the compressive strength measurement, strength: Concrete compressive strength measurement in MPa unit.

Dari struktur data di atas, semua tipe data yang diperlukan sudah bertipe numerik (tidak perlu dilakukan coercion).

Cek summary data.

summary(concrete_train)
##       id                cement           slag            flyash      
##  Length:825         Min.   :102.0   Min.   :  0.00   Min.   :  0.00  
##  Class :character   1st Qu.:194.7   1st Qu.:  0.00   1st Qu.:  0.00  
##  Mode  :character   Median :275.1   Median : 20.00   Median :  0.00  
##                     Mean   :280.9   Mean   : 73.18   Mean   : 54.03  
##                     3rd Qu.:350.0   3rd Qu.:141.30   3rd Qu.:118.20  
##                     Max.   :540.0   Max.   :359.40   Max.   :200.10  
##      water        super_plast       coarse_agg        fine_agg    
##  Min.   :121.8   Min.   : 0.000   Min.   : 801.0   Min.   :594.0  
##  1st Qu.:164.9   1st Qu.: 0.000   1st Qu.: 932.0   1st Qu.:734.0  
##  Median :184.0   Median : 6.500   Median : 968.0   Median :780.1  
##  Mean   :181.1   Mean   : 6.266   Mean   : 972.8   Mean   :775.6  
##  3rd Qu.:192.0   3rd Qu.:10.100   3rd Qu.:1028.4   3rd Qu.:826.8  
##  Max.   :247.0   Max.   :32.200   Max.   :1145.0   Max.   :992.6  
##       age            strength    
##  Min.   :  1.00   Min.   : 2.33  
##  1st Qu.:  7.00   1st Qu.:23.64  
##  Median : 28.00   Median :34.57  
##  Mean   : 45.14   Mean   :35.79  
##  3rd Qu.: 56.00   3rd Qu.:45.94  
##  Max.   :365.00   Max.   :82.60

Untuk membuat model dengan lm belum perlu dilakukan scaling.

Cek outliers

Cek outliers pada setiap variabel.

par(mfrow=c(2,3))

p1 <- boxplot(concrete_train$cement, xlab = "cement", cex.lab = 1.2)
p2 <- boxplot(concrete_train$slag, xlab = "slag (blast furnace)", cex.lab = 1.2 )
p3 <- boxplot(concrete_train$flyash,  xlab = "fly ash", cex.lab = 1.2)
p4 <- boxplot(concrete_train$water,  xlab = "water", cex.lab = 1.2)
p5 <- boxplot(concrete_train$super_plast,  xlab = "super plasticizer ", cex.lab = 1.2)
p6 <- boxplot(concrete_train$coarse_agg,  xlab = "coarse aggregate", cex.lab = 1.2)

p7 <- boxplot(concrete_train$fine_agg,xlab = "fine aggregate", cex.lab = 1.2 )
p8 <- boxplot(concrete_train$age,xlab = "age", cex.lab = 1.2)
p9 <- boxplot(concrete_train$strength,xlab = "strength", cex.lab = 1.2)

Dari visualisasi boxplot di atas, terlihat ada outliers pada variabel: slag, water, super_plast, fine_agg, dan age.

Influence Outliers?

Karena ada outliers, kita akan cek setiap outliers apakah high atau low influence outliers. Bila merupakan low influence maka kita tidak perlu membuang data tersebut karena walaupun nilainya ekstrim tetapi tidak banyak berpengaruh terhadap nilai koefisien.

  • Sebuah titik data dikatakan high influence ketika data tersebut sangat mempengaruhi bentuk garis linear (nilai intercept, slope) maupun nilai r-squared dari model.

Kita coba membuat plot hubungan linear antara prediktor dan target. Garis warna merah adalah garis linear dari model yang menyertakan outlier pada prediktirnya. Garis biru adalah garis linear dari model yang tidak menyertakan outliernya.

==slag==

ols_outlier_slag <- lm(strength ~ slag, concrete_train) #model dengan outlier slag
ols_outlier_slag
## 
## Call:
## lm(formula = strength ~ slag, data = concrete_train)
## 
## Coefficients:
## (Intercept)         slag  
##    33.96946      0.02485
# buang outliers
conc_train_clean_nooslag <- concrete_train %>% 
  filter (slag < 340)

ols_no_slag <- lm(strength ~ slag, conc_train_clean_nooslag) #model tanpa outliers slag

ols_no_slag
## 
## Call:
## lm(formula = strength ~ slag, data = conc_train_clean_nooslag)
## 
## Coefficients:
## (Intercept)         slag  
##    33.87073      0.02698
plot(concrete_train$slag, concrete_train$strength)
abline(ols_outlier_slag$coefficients[1], ols_outlier_slag$coefficients[2], col = "red")
abline(ols_no_slag$coefficients[1], ols_no_slag$coefficients[2], col = "blue")

==water==

ols_outlier_water <- lm(strength ~ water, concrete_train) #model dengan outliers water
ols_outlier_water
## 
## Call:
## lm(formula = strength ~ water, data = concrete_train)
## 
## Coefficients:
## (Intercept)        water  
##     76.3040      -0.2237
# buang outliers
conc_train_clean_noowater <- concrete_train %>% 
  filter (water<230)

ols_no_water <- lm(strength ~ water, conc_train_clean_noowater ) #model tanpa outliers water

ols_no_water
## 
## Call:
## lm(formula = strength ~ water, data = conc_train_clean_noowater)
## 
## Coefficients:
## (Intercept)        water  
##     77.4496      -0.2303
plot(concrete_train$water, concrete_train$strength)
abline(ols_outlier_water$coefficients[1], ols_outlier_water$coefficients[2], col = "red")
abline(ols_no_water$coefficients[1], ols_no_water$coefficients[2], col = "blue")

==super_plast==

ols_outlier_sp <- lm(strength ~ super_plast, concrete_train) #model dengan outliers super_plast
                    
ols_outlier_sp
## 
## Call:
## lm(formula = strength ~ super_plast, data = concrete_train)
## 
## Coefficients:
## (Intercept)  super_plast  
##      29.449        1.012
# buang outliers
conc_train_clean_nosp <- concrete_train %>% 
  filter (super_plast<25)

ols_no_sp <- lm(strength ~ super_plast, conc_train_clean_nosp) #model tanpa outliers super_plast

ols_no_sp
## 
## Call:
## lm(formula = strength ~ super_plast, data = conc_train_clean_nosp)
## 
## Coefficients:
## (Intercept)  super_plast  
##      29.303        1.043
plot(concrete_train$super_plast, concrete_train$strength)
abline(ols_outlier_sp$coefficients[1], ols_outlier_sp$coefficients[2], col = "red")
abline(ols_no_sp$coefficients[1], ols_no_sp$coefficients[2], col = "blue")

==fine_agg==

ols_outlier_fagg <- lm(strength ~ fine_agg, concrete_train) #model dengan outliers fine_agg
ols_outlier_fagg
## 
## Call:
## lm(formula = strength ~ fine_agg, data = concrete_train)
## 
## Coefficients:
## (Intercept)     fine_agg  
##    65.66147     -0.03852
# buang outliers
conc_train_clean_nofagg <- concrete_train %>% 
  filter (fine_agg< 960)

ols_no_fagg <- lm(strength ~ fine_agg, conc_train_clean_nofagg ) #model tanpa outliersfine_agg

ols_no_fagg
## 
## Call:
## lm(formula = strength ~ fine_agg, data = conc_train_clean_nofagg)
## 
## Coefficients:
## (Intercept)     fine_agg  
##    69.33996     -0.04345
plot(concrete_train$fine_agg, concrete_train$strength)
abline(ols_outlier_fagg$coefficients[1], ols_outlier_fagg$coefficients[2], col = "red")
abline(ols_no_fagg$coefficients[1], ols_no_fagg$coefficients[2], col = "blue")

==age==

ols_outlier_age <- lm(strength ~ age, concrete_train) #model dengan outliers age
ols_outlier_age
## 
## Call:
## lm(formula = strength ~ age, data = concrete_train)
## 
## Coefficients:
## (Intercept)          age  
##    31.63304      0.09204
# buang outliers
conc_train_clean_noage <- concrete_train %>% 
  filter (age< 110)

ols_no_age <- lm(strength ~ age, conc_train_clean_noage ) #model tanpa outliers age

ols_no_age
## 
## Call:
## lm(formula = strength ~ age, data = conc_train_clean_noage)
## 
## Coefficients:
## (Intercept)          age  
##     25.4506       0.3031
plot(concrete_train$age, concrete_train$strength)
abline(ols_outlier_age$coefficients[1], ols_outlier_age$coefficients[2], col = "red")
abline(ols_no_age$coefficients[1], ols_no_age$coefficients[2], col = "blue")

Dari plot di atas, tampak yang merupakan high influence outliers adalah outliers dari variableage. Kita coba lihat pengaruhnya pada R-squared:

summary(ols_outlier_age)$r.squared #model dengan outliers age
## [1] 0.1167578
summary(ols_no_age)$r.squared #model tanpa outliers age
## [1] 0.2624483

Nilai R-squared lebih tinggi pada model tanpa outliers age.

Maka kita filter data yang tidak menyertakan outliers age:

concrete_train_f <- concrete_train %>% 
  filter(age < 110)

Variabel id tidak kita gunakan dalam analisis karena merupakan data yang sangat unik dan tidak mempunyai peran penting dalam analisa. Kita gunakan fungsi select untuk membuang kolom id pada data concrete_train.

concrete_train_clean <- concrete_train_f %>% 
  select(-id)

MODEL FITTING & EVALUASI

Cross Validation

Karena data concrete_test variabel targetnya masih bernilai NA, maka kita lakukan cross validasi pada data train dengan proporsi sebesar 80% untuk train model dan 20% untuk validasi model.

library(rsample)

RNGkind(sample.kind = "Rounding")
set.seed(100)

index <- sample(nrow(concrete_train_clean),
                nrow(concrete_train_clean)*0.7)

data_train <- concrete_train_clean[index,]  # data train: digunakan untuk melatih model
data_val <- concrete_train_clean[-index,] # data validation : digunakan untuk validasi model

Cek missing values pada data train.

is.na(data_train) %>% colSums()
##      cement        slag      flyash       water super_plast  coarse_agg 
##           0           0           0           0           0           0 
##    fine_agg         age    strength 
##           0           0           0

Tidak terdapat missing values pada data concrete_train, maka kita bisa lanjut ke proses berikutnya.

MODEL

Model Neural Network

Kita akan coba menggunakan model Neural Network dengan keras mengikuti panduan langkah-langkah sebagai berikut:

  • Jumlah variabel target : 1

  • Jumlah variabel prediktor : 8

  • Berapakah dimensi data_train?

dim(data_train)
## [1] 543   9

Further Data Pre-Processing

Matrix -> Array + Scalling

Sebelum membuat model dengan Keras, ada beberapa hal yang perlu dilakukan untuk mempersiapkan data:

  1. Memisahkan prediktor dengan target variabel
  2. Mengubah prediktor menjadi matrix/array
  3. Untuk prediktor range data diubah menjadi 0-1, supaya mengurangi beban komputasi dan supaya model lebih cepat konvergen
  4. Melakukan one-hot encoding apabila target variabel adalah kategori. Namun karena pada kasus ini variabelnya adalah numerik semua makan tidak perlu dilakukan one-hot encoding
#scaling

train_x <- data_train %>% 
  select(-strength) %>% # seleksi untuk variabel prediktor
  as.matrix() %>% # diubah menjadi matrix 
  scale()

train_y <- data_train$strength # seleksi untuk variabel target

val_x <- data_val %>% 
  select(-strength) %>% 
  as.matrix() %>% 
  scale()

val_y <- data_val$strength

Untuk mengubah feature agar dapat diterima oleh keras dan python, kita harus ubah ke format array. menggunakan fungsi array_reshape(x,dim)

library(keras)

train_x_array <- array_reshape(train_x, dim=dim(train_x))
                                          
val_x_array <- array_reshape(val_x, dim=dim(val_x))

Modeling

Architecture

Mengatur set_seed(100)

Langkah selanjutnya adalah membangun arsitektur Neural Network. Beberapa ketentuan ketika membuat arsitektur Neural Network:

  1. Selalu diawali keras_model_sequential()
  2. Layer pertama yang dibuat akan menjadi hidden layer pertama
  3. Input layer dibuat dengan memasukkan parameter input_shape pada layer pertama
  4. Layer terakhir yang dibuat akan menjadi output layer

Pertama, kita buat dulu object untuk menyimpan informasi jumlah kolom dari prediktor dan jumlah kategori dari target variabel.

RNGkind(sample.kind = "Rounding")
set.seed(100)
initializer <- initializer_random_normal(seed = 100)


input_dim <- ncol(train_x)

Model base

Model base kita buat dengan parameter sebagai berikut: - layer pertama berisi 128 nodes, fungsi aktivasi relu, 8 input shape - layer kedua berisi 64 nodes, fungsi aktivasi relu - layer ketiga berisi 1 nodes, fungsi aktivasi linear

Untuk melatih model, kita susun model dengan menentukan loss function, jenis optimizer, dan metrik evaluasi. Lalu compile model dengan menggunakan parameter berikut: - MAPE sebagai loss function - optimizer_adam sebagai optimizer dengan learning rate 0.001 - gunakan mean-absolute-error sebagai metrik evaluasi karena kasus ini adalah kasus regresi di mana kita memprediksi suatu nilai numerik (strength)

Kemudian kita lakukan model fitting dengan menggunakan - epoch = 20 - batch_size = 150 dan menambahkan - parameter shuffle = F agar sampel pada tiap batch tidak diambil secara random melainkan terurut (sequence).

model_base <- keras_model_sequential(name = "model_conc1") %>% 
  
  # input layer + hidden layer 1
  layer_dense(input_shape = input_dim, # dimensi dari prediktor / jumlah pixel
              units = 128, # jumlah nodes pada hidden layer 1
              activation = "relu", # activation function
              name = "hidden_1") %>% 
  
  # hidden layer 2
  layer_dense(units = 64, # jumlah nodes pada hidden 2
              activation = "relu",
              name = "hidden_2") %>% 
  
  
  # output layer
  layer_dense(units = 1, # jumlah nodes pada output / jumlah kelas target
              activation = "linear",
              name = "output")

model_base
## Model: "model_conc1"
## ________________________________________________________________________________
## Layer (type)                        Output Shape                    Param #     
## ================================================================================
## hidden_1 (Dense)                    (None, 128)                     1152        
## ________________________________________________________________________________
## hidden_2 (Dense)                    (None, 64)                      8256        
## ________________________________________________________________________________
## output (Dense)                      (None, 1)                       65          
## ================================================================================
## Total params: 9,473
## Trainable params: 9,473
## Non-trainable params: 0
## ________________________________________________________________________________

Compile Model

model_base %>% 
  compile(loss = "MAPE",
          optimizer = optimizer_adam(lr=0.01),
          metrics = "mean_absolute_error")

Train Model

history_base <- model_base %>% fit(x = train_x_array, # data prediktor
                          y = train_y, # data target
                          epoch = 20, # jumlah step
                          validation_data = list(val_x_array, val_y), # data test (masukkan prediktor dan targetnya)
                          batch_size = 100, # jumlah data per batch
                          verbose = 1,# untuk tracing tiap epoch
                          shuffle = F) 
plot(history_base)

Model Prediction & Evaluation Data Validation

Predict Data Validation

pred_base <- predict(model_base, val_x_array)

head(as.data.frame(pred_base))
##         V1
## 1 39.01366
## 2 43.88069
## 3 41.79082
## 4 27.58593
## 5 62.16134
## 6 67.11030
data_val <- data_val %>% cbind(as.data.frame(pred_base))

Evaluation Data Validation

library(caret)
R2_keras <- R2(data_val$strength, data_val$V1)
R2_keras
## [1] 0.6962679

Masukkan nilai strength untuk submission data_test.

data_test <- concrete_test

test_x <- data_test %>% 
  select(-strength, -id) %>% 
  as.matrix() %>% 
  scale()


test_x_array <- array_reshape(test_x, dim=dim(test_x))

# predict target using your model
pred_test3 <- predict(model_base, test_x_array)

as.data.frame(pred_test3)
##           V1
## 1   32.57458
## 2   24.22050
## 3   75.28716
## 4   29.94056
## 5   64.83277
## 6   38.93812
## 7   39.85194
## 8   43.47614
## 9   53.38425
## 10  90.64542
## 11  82.90844
## 12  58.98279
## 13  42.75138
## 14  92.69809
## 15  21.58995
## 16  51.54075
## 17  46.67363
## 18  68.85043
## 19  47.05177
## 20  43.65837
## 21  47.03204
## 22  47.03204
## 23  44.24274
## 24  64.84210
## 25  50.24693
## 26  57.54071
## 27  76.36779
## 28  57.39293
## 29  60.27695
## 30  55.03392
## 31  49.75308
## 32  54.04637
## 33  60.81905
## 34  60.81905
## 35  66.67397
## 36  17.62705
## 37  22.65858
## 38  32.42190
## 39  18.66640
## 40  28.91991
## 41  16.86837
## 42  19.93578
## 43  20.89763
## 44  24.62553
## 45  21.23414
## 46  22.74746
## 47  18.42588
## 48  19.72963
## 49  17.74524
## 50  21.53815
## 51  27.33511
## 52  26.08749
## 53  33.81071
## 54  23.96517
## 55  25.51530
## 56  41.87719
## 57  26.44429
## 58  36.55034
## 59  32.16727
## 60  27.76920
## 61  35.71074
## 62  42.28911
## 63  30.11659
## 64  36.28917
## 65  22.91234
## 66  26.35909
## 67  48.92063
## 68  42.51751
## 69  40.45080
## 70  43.89380
## 71  45.05704
## 72  19.53068
## 73  19.75510
## 74  21.77291
## 75  20.15277
## 76  21.25337
## 77  21.09168
## 78  24.00453
## 79  26.07332
## 80  23.06908
## 81  28.09646
## 82  26.00266
## 83  23.03241
## 84  32.16941
## 85  32.04466
## 86  29.05465
## 87  45.39694
## 88  46.07386
## 89  54.28979
## 90  53.51581
## 91  43.41455
## 92  43.19101
## 93  41.87725
## 94  52.56909
## 95  46.56916
## 96  54.42093
## 97  42.98298
## 98  15.96926
## 99  13.66328
## 100 19.35709
## 101 31.36906
## 102 44.97278
## 103 35.62719
## 104 14.81699
## 105 19.36644
## 106 18.61287
## 107 26.01368
## 108 31.43764
## 109 21.72281
## 110 27.93874
## 111 13.69605
## 112 25.45529
## 113 33.96722
## 114 74.60609
## 115 13.64897
## 116 61.04871
## 117 17.02780
## 118 20.84869
## 119 20.14068
## 120 25.21309
## 121 30.66693
## 122 27.51463
## 123 15.53213
## 124 22.85517
## 125 33.98450
## 126 16.83189
## 127 13.60995
## 128 22.38339
## 129 11.44278
## 130 30.89508
## 131 18.80628
## 132 13.13338
## 133 40.98435
## 134 15.63742
## 135 19.38500
## 136 32.84353
## 137 28.92084
## 138 42.32024
## 139 13.28489
## 140 10.26123
## 141 24.16053
## 142 10.80512
## 143 17.55909
## 144 29.48747
## 145 23.13919
## 146 36.35224
## 147 18.54743
## 148 40.47062
## 149 49.28079
## 150 21.78569
## 151 26.31833
## 152 17.81971
## 153 30.69268
## 154 25.47974
## 155 81.55798
## 156 15.29785
## 157 15.20171
## 158 18.33035
## 159 28.40182
## 160 44.33707
## 161 57.52894
## 162 42.02027
## 163 46.58547
## 164 29.91383
## 165 17.17127
## 166 28.47148
## 167 46.49734
## 168 37.17350
## 169 30.77803
## 170 30.75510
## 171 29.08006
## 172 37.13516
## 173 26.09907
## 174 41.54782
## 175 43.70145
## 176 30.83555
## 177 14.70391
## 178 26.30681
## 179 27.11364
## 180 25.98625
## 181 26.76139
## 182 23.11620
## 183 26.05705
## 184 23.99394
## 185 28.89725
## 186 17.40555
## 187 20.80378
## 188 18.34368
## 189 32.14690
## 190 18.46708
## 191 37.92815
## 192 26.57012
## 193 12.71867
## 194 32.04294
## 195 23.89338
## 196 25.08780
## 197 38.79372
## 198 36.74041
## 199 33.06232
## 200 37.11152
## 201 30.51226
## 202 31.01650
## 203 28.36112
## 204 30.85957
## 205 27.14825
data_test <- data_test %>% cbind(as.data.frame(pred_test3))


# Create submission data
submission <- data.frame(id = data_test$id,
                         strength = pred_test3
                         )

# save data
write.csv(submission, "submission-bidhari_nn.csv", row.names = F)

# check first 3 data
head(submission, 3)
##     id strength
## 1 S826 32.57458
## 2 S827 24.22050
## 3 S828 75.28716

Dengan model neural network yang dibuat, performa yang diharapkan tercapai yaitu:

  • MAE pada data validation mencapai < 7.5.
  • R-squared pada data validation mencapai > 65%.
  • MAE pada data test mencapai < 7.5.
  • R-squared pada data test mencapai > 65%.

Dengan pemodelan ini kita dapat memprediksi kekuatan kompresi beton berdasarkan komposisi campurannya.

Apabila performa belum tercapai, kita dapat melakukan tuning model dengan mengubah jumlah layer, jumlah nodes pada tiap hidden layer, menggunakan optimizer lain, tuning jumlah epoch dan batch size.