martedì 6 aprile 2010

Repeated measures ANOVA: confronto tra gruppi di misure ripetute

Supponiamo di voler analizzare una certa variabile in n individui, e che questa variabile venga misurata più di 2 volte per ogni soggetto (ad esempio la pressione arteriosa prima del trattamento, durante il trattamento e al termine del trattamento farmacologico). E' intuitivo che i valori misurati siano tra loro appaiati, essendo stati misurati sugli stessi soggetti. Se avessimo effettuato solo 2 misurazione, avremmo potuto effettuare il confronto tra i due gruppi con un paired t-test. Avendo qui più di due gruppi, dobbiamo ricorrere invece alla ANOVA per misure ripetute.

Il modello più semplice che ora vedremo per primo, considera la presenza di un solo fattore, entro cui ciascun soggetto viene misurato più volte. Si parla di modelli within, per distinguerli dai modelli between a 1 solo fattore (che si analizzano con la classica ANOVA per campioni indipendenti).
Dire one factor within subject ANOVA (o within subject one-way ANOVA) significa proprio eseguire una ANOVA a 1 fattore di classificazione, con misure ripetute all'interno (within) dello stesso fattore.

La formula generale del modello between è la seguente:

In altri termini, ciascuna misurazione può essere spiegata considerando la media generale (Grand Mean, in inglese), l'effetto del trattamento e l'errore individuale (che dipende dal trattamento e dal soggetto).

Mentre la formula generale del modello within one-way é:

In questo modello è stato aggiunto il parametro che tiene conto delle caratteristiche intrinseche del soggetto e cioè della variabilità individuale.

In termini di varianza, nel modello a campioni indipendenti abbiamo:

Ossia la varianza totale è pari alla somma della varianza che dipende dal trattamento, più la varianza d'errore (o residua), che non viene spiegata dal modello.

Nel modello within invece abbiamo:

Quindi la varianza totale dipende dai parametri di prima, più la varianza dei soggetti (ossia la variabilità individuale che dicevamo prima).
La somma è la varianza within-subjects.

Considerando i gradi di libertà di ciascuna varianza (per i quali c'è una dimostrazione, che però non riporto qui), possiamo scrivere la tabella ANOVA per il modello a misure ripetute:





Cominciamo adesso l'esercizio, e cerchiamo dapprima di svolgerlo by-hand; in seguito passeremo alle funzioni di R. Supponiamo di avere la seguente tabella


Subject time1 time2 time3
1 45 20 30
2 56 18 29
3 59 10 24
4 49 15 25


Un certo parametro è stato misurato in 4 soggetti, e ciascun soggetto è stato sottoposto a 3 misurazioni, in 3 tempi successivi. Ci chiediamo se la variabile tempo influisce sul parametro in considerazione.

Cominciamo a calcolare i vari parametri che ci servono:













Abbiamo ora completato i calcoli (la SS residuals non l'ho calcolata manualmente, perchè i calcoli sono troppo lunghi e noiosi). I valori F-calcolati vanno confrontati con gli F-tabulati ai rispettivi gradi di libertà. Risulta comunque immediato che
- F_soggetti è minore di F_calcolato: quindi non c'è differenza significativa tra i soggetti
- F_fattore tempo è maggiore di F_calcolato: quindi il fattore tempo influenza il parametro che abbiamo calcolato (in pratica c'è differenza tra le colonne della tabella esaminata).

Cerchiamo adesso di ottenere la stessa tabella ANOVA in R.
Innanzitutto scriviamo i dati, sottoforma di dataframe:


subj <- rep(1:4, each=3)
time <- rep(c("time1", "time2", "time3"), 4)
weights <- c(45, 20, 30, 56, 18, 29, 59, 10, 24, 49, 15, 25)
mydata <- data.frame(factor(subj), factor(time), weights)
names(mydata) <- c("subj", "time", "weights")

> mydata
subj time weights
1 1 time1 45
2 1 time2 20
3 1 time3 30
4 2 time1 56
5 2 time2 18
6 2 time3 29
7 3 time1 59
8 3 time2 10
9 3 time3 24
10 4 time1 49
11 4 time2 15
12 4 time3 25


Adesso utilizziamo la funzion aov() (aov sta per Analysis Of Variance), facendo attenzione alla formula da specificare:


myANOVA <- aov(weights ~ time + Error(subj/time), mydata)


Possiamo "tradurre" la formula specificata, in questo modo:
Effettuare il confronto dei parametri weights ordinati nelle colonne relative ai vari time, considerando che le misurazioni dei vari soggetti sono ripetute nella colonna time.

L'output è il seguente:


summary(myANOVA)

Error: subj
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 34.667 11.556

Error: subj:time
Df Sum Sq Mean Sq F value Pr(>F)
time 2 2795.17 1397.58 49.086 0.0001911 ***
Residuals 6 170.83 28.47
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Abbiamo ottenuto il calcolo si SS, MS e del valore F calcolato sulle colonne time. A meno di approssimazioni, i calcoli risultano tutti corretti, e i 3 asterischi accanto al p-value indicano che la differenza tra i gruppi è altamente significativa.

Una rappresentazione grafica della situazione, può essere ottenuta con la funzione interaction.plot(). Osserviamo i grafici che otteniamo:


par(mfrow=c(1,2))
interaction.plot(time, subj, weights)
interaction.plot(subj, time, weights)




Come possiamo osservare dai grafici, le differenze tra soggetti sono irrilevanti (come abbiamo ottenuto dalla nostra analisi numerica), mentre le differenze within-subjects sono significative (le linee sono separate e non parallele tra loro).




Ecco ora un pò di dati per esercitarvi:


levels
subj high low moderate
1 9 12 8
2 3 9 6
3 11 13 18
4 8 6 8
5 3 19 12
6 1 7 5
7 5 4 5
8 7 6 7
9 9 14 18
10 13 19 10


Un certo parametro è stato misurato su 10 soggetti, in 3 condizioni diffenti. Viene chiesto di valutare se le 3 condizioni influiscono sul parametro considerato.

RISPOSTA: la variabile considerata influisce sulle misurazione, essendo p-value = 0.0436 *.




Introduco brevemente il capitolo di disegni più complessi, a 2 fattori. Si parla di two-way within subject ANOVA. I possibili disegni sperimentali che possono essere fatti sono:
- two factor ANOVA with repeated measures on one factor (= two way mixed model ANOVA / between-within design)
- two factor ANOVA with repeated measures on both factor (=two-way within subject ANOVA)

Se outcome è il valore del parametro che viene misurato, mentre factor1 e factor2 sono i due fattori di classificazione, tra i quali si vuole stimare se vi è differenza, le formule da applicare in R sono rispettivamente:

per il mixed model:
aov(outcome ~ factor1*factor2 + Error(subject/factor1))
in questo caso abbiamo factor1=within, factor2=between
o in alternativa
aov(outcome ~ factor1*factor2 + Error(subject/factor2))
e qui abbiamo factor1=between, factor2=within

mentre per la ripetizione in entrambi i fattori abbiamo:
aov(outcome ~ factor1*factor2 + Error(subject/(factor1*factor2)))

Possiamo ottenere solo i valori F usando la funzione lme della library(nlme):
anova(lme(outcome ~ factor1*factor2, random= ~1|subject/factor1)) #oppure subject/factor2
e per il mixed model:
anova(lme(outcome ~ factor1*factor2, random= ~1|subject))