class: center, middle, inverse, title-slide # Introduzione alla statistica Inferenziale con R ### Stefano Bussolon ### 18 marzo 2019 --- 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. <table class="table table-striped" style="margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:left;"> </th> <th style="text-align:left;"> dipendente nominale </th> <th style="text-align:left;"> dipendente numerica </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> indipendente nominale </td> <td style="text-align:left;"> chi quadro </td> <td style="text-align:left;"> t-test, ANOVA </td> </tr> <tr> <td style="text-align:left;"> indipendente numerica </td> <td style="text-align:left;"> analisi discriminante, regressione logistica </td> <td style="text-align:left;"> correlazione, regressione </td> </tr> </tbody> </table> 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)) ``` ``` ## varNominale2 ## varNominale1 C D ## A 128 232 ## B 82 158 ``` --- ```r plot(tabella_NN) ``` <!-- --> --- ```r marginali_riga <- apply(tabella_NN,1,sum) marginali_colonna <- apply(tabella_NN,2,sum) (prob_attese <-(marginali_riga/osservazioni) %*% t(marginali_colonna/osservazioni)) ``` ``` ## C D ## [1,] 0.21 0.39 ## [2,] 0.14 0.26 ``` ```r (frequenze_attese <- prob_attese * osservazioni) ``` ``` ## C D ## [1,] 126 234 ## [2,] 84 156 ``` ```r tabella_NN - frequenze_attese ``` ``` ## varNominale2 ## varNominale1 C D ## A 2 -2 ## B -2 2 ``` --- ## Il calcolo Utilizziamo il `\(\chi^2\)` di Pearson ```r (temp <- ((tabella_NN - frequenze_attese)^2)/frequenze_attese) ``` ``` ## varNominale2 ## varNominale1 C D ## A 0.03174603 0.01709402 ## B 0.04761905 0.02564103 ``` ```r (chi_quadro <- sum(temp)) ``` ``` ## [1] 0.1221001 ``` --- ## 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 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 ``` ``` ## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: varNominale1 and varNominale2 ## X-squared = 0.068681, df = 1, p-value = 0.7933 ``` La funzione ci dice che ha applicato la statistica *Pearson's Chi-squared test with Yates' continuity correction*. Che i gradi di libertà sono 1 ((2-1)*(2-1)). Il valore della statistica è 0.0686813 ed il p-value è 0.793. 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 ``` ``` ## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: varNominale1 and varNominale3 ## X-squared = 28.061, df = 1, p-value = 1.175e-07 ``` Anche in questo caso la statistica applicata è la Pearson's Chi-squared test with Yates' continuity correction, ed anche in questo caso i gradi di libertà sono 1 ((2-1)*(2-1)). Il valore della statistica in questo caso è 28.0613688 ed il p-value è 1.18e-07. 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 ``` ``` ## ## Shapiro-Wilk normality test ## ## data: varIntervalli1 ## W = 0.99702, p-value = 0.344 ``` La funzione ci dice che ha applicato la statistica Shapiro-Wilk normality test. Il valore della statistica è 0.997024 ed il p-value è 0.344. 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 ``` ``` ## ## One-sample Kolmogorov-Smirnov test ## ## data: varIntervalli1 ## D = 0.027046, p-value = 0.7724 ## alternative hypothesis: two-sided ``` La funzione ci dice che ha applicato la statistica One-sample Kolmogorov-Smirnov test. Il valore della statistica è 0.0270456 ed il p-value è 0.772. 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 ``` ``` ## ## Pearson's product-moment correlation ## ## data: varIntervalli1 and varIntervalli2 ## t = 1.2082, df = 598, p-value = 0.2274 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## -0.03081874 0.12888266 ## sample estimates: ## cor ## 0.04934737 ``` La funzione ci dice che ha applicato la statistica Pearson's product-moment correlation. Che i gradi di libertà sono 598 (n-2). Il valore della statistica è 1.2082145 ed il p-value è 0.227. La correlazione è pari a 0.0493474. 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 ``` ``` ## ## Pearson's product-moment correlation ## ## data: varIntervalli1 and varIntervalli3 ## t = 53.54, df = 598, p-value < 2.2e-16 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.8947110 0.9224905 ## sample estimates: ## cor ## 0.9096121 ``` In questo caso il valore della statistica è 53.5400834 ed il p-value è <2e-16. La correlazione è pari a 0.9096121. 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 # install.packages("lmtest") library(lmtest) ``` ``` ## Caricamento del pacchetto richiesto: zoo ``` ``` ## ## Caricamento pacchetto: 'zoo' ``` ``` ## I seguenti oggetti sono mascherati da 'package:base': ## ## as.Date, as.Date.numeric ``` ```r htest_1 <- harvtest(varIntervalli3 ~ varIntervalli1, order.by = ~ varIntervalli1) htest_1 ``` ``` ## ## Harvey-Collier test ## ## data: varIntervalli3 ~ varIntervalli1 ## HC = 0.46326, df = 597, p-value = 0.6433 ``` La funzione ci dice che ha applicato la statistica Harvey-Collier test. Il valore della statistica è 0.4632605 ed il p-value è 0.643. 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 ``` ``` ## ## Harvey-Collier test ## ## data: varIntervalli5 ~ varIntervalli4 ## HC = 14.138, df = 597, p-value < 2.2e-16 ``` Il valore della statistica è 14.1379485 ed il p-value è <2e-16; 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 ``` ``` ## ## Spearman's rank correlation rho ## ## data: varIntervalli4 and varIntervalli5 ## S = 834212, p-value < 2.2e-16 ## alternative hypothesis: true rho is not equal to 0 ## sample estimates: ## rho ## 0.9768274 ``` La funzione ci dice che ha applicato la statistica Spearman's rank correlation rho. Il valore della statistica è 8.34212\times 10^{5} ed il p-value è <2e-16. La correlazione è pari a 0.9768274. --- ## 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) ``` ``` ## [1] 1.179897 ``` --- ### 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 ``` ``` ## ## Bartlett test of homogeneity of variances ## ## data: varIntervalli1 by varNominale4 ## Bartlett's K-squared = 1.9795, df = 1, p-value = 0.1594 ``` La funzione ci dice che ha applicato la statistica Bartlett test of homogeneity of variances. Il p-value è 0.159. --- #### 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 ``` ``` ## ## Welch Two Sample t-test ## ## data: intervalli1A and intervalli1B ## t = 0.47351, df = 482.77, p-value = 0.6361 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -2.001797 3.272931 ## sample estimates: ## mean of x mean of y ## 100.8809 100.2453 ``` La funzione ci dice che ha applicato la statistica Welch Two Sample t-test. Che i gradi di libertà sono 482.7736812. Il valore della statistica è 0.4735102 ed il p-value è 0.6360633. 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 ``` ``` ## ## Welch Two Sample t-test ## ## data: varIntervalli5[varNominale4 == "A"] and varIntervalli5[varNominale4 == "B"] ## t = 8.734, df = 503.03, p-value < 2.2e-16 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 8.139698 12.864571 ## sample estimates: ## mean of x mean of y ## 104.79119 94.28906 ``` In questo caso il valore della statistica è 8.7339715 ed il p-value è 3.656523\times 10^{-17}, 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) ``` ``` ## [1] 4.359855 ``` --- 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 ``` ``` ## ## Bartlett test of homogeneity of variances ## ## data: varIntervalli6 by varNominale4 ## Bartlett's K-squared = 148.99, df = 1, p-value < 2.2e-16 ``` Cosa ci dice il test? -- La funzione ci dice che il p-value è <2e-16. --- ### 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 ``` ``` ## ## Wilcoxon rank sum test with continuity correction ## ## data: varIntervalli6[1:(osservazioni/2)] and varIntervalli6[(osservazioni/2 + 1):osservazioni] ## W = 11466, p-value < 2.2e-16 ## alternative hypothesis: true location shift is not equal to 0 ``` La funzione ci dice che ha applicato la statistica Wilcoxon rank sum test with continuity correction. Il valore della statistica è 1.1466\times 10^{4} ed il p-value è <2e-16. 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) ``` ``` ## Df Sum Sq Mean Sq F value Pr(>F) ## varFattore1 2 330 165 0.657 0.519 ## Residuals 597 149823 251 ``` ```r # 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=0.6574382\)`; infine, calcola il p-value: `\(p= 0.52\)`. --- ## 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) ``` ``` ## Df Sum Sq Mean Sq F value Pr(>F) ## varFattore1 2 7285 3643 14.52 6.99e-07 *** ## Residuals 597 149823 251 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r m_aov_2 <- matrix(unlist(summary(aov_2)),ncol = 5) ``` In questo caso, poiché `varIntervalli7` è stata costruita su `varFattore1`, la statistica `\(F=14.5152177\)` ed il p-value: `\(p= <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.test(aov_2$residuals) ``` ``` ## ## Shapiro-Wilk normality test ## ## data: aov_2$residuals ## W = 0.99755, p-value = 0.5278 ``` --- Per testare l'ipotesi di omoschedasticità, si può usare il test di Bartlett: ```r bartlett.test(varIntervalli7~varFattore1) ``` ``` ## ## Bartlett test of homogeneity of variances ## ## data: varIntervalli7 by varFattore1 ## Bartlett's K-squared = 8.5404, df = 2, p-value = 0.01398 ``` --- ## 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)) ``` ``` ## Tukey multiple comparisons of means ## 95% family-wise confidence level ## factor levels have been ordered ## ## Fit: aov(formula = varIntervalli7 ~ varFattore1) ## ## $varFattore1 ## diff lwr upr p adj ## Z-X 2.001363 -1.720758 5.723484 0.4165472 ## Y-X 8.186566 4.464445 11.908687 0.0000010 ## Y-Z 6.185203 2.463082 9.907324 0.0003096 ``` --- 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 8.1865656, l'intervallo di confidenza va da un minimo di 4.4644446 ad un massimo di 11.9086866. *p adj* è il p-value aggiustato, che nel secondo confronto è pari a <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 (kt_1 <- kruskal.test(varIntervalli7~varFattore1)) ``` ``` ## ## Kruskal-Wallis rank sum test ## ## data: varIntervalli7 by varFattore1 ## Kruskal-Wallis chi-squared = 29.173, df = 2, p-value = 4.626e-07 ``` ### Leggere i risultati La funzione restituisce il metodo, Kruskal-Wallis rank sum test; la statistica: 29.1726093; i gradi di libertà: 2; il p-value = <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) ``` ``` ## Df Sum Sq Mean Sq F value Pr(>F) ## varFattore1 2 7285 3643 14.696 5.89e-07 *** ## varFattore2 1 3098 3098 12.499 0.000439 *** ## varFattore1:varFattore2 2 2560 1280 5.164 0.005977 ** ## Residuals 594 147232 248 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r 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 <0.001 per il primo fattore, <0.001 per il secondo, 0.006 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