Archivio tag: Anova

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.

Rimodellare i dati

Quando dobbiamo organizzare una tabella per raccogliere dati provenienti da disegni sperimentali, la prima domanda che dovremmo porci è come organizzarli. All’università mi hanno sempre ripetuto che a ogni riga deve corrispondere un “soggetto”, mentre a ogni colonna una variabile. Questa regola, ampiamente condivisa e sulla base della quale sono programmati i software statistici, in alcuni casi può diventare ambigua. Nello specifico, questo accade quando un soggetto (unità statistica) viene misurato più volte.

Per chiarire il concetto, passiamo subito a un caso pratico. Si tratta di un esempio inventato e abbastanza semplicistico, che nella sua semplicità ci aiuterà a definire meglio il problema.

In un laboratorio di psicofisiologia del sonno si stanno testando gli effetti di una nuova psicoterapia di tipo cognitivo-comportamentale e di un trattamento di tipo farmacologico per la cura dell’insonnia. Entrambi i trattamenti hanno la durata di due mesi.
Sono state selezionate 10 persone insonni con problemi nella fase di addormentamento: 5 vengono trattate con psicoterapia (PT) e 5 con terapia farmacologica (FT).
Il giorno prima dell’inizio del trattamento, a 30 e a 60 giorni viene registrato il numero di ore dormite da ogni individuo.
Per verificare la loro ipotesi, i ricercatori vorrebbero eseguire un’analisi statistica per verificare se le persone rispondono diversamente ai due trattamenti, valutando se questi portano a delle differenze nel numero di ore di sonno.

Dovendo organizzare i dati di questo disegno sperimentale, o si inseriscono i soggetti per riga oppure le variabili per colonna, ma non è possibile soddisfare i due criteri simultaneamente. Vediamo perché.

Nella tabella qui sotto, sono mostrati i dati raccolti nel cosiddetto formato wide.

sogg tratt tempo.0 tempo.30 tempo.60
1 PT 5.60 3.30 8.90
2 PT 4.80 4.70 7.90
3 PT 3.90 5.30 8.50
4 PT 5.80 4.50 8.20
5 PT 3.50 3.70 6.90
6 FT 3.90 4.70 7.00
7 FT 4.10 4.90 5.60
8 FT 5.00 3.80 5.30
9 FT 4.00 4.70 5.30
10 FT 4.00 4.80 6.30

Come possiamo osservare, è vero che a ogni riga corrisponde un soggetto, ma non è vero che a ogni colonna corrisponde una variabile. Infatti, il fattore “tempo” è stato suddiviso in tre colonne, una per ogni livello: 0 giorni (inizio), 30 giorni e 60 giorni. “tempo.0”, “tempo.30” e “tempo.60” non sono delle variabili, ma sono tutti livelli di un’unica variabile. Software come SPSS o STATISTICA richiedono i dati disposti proprio in questo modo.

Un altro modo in cui i dati possono essere organizzati è il formato long, che possiamo vedere qui di seguito.

sogg tempo tratt ore
1 0 PT 5.60
1 30 PT 3.30
1 60 PT 8.90
2 0 PT 4.80
2 30 PT 4.70
2 60 PT 7.90
3 0 PT 3.90
3 30 PT 5.30
3 60 PT 8.50
4 0 PT 5.80
4 30 PT 4.50
4 60 PT 8.20
5 0 PT 3.50
5 30 PT 3.70
5 60 PT 6.90
6 0 FT 3.90
6 30 FT 4.70
6 60 FT 7.00
7 0 FT 4.10
7 30 FT 4.90
7 60 FT 5.60
8 0 FT 5.00
8 30 FT 3.80
8 60 FT 5.30
9 0 FT 4.00
9 30 FT 4.70
9 60 FT 5.30
10 0 FT 4.00
10 30 FT 4.80
10 60 FT 6.30

In questo formato, a ogni colonna corrisponde una variabile, ma a ogni riga corrisponde un’osservazione e non un soggetto. Infatti, dato che i soggetti sono misurati tre volte, per ognuno di essi sono presenti tre righe.

Ognuno di questi formati ha dei pro e dei contro. Il formato wide, essendo più compatto, consente di visualizzare con più facilità la tabella e permette di utilizzare ognuno dei tre tempi come se fosse una variabile diversa (in alcuni casi può tornare utile). Tuttavia, trovo che il formato long sia migliore da un punto di vista concettuale, perché consente di individuare immediatamente quali sono le variabili di uno studio, distinguendo chiaramente la dipendente “ore di sonno” dalla indipendente “tempo” (questo aspetto non è da sottovalutare nell’insegnamento della statistica).

Eccetto casi particolari, le funzioni R richiedono che i dati siano disposti in formato long. Però, capita spesso che i dati siano stati registrati in formato wide. Talvolta ci si può trovare nella situazione di dover passare i propri dati in formato long a un utente di SPSS che, comprensibilmente, li gradirebbe in formato wide. È quindi importante capire come poter passare da un formato all’altro utilizzando R.

Da wide a long


Costruiamo in R la tabella di dati in formato wide:

sonno.wide <- data.frame(
    sogg = factor(c(1,2,3,4,5,6,7,8,9,10)),
    tratt = c("PT","PT","PT","PT","PT","FT","FT","FT","FT","FT"), 
    tempo.0 = c(5.6,4.8,3.9,5.8,3.5,3.9,4.1,5,4,4),
    tempo.30 = c(3.3,4.7,5.3,4.5,3.7,4.7,4.9,3.8,4.7,4.8),
    tempo.60 = c(8.9,7.9,8.5,8.2,6.9,7,5.6,5.3,5.3,6.3)
)
&#91;/code&#93;

Per rimodellare i dati, trasformandoli in formato long, possiamo utilizzare il comando <b>reshape</b>, in questo modo:

[code]
sonno.long <- reshape(sonno.wide, varying=3:5, v.names="ore",
                      timevar="tempo", direction="long")
&#91;/code&#93;

La funzione reshape ha tanti argomenti che possono essere manipolati, io qui ho sfruttato giusto i principali. L'argomento <b>varying</b> richiede di specificare quali sono le colonne che dovranno essere impilate, cioè quelle colonne che contengono i dati (le ore di sonno): tempo.0, tempo.30 e tempo.60. Io qui ho le ho indicate specificando il loro indice di colonna, ma avrei anche potuto passare, invece degli indici, un vettore contenente i nomi delle variabili. Come visto in precedenza, queste tre colonne sono in realtà un'unica variabile, che nel formato wide era stata suddivisa. Ora, nel formato long, a partire da queste tre colonne se ne dovrà ottenere una nuova che contiene i dati. Questa nuova colonna assumerà il nome specificato nell'argomento <b>v.names</b>. Inoltre, al dataset dovrà essere aggiunta una colonna che indica in quale tempo ogni dato è stato registrato; il nome che questa colonna dovrà assumere è stato indicato nell'argomento <b>timevar</b>. Infine, bisogna esplicitare che il dataset dovrà essere convertito in formato long e questo viene dichiarato attraverso l'argomento <b>direction</b>.

Possiamo vedere di seguito un'anteprima del risultato:

[code]
> sonno.long[1:10,]
     sogg tratt tempo ore id
1.1     1    PT     1 5.6  1
2.1     2    PT     1 4.8  2
3.1     3    PT     1 3.9  3
4.1     4    PT     1 5.8  4
5.1     5    PT     1 3.5  5
6.1     6    FT     1 3.9  6
7.1     7    FT     1 4.1  7
8.1     8    FT     1 5.0  8
9.1     9    FT     1 4.0  9
10.1   10    FT     1 4.0 10

Da long a wide


Costruiamo in R la tabella di dati in formato long:

sonno.long <- data.frame(
    sogg = factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,10,10)),
    tempo = factor(c(0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60)),
    tratt = c("PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT",
        "FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT"),
    ore = c(5.6,3.3,8.9,4.8,4.7,7.9,3.9,5.3,8.5,5.8,4.5,8.2,3.5,3.7,6.9,
        3.9,4.7,7,4.1,4.9,5.6,5,3.8,5.3,4,4.7,5.3,4,4.8,6.3)
)
&#91;/code&#93;

Usiamo ora la funzione reshape:

&#91;code&#93;
sonno.wide <- reshape(sonno.long, v.names="ore", timevar="tempo",
                      idvar="sogg", direction="wide")
&#91;/code&#93;

Stavolta, con <b>v.names</b> ho specificato al comando reshape qual è la variabile da suddividere nelle colonne, e, attraverso <b>timevar</b>, ho indicato sulla base di quale variabile bisogna creare queste colonne; inoltre, ho dovuto indicare a <b>idvar</b> qual è la variabile che codifica per i soggetti. Infine, non dobbiamo dimenticarci di indicare a <b>direction</b> che dobbiamo trasformare i dati in formato wide.

Visualizziamo il risultato:

[code]
> 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

Suddivisione dei dati


Quando i dati sono in formato long, talvolta può essere comodo spezzare il dataset in più oggetti a seconda del livello di una variabile. Per esempio, potremmo avere l'esigenza di suddividere i dati delle persone trattate con psicoterapia da quelli delle persone trattate con terapia farmacologica, magari per creare dei grafici separati per ogni gruppo. Questa operazione può essere realizzata molto facilmente attraverso il comando split:

sonno.split <- split(sonno.long, f=sonno.long$tratt)
&#91;/code&#93;

Il comando richiede, oltre al dataset, la variabile da utilizzare come criterio per la suddivisione (che nell'esempio è la colonna <i>tratt</i>). Come risultato, avremo un oggetto di tipo lista con due slot, uno per ogni livello della variabile <i>tratt</i>.

Per estrarre i dati di ogni slot, possiamo usare il dollaro, per esempio:

[code]
> sonno.split$FT
   sogg tempo tratt ore
16    6     0    FT 3.9
17    6    30    FT 4.7
18    6    60    FT 7.0
19    7     0    FT 4.1
20    7    30    FT 4.9
21    7    60    FT 5.6
22    8     0    FT 5.0
23    8    30    FT 3.8
24    8    60    FT 5.3
25    9     0    FT 4.0
26    9    30    FT 4.7
27    9    60    FT 5.3
28   10     0    FT 4.0
29   10    30    FT 4.8
30   10    60    FT 6.3

Per chi volesse divertirsi, ci sono due pacchetti, reshape e reshape2, che permettono di andare ben oltre quanto presentato in questo post. Qui l'articolo di presentazione dell'autore dei pacchetti.

La dimensione temporale negli esperimenti

Quello che affronteremo adesso è un argomento abbastanza ampio, che, per essere trattato come si dovrebbe, meriterebbe la stesura non di un post, ma di un libro (almeno uno). Non sviscererò quindi tutti gli aspetti, ma proverò a lanciare degli input per eventuali approfondimenti.

Clessidra

Il tempo è una dimensione che non può essere trascurata in tutti quei disegni di ricerca che fanno uso di fattori within subjects, nei quali, in un modo o nell’altro, esso andrebbe tenuto in considerazione in fase di analisi dei dati (anche se, a onor del vero, non sempre questo viene fatto). Talvolta il tempo è considerato come una vera e propria variabile indipendente senza la quale lo studio non avrebbe motivo di esistere; altre volte, invece, è solo un impiccio che i ricercatori, se potessero, eliminerebbero come un qualsiasi altro fattore di disturbo.

Di seguito analizzeremo tre tipi di disegni di ricerca nei quali il tempo entra prepotentemente in gioco. Si tratta di tre disegni dalla struttura molto semplice: comunemente, una ricerca ha un assetto decisamente più articolato, ma per comprendere casi più complessi è bene partire da qualcosa di semplice (l’auspicio è che, alla fine, generalizzare sarà cosa semplice). Nel corso della trattazione farò riferimento prevalentemente alla ricerca sperimentale in campo psicologico, perché questo è il mio campo. Penso comunque che il lettore proveniente da altri albiti disciplinari avrà poca difficoltà a proiettare gli esempi nel proprio settore di interesse.

In questo post non entrerò nel merito di questioni legate all’elaborazione statistica dei dati, che saranno oggetto di prossimi interventi sul blog.

Caso #1

Immaginiamo di stare studiano l’evoluzione della condizione clinica di uno gruppo di pazienti in conseguenza a un trattamento terapeutico. I pazienti vengono valutati in quattro momenti distinti: prima dell’inizio del trattamento (t0), durante il trattamento (t1) e al termine del trattamento (t2); inoltre, per verificare il mantenimento di un eventuale miglioramento, i pazienti vengono misurati anche qualche tempo dopo la conclusione della terapia (t3).

Ho provato a schematizzare questo disegno nella figura sottostante, dove è rappresentata la linea del tempo con indicati i quattro momenti di misurazione. Mi si conceda l’assenza di un gruppo di controllo, che certamente dovrebbe essere contemplato in uno studio di questo tipo, ma che qui ho trascurato a favore della semplicità espositiva.

Schema caso 1

In uno studio come questo, la variabile “tempo” è un fattore fisso a quattro livelli, perché quattro sono le misurazioni. Inoltre, il tempo è anche un fattore within subjects, perché, per ogni livello, le misurazioni vengono effettuate sugli stessi pazienti.

In questo caso #1, il tempo è una variabile fondamentale e imprescindibile nella ricerca: si sta studiando l’evoluzione dei pazienti e il concetto di “evoluzione” presuppone il concetto di “tempo”.

Questo tipo di disegno spesso viene detto “a misure ripetute”. Come già avevo sottolineato in un precedente post, certamente su ogni paziente viene ripetuta la misurazione, ma in condizioni diverse: si suppone che, tra una misurazione e l’altra, i pazienti siano cambiati a seguito del trattamento. Le unità su cui vengono rilevati i dati sono di volta in volta le stesse (i pazienti), ma ogni misurazione è di per sé diversa, perché effettuata a seguito di un presunto cambiamento. Se ogni paziente fosse misurato più volte a t0, più volte a t1, ecc., quella sì che sarebbe una vera misura ripetuta!

Caso #2

Comunemente, in un esperimento di laboratorio si ha un fattore “condizione”, i cui livelli rappresentano le singole condizioni sperimentali. Quando questo fattore è di tipo within subjects, ogni partecipante alla ricerca affronta tutte le condizioni sperimentali, per cui viene testato per tutti i livelli del fattore, che nella figura qui sotto ho denominato A, B e C.

Schema caso 2

Questo disegno sperimentale è molto simile al precedente nella struttura, ma profondamente diverso nella filosofia. Infatti, di solito l’ordine di presentazione delle condizioni può essere alterato senza generare alcun problema. Un partecipante può essere testato prima nella condizione A, poi B e infine C, oppure prima in B, poi in C e infine in A, e così via. Addirittura, alcuni partecipanti possono affrontare le condizioni in un senso, altri in un altro senso, altri in un altro senso ancora, ecc. Il tempo non dovrebbe giocare alcun ruolo in in uno studio del genere, per cui l’ordine di presentazione delle condizioni può essere variato. Nello schema più sopra, infatti, non è rappresentata la linea del tempo.

Alterare l’ordine dei livelli del fattore non sarebbe possibile nel caso #1; un mondo nel quale è possibile effettuare una misurazione pre-trattamento anche dopo aver iniziato il trattamento, contemplerebbe l’esistenza di una macchina del tempo che il rasoio di Occam ci sconsiglia di contemplare.

Tra parentesi, il software statistico R è in grado di distinguere i fattori dai livelli ordinati come il fattore “tempo” del caso #1, dai fattori dai livelli non ordinati, come il fattore “condizione” di questo secondo caso. I primi vengono definiti con la funzione ordered, mentre i secondi con la funzione factor.

In un disegno di ricerca come questo del caso #2, il tempo è solo un impiccio. È infatti possibile che le condizioni si influenzino vicendevolmente, per cui, per esempio, il fatto di aver affrontato prima A, potrebbe alterare l’esito della successiva condizione B se non addirittura di C.

Si tratta di un problema molto sentito negli studi sperimentali in psicologia. Il caso più comune è quello dell’effetto stanchezza: nel corso dell’esperimento i partecipanti rischiano di affaticarsi. Questo affaticamento potrebbe portare a scoprire una differenza tra le condizioni da imputare non a una reale diversità tra esse, ma al semplice fatto che le prestazioni dei partecipanti sono calate a causa di un fattore (l’affaticamento), che nulla ha a che vedere con lo studio e che non è contemplato nel disegno della ricerca. L’affaticamento implica che se ho affrontato prima A, la mia performance sarà migliore in A, se ho affrontato prima B, la mia performance sarà migliore in B, e così via. Allo stesso modo, i partecipanti potrebbero impratichirsi nel corso della prova, con un effetto apprendimento che, di fatto, comporterebbe problemi analoghi. Strategie come il controbilanciamento e l’uso di sessioni di allenamento servono proprio a limitare questi effetti.

In un esperimento, effetti indesiderati legati alla pratica, all’apprendimento, alla fatica, sono raggruppati sotto l’etichetta di effetto ordine. Effetti invece legati all’influenza di una condizione sull’altra, dovuta sempre all’ordine di presentazione, sono chiamati effetto sequenza. Per esempio, in un esperimento relativo alla percezione di pesantezza, dopo aver sollevato un oggetto molto pesante (condizione A), potrei percepire molto leggero l’oggetto successivo (condizione B); se invece l’oggetto della condizione B fosse presentato di seguito a un oggetto molto leggero (condizione C), forse lo percepirei molto pesante.

Dunque, se mai fosse possibile, il ricercatore testerebbe ogni soggetto contemporaneamente su tutte e tre le condizioni, onde evitare i problemi descritti. Il tempo non è contemplato nello studio e costituisce solo un fattore di disturbo, il cui effetto deve essere controllato e, per quanto possibile, limitato.

Caso #3

Il terzo caso è un’estensione del secondo e forse si presenta anche più di frequente. Si ha sempre un fattore “condizione”, ma, a differenza di prima, ogni soggetto sperimentale viene testato più volte in una medesima condizione. Nello schema riportato qui sotto, ho immaginato un fattore condizione a tre livelli e quattro misurazioni ripetute all’interno di ogni condizione.

Schema caso 3

Le quattro misurazioni rappresentate qui sopra sono delle vere misure ripetute: ogni unità viene ripetutamente misurata nelle medesime condizioni. Bisogna considerare che, per forza di cose, tra una misura e l’altra passa un certo intervallo di tempo: la prova può essere ripetuta in intervalli anche molto ravvicinati, ma tra una rilevazione e l’altra sarà sempre trascorso un certo lasso di tempo.

La ”misurazione ripetuta“ è una pratica molto comune in psicologia (ma non solo). In psicologia si devono fare i conti con variabilità molto ampie, anche inter-individuali. Misurando un’abilità mentale, un’opinione, un fenomeno percettivo, difficilmente una stessa persona darà due volte la stessa identica risposta. Rispondendo a uno stimolo, talvolta saremo più rapidi e altre volte più lenti, anche se lo stimolo non è cambiato. Eseguendo uno stesso compito, talvolta saremo più accurati, altre volte meno. Per questo motivo, in psicologia si diffida della misurazione secca e si preferisce testare e ritestare una persona nelle medesime condizioni, in modo da poter misurare anche la variabilità delle risposte.

Ripetendo la misurazione in una stessa condizione, stiamo inserendo nel disegno un ulteriore fattore, che comunemente viene denominato trial (prova). Difatti, tecnicamente lo schema raffigurato qui sopra rappresenta un disegno sperimentale 3×4: tre condizioni per quattro misurazioni. Il trial può essere considerato come un fattore fisso (di tipo within subjects) ed essere inserito come variabile indipendente in un modello come quello di analisi della varianza, ma potrebbe anche essere considerato come un fattore casuale. Che cosa è un fattore casuale? Se non lo sapete, vi consiglio di andare a leggere questo articolo (che, ve ne sarete accorti, sto citando in continuazione dall’inizio del post).

Il trial ha tutte le caratteristiche di un fattore casuale. Vediamole.

1. Analizzare le differenze tra i livelli del fattore trial non è oggetto dello studio (perlomeno, di solito è così). Ciò che interessa è studiare le differenze tre le condizioni sperimentali: le differenze tra le varie misurazioni costituiscono “rumore“ del quale lo sperimentatore farebbe a meno.

2. Posso effettuare una, due, tre… cento, mille misurazioni, senza che cambi il senso dello studio: aumentando o diminuendo il numero di trial, la logica sottostante lo studio non cambia. Diversamente, aggiungendo o togliendo condizioni sperimentali lo studio cambierebbe eccome, perché la struttura dell’esperimento risulterebbe alterata: non posso rimuovere condizioni sperimentali senza intaccare il fondamento teorico che ha mosso la ricerca. Bisogna comunque andarci cauti: il numero di trial deve essere scelto accuratamente, nell’ottica sì di misurare ripetutamente, ma evitando di affaticare troppo i partecipanti.

Si presti attenzione a una questione importante, molto rilevante in psicologia ma anche in altre discipline. Il fatto di aver scelto di testare un soggetto oggi, a una certa ora, implica una certa risposta; se invece il soggetto fosse stato testato ieri, magari a un’ora diversa, molto probabilmente esso avrebbe fornito un’altra risposta. Ecco che la “misura 1” di oggi risulterà diversa dalla “misura 1” che avrei ottenuto ieri. Questo fenomeno è implicito nell’esistenza del tempo e non lo possiamo contrastare: possiamo solo adeguatamente tenerlo in considerazione in fase di analisi dei dati (ma questa è un’altra storia). Ricordiamoci che, in un certo senso, un livello del fattore trial è un qualcosa di irripetibile e che mai più si ripresenterà.

Possiamo vedere il caso #2 come un caso particolare del #3 in cui il fattore trial ha un solo livello.

A questo punto, mi piacerebbe raccontarvi come si gestiscono le “misure ripetute” in R… ma, per questo, datemi ancora un po’ di tempo!

Un’unica funzione per eseguire t-test in sequenza

L’analisi della varianza (ANOVA) consente di confrontare le medie di più campioni per verificare se questi possono essere considerati come estratti dalla medesima popolazione. Spesso, l’ANOVA mette a confronto molti gruppi e chi analizza i dati può avere l’esigenza di mettere a confronto i gruppi a due a due, eseguendo un’analisi post-hoc tramite una serie di test t di Student.

Il principale comando di R che esegue dei confronti multipli tramite test t è pairwise.t.test, ma l’output che questa funzione fornisce risulta piuttosto scarno e riporta solo i valori p; per questo motivo, ho elaborato una nuova funzione che fornisce un output più completo. La funzione all_t_tests permette di realizzare velocemente tutti i possibili t-test derivati dalla combinazione dei livelli di una data variabile indipendente. Di seguito è riportato il codice che definisce la funzione, con due esempi sul suo utilizzo.

all_t_tests = function (x, g, p.adjust.method = p.adjust.methods, paired = paired,
    alternative = c("two.sided", "less", "greater"), var.equal = NULL) {
    g <- factor(g)
    p.adjust.method <- match.arg(p.adjust.method)
    alternative <- match.arg(alternative)
    comb_g = combn(levels(g),2)
    t.val = NA
    dof = NA
    pval = NA
    adj.pval = NA
    comparison = NA
    alternat = NA
    sigp = NA
    sigadjp = NA
    pair = NA
    for (i in 1:dim(comb_g)&#91;2&#93;){
        v1 = x&#91;g == comb_g&#91;1,i&#93;&#93;
        v2 = x&#91;g == comb_g&#91;2,i&#93;&#93;
        comparison&#91;i&#93; = paste(comb_g&#91;1,i&#93;, "vs", comb_g&#91;2,i&#93;)
        t.val&#91;i&#93; = t.test(v1, v2, alternative = alternative, paired = paired,
        var.equal = var.equal)$statistic
        dof&#91;i&#93; = t.test(v1, v2, alternative = alternative, paired = paired,
        var.equal = var.equal)$parameter
        pval&#91;i&#93; = t.test(v1, v2, alternative = alternative, paired = paired,
        var.equal = var.equal)$p.value
        sigp&#91;i&#93; =
            if(pval&#91;i&#93; <= 0.001){"***"} else{
                if(pval&#91;i&#93; <= 0.01){"**"} else {
                    if(pval&#91;i&#93; <= 0.05){"*"} else{""}
                }
            }
        adj.pval&#91;i&#93; = p.adjust(pval&#91;i&#93;, method = p.adjust.method, dim(comb_g)&#91;2&#93;)
        sigadjp&#91;i&#93; =
            if(adj.pval&#91;i&#93; <= 0.001){"***"} else{
                if(adj.pval&#91;i&#93; <= 0.01){"**"} else {
                    if(adj.pval&#91;i&#93; <= 0.05){"*"} else{""}
                }
            }
        alternat&#91;i&#93; = alternative
        pair&#91;i&#93; = paired
    }
    ans = data.frame(comparison, t.val, dof, pval, sigp, adj.pval, sigadjp, alternat, pair)
    names(ans) = c("Comparison", "t", "df", "pvalue", "sig.",
    paste("adj.pvalue_", p.adjust.method, sep=""), "sig_adj.", "alternative", "paired")
    return(ans)
}
&#91;/code&#93;
<h3 style="margin-top: 2em">Esempio 1: t-test per campioni indipendenti</h3>
<p style="text-align: justify">Per fornire un esempio d’uso della funzione mi servirò del dataset PlantGrowth. Il dataset riporta la variabile dipendente <i>weight</i> che identifica il peso di tre gruppi di piante: due gruppi sono stati sottoposti a due diversi trattamenti (trt1 e trt2) mentre un gruppo di piante funge da controllo (ctrl).</p>
[code language="r"]
&gt; data(PlantGrowth)
&gt; str(PlantGrowth)
'data.frame': 30 obs. of 2 variables:
$ weight: num 4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
$ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

Al data frame originale, aggiungo la variabile subj che mi permette di identificare i diversi soggetti (le piante) nei tre gruppi:

PlantGrowth$subj = as.factor(1:30)

Utilizzando la libreria sciplot, visualizzo i dati. Nella figura sottostante è mostrato il peso medio delle piante nei tre gruppi; le barre d'errore rappresentano gli errori standard.

library(sciplot)
lineplot.CI(group,weight,data=PlantGrowth,type="p",xlab="Group",ylab="Weight")

plot1

Adesso utilizzo il paccheto ez per calcolare una ANOVA unviariata con group come un fattore tra i soggetti (3 livelli: ctrl, trt1, trt2).

library(ez)
between_anova = ezANOVA(
data = PlantGrowth
, dv = weight
, wid = subj
, between = .(group)
)

L'output dell’ANOVA indica che l'effetto del gruppo risulta significativo:

> print(between_anova)
$ANOVA
  Effect DFn DFd        F          p p<.05       ges
1  group   2  27 4.846088 0.01590996     * 0.2641483

$`Levene's Test for Homogeneity of Variance`
  DFn DFd       SSn     SSd        F         p p<.05
1   2  27 0.3495267 4.21611 1.119186 0.3412266
&#91;/code&#93;
<p style="text-align: justify">Utilizzo quindi la funzione <strong>all_t_tests</strong> per eseguire tutti i confronti multipli con correzione di Bonferroni:</p>
[code language="r"]
ttest_between = all_t_tests(
    x = PlantGrowth$weight, g = PlantGrowth$group,
    paired = FALSE, p.adjust.method = "bonferroni", var.equal = TRUE
)

Prima di tutto, bisogna specificare gli argomenti x e g, che rappresentano rispettivamente la variabile dipendente e la variabile che identifica i diversi gruppi/condizioni. L'argomento paired, impostato a TRUE, indica che devono essere eseuiti dei t-test per campioni appaiati (sarebbe FALSE per campioni indipendenti); l'argomento p.adjust.method specifica quale metodo di correzione del p-value deve essere usato (si veda ?p.adjust per un elenco dei metodi disponibili) e var.equal specifica se deve essere assunta (= TRUE) o meno (= FALSE) l'omogeneità delle varianze fra i campioni. Nel caso di campioni indipendenti (paired = FALSE), è necessario specificare l'argomento var.equal. Un ulteriore argomento che è possibile specificare è alternative, che per default è impostato come "two.sided" (si veda ?t.test per un elenco delle opzioni disponibili).

In output la funzione restituisce un comodo data frame:

> print(ttest_between)
    Comparison         t df      pvalue sig. adj.pvalue_bonferroni sig_adj. alternative paired
1 ctrl vs trt1  1.191260 18 0.249023166                 0.74706950            two.sided  FALSE
2 ctrl vs trt2 -2.134020 18 0.046851385    *            0.14055415            two.sided  FALSE
3 trt1 vs trt2 -3.010099 18 0.007518426   **            0.02255528        *   two.sided  FALSE
  • Comparison: confronto effettuato.
  • t: valore della statistica t di Student.
  • df: gradi di libertà della statistica t.
  • pvalue: p-value NON corretto.
  • sig.: * per p-value ≤ 0.05, ** per pvalue ≤ 0.01, *** per pvalue ≤ 0.001.
  • adj.pvalue_[method]: metodo di aggiustamento dei p-value.
  • sig_adj.: * per p-value ≤ 0.05, ** per pvalue ≤ 0.01, *** per pvalue ≤ 0.001.
  • alternative: c("two.sided", "less", "greater").
  • paired: TRUE o FALSE.

Esempio 2: t-test per campioni appaiati

Per realizzare dei t-test per campioni appaiti, transformo la variabile group del data frame PlantGrowth nella variabile time facendo finta che le la variabile weigth sia stata misurata sullo stesso campione di piante più volte nel tempo; inoltre, come nell’esempio precedente, aggiungo la variabile subj di modo che questa identifichi i medesimi soggetti per i 3 tempi. Nell’utilizzare la funzione all_t_tests imposteremo paired = TRUE ed useremo il metodo di correzione “holm”.

# Carico il dataset PlantGrowth:
data(PlantGrowth)
# Aggiungo la colonna subj:
PlantGrowth$subj = as.factor(as.character(rep(1:10,3)))
# Trasformo la variabile group nella variabile time e ne cambio l'ordine:
levels(PlantGrowth$group) = c("t1", "t0", "t2")
names(PlantGrowth)[2] = "time"
PlantGrowth$time = factor(PlantGrowth$time, levels = c("t0", "t1", "t2"))

Osserviamo la struttura del nuovo dataset:

> str(PlantGrowth)
'data.frame': 30 obs. of 3 variables:
$ weight: num 4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
$ time : Factor w/ 3 levels "t0","t1","t2": 2 2 2 2 2 2 2 2 2 2 ...
$ subj : Factor w/ 10 levels "1","10","2","3",..: 1 3 4 5 6 7 8 9 10 2 ...

Utilizzando sempre la libreria sciplot, visualizziamo i dati:

library(sciplot)
lineplot.CI(time,weight,data=PlantGrowth,type="p",xlab="Time",ylab="Weight")

plot2

Utiliziamo la funzione ezANOVA per fare una ANOVA univariata con time come unico fattore (3 livelli: t0, t1, t2) entro i soggetti:

library(ez)
within_anova = ezANOVA(
data = PlantGrowth
, dv = weight
, wid = subj
, within = .(time)
)

Output dell’ANOVA: l’effetto del gruppo risulta significativo.

&gt; print(within_anova)
$ANOVA
  Effect DFn DFd        F          p p<.05       ges
2   time   2  18 3.651706 0.04664882     * 0.2641483

$`Mauchly's Test for Sphericity`
  Effect         W        p p<.05
2   time 0.8700054 0.572912      

$`Sphericity Corrections`
  Effect     GGe     p&#91;GG&#93; p&#91;GG&#93;<.05      HFe      p&#91;HF&#93; p&#91;HF&#93;<.05
2   time 0.88496 0.0540233           1.085686 0.04664882         *
&#91;/code&#93;
<p style="text-align: justify">Usiamo adesso la funzione <b>all_t_tests</b> per eseguire tutti i confronti multipli con correzione di Holm:</p>
[code language="r"]
ttest_within = all_t_tests(
    x = PlantGrowth$weight, g = PlantGrowth$time,
    paired = TRUE, p.adjust.method = "holm"
)
> print(ttest_within)
  Comparison          t df     pvalue sig. adj.pvalue_holm sig_adj. alternative paired
1   t0 vs t1 -0.9938415  9 0.34626729           1.00000000            two.sided   TRUE
2   t0 vs t2 -2.8463514  9 0.01920314    *      0.05760942            two.sided   TRUE
3   t1 vs t2 -1.7720834  9 0.11014394           0.33043183            two.sided   TRUE

Rapide statistiche descrittive con tapply

Quando abbiamo dei nuovi dati da elaborare, la prima cosa da fare è prendere confidenza con loro. Prima di lanciarci in qualsiasi analisi di tipo inferenziale, infatti, è sempre una buona idea osservare un po’ i numeri, sia graficamente che attraverso dei sintetici indicatori statistici. Approcciarsi ai dati con un po’ di statistiche descrittive e visualizzazioni aiuta molto a decidere come muoversi quando ci sarà da adattare un qualsiasi modello.

Nella pratica quotidiana, una cosa che mi capita molto spesso di fare è analizzare dati di studi sperimentali, dove, tipicamente, si hanno una o più una variabili dipendenti di tipo quantitativo e una o più variabili indipendenti di tipo qualitativo (fattori). In questi casi, una delle prime cose che faccio è calcolare, per ogni livello dei fattori, un indice descrittivo sulla variabile dipendente: mediana, media, deviazione standard, ecc.

Non c’è nulla di nuovo in questo; tuttavia, quando mi trovo a mostrare come eseguire questa operazione in R, noto sempre molta difficoltà da parte di chi sta muovendo i primi passi con il linguaggio, ed ecco perché mi trovo a scrivere questo post. Impariamo insieme, una volta per tutte, a utilizzare la funzione tapply per calcolare un po’ di statistiche descrittive su dati di disegni sperimentali.

Prima di cominciare, vediamo di procurarci un dataset da analizzare. R contiene al suo interno dei dataset già caricati; uno di questi è CO2, dove sono registrati i dati di un esperimento sulla tolleranza al freddo di piante della specie echinochloa crus-galli, che fanno parte della famiglia delle graminacee. Carichiamo e osserviamo il dataset:

data(CO2)
str(CO2)
View(CO2)

Nel corso dell’esperimento è stato misurato il livello di assorbimento di CO2 di sei piante provenienti dal Quebec e di altre sei piante provenienti dal Mississippi; le misurazioni sono state effettuate a sette diversi livelli di concentrazione ambientale di CO2. Metà delle piante sono state Trattate con refrigeramento la notte prima dell’esperimento.

Abbiamo quindi una variabile dipendente (uptake) rappresentata dall’assorbimento di CO2 delle piante. Abbiamo poi tre fattori sperimentali: la provenienza delle piante (Type), il tipo di trattamento a cui queste sono state sottoposte (Treatment) e i diversi livelli di contrazione di CO2 in cui sono state effettuate le misurazioni (conc).

Quello che dobbiamo fare adesso è calcolare qualche statistica descrittiva, per studiare come il grado di assorbimento di CO2 varia fra le diverse condizioni sperimentali. Per fare questto, utilizzeremo la funzione tapply. La funzione tapply suddivide i valori di una variabile quantitativa per i livelli di una variabile qualitativa e, per ogni livello, applica una particolare funzione sui dati.

Immaginiamo per esempio di voler calcolare la media di assorbimento di CO2 suddividendo le piante a seconda del tipo di trattamento: quelle che hanno subito refrigerazione e quelle che non l’hanno subita. Alla funzione tapply dovremo passare, in ordine: 1) La variabile quantitativa, 2) la variabile qualitativa sulla base della quale suddividere i dati, 3) la funzione da applicare sui dati. Il codice da usare sarà quindi così strutturato:

tapply(CO2$uptake, CO2$Treatment, mean)

Come output, otterremo la concentrazione media di CO2 per le piante non trattate (nonchilled) e trattate (chilled):

nonchilled    chilled 
  30.64286   23.78333 

Nel caso in cui la funzione da applicare richieda la specificazione di ulteriori argomenti, questi possono esse aggiunti di seguito a quelli obbligatori. Per esempio, la funzione mean mette a disposizione un utilissimo argomento che si chiama na.rm, che, se specificato come TRUE, eliminerà dal computo i dati mancanti. Per default, na.rm è FALSE; se vogliamo modificarne il valore specificandolo come TRUE, dovrà essere aggiunto in coda agli altri:

tapply(CO2$uptake, CO2$Treatment, mean, na.rm=TRUE)

Nel nostro caso non erano presenti dati mancanti, per cui la specificazione dell’argomento na.rm non sortirà alcun effetto.

La funzione tapply consente di applicare una funzione suddividendo i dati non solo per i livelli di un fattore, ma anche per la combinazione dei livelli di più fattori. Per esempio, possiamo calcolare le medie di assorbimento di piante trattate e non trattate, suddividendo gli esemplari a seconda della loro provenienza:

tapply(CO2$uptake, list(CO2$Treatment, CO2$Type), mean)

Per svolgere questa operazione abbiamo dovuto inserire i due fattori all’interno di una lista. La funzione si occupa di incrociare i livelli delle variabili qualitative in lista e di calcolare la media. Questo sarà l’output:

             Quebec Mississippi
nonchilled 35.33333    25.95238
chilled    31.75238    15.81429

Per comodità, possiamo combinare la funzione tapply con la funzione with. La funzione with permette di esplicitare all’inizio del comando da quale data frame estrapolare le variabili, evitando di dover continuamente ripetere il nome del dataset prima di ogni variabile. Si usa in questo modo:

with(CO2, tapply(uptake, list(Treatment, Type), mean))

Utilizzare la funzione with equivale a dire a R: «Con il dataset indicato come primo argomento (CO2 nel nostro caso), esegui le operazioni descritte nel secondo argomento». La funzione with risulta comoda soprattutto quando dobbiamo ripetere molte volte il nome del dataset. Per esempio, se dobbiamo calcolare i valori di assorbimento medio per ogni incrocio dei tre fattori dell’esperimento, utilizzando la funzione with il codice sarà:

with(CO2, tapply(uptake, list(conc, Treatment, Type), mean))

Sia per le piante “Quebec” che per quelle “Mississipi”, potremo visualizzare l’assorbimento medio di CO2 a diversi gradi di concentrazione, per le piante non trattate e per quelle trattate.

, , Quebec

     nonchilled  chilled
95     15.26667 12.86667
175    30.03333 24.13333
250    37.40000 34.46667
350    40.36667 35.80000
500    39.60000 36.66667
675    41.50000 37.50000
1000   43.16667 40.83333

, , Mississippi

     nonchilled  chilled
95     11.30000  9.60000
175    20.20000 14.76667
250    27.53333 16.10000
350    29.90000 16.60000
500    30.60000 16.63333
675    30.53333 18.26667
1000   31.60000 18.73333

Chiaramente, la funzione tapply può essere utilizzata non solo per calcolare le medie, ma qualsiasi altra statistica descrittiva: deviazione standard, mediana, quantili, ecc. Con una breve stringa di codice possiamo avere un indice descrittivo per ogni cella del disegno fattoriale. Comodo, no?!