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.
# 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.
# 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.
# 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
che conferma l'impressione che la regressione isotòna corrisponda nella zona centrale alla regressione lineare, dalla quale tende però 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).
[8] Vedere il post Regressione con una funzione monotòna (antitòna) * in preparazione *
