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.

Stima della sovrapposizione tra due distribuzioni empiriche con il pacchetto ‘overlapping’

In questo post vogliamo illustrare con alcuni semplici esempi come utilizzare la funzione overlap() del pacchetto overlapping per stimare il grado di sovrapposizione tra due distribuzioni empiriche.

Esempio 1


Supponiamo di avere raccolto dei dati in due gruppi di 100 soggetti ciascuno rispetto ad una generica variabile Y, espressa da punteggi teoricamente compresi tra 0 e 30, e di essere interessati a valutare se i due gruppi possano considerarsi campioni provenienti da popolazioni con la stessa media.

Di seguito si riporta il codice per la costruzione di un data frame che include la variabile Y, che contiene i punteggi, e la variabile group, che definisce a quale dei due gruppi fa riferimento ogni osservazione.

Y <- data.frame(
    Y = c(10,10,13,9,9,11,9,10,10,9,8,10,10,10,7,9,10,9,10,12,8,9,9,11,11,10,10,
        9,10,9,7,9,12,10,11,9,10,9,10,10,9,11,10,11,9,11,11,10,12,11,10,10,11,10,
        9,11,9,9,11,11,11,11,8,10,9,9,9,10,9,10,9,10,9,11,11,8,10,8,9,11,12,9,11,
        9,9,11,11,11,11,10,8,11,10,10,11,9,8,11,10,10,8,6,10,3,9,26,13,5,9,5,9,20,
        10,5,9,16,7,4,17,13,8,9,8,8,8,5,6,13,9,5,15,5,7,8,6,11,13,10,12,5,6,5,10,
        11,16,10,7,7,12,5,12,7,7,8,11,8,17,7,10,9,7,3,11,10,12,15,11,6,5,6,8,12,8,
        5,10,12,12,7,23,8,14,8,10,11,13,6,6,10,9,11,27,9,8,10,6,12,4,6,10,9),
    group = gl(2,100,labels=c("G1","G2"))
)

Nel primo campione i punteggi variano tra 7 e 13, con media 9.88 e dev. st. 1.12, nel secondo campione variano tra 3 e 27, con media 9.5 e dev. st. 4.33. Le medie sono effettivamente molto simili, ma notiamo che la variabilità associata ai due gruppi di punteggi è molto diversa, nel secondo gruppo è circa 4 volte quella del primo.

Rappresentando graficamente le distribuzioni otteniamo la seguente figura:

histogram(~Y|group,data=Y,n=50,type="count",xlab="punteggi",ylab="frequenze",border="#0080ff",xlim=c(-2,32))

Da tale figura abbiamo la conferma diretta della disomogeneità delle varianze nei punteggi dei due gruppi; inoltre, possiamo constatare che nel secondo gruppo la distribuzione dei punteggi è piuttosto asimmetrica con una evidente coda nella parte destra. In un caso di questo tipo, il confronto statistico tra le medie può essere improprio e poco informativo; ad esempio con un t-test, corretto per la non-omogeneità, otteniamo il seguente risultato: t(112.22)=0.85, p=0.397, da cui non possiamo trarre alcuna conclusione.

A questo punto proviamo a cambiare i termini del problema: piuttosto che valutare la somiglianza dei due gruppi sulla base delle sole medie (e deviazioni standard) vogliamo utilizzare tutta l’informazione disponibile nei dati. Si tratta in pratica di stimare quanta parte delle distribuzioni di punteggi nei due gruppi siano sovrapponibili; ci aspettiamo che 0% indichi l’assenza di sovrapposizione e 100% la perfetta sovrapposizione tra le due distribuzioni. Possiamo utilizzare il pacchetto overlapping nel seguente modo:

> library(overlapping)
> dataList <- list(G1=Y1,G2=Y2)
> dataList <- overlap(dataList)$OV

    G1-G2 
0.4732179 

Y1 e Y2 sono i vettori contenenti i punteggi osservati, creiamo una lista (dataList) e la passiamo alla funzione overlap(); il valore che otteniamo in output (0.4732179) è una stima della percentuale di sovrapposizione tra le distribuzioni dei punteggi Y1 e Y2. Se vogliamo avere una rappresentazione grafica delle densità stimate delle due distribuzioni con relativa area di sovrapposizione aggiungiamo l’opzione plot=TRUE:

> overlap(dataList,plot=TRUE)

La linea azzurra è la densità stimata della distribuzione di punteggi del primo gruppo, la linea magenta la densità stimata per il secondo gruppo, l’area in verde indica la porzione di sovrapposizione tra le due densità, pari a circa il 47.32%.

Con la stessa funzione si può ottenere anche un altro tipo di grafico:

> overlap(dataList,partial.plot=TRUE)

Questa seconda rappresentazione, più grezza, disegna con due linee continue (rossa e nera) le densità stimate e con una linea verde (più spessa) la parte di densità che rappresenta la sovrapposizione; le linee verticali tratteggiate rossa e nera sono poste in corrispondenza del valore massimo di densità di ciascuna curva, le linee grigie rappresentano i punti in cui la densità è la stessa per le due curve.

Esempio 2


La funzione overlap() consente di confrontare anche più coppie di distribuzioni, indipendentemente dal numero di osservazioni di ciascuna di esse. Proviamo a generare quattro distribuzioni di dati campionate dalle seguenti densità: normale (rnorm), χ2 (rchisq), esponenziale (rexp) e uniforme (runif):

> n <- 100
> set.seed(20160108)
> Y1 <- rnorm(n,10)
> Y2 <- rchisq(n,10)
> Y3 <- rexp(n,1/10)
> Y4 <- c(runif(n,0,20))
> dataList <- list(Y1,Y2,Y3,Y4)
> overlap(dataList,plot=TRUE)

Dato che ci sono 4 diverse distribuzioni, avremo 6 possibili aree di sovrapposizione, come si vede nella figura seguente:

Ciascun pannello è relativo al confronto tra le densità stimate su una coppia di distribuzioni (Yi vs Yj) rappresentate con le linee continue mentre l’area in verde indica la sovrapposizione, il cui valore stimato è riportato nel pannello stesso.

In sintesi, la funzione overlap() riceve in input una lista composta da k vettori di valori (numerici), stima la densità associata alla distribuzione di valori di ciascun vettore e quindi confronta a coppie tali densità stimando per ciascuna coppia il grado di sovrapposizione. In aggiunta permette di rappresentare graficamente le densità messe a confronto. Un primo obiettivo per le prossime versioni del pacchetto è di aggiungere opzioni che consentano all’utente una maggiore flessibilità per la rappresentazione grafica.

Visualizzare le strutture di correlazione con “qgraph”

Non so voi ma io ho una predilezione per le rappresentazioni grafiche, che si manifesta soprattutto in fase di ricognizione dei dati. Già Davide in uno dei suoi ultimi post aveva rimarcato l’importanza di visualizzare i dati, importanza che non posso che condividere! Nello specifico, il problema descritto da Davide consisteva nel rappresentare graficamente una relazione bivariata in cui le variabili considerate erano misurate su scala Likert. Tali relazioni apparivano difficilmente interpretabili con il classico grafico di dispersione (o scatterplot). La soluzione proposta – ovvero rappresentare la densità delle variabili sul diagramma – è sicuramente ottima, ma in alcuni casi può rivelarsi poco pratica, ad esempio nel caso in cui:

  1. si devono visualizzare più relazioni bivariate (ad esempio, tra i diversi item di un questionario);
  2. si vuole avere una prima idea delle relazioni che legano tra loro le variabili indagate (ad esempio, come accade quando si effettua un’indagine esplorativa dei dati).

Per l’analisi dei dati di un questionario che ho utilizzato per la mia tesi di laurea cercavo un metodo semplice e veloce per visualizzare le relazioni tra gli item. Ero già rassegnata a dover scrutare le matrici di correlazione in cerca di qualche relazione interessante, quando mi sono imbattuta per caso in un pacchetto di facile utilizzo e comprensione: qgraph. Questa libreria permette di creare una struttura reticolare in cui ogni nodo rappresenta una variabile e ogni connessione rappresenta una relazione. In particolare è possibile visualizzare:

  1. la forza della relazione attraverso lo spessore delle linee di connessione (linee più spesse indicano relazioni più forti);
  2. la direzione della relazione attraverso il colore delle linee di connessione (linee verdi indicano una correlazione positiva mentre linee rosse indicano una correlazione negativa);
  3. la significatività della correlazione in termini di p-value.

In realtà la libreria permette di fare tante altre cose interessanti, ma ritengo sia utile iniziare dalle basi ed eventualmente approfondire altre forme di visualizzazione in un diverso post. Per mostrare l’utilizzo della libreria ricorrerò ad un caso emblematico che sicuramente sarà il pane quotidiano di molti studenti – ma non solo – di psicologia: le relazioni tra gli item di un questionario.

La libreria qgraph è ampiamente utilizzata sia per l’analisi di variabili latenti che rappresentano costrutti di personalità (come i Big Five), sia per la visualizzazione delle relazioni dirette tra variabili osservate. Prima di iniziare occorre però fare una precisazione.

Bisogna ricordare che le risposte di un questionario su scala Likert sono variabili di tipo categoriale ordinale e, come tali, godono di proprietà diverse da quelle possedute dalle variabili quantitative (sareste capaci di fare la media tra “per niente d’accordo” e “assolutamente d’accordo”?). In alcuni casi la natura categoriale delle variabili ordinali viene forzata assegnando loro un numero e trattandole di fatto come numeriche. Questa forzatura ha dei pro (ad esempio, facilitare la codifica e l’interpretazione dei risultati) e dei contro. Rimando ad una lettura interessante che tratta alcuni dei “contro”.

Ora vediamo la libreria in azione!

Il file cor.item.rda è un’area di lavoro di R che contiene una matrice di correlazione (cor.item) tra i 24 item di un questionario proposto a 100 ragazzi di 14 anni. Per ogni item si chiedeva ai partecipanti alla ricerca di esprimere un giudizio su una scala Likert da 1 (“per nulla”) a 5 (“moltissimo”). Tenendo conto della natura ordinale degli item è stata calcolata la matrice di correlazione policorica attraverso la funzione auto_cor della stessa libreria qgraph (tale funzione permette di calcolare in modo statisticamente appropriato le correlazioni tra variabili quantitative, ordinali e anche tra variabili ordinali e quantitative, vedi qui).

View(cor.item)

La visione della matrice di correlazione tra i 24 item non è certo di immediata interpretazione.

Installiamo e carichiamo la libreria qgraph:

install.packages("qgraph")
library(qgraph)

Quindi, utilizziamo la funzione qgraph per visualizzare la relazione fra gli item, utilizzando l’argomento vsize per definire manualmente la dimensione dei nodi:

qgraph(cor.item, vsize=5)

Grafico 1

Il grafico mostra gli item del questionario in relazione di tipo “tutto con tutto”. Lo spessore delle linee di collegamento indica la forza della relazione: appare evidente come siano presenti correlazioni decisamente più forti di altre. Ho ritenuto opportuno tralasciare le correlazioni che avevano un indice di correlazione inferiore a .50 con l’argomento minimum=0.5.

qgraph(cor.item, vsize=5, minimum=0.5)

Grafico 2

Consideriamo la relazione fra tre item, ad esempio l’item8.b, l’item3.b e l’item7.b: il primo correla negativamente con il secondo (il coefficiente di correlazione è infatti -0.51), il quale invece correla positivamente con il terzo (il suo coefficiente di correlazione è infatti 0.59).

Una funzione utile è quella che permette di visualizzare la struttura di correlazione tra gli item del questionario raggruppandoli secondo un determinato criterio. Lo strumento che stavo studiando io presentava tre categorie di item; se notate, infatti, il nome dell’item è composto da un numero e da una lettera, che può essere a, b oppure c: è proprio questa lettera a identificare a quale raggruppamento l’item appartiene.

Per raggruppare gli item nel grafico è prima di tutto necessario creare una lista che definisca i gruppi cui fanno riferimento le variabili.

list.cor = list(
    a = c(1,4,7,10,13,16,19,22),
    b = c(2,5,8,11,14,17,20,23),
    c = c(3,6,9,12,15,18,21,24))

La lista list.cor contiene al suo interno tre vettori e ogni vettore riporta gli indici di colonna degli item di ciascun raggruppamento.

A questo punto creiamo un vettore di nome lab (labels) che contiene dei nuovi nomi da attribuire alle variabili; si tratta di un passaggio opzionale, ma che renderà il grafico più leggibile e le variabili immediatamente identificabili. Il seguente codice estrae i nomi delle variabili ed elimina la stringa “item” da ogni nome. In questo modo, la stringa “item1.a” diventerà “1.a”, la stringa “item1.b” diventerà “1.b”, e così via.

lab = gsub("item","",rownames(cor.item))

Un altro aspetto opzionale, ma sicuramente utilissimo, è l’uso dei colori. Nel vettore hue salviamo un vettore di colori che dovranno essere usati per marcare in maniera diversa i nodi che rappresentano le variabili afferenti ai tre gruppi:

hue = c("dodgerblue3","darkgoldenrod1","snow")

Infine, possiamo costruire il grafico:

qgraph(cor.item, vsize=5, minimum=0.5, groups=list.cor, color=hue, labels=lab, edge.label.cex=3)

Grafico 3

Molto più ordinato, no?

Al posto dei valori di correlazione è possibile visualizzare i p-value che indicano se gli indici possono essere considerati significativamente diversi da zero. Anche in questo caso la libreria si rivela molto comoda, perché permette di visualizzare la significatività dei parametri di regressione aggiungendo l’argomento graph = “sig2″ come segue:

qgraph(cor.item, minimum=0.5, vsize=5, groups=list.cor, graph="sig2", color=hue, labels=lab, alpha=c(0.01,0.05))

Grafico 4

Le linee blu indicano la significatività (al 5% e all’1%) delle correlazioni positive, mentre le linee arancioni indicano la significatività delle correlazioni negative (sempre al 5% e all’1%).

La mia indagine sulle relazioni tra gli item del questionario non finisce certo qui, questa libreria infatti mi sta aiutando a pianificare le indagini successive!

Spero di riuscire presto a mostrarvi altri grafici “in azione” con le funzioni sem e factanal (molto utili in psicologia).

Info e fonti:

  1. qgraph: Network Visualizations of Relationships in Psychometric Data
  2. State of the aRt personality research: A tutorial on network analysis of personality data in R
  3. Network Model Selection Using qgraph 1.3
  4. qgraph examples

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.