Archivio tag: regressione

L’analisi di covarianza

Una delle cose che spesso mi trovo a spiegare, è come si interpretano i parametri di un modello lineare quando vengono messi insieme predittori di natura quantitativa e qualitativa, ovvero variabili di tipo numerico e di tipo fattore. Se le variabili non sono in interazione fra loro, questo tipo di modello viene chiamato “analisi di covarianza”; se invece le due variabili sono fra loro in interazione, al modello non viene dato alcun nome specifico (si tratta di una terminologia secondo me un po’ fuorviante, che personalmente non condivido).

Mi è capitato diverse volte di veder usare questo modello senza una reale consapevolezza del significato dell’analisi. Spesso un ricercatore si limita a valutare la significatività statistica di ogni singola variabile indipendente, senza pensare che i parametri stimati hanno una vera e propria interpretazione geometrica. La lettura dei parametri di un modello lineare fornisce molte informazioni, che, se sapute leggere e rappresentare adeguatamente, avrebbero tanto da raccontare. Non tutti infatti sanno che, quando nel modello sono presenti due predittori, uno numerico e uno fattore, il modello può ancora essere rappresentato bidimensionalmente, al pari di una regressione semplice.

Se vogliamo capire, è bene partire da un esempio pratico. Per il mio esempio mi sono ispirato a questo articolo, dove viene studiato come la relazione tra la prestazione a un test di matematica e l’ansia provata per il test varia fra maschi e femmine.

Presso questo link è disponibile un dataset che contiene dei dati simulati che riproducono, almeno nella logica, i dati raccolti dagli autori dell’articolo. I numeri sono stati inventati di sana pianta da me, per cui non prendete per oro colato i risultati, che sono un po’ esagerati: quello che ci interessa è il metodo.

Iniziamo scaricando i dati in locale e importandoli in R:

> math <- read.csv2("mathematics.csv")
> str(math)
'data.frame':	140 obs. of  4 variables:
 $ math.perf: int  9 4 2 8 2 4 13 10 6 5 ...
 $ math.anx : int  17 10 12 12 14 6 11 5 10 14 ...
 $ test.anx : int  19 18 24 9 26 17 7 10 15 15 ...
 $ gender   : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...

Il dataset contiene quattro variabili: la prestazione al test in matematica (math.perf), una misura dell’ansia provata verso i test di valutazione delle abilità matematiche (math.anx), una misura dell’ansia provata verso i test di valutazione di abilità in generale (test.anx) e il sesso del partecipante alla ricerca (gender).

Ci interessa studiare la relazione tra la prestazione al test di matematica (math.perf) e l’ansia provata verso le specifiche situazioni in cui le abilità matematiche vengono valutate attraverso un test (math.anx).

Attenzione: nel seguito della trattazione mi appoggerò a dei grafici per illustrare i risultati, ma, per evitare di appesantire troppo il discorso, non mostrerò il codice R utilizzato. Magari la costruzione di questi grafici sarà oggetto di altri post; per ora, vi basti sapere sono stati realizzati sfruttando il pacchetto ggplot2.

Analisi preliminari


La relazione tra math.perf e math.anx sembra esistere veramente. Disponendo math.anx sull’asse orizzontale di uno scatterplot e math.perf sull’asse verticale, notiamo infatti che all’aumentare dell’ansia diminuisce la prestazione, in maniera approssimativamente lineare.

Scatterplot modello additivo

Proviamo ad adattare un modello lineare per verificare se è possibile predire la prestazione a partire dall’ansia:

> mod1 <- lm(math.perf ~ math.anx, data=math)
> summary(mod1)

Call:
lm(formula = math.perf ~ math.anx, data = math)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.1997  -6.2352  -0.8446   5.3344  23.6583 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.57838    2.06327   6.096 1.03e-08 ***
math.anx    -0.04734    0.13835  -0.342    0.733    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.283 on 138 degrees of freedom
Multiple R-squared:  0.0008476,	Adjusted R-squared:  -0.006393 
F-statistic: 0.1171 on 1 and 138 DF,  p-value: 0.7328

La tabella dei coefficienti mostra che il parametro associato a math.anx, stimato pari a -0.047, non è significativamente diverso da zero (p = 0.773). Quindi, l’ipotesi non sembra essere valida.

Abbiamo però un’altra variabile nel dataset, gender. Proviamo a distinguere nel grafico i dati dei maschi da quelli delle femmine, colorandoli in maniera differenziata.

Scatterplot modello additivo

Guarda un po’ che coincidenza: sembra che i dati delle femmine si dispongano tutti “sopra” quelli dei maschi! In pratica, sembrerebbe che la relazione tra math.anx e math.perf esista sia nei maschi che nelle femmine, ma le femmine abbiamo delle prestazioni in matematica superiori, per cui avrebbero bisogno di una retta di regressione diversa da quella dei maschi. È infatti come se con il nostro modello dovessimo costruire non una retta di regressione, ma due: una per i maschi, con un’intercetta più bassa, e una per le femmine, con un’intercetta più elevata.

Il modello additivo


Abbiamo appena visto che, per adattare un buon modello, dovremmo costruire una retta di regressione per i maschi diversa da quella delle femmine. La retta delle femmine dovrebbe infatti essere “più alta”, ovvero dovrebbe avere un’intercetta maggiore.

Quindi, che si fa? Si divide in due il dataset e si fanno due analisi separate? Ovviamente no: si aggiunge nel modello creato precedentemente anche la variabile gender.

> mod2 <- lm(math.perf ~ math.anx + gender, data=math)
> summary(mod2)

Call:
lm(formula = math.perf ~ math.anx + gender, data = math)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1216 -3.2256  0.2374  2.8479 10.6495 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  32.22762    1.45988   22.08   <2e-16 ***
math.anx     -0.86471    0.08249  -10.48   <2e-16 ***
genderM     -16.36525    0.83473  -19.61   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.261 on 137 degrees of freedom
Multiple R-squared:  0.7375,	Adjusted R-squared:  0.7336 
F-statistic: 192.4 on 2 and 137 DF,  p-value: < 2.2e-16
&#91;/code&#93;

Adesso cominciamo a ragionare: non solo i parametri sono tutti significativamente diversi da zero per <i>p</i> < 0.05, ma anche il valore R<sup>2</sup> è molto elevato (0.7375). La prestazione in matematica, quindi, dipende sia dall'ansia verso i test di valutazione delle abilità matematiche che dal sesso.

Il modello che abbiamo appena creato può essere visualizzato, ma prima di farlo è importante capire cosa rappresentino esattamente i parametri stimati.

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">Il significato dei parametri</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

Il modello mod2 comprende tre parametri: Intercept, math.anx e genderM. Analizziamoli uno per uno.

<b>Intercept</b> è l'intercetta della retta di regressione delle femmine. Vi starete chiedendo come io faccia a saperlo... ebbene, questa informazione ci viene fornita dal comando <b>contrasts</b> applicato alla variabile gender: il livello a cui è associato il valore 0 sarà quello di <b>riferimento</b>, mentre il livello a cui è associato il valore 1 sarà quello di <b>confronto</b>. Il parametro denominato <i>Intercept</i> rappresenta l'intercetta del gruppo di riferimento, quindi delle femmine.

[code language="R"]
> contrasts(math$gender)
  M
F 0
M 1

Se non avete dimestichezza con le matrici di contrasto, credo che il modo migliore di affrontare questo passaggio sia di prenderlo come un dogma. Ci torneremo nei prossimi post.

Il parametro genderM è lo scarto fra l’intercetta dei maschi (M) e quella delle femmine. Ovvero, genderM rappresenta la quantità che dobbiamo aggiungere a Intercept per ottenere l’intercetta dei maschi. Si noti che genderM è negativo, per cui all’intercetta delle femmine dobbiamo togliere una quantità: ciò significa che l’intercetta dei maschi, come previsto, sarà inferiore a quella delle femmine.

Dato che l’intercetta delle femmine è pari a 32.23, l’intercetta dei maschi sarà pari a 32.23 + (-16.37), ovvero 15.86.

Il parametro math.anx è invece il coefficiente angolare delle due rette. Avendo in comune lo stesso coefficiente angolare, le due rette saranno parallele.

L’aver inserito nel modello il fattore genere, quindi, ha permesso di identificare due rette, parallele ma con diversa intercetta. Dato che i parametri sono tutti significativi, possiamo dire che la prestazione in matematica decresce all’aumento dell’ansia provata verso le situazioni in cui le ablità matematiche vengono sottoposte a valutazione, ma le femmine hanno di base un punteggio nelle abilità matematiche più elevato rispetto a quello dei maschi.

Di seguito è riportata una rappresentazione grafica del modello mod2, con le due rette che incrociano l’asse delle ordinate nei punti 32.23 (Intercept) e 15.86 (Intercept + genderM).

Scatterplot modello additivo

Interazione fra i predittori


Inserire nel modello l’interazione fra i due predittori significa permettere alla due rette di avere non solo diversa intercetta, ma anche diverso coefficiente angolare, per cui le rette presenteranno una pendenza diversa. Proviamo quindi a inserire nel modello anche l’effetto d’interazione.


> mod3 <- lm(math.perf ~ math.anx * gender, data=math) > summary(mod3)

Call:
lm(formula = math.perf ~ math.anx * gender, data = math)

Residuals:
Min 1Q Median 3Q Max
-10.0022 -3.1375 0.0373 2.6341 10.8452

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.8760 1.8000 21.043 < 2e-16 *** math.anx -1.2053 0.1047 -11.511 < 2e-16 *** genderM -26.5177 2.2643 -11.711 < 2e-16 *** math.anx:genderM 0.7332 0.1536 4.772 4.64e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.958 on 136 degrees of freedom Multiple R-squared: 0.7751, Adjusted R-squared: 0.7702 F-statistic: 156.3 on 3 and 136 DF, p-value: < 2.2e-16 [/code] Per prima cosa bisogna segnalare che, in questo nuovo modello, le intercette sono cambiate. Ora l'intercetta delle femmine è pari a 37.88, mentre l'intercetta dei maschi è pari a 37.88 + (-26.52) = 11.36. Inoltre, il parametro math.anx, che in precedenza era il coefficiente angolare comune a entrambe le rette, ora rappresenta il coefficiente angolare della retta del gruppo di riferimento (le femmine).

Oltre a questo, nel modello compare un nuovo parametro: math.anx:genderM. Questo parametro rappresenta la quantità da aggiungere al coefficiente angolare del gruppo di riferimento per ottenere il coefficiente angolare del gruppo di confronto (i maschi).

Quindi, dato che il coefficiente angolare delle femmine è pari a -1.21, il coefficiente angolare dei maschi sarà dato da -1.21 + 0.73, ovvero -0.48. I maschi, dunque, avranno non solo un’intercetta più bassa, ma anche un coefficiente angolare più vicino a zero.

Graficamente, la situazione è questa:

Scatterplot modello additivo

Insomma: rispetto alle femmine, i maschi presentano un’intercetta significativamente inferiore e un coefficiente angolare significativamente più vicino a zero. La relazione tra math.anx e math.perf, quindi, varia a seconda che l’individuo in questione sia un maschio oppure una femmina. Le femmine, rispetto ai maschi, hanno delle prestazioni in matematica più elevate quando l’ansia è pari a zero, ma subiscono maggiormente l’effetto negativo dell’ansia, che fa ne calare maggiormente la prestazione. L’ansia, quindi, influenza maggiormente la prestazione delle femmine.

Stepwise regRession: Parte seconda

Riassunto delle puntate precedenti. Il mondo è tutto ciò che accade (Wittgenstein, il tizio strano Viennese). Tra le cose che accadono, alcune sembrano collegate tra loro. Molte di queste relazioni tra gli accadimenti ce le mettiamo noi (tutte, a dire il vero). Ma alcune di queste relazioni non sono mere proiezioni di un pensiero magico o superstizioso, attratto dalle coincidenze. Alcune di queste relazioni mostrano un marcata ricorrenza e coerenza.

La misurazione degli elementi che compartecipano di queste relazioni esita spesso in regolarità che si è tentati di estrapolare come “legge”.  A questo servono le regressioni: a consentire la verifica delle relazioni tra le misurazioni di elementi che appaiono, all’ispezione, collegati tra loro. Alla base di una regressione c’è un’ipotesi della relazione sottesa tra gli elementi in valutazione. Questa relazione può essere espressa in formule, o algoritmi (regole di calcolo), la cui adeguatezza (fit) è verificata dalla soluzione dell’equazione stessa in base ai dati ricavati dall’osservazione.

Quando si rilevano delle regolarità nelle associazioni tra variabili, solitamente si tratta di fenomeni che hanno base biologica, meccanica o “naturale”. Fatevi un giro nei dataframe disponibili in R per farvi un’idea.

Queste relazioni sono ricorrenti in dipendenza da dinamiche che sono intrinseche al sistema misurato.

È qui che nasce il concetto di causalità (dare un’occhiata quiqui, tanto per gradire). Insomma, quando un sistema è integrato, la variazione in una sua parte implica solitamente variazioni in parti alla prima collegate. I petali e i sepali in un fiore sono dipendenti dal “metabolismo”, una cosa magica che sostiene la vita dei fiori (e di ogni essere animato), l’insieme delle reazioni chimiche che si realizzano in un dato organismo, reazioni che implicano costrizioni e costi, che si traducono in variazioni concordi tra le componenti dell’organismo stesso. A conti fatti – e mi si passi il calembour – risulta vera e verificabile l’aforisma di Lavoisier: “Nulla si crea, nulla si distrugge, tutto si trasforma”.

La regressione stepwise, dunque. Si tratta di un sistema per semplificare una regressione multipla. La regressione stepwise è un metodo di selezione delle variabili indipendenti allo scopo di selezionare un set di predittori che abbiano la migliore relazione con la variabile dipendente.

Esistono vari metodi di selezione delle variabili.

Il metodo forward (in avanti) inizia con un modello ‘nullo’ nel quale nessuna variabile tra i predittori è selezionata; nel primo step viene aggiunta la variabile con l’associazione maggiormente significativa sul piano statistico. Ad ogni step successivo è aggiunta la variabile con la maggiore associazione statisticamente significativa tra quelle non ancora incluse nel modello, ed il processo prosegue sino a quando non vi è più variabile con associazione statisticamente significativa con la variabile dipendente.

Il metodo backward (all’indietro) inizia con un modello che comprende tutte le variabili e procede, step by step, ad eliminare le variabili partendo da quella con l’associazione con la variabile dipendente meno significativa sul piano statistico.

Il processo stepwise fa avanti e indietro tra i due processi, aggiungendo e rimuovendo le variabili che, nei vari aggiustamenti del modello (con aggiunta o re-inserimento di una variabile) guadagnano o perdono in termini di significatività.

L’uso del metodo non è senza conseguenze (vedere oltre), ma può consentire una esplorazione di relazioni che a una prima ispezione possono non essere evidenti.

Esempio, just for fun.

# carichiamo un dataset disponibile in R:

data(airquality)
attach(airquality)
dim(airquality)
head(airquality)

# Qual è la relazione tra un insieme di variabili
# climatiche e il livello di ozono nell’aria?
# Consideriamo un modello ‘additivo’, nel quale non
# si realizzano interazioni tra le variabili:

model1 <- lm(Ozone ~ Temp + Wind + Solar.R)
summary(model1)

# Consideriamo ora un modello nel quale è ammessa una
# interazione tra le variabili (molto più plausibile):

model2 <- lm(Ozone ~ Temp * Wind * Solar.R)
summary(model2)

Nell’esempio qui sopra, il cui output può essere visualizzato eseguendo il codice nel terminale di R, vengono costruiti due modelli: uno additivo (model1), che non considera le interazioni fra i tre predittori Temp, Wind e Solar.R, e uno che considera queste interazioni (model2).

Come appare immediatamente evidente, nel modello2, in cui si tiene conto delle interazioni tra le variabili (modello più che plausibile), la significatività delle associazioni tra i predittori (le variabili indipendenti) e la variabile dipendente (nel nostro caso, livelli di ozono nell’aria) scompare. Verosimilmente, questo fenomeno è causato da problemi di multicollinearità che si presentano inserendo nel modello i termini d’interazione.

Inserendo le interazioni, la varianza spiegata dal modello aumenta: l’indice R-squared passa da 0.6059 a 0.6883; non si tratta di un grande incremento, ma è pur sempre un incremento (R² aumenta dell’8%). In questo secondo modello, nonostante in apparenza nessuna delle singole variabili indipendenti sembri avere alcuna relazione con la variabile dipendente, la statistica F complessiva resta significativa:

nel modello 1 abbiamo: F(3,107) = 54.83,  p-value: < 2.2e-16
nel modello 2 abbiamo: F(7,103) = 32.49,  p-value: < 2.2e-16

Ci fidiamo del risultato negativo? No.

Proviamo ad applicare un modello stepwise.

Come si fa? Si applica il comando step al modello.

# La funzione testa le interazioni ed elimina quelle
# che incrementano l'AIC:

model3 <- step(model2)
summary(model3)

# La varianza spiegata è rimasta più alta che nel modello
# semplicemente ‘additivo’ (R² =  0.6677)

Bella differenza, vero? Ora il modello restituisce alcune interazioni come significative, ma soprattutto tornano ad essere significative le relazioni già osservate nel più semplice modello additivo.

Che cosa è l’AIC? L?AIC è l’Akaike Information Criterion, un indice che consente di comparare due modelli e che consente di valutare variazioni nell‘adattamento di un modello per modifiche dello stesso (vedere qui per dettagli).

L‘AIC può essere interpretato solo in ottica comparativa: il modello migliore è sempre quello con l’AIC più basso. Quindi, quello a cui si deve badare non è il valore dell’AIC in sé, ma la differenza in AIC tra i modelli, che può essere valutata in questo modo (adattabile ad ogni successione di modelli):

# Create un vettore che comprenda tutti i valori AIC
# dei modelli (qui esempio con tre modelli):

x <- c(AIC(model1), AIC(model2), AIC(model3))

# Calcolate la differenza:

delta <- x - min(x)
delta

# Vedete da voi che il modello 3 è quello con differenza
# pari a 0, cioè, è il modello migliore.

# Potete anche derivare il ‘peso’ (weight) dell’AIC nei vari
# modelli, essendo il ‘peso’ relativo dell’AIC  di un modello
# rispetto agli altri una evidenza in suo favore:

# Calcolate la verosimiglianza (likelihood) dei modelli

L <- exp(-0.5 * delta)

# Akaike weights:

w <- L/sum(L)
w

# Si può anche utilizzare un altro indice, il Bayesian
# information criterion (BIC), che, a differenza dell'AIC,
# penalizza maggiormente modelli più complessi:

x_BIC <- c(BIC(model1), BIC(model2), BIC(model3))
delta_BIC <- x_BIC - min(x_BIC)
delta_BIC
L_BIC <- exp(-0.5 * delta_BIC)

# BIC weights:
w_BIC <- L/sum(L_BIC)
w_BIC

Ricca messe di informazioni sul tema qui. Dal sito, curato da Dolph Schluter, ho tratto il metodo rapido per calcolare l’AIC e il BIC.

Si diceva dei limiti di questa procedura. Il critico più agguerrito dell’uso della procedura stepwise è Frank Harrell, un Big Kahuna della statistica delle regressioni (non avete visto il film? Male! Monologo finale).

Le obbiezioni di Harrell sono da prendere con beneficio di inventario, ma non andrebbero trascurate:

1. R² values are biased high.
2. The F and χ² test statistics do not have the claimed distribution.
3. The standard errors of the parameter estimates are too small.
4. Consequently, the confidence intervals around the parameter estimates are too narrow.
5. p-values are too low, due to multiple comparisons, and are difficult to correct.
6. Parameter estimates are biased high in absolute value.
7. Collinearity problems are exacerbated.

Tuttavia, i benefici della stepwise regression sono indubitabili. Come nell’esempio riportato sopra, consente di identificare relazioni oscurate dall’inserimento di tutte le variabili nel modello.

Va ricordato che una regressione è un algoritmo che consente di stimare una variabile dipendente a partire da un insieme di predittori. Talvolta la stima di tale variabile può essere complicata e richiedere molte risorse; più risorse sono richieste maggiore sarà il tempo impiegato dall‘algoritmo. Inoltre, è necessario conoscere in anticipo la relazione matematica che lega i predittori alla variabile dipendente.

L’esito di una regressione sono dei parametri in base ai quali ‘correggere’ i punteggi dei predittori per ottenere una approssimazione della variabile dipendente, approssimazione che è tanto più precisa quanto migliore il modello (vedere post precedenti sulla valutazione della adeguatezza dei modelli – qui e qui).

Esempio. Poniamo che voi vogliate stimare il valore di una sostanza tossica liberata da una fabbrica. Quando i livelli di quella sostanza tossica sono elevati, c’è rischio di contaminazione. È costoso dosare la sostanza, e inoltre vi serve sapere in anticipo se la concentrazione nell’aria sta salendo verso il livello di guardia. Sapete che alcuni indici sono legati al rischio di innalzamento delle concentrazioni della sostanza tossica: temperatura, pressione nelle tubature, grado di umidità dell’aria. Queste variabili sono state stimate per un periodo di otto mesi, e una regressione stepwise le ha estratte come predittori efficienti della concentrazione della sostanza tossica.

I coefficienti di questi predittori sono:
– Temperatura: βT = 2.45
– Pressione: βP = 11.11
– Umidità: βU = –3.53
Mentre il valore dell‘intercetta è: α = –13.81.

Il modello di stima dei livelli di concentrazione della sostanza tossica è:

Sostanza_tossica = α + βT * Temperatura + βP * Pressione + βU * Umidità

Con un sistema di sensori ottenete, ogni dieci minuti, i valori di Temperatura, Pressione e Umidità.

In base ai coefficienti calcolati dal vostro modello, potete calcolare i valori presunti della sostanza tossica in base ai predittori semplicemente applicando la vostra formula:

Sostanza_tossica = – 13.81 + (2.45 x [valore misurato di Temperatura]) + (11.11 x [valore misurato di Pressione]) + (–3.53 x [valore stimato di Umidità])

Quando il valore predetto raggiunge una soglia di rischio, a quel punto effettuate la misurazione dei livelli della sostanza tossica. Poiché lo fate solo quando si raggiunge una soglia di rischio, risparmiate dei bei soldini…

Pensate a quello che si può fare con un buon modello matematico di predizione (che significa, buone analisi ma anche buoni dati), un pacchetto di sensori, e un processore Arduino

Alternative alla stepwise regression? Nessuna facile da capire. Ne parleremo in un prossimo post.

Intanto, per cominciare, date un’occhiata qui.

That’s all, folks!

Antonello Preti

Stay Tuned for our next episode!

Should I stay or should I go? Stepwise regRession – parte prima

Ah, i favolosi anni Ottanta, con le minigonne a palloncino, i capelli cotonati, la pelle liscia senza (ancora) tatuaggi… e i Clash, appunto…

L’argomento della puntata è la regressione a passi cauti (stepwise… OK, I’m joking…). L’argomento è controverso, e lo prenderemo alla lontana.

A dire il vero, R è un sistema nel quale la regressione è “cucinata” in tutte le salse. Se esiste una qualche incertezza sul perché R si chiami proprio così (forse dalle iniziali dei nomi dei primi programmatori, Ross Ihaka e Robert Gentleman), si può avanzare la ragionevole ipotesi che R stia per Regression. In effetti, una larga fetta di pacchetti statistici in R sviluppa, in via diretta o indiretta, una variazione sul tema della regressione: lineare, non-lineare, multivariata, robusta o quello che volete voi.

Non è facilissimo definire il concetto di regressione. Certo, su Wikipedia si trova che “La regressione formalizza e risolve il problema di una relazione funzionale tra variabili misurate sulla base di dati campionari estratti da un’ipotetica popolazione infinita” (vedere per credere). Tuttavia il concetto di regressione definisce più una filosofia di vita che un modello matematico. Al fondo delle credenze incentrate sui modelli di regressione è l’idea che il mondo là fuori abbia una razionalità intrinseca, che esista la possibilità di formulare previsioni sull’andamento degli eventi, e, insomma, la fede incrollabile che “God does not play dice” (ma, a sentire Stephen Hawking, “God still has a few tricks up his sleeve“, insomma, si tiene qualche asso nella manica…).

Al dunque, “l’analisi della regressione è una tecnica usata per analizzare una serie di dati che consistono in una variabile dipendente e una o più variabili indipendenti” (Wikipedia, of course). Sfortunatamente, questo nucleo incontrovertibile di ogni modello di regressione che si rispetti è quello sul quale si appunta la critica – per nulla prescindibile – di David Hume: “Che domani il sole non sorgerà, è una proposizione non meno intelligibile, e non implica più contraddizione, dell’affermazione secondo cui esso sorgerà; tenteremmo invano, perciò, di dimostrarne la falsità. Se fosse dimostrativamente falsa, essa implicherebbe una contraddizione e non potrebbe mai essere concepita dalla mente in modo distinto”.

La relazione tra le variabili, direbbe il buon Jacques, ce la mettiamo noi, è una relazione simbolica. I “dati” sono astrazioni, l’esito di operazioni di misurazione, che risentono non solo dell’errore implicito in ogni misurazione (e di ciò si tiene conto nei modelli matematici di regressione), ma anche delle “credenze” implicate nell’atto stesso della misurazione. Se crediamo di poter predire l’andamento dei raccolti dalle condizioni di semina e da quelle climatiche è perché coltiviamo (mi si perdoni il calembour) l’idea che gli eventi/fatti “semina”, “condizioni metereologiche” e “raccolto” siano tra loro collegati.

Il concetto di regressione implica – infatti – una qualche relazione tra una “causa” e un “effetto”, che il modello in valutazione ha il compito di stimare. Nella sua forma più semplice, la regressione è un confronto tra le misure di un fatto/evento lungo una dimensione e le misure del medesimo fatto/evento lungo una dimensione collegabile alla prima.

Esempio, esempio, reclama il mio fedele pubblico di 25 lettori

La regressione stepwise è applicata alla regressione multipla. Partiamo dunque da qui.

Conciossiacosaché si voglia misurare la lunghezza e la larghezza dei petali in un fiore e valutare se dette misure predicano la lunghezza dei sepali dello stesso fiore.

# caricate il dataframe (già disponibile in R)
data(iris)

# check preliminare sul dataset (approfondisci: www.insular.it/?p=1698)
dim(iris)
head(iris)

# per accedere direttamente alle variabili (approfondisci: www.insular.it/?p=2245)
attach(iris)

L’ipotesi è che in un insieme di fiori di varia specie la lunghezza e la larghezza dei petali predica la lunghezza dei sepali dei medesimi fiori, cioè, detta in regressionese, che “lunghezza e larghezza dei petali” siano le variabili indipendenti che predicono la distribuzione della variabile dipendente “lunghezza dei sepali”:

In R, la funzione che specifica una relazione lineare tra due variabili e che di conseguenza adatta un modello di lineare è lm. La relazione tra variabile dipendente (la y nei grafici cartesiani) e la o le variabili indipendenti (la x nei grafici cartesiani) è specificata dalla tilde “~”:

fit1 <- lm(Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
summary(fit1)

Come output, ci verranno restituite le principali statistiche descrittive sui residui, il test di significatività sui parametri e infine gli indici di adattamento e di significatività dell’intero modello. La tabellina “Coefficients” riporta, per ogni parametro, il valore (Estimate), l’errore standard (Std. Error) e i valore della statistica t (t value) e il p-value (Pr(>|t|)).

Call:
lm(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.18534 -0.29838 -0.02763  0.28925  1.02320 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.19058    0.09705  43.181  < 2e-16 ***
Petal.Length  0.54178    0.06928   7.820 9.41e-13 ***
Petal.Width  -0.31955    0.16045  -1.992   0.0483 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4031 on 147 degrees of freedom
Multiple R-squared:  0.7663,	Adjusted R-squared:  0.7631 
F-statistic:   241 on 2 and 147 DF,  p-value: < 2.2e-16

Il modello sopra illustrato è additivo, si suppone cioè che le due variabili indipendenti si sommino nel loro impatto ma non interagiscano. Vogliamo formulare l’ipotesi che le due variabili indipendenti interagiscano in modo moltiplicativo? Usiamo il segno “*” anziché il +:

fit2 <- lm(Sepal.Length ~ Petal.Length * Petal.Width, data = iris)
summary(fit2)

L’output sarà il seguente:

Call:
lm(formula = Sepal.Length ~ Petal.Length * Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.00058 -0.25209  0.00766  0.21640  0.89542 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               4.57717    0.11195  40.885  < 2e-16 ***
Petal.Length              0.44168    0.06551   6.742 3.38e-10 ***
Petal.Width              -1.23932    0.21937  -5.649 8.16e-08 ***
Petal.Length:Petal.Width  0.18859    0.03357   5.617 9.50e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3667 on 146 degrees of freedom
Multiple R-squared:  0.8078,	Adjusted R-squared:  0.8039 
F-statistic: 204.5 on 3 and 146 DF,  p-value: < 2.2e-16

Interessante: la varianza spiegata (Adjusted R-squared) aumenta, passando dal 76.6% all’80.8% e cambia anche l’impatto delle singole variabili indipendenti (i valori stimati per i coefficienti).

Controlliamo la validità del modello con le diagnostiche di regressione:

# regression diagnostics
par(mfcol = c(2,4))
plot(fit1,main="modello additivo")
plot(fit2, main="modello con interazione")

Regression diagnostics

A occhio, sembrerebbe che il modello moltiplicativo, che prevede l’interazione delle due variabili predittive, vada meglio, soprattutto la curva dei residui cade maggiormente sulla linea orizzontale.

Possiamo anche verificare che il secondo modello sia effettivamente migliore del primo facendo uso della funzione anova:

anova(fit2,fit1)

L’output evidenzia una differenza significativa fra i due modelli, di conseguenza sembra davvero che il secondo sia migliore del primo.

Analysis of Variance Table

Model 1: Sepal.Length ~ Petal.Length * Petal.Width
Model 2: Sepal.Length ~ Petal.Length + Petal.Width
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    146 19.637                                  
2    147 23.881 -1   -4.2441 31.556 9.502e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Siamo però sicuri che vada tutto bene? Prima di proseguire dovremmo verificare che siano rispettati alcuni requirements della regressione lineare multipla. Per esempio, bisogna verificare la normalità dei residui; per effettuare questa operazione utilizziamo la funzione qqPlot della libreria car che genera un quantile-quantile plot per i residui studentizzati.

library(car)
# qq-plot for studentized resid
par(mfcol = c(1,2))
qqPlot(fit1, main="QQ Plot – modello additivo")
qqPlot(fit2, main="QQ Plot – modello con interazione")

QQplot for studentiized residuals

Fin qui, tutto bene. Vi sono violazioni della omoscedasticità dei residui? Anche per queste analisi serve il pacchetto “car”:

library(car)
# non-constant error variance test
ncvTest(fit1)
ncvTest(fit2)
# plot studentized residuals vs. fitted values
par(mfcol = c(1,2))
spreadLevelPlot(fit1, main=paste("Spread-Level - modello additivo"))
spreadLevelPlot(fit2, main=paste("Spread-Level - modello con interazione"))

Ecco l’output:

> ncvTest(fit1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.121054    Df = 1     p = 0.2896917

> ncvTest(fit2)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6024558    Df = 1     p = 0.4376425

QQplot for studentiized residuals

Sembra che il modello si adatti con con pochi problemi. Ci possiamo accontentare? Non ancora… dobbiamo valutare l’eventuale presenza di multicollinearità. Il metodo più semplice è sicuramente quello di valutare l’entità della correlazione tra i predittori:

cor(Petal.Length, Petal.Width)

L’indice di Bravais-Pearson risulta pari a 0.96, decisamente troppo elevato. I due predittori dovrebbero infatti essere poco correlati, altrimenti le stime dei parametri risultano molto instabili. Per individuare la presenza di multicollinearità è spesso utilizzato uno stimatore noto come variance inflation factors (VIF), cioè “fattori di crescita della varianza”. Valori superiori a 1 indicano rischio di multicollinearità. Valori superiori a 10 suggerisco seri problemi di multicollinearità.

# Evaluate Collinearity
# modello additivo
vif(fit1) # variance inflation factors
sqrt(vif(fit1)) > 2 # problem?

# modello con interazione
vif(fit2) # variance inflation factors
sqrt(vif(fit2)) > 2 # problem?

Turista fai da te? No Pre-regression test? Ahi! Ahi! Ahi! Ahi!

Prima di programmare una regressione lineare è utile verificare se sono soddisfatte le assunzioni a priori per una regressione lineare. In R è disponibile un pacchetto che esegue l’intero ventaglio di valutazioni. Il pacchetto si chiama gvlma, che va preventivamente installato. Questo è il codice di valutazione:

library(gvlma)

# modello additivo
gvmodel1 <- gvlma(fit1)
summary(gvmodel1)

# modello con interazione
gvmodel2 <- gvlma(fit2)
summary(gvmodel2)

Il pacchetto produce anche un insieme di plot che definiscono la validità del modello.

# esempio per il modello addittivo
plot(gvmodel1)

Regression diagnostics with gvlma

Alcuni di questi plot sono già riportati dalle diagnostiche della regressione già viste, e queste diagnostiche sono già state spiegate in un post precedente. Altri sono nuovi.

Cosa significano i nuovi plot? Come direbbe un pirata in uno di quei magnifici libri ottocenteschi che nessuno scrive più (Treasure Island, ad esempio, da leggere rigorosamente nell’edizione integrale, senza tagli): “Ch’io possa essere dannato se lo so!...”

That's all, folks!

Antonello Preti

Stay Tuned for our next episode!