Optimisasi Seleksi Variabel untuk Indeks Ketahanan Pangan di Pulau Jawa: Pendekatan Lasso dan Quadratic Programming
Kelompok 7
2023-05-23
library(quadprog)
library(readxl)
library(ggplot2)
library(plotly)data_os<- read_excel("C:/Users/Dicky Girsang/Desktop/Semester 6/Optimisasi Statistika/Optimasi_K7/data_optimisasistk.xlsx")Berikut juga saya berikan link untuk download data di atas: https://github.com/lihardodicky/lasso-qp-fsi
A. Pre-processing Data
1. Normalisasi Data
databaku = scale(data_os)
X <- as.matrix(databaku[,2:8])
Y <- as.matrix(databaku[,1])2. Splitting Data
set.seed(12345)
indices <- sample(nrow(data_os))
train_prop <- 0.7
train_rows <- round(train_prop * nrow(data_os))
train_data <- data_os[indices[1:train_rows], ]
test_data <- data_os[indices[(train_rows + 1):nrow(data_os)], ]X <- as.matrix(train_data[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")])
Y <- train_data$Y
X_test <- as.matrix(test_data[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")])
Y_test <- test_data$YB. Seleksi Lasso dengan QP
Inisialisasi variabel n dengan jumlah baris matriks X, p dengan jumlah kolom matriks X, dan rmse_list sebagai vektor kosong yang akan digunakan untuk menyimpan RMSE untuk setiap nilai t.
Dibentuk model dengan koefisien regresi dari model OLS (ordinary least squares) dengan menggunakan formula Y ~ X - 1 (tanpa intercept). Tahap selanjutnya adalah mencari solusi beta dari masalah optimisasi, dengan matriks D, vektor d, matriks A, dan vektor b0 sebagai argumen.
calculate_RMSE <- function(X, Y, t_values){
n <- dim(X)[1]
p <- dim(X)[2]
rmse_list <- vector()
for (t in t_values) {
ols.beta <- lm(Y ~ X - 1)
t.max <- sum(abs(ols.beta$coef))
D <- 2 * t(X) %*% X
d <- 2 * t(Y) %*% X
A <- as.matrix(expand.grid(rep(list(c(-1,1)), p)))
b0 <- rep(-t, 2^p)
beta.lasso <- solve.QP(Dmat = D, dvec = d, Amat = t(A), bvec = b0)
predicted <- X_test %*% beta.lasso$solution
rmse <- sqrt(sum((Y_test - predicted)^2) / n)
rmse_list <- c(rmse_list, rmse)
}
results <- data.frame(t_values, RMSE = rmse_list)
best_model <- results[which.min(results$RMSE), ]
return(list(results = results, best_model = best_model))
}- Dengan parameter tuning (t) atau keketatan yang berbeda, dilakukan loop function untuk mendapatkan model terbaik. Semakin kecil nilai RMSE, maka model tersebut akan lebih baik dari model lainnya.
# Set the range of t values
t_values <- seq(0.05, 0.99, by = 0.01)result <- calculate_RMSE(X, Y, t_values)
print(result$results) # Menampilkan hasil RMSE untuk setiap nilai t## t_values RMSE
## 1 0.05 26.507841
## 2 0.06 26.164693
## 3 0.07 25.821647
## 4 0.08 25.478716
## 5 0.09 25.135897
## 6 0.10 24.793248
## 7 0.11 24.450664
## 8 0.12 24.108210
## 9 0.13 23.765891
## 10 0.14 23.423711
## 11 0.15 23.081679
## 12 0.16 22.739800
## 13 0.17 22.398079
## 14 0.18 22.056528
## 15 0.19 21.715155
## 16 0.20 21.373963
## 17 0.21 21.032965
## 18 0.22 20.692170
## 19 0.23 20.351587
## 20 0.24 20.011226
## 21 0.25 19.671104
## 22 0.26 19.331228
## 23 0.27 18.991613
## 24 0.28 18.652273
## 25 0.29 18.313224
## 26 0.30 17.974481
## 27 0.31 17.636064
## 28 0.32 17.297990
## 29 0.33 16.960280
## 30 0.34 16.622957
## 31 0.35 16.286045
## 32 0.36 15.949568
## 33 0.37 15.613557
## 34 0.38 15.278041
## 35 0.39 14.943054
## 36 0.40 14.608631
## 37 0.41 14.274815
## 38 0.42 13.941645
## 39 0.43 13.609174
## 40 0.44 13.277449
## 41 0.45 12.946530
## 42 0.46 12.616481
## 43 0.47 12.287370
## 44 0.48 11.959277
## 45 0.49 11.632285
## 46 0.50 11.306492
## 47 0.51 10.982005
## 48 0.52 10.658942
## 49 0.53 10.337437
## 50 0.54 10.017639
## 51 0.55 9.699718
## 52 0.56 9.383865
## 53 0.57 9.070295
## 54 0.58 8.759285
## 55 0.59 8.451049
## 56 0.60 8.145939
## 57 0.61 7.844321
## 58 0.62 7.546613
## 59 0.63 7.253296
## 60 0.64 6.964927
## 61 0.65 6.682144
## 62 0.66 6.405688
## 63 0.67 6.136227
## 64 0.68 5.875068
## 65 0.69 5.717320
## 66 0.70 5.578214
## 67 0.71 5.440025
## 68 0.72 5.302069
## 69 0.73 5.165832
## 70 0.74 5.030750
## 71 0.75 4.896919
## 72 0.76 4.764443
## 73 0.77 4.633439
## 74 0.78 4.504035
## 75 0.79 4.376373
## 76 0.80 4.251017
## 77 0.81 4.127360
## 78 0.82 4.005969
## 79 0.83 3.887055
## 80 0.84 3.770853
## 81 0.85 3.657621
## 82 0.86 3.547645
## 83 0.87 3.441236
## 84 0.88 3.338734
## 85 0.89 3.240258
## 86 0.90 3.144288
## 87 0.91 3.056293
## 88 0.92 2.973843
## 89 0.93 2.897411
## 90 0.94 2.827484
## 91 0.95 2.764557
## 92 0.96 2.709117
## 93 0.97 2.661632
## 94 0.98 2.622535
## 95 0.99 2.592205print(result$best_model) # Menampilkan model dengan RMSE terendah dan nilai t terkait## t_values RMSE
## 95 0.99 2.592205Dari rata-rata 95 parameter tuning (t), diperoleh besar RMSE model seperti dataframe di atas. Diketahui bahwa parameter t yang menghasilkan model terbaik adalah sebesar 0.99 atau mendekati 1. Selanjutnya, dilakukan pemodelan terhadap keseluruhan data dengan t = 0.99.
# Create a custom plot with advanced features and decorations
ggplot(result$results, aes(x = t_values, y = RMSE)) +
geom_line(color = "blue") +
geom_point(color = "darkred", size = 1) +
geom_vline(xintercept = result$best_model$t_values, linetype = "dashed", color = "darkred") +
labs(x = "t parameter", y = "RMSE") +
ggtitle("Relationship between t parameter and RMSE") +
theme_minimal() +
theme(plot.caption = element_text(color = "gray", size = 10, hjust = 1),
plot.margin = margin(20, 20, 50, 20)) +
labs(caption = "t is directly proportional to the RMSE value consistently")
C. Model Terbaik
# gunakan data keseluruhan
X <- as.matrix(data_os[, c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")])
Y <- data_os$Y# parameter t_best = 0.99
ols.beta = lm(Y ~ X -1)
n <- dim(X)[1]
p <- dim(X)[2]
D <- 2 * t(X) %*% X
d <- 2 * t(Y) %*% X
A <- as.matrix(expand.grid(rep(list(c(-1,1)),p)))
t = 0.99
b0 <- rep(-t,2^p)
beta.lasso<-solve.QP(Dmat=D,dvec=d,Amat=t(A),bvec=b0)
t(beta.lasso$solution)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.7891407 -0.003047194 0.1976357 7.414401e-17 2.081668e-17 8.404909e-06
## [,7] [,8]
## [1,] -0.0001623897 5.567142e-06Dari output di atas, diperoleh kesimpulan bahwa model regresi Lasso terbaik dapat dituliskan sebagai berikut:
\[ IKP = 0.7891X1 - 0.0031X2 + 0.1976X3 - 0.0002X7 \]
Interpretasi dari model di atas, dapat dituliskan sebagai
berikut:
1. Jika faktor lain dianggap konstan, ketika Angka Harapan
Hidup (AHH) meningkat 1%, maka persentase Indeks Ketahanan Pangan (IKP)
wilayah di Pulau Jawa akan meningkat sebesar 0.0.7891. Semakin tinggi
AHH suatu wilayah, maka skor komposit kondisi ketahanan pangan di
wilayah tersebut juga semakin tinggi.
2. Jika faktor lain dalam
posisi konstan, ketika Jumlah Penduduk Miskin (JPM) meningkat 1%, maka
persentase IKP wilayah di Pulau Jawa akan menurun sebesar 0.0031.
Artinya semakin banyak penduduk miskin di suatu wilayah, maka IKP
masyarakat akan menurun dalam skala yang tidak terlalu signifikan
(<0.05)
3. Jika faktor lain dianggap konstan, ketika Persentase
Rumah Tangga dengan Akses Air Minum Layak (AML) mengalami peningkatan
sebesar 1%, maka IKP masyarakat juga akan meningkat sebesar 0.1976
4. Jika faktor lain dianggap konstan, ketika Jumlah Industri di suatu
wilayah meningkat sebesar 1%, maka IKP masyarakat akan menurun secara
tidak signifikan sebesar 0.0002
D. Kesimpulan
Model regresi Least Absolute Shrinkage and Selection Operator (LASSO) menghasilkan 4 peubah signifikan, yaitu Angka Harapan Hidup (AHH), Jumlah Penduduk Miskin (JPM), Rumah Tangga dengan Akses Air Minum Layak (AML), serta Jumlah Industri di wilayah terkait.