1 Analisi di Due Variabili (Bivariata)

L’analisi congiunta di due variabili permette di valutare la presenza di relazioni (associazioni): coppie di valori delle due variabili tendono a presentarsi più o meno frequentemente di quanto ci si aspetterebbe se le variabili fossero indipendenti.

La presenza di una relazione è di grande interesse: la conoscenza di una variabile consente di prevedere con minore incertezza l’altra. L’assenza di relazione — indipendenza — implica che l’analisi congiunta non aggiunge nulla rispetto all’analisi delle singole variabili.

Gli strumenti dipendono dalla natura delle variabili:

Tipo di variabili Strumenti principali
Due categoriali Tabella a doppia entrata, statistica \(X^2\)
Quantitativa + fattore Boxplot multipli, scomposizione devianza, ANOVA
Due quantitative Scatterplot, covarianza, correlazione, regressione

2 Due Variabili Categoriali

2.1 Tabella a Doppia Entrata e Distribuzioni Condizionate

Nell’analisi bivariata è spesso utile identificare una variabile risposta (dipendente, \(Y\)) e una variabile esplicativa (indipendente, \(X\)), e studiare la distribuzione della risposta condizionatamente all’esplicativa.

Riprendiamo il dataframe AutoBi del pacchetto insuranceData, con le modifiche fatte nelle precedenti lezioni.

library(insuranceData)
data(AutoBi)

AutoBi$MARITAL   <- factor(AutoBi$MARITAL)
levels(AutoBi$MARITAL) <- c("married", "single",
                             "previouslymarried", "previouslymarried")
AutoBi$LOSSclass <- cut(AutoBi$LOSS, breaks = c(0, 0.5, 2, 4, 8, 1100))
AutoBi$ATTORNEY  <- factor(AutoBi$ATTORNEY)
levels(AutoBi$ATTORNEY) <- c("yes", "no")
AutoBi$CLMSEX    <- factor(AutoBi$CLMSEX)
levels(AutoBi$CLMSEX)   <- c("M", "F")

Analizziamo le frequenze congiunte, le frequenze di ATTORNEY condizionatamente a CLMSEX, e quella di CLMSEX, condizionatamente ad ATTORNEY.

tab1 <- table(AutoBi$ATTORNEY, AutoBi$CLMSEX)
tab1
#>      
#>         M   F
#>   yes 325 352
#>   no  261 390
rtab <- prop.table(tab1, margin = 1)  
ctab <- prop.table(tab1, margin = 2)  
rtab; ctab
#>      
#>               M         F
#>   yes 0.4800591 0.5199409
#>   no  0.4009217 0.5990783
#>      
#>               M         F
#>   yes 0.5546075 0.4743935
#>   no  0.4453925 0.5256065

Quindi produciamo un diagramma a barre per:

  • Distribuzione del ricorso all’avvocato condizionata al sesso
  • Distribuzione del ricorso all’avvocato condizionata alla classe di danno (LOSSclass)
par(mfrow = c(1, 2))

barplot(ctab, legend = TRUE, beside = TRUE, ylim = c(0, 1),
        main = "Avvocato | Sesso")

barplot(prop.table(table(AutoBi$ATTORNEY, AutoBi$LOSSclass), margin = 2),
        beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(horiz = T, x = "top"),
        main = "Avvocato | Classe di danno")

par(mfrow = c(1, 1))

📌 Il secondo grafico mostra chiaramente che il ricorso all’avvocato aumenta al crescere del danno: le distribuzioni condizionate sono molto diverse, suggerendo una forte relazione fra le due variabili.

2.2 La Statistica \(X^2\)

Una tabella a doppia entrata per una coppia \((X,Y)\) di variabili le cui modalità sono indicate per riga (la \(X\)) e per colonna (la \(Y\)), ha la seguente forma generale

\[ \begin{array}{c|cccccc|l} X\,/\,Y& y_1 & y_2 & \ldots & y_k & \ldots & y_K &{\it Tot.}\\ \hline x_1 &n_{11} & n_{12} & \ldots & n_{1k} & \ldots & n_{1K}&n_{1.}\\ x_2 &n_{21} & n_{22} & \ldots & n_{2k} & \ldots & n_{2K}&n_{2.}\\ \vdots & \vdots &\vdots &\,\cdots &\vdots &\cdots & \cdots\\ x_h &n_{h1} & n_{h2} & \ldots & n_{hk} & \ldots & n_{hK}&n_{h.}\\ \vdots& \vdots &\vdots &\,\cdots &\vdots &\cdots & \cdots\\ x_H &n_{H1} & n_{H2} & \ldots & n_{Hk} & \ldots & n_{HK}&n_{H.}\\ \hline\\[-6pt] {\it Tot.}&n_{.1} & n_{.2} & \ldots & n_{.k} & \,\ldots & n_{. K}& n_{}\\ \hline \end{array} \]

💡 Il principio fondamentale è: se le distribuzioni condizionate sono uguali le due variabili sono indipendenti. In caso di indipendenza la frequenza teorica attesa in ogni cella è:

\[\hat{n}_{hk} = \frac{n_{h\cdot} \cdot n_{\cdot k}}{n}\]

La statistica \(X^2\) misura la distanza fra frequenze osservate e attese:

\[X^2 = \sum_{h=1}^{H}\sum_{k=1}^{K} \frac{(n_{hk} - \hat{n}_{hk})^2}{\hat{n}_{hk}}\]

Sotto \(H_0: P(X=x_h, Y = y_k) = P(X=x_h) P(Y=y_k)\) (indipendenza), la statistica si distribuisce asintoticamente come \(X^2 \stackrel{\cdot}\sim \chi^2_{(H-1)(K-1)}\). Un p-value piccolo indica che i dati non supportano l’ipotesi di indipendenza.

⚠️ Valori “piccoli”, ovvero quelli ottenuti quando le frequenze osservate sono molto simili a quelle attese sotto ipotesi di indipendenza, sono una chiara indicazione di indipendenza tra la due variabili. Valori “grandi”, quindi casi in cui le frequenze osservate si discostano molto da quelle ottenute sotto ipotesi di indipendenza, evidenziano una associazione tra le due variabili (qui si analizza un tipo di relazione simmetrica quindi nessuna delle due variabili ha un ruolo sbilanciato: esplicativa vs risposta). La valutazione di “piccolo” o “grande” deve essere effettuata confrontando la statistica con la distribuzione \(\chi^2\) con opportuni gradi di libertà (oppure tramite il p-value).

Il valore della statistica \(X^2\) non può essere direttamente utilizzato per valutare la forza dell’associazione. Per questo esiste la sua versione normalizzata

\[\tilde{X}^2 = \frac{X^2}{n \min(K-1, H-1)}\] Questo indice normalizzato varia tra 0 e 1 e più è grande più indica forza nell’associazione.

💡 In R la funzione per ottenere la statistica \(X^2\) è chisq.test. Essa richiede come argomenti o una tabella a doppia entrata (specifica solo argomento x) o due vettori di factor (specfica argomenti x e y).

⚠️ Per ottenere la statistica \(X^2\) come da formula sopra, occorre disabilitare la correzione per continuità di Yates (di default in R), impostando correct = FALSE.

Vediamo sia il calcolo manuale che quello ottenuto con la funzione chisq.test.

tot  <- sum(tab1)
mr   <- margin.table(tab1, 1)        # marginale di riga
mc   <- margin.table(tab1, 2)        # marginale di colonna
hatn <- mr %*% t(mc) / tot              # frequenze teoriche di indipendenza 
hatn
#>      
#>              M        F
#>   yes 298.7364 378.2636
#>   no  287.2636 363.7364
X2 <- sum((tab1 - hatn)^2 / hatn)
X2
#> [1] 8.430048
X2.norm <- X2/(tot*min(length(mr)-1, length(mc)-1))
X2.norm
#> [1] 0.006347928
# Equivalente con chisq.test()
chisq.test(tab1, correct = FALSE)
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  tab1
#> X-squared = 8.43, df = 1, p-value = 0.003691
# Test per avvocato vs classe di danno
tab2 <- table(AutoBi$ATTORNEY, AutoBi$LOSSclass)
chisq.test(tab2, correct = FALSE)
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  tab2
#> X-squared = 342.89, df = 4, p-value < 2.2e-16
X2.LOSS <- chisq.test(tab2, correct = FALSE)$statistic
X2.norm <- X2.LOSS/(tot*min(nrow(tab2)-1, ncol(tab2)-1))
X2.norm
#> X-squared 
#> 0.2581973

📌 Per avvocato–sesso il p-value è ~0.004: c’è evidenza di associazione ma molto piccola (dal X2 normalizzato). Per avvocato–danno il p-value è \(< 2 \times 10^{-16}\): la relazione è forte e i dati non supportano l’ipotesi di indipendenza.

Esercizio 1 Il dataset Titanic contiene una tabella a 4 entrate (array 4-dimensionale) che riporta il numero di passeggeri sopravvisuti o meno al disastro del Titanic, per classe di viaggio, età e genere. Si vuole investigare se via sia una associazione tra la sopravvivenza e la classe di viaggio. A tal proposito, si analizzino graficamente le distribuzioni condizionate e le si confrontino con le distribuzioni marginali. Quindi si utilizzi la statistica \(X^2\) per valutare se emerge una relazione tra queste due variabili.

data("Titanic")
tab <- margin.table(Titanic, margin = c(1, 4))

3 Variabile Quantitativa Condizionata a un Fattore

Quando la variabile risposta è quantitativa e la variabile esplicativa è un fattore, si vuole studiare come la distribuzione di \(Y\) cambia nei gruppi definiti dal fattore. Se le distribuzioni sono uguali le due variabili sono indipendenti.

Le possibili situazioni sono:

  1. Le distribuzioni nei gruppi differiscono solo per la tendenza centrale (media o mediana), avendo la stessa forma
  2. Differiscono per tendenza centrale e variabilità, mantenendo la stessa forma
  3. Differiscono anche nella forma (simmetria, curtosi)

I boxplot multipli, il confronto delle ECDF e il QQ-plot sono gli strumenti grafici più utili per capire quale situazione si presenta.

3.1 Scomposizione della Devianza

💡 La devianza totale si decompone in due componenti:

\[DEV_{tot} = \underbrace{\sum_{g=1}^{G}\sum_{i=1}^{n_g}(x_i^{(g)} - \bar x_g)^2}_{DEV_{within}} + \underbrace{\sum_{g=1}^{G} n_g(\bar x_g - \bar x)^2}_{DEV_{between}}\]

  • \(DEV_{within}\): variabilità entro i gruppi

  • \(DEV_{between}\): variabilità fra i gruppi (scostamento delle medie di gruppo dalla media generale)

  • \(G\): numero gruppi

  • \(n_g\): numero osservazioni per il gruppo \(g\), \(g \in \{1, \ldots, G\}\)

  • \(x_i^{(g)}\): osservazioni appartenenti al gruppo \(g\)

  • \(\bar x_g\): media del gruppo \(g\)

  • \(\bar x\): media globale

💡 Il rapporto \[\eta^2 = DEV_{between}/DEV_{tot} \in [0,1]\] misura quanto la variabile di gruppo spiega della variabilità complessiva. Vicino a 1 indica gruppi ben separati.

3.2 ANOVA — Analisi della Varianza

💡 L’ ANOVA verifica formalmente \(H_0: \mu_1 = \mu_2 = \cdots = \mu_K\) (uguaglianza delle medie di \(K\) gruppi), assumendo che i dati provengano da gaussiane indipendenti con la stessa varianza (ipotesi di omoschedasticità). La statistica test è:

\[F = \frac{DEV_{between}/(K-1)}{DEV_{within}/(n-K)} \sim F_{K-1,\, n-K} \quad \text{sotto } H_0\]

Valori elevati di \(F\) indicano che i dati non supportano \(H_0\).

💡 In R il comando principale per ottenere l’analisi della varianza è aov al quale si passa in argomento la formula continua ~ factor e il data frame. Una volta ottenuto il risultato mediante la funzione summary applicata all’oggetto risultante da aov si ottiene la tabella ANOVA che scompone la variabilità totale della risposta in una componente dovuta al fattore (tra gruppi) e una componente residua (entro gruppi). Vengono riportati i gradi di libertà (Df), le devianze (Sum Sq) e le devianze medie (Mean Sq). Il rapporto tra la variabilità tra gruppi e quella residua produce la statistica F, utilizzata per testare l’ipotesi di uguaglianza delle medie. Il p-value associato indica se il fattore ha un effetto statisticamente significativo sulla variabile risposta. I gradi di libertà riportati nella prima riga sono pari al numero di livelli del fattore meno uno, mentre quelli riportati nella seconda riga sono pari al numero di osservazioni meno il numero di livelli del fattore.

3.3 Esempio IRIS

Consideriamo il dataframe iris e verifichiamo se vi sono differenze tra le medie di Sepal.Length nei gruppi definiti da Species

data(iris)
boxplot(Sepal.Length ~ Species, data = iris,
        col = c("lightblue", "blue", "darkblue"))
points(1:3, 
       tapply(iris$Sepal.Length,iris$Species, mean), pch = 7, col = "red")

testanova <- aov(Sepal.Length ~ Species, data = iris)
summary(testanova)
#>              Df Sum Sq Mean Sq F value Pr(>F)    
#> Species       2  63.21  31.606   119.3 <2e-16 ***
#> Residuals   147  38.96   0.265                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Verifica: DEVbetween + DEVwithin = DEVtot
n      <- nrow(iris)
Devtot <- var(iris$Sepal.Length) * (n - 1)
Devtot   # deve coincidere con la somma delle due Sum Sq
#> [1] 102.1683
DevBeetw <- summary(testanova)[[1]]["Species", "Sum Sq"] # Devianza between

eta2 <- DevBeetw/Devtot  # eta quadro
eta2 
#> [1] 0.6187057

📌 Il p-value è \(< 2 \times 10^{-16}\): evidenza a supporto dell’ipotesi che le medie della lunghezza del sepalo differiscano fra le tre specie. La colonna Sum Sq riporta \(DEV_{between}\) (riga Species) e \(DEV_{within}\) (riga Residuals).

📌 Il coefficiente \(\eta^2\) pari a 0.619 indica che oltre il 61% della varianza della lunghezza dei sepali è riconducibile alla diversità tra le specie. Questo suggerisce che la specie di appartenenza è un fattore determinante nel definire le dimensioni del sepalo.

📌 Nonostante il test ANOVA confermi con estrema precisione (p-value quasi nullo) l’esistenza di differenze tra le medie, è il calcolo dell’\(\eta^2\) a confermare l’importanza sostanziale del fenomeno: la specie spiega infatti la maggioranza della variabilità osservata (\(62\%\)), lasciando solo il \(38\%\) alla variabilità naturale intra-specifica (residua).

Esercizio 2 Il dataset chickwts contiene i pesi di pulcini di sei settimane che sono stati sottoposti a sei diversi tipi di alimentazione. L’obiettivo dell’analisi è determinare se il tipo di mangime influenzi in modo significativo la crescita dei pulcini.

  1. In fase esplorativa si produca un grafico che mostri la distribuzione del peso (weight) in funzione del tipo di alimentazione (feed) e si commenti se sembrano esserci differenze tra i gruppi.
  2. Si esegua un’Analisi della Varianza per testare l’ipotesi nulla che le medie dei pesi sono uguali per tutti i tipi di mangime. Cosa possiamo concludere?
  3. Si calcoli l’indice Eta-Quadro (\(\eta^2\)) e si commenti il risultato ottenuto.

4 Due Variabili Quantitative

4.1 Il Diagramma di Dispersione

💡 Con due variabili quantitative i dati si presentano come coppie \((x_i, y_i)\) e si rappresentano come punti nel piano — il diagramma di dispersione (scatterplot).

Consideriamo il dataframe whiteside del pacchetto MASS. Visulizziamo mediante la funzione plot() il diagramma di dispersione di Temp (temperatura esterna) e Gas (Consumo di gas).

library(MASS)
data(whiteside)
attach(whiteside)

plot(Temp, Gas, pch = 19,
     xlab = "Temperatura (°C)", ylab = "Consumo di gas (mc)",
     main = "Consumo di gas vs temperatura esterna")

📌 La nuvola è orientata chiaramente verso il basso: al crescere della temperatura il consumo di gas diminuisce — relazione lineare negativa, attesa.

4.2 Covarianza e Correlazione

💡 La covarianza misura la tendenza delle due variabili ad associarsi linearmente:

\[cov(X, Y) = \sum_{i=1}^{n} \frac{(x_i - \bar x)(y_i - \bar y)}{n} = \frac{1}{n}\sum_{i=1}^{n} x_iy_i - \bar x\bar y\]

Il suo valore assoluto dipende dall’unità di misura. Il coefficiente di correlazione lineare standardizza la covarianza rendendola adimensionale e compresa in \([-1, 1]\):

\[r(X, Y) = \frac{cov(X, Y)}{sd(X) \cdot sd(Y)}\]

  • \(r = 1\): perfetta relazione lineare positiva
  • \(r = -1\): perfetta relazione lineare negativa
  • \(r = 0\): assenza di relazione lineare

⚠️ In R cov() e var() dividono per \(n-1\). La covarianza descrittiva (divisa per \(n\)) va calcolata manualmente. Il problema non si pone con la correlazione.

n <- length(Temp)

# Covarianza descrittiva (÷n) — calcolo manuale
codev  <- sum((Temp - mean(Temp)) * (Gas - mean(Gas)))
covar  <- codev / n
print(paste("Covarianza (÷n):", round(covar, 4)))
#> [1] "Covarianza (÷n): -2.1548"
# Covarianza in R (÷n-1)
cov(Temp, Gas)
#> [1] -2.194
codev / (n - 1)   # verifica
#> [1] -2.194
# Correlazione
cor(Temp, Gas)
#> [1] -0.6832545

Esercizio 3: Si carichi il data frame fev, incluso nel file fev.csv. Con riferimento alle variabili fev (forced expiratory volume, misura di capacità polmonare) e height (altezza), ci aspettiamo ,senza nessuno strumento grafico o indice, una relazione positiva/negativa tra queste due variabili? Quindi si produca un diagramma di dispersione e si calcolino la covarianza e la correlazione tra queste due variabili. I risultati confermano l’intuizione? Dall’analisi del diagramma di dispersione ci aspettiamo di ottenere una correlazione maggiore o minore tra fev e height se considero solo gli individui più alti di 65 pollici. Si verifichi la propria risposta empiricamente.

5 La Regressione Lineare Semplice

5.1 La Funzione di Regressione

💡 Quando le variabili hanno ruoli diversi (\(Y\) risposta, \(X\) esplicativa), l’obiettivo è descrivere come varia la media condizionata \(\mathbb{E}(Y|x)\) al variare di \(x\). Questa funzione è detta funzione di regressione.

⚠️ Con dati reali non è possibile calcolare \(\mathbb{E}(Y|x)\) direttamente (non si osservano più valori di \(Y\) per lo stesso \(x\)). La strategia più semplice è assumere che la funzione di regressione sia lineare:

\[\mathbb{E}(Y|x) = \beta_0 + \beta_1 x\] dove \(\beta_0\) e \(\beta_1\) sono i due parametri che definiscono la retta (ignoti, quindi da stimare)

5.2 I Minimi Quadrati

💡 Il criterio dei minimi quadrati sceglie \(\hat{\beta}_0\) e \(\hat{\beta}_1\) che minimizzano la somma dei quadrati dei residui \(r_i = y_i - \hat{y}_i\):

\[\min_{\beta_0, \beta_1} L(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2\]

Le soluzioni possono essere reperite facilmente considerando la soluzione del seguente sistema di equazioni (equazioni normali): \[\begin{align} \frac{\partial L(\beta_0,\beta_1)}{\partial \beta_0} &= -2\sum_{i=1}^n (y_i-\beta_0-\beta_1 x_i)=0\\ \frac{\partial L(\beta_0,\beta_1)}{\partial \beta_1} &= -2 \sum_{i=1}^n (y_i-\beta_0-\beta_1 x_i)x_i= 0 \end{align}\]

La soluzione del sistema delle equazioni normali dà:

\[\hat{\beta}_1 = \frac{cov(X,Y)}{var(X)} = r(X,Y)\frac{sd(Y)}{sd(X)}, \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\]

Per la prima equazione normale vale anche \(\sum_{i=1}^n r_i = 0\).

b1 <- cov(Gas, Temp) / var(Temp)
b0 <- mean(Gas) - b1 * mean(Temp)
print(paste("Intercetta:", round(b0, 4)))
#> [1] "Intercetta: 5.4862"
print(paste("Coefficiente angolare:", round(b1, 4)))
#> [1] "Coefficiente angolare: -0.2902"
plot(Temp, Gas, pch = 19,
     xlab = "Temperatura (°C)", ylab = "Gas (mc)",
     main = "Retta dei minimi quadrati")
abline(b0, b1, col = 2, lwd = 2)

Otteniamo anche i valori predetti (sulla retta) e la somma dei quadrati dei residui.

Gasprev <- b0 + b1 * Temp
res      <- Gas - Gasprev
SQE      <- sum(res^2)
print(paste("Devianza residui (SQE):", round(SQE, 4)))
#> [1] "Devianza residui (SQE): 39.9949"
plot(Temp, Gas, pch = 19,
     xlab = "Temperatura (°C)", ylab = "Gas (mc)",
     main = "Retta dei minimi quadrati")
segments(Temp[3], Gas[3], Temp[3], Gasprev[3], col = "deepskyblue3", lty=2)
points(Temp[3], Gasprev[3], pch = 19, col = "red")
abline(b0, b1, col = 2, lwd = 2)

⚠️ Qualsiasi altra retta produce una SQE maggiore:

# Retta basata sui quartili
b1q <- sign(cov(Gas, Temp))*(quantile(Gas, 0.75) - quantile(Gas, 0.25)) /
         (quantile(Temp, 0.75) - quantile(Temp, 0.25))
b0q <- mean(Gas) - b1q * mean(Temp)

Gasprevq <- b0q + b1q * Temp
SQE2     <- sum((Gas - Gasprevq)^2)
print(paste("SQE minimi quadrati:", round(SQE,  4)))
#> [1] "SQE minimi quadrati: 39.9949"
print(paste("SQE retta quartili: ", round(SQE2, 4)))
#> [1] "SQE retta quartili:  40.0779"
plot(Temp, Gas, pch = 19,
     xlab = "Temperatura", ylab = "Gas",
     main = "Minimi quadrati (rosso) vs retta quartili (verde)")
abline(b0,  b1,  col = 2, lwd = 2)
abline(b0q, b1q, col = 3, lwd = 2)
legend("topright", legend = c("Minimi quadrati", "Quartili"),
       col = c(2, 3), lwd = 2)

5.3 Il Coefficiente di Determinazione \(R^2\)

Per misurare la qualità della regressione si può utilizzare il coefficiente di determinazione lineare \(R^2\).

La devianza totale di \(Y\) si scompone in:

\[\underbrace{\sum(y_i - \bar{y})^2}_{SQT} = \underbrace{\sum(y_i - \hat{y}_i)^2}_{SQE} + \underbrace{\sum(\hat{y}_i - \bar{y})^2}_{SQR}\]

  • SQT: devianza totale - variabilità totale
  • SQE: devianza dei residui- variabilità dell’errore
  • SQR: devianza di regressione - variabilità spiegata dal modello

💡 Il coefficiente di determinazione misura la quota di variabilità spiegata dalla regressione:

\[R^2 = \frac{SQR}{SQT} = 1 - \frac{SQE}{SQT} \in [0,1]\]

Si può dimostrare che \(R^2 = r(X,Y)^2\) ma solo nel caso di regressione lineare semplice (ovvero con una singola variabile esplicativa).

R2 <- 1 - (SQE/(var(Gas)*(length(Gas)-1)))
R2
#> [1] 0.4668366
cor(Gas, Temp)^2
#> [1] 0.4668366

Esercizio 4: Si carichi il data frame gatti, incluso nel file gatti.csv. Con riferimento alle variabili peso (in kg) e cibo (in g). Si produca un diagramma di dispersione e si calcoli il coefficiente di correlazione. Si calcolino i coefficienti della retta di regressione dei minimi quadrati e la si aggiunga al grafico. Si ottenga il coefficiente di determinazione \(R^2\).

5.4 Il Modello Statistico e la Funzione lm()

💡 Il modello lineare specifica:

\[Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma^2) \text{ i.i.d.}\]

Con questa assunzione si può fare inferenza sui parametri (t-test, intervalli di confidenza). La funzione lm() fornisce tutte le informazioni rilevanti:

mod1 <- lm(Gas ~ Temp, data = whiteside)
summary(mod1)
#> 
#> Call:
#> lm(formula = Gas ~ Temp, data = whiteside)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -1.6324 -0.7119 -0.2047  0.8187  1.5327 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   5.4862     0.2357  23.275  < 2e-16 ***
#> Temp         -0.2902     0.0422  -6.876 6.55e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.8606 on 54 degrees of freedom
#> Multiple R-squared:  0.4668, Adjusted R-squared:  0.457 
#> F-statistic: 47.28 on 1 and 54 DF,  p-value: 6.545e-09

L’output riporta: stime dei coefficienti con errori standard e p-value, \(R^2\) e \(R^2\) corretto, statistica \(F\) per il test globale del modello.

Soffermiamoci sull’interpretazione di alcune quantità:

  • Intercetta (\(\hat \beta_0 \approx 5.49\)): consumo medio stimato di gas quando la temperatura esterna è 0°C. Poiché 0°C rientra nell’intervallo dei valori osservati, questo valore è direttamente interpretabile.
  • Coefficiente angolare (\(\hat \beta_1 \approx -0.29\)): per ogni grado Celsius aggiuntivo di temperatura esterna, il consumo atteso di gas diminuisce di circa 0.29 m³. Il segno negativo conferma che il freddo porta a un maggiore consumo.
  • \(R^2\approx0.47\): circa il 47% della variabilità nel consumo di gas è spiegata dalla temperatura.
detach(whiteside)

Esercizio 5: Riprendendo l’esercizio precedente, si ottengano le stesse quantità richieste utilizzando il comando lm ed interpretare i coefficienti e il valore dell’\(R^2\). Infine, qual è il peso atteso per un gatto che mangia 100 grammi di crocchetti al giorno?

5.5 Lisciamento della Curva di Regressione

Assumere la linearità può essere restrittivo. Un approccio alternativo è stimare la curva di regressione non parametricamente, lasciando che i dati suggeriscano la forma. L’idea è calcolare in corrispondenza di ogni \(x\) una media locale ponderata:

\[m(x) = \frac{\sum_{i=1}^n y_i K\!\left(\frac{x-x_i}{h}\right)}{\sum_{i=1}^n K\!\left(\frac{x-x_i}{h}\right)}\]

dove \(K(\cdot)\) è una funzione nucleo (tipicamente gaussiana) e \(h\) è il parametro di lisciamento (bandwidth). A differenza della regressione lineare, non si impone alcuna forma funzionale a priori: la curva viene costruita punto per punto, combinando le osservazioni vicine a ciascun \(x\).

5.6 ksmooth() e lowess(): due approcci a confronto

Le due funzioni principali in R implementano strategie diverse per il lisciamento locale.

5.6.1 ksmooth() — Medie locali con nucleo

ksmooth() implementa direttamente la formula della media locale ponderata. Per ogni punto \(x\) di valutazione:

  1. assegna a ciascuna osservazione \((x_i, y_i)\) un peso proporzionale a \(K\!\left(\frac{x - x_i}{h}\right)\), dove \(K\) è una gaussiana (default) o una funzione box;
  2. calcola la media pesata dei valori \(y_i\).

Il parametro di controllo è bandwidth, espresso nella stessa unità di misura di \(x\): rappresenta l’ampiezza della finestra intorno a ogni punto di valutazione. Osservazioni a distanza maggiore di bandwidth/2 ricevono peso zero (per il nucleo box) o peso trascurabile (per il nucleo gaussiano).

5.6.2 lowess() — Regressione locale robusta

lowess() (Locally Weighted Scatterplot Smoothing) adotta un approccio più sofisticato: per ogni punto \(x_i\) seleziona i vicini più prossimi, assegna loro pesi decrescenti con la distanza, adatta una regressione lineare locale pesata (non una semplice media) e applica un ciclo di ri-pesatura robusta che riduce l’influenza dei valori anomali.

Il parametro di controllo è f — detto anche span — che rappresenta la proporzione del numero totale di osservazioni usate per ogni stima locale, e non un’ampiezza in unità di \(x\) come bandwidth di ksmooth(). Valori tipici sono compresi tra 0.2 e 0.8.

5.6.3 Esempio

Consideriamo il dataframe case da caricare da file esterno. Il dataset contiene, tra le altre variabili, il prezzo delle case vendute in un certo periodo ad Albuquerque (New Mexico) e la loro dimensione (in piedi quadrati). Che fra le due variabili ci sia una relazione è scontato, che il prezzo medio delle case vari linearmente al variare della dimensione lo è un po’ meno.

case <- read.csv("case.csv")
plot(case$SquareFeet, case$Price, pch = 19, cex = 0.5,
     col = "gray50",
     main = "Regressione: lineare vs lisciamento",
     xlab = "Superficie (piedi²)", ylab = "Prezzo ($)")

modl <- lm(Price ~ SquareFeet, data = case)
abline(modl, col = 2, lwd = 2)

# ksmooth(): medie locali con nucleo gaussiano
lines(ksmooth(case$SquareFeet, case$Price, bandwidth = 1400),
      col = 4, lty = 2, lwd = 2)

# lowess(): regressione locale robusta, f = proporzione di punti usati
lines(lowess(case$Price ~ case$SquareFeet, f = 0.8),
      col = 5, lwd = 2)

legend("topleft",
       legend = c("Minimi quadrati", "ksmooth (bw=1400)", "lowess (f=0.8)"),
       col = c(2, 4, 5), lty = c(1, 2, 1), lwd = 2)

📌 Le due curve di lisciamento rivelano una relazione non perfettamente lineare tra superficie e prezzo: nella fascia delle abitazioni più grandi la crescita del prezzo tende ad appiattirsi, un dettaglio che la retta dei minimi quadrati non cattura.

5.6.4 Il Parametro di Lisciamento \(h\)

Il parametro \(h\) è il principale strumento di controllo della curva stimata:

💡 Effetto di \(h\):

  • \(h\) piccolo → la finestra locale è stretta, si usano pochi punti vicini a \(x\). La curva segue da vicino i dati ma risulta molto irregolare e sensibile al rumore campionario.
  • \(h\) grande → la finestra locale è ampia, si mediano molte osservazioni anche lontane da \(x\). La curva è più liscia ma può perdere strutture locali reali, appiattendo picchi e curve.

⚠️ Non esiste un valore di \(h\) universalmente ottimale: la scelta dipende dalla struttura dei dati e dallo scopo dell’analisi. In pratica si prova più valori e si valuta il risultato graficamente.

La figura seguente illustra il trade-off al variare di \(h\):

par(mfrow = c(1, 3))
for (bw in c(300, 1400, 5000)) {
  plot(case$SquareFeet, case$Price, pch = 19, cex = 0.4,
       col = "gray50",
       main = paste("bandwidth =", bw),
       xlab = "Superficie (piedi²)", ylab = "Prezzo ($)")
  lines(ksmooth(case$SquareFeet, case$Price, bandwidth = bw),
        col = 4, lwd = 2)
}

par(mfrow = c(1, 1))

📌 Con bandwidth = 300 la curva è frastagliata e sovra-adattata; con bandwidth = 5000 è quasi piatta e perde la struttura locale; bandwidth = 1400 rappresenta un compromesso ragionevole.

5.6.5 Curve di lisciamento con ggplot()

Con ggplot2 le curve smoothed si ottengono con geom_smooth():

library(ggplot2)

ggplot(case, aes(x = SquareFeet, y = Price)) +
  geom_point(alpha = 0.4, color = "gray50") +
  geom_smooth(se = FALSE, color = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, color = "firebrick", linetype = "dashed") +
  labs(title = "Prezzo vs Superficie: loess (blu) vs lineare (rosso tratteggiato)",
       x = "Superficie (piedi²)", y = "Prezzo ($)")

ggplot(iris, aes(x = Sepal.Length, y = Petal.Length)) +
  geom_point(aes(colour = Species), alpha = 0.7) +
  geom_smooth(aes(colour = Species), se = FALSE) +
  labs(title = "Curva loess per specie",
       x = "Lunghezza sepalo (cm)", y = "Lunghezza petalo (cm)")

💡 geom_smooth() usa loess per default (\(n < 1000\)). Con method = "lm" forza la retta dei minimi quadrati; con se = FALSE nasconde la banda di confidenza. Il parametro di span si può controllare con span = (equivalente a f di lowess()).

Esercizio 6 Si utilizzi il dataset gatti già caricato in precedenza. L’obiettivo è esplorare la relazione tra peso e cibo confrontando la retta dei minimi quadrati con due curve di lisciamento.

  1. Si sovrapponga al diagramma di dispersione: la retta dei minimi quadrati (in rosso), una curva ksmooth() con bandwidth = 0.5 e una con bandwidth = 2.5. Come cambia la curva al variare del parametro? Quale dei due valori produce un lisciamento eccessivo? Quale sotto-adatta?

  2. Si sovrapponga invece una curva lowess() con f = 0.3 e una con f = 0.9. Si commenti la differenza tra le due curve e si confronti il loro comportamento con le curve ksmooth() del punto precedente.

  3. Dato che la relazione tra peso e consumo di cibo è per costruzione quasi lineare, ci si aspetta che le curve di lisciamento si discostino molto dalla retta dei minimi quadrati? Si verifichi graficamente e si commenti.

6 Il Controllo del Modello lineare: l’Analisi dei Residui

I residui \(r_i = y_i - \hat{y}_i\) forniscono informazioni preziose sulla bontà del modello. Se il modello è adeguato, essi dovrebbero:

  • essere vicini a 0 e distribuirsi in modo casuale attorno a 0 (la loro somma è sempre esattamente 0)
  • non mostrare pattern sistematici

Il grafico principale è il plot dei residui vs valori predetti (\(r_i\) contro \(\hat{y}_i\)): i residui dovrebbero distribuirsi casualmente attorno allo zero senza struttura.

Pattern preoccupanti:

  • forma a imbuto (variabilità crescente) → variabilità non constante (eteroschedasticità)
  • andamento curvo → non linearità non catturata dal modello
  • due nuvole separate → presenza gruppi che non si stanno considerando

Altri grafici utili:

  • \(r_i\) contro \(x_j\) (per ciascun predittore)
  • QQ-plot dei residui (per valutare la normalità)
  • residui standardizzati, leverage e distanza di Cook (per identificare osservazioni anomale o influenti)

Consideriamo l’esempio precedente su temperatura (predittore) e consumo di gas (variabile risposta) e analizziamo i residui del modello lineare semplice. In R, le funzioni per ottenere i residui e i valori stimati di un oggetto lm sono, rispettivamente, residuals() e fitted()

res  <- residuals(mod1)
yhat <- fitted(mod1)

plot(yhat, res, pch = 20, col = 4,
     xlab = expression(hat(y)), ylab = "residui",
     main = "Residui vs valori stimati — modello semplice")
abline(h = 0, lty = 2)

📌 Il grafico dei residui segnala due problemi: (i) la variabilità non è costante e (ii) sembrano esistere due gruppi di residui (positivi e negativi che divergono all’aumentare di \(\hat{y}\). Questo potebbe essere dovuto a una variabile omessa.

Il secondo grafico di interesse è il QQ-plot dei residui. I residui dovrebbero distribuirsi secondo una normale.

qqnorm(residuals(mod1))              # normalità dei residui
qqline(residuals(mod1))

> 📌 Il grafico evidenzia uno scarso adattamento alla distribuzione Normale.

6.1 Visualizzazione dei residui in R

💡 La funzione plot applicata a un oggetto di tipo lm produce quattro grafici diagnostici dei residui

  1. Valori predetti (\(\hat y_i\)) vs residui (\(r_i\))
  2. QQ-plot: grafico quantile quantile
  3. Valori predetti vs (radice) residui studentizzati
  4. Leverage vs residui studentizzati
par(mfrow=c(2,2))
plot(mod1)

par(mfrow=c(1,1))

📌 I grafici ottenuti a mano in precedenza corrispondono ai pannelli in alto. Gli altri due pannelli sono informativi di altri aspetti che il modello potrebbe non cogliere ma non verranno ulteriormente approfonditi.