---
title: "Regressione II (notebook 7)"
author: "Dario Trevisan"
date: "5/11/2025"
output:
  html_notebook:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: readable
    df_print: paged
    download_handler: true
  html_document:
    toc: true
    toc_depth: '3'
    df_print: paged
subtitle: "Statistica II - 750AA"
---




```{r setup, include=FALSE}
#knitr::opts_chunk$set(warning = FALSE, message = FALSE, echo=FALSE) 
```


# 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.

```{r}
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.

```{r synthetic-data}
#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

```{r}
# Scatterplot matrix
pairs(df_syn, main = "Scatterplot Matrix dei Dati Sintetici")

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

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

```{r}
library(GGally)
ggpairs(df_syn)

```

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

```{r}

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

summary(fit_syn)

```

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

```{r coef-comparison}
print("Coefficienti stimati:")
print(coef(fit_syn))
print("Coefficienti veri:")
print(c(0.5, beta))  # includiamo l'intercetta

```

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

```{r synthetic-plot}
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.

```{r}
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à.

```{r}

# 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.

```{r}

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:")
print(as.vector(beta_hat_manual))
print("Coefficienti stimati da lm():")
print(coef(fit_syn))
```

## 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.

```{r mse-rse}

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)))
print(paste("RSE calcolato:", round(rse, 4)))
print(paste("RSE da summary():", round(rse_summary, 4)))
```


## 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.


```{r ci-classic}
coefs <- summary(fit_syn)$coefficients
coefs[, "Std. Error"]  # errori standard

confint(fit_syn, level = 0.95)
```

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.

```{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)))
# 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:")
print(round(ci_manual, 4))
print("Intervalli di confidenza da confint():")
print(round(confint(fit_syn), 4))

```

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

```{r}

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

print(paste("Quantile normale per 95% CI:", round(t_quantile, 4)))
# 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):")
print(round(ci_manual_gaussian, 4))
print("Intervalli di confidenza da confint():")
print(round(confint(fit_syn), 4))
```

**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.

```{r boston-lm}

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)
```
Plottiamo i residui vs valori stimati.

```{r}

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à.

```{r}

# 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).

```{r}
# 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).

```{r boston-ci}
confint(fit_boston, level = 0.95)
```

Aggiungiamo un plot degli intervalli di confidenza per i coefficienti.
```{r}
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.

```{r}

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.

```{r}
summary(fit_boston)
```
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"``.

```{r boston-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)
head(pred_pred)

# 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")


```

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.


```{r bootstrap}

# 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
```

Plottiamo gli intervalli di confidenza bootstrap per i coefficienti.

```{r}

# 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.

```{r}

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)")
print(round(ci_t_boston, 4))
print("CI bootstrap percentile (Boston Housing)")
print(round(ci_boot_boston, 4))
```
Confrontiamoli anche graficamente.

```{r}

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.
```{r}

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.

```{r bootstrap-prediction-ci}
# 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.

```{r compare-prediction-ci}

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.

```{r}

# 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.

```{r}

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).

```{r}

# 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).


