Support Vector Regression (SVR) Analysis - Boston Housing Dataset

Analisis Komprehensif Prediksi Harga Rumah dengan SVR

1. Pendahuluan dan Teori SVR

1.1 Konsep Dasar SVR

Support Vector Regression (SVR) adalah teknik machine learning yang menggunakan konsep Support Vector Machine (SVM) untuk masalah regresi. SVR berusaha menemukan fungsi yang dapat memprediksi nilai target dengan toleransi error tertentu (epsilon-tube).

Prinsip Kerja SVR:

- Epsilon-tube (ε-tube): Zona toleransi di sekitar fungsi regresi dimana kesalahan dianggap nol

- Support Vectors: Titik data yang berada di luar epsilon-tube atau tepat di batasnya

- Kernel Trick: Memungkinkan SVR menangani hubungan nonlinear melalui transformasi ke ruang dimensi tinggi

Keunggulan SVR:

- Robust terhadap outlier karena hanya support vectors yang mempengaruhi model

- Efektif untuk data dimensi tinggi

- Menghindari overfitting melalui regularisasi

- Memberikan solusi global optimum

2. Setup dan Load Dataset

# Load required libraries
library(MASS)        # untuk dataset Boston
library(e1071)       # untuk SVR
library(caret)       # untuk preprocessing dan evaluasi
library(ggplot2)     # untuk visualisasi
library(corrplot)    # untuk correlation plot
library(gridExtra)   # untuk multiple plots
library(dplyr)       # untuk data manipulation
library(knitr)       # untuk tabel
library(kableExtra)  # untuk styling tabel
library(reshape2)    # untuk data reshaping

# Set seed untuk reproducibility
set.seed(123)

# Load Boston Housing dataset
data(Boston)
df <- Boston

# Tampilkan struktur data
str(df)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

# Informasi dataset
cat("Dimensi dataset:", dim(df)[1], "baris dan", dim(df)[2], "kolom\n")

## Dimensi dataset: 506 baris dan 14 kolom

cat("Variabel target: medv (median value of homes)\n\n")

## Variabel target: medv (median value of homes)

# Tampilkan beberapa baris pertama
kable(head(df), caption = "Preview Dataset Boston Housing") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Preview Dataset Boston Housing

crim	zn	indus	chas	nox	rm	age	dis	rad	tax	ptratio	black	lstat	medv
0.00632	18	2.31	0	0.538	6.575	65.2	4.0900	1	296	15.3	396.90	4.98	24.0
0.02731	0	7.07	0	0.469	6.421	78.9	4.9671	2	242	17.8	396.90	9.14	21.6
0.02729	0	7.07	0	0.469	7.185	61.1	4.9671	2	242	17.8	392.83	4.03	34.7
0.03237	0	2.18	0	0.458	6.998	45.8	6.0622	3	222	18.7	394.63	2.94	33.4
0.06905	0	2.18	0	0.458	7.147	54.2	6.0622	3	222	18.7	396.90	5.33	36.2
0.02985	0	2.18	0	0.458	6.430	58.7	6.0622	3	222	18.7	394.12	5.21	28.7

# Summary statistik
summary(df)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Interpretasi:

Dataset Boston Housing memiliki 506 observasi dengan 14 variabel, dimana medv adalah variabel target yang mewakili median nilai rumah dalam ribuan dollar. Variabel prediktor mencakup karakteristik sosial-ekonomi dan lingkungan seperti tingkat kriminalitas (crim), proporsi lahan industri (indus), kualitas udara (nox), dll.

3. Exploratory Data Analysis (EDA)

# Cek missing values
cat("Missing values per kolom:\n")

## Missing values per kolom:

missing_values <- sapply(df, function(x) sum(is.na(x)))
print(missing_values)

##    crim      zn   indus    chas     nox      rm     age     dis     rad     tax 
##       0       0       0       0       0       0       0       0       0       0 
## ptratio   black   lstat    medv 
##       0       0       0       0

cat("\nTotal missing values:", sum(missing_values), "\n")

## 
## Total missing values: 0

# Visualisasi distribusi variabel target
p1 <- ggplot(df, aes(x = medv)) +
  geom_histogram(bins = 30, fill = "skyblue", alpha = 0.7, color = "black") +
  labs(title = "Distribusi Harga Rumah (medv)",
       x = "Median Value (dalam $1000s)",
       y = "Frekuensi") +
  theme_minimal()

p2 <- ggplot(df, aes(y = medv)) +
  geom_boxplot(fill = "lightcoral", alpha = 0.7) +
  labs(title = "Boxplot Harga Rumah",
       y = "Median Value (dalam $1000s)") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

# Correlation matrix
cor_matrix <- cor(df)
corrplot(cor_matrix, method = "color", type = "upper", 
         order = "hclust", tl.cex = 0.8, tl.col = "black",
         title = "Correlation Matrix - Boston Housing Dataset")

# Korelasi tertinggi dengan target
cor_with_target <- cor(df[,-14], df$medv)
cor_df <- data.frame(
  Variable = rownames(cor_with_target),
  Correlation = cor_with_target[,1]
)
cor_df <- cor_df %>% arrange(desc(abs(Correlation)))

kable(cor_df, caption = "Korelasi Variabel dengan MEDV", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Korelasi Variabel dengan MEDV

	Variable	Correlation
lstat	lstat	-0.738
rm	rm	0.695
ptratio	ptratio	-0.508
indus	indus	-0.484
tax	tax	-0.469
nox	nox	-0.427
crim	crim	-0.388
rad	rad	-0.382
age	age	-0.377
zn	zn	0.360
black	black	0.333
dis	dis	0.250
chas	chas	0.175

# Scatter plot beberapa variabel penting vs medv
p1 <- ggplot(df, aes(x = lstat, y = medv)) +
  geom_point(alpha = 0.6, color = "blue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "LSTAT vs MEDV", x = "% Lower Status Population", y = "Median Value") +
  theme_bw()

p2 <- ggplot(df, aes(x = rm, y = medv)) +
  geom_point(alpha = 0.6, color = "green") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "RM vs MEDV", x = "Average Rooms per Dwelling", y = "Median Value") +
  theme_bw()

p3 <- ggplot(df, aes(x = crim, y = medv)) +
  geom_point(alpha = 0.6, color = "orange") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "CRIM vs MEDV", x = "Crime Rate", y = "Median Value") +
  theme_bw()

p4 <- ggplot(df, aes(x = dis, y = medv)) +
  geom_point(alpha = 0.6, color = "purple") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "DIS vs MEDV", x = "Distance to Employment Centers", y = "Median Value") +
  theme_bw()

grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)

Interpretasi:

- Dataset tidak memiliki missing values

- Distribusi harga rumah (medv) menunjukkan distribusi yang sedikit right-skewed dengan beberapa outliers

- Variabel lstat (% lower status population) memiliki korelasi negatif terkuat (-0.738) dengan harga rumah

- Variabel rm (average number of rooms) memiliki korelasi positif terkuat (0.696) dengan harga rumah

4. Preprocessing Data

# Normalisasi data (feature scaling)
# SVR sensitif terhadap skala data
preprocess_params <- preProcess(df[,-14], method = c("center", "scale"))
df_scaled <- predict(preprocess_params, df[,-14])
df_scaled$medv <- df$medv

# Tampilkan perbandingan sebelum dan sesudah scaling
before_after <- data.frame(
  Variable = names(df[1:6]),
  Mean_Before = sapply(df[1:6], mean),
  SD_Before = sapply(df[1:6], sd),
  Mean_After = sapply(df_scaled[1:6], mean),
  SD_After = sapply(df_scaled[1:6], sd)
)

kable(before_after, caption = "Perbandingan Sebelum dan Sesudah Scaling", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Perbandingan Sebelum dan Sesudah Scaling

	Variable	Mean_Before	SD_Before	Mean_After	SD_After
crim	crim	3.614	8.602	0	1
zn	zn	11.364	23.322	0	1
indus	indus	11.137	6.860	0	1
chas	chas	0.069	0.254	0	1
nox	nox	0.555	0.116	0	1
rm	rm	6.285	0.703	0	1

# Split data menjadi training dan testing (75:25)
trainIndex <- createDataPartition(df_scaled$medv, p = 0.75, list = FALSE)
train_data <- df_scaled[trainIndex, ]
test_data <- df_scaled[-trainIndex, ]

cat("Training set:", nrow(train_data), "samples\n")

## Training set: 381 samples

cat("Testing set:", nrow(test_data), "samples\n")

## Testing set: 125 samples

Interpretasi:

Data dinormalisasi menggunakan center dan scale untuk memastikan semua fitur memiliki mean=0 dan standard deviation=1. Hal ini penting untuk SVR karena algoritma ini sensitif terhadap skala data. Data dibagi menjadi 75% training dan 25% testing dengan stratified sampling.

5. Model Building

5.1 Fungsi Evaluasi

# Fungsi untuk menghitung metrik evaluasi
calculate_metrics <- function(actual, predicted) {
  rmse <- sqrt(mean((actual - predicted)^2))
  mae <- mean(abs(actual - predicted))
  r2 <- 1 - sum((actual - predicted)^2) / sum((actual - mean(actual))^2)
  
  return(c(RMSE = rmse, MAE = mae, R2 = r2))
}

5.2 SVR Linear

# SVR Linear
svr_linear <- svm(medv ~ ., data = train_data, 
                  type = "eps-regression", 
                  kernel = "linear",
                  cost = 1, 
                  epsilon = 0.1)

# Prediksi
pred_linear_train <- predict(svr_linear, train_data)
pred_linear_test <- predict(svr_linear, test_data)

# Evaluasi
metrics_linear_train <- calculate_metrics(train_data$medv, pred_linear_train)
metrics_linear_test <- calculate_metrics(test_data$medv, pred_linear_test)

cat("=== SVR LINEAR RESULTS ===\n")

## === SVR LINEAR RESULTS ===

cat("Training - RMSE:", round(metrics_linear_train[1], 3), 
    "MAE:", round(metrics_linear_train[2], 3), 
    "R²:", round(metrics_linear_train[3], 3), "\n")

## Training - RMSE: 4.238 MAE: 2.789 R²: 0.754

cat("Testing  - RMSE:", round(metrics_linear_test[1], 3), 
    "MAE:", round(metrics_linear_test[2], 3), 
    "R²:", round(metrics_linear_test[3], 3), "\n")

## Testing  - RMSE: 6.867 MAE: 4.071 R²: 0.598

cat("Number of Support Vectors:", svr_linear$tot.nSV, "\n\n")

## Number of Support Vectors: 294

print(svr_linear)

## 
## Call:
## svm(formula = medv ~ ., data = train_data, type = "eps-regression", 
##     kernel = "linear", cost = 1, epsilon = 0.1)
## 
## 
## Parameters:
##    SVM-Type:  eps-regression 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.07692308 
##     epsilon:  0.1 
## 
## 
## Number of Support Vectors:  294

5.3 SVR Nonlinear (RBF Kernel)

# SVR dengan RBF kernel
svr_rbf <- svm(medv ~ ., data = train_data, 
               type = "eps-regression", 
               kernel = "radial",
               cost = 10, 
               gamma = 0.1,
               epsilon = 0.1)

# Prediksi
pred_rbf_train <- predict(svr_rbf, train_data)
pred_rbf_test <- predict(svr_rbf, test_data)

# Evaluasi
metrics_rbf_train <- calculate_metrics(train_data$medv, pred_rbf_train)
metrics_rbf_test <- calculate_metrics(test_data$medv, pred_rbf_test)

cat("=== SVR RBF RESULTS ===\n")

## === SVR RBF RESULTS ===

cat("Training - RMSE:", round(metrics_rbf_train[1], 3), 
    "MAE:", round(metrics_rbf_train[2], 3), 
    "R²:", round(metrics_rbf_train[3], 3), "\n")

## Training - RMSE: 1.329 MAE: 1.001 R²: 0.976

cat("Testing  - RMSE:", round(metrics_rbf_test[1], 3), 
    "MAE:", round(metrics_rbf_test[2], 3), 
    "R²:", round(metrics_rbf_test[3], 3), "\n")

## Testing  - RMSE: 4.34 MAE: 2.388 R²: 0.84

cat("Number of Support Vectors:", svr_rbf$tot.nSV, "\n\n")

## Number of Support Vectors: 267

print(svr_rbf)

## 
## Call:
## svm(formula = medv ~ ., data = train_data, type = "eps-regression", 
##     kernel = "radial", cost = 10, gamma = 0.1, epsilon = 0.1)
## 
## 
## Parameters:
##    SVM-Type:  eps-regression 
##  SVM-Kernel:  radial 
##        cost:  10 
##       gamma:  0.1 
##     epsilon:  0.1 
## 
## 
## Number of Support Vectors:  267

5.4 Model OLS (Ordinary Least Squares)

# Ordinary Least Squares Regression untuk perbandingan
ols_model <- lm(medv ~ ., data = train_data)
pred_ols_train <- predict(ols_model, train_data)
pred_ols_test <- predict(ols_model, test_data)

# Evaluasi OLS
metrics_ols_train <- calculate_metrics(train_data$medv, pred_ols_train)
metrics_ols_test <- calculate_metrics(test_data$medv, pred_ols_test)

cat("=== OLS REGRESSION RESULTS ===\n")

## === OLS REGRESSION RESULTS ===

cat("Training - RMSE:", round(metrics_ols_train[1], 3), 
    "MAE:", round(metrics_ols_train[2], 3), 
    "R²:", round(metrics_ols_train[3], 3), "\n")

## Training - RMSE: 4.085 MAE: 2.906 R²: 0.772

cat("Testing  - RMSE:", round(metrics_ols_test[1], 3), 
    "MAE:", round(metrics_ols_test[2], 3), 
    "R²:", round(metrics_ols_test[3], 3), "\n\n")

## Testing  - RMSE: 6.421 MAE: 4.029 R²: 0.649

summary(ols_model)

## 
## Call:
## lm(formula = medv ~ ., data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.692  -2.359  -0.523   1.624  21.168 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 22.19880    0.21395 103.756  < 2e-16 ***
## crim        -0.68720    0.29453  -2.333 0.020177 *  
## zn           0.64196    0.31927   2.011 0.045088 *  
## indus       -0.06547    0.42843  -0.153 0.878633    
## chas         0.37752    0.23566   1.602 0.110028    
## nox         -1.86040    0.45475  -4.091 5.28e-05 ***
## rm           2.82595    0.29767   9.493  < 2e-16 ***
## age          0.01096    0.38360   0.029 0.977220    
## dis         -2.35504    0.40910  -5.757 1.81e-08 ***
## rad          2.27022    0.58923   3.853 0.000138 ***
## tax         -2.22790    0.63866  -3.488 0.000545 ***
## ptratio     -2.02792    0.28297  -7.167 4.25e-12 ***
## black        0.74631    0.23641   3.157 0.001727 ** 
## lstat       -3.04032    0.40049  -7.592 2.64e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.162 on 367 degrees of freedom
## Multiple R-squared:  0.7718, Adjusted R-squared:  0.7638 
## F-statistic:  95.5 on 13 and 367 DF,  p-value: < 2.2e-16

# Tabel perbandingan performa model
train_comparison <- data.frame(
  Model = c("SVR Linear", "SVR RBF", "OLS Regression"),
  RMSE = round(c(metrics_linear_train[1], metrics_rbf_train[1], metrics_ols_train[1]), 3),
  MAE = round(c(metrics_linear_train[2], metrics_rbf_train[2], metrics_ols_train[2]), 3),
  R_squared = round(c(metrics_linear_train[3], metrics_rbf_train[3], metrics_ols_train[3]), 3),
  Support_Vectors = c(svr_linear$tot.nSV, svr_rbf$tot.nSV, "-")
)

test_comparison <- data.frame(
  Model = c("SVR Linear", "SVR RBF", "OLS Regression"),
  RMSE = round(c(metrics_linear_test[1], metrics_rbf_test[1], metrics_ols_test[1]), 3),
  MAE = round(c(metrics_linear_test[2], metrics_rbf_test[2], metrics_ols_test[2]), 3),
  R_squared = round(c(metrics_linear_test[3], metrics_rbf_test[3], metrics_ols_test[3]), 3),
  Support_Vectors = c(svr_linear$tot.nSV, svr_rbf$tot.nSV, "-")
)

kable(train_comparison, caption = "Perbandingan Performa Model - Data Training") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(which.min(train_comparison$RMSE), bold = T, color = "white", background = "#D7261E")

Perbandingan Performa Model - Data Training

Model	RMSE	MAE	R_squared	Support_Vectors
SVR Linear	4.238	2.789	0.754	294
SVR RBF	1.329	1.001	0.976	267
OLS Regression	4.085	2.906	0.772

kable(test_comparison, caption = "Perbandingan Performa Model - Data Testing") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(which.min(test_comparison$RMSE), bold = T, color = "white", background = "#D7261E")

Perbandingan Performa Model - Data Testing

Model	RMSE	MAE	R_squared	Support_Vectors
SVR Linear	6.867	4.071	0.598	294
SVR RBF	4.340	2.388	0.840	267
OLS Regression	6.421	4.029	0.649

Interpretasi:

Berdasarkan hasil perbandingan, SVR dengan kernel RBF umumnya menunjukkan performa superior untuk data dengan pola nonlinear karena kemampuannya menangkap kompleksitas yang tidak dapat ditangkap oleh model linear.

6. Visualisasi Epsilon-tube dan Support Vectors

6.1 Visualisasi 1D untuk Konsep Epsilon-tube

# Untuk visualisasi 1D, gunakan variabel dengan korelasi tertinggi (lstat)
train_1d <- data.frame(lstat = train_data$lstat, medv = train_data$medv)
test_1d <- data.frame(lstat = test_data$lstat, medv = test_data$medv)

# Model SVR 1D untuk visualisasi
svr_1d_linear <- svm(medv ~ lstat, data = train_1d, 
                     type = "eps-regression", 
                     kernel = "linear",
                     cost = 1, 
                     epsilon = 0.5)

svr_1d_rbf <- svm(medv ~ lstat, data = train_1d, 
                  type = "eps-regression", 
                  kernel = "radial",
                  cost = 10, 
                  gamma = 0.1,
                  epsilon = 0.5)

# Buat data untuk plotting
lstat_range <- seq(min(train_1d$lstat), max(train_1d$lstat), length.out = 100)
pred_1d_linear <- predict(svr_1d_linear, data.frame(lstat = lstat_range))
pred_1d_rbf <- predict(svr_1d_rbf, data.frame(lstat = lstat_range))

# Plot SVR Linear dengan epsilon-tube
p1 <- ggplot() +
  geom_point(data = train_1d, aes(x = lstat, y = medv), alpha = 0.6, color = "blue", size = 1.5) +
  geom_line(data = data.frame(lstat = lstat_range, pred = pred_1d_linear), 
            aes(x = lstat, y = pred), color = "red", size = 1.2) +
  geom_ribbon(data = data.frame(lstat = lstat_range, pred = pred_1d_linear), 
              aes(x = lstat, ymin = pred - 0.5, ymax = pred + 0.5), 
              alpha = 0.3, fill = "yellow") +
  labs(title = "SVR Linear dengan Epsilon-tube (ε = 0.5)",
       subtitle = "Garis merah: prediksi SVR, Area kuning: epsilon-tube, Titik biru: data training",
       x = "LSTAT (% lower status population) - Normalized", 
       y = "MEDV (Median home value in $1000s)") +
  theme_minimal()

# Plot SVR RBF dengan epsilon-tube
p2 <- ggplot() +
  geom_point(data = train_1d, aes(x = lstat, y = medv), alpha = 0.6, color = "blue", size = 1.5) +
  geom_line(data = data.frame(lstat = lstat_range, pred = pred_1d_rbf), 
            aes(x = lstat, y = pred), color = "red", size = 1.2) +
  geom_ribbon(data = data.frame(lstat = lstat_range, pred = pred_1d_rbf), 
              aes(x = lstat, ymin = pred - 0.5, ymax = pred + 0.5), 
              alpha = 0.3, fill = "yellow") +
  labs(title = "SVR RBF dengan Epsilon-tube (ε = 0.5)",
       subtitle = "Garis merah: prediksi SVR, Area kuning: epsilon-tube, Titik biru: data training",
       x = "LSTAT (% lower status population) - Normalized", 
       y = "MEDV (Median home value in $1000s)") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

6.2 Identifikasi Support Vectors

# Informasi support vectors untuk model multivariat
cat("=== SUPPORT VECTORS INFORMATION ===\n")

## === SUPPORT VECTORS INFORMATION ===

cat("SVR Linear - Total Support Vectors:", svr_linear$tot.nSV, "\n")

## SVR Linear - Total Support Vectors: 294

cat("SVR RBF - Total Support Vectors:", svr_rbf$tot.nSV, "\n")

## SVR RBF - Total Support Vectors: 267

cat("SVR 1D Linear - Total Support Vectors:", svr_1d_linear$tot.nSV, "\n")

## SVR 1D Linear - Total Support Vectors: 131

cat("SVR 1D RBF - Total Support Vectors:", svr_1d_rbf$tot.nSV, "\n")

## SVR 1D RBF - Total Support Vectors: 112

# Identifikasi support vectors dalam model 1D RBF
sv_indices <- svr_1d_rbf$index
support_vectors_1d <- train_1d[sv_indices, ]

cat("\nFirst 10 Support Vectors (1D RBF model):\n")

## 
## First 10 Support Vectors (1D RBF model):

print(head(support_vectors_1d, 10))

##         lstat medv
## 1  -1.0744990 24.0
## 4  -1.0254866 36.2
## 6   0.9097999 27.1
## 9   0.4280788 21.7
## 10 -0.5857761 19.9
## 11 -0.8504426 23.1
## 31 -1.0941039 26.6
## 40 -1.0324884 25.0
## 43 -0.9638712 24.7
## 45 -0.8112328 23.3

# Visualisasi support vectors
ggplot() +
  geom_point(data = train_1d, aes(x = lstat, y = medv), alpha = 0.4, color = "lightblue", size = 1) +
  geom_point(data = support_vectors_1d, aes(x = lstat, y = medv), color = "red", size = 2, shape = 1, stroke = 2) +
  geom_line(data = data.frame(lstat = lstat_range, pred = pred_1d_rbf), 
            aes(x = lstat, y = pred), color = "darkred", size = 1.2) +
  labs(title = "Support Vectors dalam SVR RBF",
       subtitle = "Lingkaran merah: Support Vectors, Titik biru: Data training lainnya",
       x = "LSTAT (% lower status population) - Normalized", 
       y = "MEDV (Median home value in $1000s)") +
  theme_minimal()

Interpretasi:

- Epsilon-tube (area kuning) menunjukkan zona toleransi error sebesar ε=0.5

- Data points yang berada di dalam epsilon-tube tidak berkontribusi pada fungsi SVR

- Support vectors (lingkaran merah) adalah data points yang berada tepat pada batas atau di luar epsilon-tube

- Semakin kecil ε, semakin ketat toleransi error dan semakin banyak support vectors - SVR RBF dapat menangkap pola nonlinear lebih baik dibandingkan SVR Linear

7. Visualisasi Model Performance

7.1 Actual vs Predicted

# Buat dataframe untuk plotting
plot_data_test <- data.frame(
  Actual = test_data$medv,
  SVR_Linear = pred_linear_test,
  SVR_RBF = pred_rbf_test,
  OLS = pred_ols_test
)

# Reshape data untuk ggplot
plot_data_melted <- melt(plot_data_test, id.vars = "Actual", 
                        variable.name = "Model", value.name = "Predicted")

# Plot Actual vs Predicted
ggplot(plot_data_melted, aes(x = Actual, y = Predicted, color = Model)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_abline(slope = 1, intercept = 0, color = "black", linetype = "dashed", size = 1) +
  facet_wrap(~Model, scales = "free") +
  labs(title = "Predicted vs Actual Values - Test Data Comparison",
       subtitle = "Garis putus-putus: Perfect prediction (y = x)",
       x = "Actual MEDV", 
       y = "Predicted MEDV") +
  theme_minimal() +
  theme(legend.position = "none")

7.2 Residual Analysis

# Hitung residual untuk semua model
residuals_linear <- test_data$medv - pred_linear_test
residuals_rbf <- test_data$medv - pred_rbf_test
residuals_ols <- test_data$medv - pred_ols_test

# Residual plots
residuals_data <- data.frame(
  Predicted_SVR_Linear = pred_linear_test,
  Residuals_SVR_Linear = residuals_linear,
  Predicted_SVR_RBF = pred_rbf_test,
  Residuals_SVR_RBF = residuals_rbf,
  Predicted_OLS = pred_ols_test,
  Residuals_OLS = residuals_ols
)

p1 <- ggplot(residuals_data, aes(x = Predicted_SVR_Linear, y = Residuals_SVR_Linear)) +
  geom_point(alpha = 0.6, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "SVR Linear - Residuals", x = "Predicted Values", y = "Residuals") +
  theme_bw()

p2 <- ggplot(residuals_data, aes(x = Predicted_SVR_RBF, y = Residuals_SVR_RBF)) +
  geom_point(alpha = 0.6, color = "green") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "SVR RBF - Residuals", x = "Predicted Values", y = "Residuals") +
  theme_bw()

p3 <- ggplot(residuals_data, aes(x = Predicted_OLS, y = Residuals_OLS)) +
  geom_point(alpha = 0.6, color = "purple") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "OLS - Residuals", x = "Predicted Values", y = "Residuals") +
  theme_bw()

# Q-Q plot untuk normalitas residual
p4 <- ggplot(data.frame(residuals = residuals_rbf), aes(sample = residuals)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Q-Q Plot SVR RBF Residuals") +
  theme_bw()

grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)

Interpretasi:

- Plot residual membantu mengidentifikasi pola error dan validasi asumsi model

- Residual yang tersebar acak di sekitar garis y=0 menunjukkan model yang baik

- Q-Q plot menunjukkan distribusi residual. Jika mengikuti garis diagonal, residual berdistribusi normal

- Predicted vs Actual plot dengan titik-titik yang dekat dengan garis y=x menunjukkan prediksi yang akurat

8. Eksplorasi Parameter Lengkap

8.1 Pengaruh Parameter Epsilon

# Test berbagai nilai epsilon
epsilon_values <- c(0.01, 0.1, 0.5, 1.0, 2.0)
epsilon_results <- data.frame()

for (eps in epsilon_values) {
  svr_temp <- svm(medv ~ ., data = train_data, 
                  type = "eps-regression", 
                  kernel = "radial", 
                  cost = 10, 
                  gamma = 0.1,
                  epsilon = eps)
  
  pred_temp <- predict(svr_temp, test_data)
  metrics_temp <- calculate_metrics(test_data$medv, pred_temp)
  
  epsilon_results <- rbind(epsilon_results, 
                          data.frame(Epsilon = eps, 
                                   RMSE = metrics_temp[1], 
                                   MAE = metrics_temp[2], 
                                   R2 = metrics_temp[3],
                                   Support_Vectors = svr_temp$tot.nSV))
}

kable(epsilon_results, digits = 4, caption = "Pengaruh Parameter Epsilon terhadap Performa Model") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Pengaruh Parameter Epsilon terhadap Performa Model

	Epsilon	RMSE	MAE	R2	Support_Vectors
RMSE	0.01	4.3893	2.4540	0.8359	368
RMSE1	0.10	4.3397	2.3883	0.8396	267
RMSE2	0.50	4.5458	2.8252	0.8240	63
RMSE3	1.00	6.2257	4.6695	0.6699	30
RMSE4	2.00	10.7297	9.6696	0.0195	17

# Plot pengaruh epsilon
p1 <- ggplot(epsilon_results, aes(x = Epsilon, y = RMSE)) +
  geom_line(color = "blue", size = 1) + geom_point(color = "blue", size = 2) +
  labs(title = "Pengaruh Epsilon terhadap RMSE", y = "RMSE") + theme_bw()

p2 <- ggplot(epsilon_results, aes(x = Epsilon, y = R2)) +
  geom_line(color = "red", size = 1) + geom_point(color = "red", size = 2) +
  labs(title = "Pengaruh Epsilon terhadap R²", y = "R²") + theme_bw()

p3 <- ggplot(epsilon_results, aes(x = Epsilon, y = Support_Vectors)) +
  geom_line(color = "green", size = 1) + geom_point(color = "green", size = 2) +
  labs(title = "Pengaruh Epsilon terhadap Support Vectors", y = "Number of Support Vectors") + theme_bw()

grid.arrange(p1, p2, p3, ncol = 3)

8.2 Pengaruh Parameter Cost (C)

# Test berbagai nilai cost
cost_values <- c(0.1, 1, 10, 100, 1000)
cost_results <- data.frame()

for (c_val in cost_values) {
  svr_temp <- svm(medv ~ ., data = train_data, 
                  type = "eps-regression", 
                  kernel = "radial", 
                  cost = c_val, 
                  gamma = 0.1,
                  epsilon = 0.1)
  
  pred_temp <- predict(svr_temp, test_data)
  metrics_temp <- calculate_metrics(test_data$medv, pred_temp)
  
  cost_results <- rbind(cost_results, 
                       data.frame(Cost = c_val, 
                                RMSE = metrics_temp[1], 
                                MAE = metrics_temp[2], 
                                R2 = metrics_temp[3],
                                Support_Vectors = svr_temp$tot.nSV))
}

kable(cost_results, digits = 4, caption = "Pengaruh Parameter Cost")

Pengaruh Parameter Cost

	Cost	RMSE	MAE	R2	Support_Vectors
RMSE	1e-01	8.0947	4.6471	0.4420	283
RMSE1	1e+00	5.8094	3.0252	0.7126	256
RMSE2	1e+01	4.3397	2.3883	0.8396	267
RMSE3	1e+02	4.3466	2.6857	0.8391	268
RMSE4	1e+03	4.2066	2.7079	0.8493	279

# Plot pengaruh cost dengan dual y-axis visualization
p1 <- ggplot(cost_results, aes(x = log10(Cost), y = RMSE)) +
  geom_line(color = "blue", size = 1.2) + 
  geom_point(color = "red", size = 3) +
  labs(title = "Pengaruh Cost terhadap RMSE", x = "log10(Cost)") + 
  theme_bw()

p2 <- ggplot(cost_results, aes(x = log10(Cost), y = R2)) +
  geom_line(color = "red", size = 1.2) + 
  geom_point(color = "red", size = 3) +
  labs(title = "Pengaruh Cost terhadap R²", x = "log10(Cost)") + 
  theme_bw()

p3 <- ggplot(cost_results, aes(x = log10(Cost), y = Support_Vectors)) +
  geom_line(color = "green", size = 1.2) + 
  geom_point(color = "red", size = 3) +
  labs(title = "Pengaruh Cost terhadap Support Vectors", x = "log10(Cost)") + 
  theme_bw()

grid.arrange(p1, p2, p3, ncol = 3)

print(cost_results)

##        Cost     RMSE      MAE        R2 Support_Vectors
## RMSE  1e-01 8.094723 4.647108 0.4419590             283
## RMSE1 1e+00 5.809419 3.025219 0.7125729             256
## RMSE2 1e+01 4.339737 2.388275 0.8396057             267
## RMSE3 1e+02 4.346568 2.685743 0.8391003             268
## RMSE4 1e+03 4.206558 2.707855 0.8492991             279

8.3 Pengaruh Parameter Gamma

# Test berbagai nilai gamma
gamma_values <- c(0.001, 0.01, 0.1, 1, 10)
gamma_results <- data.frame()

for (g_val in gamma_values) {
  svr_temp <- svm(medv ~ ., data = train_data, 
                  type = "eps-regression", 
                  kernel = "radial", 
                  cost = 10, 
                  gamma = g_val,
                  epsilon = 0.1)
  
  pred_temp <- predict(svr_temp, test_data)
  metrics_temp <- calculate_metrics(test_data$medv, pred_temp)
  
  gamma_results <- rbind(gamma_results, 
                        data.frame(Gamma = g_val, 
                                 RMSE = metrics_temp[1], 
                                 MAE = metrics_temp[2], 
                                 R2 = metrics_temp[3]))
}

kable(gamma_results, digits = 4, caption = "Pengaruh Parameter Gamma")

Pengaruh Parameter Gamma

	Gamma	RMSE	MAE	R2
RMSE	1e-03	6.9072	4.0394	0.5937
RMSE1	1e-02	5.8137	3.0105	0.7121
RMSE2	1e-01	4.3397	2.3883	0.8396
RMSE3	1e+00	7.1545	4.3952	0.5641
RMSE4	1e+01	10.3316	7.0158	0.0909

# Plot pengaruh gamma
p1 <- ggplot(gamma_results, aes(x = log10(Gamma), y = RMSE)) +
  geom_line(color = "blue", size = 1.2) + 
  geom_point(color = "blue", size = 3) +
  labs(title = "Pengaruh Gamma terhadap RMSE", x = "log10(Gamma)") + 
  theme_bw()

p2 <- ggplot(gamma_results, aes(x = log10(Gamma), y = R2)) +
  geom_line(color = "red", size = 1.2) + 
  geom_point(color = "red", size = 3) +
  labs(title = "Pengaruh Gamma terhadap R²", x = "log10(Gamma)") + 
  theme_bw()

grid.arrange(p1, p2, ncol = 2)

9. Grid Search untuk Parameter Optimal

# Load required libraries
library(e1071)

# Comprehensive Grid search untuk parameter optimal
cat("Melakukan Comprehensive Grid Search untuk optimasi parameter...\n")

## Melakukan Comprehensive Grid Search untuk optimasi parameter...

tune_result <- tune(svm, medv ~ ., data = train_data,
                    type = "eps-regression",
                    kernel = "radial",
                    ranges = list(cost = c(0.1, 1, 10, 100),
                                  gamma = c(0.01, 0.1, 1, 10),
                                  epsilon = c(0.01, 0.1, 0.2)),
                    tunecontrol = tune.control(cross = 5))

# Best parameters
cat("=== GRID SEARCH RESULTS ===\n")

## === GRID SEARCH RESULTS ===

cat("Best Parameters:\n")

## Best Parameters:

print(tune_result$best.parameters)

##    cost gamma epsilon
## 19   10  0.01     0.1

cat("Best Performance (MSE):", round(tune_result$best.performance, 4), "\n")

## Best Performance (MSE): 10.3878

# Get the best model
best_svr <- tune_result$best.model
pred_best <- predict(best_svr, test_data)
rmse_best <- sqrt(mean((test_data$medv - pred_best)^2))
mae_best <- mean(abs(test_data$medv - pred_best))
r2_best <- cor(test_data$medv, pred_best)^2

cat("=== OPTIMIZED SVR RESULTS ===\n")

## === OPTIMIZED SVR RESULTS ===

cat("RMSE:", round(rmse_best, 3), "\n")

## RMSE: 5.814

cat("MAE:", round(mae_best, 3), "\n")

## MAE: 3.011

cat("R²:", round(r2_best, 3), "\n")

## R²: 0.739

cat("Number of Support Vectors:", best_svr$tot.nSV, "\n")

## Number of Support Vectors: 267

# Store optimal metrics for final comparison
metrics_optimal <- c(rmse_best, mae_best, r2_best)

# Create baseline models for comparison
cat("\nCreating baseline models for comparison...\n")

## 
## Creating baseline models for comparison...

# SVR Linear
svr_linear <- svm(medv ~ ., data = train_data, 
                  type = "eps-regression", 
                  kernel = "linear")
pred_linear <- predict(svr_linear, test_data)
rmse_linear <- sqrt(mean((test_data$medv - pred_linear)^2))
mae_linear <- mean(abs(test_data$medv - pred_linear))
r2_linear <- cor(test_data$medv, pred_linear)^2

# SVR RBF (default parameters)
svr_rbf <- svm(medv ~ ., data = train_data, 
               type = "eps-regression", 
               kernel = "radial")
pred_rbf <- predict(svr_rbf, test_data)
rmse_rbf <- sqrt(mean((test_data$medv - pred_rbf)^2))
mae_rbf <- mean(abs(test_data$medv - pred_rbf))
r2_rbf <- cor(test_data$medv, pred_rbf)^2

# OLS Regression
ols_model <- lm(medv ~ ., data = train_data)
pred_ols <- predict(ols_model, test_data)
rmse_ols <- sqrt(mean((test_data$medv - pred_ols)^2))
mae_ols <- mean(abs(test_data$medv - pred_ols))
r2_ols <- cor(test_data$medv, pred_ols)^2

cat("Baseline models created successfully.\n")

## Baseline models created successfully.

9.1 Visualisasi Grid Search Results

# Create a data frame from the results
tune_df <- data.frame(tune_result$performances)

# Load required libraries for visualization
library(reshape2)
library(dplyr)

# Prepare data for visualization
tune_df_plot <- tune_df %>%
  mutate(cost_factor = as.factor(cost),
         gamma_factor = as.factor(gamma),
         epsilon_factor = as.factor(epsilon))

# 1. Box plot showing error distribution by parameter
p1 <- ggplot(tune_df_plot, aes(x = cost_factor, y = error, fill = cost_factor)) +
  geom_boxplot() +
  labs(title = "Error Distribution by Cost Parameter",
       x = "Cost", y = "Cross-Validation Error") +
  theme_minimal() +
  theme(legend.position = "none")

p2 <- ggplot(tune_df_plot, aes(x = gamma_factor, y = error, fill = gamma_factor)) +
  geom_boxplot() +
  labs(title = "Error Distribution by Gamma Parameter",
       x = "Gamma", y = "Cross-Validation Error") +
  theme_minimal() +
  theme(legend.position = "none")

p3 <- ggplot(tune_df_plot, aes(x = epsilon_factor, y = error, fill = epsilon_factor)) +
  geom_boxplot() +
  labs(title = "Error Distribution by Epsilon Parameter",
       x = "Epsilon", y = "Cross-Validation Error") +
  theme_minimal() +
  theme(legend.position = "none")

grid.arrange(p1, p2, p3, ncol = 3)

# 2. Heatmap visualization for each epsilon value
epsilon_values_grid <- unique(tune_df$epsilon)

for(eps in epsilon_values_grid) {
  subset_data <- tune_df[tune_df$epsilon == eps, ]
  
  p <- ggplot(subset_data, aes(x = factor(cost), y = factor(gamma), fill = error)) +
    geom_tile() +
    scale_fill_gradient(low = "blue", high = "red", name = "CV Error") +
    labs(title = paste("Parameter Tuning Heatmap (ε =", eps, ")"),
         x = "Cost", y = "Gamma") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  print(p)
}

# 3. Comprehensive 3D-like visualization using facets
ggplot(tune_df_plot, aes(x = cost, y = gamma, fill = error)) +
  geom_tile() +
  facet_wrap(~epsilon_factor, labeller = label_both) +
  scale_fill_gradient(low = "blue", high = "red", name = "CV Error") +
  scale_x_log10() +
  scale_y_log10() +
  labs(title = "Grid Search Results - All Parameter Combinations",
       x = "Cost (log scale)", y = "Gamma (log scale)") +
  theme_minimal()

# Show top 10 parameter combinations
cat("\n=== TOP 10 PARAMETER COMBINATIONS ===\n")

## 
## === TOP 10 PARAMETER COMBINATIONS ===

top_results <- tune_df[order(tune_df$error), ][1:10, ]
print(top_results)

##    cost gamma epsilon    error dispersion
## 19   10  0.01    0.10 10.38781   7.542626
## 35   10  0.01    0.20 10.62208   7.690912
## 3    10  0.01    0.01 10.87840   7.534703
## 20  100  0.01    0.10 11.11934   6.885117
## 4   100  0.01    0.01 11.20698   7.160629
## 36  100  0.01    0.20 11.68139   7.120058
## 2     1  0.01    0.01 13.96195   9.384293
## 18    1  0.01    0.10 14.16763   9.529190
## 34    1  0.01    0.20 14.78303   9.787897
## 23   10  0.10    0.10 15.99168   7.515147

9.2 Traditional Heatmap dengan Base R

# Traditional heatmap visualization dengan base R
par(mfrow = c(1, length(epsilon_values_grid)))
for(eps in epsilon_values_grid) {
  subset_data <- tune_df[tune_df$epsilon == eps, ]
  
  # Create matrix for heatmap
  cost_levels <- sort(unique(subset_data$cost))
  gamma_levels <- sort(unique(subset_data$gamma))
  
  error_matrix <- matrix(NA, nrow = length(cost_levels), ncol = length(gamma_levels))
  rownames(error_matrix) <- cost_levels
  colnames(error_matrix) <- gamma_levels
  
  for(i in 1:nrow(subset_data)) {
    row_idx <- which(cost_levels == subset_data$cost[i])
    col_idx <- which(gamma_levels == subset_data$gamma[i])
    error_matrix[row_idx, col_idx] <- subset_data$error[i]
  }
  
  # Plot heatmap
  image(log10(cost_levels), log10(gamma_levels), error_matrix,
        main = paste("Grid Search Results (ε =", eps, ")"),
        xlab = "log10(Cost)", ylab = "log10(Gamma)",
        col = heat.colors(20))
  contour(log10(cost_levels), log10(gamma_levels), error_matrix, add = TRUE)
}

par(mfrow = c(1, 1))

10. Interpretasi Parameter dan Model Analysis

10.1 Detailed Parameter Analysis

Parameter Cost (C):

- Cost tinggi: Model lebih kompleks, lebih sedikit support vectors, risiko overfitting.

- Cost rendah: Model lebih sederhana, lebih banyak support vectors, risiko underfitting.

- Optimal range: Biasanya antara 1-100 untuk kebanyakan dataset.

Parameter Gamma (γ):

- Gamma tinggi: Decision boundary lebih kompleks, setiap training point punya influence kecil.

- Gamma rendah: Decision boundary lebih smooth, influence training point lebih luas.

- Optimal range: Biasanya antara 0.01-1 untuk dataset normalized.

Parameter Epsilon (ε):

- Epsilon besar: Toleransi error lebih besar, lebih sedikit support vectors.

- Epsilon kecil: Toleransi error lebih kecil, lebih banyak support vectors.

- Optimal range: Biasanya antara 0.01-0.2 untuk data yang di-scale.

# Detailed analysis of best parameters
best_params <- tune_result$best.parameters
cat("=== ANALISIS PARAMETER OPTIMAL ===\n")

## === ANALISIS PARAMETER OPTIMAL ===

cat("Best Cost:", best_params$cost, "- Memberikan balance optimal antara bias dan variance\n")

## Best Cost: 10 - Memberikan balance optimal antara bias dan variance

cat("Best Gamma:", best_params$gamma, "- Memberikan complexity yang tepat untuk decision boundary\n")

## Best Gamma: 0.01 - Memberikan complexity yang tepat untuk decision boundary

cat("Best Epsilon:", best_params$epsilon, "- Memberikan toleransi error yang optimal\n")

## Best Epsilon: 0.1 - Memberikan toleransi error yang optimal

10.2 Model Performance Comparison

# Ensure all required variables are defined
# Create baseline models for comparison if they don't exist

# Initialize variables to avoid 'object not found' errors
if (!exists("rmse_linear") || !exists("svr_linear")) {
  cat("Creating SVR Linear model...\n")
  svr_linear <- svm(medv ~ ., data = train_data, 
                    type = "eps-regression", 
                    kernel = "linear")
  pred_linear <- predict(svr_linear, test_data)
  rmse_linear <- sqrt(mean((test_data$medv - pred_linear)^2))
  mae_linear <- mean(abs(test_data$medv - pred_linear))
  r2_linear <- cor(test_data$medv, pred_linear)^2
}

if (!exists("rmse_rbf") || !exists("svr_rbf")) {
  cat("Creating SVR RBF model...\n")
  svr_rbf <- svm(medv ~ ., data = train_data, 
                 type = "eps-regression", 
                 kernel = "radial")
  pred_rbf <- predict(svr_rbf, test_data)
  rmse_rbf <- sqrt(mean((test_data$medv - pred_rbf)^2))
  mae_rbf <- mean(abs(test_data$medv - pred_rbf))
  r2_rbf <- cor(test_data$medv, pred_rbf)^2
}

if (!exists("rmse_ols")) {
  cat("Creating OLS model...\n")
  ols_model <- lm(medv ~ ., data = train_data)
  pred_ols <- predict(ols_model, test_data)
  rmse_ols <- sqrt(mean((test_data$medv - pred_ols)^2))
  mae_ols <- mean(abs(test_data$medv - pred_ols))
  r2_ols <- cor(test_data$medv, pred_ols)^2
}

# Print the metrics to verify they exist
cat("Metrics Summary:\n")

## Metrics Summary:

cat("RMSE Linear:", rmse_linear, "\n")

## RMSE Linear: 6.867489

cat("RMSE RBF:", rmse_rbf, "\n")

## RMSE RBF: 5.841939

cat("RMSE Best:", rmse_best, "\n")

## RMSE Best: 5.813745

cat("RMSE OLS:", rmse_ols, "\n")

## RMSE OLS: 6.420796

# Comprehensive final comparison table including all models
final_comparison <- data.frame(
  Model = c("SVR Linear", "SVR RBF (Default)", "SVR RBF (Optimized)", "OLS Regression"),
  RMSE = round(c(rmse_linear, rmse_rbf, rmse_best, rmse_ols), 4),
  MAE = round(c(mae_linear, mae_rbf, mae_best, mae_ols), 4),
  R_squared = round(c(r2_linear, r2_rbf, r2_best, r2_ols), 4),
  Support_Vectors = c(svr_linear$tot.nSV, svr_rbf$tot.nSV, best_svr$tot.nSV, "-")
)

# Display table with highlighting (using base kable if kableExtra not available)
if (require(kableExtra, quietly = TRUE)) {
  kable(final_comparison, caption = "Comprehensive Model Performance Comparison") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    row_spec(which.min(final_comparison$RMSE), bold = T, color = "white", background = "#D7261E")
} else {
  kable(final_comparison, digits = 4, caption = "Comprehensive Model Performance Comparison")
}

Comprehensive Model Performance Comparison

Model	RMSE	MAE	R_squared	Support_Vectors
SVR Linear	6.8675	4.0711	0.6441	294
SVR RBF (Default)	5.8419	3.0200	0.7458	263
SVR RBF (Optimized)	5.8137	3.0105	0.7394	267
OLS Regression	6.4208	4.0287	0.6790

# Identify best model
best_model_idx <- which.min(final_comparison$RMSE)
best_model_name <- final_comparison$Model[best_model_idx]

cat("=== RINGKASAN HASIL ANALISIS ===\n")

## === RINGKASAN HASIL ANALISIS ===

cat("Model Terbaik:", best_model_name, "\n")

## Model Terbaik: SVR RBF (Optimized)

cat("RMSE Terbaik:", final_comparison$RMSE[best_model_idx], "\n")

## RMSE Terbaik: 5.8137

cat("MAE Terbaik:", final_comparison$MAE[best_model_idx], "\n")

## MAE Terbaik: 3.0105

cat("R² Terbaik:", final_comparison$R_squared[best_model_idx], "\n")

## R² Terbaik: 0.7394

# Calculate improvement from baseline
baseline_rmse <- rmse_ols
best_rmse <- final_comparison$RMSE[best_model_idx]
improvement <- ((baseline_rmse - best_rmse) / baseline_rmse) * 100

cat("Improvement dari OLS Regression:", round(improvement, 2), "%\n")

## Improvement dari OLS Regression: 9.46 %

cat("Support Vectors dalam model optimal:", final_comparison$Support_Vectors[best_model_idx], "\n")

## Support Vectors dalam model optimal: 267

10.3 Performance Visualization

# Create performance comparison plots
if (require(tidyr, quietly = TRUE)) {
  # Reshape data for plotting
  performance_long <- final_comparison %>%
    select(-Support_Vectors) %>%
    pivot_longer(cols = c(RMSE, MAE, R_squared), 
                 names_to = "Metric", 
                 values_to = "Value")
  
  # Create grouped bar plot
  ggplot(performance_long, aes(x = Model, y = Value, fill = Metric)) +
    geom_bar(stat = "identity", position = "dodge") +
    facet_wrap(~Metric, scales = "free_y") +
    labs(title = "Model Performance Comparison Across All Metrics",
         x = "Model", y = "Performance Value") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none")
} else {
  # Fallback visualization without tidyr
  par(mfrow = c(1, 3))
  barplot(final_comparison$RMSE, names.arg = final_comparison$Model, 
          main = "RMSE Comparison", las = 2, cex.names = 0.8)
  barplot(final_comparison$MAE, names.arg = final_comparison$Model, 
          main = "MAE Comparison", las = 2, cex.names = 0.8)
  barplot(final_comparison$R_squared, names.arg = final_comparison$Model, 
          main = "R² Comparison", las = 2, cex.names = 0.8)
  par(mfrow = c(1, 1))
}

# Support Vectors comparison
sv_data <- final_comparison[final_comparison$Support_Vectors != "-", ]
sv_data$Support_Vectors <- as.numeric(sv_data$Support_Vectors)

ggplot(sv_data, aes(x = Model, y = Support_Vectors, fill = Model)) +
  geom_bar(stat = "identity") +
  labs(title = "Number of Support Vectors by Model",
       x = "Model", y = "Number of Support Vectors") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

11. Kesimpulan

# Final model summary
cat("Parameters:\n")

## Parameters:

cat("  - Cost (C):", best_params$cost, "\n")

##   - Cost (C): 10

cat("  - Gamma (γ):", best_params$gamma, "\n")

##   - Gamma (γ): 0.01

cat("  - Epsilon (ε):", best_params$epsilon, "\n")

##   - Epsilon (ε): 0.1

cat("Performance:\n")

## Performance:

cat("  - RMSE:", round(rmse_best, 4), "\n")

##   - RMSE: 5.8137

cat("  - MAE:", round(mae_best, 4), "\n")

##   - MAE: 3.0105

cat("  - R²:", round(r2_best, 4), "\n")

##   - R²: 0.7394

cat("  - Support Vectors:", best_svr$tot.nSV, "dari", nrow(train_data), "training points\n")

##   - Support Vectors: 267 dari 381 training points

1. Performa Model:

- SVR dengan parameter ter-optimasi menunjukkan performa terbaik dibandingkan SVR standar dan OLS regression.

• Kernel RBF lebih efektif daripada kernel linear untuk dataset Boston Housing karena mampu menangkap pola nonlinear.

- Grid search berhasil meningkatkan akurasi prediksi secara signifikan.

2. Karakteristik SVR:

- Epsilon-tube memberikan toleransi error yang fleksibel, memungkinkan model robust terhadap noise.

- Support vectors menentukan model, membuat SVR efisien dan robust terhadap outliers.

- Jumlah support vectors dapat dikontrol melalui parameter C dan epsilon.

3. Keunggulan SVR:

- Robust terhadap outliers karena hanya support vectors yang mempengaruhi model.

- Dapat menangani data nonlinear dengan kernel yang sesuai.

- Memberikan solusi global optimum (tidak terjebak di local minimum).

- Efektif untuk dataset dengan dimensi tinggi.

4. Aplikasi Praktis:

- Model SVR dapat digunakan untuk prediksi harga rumah dengan akurasi tinggi.

- Parameter tuning sangat penting untuk mendapatkan performa optimal.

- Preprocessing data (normalisasi) krusial untuk performa SVR.