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.

Visualizzare la relazione fra due variabili likert

Lavorando nel campo della psicologia, spesso mi trovo ad avere a che fare con variabili che derivano da riposte a questionari fornite utilizzando scale di tipo Likert. Si tratta di variabili che possono assumere un numero molto limitato di modalità, comunemente da tre a cinque. Quando le categorie di riposta sono almeno cinque, spesso queste variabili vengono considerate come se fossero continue e di conseguenza vengono analizzate usando indicatori e modelli statistici pensati appunto per variabili continue.

Con questo post non voglio entrare nel merito della correttezza di queste scelte (condivisibili o meno, a seconda dei casi), ma concentrarmi sul modo di utilizzare al meglio gli strumenti statistici, principalmente i grafici.

Di seguito viene costruito il data frame likert che contiene due variabili: item1 e item2, che riportano le risposte (espresse su scala Likert a cinque punti) di sessanta ipotetiche persone a due ipotetiche domande di un ipotetico questionario.

likert <- data.frame(
    item1 = c(5,3,5,5,4,1,2,3,3,5,5,2,1,3,3,3,4,2,4,1,
              3,4,4,5,3,4,5,2,1,4,3,2,3,2,4,5,5,2,5,3,4,
              4,3,2,1,1,1,5,2,3,1,1,2,2,2,3,4,2,4,3),
    item2 = c(5,2,5,4,4,1,3,3,1,3,3,2,1,3,1,3,1,2,4,1,
              3,4,2,4,4,5,5,2,1,1,1,4,4,2,4,4,5,2,5,4,5,
              3,3,1,2,2,1,5,5,5,3,4,1,1,3,2,4,2,3,2)
)
&#91;/code&#93;

Per studiare la relazione tra due variabili di questo tipo, quello che comunemente viene fatto è calcolare l'indice di correlazione lineare. Possiamo usare l'indice di Pearson oppure l'indice di Spearman, basato sui ranghi.

Attraverso il comando <b>cor</b> di R possiamo calcolare entrambi gli indici; se non viene specificato nulla nell'argomento <i>method</i>, verrà calcolata la correlazione di Pearson:

[code language="R"]
> with(likert, cor(item1, item2))
[1] 0.5901713

Se invece vogliamo calcolare l'indice di Spearman, dobbiamo esplicitare questa richiesta nell'argomento method:

> with(likert, cor(item1, item2, method="spearman"))
[1] 0.5864863

Entrambi gli indici evidenziano un'ottima correlazione fra le due variabili (“ottima” almeno per gli standard in psicologia): si sfiora lo 0.6, che è un valore abbastanza elevato.

Quello che - ahimè - non fa quasi nessuno, è visualizzare la relazione tra le due variabili. I grafici hanno sempre tanto da raccontare e talvolta è proprio dalle visualizzazioni che emergono gli aspetti più interessanti. Proviamo quindi a creare uno scatterplot per visualizzare la relazione fra item1 e item2. Utilizziamo il comando plot, aumentando la dimensione dei punti sfruttando l'argomento cex:

with(likert, plot(item1, item2, cex=2))

Scatter-plot fra item1 e item2

Beh, che ve ne pare? Si tratta di una buona visualizzazione? Secondo me, no.

A vedere questo grafico a me sorgono molte perplessità. Abbiamo appena detto che la relazione lineare tra le due variabili è buona, ma dal grafico proprio non si direbbe: i punti sono sparpagliati un po' ovunque e le due variabili sembrano tutto fuorché correlate.

Il problema di questo grafico è che le variabili possono assumere pochi valori (da 1 a 5), per cui moltissime risposte si sovrappongono. Ognuno di quei pallini in realtà ha una densità, perché su ognuno di essi si sovrappongono le risposte di più individui. Osservando bene l'immagine, infatti, possiamo notare che ci sono pallini il cui contorno è più scuro di altri; ebbene, nelle coordinate più scure si concentrano le risposte di più persone.

Il numero di osservazioni presente in ogni coordinata è un dato fondamentale per comprendere la relazione tra due variabili che assumono un numero ridotto di modalità, ma nell'immagine qui sopra questa informazione non è ben rappresentata.

La “densità” può essere calcolata semplicemente contando il numero di risposte che occorrono per ognuno degli incroci dei valori delle due variabili, ovvero costruendo una tabella di frequenza a doppia entrata:

> tab <- with(likert, table(item1, item2))
> tab
     item2
item1 1 2 3 4 5
    1 5 2 1 1 0
    2 3 6 2 1 1
    3 3 3 5 3 1
    4 2 1 2 5 2
    5 0 0 2 3 6

Osservando la tabella qui sopra possiamo notare come le frequenze maggiori siano collocate sulla diagonale, fenomeno che supporta la presenza di una relazione lineare e che giustifica valori di correlazione così elevati. Ma come fare per considerare questa informazione nel grafico?

Adesso vi proporrò due alternative; entrambe richiedono che la tabella di frequenza che abbiamo appena costruito venga convertita in un oggetto di tipo data.frame:

> tab <- as.data.frame(tab)
> head(tab)
  item1 item2 Freq
1     1     1    5
2     2     1    3
3     3     1    3
4     4     1    2
5     5     1    0
6     1     2    2

Grafico a bolle (bubble chart)


Quello che manca al grafico costruito poco sopra è l'informazione sul numero di osservazioni in ogni coordinata. Il modo più semplice di considerare quella che è a tutti gli effetti una terza variabile è fare in modo che il diametro di ogni punto dipenda dalla frequenza. Verrà creato così un grafico nel quale saranno presenti punti più grandi e punti più piccoli: più grande sarà il punto, maggiore sarà la concentrazione di dati.

Per realizzare questo grafico possiamo usare il comando symbols, specificando nell'argomento inches un'unità di misura per calibrare la dimensione dei punti.

with(tab, symbols(item1, item2, Freq, inches=0.6))

Bubble-plot fra item1 e item2

Il risultato conferma che la densità è maggiore sulla diagonale e va diminuendo con l'allontanarsi da questa. I dati quindi non sono sparpagliati in maniera casuale come poteva sembrare nel primo grafico, ma la relazione fra item1 e item2 segue un andamento ben preciso (nello specifico lineare).

Grafico a mattonelle (tile plot)


Un altro tipo di grafico che ci viene in aiuto è il tile plot. In questo tipo di visualizzazione, i dati sono rappresentati attraverso delle mattonelle che vengono colorate a seconda del valore assunto da una terza variabile. Per costruire il tileplot dobbiamo installare e attivare il pacchetto ggplot2 (del quale abbiamo già parlato).

Allo strato di base costruito con il comando ggplot dobbiamo aggiungere un livello creato con geom_tile e un gradiente di colore con scale_fill_gradient; infine, possiamo specificare un tema (io ho scelto theme_bw):

library(ggplot2)
ggplot(data=tab, aes(x=item1, y=item2, fill=Freq)) +
    geom_tile(colour="white") +
    scale_fill_gradient(low="white", high="red3") +
    theme_bw()

Tile-plot fra item1 e item2

Il grafico a mattonelle è forse quello più accattivante, ma probabilmente anche il più complicato da costruire, visto che richiede l'uso di ggplot2.

E voi, quali soluzioni prediligete in questi casi?

L’analisi di covarianza

Una delle cose che spesso mi trovo a spiegare, è come si interpretano i parametri di un modello lineare quando vengono messi insieme predittori di natura quantitativa e qualitativa, ovvero variabili di tipo numerico e di tipo fattore. Se le variabili non sono in interazione fra loro, questo tipo di modello viene chiamato “analisi di covarianza”; se invece le due variabili sono fra loro in interazione, al modello non viene dato alcun nome specifico (si tratta di una terminologia secondo me un po’ fuorviante, che personalmente non condivido).

Mi è capitato diverse volte di veder usare questo modello senza una reale consapevolezza del significato dell’analisi. Spesso un ricercatore si limita a valutare la significatività statistica di ogni singola variabile indipendente, senza pensare che i parametri stimati hanno una vera e propria interpretazione geometrica. La lettura dei parametri di un modello lineare fornisce molte informazioni, che, se sapute leggere e rappresentare adeguatamente, avrebbero tanto da raccontare. Non tutti infatti sanno che, quando nel modello sono presenti due predittori, uno numerico e uno fattore, il modello può ancora essere rappresentato bidimensionalmente, al pari di una regressione semplice.

Se vogliamo capire, è bene partire da un esempio pratico. Per il mio esempio mi sono ispirato a questo articolo, dove viene studiato come la relazione tra la prestazione a un test di matematica e l’ansia provata per il test varia fra maschi e femmine.

Presso questo link è disponibile un dataset che contiene dei dati simulati che riproducono, almeno nella logica, i dati raccolti dagli autori dell’articolo. I numeri sono stati inventati di sana pianta da me, per cui non prendete per oro colato i risultati, che sono un po’ esagerati: quello che ci interessa è il metodo.

Iniziamo scaricando i dati in locale e importandoli in R:

> math <- read.csv2("mathematics.csv")
> str(math)
'data.frame':	140 obs. of  4 variables:
 $ math.perf: int  9 4 2 8 2 4 13 10 6 5 ...
 $ math.anx : int  17 10 12 12 14 6 11 5 10 14 ...
 $ test.anx : int  19 18 24 9 26 17 7 10 15 15 ...
 $ gender   : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...

Il dataset contiene quattro variabili: la prestazione al test in matematica (math.perf), una misura dell’ansia provata verso i test di valutazione delle abilità matematiche (math.anx), una misura dell’ansia provata verso i test di valutazione di abilità in generale (test.anx) e il sesso del partecipante alla ricerca (gender).

Ci interessa studiare la relazione tra la prestazione al test di matematica (math.perf) e l’ansia provata verso le specifiche situazioni in cui le abilità matematiche vengono valutate attraverso un test (math.anx).

Attenzione: nel seguito della trattazione mi appoggerò a dei grafici per illustrare i risultati, ma, per evitare di appesantire troppo il discorso, non mostrerò il codice R utilizzato. Magari la costruzione di questi grafici sarà oggetto di altri post; per ora, vi basti sapere sono stati realizzati sfruttando il pacchetto ggplot2.

Analisi preliminari


La relazione tra math.perf e math.anx sembra esistere veramente. Disponendo math.anx sull’asse orizzontale di uno scatterplot e math.perf sull’asse verticale, notiamo infatti che all’aumentare dell’ansia diminuisce la prestazione, in maniera approssimativamente lineare.

Scatterplot modello additivo

Proviamo ad adattare un modello lineare per verificare se è possibile predire la prestazione a partire dall’ansia:

> mod1 <- lm(math.perf ~ math.anx, data=math)
> summary(mod1)

Call:
lm(formula = math.perf ~ math.anx, data = math)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.1997  -6.2352  -0.8446   5.3344  23.6583 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.57838    2.06327   6.096 1.03e-08 ***
math.anx    -0.04734    0.13835  -0.342    0.733    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.283 on 138 degrees of freedom
Multiple R-squared:  0.0008476,	Adjusted R-squared:  -0.006393 
F-statistic: 0.1171 on 1 and 138 DF,  p-value: 0.7328

La tabella dei coefficienti mostra che il parametro associato a math.anx, stimato pari a -0.047, non è significativamente diverso da zero (p = 0.773). Quindi, l’ipotesi non sembra essere valida.

Abbiamo però un’altra variabile nel dataset, gender. Proviamo a distinguere nel grafico i dati dei maschi da quelli delle femmine, colorandoli in maniera differenziata.

Scatterplot modello additivo

Guarda un po’ che coincidenza: sembra che i dati delle femmine si dispongano tutti “sopra” quelli dei maschi! In pratica, sembrerebbe che la relazione tra math.anx e math.perf esista sia nei maschi che nelle femmine, ma le femmine abbiamo delle prestazioni in matematica superiori, per cui avrebbero bisogno di una retta di regressione diversa da quella dei maschi. È infatti come se con il nostro modello dovessimo costruire non una retta di regressione, ma due: una per i maschi, con un’intercetta più bassa, e una per le femmine, con un’intercetta più elevata.

Il modello additivo


Abbiamo appena visto che, per adattare un buon modello, dovremmo costruire una retta di regressione per i maschi diversa da quella delle femmine. La retta delle femmine dovrebbe infatti essere “più alta”, ovvero dovrebbe avere un’intercetta maggiore.

Quindi, che si fa? Si divide in due il dataset e si fanno due analisi separate? Ovviamente no: si aggiunge nel modello creato precedentemente anche la variabile gender.

> mod2 <- lm(math.perf ~ math.anx + gender, data=math)
> summary(mod2)

Call:
lm(formula = math.perf ~ math.anx + gender, data = math)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1216 -3.2256  0.2374  2.8479 10.6495 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  32.22762    1.45988   22.08   <2e-16 ***
math.anx     -0.86471    0.08249  -10.48   <2e-16 ***
genderM     -16.36525    0.83473  -19.61   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.261 on 137 degrees of freedom
Multiple R-squared:  0.7375,	Adjusted R-squared:  0.7336 
F-statistic: 192.4 on 2 and 137 DF,  p-value: < 2.2e-16
&#91;/code&#93;

Adesso cominciamo a ragionare: non solo i parametri sono tutti significativamente diversi da zero per <i>p</i> < 0.05, ma anche il valore R<sup>2</sup> è molto elevato (0.7375). La prestazione in matematica, quindi, dipende sia dall'ansia verso i test di valutazione delle abilità matematiche che dal sesso.

Il modello che abbiamo appena creato può essere visualizzato, ma prima di farlo è importante capire cosa rappresentino esattamente i parametri stimati.

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">Il significato dei parametri</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

Il modello mod2 comprende tre parametri: Intercept, math.anx e genderM. Analizziamoli uno per uno.

<b>Intercept</b> è l'intercetta della retta di regressione delle femmine. Vi starete chiedendo come io faccia a saperlo... ebbene, questa informazione ci viene fornita dal comando <b>contrasts</b> applicato alla variabile gender: il livello a cui è associato il valore 0 sarà quello di <b>riferimento</b>, mentre il livello a cui è associato il valore 1 sarà quello di <b>confronto</b>. Il parametro denominato <i>Intercept</i> rappresenta l'intercetta del gruppo di riferimento, quindi delle femmine.

[code language="R"]
> contrasts(math$gender)
  M
F 0
M 1

Se non avete dimestichezza con le matrici di contrasto, credo che il modo migliore di affrontare questo passaggio sia di prenderlo come un dogma. Ci torneremo nei prossimi post.

Il parametro genderM è lo scarto fra l'intercetta dei maschi (M) e quella delle femmine. Ovvero, genderM rappresenta la quantità che dobbiamo aggiungere a Intercept per ottenere l'intercetta dei maschi. Si noti che genderM è negativo, per cui all'intercetta delle femmine dobbiamo togliere una quantità: ciò significa che l'intercetta dei maschi, come previsto, sarà inferiore a quella delle femmine.

Dato che l'intercetta delle femmine è pari a 32.23, l'intercetta dei maschi sarà pari a 32.23 + (-16.37), ovvero 15.86.

Il parametro math.anx è invece il coefficiente angolare delle due rette. Avendo in comune lo stesso coefficiente angolare, le due rette saranno parallele.

L'aver inserito nel modello il fattore genere, quindi, ha permesso di identificare due rette, parallele ma con diversa intercetta. Dato che i parametri sono tutti significativi, possiamo dire che la prestazione in matematica decresce all'aumento dell'ansia provata verso le situazioni in cui le ablità matematiche vengono sottoposte a valutazione, ma le femmine hanno di base un punteggio nelle abilità matematiche più elevato rispetto a quello dei maschi.

Di seguito è riportata una rappresentazione grafica del modello mod2, con le due rette che incrociano l'asse delle ordinate nei punti 32.23 (Intercept) e 15.86 (Intercept + genderM).

Scatterplot modello additivo

Interazione fra i predittori


Inserire nel modello l'interazione fra i due predittori significa permettere alla due rette di avere non solo diversa intercetta, ma anche diverso coefficiente angolare, per cui le rette presenteranno una pendenza diversa. Proviamo quindi a inserire nel modello anche l'effetto d'interazione.


> mod3 <- lm(math.perf ~ math.anx * gender, data=math) > summary(mod3)

Call:
lm(formula = math.perf ~ math.anx * gender, data = math)

Residuals:
Min 1Q Median 3Q Max
-10.0022 -3.1375 0.0373 2.6341 10.8452

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.8760 1.8000 21.043 < 2e-16 *** math.anx -1.2053 0.1047 -11.511 < 2e-16 *** genderM -26.5177 2.2643 -11.711 < 2e-16 *** math.anx:genderM 0.7332 0.1536 4.772 4.64e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.958 on 136 degrees of freedom Multiple R-squared: 0.7751, Adjusted R-squared: 0.7702 F-statistic: 156.3 on 3 and 136 DF, p-value: < 2.2e-16 [/code] Per prima cosa bisogna segnalare che, in questo nuovo modello, le intercette sono cambiate. Ora l'intercetta delle femmine è pari a 37.88, mentre l'intercetta dei maschi è pari a 37.88 + (-26.52) = 11.36. Inoltre, il parametro math.anx, che in precedenza era il coefficiente angolare comune a entrambe le rette, ora rappresenta il coefficiente angolare della retta del gruppo di riferimento (le femmine).

Oltre a questo, nel modello compare un nuovo parametro: math.anx:genderM. Questo parametro rappresenta la quantità da aggiungere al coefficiente angolare del gruppo di riferimento per ottenere il coefficiente angolare del gruppo di confronto (i maschi).

Quindi, dato che il coefficiente angolare delle femmine è pari a -1.21, il coefficiente angolare dei maschi sarà dato da -1.21 + 0.73, ovvero -0.48. I maschi, dunque, avranno non solo un'intercetta più bassa, ma anche un coefficiente angolare più vicino a zero.

Graficamente, la situazione è questa:

Scatterplot modello additivo

Insomma: rispetto alle femmine, i maschi presentano un'intercetta significativamente inferiore e un coefficiente angolare significativamente più vicino a zero. La relazione tra math.anx e math.perf, quindi, varia a seconda che l'individuo in questione sia un maschio oppure una femmina. Le femmine, rispetto ai maschi, hanno delle prestazioni in matematica più elevate quando l'ansia è pari a zero, ma subiscono maggiormente l'effetto negativo dell'ansia, che fa ne calare maggiormente la prestazione. L'ansia, quindi, influenza maggiormente la prestazione delle femmine.

Gestione avanzata dei fattori

Nell’ultimo post abbiamo parlato delle variabili di tipo fattore, per cercare di capire in cosa queste si differenzino dalle più semplici variabili di tipo carattere o numerico. Ora vedremo qualche funzione per manipolare i fattori, per ricodificarli e adattarli alle nostre esigenze.

Procuriamoci subito dei dati su cui lavorare. Sul portale open data del Comune di Cagliari è disponibile un dataset che contiene tutti gli incidenti stradali del 2014 avvenuti nel territorio comunale in cui sia intervenuta una pattuglia del corpo di Polizia Locale. Qui c’è la pagina di descrizione del dataset, mentre questo è il link diretto al file csv.

incidenti <- read.csv2("http://www.comune.cagliari.it/resources/cms/documents/Cagliari_PoliziaLocale_SinistriStradali_2014.csv")
&#91;/code&#93;

Il dataset contiene un bel po' di variabili (descritte sommariamente <a href="http://www.comune.cagliari.it/portale/it/opendata_visualizza_contenuto.page?modelId=45&contentId=SCH80418" target="_blank">qui</a>); noi ci concentreremo sul fattore <b>tipologia</b>, che classifica gli incidenti per tipo. Nel seguito di questo post eseguiremo una serie di operazioni su questo fattore; per evitare di alterare la variabile originale, salviamo la colonna &ldquo;tipologia&rdquo; del dataset &ldquo;incidenti&rdquo; in un nuovo oggetto, esterno al data frame:

[code language="R"]
tipoIncidenti <- incidenti$tipologia
&#91;/code&#93;

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">I livelli del fattore</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

<b>Prima domanda:</b> quanti livelli presenta il fattore tipologia? La risposta ce la fornisce il comando <b>nlevels</b>:

[code language="R"]
> nlevels(tipoIncidenti)
[1] 5

Seconda domanda: quali sono i livelli del fattore tipologia? La risposta ce la fornisce il comando levels, che estrae il vettore di etichette:

> levels(tipoIncidenti)
[1] "Danneggiamento"                         
[2] "Incidente stradale con feriti"          
[3] "Incidente stradale mortale"             
[4] "Incidente stradale per caduta di pedone"
[5] "Incidente stradale senza feriti"

Possiamo anche estrarre ogni singolo livello utilizzando le parentesi quadre. Per esempio, per estrarre il quarto livello:

> levels(tipoIncidenti)[4]
[1] "Incidente stradale per caduta di pedone"

Il fattore raggruppa gli incidenti sulla base di cinque categorie: la prima comprende generici danneggiamenti derivanti da buche sul manto stradale, caduta di rami dagli alberi, urto contro cassonetti, ecc., mentre le altre categorie comprendono incidenti stradali veri e propri.

Terza domanda: quante osservazioni sono presenti per ogni livello del fattore? La risposta ce la fornisce il comando table:

> table(tipoIncidenti)
tipoIncidenti
                         Danneggiamento           Incidente stradale con feriti 
                                    170                                     571 
             Incidente stradale mortale Incidente stradale per caduta di pedone 
                                      5                                      66 
        Incidente stradale senza feriti 
                                    634

Nel caso in cui le etichette che descrivono i nomi dei fattori siano lunghe (come in questo caso), possiamo associare il comando cbind al comando table, in modo da avere un output più ordinato:

> cbind(table(tipoIncidenti))
                                        [,1]
Danneggiamento                           170
Incidente stradale con feriti            571
Incidente stradale mortale                 5
Incidente stradale per caduta di pedone   66
Incidente stradale senza feriti          634

Ricodificare i livelli


Le etichette che costituiscono i nomi dei livelli sono piuttosto lunghe. Soprattutto quando si vogliono visualizzare i dati in un grafico, è bene che le etichette siano quanto più concise possibile. Possiamo rinominare i livelli sfruttando sempre il comando levels, utilizzandolo per attribuire al fattore un nuovo vettore di nomi:

newlev <- <- c("Danni", "con feriti", "con morti", "caduta pedone", "senza feriti")
levels(tipoIncidenti) <- newlev
&#91;/code&#93;

Attenzione: all'interno del vettore <i>newlev</i>, i nomi devono essere disposti nello stesso ordine in cui sono disposti i livelli nel fattore, pena l'attribuire ai livelli le etichette sbagliate.

Verifichiamo il risultato:

[code language="R"]
> levels(tipoIncidenti)
[1] "Danni"         "con feriti"    "con morti"     "caduta pedone"
[5] "senza feriti"

> table(tipoIncidenti)
tipoIncidenti
        Danni    con feriti     con morti caduta pedone  senza feriti 
          170           571             5            66           634

Possiamo ricodificare i livelli anche singolarmente. Per esempio, per trasformare l'etichetta “Danni” in “danni”:

> levels(tipoIncidenti)[1] <- "danni"
> levels(tipoIncidenti)
[1] "danni"         "con feriti"    "con morti"     "caduta pedone"
[5] "senza feriti"

Accorpare i livelli


Immaginiamo di voler raggruppare i sinistri in due categorie: quelli con feriti o morti (“tipo 1”) e quelli senza alcuna conseguenza fisica per le persone coinvolte (“tipo 2”). Per eseguire questa operazione possiamo sfruttare ancora una volta la funzione levels, attribuendo la medesima etichetta ai livelli che vogliamo accorpare:

levels(tipoIncidenti) <- c("tipo 2", "tipo 1", "tipo 1", "tipo 2", "tipo 2")
&#91;/code&#93;

Verifichiamo il risultato:

&#91;code language="R"&#93;
> levels(tipoIncidenti)
[1] "tipo 2" "tipo 1"

> table(tipoIncidenti)
tipoIncidenti
tipo 2 tipo 1 
   870    576

Come si vede, ci sono 576 incidenti di tipo 1 (somma di “con feriti” e “con morti”) e 870 incidenti di tipo 2 (somma di “danni”, “caduta pedone” e “senza feriti”).

Riordinare i livelli


La precedente aggregazione dei livelli ha prodotto un risultato per alcuni versi sconveniente: l'etichetta “tipo 2” costituisce il primo livello, mentre l'etichetta “tipo 1” il secondo livello. Ecco cosa succede se proviamo a visualizzare in un grafico a barre il numero di incidenti per ogni tipologia:

> barplot(table(tipoIncidenti))

barplot incidenti

I livelli vengono visualizzati nell'ordine con cui sono stati codificati, per cui il tipo 2 compare prima del tipo 1. Per ovviare a questo problema possiamo usare la funzione factor, che consente di rifattorizzare una variabile:

> tipoIncidenti <- factor(tipoIncidenti, levels=c("tipo 1","tipo 2"))
&#91;/code&#93;

Adesso i livelli sono stati riordinati, per cui sarà &ldquo;tipo 1&rdquo; a comparire per primo:

&#91;code language="R"&#93;
> table(tipoIncidenti)
tipoIncidenti
tipo 1 tipo 2 
   576    870
> barplot(table(tipoIncidenti))

barplot incidenti

Si noti che la stessa funzione factor, al pari della funzione levels, può essere usata per modificare i nomi dei livelli:

> tipoIncidenti <- factor(tipoIncidenti, labels=c("Tipo 1","Tipo 2"))
&#91;/code&#93;

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">Fattori con livelli ordinati</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

Per come l'abbiamo concepito fino a questo momento, un fattore non è altro che una variabile qualitativa su <b>scala nominale</b>. In realtà, il tipo fattore è pensato anche per registrare dati misurati su <b>scala ordinale</b>. Ripartiamo dall'inizio e riprendiamo l'originale vettore &ldquo;tipologia&rdquo; del dataset &ldquo;incidenti&rdquo;, salviamolo nell'oggetto tipoIncidenti abbreviando i nomi dei livelli:

[code language="R"]
tipoIncidenti <- factor(incidenti$tipologia, labels=c("danni", "con feriti", "con morti", "caduta pedone", "senza feriti"))
&#91;/code&#93;

Ammettiamo che sia possibile dare un ordinamento alle categorie del fattore sulla base della gravità: al livello più basso ci sarà la categoria "Danneggiamento", poi "Incidente stradale per caduta di pedone", poi "Incidente stradale senza feriti", quindi "Incidente stradale con feriti" e infine "Incidente stradale mortale".

Per prima cosa, creiamo un vettore che contiene i nomi dei livelli disposti nell'ordine desiderato:
&#91;code language="R"&#93;
lev <- c("danni", "caduta pedone", "senza feriti", "con feriti", "con morti")
&#91;/code&#93;

Si noti che avremmo potuto ottenere lo stesso risultato anche utilizzando il comando levels:

&#91;code language="R"&#93;
> lev <- levels(tipoIncidenti)&#91;c(1,4,5,2,3)&#93;
> lev
[1] "danni"         "caduta pedone" "senza feriti"  "con feriti"   
[5] "con morti"

A questo punto, rifattorizziamo il vettore utilizzando il comando factor, impostando a TRUE l'argomento ordered:


tipoIncidenti <- factor(tipoIncidenti, levels=lev, ordered=TRUE) [/code] Il comando str, adesso, indicherà che tipoIncidenti è una variabile di tipo fattore, ma i cui livelli sono ordinati, dove danni < caduta pedone < senza feriti < con feriti < con morti. [code language="R"] > str(tipoIncidenti)
Ord.factor w/ 5 levels "danni"<"caduta pedone"<..: 3 3 3 3 3 1 2 4 3 4 .. [/code] Si noti che lo stesso risultato avremmo potuto ottenerlo con la funzione ordered:


> tipoIncidenti <- ordered(tipoIncidenti, levels=lev) [/code] L'ordinamento dei livelli di un fattore è una sottigliezza che assume una sua utilità esclusivamente nell'adattamento di un modello statistico. È comunque importante che una variabile misurata su scala ordinale sia codificata nel modo corretto, soprattutto quando ci si passa i dati tra colleghi e certi dettagli rischiano di essere omessi e talvolta male interpretati.

Alcune cose da sapere sui fattori

Quando importiamo un file di dati in R, tutti i valori che il software identifica come di tipo carattere vengono automaticamente convertiti in tipo fattore. Alcune funzioni di lettura dati, come read.table o read.csv, possiedono un argomento chiamato stringsAsFactors, che, se impostato a FALSE, evita questa conversione.

Una domanda che mi viene posta abbastanza di frequente da chi non è solito usare i modelli statistici è: «che cosa sono i fattori?». La prima volta che ho dovuto a rispondere a tale quesito, devo ammetterlo, mi sono trovato un po’ in difficoltà. Sembra assurdo, ma per chi si occupa di statistica inferenziale (a qualsiasi livello) il concetto è così tanto dato per scontato da rendere complicato trovare le parole giuste per descriverlo a chi non ha mai fatto statistica o non ha mai avuto a che fare con il disegno di un esperimento.

La cosa più difficile da spiagare è quale sia la differenza tra variabili di tipo carattere e variabili di tipo fattore. Infatti, quando i dati vengono importati in R, i caratteri – e solo i caratteri – vengono convertiti in “factor”. Ma che cosa cambia? Non avremmo potuto lasciare tutto com’era e lavorare con le variabili “character” invece che “factor”? Tanto, sempre di stringhe si tratta!

Invece no, non si tratta sempre di stringhe, e adesso proverò a spiegarvi perché.

La differenza tra caratteri e fattori


Il fattore è un tipo di dato che permette di etichettare un gruppo di osservazioni utilizzando delle categorie predefinite. Immaginiamo di aver registrato in una tabella alcuni dati relativi a sette personaggi Disney residenti a Paperopoli e dintorni:

Nome Indirizzo Genere
Paperina Vico II dei ciclamini, 2 F
Paperoga Piazza C. Barks, 33 M
Nonna Papera Viale degli uliveti, s.n. F
Gastone Via della Dea Bendata, 1 M
Paperino Via della marineria, 8 M
Brigitta Via dei ciclamini, 21 F
Zio Paperone Colle del deposito, s.n. M

Questa tabella di dati contiene tre variabili: il nome del personaggio, l’indirizzo di residenza e il sesso. In tutti e tre i casi, le osservazioni sono codificate attraverso dei caratteri, ma solo una di queste variabili potrà essere utilizzata come fattore: il genere. Il nome e l’indirizzo, infatti, rappresentano dei dati univoci e caratteristici di ogni singolo personaggio. Se dovessimo raggruppare i personaggi, di certo non potremo usare l’indirizzo come criterio di raggruppamento, visto che ogni personaggio abita in un luogo diverso, né tantomeno il nome, univoco per definizione. Il genere, al contrario, è una variabile che può assumere solo due valori (F: femmina, M: maschio) e consente di suddividere i personaggi in due gruppi, uno di tre unità (le femmine) e uno di quattro unità (i maschi).

La prima differenza tra il tipo fattore e il tipo carattere consiste quindi nel fatto che il fattore serve per raggruppare le unità sottoposte a osservazione, mentre la semplice variabile di tipo carattere no. Il fattore assume un numero limitato di modalità, dette livelli, che specificano il gruppo di appartenenza di un’unità/individuo.

La seconda differenza tra il tipo fattore e il tipo carattere è che i dati che andranno a rappresentare una variabile fattore non necessariamente devono essere registrati utilizzando dei caratteri. Nell’esempio dei personaggi Disney, per codificare il genere ho usato delle etichette, ovvero F e M, ma queste etichette sono state scelte esclusivamente a mia personalissima discrezione. Nessuno mi avrebbe impedito di scrivere per esteso “femmina” o “maschio” o utilizzare addirittura dei numeri, magari 1 per indicare le femmine e 0 i maschi. Sarebbe invece stato più complicato indicare l’indirizzo attraverso dei numeri e sicuramente i caratteri erano d’obbligo per specificare il nome del personaggio.

La gestione dei fattori in R


Il tipo fattore assume un senso e gioca un ruolo fondamentale quando si deve adattare un modello statistico. Si immagini di voler eseguire un confronto tra gruppi di individui, suddivisi sulla base di una qualche caratteristica. R deve sapere ogni individuo a quale gruppo appartiene, per cui in fase di analisi dei dati sarà necessario specificare una variabile (fattore, ovviamente) che indica al software a quale raggruppamento appartiene ogni unità sottoposta a misurazione.

Da questo indirizzo è possibile scaricare in formato csv il dataset descritto nella precedente tabella. Per importare i dati in R dobbiamo prima di tutto scaricarli in locale e quindi leggerli usando la funzione read.csv2, con l’accorgimento di impostare l’argomento stringsAsFactors come FALSE; saremo noi, in seguito all’importazione, a decidere cosa convertire in fattore e cosa no.

disney <- read.csv2("disney.csv", stringsAsFactors=FALSE)
&#91;/code&#93;

Visualizziamo la struttura del dataset:

&#91;code language="R"&#93;
> str(disney)
'data.frame':	7 obs. of  4 variables:
 $ Nome     : chr  "Paperina" "Paperoga" "Nonna Papera" "Gastone" ...
 $ Indirizzo: chr  "Vico II dei ciclamini, 2" "Piazza C. Barks, 33" "Viale degli uliveti, s.n." "Via della Dea Bendata, 1" ...
 $ Genere   : chr  "F" "M" "F" "M" ...
 $ Sesso    : int  1 0 1 0 0 1 0

Rispetto alla tabella, in questo dataset è presente un’ulteriore variabile, Sesso, che codifica il genere utilizzando dei numeri: 1 per le femmine e 0 per i maschi.

La funzione per ricodificare una variabile come fattore si chiama as.factor. Iniziamo con il ricodificare la variabile di tipo carattere Genere:

disney$Genere <- as.factor(disney$Genere)
&#91;/code&#93;

Vediamo quindi l'esito della ricodifica:

&#91;code language="R"&#93;
> str(disney$Genere)
 Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2

> disney$Genere
[1] F M F M M F M
Levels: F M

Il comando str, che descrive la struttura della variabile, ci informa che Genere adesso è un fattore a due livelli: “F” e “M”. Anche quando visualizziamo tutti i valori contenuti nel vettore, di seguito alla sequenza di valori (alla voce Levels), R stampa l’elenco delle etichette che costituiscono i livelli dei fattori.

Soffermiamoci un attimo sull’esito dell’applicazione del comando str sulla variabile fattore. Di seguito all’elenco dei livelli, “F” e”M”, compare una sequenza di numeri: 1 2 1 2 2 1 2. R codifica i livelli di un fattore numerandoli; il numero 1 viene assegnato al primo livello, il numero 2 al secondo livello, il numero 3 al terzo livello, e così via.

La scelta dei numeri è strettamente legata all’ordine alfabetico delle etichette utilizzate per rappresentare i livelli. In ordine alfabetico, l’etichetta “F” viene prima dell’etichetta “M”, per cui a “F” sarà assegnato il numero 1 e a “M” il numero 2. Per questo motivo, trasformare una variabile di tipo fattore in una variabile di tipo numerico è un’operazione concessa in R:

> as.numeric(disney$Genere)
[1] 1 2 1 2 2 1 2

Ovviamente, nell’operazione di conversione da tipo fattore a tipo numerico, le etichette “F” verranno rimpiazzate dal numero 1 (perché F è il primo livello) e le etichette “M” verranno rimpiazzate dal numero 2 (perché M è il secondo livello).

Proviamo ora a convertire in fattore la variabile numerica Sesso:

> disney$Sesso <- as.factor(disney$Sesso)
> disney$Sesso
[1] 1 0 1 0 0 1 0
Levels: 0 1

Nella variabile fattore Sesso, 0 (che codifica i maschi) è il primo livello, mentre 1 (che codifica le femmine) è il secondo livello. Questo accade semplicemente perché, nella sequenza numerica, 0 viene prima di 1. Ora, quindi, è la categoria “maschio” che costituisce il primo livello, mentre la categoria “femmina” il secondo.

Si noti che, nonostante i dati siano stati etichettati con dei numeri, sui fattori non è consentito applicare alcuna funzione matematica:

> mean(disney$Sesso)
[1] NA
Warning message:
In mean.default(disney$Sesso) :
  argument is not numeric or logical: returning NA

Riconvertire il tipo fattore in tipo numerico


La rappresentazione numerica che R utilizza per codificare i fattori rischia di creare molta confusione. La categoria che noi abbiamo etichettato con 0 viene codificata da R con il numero 1, mentre la categoria che noi abbiamo etichettato con 1 viene codificata da R con il numero 2. Se volessimo riconvertire il fattore Sesso in una variabile numerica, rischieremmo di combinare un vero disastro. Questo è infatti l’esito della conversione:

> as.numeric(disney$Sesso)
[1] 2 1 2 1 1 2 1

Nel riconvertire in formato numerico la variabile, R non utilizza le etichette che noi abbiamo attribuito ai livelli, ma la codifica numerica che lui ha dato ai livelli. Come fare allora?

La soluzione è semplice: per riconvertire una variabile fattore in tipo numerico, prima la si trasforma in carattere e solo succcessivamente le stringhe potranno essere convertite in numeri. Ovvero:

> disney$Sesso <- as.character(disney$Sesso)
> disney$Sesso
[1] "1" "0" "1" "0" "0" "1" "0"
> disney$Sesso <- as.numeric(disney$Sesso)
> disney$Sesso
[1] 1 0 1 0 0 1 0

Oppure, in un unico comando:


disney$Sesso <- as.numeric(as.character(disney$Sesso)) [/code] Per ora concludiamo qui: direi che di roba da digerire ce n'è a sufficienza!