Correlazione e regressione lineare
Introduzione
Confronto fra variabili quantitative
Dopo aver visto il confronto fra due variabili categoriali ed il confronto fra una variabile categoriale ed una variabile a intervalli (limitatamente al caso del confronto fra due soli gruppi), analizziamo ora il confronto fra due variabili quantitative.
Anche in questo caso, la statistica inferenziale si propone di valutare se esiste una relazione fra le variabili, e se questa relazione è significativa.
Andamento congiunto
Nei casi visti in precedenza, nei capitoli ch:chiquadro e ch:t_test, ci si chiedeva se l'appartenenza di una osservazione ad una categoria od un'altra della variabile A influiva sulla distribuzione della variabile B.
Sia nel caso del test del $\chi^2$ che del test t di Student, che (come vedremo) nell'analisi della varianza, ci si concentra sulla significatività delle differenze.
Nel caso della relazione fra due variabili numeriche, invece, ci si chiede se le due variabili si muovono assieme: se al crescere di una cresce anche l'altra (correlazione positiva), o se al crescere di una l'altra cala (correlazione negativa), e ci si chiede se questo andamento congiunto sia significativo o sia dovuto al caso.
Ipotesi nulla e ipotesi estrema
Anche in questo caso, è opportuno partire dall'ipotesi nulla, ovvero dall'ipotesi che non vi sia alcun legame fra le due variabili.
Anche in questo caso, è necessario identificare una statistica, ovvero una misura dell'aspetto rilevante.
Da un punto di vista didattico, però, può essere utile focalizzarci sulla situazione estrema di una totale correlazione fra le due variabili.
Nell'esempio che segue, ci limiteremo ad analizzare la circostanza di una correlazione positiva, ma il principio può essere generalizzato alle correlazioni negative.
Esempio: misurare le precipitazioni
Immaginiamo che un osservatorio metereologico debba misurare le precipitazioni atmosferiche (pioggia) nel corso dell'anno. Per farlo, viene raccolto in un bacino di un metro quadrato l'acqua piovana, e ad ogni pioggia l'acqua raccolta viene misurata.
Immaginiamo che il responsabile della misurazione sia molto pignolo, e che decida di misurare sia il volume dell'acqua in litri, che il suo peso in kilogrammi. Naturalmente, ci si aspetta che vi sia una corrispondenza perfetta (salvo errori di misurazione) fra il peso e il volume. L'esempio è volutamente banale, ma ci serve per immaginare la situazione estrema della totale correlazione.
Rapporto fra litri e kg
# generiamo delle osservazioni con distribuzione normale
x <- abs(rnorm (50,8,5))+0.5
# rappresentiamo la correlazione "perfetta" fra litri e kg
plot(x,x,xlab="litri",ylab="kg", main="Correlazione perfetta")
La linea retta
L'aspetto più saliente del grafico fig:correlazione--litri_kg_1 è che le misure relative a peso e a volume si dispongono lungo una linea retta.
Il secondo aspetto è che, grazie alla linea, conoscendo il valore in litri, possiamo dedurre il peso in kg, e viceversa.
Diverse distribuzioni di esempio
sequenza <- seq(0,100, by=20)
k <- length(sequenza)
ycoeffs <- vector(mode = "numeric", length = k)
xcoeffs <- vector(mode = "numeric", length = k)
estimates <- vector(mode = "numeric", length = k)
x1 <- vector(mode = "numeric", length = 100)
y1 <- vector(mode = "numeric", length = 100)
x1 <- runif (100,0,1)
conta <- 1
par(mfrow = c(2,3))
for (b in sequenza) {
a <- 100-b
uni2 <- runif (100,0,1)
y1 <- x1*a+uni2*b
# x1 <- (x1-mean(x1))/sd(x1)
# y1 <- (y1-mean(y1))/sd(y1)
linmod <- lm (y1~x1)
ycoeff <- linmod$coefficients[2]
xcoeff <- lm (x1~y1)$coefficients[2]
cortest <- cor.test (x1,y1)
estimate <- cortest$estimate
pvalue <- cortest$p.value
ycoeffs[conta] <- ycoeff
xcoeffs[conta] <- xcoeff
estimates [conta] <- estimate
conta <- conta + 1
testo <-paste("corr.=",round(estimate,3),"p=",round(pvalue,3)) # "a=", a, "b=", b,
testo
plot (x1,y1, main=testo)
abline(linmod)
abline(h=mean(y1))
abline(v=mean(x1))
}
Più realistica la simulazione della figura fig:correlazione--distribuzioni_esempi, in cui sono rappresentate 6 figure. I valori sull'asse x sono gli stessi per ogni figura (x1 <- runif (100,0,1), runif in quanto la distribuzione è uniforme). I valori sull'asse y sono pari alla somma dei valori di x e dei valori generati casualmente, di nuovo con distribuzione per ogni grafico (uni2 <- runif (100,0,1)). La differenza è nel peso della somma. Nel primo grafico uni2 pesa 0, e dunque y=x. Nell'ultimo caso x1 0, e dunque le due variabili sono indipendenti. Nei grafici intermedi il peso di uni2 cresce al calare del peso di x1.
La retta di regressione
In questi grafici abbiamo implicitamente introdotto la retta di regressione, che verrà discussa nelle prossime pagine. Per ora, ci basti sapere che la retta di regressione ci permette di fare una stima del valore di y conoscendo il valore di x.
Nel primo grafico, proprio come nell'esempio precedente, conoscere il valore di x ci permette di inferire perfettamente il valore di y, in quanto le osservazioni di y cadono perfettamente sulla retta.
Già nel secondo grafico, però, questa previsione non è più perfetta: conoscendo x, possiamo soltanto stimare il valore di y.
Mano a mano che il legame diventa meno importante, la pendenza della linea di regressione diminuisce, fino a quasi sovrapporsi alla linea della media di y nell'ultimo grafico, dove le due variabili sono indipendenti.
Ciò che questo andamento ci dice è che conoscere il valore di x contribuisce sempre meno alla nostra conoscenza di y.
Dai grafici della seconda figura possiamo notare che la linea di regressione si incrocia, in tutti i casi, nella linea della media di x e nella linea della media di y nello stesso punto.
Questo significa che il valore stimato di y quando x è pari a $\bar$ è $\bar$.
Analisi grafica
Come abbiamo visto nelle analisi univariate e nelle altre statistiche bivariate, la visualizzazione dei dati e l'analisi visiva, qualitativa delle distribuzioni è parte integrante del processo descrittivo e inferenziale.
Il grafico utilizzato nella visualizzazione di due variabili quantitative è il grafico di dispersione, scatterplot. In R, si ottiene usando la funzione plot (x,y).
Grafico di dispersione: cosa guardare
Quali sono gli aspetti più salienti che bisogna osservare in un grafico a dispersione?
La forma generale della distribuzione.
La direzione dell'associazione: è positiva o negativa?
La forma della la relazione? È una linea retta oppure no? È una curva?
La forza dell'associazione: i punti osservati sono vicini o lontani dalla linea di regressione?
La presenza di outliers: vi sono osservazioni molto lontane dalla retta di regressione? È possibile che queste osservazioni insolite siano dovute ad errori?
È ipotizzabile che l'andamento del grafico lasci intendere l'influenza di una variabile terza?
Analisi inferenziale
Coefficiente di correlazione
La correlazione è una relazione lineare fra due variabili a intervallo o a rapporti.
È importante sottolineare che la correlazione non distingue fra variabile indipendente e dipendente, e tratta le due variabili simmetricamente: la correlazione fra Y e X equivale alla correlazione fra X e Y.
L'analisi inferenziale è finalizzata a calcolare se esiste una relazione fra le due variabili, e se la relazione è statisticamente significativa.
Correlazione: cautele
Prima di applicare la statistica, è necessario tener conto di alcuni possibili problemi:
La correlazione misura relazioni di tipo lineare. Se la relazione non è di tipo lineare, la correlazione non è appropriata.
Soprattutto se l'insieme di osservazioni ha una bassa numerosità, è possibile che degli outliers condizionino fortemente il risultato.
Se una terza variabile, anche categoriale, ha una influenza significativa su una o entrambe le variabili misurate, è possibile che, non tenendone conto, si calcolino delle correlazioni non appropriate.
Correlazione e causazione
La correlazione non implica causazione. La correlazione, infatti, può essere attribuibile a
Causazione diretta: A causa B.
Causa comune: C causa sia A che B.
Fattore confounding: l'andamento della variabile dipendente può essere condizionato da un fattore esterno che non ha nulla a che fare con la variabile indipendente.
Semplice coincidenza.
Requisiti per inferire causazione
Affinché si possa inferire causazione è necessario che:
L'associazione sia abbastanza forte.
Vi sia la possibilità di manipolare la variabile indipendente, e che il valore della variabile dipendente cambi di conseguenza.
Vi sia un chiaro rapporto temporale: la causa deve precedere l'effetto.
Le misure devono essere consistenti, e dunque replicabili.
I risultati siano teoreticamente plausibili e coerenti con altre evidenze empiriche.
L'associazione sia specifica (e dunque non possa essere attribuita a cause comuni o altri confounding).
Modelli Lineari Generalizzati
Sia l'analisi della varianza che la regressione lineare sono casi particolari della metodologia nota come Modelli Lineari Generalizzati.
La differenza più importante fra la regressione e l'ANOVA è che in un caso la variabile indipendente è a intervalli, nell'altro categoriale. Nell'ANOVA, dunque, non si fanno assuzioni sulla linearità della relazione.
Approccio intuitivo
Per misurare la correlazione lineare, abbiamo bisogno di una statistica che abbia alcune caratteristiche:
sia pari a 0 in assenza di correlazione
sia positiva quando la correlazione è positiva, e negativa quando è negativa
che abbia un valore assoluto massimo, che identifica la circostanza in cui la correlazione è perfetta
che sia standardizzata, ovvero che non dipenda dai valori assoluti delle variabili.
In termini formali: $$-1 \leq r \leq 1$$
Linea di regressione e medie
Come abbiamo osservato, la linea di regressione si incrocia sempre con le due medie. Usando le linee che identificano le medie di x e di y possiamo dividere il grafico di dispersione in 4 quadranti. Se le osservazioni si distribuiscono principalmente nel quadrante in alto a destra e in quello in basso a sinistra, la correlazione è positiva. Se sono più frequenti nei quadranti in alto a sinistra o in basso a destra, la correlazione è negativa.
Per ottenere una misura standardizzata del rapporto fra le due variabili, possiamo decidere di trasformare sia la variabile x che la y in punteggi zeta.
Trasformazione in punti z
Con la trasformazione, otteniamo due variabili con media pari a zero e deviazione standard pari ad uno. La retta di regressione, a questo punto, incrocia le due medie nella posizione 0,0 del grafico.
Il quadrante in basso a sinistra raccoglie le osservazioni in cui sia x che y sono negative.
Il quadrante in alto a destra raccoglie le osservazioni in cui sia x che y sono positive.
Il quadrante in basso a destra raccoglie le osservazioni in cui x è positiva e y è negativa.
Il quadrante in alto a sinistra raccoglie le osservazioni in cui y è positiva e x è negativa.
A questo punto appare evidente che, moltiplicando x per y, otterrò valori positivi nei due quadranti alto sinistra e basso destra, e valori negativi nei quadranti basso destra, alto sinistra.
Cosa succede se sommo la moltiplicazione x * y di tutte le osservazioni, e divido per il numero di osservazioni (per la precisione, per n-1)? Otterrò un valore che sarà positivo in caso di correlazione positiva, negativo in caso di correlazione negativa, e sarà prossimo allo zero in caso di assenza di correlazione.
Inoltre, vedremo che il valore più alto che questa misura può raggiungere è pari ad 1, e di conseguenza il valore più basso è pari a -1.
Correlazione lineare: la formula
$$\label r= \frac{1}\sum_^n [(\frac{x_i-\bar}{\sigma_x}) (\frac{y_i-\bar}{\sigma_y})]$$
Ricordando che $\bar$ e $\bar$ sono le medie di X e Y e che $\sigma_X$ e $\sigma_Y$ sono le deviazioni standard di X e Y, e ricordando che il calcolo del punteggio z è pari a $z= \frac{x_i-\bar}{\sigma_x}$
possiamo riscrivere r come
$r= \frac{1}\sum_^n (z(x_i) z(x_i))$
Correlazione lineare: assunti
Poiché il calcolo della correlazione utilizza la trasformazione dei punteggi grezzi in punti zeta, si assume che entrambe le variabili siano almeno a livello di scala di intervallo e abbiano una distribuzione normale.
Le osservazioni su entrambe le variabili devono essere stocasticamente indipendenti (ovvero, il valore di di una osservazione non deve condizionare il valore di un'altra osservazione)
Infine, il rapporto fra le due variabili sia di tipo lineare. Vedremo nelle prossime sezioni quali alternative esistono in caso di non linearità del rapporto.
Distribuzione dell'errore
Anche in questo caso utilizziamo la consueta sequenza logica:
- identificazione di una misura della relazione
- identificazione di una popolazione virtuale
- estrazione casuale di k campioni
- calcolo del vettore delle k misure
- osservazione della distribuzione delle misure
- calcolo del p-value attraverso il confronto con il vettore delle misure
- identificazione di una distribuzione teorica che, previo opportuna trasformaizone, mappa quella osservata
- calcolo del p-value utilizzando la probabilità della distribuzione teorica identificata
- utilizzo della funzione di R
La simulazione
- La misura della relazione è il coefficiente r identificato sopra.
- Attraverso la funzione
rnorm ()possiamo estrarre n campioni di osservazioni casuali da una popolazione con specifica media e deviazione standard. Poiché siamo interessati a misurare la relazione fra due variabili, estraiamo due misure per ogni campione - Utilizziamo la funzione
scale()per trasformare le osservazioni in punti z. - Calcoliamo la statistica r, e la salviamo nel vettore delle misure.
R: la simulazione
m <- 100 # numerosità dei campioni
k <- 10000 # numero di campioni generati
relazioni <- vector(mode = "numeric", length = k)
for (i in 1:length(relazioni)) {
x1 <- rnorm (m,20,2)
x2 <- rnorm (m,50,6)
x1 <- scale(x1)
x2 <- scale (x2)
erre <- sum(x1*x2)/(length(x1)-1)
relazioni[i] <- erre
}
Grafico della distribuzione dell'errore
Visualizziamo la distribuzione dell'errore di r nella figura
hist (relazioni)
Rapporto fra distribuzione dell'errore e t
Poiché il valore possibile di $r$ varia nel range $-1 \leq r \leq 1$, anche la distribuzione dell'errore varia nello stesso range, e si concentra attorno ai valori -0.1, 0.1.
La forma della distribuzione approssima la t di Student, previa opportuna trasformazione. Per arrivare alla distribuzione t va applicata la trasformazione $$t = \frac{\sqrt{\frac{1-r^2}}} \label$$
P-value calcolato sulla distribuzione t
Per calcolare il p-value basandoci sulla distribuzione t, dovremmo dunque trasformare r in t, e poi calcolare la probabilità di t.
Alcuni esempi
Facciamo ora alcuni esempi, generando differenti casi di variabili bivariate.
Primo esempio: variabili indipendenti
In questo caso, generiamo due variabili casuali indipendenti. Ci aspettiamo un r prossimo allo 0.
x1 <- rnorm (100,20,2)
y1 <- rnorm (100,50,3)
sx1 <- scale(x1)
sy1 <- scale (y1)
erre1 <- sum(sx1*sy1)/(length(sx1)-1)
erre1
## [1] 0.006969092
La statistica r, dunque, è pari a $0.0069691$.
colori <- c(2,2,1)
colore <- colori[sign(sx1*sy1)+2]
plot (sx1,sy1,col=colore)
abline (a=0,b=erre1)
abline (v=0);
abline(h=0)
Calcolo p-value su simulazione
Calcoliamo il p-value usando la distribuzione osservata. Poiché la nostra ipotesi è a due code, moltiplichiamo p per 2. Poiché le due variabili sono indipendenti, la nostra previsione è che p sia superiore a 0.05.
p_value_simulazione1 <- rank (c(-abs(erre1),relazioni))[1]/(length(relazioni)+1)
p_value_simulazione1 * 2
## [1] 0.9509049
Il p-value è dunque pari a 0.9509049.
Secondo esempio: variabili correlate
In questo secondo esempio, la variabile y è creata in modo da correlare con x. Ci aspettiamo un r relativamente alto, e un p basso.
y2 <- x1 + rnorm (100,0,2)
sy2 <- scale (y2)
erre2 <- sum(sx1*sy2)/(length(sx1)-1)
erre2
## [1] 0.7841576
Disegnamo il grafico.
colori <- c(2,2,1)
colore <- colori[sign(sx1*sy2)+2]
plot (sx1,sy2,col=colore)
abline (a=0,b=erre2)
abline (v=0); abline(h=0)
p_value_simulazione2 <- rank (c(-abs(erre2),relazioni))[1]/(length(relazioni)+1)
p_value_simulazione2 * 2 # due code
## [1] 0.00019998
Come previsto, il p-value -- calcolato sulla distribuzione osservata dalla simulazione -- è basso: $1.9998 × 10-4$.
Uso della distribuzione teorica
Dopo aver visto il calcolo del p-value usando la distribuzione della simulazione, utilizziamo la probabilità calcolata a partire dalla distribuzione teorica, t.
In primo luogo confrontiamo la distribuzione generata dalla simulazione con la distribuzione teorica, sovrapponendo le due curve.
relazioni_t <- relazioni / (sqrt((1-relazioni^2)/(100-2) ))
x <- seq.int(-4, 4, by=0.05)
y <- dt(x, 100-2)
plot (density (relazioni_t),col=4, main="Sovrapposizione delle distribuzioni")
lines (x,y, type="l",col=3)
Calcoliamo il p-value usando la distribuzione t
Dal grafico appare che le due distribuzioni sono estremamente simili. Possiamo dunque calcolare il punteggio t, usando la funzione eq:rstimat. Poi, calcoliamo il p-value usando la funzione pt
Primo esempio
Mostriamo il calcolo, ed il risultato, del primo esempio.
t1 <- abs(erre1) / (sqrt((1-erre1^2)/(100-2) ))
t1
## [1] 0.06899216
p_value_t1 <- (1 - pt(t1, df = (100-2)))
p_value_t1 * 2
## [1] 0.9451364
Dunque, t è pari a $0.0689922$, $p= 0.9451364$.
Secondo esempio
Ripetiamo il calcolo, con il secondo esempio.
t2 <- abs(erre2)/(sqrt((1 - erre2^2)/(100 - 2)))
p_value_t2 <- (1 - pt(t2, df = (100 - 2)))
p_value_t2 * 2
## [1] 0
t è pari a $12.5092524$, $p= 0$.
R: uso di cor.test
E come di consueto, terminiamo mostrando l'uso della funzione cor.test (x,y). Iniziamo con il primo esempio.
(correlazione1 <- cor.test (x1,y1))
##
## Pearson's product-moment correlation
##
## data: x1 and y1
## t = 0.068992, df = 98, p-value = 0.9451
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1897087 0.2031092
## sample estimates:
## cor
## 0.006969092
Leggere i risultati
La funzione ci dice che ha applicato la statistica Pearson's product-moment correlation. Che i gradi di libertà sono 98 (n-2). Il valore della statistica è 0.0689922 ed il p-value è 0.945. La correlazione è pari a $0.0069691$. Come prevedibile il p-value è superiore a 0.05, e dunque non si può rifiutare l'ipotesi nulla di indipendenza fra le due variabili.
Secondo esempio
(correlazione2 <- cor.test (x1,y2))
##
## Pearson's product-moment correlation
##
## data: x1 and y2
## t = 12.509, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6947462 0.8497022
## sample estimates:
## cor
## 0.7841576
In questo caso, l'ipotesi nulla va rifiutata, in quanto $p-value < <2e-16$.
Regressione lineare
Precipitazioni: dal volume al peso
Torniamo all'esempio, banale, delle precipitazioni. Grazie alla correlazione perfetta fra volume e peso, data una osservazione, conoscendo il volume, possiamo calcolare il peso.
Se guardiamo al grafico, possiamo notare che il peso stimato incrocia il volume proprio lungo la linea di regressione.
La linea di regressione, dunque, stima il valore di y a partire da x (e viceversa, il valore di x conoscendo y).
Esempio 2A: noleggio automobili
Immaginiamo che una agenzia di noleggio auto applichi, per un modello, la seguente tariffa:
- importo fisso di 40 euro al giorno
- 0.08 euro a km percorso
Conoscendo questi due valori, possiamo prevedere esattamente quanto spenderemo. Ad esempio, se ho percorso 140km, spenderò 40 + 0.08*140 = 51.20 euro.
Esempio 2B: noleggio automobili
Immaginiamo che un'altra agenzia applichi, invece, il fisso di 40 euro più il costo della benzina consumata. Immaginiamo che l'auto in questione consumi, in media, 0.06 euro a km.
In questo caso, quando spenderemo, dopo aver fatto 140 km? Il calcolo è lo stesso: 40 + 0.06*140 = 48.40. In questo caso, però, questo valore è solo una stima di quello che ci aspettiamo di spendere, in quanto non possiamo essere sicuri che la benzina consumata sia esattamente pari a 8.40 euro.
I valori di 0.06 euro a km, infatti, sono una statistica, calcolata in seguito ad una serie di osservazioni, dove si sono misurati i km effettuati e la benzina consumata.
Il numero di km fatti è il miglior predittore della benzina consumata, ma non è l'unico. Il tipo di percorso e il tipo di guida, fra gli altri, influenzano il consumo.
Quando pagheremo per il noleggio della seconda agenzia, noi possiamo aspettarci un conto di circa 48 euro e 40, ma sappiamo che in quella stima ci sarà un errore, legato a quei fattori che incidono sul consumo ma che non sono annoverati nel calcolo.
Regressione lineare: il modello
Generalizzando dall'esempio precedente, nel modello di regressione lineare bivariato (X,Y), la variabile Y può essere rappresentata tramite la relazione lineare $$Y = \beta_0 + \beta_1 X + \epsilon$$ ed il valore $$y_i = \beta_0 + \beta_1 x_i + \epsilon_i (i = 1 ... n)$$ Infine, $$\hat = \beta_0 + \beta_1 x_i (i = 1 ... n)$$ dove $\hat$ è il valore stimato di y.
Spesso si usa la forma $\alpha + \beta x_i$.
Le componenti
$\beta_0$ rappresenta l'intercetta, ovvero il valore di y quando x = 0.
$\beta_1$ (o, nella regressione bivariata, semplicemente $\beta$) rappresenta la pendenza della linea, ed è pari alla differenza fra $\hat = f(x)$ e $\hat = f(x+1)$. Nell'esempio, rappresenta il costo supplementare per ogni km percorso.
$\epsilon$ rappresenta la variabile di errore, ovvero quella varianza in y che non può essere spiegata da x.
Per massimizzare la predittività della regressione è dunque necessario scegliere due parametri $\beta_0$ e $\beta_1$ capaci di minimizzare l'errore e, possibilmente, di escludere dei bias (errori sistematici).
Somma dei quadrati degli errori
Nella regressione lineare semplice, la misura dell'errore $\epsilon$ si basa sulla somma dei quadrati degli errori (in inglese sum of squared residuals, SSR): $$SSR = \sum_^n e^2_i = \sum_^n (\hat - y_i)^2 = \sum_^n ((\beta_0 + \beta_1 x_i) - y_i)^2$$ dove $\hat$ è il valore stimato di y.
Si tratta dunque di identificare i parametri $\beta_0$ e $\beta_1$ capaci di minimizzare SSR.
Stime di $\beta_0$ e $\beta_1$
Le due stime che minimizzano la somma dei quadrati degli errori (SSR) sono $$\label \hat{\beta_1} = \frac{\sum_^n (x_i-\bar)(y_i-\bar)}{\sum_^n(x_i-\bar)^2}$$ $$\label \hat{\beta_0} = \bar - \hat{\beta_1} \bar$$
Proporzione di varianza spiegata e residua
[sec:reglinspiegata] La varianza di Y può essere divisa fra la varianza spiegata da $\beta_0 + \beta_1 X$ e la varianza residua, o di errore. La proporzione di varianza spiegata è pari a $R^2=r^2$. $$R^2 = 1- \frac{SSR / (n-1)}{var(Y)} = 1- \frac{\sum_^n e_i^2}{\sum_^n (y_i-\bar)_i^2}$$ $0 \leq R^2 \leq 1$. $R^2$ è pari al quadrato di r.
Va notato che quanto SSR è 0 $R^2$ è 1, e quando $SSR = var(Y)$ $R^2 = 0$.
Assunti della regressione lineare
- La relazione fra le variabili dev'essere lineare
- L'errore è una variabile casuale con media zero e distribuzione normale
- Gli errori non sono fra loro correlati
- La varianza dell'errore è costante
Distribuzione di $y_i$
Tenuto conto che $$y_i = \beta_0 + \beta_1 x_i + \epsilon_i (i = 1 ... n)$$ se $\epsilon_i$ ha distribuzione normale, media 0, varianza $\sigma^2$, la distribuzione di $y_i$ sarà $$y_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2) (\forall i = 1 ... n)$$ dove $N()$ è la distribuzione normale.
R: la funzione lm ()
In R, la regressione lineare si calcola utilizzando la funzione lm(). La sintassi usata è $lm (y \sim x)$ dove y è la variabile dipendente e x la variabile indipendente.
Per disegnare la retta di regressione, si passa il risultato di lm() alla funzione abline, che disegna una linea con parametri a e b.
modello1 <- lm (y1~x1)
summary (modello1)
##
## Call:
## lm(formula = y1 ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.987 -2.071 -0.001 2.186 6.431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.05800 2.97510 16.490 <2e-16 ***
## x1 0.01031 0.14950 0.069 0.945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.248 on 98 degrees of freedom
## Multiple R-squared: 4.857e-05, Adjusted R-squared: -0.01016
## F-statistic: 0.00476 on 1 and 98 DF, p-value: 0.9451
summary(lm()) ci restituisce molte informazioni. Intercept, ad esempio, ci dice se l'intercetta è significativamente diversa da 0 (informazione generalmente poco interessante).
Molto più importante la seconda linea, che calcola t e il p-value di $\beta_1$. R-squared è $R^2$. Il p-value è calcolato anche attraverso la statistica F (che non abbiamo affrontato). Il risultato, nel caso di correlazione bivariata, è lo stesso: $p-value: 0.9451364$.
Il grafico e la retta di regressione
Abbiamo già introdotto la retta di regressione, abline, che utilizza proprio il risultato del modello lineare lm().
plot (x1,y1)
abline (modello1)
Varianza dei residui, $R^2$
Nel paragrafo sec:reglinspiegata abbiamo introdotto il concetto di rapporto fra varianza spiegata e varianza residua. la funzione lm() restituisce, fra le altre cose, i residui: $residuals. Sappiamo che la varianza totale è pari alla varianza spiegata più la varianza residua. Sappiamo dunque la varianza spiegata è pari a $R^2 = 1- \frac{var(residui)}{var(Y)}$. Nelle prossime righe di R calcoliamo $R^2$, e lo confrontiamo con il $r^2$, per mostrare che sono uguali.
residui1 <- modello1$residuals
R2_1 <- 1 - var(residui1) / var (y1)
R2_1
## [1] 4.856824e-05
erre1^2
## [1] 4.856824e-05
Come vediamo, $R^2$ è pari a $4.8568238 × 10-5$, e $r^2$ è pari a $4.8568238 × 10-5$.
Grafico dei residui
Il grafico dei residui ci permette di visualizzare la distribuzione dell'errore. È importante per verificare gli assunti del modello.
plot (modello1$fitted.values, modello1$residuals)
abline (lm(modello1$residuals~modello1$fitted.values))
lines(smooth.spline(modello1$fitted.values, modello1$residuals),
col = "red", lwd = 2)
Il secondo esempio
Usiamo lm() sulle variabili del secondo esempio.
modello2 <- lm (y2~x1)
summary (modello2)
##
## Call:
## lm(formula = y2 ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9260 -1.0558 -0.0674 0.9119 5.2574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.26498 1.62812 -0.163 0.871
## x1 1.02343 0.08181 12.509 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.777 on 98 degrees of freedom
## Multiple R-squared: 0.6149, Adjusted R-squared: 0.611
## F-statistic: 156.5 on 1 and 98 DF, p-value: < 2.2e-16
In questo caso, il modello è significativo, in quanto il p-value di x1 è $5.0242137 × 10-22$.
Visualizziamo il grafico.
plot (x1,y2)
abline (modello2)
Usiamo i residui, calcoliamo $R^2$.
residui2 <- modello2$residuals
R2_2 <- 1 - var(residui2) / var (y2)
R2_2
## [1] 0.6149031
erre2^2
## [1] 0.6149031
Residui, $R^2$
Visualizziamo il grafico dei residui su x.
plot (modello2$fitted.values, modello2$residuals)
abline (lm(modello2$residuals~modello2$fitted.values))
lines(smooth.spline(modello2$fitted.values, modello2$residuals), col = "red", lwd = 2)
In questo caso, alla linea di regressione abbiamo aggiunto anche una smooth.spline, ovvero una curva che segue l'andamento della relazione fra punteggi stimati e residui. Questa curva ci può aiutare a capire se l'assunto di linearità è violato.
Violazione degli assunti
Gli assunti della regressione lineare
Ricordiamo gli assunti della regressione lineare:
- La relazione fra le variabili dev'essere lineare
- L'errore è una variabile casuale con media zero e distribuzione normale
- Gli errori non sono fra loro correlati
- La varianza dell'errore è costante (omoskedasticità)
Quali sono le conseguenze della violazione degli assunti?
La violazione di linearità è il caso più importante, di cui discuteremo estesamente nelle prossime slide.
Correlazione fra gli errori
La correlazione fra gli errori è sintomo di non indipendenza fra le misure.
La non indipendenza fra le misure è un problema di cui va tenuto conto nelle misure ripetute.
Nel caso di misure su di un campione estratto casualmente, la non indipendenza delle misure e dell'errore è meno probabile.
La violazione di questo assunto può essere diagnosticata attraverso un test di autocorrelazione dei residui
Omoskedasticità: varianza dell'errore costante
Se la varianza dell'errore non è costante, l'intervallo di confidenza della distribuzione di Y non sarà correttamente predittivo, in quanto sovrastimato nella parte del grafico in cui la varianza dell'errore è minore, e sottostimato nelle parti dove è maggiore.
La violazione dell'omoskedasticità può essere diagnosticata plottando i residui sui valori attesi: se la dispersione degli errori non è omogenea, possiamo sospettare una violazione della costanza della varianza dell'errore.
Normalità dell'errore
La violazione di questo assunto comporta la compromissione della stima sia dei coefficienti $\beta_1$ e $\beta_0$ che dei valori di confidenza della distribuzione di Y su X.
Per verificare graficamente la normalità della distribuzione dell'errore, si può utilizzare il grafico qqnorm e qqline.
Per verificarla inferenzialmente, si possono usare il test di Kolmogorov-Smirnov, ks.test, o il test di normalità Shapiro-Wilk: shapiro.test
Violazione della linearità
È uno degli aspetti più delicati, in quanto può indurre ad inferenze scorrette.
La non linearità può essere diagnosticata attraverso la visualizzazione del grafico di dispersione delle due variabili, il grafico di dispersione dei residui sui valori attesi, o sulla variabile X.
Da un punto di vista inferenziale, è possibile applicare la statistica Harvey-Collier: harvtest (va caricata la libreria lmtest: library(lmtest::harvtest).
In alcune circostanze, è possibile applicare una trasformazione non lineare ad una o entrambe le variabili, per rendere lineare la relazione.
Secondo esempio: testiamo gli assunti
Focalizziamoci sul secondo esempio, relativo a due variabili correlate, e testiamo gli assunti di normalità e di linearità.
Verifico la normalità della distribuzione dell'errore
Uso il test di Shapiro-Wilk per testare la normalità della distribuzione dell'errore.
shapiro2 <- shapiro.test(modello2$residuals)
shapiro2
##
## Shapiro-Wilk normality test
##
## data: modello2$residuals
## W = 0.9771, p-value = 0.07879
Poiché $ p= 0.0787876 $, non rifiuto l'ipotesi nulla di normalità del modello. L'assunto, dunque, non è violato.
Valuto la linearità del modello
Utilizziamo ora il test Harvey-Collier per testare la linearità del modello. Attenzione: per usare la funzione harvtest è necessario, prima, importare la libreria lmtest, con il comando library(lmtest).
library(lmtest)
harvtest2 <- harvtest(y2 ~ x1, order.by = ~ x1)
harvtest2
##
## Harvey-Collier test
##
## data: y2 ~ x1
## HC = 0.3979, df = 97, p-value = 0.6916
Poiché $p= 0.6915768$, non rifiuto l'ipotesi nulla di linearità del modello. L'assunto, dunque, non è violato.
Terzo esempio, non lineare
Infine, un ultimo esempio.
x3 <- runif (100,-2,6)
y3 <- x3^2+ rnorm(100,0,1)
x3 <- x3 + 10
modello3 <- lm (y3~x3)
summary (modello3)
##
## Call:
## lm(formula = y3 ~ x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1393 -4.4576 -0.5511 4.4691 11.9022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -34.9816 2.4777 -14.12 <2e-16 ***
## x3 3.7411 0.2034 18.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.013 on 98 degrees of freedom
## Multiple R-squared: 0.7753, Adjusted R-squared: 0.773
## F-statistic: 338.2 on 1 and 98 DF, p-value: < 2.2e-16
Come vediamo dal codice, la relazione fra y3 e x3, al netto della varianza non spiegata, è data da $x3^2$. Applichiamo il modello lineare, e leggiamo i risultati.
Il grafico
Analizziamo, ora, il grafico.
par(mfrow=c(1,2))
plot (x3,y3,main="grafico di dispersione")
abline (modello3)
plot (modello3$fitted.values,modello3$residuals,main="punteggi attesi vs residui", xlab="attesi", ylab="residui")
abline (lm(modello3$residuals~modello3$fitted.values))
lines(smooth.spline(modello3$fitted.values,modello3$residuals), col='red', lwd=2)
Risulta evidente, dal grafico a sinistra, che la relazione fra x3 e y3 non è lineare. La non linearità è ancor più evidente nel grafico dei residui, a destra.
Usiamo il test Harvey-Collier per valutare la linearità.
harvtest3 <- harvtest(y3 ~ x3, order.by = ~ x3)
harvtest3
##
## Harvey-Collier test
##
## data: y3 ~ x3
## HC = 10.157, df = 97, p-value < 2.2e-16
$p= 6.1566007 × 10-17 < 0.05$, e dunque rifiuto l'ipotesi nulla di linearità del modello. L'assunto di linearità del modello è violato, come previsto dall'analisi qualitativa del grafico.
I 4 data-set di Anscombe
In letteratura, sono noti i 4 insiemi di dati pubblicati da Anscombe. I quattro insieme di dati sono particolari: le quattro y hanno la stessa media (7.5), deviazione standard (4.12), correlazione (0.816) e linea di regressione. Com'è possibile notare dal grafico, però, i quattro insiemi sono qualitativamente molto diversi. L'esempio, è finalizzato a ricordarci che calcolare il modello lineare non basta, e che una attenta analisi dei grafici di dispersione è indispensabile, per evitare di trarre conclusioni inferenziali indebite.
Coefficiente di Spearman
Dipendenza monotona
Nelle circostanze in cui la relazione fra le due variabili non sia lineare, ma tenda ad essere comunque monotona, è possibile utilizzare il modello non-parametrico della correlazione: $\rho$ di Spearman.
In questo modello, il calcolo della relazione si effettua non sui valori delle due variabili, ma sulla loro posizione ordinale.
Questa statistica, pertanto, può essere applicata anche nella circostanza in cui una o entrambe le variabili siano di tipo ordinale.
Assunti del modello di Spearman
Gli assunti sono i seguenti:
- le due variabili devono essere almeno ordinali, non necessariamente ad intervalli;
- su entrambe le variabili, le diverse osservazioni devono essere fra loro indipendenti;
- si assume che vi sia, fra le variabili, una relazione di tipo monotono;
Quarto esempio, sigmoide
Introduciamo un quarto esempio, dove la curva fra x4 e y4 è una sigmoide.
x4 <- rnorm (200,0,4)
y4 <- (1 / (1 + exp(-x4)))*10+rnorm(200,0,0.3)
x4 <- x4+8+rnorm(200,0,0.3)
Il grafico
Prima di ogni calcolo, mostriamo il grafico.
par(mfrow=c(1,1))
plot (x4,y4)
abline(lm(y4~x4))
lines(smooth.spline(x4,y4), col='red', lwd=2)
Dal grafico appare chiaro che una relazione fra x4 e y4 esiste, ma che la relazione non è lineare. Il modello lineare riesce comunque a cogliere la relazione, ma il modello predittivo risulta sostanzialmente scorretto.
Coefficiente di Spearman
Calcoliamo, ora, a mano, il coefficiente $\rho$ di Spearman. Come abbiamo detto, la statistica non parametrica utilizza, al posto dei valori numerici di x e y, le loro posizioni ordinali. Il primo passaggio, dunque, è quello di trasformare i punteggi nei rispettivi ranking (e, nel nostro algoritmo, di scalarli). Poi, utilizziamo la consueta formula per calcolare il coefficiente.
rankx4 <- scale(rank(x4))
ranky4 <- scale(rank(y4))
spearman4 <- sum(rankx4*ranky4)/(length(rankx4)-1)
spearman4
## [1] 0.980393
plot (rankx4,ranky4)
Per capire meglio il meccanismo, disegnamo il grafico di dispersione dei ranking delle due variabili. Come possiamo vedere, la trasformazione rende lineare la relazione monotona.
Coefficiente di Spearman, cor.test
Infine, come di consueto, utilizziamo la funzione di R: cor.test(x4, y4, method = "spearman"). Se si vuole usare la $\rho$ di Spearman, è necessario specificare method = "spearman".
cor_test_S4 <- cor.test (x4,y4,method='spearman')
cor_test_S4
##
## Spearman's rank correlation rho
##
## data: x4 and y4
## S = 26142, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.980393
Leggere i risultati
La funzione calcola la statistica $S= 2.6142 × 104$, $rho= 0.980393$, $p= 0$.