1 Analisi Multivariata

L’analisi di coppie di variabili permette già di cogliere aspetti di grande interesse, ma l’analisi congiunta di un insieme più ampio di variabili \(x_1, x_2, \dots, x_p\) (con \(p \geq 3\)) consente di individuare relazioni complesse che l’analisi bivariata non sarebbe in grado di rivelare.

I dati multivariati contengono spesso variabili di natura mista (quantitative e categoriali). Vale sempre la regola aurea: definire sempre come factor() le variabili categoriali prima di qualsiasi analisi.

1.1 Più Variabili Quantitative

Quando si analizzano più variabili quantitative, si può partire dall’esame delle relazioni fra coppie di variabili attraverso la matrice di correlazione.

Con \(p\) variabili quantitative si possono calcolare \(\frac{p(p-1)}{2}\) correlazioni distinte. Per organizzarle si definiscono:

  • Matrice di varianza-covarianza \(\Sigma\): elemento generico \(s_{ij} = \mathrm{cov}(x_i, x_j)\). Simmetrica, con le varianze \(\mathrm{var}(x_i)\) sulla diagonale.
  • Matrice di correlazione \(R\): elemento generico \(r_{ij} = \mathrm{cor}(x_i, x_j)\). Simmetrica con 1 sulla diagonale. È preferibile nella fase esplorativa perché i valori sono adimensionali e comparabili.

Le funzioni cov() e cor() applicate a un data frame di variabili numeriche restituiscono direttamente queste matrici. Utilizziamo il dataset iris:

irisnum <- iris[, -5]   # rimuovo Species (fattore)

round(cov(irisnum), 2)
#>              Sepal.Length Sepal.Width Petal.Length Petal.Width
#> Sepal.Length         0.69       -0.04         1.27        0.52
#> Sepal.Width         -0.04        0.19        -0.33       -0.12
#> Petal.Length         1.27       -0.33         3.12        1.30
#> Petal.Width          0.52       -0.12         1.30        0.58
round(cor(irisnum), 2)
#>              Sepal.Length Sepal.Width Petal.Length Petal.Width
#> Sepal.Length         1.00       -0.12         0.87        0.82
#> Sepal.Width         -0.12        1.00        -0.43       -0.37
#> Petal.Length         0.87       -0.43         1.00        0.96
#> Petal.Width          0.82       -0.37         0.96        1.00

📌 Nota: in presenza di valori mancanti usare l’argomento use = "complete.obs" per calcolare la correlazione dopo aver rimosso i NA variabile per variabile.

x <- rnorm(10, 0, 1)
y <- rnorm(10, 0, 1)
x[1] <- NA
y[10] <- NA
cor(x,y)
#> [1] NA
cor(x,y, use = "complete.obs")  
#> [1] -0.4519229

Nel caso di un numero molto elevato di variabili, può essere conveniente rappresentare graficamente la matrice di correlazione, ad esempio, attraverso la funzione ggcorrplot che richiede in input una matrice di correlazione.

# install.packages("ggcorrplot")
library(ggcorrplot)
ggcorrplot(cor(irisnum),
           type   = "lower",
           method = "circle",
           lab    = TRUE,
           title  = "Matrice di correlazione — iris")

💡 L’area dei cerchi è proporzionale al valore assoluto della correlazione; il colore indica il segno. Con type = "lower" si mostra solo il triangolo inferiore (la matrice è simmetrica). Di default type = "square" (lo si provi per esercizio).

Sempre considerando coppie di variabili, è possibile organizzare gli scatterplot in un unico grafico.

pairs(irisnum)

E attraverso l’uso di colori o simboli diversi esplorare la presenza di gruppi.

pairs(irisnum, col = iris$Species, 
      pch = 20)

1.2 Variabili di natura mista

Quando il data frame contiene variabili di natura mista, ggpairs() del pacchetto GGally sceglie automaticamente la rappresentazione più appropriata per ogni coppia:

ggpairs(iris)

> 💡 Sulla diagonale compaiono le distribuzioni marginali di ciascuna variabile (density plot per le continue, barplot per i fattori). Fuori dalla diagonale: scatterplot (coppie quantitative), boxplot condizionato (quantitativa vs fattore). Sono anche riportate le correlazioni marginali e quelle condizionate ai livelli del fattore. È lo strumento esplorativo più completo per dati misti.

⚠️ Dal grafico emerge subito che Petal.Length e Petal.Width sono le variabili più correlate con Sepal.Length (guardare la prima riga o la prima colonna). Le due variabili dei petali, però, sono anche fortemente correlate fra loro (circa 0.96): questo anticipa un problema che incontreremo nella regressione multipla.

Definendo i gruppi tramite il colore in questo caso si ottengono non solo scatterplot e boxplot colorati per gruppi ma anche densità e correlazioni condizionate.

ggpairs(iris,
        aes(colour = Species, alpha = 0.5))

Esercizio 1 — Il dataset cancerdata.csv contiene variabili relative al casi di tumore al colon retto e al consumo di carne in diversi paesi americani per l’anno 2021 dati, fra cui: - meat: consumo di carne pro-capite - arate: tasso di tumori al colon-retto (aggiustato per età) - GDP: PIL pro-capite

  1. Caricare i dati e ottenere una rappresentazione grafica appropriata per l’analisi esplorativa dei dati. Prestare attenzione alle variabili da inserire nell’analisi.
  2. Commentare il grafico ottenuto, interpretando in particolare la relazione tra arate, meat e GPD.
  3. C’è associazione tra il consumo di carne e il numero di tumori? Possiamo concludere che il consumo di carne possa essere uno dei fattori che causano i tumori al colon retto?
cancerdata <- read.csv("cancerdata.csv")
ggpairs(cancerdata[,-c(1:2,6)],
        aes(colour = region, alpha = 0.5))

1.3 La Correlazione Parziale

Sia per le variabili qualitative che quantitative, l’associazione marginale può essere fuorviante. Ad esempio, tra due variabili quantitative un’elevata correlazione non prova l’esistenza di un legame causale: potrebbe essere interamente determinata dall’influenza di una terza variabile comune.

📌 “Association is not causation” — questa distinzione è fondamentale, in particolare quando si lavora con dati osservazionali (non sperimentali) nei quali non si controllano tutte le variabili in gioco.

Riprendiamo l’esempio sui dati relativi l’incidenza di tumori al colo retto dell’esercizio precedente.

La correlazione marginale fra meat e arate è circa \(r = 0.68\): elevata e positiva. Il diagramma di dispersione sembra confermare la relazione:

ggplot(cancerdata, aes(x = meat, y = arate)) +
  geom_point(colour = "red", aes( size = GDP)) +
  geom_text(label   = cancerdata$Code, color = "blue", size = 3,
            nudge_x = 3, nudge_y = 1.75, check_overlap = FALSE) +
  labs(title = "Consumo di carne vs tasso di tumori al colon-retto",
       x = "Consumo di carne pro-capite",
       y = "Tasso tumori colon-retto")

Vi è quindi una associazione fra il consumo di carne e il numero di tumori all’apparato intestinale. Ma è questo un indizio che il consumo di carne possa essere uno dei fattori che causano i tumori al colon-retto? Si osservi anche che sia meat che arate sono fortemente correlate con GDP: i paesi più ricchi hanno sia un consumo di carne più elevato che una maggiore incidenza di tumori (perché la popolazione è mediamente più anziana e meno esposta ad altre malattie fatali). La ricchezza è una variabile confondente.

Introduciamo una misura che sia in grado di misurare la correlazione fra meat e arate avendo eliminato l’influenza della terza variabile che è un fattore comune che determina sia alti livelli di casi di tumore che alti livelli di consumo della carne: si tratterebbe di valutare la relazione fra le due variabili al netto della ricchezza, cioè della variabile GDP.

Questo è possibile introducendo la correlazione parziale fra due variabili al netto di una terza. L’ideea è di calcolare la correlazione fra le variabile avendo prima eliminato l’influenza lineare della terza.

💡 La correlazione parziale fra \(x\) e \(y\) al netto di \(z\) misura la relazione lineare fra le due variabili dopo aver rimosso l’influenza di \(z\) da entrambe.

L’idea è semplice:

  1. Si regredisce meat su GDP e si prendono i residui → consumo di carne al netto del reddito
  2. Si regredisce arate su GDP e si prendono i residui → tasso di tumori al netto del reddito
  3. Si calcola la correlazione fra i due vettori di residui
resarate <- residuals(lm(arate ~ GDP, data = cancerdata))
resmeat  <- residuals(lm(meat  ~ GDP, data = cancerdata))

cor(resarate, resmeat)   # coefficiente di correlazione parziale
#> [1] 0.004255976
ggplot(cancerdata, aes(x = resmeat, y = resarate)) +
  geom_point(colour = "red", size = 2) +
  geom_text(label   = cancerdata$Code, color = "blue", size = 3,
            nudge_x = 3, nudge_y = 1.75, check_overlap = FALSE) +
  labs(title = "Correlazione parziale: meat e arate al netto di GDP",
       x = "Residui: meat | GDP",
       y = "Residui: arate | GDP")

📌 La correlazione parziale è vicina a zero: una volta rimossa l’influenza del PIL, la relazione fra consumo di carne e tasso di tumori scompare. La correlazione marginale osservata era spuria — interamente spiegata dalla variabile confondente GDP.

La correlazione parziale \(r_{xy.z}\) si può calcolare anche direttamente dalle correlazioni semplici:

\[r_{xy.z} = \frac{r_{xy} - r_{xz} \, r_{yz}}{\sqrt{1 - r_{xz}^2}\,\sqrt{1 - r_{yz}^2}}\]

dove \(r_{xy}\), \(r_{xz}\), \(r_{yz}\) sono i coefficienti di correlazione lineare semplice (o totale) fra le rispettive coppie di variabili. Il risultato coincide con quello ottenuto tramite i residui della regressione.

# Calcolo diretto via formula
r_xy <- cor(cancerdata$meat,  cancerdata$arate)
r_xz <- cor(cancerdata$meat,  cancerdata$GDP)
r_yz <- cor(cancerdata$arate, cancerdata$GDP)

r_parziale <- (r_xy - r_xz * r_yz) / (sqrt(1 - r_xz^2) * sqrt(1 - r_yz^2))
r_parziale
#> [1] 0.004255976

2 Dalla Regressione Semplice alla Regressione Multipla

Regressione con una sola esplicativa

Nelle lezioni precedenti abbiamo visto come modellare la media condizionata della variabile risposta \(Y\) data un’unica variabile esplicativa \(x\):

\[E[Y_i|x_i] = \beta_0 + \beta_1 x_i\]

Il coefficiente \(\hat\beta_1\) rappresenta la variazione media di \(Y\) per un aumento unitario di \(x\).

Regressione con due esplicative, una continua e una dicotomica

\[E[Y_i|x_i] = \beta_0 + \beta_1 x_i + \beta_2 d_i \] con \(d_i \in \{0,1\}\) una variabile dicotomica. Il modello ha un’interpretazione geometrica immediata: separando i due gruppi si ottengono due rette parallele con la stessa pendenza \(\beta_1\) ma intercette diverse. Il coefficiente \(\beta_2\) è quindi la differenza media fra i due gruppi, a parità di x.

Inserendo poi l’interazione tra il predittore continuo e quello dicotomico, si ottengono due rette con intercette e pendenze diverse.

E se avessi più predittori continui? e se il gruppo avesse \(k>2\) categorie?

2.1 Regressione lineare multipla

La regressione lineare multipla (di fatto quella già affrontata nel modello con una esplicativa continua e una categoriale) generalizza il modello includendo \(p\) variabili esplicative \(x_1, x_2, \ldots, x_p\).

La funzione di regressione multipla ha la forma:

\[E[Y \mid x_1, x_2, \ldots, x_p] = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\]

Le medie condizionate di \(Y\) giacciono su un iperpiano nello spazio \((p+1)\)-dimensionale. Le variabili \(x_j\) possono essere trasformazioni delle variabili originali (ad esempio il logaritmo o una potenza), il che rende il modello molto più flessibile di quanto possa sembrare. Ad esempio, specificando \(x_2 = x_1^2\) si ottiene la regressione polinomiale.

💡 Le variabili esplicative \(x_j\) possono essere sia quantitative che categoriali (fattori).

Come nel caso della regressione semplice, i parametri \(\beta_0, \beta_1, \ldots, \beta_p\) si stimano con il criterio dei minimi quadrati: si cercano i valori che minimizzano la somma dei quadrati dei residui

\[ L(\beta_0, \beta_1, \ldots, \beta_p) = \sum_{i=1}^{n} \left(Y_i - \beta_0 - \beta_1 x_{i1} - \beta_2 x_{i2} - \ldots - \beta_p x_{ip}\right)^2\]

In notazione matriciale, definendo:

  • \(\mathbf{y}\): vettore \((n \times 1)\) della variabile risposta
  • \(\mathbf{X}\): matrice \((n \times (p+1))\) di regressione (matrice disegno), con una colonna di uni per l’intercetta
  • \(\boldsymbol{\beta}\): vettore \(((p+1) \times 1)\) dei parametri

\[E[\mathbf{y}] = \mathbf{X}\boldsymbol{\beta}\]

La funzione da minimizzare diventa \(L(\boldsymbol{\beta}) = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\).

Derivando rispetto a \(\boldsymbol{\beta}\) e uguagliando a zero si ottiene il sistema delle equazioni normali:

\[-\mathbf{X}^T\mathbf{y} + \mathbf{X}^T\mathbf{X}\boldsymbol{\beta} = \mathbf{0}\]

la cui soluzione è:

\[\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\]

⚠️ La soluzione esiste se e solo se la matrice \(\mathbf{X}^T\mathbf{X}\) è invertibile, cioè se le colonne di \(\mathbf{X}\) sono linearmente indipendenti. Problemi emergono quando alcune variabili esplicative sono perfettamente correlate fra loro (multicollinearità). È buona pratica evitare di inserire nel modello variabili con correlazione molto elevata.

💡 Interpretazione dei coefficienti: \(\hat \beta_j\) rappresenta la variazione media di \(Y\) per un aumento unitario di \(x_j\), tenendo fisse tutte le altre variabili (ceteris paribus). Non si tratta più della pendenza della retta in un grafico bivariato, ma dell’effetto “parziale” di \(x_j\) al netto delle altre esplicative.

2.2 Il Coefficiente di Determinazione \(R^2\) corretto

Come nella regressione semplice, la devianza totale si decompone in:

\[Dev(Y) = DEV_{residua} + DEV_{spiegata} = \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{DEV_{spiegata}}{Dev(Y)} = 1 - \frac{DEV_{residua}}{Dev(Y)} \in [0,1]\]

misura la quota di variabilità spiegata dal modello. Vicino a 1 indica un buon adattamento.

⚠️ \(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à.

2.3 Costruire il Modello di regressione multipla per iris: passo dopo passo

Utilizziamo iris con Sepal.Length come variabile risposta. L’obiettivo è costruire il modello in modo ragionato, partendo dall’esplorazione delle correlazioni.

Quale predittore entra per primo?

Una strategia naturale è inserire prima la variabile più correlata con la risposta:

round(cor(irisnum), 2)["Sepal.Length", ]
#> Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
#>         1.00        -0.12         0.87         0.82

📌 Petal.Length ha la correlazione più alta con Sepal.Length (circa 0.87). È il candidato naturale per entrare per primo nel modello.

mod1 <- lm(Sepal.Length ~ Petal.Length, data = iris)
summary(mod1)
#> 
#> Call:
#> lm(formula = Sepal.Length ~ Petal.Length, data = iris)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.24675 -0.29657 -0.01515  0.27676  1.00269 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   4.30660    0.07839   54.94   <2e-16 ***
#> Petal.Length  0.40892    0.01889   21.65   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4071 on 148 degrees of freedom
#> Multiple R-squared:   0.76,  Adjusted R-squared:  0.7583 
#> F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

💡 Il coefficiente di Petal.Length è positivo e altamente significativo: mediamente, ogni centimetro in più di petalo si associa a un aumento di circa 0.4 cm nel sepalo. L’\(R^2 \approx 0.76\): la lunghezza del petalo da sola spiega il 76% della variabilità di Sepal.Length.

Aggiungiamo Petal.Width

Petal.Width è molto correlata con Sepal.Length (circa 0.82), ma è anche molto correlata con Petal.Length (circa 0.96). Proviamo ad aggiungerla al modello:

mod2 <- lm(Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
summary(mod2)
#> 
#> Call:
#> lm(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.18534 -0.29838 -0.02763  0.28925  1.02320 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   4.19058    0.09705  43.181  < 2e-16 ***
#> Petal.Length  0.54178    0.06928   7.820 9.41e-13 ***
#> Petal.Width  -0.31955    0.16045  -1.992   0.0483 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4031 on 147 degrees of freedom
#> Multiple R-squared:  0.7663, Adjusted R-squared:  0.7631 
#> F-statistic:   241 on 2 and 147 DF,  p-value: < 2.2e-16

⚠️ Confrontando i risultati di mod1 e mod2 è importante notare che: i) il coefficiente di Petal.Length cambia passando da circa 0.47 a circa 0.54, e ii) Petal.Width entra con un coefficiente negativo (circa −0.31), nonostante la sua correlazione marginale con Sepal.Length sia positiva (circa 0.82). Entrambi i coefficienti sono significativi, il che può sembrare rassicurante, ma è in realtà un sintomo del problema: quando due predittori sono fortemente correlati fra loro, il modello stima effetti parziali che compensano la collinearità, producendo coefficienti instabili e difficili da interpretare. Il VIF è lo strumento che ci permette di diagnosticare questo problema in modo sistematico.

Questo fenomeno si chiama multicollinearità: quando due o più predittori sono fortemente correlati fra loro, il modello fatica a distinguere i loro effetti separati.

2.3.1 Il VIF: Variance Inflation Factor

Un modo sistematico per diagnosticare la multicollinearità è il Variance Inflation Factor (VIF), che misura di quanto la varianza di ciascun coefficiente stimato viene “gonfiata” dalla correlazione con gli altri predittori:

\[\mathrm{VIF}_j = \frac{1}{1 - R^2_j}\]

dove \(R^2_j\) è il coefficiente di determinazione della regressione di \(x_j\) su tutte le altre esplicative.

💡 Regola empirica: VIF \(> 5\) (o \(> 10\) secondo criteri più severi) segnala un problema di multicollinearità. VIF \(= 1\) indica assenza completa di correlazione con gli altri predittori.

library(car)
vif(mod2)
#> Petal.Length  Petal.Width 
#>     13.71927     13.71927

📌 In mod2 i VIF di Petal.Length e Petal.Width sono molto elevati: la correlazione quasi perfetta (0.96) fra le due variabili rende le stime instabili e difficilmente interpretabili.

La strategia più semplice ed efficace in questo caso è la selezione parsimoniosa: tenere una sola delle due variabili collineari. Dal momento che Petal.Length è già nel modello e ha una correlazione marginale leggermente più alta con la risposta, si esclude Petal.Width.

⚠️ La multicollinearità non compromette la capacità predittiva del modello (l’\(R^2\) di mod2 non è molto diverso da quello di mod1), ma rende i coefficienti instabili e difficili da interpretare. Un modello più semplice, con predittori poco correlati fra loro, è sempre preferibile.

Aggiungere Sepal.Width

Sepal.Width ha una correlazione marginale bassa con Sepal.Length (circa −0.12), ma è poco correlata con Petal.Length (circa −0.43): aggiunge quindi informazione “ortogonale” rispetto al primo predittore. Questo è il caso ideale per la regressione multipla.

mod3 <- lm(Sepal.Length ~ Petal.Length + Sepal.Width, data = iris)
summary(mod3)
#> 
#> Call:
#> lm(formula = Sepal.Length ~ Petal.Length + Sepal.Width, data = iris)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.96159 -0.23489  0.00077  0.21453  0.78557 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   2.24914    0.24797    9.07 7.04e-16 ***
#> Petal.Length  0.47192    0.01712   27.57  < 2e-16 ***
#> Sepal.Width   0.59552    0.06933    8.59 1.16e-14 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3333 on 147 degrees of freedom
#> Multiple R-squared:  0.8402, Adjusted R-squared:  0.838 
#> F-statistic: 386.4 on 2 and 147 DF,  p-value: < 2.2e-16

📌 Inserendo Sepal.Width l’\(R^2\) sale da 0.76 a circa 0.84 e il coefficiente di Sepal.Width risulta positivo e significativo. La sua correlazione marginale negativa con Sepal.Length era fuorviante: una volta controllata la lunghezza del petalo, una larghezza maggiore del sepalo si associa a sepali più lunghi. Questo è un esempio di come la relazione fra due variabili possa apparire diversa — o addirittura di segno opposto — quando si tiene conto di una terza variabile.

⚠️ Questo fenomeno, in cui la relazione marginale fra due variabili si inverte o scompare una volta controllata una terza, è l’analogo quantitativo del paradosso di Simpson. L’analisi multivariata è lo strumento per andare oltre le relazioni marginali.

La variabile categoriale Species

Abbiamo già lavorato con variabili dicotomiche (0/1) nella regressione. Species estende questo scenario a tre livelli: setosa, versicolor, virginica.

R gestisce automaticamente un fattore a \(k\) livelli creando \(k-1\) variabili dummy, usando il primo livello in ordine alfabetico (setosa) come categoria di riferimento:

Osservazione Speciesversicolor Speciesvirginica
setosa 0 0
versicolor 1 0
virginica 0 1 
mod4 <- lm(Sepal.Length ~ Petal.Length + Sepal.Width + Species, data = iris)
summary(mod4)
#> 
#> Call:
#> lm(formula = Sepal.Length ~ Petal.Length + Sepal.Width + Species, 
#>     data = iris)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.82156 -0.20530  0.00638  0.22645  0.74999 
#> 
#> Coefficients:
#>                   Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        2.39039    0.26227   9.114 5.94e-16 ***
#> Petal.Length       0.77563    0.06425  12.073  < 2e-16 ***
#> Sepal.Width        0.43222    0.08139   5.310 4.03e-07 ***
#> Speciesversicolor -0.95581    0.21520  -4.442 1.76e-05 ***
#> Speciesvirginica  -1.39410    0.28566  -4.880 2.76e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3103 on 145 degrees of freedom
#> Multiple R-squared:  0.8633, Adjusted R-squared:  0.8595 
#> F-statistic: 228.9 on 4 and 145 DF,  p-value: < 2.2e-16

💡 Interpretazione dei coefficienti:

  • (Intercept): valore medio stimato di Sepal.Length per un fiore di setosa con Petal.Length = 0 e Sepal.Width = 0 (valore teorico, fuori dal range dei dati).
  • Petal.Length: a parità di specie e larghezza del sepalo, ogni cm in più di petalo aumenta mediamente il sepalo di questo valore.
  • Sepal.Width: a parità di specie e lunghezza del petalo, ogni cm in più di larghezza del sepalo aumenta mediamente il sepalo di questo valore.
  • Speciesversicolor: differenza media di Sepal.Length fra versicolor e setosa, a parità delle altre variabili.
  • Speciesvirginica: differenza media di Sepal.Length fra virginica e setosa, a parità delle altre variabili.

2.4 Regressione Multipla: Un Esempio Guidato con R

I dati provengono dalla Survey of Consumer Finances (SCF), un’indagine condotta negli USA su 275 cittadini che hanno acquistato una polizza vita temporanea in caso di morte (term life insurance). L’obiettivo è determinare quali caratteristiche influenzino l’ammontare della copertura assicurativa.

La variabile risposta è FACE: il valore che la compagnia pagherà in caso di morte dell’assicurato.

Variabile Descrizione
AGE Età del rispondente
MARSTAT Stato civile (single / not single)
EDUCATION Anni di istruzione
NUMHH Numero di componenti del nucleo familiare
INCOME Reddito annuo familiare
FACE Valore della polizza 
TL <- read.csv("TL.csv", header = TRUE, sep = ",", row.names = 1)
str(TL)
#> 'data.frame':    275 obs. of  6 variables:
#>  $ AGE      : int  30 50 39 43 34 29 72 51 58 73 ...
#>  $ MARSTAT  : chr  "not single" "not single" "not single" "not single" ...
#>  $ EDUCATION: int  16 9 16 17 11 16 17 16 14 12 ...
#>  $ NUMHH    : int  3 3 5 4 4 3 2 4 1 2 ...
#>  $ INCOME   : int  43000 12000 120000 40000 28000 100000 112000 15000 32000 25000 ...
#>  $ FACE     : int  20000 130000 1500000 50000 220000 600000 100000 2500000 250000 50000 ...
TL$MARSTAT <- factor(TL$MARSTAT)

Analisi esplorativa e primo modello

ggpairs(TL)

📌 Notiamo la forte asimmetria positiva delle variabili FACE e INCOME: la media è molto più alta della mediana, e lo scatterplot mostra la grande maggioranza dei punti concentrati in basso a sinistra con pochi valori molto elevati.

Proviamo il modello lineare semplice selezionando come variabile esplicativa INCOME, una delle variabili esplicative più correlate con la risposta e di interesse:

m0 <- lm(FACE ~ INCOME, data = TL)
summary(m0)
#> 
#> Call:
#> lm(formula = FACE ~ INCOME, data = TL)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -3069185  -628950  -551689  -167579 12976542 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 6.553e+05  1.019e+05   6.433 5.59e-10 ***
#> INCOME      4.414e-01  1.200e-01   3.677 0.000284 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1637000 on 273 degrees of freedom
#> Multiple R-squared:  0.04718,    Adjusted R-squared:  0.04369 
#> F-statistic: 13.52 on 1 and 273 DF,  p-value: 0.0002843

💡 Il coefficiente di INCOME è positivo e significativo (ogni dollaro in più di reddito si associa a un aumento medio di circa 0.44 dollari nella polizza), ma l’\(R^2\) è quasi 0: il modello ha pochissima capacità esplicativa. I grafici dei residui confermano che qualcosa non va.

par(mfrow = c(1, 2))
plot(m0, which = 1:2)

par(mfrow = c(1, 1))

⚠️ I residui mostrano un pattern sistematico (non sono distribuiti in modo casuale attorno a zero) e si discostano chiaramente dalla normale nel QQ-plot. La causa principale è la forte asimmetria delle variabili.

Trasformazione logaritmica

Una soluzione classica per variabili con forte asimmetria positiva è la trasformazione logaritmica, che riduce l’asimmetria e linearizza relazioni di tipo moltiplicativo.

TL$LFACE   <- log(TL$FACE)
TL$LINCOME <- log(TL$INCOME)

ggpairs(TL)

📌 Sulla scala logaritmica la nuvola di punti appare molto più regolare e la relazione lineare è più evidente.

Il nuovo modello è:

\[\mathbb{E}[\log(\mathtt{FACE}_i)] = \beta_0 + \beta_1 \log(\mathtt{INCOME}_i)\]

m1 <- lm(LFACE ~ LINCOME, data = TL)
summary(m1)
#> 
#> Call:
#> lm(formula = LFACE ~ LINCOME, data = TL)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.1967 -0.8032 -0.0018  0.8954  6.4711 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  4.23003    0.85985   4.920  1.5e-06 ***
#> LINCOME      0.69604    0.07661   9.086  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.642 on 273 degrees of freedom
#> Multiple R-squared:  0.2322, Adjusted R-squared:  0.2294 
#> F-statistic: 82.55 on 1 and 273 DF,  p-value: < 2.2e-16

💡 Interpretazione dei coefficienti su scala logaritmica.

  • \(\hat \beta_0\): quando \(\log(\mathtt{INCOME}) = 0\) (cioè INCOME = 1), il valore atteso di FACE è \(e^{4.23} \approx 69\) dollari. In questo caso l’interpretazione ha più senso che nel caso del modello m0.
  • \(\hat\beta_1 \approx 0.70\): un aumento unitario di \(\log(\mathtt{INCOME})\) corrisponde a un aumento di 0.70 in \(\log(\mathtt{FACE})\), ovvero:

Quindi considerando la variabile risposta nella scala originaria, ad un incremento unitario del reddito (nella scala log) si associa in media il raddoppio della polizza.

\[0.7 \simeq \log(FACE)_{x+1}-\log(FACE)_{x} = \log\left(\frac{FACE_{x+1}}{FACE_{x}}\right) \Rightarrow \left(\frac{FACE_{x+1}}{FACE_{x}}\right) \simeq \exp(0.7) \simeq 2.01. \]

beta <- coef(m1)
plot(TL$LINCOME, TL$LFACE,
     xlab = "log(INCOME)", ylab = "log(FACE)",
     main = "Modello m1: retta stimata",
     pch = 19, cex = 0.6, col = "steelblue")
abline(beta[1], beta[2], col = 2, lwd = 2)

par(mfrow = c(1, 2))
plot(m1, c(1, 2))

par(mfrow = c(1, 1))

📌 I residui sono migliorati rispetto a m0 e l’\(R^2\) è più elevato. Ci sono però ancora margini di miglioramento: il QQ-plot mostra code non del tutto normali, e la variabilità spiegata può aumentare includendo altre variabili.

Aggiungiamo le altre variabili socio-demografiche

m2 <- lm(LFACE ~ LINCOME + MARSTAT + EDUCATION + AGE + NUMHH,
         data = TL)
summary(m2)
#> 
#> Call:
#> lm(formula = LFACE ~ LINCOME + MARSTAT + EDUCATION + AGE + NUMHH, 
#>     data = TL)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5.8207 -0.8693  0.0870  0.8606  4.6046 
#> 
#> Coefficients:
#>                Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    3.365900   0.975654   3.450 0.000651 ***
#> LINCOME        0.466668   0.078772   5.924 9.56e-09 ***
#> MARSTATsingle -0.533526   0.262217  -2.035 0.042864 *  
#> EDUCATION      0.209715   0.038734   5.414 1.36e-07 ***
#> AGE           -0.004268   0.008032  -0.531 0.595659    
#> NUMHH          0.233228   0.074687   3.123 0.001987 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.519 on 269 degrees of freedom
#> Multiple R-squared:  0.3526, Adjusted R-squared:  0.3406 
#> F-statistic:  29.3 on 5 and 269 DF,  p-value: < 2.2e-16

📌 AGE non risulta significativa (p-value elevato): non c’è evidenza che l’età abbia un effetto aggiuntivo sul valore della polizza, una volta controllate le altre variabili. La escludiamo.

par(mfrow = c(1, 2))
plot(m2, c(1, 2))

par(mfrow = c(1, 1))
ggpairs(TL, aes(col = MARSTAT))

m3 <- lm(LFACE ~ LINCOME*MARSTAT + EDUCATION + NUMHH, data = TL)
summary(m3)
#> 
#> Call:
#> lm(formula = LFACE ~ LINCOME * MARSTAT + EDUCATION + NUMHH, data = TL)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5.7693 -0.7619  0.1405  0.9185  4.5721 
#> 
#> Coefficients:
#>                       Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)            3.85894    0.92104   4.190 3.79e-05 ***
#> LINCOME                0.40358    0.08111   4.976 1.16e-06 ***
#> MARSTATsingle         -7.21946    2.58389  -2.794 0.005580 ** 
#> EDUCATION              0.20273    0.03831   5.292 2.51e-07 ***
#> NUMHH                  0.26887    0.06935   3.877 0.000133 ***
#> LINCOME:MARSTATsingle  0.63641    0.24392   2.609 0.009587 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.501 on 269 degrees of freedom
#> Multiple R-squared:  0.3679, Adjusted R-squared:  0.3562 
#> F-statistic: 31.32 on 5 and 269 DF,  p-value: < 2.2e-16

💡 Equazione finale (categoria di riferimento: not single):

\[\widehat{\log(\mathtt{FACE})} = 3.86 + 0.40\,\log(\mathtt{INCOME}) + 0.20\,\mathtt{EDUCATION} + 0.27\,\mathtt{NUMHH}\]

Per i single:

\[\widehat{\log(\mathtt{FACE})} = -3.35 + 1.04\,\log(\mathtt{INCOME}) + 0.20\,\mathtt{EDUCATION} + 0.27\,\mathtt{NUMHH}\]

I coefficienti di EDUCATION e NUMHH si interpretano a parità di tutte le altre variabili: ogni anno aggiuntivo di istruzione aumenta in media il log-valore della polizza di 0.20, e ogni componente aggiuntivo del nucleo familiare lo aumenta di 0.27.

📌 I grafici dei residui di m3 sono molto migliorati rispetto ai modelli precedenti. Sono ancora presenti alcuni valori anomali, ma il modello appare nel complesso accettabile. L’\(R^2\) è aumentato significativamente rispetto al modello semplice m1.

r2_scf <- c(
  m1 = summary(m1)$r.squared,
  m2 = summary(m2)$r.squared,
  m3 = summary(m3)$r.squared
)
round(r2_scf, 3)
#>    m1    m2    m3 
#> 0.232 0.353 0.368