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\
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
