Archivio tag: modelli lineari

Gestire le misure ripetute con i modelli a effetti misti

In un post di qualche tempo fa ho parlato di come realizzare con R la cosiddetta “ANOVA a misure ripetute”, ovvero un’analisi della varianza in cui sono presenti uno o più fattori within subjects. Le soluzioni proposte prevedevano l’uso delle funzioni aov oppure in alternativa della combinata lm + Anova (dal pacchetto car).

Per concludere il capitolo “misure ripetute”, in questo post vedremo rapidamente come risolvere la faccenda utilizzando i modelli a effetti misti. Questi, seppure molto articolati dal punto di vista teorico, sono facili da implementare in R e consentono di analizzare le misure ripetute senza dover ricorrere a bizzarri e poco trasparenti giri di codice.

Di seguito non parlerò della teoria che sta dietro il modello, un po’ perché descriverla rigorosamente mi richiederebbe un tempo esagerato, un po’ perché su internet si trova già tanto materiale, certamente migliore di quello che potrei produrre io in questo piccolo spazio. L’interpretazione dei risultati, quindi, sarà vista sempre nell’ottica dell’ANOVA.

Per prima cosa, se non l’avete ancora fatto, vi consiglio di leggere il post sull’ANOVA, almeno nella sua parte introduttiva.

I dati ormai li conosciamo. Nel repository GitHub fakedata è presente il file sonno.zip, che contiene il dataset che useremo; il dataset che ci interessa è contenuto nel file “sonno-long.csv”, che riporta i dati di un ipotetico esperimento:

In un laboratorio di psicofisiologia del sonno si stanno testando gli effetti di una nuova psicoterapia (PT) e di un trattamento di tipo farmacologico (FT) per la cura dell’insonnia. Sono state selezionate 10 persone insonni con problemi nella fase di addormentamento: 5 vengono trattate con PT e 5 con FT. Il giorno prima dell’inizio del trattamento (tempo 0), a 30 giorni e infine a 60 giorni viene registrato il numero di ore dormite da ogni individuo.

Importiamo i dati in R:

> sonno <- read.csv2("sonno-long.csv")
> str(sonno)
'data.frame':	30 obs. of  4 variables:
 $ sogg : int  1 1 1 2 2 2 3 3 3 4 ...
 $ tempo: int  0 30 60 0 30 60 0 30 60 0 ...
 $ tratt: Factor w/ 2 levels "FT","PT": 2 2 2 2 2 2 2 2 2 2 ...
 $ ore  : num  5.6 3.3 8.9 4.8 4.7 7.9 3.9 5.3 8.5 5.8 ...

Prima di partire, dobbiamo convertire le variabili sogg e tempo in tipo fattore:


> sonno$sogg <- as.factor(sonno$sogg) > sonno$tempo <- as.factor(sonno$tempo) [/code]

Gli effetti misti


Un modello lineare a effetti misti è un modello lineare che include al suo interno contemporaneamente sia effetti fissi che effetti casuali. Per “effetto” si intende una variabile che si suppone abbia una qualunque influenza sulla variabile dipendente. Nell’esperimento descritto poco sopra, si suppone che il numero di ore dormite vari in funzione del tipo di trattamento, del momento in cui è stata fatta la rilevazione e anche da una componente individuale che rende ogni soggetto unico, con delle peculiarità che lo distinguono dagli altri e che lo portano a dormire di più o di meno a prescindere dal trattamento ricevuto e dal momento della rilevazione.

Questi tre effetti, ovvero il trattamento, il tempo e il soggetto, sono rappresentati da variabili di tipo fattore; questi fattori, però, hanno una diversa natura: se trattamento e tempo sono fissi, il soggetto è casuale. Ho già parlato di questa distinzione in un vecchio post, ma forse è meglio rinfrescarsi un po’ le idee.

Un fattore è detto fisso (fixed) se i suoi livelli sono determinati a tavolino e devono essere per forza quelli e non altri, altrimenti il fattore verrebbe snaturato (es. trattamento). Tra i fattori fissi rientrano anche tutti quei fattori i cui livelli considerati sono tutti quelli possibili presenti nella popolazione (es. sesso). Un fattore è fisso a prescindere dall’essere between o within.

Trattamento e tempo sono fattori fissi, il primo between e il secondo within.

Un fattore è detto casuale (random) se i suoi livelli sono solo un campione casuale di tutti i possibili livelli presenti nella popolazione (es. soggetto). In questo caso non si è interessati a valutare le differenze fra gli specifici livelli, ma la variabilità fra essi. A un fattore casuale possono essere aggiunti livelli al fine di aumentare la potenza dell’analisi statistica, mentre quest’operazione non avrebbe senso su un fattore fisso, i cui livelli sono, appunto, in numero fisso.

L’ANOVA considera tutti i fattori come se fossero fissi, anche il soggetto. Diversamente, i modelli a effetti misti sono in grado di distinguere i fattori a seconda della loro natura. Si tratta di un aspetto che, per la verità, adesso non approfondiremo più di tanto: pragmaticamente, ci limiteremo a mimare i risultati dell’ANOVA con i modelli a effetti misti.

Già: mimare l’ANOVA, avete capito bene. I modelli a effetti misti, infatti, rappresentano una cornice procedurale molto ampia e l’ANOVA può essere vista come un suo caso particolare. Pur trattando i dati in maniera diversa, di fatto ANOVA e modelli a effetti misti possono portare a risultati analoghi. Tuttavia, i modelli a effetti misti offrono più libertà allo statistico, libertà che l’ANOVA non offre.

La funzione lme


Principalmente, esistono in R due pacchetti per adattare i modelli a effetti misti: nlme e lme4. Il primo è già distribuito con la versione base di R, il secondo no. Per alcuni versi i due si sovrappongono, ma per altri versi questi sono complementari. Qui utilizzeremo il primo.

Il pacchetto nlme (NonLinear Mixed Effects) mette a disposizione la funzione lme, che consente di adattare un modello lineare a effetti misti.


> library(nlme)
> mod <- lme(ore ~ tratt * tempo, random=~1|sogg, data=sonno) [/code] La sintassi è molto simile a lm, con la differenza che bisogna specificare nell'argomento random il fattore casuale. La notazione 1|sogg indica che l’intercetta deve variare tra i soggetti; questo aspetto, che mi rendo conto detto così potrebbe apparire come misterioso, risulterà più chiaro approfondendo, anche solo a livello superficiale, gli aspetti teorici del modello.

La funzione anova può essere utilizzata per ottenere i test di significatività sui fattori fissi:


> anova(mod)
numDF denDF F-value p-value
(Intercept) 1 16 1471.3992 <.0001 tratt 1 8 8.5321 0.0193 tempo 2 16 42.1280 <.0001 tratt:tempo 2 16 7.7107 0.0045 [/code] Questi risultati sono identici a quelli dell'ANOVA, con la differenza fondamentale che... è stato decisamente più facile ottenerli! (basta guardare il codice utilizzato nel precedente post).

Vorrei chiudere con una nota. Non c’è un accordo unanime su come vadano calcolati i p-value nei modelli a effetti misti, tanto che il prof. Douglas Bates, autore sia di nlme che di lme4, per lme4 ha preferito non fornirli in output (qui spiega il perché). Si tenga quindi in considerazione che i p vanno presi un po’ con le pinze e che l’approccio migliore per operare sarebbe l’uso del confronto fra modelli, che chi ha seguito l’ultimo corso InsulaR conosce bene!

Analisi della varianza a misure ripetute

Diverso tempo fa ho scritto due post riguardanti le tipologie di fattori nei disegni sperimentali e la dimensione temporale negli esperimenti. In entrambi i casi ho citato l’analisi della varianza (ANOVA) a “misure ripetute”, che altro non è che un’analisi della varianza in cui, a causa della presenza di uno o più fattori entro i soggetti (within subjects), si utilizzano differenti componenti d’errore per valutare la significatività di ogni fattore.

In questo post presenterò due strategie per eseguire con R l’ANOVA in presenza di fattori di tipo within. Prima di proseguire nella lettura, però, è bene tenere a mente tre cose.

Punto 1. Non è scopo di questo post spiegare la teoria che sta dietro l’ANOVA a misure ripetute, peraltro ormai superata dai modelli a effetti misti. Si assume quindi che, chi legge, già conosca la teoria sottostante e sia semplicemente interessato a mettere in pratica l’analisi utilizzando R.

Punto 2. Come anticipato, l’ANOVA a misure ripetute è ormai considerata una tecnica obsoleta. Nonostante questo, continua a comparire nelle pubblicazioni scientifiche; personalmente, ritengo che si tratti semplicemente di una fisiologica e più che comprensibile difficoltà di molte discipline nell’aggiornare il proprio bagaglio di tecniche statistiche, ma è solo questione di tempo prima che scompaia definitivamente. Ritengo che quella dell’ANOVA a misure ripetute sia solo una lenta agonia che la condurrà ben presto all’oblio. C’è già chi si interroga se non sia il caso di seppellirla definitivamente (vedi qui). A prescindere da quelle che sono le mie idee, dato che continua a essere utilizzata e talvolta mi viene richiesta, ho pensato di scrivere questo tutorial, anche perché realizzarla con R è tutt’altro che banale.

Punto 3. L’ho già scritto nei post precedenti e qui lo ribadisco: disapprovo il termine “misure ripetute” come utilizzato in questo contesto. Vi rimando sempre ai due post citati all’inizio (qui e qui) per una breve spiegazione.

Bene, ora possiamo cominciare.

I dati


Nel repository GitHub fakedata è presente il file sonno.zip, che contiene il dataset che useremo, del quale avevamo già fatto la conoscenza in passato. Il file che ci interessa è “sonno-long.csv”, che riporta i dati di un ipotetico esperimento:

In un laboratorio di psicofisiologia del sonno si stanno testando gli effetti di una nuova psicoterapia (PT) e di un trattamento di tipo farmacologico (FT) per la cura dell’insonnia. Sono state selezionate 10 persone insonni con problemi nella fase di addormentamento: 5 vengono trattate con PT e 5 con FT. Il giorno prima dell’inizio del trattamento (tempo 0), a 30 giorni e infine a 60 giorni viene registrato il numero di ore dormite da ogni individuo.

Una volta decompresso l’archivio e posizionato il file nella cartella di lavoro, possiamo importare il dataset e visualizzarne la struttura.

> sonno <- read.csv2("sonno-long.csv")
> str(sonno)
'data.frame':	30 obs. of  4 variables:
 $ sogg : int  1 1 1 2 2 2 3 3 3 4 ...
 $ tempo: int  0 30 60 0 30 60 0 30 60 0 ...
 $ tratt: Factor w/ 2 levels "FT","PT": 2 2 2 2 2 2 2 2 2 2 ...
 $ ore  : num  5.6 3.3 8.9 4.8 4.7 7.9 3.9 5.3 8.5 5.8 ...

La variabile dipendente è il numero di ore dormite (ore), mentre i predittori sono tempo e tratt; inoltre, la variabile sogg codifica i soggetti.

Prima di partire, dobbiamo assicurarci che le variabili da includere come predittori nel modello ANOVA (soggetti compresi) siano di tipo fattore. Come possiamo notare, sogg e tempo sono codificate con dei numeri e per questo, nell’operazione di lettura dati, R gli ha assegnato il tipo numerico. Queste due variabili vanno quindi convertite in fattore:

> sonno$sogg <- as.factor(sonno$sogg)
> sonno$tempo <- as.factor(sonno$tempo)
&#91;/code&#93;

Giusto a scopo descrittivo, utilizzando il <a href="http://www.insular.it/?p=2512">comando tapply</a>, calcoliamo la media di ore dormite dai soggetti sottoposti ai due trattamenti nei tre momenti di rilevazione:

[code language="R"]
> with(sonno, tapply(ore, list(tratt,tempo), mean))
      0   30   60
FT 4.20 4.58 5.90
PT 4.72 4.30 8.08

Almeno a livello medio, tra l'inizio e il termine del periodo di monitoraggio il numero di ore dormite cresce per tutti i soggetti, ma per PT si osserva un clamoroso picco a 60 giorni che, invece, nel caso FT appare molto meno marcato.

La funzione aov


Il metodo più comune per svolgere un'ANOVA a misure ripetute con R è usare la funzione aov. Questa è molto simile alla più comune funzione lm, ma, a differenza di questa, è stata specificatamente studiata per l'ANOVA e consente di ripartire la varianza d'errore a seconda dell'effetto considerato (al contrario, la funzione lm considera un'unica varianza d'errore comune a tutti gli effetti).

L'istruzione aov, come la maggior parte delle funzioni statistiche di R, richiede dati disposti in formato long e il dataset sonno è già disposto correttamente.

Per realizzare l'ANOVA a misure ripetute bisogna passare ad aov una formula dove viene esplicitata la struttura di dipendenza delle variabili (ovvero che “ore” dipende da “tratt”, “tempo” e dalla loro interazione). Inoltre, nella formula bisogna specificare come ripartire la varianza d'errore sfruttando l'opzione Error; dentro Error bisogerà specificare che, per ogni livello del fattore “tempo”, le osservazioni sono effettuate sempre sugli stessi soggetti.

> fit <- aov(ore ~ tratt * tempo + Error(sogg/tempo), data=sonno)
&#91;/code&#93;

Quello qui presentato è un caso abbastanza semplice, in cui è presente un unico fattore within. Immaginiamo per un attimo di avere un caso decisamente più complesso, con due fattori between (bet1 e bet2), due fattori within (with1 e with2) e il fattore che identifica i soggetti (subj). In questo caso, la formula da passare al comando aov dovrà essere strutturata nel seguente modo:

<p style="background-color: white; padding: 0.2em; margin-top: 1em; margin-bottom: 1em; font-weight: bold; text-align: center;">
y ~ bet1 * bet2 * wit1 * with2 + Error( subj / (wit1 * wit2) )
</p>

Per visualizzare il risultato dell'ANOVA bisogna passare al comando <b>summary</b> l'oggetto costruito tramite aov. Gli effetti saranno suddivisi in tante tabelle a seconda della relativa sorgente d'errore. Nella prima tabella in output, la componente residua (Residuals) è data dalla varianza entro i soggetti (Error: sogg), mentre nella seconda tabella è data dall'interazione tra soggetti e tempo (Error: sogg:tempo).

[code language="R"]
> summary(fit)

Error: sogg
          Df Sum Sq Mean Sq F value Pr(>F)  
tratt      1  4.880   4.880   8.532 0.0193 *
Residuals  8  4.576   0.572                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: sogg:tempo
            Df Sum Sq Mean Sq F value   Pr(>F)    
tempo        2  43.01  21.506  42.128 4.21e-07 ***
tratt:tempo  2   7.87   3.936   7.711  0.00452 ** 
Residuals   16   8.17   0.510                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Il comando aov presenta un grosso limite al quale bisogna prestare molta attenzione: non lavora con disegni sbilanciati. La cosa veramente pericolosa è che, nonostante i messaggi di avvertimento che in questo caso potrebbero essere stampati a video, verrà comunque restituito un risultato, ma le tabelle risulteranno di fatto ininterpretabili. In questi casi, è bene procedere utilizzando il secondo metodo.

Il pacchetto car


Un altro modo di svolgere l'ANOVA a misure ripetute è combinare la funzione lm con la funzione Anova presente nel pacchetto car (Companion to Applied Regression). Si presti attenzione che in R è già presente una funzione che si chiama “anova” con la “a” minuscola; questo comando formito dal pacchetto car ha invece l'iniziale maiuscola. Il pacchetto car non è distribuito con la versione base di R, ma andrà installato a parte.

Prima di cominciare, dobbiamo ristrutturare i dati disponendoli in formato wide, disponendo le misurazioni dei tre tempi su tre diverse colonne, come se fossero variabili differenti:

> sonno.wide <- reshape(sonno, v.names="ore", timevar="tempo", idvar="sogg", direction="wide")
> sonno.wide
   sogg tratt ore.0 ore.30 ore.60
1     1    PT   5.6    3.3    8.9
4     2    PT   4.8    4.7    7.9
7     3    PT   3.9    5.3    8.5
10    4    PT   5.8    4.5    8.2
13    5    PT   3.5    3.7    6.9
16    6    FT   3.9    4.7    7.0
19    7    FT   4.1    4.9    5.6
22    8    FT   5.0    3.8    5.3
25    9    FT   4.0    4.7    5.3
28   10    FT   4.0    4.8    6.3

A questo punto, dobbiamo adattare un modello di analisi della varianza multivariata (MANOVA), considerando ore.0, ore.30 e ore.60 come tre distinte variabili dipendenti e tutti i fattori between come variabili indipendenti:

> mod <- lm(cbind(ore.0, ore.30, ore.60) ~ tratt, data=sonno.wide)
&#91;/code&#93;

Qualora nel disegno non fossero presenti fattori between, allora la formula dovrà essere scritta ponendo il valore 1 subito dopo la tilde, ovvero:

<p style="background-color: white; padding: 0.2em; margin-top: 1em; margin-bottom: 1em; font-weight: bold; text-align: center;">
cbind(y1, y2, ..., yn) ~ 1
</p>

Prima di utilizzare la funzione Anova dobbiamo costruire un nuovo data frame che specifica la struttura del piano fattoriale entro i soggetti. In pratica, questo data frame deve presentare tante colonne quanti sono i fattori between e tante righe quante sono le combinazioni dei fattori between. Nel nostro caso, essendoci un unico fattore between, la cosa è molto semplice:

[code language="R"]
> within.design <- data.frame(tempo = levels(sonno$tempo))
> within.design
  tempo
1     0
2    30
3    60

Ora possiamo eseguire l'analisi, passando al comando Anova il modello adattato attraverso la funzione lm e specificando due argomenti aggiuntivi: idata, al quale dobbiamo passare il data frame within.design appena costruito, e idesign: una formula che, sfruttando le variabili specificate nel data frame within.design, specifica il disegno entro soggetti. Inoltre, attraverso l'argomento type, possiamo indicare se vogliamo utilizzare una scomposizione dei quadrati di tipo II oppure di tipo III.

> library(car)
> av <- Anova(mod, idata = within.design, idesign = ~tempo, type = "II")
&#91;/code&#93;

Sempre utilizzando il comando <b>summary</b>, ma avendo l'accortezza di specificare che devono essere mostrati i risultati dell'ANOVA e non della MANOVA, possiamo visualizzare i risultati.

[code language="R"]
> summary(av, multivariate=FALSE)

Univariate Type II Repeated-Measures ANOVA Assuming Sphericity

                SS num Df Error SS den Df         F    Pr(>F)    
(Intercept) 841.64      1    4.576      8 1471.3992 2.343e-10 ***
tratt         4.88      1    4.576      8    8.5321   0.01926 *  
tempo        43.01      2    8.168     16   42.1280 4.208e-07 ***
tratt:tempo   7.87      2    8.168     16    7.7107   0.00452 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

            Test statistic p-value
tempo              0.85606 0.58045
tratt:tempo        0.85606 0.58045


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

             GG eps Pr(>F[GG])    
tempo       0.87417  1.938e-06 ***
tratt:tempo 0.87417   0.006847 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              HF eps   Pr(>F[HF])
tempo       1.098516 4.208025e-07
tratt:tempo 1.098516 4.520084e-03
Warning message:
In summary.Anova.mlm(av, multivariate = FALSE) : HF eps > 1 treated as 1

Come possiamo notare, il comando Anova restituisce un output più ricco rispetto ad aov, indicando anche i risultati del test di sfericità di Mauchly e i risultati corretti con il metodo di Greenhouse-Geisser. Inoltre, questa procedura non crea problemi quando si ha a che fare con disegni sbilanciati.

Qualche approfondimento sull'utilizzo del pacchetto car per l'ANOVA con fattori within è disponibile qui per il caso a una via e qui per il caso a due vie.

Infine, un altro modo di realizzare l'analisi è utilizzare la funzione ezANOVA del pacchetto ez, ma questo, magari, sarà argomento di qualche futuro post.

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.

Risolvere pRoblemi

Risolvere problemi è il metodo migliore per imparare qualcosa. Ricordate la scuola?

Con la programmazione in R è un po’ la stessa cosa. Risolvere problemi favorisce la dimestichezza con il programma e facilita l’apprendimento nell’uso dei vari pacchetti. Il suggerimento è di iniziare da cose semplici, come quelle che potete trovare in un manuale o nei siti di discussione.

Di recente l’utente Umesh Acharya ha postato un quesito su CrossValidated, un sito che si occupa di problemi di statistica. Il quesito in realtà è di natura informatica, poiché riguarda la codifica di un plot. Fino a poco tempo fa, il quesito era reperibile qui, ma è stato ormai rimosso, essendo quel genere di domande più adatte a StackOverflow, il sito gemello di CrossValidated, ma più genericamente dedicato a problemi informatici.

Domanda Umesh: “Come posso mettere tre plot in un’unica finestra?”. La domanda non è banale. Umesh ha creato una funzione che consente di “stampare” nel plot la formula di una regressione usando i valori stimati per i parametri, e per farlo ha usato ggplot2, il programma di grafica creato da Hadley Wickham, mago della programmazione in R (su InsulaR se n’era già parlato qui).

A differenza del sistema grafico base, in ggplot2 non è possibile ripartire le finestre grafiche con il comando par (che avete visto usare in precedenti post). Bisogna agire sulle griglie. E bisogna saperlo.

Dunque, abbiamo un dataset con 6 variabili, che vogliamo rappresentare a coppie di due, coordinate in una regressione. Come si rappresentano i tre grafici in un’unica finestra? Iniziamo a costruire il nostro dataset, che chiameremo ‘tom’ (toy object model):

a1 <- seq(from=1, to=10, by=1)
a2 <- seq (from=2, to=20, by=2)
b1 <- seq(from=3, to=30, by=3)
b2 <- seq(from=4, to=40, by=4)
c1 <- seq(from=5, to=50, by=5)
c2 <- seq(from=2, to=29, by=3)

# combiniamo i vettori in una matrice con il comando cbind
dataframe <- cbind(a1,a2,b1,b2,c1,c2)

# trasformiamo la matrice in un dataframe con etichette per l’intestazione delle variabili
tom <- as.data.frame(dataframe)

E adesso la bella funzione creata dal nostro Umesh per inserire il risultato della regressione nel grafico:

## function to create equation expression
lm_eqn = function(x, y, df){
    m <- lm(y ~ x, df);
    eq <- substitute(italic(y) == b %.% italic(x) + a,
    list(a = format(coef(m)[1], digits = 2),
    b = format(coef(m)[2], digits = 2)))
    as.character(as.expression(eq));
}

Chiamiamo il pacchetto:

library (ggplot2)

Creiamo i tre grafici in ggplot2:

################ a1 and a2 couple of variables
P1<-ggplot(tom, aes(x=a1, y=a2)) +
geom_point(shape=1) +      # Use hollow circles

geom_smooth(method=lm,    # Add linear regression line

se=FALSE)     # Don't add shaded confidence region

P2<- P1 + labs (x= "a1 value", y = "a2 value")## add x and y labels

P3 <- P2 + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="a1 value")                            #label and change font size

P4 <- P3 + scale_x_continuous("a1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5a <- P4 + annotate("text", x = 30, y = 10, label = lm_eqn(tom$a1, tom$a2, tom), color="black", size = 5, parse=TRUE)

P5a

Variabili a1 e a2

################ b1 and b2 couple of variables
P1b<-ggplot(tom, aes(x=b1, y=b2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2b<- P1b + labs (x= "b1 value", y = "b2 value")## add x and y labels

P3b <- P2b + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="b1 value")                            #label and change font size
P4b <- P3b + scale_x_continuous("b1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5b <- P4b + annotate("text", x = 40, y = 20, label = lm_eqn(tom$b1, tom$b2, tom), color="black", size = 5, parse=TRUE)

P5b

Variabili b1 e b2

################ c1 and c2 couple of variables
P1c<-ggplot(tom, aes(x=c1, y=c2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2c<- P1c + labs (x= "c1 value", y = "c2 value")## add x and y labels

P3c <- P2c + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="c1 value")                            #label and change font size

P4c <- P3c + scale_x_continuous("c1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5c <- P4c + annotate("text", x = 40, y = 15, label = lm_eqn(tom$c1, tom$c2, tom), color="black", size = 5, parse=TRUE)

P5c

Variabili c1 e c2

E adesso il gran finale. Come rappresentare i tre grafici in un’unica finestra?

Se il nostro Umesh, peraltro bravo programmatore, fosse stato meno pigro, avrebbe scoperto subito che CrossValidated non era il sito che faceva per lui, e avrebbe scoperto in fretta almeno due soluzioni al suo problema.

Le soluzioni che qui propongo sono tratte dal sito Cookbook for R, curato da Wiston Chang, una autentica miniera di gemme.

La prima soluzione utilizza il sistema di gestione delle griglie, e fa uso di due librerie aggiuntive a ggplot2 (librerie che dovrete scaricare se volete usarle).

library(grid)
library(gridExtra)

# semplice ed elegante (il titolo lo scegliete voi)
grid.arrange(P5a, P5b, P5c, ncol = 3, main = "Your main title")
# after saving, dev.off()
# se volete cambiare le dimensioni del titolo, seguite le istruzioni

### changing dimensions of the title
grid.arrange(P5a, P5b, P5c, ncol = 3,
main=textGrob("Your main title", gp=gpar(fontsize=20,font=3)))

Combinazione grafici

La seconda soluzione fa ricorso ad una funzione creata da Chang, e messa cortesemente a disposizione degli utenti dal creatore (il risultato è analogo a quello riportato sopra).


# Multiple plot function
# ggplot objects can be passed in …, or to plotlist (as a list of ggplot objects)
# – cols: Number of columns in layout
# – layout: A matrix specifying the layout. If present, ‘cols’ is ignored.
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
# from: http://www.cookbook-r.com/

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { require(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } # Una volta che avete caricata la funzione nell’area di lavoro, l’uso è semplicissimo # (Il numero di colonne (in questo caso 3, ma può variare a piacimento): multiplot(P5a, P5b, P5c, cols=3) [/code] That's all, folks! Antonello Preti Stay Tuned for our next episode!

Visualizzare relazioni curvilinee

Un modello statistico ha lo scopo di descrivere e prevedere uno o più fenomeni, semplificandoli, per renderli comprensibili all’uomo. Sotto questo punto di vista, i modelli statistici lineari costituiscono una soluzione ottimale per descrivere la realtà: sono semplici, facilmente interpretabili e adattabili a molte situazioni.

Nella sua formulazione più semplice, un modello lineare descrive la relazione tra due variabili attraverso una retta. Ci sono dei casi, però, in cui la relazione tra due variabili non può essere descritta attraverso una retta, ma è necessario interpolare i dati servendosi di una curva.

Descrivere una relazione di tipo curvilineo non implica necessariamente utilizzare un modello non lineare. La “linearità” di un modello, infatti, è determinata dalla forma matematica che questo assume e non dalla forma della relazione tra le variabili. Mi rendo conto che la cosa può apparire come un paradosso, ma è così: una relazione non lineare potrebbe benissimo essere descritta attraverso un modello lineare, purché questo modello sia lineare nei parametri.

In questo post vedremo un esempio di relazione non lineare e vedremo come visualizzare un modello di questo tipo utilizzando le funzioni basilari di R.

Un esempio: il compito di linea numerica


Una relazione non lineare è quella che si presenta quando un bambino piccolo, intorno ai 6/7 anni, è sottoposto al cosiddetto compito di linea numerica. Il compito consiste nel chiedere al bambino di posizionare una serie di numeri compresi entro un certo intervallo lungo un segmento. La versione più comune del compito richiede di posizionare lungo il segmento dei numeri che possono andare da 0 a 100.

Immaginiamo di presentare a un bambino trenta numeri compresi tra 0 e 100, che costituiscono gli stimoli, e di chiedergli di indicare in quale punto della linea dovrebbe essere posizionato ogni numero. Quindi, convertiamo ogni posizione indicata dal bambino nel corrispettivo numero che realmente si trova in quel punto della linea.

Registriamo nella variabile stimulus i numeri presentati e nella variable response le stime del bambino.

# Numeri presentati:
stimulus <- c(5,7,10,12,14,17,19,22,24,27,29,31,34,36,39,
            41,44,46,49,51,53,56,58,61,63,66,68,70,73,75)

# Risposte del bambino:
response <- c(38.22,42.67,40.63,45.81,44.09,44.15,45.7,52.6,
            49.51,48.13,49.68,51.83,49.4,53.48,50.06,48.76,
            53.22,51.54,53.04,56.08,53.24,50.9,52.12,53.11,
            55.95,54.81,53.96,54.62,50.92,55.61)
&#91;/code&#93;

Come ci ha insegnato Antonello Preti (<a href="">qui</a>), la primissima cosa che dovremmo fare è visualizzare la relazione tra le due variabili.

[code language="r"]
plot(stimulus, response)

Scatterplot tra stimolo e risposta nel compito di linea numerica

La relazione non sembra propriamente lineare, infatti i punti assumono un andamento leggermente curvilineo e tale curva sembra accentuata soprattutto per stimoli bassi piuttosto che elevati.

In ogni caso, se il bambino non ha risposto in maniera totalmente casuale, le sue risposte dovrebbero in qualche modo dipendere dagli stimoli presentati. Proviamo ad adattare un modello lineare:

model1 <- lm(response ~ stimulus)
summary(model1)
&#91;/code&#93;

Il modello presenta un indice R² pari a 0.6965. Non male. Vediamo allora come la retta di regressione interpola i punti:

&#91;code language="r"&#93;
plot(stimulus, response)
abline(model1, col="red")
&#91;/code&#93;

<img src="http://www.insular.it/wp-content/gallery/post-images/nlt-plot1.png" alt="Modello lineare risposta = a+b*stimolo" width="60%">

La retta non sembra interpolare benissimo i punti. Per carità, il modello non è malaccio, ma se diamo un'occhiata alla letteratura scientifica scopriremo che la relazione tra le due variabili non è lineare ma sembra essere <b>logaritmica</b>.

Quindi, a questo punto sorge un problema: come possiamo applicare un modello lineare se la relazione è non lineare (nel caso logaritmica)? Semplice: si trasforma <i>stimulus</i>, in modo da rendere lineare la sua relazione con <i>response</i>. L'assunzione è che la relazione tra il logaritmo naturale di <i>stimulus</i> e <i>response</i> sia lineare.

Adattare questo modello con R è molto semplice:

[code language="r"]
model2 <- lm(response ~ log(stimulus))
summary(model2)
&#91;/code&#93;

L'indice R² è decisamente superiore rispetto al precedente: 0.8275. La relazione potebbe davvero essere di tipo logaritmico. Come possiamo visualizzare questo modello? La funzione <b>abline</b> disegna rette, ma qui noi dobbiamo rappresentare una curva.

Per rappresentare un modello del genere si utilizza la funzione <b>lines</b>, in questo modo:

[code language="r"]
plot(stimulus, response)
lines(stimulus, predict(model2), col="blue")

Modello lineare risposta = a+b*log(stimolo)

La funzione predict estrapola i valori previsti dal modello. Alla funzione lines bisogna passare la variabile già disposta sull'asse orizzontale (stimulus) e i valori previsti dal modello, estrapolati proprio attraverso la funzione predict. Il risultato è la curva visualizzata in azzurro, che sembra interpolare i punti decisamente meglio rispetto alla retta del modello precedente.