---
title: "L'approccio simulativo"
author: "Stefano Bussolon"
date: "13 aprile 2019"
nooutput: html_document
output:
  xaringan::moon_reader:
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## 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 parametri}
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 popolazioni}
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)
mean(altezzaPopF)

```

Come prevedibile, la media delle due popolazioni si aggira intorno a `r mediaM` (`r mean(altezzaPopM)`) e `r mediaF` (`r mean(altezzaPopF)`).

---

## 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 campione}
# estraggo un campione di 50 maschi
campioneM50 <- sample(altezzaPopM,50)
mean(campioneM50)
campioneF50 <- sample(altezzaPopF,50)
mean(campioneF50)

```
---

## Errore di campionamento

La media dei due campioni si avvicina alla media delle due popolazioni, ma con uno scostamento.

```{r scostamento}
mean(campioneM50)-mean(altezzaPopM)
mean(campioneF50)-mean(altezzaPopF)
```

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))
mean(sample(altezzaPopF,50))
```

---

## 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 `r mean(mCampioniM20)`; la media delle medie dei campioni da 100 è `r mean(mCampioniM20)`.  
A causa dell'errore di campionamento, le 1000 medie saranno diverse fra loro. 

```{r headMedie}
head(mCampioniM20)
```

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 distErr}
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 due campioni}
(altezzaUominiBiondi <- mean(sample(altezzaPopM,50)))
(altezzaUominiMori <- mean(sample(altezzaPopM,50)))
```

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

```

---

Utilizziamo ora il test di Wilcoxon

```{r}
# l'approccio non parametrico *classico*
wilTest_1 <- wilcox.test(varIntervalliA, varIntervalliB)
wilTest_1
```

---

### 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 simula_simulaTTest}
# 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 simula_simulaTTest_2}
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
confidenza_simulazione
p_value_simulazione
```

---

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
risultato$nonParametrico$p.value
risultato$simP_value
```

---

### 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
risultato$nonParametrico$p.value
risultato$simP_value
```

---

## 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
risultato$nonParametrico$p.value
risultato$simP_value

```

---

### 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
risultato$nonParametrico$p.value
risultato$simP_value



```

---

### Esercizio

Creare il grafico della distribuzione delle permutazioni della correlazione, con le linee che rappresentano l'intervallo di confidenza e la correlazione osservata.
