Archivio dell'autore: AntopRet1

AntopRet1

Informazioni su AntopRet1

Psichiatra, impegnato in attività di ricerca nel campo della psicometria, dell'epidemiologia e dell'efficacia dei trattamenti

Risolvere pRoblemi

Risolvere problemi è il metodo migliore per imparare qualcosa. Ricordate la scuola?

Con la programmazione in R è un po’ la stessa cosa. Risolvere problemi favorisce la dimestichezza con il programma e facilita l’apprendimento nell’uso dei vari pacchetti. Il suggerimento è di iniziare da cose semplici, come quelle che potete trovare in un manuale o nei siti di discussione.

Di recente l’utente Umesh Acharya ha postato un quesito su CrossValidated, un sito che si occupa di problemi di statistica. Il quesito in realtà è di natura informatica, poiché riguarda la codifica di un plot. Fino a poco tempo fa, il quesito era reperibile qui, ma è stato ormai rimosso, essendo quel genere di domande più adatte a StackOverflow, il sito gemello di CrossValidated, ma più genericamente dedicato a problemi informatici.

Domanda Umesh: “Come posso mettere tre plot in un’unica finestra?”. La domanda non è banale. Umesh ha creato una funzione che consente di “stampare” nel plot la formula di una regressione usando i valori stimati per i parametri, e per farlo ha usato ggplot2, il programma di grafica creato da Hadley Wickham, mago della programmazione in R (su InsulaR se n’era già parlato qui).

A differenza del sistema grafico base, in ggplot2 non è possibile ripartire le finestre grafiche con il comando par (che avete visto usare in precedenti post). Bisogna agire sulle griglie. E bisogna saperlo.

Dunque, abbiamo un dataset con 6 variabili, che vogliamo rappresentare a coppie di due, coordinate in una regressione. Come si rappresentano i tre grafici in un’unica finestra? Iniziamo a costruire il nostro dataset, che chiameremo ‘tom’ (toy object model):

a1 <- seq(from=1, to=10, by=1)
a2 <- seq (from=2, to=20, by=2)
b1 <- seq(from=3, to=30, by=3)
b2 <- seq(from=4, to=40, by=4)
c1 <- seq(from=5, to=50, by=5)
c2 <- seq(from=2, to=29, by=3)

# combiniamo i vettori in una matrice con il comando cbind
dataframe <- cbind(a1,a2,b1,b2,c1,c2)

# trasformiamo la matrice in un dataframe con etichette per l’intestazione delle variabili
tom <- as.data.frame(dataframe)

E adesso la bella funzione creata dal nostro Umesh per inserire il risultato della regressione nel grafico:

## function to create equation expression
lm_eqn = function(x, y, df){
    m <- lm(y ~ x, df);
    eq <- substitute(italic(y) == b %.% italic(x) + a,
    list(a = format(coef(m)[1], digits = 2),
    b = format(coef(m)[2], digits = 2)))
    as.character(as.expression(eq));
}

Chiamiamo il pacchetto:

library (ggplot2)

Creiamo i tre grafici in ggplot2:

################ a1 and a2 couple of variables
P1<-ggplot(tom, aes(x=a1, y=a2)) +
geom_point(shape=1) +      # Use hollow circles

geom_smooth(method=lm,    # Add linear regression line

se=FALSE)     # Don't add shaded confidence region

P2<- P1 + labs (x= "a1 value", y = "a2 value")## add x and y labels

P3 <- P2 + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="a1 value")                            #label and change font size

P4 <- P3 + scale_x_continuous("a1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5a <- P4 + annotate("text", x = 30, y = 10, label = lm_eqn(tom$a1, tom$a2, tom), color="black", size = 5, parse=TRUE)

P5a

Variabili a1 e a2

################ b1 and b2 couple of variables
P1b<-ggplot(tom, aes(x=b1, y=b2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2b<- P1b + labs (x= "b1 value", y = "b2 value")## add x and y labels

P3b <- P2b + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="b1 value")                            #label and change font size
P4b <- P3b + scale_x_continuous("b1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5b <- P4b + annotate("text", x = 40, y = 20, label = lm_eqn(tom$b1, tom$b2, tom), color="black", size = 5, parse=TRUE)

P5b

Variabili b1 e b2

################ c1 and c2 couple of variables
P1c<-ggplot(tom, aes(x=c1, y=c2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2c<- P1c + labs (x= "c1 value", y = "c2 value")## add x and y labels

P3c <- P2c + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="c1 value")                            #label and change font size

P4c <- P3c + scale_x_continuous("c1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5c <- P4c + annotate("text", x = 40, y = 15, label = lm_eqn(tom$c1, tom$c2, tom), color="black", size = 5, parse=TRUE)

P5c

Variabili c1 e c2

E adesso il gran finale. Come rappresentare i tre grafici in un’unica finestra?

Se il nostro Umesh, peraltro bravo programmatore, fosse stato meno pigro, avrebbe scoperto subito che CrossValidated non era il sito che faceva per lui, e avrebbe scoperto in fretta almeno due soluzioni al suo problema.

Le soluzioni che qui propongo sono tratte dal sito Cookbook for R, curato da Wiston Chang, una autentica miniera di gemme.

La prima soluzione utilizza il sistema di gestione delle griglie, e fa uso di due librerie aggiuntive a ggplot2 (librerie che dovrete scaricare se volete usarle).

library(grid)
library(gridExtra)

# semplice ed elegante (il titolo lo scegliete voi)
grid.arrange(P5a, P5b, P5c, ncol = 3, main = "Your main title")
# after saving, dev.off()
# se volete cambiare le dimensioni del titolo, seguite le istruzioni

### changing dimensions of the title
grid.arrange(P5a, P5b, P5c, ncol = 3,
main=textGrob("Your main title", gp=gpar(fontsize=20,font=3)))

Combinazione grafici

La seconda soluzione fa ricorso ad una funzione creata da Chang, e messa cortesemente a disposizione degli utenti dal creatore (il risultato è analogo a quello riportato sopra).


# Multiple plot function
# ggplot objects can be passed in …, or to plotlist (as a list of ggplot objects)
# – cols: Number of columns in layout
# – layout: A matrix specifying the layout. If present, ‘cols’ is ignored.
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
# from: http://www.cookbook-r.com/

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { require(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } # Una volta che avete caricata la funzione nell’area di lavoro, l’uso è semplicissimo # (Il numero di colonne (in questo caso 3, ma può variare a piacimento): multiplot(P5a, P5b, P5c, cols=3) [/code] That's all, folks! Antonello Preti Stay Tuned for our next episode!

Stepwise regRession: Parte seconda

Riassunto delle puntate precedenti. Il mondo è tutto ciò che accade (Wittgenstein, il tizio strano Viennese). Tra le cose che accadono, alcune sembrano collegate tra loro. Molte di queste relazioni tra gli accadimenti ce le mettiamo noi (tutte, a dire il vero). Ma alcune di queste relazioni non sono mere proiezioni di un pensiero magico o superstizioso, attratto dalle coincidenze. Alcune di queste relazioni mostrano un marcata ricorrenza e coerenza.

La misurazione degli elementi che compartecipano di queste relazioni esita spesso in regolarità che si è tentati di estrapolare come “legge”.  A questo servono le regressioni: a consentire la verifica delle relazioni tra le misurazioni di elementi che appaiono, all’ispezione, collegati tra loro. Alla base di una regressione c’è un’ipotesi della relazione sottesa tra gli elementi in valutazione. Questa relazione può essere espressa in formule, o algoritmi (regole di calcolo), la cui adeguatezza (fit) è verificata dalla soluzione dell’equazione stessa in base ai dati ricavati dall’osservazione.

Quando si rilevano delle regolarità nelle associazioni tra variabili, solitamente si tratta di fenomeni che hanno base biologica, meccanica o “naturale”. Fatevi un giro nei dataframe disponibili in R per farvi un’idea.

Queste relazioni sono ricorrenti in dipendenza da dinamiche che sono intrinseche al sistema misurato.

È qui che nasce il concetto di causalità (dare un’occhiata quiqui, tanto per gradire). Insomma, quando un sistema è integrato, la variazione in una sua parte implica solitamente variazioni in parti alla prima collegate. I petali e i sepali in un fiore sono dipendenti dal “metabolismo”, una cosa magica che sostiene la vita dei fiori (e di ogni essere animato), l’insieme delle reazioni chimiche che si realizzano in un dato organismo, reazioni che implicano costrizioni e costi, che si traducono in variazioni concordi tra le componenti dell’organismo stesso. A conti fatti – e mi si passi il calembour – risulta vera e verificabile l’aforisma di Lavoisier: “Nulla si crea, nulla si distrugge, tutto si trasforma”.

La regressione stepwise, dunque. Si tratta di un sistema per semplificare una regressione multipla. La regressione stepwise è un metodo di selezione delle variabili indipendenti allo scopo di selezionare un set di predittori che abbiano la migliore relazione con la variabile dipendente.

Esistono vari metodi di selezione delle variabili.

Il metodo forward (in avanti) inizia con un modello ‘nullo’ nel quale nessuna variabile tra i predittori è selezionata; nel primo step viene aggiunta la variabile con l’associazione maggiormente significativa sul piano statistico. Ad ogni step successivo è aggiunta la variabile con la maggiore associazione statisticamente significativa tra quelle non ancora incluse nel modello, ed il processo prosegue sino a quando non vi è più variabile con associazione statisticamente significativa con la variabile dipendente.

Il metodo backward (all’indietro) inizia con un modello che comprende tutte le variabili e procede, step by step, ad eliminare le variabili partendo da quella con l’associazione con la variabile dipendente meno significativa sul piano statistico.

Il processo stepwise fa avanti e indietro tra i due processi, aggiungendo e rimuovendo le variabili che, nei vari aggiustamenti del modello (con aggiunta o re-inserimento di una variabile) guadagnano o perdono in termini di significatività.

L’uso del metodo non è senza conseguenze (vedere oltre), ma può consentire una esplorazione di relazioni che a una prima ispezione possono non essere evidenti.

Esempio, just for fun.

# carichiamo un dataset disponibile in R:

data(airquality)
attach(airquality)
dim(airquality)
head(airquality)

# Qual è la relazione tra un insieme di variabili
# climatiche e il livello di ozono nell’aria?
# Consideriamo un modello ‘additivo’, nel quale non
# si realizzano interazioni tra le variabili:

model1 <- lm(Ozone ~ Temp + Wind + Solar.R)
summary(model1)

# Consideriamo ora un modello nel quale è ammessa una
# interazione tra le variabili (molto più plausibile):

model2 <- lm(Ozone ~ Temp * Wind * Solar.R)
summary(model2)

Nell’esempio qui sopra, il cui output può essere visualizzato eseguendo il codice nel terminale di R, vengono costruiti due modelli: uno additivo (model1), che non considera le interazioni fra i tre predittori Temp, Wind e Solar.R, e uno che considera queste interazioni (model2).

Come appare immediatamente evidente, nel modello2, in cui si tiene conto delle interazioni tra le variabili (modello più che plausibile), la significatività delle associazioni tra i predittori (le variabili indipendenti) e la variabile dipendente (nel nostro caso, livelli di ozono nell’aria) scompare. Verosimilmente, questo fenomeno è causato da problemi di multicollinearità che si presentano inserendo nel modello i termini d’interazione.

Inserendo le interazioni, la varianza spiegata dal modello aumenta: l’indice R-squared passa da 0.6059 a 0.6883; non si tratta di un grande incremento, ma è pur sempre un incremento (R² aumenta dell’8%). In questo secondo modello, nonostante in apparenza nessuna delle singole variabili indipendenti sembri avere alcuna relazione con la variabile dipendente, la statistica F complessiva resta significativa:

nel modello 1 abbiamo: F(3,107) = 54.83,  p-value: < 2.2e-16
nel modello 2 abbiamo: F(7,103) = 32.49,  p-value: < 2.2e-16

Ci fidiamo del risultato negativo? No.

Proviamo ad applicare un modello stepwise.

Come si fa? Si applica il comando step al modello.

# La funzione testa le interazioni ed elimina quelle
# che incrementano l'AIC:

model3 <- step(model2)
summary(model3)

# La varianza spiegata è rimasta più alta che nel modello
# semplicemente ‘additivo’ (R² =  0.6677)

Bella differenza, vero? Ora il modello restituisce alcune interazioni come significative, ma soprattutto tornano ad essere significative le relazioni già osservate nel più semplice modello additivo.

Che cosa è l’AIC? L?AIC è l’Akaike Information Criterion, un indice che consente di comparare due modelli e che consente di valutare variazioni nell‘adattamento di un modello per modifiche dello stesso (vedere qui per dettagli).

L‘AIC può essere interpretato solo in ottica comparativa: il modello migliore è sempre quello con l’AIC più basso. Quindi, quello a cui si deve badare non è il valore dell’AIC in sé, ma la differenza in AIC tra i modelli, che può essere valutata in questo modo (adattabile ad ogni successione di modelli):

# Create un vettore che comprenda tutti i valori AIC
# dei modelli (qui esempio con tre modelli):

x <- c(AIC(model1), AIC(model2), AIC(model3))

# Calcolate la differenza:

delta <- x - min(x)
delta

# Vedete da voi che il modello 3 è quello con differenza
# pari a 0, cioè, è il modello migliore.

# Potete anche derivare il ‘peso’ (weight) dell’AIC nei vari
# modelli, essendo il ‘peso’ relativo dell’AIC  di un modello
# rispetto agli altri una evidenza in suo favore:

# Calcolate la verosimiglianza (likelihood) dei modelli

L <- exp(-0.5 * delta)

# Akaike weights:

w <- L/sum(L)
w

# Si può anche utilizzare un altro indice, il Bayesian
# information criterion (BIC), che, a differenza dell'AIC,
# penalizza maggiormente modelli più complessi:

x_BIC <- c(BIC(model1), BIC(model2), BIC(model3))
delta_BIC <- x_BIC - min(x_BIC)
delta_BIC
L_BIC <- exp(-0.5 * delta_BIC)

# BIC weights:
w_BIC <- L/sum(L_BIC)
w_BIC

Ricca messe di informazioni sul tema qui. Dal sito, curato da Dolph Schluter, ho tratto il metodo rapido per calcolare l’AIC e il BIC.

Si diceva dei limiti di questa procedura. Il critico più agguerrito dell’uso della procedura stepwise è Frank Harrell, un Big Kahuna della statistica delle regressioni (non avete visto il film? Male! Monologo finale).

Le obbiezioni di Harrell sono da prendere con beneficio di inventario, ma non andrebbero trascurate:

1. R² values are biased high.
2. The F and χ² test statistics do not have the claimed distribution.
3. The standard errors of the parameter estimates are too small.
4. Consequently, the confidence intervals around the parameter estimates are too narrow.
5. p-values are too low, due to multiple comparisons, and are difficult to correct.
6. Parameter estimates are biased high in absolute value.
7. Collinearity problems are exacerbated.

Tuttavia, i benefici della stepwise regression sono indubitabili. Come nell’esempio riportato sopra, consente di identificare relazioni oscurate dall’inserimento di tutte le variabili nel modello.

Va ricordato che una regressione è un algoritmo che consente di stimare una variabile dipendente a partire da un insieme di predittori. Talvolta la stima di tale variabile può essere complicata e richiedere molte risorse; più risorse sono richieste maggiore sarà il tempo impiegato dall‘algoritmo. Inoltre, è necessario conoscere in anticipo la relazione matematica che lega i predittori alla variabile dipendente.

L’esito di una regressione sono dei parametri in base ai quali ‘correggere’ i punteggi dei predittori per ottenere una approssimazione della variabile dipendente, approssimazione che è tanto più precisa quanto migliore il modello (vedere post precedenti sulla valutazione della adeguatezza dei modelli – qui e qui).

Esempio. Poniamo che voi vogliate stimare il valore di una sostanza tossica liberata da una fabbrica. Quando i livelli di quella sostanza tossica sono elevati, c’è rischio di contaminazione. È costoso dosare la sostanza, e inoltre vi serve sapere in anticipo se la concentrazione nell’aria sta salendo verso il livello di guardia. Sapete che alcuni indici sono legati al rischio di innalzamento delle concentrazioni della sostanza tossica: temperatura, pressione nelle tubature, grado di umidità dell’aria. Queste variabili sono state stimate per un periodo di otto mesi, e una regressione stepwise le ha estratte come predittori efficienti della concentrazione della sostanza tossica.

I coefficienti di questi predittori sono:
– Temperatura: βT = 2.45
– Pressione: βP = 11.11
– Umidità: βU = –3.53
Mentre il valore dell‘intercetta è: α = –13.81.

Il modello di stima dei livelli di concentrazione della sostanza tossica è:

Sostanza_tossica = α + βT * Temperatura + βP * Pressione + βU * Umidità

Con un sistema di sensori ottenete, ogni dieci minuti, i valori di Temperatura, Pressione e Umidità.

In base ai coefficienti calcolati dal vostro modello, potete calcolare i valori presunti della sostanza tossica in base ai predittori semplicemente applicando la vostra formula:

Sostanza_tossica = – 13.81 + (2.45 x [valore misurato di Temperatura]) + (11.11 x [valore misurato di Pressione]) + (–3.53 x [valore stimato di Umidità])

Quando il valore predetto raggiunge una soglia di rischio, a quel punto effettuate la misurazione dei livelli della sostanza tossica. Poiché lo fate solo quando si raggiunge una soglia di rischio, risparmiate dei bei soldini…

Pensate a quello che si può fare con un buon modello matematico di predizione (che significa, buone analisi ma anche buoni dati), un pacchetto di sensori, e un processore Arduino

Alternative alla stepwise regression? Nessuna facile da capire. Ne parleremo in un prossimo post.

Intanto, per cominciare, date un’occhiata qui.

That’s all, folks!

Antonello Preti

Stay Tuned for our next episode!

Should I stay or should I go? Stepwise regRession – parte prima

Ah, i favolosi anni Ottanta, con le minigonne a palloncino, i capelli cotonati, la pelle liscia senza (ancora) tatuaggi… e i Clash, appunto…

L’argomento della puntata è la regressione a passi cauti (stepwise… OK, I’m joking…). L’argomento è controverso, e lo prenderemo alla lontana.

A dire il vero, R è un sistema nel quale la regressione è “cucinata” in tutte le salse. Se esiste una qualche incertezza sul perché R si chiami proprio così (forse dalle iniziali dei nomi dei primi programmatori, Ross Ihaka e Robert Gentleman), si può avanzare la ragionevole ipotesi che R stia per Regression. In effetti, una larga fetta di pacchetti statistici in R sviluppa, in via diretta o indiretta, una variazione sul tema della regressione: lineare, non-lineare, multivariata, robusta o quello che volete voi.

Non è facilissimo definire il concetto di regressione. Certo, su Wikipedia si trova che “La regressione formalizza e risolve il problema di una relazione funzionale tra variabili misurate sulla base di dati campionari estratti da un’ipotetica popolazione infinita” (vedere per credere). Tuttavia il concetto di regressione definisce più una filosofia di vita che un modello matematico. Al fondo delle credenze incentrate sui modelli di regressione è l’idea che il mondo là fuori abbia una razionalità intrinseca, che esista la possibilità di formulare previsioni sull’andamento degli eventi, e, insomma, la fede incrollabile che “God does not play dice” (ma, a sentire Stephen Hawking, “God still has a few tricks up his sleeve“, insomma, si tiene qualche asso nella manica…).

Al dunque, “l’analisi della regressione è una tecnica usata per analizzare una serie di dati che consistono in una variabile dipendente e una o più variabili indipendenti” (Wikipedia, of course). Sfortunatamente, questo nucleo incontrovertibile di ogni modello di regressione che si rispetti è quello sul quale si appunta la critica – per nulla prescindibile – di David Hume: “Che domani il sole non sorgerà, è una proposizione non meno intelligibile, e non implica più contraddizione, dell’affermazione secondo cui esso sorgerà; tenteremmo invano, perciò, di dimostrarne la falsità. Se fosse dimostrativamente falsa, essa implicherebbe una contraddizione e non potrebbe mai essere concepita dalla mente in modo distinto”.

La relazione tra le variabili, direbbe il buon Jacques, ce la mettiamo noi, è una relazione simbolica. I “dati” sono astrazioni, l’esito di operazioni di misurazione, che risentono non solo dell’errore implicito in ogni misurazione (e di ciò si tiene conto nei modelli matematici di regressione), ma anche delle “credenze” implicate nell’atto stesso della misurazione. Se crediamo di poter predire l’andamento dei raccolti dalle condizioni di semina e da quelle climatiche è perché coltiviamo (mi si perdoni il calembour) l’idea che gli eventi/fatti “semina”, “condizioni metereologiche” e “raccolto” siano tra loro collegati.

Il concetto di regressione implica – infatti – una qualche relazione tra una “causa” e un “effetto”, che il modello in valutazione ha il compito di stimare. Nella sua forma più semplice, la regressione è un confronto tra le misure di un fatto/evento lungo una dimensione e le misure del medesimo fatto/evento lungo una dimensione collegabile alla prima.

Esempio, esempio, reclama il mio fedele pubblico di 25 lettori

La regressione stepwise è applicata alla regressione multipla. Partiamo dunque da qui.

Conciossiacosaché si voglia misurare la lunghezza e la larghezza dei petali in un fiore e valutare se dette misure predicano la lunghezza dei sepali dello stesso fiore.

# caricate il dataframe (già disponibile in R)
data(iris)

# check preliminare sul dataset (approfondisci: www.insular.it/?p=1698)
dim(iris)
head(iris)

# per accedere direttamente alle variabili (approfondisci: www.insular.it/?p=2245)
attach(iris)

L’ipotesi è che in un insieme di fiori di varia specie la lunghezza e la larghezza dei petali predica la lunghezza dei sepali dei medesimi fiori, cioè, detta in regressionese, che “lunghezza e larghezza dei petali” siano le variabili indipendenti che predicono la distribuzione della variabile dipendente “lunghezza dei sepali”:

In R, la funzione che specifica una relazione lineare tra due variabili e che di conseguenza adatta un modello di lineare è lm. La relazione tra variabile dipendente (la y nei grafici cartesiani) e la o le variabili indipendenti (la x nei grafici cartesiani) è specificata dalla tilde “~”:

fit1 <- lm(Sepal.Length ~ Petal.Length + Petal.Width, data = iris)
summary(fit1)

Come output, ci verranno restituite le principali statistiche descrittive sui residui, il test di significatività sui parametri e infine gli indici di adattamento e di significatività dell’intero modello. La tabellina “Coefficients” riporta, per ogni parametro, il valore (Estimate), l’errore standard (Std. Error) e i valore della statistica t (t value) e il p-value (Pr(>|t|)).

Call:
lm(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.18534 -0.29838 -0.02763  0.28925  1.02320 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.19058    0.09705  43.181  < 2e-16 ***
Petal.Length  0.54178    0.06928   7.820 9.41e-13 ***
Petal.Width  -0.31955    0.16045  -1.992   0.0483 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4031 on 147 degrees of freedom
Multiple R-squared:  0.7663,	Adjusted R-squared:  0.7631 
F-statistic:   241 on 2 and 147 DF,  p-value: < 2.2e-16

Il modello sopra illustrato è additivo, si suppone cioè che le due variabili indipendenti si sommino nel loro impatto ma non interagiscano. Vogliamo formulare l’ipotesi che le due variabili indipendenti interagiscano in modo moltiplicativo? Usiamo il segno “*” anziché il +:

fit2 <- lm(Sepal.Length ~ Petal.Length * Petal.Width, data = iris)
summary(fit2)

L’output sarà il seguente:

Call:
lm(formula = Sepal.Length ~ Petal.Length * Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.00058 -0.25209  0.00766  0.21640  0.89542 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               4.57717    0.11195  40.885  < 2e-16 ***
Petal.Length              0.44168    0.06551   6.742 3.38e-10 ***
Petal.Width              -1.23932    0.21937  -5.649 8.16e-08 ***
Petal.Length:Petal.Width  0.18859    0.03357   5.617 9.50e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3667 on 146 degrees of freedom
Multiple R-squared:  0.8078,	Adjusted R-squared:  0.8039 
F-statistic: 204.5 on 3 and 146 DF,  p-value: < 2.2e-16

Interessante: la varianza spiegata (Adjusted R-squared) aumenta, passando dal 76.6% all’80.8% e cambia anche l’impatto delle singole variabili indipendenti (i valori stimati per i coefficienti).

Controlliamo la validità del modello con le diagnostiche di regressione:

# regression diagnostics
par(mfcol = c(2,4))
plot(fit1,main="modello additivo")
plot(fit2, main="modello con interazione")

Regression diagnostics

A occhio, sembrerebbe che il modello moltiplicativo, che prevede l’interazione delle due variabili predittive, vada meglio, soprattutto la curva dei residui cade maggiormente sulla linea orizzontale.

Possiamo anche verificare che il secondo modello sia effettivamente migliore del primo facendo uso della funzione anova:

anova(fit2,fit1)

L’output evidenzia una differenza significativa fra i due modelli, di conseguenza sembra davvero che il secondo sia migliore del primo.

Analysis of Variance Table

Model 1: Sepal.Length ~ Petal.Length * Petal.Width
Model 2: Sepal.Length ~ Petal.Length + Petal.Width
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    146 19.637                                  
2    147 23.881 -1   -4.2441 31.556 9.502e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Siamo però sicuri che vada tutto bene? Prima di proseguire dovremmo verificare che siano rispettati alcuni requirements della regressione lineare multipla. Per esempio, bisogna verificare la normalità dei residui; per effettuare questa operazione utilizziamo la funzione qqPlot della libreria car che genera un quantile-quantile plot per i residui studentizzati.

library(car)
# qq-plot for studentized resid
par(mfcol = c(1,2))
qqPlot(fit1, main="QQ Plot – modello additivo")
qqPlot(fit2, main="QQ Plot – modello con interazione")

QQplot for studentiized residuals

Fin qui, tutto bene. Vi sono violazioni della omoscedasticità dei residui? Anche per queste analisi serve il pacchetto “car”:

library(car)
# non-constant error variance test
ncvTest(fit1)
ncvTest(fit2)
# plot studentized residuals vs. fitted values
par(mfcol = c(1,2))
spreadLevelPlot(fit1, main=paste("Spread-Level - modello additivo"))
spreadLevelPlot(fit2, main=paste("Spread-Level - modello con interazione"))

Ecco l’output:

> ncvTest(fit1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.121054    Df = 1     p = 0.2896917

> ncvTest(fit2)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6024558    Df = 1     p = 0.4376425

QQplot for studentiized residuals

Sembra che il modello si adatti con con pochi problemi. Ci possiamo accontentare? Non ancora… dobbiamo valutare l’eventuale presenza di multicollinearità. Il metodo più semplice è sicuramente quello di valutare l’entità della correlazione tra i predittori:

cor(Petal.Length, Petal.Width)

L’indice di Bravais-Pearson risulta pari a 0.96, decisamente troppo elevato. I due predittori dovrebbero infatti essere poco correlati, altrimenti le stime dei parametri risultano molto instabili. Per individuare la presenza di multicollinearità è spesso utilizzato uno stimatore noto come variance inflation factors (VIF), cioè “fattori di crescita della varianza”. Valori superiori a 1 indicano rischio di multicollinearità. Valori superiori a 10 suggerisco seri problemi di multicollinearità.

# Evaluate Collinearity
# modello additivo
vif(fit1) # variance inflation factors
sqrt(vif(fit1)) > 2 # problem?

# modello con interazione
vif(fit2) # variance inflation factors
sqrt(vif(fit2)) > 2 # problem?

Turista fai da te? No Pre-regression test? Ahi! Ahi! Ahi! Ahi!

Prima di programmare una regressione lineare è utile verificare se sono soddisfatte le assunzioni a priori per una regressione lineare. In R è disponibile un pacchetto che esegue l’intero ventaglio di valutazioni. Il pacchetto si chiama gvlma, che va preventivamente installato. Questo è il codice di valutazione:

library(gvlma)

# modello additivo
gvmodel1 <- gvlma(fit1)
summary(gvmodel1)

# modello con interazione
gvmodel2 <- gvlma(fit2)
summary(gvmodel2)

Il pacchetto produce anche un insieme di plot che definiscono la validità del modello.

# esempio per il modello addittivo
plot(gvmodel1)

Regression diagnostics with gvlma

Alcuni di questi plot sono già riportati dalle diagnostiche della regressione già viste, e queste diagnostiche sono già state spiegate in un post precedente. Altri sono nuovi.

Cosa significano i nuovi plot? Come direbbe un pirata in uno di quei magnifici libri ottocenteschi che nessuno scrive più (Treasure Island, ad esempio, da leggere rigorosamente nell’edizione integrale, senza tagli): “Ch’io possa essere dannato se lo so!...”

That's all, folks!

Antonello Preti

Stay Tuned for our next episode!

Usare i gRafici in R

Tavolozza

Predicare ai convertiti. Così dicono gli Americani.

Noi diremmo “sfondare una porta aperta”.

L’argomento è stato già trattato più volte. Comunque, repetita iuvant.

Uno dei vantaggi di R, e non dei minori, è la sua estrema ricchezza e flessibilità sul piano grafico. Nel suo formato base R ha un sofisticato sistema di rappresentazione grafica dei dati, che consente di avere una rappresentazione dei dati utile alla comprensione delle relazioni sottese.

Questa disponibilità grafica è ulteriormente arricchita dagli sviluppatori, che non mancano di dedicare degli script alla rappresentazione grafica dei dati secondo l’analisi proposta nel pacchetto cui contribuiscono.

Alcuni pacchetti sono stati sviluppati in modo specifico per rendere ancora più versatili le capacità grafiche di R. Alcuni di questi pacchetti gestiscono le forme geometriche, altri sono dedicati alle palettes, le disposizioni dei colori in un grafico, una utilità che chi deve creare grafici per lavoro è in grado di apprezzare con gratitudine.

L’importanza della rappresentazione grafica dei dati è testimoniata dalla frequenza con la quale sono scaricati i pacchetti specificamente dedicati.

Tra i pacchetti di R maggiormente scaricati figurano: ggplot2, colorspace, RColorBrewer, Scales, lattice.

Un passo indietro.

La rappresentazione grafica dei dati è importante, e permette di raggiungere una intuitiva e più immediata consapevolezza delle relazioni statistiche tra le variabili.

Usare R facilmente conduce alla scoperta di questa ovvietà.

Appena raggiunto il satori la prima reazione è quella di incazzarsi come una iena del Serengeti (copyleft, Paolo Villaggio, per i più distratti).

Se il programma statistico su cui avete studiato avesse avuto una migliore interfaccia grafica – vi verrà spontaneo pensare – avreste compreso meglio e più in fretta le basi della statistica. Ma la maggior parte dei programmi statistici commerciali concentra il focus sull’ottenimento della preziosa p < .05, e molte cose sono ad essa sacrificate.

Di seguito, facendo riferimento a dataset disponibili in R, riassumo qualche giochino utile; potrete copiare e incollare i comandi direttamente in R.

Caricate il pacchetto (già disponibile in R) e usate attach per rendere il dataset immediatamente accessibile. Dopo aver usato attach, potrete chiamare le variabili direttamente con il loro nome, risparmiandovi la seccatura di dovere specificare di continuo iris$variabile.

data(iris)
attach(iris)

Ora ottenete una descrizione del dataset:

dim(iris)
head(iris)

La funzione pairs vi consente di ottenere una matrice di correlazione delle variabili del dataset. Ogni variabile è posta in correlazione con le altre: intuitiva la comprensione di come la matrice deve essere letta. Sconsigliabile usarla con dataset con più di 10 variabili; nel caso, scomponete il dataset in vari insiemi. Come? Ripassate la lezione introduttiva su R… è possibile “forzare” pairs a fare cose parecchio interessanti date un’occhiata al manuale online.

pairs(iris)

pairs

Volete sapere se esiste una relazione tra lunghezza dei petali e lunghezza dei sepali in un fiore? Usate plot, che utilizza il formato plot(x,y): prima la variabile che va sull’asse delle x, poi quella che va sull’asse delle y.

plot(Petal.Length, Sepal.Length)

Un metodo alternativo per specificare l’ordine delle variabili in un grafico è usare il formato plot(y~x): la variabile che precede la tilde (il simbolo ~) è la variabile dipendente, che verrà posta sull’asse y, mentre la variabile inserita dopo la tilde è la variabile indipendente, che verrà posta sull’asse delle x.

plot(Sepal.Length ~ Petal.Length, data = iris)

Sembra proprio che esista una relazione lineare, con varianza spiegata pari al 75.8%. Volete esserne sicuri? Provate una regressione lineare:

fit  <- lm(Sepal.Length ~ Petal.Length, data = iris)
summary(fit)

Aggiungete la retta di regressione al grafico:

abline(fit)

Aggiungete la curva della local regression, una regressione non-parametrica. Cosa è lowess? Studiatevelo, vi tornerà utile!

fit2 <- lowess(Sepal.Length ~ Petal.Length)
lines(fit2, col = "red")

plot

Non male, la curva della local regression non si discosta troppo dalla retta della regressione lineare.

Volete l’intervallo di confidenza della retta della regressione lineare?

IntConf <- predict(fit, interval="confidence")

Nella matrice risultante, nelle colonne 2 e 3 si trovano i valori rispettivamente lower e upper dell’intervallo di confidenza della predizione. Le due linee dell’intervallo di confidenza vanno centrate sull’asse delle x (dove è il Petal.Length).

plot(Sepal.Length ~ Petal.Length)
abline(fit)
lines(Petal.Length, IntConf[,2], lty=3, col="blue")
lines(Petal.Length, IntConf[,3], lty=3, col="blue")

lines

Volete essere sicuri che la vostra regressione sia valida? Controllate le diagnostiche del modello. Sono quattro: nelle prime due (in alto) i residui dovrebbero avere una distribuzione casuale e la curva che li attraversa dovrebbe essere più o meno orizzontale; nel QQplot i vostri residui dovrebbero cadere sulla retta o non discostarsene troppo; infine, il quarto plot mostra i casi che violano la distanza di Cook, che dovrebbero essere pochi (sono numerati).

par(mfcol = c (2,2))
plot(fit)

La funzione par(mfcol = c(2,2)) serve a ripartire la lavagna grafica in quattro riquadri; se vi servono più riquadri appaiati, semplicemente variate gli input in c(num righe, num colonne).

plot(fit)

Nel complesso, non male.

Una migliore descrizione della violazione della distanza di Cook si ottiene con un plot del pacchetto car (da scaricare, se non lo avete), dove car sta per Companion to Applied Regression; creiamo un influence plot con punti proporzionali alla distanza di Cook (casi rilevanti identificati in rosso)

library(car)
layout(1)
influencePlot(fit, id.n=3,id.col="red")

Influence plot

Dato che prima, con il comando par(), avevamo ripartito la lavagna grafica in quattro riquadri, con layout(1) la riportiamo a un solo riquadro.

Sull’argomento delle regression diagnostics, utile riassunto qui, qui e qui. John Fox è il masterchef dell’argomento.

Volete avere una distribuzione della lunghezza dei petali nelle tre specie di fiori del dataframe? Perché non usare un bel density plot, sempre dal pacchetto car? La variabile di raggruppamento deve essere trattata come fattore e Species è già di classe “factor”. Non ci credete? Chiedetelo:

is.factor(Species)

La risposta è TRUE, mentre per is.numeric(Species) la risposta è FALSE. Se la vostra variabile non è un fattore, rendetela tale così (meglio assegnare un nuovo nome alla variabile, non si sa mai):

nuovo_nome_vostra_variabile <- as.factor(nome_vostra_variabile)

Ora possiamo andare sicuri a costruire i density plot:

par(mfcol = c(1,2))
densityPlot(Petal.Length~Species,data=iris)
densityPlot(Sepal.Length~Species,data=iris)

density plot

Bella differenza, eh?

Altra rappresentazione, ottenibile con un altro pacchetto, graphics (già installato e attivo nella distribuzione base di R) è il conditioning plot:

coplot(Sepal.Length ~ Petal.Length | Species, data = iris, panel = panel.smooth, rows = 1)

coplot

Il risultato suggerisce che per le tre specie esiste una relazione tra lunghezza dei petali e dei sepali, meno marcata per la specie setosa.

E adesso provate voi: caricare il dataframe mtcars:

data(mtcars)

Rapido check alle variabili:

head(mtcars)

E ora create un plot della relazione tra la variabile mpg (Miles/(US) gallon, miglia con un litro) e quella wt (weight, peso della macchina):

plot(mpg,wt)

Aha! messaggio di errore:

Errore in plot(disp, wt) : oggetto “mpg” non trovato

Capito perché è utile usare attach? Con attach potete nominare direttamente le variabili:

attach(mtcars)
plot(mpg,wt)

Ora funziona. Bene, proseguite con l’esercizio oppure variate a piacimento. Sbizzarritevi chiamando la lista dei dataset disponibili in R, digitando:

library(help = "datasets")

Per caricare il dataset, basta chiamarlo con: data(nome_del_dataset).

That’s all, folks!

Antonello Preti

Stay Tuned for our next episode!

ApRire un file SPSS oppure Excel in R

R è un sistema potente per l’analisi statistica e per la rappresentazione grafica. Tuttavia, non è esattamente user-friendly per l’archiviazione dei dati. Per un bel po’ ancora i vostri dati saranno archiviati in Excel o in SPSS o simili.

Da SPSS a R

Come aprire in R un file archiviato in altro sistema? Ci sono librerie (packages) come foreign che consentono di effettuare l’operazione. La libreria foreign è già installata nella distribuzione base del sistema: è sufficiente attivarla con la funzione library.

> library(foreign)

A questo punto potete aprire il vostro file, se sapete dove si trova…

Il metodo più semplice per localizzare un file (lo so, lo so, settare la directory, ma io sono uno di modi spicci) è digitare:

> file.choose()

Si aprirà una finestra come quella di Windows per l’accesso ai file, cercate il vostro file nelle cartella dove lo avete archiviato, e via. R vi restituisce la path, cioè il percorso:

> "C:\\Users\\Antonio\\Desktop\\Traferire in MAC\\InsulaR\\DatabaseEsempio.sav"

A questo punto potete dare il comando di lettura file in SPSS da foreign, specificando la path che conduce al file (sì, avete capito, copia e incolla tra virgolette):

dataset = read.spss("C:\\Users\\Antonio\\Desktop\\Traferire in MAC\\InsulaR\\DatabaseEsempio.sav", to.data.frame=TRUE)

Volete risparmiarvi il copia e incolla? Assegnate l’operazione file.choose() a un oggetto, chiamiamolo db (abbreviazione per database):

> db = file.choose()

Come prima, avete ottenuto il percorso al file, solo che questa volta il sistema non ve lo visualizza perché voi l’avete assegnato all’oggetto db. Quindi, ora l’oggetto db contiene una stringa di caratteri che indentifica il percorso che R dovrà seguire per recuperare il file.

Pronti via?

dataset = read.spss(db, to.data.frame=TRUE)

Il comando read.spss() legge il dataset in formato sav; bisogna stare attenti però a specificare come TRUE l’argomento to.data.frame, il quale richiede alla funzione di disporre i dati all’interno di un data frame (ovvero una tabella di dati).

Yolo, man. Un altro metodo molto semplice di aprire un file SPSS in R è salvare il file in un formato che R gestisce molto bene: il .dat (delimitato da tabulazione). Salvate il vostro file SPSS in .dat e operate come prima, andando alla ricerca del file con la funzione file.choose() e assegnando la stringa risultante a un oggetto.

La funzione per leggere il file stavolta è read.table(). Prestate molta attenzione ai dati mancanti: se ce ne sono, bisogna dire a R quale codice li identifica (es. 99), specificando un valore per l’argomento na.strings.

Avete il vostro file in .dat?

db = file.choose()
dataset = read.table(db, header = TRUE)

l’argomento header = TRUE serve per specificare che la prima riga del file contiene i nomi delle variabili per cui questi valori non sono da interpretare come dati.

Da Excel a R

Per importare in R un file da Excel ci sono molti metodi diversi e purtroppo nessuno è esente da problemi. Uno l’abbiamo già visto in questo post, ma ce ne sono di più semplici. Il metodo più immediato richiede due passaggi: 1) esportare il foglio di dati Excel in file formato testuale, 2) leggere il nuovo file attraverso le funzioni R. Per la verità, questo metodo è valido non solo per il programma Excel di Microsoft ma anche per l’applicativo Calc di LibreOffice/OpenOffice

Per prima cosa, utilizzando il menù di Excel, si deve salvare il file in formato csv. Excel fornisce diverse opzioni: CSV (delimitato da separatore di elenco), CSV (Macintosh) e CSV (MS-DOS). Nella pratica, un’opzione vale l’altra.

Il risultato sarà un file testuale che contiene tante righe quante erano le righe dell’originale Excel, con i dati separati da un punto e virgola all’interno di ogni riga. La virgola, invece, viene utilizzata come separatore decimale. Questo formato è quello che R chiama csv2, che è lo standard più utilizzato nell’Europa mediterranea. Il resto del mondo, invece, utilizza non il punto e virgola ma la virgola come separatore di campo e il punto come separatore decimale. Utilizzare il csv è molto comodo, non solo perché Excel lo riconosce e lo apre automaticamente ma anche perché è il formato standard con il quale vengono condivisi il dati.

Una volta esportato il file in csv, lo potete importare in R utilizzando i comandi:

db = file.choose()
dataset = read.csv2(db)

Nel lanciare read.csv2(), bisogna prendere alcune precauzioni. Prima di tutto, per default la funzione assume che la prima riga contiene i nomi che dovranno essere assegnati alle colonne (header = TRUE). Se invece la prima riga non contiene le etichette da utilizzare come nomi di colonna, bisognerà specificare header = FALSE. In secondo luogo, possiamo anche indicare alla funzione qual è il codice usato per codificare i dati mancanti; useremo na.strings = "" se le celle vuote corrispondono a dati mancanti, ma potremo anche utilizzare altri valori come ad esempio na.strings = 99.

Se il separatore di campo del file csv non è il punto e virgola ma la virgola e il separatore decimale è il punto, il comando da usare sarà read.csv(). Da notare che le due funzioni read.csv() e read.csv2() altro non sono che dei casi specifici di read.table() e, a dirla tutta, si potrebbe anche usare la stessa funzione read.table() per leggere un file csv:

dataset = read.table(db, sep=";", dec=",", header=TRUE)

Verificare che la lettura sia andata a buon fine

Una volta importato un file, è sempre bene fare dei controlli per verificare che la lettura sia andata a buon fine.

Per controllare le dimensioni del vostro database, usate la funzione dim. Vi restituirà due numeri, il primo si riferisce ai casi (le righe del vostro database), mentre il secondo numero è quello delle variabili (le colonne del vostro database). Facciamo un esempio con un database che avete chiamato MieiDati:

> dim(MieiDati)

Inoltre, può essere utile visualizzare un’anteprima dei dati. Per ispezionare le prime 6 righe del vostro database, usate la funzione head:

> head(MieiDati)

Per ispezionare le ultime 6 righe del vostro database, usate la funzione tail:

> tail(MieiDati)

Per ispezionare la struttura del database, usate la funzione str:

> str(MieiDati)

Volete visualizzare tutta la matrice di dati? Se la tabella di dati è grossa, vi conviene usare la funzione View, oppure fix che vi consentirà anche di modificare manualmente il contenuto delle celle:

> View(MieiDati)
> fix(MieiDati)

Antonello Preti