Archivio tag: plot

Visualizzare la relazione fra due variabili likert

Lavorando nel campo della psicologia, spesso mi trovo ad avere a che fare con variabili che derivano da riposte a questionari fornite utilizzando scale di tipo Likert. Si tratta di variabili che possono assumere un numero molto limitato di modalità, comunemente da tre a cinque. Quando le categorie di riposta sono almeno cinque, spesso queste variabili vengono considerate come se fossero continue e di conseguenza vengono analizzate usando indicatori e modelli statistici pensati appunto per variabili continue.

Con questo post non voglio entrare nel merito della correttezza di queste scelte (condivisibili o meno, a seconda dei casi), ma concentrarmi sul modo di utilizzare al meglio gli strumenti statistici, principalmente i grafici.

Di seguito viene costruito il data frame likert che contiene due variabili: item1 e item2, che riportano le risposte (espresse su scala Likert a cinque punti) di sessanta ipotetiche persone a due ipotetiche domande di un ipotetico questionario.

likert <- data.frame(
    item1 = c(5,3,5,5,4,1,2,3,3,5,5,2,1,3,3,3,4,2,4,1,
              3,4,4,5,3,4,5,2,1,4,3,2,3,2,4,5,5,2,5,3,4,
              4,3,2,1,1,1,5,2,3,1,1,2,2,2,3,4,2,4,3),
    item2 = c(5,2,5,4,4,1,3,3,1,3,3,2,1,3,1,3,1,2,4,1,
              3,4,2,4,4,5,5,2,1,1,1,4,4,2,4,4,5,2,5,4,5,
              3,3,1,2,2,1,5,5,5,3,4,1,1,3,2,4,2,3,2)
)
&#91;/code&#93;

Per studiare la relazione tra due variabili di questo tipo, quello che comunemente viene fatto è calcolare l'indice di correlazione lineare. Possiamo usare l'indice di Pearson oppure l'indice di Spearman, basato sui ranghi.

Attraverso il comando <b>cor</b> di R possiamo calcolare entrambi gli indici; se non viene specificato nulla nell'argomento <i>method</i>, verrà calcolata la correlazione di Pearson:

[code language="R"]
> with(likert, cor(item1, item2))
[1] 0.5901713

Se invece vogliamo calcolare l’indice di Spearman, dobbiamo esplicitare questa richiesta nell’argomento method:

> with(likert, cor(item1, item2, method="spearman"))
[1] 0.5864863

Entrambi gli indici evidenziano un’ottima correlazione fra le due variabili (“ottima” almeno per gli standard in psicologia): si sfiora lo 0.6, che è un valore abbastanza elevato.

Quello che – ahimè – non fa quasi nessuno, è visualizzare la relazione tra le due variabili. I grafici hanno sempre tanto da raccontare e talvolta è proprio dalle visualizzazioni che emergono gli aspetti più interessanti. Proviamo quindi a creare uno scatterplot per visualizzare la relazione fra item1 e item2. Utilizziamo il comando plot, aumentando la dimensione dei punti sfruttando l’argomento cex:

with(likert, plot(item1, item2, cex=2))

Scatter-plot fra item1 e item2

Beh, che ve ne pare? Si tratta di una buona visualizzazione? Secondo me, no.

A vedere questo grafico a me sorgono molte perplessità. Abbiamo appena detto che la relazione lineare tra le due variabili è buona, ma dal grafico proprio non si direbbe: i punti sono sparpagliati un po’ ovunque e le due variabili sembrano tutto fuorché correlate.

Il problema di questo grafico è che le variabili possono assumere pochi valori (da 1 a 5), per cui moltissime risposte si sovrappongono. Ognuno di quei pallini in realtà ha una densità, perché su ognuno di essi si sovrappongono le risposte di più individui. Osservando bene l’immagine, infatti, possiamo notare che ci sono pallini il cui contorno è più scuro di altri; ebbene, nelle coordinate più scure si concentrano le risposte di più persone.

Il numero di osservazioni presente in ogni coordinata è un dato fondamentale per comprendere la relazione tra due variabili che assumono un numero ridotto di modalità, ma nell’immagine qui sopra questa informazione non è ben rappresentata.

La “densità” può essere calcolata semplicemente contando il numero di risposte che occorrono per ognuno degli incroci dei valori delle due variabili, ovvero costruendo una tabella di frequenza a doppia entrata:

> tab <- with(likert, table(item1, item2))
> tab
     item2
item1 1 2 3 4 5
    1 5 2 1 1 0
    2 3 6 2 1 1
    3 3 3 5 3 1
    4 2 1 2 5 2
    5 0 0 2 3 6

Osservando la tabella qui sopra possiamo notare come le frequenze maggiori siano collocate sulla diagonale, fenomeno che supporta la presenza di una relazione lineare e che giustifica valori di correlazione così elevati. Ma come fare per considerare questa informazione nel grafico?

Adesso vi proporrò due alternative; entrambe richiedono che la tabella di frequenza che abbiamo appena costruito venga convertita in un oggetto di tipo data.frame:

> tab <- as.data.frame(tab)
> head(tab)
  item1 item2 Freq
1     1     1    5
2     2     1    3
3     3     1    3
4     4     1    2
5     5     1    0
6     1     2    2

Grafico a bolle (bubble chart)


Quello che manca al grafico costruito poco sopra è l’informazione sul numero di osservazioni in ogni coordinata. Il modo più semplice di considerare quella che è a tutti gli effetti una terza variabile è fare in modo che il diametro di ogni punto dipenda dalla frequenza. Verrà creato così un grafico nel quale saranno presenti punti più grandi e punti più piccoli: più grande sarà il punto, maggiore sarà la concentrazione di dati.

Per realizzare questo grafico possiamo usare il comando symbols, specificando nell’argomento inches un’unità di misura per calibrare la dimensione dei punti.

with(tab, symbols(item1, item2, Freq, inches=0.6))

Bubble-plot fra item1 e item2

Il risultato conferma che la densità è maggiore sulla diagonale e va diminuendo con l’allontanarsi da questa. I dati quindi non sono sparpagliati in maniera casuale come poteva sembrare nel primo grafico, ma la relazione fra item1 e item2 segue un andamento ben preciso (nello specifico lineare).

Grafico a mattonelle (tile plot)


Un altro tipo di grafico che ci viene in aiuto è il tile plot. In questo tipo di visualizzazione, i dati sono rappresentati attraverso delle mattonelle che vengono colorate a seconda del valore assunto da una terza variabile. Per costruire il tileplot dobbiamo installare e attivare il pacchetto ggplot2 (del quale abbiamo già parlato).

Allo strato di base costruito con il comando ggplot dobbiamo aggiungere un livello creato con geom_tile e un gradiente di colore con scale_fill_gradient; infine, possiamo specificare un tema (io ho scelto theme_bw):

library(ggplot2)
ggplot(data=tab, aes(x=item1, y=item2, fill=Freq)) +
    geom_tile(colour="white") +
    scale_fill_gradient(low="white", high="red3") +
    theme_bw()

Tile-plot fra item1 e item2

Il grafico a mattonelle è forse quello più accattivante, ma probabilmente anche il più complicato da costruire, visto che richiede l’uso di ggplot2.

E voi, quali soluzioni prediligete in questi casi?

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!