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

Come riportato anche altrove [1], in assenza di un razionale che consenta la scelta di una funzione specifica (retta o 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 diminuisca o al più resti uguale, siamo di fronte per definizione a una funzione monotòna che è anche antitòna (in quanto sempre decrescente o sempre non crescente) ovvero:

→ sempre decrescente quando y₁ > y₂ > ... > yₙ 

→ sempre non crescente quando y₁ ≥ y₂ ≥ ... ≥ yₙ

Vediamo ora come approssimare i dati con una regressione antitòna [2] impiegando sia il pacchetto fdrtool sia il pacchetto monreg che dovete scaricare dal vostro CRAN preferito. Inoltre dovete scaricare anche il pacchetto car che viene utilizzato per l'analisi preliminare dei dati.

Come dati impieghiamo i valori di BMI (indice di massa corporea) rilevati a livello europeo alcuni anni fa e pubblicati dall'Istat [3].

Per proseguire è necessario:

→ effettuare il download del file di dati bmi.csv

→ salvare il file nella cartella C:\Rdati\

Per il file di dati trovate link e modalità di download alla pagina Dati, ma potete anche semplicemente copiare i dati riportati qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarli in C:\Rdati\ in un file di testo denominato bmi.csv (assicuratevi che il file sia effettivamente salvato con l'estensione .csv).

Nazione;sottopeso;normale;sovrappeso;obeso

Austria;2.4;49.6;33.3;14.7

Belgio;2.7;48.0;35.3;14.0

Bulgaria;2.2;43.8;39.2;14.8

Cipro;3.9;47.8;33.8;14.5

Croazia;1.9;40.7;38.7;18.7

Danimarca;2.2;50.0;32.9;14.9

Estonia;2.2;43.9;33.5;20.4

Finlandia;1.2;44.1;36.4;18.3

Francia;3.2;49.6;31.9;15.3

Germania;1.8;46.1;35.2;16.9

Grecia;1.9;41.3;39.4;17.3

Irlanda;1.9;42.3;37.0;18.7

Lettonia;1.7;41.8;35.2;21.3

Lituania;1.9;42.5;38.3;17.3

Lussemburgo;2.8;49.3;32.4;15.6

Malta;2.0;37.0;35.0;26.0

Olanda;1.6;49.0;36.0;13.3

Polonia;2.4;42.9;37.5;17.2

Portogallo;1.8;44.6;36.9;16.6

Regno Unito;2.1;42.2;35.6;20.1

Repubblica Ceca;1.1;42.1;37.6;19.3

Romania;1.3;42.9;46.4;9.4

Slovacchia;2.1;43.6;38.0;16.3

Slovenia;1.6;41.8;37.4;19.2

Spagna;2.2;45.4;35.7;16.7

Svezia;1.8;48.3;35.9;14.0

Ungheria;2.9;41.9;34.0;21.2

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

# ANALISI PRELIMINARE DEI DATI 

#

library(car) # carica il pacchetto per la grafica

bmi <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

windows() # apre e inizializza una nuova finestra grafica

#

scatterplotMatrix(~sottopeso+normale+sovrappeso+obeso, 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=bmi) # mostra il grafico

#

Dopo avere caricato il pacchetto pacchetto (car) per l'analisi preliminare dei dati, questi sono importati nell'oggetto bmi, quindi 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 [4].  

L'analisi preliminare ci dice che, tranne forse per la proporzionalità inversa tra obesi e soggetti con peso normale, la regressione lineare mal si adatta a descrivere i dati. Qui ci focalizziamo sulla relazione che intercorre tra soggetti obesi (obeso) e soggetti in sovrappeso (sovrappeso) per la quale è sicuramente da sconsigliare l'impiego di una regressione lineare. 

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

# REGRESSIONE MONOTÒNA (ANTITÒNA) con il pacchetto fdrtool 

#

library(fdrtool) # carica il pacchetto per la regressione

bmi <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

windows() # apre e inizializza una nuova finestra grafica

#

fdr <- monoreg(bmi$obeso, bmi$sovrappeso, type="antitonic") # calcola la regressione antitòna

#

plot(bmi$obeso, bmi$sovrappeso, main="Regressione monotòna (antitòna)", xlab="Soggetti obesi in %", ylab="Soggetti sovrappeso in %") # riporta il grafico con la distribuzione dei dati

lines(fdr$x, fdr$yf, col="red") # sovrappone la regressione 

#

Dopo avere caricato il pacchetto e i dati, nella funzione monreg() l'argomento type consente di selezionare le due soluzioni possibili per una regressione monotòna:

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

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

I risultati della regressione calcolata con la funzione monoreg() sono salvati nell'oggetto fdr mentre compare un messaggio che avverte che alcuni valori in ascisse sono duplicati (se osservate, nella variabile obeso compare due volte il valore 14.0 e due volte compaiono anche i valori 17.3 e 18.7): ma non si tratta di un problema in quanto per questi valori la funzione provvede automaticamente a unire le y corrispondenti come riportato nel manuale del pacchetto [5].

Quindi con la funzione plot() si riporta il grafico xy dei dati aggiungendo titolo ed etichette degli assi.

Infine con lines() si sovrappone la regressione le cui ascisse (fdr$x) e ordinate (fdr$yf) sono state in precedenza salvate nell'oggetto fdr generato con la funzione monoreg(). Se nella Console di R digitate

str(fdr)

potete visualizzare tutte le variabili salvate da monoreg() nell'oggetto fdr.

Ora copiate e incollate nella Console di R questo breve addendum allo script per effettuare l'interpolazione di alcuni valori sulla regressione calcolata e premete ↵ Invio.

# interpolazione lineare sulla regressione antitòna calcolata in precedenza

reg <- approx(fdr$x, fdr$yf, xout=c(10, 12, 14, 16, 18, 20, 22, 24), ties="ordered") # interpola i valori specificati

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

#

Le percentuali (%) di soggetti sovrappeso corrispondenti al 10, 12, 14, 16, 8, 100, 20, 22, 24 % di obesi calcolate mediante interpolazione lineare sui dati della regressione antitòna (prima riga di codice) sono mostrate (seconda riga).

> round(cbind(reg$x, reg$y), digits=1) # mostra i dati e i valori interpolati

     [,1] [,2]

[1,]   10 44.8

[2,]   12 39.6

[3,]   14 36.1

[4,]   16 36.1

[5,]   18 36.1

[6,]   20 35.7

[7,]   22 34.4

[8,]   24 34.4

Se copiate e incollate nella Console di R questa riga di codice e premete ↵ Invio

# riporta sul grafico i punti interpolati

#

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

#

potete riportare sul grafico i risultati dell'interpolazione e confermare la correttezza del calcolo.

Vediamo ora in alternativa la regressione calcolata con il metodo previsto nel pacchetto monreg, copiate e incollate nella Console di R questo script e premete ↵ Invio.

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

#

library(monreg) # carica il pacchetto per la regressione

bmi <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

windows() # apre e inizializza una nuova finestra grafica

#

mon <- monreg(bmi$obeso, bmi$sovrappeso, hd=.5, hr=.5, monotonie="antiton") # calcola la regressione antitòna

#

plot(bmi$obeso, bmi$sovrappeso, main="Regressione monotòna (antitòna)", xlab="Soggetti obesi in %", ylab="Soggetti sovrappeso in %") # riporta il grafico con la distribuzione dei dati

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

#

Lo script è quasi identico al precedente, salvo il fatto che la regressione viene calcolata con la funzione monreg() e con l'argomento monotonie="antiton" salvando i risultati nell'oggetto mon dal quale la funzione lines() ricava ascisse (mon$t) e ordinate (mon$estimation) per tracciare la regressione, che ha un andamento meno spigoloso della precedente.

Questo breve addendum allo script consente di effettuare sulla regressione l'interpolazione lineare di alcuni valori a titolo di esempio, per farlo copiatelo e incollatelo nella Console di R e premete ↵ Invio.

# interpolazione lineare sulla regressione antitòna calcolata in precedenza

reg <- approx(mon$t, mon$estimation, xout=c(10, 12, 14, 16, 18, 20, 22, 24), ties="ordered") # interpola i valori specificati

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

#

Le percentuali (%) di soggetti sovrappeso corrispondenti al 10, 12, 14, 16, 8, 100, 20, 22, 24 % di obesi calcolate mediante interpolazione lineare sui dati della regressione antitòna sono ora queste.

> round(cbind(reg$x, reg$y), digits=1) # mostra i dati e i valori interpolati

     [,1] [,2]

[1,]   10 42.5

[2,]   12 38.3

[3,]   14 36.8

[4,]   16 36.3

[5,]   18 36.1

[6,]   20 35.9

[7,]   22 35.6

[8,]   24 34.6

Se si desidera riutilizzare lo script, oltre a tutto il resto anche i valori da interpolare riportati nell'argomento xout=c() devono essere ovviamente adattati al bisogno. Anche in questo caso con una semplice 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") 

#

potete verificare la correttezza del calcolo riportando sul grafico i risultati dell'interpolazione.

----------

[1] Vedere il post Regressione con una funzione monotòna (isotòna).

[2] Nella letteratura statistica anglosassone trovate "antitonic". 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".

[3] Vedere il post Indice di massa corporea (BMI).

[4] 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).

[5] A pagina 14 del Reference manual del pacchetto fdrtool è riportato: "If several identical x values are given as input, the corresponding y values and the weights w are automatically merged, and a warning is issued".

https://cran.r-project.org/web/packages/fdrtool/fdrtool.pdf

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 ed è 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).