Intervalli di previsione

Eruzioni a Yellowstone

Nel parco di Yellowstone, nello stato del Wyoming, USA, è situato uno dei geyser più famosi del mondo, chiamato “Old faithful”, le cui eruzioni di acqua bollente hanno un’altezza media di circa 44 metri. La durata delle eruzioni varia dal minuto e mezzo fino ai cinque minuti a intervalli in media di un’ora. Ecco lo scatter con i dati relativi a 272 eruzioni con
in ascissa la durata di un’eruzione e in ordinata il tempo di attesa fino all’eruzione seguente, entrambi in minuti.

data(faithful)
attach(faithful)
plot(eruptions, waiting, pch = 20, col = "RoyalBlue")

L’intervallo tra una eruzione e la seguente dipende dalla durata dell’eruzione: tipicamente l’attesa cresce con la durata dell’eruzione. La distribuzione della durata è bimodale, con due mode a 2 minuti e 4.5 minuti.

Il modello di regressione lineare tra attesa (variabile dipendente) e durata (esplicativa) appare appropriato, e le stime sono le seguenti.

mod = lm(waiting ~ eruptions, data = faithful)
summary(mod)

Call:
lm(formula = waiting ~ eruptions, data = faithful)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.080  -4.483   0.212   3.925  15.972 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   33.474      1.155    29.0   <2e-16 ***
eruptions     10.730      0.315    34.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.91 on 270 degrees of freedom
Multiple R-squared:  0.811, Adjusted R-squared:  0.811 
F-statistic: 1.16e+03 on 1 and 270 DF,  p-value: <2e-16

Pertanto, per ogni minuto in più di durata dobbiamo aspettare l’eruzione seguente circa 11 minuti in più. La variabilità residua non è bassa perché è di circa 6 minuti.

Supponiamo che un turista arrivi alla fine di un’eruzione durata \(x_* = 3.5\) minuti. Quanto si prevede che debba aspettare la prossima?

Previsione ed intervalli di previsione

Per la previsione si usa la funzione predict. Gli argomenti fondamentali sono il risultato di un modello, e un data frame che contiene i dati \(x_0\) per i quali si vogliono fare le previsioni. Per esempio per prevedere il tempo di attesa corrispondente ad una eruzione di 3 minuti, si usano le istruzioni seguenti.

x = data.frame(eruptions = 3)
predict(mod, x)
    1 
65.66

Notate che il data frame x è definito da una riga e una colonna con nome eruptions, cioè esattamente lo stesso della variabile esplicativa. Potremmo costruire le previsioni per più di un valore specificando più righe. Per esempio,

x2 = data.frame(eruptions = c(3,4))
x2
  eruptions
1         3
2         4
predict(mod, x2)
    1     2 
65.66 76.39

Infine per calcolare l’intervallo di previsione e l’errore standard di previsione si usa lo stesso comando, specificando alcuni argomenti aggiuntivi:

p = predict(mod, x, interval = "prediction", se.fit = TRUE, level = 0.95)
p
$fit
    fit   lwr   upr
1 65.66 53.99 77.33

$se.fit
[1] 0.3901

$df
[1] 270

$residual.scale
[1] 5.914

Il risultato è il valore previsto fit con i due limiti di previsione lwr e upr. Il valore residual.scale è la stima \(s\) della deviazione standard \(\sigma\) degli errori, mentre se.fit è la deviazione standard del valore previsto \(Y_0\) cioè \[ s \sqrt{\frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum (x_i - \bar x)^2}} \] Quindi l’errore standard di previsione SE\(_p\) è \[ s \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum (x_i - \bar x)^2}} \] e risulta

sqrt(p$residual.scale^2 + p$se.fit^2)
[1] 5.927

Quindi mentre l’errore standard della stima della media \(E(Y | X = x_0)\) è molto piccola, cioè 0.58 minuti, l’errore di previsione del futuro tempo di attesa è molto maggiore e cioè 5.9 minuti.

library(effects)
plot(allEffects(mod))

Cross-validation

Insieme di convalida

Stima del test error di vari modelli di regressione lineare sul data set Auto. Usiamo la funzione sample() per partizionare l’insieme delle osservazioni in due parti. Selezioniamo un sottoinsieme aleatorio di 196 osservazioni dalle originali 392 osservazioni formando il training set.

library(ISLR)
set.seed(1)
train=sample(392,196)
train
  [1] 105 146 224 354  79 348 365 255 242  24 388  68 262 391 292 188 270 372
 [19] 143 290 387 382 384  47  99 142   5 140 317 124 175 217 178  67 297 239
 [37] 283  39 257 379 289 228 275 194 185 274   9 165 252 238 164 294 149  83
 [55] 383  34 107 174 222 136 304  98 152 110 214  85 157 250  28 356 329 376
 [73] 111 336 330 323 347 123 245 301 333 334 363 101 234  63 218  38  75  44
 [91]  73  18 193 263 233 237 135 121 357 360 192 103 371 287 183  62 305 137
[109] 299 170 276 206 100 295  42   4 198  29 315 362 321 296 131 369 203 122
[127] 312  55  61 326 151  21  10 167 240 154 144 271 251 129 173 380  60  65
[145] 181 112 303 288  26 211 340 385 373 109 120  43 125 313 249  50 359 207
[163] 291 179 201  94  15  76 163 225 386 186 189  86 339 195 311 160 130 300
[181] 307  41 187 106 314  40 284 370 213 247 256 258 261 375  57 117

Usiamo l’opzione subset di lm() per adattare il modello a un sottoinsieme dei dati

lm.fit = lm(mpg ~ horsepower, data=Auto, subset = train)

Con la funzione predict si stima la risposta per tutte le 392 osservazioni e quindi si calcola il MSE per le 196 osservazioni nell’insieme di convalida (il vettore di indici -train seleziona i dati che non sono nel training set).

attach(Auto)
mean((mpg - predict(lm.fit,Auto))[-train]^2)
[1] 26.14

Quindi usiamo delle funzioni quadratiche e cubiche con poly() e stimiamo il test error.

lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
[1] 19.82
lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
[1] 19.78

Se si scegliesse un diverso training set otteremmo dei test error un po’ differenti.

set.seed(2)
train=sample(392,196)
lm.fit=lm(mpg~horsepower,subset=train)
mean((mpg-predict(lm.fit,Auto))[-train]^2)
[1] 23.3
lm.fit2=lm(mpg~poly(horsepower,2),data=Auto,subset=train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
[1] 18.9
lm.fit3=lm(mpg~poly(horsepower,3),data=Auto,subset=train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
[1] 19.26

I risultati sono in accordo con i risultati precedentemente ottenuti:un modello che prevede mpg con una funzione quadratica di horsepower predice meglio si un modello che usa una funzione lineare. Il modello cubico, più complesso` non sembra dare ulteriori miglioramenti.

Leave-One-Out Cross-Validation

La stima LLOCV si può calcolare automaticamente con le funzioni glm() e cv.glm() nel package boot. Notare che glm se non si specifica binomial adatta per default un modello con errori Gaussiani.

glm.fit = glm(mpg~horsepower, data=Auto)
coef(glm.fit)
(Intercept)  horsepower 
    39.9359     -0.1578 
lm.fit=lm(mpg~horsepower,data=Auto)
coef(lm.fit)
(Intercept)  horsepower 
    39.9359     -0.1578

Ora calcoliamo la stima del test error con glm.cv.

library(boot)
glm.fit= glm(mpg~horsepower,data=Auto)
cv.err = cv.glm(Auto,glm.fit)
cv.err$delta
[1] 24.23 24.23

I due numeri in questo caso sono uguali e corrispondono alla stima LOOCV.

Quindi adattiamo una serie di modelli polinomiali di grado crescente.

cv.error=rep(0,5)  # inizializzazione
for (i in 1:5){
 glm.fit=glm(mpg ~ poly(horsepower,i),data=Auto)
 cv.error[i] = cv.glm(Auto, glm.fit)$delta[1]
 }
cv.error
[1] 24.23 19.25 19.33 19.42 19.03

Il MSE stimato scende subito passando da un moello lineare a un modello quadratico, ma poi non sembra migliiorare molto all’aumentare del grado.

k-Fold Cross-Validation

La funzione cv.glm() permette di implementare la convalida incrociata con \(K\) blocchi. Qui usiamo \(K = 10\) e l’applichiamo a modelli per i dati Auto.

Cambiamo il random seed e definiamo un vettore iniziale che conterrà gli errori CV dei modelli polinomiali di gradi da 1 a 10.

set.seed(17)
cv.error.10 = rep(0,10)
for (i in 1:10){
 glm.fit=glm(mpg ~ poly(horsepower,i),data=Auto)
 cv.error.10[i] = cv.glm(Auto,glm.fit, K=10)$delta[1]
 }
cv.error.10
 [1] 24.21 19.19 19.31 19.34 18.88 19.02 18.90 19.71 18.95 19.50

Esercizi 9

  1. Stimare il test error del modello di regressione logistica di default da balance e income usando un insieme di convalida. Fissate set.seed all’inizio dell’analisi. La procedura è la seguente.
  1. Adattare un modello di regressione logistica che usa income e balance per prevedere default.

  2. Usando l’insieme di convalida stimare il test error del modello. Ecco i passi:

    • partizionare i dati nel training set e validation set.
    • Adattare il modello solo al training set.
    • Prevedere il default per ogni individuo nell’insieme di convalida calcolando la probablità a posteriori di default e classificando l’individuo alla classe default se la probabilità a posteriori è maggiore di 1/2.
    • Calcolare il validation set error come la proporzione di osservazioni nell’insieme di convalida che sono state classificate in modo errato.
  3. Ripetere il processo del passo b. tre volte, usando tre diverse partizioni casuali delle osservazioni. Commentate i risultati ottenuti.

  4. Quindi considerate il modello di regressione logistica che prevede la probabilità di default usando income, balance e student e stimate il test error con il metodo precedente. Verificare se l’introduzione della variabile student comporta la riduzione del test error.