STRUMENTI E MODELLI DI ANALISI DEI DATI NELLA RICERCA PSICOBIOLOGICA, Sbobinature di Statistica

Scarica STRUMENTI E MODELLI DI ANALISI DEI DATI NELLA RICERCA PSICOBIOLOGICA e più Sbobinature in PDF di Statistica solo su Docsity! 3° anno, Scienze e Tecniche Psicologiche CORSO DI STATISTICA 2: STRUMENTI E MODELLI DI ANALISI DEI DATI NELLA RICERCA PSICOBIOLOGICA Anno accademico 2020/2021 PRIMO SEMESTRE PROFESSOR GIUSEPPE PAGNONI VIDEOPILLOLE + LEZIONI + DOMANDE ESAME CORSO DI LAUREA INTERATENEO: UNIVERSITA’ DEGLI STUDI DI MODENA E REGGIO UNIVERSITA’ DEGLI STUDI DI PARMA STRUMENTI E MODELLI DI ANALISI DEI DATI NELLA RICERCA PSICOBIOLOGICA iuseppe.pagnoni@unimore.it GUARDARE | DATI Cosa facciamo quando abbiamo un database da analizzare? Per prima cosa occorre un modo per guardare i dati raccolti che suggerisca che cosa c'è in essi>+importante è la visualizzazione grafica dei dati: se grafichiamo i dati opportunamente possiamo avere un'idea delle caratteristiche essenziali del campione. Un BUON GRAFICO: - dllustra in maniera chiara gli aspetti importati dei dati (non tutto, le caratteristiche più importanti), - è adeguato al livello del pubblico (grafico semplice, divulgativo? pubblico inesperto; grafico rigoroso, più informativo? pubblico esperto), - adeguato al mezzo di presentazione (diapositive in conferenze+grafico semplice; articolo di rivista scientifica>grafico più dettagliato), - si ricorda facilmente (deve far ricordare l'informazione in maniera facile), - facilita la comprensione dei risultati, dei dati (è efficace se presenta il problema da una prospettiva che facilita la comprensione). Un - è difficile da interpretare, - può essere involontariamente fuorviante, - può contenere informazioni ridondanti o distraenti. - è quasi impossibile da interpretare, - è colmo di informazioni inutili (occorre essere essenziali: spiegazione semplice ma non troppo), - distorce intenzionalmente i risultati, - è inaccurato. — il 3D quando il 2D è sufficiente (la tridimensionalità spesso è ridondante e rende difficile il confronto delle altezze), — non bisogna inserire motivi non necessari nei bar plot, — non bisogna usare etichette mal definite sugli assi, — nonsi deve alterare la scala dell'asse delle y per far risaltare differenze trascurabili. Quando vediamo un grafico che non va conviene, se l'argomento e i dati sono per noi importanti, verificare che non ci siano errori. = http://junkcharts.typepad.com nttp://Junkcharts.ty0epad.com Itp: / nn mad arrattoss.ora f Ittp://www mediamatters.org # Esegui 41 grafico (istognanma per gruppi) mydata nist + geon histogram() + labetx = "hand length (an)", y = count") + facet_uzapl- gender) + theneltert = element tert (size = 16)) | # Inposta 41 grafico (densità) \npdeta dena <° ggplot(aydata, aes (x = hand.L, fill = gender)) (4 Esegui 41 grafico apdata dena + ‘geom denzity (alpha = 0.5) + | 3abe0x = “hand length (an)", y = "count*) + | theme (text = element text(ize = 16)) malo fem: 20 0.075 20 E E 0.050 ender 5 È male è formal 10 0.025 o ; ; ; ; 0.000 160 200 zio 160 200 160 200 260 hand length (mm) hand length (mm) | # Inposta dl grafico (box plot) © Maimumobservetion valle box <- Ggplociayaeo, seeta = gender, y © hand.1)) | Outlier (# Esegui il grafseo Upper fence (not drawn) | mpdaca box + (0 gain boxplot0 + Maximum observation | label = “hand length (an)") + 15408 below upper fence | Bleneltext © elezent tentieize = 14)) Upper quartile (Q3) T n 2 s9 Median È ni L î £ — Lower quartile (Q1) $ = _ tinimamssenaton 3 + 1.5 + IQR È T 2 Lower fence (not drawn) | # Imposta il grafico (bar plot) (npdata bar <- ggplot (nydata, aee(x = gender, y = band.L)} | (# Esegui 4L grafico mpdata bar + male female gender Mostrare affiancate Le medie di maschi e femmine per ciascun tipo di misura (piede, mano, statura) stot_sunsary(fun.y = mean, geom = "bar", fill = “unite”, colour = "black") + Ì | stat/sunmery(fun.data = mean_cl normal, geom = "errorbar*, Ì width = 0.2, colour = "red*) + | labefy = “hand length Cam") + | theme(text = element text(zize = 14)) 150° hand length (mm) so mile ferrate gender IL FORMATO OTTIMALE DEI DATI IN R: TIDY DATA È opportuno strutturare i data frame di R come segue: - ogni variabile è posta nella sua propria colonna, - ogni osservazione è posta nella sua propria riga, - ogni valore è posto nella sua propria cella, 1500 e 1000 È 500- Mel fodLL handL stalure body measure 000000 -000000 È 000000 000000 varianiee obrervatione | # Imposta 41 grafteo (scatter plot) \nydata ecatter <- ggplot(uydata, aes(x = foot.L, y = Ì | # Esegui £l grafico | nydata_scatter + (0 geom posnt0) + | feom'enooth(method = "ln") + Ì Ì pature)) labs(x = "foot length fam)", y = "stature (mm)"} + theme (text = element_tert(size = 14)) (# Imposta il grafico (seutter plot per gruppi) \nydata scatter <- geplot(uydata, aes(x = foot.L, y = stature, colour = gender}) Ì (# Esegui dl grafico mydata_scatter + (°° geom pointl) + | geom_snooth (method = "In") + | labslx = “foot length (am)", y = "atature (am)"*) + | theme(text = element text(size = 14)) E gender S mala 3 feiemala 8 stature (mm) Ò 250 270 190 210 240 260 toot length (mm) foot length (mm) zia Il comando facet_wrap (* x) esegue diverse sfaccettature del grafico a seconda del valore di una variabile. Il grafico della densità è un istogramma smussato, interpolato per dare una visione più morbida dell'andamento della distribuzione. Il comando fill dice che il colore che riempie le curve smussate dev'essere selezionato in base al livello della variabile-+ due curve nello stesso grafico. Il comando alpha indica la trasparenza (1: completamente opaco; 0: completamente trasparente). Nel box plot la linea centrale indica la mediana, le estremità della scatola indicano il quartile superiore (75%) e inferiore (25%); i dati anomali (outliers) sono al di là (sopra e sotto) di 1.5 * IQR; le estremità dei baffi indicano il valore più alto e più basso prima del limite precedentemente citato. Nel bar plot l'altezza della barra indica la media e vengono aggiunte le barre di errore (dispersione: deviazione standard...). Il comando stat summary indica una grandezza statistica (media, barre di errore...). Con il comando position=dodge si pongono le barre una di fianco all'altra. Con uno scatter plot si può osservare la relazione tra due variabili. Il comando geom smooth produce una linea che fitti i punti (banda grigia: banda di confidenza). Il comando color grafica i punti con un colore che dipende da una variabile. ESERCIZI PER CASA Temperatura corporea e frequenza cardiaca bodytemp_heartrate.dat 130 persone, con le seguenti variabili: - gender male, female - hr.bpm frequenza cardiaca (battiti/min) - temp.Ctemperatura corporea (<C) Importare i dati e realizzare i seguenti grafici: - box plot (maschi/femmine) per frequenza cardiaca, - barplot(maschi/femmine) per temperatura corporea, - scatter plot a gruppi (maschi/femmine) per frequenza cardiaca vs. temperatura corporea Suggerimenti - Aprite il file dati con un editor di testo (es. Notepad) per verificarne il formato 2 - Notate che: 1.i campi sono separati da spazi 2. la prima riga contiene il nome delle variabili -> importate i dati con la funzione read.delim() -> usando: header = TRUE poi sep = "" LA PROBABILITÀ DEI DATI OSSERVATI Variazione insita al campionamento -> Supponiamo di misurare una certa grandezza x su un campione s1 di N soggetti estratti a caso da una popolazione P, e di calcolare la media di tale variabile nel nostro campione, x—1 Se ripeto il processo con un nuovo campione s2 estratto ancora a caso da diverso. Supponendo di raccogliere un numero molto grande di campioni dalla stessa popolazione, posso costruire la distribuzione delle medie calcolate ep su ciascun campione, ovvero la distribuzione campionaria della media. La Distribuzione della popolazione i DISTRIBUZIONE CAMPIONARIA è la distribuzione dei valori medi misurati su tanti campioni. Il TEOREMA DEL LIMITE CENTRALE dice che se il campione è sufficientemente grande (N>30) la distribuzione campionaria Dettzione angina 7 Snorri risulta approssimativamente normale qualunque sia la distribuzione snnnai Seo 00 agmpioni =_= della popolazione, che la media della distribuzione campionaria è pari siae» 008, alla media della popolazione (1) e la deviazione standard è pari all'errore standard della media ( SE = SDINN). Significato dell'errore standard -> La media del campione osservato è in genere utilizzata come stima del valore nella popolazione: MA è affidabile? Ciò equivale a chiedere quanto grande sia la variabilità di tale misura da campione a campione. Questa informazione è espressa dalla deviazione standard della distribuzione campionaria, ovvero dall’errore standard. L'errore standard ci dice quanto sia affidabile la media del campione per stimare la media della popolazione. Stimare quanto affidabile è la media del campione equivale a dire quanto è grande la variabilità della media quando cambio campione > errore standard: deviazione standard della distribuzione campionaria. L'errore standard misura l’affidabilità della stima della media della popolazione. Intervalli di confidenza Sono un altro approccio per stabilire quanto la media campionaria sia accurata nello stimare la media della popolazione; indicano i confini di plausibilità intorno alla media osservata entro i quali pensiamo si situi la me osservata della popolazione. L’intervallo di confidenza (CI) è calcolato per quel campione specifico. Tale intervallo di confidenza (CI) è costruito in maniera tale che, se prelevassimo casualmente 100 campioni dalla stessa popolazione e calcolassimo per ognuno di essi il CI attorno alla media del campione, ci aspettiamo che la media della popolazione cada all’interno di questi CI in 95 casi su 100. Si può anche pensare il CI come quell'insieme di valori (ipotetici) della media della popolazione da cui la media del campione raccolto non differisce significativamente. IMPORTANTE: “confidenza” o “plausibilità” si riferiscono alla procedura di calcolo del CI e non ai valori specifici calcolati sul campione effettivamente raccolto. sample number Costruire degli intervalli di confidenza vuol dire dare un algoritmo per costruire un intervallo in cui posso essere ragionevolmente sicuro si situi la media della popolazione. In una distribuzione normale standard (media 0 e SD 1) il 95% centrale della distribuzione è delimitato da + 1.96. Se la distribuzione campionaria fosse normale standard, i limiti per una deviazione significativa (p = 0:05. della media di un campione da quella della popolazione (=0) sarebbero °1:96 ) l'intervallo ['1:96; +1:96] rappresenterebbe il 95% x-X CI. Trasformando la variabile iniziale Xin Z(2= SE) Z avrà una distribuzione campionaria normale standard#2; quindi la da formula per calcolare l'intervallo di confidenza al 95% sulla variabile X originale è da x+1.96 x SE, . Tutto questo per un campione sufficientemente grande, N>30. Un intervallo ve 95% di confidenza di percentuale arbitraria è dato invece da X + N,,, x SE. Poiché -1.96 e +1.96 vu ne so io n die corrispondono rispettivamente ai quantili del 2.5% e 97.5% della distribuzione normale standard z {NO:025 e N0:975), si può scrivere il 95% CI come: X + N(0,025) © SE » u » x+ N(0.975) * SE. Se abbiamo un campione più piccolo la distribuzione campionaria segue la distribuzione t di Student+intervalli di confidenza calcolati secondo questa distribuzione: XElyp x SE = intervallo di confidenza di percentuale arbitraria p% (es. 90%=9). In R, i quantili della distribuzione normale sono dati dalla funzione qnorm(p). ca INTERVALLI IMPORTANTI DELLA DISTRIBUZIONE NORMALE Intervalli importanti nella distribuzione normale sono + 1.96 (contiene il 95% dei valori), + 2.58 (contiene il 99% dei valori) e + 3.29 (contiene il 99.9% dei valori). Nei grafici gli intervalli di confidenza servono a capire se due campioni appartengono alla stessa popolazione (gli o intervalli si sovrappongono) o no (gli intervalli non si sovrappongono)*medie .; significativamente diverse o no. os INTERVALLI DI CONFIDENZA IN CAMPIONI PICCOLI Se il campione è piccolo, la distribuzione campionaria non è normale ma segue la distribuzione t di Student3. L'intervallo di confidenza al p% è dato da: Xtf,., x SE dove t(1+p)/2 è il quantile di una distribuzione t con df = N— 1. In R, i quantili della a distribuzione t sono dati dalla funzione gt (p, df). 3La distribuzione t tende a quella normale per N -> infinito. INTERVALLI DI CONFIDENZA NEI GRAFICI Consideriamo l’analisi della potenza statistica per un test di correlazione semplice, utilizzando lo stesso coefficiente r di Pearson come misura della dimensione dell'effetto In questo caso, possiamo usare la funzione pwr.r.test del pacchetto pwr: pwr.r.test(n =,r=,sig.level =, power =) dove specificando il valore di tre parametri, si ottiene il valore del quarto. Vogliamo calcolare il numero di soggetti necessario a rilevare una correlazione semplice di grandi dimensioni con una potenza statistica dell’80% e una soglia di significatività pari a 0.05: library(pwr) pwr.r.test(n =, r= .5, sig.level = .05, power = .8) approximate correlation power calculation (arctangh transformation) n=28.24841 0.5 sig.level = 0.05 power = 0.8 alternative = two.sided Calcoliamo la potenza di un test conoscendo & (lo stabiliamo noi), N e le dimensioni dell'effetto stimate dal nostro campione; se la potenza del test è > .8 diciamo che il test è sufficientemente potente per rilevare l’effetto. Si può quindi per esempio calcolare il numero di soggetti necessario per rilevare l’effetto: si guarda la letteratura precedente per le dimensioni dell'effetto che si spera di rilevare, poi & e 1 - 8 ; per un test di correlazione semplice, con una potenza statistica desiderata dell’80% alla soglia statistica di 0.05, abbiamo bisogno di circa: — servono 782 per un effetto piccolo, r=.1, 85 per un effetto medio, r=. 29 per un effetto grande, r=.5 In R la funzione pwr.r. test serve per l’analisi della potenza statistica. ESERCIZI PER CASA -> Intervalli di confidenza e analisi della potenza — Scrivete una funzione che calcoli un intervallo di confidenza di probabilità arbitraria per una variabile, e che abbia come argomenti un vettore contenente i dati del campione e la probabilità associata all'intervallo di confidenza desiderata (es. 90%); Usate la distribuzione t per il calcolo dei quantili, — Installate il pacchetto pwr e calcolate la potenza statistica di un test di correlazione semplice su 45 soggetti, con livello di significatività, alpha = :05 e dimensioni dell'effetto pari a r = 0:3 REQUISITI DEI TEST PARAMETRICI INTRODUZIONE: PRESUPPOSTI DI VALIDITA” DEI TEST La maggioranza dei test statistici sono parametrici, ovvero suppongono che i dati provengano da una distribuzione nota e caratterizzabile da parametri precisi. | presupposti/condizioni di validità dei test sono: — NORMALITÀ (a seconda del contesto, riferita alla distribuzione campionaria o alla distribuzione degli errori nel modello), — OMOGENEITÀ DELLA VARIANZA (con più gruppi, le varianze all'interno di ciascun gruppo devono essere simili; nelle correlazioni, la varianza di una variabile non deve mutare con il valore dell'altra variabile), — SCALA DI MISURA DATI A INTERVALLI (almeno; le variabili devono essere continue e tali che a intervalli eguali sulla scala di misura corrispondano differenze eguali nella proprietà misurata [ad es. ore della giornata]), — INDIPENDENZA (a seconda del contesto, si può riferire ai dati dei diversi partecipanti o, nella regressione, agli errori nel modello). Installare i seguenti pacchetti in Rstudio: * car * pastecs * psych Normalità In molti test statistici, il requisito di normalità non si riferisce alla distribuzione dei dati stessi, ma alla distribuzione campionaria (distribuzione delle medie di un campione estratto dalla popolazione). Tuttavia, sappiamo che: - seidati del campione sono distribuiti pressoché normalmente, anche la distribuzione campionaria lo sarà, - e che se il campione è sufficientemente grande (N > 30) la distribuzione campionaria sarà approssimativamente normale indipendentemente dalla forma distribuzionale della popolazione = teorema del limite centrale! Le analisi di regressione presuppongono che gli ERRORI nel modello siano distribuiti normalmente. VERIFICA VISIVA DELLA NORMALITÀ DEI DATI 1 Per verificare visivamente la normalità dei dati per prima cosa possiamo 18° graficare un istogramma dei dati con una curva normale modellata sugli stessi. Ei Per selezionare un sottoinsieme dei dati possiamo selezionare la funzione filter. Si i può usare l’aes density per rappresentare sull'asse y dell’istogramma la frequenza percentuale di soggetti che hanno un dato valore. Per rappresentare la curva normale usiamo la funzione stat_function (fun=dnorm, args=list (mean=, sd=)}; mean e sd devono essere quelle del campione e per rimuovere i dati mancanti usiamo la 360 funzione na.rm=TRUE. 365 s70 ars 360 body temperature (C) # Carica i pacchetti ggplot2 e dplyr library(ggplot2) library(dplyr) # Imposta la directory di lavoro gender br.bpa setwd(file.path( female:65 Min. 157.00 n "i male :65 st Qu.:69.00 ist Qu.:36.50 Sys.getenv("HOME"), Median 174.00 Median 136.80 "/Documents/Work/Didattica/reggioemilia/lezione_03")) Mean :T3.76 ja | i NI im (", n 3rd Qu.:79.00 # Carica i dati mydata <- read.delim("data/bodytemp_heartrate.dat", drd Qu.:79.00 header = TRUE, sep = "") summary(mydata) gender female :65 # Seleziona solo le femmine male : 0 mydata_f <- filter(mydata, gender == "female") Madian 196.50 summary(mydata_f) 3rd Qu. # Imposta il grafico mydata_f_plot <- ggplot(mydata_f, aes(temp.C)) + g eom_histogram(aes(y = ..density..), colour = "black", fill ="white") + labs(x = "body temperature (C)", y= "Density") + stat_function( fun = dnorm, args = list(mean = mean(mydata$temp.C, na.rm = TRUE), sd = sd(mydata$temp.C, na.rm = TRUE)), colour = "red", size = 1) # "Stampa" il grafico print(mydata_f_plot) METODO 2= QQ plot -> Un altro modo per verificare la normalità è il grafico quantile-quantile: si graficano i quantili del campione contro i quantili della distribuzione normale corrispondente; - lanormalità è rappresentata da una linea retta (ASIMMETRIA=0, CURTOSI=3), - mentre una forma a S indica asimmetria della distribuzione (SKEWNESS) -> verso l'alto si ha asimmetria (coda lunga) a destra, verso il basso si ha asimmetria (coda lunga) a sinistra; - una forma ad S indica una curtosi diversa da zero -> deflessioni verso l'alto o verso il basso indicano curtosi diverse da 0 - S invertita: leptocurtosi >3 [positiva: distribuzione appuntita e con code lunghe] -> distribuzione con valori più estremi della normale; —S dritta: platicurtosi <3 [negativa: distribuzione più piatta e con code corte] -> distribuzione con meno valori della normale; Nota= a volte per curtosi si intende la excess Kurtosis, ottenuta sottraendo 3 dalla prima. 3805 Per rappresentare il q-q plot usiamo la funzione ggplot{mydata, aes(sample=x)) + stat gg(). Per quantificare numericamente la normalità possiamo usare la funzione 315- kurt. 2SE/skew. 2SE > 1 indicano deviazioni significative dalla normalità; per campioni medi è bene alzare il criterio a > 1.29; per campioni grandi il criterio perde validità ed è bene osservare la forma della distribuzione e i valori stessi di curtosi e asimmetria. mydata_f_qplot <- ggplot(mydata_f, aes(sample = temp.C)) mydata_f_qplot + stat_qa() 204200 i 2 theoretical Per quantificare la normalità numericamente: library(pastecs) round(stat.desc(mydata_f$temp.C, basic = FALSE, norm = TRUE), digits = 3) » campioni piccoli: kurt.2SE o skew.2SE > 1 indicano deviazioni significative (p < 0:05) dalla normalità2 » campioni medi: è bene alzare il criterio a > 1:29 » campioni grandi (> 200): per il restringimento degli SE il criterio perde validità ed è bene valutare la forma della distribuzione e i valori stessi di curtosi e asimmetria 2Qui ci si riferisce alla excess Kurtosis 3. METODO 3= IL TEST DI NORMALITA” DI SHAPIRO -> confronta il campione con valori di una distribuzione normale con stessa media e deviazione standard. Possiamo anche usare il test di normalità di Shapiro-Wilk_attraverso la funzione shapiro.test(mydata$Sx); p < 0.05 indica non normalità; il test tende a essere significativo per N grande: in questo caso occorre interpretarlo congiuntamente all’istogramma, al q-q plot e ai valori di asimmetria e curtosi. Se una variabile include valori da gruppi diversi occorre esaminare ciascun gruppo separatamente; utilizzando la funzione by/data=mydata$x, INDICES=mydataSfattore, FUN=stat.desc, basic=FALSE, nori lichiamo un certa funzione (stat.desc) a sottogruj di dati identificati dai livelli di un fattore. Possiamo poi fare un istogramma per gruppi (facet_wrap) o un g-q plot per gruppi. shapiro.test(mydata_f$temp.C) Shapiro-Wilk normality test ETNIE data: mydata_f$temp.C median mean SE.nean CI.nean.0.95 var std.der _W=0.96354, p-value = 0.05259 36.80000000 36.E4153846 0.05137028 0.10262398 0.17152885 0.41016041 coef.var skevmess.-—skeu.23E lurtogia 1.25E normteat.W cadi nia conificati 0.01124167 0.12437827 0.20930949 1,37976169 1.17679s61 0.96354452 “> P < 0:05 indica non-normalità, tende ad essere significativo norstest.p 0.05258517 per N grande SR > interpretate il test congiuntamente a istogramma, Q-Q plot median mean SE.meen CI.mean.0.$5 var std.dev e Valori di asimmetria e curtosi. 26.70000000 136.6800000) 0.04879213 0,09747548 0.15475000 0,3933827 coef.var © ‘sktemess.skew.29E kurtosis | kurt.25E normtest.W? x 0.01072472 -0.15320934 -0.25782775 -0.59897144 -0.51086181 098561903 —\VALUTAZIONE DELLA NORMALITÀ PER GRUPPI normtept.p 0.64766103 Se una variabile include valori da gruppi diversi, occorre esaminare ciascun gruppo separatamente by(data = mydata$temp.C, INDICES = mydata$gender, FUN = stat.desc, ba: FALSE, norm = TRUE) ISTOGRAMMA PER GRUPPI fem aio mydata_plot <- ggplot(mydata, aes(temp.C)) mydata_plot + geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +1 abs(x = "body temperature (C)", y = "Density") + f fio acet_wrap(” gender) à vo LIL LIL io e io da mele no se Q-Q PLOT PER GRUPPI semaio malo mydata_qplot <- ggplot(mydata, aes(sample = temp.C)) sso si mydata_qplot + stat_qq() + facet_wrap(” gender) median nean SE.maan CI. mean. 0.95 var atd.dev 36.800 36.842 0.061 0.108 0.172 0414 cosf.var = akewnesa | slew.25E = kurtosis © kurt.25E normtest.W tina pa 4a i o.01t 0.124 0.208 1.380 LATT 0.964 thesretcal normtest.p 0.083 > fate un test di Levene e un test Fmax di Hartley per verificare l'omogeneità delle varianze sulla massa cerebrale per maschi e femmine > traete le vostre conclusioni CORRELAZIONE Con la correlazione si stabilisce se esiste una relazione tra due variabili (esercizio fisico e memoria a breve termine). È sempre utile fare uno scatter plot preliminare e i concetti fondamentali sono quelli di COVARIANZA e COEFFICIENTE DI CORRELAZIONE. Covarianza Varianza di una singola variabile -> La varianza di una variabile rappresenta la tipica distanza al quadrato dei punti di una distribuzione dalla loro media. Si dice che due variabili covariano quando deviazioni dalla media dell’una corrispondono a simili L.&-D0,-9) N-1 quando le due variabili cambiano nella stessa direzione e negativa quando cambiano in direzione opposta. Il valore della covarianza dipende dalle unità di misura scelte per le variabili>non possiamo confrontare covarianze diverse in maniera oggettiva>in generale, per confrontare variabili differenti occorre esprimerle in unità di deviazioni standard, rendendole così anche ‘adimensionali. -> possiamo fare la stessa cosa per la covarianza. deviazioni dalla media dell'altra > S? = sommatoria di (Xi-X)?/N-1= = COV (X, Y). La covarianza è positiva Coefficiente di correlazione di Pearson (R di Pearson) L,@&-DO0 D&G -3YZ,0,-7 è “r’è compresotra +1 è ‘“” indica direttamente le dimensioni dell'effetto con le seguenti soglie convenzionali: - |r.1 +.1 (effetto di piccole dimensione), - |r.! +.3 (effetto di medie dimensioni) - |r.1 +.5 (effetto di grandi dimensioni). > Valori anomali possono distorcere il valore di “r” PRESUPPOSTI PER IL TEST DI SIGNIFICATIVITÀ DI r - lanormalità delle variabili - l'omogeneità della varianza (la varianza di una variabile non deve cambiare significativamente al variare dell’altra variabile). Se una variabile è distribuita normalmente, conosciamo la probabilità associata a ciascun valore. La distribuzione campionaria di rnon è normale, ma possiamo normalizzarla con la formula di Fisher trasformando rin z, o trasformando rin tcon la formula Zr= 1/2ln(1+r) /1-r); dove Zr ha un errore standard pari a SÉ(Zr)= 1/radice(N-3). Ora, la distribuzione campionaria di Z, ha una deviazione standard pari a SE», quindi la variabile è così definita: Z=Z./SEx che ha deviazione standard pari a 1 -> basta cercare il valore di z osservato in una tabella della distribuzione normale standard per ottenere la significatività del risultato (p-value) Una ryN-2 Non è altro che la covarianza standardizzata: r= strategia alternativa è quella di convertire Lr intr= INTERVALLI DI CONFIDENZA PERR Ricordando che zr è distribuita normalmente con errore standard pari a SEzr, possiamo scrivere: 95% CI(zr) = [2r° (1:96 © SEzr) ; 2r + (1:96 * SEzr )] Poiché questi valori riguardano la variabile trasformata zr, dobbiamo infine applicare la trasformazione inversa per ottenere i corrispondenti valori 21/241, CORRELAZIONE NON SIGNIFICA CAUSALITA' ED È COMUNQUE A-DIREZIONALE È importante sottolineare che CORRELAZIONE NON SIGNIFICA CAUSALITÀ (le variazioni osservate nelle due variabili correlate possono essere determinate non dall'influenza diretta dell’una sull'altra, ma da una terza causa comune a entrambe); inoltre molte correlazioni sono spurie (altamente significative ma casuali); infine, anche nel caso di un'effettiva relazione causale tra le due variabili in esame, la correlazione non può stabilire se sia l'una a influenzare l’altra o viceversa. Quando si vogliono correlare diverse variabili a due a due si fa una matrice di correlazione. II COEFFICIENTE DI DETERMINAZIONE RZ: coei inte di correlazione al quadrato) esprime la quantità di variabilità condivisa da due variabili, ossia quanta _ variabilità dell'una sia potenzialmente è spiegabile dalla variabilità dell’altra. CALCOLARE LA CORRELAZIONE INR Pacchetti utili: = Hmisc - GGally - boot = polycor - ppcor è Single “è Maned vegan ona pasa 05% Deprension DI score) Soll-osto om (MSEl scor fun te i pe Mag pel Anp ep ut e De Funzioni per la correlazione: - cor) - cor.test() - Hmi orr() Graficare una matrice di scatterplot: - pairs() - GGally::ggpairs() GESTIONE DEI VALORI MANCANTI Raccogliendo misure di diverse variabili per molti soggetti è probabile che vi siano alcuni valori mancanti. Nel calcolo della matrice di correlazione tra tutte le variabili, un soggetto con dati mancanti può essere escluso: da tutte le correlazioni (listwise case deletion), solo da quelle correlazioni dove le due variabili in gioco contengono il dato mancante (pairwise case deletion), Nella funzione cor(), abbiamo l'opzione “us - "complete.obs" -> listwise case deletion - "pairwise.complete.obs" -> pairwise case deletion Posso infine eliminare del tutto i casi con valori mancanti: df_clean <- na.omit(df) STATURA, LUNGHEZZA MANO E LUNGHEZZA PIEDE: # Carica il pacchetto ggplot2 library(ggplot2) library(dplyr) # Imposta la directory di lavoro setwd(file.path( Sys.getenv("HOME"), stature hand.L foot.L " idatti icemili i » in. :1485 Min. :133.7 Min. :200.8 '/Documents/Work/Didattica/reggioemilia/lezione_04")) Het Qu.c1668 Jot Qu. (187.0 tav Qu. s230.8 # Carica i dati Median :1597 Median :189.6Median :233.0 mydata <- read.csv("data/stature_hand_foot_NA.csv") ie u tiene gu dioat lo a mydata_f <- filter(mydata, gender female") %>% dplyr::select(- Max. VATZI Max. :214.6 Max. :280.4 gender) summary(mydata_f) Malt SCATTER PLOT DEI DATI _ library(ggplot2) Fao mydata_f_plot <- ggplot(mydata_f, aes(x = hand.L, y = foot.L.)) Sa mydata_f_plot + geom_point() + geom_smooth(method = "Im")+ î Eno labs(x= "hand length (mm)", y = "foot length (mm)") + theme(text = element_text(size = 14)) 1 Éo is 200 hand length (mm) CORRELAZIONE BIVARIATA: COR.TEST() with(mydata_f, cor.test(hand.L, foot.L)) Pearson's product-moment correlation data: hand.L and foot.L t = 6.5254, df = 72, p-value = 8.189e-09 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.4427877 0.7356368 sample estimates: cor 0.6096085 Otteniamo: - significatività del test (p-value) - coefficiente di correlazione (cor! r - intervallo di confidenza Opzioni della funzione cor.test(): > alternative= "two-sided": test di significatività bilaterale (default), "less": test di significatività per r< 0 "greater": test di significatività per r > 0 conf.level = - 0.95 (default) ampiezza dell'intervallo di confidenza, method = - "pearson", - "spearman", "kendall", MATRICE DI CORRELAZIONE: COR () # Listwise case deletion cor(mydata_f{, c("stature", ‘hand.L", "foot.L")], use = "complete.obs") atature hand.L atature 1.0000000 0.7093046 06991001 foot.L hand.L 0.7093046 1.0000000 0. 6096085 # Pairwise case deletion cor(mydata_{f{, c("stature", “hand.L", "foot.L")], use = "pairwise.complete.obs") MATRICE DI CORRELAZIONE: HMISC::RCORR() library(Hmisc) rcorr(as.matrix(mydata_f[, c("stature", "hand.L", 'foot.L")])) Otteniamo: - coefficiente di correlazione (R) - numero casi validi per correlazione (n) - significatività (P) MATRICE DI SCATTER PLOT: METODO RAPIDO ri pairs(mydata_f[, c("stature", "hand.L", "foot.L")]) sure i ' A 7 48890 || hand al è,9 È } dale # too È n° È ' MATRICE DI SCATTER PLOT: METODO SOFISTICATO library(GGally) ggpairs(mydata_f, columns = c("stature", "hand.L", "foot.L"), lower = list(continuous = "smooth")) statura hand.L foor.L ateture nand.L fo01.L P steture hand L foot.L atature hand.L foot.L ori 1.00 0.61 1.00 duri 0.70 0.70 du6i 1.00 atature hand.L foot.L 76 na 78 n n n 78 TA 75 atature hand.L foot.L 0 Ò o o o o foot.L 0.6991001 0.6096085 1.0000000 ' ORDINARY NONPARAMETRIC BOOTSTRAP | Bootstrap Statistica : Call: Ì original bias std. error t1% 0.7104339 6.0793386-07 0.06896 boot(data = mydata_f, statistic = myboot_func, R = 2000) Vengono riportati: - il suo bias, cioè la differenza tra la media dei valori di fi “bootstrappati” e il valore originale, Imtervals : . - l'errore standard basato sui campioni di bootstrap, i ‘ to or des, 08607 ) # Intervalli di confidenza al 95% boot.cimyboot_obj) Level Percentile Bca 98% ( 0.6602, 0.8322 ) (0.5289, 0.8140 ) BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Calculatione and Intervals on Original Scale Based on 2000 bootstrap replicates Warning message: Im boot.ci(myboot_cbj) : CALL: boot.ci{boot.out = myboot_obj) bootstrap variances needed for studentized intervals + Abbiamo la CORRELAZIONE PUNTO-BISERIALE (misura l'associazione statistica tra una variabile continua e una variabile dicotomica [ad es. maschio/femmina]; - ilcoefficiente di correlazione punto-biseriale rp si calcola esattamente come il coefficiente di Pearson (es. usando cor.test Oi - e.ilsuo valore di significatività è lo stesso di quello di un t-test per campioni indipendenti [assumendo eguali varianze nei due gruppi]; - rà, halo stesso significato di R2, ovvero la percentuale di varianza condivisa tra la variabile continua e quella dicotomica), + la CORRELAZIONE BISERIALE (se la variabile dicotomica è stata resa tale mediante l'imposizione di una soglia, ma è intrinsecamente continua (ad es. depresso/non depresso), occorre usare il coefficiente di correlazione biseriale rs che si basa su una formula differente), # Correlazione biseriale library(polycor) with(mydata, polyserial(foot.L, gender)) [1] 0.9422766 # Per confronto, la correlazione punto-biseriale with(mydata, cor(foot.L, as.numeric(gender), use = "complete.obs")) [1] 0.7516046 footL varianza condivisa tra foot. e hand.L 3 la CORRELAZIONE PARZIALE (correlazione tra due variabili dove gli effetti di altre variabili sono tenuti costanti); Supponiamo di voler stimare la lunghezza del piede di un individuo da quella della sua mano. Sappiamo che: - r(hand.l, foot.L)=:61,R2=37:2%, - r(stature, foot.l) -> possiamo essere sicuri che la lunghezza della mano sia una tool informazione attendibile per la stima della lunghezza del piede, indipendentemente dalla statura di un individuo? Ciò che vogliamo conoscere è il valore dell’associazione statistica tra lunghezza della mano e del piede in individui della stessa statura (cioè, controllando varianza unicamente condivisa eran se tra foot. e hand per la statura) Una correlazione tra due variabili dove gli effetti di fooLL e stature altre variabili sono tenuti costanti è detta correlazione parziale. footL hand.L handL stature Per calcolare la correlazione tra x e y, controllando per z, utilizziamo la funzione pcor.test() del pacchetto ppcor: pcor.test(x, y, 2, method = c("pearson", "kendall", "spearman")). | Sostinate Ì _ _ |a Erario Nel nostro caso, la correlazione tra lunghezza della mano e del piede, controllando ——_ stature 1,0000000 0.4995456 0,4771671 mr — (R° hand.L 0.4995458 1.0000000 0.2257308 er la statura si ottiene con: librai cor] Ò (teen. O:ATTAG71 0.2257308 1.0000000 my_pcor <- with(na.omit(mydata_f), (ato stature pcor.test(hand.L, foot.L, stature, method = "pearson")) print(my_pcor) ‘statura 0.0000006+00 6.8207: cat("\nRn2 = ", round(my_pcor$estimate A 2 * 100, 1), "%\n", sep ="") Eee Sstatsstio MATRICE DI CORRELAZIONI PARZIALI e ct) Per calcolare la correlazioni per ogni coppia di variabili in un data frame, controllando S35tt 1'596128 1. assasa 0.00s000 per tutte le altre variabili, utilizziamo la funzione pcor() del pacchetto ppcor: | so earson", "kendall", "spearman" ra Nel nostro caso: library(ppcor) , pcor(na.omit(mydata_f[, 2:4]), method = "pearson") la (11) “pearson" hand ptt enel ptt » Nella correlazione parziale si controlla per l’effetto di altre variabili su entrambe le variabili poste in correlazione, » Nella correlazione semi-parziale si controlla per l’effetto di altre statue statue iobete ibi : x variabili su una soltanto delle due variabili poste in correlazione, partial correlation semi-partal correlation è e la CORRELAZIONE SEMI-PARZIALE (si controlla per l'effetto di altre variabili su una soltanto delle due variabili poste in correlazione). COEFFICIENTI DI CORRELAZIONE INDIPENDENTI Supponiamo di calcolare il coefficiente di correlazione tra lunghezze di mano e piede, per maschi e femmine: # Maschi r_male <- with(subset(mydata, gender cor(hand.L, foot.L, use = "complete.obs")) cat("Maschi, r= ", round(r_male, 2), "\n", sep Maschi, r = 0.47 # Femmine r_female <- with(subset(mydata, gendet cor(hand.L, foot.L, use = "complete.obs")) cat("Femmine, r = ", round(r_female, 2), "\n", sep = Femmine, r = 0.61 male"), "female"), È significativa questa differenza? Se due variabili x e y sono distribuite normalmente, anche la loro differenza è distribuita normalmente con: - media: Uxy=Ux-Uy varianza: 0%.y= 0°,- 0°, Può essere importante vedere se un coefficiente di correlazione è significativamente diverso da un altro coefficiente. Possiamo avere COEFFICIENTI DI CORRELAZIONE INDIPENDENTI (se due variabili X e Y sono distribuite normalmente anche la loro differenza lo sarà, con media delle differenze uguale alla differenza delle medie e varianza delle differenze uguale alla somma delle varianze; occorre poi normalizzare gli r di Pearson, poichè non è distribuito normalmente, attraverso la trasformazione di Fischer, in z, [r-> Z;=1/2In(1+r/1-r)=artanh(x)], che hanno media 0 ed errore standard _! _. Per calcolare la differenza tra due coefficienti di SN-3 r1 -> zri = artanh(r;) e poi r. -> zr. = artanh(r) e poi ne i otteniamo così una variabile con correlazione di Pearson, per prima cosa normalizziamo questi ultimi calcoliamo poi la differenza, dividendola per il suo errore standard: distribuzione campionaria normale standard), e dunque con valori di probabilità noti. Scriviamo in R una funzione che esegua questo calcolo: # Confronta due coefficienti di correlazione (Pearson) indipendenti test_r1_vs_r2_indip <- function(r1, r2, n1, n2) { zdelta <- (atanh(r1) - atanh(r2)) / sqrt(1/(n1-3)+1/(n2-3)) p_zdelta <- 2 * (1- pnorm(abs(zdelta))) print(paste("z_delta =", round(zdelta, 2))) print(paste("Two-tailed p-value =", round(p_zdelta, 3))) } Esecuzione: # Males: cor(hand.L, foot.L) = 0.47, N1 = 80 # Females: cor(hand.L, foot.L) = 0.61, N2 = 74 test_r1_vs_r2_indip(0.47, 0.61, 80, 74) [1] "2_delta = -1.21" [1] "Two-tailed p-value = 0.227" COEFFICIENTI DI CORRELAZIONE DIPENDENTI (se vogliamo per esempio confrontare il valore di due coefficienti di correlazione di Pearson provenienti dagli stessi soggetti e con una variabile condivisa [ad es. rxy e r2,]. Nel nostro esempio, potremmo voler sapere se la correlazione tra statura e lunghezza del piede (.70) è significativamente maggiore di quella tra lunghezza della mano e lunghezza del piede (.61). Occorre utilizzare una formula per ottenere una t di Student con df = N - 3, associata alla differenza tra i due t,,,-r,, = (Tzy— Tzy) ria coefficienti [serve anche il terzo coefficiente rx]). Mete ten Scriviamo la funzione in R: # Confronta due coefficienti di correlazione rxy e rzy # (Pearson) dipendenti test_my_vs_rzy_dip <- function(rxy, rzy, rxz, n) { df<-n-3 tdelta <- (my - rzy) * sart(df * (1+ rxz) / (2*(1-rxy®2-rzyA2-rmz242+2* rey *rzy* rxz))) p_tdelta <- 2 * (1 - pt(abs(tdelta), df)) print(paste("t_del , round(tdelta, 2))) print(paste("Two-tailed p-value =", round(p_tdelta, 3))) } Esecuzione: # rxy = cor(stature, foot.L) = .70 # rzy = cor(hand.L, foot.) # rxz = cor(stature, hand.L) = .90 #N=74 test_rxy_vs_rzy_dip(0.70, 0.61, 0.90, 74) [1] "t_delta = 2.38" [1] "Two-tailed p-value = 0.02" Presentazione dei risultati L’American Psychological Association suggerisce che è buona cosa: - scrivere il valore dei coefficienti di correlazione senza lo zero iniziale (ad es. r = .76), - utilizzare due cifre significative per il coefficiente, - motivare esplicitamente l'utilizzo di un test di significatività unilaterale (one-tailed), - riportare il valore di p (oggi è più comune riportare gli intervalli di confidenza). ESERCIZI PER CASA Correlazione bivariata -> Con il data set brainhead.dat già utilizzato, calcolate la correlazione tra volume cranico e massa cerebrale: selezionate il solo gruppo di età over46, graficate gli istogrammi affiancati e i Q-Q plot affiancati per maschi e femmine delle due variabili in esame, graficate gli scatter plot affiancati per maschi e femmine delle due variabili, verificate la normalità delle due variabili (test di Shapiro-Wilk) per maschi e femmine, calcolate tutti i coefficienti di correlazione appropriati (Pearson, Spearman, Kendall?) per maschi e femmine, Stabilite se i coefficienti p di Spearman siano significativamente diversi per maschi e femmine, QUARWNE Correlazione parziale -> depression.csv: 82 pazienti ricoverati per depressione, con 3 punteggi ciascuno: > simplicity: necessità di vedere il mondo in bianco e nero, 3 fatalism: tendenza al fatalismo, 3 depression: scala della depressione di Beck 1. Verificate tramite grafici la distribuzione delle variabili, - lasua significatività (p-value), I gradi di libertà (DF) sono: > peril modello: # di predittori (->1) 3 perl’errore: # di osservazioni (-> 147) - # parametri stimati (bo e bi -> 2) i l'F statistic dice se il modello preso in toto è significativo, se la varianza spiegata è maggiore di quella non spiegata, il rapporto tra varianza sistematica e varianza non sistematica; i gradi di libertà sono il numero di predittori (modello) e il numero delle osservazioni — il numero di parametri stimati (errore). INTERPRETAZIONE DEL RISULTATO: COEFFICIENTI Coefficiente: Eatimate Std. Error t value Pr(>/t]) (Intercept) 33.72471 3.31870 10.162 <2e-16 «se gup 0.06609 = 0.02984 2.245 0.0263 * In coefficients vengono riportate le stime per i coefficienti (parametri) del modello e i rispettivi valori di significatività. Ricordando l'equazione del modello: Yi = (bo + biXi) + &; si ha: - (Intercept)-> bo, - grip->b1, > (Intercept) -> bo =33:7 significa che il modello prevede per un ipotetico lavoratore con una forza di prensione uguale a 0, una valutazione da parte del supervisore uguale a 33.7, > grip -> b1 = 0:066 significa che ad un aumento di una unità nella forza di prensione corrisponde un aumento di 0.066 nella valutazione da parte del supervisore, intercept =b (valore della variabile dipendente in corrispondenza di un valore nullo della variabile esplicativa) e nome del predittore = b, (aumento del valore della variabile dipendente all'aumentare di un'unità della variabile esplicativa). Alla stima di ciascun coefficiente è associata la sua significatività (Pr(>/t/)), ossia la probabilità di ottenere da un campione casuale un valore uguale o più estremo di quello osservato, se il valore medio nella popolazione fosse zero. Ricordiamo che il valore di t si ottiene da: t(df) = b=SE», con df=N-p-1oveN=# dati, e p=#variabili esplicative nel modello. USARE IL MODELLO A SCOPO PREDITTIVO Se sostituiamo all’equazione formale del modello: Yi = (bo + b1Xi) + &;, i valori dei coefficienti stimati, otteniamo: ratings = 33:72 + 0:066 x grip da cui possiamo ricavare una previsione per la valutazione da parte di un supervisore della performance di un lavoratore che ha mostrato una data forza di prensione. RICORDATEVI DI GUARDARE | DATI library(ggplot2) mydata_scatter <- ggplot(mydata, aes(x = grip, y = ratings)) mydata_scatter + geom_point() + geom_smooth(method = "Im") + labs(x = "hand grip force", y= "performance ratings") + theme(text = element_text(size = 14)) REGRESSIONE CON 2 VARIABILI ESPLICATIVE 3 Con1solo predittore, fittare il modello equivale a trovare i coefficienti della retta che minimizza la somma delle differenze al quadrato tra i dati osservati e la retta stessa. è Con due predittori, fittare il modello Y, = (8, + 0,X,, + 8,X,;)+ €; equivale a trovare i punti del piano che minimizza la somma delle differenze al quadrato tra i dati osservati e il piano stesso. REGRESSIONE CON n VARIABILI ESPLICATIVE Con n variabili esplicative la forma generale del modello è_Y, = (8) + 8,X,; + 0,X,;+...+ 0,X,;)+ €; i fittare il modello equivale a trovare i coefficienti dell'iperpiano (piano n-dimensionale) che minimizza la somma delle differenze al quadrato tra i dati osservati e l'iperpiano stesso. Si può anche dire che ciò equivale a trovare la combinazione lineare (somma delle variabili esplicative moltiplicate per il loro coefficiente>trovare quei 8, #,. f,, P, tale che la somma sia massimamente correlata con Y) di predittori che ha la massima correlazione con la variabile dipendente. SIGNIFICATO DI MULTIPLE E ADJUSTED R? Il multiple R2 rappresenta il quadrato della correlazione tra la variabile dipendente Y; e i valori predetti dal modello È, ossia la percentuale di varianza in Y spiegata dal modello Y.= (8,+0X;+B,X,;+..+ B,X condivisa da Y e Y, spiegata dal modello). La sua effettiva utilità è tuttavia criticamente compromessa dal fatto che R? aumenta sempre all'aumentare del numero di regressori nel modello, anche quando tali variabili non avrebbero un reale potere esplicativo nella popolazione; In altre parole, il valore di R? risulta “gonfiato” (biased) rispetto a quello che si otterrebbe da un modello basato sulla intera popolazione, perché la procedura di fitting con un campione limitato tende comunque a spiegare in termini dei regressori anche quella variabilità nella Y dovuta al caso. Questa distorsione dell'A? ottenuto da un modello basato su un campione limitato è tanto maggiore quanto più piccolo è il campione e quanto più numerose sono le variabili esplicative. L'adjusted R? fornisce una misura della varianza spiegata dal modello corretta per il numero di regressori; rappresenta una stima non distorta (unbiased) per il valore di R° che si otterrebbe se il modello venisse trattato sull'intera popolazione ed è utile anche per determinare un possibile overfitting del modello (modello “sovradattato”, ossia che minimizza troppo le differenze tra i dati e il modello e quindi comincia a fittare anche le variazioni casuali perdita di capacità predittiva [modello aggiustato bene sul campione ma che non riesce a generalizzare al di fuori di esso]). )+ €; (percentuale di varianza in Y ni Modelli parsimoniosi e consigli per la costruzione dei modelli Esistono dei criteri matematici che permettono di stimare la bontà (fit) di un modello, penalizzando quelli con un maggior numero di variabili esplicative. Tra i più usati vi è il Criterio Informativo di Akaike (AIC): nin (EE) + 2k dove SSE è la sum of squared errors T per il modello, e k è il numero di regressori. -> confrontando due modelli per lo stesso campione di dati, ma con un numero diverso di regressori, è considerato migliore quello con l'AIC minore. Occorre prestare particolare cura alla scelta delle variabili esplicative, basandosi su conoscenze acquisite da ricerche precedenti e sull'importanza teorica di tali variabili. Un approccio da evitare è quello di includere nel modello tutte le variabili a disposizione e sperare in un buon risultato: nella maggioranza dei casi si diminuisce solo la capacità predittiva del modello e la sua interpretabilità. Metodi di costruzione del modello è Abbiamo il metodo SIMULTANEO o forced entry (tutti i predittori entrano nel modello senza un ordine di inserimento; è la scelta spesso migliore), 3 il metodo GERARCHICO (vengono inseriti prima i predittori basati su conoscenze teoriche pregresse e poi le nuove variabili), è il metodo STEPWISE IN AVANTI ovvero passo a passo: - si parte da un modello che contiene la sola costante Bi -_ il computercerca tra le variabili disponibili quella che ha la maggior correlazione con la variabile dipendente; - se questo regressore migliora la capacità di predire la variabile dipendente, viene aggiunto al modello; - ilcomputercerca un secondo predittore che abbia la maggior correlazione semi-parziale con la variabile dipendente, ossia che spieghi la maggior parte di varianza rimasta nella variabile dipendente dopo il fitting del primo regressore), - Rdecide di non aggiungere ulteriori predittori quando il criterio informativo di Akaike (AIC) smette di diminuire; è il metodo STEPWISE ALL'INDIETRO: - si parte da un modello che include tutti i predittori disponibili; - il computercerca la variabile la cui rimozione dal modello provoca la maggiore riduzione dell'AIC; - rimossa questa variabile, il modello viene rivalutato; - il processo si ripete fino a che la rimozione di una qualunque delle variabili rimaste provochi un aumento dell'AIC. Tale metodo è preferibile allo stepwise in avanti perché in quest'ultimo una variabile già presente nel modello può mascherare il potenziale esplicativo di un nuovo regressore). è il metodo STEPWISE IN ENTRAMBE LE DIREZIONI: - si parte come nel metodo stepwise in avanti; -. ognivolta che un predittore viene aggiunto al modello si esegue una regressione stepwise all'indietro sul nuovo modello per eliminare i predittori ridondanti; - il processo si ferma quando qualunque aggiunta o rimozione di regressori provoca un aumento dell'AIC; è il metodo SU TUTTI | SOTTOSISTEMI: il problema con i metodi stepwise è che valutano il contributo di una variabile a un modello che contiene già altre variabili. Un approccio differente consiste nel valutare e confrontare i modelli costruiti usando tutte le possibili combinazioni dei predittori disponibili -> all subject regression). Se esistono informazioni teoriche attendibili, dalla letteratura scientifica, sul ruolo predittivo di certe variabili sul fenomeno in esame è bene costruire un modello ragionevole basandosi su queste informazioni, includendo tutti i predittori rilevanti in ordine di importanza teorica; dopo questa analisi iniziale è bene ripetere la regressione escludendo le variabili che sono risultate ridondanti. In generale, minore è il numero di predittori meglio è, è bene prediligere le variabili con una solida base teorica e assicurarsi di avere un campione di dimensic ideguate. Il metodo simultaneo è quello usato dalla funzione Im() in R, ed è ritenuto da alcuni studiosi l’unico appropriato per testare una teoria. | metodi stepwise sono in genere sconsigliati perché: - levariazioni casuali dovute al campionamento rischiano di avere maggior peso nella selezione del modello rispetto alle conoscenze teoriche del ricercatore, - il modello può facilmente risultare overfitted (perde di generalità) o underfitted (non vengono identificati predittori importanti Possono essere utili a scopi esplorativi, ma è opportuno eseguire poi una convalida incrociata del modello ottenuto, ad esempio dividendo il campione in due. CONSIGLI PRATICI Se esistono informazioni teoriche attendibili, dalla letteratura scientifica, sul ruolo predittivo di certe variabili sul fenomeno in esame: - costruite un modello ragionevole basandovi su queste informazioni, includendo tutti i predittori rilevanti in ordine di importanza teorica - dopo questa analisi ale, ripetete la regressione escludendo le variabili che sono risultate ridondanti, Regole generali: - minore è il numero di predittori e meglio è, - prediligete le variabili con una solida base teorica, - assicuratevi di avere un campione di dimensioni adeguate, Valutazione del modello |: procedure diagnostiche Una volta costruito un modello basato sul campione a disposizione è importante stabilirne: lACCURATEZZA (capacità di rappresentare i dati osservati), - ela GENERALIZZABILITÀ (capacità di predire nuovi dati). Se un modello non è accurato difficilmente sarà generalizzabile, ma l'accuratezza non è garanzia di generalizzabilità (ad es. incaso di overfitting). L'accuratezza del modello dipende criticamente dalla presenza di casi anomali e casi influenti. DIAGNOSTICA DEL MODELLO: Casi anomali Nella regressione, sono quelli che si discostano dall'andamento generale degli altri casi del campione. Gli outliers possono influenzare in maniera più o meno seria la stima dei coefficienti di regressione (ad es. pendenza e intercetta della retta), ma non sono comunque dati ben rappresentati dal modello. Cam. IDENTIFICAZIONE DEI CASI ANOMALI I Possiamo cercare i casi che mostrano una grande differenza (residuo) rispetto al valore predetto dal modello; si usano i residui standardizzati, ottenuti dividendo il valore dei residui per una stima della loro deviazione standard. I residui standardizzati con valore assoluto > 3 devono destare sospetto, perché è molto improbabile che un valore così grande sia dovuto semplicemente al campionamento casuale; se più dell'1% dei casi ha un residuo standardizzato con valore assoluto > 2.5, il modello è inadeguato; se più del 5% dei casi ha un residuo standardizzato con valore assoluto > 2, il modello è inadeguato. Alcune regole utili per i residui standardizzati: — quelli con valore assoluto > 3 devono destare sospetto, perché è molto improbabile che un valore così grande sia dovuto semplicemente al campionamento casuale, — se più dell’1% dei casi ha un residuo standardizzato con valore assoluto > 2:5, il modello è probabilmente inadeguato, — se più del 5% dei casi ha un residuo standardizzato con valore assoluto > 2, il modello è probabilmente inadeguato, DIAGNOSTICA DEI MODELLI: Casi influenti Sono casi che influenzano in maniera estrema il modello, in quanto se venissero rimossi la stima dei coefficienti di regressione cambierebbe drasticamente. Un caso che ha una grande influenza può avere un residuo piccolo: è dunque bene esaminare entrambi. Diverse grandezze basate sui residui risultano utili all'identificazione dei casi influenti: deg — ADJUSTED PREDICTED VALUE (valore predetto per il caso se questo fosse escluso dall'analisi; se un caso non esercita una grossa influenza sul modello, la differenza [DFFit] tra tale valore e quello predetto dal modello originale dovrebbe essere piccola), — STUDENTIZED RESIDUAL (residuo relativo all'adjusted predicted value diviso per la sua deviazione standard; tale statistica segue la distribuzione del t di Student), — DISTANZA DI COOK (misura della influenza globale di un caso sul modello, cioè sulla sua capacità di predire tutti i casi; valori > 1 sono considerati sospetti), — Predation: rischio predatorio (da 1=min a 5=max), — Exposure: esposizione a predatori nel sonno (1-5), — Danger: esposizione complessiva al pericolo (1-5), Importiamo i dati # Importa i dati sleep0 <- read.delim("data/sleep.dat", header = TRUE, sep = "\t") # Stampa le prime 5 righe del data frame head(sleep0) # Elimina i casi con NA nelle vari sleep <- sleepO[complete.cases( sleepOl, c("TotalSleep", "Gestation", "Danger", "BrainWt")]), ] di interesse Prima versione del modello # Imposta il modello sleep_Im1 <- Im(TotalSleep * Danger, data = sleep) # Risultato summary(sleep_lm1) Seconda versione del modello sleep_lm2 <- Im(TotalSleep * Danger + Gestation + BrainWt, data = sleep) summary(sleep_lm2) Terza versione del modello sleep_lm3 <- Im(TotalSleep * Danger + Gestation, data = sleep) summary(sleep_lm3) call: Species BodyWt Brainwt NonDreaning Dreaming Totalsle na dl n fr Im(forsula = TotalSleop - Danger, data = sleep) | |a Africanelephant 6654.000 5712.0 ES (3 &fricangiantpouchedrat 1.000 6.6 63 2.0 o call: - o x = 23 (4 lm(forsula = TotalSloep - Danger 8 Î 53 Î 8 | Î Residuals: Han 19 Median 20 =.3118 caos -0.T098 2.9977 G.7989 6 Gestation 6 Brainit, data = sleep) nonna Coafficient Estimate Std. Error t value Pr(>|t1) Arr A TE (Intercept) 15.751 1.0693 14.696 < 26-16 «1 2a. PIRLO Lod Lao EP NI Danger -2:0037 © 0:3613 -5,546 9.956-07 eee RN Signif. codes: 0 «se 0.001 +» 0.01 e 0.05 da i 8 Estinate Std. Error t value Pr(>It1) (6 na) 16.8145849 0,9016814 18.648 < 20-16 #33 e deren; 3 0 Gia deere n ' -1.4224125 0.3153718 -4.510 3.936-05 ess osa ) jus: square Gestation ‘0.0211808/0/0080282 =6:212 0:000108 see ze pa i Mi Brain 0.0008795 0,0006827 1.288 0203617 Signif. codes: 0 sam 0.001 sx 0.01 «0.05 oi 1 ‘squared: 0.608, Adjusted F-statistic: 25.09 on 3 and 50 DF, p-value: 4.80! Call: | Im(formula = TotalSleep - Danger Gestation, data = sleep) AGGIORNAMENTO DEL MODELLO INR Invece della seguente sintassi nella costruzione del modello: | sleep_Im1 <- Im(TotalSleep * Danger, data = sleep) | sleep_lm2 <- Im(TotalSleep * Danger + Gestation + BrainWt, — | CIntercept) 1659845 0.89119 18.622 < 20-16 «n data = sleep) “Danger -1.50358 = 0.31100 -4.835 1.270-06 see sleep_Im3 <- Im(TotalSleep “ Danger + Gestation, data = Sestaticn -0-01607 — 0.00st1. 6.168 3.992-06 + Signif. codes: 0 em 0.001 +0 0.01 + 0.05. 0.i i sleep) Si i iù i i . Residual standard errori 3.026 on Si degresa of freedos possiamo usare più semplicemente la funzione update(): Multiple Rcaquared: 0.6876, ‘Ragustel I toe sleep_Im1 <- Im(TotalSleep * Danger, data = sleep) |Fratavistici 36.33 on 2 and'81 DF, prralue sleep_lm2 <- update(sleep_Im1, .”. + Gestation + BrainWt) sleep_Im3 <- update(sleep_Im2, .*. - BrainWt) Residuale: Min 1Q Median 39 Max -6.0169 -2.5659 -0.1704 2.0634 5.6617 Coefficiente: Estimato Std. Error + value Pr(>It1) sT14 COEFFICIENTI (PARAMETRI) DEL MODELLO Dall’ultima versione del modello (sleep_lm3) ricaviamo: TotalSleep; = bo + bi x Danger; + b. x Gestation; = 16.60 — 1.50 x Danger;— 0.02 x Gestationi I valori dei coefficienti (b) indicano quanto cambia la Y quando la variabile esplicativa corrispondente varia di una unità e gli altri predittori sono mantenuti costanti. Ogni b ha un errore standard che ne stima la variabilità da campione a campione e viene usato per determinare se b sia significativamente 6= 0 mediante t-test Se il t-test è positivo (p < 0:05) il predittore fornisce un contributo significativo al modello, tanto maggiore quanto più è piccolo il p. COEFFICIENTI STANDARDIZZATI (B) I b standardizzati5 (B) sono più semplici da interpretare perché non dipendono dalle unità di misura delle variabili In R, usiamo la funzione Im.beta() del pacchetto Im.beta: i Call: library(Im.beta) In(formula = TotalSleep - Danger Im.beta::Im.beta(sleep_lm3) MEsiUD, data = Less) 1 B indicano di quante deviazioni standard cambierà la Y quando il predittore cambia di standardized coefficients:: una deviazione standard e dunque rappresentano in maniera più efficace l’importanza (Intercept) Ca . chi csi . . io 00000000 -0.4574604 -0,4889878 di una variabile esplicativa nel modello. sCorrispondenti a un modello con variabili standardizzate. INTERVALLI DI CONFIDENZA PER | PARAMETRI Per calcolare gli intervalli di confidenza sul valore dei coefficienti di un modello: (Intercept) 14. a 416 Aa confint(sleep_lm3) Danger -2.1279454 -0.87920916 Un buon modello avrà intervalli di confidenza piuttosto stretti (e che non contengono ie ie a SIURI lo zero). REGRESSIONE GERARCHICA: CONFRONTO DI MODELLI Lo scopo è quello di vedere se un modello aumentato si adatta ai dati meglio di quello minimale, ovvero se l'R? del secondo è significativamente maggiore dell’R? del primo La significatività del fit di un modello, cioè dell’R? è calcolabile con un F-test a (k; N SSm/dfm __R°_ dfr = e SSp/dfa (4-R°)dfm dove N = # soggetti e k = # variabili esplicative. L'aggiunta di nuovi predittori al modello provoca l'aumento di R?, e possiamo calcolare il valore di F associato a questo cambiamento con la formula: F = — kaè il numero di predittori nel nuovo modello, — Ra?èilvalore di R? per il nuovo modello, — Ak=k2-kè il cambiamento nel numero di predittori, — AR2=R?.-R2è il cambiamento nel valore di R2, e la statistica F ha (‘k; N° k2 * 1) gradi di libertà. LESTOTe (EmDi AO Model 1: TotalSlsep - Danger In R usiamo la funzione anoval(): Model 2: TotalSleep - Danger anova(sleep_Im1, sleep_lm3) Rea.DÉ RSS DI Sun of Sq FP Pr6P -> sleep_lm3 ha migliorato significativamente il fit (R?) rispetto 3 Sî 164.97 1 a4a.ca 26.706 3.000-06 «ex a sleep_Im1, con F (1; 51) = 26.71, p= 3.98 x 10% Signit. codas: 0 se 0.001 +» 0.01 + 9.06. 0.1 1 DIAGNOSTICA DEI CASI ANOMALI E INFLUENTI Funzioni utili di R che hanno come argomento l’oggetto in cui è stato immagazzinato il risultato della funzione Im(): è casianomali - resid(): residui, - rstandard(): residui standardizzati, - rstudent(): residui studentizzati, fluenti - cooks.distance(): distanza di Cook, - dfbeta(): DFBeta, - dffits(): DFFit, - hatvalues(): hat values o leverages, - covratio(): covariance ratio Possiamo immagazzinare queste statistiche come variabili aggiuntive nel data frame di partenza: # Faccio una copia del data frame che contenga solo le variabili di interesse library(dplyr) sleep2 <- sleep %>% select(Species, TotalSleep, Danger, Gestation) # Aggiungo le statistiche diagnostiche per i casi anomali e influenti sleep2$residuals <- resid(sleep_Im3) sleep2$standardized.residuals <- rstandard(sleep_Im3) Gpnciea TotalSicap Dnge ontatico recite sleep2$studentized.residuals <- rstudent(sleep_Im3) + Africano ephant 8.34 6451.581090 n n — 2 Africangiantpouchedrat 8.3 3 42 -3.109730 sleep2$cooks.distance <- cooks.distance(sleep_lm3) a ArcticFox 12.6 1 60 -1.627607 . it <- . standardized residuale studentized.residuals cooks.distance sleep2$dfbeta <- dfbeta(sleep_lm3) sleep2$dffit < 4 6. 8160002 de BIAZaIO 0. d49S52068 dffits(sleep_Im3) 2 —1. 0449527 -1.0453009 0.011962313 — 8 -0. 5500626 -0. 5462659 —0.004640452 sleep2$leverage <- hatvalues(sleep_Im3) dfbeta. (Intercept) dfbeta.Danger dfbota.Cestation —1dffit leverage jance.ratio <- Ù 1 -4.376205e-02 -2.773331e-02 —1.153298e-03 0.3824224 0.28067059 sleep2$covariance. ratio covratiolsleep_Im3) 2 -4.8089350-02 -2.3117010-02 = 3.5998530-0d -O. 1896105 0,03125634 # Per comodita', aggiungo al data frame anche 8 -1.0239626-01 2.4715796-02 4,6555586-05 -0.1171744 0, 04398669 covariance.ratio fitted # i valori di Y predetti dal modello 1 1,4428588 1.T1891 < itted. 2 1.027294 11.40973 sleep2$fitted <- sleep_Im3$fitted.values È ATA # Mostra le prime 3 righe head(sleep2, 3) RESIDUI STANDARDIZZATI Se i residui sono distribuiti normalmente, ci aspettiamo che circa il 5% dei residui standardizzati abbia un valore maggiore (in modulo) di 2 # Creiamo una nuova variabile Booleana (TRUE/FALSE) che ci indichi # per ogni caso se il residuo standardizzato eccede le soglie di +/-2 # (notare l'operatore Booleano OR) sleep2$large.residual <- (sleep2$standardized.residuals > 2) | (sleep2$standardized.residuals < -2) # Numero di casi che hanno un residuo standardizzato eccessivo n_large_resid <- sum(sleep2$large.residual) # Numero di valori <-2 o >+2 che ci aspettiamo da una di: n_expected <- .05 * nrow(sleep2) ribuzione normale N residui eccessivi trovati: 1 ESAMINIAMO | CASI CON RESIDUI ECCESSIVI: sleep2[sleep2$large.residual, ] Il caso selezionato non desta preoccupazione perché: Specioa TotalSlosp Danger Costation residuale — il residuo standardizzato è < 3 15 EasternAmericannole 1 42 -6.016985 — la distanzadi Cookè<1 16 Frongardizod. residuale as «residuale cooke.distance — laleverageè<2(k+1)=n=2(2+1)=54=0.11 afbeta, (Intercept) dfbeta.Danger dfbeta.Gestation difit Leverage somi . 15 -0.3856233259 0.0877677432 —0.0002921431 -0.4579364 0.04531327 — laCVRnonè né minore di 1 + [3(k + 1) =n] = 0:83, covariance.ratio fitted large.residual né maggiore di 1 + [3(k + 1) =n] = 1.17 ti air i MULTICOLLINEARITÀ Calcoliamo la variance inflation factor: library(car) Danger Gestation # VIF per ciascun predittore Ì ART LET my_vif <- car::vif(sleep_lm3) (VIF media = 1.11 # VIF media mean_vif <- round(mean(my_vif), 2) # Stampa print(my_vif) cat(paste("\n","VIF media =", mean_vif, "\n")) Poiché nessun predittore ha una VIF > 10 e la VIF media non è sostanzialmente maggiore di 1, non vi è evidenza di seria multicollinearità. INDIPENDENZA DEGLI ERRORI PER SERIE TEMPORALI In dati con un ordinamento temporale intrinseco, la presenza di autocorrelazione (lag=1) negli errori può essere verificata con la funzione durbinWatsonTest() o dwt() del pacchetto car: library(car) dwt(my_model) La completa indipendenza degli errori corrisponde a una statistica di DW = 2: valori > 2 indicano autocorrelazione negativa e valori < 2 autocorrelazione positiva. Se DW < 1 0 > 3 la non-indipendenza è significativa. # Intercetta boot.ci (boot_resulta, type *bea", index = 1) BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 2000 bootstrap replicates CALL : boot.cifboot.out = boot_results, typa = “bca", index = i) Intervals 95% (14.53, 18.58 ) Calculations and Intervals on Original Scale Secondo le linee guide dell'American Psychological Association, # di (coefficiente di "Danger") boot.ci(boot_results, type = "bca", index # d2 (coefficiente di "Gestation") boot.ci(boot_results, type = "bca", index = 3) BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 2000 bootstrap replicates CALL = boot.ci(boot.out = boot_resulta, type = "bca", index Intervals : Level Bca 98% (-2.145, -0.923 ) Calculations and Intervals on Original Scala BOOTSTRAP_CONFIDENCE INTERVAL CALGULATIONS Based on 2000 bootstrap replicates cALL : boot.ci(boot.out = boot results, type = “bca", index = 3) Intervale Level BCa 95% (-0.0233, -0.0116 ) Calculations and Intervals on Original Scale possiamo riportare i risultati della nostra analisi di regressione come nella tabella seguente5 -> Supponendo una regressione gerarchica a 2 step: AR b SE p r Step 1 37 constant 15./2 101 < 001 Danger —2.00 035 —0.61 <.001 Step 2 22 constant 16.60 0.89 <.001 Danger 1.50 031 -—0.46 <.001 Gestation —0.02 0.003 —0.49 < .001 DUMMY CODING Uno dei requisiti della regressione è che una variabile esplicati dobbiamo scomporre la variabile categoriale in più variabili “fitti; scegliamo uno dei livelli come riferimento (ad es. valore 0 in tutte le variabi :omodo; di comodo assegniamo il valore 1 al pi nella seconda variabile di comodo assegniamo il valore 1 al tutti gli altri livelli; procediamo allo stesso modo per le rimanenti variabili di c Possiamo creare i “contrasti” che rappresentano le dummy vari: contrasti nomi che ne riflettano il significato; occorre rispettare valore del coefficiente stimato per un contrasto si riferisce alla indicati dal contrasto, con la corrispondente significatività ESEMPIO hsb2_female_small.csv: 109 studenti femmine delle scuole medie superiori americane di gruppi etnici differenti, con punteggi nelle materie di scrittura e dimatematica Variabili: id: indice dello studente race: 1=hispanic, 2=asian, 3=africanAm, 4=caucasian write: punteggio scrittura math: punteggio matematica Rinominiamo i livelli della variabile categorica: high_school$race <- factor(high_school$race, levels = 1:4, ,"caucasian")) summary(high_school) iva di tipo categoriale abbia solo 2 livelli; con più di due livelli izie” (dummy) o di comodo, ciascuna con 2 livelli soltanto (0 e 1). Se la variabile originale contiene m livelli, creiamo (m - 1) variabili di comodo; gruppo di controllo o il gruppo maggioranza) e assegniamo a esso ello di rimo livello da confrontare con ferimento e 0a tutti I secondo livello da confrontare con il livello di riferimento e 0 a comodo. ables manualmente: ciò è utile perché permette di assegnare ai l'ordine in cui sono codificati i livelli della variabile categorica. Il differenza tra i valori medi della variabile dipendente nei livelli high school <- read.csv(*data/hsb2 female small.csv", header = TRUE) summary (high school) RICODIFICA DELLA VARIABILE CATEGORICA Per la variabile race, supponiamo di prendere caucasian come livello di riferimento. Poiché abbiamo 4 livelli, creiamo le seguenti (4° 1) =3 variabili di comodo: dummyVar 1 _dummyVar 2_dummyVar 3 hispanic 1 o o asian o 1 o africanAm o o 1 caucasian o o o ASSEGNAMENTO DELL’ATTRIBUTO CONTRASTS Possiamo assegnare alla variabile categorica un attributo di tipo contrasts che ne ricodifica internamente i livelli in 123 termini di dummy variables. Si può utilizzare la funzione contr.treatment(), con la seguente sintassi: | 1100 contr.treatment(number_of_levels, base = index_of_reference_level) | g o n o dove i contrasti sono basati sul confronto di tutti i livelli con un livello che rappresenta la condizione di base | 4000 (treatment) # Esempio contr.treatment(4, base = 4) # Per il nostro data set contrasts(high school$race) <- contr.treatnent(4, base = 4) # Verifica del corretto assegnamento dell'attributo attr(high school$race, "contrasts") CODIFICA MANUALE DEI CONTRASTI Possiamo creare i “contrasti” che rappresentano le dummy variables manualment utile perché permette di assegnare ai contrasti nomi ne riflettano il significato Occorre rispettare l'ordine in cui sono codificati i livelli della variabile categorica: levels(high_school$race) [1] "hispanic" "asian" "africanAm" "caucasian" # Creazione dei contrasti hispanic_vs_caucasian <- c(1, 0, 0, 0) asian_vs_caucasian <- c(0, 1, 0, 0) africanAm_vs_caucasian <- c(0, 0, 1, 0) # Assegnazione dei contrasti alla variabile contrasts(high_school$race) <- chind( hispanic_vs_caucasian, asian_vs_caucasian, africanAm_vs_caucasian) REGRESSIONE CON VARIABILI DI COMODO high_school_Im <- Im(math * race, data = Call: Im(formla = math - race, data = high school) igh_school) summary(high_school_Im) Residuals: Min 1Q Median 39 Max -16.7792 -5.5385 0.2208 6.2208 19.4615 Coefficients: Estimate Std. Error t value Pr(>|tl) Cintercept) 53.7792 = 0.9873 Sd.A74 < 2e-16 4 racehispanic_vs caucasian -8.5065 = 2.7924 -3.046 000293 #8 raceasian_vs_caucasian 2.9708 = 3.2180 0.923 0.35804 raceafricanAm_vs_caucasian -6.2408 2.5976 -2.402 0.01804 » Signif. codes: 0 sea 0.001 »* 0.01 * 0.05. 0.1 1 Residual standard error: 8.663 on 105 degr Multiple R-squared: 0.1287, —Adjusted R-squared: 0. F-statistic: 5.169 on 3 and 105 DF, p-value: 0.002263 INTERPRETAZIONE DEI PARAMETRI STIMATI Il valore del coefficiente stimato per un contrasto si riferisce alla differenza tra i valori medi della variabile dipendente nei livelli indicati dal contrasto, con la corrispondente significatività. Nel nostro esempio, gli studenti latino-americani avevano riportato in media un voto in matematica di 8.5 punti inferiore a quello degli studenti caucasici. CODIFICA AUTOMATICA DEI CONTRASTI IN Im() # Resettiamo l'attributo "contrast" contrasts(high_school$race) <- NULL high_school_Im <- Im({math * race, data = summary(high_school_Im) igh_school) Call: Im(formla = math - race, data = high school) Residuals: Min 1Q Median 39 Max -16.7792 -5.5385 0.2208 6.2208 19.4615 Coefficients: Estimate Std. Error t value Pr>It1) (Intercept) 15.273 2.612 17.332 < 20-16 max raceasian 11.477 4.025 2.851 0.00524 raceafricanim 2.266 3.549 0.638 0.52460 racecaucasian 8.506 2.792 3.046 0.00293 # Signif. codes: 0 st 0.001 wa 0.01 a 0.058. 0.1 1 Residual standard error: 8.663 on 105 degr Multiple R-squared: 0.1287, —Adjusted R F-statistic: 5.169 on 3 and 105 DF, p-value Chamorro-Premuzic et al. (2008), Birds of a feather: Students’ preferences for lecturers’ personalities as predicted by their own personality and learning approaches. Personality and Individual Differences, 44:965-976. of freedom uared: 0.1038 0.002263 Gli autori hanno somministrato a un campione di studenti un questionario (NEO-FFI) che fornisce punteggi per i 5 tratti fondamentali della personalità (Neuroticism, Extroversion, Openness, Agreeableness, Conscientiousness). Hanno poi chiesto agli Il data set chamorro-premuzic_fixed.csv — Age:età dello studente — Gender: genere dello studente — studenth: grado di Neuroticismo dello studente — studente: grado di Estroversione dello studente — studentO: grado di Apertura all’esperienza dello studente — studentA: grado di Amicalità dello studente — studentcC: grado di Coscienziosità dello studente — lecturerN: grado di Neuroticismo desiderato in un docente — lecturerE: grado di Estroversione desiderato in un docente — lecturerO: grado di Apertura all'esperienza desiderato in un docente — lecturerA: grado di Amicalità desiderato in un docente — lecturerC: grado di Coscienziosità desiderato in un docente ANALISI Possiamo prevedere le caratteristiche di personalità desiderate in un docente dall'età, dal genere e dalla personali! studente? 1. Importatei dati in un data frame 2. Eseguite e confrontate il fitting di due modelli di regressione per il grado desiderato di Neuroticismo nel docente: -un modello semplice con i soli predittori dell’età e del genere dello studente -un modello più completo che contiene anche i tratti di personalità dello studente dello 3. Calcolate gli intervalli di confidenza dei coefficienti per il modello migliore 4. Calcolate i coefficienti standardizzati del modello 5. Verificate l'assenza di multicollinearità con la funzione car::vif() 6. Controllate la distribuzione dei residui con i seguenti grafi - residui studentizzati vs valori predetti - istogramma dei residui studentizzati - Q-Q plot dei residui studentizzati 7. Traete le opportune conclusioni e provate una regressione robusta con bootstrapping del modello 8. Ripetete il procedimento per gli altri tratti di personalità del docente Creiamo una funzione che prenda come input un data frame con i valori di Fu_indiv mesa € A) L ciascun individuo per le condizioni 1 e 2 in colonne distinte; la funzione var_nanes < nanes(nydata) ir 7 i ottificati ce — n # Nedte individuali delle due condizioni restituirà un data frame con i valori rettificati da utilizzare per il grafico a barre: individ mean <- (mydatal, 1) + mydatal, 2]) / 2 # Nedia generale grand_nean <— meanelnydata[, 1], mydatal, 21) # Differenze tra nedie individuali e media generale individ diff <- individ mean - grand nean # Valori rettificati delle variabili ydatal, 1] - indivia diff var2.adj <- nydatal, 2] — Indivia diff # Data frame rettificato come output df_out % data.frame(vari adi, var2_adj) names (df out) <= var_naneò returm(dî_out) } ESEMPIO GRAFICO setud(file.path( library (ggplot2) Sys.getenv("HOWE"), library (tidyr) "/Dacumente/Work/Didattica/reggioenilia/lezione_06")) # Conversione del data frame al formato "lunga* spider <- read .delim('dava/SpiderWide.dat", header = TRUE) spider adj long <- gather(apider adj, key = "group", value = “anziety", 1:2) head(spider, 5) # Crafica Spider bar <- ggplot(spider_adj_long, aes{x = group, y = anziety)) spider bar + atatlsunmary(fun.y = mean, geom = "bar", fill = "unite", colour = "black") + (REED Stat_simaery (fun data © sdat cl nornai, ‘gem = "linerange*, colon = "rade) è a 35 3 labelx = ‘Type of Exposure", y = “Anxiety Score") + 5 a these (text = element text (alze = 14)) 4 dò 55 8 Ca al si i | spider_adj <- rm_indiv_mean(spider) a + nead(apider adi, 5) 8 L 8» picture real O 1 48.5 To 2° 43.5 43.5 3° 41.0 46.0 è 4° 36.0 51.0 piclre Ta 5 36.0 51.0 Type of Exposure A BARRE PER MISURE RIPETUTE Viene usato per valutare se le medie di due gruppi di dati siano significativamente differenti (o anche se la media di un singolo gruppo sia significativamente diversa da un valore specificato, tipicamente zero). Vi sono due tipi di t-test, la cui rispettiva applicabilità dipende del disegno sperimentale: — independent-samples t-test (per campioni indipendenti: soggetti diversi nelle due condizioni), — epaired-samples t-test (per campioni appaiati: misure ripetute, stessi soggetti nelle due condizioni). LA LOGICA DEL T-TEST è la seguente: 1. raccogliamo due campioni di dati e calcoliamo le medie rispettive; 2. tali medie differiranno di una certa quantità; 3. sei due campioni provengono dalla stessa popolazione (ipotesi nulla) ci aspettiamo medie i, perché la probabilità di ottenere differenze grandi è bassa (se le medie del campione sono di ite normalmente i valori estremi per le differenze appartengono alle code della distribuzione e quindi hanno probabilità via via decrescente); 4. quanto simili dipende dalla deviazione standard della popolazione, stimata dall'errore standard dei dati; 5. sela differenza tra le medie osservate è grande rispetto all'errore standard, allora rifiutiamo l'ipotesi nulla e concludiamo che le medie sono significativamente differenti (cioè i due campioni provengono da popolazioni diverse). Il t-test si basa sulla T-STATISTIC, che, come la maggior parte delle test statistics, esprime il rapporto tra: — lavariabilità spiegata dal modello (effetto: differenza tra il valore osservato nei dati e il valore atteso sotto l'ipotesi nulla), — ela variabilità non spiegata dal modello (errore: entità della variazione casuale da campione a campione, stimata dall'errore standard); = più in specifico: , _ (X Xa)-(14,- Ha), dove la differenza attesa tra le medie di popolazione sotto l'ipotesi nulla, cioè pa - 48, SEx.%, è tipicamente zero. IL T-TEST COME MODELLO LINEARE Possiamo scrivere il modello statistico che rappresenta il confronto tra due medie come x = (b +b,G, )+ e, dove G è una dummy variabile con valore 0 per il gruppo A e 1 per il gruppo B, bo è la media del gruppo A e b; è la differenza tra la media del gruppo B e quella del gruppo A. Notate che: * perG=0 (gruppo A): X, = bo -> bo rappresenta la media del gruppo A * perG=1(gruppo B): Xa = (bo+ bi) ->b1=(Xg- Xa) > b1 rappresenta la differenza tra la media del gruppo B e quella del gruppo A; Poiché il t-test è essenzialmente una regressione, si applicano gli stessi requisiti: * normalità della distribuzione campionaria delle medie e, per il t-test a campioni appaiati, normalità della distribuzione campionaria delle differenze tra medie; * scala di misura dati (almeno) a intervalli. Il t-test per campio ipendet lede inoltre: * indipendenza dei valori n * e omogeneità della varianza nei due gruppi, (in pratica questo requisito non è necessario se si usa la variante di Welch del test [di default in R] che produce risultati validi anche in caso di eteroschedasticità). Il calcolo della t-statistic La t è data dal rapporto tra variabil ‘à spiegata dal modello e variabilità dovuta al caso: pa RR) () — RX SERIA: SEX A poiché tipicamente l'ipotesi nulla corrisponde a U = U» Come stimiamo l'errore standard SEX1.X,? Per la legge della somma delle varianze, la varianza della differenza tra due variabili indipendenti è uguale alla somma delle loro varianze. Se n1= n. = n, possiamo scrivere: SEX), = VVar(Xi = Ra) = VNar(xi) + var(xa) = (er +ser= (48 [+ k+SE = 7 += dove 52, e 5? sono le varianze dei campioni, e quindi: ? = XX , dove s? è la varianza e n la numerosità del campione. In questo caso si calcola prima una varianza complessiva (pooled variance) pesando le varianze dei campioni per i rispettivi gradi di libertà: (mi — 1)sî + (no — 1)s3 2 Sp= mi +na-2 oi si sostituisce tale varianza nella formula precedente: SEX} = I gradi di libertà per ilt sono (n1 + n2- 2) 1. 1Formattare idati in un data frame (long form preferibile) 2. 2Graficarei dati e calcolare statistiche descrittive (ve. 3. 3 Eseguire il test 4. 4Calcolare le dimensioni dell’effetto icando anche i req STATISTICHE DESCRITTIVE DEI DUE CAMPIONI 4 Conversione nel formato "Lungo" apider_long <- gather (spider, key = "group", valua = "anziety", 1:2) #8 Statistiche descrittive library (pasteca) by(opider_longSanziety, spider_long$group, stat.dese, ‘basic = FALSE, norm = TRUE) spider_longSgroup: picture ‘nedian sean SE.mean CI.nean.0.95 var std.dev ‘40.0000000 140.0000000 —2.GEZT168 5.9046200 86.3636364 9.2932038 Coef.var = akemess | skew.2SE lurtosis | kurt.2SE | normtest.W 0.2393301 00000000 0.0000000 -1.3939289 -0.5656047 0,9660165 notmtest.p 0,8522870 spider_long$group: real nedian mean SE.nean CI .mean.0.95 50.000000000 47.000000000 3.183765698 ‘7.007420922 121. 636363696 std.dev —’coef.var —skeness = skeu.2SE = kurtosis 11.028887688 0.254657185 -0.005590699 -0.004386224 -1.459758279 kurt.25E | normtest.W | norntest.p -0.592315868 0.948872904 0620569431 ESECUZIONE DEL TEST La sintassi del comando è la seguente: t.test(outcome * predictor, data = dataFrame) che applicata al nostro esempio diventa: ind_ttest <- t.test(anxiety * group, data = spider_long) print(ind_ttest) Welch Tuo Sanple t-test data: anciety by group © = -1.6813, df = 21.385, p-value = 0.1072 alternative hypothesis: true difference in means is not equal to 0 ‘95 percent confidence interval: -15.640641 sample estimati mean in group picture mean in group real 20 at ARGOMENTI DELLA FUNZIONE t.test() Possiamo specificare le seguenti opzioni? : - alternative = "two.sided":indica un test a due code, cioè una ipotesi alternativa X - X. #0; possiamo eseguire un test a una coda usando "less" o "greater", per le ipotesi alternative X1 - Xa<00 X1-X2>0 - mu=0:la differenza tra le medie attesa sotto - var.equal = FALSE: assume varianze diseguali nei due gruppi - conf.level=0.95: significa à dell’intervallo di confidenza 2 il valore indicato è quello di default SINTASSI ALTERNATIVA DELLA FUNZIONE t-test() Se abbiamo i dati dei due gruppi in due colonne diverse del data frame, possiamo usare il comando nella forma: t.test(mydata$x1, mydata$x2) dove x1 e x2 sono i nomi delle colonne contenenti i valori del gruppo 1 e del gruppo 2, rispettivamente. METODI ROBUSTI PER MEDIE INDIPENDENTI Nel caso di dati distribuiti asimmetricamente o con dati anomali, possiamo usare le funzioni di Rand Wilcox nella libreria WRS2 per confrontare le medie di due campioni indipendenti La prima funzione, yuen(), si basa sull'uso della media modificata o “troncata” (trimmed mean) e ha la forma generale: yuen(outcome * group, data = mydata, tr = .2, alpha = 0.05) dove sono stati indicati i valori di default per la percentuale di trimming (tr) e per la soglia di significatività (alpha). Per il nostro esempio: Library (VRS2) guen (anxiety - group, data = apider_long) call: guen (formula = aniety - group, data = spider_long) Test statistic: 1.2058 (df = 13.91), p-raluo = 0.21614 Trismed mean difference: -6.76 95 percent confidence interval: -17.9294 4.4294 Explanatory neasure of effect size: 0,38 Possiamo includere una procedura di bootstrapping nel confrontare le medie modificate: Gli autori hanno somministrato a 39 direttori commerciali e amministratori delegati il Minnesota Multiphasic Personality Inventory Scales for DSM Ill Personality Disorders (MMPI-PD), uno strumento validato per 11 disturbi della personalità: Histrionic, Narcissistic, Antisocial, Borderline, Dependent, Compulsive, Passive-aggressive, Paranoid, Schizotypal, Schizoid e Avoidant. 317 pazienti del Broadmoor Hospital classificati legalmente come psicopatici hanno costituito il gruppo di controllo. IL DATA SET: Board_Fritzon_2005.dat Non contiene i dati originali, ma solo medie, deviazioni standard e dimensioni del campione per ciascun gruppo: — disorder: tipo di disturbo della personalità — ceo.mean: media nel gruppo di CEO e direttori commerciali — psycho.mean: media nel gruppo dei pazienti psicopatici ceo.sd: deviazione standard nel gruppo di CEO e direttori commerciali — psycho.sd: deviazione standard nel gruppo dei pazienti psicopati — ceo.n: numero di soggetti nel gruppo di CEO e direttori commerciali — psycho.n: numero di soggetti nel gruppo dei pazienti psicopatici ANALISI 1. Scrivete una funzione in R che calcoli il valore del t-test a gruppi indipendenti, avendo a disposizione medie, deviazioni standard e numero di soggetti, 2. Applicate tale funzione a tutti i disturbi della personalità misurati nel vostro data set, 3. Riportate in forma scritta il risultato dei test, Volendo confrontare le medie di più di due gruppi, possiamo eseguire una serie di t-test. Se abbiamo scelto a = .05 come livello di significatività per ciascun t-test, ciò significa che la probabilità che i risultati non rappresentino un falso positivo è del 95%, la probabilità di non ottenere un risultato significativo sotto l'ipotesi nulla è del 95%. Con 3 gruppi, ad esempio, vi sono tre confronti possibili: (A vs. B), (A vs. C) e (B vs. C). La probabilità che non vi sia alcun falso positivo tra di essi è 0.95 x 0.95 x 0.95 = 0.857; la probabilità di commettere almeno un falso positivo è dunque (1 - 0.857) = 0.143, ovvero del 14.3%. -> Il rischio complessivo di errori di tipo | sull'intero esperimento (family-wise error rate) è quasi triplicato. CONFRONTI MULTIPLI E FALSI POSITIVI nl nl Con n gruppi il numero totale di confronti a coppie è dato dal coefficiente binomiale (k = 2): =—_ k) kin-h! 2-2)! ad esempio, con 5 gruppi il numero totale di t-test possibili è 10 e il family-wise error rate è pari a (1 - 0.9519) = 0.40: abbiamo il 40% di probabilità di avere almeno un falso positivo tra i risultati dei test. L'analisi della varianza (ANOVA) è una procedura che permette di eseguire confronti tra più medie controllando efficacemente il family-wise error rate. ANOVA E IPOTESI NULLA L'ipotesi nulla testata dall'ANOVA è che le medie dei vari gruppi o condizioni siano tutte uguali tra loro. La Festatistic o F-ratio prodotta dall'ANOVA è ancora una volta il rapporto tra la varianza sistematica (spiegata dal modello) e quella non sistematica (errore). Un’ANOVA significativa indica soltanto che alcune medie sono tra loro differenti, ma non ci dice quali: in altre parole, I'ANOVA è un omnibus test, un test complessivo che ci dice se, per spiegare i dati, è più vantaggioso dividerli in gruppi rispetto a non farlo, ossia se i residui lasciati dal prendere come modello la media complessiva sono maggiori dei residui lasciati prendendo medie separate per ciascun gruppo. L’ANOVA COME REGRESSIONE La logica della regressione lineare con predittori categoriali a più di due categorie si applica direttamente all'ANOVA. Con n gruppi occorre definire n - 1 variabili di comodo che rappresentino altrettanti confronti di ciascun gruppo con un gruppo scelto come riferimento. Ad esempio, in un’esperimento mirante a stabilire l'efficacia di un farmaco, possiamo avere un gruppo placebo, uno a bassa dose di farmaco e uno ad alta dose di farmaco. Definiamo le dummy variables nel seguente modo: Group dummyVar 1_dummyVar 2 placebo o o Low dose 1 o high dose o 1 NOTA SU DUMMY VARIABLES E CONTRASTI * Ilgruppoo condizione da scegliere come riferimento dipende dai confronti che intendiamo eseguire, * Nei disegni sperimentali non-bilanciati (numero di soggetti diverso in ciascun gruppo)? è importante che il gruppo di riferimento sia piuttosto grande perché le stime dei coefficienti di regressione siano attendi * Esistono schemi alternativi alle dummy variables per specificare regressori che corrispondano a confronti di interesse tra gruppi, chiamati “contrasti”, ?Per varie ragioni, è sempre meglio cercare di avere disegni bilanciati. INTERPRETAZIONE DEI COEFFICIENTI DI REGRESSIONE In un’ANOVA specificata come modello di regressione utilizzando dummy variables (DV), il modello è dato da X, = b, +b,DV.,+b,DV.,;+.;: «dove bo rappresenta la media di X per il gruppo di base (Xy), * bi rappresenta la differenza X, - Yo, * ebzrappresenta a differenza XY, - Yo. 1 valori di F e p ci dicono se l'avere usato come modello medie diverse per ciascun gruppo ha portato a previsioni significativamente migliori per il valore della variabile dipendente, ' rispetto all'uso della sola media complessiva di tutti i gruppi (l'ipotesi nulla è che le medie dei te * Frigo iena gruppi siano uguali a questa media complessiva). tr DEVIANZE (SUM OF SQUARES) : La valutazione della bontà (fit) di un modello ANOVA si basa sul calcolo di tre tipi di devianze k o scarti quadratici: dl — la TOTAL SUM OF SQUARES (SSr: devianza dei dati dalla media complessiva di tuttii‘ dati), — la RESIDUAL SUM OF SQUARES (SSx: devianza dei dati di ciascun gruppo dalle rispettive medie; è ciò che il modello non spiega, la differenza tra i dati e il modello) — e la MODELSUMIOF SQUARES (Sw: devianza delle medie di ciascun gruppo dalla media complessiva), dove SSm = St — SSa. er TOTAL SUM OF SQUARES (SST) È definita come: n SST = }(w— Karan)? = s2(N — 1) dove s2 è la varianza complessiva (grand variance) e N il numero totale di soggetti Il numero ul di gradi di libertà per SSt è: dfr=N—1. MODEL SUM OF SQUARES (SSM) È definita come: & SSm= )O TI (E; — Ryrand)” #1 dove k è il numero di gruppi, n; il numero di soggetti nel gruppo j-esimo e xj la media del gruppo j-esimo Il numero di gradi di libertà per SSM è: dfw 1 RESIDUAL SUM OF SQUARES (SSR) È definita come: SS =) Van 8 = D0-1) d=Lis1 gui = SSr — SSw dove s?;e n; sono la varianza il numero di soggetti, rispettivamente, del gruppo j-esimo Il numero di gradi di libertà per SSR è: dfr = dfr — dfmw= (N-1)—(kK-1) =N-k MEAN SQUARES E F-RATIO Le test statistics esprimono il rapporto tra variabilità sistematica dei dati (cioè spiegata dal modello) e variabilità non sistematica (non spiegata dal modello); in particolare, la F-statistic (o F-ratio) è definita come il rapporto tra le mean squares del modello e dei residui, dove queste sono le sum of squares_ normalizzate per i rispettivi gradi di libertà: MS, _ SS, 1df, F(df,,, df,)= = Qu lfu (dfy» dfx) MS, SS, /df, * REQUISITI DISTRIBUZIONALI DELL'’ANOVA Gli assunti sotto i quali la F-statistic è attendibile sono gli stessi di tutti i test parametrici basati sulla distribuzione normale: _ omogeneità della varianza tra gruppi, ipendenti tra i diversi gruppi, isurate almeno su scala a intervalli, — normalità all'interno di ciascun gruppo (normalità dei residui). Se il numero di soggetti è lo stesso in ciascun gruppo l'ANOVA è piuttosto robusta per violazioni della normalità e omogeneità della varianza. SE I REQUISITI NON SONO SODDISFATTI Ma se i requisiti non sono soddisfatti: — si possono usare i test robusti di Wilcox basati su trimmed mean, M-estimatore e bootstrapping, — in caso di eteroschedasticità usare l'F di Welch, — incaso di non normalità usare il test di Kruskal-Wallis (versione non parametrica dell’ANOVA per medie indipendenti), — trasformare i dati. CONFRONTI TRA GRUPPI E FALSI POSITIVI Se _l'ANOVA è significativa occorre individuare le differenze specifiche tra gruppi responsabili dell'effetto complessivo. Nello schema della regressione multipla otteniamo automaticamente t-test per le differenze tra gruppi corrispondenti ai coefficienti delle dummy variables. Tuttavia, eseguendo confronti statistici multipli, il rischio complessivo di falsi positivi aumenta: è necessaria una strategia per controllare il family-wise error rate di tipo |. CONTRASTIA PRIORI E POSTERIORI Come facciamo a identificare i gruppi che mostrano medie differenti in maniera più dettagliata? Il rischio complessivo di falsi positivi può essere controllato in due modi: 1. seabbiamo ipotesi specifiche sugli effetti, prima di raccogliere i dati, possiamo scomporre la varianza spiegata dal modello in componenti disgiunte tramite contrasti opportuni a priori (planned contrasts); 2. senonabbiamo ipotesi specifiche, possiamo confrontare tutti i gruppi a due a due tramite contrasti a posteriori (post-hoc contrasts), e poi applicare una correzione della significatività per il numero totale dei confronti; Occorre eseguire confronti specifici tra gruppi, tramite CONTRASTI PIANIFICATI A PRIORI LO anova (scomposizione della varianza spiegata dal modello in componenti disgiunte). Per quanto riguarda i contrasti pianificati, la prima scomposizione della varianza è quella della varianza totale, all’omnibus F-test dell’Anova, nei dati cioè al confronto tra varianza DI spiegata dal modello e varianza non spiegata dal modello. Oppure tramite CONTRASTI POST-HOC (confronto di tutti i gruppi a due a due e applicazione di una correzione della significatività per il numero totale dei confronti). La seconda scomposizione è quella della varianza spiegata dal modello e corrisponde, nel nostro caso, al confronto tra i soggetti che hanno ricevuto il farmaco e quelli che hanno ricevuto il placebo (contrasto 1: confronto tra media dei due gruppi farmaco e media del gruppo placebo); la terza scomposizione confronta i gruppi a bassa e ad alta dose di farmaco 0 Ce) L cm: (contrasto 2). dose | | ces 3 REGOLE PER I CONTRASTI PIANIFICATI In generale: — se esiste un gruppo di controllo, il primo contrasto confronta quest'ultimo con l'insieme di tutti gli altri; — ogni contrasto deve corrispondere alla divisione di un “blocco” di varianza; — una volta che un gruppo è stato isolato in un contrasto non può venire riutilizzato in un altro contrasto. Se si procede in questo modo i contrasti risultano indipendenti e ne segue che IL NUMERO MASSIMO DI CONTASTI È SEMPRE K- 1, dove k è il numero dei gruppi. + il meropo DI MioiM il METODO DI ÎNUREY -> L'hsd= honest significant difference di Tukey è una procedura che fornisce valori di significatività corretti per tutti i confronti a coppie possibili; è in genere preferibile a Bonferroni data la sua maggiore potenza statistica, specialmente quando il numero dei confronti è elevato; il metodo è esatto se il numero di soggetti è lo stesso in tutti i gruppi, ma vi è una variante [Tukey-Kramer] che fornisce valori approssimati per gruppi di dimensioni diverse; richiede l'omogeneità della varianza), > la procedura è la seguente: p-value per tutte le differenze tra gruppi tramite t-test; si ordinano i p in ordine crescente; si assegna un indice intero j a ciascun p, partendo da 1 per il p maggiore e procedendo in maniera crescente verso il p minore; io corretto per ciascun confronto è dato da pir=@/j; il metodo dà un criterio pari a quello di Bonferroni per il p più piccolo ma, procedendo via via verso p meno significativi, a si calcolan viene diviso solo per il numero di confronti rimanenti: è un metodo meno conservativo di Bonferroni; si procede in sequenza dal p più significativo in giù: una volta incontrato un test che non supera il criterio di significatività corretto, anche i rimanenti non lo supereranno il METODO DI BENJAMINI E HOCKBERG -> le tecniche precedenti hanno lo scopo di minimizzare la probabilità di commettere anche solo un errore di tipo | [falso positivo] sull'intero insieme dei confronti; Il metodo di Benjamini e Hochberg si propone invece di controllare la proporzione di ipotesi nulle erroneamente rifiutate, ovvero il rapporto tra il numero di ipotesi nulle erroneamente rifiutate e il numero totale di ipotesi nulle rifiutate [FDR (false discovery rate) = numero di falsi positivi /numero di confronti significativi]. Per calcolare il criterio di significatività corretto per i singoli test si procede nel seguente modo: - come nel metodo di Holm si ordinano i p-value per i test singoli in ordine crescente; - si assegna un indice intero ja ciascun p, partendo da 1 per il p minore e procedendo in maniera crescente verso il p maggiore; - . ilcriterio corretto per ciascun confronto è dato da poi = a(j/m), dove m è il numero totale di confronti; Si procede in sequenza dal p meno significativo verso quello più significativo: una volta incontrato il primo test che supera il criterio di significatività corretto, anche i rimanenti lo supereranno). 4Gli indici sono cioè assegnati nell'ordine opposto a Holm. BONFERRONI, HOLM E BENJAMINI-HOCHBERG Confronto Banferron Benjamini-Hochbera La tabella seguente esemplifica le tre procedure nel caso di r rem=# i i senza ld) 4 gruppi (A, B, C, D: m = 6 confronti a coppie possibili). rat nni gruppi (A, 8,G ppie possibi) Bee gaia oo * 5 ‘a Due . AGE 9127 ‘nos + 3 0250 . AGD a252 ‘nos 3 + casa . CAD rsa 008 2 5 AT Aduc sa nos 1 î “os0n PROCEDURE POST HOC E VIOLAZIONE DEI REQUISITI La maggior parte dei metodi di correzione per confronti multipli si comporta bene in caso di piccole deviazioni dalla normalità dei dati. AI contrario, disegni sperimentali non bilanciati (numero di soggetti diseguale nei gruppi) ed eteroschedasticità ne possono compromettere la validità in maniera non trascurabile; in questo caso, si possono provare metodi robusti basati su trimmed means, M-estimators e bootstrapping, anche se sono tutte procedure piuttosto nuove. RIASSUNTO SUI METODI POST-HOC Dopo un'ANOVA quindi abbiamo bisogno di ulteriori analisi per identificare le medie dei gruppi che differiscono; quando non abbiamo stabilito ipotesi precise prima dell'esperimento, eseguiamo dei test post-hoc; quando abbiamo lo stesso numero di soggetti in ogni gruppo e varianze simili, usiamo Tukey; se vogliamo un controllo rigoroso sull'errore di tipo | usiamo Bonferroni o Holm; se preferiamo controllare semplicemente la FDR usiamo Benjamini-Hochberg; se le varianze sembrano non omogenee usiamo un metodo robusto. ESEMPIO: IQ DI STUDENTI IN DIVERSE FACOLTÀ GUARDIAMO | DATI Ciuinooy | Maboveloo — Pivsico Scuse Se abbiamo ragione di pensare che le varianze della popolazione non siano uguali per tutti i gruppi, possiamo usare il test F di Welch: oneway.test(IQ * Course, data = iq_data) One-way analysis of means (not assuming equal variances) data: IQ and Course F = 34.136, num df = 2.000, denom df = 22.034, p-value = 1.781e-07 Notate che il numero di gradi di libertà dei residui è stato modificato (come nella variante di welch del t-test). Sono quelli basati sulle dummy variables con il primo gruppo in ordine alfanumerico come gruppo di riferimento CONTRASTI PIANIFICATI SPECIFICATI MANUALMENTE Supponiamo aesevistazets) di voler confrontare il gruppo di studenti di Chimica con gli studenti 1) e poi il gruppo degli studenti di Matematica (contrasto 2). Poiché i contrasti sono ortogonali, i p-value associati ai coefficienti del modello sono già corretti per confronti multipli (e a due code). isica e Matematica (contrasto ‘a con il gruppo degli studenti di Le dimensioni degli EFFETTI SPECIFICI si riferiscono ai confronti specificitra due —## Dsfinisco una funzione che Eine il calcolo — i - — — - O - rcontr <- function(my_t, my_df) gruppi. Nel caso di contrasti ortogonali possiamo anche ricavare la dimensione my_r <- sqrt(my_t © 2 / (my_t © 2 + my_df)} pr int (paste("r =", signif(my_r, 3))) degli effetti specifici dai valori di t e dai gradi di libertà: contrast } ## Confronto: Chemistry us (Maths & Physics) rcontr(6.327, 42) df=(N-p -1), e p = numero predittori nel modello. LA] "r = 0.699" COMUNICARE | RISULTATI DELL’ANOVA ESEMPIO È stato rilevato un effetto significativo del corso di laurea sull’IQ medio degli studenti, F (2; 42) = 20.02, p < .001, w = .68. | test post-hoc (Tukey) hanno identificato differenze significative tra studenti di Chimica (X +/- SE: 108.20 +/- 0.62) e studenti sia di Matematica (95.87 +/- 1.79; d= 12.3 +/- 2.26, p< .001,r = .76), che di Fisica (95.80 +/- 2.01; d= 12.4 */- 2.26, p< .001,r=.73); la differenza tra studenti di Fisica e Matematica non è risultata significativa (d = 0.07 */- 2.26). ESERCIZI DIETA E PERDITA DI PESO diet.csv: Peso corporeo in un campione di soggetti maschi e femmine, misurato prima e dopo 6 settimane di trattamento con tre tipi di dieta differenti. - . id: indice del soggetto - gender: sesso (1=maschio, 0=femmina) - age:età (anni) - height: altezza (cm) - pre.weight: peso pre-dieta (Kg) - diet: tipo di dieta - weight6weeks: peso post-dieta (Kg) ANALISI 1. 1Importate i dati e esaminateli con summary() 2. 2 Rimuovete i casi con gender = NA 3. 3 Convertite la variabile numerica gender in un fattore con etichette appropriate (0! "female", 1! "male") 4. 4Convertite la variabile numerica diet in un fattore (1! "A", 2! "B",31"C") 5. 5 Aggiungete al data frame una variabile (weight.loss) che rappresenti la differenza peso pre-post dieta 6. 6 Create un data frame con i soli soggetti femmina 7. 7 Guardate le statistiche descrittive con stat.desc() 8. 8 Graficate i dati 9. 9 Eseguite il test di Levene per l'omogeneità della varianza 10. 10 Eseguite l’ANOVA ed esaminate i grafici diagnostici 11. 11 Eseguite i test post-hoc con il metodo di Holm 12. 12 Calcolate le dimensioni dell'effetto omnibus dell’ANOVA e degli effetti specifici post-hoc 13. 13 Ripetete l’analisi per i soggetti di sesso maschile 14. 14 Riportate in forma scritta il risultato delle analisi - È l'estensione di un modello ANOVA con l'inclusione addizionale di una o più variabili continue; - Tali COVARIATE non fanno parte della manipolazione sperimentale principale ma possono contribuire significativamente alla pre ne della variabile dipendente. =. L'inclusione di covariate opportune permette di * ridurre la variabilità residua all'interno di ciascun gruppo, * edeliminare effetti confondenti (confounds). Se in una regressione gerarchica inseriamo prima le covariate e poi le dummy variables della manipolazione sperimentale, possiamo valutare gli effetti corrispondenti a quest'ultima già “puliti” dagli effetti delle covariate. | REQUISITI DELL’ANCOVA sono gli stessi dell’ANOVA, ai quali si aggiungono * l'indipendenza di covariata e manipolazione sperimentale (se ciò non avviene l'effetto associato alla covariata coincide parzialmente con quello della manipolazione sperimentale, oscurandolo o distorcendolo; per questo è importante assegnare i soggetti ai gruppi sperimentali in maniera casuale o assicurarsi che questi combacino per quanto riguarda i valori della covariata), * el’omogeneità della pendenza delle rette di regressione (rette parallele negli scatter plot per gruppi [variabile indipendente] tra variabile dipendente e covariata; per valutare l'effetto della manipolazione sperimentale sulla variabile dipendente, corretto per l'influenza della covariata, questa deve venire “scontata” allo stesso modo da tutti i gruppi [il contributo della covariata è uguale per tutti i gruppi). Quindi lo scopo dell'ANCOVA è quello di vedere se ci sono differenze significative tra le medie “tirando via” il contributo della covariata; se il contributo della covariata non è uguale in tutti i gruppi c'è un'interazione tra variabile dipendente che definisce i gruppi e covariata>in presenza di un’interazione non ha senso chiedersi se le medie sono diverse nei diversi gruppi, in quanto ciò dipende dal valore della covariata. ANCOVA INR EFFETTO DEL TIPO DI DIETA E ETÀ GUARDIAMO | DATI Variabile dipendente e covariata nei diversi gruppi setudlfile.pathl Sys.getenv("HOME"), [api i) », » weightloss_plot < lot. lata, aes(diet, weight.losa}) + /Documenta/Work/Didattica/reggioeni lia/lezione_08")) |velpitlosa plot € gqplot\apite, serCaiet, veigit.) ta <- read.cevl'data/diet_age.cev", header = TRUE) E i cao mydata$diet <- factor(mydata$diet) let <- egplot(axdata, aen(diet, age)) + om boxplot () + labaly = " a") + sumnary(aydata) | RT immane claetalte DI0) | grid.arrange (ueightloss plot, age_plot, nrov = 1) diet weight.loss 6.00 Aid Min, :-2,100 n ni 1.00 B:î4 diat Qu: 2.000 SB! d 7.00 C:15 Median : 3.600 È 4 9.12 Mean 3.893 5, So 00 3rd Qu: 6.150 È Gi 8 00 Mar. : 8.600 È. . ” n À dè A dè dist diet STATISTICHE DESCRITTIVE: VARIABILE DIPENDENTE STATISTICHE DESCRITTIVE: COVARIATA library{pasteca) \vith(nydata, by(ueight.loss, diet, stat.desc, basic = FALSE)} \vathi(apdata, DyCage, diet, otat.desc, bagic = FALSE) diet: A a senn-— Ssnenn CI mean.0.38 rn ERE ES | 2,8500000 30500000 0,5518948 1,1922963 4,2642308 2.0650014 41.4DISTIA 2.9663206 © 6.4083461 123.1968132 | 11.09805S5 coef. var 06770496 diet: B median mean SE.mean CI.mean.0.$5 var std.der ca CE RE Es © auafbtodo assortito o.eritioe © irarisler c.zseiti6 2.266625 GI saniTid/ a.oniboro = e ssbrioa sasorzie 11.680? coet.var 0,8779437 DI O median mean SE.mean CI.mean. 0.96 var std. dev ME mean SE.nean CI.nean.0.gE var std.der 136.0000000 38.0666667 13.06GY8TO 6.6204997 142.9238095 1,9550746 6.8000000 5,8800190 0.4878720 1.0463813 3,5702857 1,8895200 c0et var Î coef.var ©. 3140563 | 0.3210469 VERIFICA DELL’OMOGENEITÀ DELLA VARIANZA ## Test di Levene Levene's Test for Homogeneity of Variance (center = median) tibrary(car) aroup "a "o dani o.teso with(mydata, leveneTest(weight.loss, diet, center = median)) 40 Possiamo fare un controllo incrociato con il test dell’F-max di Hartley, valutando il rapporto tra la varianza maggiore e la varianza minore nei gruppi: Fmax <- 5.24; Fmin <- 3.57 print(round(Fmax / Fmin, 2)) [1] 1.47 Il valore è < 5 --> OK (v. Lezione 03) INDIPENDENZA DI COVARIATA E TRATTAMENTO Controlliamo che la covariata non condivida una quota rilevante di varianza con la PRA nostss TRI manipolazione sperimentale, mediante una ANOVA per la covariata in funzione dei gruppi: [ Dé Sum Sq Hean Sq F value Pr(>F) a cionificati der 2 trito se.sé 0.416 0.ess L’ANOVA non è significativa -> OK Recidualo 40 63S 1: Un modello ANOVA/ANCOVA può essere basato sul calcolo di sum of squares di » miPO | (il contributo di ogni predittore viene calcolato tenendo in considerazione i contributi dei soli predittori che lo precedono nella specificazione del modello; i valori legati alla significatività dei predittori (F, p-value...) dipendono dunque dal loro ordine nel modello, anche se il valore dei coefficienti b non cambia->le sum of squares di tipo Inon sono adatte a valutare ipotesi su main effects e interazioni, » di TIPO Il (il contributo di ogni predittore viene calcolato tenendo in considerazione i contributi di tutti gli altri predittori nel modello, eccetto le interazioni di ordine superiore in cui compare il predittore in questione; questa variante è adatta a valutare i main effects se non prevediamo un effetto delle interazioni, perché i main effects vengono calcolati senza che parte della varianza da loro spiegata possa venire “assorbita” dalle interazioni; se tuttavia un’interazione è presente, le sum of squares di tipo Il non stimano correttamente i main effects; tuttavia un vantaggio delle sum of squares di tipo Il è che il loro risultato non dipende dal tipo di contrasto utilizzato » Edi TIPO Ill (sono la variante di default in molti software statistici; il contributo di ogni predittore viene calcolato tenendo in considerazione i contributi di tutti gli altri predittori nel modello; rispetto alle sum of squares di tipo Il hanno il vantaggio di poter stimare correttamente i main effects quando sono presenti interazioni [anche se i main effects non hanno in realtà molto senso in presenza di interazioni]; sono inoltre preferibili agli altri tipi di sum of squares quando i gruppi hanno dimensioni diverse, ma richiedono contrasti ortogonali). ESECUZIONE DELL’ANCOVA library(car) ## Contrasti di Helmert (ortogonali) La stima del modello basata su Sum of Squares di Tipo Ill assicura che il contrasta (nydata$diet) <- contr.helmert (3) colnanes(contraste (nydataSdiet)) < c('_va_B", “AB_vs_c') COntributo di ogni predittore venga valutato tenendo in considerazione i Ra ear pini nere date = nydete) contributi di tutti gli altri predittori nel modello; richiede contrasti ortogonali Anova Table (Type III teate) Response: veighit.losa ‘Sun Sq Dî F value PrGF) (Intercept) 83.422 1 19.3458 B.1770-05 sv age 5.354 11/2416 0,2719755 diet 89.859 2 10.4204 0,0002368 «se Residuals 168.174 39 Signif. codes: 0 'a4s' 0.001 "ee Q.01 "e" 0.05 '.' 0.11 MEDIE CORRETTE library(effecto) Se il contributo della covariata è significativo, le medie della variabile dipendente dfect(*diet', mmodal, se = TRUE) ) nei diversi gruppi devono venire corrette per l'effetto della covariata (ADJUSTEDO gar('\aStandara Erzozo: *, cignif(adj meanssse, 3). ‘\a') —MARGINAL MEANS). Utilizziamo la funzione effect() del pacchetto effects, con la seguente sintassi: adj_means <- effect("NameOfeffect", modelObject, se = TRUE). algo: effet Possiamo poi ottenere gli errori standard delle medie corrette semplicemente con: A a ° — 3.120102 2.560658 5.ev01 adj_means$se. Lower 95 Percent Confidence Limita diet A 8 e 1.902781 1.444913 4.780602 Upper 05 Percent Confidenca Limita diet A 8 e 4.253583 3. 694198 6,9329860 Standard Errore: 0.559 0.556 0.537 aummery.1n (nynode1) I CONTRASTI PIANIFICATI NELL’ANCOVA a | coefficienti dei contrasti pianificati rappresentano differenze tra le medie corrette acu(formula = uoight_losa - age + diet, data = eydata) dove il valore di tali differenze è stato diviso per il numero di gruppi attivi nel Resttuale: 10 Nedsan 30 ax contrasto. Il valore del coefficiente della covariata rappresenta la variazione nel -9+2087 -1.0070 0.1351 1.4429. 5.0250 valore della variabile dipendente in corrispondenza della variazione di una unità °9**ficients; nella covariata. Possiamo esaminarli nel modo usuale come mostrato in immagine: inate Std. Erzor © value Pr(>It]) CIntercepe) 5.0940 1.15509 4.308 8108-05 «ev sa 0:02640 -i.114 © 0.272 aietà_va_B 0.39557 -0.700 0.488 dsetd va_O 1.00018 022201 4.505 5.894-05 see Signif. codes: 0 'ses' 0.001 "es' 0.01 "0.05 !.'0.1''1 Residual standard error: 2.077 on 39 degre: Multiple R-squarad: 0.3674, —Adjueted R-gquari Fostetistic: 7.55 on 3 and 39 DE, p-value: 0. VISUALIZZAZIONE DELL'EFFETTO DELLA COVARIATA Eseguiamo semplicemente uno scatter plot della variabile dipendente in funzione della covariata: