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 29.13 100.000
## 2   5 33.61  79.430
## 3   5 38.23  63.100
## 4   5 42.83  50.120
## 5   5 47.27  39.810
## 6   5 51.39  31.620
## 7   5 55.08  25.120
## 8   5 58.26  19.950
## 9   5 60.91  15.850
## 10  5 63.07  12.590
## 11  5 64.77  10.000
## 12  5 66.09   7.943
## 13  5 67.11   6.310
## 14  5 67.88   5.012
## 15  5 68.47   3.981
## 16  5 68.92   3.162
## 17  5 69.27   2.512
## 18  5 69.54   1.995
## 19  5 69.76   1.585
## 20  5 69.94   1.259
## 21  5 70.09   1.000
## 22  5 70.21   0.794
## 23  5 70.31   0.631
## 24  5 70.39   0.501
## 25  5 70.46   0.398
## 26  5 70.51   0.316
## 27  5 70.55   0.251
## 28  5 70.58   0.200
## 29  5 70.60   0.158
## 30  5 70.62   0.126
## 31  5 70.63   0.100
## 32  5 70.64   0.079
## 33  5 70.65   0.063
## 34  5 70.65   0.050
## 35  5 70.65   0.040
## 36  5 70.66   0.032
## 37  5 70.66   0.025
## 38  5 70.66   0.020
## 39  5 70.66   0.016
## 40  5 70.66   0.013
## 41  5 70.66   0.010
## 42  5 70.66   0.008
## 43  5 70.66   0.006
## 44  5 70.66   0.005
## 45  5 70.66   0.004
## 46  5 70.66   0.003
## 47  5 70.66   0.003
## 48  5 70.66   0.002
## 49  5 70.66   0.002
## 50  5 70.66   0.001
## 51  5 70.66   0.001

Pada cara kedua ini 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.

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.3981072

lambda minimum bernilai 0.2511886 yang mendekati output ridge_reg dengan 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.050319569
## cylinders    -0.509352541
## displacement -0.008765559
## horsepower   -0.043381131
## weight       -0.003957876
## acceleration -0.037746155

# 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.356121 0.7045943

# 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 3.945211 0.7022579

Diperoleh hasil bahwa nilai akurasi (RMSE) dan R-square dari data training dan testing yang terlihat keduanya memiliki nilai RMSE tidak berbeda jauh

Pemilihan Banyaknya Cluster

Dalam hal ini, digunakan observasi terhadap dendogram serta melihat koefisien silhouette dan koefisien WSS.

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).

Hierarchical Clustering dengan Metode Complete Linkage

Diterapkan hierarchical clustering dengan menggunakan metode yang terlihat baik dari hasil dendogramnya, yaitu complete linkage dengan banyak cluster adalah 2 (k optimum)

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

Diketahui bahwa cluster 1 berisi kendaraan dengan spesifikasi yang lebih tinggi dibandingkan dengan cluster 2. Pernyataan 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 terlampau 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)

Berdasarkan plot di atas, dapat dilihat secara visual lebih banyak kendaraan yang berada di cluster 2 dibandingkan pada cluster 1. Kedua komponen utama dapat menjelaskan sekitar 91.9% keragaman pada data.

Pemilihan Banyaknya Cluster

Dalam hal ini penentuan banyaknya cluster dengan melihat koefisien silhouette dan koefisien 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, digunakan fungsi aggregate untuk menampilkan rata-rata peubah asli.

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

Tabel di atas adalah rata-rata setiap peubah pada masing-masing cluster yang dibuat. 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, misalnya blok silinder, kekuatan mesin, dan berat. Untuk peubah mpg, cluster 2 lebih rendah daripada cluster 1 dan 3. Hal ini menunjukkan bahwa dengan banyak 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 membandingkan kedua metode yang digunakan, terjadi perbedaan yang signifikan terhadap hasil penggerombolan. Jumlah cluster yang lebih baik dalam mengelompokkan data adalah tiga (K-Means Clustering), contohnya 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. Hal ini sesuai dengan prinsip penggerombolan, di mana penggerombolan akan lebih baik jika antar gerombol heterogen. Pada kedua metode, tidak terjadi overlap sehingga penggerombolan terbilang cukup baik.

Menangani Missing Values

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

## # A tibble: 263 x 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 atas. 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, seluruh peubah kategorik memiliki lebih dari satu kategori sehingga tidak perlu penghapusan peubah kategorik. Selain itu, terdapat perubahan pada informasi peubah pada 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="maroon"))

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

plot_bar(df.hitters1)

Diagram batang di atas merupakan besaran dari masing-masing peubah kategorik, yaitu league, division, dan new league. Frekuensi kategori dari masing-masing peubah tidak berbeda signifikan atau hampir balance.

Periksa korelasi peubah yang digunakan

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 x 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 lebih besar dari 0.6.

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

## # A tibble: 94 x 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 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 x 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 x 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.
## i 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 x 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.