Analizzare graficamente la distribuzione di una variabile

Questa volta l'obiettivo è di fornire uno strumento pratico, immediatamente riadattabile a nuovi dati, che condensa in una sola immagine quattro grafici che illustrano la distribuzione di una variabile: una sintesi grafica che rappresenta un utile complemento ai test di normalità [1].

Come dati impieghiamo la concentrazione delle ferritina nel sangue rilevata in 202 atleti australiani riportata nella variabile ferr della tabella ais inclusa nel pacchetto DAAG. Potete scaricare il pacchetto dal CRAN oppure acquisire la tabella seguendo le indicazioni alternative riportate in [2].

Incollate questo script nella Console di R e premete ↵ Invio:

# ANALIZZARE GRAFICAMENTE LA DISTRIBUZIONE DI UNA VARIABILE

#

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

mydata <- ais$ferr # salva in mydata i valori della variabile ferritina

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

#

# istogramma e kernel density plot

#

hist(mydata, xlim=c(0, 250), ylim=c(0, 0.012), freq=FALSE, breaks="FD", main="Istogramma\ne kernel density plot", xlab="Ferritina in µg/L", ylab="Densità di probabilità") # traccia l'istogramma dei dati

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

plot(density(mydata, n=1024, from=0, to=250), xlim=c(0, 250), ylim=c(0, 0.012), yaxt="n", col="black") # sovrappone il kernel density plot della distribuzione campionaria

#

# distribuzione campionaria vs. gaussiana teorica

#

plot(density(mydata, n=1024, from=0, to=250), xlim=c(0, 250), ylim=c(0, 0.012)) # traccia il kernel density plot della distribuzione campionaria

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

x <- seq(0, 250, length.out=1024) # calcola i valori in ascisse della gaussiana teorica

y <- dnorm(x, mean=mean(mydata), sd=sd(mydata)) # calcola i valori in ordinate della gaussiana teorica

plot(x, y, xlim=c(0, 250), ylim=c(0, 0.012), yaxt="n", col="red", type="l") # sovrappone la distribuzione gaussiana teorica in colore rosso

title(main="Distribuzione campionaria\nvs. gaussiana teorica", xlab="Ferritina in µg/L", ylab="Densità di probabilità") # aggiunge titolo e legende degli assi

#

# distribuzione cumulativa campionaria vs. distribuzione cumulativa teorica

#

par(new=FALSE, ann=TRUE) # per mostrare titolo e legende degli assi

plot(ecdf(scale(mydata)), main="Cumulativa campionaria\nvs. cumulativa teorica", xlab="Deviata normale standardizzata z", ylab="Frequenza cumulativa", xlog = FALSE, ylog = FALSE, xlim=c(-4, 4), ylim=c(0, 1), xaxp=c(-4, 4, 5), yaxp=c(0, 1, 5)) # traccia il grafico della distribuzione cumulativa campionaria

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

plot(ecdf(rnorm(10000, mean=0, sd=1)), col="red", xlog=FALSE, ylog=FALSE, xlim=c(-4, 4), ylim=c(0, 1), xaxp=c(-4, 4, 5), yaxp=c(0, 1, 5)) # sovrappone la distribuzione cumulativa teorica in colore rosso

#

# quantili campionari vs. quantili teorici

#

par(new=FALSE, ann=TRUE) # per mostrare titolo e legende degli assi

qqnorm(scale(mydata), main="Quantili campionari\nvs. quantili teorici", xlab="Quantili teorici", ylab="Quantili campionari", xlim=c(-4, 4), ylim=c(-4, 4)) # traccia il grafico della distribuzione dei quantili campionari

abline (0, 1, col="red") # sovrappone la distribuzione dei quantili teorica in colore rosso

#

La ferritina è una proteina che, pur con alcune limitazioni, fornisce una misura dei depositi di ferro presenti nell'organismo, necessari per garantire una adeguata produzione di emoglobina. L'analisi grafica dei dati consente di confermare graficamente, per la concentrazione della ferritina nel sangue, una distribuzione non gaussiana.

Il codice dello script è illustrato dai commenti via via inseriti, ma vediamo brevemente il significato di questa rappresentazione.

Il primo grafico presenta la distribuzione dei dati sotto forma del classico istogramma [3] impiegando la funzione hist(), calcolando il numero di classi (breaks="FD") mediante la regola di Freedman-Diaconis [4], e ponendo in ordinate la densità di probabilità (freq=FALSE).

Questo consente di sovrapporre all'istogramma mediante la funzione plot(), con la stessa scala delle ordinate, il corrispondente kernel density plot, la cui densità di probabilità è calcolata mediante la funzione density() in corrispondenza di 1024 valori (n=1024) compresi tra 0 (from=0) e 250 (to=250).

Con il secondo grafico entriamo nel vivo del confronto tra la distribuzione dei dati osservata, o distribuzione empirica o distribuzione campionaria che dir si voglia, e quella che i dati dovrebbero avere se fossero distribuiti secondo una gaussiana. 

Questo viene fatto riportando prima il kernel density plot calcolato come al punto precedente, e riportando per confronto (in colore rosso) la distribuzione gaussiana teorica corrispondente, ottenuta calcolando con la funzione dnorm() la densità di probabilità (y) - quella che avremmo se dati con la media mean=mean(mydata) e la deviazione standard sd=sd(mydata) dei dati osservati fossero distribuiti in modo gaussiano - in corrispondenza di 1024 valori (length.out=1024) della x ottenuti con la funzione seq().

Il terzo grafico riporta la distribuzione cumulativa campionaria calcolata con la funzione ecdf() sui dati (mydata) del campione, e riporta per confronto (in colore rosso) la corrispondente distribuzione cumulativa teorica, calcolata sempre con la funzione ecdf() su un campione di 10 000 dati distribuiti in modo gaussiano con media 0 e deviazione standard 1, generati impiegando la funzione rnorm() con gli argomenti 10 000, mean=0, sd=1.

Il quarto grafico impiega la funzione qqplot() per calcolare i quantili della distribuzione campionaria in funzione dei quantili teorici, e ne confronta l'andamento con quello previsto (la linea retta in colore rosso) per dati che seguono una distribuzione gaussiana. Concettualmente, si tratta della trasformazione lineare del grafico della distribuzione cumulativa di cui al punto precedente. 

In definitiva vediamo come i grafici, realizzati in modo semplice, forniscono una indicazione concisa di come e di quanto una distribuzione campionaria si discosta da una distribuzione gaussiana. E, come già detto all'inizio, offrono una sintesi che rappresenta un indispensabile complemento ai test statistici di normalità (gaussianità) [1].

Aggiungo che, nonostante siano ovviamente tecnicamente ineccepibili, distribuzione cumulativa e quantili tendono a comprimere, dal punto di vista grafico, la differenza tra la distribuzione campionaria e la distribuzione attesa nel caso di dati distribuiti in modo gaussiano. Tale differenza per certi versi appare più evidente nel semplice confronto tra la distribuzione osservata, rappresentata sotto forma di kernel density plot e la corrispondente distribuzione gaussiana (grafico in alto a destra).

I dettagli delle funzioni e degli argomenti impiegati sono riportati nei post che trattano le rispettive rappresentazioni grafiche (vedere alla pagina Indice) e sono come sempre accessibili digitando help(nomedellafunzione) nella Console di R.

Si fa infine notare la comparsa nel codice R di un \n che impone un “a capo” in alcune delle righe di testo.

Per riutilizzare questo script per analizzare altri dati è sufficiente:

→ nella prima riga di codice sostituire ais con il nome che volete assegnare ai vostri dati e sostituire C:/Rdati/ais.csv con nome e posizione del file dal quale importare i dati, adeguando al bisogno il separatore dei campi (sep=";") e il separatore dei decimali (dec=",") a quelli da voi impiegati;

→ essendo (ad esempio) dati il nome dei vostri dati e var il nome della variabile che volete analizzare, sostituire nella seconda riga di codice ais$ferr con dati$var;

→ adattare opportunamente xlim, ylim, titoli e legende degli assi.

Vale infine la pena di ricordare che, impiegando l'utilità riportata nel post Salvare i grafici di R in un file, potete anche salvare le immagini realizzate con R, come quella qui riportata, sotto forma di file .bmp, .jpeg, .png, .pdf, .ps per poterle stampare, archiviare, inserire in una pubblicazione, in un post o in un sito web.

----------

[1] Vedere ad esempio la sintesi dei test statistici riportata nel post Tabulare una serie di test di normalità (gaussianità).

[2] Vedere nel post Il set di dati ais come acquisire i dati della tabella ais anche senza istallare il pacchetto DAAG. La concentrazione della ferritina è espressa in microgrammi per litro di siero (µL). 

[3] Da notare che a partire da R versione 4.0 le barre dell'istogramma di default sono riempite in colore grigio. Per una trattazione dettagliata, che include anche vari script con rappresentazoni utili nella pratica, vedere il post Istogrammi.

[4] Vedere: Freedman–Diaconis rule.

https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule

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 non parametrica

La caratteristica saliente dei metodi di regressione lineare non parametrica è costituita dalla loro robustezza. Questo significa che si tratta di metodi che risentono poco della eventuale presenza di dati apparentemente anomali, ma che non possono essere esclusi a priori dal calcolo della regressione.

Essendo n il numero delle coppie di dati/punti esistono n · (n - 1) / 2 modi di connettere due punti qualsiasi con una retta, cioè esistono n · (n - 1) / 2 possibili coefficienti angolari [1]. Nel caso più semplice il coefficiente angolare b della retta di regressione è la mediana di questi valori e l'intercetta a è la mediana degli n valori possibili calcolati mediante il coefficiente angolare trovato. Questi sono concettualmente i principi alla base della regressione lineare non parametrica, con possibili varianti da metodo a metodo.

Qui vediamo, confrontati con la regressione lineare semplice parametrica (regressione lineare ordinaria), tre modelli di regressione lineare non parametrica:

→ la regressione lineare di Thiel-Sen;

→ la regressione lineare di Siegel;

→ la regressione quantilica.

Da notare che questi tre modelli di regressione lineare non parametrica assumono la x come variabile indipendente analogamente alla regressione lineare ordinaria.

Le funzioni necessarie per il calcolo delle tre regressioni lineari non parametriche sono incluse in due pacchetti aggiuntivi, mblm e quantreg, che devono essere scaricati dal CRAN. Trovate come al solito la documentazione dei pacchetti e in particolare il loro manuale di riferimento sul repository della documentazione [2].

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

# REGRESSIONE LINEARE ORDINARIA E NON PARAMETRICA DI THEIL-SEN, DI SIEGEL, QUANTILICA

#

library(mblm) # carica il pacchetto per le regressioni di Theil-Sen e Siegel

library(quantreg) # carica il pacchetto per la regressione quantilica

#

var_x <- c(18, 46, 104, 72, 24, 63, 23, 58, 36, 82, 94) # variabile in ascisse

var_y <- c(21, 42, 109, 64, 102, 78, 19, 66, 42, 89, 91) # variabile in ordinate

#

regpar_yx <- lm(var_y ~ var_x) # regressione lineare ordinaria x variabile indipendente

a_yx <- regpar_yx$coefficients[1] # intercetta a

b_yx <- regpar_yx$coefficients[2] # coefficiente angolare b

#

reg_TS <- mblm(var_y ~ var_x, repeated=FALSE) # regressione lineare di Theil-Sen

a_TS <- reg_TS$coefficients[1] # intercetta a

b_TS <- reg_TS$coefficients[2] # coefficiente angolare b

#

reg_S = mblm(var_y ~ var_x, repeated=TRUE) # regressione lineare di Siegel

a_S <- reg_S$coefficients[1] # intercetta a

b_S <- reg_S$coefficients[2] # coefficiente angolare b

#

reg_q <- rq(var_y ~ var_x, tau=0.5) # regressione lineare quantilica

a_q <- reg_q$coefficients[1] # intercetta a

b_q <- reg_q$coefficients[2] # coefficiente angolare b

#

# traccia il grafico dei dati con le rette di regressione e aggiunge una legenda

#

windows() # apre una nuova finestra

plot(var_x, var_y, xlim = c(0,120), ylim = c(0,120), pch=1, xlab="Variabile x", ylab="Variabile y", main="Regressioni ordinaria, Theil-Sen, Siegel, quantilica", cex.main = 0.9) # grafico dei dati

abline(a_yx, b_yx, col="blue", lty=4) # retta regressione ordinaria

abline(a_TS, b_TS, col="red", lty=4) # retta regressione di Theil-Sen

abline(a_S, b_S, col="green", lty=1) # retta regressione di Siegel

abline(a_q, b_q, col="orange", lty=1) # retta regressione quantilica

legend(60, 30, legend=c("Regressione ordinaria", "Regressione di Theil-Sen", "Regressione di Siegel", "Regressione quantilica"), col=c("blue", "red", "green", "orange"), lty=c(4,4,1,1), cex=0.8) # aggiunge al grafico la legenda

#

# crea una tabella con intercetta e coefficiente angolare delle quattro rette di regressione

#

cells <- c(a_yx,b_yx,a_TS,b_TS,a_S,b_S,a_q,b_q) # genera l'array cells con i valori numerici di a e di b

rnames <- c("Regressione ordinaria", "Regressione di Theil-Sen", "Regressione di Siegel", "Regressione quantilica") # genera l'array rnames con i nomi delle righe

cnames <- c("Intercetta a", "Coefficiente angolare b") # genera l'array cnames con i nomi delle colonne

tabris <- matrix(cells, nrow=4, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la tabella con i risultati a partire dagli array cells, rnames, cnames

tabris # mostra la tabella con i risultati

#

Questo è il grafico di dispersione (xy) dei dati, con le rette calcolate:

I risultati delle quattro equazioni

y = a + b · x

sono così riportati nella Console di R:

> tabris # mostra la tabella con i risultati

                         Intercetta a Coefficiente angolare b

Regressione ordinaria       24.078677               0.7389267

Regressione di Theil-Sen     6.529412               0.9852941

Regressione di Siegel        5.873402               1.0042750

Regressione quantilica       2.581395               1.0232558

Le cose da notare sono:

→ è evidente la presenza di un dato che si discosta in modo importante da tutti gli altri;

→ i risultati dei tre metodi non parametrici sono simili tra loro;

→ i risultati delle regressione ordinaria (parametrica) si discostano da quelli delle regressioni non parametriche.

Per verificare quanto il dato, chiaramente anomalo, influenzi le regressioni, eseguiamo nuovamente lo script, sostituendo la terza e la quarta riga con queste due, nelle quali è stato eliminato il dato, il quinto, quello con le coordinate (x,y) uguali a (24, 102):

var_x <- c(18, 46, 104, 72, 63, 23, 58, 36, 82, 94) # variabile in ascisse

var_y <- c(21, 42, 109, 64, 78, 19, 66, 42, 89, 91) # variabile in ordinate

Ora i grafici delle rette sono praticamente sovrapposti (e difficilmente distinguibili tra loro)

mentre i risultati sono diventati i seguenti:

> tabris # mostra la tabella con i risultati

                         Intercetta a Coefficiente angolare b

Regressione ordinaria        1.337096                1.019512

Regressione di Theil-Sen     2.729167                1.020833

Regressione di Siegel        2.654334                1.022497

Regressione quantilica       2.581395                1.023256

Come si vede dal confronto con i risultati precedenti, con l'eliminazione del dato anomalo i risultati della regressione ordinaria sono cambiati: e questa fornisce ora risultati sovrapponibili a quelli dei tre metodi non parametrici. Poco sono cambiati, rispetti ai precedenti, i risultati della regressione di Thiel-Sen e di Siegel, mentre la regressione quantilica ha addirittura fornito identici valori nei due casi, un fatto molto interessante. I metodi non parametrici confermano la loro robustezza, mentre viene documentato quanto un singolo dato anomalo possa influire sulla regressione lineare ordinaria.

----------

[1] Per 100 coppie di dati il numero dei coefficienti angolari (sui quali calcolare la mediana) è uguale quindi a 4 950. Con 500 dati passiamo a 124 740. Con 1 000 dati passiamo a 499 500. 

[2] Available CRAN Packages By Name.

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