Il fine della statistica inferenziale bivariata è di applicare delle appropriate statistiche, e di stimare se il risultato di tali statistiche sia attribuibile al caso o ad una relazione fra le due variabili.
Nell'approccio parametrico si confronta la statistica con la distribuzione teorica di riferimento.
Nell'approccio non parametrico classico si trasformano i punteggi numerici in misure ordinali.
Nell'approccio simulativo vengono utilizzate le osservazioni per simulare la distribuzione dell'errore di campionamento, e si confronta la statistica con la distribuzione emersa dalla simulazione.
Confronto della media di due gruppi
Partiamo dalla situazione in cui la variabile indipendente sia categoriale con due livelli, e la dipendente sia numerica. Nell'approccio parametrico si usa il t test, in quello non parametrico il test di Wilcoxon.
La domanda sottostante al test è: la differenza fra le medie dei due gruppi va attribuita al caso oppure può essere attribuita all'effetto della variabile indipendente?
Semplificando, non rifiutare l'ipotesi nulla significa implicitamente assumere che i due gruppi appartengano alla stessa popolazione, e che le differenze siano attribuibili al caso.
Ma come possiamo simulare questo scenario? Semplicemente, mescolando i due gruppi.
Se si mettono le osservazioni dei due gruppi in una stessa urna, e si generano due nuovi gruppi estraendo gli elementi dall'urna, la differenza fra la media di questi due gruppi sarà - per definizione - da attribuire al caso. Questo processo è una permutazione.
Se ripetiamo questa operazione molte volte (ad esempio, mille volte) e salviamo in un vettore le differenze fra i due vettori, otterremo la distribuzione delle differenze dovute al caso (la distribuzione delle permutazioni.
A questo punto è possibile confrontare la differenza dei gruppi originali con la distribuzione ottenuta. Se la differenza originale si colloca all'esterno delle code della distribuzione (definendo preventivamente l'ipotesi alternativa e alpha) possiamo rifiutare l'ipotesi nulla, altrimenti no.
In questo caso il p-value corrisponde alla posizione della differenza osservata rispetto alle differnze emerse dalla simulazione.
La simulazione
Generiamo due gruppi di osservazioni, con media di poco diversa.
osservazioni <- 600
varIntervalliA <- rnorm(osservazioni/2,(100+1),8)
varIntervalliB <- rnorm(osservazioni/2,(100-1),8)
Calcoliamo la significatività della differenza utilizzando il t.test
# l'approccio parametrico
ttest_1 <- t.test(varIntervalliA, varIntervalliB)
ttest_1
##
## Welch Two Sample t-test
##
## data: varIntervalliA and varIntervalliB
## t = 2.4949, df = 597.12, p-value = 0.01287
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.3403928 2.8587703
## sample estimates:
## mean of x mean of y
## 101.08516 99.48558
Utilizziamo ora il test di Wilcoxon
# l'approccio non parametrico *classico*
wilTest_1 <- wilcox.test(varIntervalliA, varIntervalliB)
wilTest_1
##
## Wilcoxon rank sum test with continuity correction
##
## data: varIntervalliA and varIntervalliB
## W = 49803, p-value = 0.0237
## alternative hypothesis: true location shift is not equal to 0
Ora mostriamo l'approccio simulativo. Calcoliamo la differenza osservata fra le medie; poi mescoliamo i 2 gruppi per 10000 volte, e memorizziamo le differenze fra le coppie di gruppi così generati.
# l'approccio simulativo
differenzaOsservata <- mean(varIntervalliA)-mean(varIntervalliB)
cicli <- 10000
differenze <- rep(0,cicli)
for(c in 1:cicli) {
mischiato <- sample(c(varIntervalliA,varIntervalliB))
mischiatoA <- mischiato[1:(osservazioni/2)]
mischiatoB <- mischiato[(osservazioni/2+1):osservazioni]
differenza <- mean(mischiatoA)-mean(mischiatoB)
differenze[c] <- differenza
}
Ora, partendo dalla distribuzione delle permutazioni, calcoliamo gli intervalli di confidenza.
Inoltre, sempre usando la distribuzione delle permutazioni, calcoliamo la posizione della differenza osservata rispetto alle permutazioni. La posizione, divisa per la numerosità della distribuzione, costituisce il p-value. Se la differenza è positiva, il calcolo è 1 - posizione/numerosità
confidenza_simulazione <- quantile(differenze, probs = c(0.025, 0.975))
p_value_simulazione <- 1-(rank(c(differenzaOsservata,differenze))[1]+1)/(length(differenze)+1)
Visualizziamo i risultati: l'intervallo di confidenza della simulazione, ed il p-value (calcolando la posizione della differenza osservata rispetto alla distribuzione delle permutazioni).
differenzaOsservata
## [1] 1.599582
confidenza_simulazione
## 2.5% 97.5%
## -1.283872 1.248990
p_value_simulazione
## [1] 0.0069993
Disegnamo ora il grafico della distribuzione delle permutazioni, la posizione dell'intervallo di confidenza, e la differenza. Se la linea della differenza osservata (in verde) cade all'esterno delle due linee rosse, rifiutiamo l'ipotesi nulla.
plot(density(differenze))
abline(v = confidenza_simulazione[1], col = 2, lty = 2)
abline(v = confidenza_simulazione[2], col = 2, lty = 2)
abline(v = differenzaOsservata, col = 3, lty = 2)
Abbiamo trasformato gli script appena mostrati nella funzione confronto2gruppi che abbiamo salvato nel file fSimulazione.r. Ora carichiamo la funzione e la usiamo per confrontare due gruppi utilizzando gli approcci visti in precedenza.
source("../dispensa/fSimulazione.r")
Primo esempio
varIntervalliA <- rnorm(osservazioni/2,(100+1),10)
varIntervalliB <- rnorm(osservazioni/2,(100-1),10)
risultato <- confronto2gruppi(varIntervalliA, varIntervalliB)
risultato$parametrico$p.value
## [1] 0.01558901
risultato$nonParametrico$p.value
## [1] 0.01692803
risultato$simP_value
## [1] 0.00789921
Secondo esempio
varIntervalliA <- rnorm(osservazioni/2,(100+4),8)
varIntervalliB <- rnorm(osservazioni/2,(100-4),8)
risultato <- confronto2gruppi(varIntervalliA, varIntervalliB)
risultato$parametrico$p.value
## [1] 1.872977e-33
risultato$nonParametrico$p.value
## [1] 1.479295e-30
risultato$simP_value
## [1] -9.999e-05
Correlazione fra due variabili numeriche
L'approccio simulativo attraverso la permutazione può essere usato anche per valutare la significatività della correlazione fra due variabili numeriche.
In questo caso si permuta l'ordine delle osservazioni di una delle due variabili. In questo modo si annulla l'eventuale relazione di una variabile con l'altra. La simulazione consiste nel ripetere l'operazione (permutazione dell'ordine di una variabile) per molte volte (ad esempio 10.000) ed ogni volta calcolare l'indice di correlazione, che verrà salvato nel vettore simulazioni. Il vettore costituisce la distribuzione delle permutazioni su cui calcolare l'intervallo di confidenza ed il p-value della correlazione osservata, nella stessa maniera esposta nella parte dedicata al confronto fra due gruppi.
Lo script della simulazione è stato salvato - assieme all'uso del test parametrico (Pearson) e non parametrico classico (Spearman) nella funzione correlazione.
Nel primo esempio usiamo ancora i vettori varIntervalliA e varIntervalliB. Poiche sono stati costruiti in maniera indipendente, ci aspettiamo che la correlazione non sia significativa.
risultato <- correlazione(varIntervalliA, varIntervalliB)
risultato$parametrico$p.value
## [1] 0.9452157
risultato$nonParametrico$p.value
## [1] 0.6826092
risultato$simP_value
## [1] 0.4645535
Costruiamo ora la variabile varIntervalliC, che in questo caso si basa su varIntervalliA, e dunque ci aspettiamo che la correlazione sia significativa.
varIntervalliC <- varIntervalliA + rnorm(osservazioni/2, 0, 3)
risultato <- correlazione(varIntervalliA, varIntervalliC)
risultato$parametrico$p.value
## [1] 1.891757e-133
risultato$nonParametrico$p.value
## [1] 0
risultato$simP_value
## [1] -9.999e-05
Esercizio
Creare il grafico della distribuzione delle permutazioni della correlazione, con le linee che rappresentano l'intervallo di confidenza e la correlazione osservata.