Regressione lineare in R: lm(), summary() e predict()

L'idea: la retta dei minimi quadrati

La regressione lineare fa passare una retta attraverso una nuvola di punti: y = intercept + slope × x. Tra tutte le rette possibili, lm() sceglie quella che minimizza la somma dei residui al quadrato, dove un residuo è la distanza verticale tra un punto e la retta. Elevare al quadrato fa pesare in modo sproporzionato gli errori grandi, ed è per questo che un solo outlier estremo può inclinare l'intera stima.

Mentre la correlazione ti dà un unico numero senza unità per dire "quanto si muovono insieme", la regressione ti dà un'equazione: con le unità, una pendenza che puoi interpretare e gli strumenti per fare previsioni.

fit <- lm(mpg ~ wt, data = mtcars)
print(coef(fit))

Leggi la formula come "modella mpg in funzione di wt". I due coefficienti sono la retta stimata: mpg ≈ 37.3 − 5.3 × peso. La pendenza ha unità reali: ogni 1000 libbre in più di auto (wt è in unità da 1000 libbre) costano circa 5.3 miglia per gallone. L'intercetta (37.3 mpg a peso zero) è solo il punto in cui la retta incrocia lo zero; nessuna auto pesa niente, quindi non darle troppo peso.

Leggere summary() passo per passo

summary(fit) è l'output che ogni corso di statistica ti chiede di interpretare. Eseguilo, poi affrontalo blocco per blocco:

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

Call: ripete il modello che hai stimato. Banale ora, salvavita quando destreggi sei oggetti modello.

Residuals: il riepilogo a cinque numeri degli scarti (valore reale − valore previsto). Vuoi la mediana vicina a 0 e una simmetria approssimativa tra Min/Max e 1Q/3Q; una forte asimmetria suggerisce che al modello a retta manca qualcosa.

Coefficients: il cuore dell'output, una riga per termine:

• Estimate: il valore stimato. Per wt, −5.34: ogni 1000 libbre in più sono associate a circa 5.3 mpg in meno. Traduci sempre la pendenza in una frase con le unità; quella frase è tutto il contenuto pratico del modello.
• Std. Error: quanto oscillerebbe la stima su campioni ripetuti. Le stime che distano da zero solo un paio di errori standard sono traballanti.
• t value: Estimate ÷ Std. Error, cioè a quanti errori standard da zero si trova il coefficiente (−9.56 qui).
• Pr(>|t|): il p-value per la domanda "questo coefficiente potrebbe essere davvero zero?". Per wt è circa 1.3e-10: se il peso non avesse davvero alcuna relazione lineare con mpg, una pendenza così ripida non comparirebbe praticamente mai in un campione di 32. P-value piccolo significa prove che l'associazione esiste, non la dimostrazione che il modello sia corretto, e nemmeno una misura di importanza (anche un effetto minuscolo stimato con precisione ottiene un p-value minuscolo).
• Signif. codes / asterischi: un'abbreviazione visiva della colonna dei p-value. Comoda, ma non aggiunge informazioni.

Residual standard error: 3.05 on 30 degrees of freedom: la dimensione tipica di un errore di previsione, nelle unità della risposta: le previsioni sbagliano di solito di circa 3 mpg. Valutalo rispetto alla scala di mpg (che va più o meno da 10 a 34).

Multiple R-squared: 0.75: il peso spiega circa il 75% della varianza di mpg. L'Adjusted R-squared (0.74) ricalcola lo stesso valore con una penalità per predittore, perché la versione grezza può solo aumentare man mano che aggiungi variabili, anche rumore casuale. Quando confronti modelli con un numero diverso di predittori, quello onesto è l'adjusted. E resisti al riflesso "buon modello = R² alto": un effetto davvero utile può vivere in un modello con R² basso (risultato rumoroso, uno tra tanti fattori), mentre un R² alto può venire da un overfitting o da una variabile che fa trapelare la risposta.

F-statistic: 91.4 ... p-value: 1.29e-10: il test sull'intero modello, cioè se questo modello batte "prevedi semplicemente la media per tutti". Con un solo predittore duplica il t-test della pendenza (nota che 9.56² ≈ 91.4); con più predittori diventa il test congiunto che almeno un coefficiente sia diverso da zero. Il suo meccanismo è la stessa scomposizione della varianza dell'ANOVA.

Regressione multipla: a parità degli altri

Aggiungi predittori con +:

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

L'interpretazione cambia in un punto cruciale. Ogni Estimate ora è l'effetto di quel predittore a parità degli altri: il coefficiente di wt (circa −3.9, sceso da −5.3) è il costo in mpg del peso in più confrontando auto con la stessa potenza. Il −5.3 della regressione semplice includeva in silenzio il fatto che le auto più pesanti tendono anche a essere più potenti; la regressione multipla separa le due cose. È anche il motivo per cui i coefficienti cambiano quando aggiungi variabili: se il nuovo predittore è correlato con uno vecchio, il ruolo di quello vecchio cambia. L'R quadro sale a circa 0.83, e qui è la versione adjusted il confronto equo con il modello a un solo predittore.

Previsioni: predict()

Il modello stimato è una funzione; predict() la valuta. Costruisci un data frame newdata i cui nomi di colonna corrispondono esattamente ai predittori:

fit <- lm(mpg ~ wt, data = mtcars)

new_cars <- data.frame(wt = c(2.5, 3.5))
print(predict(fit, newdata = new_cars))

print(predict(fit, newdata = new_cars, interval = "confidence"))
print(predict(fit, newdata = new_cars, interval = "prediction"))

I due tipi di intervallo rispondono a domande diverse, e confonderli è un classico errore d'esame:

• interval = "confidence": l'incertezza sulla media, cioè "per tutte le auto da 2500 libbre, dove sta l'mpg medio?". Stretto, e si restringe man mano che crescono i dati.
• interval = "prediction": l'intervallo in cui è probabile che cada una singola nuova auto di quel peso. Molto più ampio, perché una singola auto porta con sé la propria dispersione intorno alla retta, una dispersione che nessuna quantità di dati riesce a mediare via.

Riportare un intervallo di confidenza quando la domanda riguarda una nuova osservazione sovrastima in modo drastico la tua precisione.

Diagnostica e la trappola dell'estrapolazione

summary() ti dice cosa stima il modello; i grafici dei residui ti dicono se crederci. In una sessione interattiva:

par(mfrow = c(2, 2))
plot(fit)   # four diagnostic plots

Cosa cercare: Residuals vs Fitted dovrebbe essere una nuvola senza forma; una curva significa che la relazione non è rettilinea, un imbuto (dispersione che cresce con i valori stimati) significa varianza non costante, e i tuoi errori standard sono sbagliati. Nel grafico Q-Q i punti dovrebbero aderire alla retta; code pesanti significano che gli outlier stanno distorcendo la stima. Scale-Location è di nuovo il controllo dell'imbuto. Residuals vs Leverage segnala i punti influenti, osservazioni che da sole trascinano i coefficienti (in mtcars, qui tendono a comparire auto particolari come la Chrysler Imperial). Un rapido grafico a dispersione dei dati grezzi prima della stima ti fa accorgere presto della maggior parte di questi problemi.

Infine, la trappola che nessuna diagnostica cattura: l'estrapolazione. Il modello ha imparato da auto che pesano più o meno da 1500 a 5400 libbre. Passa a predict() un wt di 8 e ti restituirà allegramente un mpg negativo: la matematica prolunga la retta all'infinito, ma le prove si fermano al bordo dei dati. Fai previsioni solo dentro (o vicino) all'intervallo su cui hai stimato il modello.

Cosa ti porti a casa

• fit <- lm(y ~ x, data = df) stima la retta dei minimi quadrati; coef(fit) è l'equazione, summary(fit) il report completo.
• Leggi le Estimate come frasi con le unità; Pr(>|t|) chiede "potrebbe essere zero?", non "conta qualcosa?".
• Il residual standard error è l'errore tipico in unità reali; l'adjusted R-squared è il numero equo per confrontare modelli.
• Nella regressione multipla ogni coefficiente significa "a parità degli altri", e i coefficienti cambiano quando entrano predittori correlati.
• predict(fit, newdata, interval = ...): "confidence" per la media, "prediction" per un singolo nuovo caso, quello ampio.
• Controlla plot(fit) per curve, imbuti e punti influenti; non fidarti mai delle previsioni fuori dall'intervallo dei dati.

Prossimo passo: quando il risultato è un sì/no invece di un numero, la regressione logistica con glm().