library(dplyr)        #wrangling
library(tidymodels)   #recipe
library(GGally)       #correlation check
library(rsample)      #initial_split
library(MLmetrics)    #error check
library(lmtest)       #heteroscedasticity
library(car)          #multicolinearity/vif

library(randomForest)
library(lime)
library(ranger)

Berangkat dari masalah yang dihadapi para insinyur teknik sipil, mereka ingin mendapatkan ramuan atau racikan terbaik untuk bisa memperoleh beton berkualitas tinggi. Nyatanya beton sendiri tersusun dari berbagai material dengan porsi yang beragam. Maka dari itu, kita akan melakukan analisis tentang resep beton terbaik dengan menggunakan dataset yang telah disediakan.

test <- read.csv("data_concreteanalysis/data-test.csv")
train <- read.csv("data_concreteanalysis/data-train.csv")
head(train)

Data terdiri dari beberapa variabel di antaranya:

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.

summary(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

Data Preprocess

Tahap persiapan data mencakup transformasi tipe data, cek nilai NA, outlier, termasuk juga melakukan scaling. Mengapa dilakukan scaling? Karena berdasarkan summary(train), range dari variabel numerik pada data tersebut punya nilai yang beragam.

Misalnya nilai max di coarse_agg hingga ribuan sedangkan super_plast hanya puluhan. Variabel yang perlu dilakukan scaling adalah seluruh variabel numerik, baik target maupun prediktor.

Sebelum melakukan scaling, kita akan drop variabel yang dianggap sama sekali tidak punya pengaruh/hubungan terhadap target yaitu kolom id dengan menggunakan select() dari library(dplyr).

# drop kolom id
train <- train %>% 
  select(-1)
test <- test %>% 
  select(-1)

# cek NA
anyNA(train)

## [1] FALSE

Cross Validation & Scaling

Kita akan melakukan scaling terhadap variabel-variabel prediktor yang sudah linear terhadap target. Namun satu dari 8 prediktor yang ada, hubungannya tidak linear dengan strength, yaitu variabel age. Maka dari itu, di tahap preprocessing data ini kita juga menggunakan fungsi log untuk mengubah data non-linear menjadi linear dengan target. Kemudian setelah itu akan discaling dengan bantuan library(recipes) dan masuk ke tahap cross validation, membagi data train menjadi datatrain dan datavalidation menggunakan initial_split. Proporsi yang digunakan adalah 80% untuk datatrain dan 20% untuk datavalidation.

pre.train_log <- train %>% mutate_if(~is.integer(.), ~log10(.))

# sampling, recipe, recipe revert
RNGkind(sample.kind = "Rounding")
set.seed(100)
idx <- initial_split(data = pre.train_log, prop = 0.8)

rec <- recipe(strength ~ ., data = training(idx)) %>% 
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>% 
  step_scale(all_numeric()) %>% 
  prep()

rec_rev <- function(x, rec){

  means <- rec$steps[[2]]$means[["strength"]]
  sds <- rec$steps[[3]]$sds[["strength"]]

  x <- (x * sds + means) ^ 2

  x
}

# Cross Validation:

# data train & validation TANPA scaling
datatrain <- training(idx)
datavalidation <- testing(idx)

# train & validation DENGAN scaling
train_scaled <- juice(rec)
validation_scaled <- bake(rec, testing(idx))

#range data train sebelum dan setelah di-scale
range(train)

## [1]    0 1145

range(train_scaled)

## [1] -4.695796  2.511334

Range dari data yang sudah di-scale berada pada kisaran -4.7 hingga 2.6. Berikutnya, untuk memastikan jumlah observasi dari tiap data yang nantinya akan kita gunakan, bisa dilihat dengan fungsi dim(). Keterangan:
- train : data historis secara keseluruhan.
- datatrain: hasil splitting 80% dari train, dipakai agar model yang nanti dibuat bisa belajar dengan baik.
- datavalidation : hasil splitting 20% dari train untuk mengukur performa/kebaikan model setelah belajar dari datatrain.
- train_scaled: datatrain yang telah dilakukan scaling agar variabel-variabel numerik berada pada range yang seimbang.
- validation_scaled: datavalidation yang telah dilakukan scaling agar variabel-variabel numerik berada pada range yang seimbang.
Kedua data yang sudah di-scaling akan digunakan untuk membuat model dan prediksi.

dim(train)

## [1] 825   9

dim(datatrain)

## [1] 661   9

dim(datavalidation)

## [1] 164   9

dim(train_scaled)

## [1] 661   9

dim(validation_scaled)

## [1] 164   9

# cek outlier
sort(boxplot(train_scaled$strength, plot = T)$out)

## [1] -2.986445

Setelah dilakukan scaling dan pengecekan outlier, kita bisa ketahui bahwa ada 1 data outlier pada target kita dengan nilai (scaled) -2.947195. Walau ada outlier, data kita bisa dikatakan cukup rapi dan keadaan ini bisa ditolerir.

Exploratory Data Analysis

Berikutnya kita akan memilah variabel mana saja yang bisa secara signifikan punya hubungan/korelasi dengan target variabel kita, strength. Metode yang bisa digunakan adalah dengan menggunakan visualisasi dari library(GGally).

ggcorr(train_scaled, label = T, color = "black", label_size = 3, hjust = 1) +
  labs(title = "Correlation Matrix",
       subtitle = "on Strength")

Hasil visualisasi di atas menjelaskan bahwa pada data train yang sudah di-scale, ada 2 variabel prediktor yang punya korelasi cukup kuat positif denganstrength, yaitu age dan cement serta 2 prediktor lain yang berkorelasi lemah positif, yakni slag dan super_plast. Sedangkan variabel-variabel prediktor lain punya nilai korelasi lemah negatif terhadap strength.

Berbeda dengan korelasi prediktor terhadap target, dari plot di atas kita bisa tahu bahwa korelasi antar variabel prediktor justru cenderung minim bahkan bernilai 0. Hal ini bisa menjelaskan bahwa tidak ada multikolinearitas antar variabel prediktor: age, cement, super_plast, dan slag. Keadaan ini akan dibuktikan melalui uji 4 asumsi regresi linear terhadap model yang dibuat.

Kita juga dapat melihat distribusi data dari tiap variabel lewat visualisasi boxplot maupun bar. Visualisasi akan mempermudah kita dalam melihat sebaran dan range dari tiap variabel, misalnya pada boxplot di bawah. Terlihat bagaimana variabel dalam data train begitu beragam dengan range yang sangat berbeda jauh. Maka dari itu tahap scaling data menjadi tahap krusial yang harus dilakukan di bagian persiapan data sebelum dieksekusi lebih jauh. Data yang sudah di-scaling dapat menghasilkan boxplot pada gambar kedua.

boxplot(datatrain, horizontal = TRUE)

# scaled data train
boxplot(train_scaled, horizontal = TRUE)

Secara visual terlihat bahwa sebaran data dari tiap variabel di train_scaled sudah berada pada range yang serupa, antara -5 hingga 5. Untuk lebih jelasnya, bisa dilihat melalui visualisasi bar plot berikut:

train_scaled2 <- as.data.frame(train_scaled)

par(mfrow=c(3, 3))
colnames <- dimnames(train_scaled2)[[2]]
for (i in 1:9) {
    hist(train_scaled2[,i], xlim=c(-10, 10),
         breaks= 2,
         main=colnames[i],
         probability=TRUE, col="gray", border="white")
}

Model Fitting and Interpretation

Kita akan membuat model regresi linear dengan metode stepwise untuk memilih variabel prediktor mana yang mungkin berpengaruh signifikan terhadap target. Metode ini akan mengembalikan model dengan beberapa variabel prediktor dengan tingkat AIC (information loss) terkecil.

Dengan mengetahui variabel mana saja yang berpengaruh secara signifikan, maka kita bisa memperoleh ‘racikan’ yang dianggap paling efisien dengan mengabaikan variabel-variabel yang dianggap tidak berkontribusi terhadap target.

model_lm <- stats::step(lm(strength ~., data = train_scaled), direction = "backward", trace = 0)
summary(model_lm)

## 
## Call:
## lm(formula = strength ~ cement + slag + flyash + water + super_plast + 
##     coarse_agg + fine_agg + age, data = train_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.28172 -0.27220  0.01624  0.27621  1.24518 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.614e-16  1.525e-02   0.000  1.00000    
## cement       7.613e-01  3.308e-02  23.012  < 2e-16 ***
## slag         5.072e-01  3.331e-02  15.225  < 2e-16 ***
## flyash       2.645e-01  3.243e-02   8.158 1.75e-15 ***
## water       -1.046e-01  3.302e-02  -3.167  0.00161 ** 
## super_plast  1.201e-01  2.893e-02   4.151 3.75e-05 ***
## coarse_agg   1.442e-01  2.760e-02   5.226 2.34e-07 ***
## fine_agg     9.023e-02  3.204e-02   2.816  0.00501 ** 
## age          6.506e-01  1.553e-02  41.902  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3921 on 652 degrees of freedom
## Multiple R-squared:  0.8482, Adjusted R-squared:  0.8463 
## F-statistic: 455.2 on 8 and 652 DF,  p-value: < 2.2e-16

Interpretasi:

Berdasarkan summary model_lm, dapat diketahui bahwa ada 8 variabel prediktor yang dianggap berkontribusi secara signifikan terhadap strength, dibuktikan dengan nilai p-value yang kurang dari \(2.2*10^{-16}\). Nilai Adjusted R-squared bisa diinterpretasikan bahwa kedelapan prediktor tersebut mampu menjelaskan 84.63% dari nilai strength. Namun jika ditelaah satu-persatu, ada 2 variabel yang signifikansinya tidak setinggi variabel prediktor lain, yaitu water dan fine_agg. Walau begitu, kedua prediktor tersebut tetap memiliki kontribusi terhadap kenaikan adj.R-squared yang menunjukkan seberapa besar variabel prediktor mampu menjelaskan variabel target.

Jika dituliskan secara matematis, kontribusi tiap variabel bisa dilihat dari nilai koefisiennya:

\[ strength = (-2.614 * 10^{-16}) + (7.613*10^{-1}.cement) + (5.072*10^{-1}.slag) + (2.645*10^{-1}.flyash) + (-1.046*10^{-1}.water) + (1.201*10^{-1}.super_plast) + (1.442*10^{-1}.coarse_agg) + (9.023*10^{-1}.fine_agg)+ (6.506*10^{-1}.age) \] Artinya nilai intercept/konstanta/variabel strength sebesar\(-2.614 * 10^{-16}\) saat semua prediktornya bernilai 0.Sedangkan koefisien pada tiap variabel menjelaskan seberapa besar kontribusinya terhadap nilai strength.

Misalnya cement memiliki kontribusi sebesar \(7.613*10^{-1}\) untuk tiap pertambahan satu satuan terhadap kenaikan strength, dengan estimasi kesalahan/standard error kurang lebih \(3.308*10^{-2}\). Begitu pula dengan variabelslagmemiliki kontribusi sebesar \(5.072*10^{-1}\) untuk tiap pertambahan satu satuan terhadap kenaikan variabel target, dengan estimasi kesalahan/standard error kurang lebih \(3.331*10^{-2}\).

Lantas bagaimana dengan koefisien yang bernilai negatif seperti water? Interpretasinya sama dengan koefisien dari dua variabel sebelumnya, bedanya prediktor ini justru berkontribusi terhadap penurunan nilai variabel strength sebesar \(1.046*10^{-1}\) untuk setiap pertambahan satu satuan water. Lalu kenapa tetap dipakai sebagai prediktor? Jika dikaitkan dengan realitanya, air merupakan bahan penting yang digunakan dalam campuran beton. Fungsinya sebagai pemicu reaksi kimiawi dan perekat/pengikat agregat. Jadi meskipun secara angka bernilai negatif, pada praktiknya bahan ini sangatlah krusial dalam proses pembuatan beton.

Prediction and Evaluation

Setelah melihat kontribusi tiap prediktor dan memahami signifikansi terhadap target secara keseluruhan, kita akan menguji model_lm yang sudah dibuat untuk melihat apakah model ini bisa memprediksi nilai strength dari data validation_scaled dengan baik.

predict_lm <- predict(model_lm, newdata = data.frame(validation_scaled))

Ada satu langkah penting yang perlu dilakukan dalam proses evaluasi kebaikan model kita. Ketika mengukur kebaikan model dengan melihat MAE, kita perlu ingat bahwa data yang dipakai haruslah data mentah yang belum dilakukan transformasi (log/scale). Jadi sangat penting untuk mengembalikan hasil prediksi yang datanya sudah ditransformasi, menjadi data semula. Untuk melakukan ini, kita bisa memakai fungsi rev_rec() yang sudah kita buat di atas.

Revert predict_lm agar rangenya sesuai dengan data asli:

predict_lm_revert <- rec_rev(predict_lm, rec)

Setelah mengembalikan nilai prediksi dari tiap variabel ke nilai aslinya, barulah kita bisa mengukur kebaikan model dari MAE (mean absolute error). MAE diperoleh dari nilai absolut error yang dirata-ratakan.

# Prediction Performance
summary(model_lm)$adj

## [1] 0.8462938

MAE(y_pred = predict_lm_revert, y_true = datavalidation$strength)

## [1] 4.973242

*R-squared in (your own) validation dataset reach > 65%. MAE in (your own) validation dataset reach < 7.5.*

Interpretasi:

Sebagaimana yang telah dipaparkan di atas, bahwa nilai adjusted R-squared pada model_lm berarti 8 variabel prediktor mampu menjelaskan sebesar 84.62% dari keseluruhan nilai strength. Evaluasi model_lm dengan MAE menghasilkan angka 4.9 dari batas maksimal 7.5. MAE sebesar 4.973 berarti nilai prediksi akan berada kurang/lebih 4.973 dari nilai aslinya.

Contoh, pada hasil prediksi di observasi pertama sebesar 34.183659. Dengan error 4.973, maka nilai aslinya adalah 34.183659 -/+ 4.973 yaitu 29.21066 atau 39.15666. Dengan nilai adj.R-squared dan MAE sebesar ini, bisa dikatakan bahwa model_lm sudah cukup baik dalam belajar dan latihan memprediksi data validation.

Model Evaluation : 4 Limitations

Ada 4 asumsi untuk mengukur kebaikan model regresi linear yang telah kita buat:

Linearity : `cor.test()`

Model yang baik adalah model dengan variabel target dan prediktor yang memiliki hubungan linear.
H0 : Tidak linear
H1 : Linear
p-value < 0.05, tolak H0
* Target : Mencari p-value < 0.05 agar tolak H0, sehingga didapatkan linear.

cor.test(datatrain$cement, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$cement and datatrain$strength
## t = 14.084, df = 659, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4201430 0.5375336
## sample estimates:
##       cor 
## 0.4809912

cor.test(datatrain$water, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$water and datatrain$strength
## t = -7.6618, df = 659, p-value = 6.575e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3545206 -0.2144107
## sample estimates:
##        cor 
## -0.2859935

cor.test(datatrain$slag, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$slag and datatrain$strength
## t = 3.9344, df = 659, p-value = 9.224e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07611567 0.22515316
## sample estimates:
##       cor 
## 0.1514953

cor.test(datatrain$flyash, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$flyash and datatrain$strength
## t = -2.635, df = 659, p-value = 0.00861
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.17699098 -0.02605404
## sample estimates:
##        cor 
## -0.1021102

cor.test(datatrain$super_plast, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$super_plast and datatrain$strength
## t = 9.7859, df = 659, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2877598 0.4210247
## sample estimates:
##       cor 
## 0.3562022

cor.test(datatrain$coarse_agg, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$coarse_agg and datatrain$strength
## t = -3.1855, df = 659, p-value = 0.001513
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.19755014 -0.04733181
## sample estimates:
##        cor 
## -0.1231463

cor.test(datatrain$age, datatrain$strength)

## 
##  Pearson's product-moment correlation
## 
## data:  datatrain$age and datatrain$strength
## t = 16.783, df = 659, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4914599 0.5984944
## sample estimates:
##       cor 
## 0.5472104

Nilai p-value < 0.05 berarti tolak H0, atau bisa dinyatakan bahwa variabel prediktor dari model_lm kita punya hubungan linear dengan target variabel.

Normality of Residuals : `shapiro.test()`

Model regresi linear menunjukkan distribusi residual/error yang normal.
H0 : Residual berdistribusi normal
H1 : Residual berdistribusi tidak normal
p-value < 0.05, tolak H0
* Target : Mencari p-value > 0.05 agar gagal tolak H0, sehingga didapatkan residual berdistribusi normal.

shapiro.test(model_lm$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_lm$residuals
## W = 0.9984, p-value = 0.8131

p-value = 0.8131, maka gagal tolak H0. Artinya residual/error pada model_lm berdistribusi normal atau range nilai error tidak terlalu jauh dan berkumpul di area distribusi normal.

Heteroskedastisitas : `bptest()`

Kondisi yang harus dihindari, dimana variansi error dari model cenderung memiliki pola. Model yang ideal adalah model yang homoskedastisitas yakni errornya tidak berpola.
H0 : Homoscedasticity
H1 : Heteroscedasticity
p-value < 0.05, tolak H0
* Target : Mencari p-value > 0.05 agar gagal tolak H0, sehingga didapatkan kondisi Homoscedasticity.

bptest(model_lm)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_lm
## BP = 45.244, df = 8, p-value = 3.308e-07

p-value nilainya jauh lebih kecil dari 0.05, maka tolak H0 dan terima H1. Berarti dapat dikatakan model_lm memiliki error yang berpola/heteroscedasticity. Keadaan ini bisa disebabkan oleh beberapa hal, misalnya karena nilai outlier, nilai error yang besar, ataupun beberapa prediktor yang tidak punya hubungan linear dengan target.

Multikolinearitas : `vif()`

Tiap variabel prediktor pada model, diharapkan tidak saling memiliki korelasi satu sama lain. Jika nilai VIF(Variance Inflation Factor) mencapai atau lebih dari 10, maka bisa dinyatakan variabel tersebut memiliki multikolinearitas.

vif(model_lm)

##      cement        slag      flyash       water super_plast  coarse_agg 
##    4.699001    4.765616    4.514625    4.681496    3.593003    3.270860 
##    fine_agg         age 
##    4.409176    1.035187

Nilai VIF dari model_lm tidak ada yang mencapai 10, berarti tidak ada multikolinearitas dari variabel-variabel prediktor dalam model.

Kesimpulannya, dari uji 4 asumsi ini menunjukkan bahwa model_lm dengan 8 prediktor berhasil memenuhi 3 asumsi: linearity, normality, dan multicolinearity. Meski tidak memenuhi standar homoscedasticity, performa model_lm sudah cukup baik dan memenuhi standar, yakni dengan adj.R-squared 84.62% dan MAE sebesar 4.973.

Model Improvement

Pada tahap ini, kita akan melakukan tuning pada model_lm untuk meningkatkan performanya dan memperkecil error. Tuning/improvement model dilakukan dengan mentransform beberapa variabel target menggunakan sqrt maupun log10, juga menghapus data outlier. Berikut transformasi variabel prediktor yang dilakukan:

pre.train_logg <- train %>%
  mutate(age = log10(age),
         slag = slag^2,
         water = sqrt(water),
         super_plast = sqrt(super_plast),
         fine_agg = log10(fine_agg),
         coarse_agg = log10(coarse_agg))

Mengapa variabel target tidak ikut ditransformasi? Karena setelah dicoba, performa yang dihasilkan malah cenderung turun Maka tranformasi variabel target tidak dilakukan.

RNGkind(sample.kind = "Rounding")
set.seed(100)
coba_idx <- initial_split(data = pre.train_logg, prop = 0.8)

coba_rec <- recipe(strength ~ ., data = training(coba_idx)) %>% 
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>% 
  step_scale(all_numeric()) %>% 
  prep()

coba_rec_rev <- function(x, coba_rec){

  means <- rec$steps[[2]]$means[["strength"]]
  sds <- rec$steps[[3]]$sds[["strength"]]

  x <- (x * sds + means) ^ 2

  x
}

# Cross Validation:

# data train & validation TANPA scaling
coba_datatrain <- training(coba_idx)
coba_datavalidation <- testing(coba_idx)

# train & validation DENGAN scaling
coba_train_scaled <- juice(coba_rec)
coba_validation_scaled <- bake(coba_rec, testing(coba_idx))

# drop outlier
coba_outliers <- boxplot(coba_train_scaled$strength, plot = F)$out
coba_train_scaled <- coba_train_scaled[-which(coba_train_scaled$strength %in% coba_outliers),]

# modeling
coba_model <- stats::step(lm(strength ~., data = coba_train_scaled), direction = "backward", trace = 0)
summary(coba_model)$adj

## [1] 0.8499894

bptest(coba_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  coba_model
## BP = 42.241, df = 8, p-value = 1.22e-06

Adj.R-squared model meningkat sekitar 0.3% dengan asumsi heteroscedasticity yang masih belum bisa terpenuhi. Berikutnya kita ukur kebaikan model dalam melakukan prediksi terhadap data validation.

#Prediction

coba_predict <- predict(coba_model, newdata = data.frame(coba_validation_scaled))

# revert
coba_predict_lm_revert <- coba_rec_rev(coba_predict, coba_rec)
MAE(y_pred = coba_predict_lm_revert, y_true = datavalidation$strength)

## [1] 4.778795

Melalui proses tuning model dengan cara transformasi beberapa variabel prediktor, diperoleh hasil akhir adj.R-squared 84.99% dan MAE sebesar 4.778.

Predicting Data Test

Trial & error yang sudah dilakukan dalam proses training model, sekarang saatnya menguji model untuk melakukan prediksi terhadap nilai strength di data test. Model yang nanti dibuat akan belajar dari data train keseluruhan (tanpa di-split). Tujuannya supaya model kita bisa melihat data/observasi lebih banyak dan belajar dengan baik.

head(test)

# preprocess/tranform variable
train_log <- train %>% mutate_if(~is.integer(.), ~log10(.)) %>% 
  mutate(slag = slag^2,
         water = sqrt(water),
         super_plast = sqrt(super_plast),
         fine_agg = log10(fine_agg),
         coarse_agg = log10(coarse_agg))
test_log <- test %>% mutate_if(~is.integer(.), ~log10(.)) %>% 
  mutate(slag = slag^2,
         water = sqrt(water),
         super_plast = sqrt(super_plast),
         fine_agg = log10(fine_agg),
         coarse_agg = log10(coarse_agg))

# recipe & revert
recipe <- recipe(strength ~ ., data = train_log) %>% 
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>% 
  step_scale(all_numeric()) %>% 
  prep()

recipe_rev <- function(x, recipe){

  means <- rec$steps[[2]]$means[["strength"]]
  sds <- rec$steps[[3]]$sds[["strength"]]

  x <- (x * sds + means) ^ 2

  x
}

# scaling
train.sc <- juice(recipe)
test.sc <- bake(recipe, test_log)

# drop outlier
outliers <- boxplot(train.sc$strength, plot = F)$out
coba_train_scaled <- train.sc[-which(train.sc$strength %in% outliers),]

# modeling
model <- stats::step(lm(strength ~., data = train.sc), direction = "backward", trace = 0)
summary(model)$adj

## [1] 0.8547752

# predicting
predict <- predict(model, newdata = data.frame(test.sc))

# revert
predict_lm_revert <- recipe_rev(predict, recipe)
test$strength <- predict_lm_revert

test_final <- read.csv("data_concreteanalysis/data-test.csv")

# create new csv
id <- test_final %>% 
  select(1)
prediction <- test %>% 
  select(9)
bind <- cbind(id, prediction)
write.csv(bind, "prediction3.csv", row.names = FALSE)

Performa modeldalam memprediksi strength di data test sudah baik dan memenuhi standar yang diminta. Kebaikannya diukur dari nilai adj.R-squared 82% dan MAE sebesar 5.57, dengan kontribusi tiap variabel/komponen dalam adonan beton sebagai berikut:

summary(model)

## 
## Call:
## lm(formula = strength ~ cement + slag + flyash + water + super_plast + 
##     coarse_agg + fine_agg + age, data = train.sc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.34981 -0.26871 -0.00234  0.26520  1.22981 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.932e-16  1.327e-02   0.000 1.000000    
## cement       8.160e-01  3.101e-02  26.313  < 2e-16 ***
## slag         5.384e-01  3.111e-02  17.307  < 2e-16 ***
## flyash       3.052e-01  3.204e-02   9.527  < 2e-16 ***
## water       -1.032e-01  2.661e-02  -3.879 0.000113 ***
## super_plast  1.299e-01  2.485e-02   5.227 2.19e-07 ***
## coarse_agg   1.408e-01  2.321e-02   6.068 1.98e-09 ***
## fine_agg     1.317e-01  2.846e-02   4.627 4.31e-06 ***
## age          6.577e-01  1.358e-02  48.426  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3811 on 816 degrees of freedom
## Multiple R-squared:  0.8562, Adjusted R-squared:  0.8548 
## F-statistic: 607.2 on 8 and 816 DF,  p-value: < 2.2e-16

\[ strength = (1.932 * 10^{-16}) + (8.160*10^{-1}.cement) + (5.384*10^{-1}.slag) + (3.052*10^{-1}.flyash) + (-1.032*10^{-1}.water) + (1.299*10^{-1}.super_plast) + (1.408*10^{-1}.coarse_agg) + (1.317*10^{-1}.fine_agg)+ (6.577*10^{-1}.age) \]

Finding the Right Material Composition

Pada tahap ini, kita akan membandingkan kekuatan beton dengan racikan yang berbeda. Titik beratnya berada variabel prediktor yang memiliki korelasi terkuat dengan target, yaitu age yang menjelaskan seberapa lama adonan beton didiamkan. Ada pula variabel flyash yang nilai berkisar dari 0 hingga 200, menjadi variabel yang menarik untuk diberi label. Apakah dengan ada/tidaknya flyash mempengaruhi strength?
Lebih dulu, kita perlu membuat label sesuai dengan nilai age, strength, dan flyash:

train2 <- train

#3 class of age (short, average, long)
train2 <- train2 %>% 
  mutate(age.label = as.factor
         (case_when(age <= quantile(train$age, 0.25) ~ "short",
                    age > quantile(train$age, 0.25) & train$age < quantile(train$age, 0.75) ~ "average",
                    TRUE ~ "long")))

#4 class of strength (weak, average, strong, very strong)
train2 <- train2 %>% 
 mutate(strength.label = as.factor
         (case_when(strength <= quantile(train$strength, 0.25) ~ "weak",
                    strength > quantile(train$strength, 0.25) & strength <= quantile(train$strength, 0.5) ~ "average",
                    strength > quantile(train$strength, 0.5) & strength <= quantile(train$strength, 0.75) ~ "strong",
                    TRUE ~ "very strong")))

#2 class of flyash (none, exist)
train2$flyash.label <- ifelse(train2$flyash == 0, "none", "exist")
head(train2)

Selanjutnya, kita akan melakukan uji statistik untuk mencari perbedaan rerata yaitu dengan ANOVA.

aov <- aov(strength ~ age.label + flyash.label, data = train2)
summary(aov)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## age.label      2  73647   36823  196.48  < 2e-16 ***
## flyash.label   1   2692    2692   14.37 0.000162 ***
## Residuals    821 153872     187                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(aov, 1)

Berdasarkan summary(aov) diperoleh bahwa kedua variabel kategorik berupa label age dan flyash signifikan berpengaruh terhadap nilai strength. Hal ini dibuktikan dengan tanda ’***’ pada kedua prediktor.

Conclusion

Model regresi linear secara statistik bisa dinyatakan tepat dan baik untuk melakukan prediksi nilai strength pada data test. Awalnya model belajar terlebih dahulu dengan memanfaatkan data train saja–yang kemudian dibagi menjadi train dan validation. Barulah model dibuat serta diuji dengan memanfaatkan 100% data train dan diukur kebaikannya saat memprediksi data test. Akurasi adj.R-squared sebesar 82% dan MAE sebesar 5.57 menandakan performa yang lebih dari cukup.

Maka dari itu bisa disimpulkan bahwa machine learning dengan model regresi linear sangat tepat digunakan untuk memprediksi variabel target numerik, apalagi saat variabel prediktornya juga numerik. Penggunaan model ini bisa diaplikasikan pada kasus bisnis lain, misalnya prediksi omset perusahaan, harga jual properti/rumah, usia mesin (dalam satuan minggu/bulan), dan sebagainya.

Extra : Random Forest

Jika regresi model unggul di bagian interpretasi, random forest adalah salah satu model machine learning yang robust. Model ini memiliki akurasi yang tinggi, tapi punya perhitungan tersendiri di baliknya sehingga agak sulit untuk diinterpretasikan. Di sisi lain, ada LIME dari library (lime) untuk membantu kita memahami variabel prediktor mana saja yang punya kontribusi paling besar dalam mempengaruhi target, namun kita tidak bisa tahu bobot kontribusi dari tiap prediktor.

# recipe & revert
recrev <- function(x){

  y <- x %>% select_if(is.numeric)

  for (i in 1:length(names(y))) {
    y[ , i] <- y[ ,i] * rec$steps[[3]]$sds[names(y)[i]]
  }

  x <- x %>% select_if(is.factor) %>% bind_cols(y)
  return(x)
}

# modeling
model_rf <- randomForest(
  formula =  strength ~ .,
  data = datatrain,
  importance=TRUE,
  ntree = 500)

# prediction
predrf_toval <- predict(model_rf, newdata = data.frame(datavalidation))

# MAE
MAE(predrf_toval, datavalidation$strength)

## [1] 4.406124

select(datavalidation, strength) %>%
  mutate(
    pred_lm = predict(model_rf, datavalidation)
    # pred_lm = rec_rev(pred_lm, rec)
  ) %>%
  summarise(
    mae = mae_vec(truth = strength, estimate = pred_lm),
    rsq = rsq_vec(truth = strength, estimate = pred_lm),
  )

Nilai adj.R-squared dan MAE dengan random forest lebih unggul beberapa persen dari model linear regresi. Namun model random forest ini tidak bisa dengan mudah diinterpretasikan seperti regresi.

Visualisasi

#define model spec
model_spec <- rand_forest(
  mode = "regression",
  mtry = 2,
  trees = 500,
  min_n = 1)

#define model engine
model_spec <- set_engine(model_spec,
                         engine = "ranger",
                         seed = 123,
                         num.threads = parallel::detectCores(),
                         importance = "impurity")

Variabel prediktor apa saja yang paling penting?

#model fitting
set.seed(123)
model <- fit_xy(
  object = model_spec,
  x = select(datatrain, -strength),
  y = select(datatrain, strength)
)

# get variable importance
var_imp <- tidy(model$fit$variable.importance) %>%
  arrange(desc(x))

# tidying
var_imp <- var_imp %>%
  head(10) %>%
  rename(variable = names, importance = x) %>%
  mutate(variable = reorder(variable, importance))

# variable importance plot
ggplot(var_imp, aes(x = variable, y = importance)) +
  geom_col(aes(fill = importance), show.legend = F) +
  geom_text(aes(label = round(importance, 2)), nudge_y = 1)+
  coord_flip() +
  labs(title = "Variables Importance (Top 10)", x = NULL, y = NULL, fill = NULL) +
  scale_y_continuous(expand = expand_scale(mult = c(0, 0.1))) +
  scale_fill_viridis_c()+
  theme_minimal()

set.seed(123)
explainer <- lime(x = recrev(datatrain) %>% select(-strength), model = model)
explanation <- explain(x = recrev(datavalidation) %>% select(-strength) %>% slice(1:4),
                      n_permutations = 500,
                       explainer = explainer,
                       n_features = 10)
plot_features(explanation)

Interpretasi: Visualisasi grafik pertama menjelaskan bahwa variabel age punya pengaruh paling besar terhadap strength. Kemudian ada cement dan water di posisi dua dan tiga, namun sayang kita tidak tahu bobot dari masing-masing prediktor. Pada visualisasi kedua, tampak lebih jelas ada prediktor yang berkontribusi negatif terhadap strength, misalnya saat super_plast kurang dari 10.2. Sedangkan saat age nilainya lebih besar dari 194.7 maka variabel ini berkontribusi positif bagi strength. Jadi model random forest tepat digunakan saat ingin menangani kasus yang mengedepankan hasil prediksi ketimbang interpretasi model.

Secret Recipe of the Great Concrete

Wenny Kezia Jo

Aug, 2020

Data Preprocess

Cross Validation & Scaling

Exploratory Data Analysis

Model Fitting and Interpretation

Prediction and Evaluation

Model Evaluation : 4 Limitations

Linearity : `cor.test()`

Normality of Residuals : `shapiro.test()`

Heteroskedastisitas : `bptest()`

Multikolinearitas : `vif()`

Model Improvement

Predicting Data Test

Finding the Right Material Composition

Conclusion

Extra : Random Forest

Visualisasi

Secret Recipe of the Great Concrete

Wenny Kezia Jo

Aug, 2020

Data Preprocess

Cross Validation & Scaling

Exploratory Data Analysis

Model Fitting and Interpretation

Prediction and Evaluation

Model Evaluation : 4 Limitations

Linearity : cor.test()

Normality of Residuals : shapiro.test()

Heteroskedastisitas : bptest()

Multikolinearitas : vif()

Model Improvement

Predicting Data Test

Finding the Right Material Composition

Conclusion

Extra : Random Forest

Visualisasi

Linearity : `cor.test()`

Normality of Residuals : `shapiro.test()`

Heteroskedastisitas : `bptest()`

Multikolinearitas : `vif()`