Grafici di dispersione (scatterplot)

I grafici di dispersione (scatter plot o scatterplot) sono lo strumento base per visualizzare su un sistema di coordinate cartesiane le coppie di valori x,y di due variabili numeriche nel caso di distribuzioni bivariate.

Prendiamo come esempio i dati ematologici relativi agli eritrociti (o globuli rossi), la cui funzione essenziale come noto è trasportare l'ossigeno dai polmoni a tutti gli altri organi e tessuti, e per i quali valgono in media i seguenti dati (arrotondati per semplicità) :

→ gli eritrociti hanno una concentrazione nel sangue all'incirca di 5·10¹²/L (ovvero ve ne sono all'incirca cinquemila miliardi in un litro di sangue);

→ un eritrocito ha un volume all'incirca di 90 fL (un femtolitro corrisponde a 10^-15 litri quindi un globulo rosso ha un volume all'incirca di 0.000 000 000 000 090 litri).

Da questi dati si ricava che la percentuale del volume del sangue occupata dagli eritrociti (detta valore ematòcrito, e uguale alla concentrazione degli eritrociti moltiplicata per il volume di un eritrocito) è il 45% (la restante parte, esclusi i globuli bianchi e le piastrine, che occupano un volume irrilevante, è liquida ed è rappresentata dal plasma). Ma non solo. Da questi dati si ricava anche che il valore ematòcrito dipende in modo lineare dalla concentrazione degli eritrociti. Possiamo quindi dare un senso al prossimo script, che traccia i grafici di dispersione di queste due grandezze.

I dati sono contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [1].

Copiate lo script e incollatelo nella Console di R e premete ↵ Invio:

# GRAFICI DI DISPERSIONE (scatterplot)

#

library(DAAG) # carica il pacchetto che include il set di dati ais

str(ais) # struttura di ais

attach(ais) # consente di impiegare direttamente i nomi delle variabili

#

windows() # apre e inizializza una nuova finestra grafica

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

plot(rcc, hc, xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Grafico di dispersione") # grafico di dispersione semplice

#

plot(rcc, hc, pch=2, col="goldenrod", xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Cambia simbolo e colore") # cambia simbolo e colore

#

plot(rcc, hc, pch=1, col="black", abline(lm(hc~rcc), col="red", lty=1, lwd=1), xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Retta di regressione") # grafico con retta di regressione

#

plot(rcc, hc, xlim=c(3,7), ylim=c(30,70), xaxp=c(3, 7, 4), yaxp=c(30, 70, 4), pch=20, col="black", abline(lm(hc~rcc), col="black", lty=1, lwd=1), xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Ridimensiona gli assi") # ridimensiona gli assi

#

detach(ais) # termina l'impiego diretto dei nomi delle variabili 

#

Dopo avere caricato il pacchetto che include il set di dati ais e dopo averne mostrato la sua struttura con str(ais), viene impiegata la funzione attach(ais) al fine di consente, nelle funzioni successive, di impiegare direttamente i nomi delle variabili incluse nel set di dati.

Per semplificare la presentazione grafica, con par(mfrow=c(2,2)) la finestra viene suddivisa in quattro quadranti, uno per ciascuno dei grafici che verranno preparati.

Il primo grafico, in alto a sinistra, lascia per la maggior parte degli argomenti della funzione plot() i valori di default, specificando solamente quelli essenziali [2]:

→ rcc la variabile da riportare sull'asse delle ascisse x;

→ hc la variabile da riportare sull'asse delle ordinate y;

→ xlab="" l'etichetta da riportare per l'asse delle x;

→ ylab="" l'etichetta riportata per l'asse delle y;

→ main="" il titolo del grafico che compare in alto.

Nel secondo grafico, in alto a destra, sono aggiunti gli argomenti:

→ pch=2 che specifica per i dati l'impiego di un simbolo (triangolo) differente da quello di default (che è il cerchio) [3];

→ col="goldenrod" che specifica per il simbolo dei dati l'impiego di un colore diverso dal nero [4].

Nel terzo grafico, in basso a sinistra, viene riportata una retta mediante la funzione abline() che impiega come argomenti:

→ l'equazione della retta, che altro non è se non la retta di regressione calcolata con la funzione lm() che a sua volta impiega come argomenti hc~rcc per specificare rispettivamente la variabile in ordinate e la variabile in ascisse;

→ il colore della retta (col="red");

→ il tipo di linea da impiegare (lty=1), uno dei sei possibili con i valori da 1 a 6 [5];

→ lo spessore della linea (lwd=1) espresso come multiplo dello spessore di default (1 uguale spessore di default, 2 spessore doppio, eccetera).

Nel quarto grafico, in basso a destra, compaiono quattro nuovi argomenti, che consentono di gestire in modo personalizzabile a piacere le scale degli assi e la loro suddivisione:

→ xlim = c(3,7) indica limite inferiore e limite superiore della scala applicata all'asse delle ascisse, trattandosi di valori di concentrazione degli eritrociti espressi in 10¹²/L questo significa che l'asse delle ascisse riporterà i valori compresi tra 3·10¹²/L e 7·10¹²/L;

→ ylim = c(30,70) indica limite inferiore e limite superiore della scala applicata all'asse delle ordinate, che, trattandosi di valori di percentuale, varieranno quindi tra 30% e 70%;

→ xaxp = c(3, 7, 4) dice che sull'asse delle ascisse la scala che va da 3·10¹²/L a 7·10¹²/L deve essere suddivisa in quattro intervalli, quindi sull'asse delle ascisse compariranno cinque tacche con i valori 3, 4, 50, 6, 7;

→ yaxp = c(30,70,4) dice che sull'asse delle ordinate la scala che va da 30% e 70% deve essere suddivisa in quattro intervalli, quindi sull'asse delle ordinate compariranno cinque tacche con i valori 30, 40, 50, 60, 70.

Lo script si chiude con la funzione detach() che pone termine all'impiego diretto dei nomi delle variabili. In altre parole la variabile rcc ora diventa ais$rcc, la variabile hc ora diventa ais$hc e così via.

I grafici di dispersione riportati qui sopra sono interessanti in quanto sono stati realizzati impiegando esclusivamente le funzioni base di R, tuttavia si può fare di meglio come nell'esempio che segue, che per tracciare lo scatterplot impiega le funzioni del pacchetto ggplot2.

Copiate lo script e incollatelo nella Console di R e premete ↵ Invio:

# GRAFICO DI DISPERSIONE con ggplot

# colori differenziati in base a un fattore

#

library(DAAG) # carica il pacchetto che include il set di dati ais

library(ggplot2) # carica il pacchetto per la grafica

#

windows() # apre e inizializza una nuova finestra grafica

ggplot(ais, aes(x=rcc, y=hc, color=sport)) + geom_point(size=4) + theme_minimal() # traccia il grafico

#

Ed ecco il grafico, nel quale è interessante notare come la funzione ggplot() consente di attribuire ai punti un differente colore (color=) sulla base dei valori assunti dal fattore (sport) – lo sport praticato – evidenziando così nei dati i corrispondenti sottoinsiemi. Altri esempi che illustrano le varie possibilità di realizzare scatterplot con ggplot2 sono riportati in un post al quale rimando [6].

Vediamo ora due grafici di dispersione multipli che aiutano ad illustrare le relazioni tra le variabili riportando tutti i possibili grafici di dispersione singoli che risultano incrociando i dati ematologici contenuti nelle seguenti variabili del set di dati ais: rcc (eritrociti), wcc (leucociti), hc (ematòcrito), hg (emoglobina) e ferr (ferritina). 

I grafici sono realizzati mediante i pacchetti aggiuntivi car e gclus: se non l'avete già fatto, dovete scaricarli e installarli dal CRAN.

Copiate e incollate nella Console di R lo script e premete ↵ Invio:

# GRAFICI DI DISPERSIONE MULTIPLI 

#

library(DAAG) # carica il pacchetto che include il set di dati ais

library(car) # carica il pacchetto per la grafica

windows() # apre e inizializza una nuova finestra grafica

#

# grafici di dispersione con istogrammi delle variabili

#

scatterplotMatrix(~rcc+wcc+hc+hg+ferr, col="black", pch=20, regLine = list(method=lm, lty=1, lwd=2, col="chartreuse3"), smooth=FALSE, diagonal=list(method ="histogram", breaks="FD"), main="Matrice di dispersione con rette di regressione", data=ais) # mostra il grafico

#

L'argomento ~rcc+wcc+hc+hg+ferr definisce le variabili da rappresentare, e può essere facilmente riscritto per cambiare la rappresentazione a piacere. Ad esempio ponendo questo argomento uguale a ~rcc+wcc+hc+hg+ferr+bmi+ssf+pcBfat+lbm+ht+wt si possono rappresentare in uno stesso grafico tutte le variabili del set di dati ais.

Per la retta di regressione sono impiegati una linea continua (lty=1), uno spessore doppio (lwd=2) e un colore (col="chartreuse3") diverso dal nero, mentre l'argomento smooth=FALSE fa si che non venga aggiunta ai grafici la rappresentazione della funzione non lineare prevista di default nel pacchetto car.

Infine l'argomento diagonal=list() indica la rappresentazione della variabile riportata nella diagonale, che qui è l'istogramma (method ="histogram") con il numero di classi (breaks="FD") calcolato mediante la regola di Freedman-Diaconis. L'istogramma può essere sostituito con altri tipi di rappresentazione, l'elenco delle espressioni che è possibile impiegare in alternativa per questo argomento (provatele) è:

→ diagonal=list(method="density",  bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE)

→ diagonal=list(method="boxplot")

→ diagonal=list(method="qqplot")

→ diagonal=list(method="oned")

Copiate e incollate nella Console di R lo script e premete ↵ Invio:

# GRAFICI DI DISPERSIONE MULTIPLI ordinati in base al valore di r

#

library(DAAG) # carica il pacchetto incluso il set di dati ais

library(gclus) # carica il pacchetto per la grafica

windows() # apre e inizializza una nuova finestra grafica

mydata <- ais[c(1,2,3,4,5)] # tabella con i dati delle variabili delle colonne 1,2,3,4,5

mydata.r <- abs(cor(mydata)) # calcola la correlazione

mydata.col <- dmat.color(mydata.r) # applica i colori

mydata.o <- order.single(mydata.r) # ordina le variabili, le meglio correlate più vicine alla diagonale

cpairs(mydata, mydata.o, panel.colors=mydata.col, gap=.5, main="Variabili meglio correlate vicine alla diagonale") # mostra il grafico

#

Viene realizzato un grafico di dispersione multiplo che incrocia ancora i dati ematologici degli atleti contenuti nelle cinque variabili rcc (eritrociti), wcc (leucociti), hc (ematòcrito), hg (emoglobina) e ferr (ferritina) ma con due novità.

La prima è che le variabili ora sono individuate non mediante il loro nome bensì indicando, con l'espressione ais[c(1,2,3,4,5)], il set di dati che le contiene e il numero della colonna corrispondente a ciascuna variabile

La seconda novità è che impiegando la funzione cpairs() [7] i singoli grafici di dispersione vengono ordinati in base al valore del coefficiente di correlazione r di Pearson, in modo che le variabili meglio correlate, quelle per le quali la retta di regressione pare essere un'approssimazione ragionevole, siano le più vicine alla diagonale ed evidenziate in colore, aumentando così la chiarezza della rappresentazione grafica.

Potete riutilizzare facilmente lo script sostituendo all'oggetto ais l'oggetto contenente i vostri dati, opportunamente strutturati, e modificando conseguentemente i nomi della variabili. Per una guida rapida all'importazione dei dati potete consultare i post:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestione dei dati mancanti

----------

[1] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG.

[2] Per la documentazione sugli argomenti della funzione plot() digitare help(plot) nella Console di R, e in aggiunta consultare la documentazione sugli argomenti della funzione par() con help(par).

[3] Per i diversi simboli che è possibile impiegare per tracciare i punti vedere il post I simboli dei punti di R.

[4] Per i colori che è possibile impiegare vedere il post La tabella dei colori di R.

[5] Per le linee che è possibile impiegare vedere il post Gli stili delle linee di R.

[6] Vedere il post Grafici di dispersione (scatterplot) con ggplot.

[7] Digitate help(cpairs) nella Console di R per la documentazione della funzione cpairs() ovvero consultate il manuale di riferimento del pacchetto gclus su: Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html

Regressione lineare semplice parametrica

La regressione lineare tout court, quella illustrata in tutti i testi di statistica e qui sviluppata con R, è:

→ una regressione lineare semplice (in contrapposizione alla regressione lineare multipla);

→ parametrica (in contrapposizione alla regressione lineare non parametrica);

→ x variabile indipendente (in contrapposizione ad alternative che non prevedono questo assunti).

Dato che la denominazione razionale, quella completa, diverrebbe chilometrica, e dato che l'espressione regressione lineare è però troppo generica, laddove è opportuno semplificare impiegherò la denominazione regressione lineare ordinaria, con riferimento al fatto che è quella non solo dovunque illustrata ma anche più ampiamente impiegata.

Gli assunti alla base del modello matematico-statistico implicano una serie di requisiti che devono essere soddisfatti dai dati, che sono ben chiariti ad esempio in Marubini [1] e in Snedecor [2], ma anche online [3]. Per il caso speciale, ma non infrequente, nel quale nessuna delle due variabili confrontate abbia i requisiti richiesti per essere considerata come variabile indipendente, vedere anche [4].

Qualora invece sia necessario l'impiego di "metodi robusti" di regressione - cioè di metodi che risentono poco della eventuale presenza di dati apparentemente anomali ma che non possono essere esclusi a priori dal calcolo della regressione - è possibile impiegare un metodo non parametrico [5].

Tornando alla nostra regressione lineare semplice parametrica x variabile indipendente, essendo x la variabile indipendente posta sull'asse delle ascisse, y la variabile dipendente posta sull'asse delle ordinate, ed essendo soddisfatti i requisiti previsti per i dati, il metodo dei minimi quadrati consente di calcolare l'intercetta a e il coefficiente angolare b della retta di regressione di equazione

y = a + b∙x

che meglio approssima la distribuzione dei dati sperimentali.

Il calcolo dell'equazione della retta di regressione viene effettuato mediante la funzione lm(), che può essere applicata anche al caso di più variabili indipendenti, consentendo quindi il calcolo della regressione lineare multipla [4]. Ma la cosa più interessante è che R, oltre al calcolo dell'equazione della retta di regressione, un calcolo di per sé semplice, fornisce una serie molto interessante di strumenti che consentono di valutare quanto i dati soddisfano i requisiti richiesti, o in altre parole di valutare se la regressione lineare descrive in modo adeguato la relazione tra le due variabili.

Qui impieghiamo come esempi due diversi set di dati, per uno solo dei quali, come vedremo, la regressione lineare fornisce risultati adeguati.

Il primo è il set di dati ais, nel quale prendiamo in considerazione la concentrazione degli eritrociti (espressa in 10¹²/L) e il valore ematòcrito (espresso in %), due variabili che, analizzate mediante i coefficienti di correlazione parametrici e non parametrici, hanno mostrato di essere correlate, e che, all'ispezione visiva di una serie di grafici di dispersione (xy) realizzati mediante i correlogrammi, hanno mostrato valori ben allineati su una possibile retta.

Per procedere dovete scaricare e installare dal CRAN i pacchetti aggiuntivi gvlma e car, oltre al pacchetto aggiuntivo DAAG [6].

Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x

#

library(DAAG) # carica il pacchetto che include il set di dati ais

library(gvlma) # carica il pacchetto per la funzione gvlma()

library(car) # carica il pacchetto per analisi ouliers e analisi grafica

str(ais) # mostra la struttura dei dati

#

var_x <- ais$rcc # variabile eritrociti x in ascisse

var_y <- ais$hc # variabile ematòcrito y in ordinate

reglin <- lm(var_y ~ var_x) # calcola intercetta (a) e coefficiente angolare (b)

coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc

confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare

#

# analisi statistica della adeguatezza della regressione lineare

#

summary(reglin) # mostra un riepilogo dei risultati

t.test(residuals(reglin)) # verifica che la media degli errori non sia significativamente diversa da zero

shapiro.test(residuals(reglin)) # verifica la normalità della distribuzione degli errori

summary(gvlma(reglin)) # test globale per l'assunto di linearità

outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)

#

# analisi grafica della adeguatezza della regressione lineare

#

windows() # apre una nuova finestra

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

newx = seq(min(var_x), max(var_x), by = 0.01) # valori della x per i quali calcolare l'intervallo di confidenza

conf_interval <- predict(reglin, newdata=data.frame(var_x=newx), interval="confidence", level = 0.95) # calcola gli intervalli di confidenza

plot(var_x, var_y, xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Regressione lineare y = a + b·x") # grafico dei dati

abline(reglin, col="lightblue") # retta di regressione

lines(newx, conf_interval[,2], col="blue", lty=2) # limite di confidenza inferiore

lines(newx, conf_interval[,3], col="blue", lty=2) # limite di confidenza superiore

#

plot(var_y, residuals(reglin), xlab="Ematòcrito (%)", ylab="Ematòcrito osservato - calcolato (%)", main="Analisi delle differenza residue") # grafico delle differenza tra ematòcrito osservato e ematòcrito calcolato con l'equazione della retta

#

influencePlot(reglin, fill=FALSE, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

#

qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

#

Dopo avere caricato i pacchetti aggiuntivi e i dati, e dopo avere mostrato con str(ais) la struttura di questi ultimi, con var_x <- ais$rcc e con var_y <- ais$hc sono memorizzate negli oggetti var_x e var_y rispettivamente la variabile indipendente x, posta in ascisse, e la variabile dipendente y, posta in ordinate. In questo modo sarà possibile riutilizzare per intero lo script semplicemente sostituendo ad ais$rcct e ad ais$hc i vettori contenenti i propri dati.

Con reglin <- lm(var_y ~ var_x) l'equazione della retta di regressione che esprime la y in funzione della x viene calcolata e memorizzata nell'oggetto reglin, che a questo punto diventa l'argomento chiave, l'argomento contenente i dati in ingresso impiegati nelle funzioni successive [7].

Con le funzioni coefficients() e confint() sono infine mostrati i coefficienti della retta di regressione e i loro intervalli di confidenza al 95%

> coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc 

(Intercept)       var_x 

   8.183033    7.398052 

> confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare 

               2.5 %    97.5 %

(Intercept) 6.173717 10.192350

var_x       6.974206  7.821898

e pertanto questa risulta essere l'equazione della retta di regressione:

hc = 8.183033 + 7.398052 · rcc

A questo punto seguono due blocchi di codice, il primo per effettuare una analisi statistica, e il secondo per effettuare una analisi grafica della regressione, entrambe allo scopo di valutare, come già detto, se la regressione lineare descrive in modo adeguato la relazione tra le due variabili.

Per quanto concerne l'analisi statistica, le principali conclusioni sono quelle tratte con il test globale per l'assunto di linearità effettuato mediante la funzione gvlma(), che conferma il fatto che gli assunti che stanno alla base del modello di regressione lineare sono tutti soddisfatti:

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                      Value p-value                Decision

Global Stat        3.479904  0.4809 Assumptions acceptable.

Skewness           2.233166  0.1351 Assumptions acceptable.

Kurtosis           0.958154  0.3277 Assumptions acceptable.

Link Function      0.005239  0.9423 Assumptions acceptable.

Heteroscedasticity 0.283345  0.5945 Assumptions acceptable.

Il test di Bonferroni non evidenzia dati anomali, ma segnala il dato numero 68 come quello che maggiormente si discosta dai rimanenti:

> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers) 

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

    rstudent unadjusted p-value Bonferonni p

68 -3.569525         0.00044828     0.090552

Le conclusioni dell'analisi grafica confermano, con i valori delle differenze residue dispersi in modo casuale, che la regressione lineare descrive in modo adeguato la relazione tra le due variabili

e confermano anche la presenza di possibili dati anomali con il grafico delle distanze di Cook e con il grafico dei quantili, dati che sono evidenziati sia direttamente nei grafici, sia nei rispettivi riepiloghi riportati nella Console di R:

> influencePlot(reglin, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni 

      StudRes         Hat      CookD

68  -3.569525 0.006301040 0.03815682

78   3.186317 0.009724284 0.04766681

161 -2.179789 0.039758811 0.09655642

166  1.363972 0.099962743 0.10287127

> #  

> qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati 

[1] 68 78

I dati per i quali i metodi di identificazione meglio concordano sono il numero 68 e il numero 78, ma anche i dati numero 161 e numero 166 meritano di essere valutati. Si tratta di casi dei quali sarebbe importante controllare la validità, o che potrebbero essere situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati.

Vediamo ora la stessa identica analisi applicata al set di dati galton [8]. I pacchetti aggiuntivi psychTools, gvlma e car, se non l'avete già fatto, devono essere preventivamente scaricati e installati dal CRAN [9].

Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x

#

library(psychTools) # carica il pacchetto che include il set di dati galton

library(gvlma) # carica il pacchetto per la funzione gvlma()

library(car) # carica il pacchetto per analisi outliers e analisi grafica

str(galton) # mostra la struttura dei dati

#

var_x <- galton$parent # variabile altezza dei padri x in ascisse

var_y <- galton$child # variabile altezza dei figli y in ordinate

reglin <- lm(var_y ~ var_x) # calcola intercetta (a) e coefficiente angolare (b)

coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc

confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare

#

# analisi statistica della adeguatezza della regressione lineare

#

summary(reglin) # mostra un riepilogo dei risultati

t.test(residuals(reglin)) # verifica che la media degli errori non sia significativamente diversa da zero

shapiro.test(residuals(reglin)) # verifica la normalità della distribuzione degli errori

summary(gvlma(reglin)) # test globale per l'assunto di linearità

outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)

#

# analisi grafica della adeguatezza della regressione lineare

#

windows() # apre una nuova finestra

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

newx = seq(min(var_x), max(var_x), by = 0.01) # valori della x per i quali calcolare l'intervallo di confidenza

conf_interval <- predict(reglin, newdata=data.frame(var_x=newx), interval="confidence", level = 0.95) # calcola gli intervalli di confidenza

plot(var_x, var_y, xlab="Altezza dei padri (pollici)", ylab="Altezza dei figli (pollici)", main="Regressione lineare y = a + b·x") # grafico dei dati

abline(reglin, col="lightblue") # retta di regressione

lines(newx, conf_interval[,2], col="blue", lty=2) # limite di confidenza inferiore

lines(newx, conf_interval[,3], col="blue", lty=2) # limite di confidenza superiore

#

plot(var_y, residuals(reglin), xlab="Altezza dei figli (pollici)", ylab="Altezza osservata - calcolata (pollici)", main="Analisi delle differenza residue") # grafico delle differenza tra altezza osservata e altezza calcolata con l'equazione della retta

#

influencePlot(reglin, fill=FALSE, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

#

qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

#

Anche in questo caso, per quanto concerne l'analisi statistica, le principali conclusioni sono quelle tratte con il test globale per l'assunto di linearità effettuato mediante la funzione gvlma(), che però questa volta ci dice che gli assunti che stanno alla base del modello di regressione lineare, a parte la curtosi, non sono per niente soddisfatti:

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                    Value   p-value                   Decision

Global Stat        25.489 4.011e-05 Assumptions NOT satisfied!

Skewness            8.979 2.731e-03 Assumptions NOT satisfied!

Kurtosis            1.965 1.610e-01    Assumptions acceptable.

Link Function       4.632 3.138e-02 Assumptions NOT satisfied!

Heteroscedasticity  9.913 1.641e-03 Assumptions NOT satisfied!

Il test di Bonferroni non evidenzia dati anomali, ma segnala il dato numero 1 come quello che maggiormente si discosta dai rimanenti:

> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers) 

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

   rstudent unadjusted p-value Bonferonni p

1 -3.512671         0.00046509      0.43161

Le conclusioni dell'analisi grafica confermano, con i valori delle differenze residue che mostrano una proporzionalità diretta con l'altezza, che i dati in questione violano uno degli assunti del modello di regressione lineare 

e confermano anche la possibile presenza di dati anomali con il grafico delle distanze di Cook e con il grafico dei quantili, dati che sono evidenziati sia direttamente nei grafici, sia nei rispettivi riepiloghi riportati nella Console di R:

> influencePlot(lm(galton$child ~ galton$parent), xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

       StudRes         Hat       CookD

1   -3.5126706 0.002699826 0.016499436

2   -2.9226403 0.001090010 0.004622768

857  0.4839886 0.008511029 0.001006222

897  2.6611087 0.003740530 0.013207269

898  0.9327550 0.008511029 0.003734739

> #

> qqPlot(lm(galton$child ~ galton$parent), xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

[1] 1 2

I dati numero 1, 2, 857, 897, 898 forniti dalla funzione influencePlot() stanno di nuovo a indicare i casi che influenzano in modo importante la regressione, casi dei quali sarebbe importante controllare la validità, o che potrebbero essere situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati. Anche la funzione qqPlot() fornisce, oltre al grafico, l'indicazione di due punti da controllare, che sono punti 1, 2 già indicati dalla funzione precedente.

Il set di dati galton, oltre a non soddisfare i requisiti fondamentali che si richiedono ai dati per l'applicazione della regressione lineare ordinaria, è anche un esempio di dati nei quali non è chiaro, a dispetto delle conclusioni tratte da Francis Galton di una "regressione verso la media" delle altezze dei figli rispetto a quelle dei padri, quale delle due variabili debba essere considerata la variabile indipendente. Le implicazioni di questo fatto, le conseguenze che esso determina nelle conclusioni tratte dalla regressione lineare, e il calcolo della regressione lineare con modelli alternativi, sono discussi a parte nel post La regressione lineare: assunti e modelli.

Come potete notare entrambi gli script sono stato realizzati in modo da rendere immediato il loro riutilizzo: è sufficiente assegnare alla variabile x (var_x <-) e alla variabile y (var_y <-) i vostri nuovi dati (e personalizzare opportunamente titoli e legende).

Per una guida rapida all'importazione dei dati potete consultare i link:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestione dei dati mancanti

----------

[1] Bossi A, Cortinovis I, Duca PG, Marubini E. Introduzione alla statistica medica. La Nuova Italia Scientifica, Roma, 1994, ISBN 88-430-0284-8. Il modello statistico nella regressione lineare, pp-305-308.

[2] Snedecor GW, Cochran WG. Statistical Methods. The Iowa State University Press, 1980, ISBN 0-8138-1560-6. The matematical model in linear regression, pp.153-157.

[3] Regressione lineare.

https://it.wikipedia.org/wiki/Regressione_lineare

[4] Vedere il post La regressione lineare: assunti e modelli. 

[5]  Vedere il post Regressione lineare semplice non parametrica.

[6] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[7] Digitate help(lm) nella Console di R per la documentazione della funzione lm().

[8] Vedere il post Il set di dati galton.

[9] Informazioni esaurienti sono contenute nei manuali di riferimento dei pacchetti aggiuntivi qui impiegati, che trovate alla pagina Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html