Tutorial sulle tabelle di contingenza in R

Come creare una tabella

Per prima cosa, procuriamoci dei dati. Il pacchetto MASS contiene dati su 93 auto in vendita negli USA nel 1993. Sono memorizzati nell'oggetto Cars93 e includono 27 caratteristiche per ogni auto, alcune delle quali sono categoriche. Quindi carichiamo il pacchetto MASS e guardiamo il tipo di veicoli inclusi in cars93:

library(MASS)
Cars93$Type

##  [1] Small   Midsize Compact Midsize Midsize Midsize Large   Large  
##  [9] Midsize Large   Midsize Compact Compact Sporty  Midsize Van    
## [17] Van     Large   Sporty  Large   Compact Large   Small   Small  
## [25] Compact Van     Midsize Sporty  Small   Large   Small   Small  
## [33] Compact Sporty  Sporty  Van     Midsize Large   Small   Sporty
## [41] Sporty  Small   Compact Small   Small   Sporty  Midsize Midsize
## [49] Midsize Midsize Midsize Large   Small   Small   Compact Van    
## [57] Sporty  Compact Midsize Sporty  Midsize Small   Midsize Small  
## [65] Compact Van     Midsize Compact Midsize Van     Large   Sporty
## [73] Small   Compact Sporty  Midsize Large   Compact Small   Small  
## [81] Small   Compact Small   Small   Sporty  Midsize Van     Small  
## [89] Van     Compact Sporty  Compact Midsize
## Levels: Compact Large Midsize Small Sporty Van

Abbiamo 6 tipi di auto. La funzione table indica quante ne abbiamo per ciascun tipo:

table(Cars93$Type)

##
## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

prop.table lo converte in frazioni:

prop.table(table(Cars93$Type))

##
##    Compact      Large    Midsize      Small     Sporty        Van
## 0.17204301 0.11827957 0.23655914 0.22580645 0.15053763 0.09677419

Lo stesso con l'origine delle auto:

table(Cars93$Origin)

##
##     USA non-USA
##      48      45

prop.table(table(Cars93$Origin))

##
##      USA  non-USA
## 0.516129 0.483871

Come creare una tabella di contingenza

Bene, abbiamo visto che il nostro dataset contiene un numero simile di auto USA e non USA e che i tipi più diffusi sono Midsize e Small. Tuttavia, potrebbe esserci una differenza di tipologia tra auto USA e non USA?

Guardiamo i tipi di auto in relazione alla loro origine. Possiamo usare di nuovo table, ma ora con due argomenti. Il primo diventerà la variabile di riga e il secondo la variabile di colonna:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

Ora abbiamo visto ciò che tutti sanno: gli americani adorano i veicoli grandi! La tabella sopra mostra la distribuzione congiunta di due variabili categoriche (Type e Origin). Tali tabelle si chiamano tabelle di contingenza.

Come ottenere i marginali da una tabella di contingenza

Le funzioni rowSums e colSums si spiegano da sole

(tab1<-table(Cars93$Type, Cars93$Origin))

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

rowSums(tab1)

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

colSums(tab1)

##     USA non-USA
##      48      45

Come ottenere le percentuali da una tabella di contingenza

prop.table annidato con table restituisce le frequenze

prop.table(table(Cars93$Type, Cars93$Origin))

##          
##                  USA    non-USA
##   Compact 0.07526882 0.09677419
##   Large   0.11827957 0.00000000
##   Midsize 0.10752688 0.12903226
##   Small   0.07526882 0.15053763
##   Sporty  0.08602151 0.06451613
##   Van     0.05376344 0.04301075

La conversione in percentuali è semplicemente una moltiplicazione per 100:

prop.table(table(Cars93$Type, Cars93$Origin))*100

##          
##                 USA   non-USA
##   Compact  7.526882  9.677419
##   Large   11.827957  0.000000
##   Midsize 10.752688 12.903226
##   Small    7.526882 15.053763
##   Sporty   8.602151  6.451613
##   Van      5.376344  4.301075

Nota che questa è una distribuzione di probabilità congiunta, dalla quale possiamo vedere che, ad esempio, circa il 7,5% delle auto è di piccole dimensioni e di origine americana.

Più spesso, ci interessa la distribuzione di una variabile all'interno dei gruppi creati da un'altra. Qui sembra interessante la distribuzione dei tipi di auto tra le auto USA e, separatamente, quelle non USA. Per ottenerla, usiamo l'argomento margin della funzione prop.table. Indica dove si trova la variabile di raggruppamento: nelle righe (margin=1) o nelle colonne (margin=2):

prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

Ora possiamo vedere facilmente che le auto piccole sono due volte più frequenti tra le non USA che tra le USA nella nostra base dati.

Nota anche che le percentuali sommano a 100 nelle colonne, mentre nella tabella della distribuzione congiunta (quella senza l'argomento margin) 100 era la somma dell'intera tabella.

(tab2<-prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100)

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

colSums(tab2)

##     USA non-USA
##     100     100

Test chi-quadrato

La domanda più comune che nasce dalle tabelle di contingenza è se le variabili di riga e di colonna siano indipendenti. Il modo più basilare per rispondervi è eseguire un test chi-quadrato. È trattato in grande dettaglio in questo tutorial. Verifichiamo se Type e Origin sono indipendenti:

chisq.test(Cars93$Type, Cars93$Origin)

## Warning in chisq.test(Cars93$Type, Cars93$Origin): Chi-squared
## approximation may be incorrect

##
##  Pearson's Chi-squared test
##
## data:  Cars93$Type and Cars93$Origin
## X-squared = 14.08, df = 5, p-value = 0.01511

A quanto pare non lo sono, ma abbiamo anche ottenuto l'avviso Chi-squared approximation may be incorrect. Questo perché la statistica chi-quadrato segue la distribuzione chi-quadrato solo approssimativamente. Più osservazioni abbiamo, migliore è l'approssimazione. La funzione chisq.test mostra l'avviso sopra ogni volta che uno dei conteggi attesi è inferiore a 5 (per cosa sia un "conteggio atteso", vedi il tutorial linkato sopra).

Test esatto di Fisher

Il test esatto di Fisher è un'alternativa al test chi-quadrato usata soprattutto quando l'approssimazione chi-quadrato non è soddisfacente. Eseguiamolo:

fisher.test(Cars93$Type, Cars93$Origin)

##
##  Fisher's Exact Test for Count Data
##
## data:  Cars93$Type and Cars93$Origin
## p-value = 0.007248
## alternative hypothesis: two.sided

I risultati sono piuttosto simili a quelli del chi-quadrato, ma non è sempre così. Un grande svantaggio del test esatto di Fisher è che per tabelle grandi (o campioni grandi) diventa computazionalmente inefficiente.

G-test

Un'altra alternativa è il cosiddetto G-test. La sua statistica è anch'essa approssimativamente distribuita come chi-quadrato, ma per campioni piccoli questa approssimazione è più accurata di quella usata dal test chi-quadrato. Per il G-test possiamo usare la funzione GTest del pacchetto DescTools. I risultati sono di nuovo abbastanza simili ai due test precedenti: Type e Origin non sono indipendenti.

library(DescTools)
GTest(Cars93$Type, Cars93$Origin)

##
##  Log likelihood ratio (G-test) test of independence without
##  correction
##
## data:  Cars93$Type and Cars93$Origin
## G = 18.362, X-squared df = 5, p-value = 0.002526

Correzione di Yates

Nelle tabelle di contingenza 2x2, il test chi-quadrato può essere migliorato con la correzione di continuità di Yates. Semplicemente sottrae 0,5 da ciascun termine | Osservato - Atteso | nella statistica chi-quadrato. Fai di nuovo riferimento a questo tutorial se ora ti senti spaesato. Inoltre, R applica automaticamente la correzione di Yates quando necessario. Guardiamo la disponibilità di versioni con cambio manuale tra auto USA e non USA:

(tab3<-table(Cars93$Man.trans.avail, Cars93$Origin))

##      
##       USA non-USA
##   No   26       6
##   Yes  22      39

chisq.test(tab3)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  tab3
## X-squared = 15.397, df = 1, p-value = 8.712e-05

La correzione di Yates è stata applicata automaticamente, e R ce l'ha comunicato.

Tabella a 3 (o più) dimensioni

Guardiamo insieme Man.trans.avail, Origin e Type. table suddivide il dataset in base alla variabile fornita come terza:

table(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

## , ,  = Compact
##
##      
##       USA non-USA
##   No    2       0
##   Yes   5       9
##
## , ,  = Large
##
##      
##       USA non-USA
##   No   11       0
##   Yes   0       0
##
## , ,  = Midsize
##
##      
##       USA non-USA
##   No    9       4
##   Yes   1       8
##
## , ,  = Small
##
##      
##       USA non-USA
##   No    0       0
##   Yes   7      14
##
## , ,  = Sporty
##
##      
##       USA non-USA
##   No    0       0
##   Yes   8       6
##
## , ,  = Van
##
##      
##       USA non-USA
##   No    4       2
##   Yes   1       2

ftable offre una vista più compatta:

ftable(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##              Compact Large Midsize Small Sporty Van
##                                                    
## No  USA            2    11       9     0      0   4
##     non-USA        0     0       4     0      0   2
## Yes USA            5     0       1     7      8   1
##     non-USA        9     0       8    14      6   2

Test di Cochran-Mantel-Haenszel

Dai risultati di table sopra, dovrebbe essere chiaro che la relazione tra Origin e Man.trans.avail varia con Type. Per le auto small e sporty, ad esempio, non c'è associazione: ogni auto small o sporty USA, così come ogni non USA, ha una versione manuale. D'altra parte, la maggior parte delle auto Midsize USA non ha la versione manuale, mentre la maggior parte delle non USA di questo tipo sì. Per tenere conto di possibili relazioni diverse tra strati, si può usare il test di Cochran-Mantel-Haenszel:

mantelhaen.test(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##
##  Mantel-Haenszel chi-squared test with continuity correction
##
## data:  Cars93$Man.trans.avail and Cars93$Origin and Cars93$Type
## Mantel-Haenszel X-squared = 8.0153, df = 1, p-value = 0.004638
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##   2.226531 76.891307
## sample estimates:
## common odds ratio
##          13.08438

Il terzo argomento passato a mantelhaen.test identifica gli strati. Confronta i risultati sopra con quelli senza stratificazione (l'esempio sulla correzione di Yates sopra). L'associazione è ancora presente, ma l'evidenza a suo favore è molto più debole.

Misure di associazione

Una volta scoperte alcune associazioni tra variabili, è il momento di misurarne l'intensità. Esiste una pletora di possibili misure. Molte di esse sono descritte qui. Ora mi concentrerò sulle due più utilizzate.

V di Cramer

V si basa sulla statistica chi-quadrato:

$$ V = \sqrt{\frac{\chi^2/N}{\min(C-1, R-1)}}, $$ dove:

• N è il totale complessivo della tabella di contingenza (somma di tutte le sue celle),
• C è il numero di colonne
• R è il numero di righe.

V ∈ [0; 1]. Quanto più grande è V, tanto più forte è la relazione tra le variabili. V = 0 può essere interpretato come indipendenza (poiché V = 0 se e solo se χ² = 0). Il principale svantaggio di V è la mancanza di un'interpretazione precisa. V = 0,6 è un'associazione forte, moderata o debole?

La funzione CramerV di DescTools può calcolarlo per noi:

CramerV(Cars93$Type, Cars93$Origin)

## [1] 0.3890967

Ancora una volta: si tratta di un'associazione forte, moderata o debole?

Lambda di Goodman e Kruskal

La lambda di Goodman e Kruskal è un esempio di misura basata sulla riduzione proporzionale della variazione. Questi tipi di misure mirano a imitare R² - il coefficiente di determinazione della regressione lineare:

• assumono valori in [0, 1]
• esprimono la frazione di variazione spiegata dalla variabile indipendente.

Usiamo la variabile di colonna come indipendente. La formula della lambda di Goodman e Kruskal è quindi: $$\lambda = \frac{L - \sum_j L_j}{L},$$ dove:

• L_j è la somma delle frequenze non modali nella j-esima colonna,
• L è la somma delle frequenze non modali nella colonna "totale"

Questo si capisce meglio con un esempio. Guardiamo Type e Origin, la seconda come indipendente:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

La moda della colonna USA è 11, quindi L₁ = 7 + 10 + 7 + 8 + 5 = 37. La moda della colonna non USA è 14, quindi L₂ = 9 + 0 + 12 + 6 + 4 = 31. La colonna "Totale" è

rowSums(table(Cars93$Type, Cars93$Origin))

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

La moda è 22, quindi L = 16 + 11 + 21 + 14 + 9 = 71 e $$\lambda = \frac{71 - (37+31)}{71}=0.042$$

Possiamo invece usare la funzione Lambda del pacchetto DescTools:

Lambda(Cars93$Type, Cars93$Origin, direction='row')

## [1] 0.04225352

Il parametro direction determina dove si trova la variabile dipendente (in row o in column).

Origin spiega solo circa il 4% della variazione in Type. Nota dalla formula che la lambda definisce la variazione come dicotomia tra appartenere e non appartenere al gruppo più numeroso.

Vale la pena menzionare che la lambda è zero ogni volta che la categoria modale di ciascuna colonna è la stessa. Dai un'occhiata a questa tabella:

(lambdaTab<-cbind(c(0,100), c(49,51)))

##      [,1] [,2]
## [1,]    0   49
## [2,]  100   51

chisq.test(lambdaTab)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  lambdaTab
## X-squared = 62.279, df = 1, p-value = 2.981e-15

La categoria modale è la stessa per ciascuna colonna (è la seconda riga). Quindi la lambda deve essere zero nonostante un'associazione significativa e visibile:

Lambda(lambdaTab, direction='row')

## [1] 0

Se vuoi saperne di più sulla statistica in R, segui il corso Statistical Modeling in R (Parte 1) di DataCamp.