Regressione II (notebook 7)

Statistica II - 750AA

Introduzione

In questo notebook mostriamo come costruire e interpretare un modello di regressione lineare multipla, stimare l’incertezza dei coefficienti e delle previsioni (bootstrap e parametrici), eseguire test di significatività ed eventualmente applicarli per la riduzione dei parametri di un modello.

Carichiamo delle librerie che saranno utili.

library(ggplot2)
library(caret)    # per la suddivisione train/test
library(boot)     # per il bootstrap
library(MASS)     # per esempi (se necessario)

Dati sintetici

Ricordiamo che partendo da un modello di regressione lineare multipla \[ y = x \beta + \varepsilon,\qquad \varepsilon\sim \mathcal{N}(0,\sigma^2 I) \]

lo stimatore OLS è \[ \beta_{OLS} = ({\bf{x}}^T\bf{x})^{-1}{\bf{x}}^T\bf{y}. \] dove \(\bf{x}\) indica la matrice dei fattori di ingresso (con una colonna di 1 per l’intercetta), \(\bf{y}\) il vettore delle risposte osservate, e \(\beta\) il vettore dei coefficienti da stimare,

Verifichiamo questa formula generando dei dati con tre fattori di ingresso correlati tra loro e un termine di rumore aggiuntivo.

#set.seed(42)
n <-50
x1 <- rpois(n, lambda=2)
x2 <- 0.8*x1 + 0.4 * rnorm(n)   # correlato con x1
x3 <- rt(n, df=4 )
beta <- c(1.5, 2.0, -1.0)          # coefficienti veri (senza intercept)
y <- 0.5 + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + rnorm(n, sd = 0.6)
df_syn <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3)

Visualizziamo i dati con uno scatterplot e la matrice di correlazione

# Scatterplot matrix
pairs(df_syn, main = "Scatterplot Matrix dei Dati Sintetici")

# Correlation matrix
cor_matrix <- cor(df_syn)
print("Matrice di correlazione:")

[1] "Matrice di correlazione:"

print(round(cor_matrix, 2))

       y   x1   x2    x3
y   1.00 0.90 0.89 -0.31
x1  0.90 1.00 0.89  0.06
x2  0.89 0.89 1.00  0.05
x3 -0.31 0.06 0.05  1.00

Il comando ggpairs (dal pacchetto GGally) può essere usato per una visualizzazione più avanzata.

library(GGally)
ggpairs(df_syn)

Eseguiamo quindi la regressione lineare multipla sui dati generati con il comando lm().

fit_syn <- lm(y ~ ., data = df_syn)

summary(fit_syn)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-1.10479 -0.37649 -0.00144  0.33161  1.18882 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.23102    0.12779   1.808   0.0772 .  
x1           1.85919    0.16939  10.976 1.94e-14 ***
x2           1.63633    0.20344   8.043 2.55e-10 ***
x3          -0.93915    0.05997 -15.661  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5287 on 46 degrees of freedom
Multiple R-squared:  0.9945,    Adjusted R-squared:  0.9941 
F-statistic:  2772 on 3 and 46 DF,  p-value: < 2.2e-16

Confrontiamo i coefficienti stimati con quelli veri impostati da noi. Cosa notiamo?

print("Coefficienti stimati:")

[1] "Coefficienti stimati:"

print(coef(fit_syn))

(Intercept)          x1          x2          x3 
  0.3486526   1.8305423   1.6995324  -1.0507698 

print("Coefficienti veri:")

[1] "Coefficienti veri:"

print(c(0.5, beta))  # includiamo l'intercetta

[1]  0.5  1.5  2.0 -1.0

Trattandosi di una funzione da più variabili, non possiamo visualizzare direttamente la retta di regressione, però possiamo plottare le predizioni vs i valori reali.

df_syn$yhat <- fit_syn$fitted.values #per ottenere i valori predetti dal modello

ggplot(df_syn, aes(x = yhat, y = y)) +
  geom_point(alpha = 0.6) + geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni vs valori reali (dati sintetici)", x="Predetto", y="Osservato")

Un altro modo di visualizzare graficamente la qualità del modello è plottare i residui.

df_syn$residuals <- fit_syn$residuals
ggplot(df_syn, aes(x = yhat, y = residuals)) +
  geom_point(alpha = 0.6) + geom_hline(yintercept=0, linetype=2) +
  labs(title="Residui vs valori stimati", x="Valori stimati", y="Residui")

Anche un istogramma dei residui può essere utile per valutare la normalità.

# Istogramma dei residui

hist(df_syn$residuals, breaks=10, main="Istogramma dei residui", xlab="Residui", col="lightblue", border="black")

# oppure con ggplot2
ggplot(df_syn, aes(x = residuals)) + 
  geom_histogram(aes(y=..density..), bins=10, fill="lightblue", color="black", alpha=0.7) +
  geom_density(color="red", size=1) +
  labs(title="Istogramma dei residui con densità stimata", x="Residui", y="Densità")

Consideriamo ora la matrice di Gram e verifichiamo numericamente che la formula dell’OLS sia rispettata.

Xmat <- as.matrix(cbind( rep(1,5), df_syn[2:4]))

gram_mat <- t(Xmat) %*% Xmat # matrice di Gram

# il comando solve() permette di invertire una matrice
beta_hat_manual <- solve(gram_mat) %*% t(Xmat) %*% df_syn$y

print("Coefficienti stimati manualmente:")

[1] "Coefficienti stimati manualmente:"

print(as.vector(beta_hat_manual))

[1]  0.2310203  1.8591892  1.6363304 -0.9391484

print("Coefficienti stimati da lm():")

[1] "Coefficienti stimati da lm():"

print(coef(fit_syn))

(Intercept)          x1          x2          x3 
  0.2310203   1.8591892   1.6363304  -0.9391484 

Errore quadratico medio e RSE

Ricordiamo che il MSE (per il modello stimato con OLS) è \[ MSE(h(;\beta_{OLS})) = \frac 1 n \sum_{i=1}^n (y_i - x_i \beta_{OLS})^2. \] L’RSE (Residual Standard Error) è invece definito come \[ RSE = \sqrt{\frac{1}{n-p} \sum_{i=1}^n (y_i - x_i \beta_{OLS})^2}, \] dove \(p\) è il numero di parametri stimati (inclusa l’intercetta). Verifichiamo numericamente questi valori. Partiamo direttamente dai residui.

mse <- mean(fit_syn$residuals^2)
rse <- sqrt(sum(fit_syn$residuals^2) / (n - length(coef(fit_syn))))

# la summary di lm fornisce pure l'RSE. confrontiamoli

rse_summary <- summary(fit_syn)$sigma
print(paste("MSE calcolato:", round(mse, 4)))

[1] "MSE calcolato: 0.2571"

print(paste("RSE calcolato:", round(rse, 4)))

[1] "RSE calcolato: 0.5287"

print(paste("RSE da summary():", round(rse_summary, 4)))

[1] "RSE da summary(): 0.5287"

Errori standard e intervalli di confidenza t-based

R fornisce anche gli errori standard delle stime dei coefficienti e gli intervalli di confidenza basati sulla distribuzione t di Student. Il comando summary() mostra gli errori standard, mentre confint() calcola gli intervalli di confidenza.

coefs <- summary(fit_syn)$coefficients
coefs[, "Std. Error"]  # errori standard

(Intercept)          x1          x2          x3 
 0.12779494  0.16939401  0.20343752  0.05996639 

confint(fit_syn, level = 0.95)

                  2.5 %     97.5 %
(Intercept) -0.02621757  0.4882582
x1           1.51821676  2.2001617
x2           1.22683194  2.0458289
x3          -1.05985444 -0.8184423

Il quantile della distribuzione t può essere calcolato con il comando qt() di R. Proviamo a confrontare le formule teoriche con i risultati di R.

alpha <- 0.05
t_quantile <- qt(1 - alpha/2, df = n - length(coef(fit_syn)))

print(paste("Quantile t per 95% CI:", round(t_quantile, 4)))

[1] "Quantile t per 95% CI: 2.0129"

# Calcolo manuale degli intervalli di confidenza
se <- coefs[, "Std. Error"]
ci_manual <- cbind(coef(fit_syn) - t_quantile * se,
                   coef(fit_syn) + t_quantile * se)

colnames(ci_manual) <- c("2.5 %", "97.5 %")
print("Intervalli di confidenza calcolati manualmente:")

[1] "Intervalli di confidenza calcolati manualmente:"

print(round(ci_manual, 4))

              2.5 %  97.5 %
(Intercept) -0.0262  0.4883
x1           1.5182  2.2002
x2           1.2268  2.0458
x3          -1.0599 -0.8184

print("Intervalli di confidenza da confint():")

[1] "Intervalli di confidenza da confint():"

print(round(confint(fit_syn), 4))

              2.5 %  97.5 %
(Intercept) -0.0262  0.4883
x1           1.5182  2.2002
x2           1.2268  2.0458
x3          -1.0599 -0.8184

Sappiamo che per \(n\) grande gli intervalli t-based si avvicinano a quelli basati sulla distribuzione normale.

# confronto con intervalli gaussiani
norm_quantile <- qnorm(1 - alpha/2)

print(paste("Quantile normale per 95% CI:", round(t_quantile, 4)))

[1] "Quantile normale per 95% CI: 2.0129"

# Calcolo manuale degli intervalli di confidenza
se <- coefs[, "Std. Error"]
ci_manual_gaussian <- cbind(coef(fit_syn) - norm_quantile * se,
                   coef(fit_syn) + norm_quantile * se)

colnames(ci_manual_gaussian) <- c("2.5 %", "97.5 %")
print("Intervalli di confidenza calcolati manualmente (normale):")

[1] "Intervalli di confidenza calcolati manualmente (normale):"

print(round(ci_manual_gaussian, 4))

              2.5 %  97.5 %
(Intercept) -0.0195  0.4815
x1           1.5272  2.1912
x2           1.2376  2.0351
x3          -1.0567 -0.8216

print("Intervalli di confidenza da confint():")

[1] "Intervalli di confidenza da confint():"

print(round(confint(fit_syn), 4))

              2.5 %  97.5 %
(Intercept) -0.0262  0.4883
x1           1.5182  2.2002
x2           1.2268  2.0458
x3          -1.0599 -0.8184

Esercizio: verificare che al crescere di \(n\) gli intervalli t-based e quelli normali si avvicinano.

Dati reali: esempio con Boston Housing

Usiamo ora il dataset Boston Housing per mostrare come eseguire una regressione lineare multipla su dati reali.

df_boston <- Boston

# dividiamo in train e test set
set.seed(42)
train_index <- createDataPartition(df_boston$medv, p = 0.8, list = FALSE)
train_set <- df_boston[train_index, ]
test_set <- df_boston[-train_index, ]

# regressione lineare multipla
fit_boston <- lm(medv ~ ., data = train_set)
summary(fit_boston)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6657  -2.6321  -0.5241   1.5665  26.4389 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  33.998764   5.791932   5.870 9.27e-09 ***
crim         -0.122347   0.033662  -3.635 0.000315 ***
zn            0.043296   0.015318   2.826 0.004947 ** 
indus         0.003389   0.073094   0.046 0.963038    
chas          1.554178   0.917885   1.693 0.091207 .  
nox         -15.555219   4.177618  -3.723 0.000225 ***
rm            3.859758   0.486968   7.926 2.36e-14 ***
age           0.004788   0.014694   0.326 0.744717    
dis          -1.305775   0.220032  -5.934 6.48e-09 ***
rad           0.319482   0.075705   4.220 3.04e-05 ***
tax          -0.012468   0.004393  -2.838 0.004770 ** 
ptratio      -0.953531   0.148188  -6.435 3.61e-10 ***
black         0.009516   0.002965   3.209 0.001440 ** 
lstat        -0.516037   0.057726  -8.939  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.706 on 393 degrees of freedom
Multiple R-squared:  0.7382,    Adjusted R-squared:  0.7295 
F-statistic: 85.22 on 13 and 393 DF,  p-value: < 2.2e-16

Plottiamo i residui vs valori stimati.

boston_fit_residuals_fitted <- data.frame( "residuals"= fit_boston$residuals, "yhat" = fit_boston$fitted.values)

ggplot(boston_fit_residuals_fitted, aes(x = yhat, y = residuals)) +
  geom_point(alpha = 0.6) + geom_hline(yintercept=0, linetype=2) +
  labs(title="Residui vs valori stimati (Boston Housing)", x="Valori stimati", y="Residui")

Aggiungiamo un istogramma dei residui per valutarne la gaussianità.

# Istogramma dei residui

hist(fit_boston$residuals, breaks=15, main="Istogramma dei residui (Boston Housing)", xlab="Residui", col="lightblue", border="black")

# usiamo ggplot2 e aggiungiamo una stima della densità
ggplot(boston_fit_residuals_fitted, aes(x = residuals)) + 
  geom_histogram(aes(y=..density..), bins=15, fill="lightblue", color="black", alpha=0.7) +
  geom_density(color="red", size=1) +
  labs(title="Istogramma dei residui con densità stimata (Boston Housing)", x="Residui", y="Densità")

Un altro metodo per valutare la normalità dei residui è il QQ-plot gaussiano. Si tratta di un grafico che confronta i quantili campionari dei residui con i quantili di una distribuzione normale teorica con gli stessi primo e terzo quartile (se i residui sono normali, i punti dovrebbero disporsi lungo una retta).

# per qqplot generali il comando è qqplot
qqnorm(fit_boston$residuals)
qqline(fit_boston$residuals, col = "red")

# possiamo anche usare ggplot2
ggplot(boston_fit_residuals_fitted, aes(sample = residuals)) +
  stat_qq() +
  stat_qq_line(colour = "red") +
  labs(title = "QQ-Plot dei residui (Boston Housing)")

Vediamo che i residui si discostano dalla normalità, specialmente nelle code.

Intervalli di confidenza t-based per i coefficienti

In presenza di non-normalità dei residui, gli intervalli di confidenza t-based potrebbero non essere affidabili. Calcoliamo comunque gli intervalli di confidenza per i coefficienti usando il metodo parametrico (basato sulla t di Student o la gaussiana essendo la taglia molto grande).

confint(fit_boston, level = 0.95)

                    2.5 %       97.5 %
(Intercept)  22.611717156 45.385810568
crim         -0.188527142 -0.056167113
zn            0.013180515  0.073410802
indus        -0.140313846  0.147092644
chas         -0.250401706  3.358758099
nox         -23.768494627 -7.341943227
rm            2.902368695  4.817146589
age          -0.024100269  0.033675953
dis          -1.738363063 -0.873187376
rad           0.170643789  0.468320476
tax          -0.021103675 -0.003831841
ptratio      -1.244872542 -0.662190305
black         0.003686281  0.015344989
lstat        -0.629526620 -0.402547632

Aggiungiamo un plot degli intervalli di confidenza per i coefficienti.

coefs_boston <- summary(fit_boston)$coefficients
se_boston <- coefs_boston[, "Std. Error"]
beta_hat_boston <- coef(fit_boston)
t_quantile_boston <- qt(1 - 0.05/2, df = nrow(train_set) - length(beta_hat_boston))
ci_boston <- cbind(beta_hat_boston - t_quantile_boston * se_boston,
                    beta_hat_boston + t_quantile_boston * se_boston)
ci_df <- data.frame(
  Term = names(beta_hat_boston),
  Estimate = beta_hat_boston,
  Lower = ci_boston[,1],
  Upper = ci_boston[,2]
)

ggplot(ci_df, aes(x = Term, y = Estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2) +
  coord_flip() +
  labs(title = "Intervalli di confidenza per i coefficienti (Boston Housing)", y = "Stima del coefficiente", x = "Termine")

Rimuoviamo l’intercetta, e le variabili, nox, rm, e chas per una visualizzazione più chiara dei coefficienti significativi.

ci_df_filtered <- ci_df[ci_df$Term != "(Intercept)" & ci_df$Term != "nox" & ci_df$Term != "rm" & ci_df$Term != "chas", ]
ggplot(ci_df_filtered, aes(x = Term, y = Estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2) +
  coord_flip() +
  labs(title = "Intervalli di confidenza per i coefficienti (Boston Housing)", y = "Stima del coefficiente", x = "Termine")

Confrontiamo con il risultato di summary al modello.

summary(fit_boston)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6657  -2.6321  -0.5241   1.5665  26.4389 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  33.998764   5.791932   5.870 9.27e-09 ***
crim         -0.122347   0.033662  -3.635 0.000315 ***
zn            0.043296   0.015318   2.826 0.004947 ** 
indus         0.003389   0.073094   0.046 0.963038    
chas          1.554178   0.917885   1.693 0.091207 .  
nox         -15.555219   4.177618  -3.723 0.000225 ***
rm            3.859758   0.486968   7.926 2.36e-14 ***
age           0.004788   0.014694   0.326 0.744717    
dis          -1.305775   0.220032  -5.934 6.48e-09 ***
rad           0.319482   0.075705   4.220 3.04e-05 ***
tax          -0.012468   0.004393  -2.838 0.004770 ** 
ptratio      -0.953531   0.148188  -6.435 3.61e-10 ***
black         0.009516   0.002965   3.209 0.001440 ** 
lstat        -0.516037   0.057726  -8.939  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.706 on 393 degrees of freedom
Multiple R-squared:  0.7382,    Adjusted R-squared:  0.7295 
F-statistic: 85.22 on 13 and 393 DF,  p-value: < 2.2e-16

Vediamo quindi che molte variabili non sono significative (p-value alto): in effetti i loro intervalli di confidenza includono lo zero (questo è il criterio di significatività basato sugli intervalli di confidenza).

Intervalli di previsione

Possiamo calcolare gli intervalli di confidenza per le previsioni su nuovi dati usando il comando predict() con l’opzione interval = "confidence" o interval = "prediction".

# intervalli di confidenza e predizione
pred_conf <- predict(fit_boston, test_set, interval = "confidence", level = 0.99)
pred_pred <- predict(fit_boston, test_set, interval = "prediction", level = 0.99)

# mostriamo i primi risultati
head(pred_conf)

         fit       lwr      upr
1  30.012773 28.191778 31.83377
9  12.121378  8.990245 15.25251
15 19.434768 17.790500 21.07904
19 15.979750 13.853413 18.10609
26 13.552831 11.791572 15.31409
33  8.976062  6.316467 11.63566

head(pred_pred)

         fit        lwr      upr
1  30.012773 17.6960224 42.32952
9  12.121378 -0.4559953 24.69875
15 19.434768  7.1429029 31.72663
19 15.979750  3.6141675 28.34533
26 13.552831  1.2447704 25.86089
33  8.976062 -3.4922888 21.44441

# plot delle predizioni con intervalli
pred_df <- data.frame(
  Actual = test_set$medv,
  Predicted = pred_conf[, "fit"],
  Conf_Lower = pred_conf[, "lwr"],
  Conf_Upper = pred_conf[, "upr"],
  Pred_Lower = pred_pred[, "lwr"],
  Pred_Upper = pred_pred[, "upr"]
)

ggplot(pred_df, aes(x = Predicted, y = Actual)) +
  geom_point(alpha = 0.6) +
  geom_errorbar(aes(ymin = Conf_Lower, ymax = Conf_Upper), width = 0.2, color = "blue", alpha = 0.5, linewidth = 2) +
  geom_errorbar(aes(ymin = Pred_Lower, ymax = Pred_Upper), width = 0.2, color = "red", alpha = 0.5) +
  geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni con intervalli di confidenza (blu) e predizione (rosso)", x="Predetto", y="Osservato")

NA
NA

Notiamo degli intervalli di predizione (rossi) più ampi rispetto a quelli di confidenza (blu), come previsto teoricamente. Tuttavia nonostante la grandezza del dataset, gli intervalli di confidenza non sono molto stretti, indicando una certa incertezza nelle stime dei coefficienti.

Stime bootstrap per i parametri

Il bootstrap è utile quando le ipotesi di normalità e omoschedasticità dei residui \(\varepsilon_i \approx \mathcal{N}(0, \sigma^2)\) non sono così evidenti (come nel caso sopra), o vogliamo comunque una stima empirica degli errori. Lavoriamo sul dataset Boston.

Usiamo la funzione boot() del pacchetto boot per eseguire il bootstrap sui coefficienti del modello di regressione lineare.

# usiamo solamente il train set

head(train_set)

# Funzione che ritorna i coefficienti sui dati ri-campionati
boot_coef_fun <- function(data, indices) {
  d <- data[indices, ]  # campionamento con rimpiazzo
  lm( medv ~., data = d)$coefficients
}

boot_out <- boot(train_set, statistic = boot_coef_fun, R = 2000)
summary(boot_out)  # riassunto del bootstrap

          Length Class      Mode     
t0           14  -none-     numeric  
t         28000  -none-     numeric  
R             1  -none-     numeric  
data         14  data.frame list     
seed        626  -none-     numeric  
statistic     1  -none-     function 
sim           1  -none-     character
call          4  -none-     call     
stype         1  -none-     character
strata      407  -none-     numeric  
weights     407  -none-     numeric  

Plottiamo gli intervalli di confidenza bootstrap per i coefficienti.

# Funzione per plottare gli intervalli di confidenza bootstrap
plot_boot_ci <- function(boot_out, level = 0.95) {
  alpha <- 1 - level
  ci_lower <- apply(boot_out$t, 2, function(tcol) quantile(tcol, probs = alpha/2, na.rm=TRUE))
  ci_upper <- apply(boot_out$t, 2, function(tcol) quantile(tcol, probs = 1 - alpha/2, na.rm=TRUE))
  
  ci_df <- data.frame(
    Term = names(coef(fit_boston)),
    Estimate = coef(fit_boston),
    Lower = ci_lower,
    Upper = ci_upper
  )
  
  ggplot(ci_df, aes(x = Term, y = Estimate)) +
    geom_point() +
    geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2) +
    coord_flip() +
    labs(title = paste("Intervalli di confidenza bootstrap per i coefficienti (level =", level, ")"),
         y = "Stima del coefficiente", x = "Termine")
}

plot_boot_ci(boot_out, level = 0.95)

Confrontiamo gli intervalli di bootstrap con quelli parametrici t-based.

ci_t_boston <- confint(fit_boston)
ci_boot_boston <- t(apply(boot_out$t, 2, function(tcol) quantile(tcol, probs = c(0.025, 0.975))))
print("CI t-based (Boston Housing)")

[1] "CI t-based (Boston Housing)"

print(round(ci_t_boston, 4))

               2.5 %  97.5 %
(Intercept)  22.6117 45.3858
crim         -0.1885 -0.0562
zn            0.0132  0.0734
indus        -0.1403  0.1471
chas         -0.2504  3.3588
nox         -23.7685 -7.3419
rm            2.9024  4.8171
age          -0.0241  0.0337
dis          -1.7384 -0.8732
rad           0.1706  0.4683
tax          -0.0211 -0.0038
ptratio      -1.2449 -0.6622
black         0.0037  0.0153
lstat        -0.6295 -0.4025

print("CI bootstrap percentile (Boston Housing)")

[1] "CI bootstrap percentile (Boston Housing)"

print(round(ci_boot_boston, 4))

          2.5%   97.5%
 [1,]  14.1805 51.8050
 [2,]  -0.1796 -0.0425
 [3,]   0.0133  0.0736
 [4,]  -0.1015  0.1149
 [5,]  -0.5593  3.9884
 [6,] -23.7595 -6.9867
 [7,]   2.0189  6.0812
 [8,]  -0.0299  0.0386
 [9,]  -1.7991 -0.8593
[10,]   0.1638  0.4792
[11,]  -0.0198 -0.0064
[12,]  -1.1883 -0.6993
[13,]   0.0030  0.0158
[14,]  -0.7391 -0.2634

Confrontiamoli anche graficamente.

ci_df_boston <- data.frame(
  Term = names(coef(fit_boston)),
  Estimate = coef(fit_boston),
  T_Lower = ci_t_boston[,1],
  T_Upper = ci_t_boston[,2],
  Boot_Lower = ci_boot_boston[,1],
  Boot_Upper = ci_boot_boston[,2]
)
ggplot(ci_df_boston, aes(x = Term, y = Estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_Lower, ymax = T_Upper), width = 1, color = "blue", alpha = 0.5, linewidth = 1) +
  geom_errorbar(aes(ymin = Boot_Lower, ymax = Boot_Upper), width = 0.2, color = "red", alpha = 1, linewidth =  0.5) +
  coord_flip() +
  labs(title = "Confronto tra intervalli di confidenza t-based (blu) e bootstrap (rosso)",
       y = "Stima del coefficiente", x = "Termine")

Eliminiamo al solito l’intercetta, rm, nox e chas per una visualizzazione più chiara degli intervalli più piccoli.

ci_df_boston_filtered <- ci_df_boston[ci_df_boston$Term != "(Intercept)" & ci_df_boston$Term != "nox" & ci_df_boston$Term != "rm" & ci_df_boston$Term != "chas", ]

ggplot(ci_df_boston_filtered, aes(x = Term, y = Estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = T_Lower, ymax = T_Upper), width = 0.5, color = "blue", alpha = 0.5) +
  geom_errorbar(aes(ymin = Boot_Lower, ymax = Boot_Upper), width = 0.5, color = "red", alpha = 0.5) +
  coord_flip() +
  labs(title = "Confronto tra intervalli di confidenza t-based (blu) e bootstrap (rosso)",
       y = "Stima del coefficiente", x = "Feature")

Vediamo che gli intervalli di confidenza bootstrap e t-based sono abbastanza simili ma non perfettamente identici, indicando che le ipotesi di normalità non sono troppo violate ma neppure perfettamente rispettate (come era evidente anche dal qqplot).

Stime bootstrap per le previsioni

Possiamo anche usare il bootstrap per stimare gli intervalli di confidenza e di previsione per le predizioni su nuovi dati. Usiamo sempre il dataset Boston.

Cominciamo con gli intervalli di confidenza per le predizioni. Qui si tratta semplicemente di calcolare il valore predetto per ogni campione bootstrap e poi calcolare i quantili desiderati.

# Funzione per calcolare le predizioni bootstrap

boot_pred_fun <- function(data, indices) {
  d <- data[indices, ]  # campionamento con rimpiazzo
  model <- lm(medv ~ ., data = d)
  predict(model, test_set)
}

boot_pred_out <- boot(train_set, statistic = boot_pred_fun, R = 2000)

# Calcolo degli intervalli di confidenza bootstrap per le predizioni
pred_ci_lower <- apply(boot_pred_out$t, 2, function(tcol) quantile(tcol, probs = 0.025))
pred_ci_upper <- apply(boot_pred_out$t, 2, function(tcol) quantile(tcol, probs = 0.975))
pred_df$Boot_Conf_Lower <- pred_ci_lower
pred_df$Boot_Conf_Upper <- pred_ci_upper
ggplot(pred_df, aes(x = Predicted, y = Actual)) +
  geom_point(alpha = 0.6) +
  geom_errorbar(aes(ymin = Boot_Conf_Lower, ymax = Boot_Conf_Upper), width = 0.2, color = "blue", alpha = 0.5, linewidth=2) +
  geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni con intervalli di confidenza bootstrap (blu)", x="Predetto", y="Osservato")

Confrontiamo con gli intervalli di confidenza t-based calcolati in precedenza.

ggplot(pred_df, aes(x = Predicted, y = Actual)) +
  geom_point(alpha = 0.6) +
  geom_errorbar(aes(ymin = Conf_Lower, ymax = Conf_Upper), width = 0.2, color = "red", alpha = 0.5, linewidth=2) +
  geom_errorbar(aes(ymin = Boot_Conf_Lower, ymax = Boot_Conf_Upper), width = 0.2, color = "blue", alpha = 0.5, linewidth=2) +
  geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni con intervalli di confidenza t-based (rosso) e bootstrap (blu)", x="Predetto", y="Osservato")

Vediamo che nonostante la mancanza di gaussianità, gli intervalli di confidenza t-based e bootstrap per le predizioni sono dopotutto abbastanza simili.

Passiamo infine agli intervalli di previsione con bootstrap. Qui la procedura è simile, ma dobbiamo aggiungere un termine di rumore simulato per ogni predizione.

# Funzione per calcolare le predizioni bootstrap con rumore
boot_pred_fun_pred <- function(data, indices, newdata) {
  d <- data[indices, ]  # campionamento con rimpiazzo
  model <- lm(medv ~ ., data = d)
  preds <- predict(model, newdata)
  # aggiungiamo rumore ottenuto dai residui sui punti campionati
  residuals <- d$medv - predict(model, d)
  noise <- sample(residuals, size = nrow(newdata), replace = TRUE)
  return(preds + noise)
}

boot_pred_out_pred <- boot(train_set, statistic = function(data, indices) boot_pred_fun_pred(data, indices, test_set), R = 200)

# Calcolo degli intervalli di previsione bootstrap
pred_pred_lower <- apply(boot_pred_out_pred$t, 2, function(tcol) quantile(tcol, probs = 0.025))
pred_pred_upper <- apply(boot_pred_out_pred$t, 2, function(tcol) quantile(tcol, probs = 0.975))
pred_df$Boot_Pred_Lower <- pred_pred_lower
pred_df$Boot_Pred_Upper <- pred_pred_upper
ggplot(pred_df, aes(x = Predicted, y = Actual)) +
  geom_point(alpha = 0.6) +
  geom_errorbar(aes(ymin = Boot_Pred_Lower, ymax = Boot_Pred_Upper), width = 0.2, color = "red", alpha = 0.5, linewidth=2) +
  geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni con intervalli di previsione bootstrap (rosso)", x="Predetto", y="Osservato")

Confrontiamo con gli intervalli di previsione t-based calcolati in precedenza.

ggplot(pred_df, aes(x = Predicted, y = Actual)) +
  geom_point(alpha = 0.6) +
  geom_errorbar(aes(ymin = Pred_Lower, ymax = Pred_Upper), width = 0.2, color = "blue", alpha = 0.5, linewidth=2) +
  geom_errorbar(aes(ymin = Boot_Pred_Lower, ymax = Boot_Pred_Upper), width = 0.2, color = "red", alpha = 0.5, linewidth=2) +
  geom_abline(slope=1, intercept=0, linetype=2) +
  labs(title="Predizioni con intervalli di previsione t-based (blu) e bootstrap (rosso)", x="Predetto", y="Osservato")

Come mai gli intervalli di previsione bootstrap sono più stretti di quelli t-based? Questo può accadere per diversi motivi, ma in particolare notiamo che abbiamo stimato i residui sul training set (ossia i punti selezionati dal bootstrap). Tendenzialmente i residui sul training set sono più piccoli di quelli sui nuovi dati (test set), portando a una sottostima del rumore e quindi a intervalli di previsione troppo stretti. Un modo per ovviare a questo problema è di usare i punti out-of-bag (OOB) per stimare i residui (in realtà un modo ancora migliore è una media pesata tra i residui in-bag e OOB, ma per semplicità usiamo solo gli OOB).

# Funzione per calcolare le predizioni bootstrap con rumore OOB
boot_pred_fun_pred_oob <- function(data, indices, newdata) {
  d <- data[indices, ]  # campionamento con rimpiazzo
  model <- lm(medv ~ ., data = d)
  preds <- predict(model, newdata)
  # calcoliamo i residui OOB
  oob_indices <- setdiff(1:nrow(data), unique(indices))
  if (length(oob_indices) == 0) {
    residuals <- rep(0, nrow(newdata))  # nessun OOB, nessun rumore
  } else {
    oob_data <- data[oob_indices, ]
    oob_preds <- predict(model, oob_data)
    residuals <- oob_data$medv - oob_preds
  }
  noise <- sample(residuals, size = nrow(newdata), replace = TRUE)
  return(preds + noise)
}

boot_pred_out_pred_oob <- boot(train_set, statistic = function(data, indices) boot_pred_fun_pred_oob(data, indices, test_set), R = 200) 
# Calcolo degli intervalli di previsione bootstrap OOB 
pred_pred_oob_lower <- apply(boot_pred_out_pred_oob$t, 2, function(tcol) quantile(tcol, probs = 0.025)) 
pred_pred_oob_upper <- apply(boot_pred_out_pred_oob$t, 2, function(tcol) quantile(tcol, probs = 0.975)) 
pred_df$Boot_Pred_OOB_Lower <- pred_pred_oob_lower 
pred_df$Boot_Pred_OOB_Upper <- pred_pred_oob_upper 

# plot degli intervalli e confronto con i true metodi (bootstrap con residui sul training set, OOB e t-based)
ggplot(pred_df, aes(x = Predicted, y = Actual)) + 
  geom_point(alpha = 0.6) + 
  geom_errorbar(aes(ymin = Boot_Pred_OOB_Lower, ymax = Boot_Pred_OOB_Upper), width = 0.2, color = "green", alpha = 0.5, linewidth=2) + 
  geom_errorbar(aes(ymin = Pred_Lower, ymax = Pred_Upper), width = 0.2, color = "blue", alpha = 0.5, linewidth=0.5) +
  geom_errorbar(aes(ymin = Boot_Pred_Lower, ymax = Boot_Pred_Upper), width = 0.2, color = "red", alpha = 0.5, linewidth=0.5) +
  geom_abline(slope=1, intercept=0, linetype=2) + 
  labs(title="Predizioni con intervalli di previsione bootstrap residui training set (rosso), OOB (verde), t-based (blu)", x="Predetto", y="Osservato")

Esercizi suggeriti

1. Ripetere l’esempio generato variando la correlazione tra \(x_1\) e \(x_2\) e osservare l’effetto sugli errori standard.

2. Confrontare gli intervalli t-based e bootstrap per i coefficienti in iris.

3. Usare il dataset Boston per eseguire una regressione lineare multipla con selezione di alcune caratteristiche e confrontare gli errori di test e gli intervalli di previsione.

4. Esplorare l’effetto di una trasformazione non lineare variabili (\(\log\), \(\sqrt{\cdot}\)) sui residui e sugli intervalli di confidenza (ad esempio nel dataset trees).

5. Implementare una procedura di bootstrap per stimare empiricamente l’errore di confidenza e di previsione per un modello KNN (anche su dati generati).