Salvare i grafici di R in un file

Una rappresentazione grafica dei dati è importante non solo visualizzarla, ma anche salvarla sotto forma di file per archiviarla o inserirla in una pubblicazione, in un post o in un sito web.

Tecnicamente vi sono due modi per salvare un'immagine sotto forma di file:

→ salvarla in formato raster (detto anche bitmap);

→ salvarla in formato vettoriale.

Un'immagine in formato raster/bitmap è costituita da una matrice di punti (pixel) opportunamente colorati. Quando viene ingrandita la qualità dell'immagine peggiora. In compenso è possibile agire sui pixel modificando luminosità, contrasto, saturazione dei colori dell'immagine. Il formato raster/bitmap è tipico della fotografia.

Un'immagine in formato vettoriale è invece costruita con una tecnica che consente, mediante trasformazioni matematiche (vettori, onde il nome, e matrici), di ingrandirla senza perdita della qualità.

In questo esempio nel quale i caratteri tipografici sono stati ingranditi di circa sei volte è evidente la perdita di qualità del formato raster/bitmap (a sinistra) rispetto al formato vettoriale (a destra):

Nel primo script viene salvato sotto forma di file nei tre formati grafici raster/bitmap più comuni e cioè

→ .bmp (Windows bitmap)

→ .jpeg (Joint Photographic Experts Group)

→ .png (Portable Network Graphics)

questo grafico a scatola con i baffi (boxplot) che illustra la distribuzione della concentrazione nel sangue degli eritrociti (globuli rossi) rcc per sport praticato ricavata dal set di dati ais [1]:

Prima di eseguire lo script è necessario creare (se non l'avete ancora fatto) sul vostro PC o notebook la cartella C:\Rdati\ nella quale verranno salvati i tre file, denominati rispettivamente boxplot.bmp, boxplot.jpg e boxplot.png. 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 per intero lo script riportato qui sotto quindi incollatelo nella Console di R e premete ↵ Invio:

# SALVA I GRAFICI SU DISCO SOTTO FORMA DI FILE come immagini in formato raster (bitmap)

#

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

#

# salva l'immagine in formato .bmp (Windows bitmap)

#

bmp("C:/Rdati/boxplot.bmp", units = "px", width = 1000, height = 1000, pointsize = 24, bg = "white") # predispone nome e formato del file

boxplot(rcc~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # traccia boxplot

dev.off() # salva il file

#

# salva l'immagine in formato .jpeg (Joint Photographic Experts Group)

#

jpeg("C:/Rdati/boxplot.jpg", units = "px", width = 1000, height = 1000, pointsize = 24, bg = "white") # predispone nome e formato del file

boxplot(rcc~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # traccia boxplot

dev.off() # salva il file

#

# salva l'immagine in formato .png (Portable Network Graphics)

#

png("C:/Rdati/boxplot.png", units = "px", width = 1000, height = 1000, pointsize = 24, bg = "white") # predispone nome e formato del file

boxplot(rcc~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # traccia boxplot

dev.off() # salva il file

#

Gli argomenti width=1000 e height=1000 specificano larghezza e altezza dell'immagine in pixel. Questa unità di misura può essere cambiata con l'argomento units che qui è posto uguale a "px" ma che può essere espresso in alternativa come "in" (pollici), "cm" (centimetri) o "mm" (millimetri).

L'argomento pointsize = 24 specifica la dimensione dei caratteri, mentre l'argomento bg = "white" specifica il colore dello sfondo.

Da notare che nella rappresentazione del grafico a scatola con i baffi nella funzione boxplot() compaiono come argomenti:

→ l'argomento boxwex = 0.4 che gestisce la larghezza dei boxplot;

→ l'argomento cex.axis = 0.8 che gestisce la dimensione dei caratteri;

→ l'argomento las = 2 che ruota verticalmente le etichette.

In questo secondo script lo stesso grafico viene salvato nei formati vettoriali

→ .pdf (Portable Document Format)

→ .ps (Postscript)

nella cartella C:\Rdati\ nei file denominati rispettivamente boxplot.pdf e boxplot.ps.

Copiate per intero lo script riportato qui sotto quindi incollatelo nella Console di R e premete ↵ Invio:

# SALVA I GRAFICI SU DISCO SOTTO FORMA DI FILE come immagini in formato vettoriale

#

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

#

# salva l'immagine in formato .pdf (Portable Document Format)

#

pdf("C:/Rdati/boxplot.pdf", width = 7, height = 7) # predispone nome e formato del file

boxplot(rcc~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # traccia boxplot

dev.off() # salva il file

#

# salva l'immagine in formato .ps (postscript)

#

postscript("C:/Rdati/boxplot.ps", width = 7, height = 7) # predispone nome e formato del file

boxplot(rcc~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # traccia boxplot

dev.off() # salva il file

#

Da notare che i formati .pdf e .ps in quanto formati vettoriali prevedono che le dimensioni dell'immagine (width = 7, height = 7) siano espresse in pollici (non avrebbe senso esprimerle in pixel). Inoltre quando estraete dal file .ps il file boxplot.pdf il file preesistente con lo stesso nome generato dalla funzione pdf() verrà sovrascritto.

In entrambi gli script la funzione dev(off) salva il file e chiude il dispositivo (virtuale) che ha generato la stampa, e a questo punto il dispositivo al quale sono inviati i grafici torna ad essere la finestra che viene aperta con la funzione windows() sul monitor/schermo del PC o notebook.

----------

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

Aggiungere la legenda a un grafico

I modi per aggiungere una legenda a un grafico sono dispersi qua e la nei post di questo blog. Vale quindi la pena di farne una sintesi riportando gli esempi più significativi, basati sull'impiego della funzione legend() e realizzati senza la necessità di installare pacchetti aggiuntivi ma impiegando solamente le funzioni statistiche e grafiche di base di R.

Nel primo esempio realizziamo la legenda che combina simboli e testo per identificare le tre rette di regressione riportate nel grafico, calcolate con altrettanti modelli basati su differenti assunti [1].

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

# LEGENDA con simboli e testo

#

x <- c(9,10,11,9,10,10,11) # variabile in ascisse

y <- c(9,10,11,10,9,11,10) # variabile in ordinate

windows() # apre e inizializza una nuova finestra grafica

#

plot(x, y, xlim=c(0,15), ylim=c(0,15), main="LEGENDA con simboli e testo") # grafico dei dati

abline(5, 0.5, lty=2, col="blue") # retta x variabile indipendente

abline(-10, 2, lty=3, col="red") # retta y variabile indipendente

abline(0, 1, lty=4, col="black") # retta media geometrica

#

legend(0, 15, legend=c("regressione x variabile indipendente", "regressione y variabile indipendente", "regressione media geometrica"), lty=c(2,3,4), col=c("blue", "red", "black"), cex=0.8, title="Regressione lineare con tre metodi") # aggiunge al grafico la legenda

#

La legenda necessaria per identificare le tre rette riportate nel grafico viene realizzata mediante la funzione legend() impiegando come argomenti:

→ 0 che specifica il valore della x al quale posizionare l'angolo superiore sinistro della legenda;

→ 15 che specifica il valore della y al quale posizionare l'angolo superiore sinistro della legenda;

→ legend= che specifica con la funzione c( , , ) le tre righe di testo, uno per ciascuna retta, che compaiono nella legenda;

→ lty= che specifica con la funzione c( , , ) i tre stili delle linee corrispondenti a quelli delle rette tracciate nel grafico - da notare che specificando lty (o lwd cioè lo spessore della linea) la funzione legend() traccia automaticamente accanto alle righe di testo della legenda il simbolo rappresentato da un segmento di retta;

→ col= che specifica con la funzione c( , , ) i tre colori corrispondenti e quelli delle rette tracciate nel grafico [2];

→ cex=0.8 che gestisce la dimensione della legenda (la rimpicciolisce lievemente rispetto al valore 1 di default);

→ title= che mostra sulla prima riga della legenda un titolo (facoltativo).

La posizione della legenda può essere definita anche sostituendo i due primi argomenti con "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" e "center", con una personalizzazione meno spinta di quella possibile con l'impiego delle coordinate cartesiane, ma che può essere adeguata, come vediamo nel prossimo script.

Il secondo esempio illustra come generare una legenda che combina ancora simboli e testo, ma in modo leggermente diverso.

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

# LEGENDA con simboli e testo

#

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

f <-as.matrix(subset(ais, sex=="f", select=c(rcc))) # estrae il sottoinsieme dei dati con sesso=f

m <-as.matrix(subset(ais, sex=="m", select=c(rcc))) # estrae il sottoinsieme dei dati con sesso=m

windows() # apre e inizializza una nuova finestra grafica

#

hist(f, xlim=c(3,7), ylim=c(0,30), breaks=sqrt(length(f)), border="red", col="red", density=20, angle=45, main="LEGENDA con simboli e testo", xlab="Eritrociti in 10^12/L", ylab="Frequenza (numero dei casi)") # istogramma di f

hist(m, add=TRUE, xlim=c(3,7), ylim=c(0,30), breaks=sqrt(length(m)), border="green4", col="green4", density=20, angle=-45) # istogramma di m

#

legend("topright", legend=c("sesso f","sesso m"), fill=c("red","green4"), border=c("red","green4"), density=c(20,20), angle=c(45,-45)) # riporta la legenda

#

La legenda che identifica i due istogrammi è realizzata con la funzione legend() impiegando come argomenti:

→ "topright" che posiziona la legenda in alto a destra;

→ legend= che specifica le due righe di testo che compaiono nella legenda;

→ fill= che determina la comparsa accanto alle righe di testo di due quadrati che avranno all'interno i colori specificati;

→ border= che definisce i colori del bordo dei quadrati di cui al punto precedente;

→ density= che specifica la densità del colore indicato dall'argomento fill= e che con un valore uguale a 20 (linee per pollice) genera in entrambi i quadrati un tratteggio con questa densità, identica a  quella dei corrispondenti istogrammi;

→ angle=... che indicare l'angolo di pendenza del tratteggio che con un valore positivo dell'angolo (45) riprende lo stile del primo istogramma, inclinato da sinistra a destra dal basso in alto, e con un valore negativo dell'angolo (-45) riprende lo stile del secondo istogramma, inclinato da sinistra a destra dall'alto in basso.

Anche in questo caso la funzione c() consente di specificare uno per uno i valori degli argomenti per riportarli ai corrispondenti elementi del grafico.

Il terzo esempio illustra come impiegare la funzione paste() per riportare all'interno di una legenda righe che combinano testo e numeri.

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

# LEGENDA con testo e valori numerici

#  

x <- c(10376,9159,8290,7924,6872) # valori in ascisse 

y <- c(2441,2135,1948,1812,1756) # valori in ordinate 

windows() # apre e inizializza una nuova finestra grafica

#

reglin <- lm(y ~ x) # calcola la regressione lineare 

a <- reglin$coefficients[1] # intercetta a

b <- reglin$coefficients[2] # coefficiente angolare b

r <- cor(x, y, method="pearson") # coefficiente di correlazione r 

#  

plot(x, y, xlim=c(6000,11000), ylim=c(1500,2500), main="LEGENDA con testo e valori numerici") # grafico dei dati 

lines(c(min(x), max(x)), c(a+b*min(x), a+b*max(x)), col="black", lty=2, lwd=1) # retta di regressione  

#

legend(6000, 2500, legend=c(paste("Intercetta a =", round(a, digits=0)), paste("Coefficiente angolare b =", round(b, digits=3)), paste("Equazione: y =", round(a, digits=0), ifelse(sign(b) == 1, "+", "-"), round(b, digits=3), "· x"), paste("Coefficiente di correlazione r =", round(r, digits=3))), cex=0.8, title="Regressione lineare") # aggiunge al grafico la legenda

#

In questo caso la legenda, che riporta nel grafico l'equazione della retta di regressione e il coefficiente di correlazione lineare, viene realizzata impiegando come argomenti della funzione legend():

→ 6000 che specifica il valore della x al quale posizionare l'angolo superiore sinistro della legenda;

→ 2500 che specifica il valore della y al quale posizionare l'angolo superiore sinistro della legenda;

→ legend= che specifica mediante la funzione c() le quattro righe che compaiono nella legenda, ciascuna delle quali è generata mediante la funzione paste() che combina uno o più testi inclusi nelle virgolette "...." con uno o più valori numerici arrotondati mediante la funzione round(), mentre la funzione ifelse() riporta il segno corretto nell'equazione della retta di regressione;

→ cex=0.8 che gestisce la dimensione della legenda (la rimpicciolisce lievemente rispetto al valore 1 di default);

→ title= che mostra sulla prima riga della legenda un titolo (facoltativo).

Il quarto esempio impiega i dati della tabella ais inclusa nel pacchetto DAAG. Potete scaricare il pacchetto dal CRAN oppure acquisire la tabella seguendo le indicazioni alternative riportate in [3].

L'esempio illustra come generare una legenda liberamente posizionabile: dopo avere tracciato i grafici separati per sport praticato e per sesso, lo script rimane in attesa. A questo punto bisogna posizionare il mouse dove si vuole che compaia la legenda e fare click con il tasto sinistro del mouse per farla comparire e terminare lo script.

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

# LEGENDA liberamente posizionabile

#

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

attach(ais) # impiega direttamente i nomi delle variabili senza specificare la tabella

windows() # apre e inizializza una nuova finestra grafica

#

boxplot(rcc~sex+sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="LEGENDA liberamente posizionabile", xlab="", ylab="Eritrociti in 10^12/L", notch=FALSE, col=c("green", "yellow")) # eritrociti per ciascuno sport praticato

#

legend(locator(1), legend=c("sesso f","sesso m"), fill=c("green", "yellow")) # posiziona la legenda

#

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

#

La legenda viene realizzata mediante la funzione legend() impiegando come argomenti:

→ la funzione locator() che resta in attesa, legge la posizione del cursore grafico quando viene fatto click con il tasto sinistro del mouse, e posiziona la legenda;

→ legend= che specifica le due righe di testo che compaiono nella legenda;

→ fill= che determina la comparsa accanto alle righe di testo di due quadrati dei colori specificati.

Il codice riportato non consente di spostare la legenda. Se non si è soddisfatti della sua posizione, è necessario rieseguire l’intero script e fare nuovamente click con il tasto sinistro del mouse nel punto in cui si vuole posizionare la legenda.

Per le semplici legende che si possono inserire nei grafici a torta e nei grafici a barre si rimanda ai relativi post [4, 5].

Gli esempi forniti consentono di realizzare facilmente il codice necessario per inserire le legende nei propri grafici. Si consiglia di consultare il Manuale di Riferimento di R ovvero di digitare help(nomedellafunzione) nella Console di R per la documentazione degli argomenti che consentono di adattare e personalizzare le funzioni impiegate.

Infine se al vostro grafico avete aggiunto una legenda probabilmente volete non solo visualizzarlo, ma anche salvarlo sotto forma di file per archiviarlo o inserirlo in una pubblicazione, in un post o in un sito web. Trovate nel post Salvare i grafici di R in un file come farlo nei formati più diffusi.

----------

[1] L'argomento è trattato nel post La regressione lineare: assunti e modelli. 

[2] Vedere il post Visualizzare i colori disponibili in R.

[3] Vedere nel post Il set di dati ais come acquisire i dati della tabella ais anche senza istallare il pacchetto DAAG.

[4] Vedere il post Grafici a torta.

[5] Vedere il post Grafici a barre [1] e il post Grafici a barre [2].

La regressione lineare: assunti e modelli

La regressione lineare (x variabile indipendente) viene trattata estesamente in tutti i testi di statica, ai quali si rimanda per gli aspetti tecnici [1]. Tuttavia esiste un aspetto relativo al calcolo della regressione lineare che viene spesso sottovalutato, ed è il fatto che fornisce risultati che dipendono da quale delle due variabili in esame sia assunta come variabile indipendente, che per definizione, per le esigenze imposte dal modello matematico impiegato, deve essere posta sull'asse delle ascisse (x). L'altra variabile, la variabile dipendente, viene posta sull'asse delle ordinate (y).

Si considerino i seguenti dati (simulati):

Variabile 1	Variabile 2
9	9
10	10
11	11
9	10
10	9
10	11
11	10

Non avendo alcuna ragione per assumere l'una o l'altra delle due variabili come la variabile indipendente:

→ prima calcoliamo l'equazione della retta di regressione x variabile indipendente ponendo sull'asse delle ascisse (x) la Variabile 1;

→ poi ricalcoliamo l'equazione della retta di regressione x variabile indipendente ponendo sull'asse delle ascisse (x) la Variabile 2;

→ ricaviamo da quest'ultima l'equazione della retta di regressione y variabile indipendente;

→ infine calcoliamo l'equazione della retta di regressione con un metodo che non impiega nessuna delle due come variabile indipendente, ma calcola la media geometrica dei coefficienti angolari della regressione x variabile indipendente e della regressione y variabile indipendente (regressione media geometrica) [2, 3].

Copiate la prima parte dello script [4], quindi incollatela nella Console di R e premete ↵ Invio:

# MODELLI DI REGRESSIONE PARAMETRICA retta y = a + b·x con dati simulati (1/4)

#

windows() # apre una nuova finestra

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

#

var_1 <- c(9,10,11,9,10,10,11) # variabile 1

var_2 <- c(9,10,11,10,9,11,10) # variabile 2

#

x <- var_1 # variabile 1 in ascisse

y <- var_2 # variabile 2 in ordinate

plot(x, y, xlim = c(0,15), ylim = c(0,15), xlab="Variabile 1", ylab="Variabile 2", main="Regressione x variabile indipendente", cex.main = 0.9) # grafico dei dati

regpar <- lm(y ~ x) # regressione lineare standard

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

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

abline(a_yx, b_yx, col="blue") # retta di regressione x variabile indipendente

c(a_yx, b_yx) # mostra a e b

#

Il risultato di questa prima parte dello script, nel quale l'equazione della retta di regressione x variabile indipendente è calcolata ponendo sull'asse delle ascisse (x) la Variabile 1, è il primo grafico in alto a sinistra:

A questo punto scambiamo tra di loro le due variabili e ricalcoliamo l'equazione della retta di regressione con questa seconda parte dello script. Copiatela e incollatela nella Console di R e premete ↵ Invio:

x <- var_2 # variabile 2 in ascisse (2/4)

y <- var_1 # variabile 1 in ordinate

plot(x, y, xlim = c(0,15), ylim = c(0,15), xlab="Variabile 2", ylab="Variabile 1", main="Variabili invertite", cex.main = 0.9) # grafico dei dati

regpar_1 <- lm(y ~ x) # regressione lineare standard

a_yx_1 <- regpar_1$coefficients[1] # intercetta a

b_yx_1 <- regpar_1$coefficients[2] # coefficiente angolare b

abline(a_yx_1, b_yx_1, col="blue") # retta di regressione con variabili scambiate

c(a_yx_1, b_yx_1) # mostra a e b

#

Ora la Variabile 2 è riportata in ascisse (e diventa quindi la variabile indipendente) e la Variabile 1 è riportata in ordinate (grafico in alto a destra). L'equazione della retta di regressione è identica alla precedente (questo è atteso visto che Variabile 1 e Variabile 2 comprendono identici valori). Per confrontare questo grafico e questa retta di regressione con quelli del primo grafico dobbiamo capovolgere il grafico orizzontalmente, quindi ruotarlo a destra di 90°, riportando così la Variabile 2 in ordinate e riportando la Variabile 1 in ascisse, cosa che facciamo con questa terza parte dello script. Copiatela e incollatela nella Console di R e premete ↵ Invio:

x <- var_1 # variabile 1 in ascisse (3/4)

y <- var_2 # variabile 2 in ordinate

a_xy <- - a_yx_1 / b_yx_1 # intercetta a

b_xy <- 1 / b_yx_1 # coefficiente angolare b

plot(x, y, xlim = c(0,15), ylim = c(0,15), xlab="Variabile 1", ylab="Variabile 2", main="Regressione y variabile indipendente", cex.main = 0.9) # grafico dei dati

abline(a_xy, b_xy, col="red") # retta di regressione y variabile indipendente

c(a_xy, b_xy) # mostra a e b

#

Abbiamo così ottenuta la regressione y variabile indipendente, illustrata nel grafico in basso a sinistra.

Infine calcoliamo l'equazione della retta di regressione espressa come regressione media geometrica, copiate la quarta e ultima parte dello script e incollatela nella Console di R e premete ↵ Invio:

x <- var_1 # variabile 1 in ascisse (4/4)

y <- var_2 # variabile 2 in ordinate

b_mg <- sqrt(b_yx*b_xy) # coefficiente angolare b

a_mg <- mean(y) - (b_mg * (mean(x))) # intercetta a

plot(x, y, xlim = c(0,15), ylim = c(0,15), xlab="Variabile 1", ylab="Variabile 2", main="Regressione media geometrica", cex.main = 0.9) # grafico dei dati

abline(a_mg, b_mg, col="green") # retta di regressione media geometrica

c(a_mg, b_mg) # mostra a e b

#

Questa regressione (grafico in basso a destra) assume come coefficiente angolare la media geometrica dei coefficienti angolari della regressione x variabile indipendente e della regressione y variabile indipendente. E con la salomonica conclusione che y = x fornisce una soluzione intermedia, che appare più convincente rispetto all'una e all'altra regressione lineare, per il fatto che la retta calcolata risulta essere allineata con l'asse maggiore dell'elissoide che include i dati.

Nel caso del set di dati galton [5] incluso nel pacchetto psychTools ci troviamo proprio in una situazione di questo genere. Nell'analizzare mediante la regressione lineare le altezze dei genitori e le altezze dei figli non vi sono ragioni per assumere le prime come variabile indipendente. Ma non vi sono neppure ragioni per assumere le altezze dei figli come variabile indipendente. Possiamo allora effettuare un'analisi di questi dati provando a calcolare anche per essi la regressione lineare nei tre modi possibili qui riportati:

→ regressione x variabile indipendente;

→ regressione y variabile indipendente;

→ regressione media geometrica.

Se non l'avete già fatto, prima di eseguire lo script dovete scaricare e installare il pacchetto aggiuntivo psychTools. Copiate per intero lo script riportato qui sotto quindi incollatelo nella Console di R e premete ↵ Invio:

# MODELLI DI REGRESSIONE PARAMETRICA retta y = a + b·x set di dati galton

#

library(psychTools) # carica il pacchetto psych incluso il set di dati galton

#

windows() # apre una nuova finestra

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

#

var_1 <- galton$parent # variabile 1

var_2 <- galton$child # variabile 2

#

x <- var_1 # variabile 1 in ascisse

y <- var_2 # variabile 2 in ordinate

plot(x, y, xlim = c(50,80), ylim = c(50,80), xlab="Altezza dei padri", ylab="Altezza dei figli", main="Regressione x variabile indipendente", cex = 0.5, cex.main = 0.9) # grafico dei dati

regpar <- lm(y ~ x) # regressione lineare standard

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

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

abline(a_yx, b_yx, col="blue") # retta di regressione x variabile indipendente

c(a_yx, b_yx) # visualizza a e b

#

x <- var_2 # variabile 2 in ascisse

y <- var_1 # variabile 1 in ordinate

plot(x, y, xlim = c(50,80), ylim = c(50,80), xlab="Altezza dei figli", ylab="Altezza dei padri", main="Variabili invertite", cex = 0.5, cex.main = 0.9) # grafico dei dati

regpar_1 <- lm(y ~ x) # regressione lineare standard

a_yx_1 <- regpar_1$coefficients[1] # intercetta a

b_yx_1 <- regpar_1$coefficients[2] # coefficiente angolare b

abline(a_yx_1, b_yx_1, col="blue") # retta di regressione con variabili scambiate

c(a_yx_1, b_yx_1) # visualizza a e b

#

x <- var_1 # variabile 1 in ascisse

y <- var_2 # variabile 2 in ordinate

a_xy <- - a_yx_1 / b_yx_1 # intercetta a

b_xy <- 1 / b_yx_1 # coefficiente angolare b

plot(x, y, xlim = c(50,80), ylim = c(50,80), xlab="Altezza dei padri", ylab="Altezza dei figli", main="Regressione y variabile indipendente", cex = 0.5, cex.main = 0.9) # grafico dei dati

abline(a_xy, b_xy, col="red") # retta di regressione y variabile indipendente

c(a_xy, b_xy) # visualizza a e b

#

x <- var_1 # variabile 1 in ascisse

y <- var_2 # variabile 2 in ordinate

b_mg <- sqrt(b_yx * b_xy) # coefficiente angolare b

a_mg <- mean(y) - (b_mg * (mean(x))) # intercetta a

plot(x, y, xlim = c(50,80), ylim = c(50,80), xlab="Altezza dei padri", ylab="Altezza dei figli", main="Regressione media geometrica", cex = 0.5, cex.main = 0.9) # grafico dei dati

abline(a_mg, b_mg, col="green") # retta di regressione media geometrica

c(a_mg, b_mg) # visualizza a e b

#

Questo è il grafico risultante.

Come si vede chiaramente assumendo come variabile indipendente i padri o i figli i risultati cambiano. Il problema, generato a scopo didattico nei dati simulati analizzati con il primo script e riportati nella prima figura e che si presenta spontaneamente nei dati di Francis Galton riportati in quest'ultima figura [6], nasce dal fatto che l'equazione della retta di regressione è determinata dall'interazione di due fattori:

→ gli assunti alla base del modello di regressione lineare impiegato;

→ l'informazione che il modello di regressione lineare può ricavare dai dati forniti.

Più i dati sono allineati su una retta:

→ più consistente è l'informazione in essi contenuta;

→ più r si avvicina a 1;

→ meno rilevante è l'influenza degli assunti alla base del modello di regressione sulle conclusioni.
Tanto che nel caso limite di dati perfettamente allineati su una retta, quando r = 1 e l'informazione contenuta nei dati è univoca, tutti i modelli di regressione forniscono lo stesso identico risultato, indipendentemente dagli assunti alla base del modello di regressione.

Meno i dati sono allineati su una retta:

→ meno consistente è l'informazione in essi contenuta;

→ meno r è vicino a 1 (più r si avvicina a 0);

→ più rilevante è l'influenza degli assunti alla base del modello di regressione sulle conclusioni.

Nel caso limite di dati distribuiti in una palla rotonda, quando r = 0 e l'informazione contenuta nei dati è nulla (nessuno sa come far passare una retta da punti così distribuiti), i modelli di regressione manifestano il massimo di divergenza nei risultati, che a questo punto riflettono esclusivamente gli assunti alla base del modello di regressione.

I dati di Galton presentano un coefficiente di correlazione r uguale a 0.4587624 e un coefficiente di determinazione R² = 0.21 che indica che solamente il 21% della variabilità osservata viene spiegato dalla retta. Ci troviamo di fronte a un caso tipico di una situazione nella quale le conclusioni dipendono più dagli assunti alla base del modello impiegato per interpretare i dati che dalla oggettività fornita dai dati stessi. Non potendo assumere né le altezze dei padri né le altezze dei figli come variabile indipendente, è opportuno accettare le conclusioni tratte dal calcolo della (retta di regressione espressa come) regressione media geometrica, la cui caratteristica più convincente è di essere allineata con l'asse maggiore dell'ellissoide che include i dati (Fig. 2 grafico in basso a destra).

----------

[1] Vedere i riferimenti nel post: Regressione lineare semplice parametrica.

[2] Ling Leng et al. Ordinary least square regression, orthogonal regression, geometric mean regression and their applications in aerosol science. J. Phys. 2007;Conf. Ser. 78 012084. DOI: 10.1088/1742-6596/78/1/012084

[3] Shaoji Xu. A Property of Geometric Mean Regression. The American Statistician 2014;68(4);277-281. DOI: 10.1080/00031305.2014.962763 

[4] Nei link che seguono trovate i dettagli:

→ sui simboli dei punti di R che è possibile impiegare nella funzione plot(); 

→ gli stili delle linee di R che è possibile impiegare nella funzione abline(); 

→ i colori di R che potete impiegare;

→ la possibilità di aggiungere la legenda a un grafico.

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

[6] Il codice per riportare le tre rette su di un unico grafico è riportato nel post: Aggiungere la legenda a un grafico.