I dati whiteside del pacchetto MASS
registrano il consumo di gas (Gas) e la temperatura esterna
(Temp) in un’abitazione, prima e dopo un intervento di
isolamento termico (Insul: “Before” / “After”).
#> 'data.frame': 56 obs. of 3 variables:
#> $ Insul: Factor w/ 2 levels "Before","After": 1 1 1 1 1 1 1 1 1 1 ...
#> $ Temp : num -0.8 -0.7 0.4 2.5 2.9 3.2 3.6 3.9 4.2 4.3 ...
#> $ Gas : num 7.2 6.9 6.4 6 5.8 5.8 5.6 4.7 5.8 5.2 ...
📌 L’analisi dei residui del modello di regressione lineare semplice (avente solo temperatura come predittore) ha evidenziato alcune problematiche. Una di esse potrebbe essere connessa con la mancata inclusione di una variabile nel modello. Abbiamo a disposizione una variabile dicotomica
Insul, proviamo a valutarne l’impatto sul consumo medio di gas.
Iniziamo con un’analisi esplorativa che metta in luce la relazione tra le tre variabili: coloriamo il diagramma di dispersione distinguendo le due condizioni:
📌 È evidente che i punti “B” (Before) tendono a essere più alti dei punti “A” (After) a parità di temperatura: l’isolamento riduce il consumo di gas. I due gruppi hanno chiaramente andamenti paralleli — un chiaro segnale che
Insulva inclusa nel modello.
Insul è una variabile categoriale con due modalità. Per
includerla nel modello lineare si introduce una variabile
indicatrice (dummy):
\[d_i = \begin{cases} 0 & \text{se } Insul_i = \text{"Before"} \\ 1 & \text{se } Insul_i = \text{"After"} \end{cases}\]
Denotando con \(Y_i = Gas_i\) e con \(X_i = Temp_i\), il modello diventa con l’inclusione di \(d_i\):
\[\mathbb{E}(Gas_i|Temp_i, d_i) = \beta_0 + \beta_1 Temp_i + \beta_2 d_i\]
che equivale a:
\[\mathbb{E}(Gas_i|Temp_i, d_i) = \begin{cases} \beta_0 + \beta_1 Temp_i & \text{prima dell'intervento} \\ (\beta_0 + \beta_2) + \beta_1 Temp_i & \text{dopo l'intervento} \end{cases}\]
L’introduzione di \(d_i\) produce due rette parallele — una per gruppo — con la stessa pendenza \(\beta_1\) ma intercette diverse: \(\beta_0\) (Before) e \(\beta_0 + \beta_2\) (After).
💡 In R non è necessario creare manualmente la variabile indicatrice:
lm()converte automaticamente i fattori in variabili dummy, usando il primo livello come livello di riferimento. Per controllare quale livello è il riferimento si usarelevel()o in alternativa si controlla il primo livello della variabile.
#>
#> Call:
#> lm(formula = Gas ~ Temp + Insul, data = whiteside)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.74236 -0.22291 0.04338 0.24377 0.74314
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.55133 0.11809 55.48 <2e-16 ***
#> Temp -0.33670 0.01776 -18.95 <2e-16 ***
#> InsulAfter -1.56520 0.09705 -16.13 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.3574 on 53 degrees of freedom
#> Multiple R-squared: 0.9097, Adjusted R-squared: 0.9063
#> F-statistic: 267.1 on 2 and 53 DF, p-value: < 2.2e-16
💡 Come nella regressione semplice, la devianza totale si scompone in:
\[ \sum_{i=1}^n (Y_i - \bar{Y})^2 = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 + \sum_{i=1}^n (\hat{Y}_i - \bar{Y})^2\]
Il coefficiente di determinazione:
\[R^2 = \frac{SQR}{SQT} = 1 - \frac{SQE}{SQT} \in [0,1]\]
misura la quota di variabilità spiegata dal modello. Vicino a 1 indica un buon adattamento.
📌 A differenza della regressione lineare semplice il coefficiente \(R^2\) non è più dato dal quadrato del coefficiente di correlazione tra \(X\) e \(Y\).
📌 L’\(R^2\) è molto maggiore rispetto a quello ottenuto con il modello senza la variabile
Insul.
I coefficienti nella regressione sono parziali: misurano l’effetto di un predittore su \(Y\) tenendo ferme l’altra variabile.
La funzione di regressione stimata consente di prevedere il consumo medio di gas per qualsiasi combinazione di temperatura e condizione:
\[\widehat{Gas} = \hat{\beta}_0 + \hat{\beta}_1 \cdot Temp + \hat{\beta}_2 \cdot d\]
#> [1] "Before, Temp = 5: 4.868 mc"
#> [1] "After, Temp = 5: 3.303 mc"
#> 1 2
#> 4.867844 3.302639
📌 L’analisi dei residui evidenzia il migliore adattamento del modello ai dati. Nel grafico di sinistra non sono più evidenti le due nuvole di punti e nel grafico di destra i residui mostrano un notevole adattamento con la distribuzione normale.
Il modello adattato in precedenza impone che le due rette abbiano la stessa pendenza: l’effetto della temperatura sul consumo di gas è identico prima e dopo l’isolamento. Questa assunzione di parallelismo può essere troppo restrittiva — è plausibile che l’isolamento modifichi non solo il livello medio di consumo (intercetta), ma anche la sua sensibilità alla temperatura (pendenza).
Per catturare questo fenomeno si aggiunge al modello il termine di
interazione tra Temp e
Insul:
\[\mathbb{E}(Gas_i \mid Temp_i, d_i) = \beta_0 + \beta_1 Temp_i + \beta_2 d_i + \beta_3 (Temp_i \cdot d_i)\]
che equivale a due rette non parallele:
\[\mathbb{E}(Gas_i \mid Temp_i, d_i) = \begin{cases} \beta_0 + \beta_1 \, Temp_i & \text{prima dell'isolamento} \\ (\beta_0 + \beta_2) + (\beta_1 + \beta_3) \, Temp_i & \text{dopo l'isolamento} \end{cases}\]
Il coefficiente \(\beta_3\) misura la differenza di pendenza tra i due gruppi: se \(\beta_3 \neq 0\) le rette non sono parallele e l’effetto della temperatura dipende dalla condizione di isolamento.
In R il termine di interazione si specifica con
Temp:Insul (solo l’interazione) oppure con
Temp * Insul (effetti principali più interazione,
equivalente a Temp + Insul + Temp:Insul):
#>
#> Call:
#> lm(formula = Gas ~ Temp * Insul, data = whiteside)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.97802 -0.18011 0.03757 0.20930 0.63803
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.85383 0.13596 50.409 < 2e-16 ***
#> Temp -0.39324 0.02249 -17.487 < 2e-16 ***
#> InsulAfter -2.12998 0.18009 -11.827 2.32e-16 ***
#> Temp:InsulAfter 0.11530 0.03211 3.591 0.000731 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.323 on 52 degrees of freedom
#> Multiple R-squared: 0.9277, Adjusted R-squared: 0.9235
#> F-statistic: 222.3 on 3 and 52 DF, p-value: < 2.2e-16
💡 Il modello con interazione sembra adattarsi meglio ai dati (\(R^2\) leggermente più alto, errore standard residuo più basso e il termine aggiuntivo è statisticamente significativo).
Con i risultati ottenuti:
📌 Tutti e quattro i coefficienti sono significativi. Riferendosi al termine di interazione si può dire l’isolamento i) riduce il consumo medio e ii) attenua la relazione tra consumo di gas e temperatura.
Calcoliamo il consumo atteso quando la temperatura è pari a 5 per due case con e senza intervento.
#> 1 2
#> 4.888 3.334
📌 L’analisi dei residui non evidenzia particolari problematiche.
La tabella riporta un confronto tra modelli in termini di \(R^2\) e Residual Standard error. Il
confronto evidenzia che l’inclusione della variabile Insul
nel modello comporta un netto miglioramento (un marcato abbassamento
dell’RSE e un marcato aumento dell’\(R^2\)). L’inclusione del termine di
interazione determina un ulteriore miglioramento, ma meno marcato.
| Modello Temp | Modello Temp + Insul | Modello Temp * Insul | |
|---|---|---|---|
| \(R^2\) | 0.467 | 0.900 | 0.928 |
| \(R^2\) adj | 0.457 | 0.896 | 0.924 |
| RSE | 0.8606 | 0.367 | 0.323 |
| Parametri | 2 | 3 | 4 |
⚠️ L’\(R^2\) aumenta sempre quando si aggiunge una nuova variabile, anche se non significativa. Per penalizzare la complessità del modello si usa l’\(R^2\) corretto (adjusted):
\[R^2_C = 1 - (1 - R^2) \cdot \frac{n-1}{n-p}\]
\(R^2_C\) può diminuire se una nuova variabile aggiunta non contribuisce abbastanza. Raggiunge il massimo quando il modello realizza il compromesso ottimale fra capacità esplicativa e complessità.
Esercizio 1 Il dataset
neonati(disponibile su moodle) registra il peso (in grammi) (Peso), le settimane di gestazione (settimane) e l’abitudine al fumo della madre (fuma) di 32 neonati. Si vuole studiare se il peso dei neonati dipende dalle settimane di gestazione e dall’abitudine al fumo della madre.
Si analizzi la relazione tra le tre variabili con un opportuno strumento grafico. Dall’analisi emergono pattern o relazioni particolari?
Si stimi il modello di regressione lineare semplice considerando solo la variabile
settimanecome predittore. Si interpretino i risultati e si analizzino i residui.Valutare l’inclusione della variabile dicotomica
fuma. Come differiscono i due modelli stimati (coefficienti, \(R^2\), residui)? Ottenere una rappresentazione grafica delle due rette.Si aggiunga il termine di interazione tra settimane e fumo e confrontare questo modello con i precenti. Quale sembra più appropriato per descrivere il peso dei neonati? Perchè?
Con il modello selezionato, calcolare quanto vale la differenza media di peso tra due neonati con madri fumatrici e non, scegliendo un valore delle settimane.
#>
#> Call:
#> lm(formula = Peso ~ settimane, data = neo)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -354.03 -115.09 18.07 100.22 263.34
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2037.00 498.11 -4.089 0.000298 ***
#> settimane 130.82 12.86 10.170 3.09e-11 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 167.3 on 30 degrees of freedom
#> Multiple R-squared: 0.7752, Adjusted R-squared: 0.7677
#> F-statistic: 103.4 on 1 and 30 DF, p-value: 3.085e-11
#>
#> Call:
#> lm(formula = Peso ~ settimane + fuma, data = neo)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -223.693 -92.063 -9.365 79.663 197.507
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2389.573 349.206 -6.843 1.63e-07 ***
#> settimane 143.100 9.128 15.677 1.07e-15 ***
#> fumaS -244.544 41.982 -5.825 2.58e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 115.5 on 29 degrees of freedom
#> Multiple R-squared: 0.8964, Adjusted R-squared: 0.8892
#> F-statistic: 125.4 on 2 and 29 DF, p-value: 5.289e-15
#>
#> Call:
#> lm(formula = Peso ~ settimane * fuma, data = neo)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -228.528 -89.560 0.273 83.629 184.529
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -2546.138 501.067 -5.081 2.22e-05 ***
#> settimane 147.207 13.120 11.220 7.15e-12 ***
#> fumaS 71.574 716.950 0.100 0.921
#> settimane:fumaS -8.178 18.515 -0.442 0.662
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 117.2 on 28 degrees of freedom
#> Multiple R-squared: 0.8971, Adjusted R-squared: 0.8861
#> F-statistic: 81.37 on 3 and 28 DF, p-value: 6.144e-14
# FATTORI (variabili indicatrici automatiche)
df$var <- relevel(df$var, ref = "livello_riferimento")
mod <- lm(Y ~ x + var, data = df)
summary(mod)
coef(mod) # vettore dei coefficienti stimati
fitted(mod) # valori stimati ŷ
residuals(mod) # residui r = y - ŷ
# PREVISIONE
# valori predetti basati sulle x osservate
predict(mod)
# valori predetti basati su x non osservate
predict(mod, newdata = data.frame(x = 5, var = "level"))
# ANALISI DEI RESIDUI
plot(fitted(mod), residuals(mod)) # residui vs ŷ
abline(h = 0, lty = 2)
qqnorm(residuals(mod)) # normalità dei residui
qqline(residuals(mod))
plot(mod) # 4 grafici diagnostici automatici
# Modello con interazione
mod <- lm(Y ~ x * var, data = df)