Archivio tag: Modelli a effetti misti

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!