---
title: "Introduzione alla statistica Inferenziale con R"
author: "Stefano Bussolon"
date: "18 marzo 2019"
output:
  xaringan::moon_reader:
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false

---

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



Il fine dell'analisi inferenziale è di fare delle inferenze su di una popolazione a partire dalle osservazioni di un campione.

Il fine dell'analisi inferenziale univariata è di stimare il valore di un parametro della popolazione a partire da una statistica calcolata sul campione.  
Il fine dell'analisi inferenziale bivariata è quello di stimare la significatività di una relazione fra due variabili.  
Le analisi multivariate sono concettualmente un'estensione dell'analisi bivariata.

---

## Le analisi bivariate

### Tipi di variabili e statistica

La tipologia di statistica inferenziale da applicare si basa sulla tipologia di variabili. Possiamo distinguere fra variabili categoriali, ordinali, ad intervalli e a rapporti.

Queste quattro tipologie possono essere raggruppate in variabili nominali (categoriali e, generalmente, ordinali) e variabili numeriche (a intervalli, a rapporti).

---

La tipologia di statistica che può essere applicata si basa sulla tipologia delle variabili indipendenti e dipendenti.

```{r tabella_stat, echo=FALSE}
library(kableExtra)
tabella_statistiche <- matrix(c(
																"chi quadro",
																"t-test, ANOVA",
																"analisi discriminante, regressione logistica",
																"correlazione, regressione"
																), 
															ncol = 2, byrow=TRUE)
colnames(tabella_statistiche)<- c("dipendente nominale", "dipendente numerica")
rownames(tabella_statistiche)<- c("indipendente nominale", "indipendente numerica")

kable(tabella_statistiche) %>%
  kable_styling("striped")
```

La distinzione fra variabile indipendente e dipendente è legata al disegno sperimentale.

---

## La logica

La logica sottostante è quella di identificare una statistica che possa misurare la relazione fra le due (o più) variabili.

* nel caso di due variabili numeriche la misura è la correlazione
* nel caso di due variabili categoriali si misura la differenza fra le frequenze osservate e le frequenze attese
* nel caso del confronto fra due gruppi la misura è la differenza fra le medie
* nel caso del confronto fra più gruppi si confronta la varianza attribuita al gruppo con la varianza interna ai gruppi

---

## Distribuzione dell'errore

Come abbiamo visto nell'approccio simulativo, è necessario stimare se la misura (correlazione, differenza ...) va attribuita all'errore di campionamento (H0) oppure no.

La distribuzione dell'errore di campionamento, negli scenari citati, è riconducibile a delle distribuzioni teoriche:

* nella correlazione (metodo Pearson) la distribuzione t di student
* nella differenza della media di due campioni, di nuovo la t di student
* nel caso di due variabili categoriali, è la distribuzione $\chi^{2}$
* nel caso dell'analisi della varianza, la distribuzione F di Fisher-Snedecor

---

# Due variabili nominali

---

## Il processo

* si crea 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.
* si calcola la tabella delle frequenze attese.
* si calcola la *distanza* fra le frequenze attese e quelle osservate.
* si calcola il p-value

---

## Le frequenze attese

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).

---

## Un esempio

Per comprendere il processo, introduciamo un esempio: generiamo due variabili categoriali indipendenti, visualizziamo il grafico delle variabili, e calcoliamo la statistica.

---

```{r }
osservazioni <- 600 # è più facile fare i calcoli
varNominale1 <- sample(c(rep("A",osservazioni*.6),rep("B",osservazioni*.4)))
varNominale2 <- sample(c(rep("C",osservazioni*.35),rep("D",osservazioni*.65)))
#df_Nominale <- data.frame(varNominale1,varNominale2)
(tabella_NN <- table(varNominale1,varNominale2))

```
---

```{r }
plot(tabella_NN)

```

---

```{r chi2--prob-attese-1}
marginali_riga <- apply(tabella_NN,1,sum)
marginali_colonna <- apply(tabella_NN,2,sum)
(prob_attese <-(marginali_riga/osservazioni) %*% t(marginali_colonna/osservazioni))
(frequenze_attese <- prob_attese * osservazioni)
tabella_NN - frequenze_attese
```


---

## Il calcolo

Utilizziamo il $\chi^2$ di Pearson

```{r}

(temp <- ((tabella_NN - frequenze_attese)^2)/frequenze_attese)
(chi_quadro <- sum(temp))

```

---

## La simulazione

```{r}

distanza_oss <-vector(mode = "numeric", length = 10000)

for (ciclo in 1:length(distanza_oss)) {
	tabella_NN <- table(varNominale1,sample(varNominale2))
	chi_quadro <- sum(((tabella_NN - frequenze_attese)^2)/frequenze_attese)
	distanza_oss [ciclo] <- chi_quadro
}
```

---

```{r chi2--hist_distanza_oss}
plot_range <- seq(0,15,by=.10)
prob_dist <- dchisq(plot_range,1)
plot(density(distanza_oss))
lines(plot_range,prob_dist,type="l", col=2)


```


---
## La funzione chisq.test

In `R`, per calcolare il test di $\chi^{2}$ si usa la funzione `chisq.test`.

[The Chi-square test of independence](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3900058/) 

---

```{r}

chi2_a <- chisq.test(varNominale1,varNominale2)
chi2_a
```

La funzione ci dice che ha applicato la statistica *`r chi2_a$method`*. Che i gradi di libertà sono `r chi2_a$parameter` ((2-1)*(2-1)). Il valore della statistica è `r chi2_a$statistic` ed il p-value è `r format.pval(chi2_a$p.value, 3)`. Come prevedibile il p-value è superiore a `0.05`, e dunque non si può rifiutare l'ipotesi nulla di indipendenza fra le due variabili.

---

Nell'esempio sucessivo, creiamo `varNominale3`, che invece non è indipendente da `varNominale1`. Disegnamo il grafico e calcoliamo la statistica.

---

```{r}
isNominale1A <- as.integer(varNominale1=="A")
# genero una variabile nominale che "dipende" da varNominale1
varNominale3 <- rbinom(length(varNominale1),1,.35+isNominale1A*.2)
plot(table(varNominale1,varNominale3))

```
---

```{r}
chi2_b <- chisq.test(varNominale1,varNominale3)
chi2_b

```

Anche in questo caso la statistica applicata è la `r chi2_b$method`, ed anche in questo caso i gradi di libertà sono `r chi2_b$parameter` ((2-1)*(2-1)). Il valore della statistica in questo caso è `r chi2_b$statistic` ed il p-value è `r format.pval(chi2_b$p.value, 3) `. In questo caso il p-value è inferiore a `0.05`, e dunque dobbiamo rifiutare l'ipotesi nulla di indipendenza fra le due variabili.

---

# Due variabili numeriche

---

## Il test di correlazione



In caso di due variabili numeriche, si utilizza la funzione `cor.test(x,y)`.

Anche in questo caso, generiamo due variabili (numeriche) indipendenti e creiamo lo scatterplot.
Valutiamo se le variabili superano il test di normalità, e calcoliamo il test di correlazione.

```{r}
#osservazioni <- 300
varIntervalli1 <- rnorm(osservazioni, mean=100, sd=15)
varIntervalli2 <- rnorm(osservazioni, mean=30, sd=5)

```

---

```{r}

plot(varIntervalli1,varIntervalli2)
lineare <- lm(varIntervalli2 ~ varIntervalli1)
abline(lineare)
```
---

### Normalità delle variabili numeriche

Poiché le statistiche inferenziali parametriche assumono una distribuzione delle osservazioni di tipo normale, è generalmente opportuno valutare la distribuzione osservata di una variabile non soltanto attraverso metodi grafici e descrittivi, ma anche attraverso dei test di normalità.  
Utilizzeremo due di questi test:

* Il test di *Kolmogorov-Smirnov* permette di confrontare due distribuzioni arbitrarie, e può essere usato per il confronto fra la distribuzione osservata e la distribuzione normale;
* Il test di normalità *Shapiro-Wilk* è finalizzato a valutare la normalità della distribuzione osservata.

---

Le due misure possono dare risultati differenti. Risulta pertanto necessario un processo di valutazione che tenga conto sia dei risultati dei test che dell’analisi grafica della distribuzione.

Questa regola pratica vale in ogni ambito della ricerca e dell’analisi dei dati: la metodologia ci indica delle procedure che è opportuno seguire, per minimizzare il rischio di errori che mettano a repentaglio affidabilità e validità della ricerca.

Le procedure, però, non vanno seguite pedissequamente. Conoscere i princîpi e gli assunti dell’analisi dei dati ci permette di fare delle inferenze ragionevolmente robuste anche nei casi, e sono molti, in cui non è possibile una applicazione meccanica della procedura.

---

### Il processo

In genere si segue un processo che prevede:

* la visualizzazione del grafico di dispersione e la linea di regressione
* la valutazione della normalità, sia graficamente che attraverso gli opportuni test
* la valutazione grafica della linearità dell'eventuale correlazione
* il test di linearità
* se gli assunti di normalità o di linearità sono violati, l'utilizzo del test inferenziale non parametrico; altrimenti, del test inferenziale parametrico

---

### QQnorm e qqline

le funzioni `qqnorm()` e `qqline()` ci permettono di visualizzare graficamente la normalità di una distribuzione.

---

```{r}

qqnorm(varIntervalli1)
# disegno una linea blu, di spessore 2
qqline(varIntervalli1, col = "steelblue", lwd = 2)

```

---

Se la distribuzione delle osservazioni nel grafico è lineare e sovrapponibile alla linea, si può assumere che la distribuzione sia normale.

---

### Test di Shapiro-Wilk

Applichiamo il test di Shapiro-Wilk per valutare la normalità della distribuzione `varIntervalli1`.

```{r}

stest_1 <- shapiro.test(varIntervalli1)
stest_1

```
La funzione ci dice che ha applicato la statistica `r stest_1$method`. Il valore della statistica è `r stest_1$statistic` ed il p-value è `r format.pval(stest_1$p.value, 3) `. Il p-value è superiore a `0.05`, e dunque non si rifiuta l'ipotesi nulla di normalità della variabile.

---

### Test di Kolmogorov-Smirnov

Confrontiamo la distribuzione di `varIntervalli1` con la distribuzione normale teorica, attraverso la funzione `ks.test`.

```{r}

kstest_1 <- ks.test(varIntervalli1, "pnorm", 
						mean = mean(varIntervalli1), 
						sd=sd(varIntervalli1))
kstest_1

```

La funzione ci dice che ha applicato la statistica `r kstest_1$method`. Il valore della statistica è `r kstest_1$statistic` ed il p-value è `r format.pval(kstest_1$p.value, 3) `. Il p-value è superiore a `0.05`, e dunque non si può rifiutare l'ipotesi nulla di normalità della variabile.

---

#### Esercizio

testare la normalità della variabile `varIntervalli2`.

---

### Correlazione: test parametrico

Per valutare la correlazione fra due variabili numeriche si utilizza la funzione `cor.test()`. La funzione permette di applicare tre metodi diversi: Pearson (di default), Kendall e Spearman. In caso di non violazione degli assunti (normalità, linearità) si usa il metodo parametrico, di Pearson.

---

```{r}

# non è necessario specificare method="pearson"
corTest_1 <- cor.test(varIntervalli1,varIntervalli2)
corTest_1

```

La funzione ci dice che ha applicato la statistica `r corTest_1$method`. Che i gradi di libertà sono `r corTest_1$parameter` (n-2). Il valore della statistica è `r corTest_1$statistic` ed il p-value è `r format.pval(corTest_1$p.value, 3) `. La correlazione è pari a `r corTest_1$estimate`. Come prevedibile il p-value è superiore a `0.05`, e dunque non si può rifiutare l'ipotesi nulla di indipendenza fra le due variabili.

---

Nel prossimo esempio, generiamo `varIntervalli3`, una variabile ad intervalli *dipendente* da `varIntervalli1`

```{r}

varIntervalli3 <- varIntervalli1 + rnorm(100, mean=0, sd=8)
plot(varIntervalli1,varIntervalli3)
lineare_1 <- lm(varIntervalli3 ~ varIntervalli1)
abline(lineare_1)
```
---
```{r}

corTest_2 <- cor.test(varIntervalli1,varIntervalli3)
corTest_2

```

In questo caso il valore della statistica è `r corTest_2$statistic` ed il p-value è `r format.pval(corTest_2$p.value, 3) `. La correlazione è pari a `r corTest_2$estimate`. Come prevedibile il p-value è inferiore a `0.05`, e dunque si rifiuta l'ipotesi nulla di indipendenza fra le due variabili.

---

## Testare la linearità

Il test parametrico assume che la correlazione sia lineare. 

La non linearità può essere diagnosticata sia attraverso la visualizzazione del grafico di dispersione delle due variabili che il grafico di dispersione dei residui sui valori attesi.

Da un punto di vista inferenziale, è possibile applicare la statistica *Harvey-Collier*, che valuta la linearità della correlazione: `harvtest` (va caricata la libreria `lmtest`).

---

```{r warning=FALSE}
# install.packages("lmtest")
library(lmtest)
htest_1 <- harvtest(varIntervalli3 ~ varIntervalli1, order.by = ~ varIntervalli1)
htest_1

```

La funzione ci dice che ha applicato la statistica `r htest_1$method`. Il valore della statistica è `r htest_1$statistic` ed il p-value è `r format.pval(htest_1$p.value, 3) `. Il p-value è superiore a `0.05`, e dunque non si può rifiutare l'ipotesi nulla di linearità della regressione.

---

```{r}
plot(lineare_1$fitted.values,lineare_1$residuals)
```

---

### Esempio di correlazione non lineare

```{r}

varIntervalli4 <- rnorm (osservazioni,0,4)
varIntervalli5 <- (1 / (1 + exp(-varIntervalli4)))*10+rnorm(osservazioni,0,0.3)
varIntervalli4 <- varIntervalli4+8+rnorm(osservazioni,0,0.3)
```

---

```{r}


plot(varIntervalli4,varIntervalli5)
lineare_2 <- lm(varIntervalli5 ~ varIntervalli4)
abline(lineare_2)
# disegno una curva che segue i punti
lines(smooth.spline(varIntervalli4,varIntervalli5), col='red', lwd=2)

```

---

```{r}

htest_2 <- harvtest(varIntervalli5 ~ varIntervalli4, order.by = ~ varIntervalli4)
htest_2


```

Il valore della statistica è `r htest_2$statistic` ed il p-value è `r format.pval(htest_2$p.value, 3) `; essendo inferiore a `0.05`, dobbiamo rifiutare l'ipotesi nulla di linearità della regressione.

---

La non linearità appare evidente anche dalla curva `smooth.spline` e dal grafico di dispersione dei residui.

```{r}

plot(lineare_2$fitted.values,lineare_2$residuals)

```
---

## Coefficiente di Spearman


Nelle circostanze in cui la relazione fra le due variabili non sia lineare, ma tenda ad essere comunque monotona, è possibile utilizzare il modello non-parametrico della correlazione: $\rho$ di Spearman.

In questo modello, il calcolo della relazione si effettua non sui valori delle due variabili, ma sulla loro posizione ordinale.

Questa statistica, pertanto, può essere applicata anche nella circostanza in cui una o entrambe le variabili siano di tipo ordinale.

---

```{r}
# specifico il metodo: Spearman
corTest_3 <- cor.test (varIntervalli4,
											 varIntervalli5, 
											 method='spearman')
corTest_3

```
La funzione ci dice che ha applicato la statistica `r corTest_3$method`. Il valore della statistica è `r corTest_3$statistic` ed il p-value è `r format.pval(corTest_3$p.value, 3) `. La correlazione è pari a `r corTest_3$estimate`. 

---

## Indipendente categoriale, dipendente numerica

In caso di indipendente categoriale (a 2 livelli) e dipendente numerica si usa la funzione `t.test`.

Anche in questo caso, prima di optare per il test parametrico, è necessario valutare la normalità delle variabili dipendenti. In secondo luogo, è necessario valutare l'omogeneità delle varianze delle due variabili.

Nel primo esempio, applichiamo il test alla variabile dipendente `varNominale1` e alla indipendente `varIntervalli1`.

`t.test` si aspetta i due gruppi di osservazione, che possono essere estratti attraverso un filtro. In alternativa si può usare la sintassi `t.test(dipendente ~ indipendente)`

---

```{r}

varNominale4 <- sample(c(rep("A",osservazioni*.6),rep("B",osservazioni*.4)))


boxplot(varIntervalli1~varNominale4)

```
---

```{r}
# per comodità, creiamo le due variabili distinte
intervalli1A <- varIntervalli1[varNominale4=="A"]
intervalli1B <- varIntervalli1[varNominale4=="B"]
# calcoliamo la varianza
var_1A <- var(intervalli1A)
var_1B <- var(intervalli1B)
# calcolo il rapporto fra la maggiore e la minore
max(var_1A,var_1B)/min(var_1A,var_1B)


```

---

### Testare l'omogeneità della varianza (omoschedasticità)

Uno dei test per valutare che le varianze non siano significativamente diverse è il test di Bartlett, attraverso la funzione `bartlett.test`. Se il test risulta non significativo (p-value > 0.05) non si rifiuta l'ipotesi nulla di omoschedasticità.

```{r}
bartlett_1 <- bartlett.test(varIntervalli1 ~ varNominale4)
bartlett_1

```
La funzione ci dice che ha applicato la statistica `r bartlett_1$method`. Il p-value è `r format.pval(bartlett_1$p.value, 3) `.

---

#### Esercizio

Valutare la normalità di intervalli1A e intervalli1B

---

Se gli assunti non sono violati, si applica il Test di Student con la funzione `t.test()`.

```{r}

ttest_1 <- t.test(intervalli1A,intervalli1B)
# sintassi alternativa: t.test(varIntervalli1 ~ varNominale4)
ttest_1									

```

La funzione ci dice che ha applicato la statistica `r ttest_1$method`. Che i gradi di libertà sono `r ttest_1$parameter`. Il valore della statistica è `r ttest_1$statistic` ed il p-value è `r ttest_1$p.value `. Come prevedibile il p-value è superiore a `0.05`, e dunque non si può rifiutare l'ipotesi nulla di indipendenza fra le due variabili.

---

Nel prossimo esempio, generiamo `varIntervalli5`, una variabile ad intervalli *dipendente* da `varNominale4`

```{r}
isNominale4A <- as.integer(varNominale4=="A")
varIntervalli5 <- rnorm(osservazioni, mean=95, sd=15)+isNominale4A*10
```
---

```{r}

boxplot(varIntervalli5~varNominale4)
```
---

```{r}

ttest_2 <- t.test(varIntervalli5[varNominale4=="A"],varIntervalli5[varNominale4=="B"])
ttest_2									

```

In questo caso il valore della statistica è `r ttest_2$statistic` ed il p-value è `r ttest_2$p.value `, inferiore a `0.05`, e dunque si deve rifiutare l'ipotesi nulla di indipendenza fra le due variabili.

---

Confrontiamo ora due gruppi di osservazioni con varianza diversa.

```{r}

varNominale4 <- c(rep("A",osservazioni/2),rep("B",osservazioni/2))
varIntervalli6 <- c(rnorm(osservazioni/2,10,2),rnorm(osservazioni/2,15,4))
var_6A <- var(varIntervalli6[1:(osservazioni/2)])
var_6B <- var(varIntervalli6[(osservazioni/2+1):osservazioni])
# calcolo il rapporto fra la maggiore e la minore
max(var_6A,var_6B)/min(var_6A,var_6B)
```
---

Confrontiamo ora due gruppi di osservazioni con varianza diversa.

```{r}

boxplot(varIntervalli6 ~ varNominale4)
```
---

Applichiamo il test di Bartlett per verificare l'omoschedasticità delle due variabili.

```{r}

bartlett_2 <- bartlett.test(varIntervalli6 ~ varNominale4)
bartlett_2

```

Cosa ci dice il test?

--

La funzione ci dice che il p-value è `r format.pval(bartlett_2$p.value, 3) `.


---

### Test non parametrico: test di Wilcoxon

In caso di violazione degli assunti di normalità o di omoschedasticità, si può usare il test non parametrico di Wilcoxon, con la funzione `wilcox.test`.

```{r}
wilcox_1 <- wilcox.test(varIntervalli6[1:(osservazioni/2)],
			varIntervalli6[(osservazioni/2+1):osservazioni])
wilcox_1
```

La funzione ci dice che ha applicato la statistica `r wilcox_1$method`. Il valore della statistica è `r wilcox_1$statistic` ed il p-value è `r format.pval(corTest_3$p.value, 3) `.
Essendo il p-value inferiore a 0.05, rifiuto l'ipotesi nulla.

---



# L'analisi della varianza

Nelle circostanze in cui la variabile dipendente abbia più di due livelli, oppure nel caso di più di una variabile indipendente, è necessario applicare l'Analisi della Varianza (anova).

Anche in questo caso, prima di optare per il test parametrico, è necessario valutare sia la normalità dei residui della variabile indipendente che l'omogeneità delle varianze delle due variabili.

---

## La funzione aov

R mette a disposizione, per il calcolo dell'analisi della varianza, la funzione `aov(y~x)`, dove y è la variabile dipendente, numerica, e x è il fattore.

Per mostrare l'uso di aov creiamo una variabile nominale (`varFattore1`) con 3 livelli (`X`,`Y`,`Z`) e calcoliamo l'analisi della varianza utilizzando la variabile ad intervalli `varIntervalli1`. Poiché le due variabili sono del tutto indipendenti, ci aspettiamo che l'ipotesi nulla non sia falsificata.

---

```{r}
varFattore1 <- factor (c(rep('X',osservazioni/3),
												 rep('Y',osservazioni/3),
												 rep('Z',osservazioni/3))) 

aov_1 <- aov(varIntervalli1~varFattore1)
summary(aov_1)
# la variabile m_aov_1 mi serve per riportare i risultati
m_aov_1 <- matrix(unlist(summary(aov_1)),ncol = 5)
```


Utilizzando la funzione `summary` su `aov` è possibile avere il dettaglio dei risultati dell'analisi. Nel caso di una analisi ad una via, avremo una tabella con due righe. La seconda riga calcola i gradi di libertà, la somma dei quadrati, e la media dei quadrati dei residui. La prima riga calcola i gradi di libertà, la somma dei quadrati, e la media dei quadrati del modello; inoltre, calcola la statistica $F=`r m_aov_1[1,4]`$; infine, calcola il p-value: $p= `r format.pval(m_aov_1[1,5], 2, 0.001)`$.

---

## Secondo esempio

Creiamo la variabile numerica `varIntervalli7`, che non è indipendente da `varFattore1`. In questo caso ci aspettiamo che l'analisi della varianza risulti significativa.

---

```{r}
varIntervalli7 <- varIntervalli1
# rendiamo varIntervalli7 non indipendente da varFattore1
varIntervalli7[varFattore1=="X"]<- varIntervalli7[varFattore1=="X"]-3
varIntervalli7[varFattore1=="Y"]<- varIntervalli7[varFattore1=="Y"]+7
aov_2 <- aov(varIntervalli7~varFattore1)
summary(aov_2)
m_aov_2 <- matrix(unlist(summary(aov_2)),ncol = 5)

```

In questo caso, poiché `varIntervalli7` è stata costruita su `varFattore1`, la statistica $F=`r m_aov_2[1,4]`$ ed il p-value: $p= `r format.pval(m_aov_2[1,5], 2, 0.001)`$; l'ipotesi nulla va dunque rifiutata.

---

## Verifica degli assunti

Come ogni approccio parametrico, anche l'analisi della varianza fa delle assunti:

* indipendenza delle osservazioni
* distribuzione normale dei *residui*
* omoschedasticità: la varianza dell'errore è costante
* gli errori sono fra loro indipendenti

---

## Distribuzione dei *residui*

Si assume che i residui abbiano una distribuzione normale, con media pari a 0, e varianza costante fra i gruppi. Per testare l'ipotesi di normalità, è possibile usare il test di Shapiro-Wilk sui residui del modello dell'analisi della varianza:

```{r shapiro}
shapiro.test(aov_2$residuals)
```

---

Per testare l'ipotesi di omoschedasticità, si può usare il test di Bartlett:

```{r bartlett}
bartlett.test(varIntervalli7~varFattore1)
```

---

## Confronti multipli

---

### Confronti multipli ed errore

L'analisi della varianza ci permette di verificare se le differenze fra le medie di tre o più campioni sono da attribuire all'errore campionario, o se sono significative.

Una volta rifiutata l'ipotesi nulla, però, resta da determinare quali differenze sono significative. L'analisi della varianza, infatti, ci dice se vi è almeno una coppia di gruppi la cui differenza è significativa, ma non ci dice quali differenze lo sono.

---

Per poter determinare quali differenze sono significative, diventa necessario confrontare i gruppi a due a due.

Si potrebbe decidere di utilizzare, per confrontare a due a due i diversi gruppi, il t-test. Ma applicare ripetutamente il t-test aumenta la probabilità di incorrere in un errore del primo tipo.

Diventa dunque necessario adottare dei test di confronti multipli capaci di mantenere sotto controllo l'errore del I tipo.

---

### Il test di Tukey

Attraverso il test di Tukey è possibile mantenere l'errore di tipo I entro un predeterminato valore di $\alpha$ (generalmente pari a 0.05).

Il test di Tukey permette di *correggere* il p-value in base al numero di confronti che vengono effettuati nel confronto multiplo, senza però penalizzare eccessivamente la statistica.

---

### La funzione R TukeyHSD

La funzione di R per il calcolo del confronto con il metodo Tukey è `TukeyHSD`.
La funzione si applica sul risultato della corrispondente analisi della varianza.

```{r}
(confronti1 <- TukeyHSD(aov_2, ordered = TRUE))

```

---

La funzione ritorna una tabella, con una riga per ogni confronto, dove vengono mostrate:

* la coppia confrontata (es, il confronto fra il gruppo Y ed il gruppo X); l'ordine è tale che il gruppo con media più alta è davanti all'altro;
* la differenza (positiva) fra i due gruppi;
* l'intervallo di confidenza della differenza; ad esempio, nel secondo confronto (Y-X), la differenza è di `r confronti1$varFattore1[2,1]`, l'intervallo di confidenza va da un minimo di `r confronti1$varFattore1[2,2]` ad un massimo di `r confronti1$varFattore1[2,3]`. *p adj* è il p-value aggiustato, che nel secondo confronto è pari a `r format.pval(confronti1$varFattore1[2,4], 2, 0.001)`.


---

## Test non parametrico

Vi sono circostanze in cui l'analisi della varianza non può essere applicata, in quanto vengono meno alcuni assunti o condizioni:

* non si può assumere la normalità della distribuzione degli errori
* il numero di osservazioni per ogni gruppo è minore di 10
* la variabile dipendente non è ad intervalli, ma ordinale

In questi casi è possibile applicare il test non parametrico di Kruskal-Wallis

---

### R: la funzione kruskal.test

Applichiamo il test di *Kruskal-Wallis* al nostro secondo esempio.

```{r kruskal}
(kt_1 <- kruskal.test(varIntervalli7~varFattore1))
```


### Leggere i risultati

La funzione restituisce il metodo, `r kt_1$method`; la statistica: `r kt_1$statistic`; i gradi di libertà: `r kt_1$parameter`; il p-value = `r  format.pval(kt_1$p.value, 2, 0.001)`.

---

## Anova a due vie

### Due variabili indipendenti

L'analisi della varianza che abbiamo introdotto, può essere estesa anche ai casi in cui le variabili indipendenti sono più di una.

Nell'analisi della varianza a due vie, si indaga la relazione fra due variabili indipendenti, entrambe categoriali, ed una variabile dipendente, quantitativa.

In questa sezione analizziamo la circostanza in cui le variabili indipendenti sono due, ma la logica rimane la stessa anche nelle circostanze in cui le variabili indipendenti sono più di due.

---

Creiamo la variabile `varFattore2` con due livelli.

```{r}
varFattore2 <- factor (rep(c(rep('E',osservazioni/6),
														 rep('F',osservazioni/6)),3))
# modifichiamo varIntervalli7 in modo da creare un effetto di varFattore2
varIntervalli7[varFattore2=="E"] <- 
	varIntervalli7[varFattore2=="E"]+5
```

---
Rappresentiamo l'effetto delle due variabili indipendenti sulla variabile dipendente usando la funzione `interaction.plot`.

```{r}

interaction.plot (varFattore1,varFattore2,varIntervalli7)
```

---

Calcoliamo l'analisi della varianza, utilizzando la funzione `aov(dipendente~indipendente1+indipendente2+indipendente1:indipendente2)`, dove `indipendente1:indipendente2` rappresenta l'interazione delle due variabili indipendenti.  
Una sintassi del tutto equivalente è  
`aov(dipendente~indipendente1*indipendente2)`.


---

```{r}
aov_3 <- aov(varIntervalli7~varFattore1+varFattore2+varFattore1:varFattore2)
summary(aov_3)
m_aov_3 <- matrix(unlist(summary(aov_3)),ncol = 5)

```

`summary(aov_3)` ha 4 righe: i residui (ultima riga), `varFattore1`, `varFattore2` e l'interazione. I p-value sono `r format.pval(m_aov_3[1,5], 2, 0.001)` per il primo fattore, `r format.pval(m_aov_3[2,5], 2, 0.001)` per il secondo, `r format.pval(m_aov_3[3,5], 2, 0.001)` per l'interazione.

---

## Misure ripetute

Nell'analisi della varianza ad una via con un disegno a misure ripetute, la sintassi è la seguente:  
`aov(dipendenteNumerica~IndipendenteCategoriale
	+Error(IdPartecipante/IndipendenteCategoriale)`

[How to do Repeated Measures ANOVAs in R | R-bloggers](https://www.r-bloggers.com/how-to-do-repeated-measures-anovas-in-r/)  


---

## Esercizi

Caricare il file <https://s3.eu-central-1.amazonaws.com/bussolon/dati/df_simulato_3.RDS> nel dataframe df_simulato_2.  
Per ognuna delle seguenti coppie di variabili, generare l'appropriata rappresentazione grafica bivariata, valutare quando opportuno gli assunti (normalità, omoschedasticità) e le appropriate statistiche inferenziali.

* nominale1 e nominale2
* nominale2 e nominale3
* nominale2 e intervalli1
* nominale1 e intervalli3
* intervalli1 e intervalli2
* intervalli1 e intervalli4
* nominale1, nominale2 e intervalli5


```{r echo=FALSE, eval=FALSE}
df_simulato_2 <- readRDS("/home/bussolon/documenti/didattica/r_dottorato/dati/df_simulato_2.RDS")

# nominale1, nominale2

plot(df_simulato_2$nominale1, df_simulato_2$nominale2)
chisq.test(df_simulato_2$nominale1, df_simulato_2$nominale2)

# nominale2, nominale3

plot(df_simulato_2$nominale2, df_simulato_2$nominale3)
chisq.test(df_simulato_2$nominale2, df_simulato_2$nominale3)

# nominale2, intervalli1

boxplot(df_simulato_2$intervalli1 ~ df_simulato_2$nominale2)

# devo farlo per entrambi i gruppi
s2Intervalli1A <- df_simulato_2$intervalli1[df_simulato_2$nominale2=="catA"]
s2Intervalli1B <- df_simulato_2$intervalli1[df_simulato_2$nominale2=="catB"]

qqnorm(s2Intervalli1A)
# disegno una linea blu, di spessore 2
qqline(s2Intervalli1A, col = "steelblue", lwd = 2)
shapiro.test(s2Intervalli1A)

qqnorm(s2Intervalli1B)
qqline(s2Intervalli1B, col = "steelblue", lwd = 2)
shapiro.test(s2Intervalli1B)

bartlett.test(df_simulato_2$intervalli1 ~ df_simulato_2$nominale2)

t.test(df_simulato_2$intervalli1 ~ df_simulato_2$nominale2)

# nominale1 e intervalli3 (sì)

boxplot(df_simulato_2$intervalli3 ~ df_simulato_2$nominale1)

# devo farlo per entrambi i gruppi
s2Intervalli3X <- df_simulato_2$intervalli3[df_simulato_2$nominale1=="catX"]
s2Intervalli3Y <- df_simulato_2$intervalli3[df_simulato_2$nominale1=="catY"]


qqnorm(s2Intervalli3X)
qqline(s2Intervalli3X, col = "steelblue", lwd = 2)
shapiro.test(s2Intervalli3X)

qqnorm(s2Intervalli3Y)
qqline(s2Intervalli3Y, col = "steelblue", lwd = 2)
shapiro.test(s2Intervalli3Y)

bartlett.test(df_simulato_2$intervalli3 ~ df_simulato_2$nominale2)


t.test(df_simulato_2$intervalli3 ~ df_simulato_2$nominale1)

# intervalli1 e intervalli2 (sì)
plot(df_simulato_2$intervalli1,df_simulato_2$intervalli2)
cor.test(df_simulato_2$intervalli1,df_simulato_2$intervalli2)

# intervalli1 e intervalli4 (sì, non lineare)

plot(df_simulato_2$intervalli1,df_simulato_2$intervalli4)
abline(lm(df_simulato_2$intervalli4 ~ df_simulato_2$intervalli1))
cor.test(df_simulato_2$intervalli1,df_simulato_2$intervalli4, method='spearman')

```

