Regressione lineare con R

Strumenti quantitativi per la gestione

Emanuele Taufer

Esempio: Diamond

Il dataset diamond dalla libreria UsingR contiene i prezzi di una serie di \(48\) diamanti (in dollari di Singapore) ed il loro peso in carati.

Per caricare i dati è necessario installare prima il pacchetto UsingR.

library(UsingR)
data(diamond)

Con i comandi names e head è possibile vedere, rispettivamente, i nomi delle variabili del dataset e le prime righe di dati

names(diamond)
## [1] "carat" "price"
head(diamond)
##   carat price
## 1  0.17   355
## 2  0.16   328
## 3  0.17   350
## 4  0.18   325
## 5  0.25   642
## 6  0.16   342

Alcune statistiche

summary(diamond)
     carat            price       
 Min.   :0.1200   Min.   : 223.0  
 1st Qu.:0.1600   1st Qu.: 337.5  
 Median :0.1800   Median : 428.5  
 Mean   :0.2042   Mean   : 500.1  
 3rd Qu.:0.2500   3rd Qu.: 657.0  
 Max.   :0.3500   Max.   :1086.0  

Plot

plot(diamond$price,diamond$carat,
  xlab = "Massa (carati)", ylab = "Prezzo (SIN $)",
  bg = "red", col = "black", cex = 1.5, pch = 21,
  frame = FALSE)

RL con R

La stima di un modello di regressione lineare in R viene fatta utilizzando la funzione lm()

La funzione lm() ha due argomenti base: l’equazione del modello che si vuole stimare (formula) e il nome del dataset dove trovare i dati.

La formula è espressa come \[Y \sim X.\] Vedremo che nel caso della regressione multipla sarà semplicemente estesa con \[Y \sim X_1+X_2 + \dots + X_q.\]

Nell’esempio Diamond la sintassi è

reg<-lm(price~carat,data=diamond)

Il codice sopra produce un oggetto R (reg) che contiene i risultati della stima. Per accedere ai risultati principali si può usare

summary(reg)

RLS con R e output

reg<-lm(price ~ carat, data = diamond)
summary(reg)

Call:
lm(formula = price ~ carat, data = diamond)

Residuals:
    Min      1Q  Median      3Q     Max 
-85.159 -21.448  -0.869  18.972  79.370 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -259.63      17.32  -14.99   <2e-16 ***
carat        3721.02      81.79   45.50   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31.84 on 46 degrees of freedom
Multiple R-squared:  0.9783,    Adjusted R-squared:  0.9778 
F-statistic:  2070 on 1 and 46 DF,  p-value: < 2.2e-16

Nell’output troviamo:

Interpretazione dei risultati

Retta stimata

La funzione abline(), con l’oggetto creato da lm() come argomento, produce la retta stimata nel grafico

plot(diamond$carat, diamond$price,
  xlab = "Massa (carati)",
  ylab = "Prezzo (SIN $)",
  bg = "red",
  col = "black", cex = 1.1, pch = 21,frame = FALSE)
abline(reg, lwd = 2)

IC per la retta

Per ottenere gli IC al 95% per i parametri della retta si usi (l’argomento è sempre l’oggetto reg)

confint(reg)
               2.5 %    97.5 %
(Intercept) -294.487 -224.7649
carat       3556.398 3885.6513

Come ottenere un’intercetta più interpretabile

Nel codice sotto, la funzione I() è definita wrapper e obbiga R a effettuare i calcoli definiti in I() (in questo caso calcolo delle deviazioni dalla media) prima di stimare il modello.

reg2 <- lm(price ~ I(carat - mean(carat)), data = diamond)
summary(reg2)
## 
## Call:
## lm(formula = price ~ I(carat - mean(carat)), data = diamond)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -85.159 -21.448  -0.869  18.972  79.370 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             500.083      4.596   108.8   <2e-16 ***
## I(carat - mean(carat)) 3721.025     81.786    45.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.84 on 46 degrees of freedom
## Multiple R-squared:  0.9783, Adjusted R-squared:  0.9778 
## F-statistic:  2070 on 1 and 46 DF,  p-value: < 2.2e-16

Interpolazione

Per calcolare i valori interpolati dalla retta stimata per i training data, ci basta semplicemente (si ricordi che abbiamo chiamato reg l’oggetto che contiene i risultati della procedura lm)

interp<-predict(reg)

Il codice sotto costruisce un piccolo data.frame, df, in cui sono inseriti i dati e interpolazioni

df<-data.frame("Prezzo"=diamond$price, "Carati"=diamond$carat,"Prezzo modello"=interp)
head(df)
##   Prezzo Carati Prezzo.modello
## 1    355   0.17       372.9483
## 2    328   0.16       335.7381
## 3    350   0.17       372.9483
## 4    325   0.18       410.1586
## 5    642   0.25       670.6303
## 6    342   0.16       335.7381

Previsione

Per fare le previsioni sulla base di un predittore \(x_0\) non presente nei training data è necessario creare un data.frame che contiene i valori dei predittori

Proviamo a prevedere il prezzo per dei diamanti di 1, 2 e 3 carati

dfP<-data.frame("carat"=c(1,2,3))
predict(reg,dfP)
        1         2         3 
 3461.399  7182.424 10903.449 

E’ possibile ottenere gli IC per il valor medio previsto con

predict(reg,dfP,interval="confidence")
        fit       lwr      upr
1  3461.399  3330.058  3592.74
2  7182.424  6886.637  7478.21
3 10903.449 10443.088 11363.81

In questo caso gli intervalli (lwr, upr) - di livello 0.95 - forniscono un range di valori per il prezzo medio, dato il valore del predittore carat

Se nell’intervallo vogliamo considerare anche la variabilità data dal termine di errore \(\varepsilon\) (ossia la parte irriducibile) è necessario specificare l’opzione prediction

predict(reg,dfP,interval="prediction")
        fit       lwr       upr
1  3461.399  3315.254  3607.544
2  7182.424  6879.773  7485.074
3 10903.449 10438.648 11368.250

Si noti che gli intervalli così ottenuti sono più ampi rispetto ai precedenti