Regressione I (notebook 6)

Statistica II - 750AA

Introduzione

In questo notebook esploreremo diverse tecniche di regressione: \(k\)-NN per la regressione, regressione lineare semplice, regressione polinomiale e un esempio di regressione (apparentemente) non lineare. Utilizzeremo vari dataset per illustrare i concetti.

Carichiamo delle librerie che saranno utili.

library(caret) 
library(FNN) # pacchetto con knn per regressione (anche per classificazione)
library(ggplot2)
library(MASS)     # Per dataset reali

Dataset

Come al solito consideriamo tre datasets, uno generato e due reali (uno precaricato).

Dataset generato

Nel primo esempio generiamo un data frame con due colonne, legate da una relazione di associazione lineare cui aggiungiamo del rumore.

set.seed(42)

# scegliamo coefficiente angolare e intercetta
a = 3
b = 2
n <- 100
x <- rt(n, df=6)
y <- sin(x) + 0.1*rt(n, df=4)

df_generated <- data.frame(x, y)

# Visualizza il dataset generato
ggplot(df_generated, aes(x = x, y = y)) +
  geom_point() +
  ggtitle("Dataset Generato")

Dataset precaricato (Boston Housing)

Il dataset Boston dal pacchetto MASS riguarda i prezzi delle case in 506 quartieri di Boston e include varie caratteristiche delle abitazioni.

df_boston <- Boston

# Visualizza il dataset Boston
head(df_boston)

Selezioniamo solo alcune colonne per la visualizzazione (esplorate anche le altre features).

# Visualizza il dataset
ggplot(Boston, aes(x = lstat, y = medv, colour=rm)) +
  geom_point() +
  ggtitle("Boston Housing Dataset")

NA
NA

È già evidente una correlazione negativa tra lstat (stato della popolazione nel quartiere) e medv (valore mediano delle case occupate da proprietari).

Dataset Reale (Ricoveri Careggi 2012)

Usiamo i dati reali raccolti nel 2012 circa i ricoveri presso l’ospedale Careggi di Firenze scaricabili dal sito OpenToscano.

df_careggi <- read.csv("datasets/ricoverati-ammessi-dimessi-e-permanenza-in-reparti-ordinari-per-area-di-attivita.csv")

head(df_careggi)

#usiamo ggpairs per un grafico che combina nuvola di punti, densità stimata e correlogramma (in alternativa usare solo il comando plot)

library(GGally)
ggpairs(df_careggi[2:4])

NA
NA

Vediamo una perfetta correlazione tra Ammessi e Dimessi, mentre una correlazione positiva tra Ammessi (o Dimessi) e giorni accumulati di degenza presso il reparto.

KNN per la Regressione

Per la regressione usiamo inizialmente il modello \(k\)-NN che si basa sulla media dei fattori di uscita dei \(k\) punti più vicini a un punto test in termini di distanza euclidea delle caratteristiche di ingresso. Possiamo usare comandi da due pacchetti: knn.reg() da FNN (che funziona grossomodo come knn del pacchetto class) oppure train specificando il metodo knn e predict di caret.

Dataset generato

library(FNN)

# scegliamo intanto un valore di k a piacere

k <- 5

knn_generated <- knn.reg(train = as.matrix(df_generated$x), test = as.matrix(x), y = df_generated$y, k = k)

df_generated_pred <- data.frame(x = x, y = knn_generated$pred)

ggplot() +
  geom_point(data = df_generated, aes(x = x, y = y)) +
  geom_line(data = df_generated_pred, aes(x = x, y = y), colour='red')+
  ggtitle(paste("Regressione KNN sul dataset generato con k =", k))

Per l’allenamento del modello (ossia la scelta del parametro \(k\)) conviene usare il pacchetto caret.

# Suddivisione dei dati
train_index <- createDataPartition(df_generated$y, p = 0.7, list = FALSE)
train_set <- df_generated[train_index, ]
test_set <- df_generated[-train_index, ]

# Allenamento del modello KNN
knn_model <- train(y ~ x, data = train_set, method = "knn", tuneGrid = data.frame("k"= 1:50))
# per specificare i valori di k dare l'opzione tuneGrid=data.frame("k"=...)

plot(knn_model)

print(knn_model)

k-Nearest Neighbors 

72 samples
 1 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 72, 72, 72, 72, 72, 72, ... 
Resampling results across tuning parameters:

  k   RMSE       Rsquared   MAE      
   1  0.2285392  0.8926993  0.1643142
   2  0.2186904  0.9005930  0.1597419
   3  0.2123933  0.9054767  0.1577372
   4  0.2076258  0.9098265  0.1553781
   5  0.1997857  0.9154713  0.1489045
   6  0.1995335  0.9147411  0.1473896
   7  0.1978503  0.9168041  0.1439567
   8  0.1983277  0.9161153  0.1443514
   9  0.2000995  0.9146540  0.1459419
  10  0.2010859  0.9141720  0.1477840
  11  0.2007175  0.9148232  0.1481636
  12  0.2025846  0.9150106  0.1506338
  13  0.2054951  0.9130907  0.1529180
  14  0.2070120  0.9132680  0.1560020
  15  0.2100132  0.9130388  0.1588261
  16  0.2116558  0.9133411  0.1606263
  17  0.2135414  0.9150390  0.1623269
  18  0.2169669  0.9150046  0.1662987
  19  0.2201531  0.9160524  0.1692551
  20  0.2232591  0.9159149  0.1728527
  21  0.2297069  0.9132365  0.1789611
  22  0.2331192  0.9136055  0.1824157
  23  0.2395994  0.9116681  0.1883441
  24  0.2440008  0.9105096  0.1932502
  25  0.2502260  0.9107947  0.2004753
  26  0.2580332  0.9102989  0.2081613
  27  0.2641159  0.9099707  0.2141989
  28  0.2718084  0.9094414  0.2215129
  29  0.2797407  0.9083466  0.2297738
  30  0.2876010  0.9066373  0.2371173
  31  0.2943189  0.9061321  0.2430919
  32  0.3024447  0.9051786  0.2503504
  33  0.3078908  0.9061798  0.2559810
  34  0.3150735  0.9066025  0.2628268
  35  0.3227413  0.9071477  0.2692529
  36  0.3317649  0.9083220  0.2775795
  37  0.3384199  0.9084589  0.2835297
  38  0.3487506  0.9057239  0.2924154
  39  0.3583677  0.9064073  0.3009588
  40  0.3660420  0.9067228  0.3087697
  41  0.3741855  0.9074475  0.3165264
  42  0.3825200  0.9076738  0.3245892
  43  0.3919679  0.9073862  0.3330232
  44  0.4014234  0.9081728  0.3421511
  45  0.4107612  0.9095466  0.3504502
  46  0.4187719  0.9135308  0.3571954
  47  0.4298985  0.9119670  0.3672567
  48  0.4388752  0.9107123  0.3740874
  49  0.4482458  0.9121328  0.3815309
  50  0.4581600  0.9116745  0.3899834

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was k = 7.

Per la previsione usiamo la funzione predict specificando modello e test set.

# Previsione sul test set
knn_pred <- predict(knn_model, test_set)

# Calcolo dell'errore
knn_rmse <- sqrt(mean((knn_pred - test_set$y)^2))
print(paste("RMSE KNN:", knn_rmse))

[1] "RMSE KNN: 0.417918495369135"

# plot della previsione vs fattori di uscita effettivi

df_test_pred <- data.frame(x = test_set$x, y_true = test_set$y, y_pred = knn_pred)
ggplot(df_test_pred, aes(x = x)) +
  geom_point(aes(y = y_true), colour = "black", shape=1) +
  geom_point(aes(y = y_pred), colour = "red", shape=2) +
  ggtitle("Previsioni KNN vs Valori Reali sul Test Set") +
  ylab("Valore di uscita") +
  xlab("Caratteristica di ingresso")

Possiamo anche considerare altri indicatori di performance sul test set (confrontiamoli con quelli della cross validation usata per determinare \(k\)).

# calcola il MAE sul test set

knn_mae <- mean(abs(knn_pred - test_set$y))

# calcola il coefficiente di determinazione R^2 sul test set

ss_total <- sum((test_set$y - mean(test_set$y))^2)
ss_res <- sum((knn_pred - test_set$y)^2)
knn_r2 <- 1 - (ss_res / ss_total)

L`erore medio assoluto (MAE) vale 0.1923258 e il coefficiente di determinazione \(R^2\) vale 0.5849585 (calcolati sul test set).

Dataset Boston

Passiamo ora al dataset sulle case dei quartieri di Boston. Esso contiene molte caratteristiche, quindi possiamo valutare di effettuare una riduzione dimensionale (ad esempio tramite PCA). Ricordate anche che con i metodi basati sulla distanza è consigliato standardizzare le caratteristiche di input (non è necessario standardizzare quelle di output).

Usiamo la caratteristica medv (valore della casa mediano nel quartiere) come caratteristica di uscita.

# usiamo il pacchetto caret per knn sul dataset boston

# Suddivisione dei dati
set.seed(42)
train_index <- createDataPartition(df_boston$medv, p = 0.7, list = FALSE)
train_set <- df_boston[train_index, ]
test_set <- df_boston[-train_index, ]

# Allenamento del modello KNN specifichiamo di centrare e scalare i dati
knn_model_boston <- train(medv ~ ., data = train_set,   preProcess = c("center", "scale"), method = "knn", tuneGrid = data.frame("k"= 1:20))

plot(knn_model_boston)

print(knn_model_boston)

k-Nearest Neighbors 

356 samples
 13 predictor

Pre-processing: centered (13), scaled (13) 
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 356, 356, 356, 356, 356, 356, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   1  5.286680  0.6936377  3.432798
   2  5.203928  0.6912242  3.388151
   3  5.119846  0.7014548  3.353115
   4  5.169694  0.6963904  3.352357
   5  5.210683  0.6921727  3.361788
   6  5.213923  0.6941756  3.358467
   7  5.182495  0.6997888  3.332981
   8  5.140886  0.7095845  3.320832
   9  5.195486  0.7045500  3.357256
  10  5.157155  0.7115381  3.346632
  11  5.172255  0.7122806  3.377432
  12  5.150771  0.7171045  3.375410
  13  5.174643  0.7159958  3.412493
  14  5.190267  0.7162546  3.424034
  15  5.227784  0.7129960  3.455417
  16  5.272501  0.7083041  3.490901
  17  5.312728  0.7042011  3.514802
  18  5.334650  0.7025041  3.534245
  19  5.373924  0.6991709  3.561547
  20  5.403465  0.6977176  3.572611

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was k = 3.

Possiamo ora usare il modello per prevedere il valore sui quartieri di test.

# Previsione sul test set
knn_pred_boston <- predict(knn_model_boston, test_set)
# Calcolo dell'errore di test
knn_rmse_boston <- sqrt(mean((knn_pred_boston - test_set$medv)^2))
print(paste("RMSE KNN Boston:", knn_rmse_boston))

[1] "RMSE KNN Boston: 3.8008312930882"

# Coefficiente R^2 di test
ss_total_boston <- sum((test_set$medv - mean(test_set$medv))^2)
ss_res_boston <- sum((knn_pred_boston - test_set$medv)^2)
knn_r2_boston <- 1 - (ss_res_boston / ss_total_boston)
print(paste("R^2 KNN Boston:", knn_r2_boston))

[1] "R^2 KNN Boston: 0.815582465902293"

Esercizio: ripetere l’analisi selezionando solamente alcune caratteristiche del dataset. Come cambia il RMSE? come cambia il coefficiente di determinazione \(R^2\)? Come cambia la performance del metodo senza standardizzare i dati?

Confrontiamo l’errore di training con l’errore di test al variare di \(k\).

# calcola train e test error al variare di k
train_error_boston <- numeric(20)
for (k in 1:20){
  # calcola le previsioni sul train set
  knn_pred_train <- predict(train(medv ~ ., data = train_set,   preProcess = c("center", "scale"), method = "knn", tuneGrid = data.frame("k"= k)), train_set)
  # calcola il RMSE
  train_error_boston[k] <- sqrt(mean((knn_pred_train - train_set$medv)^2))
}

# plot dei risultati della cross-validation
df_error_boston <- data.frame(k = 1:20, train_error = train_error_boston, test_error = knn_model_boston$results$RMSE)
ggplot(df_error_boston, aes(x = k)) +
  geom_line(aes(y = train_error, colour = "Training Error")) +
  geom_point(aes(y = train_error, colour = "Training Error"), shape=1) +
  geom_line(aes(y = test_error, colour = "Test Error")) +
  geom_point(aes(y = test_error, colour = "Test Error"), shape=1) +
  ylab("RMSE") +
  ggtitle("Training e Test Error al variare di k (Boston Housing)") +
  scale_colour_manual("", 
                      breaks = c("Training Error", "Test Error"),
                      values = c("Training Error"="violet", "Test Error"="blue"))+
  theme_minimal()

Dataset Careggi

Consideriamo infine il dataset dei ricoveri al Careggi. L’obiettivo della regressione è prevedere il numero di Giorni Accumulati in funzione del numero di Ammessi (o Dimessi).

# Suddivisione dei dati
set.seed(42)
train_index <- createDataPartition(df_careggi$GiorniAccumulati, p = 0.7, list = FALSE)
train_set <- df_careggi[train_index, ]
test_set <- df_careggi[-train_index, ]

# Allenamento del modello KNN
knn_model_careggi <- train(GiorniAccumulati ~ Ammessi, data = train_set, method = "knn", tuneGrid = data.frame("k"= 1:40))
plot(knn_model_careggi)

print(knn_model_careggi)

k-Nearest Neighbors 

51 samples
 1 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 51, 51, 51, 51, 51, 51, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   1  4312.127  0.2244879  3235.359
   2  3852.788  0.2679383  2870.393
   3  3751.619  0.2902639  2754.002
   4  3564.636  0.3179102  2643.992
   5  3495.508  0.3222912  2580.224
   6  3391.071  0.3391636  2520.578
   7  3318.963  0.3483867  2467.923
   8  3325.561  0.3417332  2467.559
   9  3312.029  0.3372959  2461.531
  10  3265.029  0.3526664  2426.125
  11  3226.596  0.3582410  2410.690
  12  3217.623  0.3646161  2412.532
  13  3182.658  0.3751188  2397.457
  14  3203.922  0.3654594  2412.722
  15  3212.363  0.3637422  2412.042
  16  3198.851  0.3657141  2402.494
  17  3187.655  0.3697047  2392.078
  18  3181.200  0.3756738  2392.754
  19  3194.256  0.3726998  2390.529
  20  3194.230  0.3793524  2385.257
  21  3205.937  0.3777413  2396.998
  22  3219.338  0.3754725  2401.059
  23  3238.105  0.3723319  2412.547
  24  3249.707  0.3723834  2414.864
  25  3260.388  0.3659392  2421.295
  26  3281.801  0.3644368  2439.896
  27  3310.523  0.3432603  2450.157
  28  3318.563  0.3459337  2449.805
  29  3327.317  0.3406555  2456.821
  30  3337.274  0.3414789  2469.056
  31  3351.196  0.3378359  2475.937
  32  3373.304  0.3348179  2498.442
  33  3375.098  0.3378106  2495.380
  34  3389.458  0.3301105  2512.359
  35  3410.479  0.3210328  2525.236
  36  3436.087  0.3091044  2544.050
  37  3468.837  0.2984539  2578.627
  38  3499.223  0.2832781  2609.325
  39  3523.468  0.2810815  2646.532
  40  3542.763  0.2869567  2671.431

RMSE was used to select the optimal model using
 the smallest value.
The final value used for the model was k = 18.

Il modello con \(k=\) rknn_model_careggi$bestTune$k minimizza l’errore di cross-validation. Tuttavia il coefficiente \(R^2=\) 0.3756738 indica che il modello spiega solo una piccola parte della varianza.

Regressione Lineare

Passiamo ora alla regressione lineare semplice (una variabile di input). Il comando di base è lm(), ma anche caret può essere usato per l’allenamento del modello.

Dataset generato

Per il dataset generato possiamo aspettarci un’ottima aderenza al modello sottostante.

# divisione train e test set
train_index <- createDataPartition(df_generated$y, p = 0.7, list = FALSE)
train_set <- df_generated[train_index, ]
test_set <- df_generated[-train_index, ]

# Allenamento del modello di regressione lineare semplice
lm_model_generated <- lm(y ~ ., data = train_set)

summary(lm_model_generated)

Call:
lm(formula = y ~ ., data = train_set)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.86328 -0.24552 -0.04007  0.33662  1.12190 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.04895    0.05818  -0.841    0.403    
x            0.41069    0.05225   7.860 3.31e-11 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4931 on 70 degrees of freedom
Multiple R-squared:  0.4688,    Adjusted R-squared:  0.4612 
F-statistic: 61.78 on 1 and 70 DF,  p-value: 3.307e-11

Dalla teoria sappiamo che il coefficiente angolare è proporzionale alla correlazione. Verifichiamolo numericamente.

lambda <- cor(train_set$y, train_set$x)*sd(train_set$y)/sd(train_set$x)


print(paste("Coefficiente angolare stimato:", coef(lm_model_generated)[2]))

[1] "Coefficiente angolare stimato: 0.41068744537476"

print(paste("Coefficiente angolare teorico:", lambda))

[1] "Coefficiente angolare teorico: 0.410687445374759"

Possiamo effettuare ora una previsione sul test set e calcolare l’errore (RMSE).

# Previsione sul test set
lm_pred_generated <- predict(lm_model_generated, test_set)


# Calcolo dell'errore
lm_rmse_generated <- sqrt(mean((lm_pred_generated - test_set$y)^2))
print(paste("RMSE Regressione Lineare (generato):", lm_rmse_generated))

[1] "RMSE Regressione Lineare (generato): 0.309452490992269"

Visualizziamo anche la retta di regressione.

ggplot(df_generated, aes(x = x, y = y)) +
  geom_point() +
geom_abline(slope=lm_model_generated$coefficients[2], intercept=lm_model_generated$coefficients[1], col='red') +  ggtitle("Regressione Lineare Semplice sul Dataset Generato")

Come possiamo aggiungere dell’informazione riguardo all’incertezza della previsione? Vedremo la teoria nella prossima lezione, per il momento usiamo il comando predict con l’opzione interval="prediction".

# Intervalli di previsione
lm_pred_conf <- predict(lm_model_generated, test_set, interval="prediction", level=0.95)

# Aggiungi delle barre di previsione al plot
df_test_pred_conf <- data.frame(x = test_set$x, y_true = test_set$y, y_pred = lm_pred_conf[,1], lwr = lm_pred_conf[,2], upr = lm_pred_conf[,3])
ggplot(df_test_pred_conf, aes(x = x)) +
  geom_point(aes(y = y_true), colour = "blue", shape=1) +
  geom_line(aes(y = y_pred), colour = "red") +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2) +
  ggtitle("Regressione Lineare con Intervalli di previsione") +
  ylab("Valore di uscita") +
  xlab("Caratteristica di ingresso")

Dataset Boston housing

Consideriamo il dataset Boston Housing e la relazione tra lstat (percentuale della popolazione a basso reddito) e medv (valore mediano delle case occupate da proprietari). Usiamo stavolta il pacchetto caret.

# Allenamento del modello di regressione lineare semplice
set.seed(42)
train_index <- createDataPartition(Boston$medv, p = 0.7, list = FALSE)
train_set <- Boston[train_index, ]
test_set <- Boston[-train_index, ]

# usiamo caret
lm_model_boston <- train(medv ~ lstat, data = train_set, method = "lm")

print(lm_model_boston)

Linear Regression 

356 samples
  1 predictor

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 356, 356, 356, 356, 356, 356, ... 
Resampling results:

  RMSE      Rsquared   MAE     
  6.492122  0.5251329  4.840703

Tuning parameter 'intercept' was held constant at
 a value of TRUE

Confrontiamo i parametri trovati con i due metodi.

summary(lm(medv ~ lstat, data = train_set))

Call:
lm(formula = medv ~ lstat, data = train_set)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.293  -4.114  -1.476   2.286  24.351 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.56528    0.69568   49.69   <2e-16 ***
lstat       -0.93555    0.04758  -19.66   <2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.464 on 354 degrees of freedom
Multiple R-squared:  0.522, Adjusted R-squared:  0.5207 
F-statistic: 386.7 on 1 and 354 DF,  p-value: < 2.2e-16

lm_model_boston$finalModel

Call:
lm(formula = .outcome ~ ., data = dat)

Coefficients:
(Intercept)        lstat  
    34.5653      -0.9356  

Notiamo che il coefficiente \(R^2\) è peggiore di quello trovato con \(k\)-NN (questo è dovuto al bias del modello lineare oltre al fatto che stiamo usando una sola caratteristica: provare ad applicare \(k\)-NN con il solo fattore di ingresso lstat). Confrontiamo quindi l’errore di train con l’errore di test.

# calcolar RMSE e R**2 sul test set
train_pred_boston <- predict(lm_model_boston, train_set)
train_rmse_boston <- sqrt(mean((train_pred_boston - train_set$medv)^2))
ss_total_train_boston <- sum((train_set$medv - mean(train_set$medv))^2)
ss_res_train_boston <- sum((train_pred_boston - train_set$medv)^2)
train_r2_boston <- 1 - (ss_res_train_boston / ss_total_train_boston)
test_pred_boston <- predict(lm_model_boston, test_set)
test_rmse_boston <- sqrt(mean((test_pred_boston - test_set$medv)^2))
ss_total_test_boston <- sum((test_set$medv - mean(test_set$medv))^2)
ss_res_test_boston <- sum((test_pred_boston - test_set$medv)^2)
test_r2_boston <- 1 - (ss_res_test_boston / ss_total_test_boston)
print(paste("RMSE Train Boston:", train_rmse_boston))

[1] "RMSE Train Boston: 6.44581268369236"

print(paste("R^2 Train Boston:", train_r2_boston))

[1] "R^2 Train Boston: 0.522041789419562"

print(paste("RMSE Test Boston:", test_rmse_boston))

[1] "RMSE Test Boston: 5.60103528132736"

print(paste("R^2 Test Boston:", test_r2_boston))

[1] "R^2 Test Boston: 0.599519254079983"

Procediamo comunque con la previsione sul test set e plottiamo la retta di regressione.

# Previsione sul test set
lm_pred <- predict(lm_model_boston, test_set)

# Calcolo dell'errore
lm_rmse <- sqrt(mean((lm_pred - test_set$medv)^2))
print(paste("RMSE Regressione Lineare:", lm_rmse))

[1] "RMSE Regressione Lineare: 5.60103528132736"

# Visualizza la regressione
ggplot(Boston, aes(x = lstat, y = medv)) +
  geom_point() +
  geom_abline(slope=lm_model_boston$finalModel$coefficients[2], intercept = lm_model_boston$finalModel$coefficients[1], col="blue")+
  ggtitle("Regressione Lineare Semplice")

Aggiungiamo l’intervallo di previsione (confrontare con intervallo di confidenza).

lm_model_boston = lm(medv ~ lstat, data = train_set)

# Intervalli di previsione
lm_pred_conf_boston <- predict(lm_model_boston, test_set, interval="prediction", level=0.95)
# Aggiungi delle barre di previsione al plot
df_test_pred_conf_boston <- data.frame(x = test_set$lstat, y_true = test_set$medv, y_pred = lm_pred_conf_boston[,1], lwr = lm_pred_conf_boston[,2], upr = lm_pred_conf_boston[,3])
ggplot(df_test_pred_conf_boston, aes(x = x)) +
  geom_point(aes(y = y_true), colour = "blue", shape=1) +
  geom_line(aes(y = y_pred), colour = "red") +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2) +
  ggtitle("Regressione Lineare Boston con Intervalli di previsione") +
  ylab("Valore di uscita") +
  xlab("Caratteristica di ingresso")

Regressione polinomiale

Proviamo ad aumentare la complessità del modello mantenendo una singola caratteristica (lstat) ma usando una funzione polinomiale. Consideriamo per confronto il dataset generato (per cui non dovrebbe essere necessiario).

# Allenamento del modello di regressione polinomiale con grado 2
poly_model_generated <- lm(y ~ x + I(x^2), data = df_generated)
summary(poly_model_generated)

Call:
lm(formula = y ~ x + I(x^2), data = df_generated)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.78027 -0.26017 -0.03059  0.26429  1.80353 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.06442    0.04705   1.369    0.174    
x            0.46085    0.03653  12.617  < 2e-16 ***
I(x^2)      -0.07120    0.01689  -4.214 5.62e-05 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4134 on 97 degrees of freedom
Multiple R-squared:  0.6261,    Adjusted R-squared:  0.6184 
F-statistic: 81.21 on 2 and 97 DF,  p-value: < 2.2e-16

Una alternativa tecnicamente migliore è di usare il comando poly() specificando il grado del modello polinomiale. Attenzione però che i coefficienti non sono direttamente quelli del polinomio naturale (il comando usa combinazioni lineari dei polinomi per migliorare la precisione).

poly_model_generated <- lm(y ~ poly(x, degree =10), data = df_generated)
summary(poly_model_generated)

Call:
lm(formula = y ~ poly(x, degree = 10), data = df_generated)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.32774 -0.08232 -0.01211  0.07392  0.54405 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -0.02860    0.01366  -2.095   0.0391 *  
poly(x, degree = 10)1   4.97241    0.13656  36.411  < 2e-16 ***
poly(x, degree = 10)2  -1.74226    0.13656 -12.758  < 2e-16 ***
poly(x, degree = 10)3  -3.68979    0.13656 -27.019  < 2e-16 ***
poly(x, degree = 10)4   0.78815    0.13656   5.771 1.13e-07 ***
poly(x, degree = 10)5   0.73409    0.13656   5.375 6.08e-07 ***
poly(x, degree = 10)6  -0.06162    0.13656  -0.451   0.6529    
poly(x, degree = 10)7  -0.02097    0.13656  -0.154   0.8783    
poly(x, degree = 10)8   0.21391    0.13656   1.566   0.1208    
poly(x, degree = 10)9  -0.07845    0.13656  -0.574   0.5671    
poly(x, degree = 10)10  0.29613    0.13656   2.168   0.0328 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1366 on 89 degrees of freedom
Multiple R-squared:  0.9626,    Adjusted R-squared:  0.9584 
F-statistic: 228.8 on 10 and 89 DF,  p-value: < 2.2e-16

Vediamo che in entrambi i casi il coefficiente del termine di grado \(0\) non è significativo (quindi può essere rimosso). Vedremo meglio la teoria di questo test la prossima lezione. Per il momento applichiamo la stessa strategia al datset Boston.

# Allenamento del modello di regressione polinomiale con grado 2
poly_model <- lm(medv ~ poly(lstat, degree = 10), data = train_set)

summary(poly_model)

Call:
lm(formula = medv ~ poly(lstat, degree = 10), data = train_set)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.7454  -3.2492  -0.8553   2.1386  26.1995 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  22.6596     0.2851  79.469  < 2e-16 ***
poly(lstat, degree = 10)1  -127.1043     5.3800 -23.625  < 2e-16 ***
poly(lstat, degree = 10)2    58.2666     5.3800  10.830  < 2e-16 ***
poly(lstat, degree = 10)3   -24.0367     5.3800  -4.468 1.07e-05 ***
poly(lstat, degree = 10)4    18.4519     5.3800   3.430 0.000677 ***
poly(lstat, degree = 10)5   -17.1279     5.3800  -3.184 0.001587 ** 
poly(lstat, degree = 10)6     5.1373     5.3800   0.955 0.340296    
poly(lstat, degree = 10)7     4.4559     5.3800   0.828 0.408106    
poly(lstat, degree = 10)8    -4.2425     5.3800  -0.789 0.430908    
poly(lstat, degree = 10)9     8.8995     5.3800   1.654 0.098999 .  
poly(lstat, degree = 10)10   -7.4488     5.3800  -1.385 0.167085    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.38 on 345 degrees of freedom
Multiple R-squared:  0.6773,    Adjusted R-squared:  0.668 
F-statistic: 72.42 on 10 and 345 DF,  p-value: < 2.2e-16

Vediamo che il coefficiente \(R^2\) è migliorato (ma attenzione all’overfitting).

# Previsione sul test set
poly_pred <- predict(poly_model, test_set)

# Calcolo dell'errore
poly_rmse <- sqrt(mean((poly_pred - test_set$medv)^2))
print(paste("RMSE Regressione Polinomiale:", poly_rmse))

[1] "RMSE Regressione Polinomiale: 4.80813588113572"

# Visualizza la regressione polinomiale
ggplot(df_boston, aes(x = lstat, y = medv)) +
  geom_point() +
  geom_line(data=test_set, aes(x = lstat, y = predict(poly_model, newdata = test_set), col="red" ))+
 # geom_smooth(method = "lm", formula = y ~ poly(x, 2), col = "red") +
  ggtitle("Regressione Polinomiale")

Studiamo al variare del grado l’errore di training (usiamo qui \(R^2\) e cerchiamo un punto di gomito).

set.seed(42)

train_r2 <- numeric(10)

for (grado in 1:10){
  # allenamento del modello polinomiale usando lm()
  poly_model_cv <- glm(medv ~ poly(lstat, degree =  grado ), data = train_set)
  
  # calcolo R^2 di train
  train_pred <- predict(poly_model_cv, train_set)
  ss_total_train <- sum((train_set$medv - mean(train_set$medv))^2)
  ss_res_train <- sum((train_pred - train_set$medv)^2)
  train_r2[grado] <- 1 - (ss_res_train / ss_total_train)
}

# plot
df_r2 <- data.frame(grado = 1:10, train_r2 = train_r2)
ggplot(df_r2, aes(x = grado, y = train_r2)) +
  geom_line() +
  geom_point() +
  ylab("R^2 di Train") +
  xlab("Grado del Polinomio") +
  ggtitle("R^2 di Train al variare del grado del polinomio")

Regressione esponenziale

Possiamo usare i modelli lineari anche per dipendenze più complesse, ad esempio possiamo cercare una relazione del tipo \(y \approx exp(ax+b)\): basta infatti passare al logaritmo e diventa una regressione lineare semplice \(\log(y) \approx ax+b\).

# Allenamento del modello di regressione esponenziale

There were 50 or more warnings (use warnings() to see the first 50)

exp_model <- lm(log(medv) ~ lstat, data = train_set)

summary(exp_model)

Call:
lm(formula = log(medv) ~ lstat, data = train_set)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.97429 -0.15022 -0.03139  0.12516  0.87385 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.605393   0.027247  132.32   <2e-16 ***
lstat       -0.044459   0.001863  -23.86   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2532 on 354 degrees of freedom
Multiple R-squared:  0.6166,    Adjusted R-squared:  0.6155 
F-statistic: 569.2 on 1 and 354 DF,  p-value: < 2.2e-16

# Previsione sul test set
exp_pred_log <- predict(exp_model, test_set)
exp_pred <- exp(exp_pred_log)
# Calcolo dell'errore
exp_rmse <- sqrt(mean((exp_pred - test_set$medv)^2))
print(paste("RMSE Regressione Esponenziale:", exp_rmse))

[1] "RMSE Regressione Esponenziale: 5.27754108639585"

# Visualizza la regressione esponenziale
ggplot() +
  geom_point(data = test_set, aes(x = lstat, y = medv), color = "blue") +
  geom_line(data = test_set, aes(x = lstat, y = exp(predict(exp_model, newdata = test_set))), color = "red") +
  labs(title = "Regressione Esponenziale", x = "lstat", y = "medv") +
  theme_minimal()

Dataset Careggi

Confrontiamo infine la performance di una regressione lineare semplice, polinomiale e infine esponenziale sul dataset dell’ospedale Careggi.

# regressione lineare

set.seed(42)
train_index <- createDataPartition(df_careggi$GiorniAccumulati, p = 0.8, list = FALSE)
train_set <- df_careggi[train_index, ]
test_set <- df_careggi[-train_index, ]


linear_model_careggi <- lm(GiorniAccumulati ~ Ammessi, data = train_set)
  
summary(linear_model_careggi)

Call:
lm(formula = GiorniAccumulati ~ Ammessi, data = train_set)

Residuals:
    Min      1Q  Median      3Q     Max 
-5895.3 -1617.4  -703.6  1387.1 13823.8 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1190.8305   566.3918   2.102   0.0399 *  
Ammessi        3.7817     0.3837   9.857 6.36e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3119 on 57 degrees of freedom
Multiple R-squared:  0.6303,    Adjusted R-squared:  0.6238 
F-statistic: 97.16 on 1 and 57 DF,  p-value: 6.359e-14

linear_pred_careggi <- predict(linear_model_careggi, test_set)
linear_rmse_values <- sqrt(mean((poly_pred_careggi - test_set$GiorniAccumulati)^2))

Consideriamo ora una regressione quadratica (polinomiale di grado \(2\)).

poly_model_careggi <- lm(GiorniAccumulati ~ poly(Ammessi, degree = 2), data = train_set)
summary(poly_model_careggi)

Call:
lm(formula = GiorniAccumulati ~ poly(Ammessi, degree = 2), data = train_set)

Residuals:
    Min      1Q  Median      3Q     Max 
-6132.5 -1652.1  -464.4  1283.0 13487.4 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    5083        408  12.459  < 2e-16 ***
poly(Ammessi, degree = 2)1    30744       3134   9.811 9.08e-14 ***
poly(Ammessi, degree = 2)2    -2151       3134  -0.686    0.495    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3134 on 56 degrees of freedom
Multiple R-squared:  0.6333,    Adjusted R-squared:  0.6202 
F-statistic: 48.36 on 2 and 56 DF,  p-value: 6.301e-13

poly_pred_careggi <- predict(poly_model_careggi, test_set)
poly_rmse_values <- sqrt(mean((poly_pred_careggi - test_set$GiorniAccumulati)^2))

Vediamo che il termine di secondo grado non è significativo. Possiamo limitarci al caso lineare. Consideriamo infine una relazione del tipo \(Giorni \approx Visite ^\alpha\), dove \(\alpha\) è da determinare. Per questo possiamo passare al logaritmo ambo i membri e diventa di nuovo una regressione lineare semplice.

log_model_careggi <- lm(log(GiorniAccumulati) ~ log(Ammessi), data = train_set)
summary(log_model_careggi)

Call:
lm(formula = log(GiorniAccumulati) ~ log(Ammessi), data = train_set)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.03909 -0.44110  0.00677  0.48194  1.88021 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.5541     0.5340    2.91  0.00515 ** 
log(Ammessi)   0.9937     0.0822   12.09  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.825 on 57 degrees of freedom
Multiple R-squared:  0.7194,    Adjusted R-squared:  0.7145 
F-statistic: 146.1 on 1 and 57 DF,  p-value: < 2.2e-16

log_pred_careggi_log <- predict(log_model_careggi, test_set)
log_pred_careggi <- exp(log_pred_careggi_log)
log_rmse_values <- sqrt(mean((log_pred_careggi - test_set$GiorniAccumulati)^2))

Il coefficiente del logaritmo è \(0.99\), quindi l’ipotesi di relazione lineare è ulteriormente confermata.

Confrontiamo infine i RMSE.

print(paste("RMSE Regressione Lineare (Careggi):", linear_rmse_values))

[1] "RMSE Regressione Lineare (Careggi): 3756.44220814634"

print(paste("RMSE Regressione Polinomiale (Careggi):", poly_rmse_values)) 

[1] "RMSE Regressione Polinomiale (Careggi): 3756.44220814634"

print(paste("RMSE Regressione Esponenziale (Careggi):", log_rmse_values)) 

[1] "RMSE Regressione Esponenziale (Careggi): 3579.85710029027"

Plottiamo i tre modelli.

ggplot() +
  geom_point(data = df_careggi, aes(x = Ammessi, y = GiorniAccumulati), color = "blue") +
  geom_smooth(data = train_set, aes(x = Ammessi, y = GiorniAccumulati), method = "lm", formula = y ~ x, color = "red", se = FALSE) +
  geom_smooth(data = train_set, aes(x = Ammessi, y = GiorniAccumulati), method = "lm", formula = y ~ poly(x, 2), color = "green", se = FALSE) +
  geom_line(data = test_set, aes(x = Ammessi, y = exp(predict(log_model_careggi, newdata = test_set))), color = "purple") +
  labs(title = "Confronto Modelli di Regressione (Careggi)", x = "Ammessi", y = "Giorni Accumulati") +
  theme_minimal() +
  scale_color_manual(name = "Modello", values = c("Lineare" = "red", "Polinomiale" = "green", "Esponenziale" = "purple"))

Esercizi

1. Utilizza il dataset mtcars per prevedere il consumo di carburante (mpg) in base al peso dell’auto (wt). Valuta le performance del modello \(k\)-NN al variare del valore di \(k\).

2. Utilizza il dataset iris per prevedere la lunghezza del petalo (Petal.Length) in base alla lunghezza del sepalo (Sepal.Length). Calcola e interpreta il coefficiente di determinazione \(R^2\).

3. Genera un dataset con una relazione quadratica (ad esempio, \(y = x^2 + \text{rumore}\)). Applica una regressione polinomiale di secondo grado e confronta il risultato con una regressione lineare semplice.

4. Crea un dataset in cui la variabile dipendente cresce esponenzialmente rispetto a una variabile indipendente. Applica una regressione esponenziale e confronta il RMSE con quello di una regressione lineare.

5. Utilizza il dataset Iris per prevedere la lunghezza del petalo come funzione della lunghezza del sepalo elevata ad un esponente \(\alpha\) da determinare tramite regressione lineare.

6. Utilizza il dataset airquality per prevedere la temperatura (Temp) in base all’ozono (Ozone). Sperimenta con KNN e regressione lineare semplice e confronta gli errori di previsione.

7. Per un modello KNN costruito con un dataset generato, calcola l’errore medio assoluto (MAE) e il root mean squared error (RMSE) e discuti le loro implicazioni.

8. Utilizza il dataset diamonds per prevedere il prezzo (price) in base alle caratteristiche del diamante (carat, cut, color). Prova a usare KNN e regressione polinomiale e discuti come la scelta delle caratteristiche influisce sulle performance.

9. Costruisci un modello di regressione lineare semplice utilizzando il dataset ChickWeight per prevedere il peso (weight) in base all’età (Time). Visualizza i risultati e gli intervalli di previsione.