Stacco estivo

 "Sai ched’è la statistica? È ’na cosa / che serve pe’ fa’ un conto in generale / de la gente che nasce, che sta male, / che more, che va in carcere e che spósa. / Ma pe’ me la statistica curiosa / è dove c’entra la percentuale, / pe’ via che, lì, la media è sempre eguale / puro co’ la persona bisognosa. / Me spiego: da li conti che se fanno / secondo le statistiche d’adesso / risurta che te tocca un pollo all’anno: / e, se nun entra ne le spese tue, / t’entra ne la statistica lo stesso / perché c’è un antro che ne magna due."
(Trilussa)

"Nei tempi antichi non c’erano le statistiche, perciò era necessario ripiegare sulle menzogne."
(Stephen Leacock)

Grafici di dispersione (scatterplot) con ggplot

I grafici di dispersione (scatter plot o scatterplot) forniscono un ottimo esempio delle potenzialità delle funzioni del pacchetto ggplot2 [1] nella rappresentazione grafica dei dati. Funzioni che possono risultare ostiche a chi le affronta per la prima volta, ma che consentono risultati che non si riescono a ottenere impiegando esclusivamente le funzioni base di R [2].

Oltre al pacchetto ggplot2 bisogna avere installato – o bisogna scaricare anche – il pacchetto DAAG [3] che contiene i dati impiegati nell'esempio e il pacchetto gridExtra che nel nostro caso è necessario per combinare più grafici in una sola figura. 

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

# SCATTERPLOT con ggplot 

# 

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

library(ggplot2) # carica il pacchetto per la grafica

library(gridExtra) # carica il pacchetto per combinare i grafici in una sola figura

#

# scatterplot 

plot1 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point()

#

# con regressione loess

plot2 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point() + geom_smooth()

#

# con regressione lineare

plot3 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point() + geom_smooth(method="lm")

#

# con la densità di distribuzione

plot4 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color="red")) + geom_density_2d() + theme_classic() + labs(title="Densità di distribuzione \n tutti i casi", x="Eritrociti in 10^6/μL", y="Emoglobina in g/dL") + theme(plot.title=element_text(hjust=0.5), legend.position="none")

#

# combina i grafici in una sola figura

grid.arrange(plot1, plot2, plot3, plot4, nrow=2, ncol=2) 

#

Lo script è semplice: dopo avere caricato i pacchetti necessari, sono generati quattro grafici che sono salvati separatamente per essere poi, al termine, combinati in una sola figura.

Per realizzare il primo grafico (plot1) viene impiegata la funzione ggplot() che inizializza l'oggetto al quale sono collegate con il segno più [+] le successive funzioni, che sviluppano la grafica e che prevede:

→ come primo argomento il nome della tabella che contiene i dati (ais) [3];

→ come secondo argomento la funzione aes() che specifica la variabile da riportare in ascisse (rcc) e la variabile da riportare in ordinate (hg), che sono rispettivamente la concentrazione degli eritrociti nel sangue espressa in milioni per microlitro di sangue (10^6/µL) e la concentrazione di emoglobina espressa in grammi per decilitro di sangue (g/dL).

Alla funzione ggplot() viene quindi collegata con il segno più [+] la funzione geom_point() che realizza il grafico xy lasciando per tutti gli argomenti i valori previsti di default.

Da notare subito che questa è la struttura base che rimane identica in tutti e quattro i grafici realizzati, che si differenziano tra loro per le componenti grafiche via via aggiunte.

Il secondo grafico (plot2) quindi riprende tal quale la precedente struttura base e aggiunge semplicemente tramite la funzione geom_smooth(), la regressione loess prevista di default, con i suoi limiti di confidenza (la regressione loess è una regressione che collega tra loro una serie di regressioni polinomiali, quelle che localmente meglio si adattano ai dati).

Nel terzo grafico (plot3) impiegando l'argomento method la regressione loess è sostituita con la usuale retta di regressione ("lm"), riportando anche in questo caso i limiti i confidenza. Se non li volete rappresentare correggete la riga di codice sostituendola completamente con la seguente o semplicemente aggiungendo quanto qui riportato in colore 

# con regressione lineare

plot3 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point() + geom_smooth(method="lm", se=FALSE)

#

Il quarto grafico (plot4) riprende la struttura base aggiungendo:

→ con aes(color="red") un colore ai punti che rappresentano i dati;

→ con la funzione geom_density_2d() la rappresentazione sul piano dei livelli di densità nella distribuzione dei dati;

→ con la funzione theme_classic() il cambiamento di stile che risulta evidente nell'ultimo grafico;

→ con la funzione labs() il titolo e le etichette degli assi, Da notare il simbolo \n che all'interno di una stringa di testo determina il ritorno a capo;

→ con la funzione theme() il codice (abbastanza arzigogolato) che determina la centratura del titolo e il valore "none" per l'argomento legend.position che determina la posizione della legenda, che pertanto non viene riportata.

Infine con la funzione grid.arrange() i quattro grafici sono combinati in un'unica figura.

Ora copiate e incollate nella Console di R questo secondo script che ha come obiettivo evidenziare la semplicità con la quale possiamo effettuare in un grafico rappresentazioni separate in base a un fattore che nel nostro caso è la variabile sesso (sex), in quanto il set di dati ais contiene dati ematologici e dati biometrici, rilevati in atleti australiani sia di sesso femminile sia di sesso maschile.

# SCATTERPLOT con ggplot (differenziaziati in base a un fattore)

# 

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

library(ggplot2) # carica il pacchetto per la grafica

library(gridExtra) # carica il pacchetto per combinare i grafici in una sola figura

#

# scatterplot per sesso

plot1 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex))

#

# con simboli dei punti differenziati

plot2 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none")

#

# con rette di regressione e limiti di confidenza

plot3 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + geom_smooth(aes(group=sex, color=sex), method="lm", se=TRUE)

#

# con rette di regressione estese

plot4 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + geom_smooth(aes(group=sex, color=sex), method="lm", se=FALSE, fullrange=TRUE)

#

# combina i grafici in una sola figura

grid.arrange(plot1, plot2, plot3, plot4, nrow=2, ncol=2) 

#

Per realizzare il primo grafico (plot1) viene impiegata la struttura base già vista nello script precedente, ma con la semplice aggiunta a geom.point della funzione aes() che con l'argomento color posto uguale a sex riporta i dati con un colore differenziato in base al sesso.

Il secondo grafico (plot2) riprende tal quale il codice precedente ma:

→ alla funzione geom_point() aggiunge l'argomento shape=sex che automaticamente riporta un simbolo differente in base al sesso;

→ con la funzione theme() esclude dal grafico la legenda.

Il terzo grafico (plot3) è identico al secondo ma con la funzione geom_smooth() aggiunge la retta di regressione ordinaria (method="lm") con i limiti di confidenza (se=TRUE) separatamente (group=sex) e con un colore diverso (color=sex) per ciascun sesso.

Il quarto grafico (plot4) riprende esattamente il codice del terzo ma con se=FALSE esclude la rappresentazione dei limiti di confidenza e con fullrange=TRUE riporta la retta di regressione a tutto campo e al di la dei limiti dei dati osservati (per inciso si ricorda che, a meno di casi particolari e adeguatamente giustificati, l'estrapolazione non andrebbe mai fatta, ma viene qui riportata per completezza).

Infine con la funzione grid.arrange() i quattro grafici sono di nuovo combinati, per semplicità, in un'unica figura.

Ora copiate e incollate nella Console di R questo script e premete ↵ Invio.

# SCATTERPLOT con ggplot (densità delle distribuzioni)

# 

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

library(ggplot2) # carica il pacchetto per la grafica

library(gridExtra) # carica il pacchetto per combinare i grafici in una sola figura

#

# scatterplot con ellissi

plot1 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + stat_ellipse(aes(color=sex))

#

# con densità di distribuzione della y

plot2 <- ggplot(ais, aes(y=hg, fill=sex)) + geom_density(alpha=0.5) + theme(legend.position="none")

#

# con densità di distribuzione della x

plot3 <- ggplot(ais, aes(x=rcc, fill=sex)) + geom_density(alpha=0.5) + theme(legend.position="none")

#

# con livelli di densità separati per sesso

plot4 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + geom_density_2d(aes(color=sex))

#

# combina i grafici in una sola figura

grid.arrange(plot1, plot2, plot3, plot4, nrow=2, ncol=2) 

#

Il primo grafico (plot1) riprende le funzioni già note, ma aggiunge, con la funzione stat_ellipse(), le ellissi che riportano i limiti di confidenza al 95% delle distribuzioni dei dati, separate per ciascun sesso e ulteriormente analizzate con il secondo e il terzo grafico.

Il secondo grafico (plot2) riporta in ordinate per la variabile emoglobina (y=hg) e separatamente per ciascun sesso (fill=sex), con la funzione geom_density(), il kernel density plot, con un colore trasparente al 50% (alpha=0.5) per consentirne la sovrapposizione. Da notare due aspetti non banali gestiti automaticamente dalle funzioni qui impiegate: gli assi sono ruotati nella posizione corretta e per l'asse delle ordinate la scala nella quale sono espressi i valori di concentrazione dell'emoglobina (hg) è allineata con quella del primo grafico.

Il terzo grafico (plot3) riporta in modo analogo ma in ascisse il kernel density plot per sesso della variabile concentrazione degli eritrociti (x=rcc). Da notare che anche in questo caso gli assi sono automaticamente ruotati nella posizione corretta e per l'asse delle ascisse la scala nella quale sono espressi i valori di concentrazione degli eritrociti (rcc) è automaticamente allineata con quella del primo grafico.

In definitiva i primi tre grafici formano un tutt'uno organico – con le variabili rappresentate contemporaneamente sia singolarmente in forma di grafico xy sia in forma di grafici che ne riportano in continuo la densità (kernel density plot) – che fornisce un'analisi molto interessante dei dati a riprova dell'utilità delle funzioni incluse nel pacchetto.

Il quarto grafico (plot4) riprende esattamente il primo aggiungendo questa volta al posto delle ellissi con la funzione geom_density_2d() la rappresentazione sul piano dei livelli di densità nella distribuzione dei dati.

Se vedendo il grafico notate che preferireste (io lo preferisco di gran lunga) che i due kernel density plot fossero riferiti al grafico in basso a destra, niente di più facile, copiate e incollate nella Console di R questo script e premete ↵ Invio.

# SCATTERPLOT con ggplot (densità delle distribuzioni)

# 

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

library(ggplot2) # carica il pacchetto per la grafica

library(gridExtra) # carica il pacchetto per combinare i grafici in una sola figura

#

# scatterplot con ellissi

plot1 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + stat_ellipse(aes(color=sex))

#

# con densità di distribuzione della y

plot2 <- ggplot(ais, aes(y=hg, fill=sex)) + geom_density(alpha=0.5) + theme(legend.position="none") + scale_x_reverse()

#

# con densità di distribuzione della x

plot3 <- ggplot(ais, aes(x=rcc, fill=sex)) + geom_density(alpha=0.5) + theme(legend.position="none")

#

# con livelli di densità separati per sesso

plot4 <- ggplot(ais, aes(x=rcc, y=hg)) + geom_point(aes(color=sex, shape=sex)) + theme(legend.position="none") + geom_density_2d(aes(color=sex))

#

# combina i grafici in una sola figura

grid.arrange(plot1, plot3, plot2, plot4, nrow=2, ncol=2) 

#

Come si vede il codice è identico al precedente, a parte il fatto che dovete:

→ scambiare tra di loro i grafici (plot3, plot2,) quando li combinate in un'unica figura (ultima riga di codice); 

→ aggiungere con (+) al plot2 (che ora diventa il terzo grafico della figura) la funzione scale_x_reverse() per invertire la scala dell'asse delle ascisse.

Ricordate, perché prima o poi vi tornerà utile, che potete anche invertire la scala dell'asse delle ordinale con scale_y_reverse() e scambiare di posizione gli assi con coord_flip().

In questo e negli altri post nei quali il pacchetto è stato (e sarà) impiegato [4] ho cercato (e cercherò) di riportare alcuni esempi dell'interessante ausilio che una analisi grafica con ggplot2 può fornire nella introspezione dei dati e nella comunicazione dell'informazione in essi contenuta. 

Ovviamente potete procedere con ulteriori approfondimenti impiegando la documentazione del pacchetto [1], ma potete anche valutare l'opportunità di ricorrere a qualcuno dei corsi sull'argomento disponibili sul web [5].

----------

[1] Vedere la documentazione e il manuale di riferimento su:

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

[2] Vedere ad esempio il post Grafici di dispersione (grafici xy).

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

[4] Fate click su ggplot2 nelle Parole chiave o, per una ricerca meno selettiva, digitate ggplot nella casella Cerca nel blog quindi fate click su Cerca.

[5] Vedasi ad esempio "Introduction to Data Visualization with ggplot2":

https://www.datacamp.com/courses/introduction-to-data-visualization-with-ggplot2

Multipli e sottomultipli delle unità di misura SI

Statistica significa analisi dei dati, cioè dei risultati di processi di misura espressi in termini di grandezze, di unità di misura e dei loro multipli e sottomultipli. 

Ho riportato altrove [1] la storia che dalle antiche misure antropomorfe, attraverso la svolta che è seguita alla nascita della scienza moderna, ha condotto anche l'Italia a riconoscere la necessità di adottare un linguaggio delle misure rigoroso e ampiamente concordato e condiviso a livello scientifico, oltre che in grado di risolvere le esigenze pratiche quotidiane, prima con l'adesione nel 1875 al Sistema metrico decimale e successivamente adottando per legge nel 1982  il Sistema internazionale di unità (SI) che ne rappresenta la moderna evoluzione [2].

Così come si trova citato anche nella Appendice 1 della Brochure del SI [3], nel 2022 la 27a CGPM [4] con una delle sue risoluzioni [5] ha aggiunto, ai precedenti, due altri multipli (ronna e quetta) e due altri sottomultipli (ronto e quecto), quindi per indicare multipli e sottomultipli delle unità di misura nel SI sono ad oggi previsti i prefissi, i simboli e i corrispondenti fattori di moltiplicazione riportati in questa tabella:

Lo scopo è evidente: se il numero di cifre diventa eccessivo e si trova ostica per l'uso corrente e non si vuole impiegare la notazione scientifica [6], ma si preferisce esprimere in forma verbale i numeri superiori al milione, i prefissi SI consentono di evitare la trappola delle denominazioni che in tutta Europa con la scala lunga [7], e negli USA e in molti altri Paesi con la scala corta, prevedono di impiegare gli stessi termini per indicare numeri completamente diversi, come qui evidenziato,

un fatto che determina confusione e possibili errori nella interpretazione dei numeri espressi in forma verbale quando questi sono riportati da una lingua all'altra [8].

Pertanto – giusto per fare un esempio – per la constatazione che "... as of April 2024 Apple has a market cap of $2.547 Trillion..." [9] è opportuno precisare che i 2.547 trilioni di dollari di capitalizzazione dell'Apple sono espressi nella scala corta dei Paesi anglosassoni, e quindi sono 2.547 ∙ 10¹² dollari, e non sono i 2.547 ∙ 10¹⁸ dollari corrispondenti ai trilioni della scala lunga.

Invece nessun dubbio vi può essere nel caso della larghezza di banda espressa in bit (b) di una rete a 100 Mb/s (100 megabit al secondo ovvero 100 ∙ 10⁶ bit al secondo) e della capacità espressa in byte (B) di una RAM di 16 GB (16 gigabyte ovvero 16 ∙ 10⁹ byte) o di un hard-disk di 4 TB (4 terabyte ovvero 4 ∙ 10¹² byte), dato il significato univoco dei prefissi SI.  

----------

[1] Grandezze e unità di misura. Breve storia dall'antichità al Sistema Internazionale di Unità (SI). Vedere il post:

https://impararfacendo.blogspot.com/2020/10/grandezze-e-unita-di-misura.html

Il testo può essere liberamente scaricato (e senza pubblicità diretta nè occulta) da questi siti

https://www.academia.edu/41041923/

https://books.google.it/books?id=GciVDwAAQBAJ

come pure dal mio sito personale

https://www.bayes.it/SI.pdf

[2] I riferimenti normativi completi altre che nel testo citato sono riportati  nel  post Scale nominali, scale ordinali e scale numeriche.

[3] SI Brochure: The International System of Units (SI) – 9a edizione, 2019.

https://www.bipm.org/en/publications/si-brochure

[4] Sigla della "Conferenza Generale dei Pesi e delle Misure", in sostanza il parlamento attraverso il quale si esprimono in tema di misure gli Stati che hanno adottato il SI.

[5] Vedere: Resolution 3 of the 27th CGPM (2022). On the extension of the range of SI prefixes. The General Conference on Weights and Measures (CGPM), at its 27th meeting.

https://www.bipm.org/en/cgpm-2022/resolution-3

[6] Vedere il post Come riportare i numeri in formato fisso.

[7] La scala lunga la si trova riportata anche in Italia nella legislazione, come ad esempio nella tabella a pagina 12 dell'Allegato al DECRETO MINISTERIALE 4 settembre 1996 Attuazione della direttiva 94/55/CE del Consiglio concernente il ravvicinamento delle legislazioni degli Stati membri relative al trasporto di merci pericolose su strada. (GU Serie Generale n.282 del 02-12-1996 - Suppl. Ordinario n. 211).

[8] Vedere: Long and short scales. 

https://en.wikipedia.org/wiki/Long_and_short_scales

 

[9] Dati tratti da: Largest Companies by Market Cap

https://companiesmarketcap.com/

Regressione con una funzione monotòna (isotòna)

In assenza di un razionale che consenta la scelta di una funzione specifica (una retta o quant'altro) per approssimare una serie di dati, una possibile alternativa è rappresentata da una funzione con una forma determinata da condizioni poco restrittive.

Se partiamo dall'assunto minimo che essendo x₁ < x₂ < ... < xₙ all'aumentare di x la y corrispondente aumenti o al più resti uguale, siamo di fronte per definizione a una funzione monotòna (in quanto non presenta nè minimi né massimi) che è anche isotòna (in quanto sempre crescente o sempre non decrescente) ovvero:

→ sempre crescente quando y₁ < y₂ < ... < yₙ 

→ sempre non decrescente quando y₁ ≤ y₂ ≤ ... ≤ yₙ

Vediamo ora come approssimare i dati con una regressione isotòna [1] impiegando le funzioni del pacchetto isotone [2] e del pacchetto monreg che dovete scaricare dal vostro CRAN preferito e dei quali è utile scaricare anche i rispettivi Reference manual [3].

Dovete inoltre accertarvi di avere installato – o dovete scaricare anche – il pacchetto DAAG [4] che contiene il set di dati ais impiegati nell'esempio e il pacchetto car che viene utilizzato per l'analisi preliminare dei dati.

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

# ANALISI PRELIMINARE DEI DATI 

#

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

library(car) # carica il pacchetto per lo scatterplot

windows() # apre e inizializza una nuova finestra grafica

#

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="Distribuzione dei dati e rette di regressione", data=ais) # mostra il grafico

#

Dopo avere caricato il pacchetto (DAAG) che contiene i dati da analizzare e il pacchetto (car), con windows() viene aperta e inizializzata una nuova finestra nella quale comparirà il grafico predisposto con la funzione scatterplotMatrix() nella quarta e ultima riga di codice [5].  

L'analisi preliminare dei dati ci dice che se in alcuni casi la regressione lineare sembra appropriata, in altri lo è molto meno. Qui ci focalizziamo su emoglobina (hg) e ferritina (ferr) [6] che mostrano distribuzioni dei dati non gaussiane, quindi in linea di principio inadatte all'impiego della regressione lineare ordinaria [7]. 

Ora copiate e incollate nella Console di R questo script e premete ↵ Invio.

# REGRESSIONE MONOTÒNA (ISOTÒNA) con il pacchetto isotone 

#

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

library(isotone) # carica il pacchetto per la regressione 

windows() # apre e inizializza una nuova finestra grafica

#

iso <- gpava(ais$ferr, ais$hg, solver=weighted.mean, ties="secondary") # calcola la regressione isotòna

#

plot(iso, main="Regressione monotòna (isotòna)", xlab="Ferritina in µg/L", ylab="Emoglobina in g/dL", col="red") # riporta il grafico con la distribuzione dei dati e sovrappone la regressione

#

Impieghiamo quindi il pacchetto isotone per calcolare la regressione isotòna con la funzione gpava che prevede come argomenti:

→ il valore in ascisse, nel nostro caso ais$ferr;

→ il valore in ordinate, nel nostro caso ais$hg;

→ l'argomento solver che consente di approssimare i dati impiegando la media (weighted.mean), la mediana (weighted.median), i frattili (weighted.fractile) pesati (qui i singoli pesi non sono specificati quindi per default sono uguali a 1);

→ l'argomento ties che consente di gestire l'ordinamento dei dati nel caso di valori identici, con "primary" l'ordinamento viene effettuato tra i blocchi di dati identici, non all'interno di questi, con "secondary" l'ordinamento viene effettuato tra i blocchi di dati identici, ma viene richiesta l'uguaglianza dei valori all'interno di questi, con "tertiary" viene effettuato solamente l'ordinamento delle medie dei blocchi di dati identici.

Nell'oggetto iso generato dalla funzione gpava sono contenute tutte le informazioni necessarie per riportare in un grafico con la funzione plot() dati e regressione, quindi è sufficiente aggiungere semplicemente titolo, legende e colore ("red") della regressione.

Nel caso di una regressione una cosa che interessa è impiegare la funzione calcolata per interpolare nuovi valori, per farlo copiate e incollate nella Console di R questo breve addendum allo script.

# interpolazione lineare sulla regressione isotòna calcolata in precedenza

#

ris <- as.data.frame(cbind(z=iso$z, y=iso$x)) # riporta in una tabella i dati della regressione

risord <- ris[order(ris$z),] # ordina la tabella

#

reg <- approx(risord$z, risord$y, xout=c(15, 20, 40, 60, 80, 100, 110, 120, 150, 170, 190, 220), ties="ordered") # interpola i valori specificati

#

round(cbind(reg$x, reg$y), digits=1) # mostra i risultati dell'interpolazione

#

Dopo avere riportato i dati della regressione in una tabella e avere ordinato i valori in ascisse in ordine crescente (prerequisito indispensabile per effettuare una interpolazione), i valori di emoglobina corrispondenti a 15, 20, 40, 60, 80, 100, 110, 120, 150, 170, 190, 220 µg/L di ferritina sono calcolati mediante interpolazione lineare sui dati della regressione isotòna con la funzione approx() (terza riga di codice) e sono quindi mostrati (quarta riga) arrotondandoli a un decimale.

> round(cbind(reg$x, reg$y), digits=1) # mostra i risultati dell'interpolazione

      [,1] [,2]

 [1,]   15 13.5

 [2,]   20 13.5

 [3,]   40 14.4

 [4,]   60 14.4

 [5,]   80 14.7

 [6,]  100 14.9

 [7,]  110 14.9

 [8,]  120 14.9

 [9,]  150 15.3

[10,]  170 15.3

[11,]  190 15.9

[12,]  220 15.9

I valori da interpolare riportati nell'argomento xout=c() vanno ovviamente adattati al bisogno. In ogni caso se volete verificare la correttezza del calcolo potete farlo riportando sul grafico i risultati dell'interpolazione con questa semplice riga di codice

# riporta sul grafico i punti interpolati

#

points(reg$x, reg$y, pch=16, col="blue") 

#

e questi sono i punti interpolati riportati sulla regressione isotòna calcolata.

Per meglio sviluppare il tema vediamo ora la regressione calcolata con un metodo diverso, previsto nel pacchetto monreg, copiate e incollate nella Console di R questo script e premete ↵ Invio.

# REGRESSIONE MONOTÒNA (ISOTÒNA) con il pacchetto monreg 

#

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

library(monreg) # carica il pacchetto per la regressione 

windows() # apre e inizializza una nuova finestra grafica

#

mon <- monreg(ais$ferr, ais$hg, hd=.5, hr=.5, monotonie="isoton") # calcola la regressione isotòna

#

plot(ais$ferr, ais$hg, main="Regressione monotòna (isotòna)", xlab="Ferritina in µg/L", ylab="Emoglobina in g/dL") # riporta il grafico con la distribuzione dei dati

lines(mon$t, mon$estimation, col="red") # sovrappone la regressione 

#

In questo caso il passo più interessante si trova nella funzione monreg() ed è l'argomento monotonie che consente di selezionare le due soluzioni possibili per una regressione monotòna:

→ con "isoton" si calcola la regressione isotòna (crescente o sempre non decrescente) come facciamo in questo caso;

→ con "antiton" si calcola la regressione antitòna (decrescente o sempre non crescente), che è trattata a parte in un altro post [8].

Per gli altri argomenti della funzione rimando alla guida specifica richiamabile digitando help(monreg) nella Console di R.

Infine con plot() si riporta il grafico xy dei dati e con lines() si sovrappone la regressione le cui ascisse (mon$t) e ordinate (mon$estimation) sono state in precedenza salvate nell'oggetto mon generato con la funzione monreg(). Digitate str(mon) nella Console di R per visualizzare tutte le variabili salvate da monreg() nell'oggetto mon.

Ora copiate e incollate nella Console di R questo breve addendum allo script che riporta nuovamente – in quanto lievemente diverso e più semplice del precedente poiché la funzione monreg() restituisce i dati con la x già ordinata in ordine crescente – il codice necessario per effettuare l'interpolazione sulla regressione calcolata.

# interpolazione lineare sulla regressione isotòna calcolata in precedenza

#

reg <- approx(mon$t, mon$estimation, xout=c(15 , 20, 40, 60, 80, 100, 110, 120, 150, 170, 190, 220), ties="ordered") # interpola i valori specificati

#

round(cbind(reg$x, reg$y), digits=1) # mostra i risultati dell'interpolazione

#

Anche in questo caso i valori di emoglobina corrispondenti a 15, 20, 40, 60, 80, 100, 110, 120, 150, 170, 190, 220 µg/L di ferritina sono calcolati mediante interpolazione lineare sui dati della regressione isotòna con la funzione approx() (prima riga di codice), quindi sono mostrati (seconda riga).

> round(cbind(reg$x, reg$y), digits=1) # mostra i risultati dell'interpolazione

      [,1] [,2]

 [1,]   15 13.8

 [2,]   20 13.9 

 [3,]   40 14.2

 [4,]   60 14.4

 [5,]   80 14.6

 [6,]  100 14.8

 [7,]  110 14.9

 [8,]  120 15.0

 [9,]  150 15.2

[10,]  170 15.3

[11,]  190 15.5

[12,]  220 15.8

Infine se volete verificare la correttezza del calcolo potete riportare sul grafico i risultati dell'interpolazione, per farlo copiate e incollate nella Console di R questa riga di codice, identica a quella già vista a conclusione dello script precedente.

# riporta sul grafico i punti interpolati

#

points(reg$x, reg$y, pch=16, col="blue") 

#

Tuttavia la cosa più interessante in questo caso è probabilmente il confronto diretto tra la regressione isotòna così calcolata e la regressione lineare standard, cosa che può essere agevolmenta fatta copiando e incollando nella Console di R queste due ultime righe di codice

# 

reglin <- lm(ais$hg ~ ais$ferr) # calcola la regressione lineare

abline(reglin, col="green3") # traccia la retta calcolata

#

con questo risultato finale (retta di regressione in colore verde)

che conferma l'impressione che, in questo specifico caso, e con questo specifico pacchetto, la regressione isotòna corrisponde in larga parte a una regressione lineare, dalla quale tende peraltro a discostarsi significativamente nelle regioni più esterne.

----------

[1] Nella letteratura statistica anglosassone trovate "isotonic" che talora viene erroneamente considerato sinonimo di "monotonic". Qui "monotòna" viene riportato con il significato matematico originario di "funzione che non ha massimi e non ha minimi" mentre "isotòna" e "antitòna" sono gli attributi aggiuntivi di una funzione monotòna quando è rispettivamente una "funzione crescente o sempre non decrescente" e quando è una "funzione decrescente o sempre non crescente".

[2] de Leeuw, J., Hornik, K., & Mair, P. (2009). Isotone Optimization in R: Pool-Adjacent-Violators Algorithm (PAVA) and Active Set Methods. Journal of Statistical Software, 32(5), 1–24. https://doi.org/10.18637/jss.v032.i05

[3] Vedere Reference manual: isotone.pdf al link:

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

e Reference manual: monoreg.pdf al link:

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

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

[5] Per i dettagli della funzione potete digitare help(scatterplotMatrix) nella Console di R ma trovate anche una breve spiegazione nel post Grafici di dispersione (grafici xy).

[6] L'emoglobina è la proteina che trasporta l'ossigeno contenuta nei globuli rossi del sangue, la sua concentrazione è espressa in grammi per decilitro di sangue (g/dL), la ferritina è la proteina collegata ai depositi di ferro presenti nell'organismo che contribuiscono a garantire una adeguata produzione di emoglobina, la sua concentrazione è espressa in microgrammi per litro di siero (µg/L).

[7] Vedere il post Regressione lineare semplice parametrica.

[8] Vedere il post Regressione con una funzione monotòna (antitòna) * in preparazione *