1 Modello di regressione con una variabile continua e una dicotomica

1.1 Framework

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.

1.2 Diagramma di dispersione con la variabile Isolamento

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 Insul va inclusa nel modello.

1.3 La Variabile Indicatrice

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 usa relevel() o in alternativa si controlla il primo livello della variabile.

1.4 Adattamento del Modello

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

1.4.1 Interpretazione dei Coefficienti

I coefficienti nella regressione sono parziali: misurano l’effetto di un predittore su \(Y\) tenendo ferme l’altra variabile.

  • \(\hat{\beta}_0 \approx 6.55\): consumo medio di gas quando la temperatura è 0°C e si è prima dell’isolamento. Ha senso perché 0°C rientra nel range dei dati.
  • \(\hat{\beta}_1 \approx -0.34\): per ogni grado in più di temperatura, il consumo medio di gas diminuisce di circa 0.34 mc, indipendentemente dal periodo (Before/After). È il coefficiente di regressione parziale di Temp al netto di Insul.
  • \(\hat{\beta}_2 \approx -1.57\): dopo l’isolamento il consumo medio di gas è inferiore di circa 1.57 mc rispetto a prima, a parità di temperatura.

1.5 Visualizzazione delle rette di regressione

1.6 Previsione

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

1.6.1 Analisi dei residui

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


1.7 Modello con interazione

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

1.7.1 Interpretazione dei Coefficienti

Con i risultati ottenuti:

  • \(\hat{\beta}_0 \approx 6.85\): consumo medio stimato prima dell’isolamento a 0°C. Il valore è direttamente interpretabile perché 0°C rientra nel range dei dati.
  • \(\hat{\beta}_1 \approx -0.39\): pendenza della retta Before — per ogni grado in più di temperatura, il consumo diminuisce di circa 0.39 mc prima dell’isolamento.
  • \(\hat{\beta}_2 \approx -2.13\): differenza di intercetta tra After e Before a 0°C — dopo l’isolamento il consumo è in media 2.13 mc inferiore rispetto a prima, a temperatura fissata.
  • \(\hat{\beta}_3 \approx 0.12\): differenza di pendenza tra After e Before — dopo l’isolamento l’effetto della temperatura sul consumo è attenuato di circa 0.12 mc per grado rispetto a prima. La pendenza della retta After è quindi \(-0.39 + 0.12 = -0.27\).

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

1.7.2 Visualizzazione delle Rette Non Parallele

1.7.3 Previsione con il Modello ad Interazione

Calcoliamo il consumo atteso quando la temperatura è pari a 5 per due case con e senza intervento.

#>     1     2 
#> 4.888 3.334

1.7.4 Analisi dei residui

📌 L’analisi dei residui non evidenzia particolari problematiche.

1.7.5 Confronto tra i modelli

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.

  1. Si analizzi la relazione tra le tre variabili con un opportuno strumento grafico. Dall’analisi emergono pattern o relazioni particolari?

  2. Si stimi il modello di regressione lineare semplice considerando solo la variabile settimane come predittore. Si interpretino i risultati e si analizzino i residui.

  3. Valutare l’inclusione della variabile dicotomica fuma. Come differiscono i due modelli stimati (coefficienti, \(R^2\), residui)? Ottenere una rappresentazione grafica delle due rette.

  4. 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è?

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