61  La regressione lineare bivariata: un approccio frequentista

In questo capitolo imparerai a:

• comprendere il funzionamento del modello lineare secondo l’approccio frequentista;
• stimare i coefficienti del modello utilizzando il metodo dei minimi quadrati e interpretarli correttamente;
• valutare la qualità del modello attraverso l’indice di determinazione (\(R^2\)).

Prerequisiti

• Leggere il capitolo Basic Regression di Statistical Inference via Data Science: A ModernDive into R and the Tidyverse (Second Edition).
• Consulta l’appendice Appendice J per un’introduzione alle funzioni lineari.
• Leggere Navigating the Bayes maze: The psychologist’s guide to Bayesian statistics, a hands-on tutorial with R code (Alter et al., 2025).
• Leggere il capitolo Linear Statistical Models (Schervish & DeGroot, 2014).

Preparazione del Notebook

here::here("code", "_common.R") |> 
  source()

# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(broom)

61.1 Introduzione

La regressione è una tecnica statistica che permette ai ricercatori di analizzare come i valori medi di una variabile risultato (o variabile dipendente) cambiano in relazione a un insieme di predittori (variabili indipendenti).

Secondo Gelman et al. (2021), i principali utilizzi della regressione includono:

1. Previsione – Modellare dati esistenti o prevedere nuove osservazioni, sia continue che categoriche.
   • Esempi: Prevedere punteggi futuri in un test, monitorare il benessere psicologico in uno studio longitudinale o classificare individui in base alla probabilità di successo in un compito cognitivo.
2. Esplorazione delle associazioni – Identificare e quantificare le relazioni tra una o più variabili indipendenti e un risultato.
   • Esempi: Analizzare i tratti di personalità legati alla resilienza allo stress, indagare la correlazione tra stili di attaccamento infantile e capacità relazionali in età adulta, o valutare l’influenza di fattori socio-economici sullo sviluppo cognitivo nei bambini.
3. Estrapolazione – Estendere i risultati osservati in un campione a una popolazione più ampia.
   • Esempi: Stimare l’efficacia di una terapia testata su studenti universitari alla popolazione generale, oppure prevedere l’impatto di un intervento scolastico su un intero distretto partendo dai dati di alcune scuole pilota.
4. Inferenza causale – Analizzare gli effetti di un trattamento o intervento, purché supportata da un disegno sperimentale adeguato.
   • Esempi: Valutare l’efficacia di un programma di mindfulness sui livelli di ansia, stimare l’impatto di una psicoterapia per il disturbo post-traumatico da stress o determinare l’effetto di un intervento educativo su una popolazione diversificata.

In tutti questi contesti, è essenziale che il modello includa tutte le variabili rilevanti per la variabile di interesse. Ad esempio, in uno studio sull’efficacia di una psicoterapia per la depressione, fattori come età, condizioni di salute preesistenti e supporto sociale possono influenzare significativamente i risultati e devono essere integrati nell’analisi. La trascuratezza di tali variabili può distorto le relazioni stimate tra predittori e risultato. L’accuratezza dei modelli di regressione dipende quindi sia dalla completezza del modello che dalla selezione mirata delle informazioni incluse.

Un problema critico legato a questa scelta è il cosiddetto “errore di specificazione”, che si verifica quando uno o più predittori rilevanti vengono omessi dal modello. Questo porta a stime sistematicamente errate, compromettendo la validità delle conclusioni tratte. Approfondiremo questo concetto in seguito, esaminando le implicazioni pratiche e strategie per evitarlo.

61.1.1 La storia dei modelli lineari

I modelli lineari hanno una lunga storia nella statistica. Come riportato da Stigler (1986), il metodo dei minimi quadrati per adattare un modello di regressione lineare bivariata fu introdotto nel XVIII secolo per l’analisi dei dati astronomici. Gli astronomi lo utilizzavano, ad esempio, per determinare il moto della Luna o per modellare i movimenti non periodici di Giove e Saturno.

L’adozione di questi metodi fu favorita dalla qualità e dall’omogeneità delle misurazioni astronomiche, raccolte con strumenti standardizzati e sotto condizioni controllate. Al contrario, nelle scienze sociali, l’elevata variabilità dei dati—derivante da differenze individuali, contesti mutevoli e strumenti di misura meno precisi—rese più complessa l’applicazione della regressione e ne ritardò la diffusione. Solo con lo sviluppo di tecniche statistiche più avanzate, capaci di gestire questa complessità, la regressione divenne uno strumento essenziale anche nelle scienze umane e sociali.

61.1.2 Tipologie di regressione

È utile distinguere tra tre principali categorie di regressione:

1. Regressione bivariata – Un solo predittore e una sola variabile dipendente.
   
2. Regressione multipla – Un solo risultato ma molteplici predittori.
   
3. Regressione multivariata – Molteplici risultati, con uno o più predittori.

Il caso bivariato, pur essendo il più semplice, è raramente sufficiente per descrivere fenomeni complessi, soprattutto in psicologia, dove molteplici fattori influenzano le variabili di interesse. Tuttavia, lo analizzeremo in dettaglio perché permette di comprendere la logica dell’analisi di regressione in un contesto semplificato, privo di complicazioni matematiche. I concetti appresi saranno poi estesi ai modelli più complessi, in cui la regressione multipla e la regressione multivariata richiederanno calcoli più articolati e un’interpretazione più sofisticata dei parametri.

61.2 Il modello lineare bivariato

Nel contesto frequentista, il modello di regressione lineare bivariata consente di predire una variabile continua \(y\) sulla base di un singolo predittore continuo \(x\). La relazione tra le due variabili è espressa dall’equazione della retta di regressione:

\[ y_i = a + b x_i + e_i, \quad i = 1, \dots, n \]

dove:

• \(a\) è l’intercetta (il valore atteso di \(y\) quando \(x = 0\)),
• \(b\) è la pendenza della retta (coefficiente di regressione, che misura il cambiamento atteso in \(y\) per ogni unità di incremento in \(x\)),
• \(e_i\) è l’errore residuo (la differenza tra il valore osservato di \(y_i\) e il valore predetto dal modello).

61.2.1 Aspetti principali

1. Stima dei coefficienti
   I coefficienti \(a\) e \(b\) vengono stimati mediante il metodo dei minimi quadrati, che minimizza la somma dei quadrati degli errori residui (\(\sum e_i^2\)).

2. Interpretazione dei coefficienti

   • \(a\): rappresenta il valore medio previsto di \(y\) quando \(x = 0\).
     
   • \(b\): indica la variazione media prevista in \(y\) per ogni unità di variazione in \(x\).

3. Valutazione del modello
   La bontà di adattamento del modello viene valutata attraverso:

   • L’indice di determinazione (\(R^2\)), che misura la proporzione della varianza di \(y\) spiegata dal modello.
   • L’analisi dei residui, utile per individuare pattern non catturati dal modello e diagnosticare eventuali problemi.

61.2.2 Verso modelli più complessi

Questo capitolo illustrerà come applicare e interpretare il modello di regressione bivariata, fornendo una base per comprendere il modello lineare multiplo e, più in generale, gli approcci avanzati per l’inferenza causale. La logica della regressione bivariata sarà estesa ai modelli con più predittori, dove sarà fondamentale distinguere tra relazioni spurie e reali effetti predittivi.

61.3 La Predizione dell’Intelligenza

Nella presente discussione, esamineremo i dati kidiq che consistono in una raccolta di dati provenienti da un sondaggio su donne adulte americane e i loro figli, selezionati da un sotto-campione del National Longitudinal Survey of Youth (Gelman et al., 2021).

Nello specifico, ci concentreremo sulla relazione tra il punteggio di intelligenza del bambino (kid_score) e quello della madre (mom_iq). Ci proponiamo di valutare se e in quale misura l’intelligenza della madre possa prevedere l’intelligenza del bambino. Per fare ciò, inizieremo ad importare i dati nell’ambiente R.

# Caricamento dei dati
kidiq <- rio::import(here::here("data", "kidiq.dta"))

# Anteprima dei dati
head(kidiq)
#>   kid_score mom_hs mom_iq mom_work mom_age
#> 1        65      1 121.12        4      27
#> 2        98      1  89.36        4      25
#> 3        85      1 115.44        4      27
#> 4        83      1  99.45        3      25
#> 5       115      1  92.75        4      27
#> 6        98      0 107.90        1      18

Un diagramma a dispersione per i dati di questo campione suggerisce la presenza di un’associazione positiva tra l’intelligenza del bambino (kid_score) e l’intelligenza della madre (mom_iq).

ggplot(kidiq, aes(x = mom_iq, y = kid_score)) +
  geom_point(alpha = 0.6) +
  labs(x = "QI della madre", y = "QI del bambino") +
  ggtitle("Diagramma a dispersione")

61.4 Stima del modello di regressione lineare

Calcoliamo i coefficienti della retta di regressione utilizzando la funzione lm.

# Modello di regressione lineare
mod <- lm(kid_score ~ mom_iq, data = kidiq)

# Coefficienti stimati
coef(mod)
#> (Intercept)      mom_iq 
#>       25.80        0.61

Ci sono però infinite rette che, in linea di principio, possono essere usate per “approssimare” la nube di punti nel diagramma a dispersione. È dunque necessario introdurre dei vincoli per selezionare una di queste possibili rette. Il vincolo che viene introdotto dal modello di regressione è quello di costringere la retta a passare per il punto \((\bar{x}, \bar{y})\).

# Calcola le medie per mom_iq e kid_score
mean_x <- mean(kidiq$mom_iq, na.rm = TRUE)
mean_y <- mean(kidiq$kid_score, na.rm = TRUE)

# Aggiungi il punto medio al grafico
ggplot(kidiq, aes(x = mom_iq, y = kid_score)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  annotate(
    "point", x = mean_x, y = mean_y, color = "red", size = 5, 
    shape = 4, stroke = 3) +
  labs(x = "QI della madre", y = "QI del bambino") +
  ggtitle("Retta di regressione")

Una retta di regressione che passa per il punto medio \((\bar{x}, \bar{y})\) (che rappresenta il centro di massa dei dati) è preferibile dal punto di vista statistico poiché minimizza la somma dei quadrati degli errori residui.

Il campione è costituito da \(n\) coppie di osservazioni (\(x, y\)).

\[ \begin{array}{cc} \hline x_1 & y_1 \\ x_2 & y_2 \\ x_3 & y_3 \\ \vdots & \vdots \\ x_n & y_n \\ \hline \end{array} \]

Per ciascuna coppia di valori \(x_i, y_i\), il modello di regressione si aspetta che il valore \(y_i\) sia associato al corrispondente valore \(x_i\) come indicato dalla seguente equazione

\[ \begin{equation} \mathbb{E}(y_i) = a + b x_i , \end{equation} \tag{61.1}\]

ovvero:

\[ \begin{array}{ccc} \hline x_i & y_i & \mathbb{E}(y_i) = a + b x_i \\ \hline x_1 & y_1 & a + b x_1 \\ x_2 & y_2 & a + b x_2 \\ x_3 & y_3 & a + b x_3 \\ \vdots & \vdots & \vdots \\ x_n & y_n & a + b x_n \\ \hline \end{array} \]

I valori \(y_i\) corrispondono, nell’esempio che stiamo discutendo, alla variabile kid_score. I primi 10 valori della variabile \(y\) sono i seguenti:

kidiq$kid_score[0:10]
#>  [1]  65  98  85  83 115  98  69 106 102  95

Per fare riferimento a ciascuna osservazione usiamo l’indice \(i\). Quindi, ad esempio, \(y_2\) è uguale a

kidiq$kid_score[2]
#> [1] 98

Il modello di regressione lineare bivariata, rappresentato dall’equazione \(y_i = a + b x_i + e_i\), descrive la relazione tra le variabili \(x\) e \(y\), dove \(y\) è la variabile dipendente (nel nostro esempio, la variabile kid_score) e \(x\) è la variabile indipendente (nel nostro esempio, la variabile mom_iq).

Il valore di \(y\) è la somma di due componenti:

• la componente deterministica, \(\hat{y}_i = a + b x_i\), rappresenta la porzione della \(y\) che è prevedibile conoscendo il valore di \(x\);
• la componente aleatoria, \(e_i\), rappresenta la porzione della \(y\) che non è prevedibile dal modello.

Il modello lineare cerca di trovare i coefficienti \(a\) e \(b\) che permettono di prevedere la componente deterministica di \(y\) conoscendo il valore di \(x\). Tuttavia, poiché la retta è solo un’approssimazione della relazione tra \(x\) e \(y\), la componente deterministica rappresenta solo una stima approssimata della vera relazione tra le due variabili.

Per valutare l’accuratezza del modello di regressione lineare, è necessario calcolare il residuo

\[ e_i = y_i - (a + b x_i) , \tag{61.2}\]

ovvero la differenza tra il valore osservato di \(y\) e il valore previsto dal modello, \(\hat{y}\). La dimensione del residuo indica quanto la componente aleatoria contribuisce al valore osservato di \(y\).

Il modello di regressione lineare ha tre obiettivi (Fox, 2015):

1. il primo è quello di trovare i coefficienti \(a\) e \(b\) che permettono di prevedere la componente deterministica di \(y\) conoscendo il valore di \(x\);
2. il secondo obiettivo è quello di valutare l’accuratezza della predizione fornita dal modello di regressione lineare;
3. infine, il terzo obiettivo è quello dell’inferenza, ovvero quello di capire quali relazioni esistono tra la relazione tra \(x\) e \(y\) osservata nel campione e la relazione tra le due variabili nella popolazione.

61.5 Stima dei coefficienti di regressione

In breve, stiamo cercando di descrivere una relazione tra due variabili, il QI della madre e il QI del bambino, utilizzando un modello di regressione lineare. L’equazione lineare che descrive la relazione tra le due variabili è della forma \(\hat{y}_i = a_i + bx_i\), dove \(\hat{y}_i\) rappresenta la previsione per il QI del bambino \(i\)-esimo, \(a_i\) e \(b\) sono i coefficienti di regressione che vogliamo trovare e \(x_i\) è il QI della madre del bambino \(i\)-esimo.

Per trovare i coefficienti di regressione, dobbiamo introdurre dei vincoli per limitare lo spazio delle possibili soluzioni. Il primo vincolo è che la retta di regressione deve passare per il baricentro del grafico a dispersione. Il secondo vincolo è che vogliamo minimizzare la somma dei quadrati dei residui, ovvero la differenza tra il valore osservato e il valore previsto dal modello. I coefficienti di regressione che soddisfano questi vincoli si chiamano coefficienti dei minimi quadrati.

Il problema di trovare i coefficienti di regressione \(a\) e \(b\) che minimizzano la somma dei quadrati dei residui ha una soluzione analitica. Questa soluzione si ottiene trovando il punto di minimo di una superficie tridimensionale che rappresenta la somma dei quadrati dei residui. Il punto di minimo è quello per cui il piano tangente alla superficie nelle due direzioni \(a\) e \(b\) è piatto, cioè le derivate parziali rispetto ad \(a\) e \(b\) sono uguali a zero. In pratica, ciò significa risolvere un sistema di equazioni lineari con due incognite \(a\) e \(b\), noto come equazioni normali.

La soluzione delle equazioni normali ci fornisce i coefficienti di regressione stimati, che minimizzano la somma dei quadrati dei residui. La formula per il coefficiente \(a\) è

\[ a = \bar{y} - b \bar{x} , \tag{61.3}\]

la formula per il coefficiente \(b\) è

\[ b = \frac{Cov(x, y)}{Var(x)}, \tag{61.4}\]

dove \(\bar{x}\) e \(\bar{y}\) sono le medie delle variabili \(x\) e \(y\), \(Cov(x,y)\) è la covarianza tra \(x\) e \(y\) e \(Var(x)\) è la varianza di \(x\).

Queste equazioni rappresentano la stima dei minimi quadrati dei coefficienti di regressione che ci permettono di trovare la retta che minimizza la somma dei quadrati dei residui.

61.5.1 Calcolo manuale dei coefficienti di regressione

Calcoliamo i coefficientii dei minimi quadrati con l’Equazione 61.3 e l’Equazione 61.4:

# Calcolo manuale dei coefficienti
cov_xy <- cov(kidiq$kid_score, kidiq$mom_iq)  # Covarianza tra le due variabili
var_x <- var(kidiq$mom_iq)  # Varianza della variabile indipendente
b <- cov_xy / var_x  # Pendenza (coefficiente b)
b
#> [1] 0.61

# Intercetta (coefficiente a)
a <- mean(kidiq$kid_score) - b * mean(kidiq$mom_iq)
a
#> [1] 25.8

I risultati replicano quelli ottenuti in precedenza con lm().

61.5.2 Interpretazione

Il coefficiente \(a\) indica l’intercetta della retta di regressione nel diagramma a dispersione. Questo valore rappresenta il punto in cui la retta di regressione interseca l’asse \(y\) del sistema di assi cartesiani. Tuttavia, in questo caso specifico, il valore di \(a\) non è di particolare interesse poiché corrisponde al valore della retta di regressione quando l’intelligenza della madre è pari a 0, il che non ha senso nella situazione reale. Successivamente, vedremo come è possibile trasformare i dati per fornire un’interpretazione utile del coefficiente \(a\).

Invece, il coefficiente \(b\) indica la pendenza della retta di regressione, ovvero di quanto aumenta (se \(b\) è positivo) o diminuisce (se \(b\) è negativo) la retta di regressione in corrispondenza di un aumento di 1 punto della variabile \(x\). Nel caso specifico del QI delle madri e dei loro figli, il coefficiente \(b\) ci indica che un aumento di 1 punto del QI delle madri è associato, in media, a un aumento di 0.61 punti del QI dei loro figli.

In pratica, il modello di regressione lineare cerca di prevedere le medie dei punteggi del QI dei figli in base al QI delle madri. Ciò significa che non è in grado di prevedere esattamente il punteggio di ciascun bambino in funzione del QI della madre, ma solo una stima della media dei punteggi dei figli quando il QI delle madri aumenta o diminuisce di un punto.

Il coefficiente \(b\) ci dice di quanto aumenta (o diminuisce) in media il QI dei figli per ogni unità di aumento (o diminuzione) del QI della madre. Nel nostro caso, se il QI della madre aumenta di un punto, il QI dei figli aumenta in media di 0.61 punti.

È importante comprendere che il modello statistico di regressione lineare non è in grado di prevedere il valore preciso di ogni singolo bambino, ma solo una stima della media dei punteggi del QI dei figli quando il QI delle madri aumenta o diminuisce. Questa stima è basata su una distribuzione di valori possibili che si chiama distribuzione condizionata \(p(y \mid x_i)\).

Una rappresentazione grafica del valore predetto dal modello di regressione, \(\hat{y}_i = a + bx_i\) è stato fornito in precedenza. Il diagramma presenta ciascun valore \(\hat{y}_i = a + b x_i\) in funzione di \(x_i\). I valori predetti dal modello di regressione sono i punti che stanno sulla retta di regressione.

61.6 Residui

In precedenza abbiamo detto che il residuo, ovvero la componente di ciascuna osservazione \(y_i\) che non viene predetta dal modello di regressione, corrisponde alla distanza verticale tra il valore \(y_i\) osservato e il valore \(\hat{y}_i\) predetto dal modello di regressione:

\[ e_i = y_i - (a + b x_i). \]

Per fare un esempio numerico, consideriamo il punteggio osservato del QI del primo bambino.

kidiq$kid_score[1]
#> [1] 65

Il QI della madre è

kidiq$mom_iq[1]
#> [1] 121.1

Per questo bambino, il valore predetto dal modello di regressione è

a + b * kidiq$mom_iq[1]
#> [1] 99.68

L’errore che compiamo per predire il QI del bambino utilizzando il modello di regressione (ovvero, il residuo) è

kidiq$kid_score[1] - (a + b * kidiq$mom_iq[1])
#> [1] -34.68

Per tutte le osservazioni abbiamo

res <- kidiq$kid_score - (a + b * kidiq$mom_iq)

È una proprietà del modello di regressione (calcolato con il metodo dei minimi quadrati) che la somma dei residui sia uguale a zero.

sum(res)
#> [1] 1.444e-11

Questo significa che ogni valore osservato \(y_i\) viene scomposto dal modello di regressione in due componenti distinte. La componente deterministica \(\hat{y}_i\), che è predicibile da \(x_i\), è data da \(\hat{y}_i = a + b x_i\). Il residuo, invece, è dato da \(e_i = y_i - \hat{y}_i\). La somma di queste due componenti, ovviamente, riproduce il valore osservato.

# Creazione di un data frame con i valori calcolati
df <- data.frame(
  kid_score = kidiq$kid_score,
  mom_iq = kidiq$mom_iq,
  y_hat = a + b * kidiq$mom_iq,
  e = kidiq$kid_score - (a + b * kidiq$mom_iq),
  y_hat_plus_e = (a + b * kidiq$mom_iq) + (kidiq$kid_score - (a + b * kidiq$mom_iq))
)

# Visualizzazione dei primi 6 valori
head(df)
#>   kid_score mom_iq y_hat       e y_hat_plus_e
#> 1        65 121.12 99.68 -34.678           65
#> 2        98  89.36 80.31  17.692           98
#> 3        85 115.44 96.22 -11.217           85
#> 4        83  99.45 86.46  -3.462           83
#> 5       115  92.75 82.37  32.628          115
#> 6        98 107.90 91.62   6.383           98

61.7 Trasformazione dei dati

In generale, per variabili a livello di scala ad intervalli, l’intercetta del modello di regressione lineare non ha un’interpretazione utile. Questo perché l’intercetta indica il valore atteso di \(y\) quando \(x = 0\), ma in caso di variabili a scala di intervalli, il valore “0” di \(x\) è arbitrario e non corrisponde ad un “assenza” della variabile \(x\). Ad esempio, un QI della madre pari a 0 non indica un’assenza di intelligenza, ma solo un valore arbitrario del test usato per misurare il QI. Quindi, sapere il valore medio del QI dei bambini quando il QI della madre è 0 non è di alcun interesse.

Per fornire all’intercetta del modello di regressione un’interpretazione più utile, dobbiamo trasformare le osservazioni di \(x\). Per esempio, esprimiamo \(x\) come differenza dalla media. Chiamiamo questa nuova variabile \(xd\):

kidiq$xd <- kidiq$mom_iq - mean(kidiq$mom_iq)

kidiq |> 
  head()
#>   kid_score mom_hs mom_iq mom_work mom_age       xd
#> 1        65      1 121.12        4      27  21.1175
#> 2        98      1  89.36        4      25 -10.6381
#> 3        85      1 115.44        4      27  15.4432
#> 4        83      1  99.45        3      25  -0.5504
#> 5       115      1  92.75        4      27  -7.2543
#> 6        98      0 107.90        1      18   7.9018

Se ora usiamo le coppie di osservazioni \((xd_i, y_i)\), il diagramma a dispersione assume la forma seguente.

# Aggiungiamo una nuova variabile centrata (scarti dalla media)
kidiq <- kidiq %>%
  mutate(xd = mom_iq - mean(mom_iq))

# Calcolo della retta di regressione
b <- cov(kidiq$xd, kidiq$kid_score) / var(kidiq$xd)
a <- mean(kidiq$kid_score) - b * mean(kidiq$xd)

# Grafico con ggplot2
ggplot(kidiq, aes(x = xd, y = kid_score)) +
  geom_point(alpha = 0.6) +  # Punti del grafico
  geom_abline(intercept = a, slope = b, color = "blue") +  # Retta di regressione
  labs(
    x = "QI della madre (scarti dalla media)", 
    y = "QI del bambino"
  ) +
  ggtitle("Retta di regressione sui dati centrati")

In pratica, abbiamo spostato tutti i punti del grafico lungo l’asse delle \(x\), in modo tale che la media dei valori di \(x\) sia uguale a 0. Questo non ha cambiato la forma dei punti nel grafico, ma ha solo spostato l’origine dell’asse \(x\). La pendenza della linea di regressione tra \(x\) e \(y\) rimane la stessa, sia per i dati originali che per quelli trasformati. L’unica cosa che cambia è il valore dell’intercetta della linea di regressione, che ora ha un’interpretazione più significativa.

fm1 <- lm(kid_score ~ xd, data = kidiq)
coef(fm1)
#> (Intercept)          xd 
#>       86.80        0.61

L’intercetta rappresenta il punto in cui la retta di regressione incontra l’asse \(y\) nel diagramma a dispersione. Nel caso dei dati trasformati, abbiamo spostato la nube di punti lungo l’asse \(x\) di una quantità pari a \(x - \bar{x}\), ma le relazioni spaziali tra i punti rimangono invariate. Pertanto, la pendenza della retta di regressione non cambia rispetto ai dati non trasformati. Tuttavia, il valore dell’intercetta viene influenzato dalla trasformazione. In particolare, poiché \(xd = 0\) corrisponde a \(x = \bar{x}\) nei dati grezzi, l’intercetta del modello di regressione lineare calcolata sui dati trasformati corrisponde al valore atteso di \(y\) quando \(x\) assume il valore medio sulla scala dei dati grezzi. In altre parole, l’intercetta del modello di regressione lineare sui dati trasformati rappresenta il valore atteso del QI dei bambini corrispondente al QI medio delle madri.

61.8 Il metodo dei minimi quadrati

Per stimare i coefficienti \(a\) e \(b\), possiamo minimizzare la somma dei quadrati dei residui tra i valori osservati \(y_i\) e quelli previsti \(a + b x_i\).

Iniziamo con il creare una griglia per i valori di \(b\). Supponiamo che il valore di \(a\) sia noto (\(a = 25.79978\)). Usiamo R per creare una griglia di valori possibili per \(b\).

# Griglia di valori per b
b_grid <- seq(0, 1, length.out = 1001)
a <- 25.79978  # Intercetta nota

Definiamo ora una funzione che calcola la somma dei quadrati dei residui (\(SSE\)) per ciascun valore di \(b\).

# Funzione per la somma dei quadrati dei residui
sse <- function(a, b, x, y) {
  sum((y - (a + b * x))^2)
}

Applichiamo la funzione sse alla griglia di valori \(b\) per calcolare la somma dei quadrati dei residui per ogni valore di \(b\).

# Calcolo di SSE per ciascun valore di b
sse_vals <- sapply(
  b_grid, 
  function(b) sse(a, b, kidiq$mom_iq, kidiq$kid_score)
)

• sapply:
  • È una funzione di R che applica una funzione ad ogni elemento di un vettore (o lista) e restituisce i risultati in un vettore.
  • Qui, applica la funzione sse a ciascun valore di \(b\) contenuto in b_grid.
• function(b):
  • È una funzione anonima definita al volo per specificare come calcolare \(SSE\) per ciascun valore di \(b\).
  • All’interno, viene chiamata la funzione sse(a, b, x, y) con i seguenti parametri:
    • a: il valore dell’intercetta (fissato in precedenza o noto).
    • b: il valore corrente nella griglia b_grid.
    • x: la variabile indipendente del dataset (kidiq$mom_iq).
    • y: la variabile dipendente del dataset (kidiq$kid_score).
• Il risultato è un vettore, sse_vals, che contiene i valori di \(SSE\) corrispondenti a ciascun valore di \(b\) in b_grid.

Tracciamo un grafico che mostra la somma dei quadrati dei residui (\(SSE\)) in funzione dei valori di \(b\), evidenziando il minimo.

# Identificazione del valore di b che minimizza SSE
b_min <- b_grid[which.min(sse_vals)]

# Creazione del dataframe per ggplot
dat <- data.frame(b_grid = b_grid, sse_vals = sse_vals)

# Genera il grafico
ggplot(dat, aes(x = b_grid, y = sse_vals)) +
  geom_line(color = "blue", linewidth = 1) +  
  annotate(
    "point", x = b_min, y = min(sse_vals),
    color = "red", size = 3
  ) +  # Punto minimo
  labs(
    x = expression(paste("Possibili valori di ", hat(beta))),
    y = "Somma dei quadrati dei residui (SSE)",
    title = "Minimizzazione dei residui quadratici"
  ) +
  annotate(
    "text", x = b_min, y = min(sse_vals), 
    label = expression(hat(beta)), color = "red", vjust = -1, hjust = 0.5
  )

Infine, identifichiamo il valore di \(b\) che minimizza la somma dei quadrati dei residui.

b_min
#> [1] 0.61

Con questa simulazione, abbiamo stimato il coefficiente \(b\) minimizzando la somma dei quadrati dei residui.

Questo approccio può essere esteso per stimare simultaneamente entrambi i coefficienti (\(a\) e \(b\)) utilizzando metodi di ottimizzazione più avanzati, come optim in R.

optim_result <- optim(
  par = c(a = 25, b = 0.5),  # Valori iniziali
  fn = function(params) {
    a <- params[1]
    b <- params[2]
    sse(a, b, kidiq$mom_iq, kidiq$kid_score)
  }
)

# Coefficienti stimati
optim_result$par
#>       a       b 
#> 25.7874  0.6101

Questa simulazione illustra come, tramite il metodo dei minimi quadrati, sia possibile stimare i parametri di un modello bivariato di regressione.

61.9 L’errore standard della regressione

Il secondo obiettivo del modello di regressione lineare è quello di misurare quanto della variabilità di \(y\) possa essere spiegata dalla variabilità di \(x\) per ogni osservazione. L’indice di bontà di adattamento del modello viene fornito dalla deviazione standard dei residui, chiamata anche errore standard della stima (o errore standard della regressione), \(s_e\).

Per calcolare \(s_e\), si sommano i quadrati dei residui \(e_i\) per ogni osservazione e si divide per \(n-2\), dove \(n\) rappresenta la numerosità del campione e \(2\) il numero di coefficienti stimati nel modello di regressione. Si prende poi la radice quadrata del risultato:

\[ \sqrt{\frac{1}{n-2} \sum_{i=1}^n \big(y_i - (\hat{a} + \hat{b}x_i)\big)^2}, \tag{61.5}\]

L’indice \(s_e\) possiede la stessa unità di misura di \(y\) ed è una stima della deviazione standard dei residui nella popolazione.

Illustriamo il calcolo di \(s_e\) con i dati a disposizione. I residui \(e\) possono essere calcolati sottraendo ai valori osservati \(y_i\) i valori predetti dal modello \(a + b x_i\):

# Calcolo dei residui
e <- kidiq$kid_score - (a + b * kidiq$mom_iq)

# Mostriamo i primi 10 residui
head(e, 10)
#>  [1] -34.678  17.692 -11.217  -3.462  32.628   6.383 -41.521   3.865  26.414
#> [10]  11.208

Calcoliamo il valore medio assoluto dei residui per avere un’indicazione della deviazione media rispetto alla retta di regressione.

# Media assoluta dei residui
mean(abs(e))
#> [1] 14.47

L’errore standard della stima \(s_e\) si calcola come la radice quadrata della somma dei quadrati dei residui divisa per \(n-2\):

# Calcolo di s_e
se <- sqrt(sum(e^2) / (length(e) - 2))
se
#> [1] 18.27

Notiamo che il valore medio assoluto dei residui e l’errore standard \(s_e\) non sono identici, ma hanno lo stesso ordine di grandezza. \(s_e\) è una misura più rigorosa della deviazione standard dei residui.

Questa analisi dimostra come \(s_e\) consenta di valutare quanto le previsioni del modello si discostino (in media) dai dati osservati.

61.9.1 Sottostima dell’Errore nel Modello di Regressione

Come discusso da Gelman et al. (2021), l’errore standard della regressione tende a sottostimare la vera deviazione standard \(\sigma\) dell’errore nel modello di regressione nella popolazione. Questa sottostima è dovuta al fenomeno del sovradimensionamento, dato che i parametri \(a\) e \(b\) sono stimati utilizzando gli stessi \(n\) punti dati su cui vengono calcolati i residui. In altre parole, i residui non sono del tutto indipendenti dal modello.

Un approccio alternativo per valutare l’errore predittivo e mitigare il problema del sovradimensionamento è la validazione incrociata. In particolare, l’approccio leave-one-out (LOOCV) offre una soluzione semplice ed efficace. Questo metodo consiste nell’adattare il modello \(n\) volte, escludendo ogni volta un punto dati, adattando il modello ai rimanenti \(n-1\) punti, e utilizzando tale modello per predire l’osservazione esclusa.

61.9.1.1 Procedura Leave-One-Out:

• Per \(i = 1, \ldots, n\):

  1. Adatta il modello \(y = a + bx + \text{errore}\) ai \(n-1\) punti dati \((x, y)_j, j \neq i\). Denomina i coefficienti stimati come \(\hat{a}_{-i}\) e \(\hat{b}_{-i}\).

  2. Calcola il residuo validato incrociato:

     \[ e_{\text{CV}} = y_i - (\hat{a}_{-i} + \hat{b}_{-i} x_i). \]

  3. Salva il residuo al quadrato per il calcolo successivo.

• Calcola infine la stima di \(\sigma_{\text{CV}}\) come:

  \[ \sigma_{\text{CV}} = \sqrt{\frac{1}{n} \sum_{i=1}^n e_{\text{CV}}^2}. \]

61.9.1.2 Applicazione Pratica:

Ecco un esempio applicato al modello che predice l’intelligenza del bambino (\(\texttt{kid\_score}\)) in funzione dell’intelligenza della madre (\(\texttt{mom\_iq}\)) utilizzando il dataset kidiq.

options(round = 5)

# Array per salvare i residui validati incrociati
residuals_cv <- numeric(nrow(kidiq))

# Loop per la validazione incrociata leave-one-out
for (i in 1:nrow(kidiq)) {
  # Dati di training escludendo l'i-esimo punto
  train_data <- kidiq[-i, ]
  test_data <- kidiq[i, ]
  
  # Addestramento del modello
  model <- lm(kid_score ~ mom_iq, data = train_data)
  
  # Predizione sull'i-esimo punto
  y_pred <- predict(model, newdata = test_data)
  
  # Calcolo del residuo validato incrociato
  residual_cv <- test_data$kid_score - y_pred
  residuals_cv[i] <- residual_cv^2
}

# Calcolo di sigma_cv
sigma_cv <- sqrt(mean(residuals_cv))

cat("Stima di σ_CV:", sigma_cv, "\n")
#> Stima di σ_CV: 18.31

Calcoliamo ora la stima tradizionale di \(\sigma\) utilizzando il modello completo e confrontiamola con \(\sigma_{\text{CV}}\).

# Modello completo
fm2 <- lm(kid_score ~ mom_iq, data = kidiq)

# Stima tradizionale dell'errore standard della regressione
res <- summary(fm2)
cat("Errore standard della regressione (tradizionale):", res$sigma, "\n")
#> Errore standard della regressione (tradizionale): 18.27

Nel caso analizzato, i valori stimati di \(\sigma_{\text{CV}}\) e \(\hat{\sigma}_e\) tradizionale possono risultare molto simili. Tuttavia, in generale, la stima di \(\sigma_{\text{CV}}\) tende a essere leggermente superiore, in quanto riflette meglio l’errore predittivo su dati non utilizzati per adattare il modello. Questo rende \(\sigma_{\text{CV}}\) una misura più robusta e conservativa dell’incertezza del modello.

In conclusione, la validazione incrociata, e in particolare l’approccio LOOCV, rappresenta uno strumento importante per valutare le performance predittive di un modello di regressione e per ottenere stime più affidabili della deviazione standard dell’errore.

61.10 Indice di determinazione

Un importante risultato dell’analisi di regressione riguarda la scomposizione della varianza della variabile dipendente \(y\) in due componenti: la varianza spiegata dal modello e la varianza residua. Questa scomposizione è descritta mediante l’indice di determinazione \(R^2\), che fornisce una misura della bontà di adattamento del modello ai dati del campione.

Per una generica osservazione \(x_i, y_i\), la deviazione di \(y_i\) rispetto alla media \(\bar{y}\) può essere espressa come la somma di due componenti: il residuo \(e_i=y_i- \hat{y}_i\) e lo scarto di \(\hat{y}_i\) rispetto alla media \(\bar{y}\):

\[ y_i - \bar{y} = (y_i- \hat{y}_i) + (\hat{y}_i - \bar{y}) = e_i + (\hat{y}_i - \bar{y}). \]

La varianza totale di \(y\) può quindi essere scritta come:

\[ \sum_{i=1}^{n}(y_i - \bar{y})^2 = \sum_{i=1}^{n}(e_i + (\hat{y}_i - \bar{y}))^2. \]

Sviluppando il quadrato e sommando, si ottiene:

\[ \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. \tag{61.6}\]

Il primo termine rappresenta la varianza residua, mentre il secondo termine rappresenta la varianza spiegata dal modello. Questa scomposizione della devianza va sotto il nome di teorema della scomposizione della devianza.

L’indice di determinazione \(R^2\) è definito come il rapporto tra la varianza spiegata e la varianza totale:

\[ R^2 = \frac{\sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2}. \tag{61.7}\]

Questo indice varia tra 0 e 1 e indica la frazione di varianza totale di \(y\) spiegata dal modello di regressione lineare. Un valore alto di \(R^2\) indica che il modello di regressione lineare si adatta bene ai dati, in quanto una grande parte della varianza di \(y\) è spiegata dalla variabile indipendente \(x\).

Per l’esempio in discussione, possiamo calcolare la devianza totale, la devianza spiegata e l’indice di determinazione \(R^2\) come segue:

La devianza totale misura la variabilità complessiva dei punteggi osservati \(y\) rispetto alla loro media:

# Devianza totale
dev_t <- sum((kidiq$kid_score - mean(kidiq$kid_score))^2)
dev_t
#> [1] 180386

La devianza spiegata misura la variabilità che il modello è in grado di spiegare, considerando i valori predetti \(a + b x\):

# Devianza spiegata
dev_r <- sum(((a + b * kidiq$mom_iq) - mean(kidiq$kid_score))^2)
dev_r
#> [1] 36249

L’indice \(R^2\) è il rapporto tra la devianza spiegata e la devianza totale, e indica la frazione della variabilità totale che è spiegata dal modello di regressione:

# Indice di determinazione
R2 <- dev_r / dev_t
round(R2, 3)
#> [1] 0.201

Per verificare i calcoli, utilizziamo il modello di regressione lineare in R e leggiamo \(R^2\) direttamente dal sommario del modello:

# Modello di regressione lineare
mod <- lm(kid_score ~ mom_iq, data = kidiq)

# Sommario del modello per leggere R^2
summary(mod)$r.squared
#> [1] 0.201

Il risultato mostra che circa il 20% della variabilità nei punteggi del QI dei bambini è spiegabile conoscendo il QI delle madri. Questo significa che il modello cattura una porzione rilevante della relazione, ma lascia anche spazio a fattori non inclusi nel modello che influenzano il QI dei bambini.

61.11 Inferenza sul modello di regressione

L’inferenza statistica sui coefficienti di regressione richiede la definizione della distribuzione campionaria dei coefficienti di regressione. Il modello di regressione bivariata (o lineare semplice) è:

\[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i,\]

dove \(Y_i\) è la variabile dipendente per l’osservazione \(i\), \(X_i\) è la variabile indipendente, \(\beta_0\) è l’intercetta (parametro ignoto), \(\beta_1\) è il coefficiente angolare (il parametro che ci interessa stimare), e \(\varepsilon_i\) è il termine di errore per l’osservazione \(i\).

Nell’inferenza, il nostro obiettivo è stimare i parametri ignoti \(\beta_0\) e \(\beta_1\) usando i dati campionari disponibili. Il metodo più comune è quello dei Minimi Quadrati Ordinari (OLS), che ci fornisce gli stimatori \(\hat{\beta}_0\) e \(\hat{\beta}_1\) (spesso indicati semplicemente con \(b_0\) e \(b_1\)). Lo stimatore per il coefficiente angolare è dato dalla formula:

\[\hat{\beta}_1 = b_1 = \frac{\sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i=1}^n (X_i - \bar{X})^2}.\]

61.12 Cos’è la Distribuzione Campionaria di \(b_1\) ?

Lo stimatore \(b_1\) è una variabile casuale. Perché? Perché il suo valore dipende dal campione casuale di dati \((X_i, Y_i)\) che abbiamo estratto dalla popolazione.

Immaginiamo di poter ripetere il processo di campionamento e stima del modello infinite volte:

1. Estraiamo un campione casuale di dimensione \(n\).
2. Calcoliamo lo stimatore \(b_1\) per quel campione.
3. Registriamo il valore di \(b_1\).
4. Estraiamo un nuovo campione casuale (indipendentemente dal primo).
5. Calcoliamo il “nuovo” \(b_1\).
6. Registriamo il valore.
7. … e così via, per infinite volte.

La distribuzione campionaria di \(b_1\) è la distribuzione di probabilità di tutti i valori di \(b_1\) che otterremmo da questi infiniti campioni casuali di dimensione \(n\) estratti dalla stessa popolazione.

61.13 Proprietà della Distribuzione Campionaria di \(b_1\)

Il modello statistico standard della regressione bivariata assume che:

• 1. Linearità nei parametri: La relazione tra la variabile dipendente (\(Y\)) e la variabile indipendente (\(X\)) è lineare nei parametri (\(\beta_0\) e \(\beta_1\)). Questo implica che il valore atteso di \(Y_i\) condizionato a \(X_i\), indicato con \(E(Y_i \mid X_i)\), sia espresso correttamente come \(E(Y_i) = \beta_0 + \beta_1 X_i\).

  • È importante notare che questa assunzione non richiede che i singoli punti dati giacciano esattamente su una retta, ma che la funzione di regressione della popolazione (cioè la relazione tra la media di \(Y\) e \(X\)) sia lineare. L’impatto medio di una variazione unitaria in \(X\) su \(Y\), quantificato attraverso il parametro \(\beta_1\), è additivo (contribuisce in modo indipendente da altre variabili) e costante (non dipende dal livello di \(X\)). Ad esempio, un aumento di 1 unità in \(X\) è associato sempre a un aumento di \(\beta_1\) unità in \(Y\), indipendentemente dal valore iniziale di \(X\).

• 2. Campionamento Casuale: I dati osservati \((X_i, Y_i)\) per \(i = 1, \dots, n\) rappresentano un campione casuale estratto dalla popolazione di interesse. Questo implica che le osservazioni sono indipendenti l’una dall’altra e provengono dalla stessa distribuzione congiunta di \((X, Y)\).

• 3. Esogeneità (o Media Condizionata degli Errori Nulla): Il valore atteso del termine di errore (\(\epsilon_i\)) è zero, condizionatamente al valore della variabile indipendente (\(X_i\)). Formalmente: \(E(\epsilon_i \mid X_i) = 0\). Questa è un’assunzione cruciale. Implica che non ci siano altri fattori (contenuti nell’errore \(\epsilon_i\)) che sono correlati con la variabile indipendente \(X_i\) e che influenzano anche \(Y_i\).

  • Se questa assunzione è violata (ad esempio, a causa di variabili omesse correlate sia con \(X\) che con \(Y\)), lo stimatore OLS \(b_1\) sarà distorto e inconsistente.

• 4. Omoschedasticità (Varianza Costante degli Errori): La varianza del termine di errore (\(\epsilon_i\)) è costante per tutti i valori della variabile indipendente (\(X_i\)). Formalmente: \(Var(\epsilon_i \mid X_i) = \sigma^2\), dove \(\sigma^2\) è una costante sconosciuta e finita.

  • Questa assunzione implica che la dispersione degli errori attorno alla retta di regressione è la stessa lungo tutto il range di valori di \(X\). Se questa assunzione è violata (eteroschedasticità), lo stimatore OLS \(b_1\) rimane non distorto, ma i suoi errori standard stimati non saranno corretti, invalidando i test di ipotesi e gli intervalli di confidenza basati su di essi.

• 5. Assenza di Collinearità Perfetta: Nel modello bivariato, questa assunzione si riduce semplicemente al fatto che la variabile indipendente \(X_i\) non deve essere una costante per tutte le osservazioni nel campione. Ovvero, \(\sum_{i=1}^n (X_i - \bar{X})^2 > 0\).

  • Se \(X_i\) fosse costante, non ci sarebbe variazione in \(X\) da cui stimare l’effetto sulla variazione in \(Y\), e la formula per \(b_1\) implicherebbe una divisione per zero.

• 6. Normalità degli Errori (Opzionale per piccoli campioni): Il termine di errore (\(\epsilon_i\)) è distribuito normalmente con media zero e varianza costante \(\sigma^2\). Formalmente: \(\epsilon_i \sim N(0, \sigma^2)\).

  • Questa assunzione non è necessaria perché gli stimatori OLS \(b_0\) e \(b_1\) siano non distorti e per le proprietà di grandi campioni (come la consistenza e la normalità asintotica). Tuttavia, è cruciale per derivare le distribuzioni di probabilità esatte (come la distribuzione t di Student per i test sui coefficienti) degli stimatori nei campioni di piccole dimensioni, permettendo un’inferenza valida in questi casi. Per grandi campioni, la distribuzione asintotica (approssimata) normale degli stimatori vale anche senza questa assunzione.

Sotto le precedenti assunzioni, chiamate assunzioni di Gauss-Markov, più l’assunzione di normalità degli errori per l’inferenza nei piccoli campioni, la distribuzione campionaria di \(b_1\) viene specificata come segue:

1. Media (o Valore Atteso): Unbiasedness (Non Distorsione) Sotto le prime tre assunzioni di Gauss-Markov (Linearità, Campionamento Casuale, Esogeneità dei regressori), lo stimatore OLS \(b_1\) è non distorto (unbiased) per il vero parametro di popolazione \(\beta_1\). Questo significa che la media della distribuzione campionaria di \(b_1\) è uguale al vero valore di \(\beta_1\):

   \[E(b_1) = \beta_1.\]

   In altre parole, se potessimo prendere infiniti campioni, la media degli stimatori \(b_1\) che otterremmo sarebbe esattamente il vero coefficiente angolare della popolazione. Il nostro singolo stimatore \(b_1\) calcolato su un campione specifico è solo una realizzazione di questa variabile casuale la cui distribuzione è centrata sul vero valore \(\beta_1\).

2. Varianza (o Errore Standard): Precisione dello Stimatore La varianza della distribuzione campionaria di \(b_1\) ci dice quanto sono dispersi i valori di \(b_1\) attorno alla loro media (che è \(\beta_1\)). Una varianza minore indica che gli stimatori ottenuti da campioni diversi tendono ad essere più vicini al vero \(\beta_1\), il che significa che lo stimatore è più preciso. Sotto le assunzioni di Gauss–Markov (inclusa l’omoschedasticità, ossia errori con varianza costante \(\sigma^2_\varepsilon\)), la varianza di \(b_1\) è data da:

   \[ \mathrm{Var}(b_1) \;=\; \frac{\sigma^2_\varepsilon}{\sum_{i=1}^n (X_i - \bar{X})^2}. \]

   L’errore standard di \(b_1\), radice quadrata della varianza, è la misura più comunemente usata per indicare la precisione dello stimatore:

   \[ SE(b_1) \;=\; \sqrt{\frac{\sigma^2_\varepsilon}{\sum_{i=1}^n (X_i - \bar{X})^2}}. \]

   Dalla formula, possiamo osservare che la precisione di \(b_1\):

   • diminuisce all’aumentare della varianza degli errori (\(\sigma^2_\varepsilon\)): più “rumore” c’è nei dati, meno preciso sarà lo stimatore;
   • aumenta all’aumentare della variabilità della variabile indipendente \(X\) (misurata da \(\sum_{i=1}^n (X_i - \bar{X})^2\)): se i valori di \(X\) sono molto concentrati, è difficile stimare la pendenza; se sono ben distribuiti, la stima è più precisa;
   • aumenta all’aumentare della dimensione del campione (\(n\)): campioni più grandi forniscono più informazioni e quindi stime più precise.

   Nella pratica, non conosciamo la vera varianza degli errori \(\sigma^2_\varepsilon\), perciò la stimiamo usando la varianza dei residui del modello (\(s^2_e\)), calcolata come:

   \[ s^2_e \;=\; \frac{\sum_{i=1}^n e_i^2}{n - 2}, \]

   dove \(n - 2\) riflette i gradi di libertà persi nella stima dei due parametri (\(b_0\) e \(b_1\)). L’errore standard stimato diventa quindi:

   \[ SE(b_1) \;=\; \sqrt{\frac{s^2_e}{\sum_{i=1}^n (X_i - \bar{X})^2}}. \]

   Questa è la formula riportata nei principali software statistici (R, Stata, Python).

3. Forma della Distribuzione:

   • Sotto le assunzioni di Gauss-Markov (ma senza assumere la normalità degli errori), per grandi campioni (in virtù del Teorema del Limite Centrale), la distribuzione campionaria di \(b_1\) è approssimativamente Normale.
   • Se aggiungiamo l’assunzione che i termini di errore \(\varepsilon_i\) si distribuiscano normalmente, \(\varepsilon_i \sim N(0, \sigma^2)\), allora la distribuzione campionaria di \(b_1\) è esattamente Normale per qualsiasi dimensione del campione \(n\). Tuttavia, quando usiamo l’errore standard stimato (che dipende dalla stima della varianza degli errori), la statistica test (t-statistic) segue una distribuzione t di Student.

Importanza della Distribuzione Campionaria di \(b_1\)

Conoscere la distribuzione campionaria di \(b_1\) permette all’approccio frequentista di:

• Testare Ipotesi: Possiamo testare se la variabile indipendente \(X\) ha un effetto statisticamente significativo sulla variabile dipendente \(Y\). L’ipotesi nulla più comune è \(H_0: \beta_1 = 0\) (non c’è relazione lineare tra \(X\) e \(Y\)). Usiamo il valore stimato \(b_1\) e il suo errore standard per calcolare una statistica test (il \(t\)-value) e confrontarla con una distribuzione teorica (Normale o \(t\) di Student) per ottenere un \(p\)-value.
• Costruire Intervalli di Confidenza: Possiamo costruire un intervallo di valori plausibili per il vero parametro di popolazione \(\beta_1\) basato sul nostro campione. Un intervallo di confidenza al 95% per \(\beta_1\) significa che se ripetessimo il campionamento e la stima molte volte, il 95% degli intervalli costruiti in questo modo conterrebbe il vero valore di \(\beta_1\).

In sintesi, la distribuzione campionaria di \(b_1\) descrive la variabilità attesa dello stimatore del coefficiente angolare attraverso diversi campioni casuali. Comprendere le sue proprietà (media, varianza, forma) è essenziale per interpretare correttamente i risultati di una regressione e trarre conclusioni affidabili sulla relazione nella popolazione.¹

61.14 Simulazione in R

Modello Simulato

Studiamo il modello di regressione lineare senza intercetta, per variabili \(X\) e \(Y\) standardizzate:

\[ Y_i = \beta X_i + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma_\varepsilon = 0.5). \] Fissiamo \(\beta = 1.5\), \(n = 30\) (dimensione del campione) e ripetiamo la simulazione \(100,\!000\) volte per approssimare la distribuzione campionaria di \(\hat{\beta}_1\).

# Parametri del modello
beta <- 1.5           # Vero valore del coefficiente
sigma_e <- 0.5        # Deviazione standard degli errori
n_sample <- 30        # Dimensione del campione per ogni simulazione
nrep <- 100000        # Numero di repliche (campioni) da generare

# Preallocazione di un vettore per salvare gli stimatori
b_hat <- rep(NA, nrep)

# Genera la variabile indipendente X una volta (fissa tra i campioni)
# Questo corrisponde a condizionare sulle X, come nelle assunzioni di Gauss-Markov
x <- rnorm(n_sample, mean = 0, sd = 1)

# Ciclo per replicare il processo di campionamento
for (i in 1:nrep) {
  # 1. Genera errori casuali ε_i ~ N(0, σ_ε)
  e <- rnorm(n_sample, mean = 0, sd = sigma_e)
  
  # 2. Genera la variabile dipendente Y_i = β * X_i + ε_i
  y <- beta * x + e
  
  # 3. Stima β_1 usando la formula OLS: Cov(X,Y)/Var(X)
  # Questa formula è valida per modelli con intercetta, ma qui X è centrata (media ≈ 0)
  b_hat[i] <- cov(x, y) / var(x)
}

# Analisi dei risultati
hist(
  b_hat, 
  breaks = 50, 
  col = "lightblue", 
  main = expression("Distribuzione Campionaria di" ~ hat(beta)[1]), 
  xlab = expression("Valore Stimato" ~ hat(beta)[1]),
  freq = FALSE
)

# Sovrappone una curva normale con media teorica e deviazione standard osservata
curve(dnorm(x, mean = mean(b_hat), sd = sd(b_hat)), add = TRUE, col = "red", lwd = 2)

# Media empirica di β̂₁ (dovrebbe essere vicina a β = 1.5)
mean_b_hat <- mean(b_hat)
cat("Media empirica di β̂₁:", round(mean_b_hat, 3), "\n")
#> Media empirica di β̂₁: 1.5

# Deviazione standard empirica di β̂₁ (confronto con teoria)
sd_b_hat <- sd(b_hat)
cat("Deviazione standard empirica di β̂₁:", round(sd_b_hat, 3), "\n")
#> Deviazione standard empirica di β̂₁: 0.074

# Calcolo dell'errore standard teorico
# Var(β̂₁) = σ_ε² / Σ(X_i - X̄)²
theoretical_se <- sqrt(sigma_e^2 / sum((x - mean(x))^2))
cat("Errore standard teorico:", round(theoretical_se, 3), "\n")
#> Errore standard teorico: 0.074

Interpretazione dei Risultati

1. Non Distorsione (Unbiasedness)

• La media empirica di \(\hat{\beta}_1\) (calcolata come mean(b_hat)) dovrebbe essere molto vicina al valore vero \(\beta = 1.5\).
  
• Teoria: Sotto le assunzioni di Gauss-Markov (linearità, esogeneità, campionamento casuale), lo stimatore OLS è non distorto: \(E(\hat{\beta}_1) = \beta_1\).
  
• Simulazione: Se la media è prossima a \(1.5\), confermiamo questa proprietà.

2. Precisione (Varianza)

• La deviazione standard empirica (sd(b_hat)) misura la dispersione dello stimatore attorno al valore vero.
  
• Teoria: La varianza di \(\hat{\beta}_1\) dipende da:
  • La varianza degli errori \(\sigma_\varepsilon^2\) (maggiore rumore → minore precisione),
  • La variabilità di \(X\) (\(\sum (X_i - \bar{X})^2\); maggiore variabilità → maggiore precisione).
• Simulazione: Confrontiamo sd(b_hat) con l’errore standard teorico calcolato come: \[ SE(\hat{\beta}_1) = \sqrt{\frac{\sigma_\varepsilon^2}{\sum_{i=1}^n (X_i - \bar{X})^2}}. \] I valori dovrebbero essere simili.

3. Normalità Asintotica

• L’istogramma mostra la distribuzione degli stimatori \(\hat{\beta}_1\) ottenuti nei \(100,\!000\) campioni.
  
• Teoria: Per grandi campioni (o se gli errori sono normali), la distribuzione campionaria di \(\hat{\beta}_1\) è normale.
  
• Simulazione: La curva rossa (densità normale) dovrebbe sovrapporsi bene all’istogramma, confermando la forma normale.

Collegamento con le Assunzioni del Modello

1. Esogeneità (\(E(\varepsilon_i \mid X_i) = 0\)):
   Nella simulazione, gli errori \(\varepsilon_i\) sono generati indipendentemente da \(X_i\), soddisfacendo questa assunzione. Se invece \(\varepsilon_i\) fosse correlato con \(X_i\) (es. \(X_i\) influenzasse \(\varepsilon_i\)), lo stimatore sarebbe distorto.

2. Omoschedasticità:
   La varianza degli errori è costante (\(\sigma_\varepsilon = 0.5\)), garantendo che la formula dell’errore standard sia corretta. In presenza di eteroschedasticità, l’errore standard teorico non sarebbe più affidabile.

3. Normalità degli Errori:
   Gli errori sono generati da una distribuzione normale, ma la normalità asintotica della distribuzione campionaria di \(\hat{\beta}_1\) sarebbe comunque approssimativamente valida anche con errori non normali (per il Teorema del Limite Centrale).

Conclusione

La simulazione conferma le proprietà teoriche della distribuzione campionaria di \(\hat{\beta}_1\):

• non distorsione: la media empirica coincide con il valore vero;
• precisione: la varianza dello stimatore diminuisce con la variabilità di \(X\) e aumenta con la dimensione del campione;
• normalità: la distribuzione è simmetrica e approssimativamente normale, anche per \(n = 30\).

Questo esempio pratico chiarisce come le assunzioni del modello influenzino la qualità delle stime e l’affidabilità dell’inferenza statistica.

61.15 Riflessioni conclusive 

Il modello lineare bivariato costituisce uno strumento basilare per analizzare la relazione tra due variabili e rappresenta una pietra miliare dell’approccio frequentista. In questo capitolo abbiamo mostrato come il modello permetta di quantificare l’associazione fra una variabile indipendente \(X\) e una variabile dipendente \(Y\) tramite una relazione lineare.

Nell’ambito frequentista, i coefficienti \(b_0\) (intercetta) e \(b_1\) (pendenza) vengono stimati con il metodo dei minimi quadrati, che minimizza la somma dei quadrati dei residui e individua la retta che meglio approssima i dati osservati. L’indice di determinazione (\(R^2\)) misura la bontà di adattamento del modello, indicando quale quota della variabilità di \(Y\) è spiegata da \(X\).

In pratica, il modello lineare bivariato consente di rispondere a domande fondamentali quali:

• Come varia la variabile dipendente al variare della variabile indipendente?
• Qual è intensità e segno dell’associazione tra le due variabili?

La semplicità del modello lo rende utile non solo per descrivere e analizzare relazioni, ma anche per formulare previsioni. Pur limitandosi a un singolo predittore, costituisce la base concettuale per modelli più complessi (ad es. la regressione multipla).

Derivare la distribuzione campionaria dei coefficienti permette di svolgere inferenza sui parametri. Tuttavia, nel paradigma frequentista l’inferenza si fonda:

1. sul test dell’ipotesi nulla (che assume, nella popolazione, pendenza esattamente pari a zero);
2. sulla costruzione degli intervalli di fiducia, i quali descrivono le proprietà a lungo termine della procedura di campionamento, ma non forniscono probabilità per il singolo intervallo.

La riflessione metodologica contemporanea ha evidenziato i limiti di questa impostazione, mostrando come l’inferenza frequentista sia vulnerabile a pratiche scorrette che hanno alimentato la crisi di riproducibilità in psicologia. Per tali motivi, risulta spesso più informativo riformulare l’inferenza in termini bayesiani, dove si ottengono distribuzioni a posteriori dei parametri che consentono interpretazioni probabilistiche dirette e maggior trasparenza nel processo inferenziale.

Problemi

1 – Definizione e scopi della regressione
Secondo Gelman et al. (2021), quali sono i quattro principali utilizzi della regressione? Fornisci una breve descrizione di ciascuno.

2 – Errore di specificazione
Cos’è l’“errore di specificazione” in un modello di regressione? Quali effetti ha sulle stime dei parametri?

3 – Stima del modello con lm()
Usando il data‑set kidiq (variabili kid_score e mom_iq):

# carica i dati
kidiq <- rio::import(here::here("data", "kidiq.dta"))

Adatta il modello kid_score ~ mom_iq e riporta intercetta e pendenza.

4 – Interpretazione della pendenza
Interpreta in parole la pendenza stimata al punto 3 nel contesto dei QI di madri e figli.

5 – Indice R²
Calcola l’R² del modello di cui al punto 3. Cosa indica il suo valore?

6 – Centratura del predittore
Crea la variabile centrata mom_iq_c = mom_iq - mean(mom_iq) e ri‑adatta il modello kid_score ~ mom_iq_c. Qual è la nuova intercetta e perché adesso è più interpretabile?

7 – Calcolo manuale di \(b\)
Calcola manualmente la pendenza con la formula
\(b = \frac{\text{Cov}(x,y)}{\text{Var}(x)}\)
e confrontala col risultato di lm().

8 – Confronto tra σ̂ (tradizionale) e σ_CV
* (a) Calcola l’errore standard della regressione (σ̂) usando il modello completo.
* (b) Esegui una validazione incrociata leave‑one‑out (LOOCV) e ottieni σ_CV.
* (c) Spiega perché, in genere, σ_CV è ≥ σ̂.

9 – Assunzioni di Gauss‑Markov
Elenca le cinque assunzioni di Gauss‑Markov per la regressione lineare semplice e indica quale, se violata, rende distorto lo stimatore OLS della pendenza.

10 – Simulazione della distribuzione campionaria di \(b\)
Simula 100 000 campioni (n = 30) dal modello
\(Y_i = 1.5\,X_i + \varepsilon_i,\; X_i \sim \mathcal N(0,1),\; \varepsilon_i \sim \mathcal N(0,0.5^2)\).
Per ogni campione calcola \(b̂\). Rappresenta l’istogramma dei 100 000 \(b̂\), riporta media e deviazione standard empiriche e confrontale con l’errore standard teorico.

Soluzioni

1 – Definizione e scopi
1. Previsione – modellare / predire nuove osservazioni.
2. Esplorazione delle associazioni – quantificare relazioni \(X \rightarrow Y\).
3. Estrapolazione – generalizzare dal campione alla popolazione.
4. Inferenza causale – stimare l’effetto di un intervento quando il design lo consente.

2 – Errore di specificazione
Omettere un predittore rilevante ⇒ i residui assorbono la sua varianza ⇒ stime di \(b\) distorte (bias) e varianze sottostimate → inferenze non valide.

3 – Stima del modello

library(rio); library(here)
kidiq <- import(here("data","kidiq.dta"))

mod <- lm(kid_score ~ mom_iq, data = kidiq)
coef(mod)
#> (Intercept)      mom_iq 
#>       25.80        0.61

Esempio di output

(Intercept)     mom_iq 
 25.800        0.610 

4 – Interpretazione della pendenza
Un punto in più di QI materno è associato, in media, a 0.61 punti di QI del figlio.

5 – Indice R²

summary(mod)$r.squared
#> [1] 0.201

Output ≈ 0.20 ⇒ il 20 % della varianza di kid_score è spiegato da mom_iq.

6 – Centratura

kidiq$mom_iq_c <- kidiq$mom_iq - mean(kidiq$mom_iq)
mod_c <- lm(kid_score ~ mom_iq_c, data = kidiq)
coef(mod_c)
#> (Intercept)    mom_iq_c 
#>       86.80        0.61

L’intercetta ora ≈ media del QI dei bambini quando il QI materno è medio.
La pendenza resta 0.61.

7 – Calcolo manuale di b

b_manual <- cov(kidiq$mom_iq, kidiq$kid_score) / var(kidiq$mom_iq)
all.equal(b_manual, coef(mod)["mom_iq"])
#> [1] "names for current but not for target"

TRUE → concordanza perfetta (salvo arrotondamenti).

8 – σ̂ vs σ_CV

# (a) σ̂
sigma_hat <- summary(mod)$sigma

# (b) LOOCV
n <- nrow(kidiq)
res_cv2 <- numeric(n)
for (i in seq_len(n)){
  fit_i <- lm(kid_score ~ mom_iq, data = kidiq[-i,])
  res_cv2[i] <- (kidiq$kid_score[i] - predict(fit_i, kidiq[i,]))^2
}
sigma_CV <- sqrt(mean(res_cv2))

c(sigma_hat = sigma_hat, sigma_CV = sigma_CV)
#> sigma_hat  sigma_CV 
#>     18.27     18.31

σ_CV tende a superare σ̂ perché ogni predizione è fatta su dati che non hanno “visto” quell’osservazione → niente sovradimensionamento.

9 – Assunzioni Gauss‑Markov
1. Linearità nei parametri
2. Campionamento casuale IID
3. Esogeneità \(E(\varepsilon_i\!\mid X_i)=0\) ← questa garantisce non‑distorsione
4. Omoschedasticità
5. Assenza di collinearità perfetta

10 – Simulazione

set.seed(123)
beta  <- 1.5; sigma_e <- 0.5; n  <- 30; nrep <- 1e5
x <- rnorm(n)
b_hat <- replicate(nrep, {
  y <- beta * x + rnorm(n, sd = sigma_e)
  cov(x,y) / var(x)
})

mean_emp <- mean(b_hat)
sd_emp   <- sd(b_hat)
se_theo  <- sqrt(sigma_e^2 / sum((x - mean(x))^2))
c(media_empirica = mean_emp, sd_empirica = sd_emp, SE_teorico = se_theo)
#> media_empirica    sd_empirica     SE_teorico 
#>        1.50006        0.09477        0.09464

L’istogramma (codice sotto) mostra la forma quasi normale centrata su 1.5.

hist(b_hat, breaks = 60, freq = FALSE, main = "Distribuzione campionaria di b̂",
     xlab = "b̂"); curve(dnorm(x, mean_emp, sd_emp), add = TRUE, lwd = 2)

- La media empirica ≈ 1.5 ⇒ unbiased.
- La sd empirica ≈ SE teorico ⇒ formula varianza confermata.
- Distribuzione simmetrica ≈ Normale ⇒ normalità asintotica verificata.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.4.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/UTF-8/C/C/C/C
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] broom_1.0.8      thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0
#>  [5] see_0.11.0       gridExtra_2.3    patchwork_1.3.0  bayesplot_1.12.0
#>  [9] psych_2.5.3      scales_1.4.0     markdown_2.0     knitr_1.50      
#> [13] lubridate_1.9.4  forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4     
#> [17] purrr_1.0.4      readr_2.1.5      tidyr_1.3.1      tibble_3.2.1    
#> [21] ggplot2_3.5.2    tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.3     stringi_1.8.7      lattice_0.22-7    
#>  [4] hms_1.1.3          digest_0.6.37      magrittr_2.0.3    
#>  [7] evaluate_1.0.3     grid_4.5.0         timechange_0.3.0  
#> [10] RColorBrewer_1.1-3 fastmap_1.2.0      Matrix_1.7-3      
#> [13] R.oo_1.27.0        rprojroot_2.0.4    jsonlite_2.0.0    
#> [16] R.utils_2.13.0     backports_1.5.0    mgcv_1.9-3        
#> [19] mnormt_2.1.1       cli_3.6.5          rlang_1.1.6       
#> [22] R.methodsS3_1.8.2  splines_4.5.0      withr_3.0.2       
#> [25] tools_4.5.0        parallel_4.5.0     tzdb_0.5.0        
#> [28] pacman_0.5.1       vctrs_0.6.5        R6_2.6.1          
#> [31] lifecycle_1.0.4    htmlwidgets_1.6.4  pkgconfig_2.0.3   
#> [34] pillar_1.10.2      gtable_0.3.6       glue_1.8.0        
#> [37] haven_2.5.4        xfun_0.52          tidyselect_1.2.1  
#> [40] rstudioapi_0.17.1  farver_2.1.2       htmltools_0.5.8.1 
#> [43] nlme_3.1-168       labeling_0.4.3     rmarkdown_2.29    
#> [46] compiler_4.5.0

Bibliografia

Alter, U., Too, M. A., & Cribbie, R. A. (2025). Navigating the Bayes maze: The psychologist’s guide to Bayesian statistics, a hands-on tutorial with R code. International Journal of Psychology, 60(1), e13271.

Caudek, C., & Luccio, R. (2001). Statistica per psicologi (III rist. 2023, Vol. 11, p. 320). Laterza.

Fox, J. (2015). Applied regression analysis and generalized linear models. Sage publications.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.

Schervish, M. J., & DeGroot, M. H. (2014). Probability and statistics (Vol. 563). Pearson Education London, UK:

Stigler, S. (1986). The History of Statistics. Belknap Harvard.

---

1. Per un approfondimento sull’approccio frequentista alla regressione, si veda, per esempio, Caudek & Luccio (2001).