Variabili nominali
Statistiche sulle variabili nominali
A partire da una variabile nominale è possibile ottenere un numero limitato di statistiche descrittive, univariate. Quello che è possibile fare, in pratica, è creare una tabella di contingenza unidimensionale, ovvero un vettore la cui lunghezza è pari al numero di livelli, e dove il valore di ogni cella è pari al numero di osservazioni che appartengono alla corrispondente categoria.
A partire da questa rappresentazione, è possibile calcolare l'indice centrale della moda.
Distribuzioni categoriali
Da un punto di vista inferenziale, è possibile utilizzare la distribuzione delle osservazioni del campione come stima della frequenza di ogni categoria nella popolazione.
Confronto di una distribuzione campionaria con una distribuzione teorica
Vi sono circostanze, in cui l'ipotesi di ricerca è finalizzata a confrontare la distribuzione categoriale di un campione, rispetto ad una distribuzione teorica riguardante una popolazione.
Un esempio: distribuzione occupati
Facciamo un esempio. Assumiamo che, in Italia, il 17% della popolazione occupata lavori nell'agricoltura, il 51% nell'industria, e il rimanente 32% nel terziario (le cifre sono assolutamente inventate).
Potremmo chiederci se, in una determinata provincia, la distribuzione categoriale si discosta o meno dalla distribuzione ipotizzata.
Per fare questo, possiamo estrarre un campione dalla popolazione provinciale degli occupati, creare la relativa tabella di contingenza, e capire se la distribuzione del campione rispetta le cosiddette frequenze attese.
Calcolo frequenze attese
Il calcolo delle frequenze attese è semplice: basta moltiplicare la probabilità attesa (ovvero la probabilità che, data la popolazione di riferimento, venga estratta una osservazione appartenente ad una determinata categoria) per la numerosità del campione.
Se immaginiamo di utilizzare un campione di 100 persone, la frequenza attesa sarà di 17 persone che lavorano nell'agricoltura, 51 nell'industria, 33 nel terziario.
Stima dell'errore
Errore di campionamento
A causa dell'errore non sistematico di campionamento, però, difficilmente le frequenze osservate saranno uguali alle frequenze attese.
Il compito della statistica inferenziale, in questo caso, sarà quello di stabilire se la differenza fra le frequenze osservate e quelle attese sono da attribuire o meno al caso. Nella prima ipotesi, non si rifiuta l'ipotesi nulla, secondo cui non vi è differenza significativa fra le frequenze osservate e quelle attese, e dunque si assume che non vi sia differenza fra la distribuzione della provincia in esame e quella nazionale, di riferimento.
Nella seconda ipotesi, si rifiuterà l'ipotesi nulla, e di conseguenza si assumerà che vi è una differenza significativa fra la distribuzione di frequenza del campione e quella della popolazione.
Nel nostro esempio, si assumerà che una differenza fra il nostro campione e la popolazione nazionale di riferimento, e pertanto che la distribuzione di frequenza nelle tre categorie è, nella popolazione provinciale, diversa da quella nazionale.
Stimare la differenza fra le fequenze
Per fare questo, abbiamo bisogno di una misura che ci permetta di calcolare la differenza fra due tabelle.
Informalmente, abbiamo bisogno di una misura che stimi la distanza fra due tabelle, e dunque che sia pari a zero se le due tabelle sono uguali, che sia positiva se le due tabelle sono differenti, e che cresca al crescere delle differenze.
In letteratura, vengono citate tre possibili misure.
il $\chi^2$ di Pearson \index{$\chi^2$!Pearson}: \begin\label \chi^{2} = \sum^r_ \sum^c_ \frac{(O_-E_)^2}{E_} \end
Il likelihood ratio chi square \index{$\chi^2$!likelihood ratio chi square}: \begin\label G^{2} = 2 \sum^r_ \sum^c_ (O_ \times log(\frac{O_}{E_})) \end
il test di Fisher.
La misura più nota, che utilizzeremo, è il $\chi^2$ di Pearson \ref.
Numerosità delle frequenze attese
Com'è intuibile, sia l'equazione \ref che la \ref sono poco adatte a circostanze in cui la frequenza attesa di una delle celle è molto bassa (generalmente, si assume che 5 sia la frequenza attesa minima).
La simulazione
Per introdurre la stima della probabilità che la distribuzione di frequenze di un campione sia significativamente diverso dalla distribuzione attesa, useremo di nuovo il metodo della simulazione.
Più in particolare, andiamo a generare k campioni di una variabile categoriale con una determinata frequenza attesa, e per ogni campione misuriamo, utilizzando l'equazione \ref, la distanza fra la distribuzione osservata e quella attesa.
La distribuzione delle distanze dei k campioni ci permette di stimare la distribuzione dell'errore di campionamento, di stabilire i valori critici che corrispondono ad un errore $\alpha$, e al calcolo del p-value di una distribuzione.
Generiamo un'urna
Generiamo un'urna, di 100 valori, da cui estrarre i campioni.
# 17 in agricoltura, 51 industria, 32 terziario
atteso <- c(rep(1,17),rep(2,51),rep(3,32))
length(atteso)
## [1] 100
t_atteso <- table(atteso)
(prob_attesa <- t_atteso/length(atteso))
## atteso
## 1 2 3
## 0.17 0.51 0.32
Di fatto, in questa simulazione utilizziamo il metodo delle permutazioni \index{Simulazione!permutazioni}.
Generiamo i campioni
Generiamo 10.000 tabelle da 100 elementi, estratti (con ripetizione) dall'urna sopra creata. Per ogni tabella, calcoliamo la distanza dalla tabella di riferimento (t_atteso) usando la formula del $\chi^{2}$, e la salviamo in un vettore, distanza_oss (distenza osservata).
R: generazione dei campioni
# creo il vettore dove salvare distanze osservate
distanza_oss <-vector(mode = "numeric", length = 10000)
for (ciclo in 1:length(distanza_oss)) {
campione <- sample(atteso,100,replace=TRUE)
t_campione <- table(campione)
chi_quadro <- sum(((t_campione-t_atteso)^2)/t_atteso)
distanza_oss [ciclo] <- chi_quadro
}
Disegno il grafico della distribuzione della distanza osservata.
hist(distanza_oss,freq=FALSE,breaks=20)
R: Valori critici
calcoliamo i valori critici per $\alpha = 0.05$ e 0.01, ovvero il $95^o$ ed il $99^o$ percentile.
quantile(distanza_oss, probs = c(0.95,0.99))
## 95% 99%
## 5.825368 9.027298
Valori critici e inferenza
Questo calcolo ci permette di inferire che possiamo rifiutare l'ipotesi nulla quando la distanza della tabella osservata da quella attesa (misurata con la formula \ref) è > 5.88 (con $\alpha = 0.05$) o > 9.47 (con $\alpha = 0.01$)
Possiamo dunque generare un nuovo campione, calcolare la distanza dalla distribuzione attesa, e confrontarla con i valori critici.
Generazione di un nuovo campione
Generiamo un nuovo campione, calcoliamo la tabella di contingenza, e calcoliamo il $\chi^{2}$. Poiché il campione è generato a partire dall'urna, ci aspettiamo che la statistica calcolata non sia significativa.
campione <- sample(atteso,100,replace=TRUE)
t_campione <- table(campione)
chi_quadro1 <- sum(((t_campione-t_atteso)^2)/t_atteso)
chi_quadro1
## [1] 0.09007353
Stima del p-value
Inoltre, possiamo stimare il p-value, ovvero la probabilità di compiere un errore di tipo I rifiutando l'ipotesi nulla. Per fare questo, basta aggiungere la distanza del nuovo campione al vettore delle 10.000 distanze, e calcolare la posizione della distanza rispetto a tutte le altre (usando la funzione rank()).
Calcolo della posizione
posizione <- rank (c(chi_quadro1,distanza_oss))[1]
p_value1 <- 1-posizione/length(distanza_oss)
p_value1
## [1] 0.9449
Non rifiuto dell'ipotesi nulla
In questo caso, dunque, non rifiutiamo l'ipotesi nulla, in quanto
- il valore del $\chi^{2}$ è inferiore al valore critico con $\alpha = 0.05$: 0.0900735 < 5.88
- il p-value è pari a 0.9449, superiore ad $\alpha = 0.05$.
Calcolo su nuovo campione
Proviamo ora a fare la stessa verifica, ma questa volta partendo da un campione da noi generato, in cui le frequenze osservate sono pari a 22, 35, 43.
È facile intuire che questa distribuzione è molto diversa da quella attesa: 17 51 33.
Il calcolo inferenziale, però, ci permette una stima più precisa.
R: Calcolo su nuovo campione
t_campione2 <- c (22, 35,43)
chi_quadro2 <- sum(((t_campione2-t_atteso)^2)/t_atteso)
chi_quadro2
## [1] 10.27145
posizione <- rank (c(chi_quadro2,distanza_oss))[1]
p_value2 <- 1-posizione/length(distanza_oss)
p_value2
## [1] 0.00545
Rifiuto dell'ipotesi nulla
In questo caso, dunque, rifiutiamo l'ipotesi nulla, in quanto
- il valore del $\chi^{2}$ è superiore ad entrambi i valori critici: con $\alpha = 0.01$: 10.2714461 > 9.47
- conseguentemente, il p-value (pari a 0.00545), è inferiore ad $\alpha = 0.01$.
La distribuzione $\chi^{2}$
Nella sezione precedente, abbiamo calcolato valori critici e p-value basandoci sulla distribuzione dell'errore di campionamento generata dalla simulazione.
La distribuzione che abbiamo ottenuto, applicando la formula del $\chi^{2}$, è una distribuzione nota con il nome distribuzione $\chi^{2}$.
L'ambiente R mette a disposizione delle funzioni che, similmente alle distribuzioni normale e t di Student, permette di calcolare alcuni valori legati alla distribuzione $\chi^{2}$.
$\chi^{2}$: funzioni in R
Con rchisq è possibile generare dei valori casuali, con distribuzione $\chi^{2}$. Con dchisq possiamo ottenere la densità della distribuzione per un determinato valore.
Così come per la distribuzione t di Student, anche la $\chi^{2}$ è una famiglia di distribuzioni, che differiscono fra loro in base ai gradi di libertà (df, degree of freedom).
Pertanto, le funzioni legate al $\chi^{2}$ si attendono, fra gli argomenti, anche i gradi di libertà.
Gradi di libertà
Nell'esempio precedente, la tabella delle distribuzioni aveva un rango pari a 3 (ovvero, avevamo 3 categorie: agricoltura, industria, terziario).
In una tabella unidimensionale, i gradi di libertà sono pari a r-1. Nel nostro caso, dunque df=2.
Nel prossimo grafico, visualizziamo nuovamente l'istogramma delle distribuzioni della distanza dai campioni generati alle frequenze attese (figura \ref). All'istogramma sovrapponiamo la distribuzione $\chi^{2}$ con 2 gradi di libertà.
Sovrapposizione fra distribuzione osservata e teorica
plot_range <- seq(0,15,by=.25)
prob_dist <- dchisq(plot_range,2)
hist(distanza_oss,freq=FALSE,breaks=20)
lines(plot_range,prob_dist,type="l", col=2)
Utilizzo della distribuzione
Come vediamo, la distribuzione della nostra simulazione si sovrappone quasi perfettamente alla distribuzione $\chi^{2}$. Appurata questa sovrapposizione, possiamo sfruttare la distribuzione $\chi^{2}$ per calcolare valori critici e p-value.
Ad esempio, grazie alla funzione pchisq possiamo calcolare il p-value dei due campioni che abbiamo utilizzato negli esperimenti precedenti.
R: Calcolo del p-value usando pchisq
la funzione pchisq, analogamente a pnorm, ci permette di calcolare l'area della distribuzione a destra di un determinato valore. Con 1 - pchisq calcoliamo l'area rimanente, che corrisponde al p-value.
1-pchisq(chi_quadro1, df=2)
## [1] 0.9559623
1-pchisq(chi_quadro2, df=2)
## [1] 0.005882796
Approssimazione dei risultati
Come possiamo notare, i p-value non sono identici a quelli ottenuti partendo dalla simulazione, in quanto la sovrapposizione fra la distribuzione ottenuta dalla simulazione non è perfettamente identica alla distribuzione teorica.
Ciononostante, i valori sono molto simili, e ci portano a trarre le stesse conclusioni: non rifiuto di $H_0$ nel primo caso, rifiuto nel secondo.
La funzione chisq.test \index{Funzioni R!chisq.test}
Il nostro excursus ha, naturalmente, una finalità esclusivamente didattica. In pratica, per calcolare il test del $\chi^{2}$, R ci mette a disposizione la funzione chisq.test, che ci restituisce il calcolo del $\chi^{2}$, i gradi di libertà, e il p-value.
Applichiamo la funzione al primo gruppo.
chisq.test (x=t_campione,p=prob_attesa)
##
## Chi-squared test for given probabilities
##
## data: t_campione
## X-squared = 0.090074, df = 2, p-value = 0.956
La funzione ci restituisce il calcolo del $\chi^{2}$, i gradi di libertà (df), il p-value.
chisq.test, secondo gruppo
Applichiamo la stessa funzione al secondo gruppo.
chisq.test (x=t_campione2,p=prob_attesa)
##
## Chi-squared test for given probabilities
##
## data: t_campione2
## X-squared = 10.271, df = 2, p-value = 0.005883
Confronto fra due variabili nominali
Statistica bivariata
Nella sezione precedente abbiamo analizzato il caso di del confronto fra una distribuzione osservata ed una attesa.
La procedura, però, può essere utilizzata anche per valutare delle ipotesi relative al rapporto che intercorre fra due variabili nominali, ovvero nel contesto di una statistica bivariata.
Esempio: categoria lavorativa e genere
Per introdurre questa statistica, modifichiamo l'esempio precedente. Immaginiamo di voler capire se, nella distribuzione della forza lavoro fra agricoltura, industria e terziario, vi sono differenze di genere.
Per fare questo, raccogliamo un campione di persone attive, e per ognuno di loro identifichiamo il genere e la categoria lavorativa (agricoltura, industria, terziario).
Il processo
- A partire da questi dati, possiamo creare una tabella di contingenza a doppia entrata, di dimensione r * c, dove r è pari al rango della prima variabile, e c a quello della seconda.
- In secondo luogo, calcoliamo la tabella delle frequenze attese.
- Calcoliamo, attraverso la formula \ref, la distanza fra le frequenze attese e quelle osservate.
- Calcoliamo il p-value, attraverso la funzione $1 - pchisq(\chi^{2}, df = (r-1)(c-1))$.
Calcolare le frequenze attese
L'unica novità di rilievo, rispetto all'algoritmo relativo ad una sola variabile, è il calcolo delle frequenze attese è più complicato (ma non troppo)
Nuovamente, le frequenze attese si basano sull'ipotesi nulla, ovvero che non vi sia alcun legame fra le due variabili.
In termini di probabilità condizionale, si assume che la probabilità che un individuo appartenga ad una delle categorie della seconda variabile non cambi a seconda che l'individuo appartenga ad una categoria della prima, e viceversa.
L'ipotesi di indipendenza
Nel nostro esempio, l'ipotesi di indipendenza assume che il fatto di essere maschio (o femmina) non influisca sulla probabilità di essere occupato nell'agricoltura, nell'industria o nel terziario, e viceversa: il fatto di essere impiegato nel terziario non incide sulla probabilità di essere maschio o femmina.
\todo{Formalizzare l'ipotesi di indipendenza}
In base a questa assunzione, la frequenza attesa, nella categoria agricoltura, maschio è pari alla probabilità associata alla categoria agricoltura, moltiplicata per la probabilità associata alla categoria maschio, moltiplicata per la numerosità del campione.
Frequenze attese
Continuiamo ad assumere probabilità pari a .17, .51 e .32 per la variabile tipo di occupazione, e assumiamo che la popolazione attiva sia per il 54% maschile e per il 46% femminile (di nuovo, sono percentuali inventate). Decidendo per un campione di 100 persone, la frequenza attesa, per la casella agricoltore maschio, sarà pari a .17 * .54 * 100 = 9.18
\begin f_{e[i,j]} = \frac{f_ f_} \end \begin f_{e[i,j]} = p_ p_ n \end
R: Generare il data-frame
Creiamo un data.frame con due colonne: il genere e l'occupazione. In primo luogo creiamo un'urna, con probabilità .54 e .46, da cui estrarre il genere. Poi, usando l'urna precedente (relativa all'occupazione), creiamo un vettore di 100 osservazioni. Creiamo un secondo vettore, nuovamente di 100 osservazioni, relative al genere. Creiamo infine il data.frame con i due vettori, campione_genere e campione_occupazione.
genere_atteso <- c(rep(1,54),rep(2,46))
length(genere_atteso)
## [1] 100
t_genere_atteso <- table(genere_atteso)
prob_genere_attesa <- t_genere_atteso/length(genere_atteso)
campione_occupazione <- sample(atteso,100,replace=FALSE)
campione_genere <- sample(genere_atteso,100,replace=FALSE)
campione <- data.frame (campione_genere,campione_occupazione)
t_campione <- table(campione)
R: calcolo delle probabilità attese
Calcoliamo le probabilità attese. Calcoliamo la somma marginale delle righe e delle colonne. Calcoliamo le probabilità attese, pari a (marginale riga/numerosità campione) * (marginale colonna/numerosità campione).
marginali_riga <- apply(t_campione,1,sum)
marginali_colonna <- apply(t_campione,2,sum)
prob_attese <-(marginali_riga/100) %*% t(marginali_colonna/100)
Frequenze attese e osservate
Le frequenze attese sono pari alle probabilità attese, moltiplicate per il numero di osservazioni.
t_atteso <- prob_attese * 100
t_atteso
## 1 2 3
## [1,] 9.18 27.54 17.28
## [2,] 7.82 23.46 14.72
t_campione
## campione_occupazione
## campione_genere 1 2 3
## 1 8 27 19
## 2 9 24 13
R: calcolo di $\chi^{2}$ e p-value
Calcoliamo il valore della statistica $\chi^{2}$. Calcoliamo poi il p-value, utilizzando la funzione pchisq.\index{Funzioni R!pchisq}
chi_quadro <- sum(((t_campione-t_atteso)^2)/t_atteso)
p_value <- 1-pchisq(chi_quadro, df=2)
chi_quadro
## [1] 0.7249337
p_value
## [1] 0.6959574
R: uso di chisq.test
Lo stesso calcolo può essere eseguito -- più agevolmente -- usando la funzione chisq.test.
(chisq_test1 <- chisq.test (t_campione))
##
## Pearson's Chi-squared test
##
## data: t_campione
## X-squared = 0.72493, df = 2, p-value = 0.696
Leggere l'output
La funzione chisq.test ci restituisce il nome del test: Pearson's Chi-squared test; il valore della statistica: X-squared = 0.7249337; i gradi di libertà: 2; il p-value = 0.6959574.
Come prevedibile -- considerata la modalità con cui abbiamo generato il campione -- dal calcolo del $\chi^{2}$ non possiamo rifiutare l'ipotesi nulla, in quanto $p-value= 0.6959574 > \alpha=0.05$.
Le frequenze osservate nel data frame generato scegliendo le due variabili in maniera indipendente non si discostano significativamente dalle frequenze attese.