class: center, middle, inverse, title-slide # L’approccio simulativo ### Stefano Bussolon ### 13 aprile 2019 --- ## Esempio: l'altezza media [Our World in Data](https://ourworldindata.org/human-height) ci dice che l'altezza media della popolazione mondiale maschile è di 178.4 centimetri, quella femminile di 164.7. Dai grafici del documento si evince che la deviazione standard è di 7.6 e 7.1cm, rispettivamente. ```r mediaM <- 178.4 dsM <- 7.6 mediaF <- 164.7 dsF <- 7.1 ``` --- ## Altezza: popolazioni Usiamo questi parametri per creare una *popolazione virtuale* di 100.000 osservazioni. ```r nPop <- 100000 # creo una *popolazione* di 100.000 *maschi* altezzaPopM <- rnorm(nPop,mediaM,dsM) # creo una *popolazione* di 100.000 *femmine* altezzaPopF <- rnorm(nPop,mediaF,dsF) mean(altezzaPopM) ``` ``` ## [1] 178.3992 ``` ```r mean(altezzaPopF) ``` ``` ## [1] 164.704 ``` Come prevedibile, la media delle due popolazioni si aggira intorno a 178.4 (178.3992402) e 164.7 (164.7040009). --- ## Estrarre un campione Immaginiamo ora di non conoscere i parametri (media e sd dell'altezza) e di volerli stimare attraverso la misura di un campione della popolazione. Estraiamo un campione di 50 maschi e uno di 50 femmine e calcoliamo la media. ```r # estraggo un campione di 50 maschi campioneM50 <- sample(altezzaPopM,50) mean(campioneM50) ``` ``` ## [1] 179.0803 ``` ```r campioneF50 <- sample(altezzaPopF,50) mean(campioneF50) ``` ``` ## [1] 164.3394 ``` --- ## Errore di campionamento La media dei due campioni si avvicina alla media delle due popolazioni, ma con uno scostamento. ```r mean(campioneM50)-mean(altezzaPopM) ``` ``` ## [1] 0.6810409 ``` ```r mean(campioneF50)-mean(altezzaPopF) ``` ``` ## [1] -0.3645841 ``` Questo scostamento è definito [errore di campionamento](https://en.wikipedia.org/wiki/Sampling_error). Estraendo un nuovo campione, otterremo medie diverse. ```r mean(sample(altezzaPopM,50)) ``` ``` ## [1] 176.7765 ``` ```r mean(sample(altezzaPopF,50)) ``` ``` ## [1] 165.1337 ``` --- ## Estrarre molti campioni Nelle slide precedenti abbiamo estratto due campioni di maschi e due campioni di femmine, e abbiamo visto che le medie sono diverse. Cosa succederebbe se provassimo ad estrarre mille campioni di maschi (o di femmine)? Inoltre, abbiamo deciso di estrarre campioni di 50 osservazioni. Cosa succederebbe se diminuissimo (o aumentassimo) la numerosità dei campioni? Per scoprirlo, possiamo fare delle simulazioni. --- ```r nCampioni <- 1000 mCampioniM20 <- rep(0,nCampioni) for (i in 1:length(mCampioniM20)) { # estraiamo 1000 campioni di numerosità 20 # e calcoliamo la media mCampioniM20[i] <- mean(sample(altezzaPopM,20)) } mCampioniM100 <- rep(0,nCampioni) for (i in 1:length(mCampioniM100)) { # estraiamo 1000 campioni di numerosità 100 # e calcoliamo la media mCampioniM100[i] <- mean(sample(altezzaPopM,100)) } ``` --- `mCampioniM20` è il vettore delle medie dei 1000 campioni da 20 osservazioni che abbiamo generato. `mCampioniM100` quello dei campioni da 100 osservazioni. La media delle medie dei campioni da 20 è pari a 178.4142246; la media delle medie dei campioni da 100 è 178.4142246. A causa dell'errore di campionamento, le 1000 medie saranno diverse fra loro. ```r head(mCampioniM20) ``` ``` ## [1] 178.4310 179.7994 177.8794 179.1827 179.9245 178.2687 ``` Se disegno l'istogramma (o il diagramma di densità) di tutte le medie, ottengo la *distribuzione dell'errore di campionamento*. --- ## La distribuzione dell'errore di campionamento ```r library(ggplot2) # creo un data.frame con il fattore numerosità e i valori medi numerosità <- gl (2,1000,labels=c("20","100")) erroreCampionario <- data.frame(numerosità, medie= c(mCampioniM20,mCampioniM100)) # creo il density plot, distinto per i due vettori # (campioni 20 e da 100) ggDensityErrCamp <- ggplot(data=erroreCampionario, aes(medie)) + geom_density(aes(fill=factor(numerosità)), alpha=0.5) + labs(title="Errore campionario", x="altezza (cm)", y="", subtitle="distribuzione dell'errore", caption="distribuzione dell'errore campionario") ``` --- ```r ggDensityErrCamp ``` <!-- --> --- ## Differenza fra due campioni ```r (altezzaUominiBiondi <- mean(sample(altezzaPopM,50))) ``` ``` ## [1] 178.8027 ``` ```r (altezzaUominiMori <- mean(sample(altezzaPopM,50))) ``` ``` ## [1] 177.8789 ``` Riprendendo l'esempio dei due campioni. Noi sappiamo che la differenza fra le due medie è dovuta all'errore di campionamento. Ma ammettiamo di non saperlo. Ammettiamo che i due campioni appartengano il primo a uomini biondi, il secondo a uomini mori. Il fatto che le due medie non siano uguali ci porta a dire che l'altezza media della popolazione dei biondi è diversa da quella dei mori? Il fine della statistica inferenziale (bivariata) è quello di rispondere a questo tipo di domande. --- 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. ```r 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` ```r # l'approccio parametrico ttest_1 <- t.test(varIntervalliA, varIntervalliB) ttest_1 ``` ``` ## ## Welch Two Sample t-test ## ## data: varIntervalliA and varIntervalliB ## t = 4.0202, df = 597.66, p-value = 6.559e-05 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 1.371784 3.992188 ## sample estimates: ## mean of x mean of y ## 101.36213 98.68015 ``` --- Utilizziamo ora il test di Wilcoxon ```r # 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 = 53698, p-value = 4.192e-05 ## alternative hypothesis: true location shift is not equal to 0 ``` --- ### La simulazione 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. ```r # 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à ```r 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). ```r differenzaOsservata ``` ``` ## [1] 2.681986 ``` ```r confidenza_simulazione ``` ``` ## 2.5% 97.5% ## -1.282811 1.333749 ``` ```r p_value_simulazione ``` ``` ## [1] -9.999e-05 ``` --- 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. --- ```r 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 sript 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. ```r source("https://www.bussolon.it/didattica/r_intro/rmd/fSimulazione.r") ``` --- ### Primo esempio ```r varIntervalliA <- rnorm(osservazioni/2,(100+1),40) varIntervalliB <- rnorm(osservazioni/2,(100-1),40) risultato <- confronto2gruppi(varIntervalliA, varIntervalliB) risultato$parametrico$p.value ``` ``` ## [1] 0.9093201 ``` ```r risultato$nonParametrico$p.value ``` ``` ## [1] 0.8644267 ``` ```r risultato$simP_value ``` ``` ## [1] 0.5364635 ``` --- ### Secondo esempio ```r varIntervalliA <- rnorm(osservazioni/2,(100+4),8) varIntervalliB <- rnorm(osservazioni/2,(100-4),8) risultato <- confronto2gruppi(varIntervalliA, varIntervalliB) risultato$parametrico$p.value ``` ``` ## [1] 4.081605e-39 ``` ```r risultato$nonParametrico$p.value ``` ``` ## [1] 3.359765e-35 ``` ```r risultato$simP_value ``` ``` ## [1] -0.000999001 ``` --- ## 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`. --- ### Primo esempio 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. ```r risultato <- correlazione(varIntervalliA, varIntervalliB) risultato$parametrico$p.value ``` ``` ## [1] 0.08202336 ``` ```r risultato$nonParametrico$p.value ``` ``` ## [1] 0.04232428 ``` ```r risultato$simP_value ``` ``` ## [1] 0.03596404 ``` --- ### Secondo esempio Costruiamo ora la variabile `varIntervalliC`, che in questo caso si basa su `varIntervalliA`, e dunque ci aspettiamo che la correlazione sia significativa. ```r varIntervalliC <- varIntervalliA + rnorm(length(varIntervalliA), 0, 3) risultato <- correlazione(varIntervalliA, varIntervalliC) risultato$parametrico$p.value ``` ``` ## [1] 3.181404e-138 ``` ```r risultato$nonParametrico$p.value ``` ``` ## [1] 0 ``` ```r risultato$simP_value ``` ``` ## [1] -0.000999001 ``` --- ### Esercizio Creare il grafico della distribuzione delle permutazioni della correlazione, con le linee che rappresentano l'intervallo di confidenza e la correlazione osservata.