Cara 1

custom <- trainControl(method = "repeatedcv",
                       number = 10,
                       repeats = 5,
                       verboseIter = F)

Langkah pertama adalah mengkustom control parameters dengan sepuluh lipatan yang diulang sebanyak lima kali. Logical verboseIter digunakan untuk mencetak sebuah training log atau mencetak iterasi.

#fitting  Ridge Regression model
set.seed(1234)
ridge <- train(mpg~.,train,
               method = "glmnet",
               tuneGrid = expand.grid(alpha = 0, 
                                      lambda = seq(0.0001, 1, length = 10)),
               trControl = custom)
ridge

## glmnet 
## 
## 313 samples
##   5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 283, 283, 282, 281, 281, 281, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0.0001  4.247618  0.6867413  3.262458
##   0.1112  4.247618  0.6867413  3.262458
##   0.2223  4.247618  0.6867413  3.262458
##   0.3334  4.247618  0.6867413  3.262458
##   0.4445  4.247618  0.6867413  3.262458
##   0.5556  4.247618  0.6867413  3.262458
##   0.6667  4.248894  0.6866187  3.262537
##   0.7778  4.251497  0.6863436  3.262971
##   0.8889  4.254314  0.6860585  3.263812
##   1.0000  4.257276  0.6857711  3.264911
## 
## Tuning parameter 'alpha' was held constant at a value of 0
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0 and lambda = 0.5556.

Model regresi Ridge dibangun dengan data training dengan alpha = 0 dalam fungsi tune.grid yang digunakan untuk mengatur nilai lambda. Lambda yang dikeluarkan diatur secara sekuensial dari 0.0001 sampai dengan 1 dengan panjang 10. Output di atas membentuk model regresi berbasis Ridge regression dari beberapa jenis lambda dengan cross validation sejumlah sepuluh lipatan yang diulang sebanyak lima kali. Lambda optimum yang didapatkan adalah 0.556.

#mean validation score
mean(ridge$resample$RMSE)

## [1] 4.247618

Mean validation score merupakan rata-rata dari 50 RMSE yang dikeluarkan, di mana angka 50 ini berasal dari sepuluh lipatan yang diulang sebanyak lima kali. Didapatkan rata-ratanya berjumlah 4.247618.

#plotting the model
plot1 <- plot(ridge, main = "Ridge Regression")

#plotting important variables
plot2 <- plot(varImp(ridge, scale = F), main = "Variables Importance")

plot_grid(plot1, plot2)

Plot pertama menunjukkan bahwa RMSE minimum terjadi pada saat regularization parameter berada di sekitar angka 0.5 sehingga diperoleh lambda optimum di sekitar angka tersebut, yaitu 0.556 (berdasarkan fitting model). Plot kedua menunjukkan variables importance dari kelima peubah prediktor. Terlihat bahwa peubah yang paling berpengaruh terhadap peubah respon mpg adalah peubah cylinders, di mana skor kepentingan peubah ini berbeda signifikan dengan keempat peubah prediktor lainnya. Silinder atau cylinder block adalah komponen yang sangat penting dalam rangkaian mesin mobil. Komponen ini merupakan penopang untuk komponen mesin lainnya. Saat silinder bekerja dengan optimal, mesin mobil akan memiliki performa terbaik. Oleh karena itu, wajar jika komponen ini merupakan peubah yang paling berpengaruh terhadap peubah mpg, yaitu jarak yang ditempuh kendaraan per galon bahan bakarnya.

Cara 2

lambdas <- 10^seq(2,-3, by = -.1)
ridge_reg = glmnet(as.matrix(train[,-1]), train[,1], nlambda = 25, alpha = 0, family = 'gaussian', lambda = lambdas)
summary(ridge_reg)

##           Length Class     Mode   
## a0         51    -none-    numeric
## beta      255    dgCMatrix S4     
## df         51    -none-    numeric
## dim         2    -none-    numeric
## lambda     51    -none-    numeric
## dev.ratio  51    -none-    numeric
## nulldev     1    -none-    numeric
## npasses     1    -none-    numeric
## jerr        1    -none-    numeric
## offset      1    -none-    logical
## call        7    -none-    call   
## nobs        1    -none-    numeric

ridge_reg

## 
## Call:  glmnet::glmnet(x = as.matrix(train[, -1]), y = train[, 1], nlambda = 25,      alpha = 0, family = "gaussian", lambda = lambdas) 
## 
##    Df  %Dev  Lambda
## 1   5 26.70 100.000
## 2   5 30.96  79.430
## 3   5 35.39  63.100
## 4   5 39.87  50.120
## 5   5 44.25  39.810
## 6   5 48.37  31.620
## 7   5 52.10  25.120
## 8   5 55.37  19.950
## 9   5 58.14  15.850
## 10  5 60.42  12.590
## 11  5 62.24  10.000
## 12  5 63.67   7.943
## 13  5 64.78   6.310
## 14  5 65.62   5.012
## 15  5 66.27   3.981
## 16  5 66.77   3.162
## 17  5 67.16   2.512
## 18  5 67.46   1.995
## 19  5 67.70   1.585
## 20  5 67.89   1.259
## 21  5 68.04   1.000
## 22  5 68.16   0.794
## 23  5 68.25   0.631
## 24  5 68.33   0.501
## 25  5 68.38   0.398
## 26  5 68.42   0.316
## 27  5 68.46   0.251
## 28  5 68.48   0.200
## 29  5 68.50   0.158
## 30  5 68.51   0.126
## 31  5 68.51   0.100
## 32  5 68.52   0.079
## 33  5 68.52   0.063
## 34  5 68.53   0.050
## 35  5 68.53   0.040
## 36  5 68.53   0.032
## 37  5 68.53   0.025
## 38  5 68.53   0.020
## 39  5 68.53   0.016
## 40  5 68.53   0.013
## 41  5 68.53   0.010
## 42  5 68.53   0.008
## 43  5 68.53   0.006
## 44  5 68.53   0.005
## 45  5 68.53   0.004
## 46  5 68.53   0.003
## 47  5 68.53   0.003
## 48  5 68.53   0.002
## 49  5 68.53   0.002
## 50  5 68.53   0.001
## 51  5 68.53   0.001

Dilakukan setting nilai lambda yang akan dikeluarkan, yaitu 10^seq(2,-3, by = -.1) yang berjumlah 51 lambda dari nilai 100 sampai dengan 0.001. Selain itu, terdapat pula persentase deviasi dari masing-masing lambda.

set.seed(100)
#tunning parameter lambda dg cv
cv_ridge <- cv.glmnet(as.matrix(train[,-1]), train[,1], alpha = 0, lambda = lambdas)
optimal_lambda <- cv_ridge$lambda.min
optimal_lambda

## [1] 0.2511886

Didapatkan lambda minimum bernilai 0.2511886, lalu dari output ridge_reg sebelumnya dilihat lambda minimum tersebut (0.2512) paling mendekati nilai lambda yang mana. Berdasarkan output ridge_reg, didapatkan lambda = 0.251 pada indeks ke-27 dengan persentase deviasi adalah 68.46%.

coef(cv_ridge, s = "lambda.min")

## 6 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept)  45.656149183
## cylinders    -0.321199471
## displacement -0.010287495
## horsepower   -0.050784436
## weight       -0.003870524
## acceleration -0.106145478

Output di atas merupakan nilai intersep dan koefisien regresi dari masing-masing peubah prediktor yang selanjutnya digunakan untuk membuat persamaan model regresi Ridge.

# Compute R^2 from true and predicted values
eval_results <- function(true, predicted, df) {
  SSE <- sum((predicted - true)^2)
  SST <- sum((true - mean(true))^2)
  R_square <- 1 - SSE / SST
  RMSE = sqrt(SSE/nrow(df))
 
  # Model performance metrics
  data.frame(RMSE = RMSE,
             Rsquare = R_square)
}

# Prediction and evaluation on train data
predictions_train <- predict(ridge_reg, s = optimal_lambda, newx = as.matrix(train[,-1]))
eval_results(train[,1], predictions_train, train)

##       RMSE   Rsquare
## 1 4.217771 0.6845617

# Prediction and evaluation on test data
predictions_test <- predict(ridge_reg, s = optimal_lambda, newx = as.matrix(test[,-1]))
eval_results(test[,1], predictions_test, test)

##       RMSE   Rsquare
## 1 4.279487 0.7520407

Didapatkan nilai akurasi (RMSE) dan R-square dari data training dan testing seperti pada tabel di bawah ini. Terlihat bahwa nilai RMSE keduanya tidak berbeda jauh, yang mengindikasikan bahwa model yang dilatih pada tahapan training sudah cukup baik saat dilakukan validasi terhadap data testing.

Pemilihan Banyaknya Cluster

Penentuan cluster pada metode hirarki dapat menggunakan dendogram, lalu diperkuat dengan koefisien silhouette atau WSS. Tujuan melihat banyaknya cluster menggunakan dendogram adalah agar data dapat dipisahkan atau dikelompokkan berdasarkan karakteristik yang jelas.

Koefisien Sihouette

# Complete Linkage
plota <- fviz_nbclust(data.stdz, FUNcluster = hcut, method = "silhouette", hc_method = "complete", hc_metric="euclidean")

# Average Linkage
plotb <- fviz_nbclust(data.stdz, FUNcluster = hcut, method = "silhouette", hc_method = "average", hc_metric="euclidean")

plot_grid(plota, plotb)

# Centroid Linkage
plotc <- fviz_nbclust(data.stdz, FUNcluster = hcut, method = "silhouette", hc_method = "centroid", hc_metric="euclidean")

# Single Linkage
plotd <- fviz_nbclust(data.stdz, FUNcluster = hcut, method = "silhouette", hc_method = "single", hc_metric="euclidean")

plot_grid(plotc, plotd)

# Ward Linkage
fviz_nbclust(data.stdz, FUNcluster = hcut, method = "silhouette", hc_method = "ward.D", hc_metric="euclidean")

Hanya satu metode linkage yang menghasilkan koefisien silhouette tertinggi saat jumlah cluster sama dengan delapan (k = 8), yaitu single linkage. Sementara itu, keempat metode lainnya menghasilkan koefisien silhouette tertinggi saat jumlah cluster sama dengan dua (k = 2).

Ringkasan dari penentuan jumlah cluster dengan menggunakan koefisien silhouette pada beberapa metode penggerombolan dapat dilihat pada tabel di bawah.

Hierarchical Clustering dengan Metode Complete Linkage

Diterapkan hierarchical clustering dengan menggunakan metode complete linkage dengan banyak cluster adalah dua (k = 2). Metode ini dipilih karena hasil dendogram terlihat baik.

hc.auto <- eclust(df.auto, stand = TRUE, FUNcluster = "hclust", k = 2, hc_method = "complete", hc_metric = "euclidean", graph = F)

aggregate(df.auto, by=list(cluster=hc.auto$cluster), FUN = mean)

##   cluster      mpg cylinders displacement horsepower   weight acceleration
## 1       1 14.68400  7.980000     345.4700  160.40000 4121.560      12.7020
## 2       2 26.44658  4.613014     142.6798   85.31507 2585.812      16.5137

Output di atas menghasilkan nilai rata-rata keenam peubah pada setiap cluster, di mana terdapat dua cluster. Dapat disimpulkan bahwa cluster 1 berisi kendaraan dengan spesifikasi yang lebih tinggi dibandingkan dengan cluster 2. Hal ini ditunjukkan dengan nilai rata-rata beberapa peubah yang lebih tinggi, seperti blok silinder, kekuatan mesin, dan lainnya. Sementara itu, peubah mpg pada cluster 1 lebih rendah daripada cluster 2. Hal ini menunjukkan bahwa dengan banyaknya bahan bakar yang sama, kendaraan pada cluster 2 dapat menempuh jarak yang lebih jauh dibandingkan dengan kendaraan pada cluster 1.

fviz_cluster(hc.auto)

Dengan melihat cluster plot di atas, lebih banyak kendaraan yang berada di cluster 2 dibandingkan pada cluster 2. Kedua komponen utama dapat menjelaskan sekitar 91.9% keragaman pada data.

Pemilihan Banyaknya Cluster

Banyaknya cluster ditentukan dengan menggunakan kriteria statistik, yaitu koefisien silhouette dan Within Sum of Square (WSS).

fviz_nbclust(data.stdz, FUNcluster = kmeans, method = "silhouette")

fviz_nbclust(data.stdz, FUNcluster = kmeans, method = "wss")

K-Means

kmeans.auto <- eclust(df.auto, stand = TRUE, FUNcluster = "kmeans", k = 3, graph = F)

kmeans.auto$centers

##          mpg  cylinders displacement  horsepower     weight acceleration
## 1 -0.4653346  0.3347858    0.2281863 -0.06918306  0.2926493    0.3077339
## 2 -1.1424784  1.4699656    1.4679162  1.48955594  1.3731827   -1.0511897
## 3  0.7628541 -0.8600086   -0.8099465 -0.68343021 -0.7941342    0.3630997

Output di atas menghasilkan centroid untuk data yang masih terstandardisasi sehingga sulit diinterpretasikan.

aggregate(df.auto, by=list(cluster=kmeans.auto$cluster), FUN = mean)

##   cluster      mpg cylinders displacement horsepower   weight acceleration
## 1       1 19.81398  6.043011     218.2903  101.80645 3226.161     16.39032
## 2       2 14.52887  7.979381     348.0206  161.80412 4143.969     12.64124
## 3       3 29.40000  4.004950     109.6559   78.16337 2303.045     16.54307

Output di atas menghasilkan nilai rata-rata keenam peubah pada setiap cluster, di mana terdapat tiga cluster. Dapat disimpulkan bahwa cluster 2 berisi kendaraan dengan spesifikasi yang lebih tinggi dibandingkan dengan cluster 1 dan 3. Hal ini ditunjukkan dengan rata-rata beberapa peubah bernilai paling tinggi, seperti blok silinder, kekuatan mesin, dan lainnya. Sementara itu, peubah mpg pada cluster 2 lebih rendah daripada cluster 1 dan 3. Hal ini menunjukkan bahwa dengan banyaknya bahan bakar yang sama, kendaraan pada cluster 2 dapat menempuh jarak yang lebih sedikit dibandingkan dengan kendaraan pada cluster 1 dan 3.

fviz_cluster(kmeans.auto)

Dengan menggunakan kedua metode, yaitu hierarchical clustering dengan dua cluster dan K-Means clustering dengan tiga cluster, terjadi perbedaan yang signifikan terhadap hasil penggerombolan. Jumlah cluster yang lebih baik dalam mengelompokkan data adalah tiga, salah satunya dapat disesuaikan dengan data aktual peubah displacement yang bernilai 4, 6, dan 8. Selain itu, nilai rata-rata peubah pada setiap cluster dengan (k = 3) sangat berbeda satu sama lain. Pada kedua metode, tetap tidak terjadi overlap atau tumpang tindih antara individu pada satu cluster dengan cluster lainnya.

Menangani Missing Values

df.hitters1 <- na.omit(df.hitters)
tibble(df.hitters1)

## # A tibble: 263 × 20
##    AtBat  Hits HmRun  Runs   RBI Walks Years CAtBat CHits CHmRun CRuns  CRBI
##    <int> <int> <int> <int> <int> <int> <int>  <int> <int>  <int> <int> <int>
##  1   315    81     7    24    38    39    14   3449   835     69   321   414
##  2   479   130    18    66    72    76     3   1624   457     63   224   266
##  3   496   141    20    65    78    37    11   5628  1575    225   828   838
##  4   321    87    10    39    42    30     2    396   101     12    48    46
##  5   594   169     4    74    51    35    11   4408  1133     19   501   336
##  6   185    37     1    23     8    21     2    214    42      1    30     9
##  7   298    73     0    24    24     7     3    509   108      0    41    37
##  8   323    81     6    26    32     8     2    341    86      6    32    34
##  9   401    92    17    49    66    65    13   5206  1332    253   784   890
## 10   574   159    21   107    75    59    10   4631  1300     90   702   504
## # … with 253 more rows, and 8 more variables: CWalks <int>, League <fct>,
## #   Division <fct>, PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>,
## #   NewLeague <fct>

Dilakukan penghapusan terhadap baris dengan missing values sehingga didapatkan data bersih sebanyak 263 amatan.

plot_intro(df.hitters1)

Dilihat lagi metrik gambaran umum dari data bersih yang sudah tidak memiliki missing value, didapatkan hasil seperti pada plot di bawah. Memori yang terpakai lebih rendah daripada penggunaan memori pada data yang belum ditangani missing values-nya. Hal ini disebabkan oleh baris data yang berkurang setelah dilakukan penghapusan. Terlihat pula seluruh baris memiliki amatan yang lengkap.

skim_without_charts(df.hitters1)

Variable type: factor

Variable type: numeric

Berdasarkan output di atas, tidak perlu dilakukan penghapusan terhadap peubah kategorik pada data karena seluruh peubah kategorik memiliki lebih dari satu kategori. Selain itu, terdapat perubahan pada informasi peubah (rata-rata, standar deviasi, dan persentil) setelah dilakukan penanganan missing values pada data.

Memeriksa Sebaran Data

plot_histogram(data = df.hitters1, nrow=3, ncol = 3,
               geom_histogram_args = list(fill="#b8b894"))

Histogram di atas menggambarkan sebaran dari 17 peubah numerik. Terlihat sebaran sebagian besar peubah cenderung menjulur ke kanan. Artinya, frekuensi tertinggi dari peubah yang digunakan berada pada rentang nilai yang rendah.

par(mfrow=c(1,3))
pl1 <- ggplot(df.hitters1) +
 aes(x = League) +
 geom_bar(fill = "#b8b894") +
 labs(title = "League") +
 coord_flip() +
 theme_gray() +
 theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5))


pl2 <- ggplot(df.hitters1) +
 aes(x = Division) +
 geom_bar(fill = "#b8b894") +
 labs(title = "Division") +
 coord_flip() +
 theme_gray() +
 theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5))


pl3 <- ggplot(df.hitters1) +
 aes(x = NewLeague) +
 geom_bar(fill = "#b8b894") +
 labs(title = "NewLeague") +
 coord_flip() +
 theme_gray() +
 theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5))

plot_grid(pl1, pl2, pl3)

Diagram batang di atas merupakan besaran dari masing-masing peubah kategorik, yaitu league, division, dan new league.

Memeriksa Korelasi Peubah

cor_mat <- cor(df.hitters1%>%
                   select_if(is.numeric),method = "spearman")
cor_mat[upper.tri(cor_mat,diag = TRUE)] <- NA 
cor_df <- cor_mat   %>%
    as.data.frame() %>% 
    rownames_to_column(var = "Var1") %>%
  pivot_longer(names_to = "Var2",
               values_to = "corr",
               -Var1) %>% na.omit
cor_df %>% filter(abs(corr)>0.6) %>% arrange(desc(abs(corr)))

## # A tibble: 42 × 3
##    Var1   Var2    corr
##    <chr>  <chr>  <dbl>
##  1 CHits  CAtBat 0.997
##  2 CRuns  CHits  0.990
##  3 CRuns  CAtBat 0.989
##  4 Hits   AtBat  0.972
##  5 CRBI   CHits  0.968
##  6 CRBI   CAtBat 0.968
##  7 CRBI   CRuns  0.962
##  8 CWalks CRuns  0.950
##  9 CWalks CAtBat 0.944
## 10 CWalks CHits  0.940
## # … with 32 more rows

Dilakukan pemeriksaan terhadap korelasi antarpeubah. Nilai korelasi di atas difilter sehingga hanya mencetak output berupa peubah dengan korelasi yang kuat (lebih besar dari 0.6).

cor_df %>% filter(abs(corr)<=0.6)  

## # A tibble: 94 × 3
##    Var1   Var2    corr
##    <chr>  <chr>  <dbl>
##  1 HmRun  AtBat 0.534 
##  2 HmRun  Hits  0.529 
##  3 Walks  HmRun 0.448 
##  4 Years  AtBat 0.0526
##  5 Years  Hits  0.0814
##  6 Years  HmRun 0.170 
##  7 Years  Runs  0.0427
##  8 Years  RBI   0.166 
##  9 Years  Walks 0.154 
## 10 CAtBat AtBat 0.320 
## # … with 84 more rows

Nilai korelasi di atas difilter sehingga hanya mencetak output berupa peubah dengan korelasi yang rendah (lebih kecil dari 0.6).

1. Regresi Linear

regresi1 <- lm(formula = Salary~.,data = df.hitters1)
broom::tidy(regresi1) %>% mutate(across(where(is.numeric),                                ~ format(round(.x,3),                                             big.mark=",",                                             scientific=F)))

## # A tibble: 20 × 5
##    term        estimate std.error statistic  p.value
##    <chr>          <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)  163.      90.8        1.80  0.0736  
##  2 AtBat         -1.98     0.634     -3.12  0.00201 
##  3 Hits           7.50     2.38       3.15  0.00181 
##  4 HmRun          4.33     6.20       0.698 0.486   
##  5 Runs          -2.38     2.98      -0.797 0.426   
##  6 RBI           -1.04     2.60      -0.402 0.688   
##  7 Walks          6.23     1.83       3.41  0.000766
##  8 Years         -3.49    12.4       -0.281 0.779   
##  9 CAtBat        -0.171    0.135     -1.27  0.206   
## 10 CHits          0.134    0.675      0.199 0.843   
## 11 CHmRun        -0.173    1.62      -0.107 0.915   
## 12 CRuns          1.45     0.750      1.94  0.0538  
## 13 CRBI           0.808    0.693      1.17  0.245   
## 14 CWalks        -0.812    0.328     -2.47  0.0141  
## 15 LeagueN       62.6     79.3        0.790 0.430   
## 16 DivisionW   -117.      40.4       -2.89  0.00414 
## 17 PutOuts        0.282    0.0774     3.64  0.000333
## 18 Assists        0.371    0.221      1.68  0.0947  
## 19 Errors        -3.36     4.39      -0.765 0.445   
## 20 NewLeagueN   -24.8     79.0       -0.313 0.754

Output di atas menghasilkan koefisien regresi dari masing-masing peubah yang selanjutnya digunakan untuk membuat persamaan model.

2. Forward Selection

best_forward <- regsubsets(Salary~.,data = df.hitters1, method = "forward", nvmax = 19)

summ_forward <- summary(best_forward)
result_gof <- data.frame(
  Adj.R2 = c(which.max(summ_forward$adjr2),max(summ_forward$adjr2)),
  CP = c(which.min(summ_forward$cp),min(summ_forward$cp)),
  BIC = c(which.min(summ_forward$bic),min(summ_forward$bic))
)
result_gof

##       Adj.R2        CP       BIC
## 1 11.0000000 10.000000    6.0000
## 2  0.5225706  5.009317 -147.9169

coef_forward <- coef(best_forward, result_gof$BIC[1])
enframe(coef_forward) %>% mutate(across(where(is.numeric),                                ~ format(round(.x,3),                                             big.mark=",",                                             scientific=F)))

## # A tibble: 7 × 2
##   name           value
##   <chr>          <dbl>
## 1 (Intercept)   91.5  
## 2 AtBat         -1.87 
## 3 Hits           7.60 
## 4 Walks          3.70 
## 5 CRBI           0.643
## 6 DivisionW   -123.   
## 7 PutOuts        0.264

plot_forward <- data.frame(subset=seq_along(summ_forward$bic),
                           BIC = summ_forward$bic)
ggplot(plot_forward, aes(subset, BIC)) +
  geom_line(size = 1, color = "darkgrey") +
  geom_point(color = "steelblue", size = 2)+
  geom_vline(aes(xintercept = result_gof$BIC[1]),
             color = "steelblue", size = 1)+
  theme_bw()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Berdasarkan hasil seleksi peubah dengan metode forward selection, subset terbaik dengan nilai BIC terendah adalah -147.9169, yang terjadi saat subset = 6. Peubah penjelas yang masuk ke dalam model adalah AtBat, Hits, Walks, CRBI, DivisionW, dan PutOuts.

3. Backward Selection

best_backward <- regsubsets(Salary~., data = df.hitters1, method = "backward", nvmax = 19)

summ_backward <- summary(best_backward)
result_gof2 <- data.frame(
  Adj.R2 = c(which.max(summ_backward$adjr2),max(summ_backward$adjr2)),
  CP = c(which.min(summ_backward$cp),min(summ_backward$cp)),
  BIC = c(which.min(summ_backward$bic),min(summ_backward$bic))
)
result_gof2

##       Adj.R2        CP      BIC
## 1 11.0000000 10.000000    8.000
## 2  0.5225706  5.009317 -147.382

coef_backward <-  coef(best_backward,result_gof2$BIC[1])
enframe(coef_backward) %>% mutate(across(where(is.numeric),                                ~ format(round(.x,3),                                             big.mark=",",                                             scientific=F)))

## # A tibble: 9 × 2
##   name           value
##   <chr>          <dbl>
## 1 (Intercept)  117.   
## 2 AtBat         -2.03 
## 3 Hits           6.85 
## 4 Walks          6.44 
## 5 CRuns          0.705
## 6 CRBI           0.527
## 7 CWalks        -0.807
## 8 DivisionW   -124.   
## 9 PutOuts        0.275

plot_backward <- data.frame(subset = seq_along(summ_backward$bic),
                           BIC = summ_backward$bic)
ggplot(plot_backward, aes(subset, BIC))+
  geom_line(size=1, color = "darkgrey")+
  geom_point(color="steelblue", size=2)+
  geom_vline(aes(xintercept = result_gof$BIC[1]),
             color = "steelblue", size = 1.2)+
  theme_bw()

Berdasarkan hasil seleksi peubah dengan metode forward selection, subset terbaik dengan nilai BIC terendah adalah -147.382, yang terjadi saat subset = 6. Peubah penjelas yang masuk ke dalam model adalah AtBat, Hits, Walks, CRuns, CRBI, CWalks, DivisionW, dan PutOuts.