Archivio dell'autore: Davide Massidda

Davide Massidda

Informazioni su Davide Massidda

Dottore di ricerca in psicologia sperimentale, mi occupo di consulenza statistica nel campo delle scienze psicologiche, sanitarie e sociali.

La variabile età

Svincolarsi da vecchi retaggi del passato, tanto consolidati da diventare prassi, non è semplice. Le abitudini ottimizzano il nostro modo di agire, rendendolo per molti aspetti più efficiente. Capita però che i tempi cambino e certe abitudini ci portino a sprecare risorse invece di risparmiale.

Qualche decennio fa, nelle nostre vite è entrato il computer. Non solo: negli ultimi anni le prestazioni di calcolo di queste macchine sono diventate tali da consentire elaborazioni quasi impensabili fino a soli dieci anni fa. Eppure, nonostante questo, per certi aspetti continuiamo a costruire i dataset in modo simile a come venivano costruiti carta e matita.

È per esempio il caso della variabile “età”.

Codificare la variabile età in un dataset è una questione di scelte. Comunemente, se parliamo di adulti o di bambini in età scolare, l’età viene registrata in anni compiuti. Se invece si parla di bambini molto piccoli, l’età viene più spesso registrata in mesi compiuti, se non addirittura in giorni se i soggetti sono neonati. Ci sono poi i ricercatori a cui piace aggregare, per cui registrano le età come fasce: “da 30 a 35”, “da 36 a 40”, eccetera eccetera. In questo caso, le fasce sono spesso definite a seconda del legame che l’età potrebbe avere con altre variabili incluse nel dataset, per cui il modo in cui questa viene codificata è contestuale alla raccolta dati e non indipendente da essa.

La codifica in fasce d’età è secondo me retaggio di un passato in cui i computer non c’erano o avevano delle prestazioni molto limitate rispetto a quello a cui siamo abituati oggi. Anni addietro era necessario economizzare e semplificare il più possibile i calcoli: aggregare faceva comodo.

Oggi, aggregare è anacronistico. Chiaro, un minimo di livello di aggregazione ci può essere e fa comunque comodo, ma spesso si esagera, perdendo della preziosissima informazione. Ricordiamoci infatti che da un’informazione micro se ne può ricavare una macro, mentre non è sempre vero il contrario. Se io vi dicessi che ho 33 anni, potrete facilmente dedurre che posso essere collocato nella fascia “30-35”. Se però vi dicessi che mi trovo nella fascia “30-35” e vi chiedessi quanti anni ho, sapreste rispondermi?

Aggregare le età in fase di codifica dati è secondo me controproducente per tre ragioni.

Primo, ad aggregare si fa sempre in tempo in un momento successivo e più consono (sempre che usiate un software abbastanza flessibile che vi consenta di farlo).

Secondo, in fase di analisi dei dati potreste rendervi conto che quel tipo di aggregazione non è adeguato al fenomeno che state studiando e potreste desiderare delle fasce d’età costruite in un altro modo.

Terzo, l’età è una variabile continua. Per comodità noi possiamo concepire il tempo come discreto, ma ovviamente non è tale. Classificando l’età, anche semplicemente esplicitandola in anni compiuti, la stiamo in un certo senso “discretizzando”. Certo, bisogna dire che “età” non è sinonimo di “tempo”, ma è un suo derivato; inoltre, l’età è un concetto che può assumere forme diverse a seconda del contesto (fasce, anni, mesi, giorni, ecc.), mentre il tempo no. A mio parere, però, codificando l’età dovremmo per quanto possibile evitare di spezzare il tempo: dovremmo invece provare a fermarlo.

Quando possibile, io chiedo sempre che nei dataset che mi vengono forniti non venga registrata l’età, o meglio, non venga registrata solo l’età. Soprattutto se le unità statistiche sono soggetti in età evolutiva, trovo molto utile che l’età sia accompagnata da altre due variabili: la data di nascita del soggetto e la data in cui sono state rilevate le misure. L’età può infatti essere tranquillamente dedotta a partire dallo scarto fra queste due variabili, con la differenza però che in questo caso sarà chi effettua le analisi a decidere su quale unità di misura esprimerla.

R dispone di un tipo di dato particolare che consente di esplicitare che una variabile è di tipo “data”. Questo tipo di dato non è in realtà semplicissimo da gestire, principalmente a causa dell’elevatissima varietà di modi che si possono usare per codificare le date. Per esempio, la data di oggi, 29 luglio 2016, può essere scritta come “29.07.16”, “29/7/2016”, “29 lug 2016”, eccetera, eccetera, eccetera… R utilizza come formato standard per la codifica delle date quello anglosassone, per cui richiede che la data sia scritta come “anno-mese-giorno”, tutto in formato numerico e facendo in modo che i giorni e i mesi siano sempre scritti con due cifre, mentre l’anno con quattro. Per intenderci, la data di oggi in questo formato sarebbe “2016-07-29”.

Se la data è già scritta in questo formato, possiamo utilizzare il comando as.Date() per esplicitare che la stringa è una data:

> x <- as.Date("2016-07-29") 
> x 
[1] "2016-07-29"
> str(x) 
 Date[1:1], format: "2016-07-29"

Se la data non è in questo formato… beh, lo deve diventare. Attraverso l’argomento format, il comando as.Date() offre la possibilità di esplicitare la struttura del formato di origine in modo che R sia in grado di comprendere come sono disposte le informazioni all’interno della stringa.

Vediamo subito un esempio:

> x <- as.Date("29/07/2016", format="%d/%m/%Y") 
> x 
[1] "2016-07-29" 
> str(x) 
 Date[1:1], format: "2016-07-29"

All’argomento format ho passato una stringa nella quale viene dichiarato che il primo numero rappresenta il giorno in formato numerico (%d), il secondo numero rappresenta il mese sempre in formato numerico (%m) e il terzo numero rappresenta l’anno in formato numerico a quattro cifre (%Y); inoltre, ogni valore è separato da uno slash e non da un trattino.

Come facevo a sapere come costruire la stringa da passare a format? Semplice: richiamando l’help ?strptime sarà possibile studiarsi tutti i dettagli.

Torniamo ora alla questione delle età. Immaginiamo di avere un vettore born che indica la data di nascita di ogni soggetto e un vettore test che indica la data in cui è stata rilevata una misura.

> born <- c("02/12/2012","12/11/2012","17/05/2013","04/07/2013","29/02/2013") 
> test <- c("21/10/2015","21/10/2015","21/10/2015","26/10/2015","26/10/2015")
&#91;/code&#93;

I dati contenuti in queste due variabili possono essere convertiti in date:

&#91;code language="R"&#93;
> born <- as.Date(born, format="%d/%m/%Y") 
> test <- as.Date(test, format="%d/%m/%Y") 
&#91;/code&#93;

Ecco il risultato:

&#91;code language="R"&#93;
> born 
[1] "2012-12-02" "2012-11-12" "2013-05-17" "2013-07-04" "2013-02-29" 
> test 
[1] "2015-10-21" "2015-10-21" "2015-10-21" "2015-10-26" "2015-10-26"

A questo punto, per ogni soggetto è possibile calcolare l’età facendo la sottrazione fra le due variabili:

> test - born 
Time differences in days 
[1] 1053 1073  887  844  970

R esprime la differenza fra due date in giorni. Per ogni soggetto, quindi, possiamo conoscere l’età espressa in giorni.

Chiaramente, non sempre l’età espressa in giorni ci può andare bene. Se desideriamo un’unità di misura più consona, come i mesi o gli anni compiuti, possiamo ricorrere alla funzione age_calc() del pacchetto kyotil.

> library(kyotil)
# Mesi:
> age_calc(born, test, units="months") 
[1] 34.61290 35.29032 29.12903 27.70968 31.90538 
# Anni:
> age_calc(born, test, units="years") 
[1] 2.884707 2.939352 2.430137 2.312329 2.657534

La stessa funzione age_calc() del pacchetto kyotil è contenuta anche nel pacchetto eeptools, ma probabilmente quella di eeptools è una versione più vecchia, perché sembra avere delle limitazioni che quella di kyotil non ha.

Infine, non dimentichiamoci che possiamo utilizzare il comando cut(), distributo con la versione base di R, per creare delle fasce d’età. Questa, però, è un’altra storia.

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.

Visualizzare la relazione fra due variabili likert

Lavorando nel campo della psicologia, spesso mi trovo ad avere a che fare con variabili che derivano da riposte a questionari fornite utilizzando scale di tipo Likert. Si tratta di variabili che possono assumere un numero molto limitato di modalità, comunemente da tre a cinque. Quando le categorie di riposta sono almeno cinque, spesso queste variabili vengono considerate come se fossero continue e di conseguenza vengono analizzate usando indicatori e modelli statistici pensati appunto per variabili continue.

Con questo post non voglio entrare nel merito della correttezza di queste scelte (condivisibili o meno, a seconda dei casi), ma concentrarmi sul modo di utilizzare al meglio gli strumenti statistici, principalmente i grafici.

Di seguito viene costruito il data frame likert che contiene due variabili: item1 e item2, che riportano le risposte (espresse su scala Likert a cinque punti) di sessanta ipotetiche persone a due ipotetiche domande di un ipotetico questionario.

likert <- data.frame(
    item1 = c(5,3,5,5,4,1,2,3,3,5,5,2,1,3,3,3,4,2,4,1,
              3,4,4,5,3,4,5,2,1,4,3,2,3,2,4,5,5,2,5,3,4,
              4,3,2,1,1,1,5,2,3,1,1,2,2,2,3,4,2,4,3),
    item2 = c(5,2,5,4,4,1,3,3,1,3,3,2,1,3,1,3,1,2,4,1,
              3,4,2,4,4,5,5,2,1,1,1,4,4,2,4,4,5,2,5,4,5,
              3,3,1,2,2,1,5,5,5,3,4,1,1,3,2,4,2,3,2)
)
&#91;/code&#93;

Per studiare la relazione tra due variabili di questo tipo, quello che comunemente viene fatto è calcolare l'indice di correlazione lineare. Possiamo usare l'indice di Pearson oppure l'indice di Spearman, basato sui ranghi.

Attraverso il comando <b>cor</b> di R possiamo calcolare entrambi gli indici; se non viene specificato nulla nell'argomento <i>method</i>, verrà calcolata la correlazione di Pearson:

[code language="R"]
> with(likert, cor(item1, item2))
[1] 0.5901713

Se invece vogliamo calcolare l'indice di Spearman, dobbiamo esplicitare questa richiesta nell'argomento method:

> with(likert, cor(item1, item2, method="spearman"))
[1] 0.5864863

Entrambi gli indici evidenziano un'ottima correlazione fra le due variabili (“ottima” almeno per gli standard in psicologia): si sfiora lo 0.6, che è un valore abbastanza elevato.

Quello che - ahimè - non fa quasi nessuno, è visualizzare la relazione tra le due variabili. I grafici hanno sempre tanto da raccontare e talvolta è proprio dalle visualizzazioni che emergono gli aspetti più interessanti. Proviamo quindi a creare uno scatterplot per visualizzare la relazione fra item1 e item2. Utilizziamo il comando plot, aumentando la dimensione dei punti sfruttando l'argomento cex:

with(likert, plot(item1, item2, cex=2))

Scatter-plot fra item1 e item2

Beh, che ve ne pare? Si tratta di una buona visualizzazione? Secondo me, no.

A vedere questo grafico a me sorgono molte perplessità. Abbiamo appena detto che la relazione lineare tra le due variabili è buona, ma dal grafico proprio non si direbbe: i punti sono sparpagliati un po' ovunque e le due variabili sembrano tutto fuorché correlate.

Il problema di questo grafico è che le variabili possono assumere pochi valori (da 1 a 5), per cui moltissime risposte si sovrappongono. Ognuno di quei pallini in realtà ha una densità, perché su ognuno di essi si sovrappongono le risposte di più individui. Osservando bene l'immagine, infatti, possiamo notare che ci sono pallini il cui contorno è più scuro di altri; ebbene, nelle coordinate più scure si concentrano le risposte di più persone.

Il numero di osservazioni presente in ogni coordinata è un dato fondamentale per comprendere la relazione tra due variabili che assumono un numero ridotto di modalità, ma nell'immagine qui sopra questa informazione non è ben rappresentata.

La “densità” può essere calcolata semplicemente contando il numero di risposte che occorrono per ognuno degli incroci dei valori delle due variabili, ovvero costruendo una tabella di frequenza a doppia entrata:

> tab <- with(likert, table(item1, item2))
> tab
     item2
item1 1 2 3 4 5
    1 5 2 1 1 0
    2 3 6 2 1 1
    3 3 3 5 3 1
    4 2 1 2 5 2
    5 0 0 2 3 6

Osservando la tabella qui sopra possiamo notare come le frequenze maggiori siano collocate sulla diagonale, fenomeno che supporta la presenza di una relazione lineare e che giustifica valori di correlazione così elevati. Ma come fare per considerare questa informazione nel grafico?

Adesso vi proporrò due alternative; entrambe richiedono che la tabella di frequenza che abbiamo appena costruito venga convertita in un oggetto di tipo data.frame:

> tab <- as.data.frame(tab)
> head(tab)
  item1 item2 Freq
1     1     1    5
2     2     1    3
3     3     1    3
4     4     1    2
5     5     1    0
6     1     2    2

Grafico a bolle (bubble chart)


Quello che manca al grafico costruito poco sopra è l'informazione sul numero di osservazioni in ogni coordinata. Il modo più semplice di considerare quella che è a tutti gli effetti una terza variabile è fare in modo che il diametro di ogni punto dipenda dalla frequenza. Verrà creato così un grafico nel quale saranno presenti punti più grandi e punti più piccoli: più grande sarà il punto, maggiore sarà la concentrazione di dati.

Per realizzare questo grafico possiamo usare il comando symbols, specificando nell'argomento inches un'unità di misura per calibrare la dimensione dei punti.

with(tab, symbols(item1, item2, Freq, inches=0.6))

Bubble-plot fra item1 e item2

Il risultato conferma che la densità è maggiore sulla diagonale e va diminuendo con l'allontanarsi da questa. I dati quindi non sono sparpagliati in maniera casuale come poteva sembrare nel primo grafico, ma la relazione fra item1 e item2 segue un andamento ben preciso (nello specifico lineare).

Grafico a mattonelle (tile plot)


Un altro tipo di grafico che ci viene in aiuto è il tile plot. In questo tipo di visualizzazione, i dati sono rappresentati attraverso delle mattonelle che vengono colorate a seconda del valore assunto da una terza variabile. Per costruire il tileplot dobbiamo installare e attivare il pacchetto ggplot2 (del quale abbiamo già parlato).

Allo strato di base costruito con il comando ggplot dobbiamo aggiungere un livello creato con geom_tile e un gradiente di colore con scale_fill_gradient; infine, possiamo specificare un tema (io ho scelto theme_bw):

library(ggplot2)
ggplot(data=tab, aes(x=item1, y=item2, fill=Freq)) +
    geom_tile(colour="white") +
    scale_fill_gradient(low="white", high="red3") +
    theme_bw()

Tile-plot fra item1 e item2

Il grafico a mattonelle è forse quello più accattivante, ma probabilmente anche il più complicato da costruire, visto che richiede l'uso di ggplot2.

E voi, quali soluzioni prediligete in questi casi?

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.