Connessione di R ad SQL Server 2012/2014

Uno dei più importanti aspetti di R è la gestione dei dati, i quali possono essere contenuti in diverse fonti, siano essi file csv, data base, file excel o altro ancora.

In questa breve guida verrà illustrato come effettuare la connessione tra R e Microsoft SQL Server in modo da poter estrarre i dati direttamente da un database usando le istruzioni SQL. L’approccio descritto in questa guida è supportato sia da SQL Server 2012 che dalla più recente versione 2014. È possibile connettersi ad SQL Server in diversi modi, uno dei quali è per mezzo dell’uso di ODBS, che è quello scelto in questa guida.

La prima volta che ci si collega ad un database bisogna effettuare delle operazioni preliminari, che verranno eseguite solo in questa prima occasione, che sono:

  • Creare un ODCB DSN per poter accedere ai dati;
  • Installare i packages necessari in R.

Nella screenshot sottostante viene mostrata una tabella contenente il data set d’interesse, il cui nome nel caso di specie è “TAddizionali”. La tabella fa parte del database nominato “DBTestRAddizionali”. Lo scopo è quello di importare tutti questi dati in R

1

Creazione del DSN

La prima cosa da fare è quella di realizzare DSN utente, da cui recuperare i dati, che fa riferimento al SQL Server usando ODBC”.

1)  Aprire “Strumenti di Amministrazione” e selezionare “ODBS Data Source (64 bit)”

2

2)  Dalla scheda “DSN utente” clickare su “Aggiungi”

3)  Dalla finestra appena apertasi selezionare “SQL Server” nella lista provider

4)  Bisogna ora attribuire un nome ed una descrizione alla sorgente dei dati. È bene tenere a mente questo nome perché sarà necessario più avanti nella procedura. In ultimo si deve forinre il nome del server in cui è installato l’SQL Server. A questo punto è possibile andare “Avanti”

3

5)  Selezionare il modo in cui si intende autenticare all’SQL Server. In questo caso si usa la Sicurezza Integrata, pertanto andiamo “Avanti”

4

 

6)  Si ha ora la possibilità di selezionare il database di default al quale accedere come sorgente dei dati; si sceglie di fare riferimento a “DBTestRAddizionali”. Selezioniamo quindi di andare “Avanti”

5

7)  Nell’ultima finestra che compare si ricordi di clickare su “Test di connessione”, in modo da verificare se effettivamente è possibile stabilire una connessione, e subito dopo si sceglie “Fine”

8)  Il DSN Utente è finalmente stato creato ed è attivo.

6

Installare e caricare RODBC

Il supporto per SQL Server non fa parte dei packages base di R, pertanto è necessario installare il pacchetto “RODBC” dai server CRAN.

1)  Aprire R e della console scrivere: install.packages(“RODBC”)

console

2)  Una volta scaricato ed installato il package “RODBC”, bisogna caricare il package in R così da poter usare tutte le sue funzioni. A tale scopo si digita nella console: library(“RODBC”);

3)  Il package è stato caricato e tutte le sue funzioni sono pronte per essere usate in R.

Connessione di R ad SQL Server

Arrivati a questo punto si è in grado di connettersi al database di SQL Server ed accedere al dataset d’interesse.

1)  Quando si vuole richiamare il database bisogna innanzitutto stabilire una connessione col database, e dopo aver eseguito tutte le operazioni necessarie la connessione va chiusa. In R i comandi, compresi nel pacchetto “RODBC” sono odbcConnetc()e odbcClose().

8

 

Nello screenshot sopra l’intero contenuto del database è stato attribuito all’oggetto R “prova” attraverso la funzione sqlFetch(), che non fa altro che leggere l’intero contenuto del database. L’assegnazione del contenuto del database alla variabile “prova” ha fatto sì che quest’ultima venisse interpretata da R come un data frame, contenente l’intera struttura del database.

9

A questo punto è possibile trattare “prova” come un’oggetto tipico dell’ambiente R, eseguendo su di esso qual si voglia funzione disponibile nell’ambiente.

Alfonso Izzo

La variabile età

Svincolarsi da vecchi retaggi del passato, tanto consolidati da diventare prassi, non è semplice. Le abitudini ottimizzano il nostro modo di agire, rendendolo per molti aspetti più efficiente. Capita però che i tempi cambino e certe abitudini ci portino a sprecare risorse invece di risparmiale.

Qualche decennio fa, nelle nostre vite è entrato il computer. Non solo: negli ultimi anni le prestazioni di calcolo di queste macchine sono diventate tali da consentire elaborazioni quasi impensabili fino a soli dieci anni fa. Eppure, nonostante questo, per certi aspetti continuiamo a costruire i dataset in modo simile a come venivano costruiti carta e matita.

È per esempio il caso della variabile “età”.

Codificare la variabile età in un dataset è una questione di scelte. Comunemente, se parliamo di adulti o di bambini in età scolare, l’età viene registrata in anni compiuti. Se invece si parla di bambini molto piccoli, l’età viene più spesso registrata in mesi compiuti, se non addirittura in giorni se i soggetti sono neonati. Ci sono poi i ricercatori a cui piace aggregare, per cui registrano le età come fasce: “da 30 a 35”, “da 36 a 40”, eccetera eccetera. In questo caso, le fasce sono spesso definite a seconda del legame che l’età potrebbe avere con altre variabili incluse nel dataset, per cui il modo in cui questa viene codificata è contestuale alla raccolta dati e non indipendente da essa.

La codifica in fasce d’età è secondo me retaggio di un passato in cui i computer non c’erano o avevano delle prestazioni molto limitate rispetto a quello a cui siamo abituati oggi. Anni addietro era necessario economizzare e semplificare il più possibile i calcoli: aggregare faceva comodo.

Oggi, aggregare è anacronistico. Chiaro, un minimo di livello di aggregazione ci può essere e fa comunque comodo, ma spesso si esagera, perdendo della preziosissima informazione. Ricordiamoci infatti che da un’informazione micro se ne può ricavare una macro, mentre non è sempre vero il contrario. Se io vi dicessi che ho 33 anni, potrete facilmente dedurre che posso essere collocato nella fascia “30-35”. Se però vi dicessi che mi trovo nella fascia “30-35” e vi chiedessi quanti anni ho, sapreste rispondermi?

Aggregare le età in fase di codifica dati è secondo me controproducente per tre ragioni.

Primo, ad aggregare si fa sempre in tempo in un momento successivo e più consono (sempre che usiate un software abbastanza flessibile che vi consenta di farlo).

Secondo, in fase di analisi dei dati potreste rendervi conto che quel tipo di aggregazione non è adeguato al fenomeno che state studiando e potreste desiderare delle fasce d’età costruite in un altro modo.

Terzo, l’età è una variabile continua. Per comodità noi possiamo concepire il tempo come discreto, ma ovviamente non è tale. Classificando l’età, anche semplicemente esplicitandola in anni compiuti, la stiamo in un certo senso “discretizzando”. Certo, bisogna dire che “età” non è sinonimo di “tempo”, ma è un suo derivato; inoltre, l’età è un concetto che può assumere forme diverse a seconda del contesto (fasce, anni, mesi, giorni, ecc.), mentre il tempo no. A mio parere, però, codificando l’età dovremmo per quanto possibile evitare di spezzare il tempo: dovremmo invece provare a fermarlo.

Quando possibile, io chiedo sempre che nei dataset che mi vengono forniti non venga registrata l’età, o meglio, non venga registrata solo l’età. Soprattutto se le unità statistiche sono soggetti in età evolutiva, trovo molto utile che l’età sia accompagnata da altre due variabili: la data di nascita del soggetto e la data in cui sono state rilevate le misure. L’età può infatti essere tranquillamente dedotta a partire dallo scarto fra queste due variabili, con la differenza però che in questo caso sarà chi effettua le analisi a decidere su quale unità di misura esprimerla.

R dispone di un tipo di dato particolare che consente di esplicitare che una variabile è di tipo “data”. Questo tipo di dato non è in realtà semplicissimo da gestire, principalmente a causa dell’elevatissima varietà di modi che si possono usare per codificare le date. Per esempio, la data di oggi, 29 luglio 2016, può essere scritta come “29.07.16”, “29/7/2016”, “29 lug 2016”, eccetera, eccetera, eccetera… R utilizza come formato standard per la codifica delle date quello anglosassone, per cui richiede che la data sia scritta come “anno-mese-giorno”, tutto in formato numerico e facendo in modo che i giorni e i mesi siano sempre scritti con due cifre, mentre l’anno con quattro. Per intenderci, la data di oggi in questo formato sarebbe “2016-07-29”.

Se la data è già scritta in questo formato, possiamo utilizzare il comando as.Date() per esplicitare che la stringa è una data:

> x <- as.Date("2016-07-29") 
> x 
[1] "2016-07-29"
> str(x) 
 Date[1:1], format: "2016-07-29"

Se la data non è in questo formato… beh, lo deve diventare. Attraverso l’argomento format, il comando as.Date() offre la possibilità di esplicitare la struttura del formato di origine in modo che R sia in grado di comprendere come sono disposte le informazioni all’interno della stringa.

Vediamo subito un esempio:

> x <- as.Date("29/07/2016", format="%d/%m/%Y") 
> x 
[1] "2016-07-29" 
> str(x) 
 Date[1:1], format: "2016-07-29"

All’argomento format ho passato una stringa nella quale viene dichiarato che il primo numero rappresenta il giorno in formato numerico (%d), il secondo numero rappresenta il mese sempre in formato numerico (%m) e il terzo numero rappresenta l’anno in formato numerico a quattro cifre (%Y); inoltre, ogni valore è separato da uno slash e non da un trattino.

Come facevo a sapere come costruire la stringa da passare a format? Semplice: richiamando l’help ?strptime sarà possibile studiarsi tutti i dettagli.

Torniamo ora alla questione delle età. Immaginiamo di avere un vettore born che indica la data di nascita di ogni soggetto e un vettore test che indica la data in cui è stata rilevata una misura.

> born <- c("02/12/2012","12/11/2012","17/05/2013","04/07/2013","29/02/2013") 
> test <- c("21/10/2015","21/10/2015","21/10/2015","26/10/2015","26/10/2015")
&#91;/code&#93;

I dati contenuti in queste due variabili possono essere convertiti in date:

&#91;code language="R"&#93;
> born <- as.Date(born, format="%d/%m/%Y") 
> test <- as.Date(test, format="%d/%m/%Y") 
&#91;/code&#93;

Ecco il risultato:

&#91;code language="R"&#93;
> born 
[1] "2012-12-02" "2012-11-12" "2013-05-17" "2013-07-04" "2013-02-29" 
> test 
[1] "2015-10-21" "2015-10-21" "2015-10-21" "2015-10-26" "2015-10-26"

A questo punto, per ogni soggetto è possibile calcolare l’età facendo la sottrazione fra le due variabili:

> test - born 
Time differences in days 
[1] 1053 1073  887  844  970

R esprime la differenza fra due date in giorni. Per ogni soggetto, quindi, possiamo conoscere l’età espressa in giorni.

Chiaramente, non sempre l’età espressa in giorni ci può andare bene. Se desideriamo un’unità di misura più consona, come i mesi o gli anni compiuti, possiamo ricorrere alla funzione age_calc() del pacchetto kyotil.

> library(kyotil)
# Mesi:
> age_calc(born, test, units="months") 
[1] 34.61290 35.29032 29.12903 27.70968 31.90538 
# Anni:
> age_calc(born, test, units="years") 
[1] 2.884707 2.939352 2.430137 2.312329 2.657534

La stessa funzione age_calc() del pacchetto kyotil è contenuta anche nel pacchetto eeptools, ma probabilmente quella di eeptools è una versione più vecchia, perché sembra avere delle limitazioni che quella di kyotil non ha.

Infine, non dimentichiamoci che possiamo utilizzare il comando cut(), distributo con la versione base di R, per creare delle fasce d’età. Questa, però, è un’altra storia.

Stima della sovrapposizione tra due distribuzioni empiriche con il pacchetto ‘overlapping’

In questo post vogliamo illustrare con alcuni semplici esempi come utilizzare la funzione overlap() del pacchetto overlapping per stimare il grado di sovrapposizione tra due distribuzioni empiriche.

Esempio 1


Supponiamo di avere raccolto dei dati in due gruppi di 100 soggetti ciascuno rispetto ad una generica variabile Y, espressa da punteggi teoricamente compresi tra 0 e 30, e di essere interessati a valutare se i due gruppi possano considerarsi campioni provenienti da popolazioni con la stessa media.

Di seguito si riporta il codice per la costruzione di un data frame che include la variabile Y, che contiene i punteggi, e la variabile group, che definisce a quale dei due gruppi fa riferimento ogni osservazione.

Y <- data.frame(
    Y = c(10,10,13,9,9,11,9,10,10,9,8,10,10,10,7,9,10,9,10,12,8,9,9,11,11,10,10,
        9,10,9,7,9,12,10,11,9,10,9,10,10,9,11,10,11,9,11,11,10,12,11,10,10,11,10,
        9,11,9,9,11,11,11,11,8,10,9,9,9,10,9,10,9,10,9,11,11,8,10,8,9,11,12,9,11,
        9,9,11,11,11,11,10,8,11,10,10,11,9,8,11,10,10,8,6,10,3,9,26,13,5,9,5,9,20,
        10,5,9,16,7,4,17,13,8,9,8,8,8,5,6,13,9,5,15,5,7,8,6,11,13,10,12,5,6,5,10,
        11,16,10,7,7,12,5,12,7,7,8,11,8,17,7,10,9,7,3,11,10,12,15,11,6,5,6,8,12,8,
        5,10,12,12,7,23,8,14,8,10,11,13,6,6,10,9,11,27,9,8,10,6,12,4,6,10,9),
    group = gl(2,100,labels=c("G1","G2"))
)

Nel primo campione i punteggi variano tra 7 e 13, con media 9.88 e dev. st. 1.12, nel secondo campione variano tra 3 e 27, con media 9.5 e dev. st. 4.33. Le medie sono effettivamente molto simili, ma notiamo che la variabilità associata ai due gruppi di punteggi è molto diversa, nel secondo gruppo è circa 4 volte quella del primo.

Rappresentando graficamente le distribuzioni otteniamo la seguente figura:

histogram(~Y|group,data=Y,n=50,type="count",xlab="punteggi",ylab="frequenze",border="#0080ff",xlim=c(-2,32))

Da tale figura abbiamo la conferma diretta della disomogeneità delle varianze nei punteggi dei due gruppi; inoltre, possiamo constatare che nel secondo gruppo la distribuzione dei punteggi è piuttosto asimmetrica con una evidente coda nella parte destra. In un caso di questo tipo, il confronto statistico tra le medie può essere improprio e poco informativo; ad esempio con un t-test, corretto per la non-omogeneità, otteniamo il seguente risultato: t(112.22)=0.85, p=0.397, da cui non possiamo trarre alcuna conclusione.

A questo punto proviamo a cambiare i termini del problema: piuttosto che valutare la somiglianza dei due gruppi sulla base delle sole medie (e deviazioni standard) vogliamo utilizzare tutta l’informazione disponibile nei dati. Si tratta in pratica di stimare quanta parte delle distribuzioni di punteggi nei due gruppi siano sovrapponibili; ci aspettiamo che 0% indichi l’assenza di sovrapposizione e 100% la perfetta sovrapposizione tra le due distribuzioni. Possiamo utilizzare il pacchetto overlapping nel seguente modo:

> library(overlapping)
> dataList <- list(G1=Y1,G2=Y2)
> dataList <- overlap(dataList)$OV

    G1-G2 
0.4732179 

Y1 e Y2 sono i vettori contenenti i punteggi osservati, creiamo una lista (dataList) e la passiamo alla funzione overlap(); il valore che otteniamo in output (0.4732179) è una stima della percentuale di sovrapposizione tra le distribuzioni dei punteggi Y1 e Y2. Se vogliamo avere una rappresentazione grafica delle densità stimate delle due distribuzioni con relativa area di sovrapposizione aggiungiamo l’opzione plot=TRUE:

> overlap(dataList,plot=TRUE)

La linea azzurra è la densità stimata della distribuzione di punteggi del primo gruppo, la linea magenta la densità stimata per il secondo gruppo, l’area in verde indica la porzione di sovrapposizione tra le due densità, pari a circa il 47.32%.

Con la stessa funzione si può ottenere anche un altro tipo di grafico:

> overlap(dataList,partial.plot=TRUE)

Questa seconda rappresentazione, più grezza, disegna con due linee continue (rossa e nera) le densità stimate e con una linea verde (più spessa) la parte di densità che rappresenta la sovrapposizione; le linee verticali tratteggiate rossa e nera sono poste in corrispondenza del valore massimo di densità di ciascuna curva, le linee grigie rappresentano i punti in cui la densità è la stessa per le due curve.

Esempio 2


La funzione overlap() consente di confrontare anche più coppie di distribuzioni, indipendentemente dal numero di osservazioni di ciascuna di esse. Proviamo a generare quattro distribuzioni di dati campionate dalle seguenti densità: normale (rnorm), χ2 (rchisq), esponenziale (rexp) e uniforme (runif):

> n <- 100
> set.seed(20160108)
> Y1 <- rnorm(n,10)
> Y2 <- rchisq(n,10)
> Y3 <- rexp(n,1/10)
> Y4 <- c(runif(n,0,20))
> dataList <- list(Y1,Y2,Y3,Y4)
> overlap(dataList,plot=TRUE)

Dato che ci sono 4 diverse distribuzioni, avremo 6 possibili aree di sovrapposizione, come si vede nella figura seguente:

Ciascun pannello è relativo al confronto tra le densità stimate su una coppia di distribuzioni (Yi vs Yj) rappresentate con le linee continue mentre l’area in verde indica la sovrapposizione, il cui valore stimato è riportato nel pannello stesso.

In sintesi, la funzione overlap() riceve in input una lista composta da k vettori di valori (numerici), stima la densità associata alla distribuzione di valori di ciascun vettore e quindi confronta a coppie tali densità stimando per ciascuna coppia il grado di sovrapposizione. In aggiunta permette di rappresentare graficamente le densità messe a confronto. Un primo obiettivo per le prossime versioni del pacchetto è di aggiungere opzioni che consentano all’utente una maggiore flessibilità per la rappresentazione grafica.

Visualizzare le strutture di correlazione con “qgraph”

Non so voi ma io ho una predilezione per le rappresentazioni grafiche, che si manifesta soprattutto in fase di ricognizione dei dati. Già Davide in uno dei suoi ultimi post aveva rimarcato l’importanza di visualizzare i dati, importanza che non posso che condividere! Nello specifico, il problema descritto da Davide consisteva nel rappresentare graficamente una relazione bivariata in cui le variabili considerate erano misurate su scala Likert. Tali relazioni apparivano difficilmente interpretabili con il classico grafico di dispersione (o scatterplot). La soluzione proposta – ovvero rappresentare la densità delle variabili sul diagramma – è sicuramente ottima, ma in alcuni casi può rivelarsi poco pratica, ad esempio nel caso in cui:

  1. si devono visualizzare più relazioni bivariate (ad esempio, tra i diversi item di un questionario);
  2. si vuole avere una prima idea delle relazioni che legano tra loro le variabili indagate (ad esempio, come accade quando si effettua un’indagine esplorativa dei dati).

Per l’analisi dei dati di un questionario che ho utilizzato per la mia tesi di laurea cercavo un metodo semplice e veloce per visualizzare le relazioni tra gli item. Ero già rassegnata a dover scrutare le matrici di correlazione in cerca di qualche relazione interessante, quando mi sono imbattuta per caso in un pacchetto di facile utilizzo e comprensione: qgraph. Questa libreria permette di creare una struttura reticolare in cui ogni nodo rappresenta una variabile e ogni connessione rappresenta una relazione. In particolare è possibile visualizzare:

  1. la forza della relazione attraverso lo spessore delle linee di connessione (linee più spesse indicano relazioni più forti);
  2. la direzione della relazione attraverso il colore delle linee di connessione (linee verdi indicano una correlazione positiva mentre linee rosse indicano una correlazione negativa);
  3. la significatività della correlazione in termini di p-value.

In realtà la libreria permette di fare tante altre cose interessanti, ma ritengo sia utile iniziare dalle basi ed eventualmente approfondire altre forme di visualizzazione in un diverso post. Per mostrare l’utilizzo della libreria ricorrerò ad un caso emblematico che sicuramente sarà il pane quotidiano di molti studenti – ma non solo – di psicologia: le relazioni tra gli item di un questionario.

La libreria qgraph è ampiamente utilizzata sia per l’analisi di variabili latenti che rappresentano costrutti di personalità (come i Big Five), sia per la visualizzazione delle relazioni dirette tra variabili osservate. Prima di iniziare occorre però fare una precisazione.

Bisogna ricordare che le risposte di un questionario su scala Likert sono variabili di tipo categoriale ordinale e, come tali, godono di proprietà diverse da quelle possedute dalle variabili quantitative (sareste capaci di fare la media tra “per niente d’accordo” e “assolutamente d’accordo”?). In alcuni casi la natura categoriale delle variabili ordinali viene forzata assegnando loro un numero e trattandole di fatto come numeriche. Questa forzatura ha dei pro (ad esempio, facilitare la codifica e l’interpretazione dei risultati) e dei contro. Rimando ad una lettura interessante che tratta alcuni dei “contro”.

Ora vediamo la libreria in azione!

Il file cor.item.rda è un’area di lavoro di R che contiene una matrice di correlazione (cor.item) tra i 24 item di un questionario proposto a 100 ragazzi di 14 anni. Per ogni item si chiedeva ai partecipanti alla ricerca di esprimere un giudizio su una scala Likert da 1 (“per nulla”) a 5 (“moltissimo”). Tenendo conto della natura ordinale degli item è stata calcolata la matrice di correlazione policorica attraverso la funzione auto_cor della stessa libreria qgraph (tale funzione permette di calcolare in modo statisticamente appropriato le correlazioni tra variabili quantitative, ordinali e anche tra variabili ordinali e quantitative, vedi qui).

View(cor.item)

La visione della matrice di correlazione tra i 24 item non è certo di immediata interpretazione.

Installiamo e carichiamo la libreria qgraph:

install.packages("qgraph")
library(qgraph)

Quindi, utilizziamo la funzione qgraph per visualizzare la relazione fra gli item, utilizzando l’argomento vsize per definire manualmente la dimensione dei nodi:

qgraph(cor.item, vsize=5)

Grafico 1

Il grafico mostra gli item del questionario in relazione di tipo “tutto con tutto”. Lo spessore delle linee di collegamento indica la forza della relazione: appare evidente come siano presenti correlazioni decisamente più forti di altre. Ho ritenuto opportuno tralasciare le correlazioni che avevano un indice di correlazione inferiore a .50 con l’argomento minimum=0.5.

qgraph(cor.item, vsize=5, minimum=0.5)

Grafico 2

Consideriamo la relazione fra tre item, ad esempio l’item8.b, l’item3.b e l’item7.b: il primo correla negativamente con il secondo (il coefficiente di correlazione è infatti -0.51), il quale invece correla positivamente con il terzo (il suo coefficiente di correlazione è infatti 0.59).

Una funzione utile è quella che permette di visualizzare la struttura di correlazione tra gli item del questionario raggruppandoli secondo un determinato criterio. Lo strumento che stavo studiando io presentava tre categorie di item; se notate, infatti, il nome dell’item è composto da un numero e da una lettera, che può essere a, b oppure c: è proprio questa lettera a identificare a quale raggruppamento l’item appartiene.

Per raggruppare gli item nel grafico è prima di tutto necessario creare una lista che definisca i gruppi cui fanno riferimento le variabili.

list.cor = list(
    a = c(1,4,7,10,13,16,19,22),
    b = c(2,5,8,11,14,17,20,23),
    c = c(3,6,9,12,15,18,21,24))

La lista list.cor contiene al suo interno tre vettori e ogni vettore riporta gli indici di colonna degli item di ciascun raggruppamento.

A questo punto creiamo un vettore di nome lab (labels) che contiene dei nuovi nomi da attribuire alle variabili; si tratta di un passaggio opzionale, ma che renderà il grafico più leggibile e le variabili immediatamente identificabili. Il seguente codice estrae i nomi delle variabili ed elimina la stringa “item” da ogni nome. In questo modo, la stringa “item1.a” diventerà “1.a”, la stringa “item1.b” diventerà “1.b”, e così via.

lab = gsub("item","",rownames(cor.item))

Un altro aspetto opzionale, ma sicuramente utilissimo, è l’uso dei colori. Nel vettore hue salviamo un vettore di colori che dovranno essere usati per marcare in maniera diversa i nodi che rappresentano le variabili afferenti ai tre gruppi:

hue = c("dodgerblue3","darkgoldenrod1","snow")

Infine, possiamo costruire il grafico:

qgraph(cor.item, vsize=5, minimum=0.5, groups=list.cor, color=hue, labels=lab, edge.label.cex=3)

Grafico 3

Molto più ordinato, no?

Al posto dei valori di correlazione è possibile visualizzare i p-value che indicano se gli indici possono essere considerati significativamente diversi da zero. Anche in questo caso la libreria si rivela molto comoda, perché permette di visualizzare la significatività dei parametri di regressione aggiungendo l’argomento graph = “sig2″ come segue:

qgraph(cor.item, minimum=0.5, vsize=5, groups=list.cor, graph="sig2", color=hue, labels=lab, alpha=c(0.01,0.05))

Grafico 4

Le linee blu indicano la significatività (al 5% e all’1%) delle correlazioni positive, mentre le linee arancioni indicano la significatività delle correlazioni negative (sempre al 5% e all’1%).

La mia indagine sulle relazioni tra gli item del questionario non finisce certo qui, questa libreria infatti mi sta aiutando a pianificare le indagini successive!

Spero di riuscire presto a mostrarvi altri grafici “in azione” con le funzioni sem e factanal (molto utili in psicologia).

Info e fonti:

  1. qgraph: Network Visualizations of Relationships in Psychometric Data
  2. State of the aRt personality research: A tutorial on network analysis of personality data in R
  3. Network Model Selection Using qgraph 1.3
  4. qgraph examples

Gestire le misure ripetute con i modelli a effetti misti

In un post di qualche tempo fa ho parlato di come realizzare con R la cosiddetta “ANOVA a misure ripetute”, ovvero un’analisi della varianza in cui sono presenti uno o più fattori within subjects. Le soluzioni proposte prevedevano l’uso delle funzioni aov oppure in alternativa della combinata lm + Anova (dal pacchetto car).

Per concludere il capitolo “misure ripetute”, in questo post vedremo rapidamente come risolvere la faccenda utilizzando i modelli a effetti misti. Questi, seppure molto articolati dal punto di vista teorico, sono facili da implementare in R e consentono di analizzare le misure ripetute senza dover ricorrere a bizzarri e poco trasparenti giri di codice.

Di seguito non parlerò della teoria che sta dietro il modello, un po’ perché descriverla rigorosamente mi richiederebbe un tempo esagerato, un po’ perché su internet si trova già tanto materiale, certamente migliore di quello che potrei produrre io in questo piccolo spazio. L’interpretazione dei risultati, quindi, sarà vista sempre nell’ottica dell’ANOVA.

Per prima cosa, se non l’avete ancora fatto, vi consiglio di leggere il post sull’ANOVA, almeno nella sua parte introduttiva.

I dati ormai li conosciamo. Nel repository GitHub fakedata è presente il file sonno.zip, che contiene il dataset che useremo; il dataset che ci interessa è contenuto nel file “sonno-long.csv”, che riporta i dati di un ipotetico esperimento:

In un laboratorio di psicofisiologia del sonno si stanno testando gli effetti di una nuova psicoterapia (PT) e di un trattamento di tipo farmacologico (FT) per la cura dell’insonnia. Sono state selezionate 10 persone insonni con problemi nella fase di addormentamento: 5 vengono trattate con PT e 5 con FT. Il giorno prima dell’inizio del trattamento (tempo 0), a 30 giorni e infine a 60 giorni viene registrato il numero di ore dormite da ogni individuo.

Importiamo i dati in R:

> sonno <- read.csv2("sonno-long.csv")
> str(sonno)
'data.frame':	30 obs. of  4 variables:
 $ sogg : int  1 1 1 2 2 2 3 3 3 4 ...
 $ tempo: int  0 30 60 0 30 60 0 30 60 0 ...
 $ tratt: Factor w/ 2 levels "FT","PT": 2 2 2 2 2 2 2 2 2 2 ...
 $ ore  : num  5.6 3.3 8.9 4.8 4.7 7.9 3.9 5.3 8.5 5.8 ...

Prima di partire, dobbiamo convertire le variabili sogg e tempo in tipo fattore:


> sonno$sogg <- as.factor(sonno$sogg) > sonno$tempo <- as.factor(sonno$tempo) [/code]

Gli effetti misti


Un modello lineare a effetti misti è un modello lineare che include al suo interno contemporaneamente sia effetti fissi che effetti casuali. Per “effetto” si intende una variabile che si suppone abbia una qualunque influenza sulla variabile dipendente. Nell’esperimento descritto poco sopra, si suppone che il numero di ore dormite vari in funzione del tipo di trattamento, del momento in cui è stata fatta la rilevazione e anche da una componente individuale che rende ogni soggetto unico, con delle peculiarità che lo distinguono dagli altri e che lo portano a dormire di più o di meno a prescindere dal trattamento ricevuto e dal momento della rilevazione.

Questi tre effetti, ovvero il trattamento, il tempo e il soggetto, sono rappresentati da variabili di tipo fattore; questi fattori, però, hanno una diversa natura: se trattamento e tempo sono fissi, il soggetto è casuale. Ho già parlato di questa distinzione in un vecchio post, ma forse è meglio rinfrescarsi un po’ le idee.

Un fattore è detto fisso (fixed) se i suoi livelli sono determinati a tavolino e devono essere per forza quelli e non altri, altrimenti il fattore verrebbe snaturato (es. trattamento). Tra i fattori fissi rientrano anche tutti quei fattori i cui livelli considerati sono tutti quelli possibili presenti nella popolazione (es. sesso). Un fattore è fisso a prescindere dall’essere between o within.

Trattamento e tempo sono fattori fissi, il primo between e il secondo within.

Un fattore è detto casuale (random) se i suoi livelli sono solo un campione casuale di tutti i possibili livelli presenti nella popolazione (es. soggetto). In questo caso non si è interessati a valutare le differenze fra gli specifici livelli, ma la variabilità fra essi. A un fattore casuale possono essere aggiunti livelli al fine di aumentare la potenza dell’analisi statistica, mentre quest’operazione non avrebbe senso su un fattore fisso, i cui livelli sono, appunto, in numero fisso.

L’ANOVA considera tutti i fattori come se fossero fissi, anche il soggetto. Diversamente, i modelli a effetti misti sono in grado di distinguere i fattori a seconda della loro natura. Si tratta di un aspetto che, per la verità, adesso non approfondiremo più di tanto: pragmaticamente, ci limiteremo a mimare i risultati dell’ANOVA con i modelli a effetti misti.

Già: mimare l’ANOVA, avete capito bene. I modelli a effetti misti, infatti, rappresentano una cornice procedurale molto ampia e l’ANOVA può essere vista come un suo caso particolare. Pur trattando i dati in maniera diversa, di fatto ANOVA e modelli a effetti misti possono portare a risultati analoghi. Tuttavia, i modelli a effetti misti offrono più libertà allo statistico, libertà che l’ANOVA non offre.

La funzione lme


Principalmente, esistono in R due pacchetti per adattare i modelli a effetti misti: nlme e lme4. Il primo è già distribuito con la versione base di R, il secondo no. Per alcuni versi i due si sovrappongono, ma per altri versi questi sono complementari. Qui utilizzeremo il primo.

Il pacchetto nlme (NonLinear Mixed Effects) mette a disposizione la funzione lme, che consente di adattare un modello lineare a effetti misti.


> library(nlme)
> mod <- lme(ore ~ tratt * tempo, random=~1|sogg, data=sonno) [/code] La sintassi è molto simile a lm, con la differenza che bisogna specificare nell'argomento random il fattore casuale. La notazione 1|sogg indica che l’intercetta deve variare tra i soggetti; questo aspetto, che mi rendo conto detto così potrebbe apparire come misterioso, risulterà più chiaro approfondendo, anche solo a livello superficiale, gli aspetti teorici del modello.

La funzione anova può essere utilizzata per ottenere i test di significatività sui fattori fissi:


> anova(mod)
numDF denDF F-value p-value
(Intercept) 1 16 1471.3992 <.0001 tratt 1 8 8.5321 0.0193 tempo 2 16 42.1280 <.0001 tratt:tempo 2 16 7.7107 0.0045 [/code] Questi risultati sono identici a quelli dell'ANOVA, con la differenza fondamentale che... è stato decisamente più facile ottenerli! (basta guardare il codice utilizzato nel precedente post).

Vorrei chiudere con una nota. Non c’è un accordo unanime su come vadano calcolati i p-value nei modelli a effetti misti, tanto che il prof. Douglas Bates, autore sia di nlme che di lme4, per lme4 ha preferito non fornirli in output (qui spiega il perché). Si tenga quindi in considerazione che i p vanno presi un po’ con le pinze e che l’approccio migliore per operare sarebbe l’uso del confronto fra modelli, che chi ha seguito l’ultimo corso InsulaR conosce bene!