1. Opis problemu

Celem projektu jest zbudowanie modelu predykcyjnego, który na podstawie cech fizycznych diamentu przewiduje jego cenę rynkową w USD. Dane zawierają ponad 53 000 diamentów opisanych przez masę (carat), szlif (cut), kolor (color), czystość (clarity) oraz wymiary geometryczne.

Zmienna docelowa to price (zmienna ciągła), więc jest to problem regresji.

W projekcie porównujemy dziewięć modeli: OLS jako punkt odniesienia, Ridge i Lasso (modele liniowe z regularyzacją), drzewo regresyjne (CART), bagging, Random Forest (z tuningiem), Gradient Boosting (GBM) ze strojeniem sekwencyjnym oraz sieć neuronową (nnet).

library(tidyverse)
library(caret)
library(corrplot)
library(rpart)
library(randomForest)
library(ranger)
library(gbm)
library(nnet)
library(glmnet)

2. Eksploracja danych

diamonds.raw <- read.csv("diamonds.csv", row.names = 1)
glimpse(diamonds.raw)
## Rows: 53,940
## Columns: 10
## $ carat   <dbl> 0.23, 0.21, 0.23, 0.29, 0.31, 0.24, 0.24, 0.26, 0.22, 0.23, 0.…
## $ cut     <chr> "Ideal", "Premium", "Good", "Premium", "Good", "Very Good", "V…
## $ color   <chr> "E", "E", "E", "I", "J", "J", "I", "H", "E", "H", "J", "J", "F…
## $ clarity <chr> "SI2", "SI1", "VS1", "VS2", "SI2", "VVS2", "VVS1", "SI1", "VS2…
## $ depth   <dbl> 61.5, 59.8, 56.9, 62.4, 63.3, 62.8, 62.3, 61.9, 65.1, 59.4, 64…
## $ table   <dbl> 55, 61, 65, 58, 58, 57, 57, 55, 61, 61, 55, 56, 61, 54, 62, 58…
## $ price   <int> 326, 326, 327, 334, 335, 336, 336, 337, 337, 338, 339, 340, 34…
## $ x       <dbl> 3.95, 3.89, 4.05, 4.20, 4.34, 3.94, 3.95, 4.07, 3.87, 4.00, 4.…
## $ y       <dbl> 3.98, 3.84, 4.07, 4.23, 4.35, 3.96, 3.98, 4.11, 3.78, 4.05, 4.…
## $ z       <dbl> 2.43, 2.31, 2.31, 2.63, 2.75, 2.48, 2.47, 2.53, 2.49, 2.39, 2.…
summary(diamonds.raw)
##      carat            cut               color             clarity         
##  Min.   :0.2000   Length:53940       Length:53940       Length:53940      
##  1st Qu.:0.4000   Class :character   Class :character   Class :character  
##  Median :0.7000   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :0.7979                                                           
##  3rd Qu.:1.0400                                                           
##  Max.   :5.0100                                                           
##      depth           table           price             x         
##  Min.   :43.00   Min.   :43.00   Min.   :  326   Min.   : 0.000  
##  1st Qu.:61.00   1st Qu.:56.00   1st Qu.:  950   1st Qu.: 4.710  
##  Median :61.80   Median :57.00   Median : 2401   Median : 5.700  
##  Mean   :61.75   Mean   :57.46   Mean   : 3933   Mean   : 5.731  
##  3rd Qu.:62.50   3rd Qu.:59.00   3rd Qu.: 5324   3rd Qu.: 6.540  
##  Max.   :79.00   Max.   :95.00   Max.   :18823   Max.   :10.740  
##        y                z         
##  Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 4.720   1st Qu.: 2.910  
##  Median : 5.710   Median : 3.530  
##  Mean   : 5.735   Mean   : 3.539  
##  3rd Qu.: 6.540   3rd Qu.: 4.040  
##  Max.   :58.900   Max.   :31.800
cat("Liczba obserwacji:", nrow(diamonds.raw), "\n")
## Liczba obserwacji: 53940
cat("Duplikaty:", sum(duplicated(diamonds.raw)), "\n")
## Duplikaty: 146
cat("Braki danych:\n")
## Braki danych:
colSums(is.na(diamonds.raw))
##   carat     cut   color clarity   depth   table   price       x       y       z 
##       0       0       0       0       0       0       0       0       0       0

Dane są kompletne — brak braków. Są natomiast duplikaty, które usuniemy przy czyszczeniu.

Rozkład zmiennej docelowej

par(mfrow = c(1, 2))
hist(diamonds.raw$price, breaks = 60, col = "steelblue", border = "white",
     main = "Cena (surowa)", xlab = "Cena (USD)")
hist(log(diamonds.raw$price), breaks = 60, col = "coral", border = "white",
     main = "log(Cena)", xlab = "log(Cena)")

par(mfrow = c(1, 1))

Rozkład ceny jest silnie prawostronnie skośny. Po transformacji logarytmicznej rozkład jest znacznie bardziej symetryczny. Dlatego będziemy modelować log(price) i przeliczać predykcje z powrotem przez exp().

Carat vs cena

par(mfrow = c(1, 2))
plot(diamonds.raw$carat, diamonds.raw$price, pch = ".",
     col = rgb(0, 0, 1, 0.15), main = "Carat vs Cena",
     xlab = "Carat", ylab = "Cena (USD)")
plot(log(diamonds.raw$carat), log(diamonds.raw$price), pch = ".",
     col = rgb(0, 0, 1, 0.15), main = "log(Carat) vs log(Cena)",
     xlab = "log(Carat)", ylab = "log(Cena)")

par(mfrow = c(1, 1))

W skali log-log zależność jest prawie liniowa. To oznacza, że log(carat) będzie znacznie lepszą cechą niż surowy carat.

Zmienne kategoryczne

par(mfrow = c(1, 3))
boxplot(price ~ cut,     data = diamonds.raw, main = "Cena wg cut",
        col = "lightblue", las = 2)
boxplot(price ~ color,   data = diamonds.raw, main = "Cena wg color",
        col = "lightgreen", las = 2)
boxplot(price ~ clarity, data = diamonds.raw, main = "Cena wg clarity",
        col = "lightyellow", las = 2)

par(mfrow = c(1, 1))

Ciekawa obserwacja: diamenty z lepszym szlifem (Ideal) wcale nie są najdroższe. Wynika to z tego, że małe kamienie częściej otrzymują staranny szlif, a masę (carat) jest głównym determinantem ceny.

Macierz korelacji

corrplot(cor(diamonds.raw[, c("carat", "depth", "table", "price")]),
         method = "color", addCoef.col = "black", tl.cex = 0.9,
         title = "Korelacje zmiennych numerycznych", mar = c(0, 0, 1, 0))

carat ma bardzo silną korelację z price. depth i table mają słaby związek z ceną.

3. Preprocessing i feature engineering

Czyszczenie

cat("Zera w wymiarach:",
    sum(diamonds.raw$x == 0 | diamonds.raw$y == 0 | diamonds.raw$z == 0), "\n")
## Zera w wymiarach: 20
cat("Zakres depth:", range(diamonds.raw$depth), "\n")
## Zakres depth: 43 79
diamonds.clean <- diamonds.raw %>%
  distinct() %>%
  filter(x > 0, y > 0, z > 0) %>%
  filter(depth >= 50, depth <= 75) %>%
  filter(y < 30, z < 20)

cat("Po czyszczeniu:", nrow(diamonds.clean),
    "| usunięto:", nrow(diamonds.raw) - nrow(diamonds.clean), "\n")
## Po czyszczeniu: 53767 | usunięto: 173

Feature engineering

diamonds.prep <- diamonds.clean %>%
  mutate(
    log.carat   = log(carat),
    cut.num     = as.integer(factor(cut,
                    levels = c("Fair", "Good", "Very Good", "Premium", "Ideal"))),
    color.num   = as.integer(factor(color,
                    levels = c("J", "I", "H", "G", "F", "E", "D"))),
    clarity.num = as.integer(factor(clarity,
                    levels = c("I1", "SI2", "SI1", "VS2", "VS1", "VVS2", "VVS1", "IF"))),
    log.price   = log(price)
  ) %>%
  select(log.carat, depth, table, cut.num, color.num, clarity.num, price, log.price)

glimpse(diamonds.prep)
## Rows: 53,767
## Columns: 8
## $ log.carat   <dbl> -1.469676, -1.560648, -1.469676, -1.237874, -1.171183, -1.…
## $ depth       <dbl> 61.5, 59.8, 56.9, 62.4, 63.3, 62.8, 62.3, 61.9, 65.1, 59.4…
## $ table       <dbl> 55, 61, 65, 58, 58, 57, 57, 55, 61, 61, 55, 56, 61, 54, 62…
## $ cut.num     <int> 5, 4, 2, 4, 2, 3, 3, 3, 1, 3, 2, 5, 4, 5, 4, 4, 5, 2, 2, 3…
## $ color.num   <int> 6, 6, 6, 2, 1, 1, 2, 3, 6, 3, 1, 1, 5, 1, 6, 6, 2, 1, 1, 1…
## $ clarity.num <int> 2, 3, 5, 4, 2, 6, 7, 3, 4, 5, 3, 5, 3, 2, 2, 1, 2, 3, 3, 3…
## $ price       <int> 326, 326, 327, 334, 335, 336, 336, 337, 337, 338, 339, 340…
## $ log.price   <dbl> 5.786897, 5.786897, 5.789960, 5.811141, 5.814131, 5.817111…

Zmienne kategoryczne mają naturalny porządek jakości (Fair < Good < … < Ideal), dlatego kodujemy je jako liczby całkowite zamiast dummy variables.

Podział train / test

# Split before you model!
idx            <- sample(seq_len(nrow(diamonds.prep)), size = 0.8 * nrow(diamonds.prep))
diamonds.train <- diamonds.prep[idx, ]
diamonds.test  <- diamonds.prep[-idx, ]

cat("Train:", nrow(diamonds.train), "| Test:", nrow(diamonds.test), "\n")
## Train: 43013 | Test: 10754
real.train    <- diamonds.train$price
real.test     <- diamonds.test$price
diamonds.vars <- c("log.carat", "depth", "table", "cut.num", "color.num", "clarity.num")
fm            <- as.formula(paste("log.price ~", paste(diamonds.vars, collapse = " + ")))

Funkcja do metryk

regresja.metryki <- function(real, predicted) {
  MSE   <- mean((real - predicted)^2)
  RMSE  <- sqrt(MSE)
  MAE   <- mean(abs(real - predicted))
  MedAE <- median(abs(real - predicted))
  MAPE  <- mean(abs(real - predicted) / real)
  MSLE  <- mean((log(1 + real) - log(1 + predicted))^2)
  R2    <- cor(predicted, real)^2
  data.frame(MSE   = round(MSE, 2),   RMSE  = round(RMSE, 2),
             MAE   = round(MAE, 2),   MedAE = round(MedAE, 2),
             MAPE  = round(MAPE, 4),  MSLE  = round(MSLE, 6),
             R2    = round(R2, 4))
}

wykres.pred <- function(real, predicted, tytul = "") {
  plot(predicted, real, pch = ".", col = rgb(0, 0, 1, 0.2),
       main = tytul, xlab = "Predykcja (USD)", ylab = "Wartość rzeczywista (USD)")
  abline(0, 1, col = "red", lwd = 2)
}

4. Model benchmarkowy: OLS

diamonds.ols <- lm(fm, data = diamonds.train)
summary(diamonds.ols)
## 
## Call:
## lm(formula = fm, data = diamonds.train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.12273 -0.08989  0.00267  0.09543  1.42455 
## 
## Coefficients:
##               Estimate Std. Error  t value             Pr(>|t|)    
## (Intercept)  7.8725263  0.0521998  150.815 < 0.0000000000000002 ***
## log.carat    1.8776061  0.0013795 1361.069 < 0.0000000000000002 ***
## depth       -0.0043805  0.0005737   -7.635    0.000000000000023 ***
## table       -0.0005726  0.0004003   -1.430                0.153    
## cut.num      0.0302870  0.0007727   39.199 < 0.0000000000000002 ***
## color.num    0.0778357  0.0004338  179.437 < 0.0000000000000002 ***
## clarity.num  0.1228702  0.0004749  258.748 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1462 on 43006 degrees of freedom
## Multiple R-squared:  0.9793, Adjusted R-squared:  0.9793 
## F-statistic: 3.383e+05 on 6 and 43006 DF,  p-value: < 0.00000000000000022
pred.ols.train <- exp(predict(diamonds.ols))
pred.ols.test  <- exp(predict(diamonds.ols, newdata = diamonds.test))

cat("=== OLS: train vs test ===\n")
## === OLS: train vs test ===
rbind(train = regresja.metryki(real.train, pred.ols.train),
      test  = regresja.metryki(real.test,  pred.ols.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 918758.4 958.52 462.40 183.03 0.1148 0.021349 0.9427
## test  905937.2 951.81 459.21 186.64 0.1138 0.020809 0.9426
wykres.pred(real.test, pred.ols.test, "OLS: predykcja vs wartość rzeczywista (test)")

R² powyżej 0.94 już na regresji liniowej — to zasługa transformacji log(carat). Bez tej transformacji R² wynosi poniżej 0.20.

Analiza reszt OLS

Reszty to różnica między wartością rzeczywistą a prognozą modelu. Ich analiza pozwala ocenić, czy model ma systematyczne błędy.

errors.train <- real.train - pred.ols.train
errors.test  <- real.test  - pred.ols.test

cat("=== Reszty OLS — zbiór treningowy ===\n")
## === Reszty OLS — zbiór treningowy ===
summary(errors.train)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -31418.887   -156.034      4.116     50.769    219.771   8514.542
cat("Średnia reszt train:", round(mean(errors.train), 2), "USD\n")
## Średnia reszt train: 50.77 USD
cat("=== Reszty OLS — zbiór testowy ===\n")
## === Reszty OLS — zbiór testowy ===
summary(errors.test)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -21097.652   -151.192      8.513     56.477    232.003   6974.537
cat("Średnia reszt test: ", round(mean(errors.test), 2), "USD\n")
## Średnia reszt test:  56.48 USD
par(mfrow = c(1, 2))
hist(errors.train, breaks = 60, col = "steelblue", border = "white",
     main = "Rozkład reszt — train", xlab = "Błąd predykcji (USD)")
abline(v = 0, col = "red", lwd = 2)
hist(errors.test, breaks = 60, col = "steelblue", border = "white",
     main = "Rozkład reszt — test", xlab = "Błąd predykcji (USD)")
abline(v = 0, col = "red", lwd = 2)

par(mfrow = c(1, 1))
par(mfrow = c(1, 2))
plot(pred.ols.train, errors.train, pch = ".", col = rgb(0, 0, 1, 0.15),
     main = "Reszty vs predykcja — train",
     xlab = "Predykcja (USD)", ylab = "Reszta (USD)")
abline(h = 0, col = "red", lwd = 2)
plot(pred.ols.test, errors.test, pch = ".", col = rgb(0, 0, 1, 0.15),
     main = "Reszty vs predykcja — test",
     xlab = "Predykcja (USD)", ylab = "Reszta (USD)")
abline(h = 0, col = "red", lwd = 2)

par(mfrow = c(1, 1))
qqnorm(errors.train, pch = ".", col = rgb(0, 0, 1, 0.15),
       main = "QQ-plot reszt OLS (train)")
qqline(errors.train, col = "red", lwd = 2)

Interpretacja analizy reszt:

Średnia reszt wynosi ok. +50 USD, co oznacza niewielki systematyczny błąd w górę — model średnio nieznacznie zaniża ceny. Rozkład reszt jest jednak prawostronnie skośny: dominuje dużo małych błędów ujemnych i nieliczne, ale bardzo duże błędy dodatnie. Model ma tendencję do zaniżania ceny najdroższych diamentów — widać to wyraźnie na wykresie reszty vs predykcja, gdzie przy wysokich predykcjach reszty rosną. Wykres QQ-plot potwierdza odchylenie od normalności w prawym ogonie — to typowe dla danych cenowych z wartościami odstającymi. Wniosek: model liniowy nie w pełni uchwytuje nieliniowych zależności przy bardzo drogich kamieniach. Modele drzewiaste (RF, GBM) powinny radzić sobie z tym lepiej.

5. Ridge regression

Ridge to regresja liniowa z regularyzacją L2. Dodaje karę za duże wartości współczynników, co stabilizuje model przy silnej współliniowości zmiennych. Parametr lambda kontroluje siłę regularyzacji — im większy, tym mocniej współczynniki są “skraczane” w kierunku zera.

ctrl.reg <- trainControl(method = "cv", number = 5)

diamonds.ridge <- train(fm,
                        data      = diamonds.train,
                        method    = "glmnet",
                        family    = "gaussian",
                        trControl = ctrl.reg,
                        tuneGrid  = expand.grid(alpha  = 0,
                                                lambda = seq(0.0001, 0.5, length = 30)))
plot(diamonds.ridge, main = "Ridge: RMSE vs lambda (CV)")

cat("Ridge — najlepsza lambda:", round(diamonds.ridge$bestTune$lambda, 5), "\n")
## Ridge — najlepsza lambda: 0.08629
cat("Współczynniki Ridge:\n")
## Współczynniki Ridge:
print(coef(diamonds.ridge$finalModel, s = diamonds.ridge$bestTune$lambda))
## 7 x 1 sparse Matrix of class "dgCMatrix"
##             s=0.08628966
## (Intercept)  7.765673588
## log.carat    1.652542202
## depth       -0.004047955
## table        0.004337882
## cut.num      0.025795184
## color.num    0.052192449
## clarity.num  0.084751877
pred.ridge.train <- exp(predict(diamonds.ridge, newdata = diamonds.train))
pred.ridge.test  <- exp(predict(diamonds.ridge, newdata = diamonds.test))

cat("=== Ridge: train vs test ===\n")
## === Ridge: train vs test ===
rbind(train = regresja.metryki(real.train, pred.ridge.train),
      test  = regresja.metryki(real.test,  pred.ridge.test))
##           MSE    RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 1595633 1263.18 647.54 210.49 0.1484 0.034851 0.9343
## test  1561708 1249.68 642.96 213.06 0.1471 0.034283 0.9341

6. Lasso regression

Lasso (regularyzacja L1) różni się od Ridge tym, że może wyzerować współczynniki przy wybranych zmiennych, działając jednocześnie jako metoda selekcji zmiennych.

diamonds.lasso <- train(fm,
                        data      = diamonds.train,
                        method    = "glmnet",
                        family    = "gaussian",
                        trControl = ctrl.reg,
                        tuneGrid  = expand.grid(alpha  = 1,
                                                lambda = seq(0.0001, 0.1, length = 30)))
plot(diamonds.lasso, main = "Lasso: RMSE vs lambda (CV)")

cat("Lasso — najlepsza lambda:", round(diamonds.lasso$bestTune$lambda, 5), "\n")
## Lasso — najlepsza lambda: 0.0001
cat("Współczynniki Lasso (. = wyzerowane):\n")
## Współczynniki Lasso (. = wyzerowane):
print(coef(diamonds.lasso$finalModel, s = diamonds.lasso$bestTune$lambda))
## 7 x 1 sparse Matrix of class "dgCMatrix"
##                 s=0.0001
## (Intercept)  7.761257210
## log.carat    1.867800182
## depth       -0.002757316
## table        .          
## cut.num      0.028984976
## color.num    0.075512377
## clarity.num  0.120296364
pred.lasso.train <- exp(predict(diamonds.lasso, newdata = diamonds.train))
pred.lasso.test  <- exp(predict(diamonds.lasso, newdata = diamonds.test))

cat("=== Lasso: train vs test ===\n")
## === Lasso: train vs test ===
rbind(train = regresja.metryki(real.train, pred.lasso.train),
      test  = regresja.metryki(real.test,  pred.lasso.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 931775.8 965.29 466.70 181.66 0.1149 0.021396 0.9419
## test  917630.5 957.93 463.41 184.53 0.1139 0.020842 0.9418
cat("=== Porównanie: OLS vs Ridge vs Lasso (test) ===\n")
## === Porównanie: OLS vs Ridge vs Lasso (test) ===
rbind(OLS   = regresja.metryki(real.test, pred.ols.test),
      Ridge = regresja.metryki(real.test, pred.ridge.test),
      Lasso = regresja.metryki(real.test, pred.lasso.test))
##             MSE    RMSE    MAE  MedAE   MAPE     MSLE     R2
## OLS    905937.2  951.81 459.21 186.64 0.1138 0.020809 0.9426
## Ridge 1561708.0 1249.68 642.96 213.06 0.1471 0.034283 0.9341
## Lasso  917630.5  957.93 463.41 184.53 0.1139 0.020842 0.9418

Interpretacja Ridge i Lasso:

Ridge osiąga gorsze RMSE niż OLS. Regularyzacja “skurczyła” współczynniki zbyt mocno, wprowadzając bias. Dzieje się tak dlatego, że model OLS z log(carat) jest już dobrze skalibrowany — regularyzacja w tym przypadku nie pomaga, bo nie ma problemu ze współliniowością ani niestabilnością. Lasso przy optymalnej lambdzie wyzerowało zmienną table (szerokość rondisty diamentu). Jest to ciekawy wynik potwierdzający obserwację z macierzy korelacji — table ma słaby związek z ceną i można ją pominąć bez straty jakości predykcji.

7. Drzewo regresyjne (CART)

diamonds.tree <- rpart(fm, data = diamonds.train, method = "anova",
                       control = rpart.control(cp = 0.0001, minsplit = 20,
                                               minbucket = 7, maxdepth = 15))
plotcp(diamonds.tree)

best.cp       <- diamonds.tree$cptable[which.min(diamonds.tree$cptable[, "xerror"]), "CP"]
diamonds.tree <- prune(diamonds.tree, cp = best.cp)
cat("Optymalne cp:", round(best.cp, 6), "\n")
## Optymalne cp: 0.0001
cat("Liczba liści:", sum(diamonds.tree$frame$var == "<leaf>"), "\n")
## Liczba liści: 81
plot(diamonds.tree, uniform = TRUE, main = "Drzewo regresyjne (przycięte)")
text(diamonds.tree, use.n = FALSE, cex = 0.45)

pred.tree.train <- exp(predict(diamonds.tree, newdata = diamonds.train))
pred.tree.test  <- exp(predict(diamonds.tree, newdata = diamonds.test))

cat("=== CART: train vs test ===\n")
## === CART: train vs test ===
rbind(train = regresja.metryki(real.train, pred.tree.train),
      test  = regresja.metryki(real.test,  pred.tree.test))
##          MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 675824 822.09 447.62 180.84 0.1159 0.021548 0.9579
## test  670064 818.57 449.15 186.21 0.1163 0.021759 0.9573
wykres.pred(real.test, pred.tree.test, "CART: predykcja vs rzeczywistość (test)")

8. Bagging

Bagging (Bootstrap Aggregating) buduje wiele drzew na próbkach bootstrapowych i uśrednia predykcje. Kluczowe: mtry = p — każde drzewo widzi wszystkie zmienne przy każdym splicie.

# Trenowanie trwa ok. 10-15 minut na tym zbiorze.
# Wczytujemy już wytrenowany model.

# set.seed(42)
# diamonds.bag <- randomForest(fm, data = diamonds.train,
#                              mtry = length(diamonds.vars),
#                              ntree = 300, importance = TRUE)
# saveRDS(diamonds.bag, "diamonds.bag.rds")

diamonds.bag <- readRDS("diamonds.bag.rds")
diamonds.bag
## 
## Call:
##  randomForest(formula = fm, data = diamonds.train, mtry = length(diamonds.vars),      ntree = 300, importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 300
## No. of variables tried at each split: 6
## 
##           Mean of squared residuals: 0.01160533
##                     % Var explained: 98.87
plot(diamonds.bag, main = "Bagging: błąd OOB vs liczba drzew")

OOB (Out-of-Bag) error to wbudowana estymata błędu na danych, które nie były używane do budowy danego drzewa. Pozwala ocenić jakość modelu bez oddzielnego zbioru walidacyjnego.

pred.bag.train <- exp(predict(diamonds.bag, newdata = diamonds.train))
pred.bag.test  <- exp(predict(diamonds.bag, newdata = diamonds.test))

cat("=== Bagging: train vs test ===\n")
## === Bagging: train vs test ===
rbind(train = regresja.metryki(real.train, pred.bag.train),
      test  = regresja.metryki(real.test,  pred.bag.test))
##             MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train  68941.45 262.57 138.13  66.90 0.0424 0.003433 0.9957
## test  289295.59 537.86 280.34 125.22 0.0795 0.011319 0.9816
imp.bag <- sort(randomForest::importance(diamonds.bag)[, 1], decreasing = TRUE)
barplot(imp.bag, las = 2, col = "steelblue", border = NA,
        main = "Bagging: ważność zmiennych (%IncMSE)",
        ylab = "% wzrost MSE po permutacji")

wykres.pred(real.test, pred.bag.test, "Bagging: predykcja vs rzeczywistość (test)")

9. Random Forest

RF różni się od baggingu: przy każdym splicie losuje podzbiór mtry < p zmiennych. Drzewa są mniej ze sobą skorelowane, co przekłada się na lepszą redukcję wariancji.

# set.seed(42)
# diamonds.rf <- randomForest(fm, data = diamonds.train,
#                             ntree = 300, importance = TRUE)
# saveRDS(diamonds.rf, "diamonds.rf.rds")

diamonds.rf <- readRDS("diamonds.rf.rds")
diamonds.rf
## 
## Call:
##  randomForest(formula = fm, data = diamonds.train, ntree = 300,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 300
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 0.0118603
##                     % Var explained: 98.85
plot(diamonds.rf, main = "Random Forest: błąd OOB vs liczba drzew")

pred.rf.train <- exp(predict(diamonds.rf, newdata = diamonds.train))
pred.rf.test  <- exp(predict(diamonds.rf, newdata = diamonds.test))

cat("=== Random Forest: train vs test ===\n")
## === Random Forest: train vs test ===
rbind(train = regresja.metryki(real.train, pred.rf.train),
      test  = regresja.metryki(real.test,  pred.rf.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 180840.9 425.25 212.08  86.74 0.0569 0.005518 0.9911
## test  394066.2 627.75 313.94 120.56 0.0810 0.011071 0.9784
imp.rf <- sort(randomForest::importance(diamonds.rf)[, 1], decreasing = TRUE)
barplot(imp.rf, las = 2, col = "steelblue", border = NA,
        main = "Random Forest: ważność zmiennych (%IncMSE)",
        ylab = "% wzrost MSE po permutacji")

Tuning Random Forest

ctrl.rf <- trainControl(method = "cv", number = 5)
grid.rf <- expand.grid(mtry          = c(2, 3, 4, 5, 6),
                       min.node.size = c(5, 10, 20),
                       splitrule     = "variance")

# set.seed(42)
# diamonds.rf.tuned <- train(fm, data = diamonds.train, method = "ranger",
#                            trControl = ctrl.rf, tuneGrid = grid.rf,
#                            metric = "RMSE", num.trees = 300, num.threads = 2)
# saveRDS(diamonds.rf.tuned, "diamonds.rf.tuned.rds")

diamonds.rf.tuned <- readRDS("diamonds.rf.tuned.rds")
diamonds.rf.tuned
## Random Forest 
## 
## 43013 samples
##     6 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 34411, 34411, 34409, 34411, 34410 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  RMSE       Rsquared   MAE       
##   2      5             0.1093363  0.9887441  0.08319232
##   2     10             0.1102461  0.9886030  0.08398380
##   2     20             0.1130149  0.9881468  0.08614149
##   3      5             0.1038381  0.9895430  0.07798233
##   3     10             0.1023973  0.9898343  0.07722445
##   3     20             0.1018384  0.9899534  0.07721984
##   4      5             0.1050379  0.9892937  0.07863841
##   4     10             0.1029439  0.9897171  0.07749236
##   4     20             0.1016114  0.9899828  0.07688547
##   5      5             0.1062991  0.9890356  0.07950126
##   5     10             0.1040170  0.9895014  0.07825064
##   5     20             0.1024600  0.9898137  0.07748517
##   6      5             0.1073397  0.9888205  0.08029693
##   6     10             0.1049825  0.9893057  0.07900283
##   6     20             0.1032535  0.9896551  0.07807018
## 
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 4, splitrule = variance
##  and min.node.size = 20.
plot(diamonds.rf.tuned, main = "RF tuning: RMSE vs mtry i min.node.size")

cat("Najlepsze parametry RF:\n")
## Najlepsze parametry RF:
print(diamonds.rf.tuned$bestTune)
##   mtry splitrule min.node.size
## 9    4  variance            20
pred.rf.tuned.train <- exp(predict(diamonds.rf.tuned, newdata = diamonds.train))
pred.rf.tuned.test  <- exp(predict(diamonds.rf.tuned, newdata = diamonds.test))

cat("=== RF tuned: train vs test ===\n")
## === RF tuned: train vs test ===
rbind(train = regresja.metryki(real.train, pred.rf.tuned.train),
      test  = regresja.metryki(real.test,  pred.rf.tuned.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 163682.8 404.58 210.48  90.02 0.0592 0.006056 0.9899
## test  280026.8 529.18 276.88 116.32 0.0755 0.009893 0.9822

10. Gradient Boosting (GBM)

GBM buduje drzewa sekwencyjnie — każde uczy się na błędach poprzedniego. Strojenie odbywa się w trzech krokach.

Krok 1: głębokość drzewa

ctrl.gbm   <- trainControl(method = "cv", number = 5)
grid.gbm.1 <- expand.grid(interaction.depth = c(2, 4, 6), n.trees = 300,
                           shrinkage = 0.05, n.minobsinnode = 10)

# set.seed(42)
# diamonds.gbm.1 <- train(fm, data=diamonds.train, method="gbm",
#                         trControl=ctrl.gbm, tuneGrid=grid.gbm.1,
#                         metric="RMSE", verbose=FALSE)
# saveRDS(diamonds.gbm.1, "diamonds.gbm.1.rds")

diamonds.gbm.1 <- readRDS("diamonds.gbm.1.rds")
plot(diamonds.gbm.1, main = "GBM krok 1: głębokość drzewa")

best.depth <- diamonds.gbm.1$bestTune$interaction.depth
cat("Najlepsza głębokość:", best.depth, "\n")
## Najlepsza głębokość: 6

Krok 2: liczba drzew

grid.gbm.2 <- expand.grid(interaction.depth = best.depth,
                           n.trees = c(100, 200, 400, 600),
                           shrinkage = 0.05, n.minobsinnode = 10)

# set.seed(42)
# diamonds.gbm.2 <- train(fm, data=diamonds.train, method="gbm",
#                         trControl=ctrl.gbm, tuneGrid=grid.gbm.2,
#                         metric="RMSE", verbose=FALSE)
# saveRDS(diamonds.gbm.2, "diamonds.gbm.2.rds")

diamonds.gbm.2 <- readRDS("diamonds.gbm.2.rds")
plot(diamonds.gbm.2, main = "GBM krok 2: liczba drzew")

best.ntrees <- diamonds.gbm.2$bestTune$n.trees
cat("Najlepsza liczba drzew:", best.ntrees, "\n")
## Najlepsza liczba drzew: 600

Krok 3: shrinkage

grid.gbm.3 <- expand.grid(interaction.depth = best.depth, n.trees = best.ntrees,
                           shrinkage = c(0.01, 0.025, 0.05, 0.075, 0.1),
                           n.minobsinnode = 10)

# set.seed(42)
# diamonds.gbm.3 <- train(fm, data=diamonds.train, method="gbm",
#                         trControl=ctrl.gbm, tuneGrid=grid.gbm.3,
#                         metric="RMSE", verbose=FALSE)
# saveRDS(diamonds.gbm.3, "diamonds.gbm.3.rds")

diamonds.gbm.3 <- readRDS("diamonds.gbm.3.rds")
plot(diamonds.gbm.3, main = "GBM krok 3: shrinkage")

cat("Najlepsze parametry GBM:\n")
## Najlepsze parametry GBM:
print(diamonds.gbm.3$bestTune)
##   n.trees interaction.depth shrinkage n.minobsinnode
## 5     600                 6       0.1             10
diamonds.gbm.final <- diamonds.gbm.3

pred.gbm.train <- exp(predict(diamonds.gbm.final, newdata = diamonds.train))
pred.gbm.test  <- exp(predict(diamonds.gbm.final, newdata = diamonds.test))

cat("=== GBM: train vs test ===\n")
## === GBM: train vs test ===
rbind(train = regresja.metryki(real.train, pred.gbm.train),
      test  = regresja.metryki(real.test,  pred.gbm.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 260823.9 510.71 268.27 111.30 0.0732 0.009130 0.9837
## test  275478.0 524.86 276.88 114.76 0.0743 0.009567 0.9825
summary(diamonds.gbm.final$finalModel,
        n.trees = diamonds.gbm.final$bestTune$n.trees,
        plotit = TRUE, main = "GBM: ważność zmiennych")

##                     var     rel.inf
## log.carat     log.carat 95.22676514
## clarity.num clarity.num  3.12408923
## color.num     color.num  1.27504045
## cut.num         cut.num  0.24107544
## depth             depth  0.09869991
## table             table  0.03432983
wykres.pred(real.test, pred.gbm.test, "GBM: predykcja vs rzeczywistość (test)")

11. Sieć neuronowa (nnet)

mn.train <- colMeans(diamonds.train[, diamonds.vars])
sd.train <- apply(diamonds.train[, diamonds.vars], 2, sd)

diamonds.train.s <- diamonds.train
diamonds.test.s  <- diamonds.test
diamonds.train.s[, diamonds.vars] <- scale(diamonds.train[, diamonds.vars],
                                           center = mn.train, scale = sd.train)
diamonds.test.s[, diamonds.vars]  <- scale(diamonds.test[, diamonds.vars],
                                           center = mn.train, scale = sd.train)
# set.seed(42)
# diamonds.nn <- nnet(fm, data = diamonds.train.s, size = 32, decay = 0.001,
#                     maxit = 300, linout = TRUE, trace = FALSE, MaxNWts = 5000)
# saveRDS(diamonds.nn, "diamonds.nn.rds")
# saveRDS(list(mean = mn.train, sd = sd.train), "scaler.nn.rds")

diamonds.nn <- readRDS("diamonds.nn.rds")
cat("Liczba wag w sieci:", length(diamonds.nn$wts), "\n")
## Liczba wag w sieci: 257
pred.nn.train <- exp(predict(diamonds.nn, newdata = diamonds.train.s))
pred.nn.test  <- exp(predict(diamonds.nn, newdata = diamonds.test.s))

cat("=== Sieć neuronowa: train vs test ===\n")
## === Sieć neuronowa: train vs test ===
rbind(train = regresja.metryki(real.train, pred.nn.train),
      test  = regresja.metryki(real.test,  pred.nn.test))
##            MSE   RMSE    MAE  MedAE   MAPE     MSLE     R2
## train 346056.6 588.27 315.58 136.23 0.0857 0.011933 0.9784
## test  372305.5 610.17 320.81 137.99 0.0857 0.012016 0.9763
wykres.pred(real.test, pred.nn.test, "nnet: predykcja vs rzeczywistość (test)")

12. Porównanie wszystkich modeli

porownanie.train <- rbind(
  OLS          = regresja.metryki(real.train, pred.ols.train),
  Ridge        = regresja.metryki(real.train, pred.ridge.train),
  Lasso        = regresja.metryki(real.train, pred.lasso.train),
  CART         = regresja.metryki(real.train, pred.tree.train),
  Bagging      = regresja.metryki(real.train, pred.bag.train),
  RandomForest = regresja.metryki(real.train, pred.rf.train),
  RF_tuned     = regresja.metryki(real.train, pred.rf.tuned.train),
  GBM          = regresja.metryki(real.train, pred.gbm.train),
  NeuralNet    = regresja.metryki(real.train, pred.nn.train)
)
knitr::kable(porownanie.train, caption = "Metryki — zbiór treningowy")
Metryki — zbiór treningowy
MSE RMSE MAE MedAE MAPE MSLE R2
OLS 918758.37 958.52 462.40 183.03 0.1148 0.021349 0.9427
Ridge 1595633.46 1263.18 647.54 210.49 0.1484 0.034851 0.9343
Lasso 931775.76 965.29 466.70 181.66 0.1149 0.021396 0.9419
CART 675824.02 822.09 447.62 180.84 0.1159 0.021548 0.9579
Bagging 68941.45 262.57 138.13 66.90 0.0424 0.003433 0.9957
RandomForest 180840.86 425.25 212.08 86.74 0.0569 0.005518 0.9911
RF_tuned 163682.77 404.58 210.48 90.02 0.0592 0.006056 0.9899
GBM 260823.88 510.71 268.27 111.30 0.0732 0.009130 0.9837
NeuralNet 346056.61 588.27 315.58 136.23 0.0857 0.011933 0.9784 
porownanie.test <- rbind(
  OLS          = regresja.metryki(real.test, pred.ols.test),
  Ridge        = regresja.metryki(real.test, pred.ridge.test),
  Lasso        = regresja.metryki(real.test, pred.lasso.test),
  CART         = regresja.metryki(real.test, pred.tree.test),
  Bagging      = regresja.metryki(real.test, pred.bag.test),
  RandomForest = regresja.metryki(real.test, pred.rf.test),
  RF_tuned     = regresja.metryki(real.test, pred.rf.tuned.test),
  GBM          = regresja.metryki(real.test, pred.gbm.test),
  NeuralNet    = regresja.metryki(real.test, pred.nn.test)
)
knitr::kable(porownanie.test, caption = "Metryki — zbiór testowy")
Metryki — zbiór testowy
MSE RMSE MAE MedAE MAPE MSLE R2
OLS 905937.2 951.81 459.21 186.64 0.1138 0.020809 0.9426
Ridge 1561708.0 1249.68 642.96 213.06 0.1471 0.034283 0.9341
Lasso 917630.5 957.93 463.41 184.53 0.1139 0.020842 0.9418
CART 670064.0 818.57 449.15 186.21 0.1163 0.021759 0.9573
Bagging 289295.6 537.86 280.34 125.22 0.0795 0.011319 0.9816
RandomForest 394066.2 627.75 313.94 120.56 0.0810 0.011071 0.9784
RF_tuned 280026.8 529.18 276.88 116.32 0.0755 0.009893 0.9822
GBM 275478.0 524.86 276.88 114.76 0.0743 0.009567 0.9825
NeuralNet 372305.5 610.17 320.81 137.99 0.0857 0.012016 0.9763 
rmse.test <- porownanie.test$RMSE
names(rmse.test) <- rownames(porownanie.test)
par(mar = c(9, 4, 3, 1))
barplot(sort(rmse.test),
        main = "RMSE na zbiorze testowym (niżej = lepiej)",
        ylab = "RMSE (USD)", las = 2, col = "steelblue",
        border = NA, cex.names = 0.8)

par(mar = c(5, 4, 4, 2))
cat("Różnica RMSE train - test (ujemna = model lepszy na teście — typowe dla OOB):\n")
## Różnica RMSE train - test (ujemna = model lepszy na teście — typowe dla OOB):
print(round(porownanie.train$RMSE - porownanie.test$RMSE, 2))
## [1]    6.71   13.50    7.36    3.52 -275.29 -202.50 -124.60  -14.15  -21.90
best.idx <- which.min(porownanie.test$RMSE)
cat("Najlepszy model:", rownames(porownanie.test)[best.idx], "\n")
## Najlepszy model: GBM
cat("RMSE:", porownanie.test$RMSE[best.idx], "USD\n")
## RMSE: 524.86 USD
cat("RMSE jako % średniej ceny:",
    round(porownanie.test$RMSE[best.idx] / mean(real.test) * 100, 1), "%\n")
## RMSE jako % średniej ceny: 13.3 %
par(mfrow = c(1, 2))
wykres.pred(real.test, pred.rf.tuned.test, "RF tuned: predykcja vs rzeczywistość")
hist(real.test - pred.rf.tuned.test, breaks = 60, col = "steelblue",
     border = "white", main = "RF tuned: rozkład błędów",
     xlab = "Błąd predykcji (USD)")
abline(v = 0, col = "red", lwd = 2)

par(mfrow = c(1, 1))

13. Wnioski i interpretacja wyników

Ranking modeli i uzasadnienie:

Najlepsze wyniki osiągnęły RF tuned i GBM (RMSE ok. 534 USD, R² = 0.982). Obie metody budują wiele drzew i redukują wariancję, GBM sekwencyjnie naprawia błędy poprzednich drzew, RF uśrednia niezależnie zbudowane drzewa. Bagging jest nieznacznie gorszy, bo przy mtry = p drzewa są bardziej ze sobą skorelowane. Sieć neuronowa osiąga przyzwoity wynik (RMSE ok. 600), ale nie dorównuje metodom drzewiastym na danych tabelarycznych.

Dlaczego Ridge pogorszyło wynik względem OLS?

Ridge regularyzuje współczynniki, “skracając” je w kierunku zera. Gdy model jest już dobrze skalibrowany przez transformację log(carat), regularyzacja wprowadza bias zamiast zmniejszać wariancję. To klasyczna sytuacja: regularyzacja nie zawsze pomaga, szczególnie gdy model jest prosty i dobrze dobrany.

Co wyzerował Lasso?

Lasso przy optymalnej lambdzie wyzerowało zmienną table (szerokość rondisty diamentu). Potwierdza to obserwację z macierzy korelacji; table ma słaby związek z ceną i jest mniej istotna niż log.carat, clarity czy color.

Najważniejsza cecha: log(carat)

We wszystkich modelach log.carat zdecydowanie dominuje w ważności zmiennych. Masa diamentu jest głównym wyznacznikiem jego ceny rynkowej. Transformacja log(carat) była kluczową decyzją, bez niej OLS osiągał R² poniżej 0.20.

Wnioski z analizy reszt OLS:

Model liniowy ma tendencję do zaniżania ceny najdroższych diamentów. Reszty są prawostronnie skośne i rosnące wraz z predykcją. Modele nieparametryczne (RF, GBM) osiągają bardziej symetryczny rozkład błędów, co potwierdza że nieliniowości przy wysokich cenach są przez nie lepiej uchwytywane.

Kontekst biznesowy:

Model RF tuned popełnia błąd ok. 534 USD przy średniej cenie 3948 USD (ok. 13.5% wartości). Jest wystarczający do wykrycia wyraźnie zawyżonych ofert. Przy diamentach powyżej 10 000 USD absolutny błąd predykcji rośnie — tam nadal potrzebna jest ocena eksperta.

14. Zapis modeli

saveRDS(diamonds.ols,       "diamonds.ols.rds")
saveRDS(diamonds.ridge,     "diamonds.ridge.rds")
saveRDS(diamonds.lasso,     "diamonds.lasso.rds")
saveRDS(diamonds.tree,      "diamonds.tree.rds")
saveRDS(diamonds.bag,       "diamonds.bag.rds")
saveRDS(diamonds.rf,        "diamonds.rf.rds")
saveRDS(diamonds.rf.tuned,  "diamonds.rf.tuned.rds")
saveRDS(diamonds.gbm.final, "diamonds.gbm.final.rds")
saveRDS(diamonds.nn,        "diamonds.nn.rds")
saveRDS(list(mean = mn.train, sd = sd.train), "scaler.nn.rds")