Archivio tag: fattori within

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.