Regressione lineare multipla

La regressione lineare semplice (regressione bivariata) costruita sul piano cartesiano

y = a + b·x

può essere estesa a uno spazio n-dimensionale: abbiamo così la regressione lineare multipla nella quale la variabile dipendente y viene fatta dipendere non più da una, bensì da k variabili indipendenti, indicate come x₁, x₂, x₃ … xₖ secondo l'equazione

y = a + b₁·x₁ + b₂·x₂ + b₃·x₃ + ... + bₖ·xₖ

Il risultato è rappresentato quindi da una intercetta a e da tanti coefficienti angolari b (indicati con b₁, b₂, b₃ ... bₖ) quante sono le variabili indipendenti (x₁, x₂, x₃ ... xₖ) che contribuiscono a determinare il valore della variabile dipendente y.

Come esempio impieghiamo il set di dati ais nel quale, tra i dati ematologici rilevati in un gruppo di atleti australiani, abbiamo tre grandezze:

→ la concentrazione degli eritrociti (globuli rossi) espressa in milioni per microlitro di sangue (10⁶/µL) e contenuta nella variabile rcc;

→ il valore ematòcrito cioè la frazione del volume del sangue occupata dagli eritrociti espressa in percentuale (%) e contenuta nella variabile hc;

→ la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue (g/dL) e contenuta nella variabile hg.

Le tre grandezze sono strettamente legate tra loro dal punto di vista biologico, in quanto gli eritrociti hanno ciascuno un certo volume e contengono ciascuno una data quantità di emoglobina: pertanto all'aumentare del numero di eritrociti aumenta la frazione percentuale del volume del sangue occupata dagli eritrociti (il valore ematòcrito) e aumenta la quantità di emoglobina. Mentre la diminuzione dell'ematòcrito e la diminuzione dell'emoglobina sono indice della diminuzione del numero degli eritrociti.

L'esempio è stato scomposto in una serie di script indipendenti - onde facilitare il riutilizzo del codice per analizzare i propri dati - che impiegano i pacchetti aggiuntivi gvlma, car e rgl, che devono essere scaricati dal CRAN. Accertatevi di avere installato inoltre il pacchetto DAAG che contiene i dati impiegati nello script, o in alternativa  procedete come indicato nel post Il set di dati ais. 

Iniziamo con una analisi esplorativa preliminare dei dati, copiate questo primo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - analisi esplorativa preliminare dei dati

#

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

library(car) # carica il pacchetto per la grafica

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#

scatterplotMatrix(~rcc+hc+hg, col="black", pch=20, regLine=list(method=lm, lty=1, lwd=2, col="chartreuse3"), smooth=FALSE, diagonal=list(method="histogram", breaks="FD"), main="Istogrammi e grafici xy con rette di regressione", data=mydata) # istogrammi e grafici xy con rette di regressione

#

Dopo avere caricato i pacchetti necessari, copiato nell'oggetto mydata le tre variabili che ci interessano, aperto e inizializzato una nuova finestra grafica, l'analisi preliminare dei dati viene effettuata con la funzione scatterplotMatrix() nella quale:

→ ~rcc+hc+hg specifica le variabili da rappresentare, che sono la concentrazione degli eritrociti rcc, il valore ematòcrito hc e la concentrazone di emoglobia hg;

→ col="..." specifica il colore dei grafici; 

→ pch=... specifica il tipo di simbolo da impiegare (un piccolo cerchio pieno);

→ per la retta di regressione regLine sono impiegati il modello lineare (lm), una linea continua (lty=1) con uno spessore doppio (lwd=2) e un colore (col="chartreuse3"), mentre l'argomento smooth=FALSE fa si che non venga aggiunta ai grafici la rappresentazione della funzione non lineare prevista di default nel pacchetto car;

→ l'argomento diagonal=list() indica la rappresentazione della variabile riportata nella diagonale, che qui è l'istogramma (method="histogram") con il numero di classi (breaks="FD") calcolato mediante la regola di Freedman-Diaconis. 

Da notare che le espressioni che è possibile impiegare per l'argomento  diagonal=... sono:

→ list(method="boxplot");

→ list(method="density", bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE);

→ list(method="histogram");

→ list(method="oned");

→ list(method="qqplot").

Vista la relazione lineare molto buona che intercorre tra le variabili sembra ragionevole continuare ed effettuare il calcolo della regressione lineare multipla, copiate questo secondo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - calcolo della regressione lineare multipla

#  

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della regressione 

#

coefficients(reglin) # intercetta e coefficienti angolari 

confint(reglin, level=0.95) # intervalli di confidenza al 95% 

#  

La funzione lm(...) viene impiegata per calcolare la regressione lineare multipla specificando con l'espressione rcc~hg+hc la variabile dipendente (rcc), le due variabili indipendenti (hg e hc) e l'equazione 

rcc = a + b₁·hg + b₂·hc

quindi possiamo ricavare con la funzione coefficients() i valori dell'intercetta (a), del coefficiente angolare di hg (b₁) e del coefficiente angolare di hc (b₂):

(Intercept)          hg          hc 

-0.24269563  0.03283747  0.10403395 

e riportare l'equazione della regressione lineare multipla come

rcc = -0.24269563 + 0.03283747·hg + 0.10403395·hc

Successivamente con la funzione confint() sono riportati gli intervalli di confidenza al 95% dell'intercetta e dei due coefficienti angolari:

                  2.5 %     97.5 %

(Intercept) -0.53163317 0.04624192

hg          -0.02459817 0.09027311

hc           0.08267073 0.12539718

Da notare che gli intervalli di confidenza forniscono un test di significatività per intercetta e coefficienti angolari:

→ dato un intervallo di confidenza al 95% che va da -0.53163317 a 0.04624192 (che quindi include il valore 0) l'intercetta non è significativamente diversa da 0, cosa molto logica in quanto al diminuire fino ad azzerarsi dei valori di emoglobina ed ematòcrito deve necessariamente diminuire e azzerarsi il numero degli eritrociti;

→ con un intervallo di confidenza al 95% che va da -0.02459817 a 0.09027311 (che quindi include il valore 0) il coefficiente angolare della concentrazione di emoglobina hg non è significativamente diversa da 0, cosa meno logica in quanto non solo dal punto di vista biologico deve esistere una proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina, ma anche nell'analisi esplorativa preliminare dei dati riportata con lo script precedente risulta evidente la proporzionalità tra numero degli eritrociti e concentrazione dell'emoglobina che è una dato di fatto dal punto di vista biologico;

→ con un intervallo di confidenza al 95% che va da 0.08267073 a 0.12539718 (che quindi non include il valore 0) il coefficiente angolare del valore ematòcrito hc risulta significativamente diverso da 0.

Come detto all'inizio, oltre ad effettuare una analisi esplorativa preliminare dei dati è opportuno sottoporre i risultati ad una valutazione critica a posteriori, che effettuiamo ora sotto forma di una diagnostica numerica della regressione lineare multipla, copiate il terzo script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica numerica della regressione

#

library(gvlma) # carica il pacchetto per la diagnostica numerica 

library(car) # carica il pacchetto per la grafica 

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

summary(gvlma(reglin)) #  test globale per l'assunto di linearità 

outlierTest(reglin) # valore p di Bonferroni per la presenza di dati aberranti (outliers) 

#

Ecco la prima parte del riepilogo dei risultati fornito dalla funzione summary():

Call:

lm(formula = rcc ~ hg + hc, data = mydata)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.55791 -0.11284 -0.00699  0.09950  0.53595 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept) -0.24270    0.14652  -1.656   0.0992 .  

hg           0.03284    0.02913   1.127   0.2609    

hc           0.10403    0.01083   9.603   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1744 on 199 degrees of freedom

Multiple R-squared:  0.8565,    Adjusted R-squared:  0.855 

F-statistic: 593.8 on 2 and 199 DF,  p-value: < 2.2e-16

Il valore di probabilità Pr(>|t|) ottenuto con il test t di Student conferma (ovviamente) i risultati della significatività ricavati dagli intervalli di confidenza riportati sopra, per i quali valgono le considerazioni già fatte. E se la significatività statistica non è coerente con una realtà biologica assodata, dovrà ovviamente essere quest'ultima a prevalere: per cui assumiamo che l'equazione calcolata sia, fino a prova contraria, una descrizione sufficientemente adeguata della relazione di funzione che intercorre tra le tre variabili esaminate.

Che questo sia vero ce lo conferma la funzione gvlma(), acronimo di global validation of linear models assumptions, che fornisce questi risultati:

→ i dati sono distribuiti in modo lineare (Global Stat), assunto accettabile;

→ i dati non presentano asimmetria (Skewness), assunto accettabile;

→ i dati non presentano curtosi (Skewness), assunto accettabile;

→ nell'ambito dei valori osservati la varianza è costante (Heteroscedasticity), assunto accettabile.

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                    Value p-value                   Decision

Global Stat        6.5801 0.15981    Assumptions acceptable.

Skewness           0.5149 0.47304    Assumptions acceptable.

Kurtosis           0.7559 0.38461    Assumptions acceptable.

Link Function      3.9844 0.04592 Assumptions NOT satisfied!

Heteroscedasticity 1.3249 0.24972    Assumptions acceptable.

L'assenza nei dati di asimmetria e di curtosi, e la loro omoschedasticità, sono condizioni sufficienti per garantire la loro adeguata descrizione mediante una regressione lineare. Il primo test (Global Stat) di valutazione globale degli assunti di linearità, che sostanzialmente comprende tutti e tre i test precedenti in forma compatta, li conferma. In questo contesto il test rimanente e l'unico discordante (Link Function) può essere ignorato.

Con la funzione outlierTest() viene  effettuato il test di Bonferroni, che non evidenzia dati anomali (outliers) ma segnala il dato 78 come quello che si discosta maggiormente dagli altri:

No Studentized residuals with Bonferroni p < 0.05

Largest |rstudent|:

   rstudent unadjusted p-value Bonferroni p

78 -3.28855          0.0011917      0.24073

In tema di valutazione critica a posteriori, aggiungiamo ora alla diagnostica numerica anche una diagnostica grafica della regressione lineare multipla, per effettuarla copiate questo quarto script, incollatelo nella Console di R e premete ↵ Invio:

# REGRESSIONE LINEARE MULTIPLA - diagnostica grafica della regressione

#

library(car) # carica il pacchetto per la grafica

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

reglin <- lm(rcc~hg+hc, data=mydata) # calcolo della retta di regressione lineare multipla

#

windows() # apre e inizializza una nuova finestra grafica

qqPlot(reglin, main="Grafico dei quantili per i residui", xlab="t-quantili", ylab="Residui studentizzati") # grafico globale per linearità e dati aberranti

#

windows() # apre e inizializza una nuova finestra grafica

ceresPlots(reglin, ask=FALSE, main="Grafico di Ceres", ylab="Residui di Ceres (rcc)") # grafico di Ceres per la valutazione separata della linearità di hg e hc vs. rcc

#

windows() # apre e inizializza una nuova finestra grafica

plot(mydata$rcc, reglin$fitted.values, xlim=c(4,7), ylim=c(4,7)) # grafico valori originari vs valori stimati

abline(0, 1, lty=2) # linea di indentità y = x

#

La funzione qqPlot() riporta una distribuzione dei residui lineare, evidenziando nel contempo i dati 78 e 161 come quelli che dei quali sarebbe utile controllare la validità, ovvero che sono situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati. 

La funzione ceresPlot() genera due grafici di Ceres separati per emoglobina ed ematòcrito, e consente di  visualizzare una relazione lineare tra ematòcrito ed eritrociti migliore di quella che intercorre tra emoglobina ed eritrociti, confermando in un certo senso i valori di significatività riportati sopra.

La diagnostica grafica della regressione lineare multipla si conclude con questo grafico minimalista, ma che entra nell'essenza dei risultati.

Il grafico, generato con la funzione plot(), riporta in ascisse i valori originari della concentrazione degli eritrociti (mydata$rcc) e in ordinate i corrispondenti valori reglin$fitted.values ricalcolati da hg e hc impiegando l'equazione della regressione lineare multipla ove

reglin$fitted.values = -0.24269563 + 0.03283747·hg + 0.10403395·hc

La linea di identità y = x che viene riportata evidenzia, senza la necessità di ulteriori statistiche, che l'informazione acquisita misurando la concentrazione dell'emoglobina e l'ematòcrito può essere compressa in un unico valore, che fornisce una stima mediamente adeguata della concentrazione degli eritrociti.

Nota: con "mediamente adeguata" si intende che, data la misura diretta della concentrazione degli eritrociti riportata in ascisse, la sua misura indiretta - ottenuta misurando emoglobina ed ematòcrito ed applicando la regressione lineare multipla - comporta una perdita di informazione [aumento dell'incertezza della misura] che si riflette nella "dispersione media" dei valori sull'asse delle ordinate. Se non ci fosse perdita di informazione i punti sarebbero perfettamente allineati sulla retta y = x.

Ora, visto che con tre variabili la cosa è ancora fattibile, non resta che aggiungere un grafico tridimensionale (3D).

Copiate quest'ultimo script, incollatelo nella Console di R e premete ↵ Invio:

# GRAFICO 3D ANIMATO con il pacchetto rgl

#  

library(rgl) # carica il pacchetto per la grafica 3D

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

mydata <- ais[,c(1,3,4)] # le colonne 1, 3, 4 contengono le variabili rcc, hc, hg

windows() # apre e inizializza una nuova finestra grafica

#  

# grafico 3D animato

#  

plot3d(mydata$rcc, y=mydata$hg, z=mydata$hc, type="p", col="red", size=3) # realizza il grafico

play3d(spin3d(axis=c(0,0,1), rpm=10), duration=10) # anima il grafico

#  

Dopo avere caricato il pacchetto necessario con il comando library(rgl) il grafico viene realizzato con la funzione plot3d().

Quindi con la funzione play3d() viene realizzata una animazione della durata di 10 secondi (duration=10) con i parametri specificati dalla funzione spin3d() e cioè l'asse sul quale fare ruotare l'animazione, nel nostro caso  l'asse delle z (0,0,1), e il numero di rotazioni al minuto (rpm=10).

Quando l'animazione che compare nella finestra "RGL device" aperta nella Console di R dallo script si ferma, se “afferrate” il grafico 3D facendo click con il tasto sinistro del mouse e tenendolo premuto senza rilasciarlo potete ruotarlo a vostro piacimento muovendo il mouse. 

Anche se non aggiunge nulla di nuovo, il grafico [1] ci conforta documentando la linearità della relazione tra le tre variabili.

Conclusione: mettere in relazione tra di loro troppe variabili contemporaneamente può portare a risultati di difficile comprensione. Tuttavia una analisi esplorativa preliminare dei dati e una attenta valutazione critica a posteriori dei risultati possono aiutare a identificare i casi nei quali l'impiego della regressione lineare multipla può essere appropriato.

----------

[1] Trovate altre modalità di rappresentazione sotto forma di grafici tridimensionali nel post Grafici 3D.

Sensibilità, specificità e valore predittivo

Nel post Teorema di Bayes e diagnosi medica abbiamo visto che è possibile impiegare il teorema per rispondere alle due domande: quale probabilità ha di essere affetto da una specifica malattia un soggetto cui un dato test per quella specifica malattia è risultato positivo? quale probabilità ha di essere sano un soggetto cui il test è risultato negativo?

Nel post Curve ROC e valore soglia abbiamo visto, dati i risultati di un test per la diagnosi di una malattia X, come impiegare il pacchetto pROC per tracciare la curva ROC [1] e individuare il valore soglia che consente di classificare correttamente il maggior numero possibile di soggetti sani (VN) e di classificare correttamente il maggior numero possibile di soggetti malati (VP), minimizzando sia la possibilità di sbagliare classificando dei soggetti malati come sani (FN) sia la possibilità di sbagliare classificando dei soggetti sani come malati (FP). Dato il valore soglia individuato, risultano conseguentemente fissate la sensibilità e la specificità del test.

Ora impieghiamo il pacchetto epiR per calcolare le probabilità che rispondono alle due domande e cioè:

→ il valore predittivo di un test positivo, ovvero la probabilità di essere malato per un soggetto che presenta un test positivo, che risponde alla prima domanda

→ il valore predittivo di un test negativo, ovvero la probabilità di essere sano per un soggetto che presenta un test negativo, che risponde alla seconda domanda.

Si stabilisce se il test è positivo o negativo impiegando il valore soglia individuato, quindi si calcola il valore predittivo del test impiegando:

→ la sensibilità del test, cioè la probabilità che il test sia positivo in un soggetto malato;

→ la specificità del test, cioè la probabilità che il test sia negativo in un soggetto sano;

→ la prevalenza, cioè la probabilità della malattia.

La sensibilità e la specificità sono caratteristiche intrinseche al test, dipendono dal valore soglia individuato e sono indipendenti dalla prevalenza, che è invece una caratteristica della popolazione alla quale il test viene applicato.

Notare che la prevalenza, cioè la probabilità della malattia nella popolazione alla quale il test viene applicato, è la probabilità a priori della malattia, ovvero la probabilità per un dato paziente di essere malato (o sano) prima di avere eseguito il test. Il teorema di Bayes impiega la probabilità a priori (nel nostro caso la prevalenza della malattia), insieme a sensibilità e specificità del test, per calcolare la probabilità a posteriori, ovvero la probabilità per un dato paziente di essere malato (o sano) dopo avere eseguito il test (il valore predittivo del test positivo o del test negativo).

Notare inoltre che nella (iii) del post Teorema di Bayes e diagnosi medica:

→ P(T+|M+)/P(T+) è una costante e rappresenta l'evidenza fornita dai fatti (nel nostro caso il risultato del test per la diagnosi della malattia)

→ P(M+) è l'informazione/probabilità a priori

→ P(M+|T+) è l'informazione/probabilità a posteriori

Pertanto il teorema può essere riscritto come

probabilità a posteriori = evidenza · probabilità a priori

espressione che sintetizza l'essenza del teorema di Bayes: che consente di aggiornare l'informazione in nostro possesso, la probabilità a priori, sulla base dell'evidenza fornita dai fatti.

Mentre i calcoli eseguiti impiegando i valori di probabilità sono riportati nell'appendice al termine del post, vediamo qui i calcoli eseguiti impiegando il numero di casi osservati, organizzati in questo modo

	Malattia +	Malattia -
Test +	VP	FP
Test -	FN	VN

e definiti come [2]:

→ veri positivi (VP), i soggetti malati con il test positivo;

→ falsi positivi (FP), i soggetti sani con il test positivo;

→ falsi negativi (FN), i soggetti malati con il test negativo;

→ veri negativi (VN), i soggetti sani con il test negativo.

Supponiamo ora (primo esempio) che il test per una data malattia sia stato eseguito in 800 soggetti sani, dipendenti di un Ospedale, impiegati come controlli, e in 800 malati, che sono stati raccolti in numero così elevato in quanto nell'Ospedale esiste un Reparto che si occupa specificamente di questa malattia. Questa è la tabella con i risultati del test

	Malattia +	Malattia -
Test +	700 (VP)	44 (FP)
Test -	100 (FN)	756 (VN)

I calcoli sono effettuati mediante il pacchetto aggiuntivo epiR [3] che deve essere preventivamente scaricato dal CRAN. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# SENSIBILITA', SPECIFICITA' E VALORE PREDITTIVO

# 

#

library(epiR) # carica il pacchetto

cells <- c(700,44,100,756) # genera l'array cells con i valori numerici contenuti nelle celle

rnames <- c("Test +", "Test -") # genera l'array rnames con i nomi delle righe

cnames <- c("Malattia +", "Malattia -") # genera l'array cnames con i nomi delle colonne

mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la matrice dati mymatrix a partire dagli array cells, rnames, cnames

mymatrix # mostra la matrice dati mymatrix

mystat <- epi.tests(mymatrix, conf.level=0.95) # calcola le statistiche

summary(mystat) # mostra le statistiche

#

Con la funzione c() [4] sono generati gli array che mediante l'operatore di assegnamento <- vengono denominati cells, rnames, cnames e che ora contengono rispettivamente i dati (cells), i nomi delle righe (rnames) e i nomi delle colonne (cnames).

Quindi mediante la funzione matrix() [5] gli array sono combinati nella matrice mymatrix di 2 righe (nrow=2) per due colonne (ncol=2) inserendo i valori per riga (byrow=TRUE) e quindi da sinistra a destra e dall'alto in basso ottenendo la matrice desiderata alla quale con l'ultimo argomento (dimnames=list(rnames, cnames)) sono assegnati i nomi alle righe e alle colonne.

Questo è il contenuto dell'oggetto mymatrix riportato da R con il comando mymatrix:

> mymatrix # mostra la matrice dati mymatrix

       Malattia + Malattia -

Test +        700         44

Test -        100        756

Impiegando la funzione epi.test() sono infine calcolate sui dati contenuti in mymatrix le statistiche, ciascuna con l'intervallo di confidenza al 95% (conf.level=0.95), e i risultati sono riportati con la funzione summary():

> summary(mystat) # mostra le statistiche

   statistic          est       lower        upper

1         ap   0.46500000  0.44033136   0.48979728

2         tp   0.50000000  0.47520704   0.52479296

3         se   0.87500000  0.85006839   0.89712668

4         sp   0.94500000  0.92686534   0.95975609

5    diag.ac   0.91000000  0.89490175   0.92357301

6    diag.or 120.27272727 83.14513765 173.97925284

7       nndx   1.21951220  1.16702078   1.28711106

8     youden   0.82000000  0.77693373   0.85688277

9     pv.pos   0.94086022  0.92141771   0.95670377

10    pv.neg   0.88317757  0.85974563   0.90393102

11    lr.pos  15.90909091 11.92293523  21.22792488

12    lr.neg   0.13227513  0.11003359   0.15901245

13    p.rout   0.53500000  0.51020272   0.55966864

14     p.rin   0.46500000  0.44033136   0.48979728

15    p.tpdn   0.05500000  0.04024391   0.07313466

16    p.tndp   0.12500000  0.10287332   0.14993161

17    p.dntp   0.05913978  0.04329623   0.07858229

18    p.dptn   0.11682243  0.09606898   0.14025437

Le prime 12 statistiche calcolate sono quelle che ci interessano, riportiamo (tra parentesi) il significato della sigla, le formule con cui sono calcolate e il risultato numerico ottenuto (est) con il limite inferiore (lower) e il limite superiore (upper) del relativo intervallo di confidenza al 95% (le statistiche da 13 a 18 riguardano l'outcome cioè il risultato, e non le consideriamo).

ap (prevalenza apparente, soggetti con il test positivo)

(VP+FP)/(VP+FP+FN+VN) = 0.4650000 (0.4403314 - 0.4897973)

tp (prevalenza reale, soggetti con la malattia)

(VP+FN)/(VP+FP+FN+VN) = 0.5000000 (0.4752070 - 0.5247930)

se (sensibilità, positività nei malati)

VP/(VP+FN) = 0.8750000 (0.8500684 - 0.8971267)

sp (specificità, negatività nei sani)

VN/(VN+FP) = 0.9450000 (0.9268653 - 0.9597561)

diag.ac (accuratezza diagnostica)

(VP+VN)/(VP+FP+FN+VN) = 0.9100000 (0.8949017 - 0.9235730)

diag.or (odd ratio)

rapporto LR+/LR- = 120.2727273 (83.1451376 - 173.9792528)

nndx (numero necessario per la diagnosi)

1/(sensibilità - (1 - specificità)) = 1.2195122 (1.1670208 - 1.2871111)

youden (indice di Youden)

sensibilità + specificità - 1 = 0.8200000 (0.7769337 - 0.8568828)

pv.pos (valore predittivo di un test positivo)

VP/(VP+FP) = 0.9408602 (0.9214177 - 0.9567038)

pv.neg (valore predittivo di un test negativo)

VN/(VN+FN) = 0.8831776 (0.8597456 - 0.9039310)

lr.pos (rapporto di verosimiglianza LR+ per un test positivo) [6]

(VP/(VP+FN))/(FP/(FP+VN)) = 15.9090909 (11.9229352 - 21.2279249)

lr.neg (rapporto di verosimiglianza LR- per un test negativo) [6]

(FN/(VP+FN))/(VN/(FP+VN)) = 0.1322751 (0.1100336 - 0.1590125)

La conclusione è che il valore predittivo di un test positivo, ovvero la probabilità di essere malato per un soggetto al quale il test è risultato positivo, è uguale al 94.1% e il valore predittivo di un test negativo, ovvero la probabilità di essere sano per un soggetto al quale il test è risultato negativo, è uguale all'88.3%.

Nel caso appena esaminato la prevalenza della malattia è del 50%, un valore che è dipeso dalle modalità di raccolta dei dati: per valutare adeguatamente il test è stato raccolto il maggior numero possibile di malati, un fatto che non riflette la reale diffusione della malattia. Ora supponiamo che il test venga applicato per valutare lo stato di salute di un individuo della popolazione generale, nella quale la reale diffusione della malattia, la sua prevalenza, è del 2%. I dati (secondo esempio) sono questi:

	Malattia +	Malattia -
Test +	14 (VP)	44 (FP)
Test -	2 (FN)	756 (VN)

Copiate questo script, identico al precedente tranne che per il fatto che contiene i nuovi dati, incollatelo nella Console di R e premete ↵ Invio:

# SENSIBILITA', SPECIFICITA' E VALORE PREDITTIVO

# 

#

library(epiR) # carica il pacchetto

cells <- c(14,44,2,756) # genera l'array cells con i valori numerici contenuti nelle celle

rnames <- c("Test +", "Test -") # genera l'array rnames con i nomi delle righe

cnames <- c("Malattia +", "Malattia -") # genera l'array cnames con i nomi delle colonne

mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la matrice dati mymatrix a partire dagli array cells, rnames, cnames

mymatrix # mostra la matrice dati mymatrix

mystat <- epi.tests(mymatrix, conf.level=0.95) # calcola le statistiche

summary(mystat) # mostra le statistiche

#

Questi sono i dati contenuti nell'oggetto mymatrix

> mymatrix # mostra la matrice dati mymatrix

            Malattia + Malattia -

Test +         14         44

Test -          2        756

e questi sono, riportati con la funzione summary(), i risultati calcolati con la funzione epi.test(), ciascuno con l'intervallo di confidenza al 95% (conf.level=0.95)

> summary(mystat) # mostra le statistiche

   statistic           est         lower         upper

1         ap   0.071078431  0.0544119039   0.090919329

2         tp   0.019607843  0.0112480876   0.031647030

3         se   0.875000000  0.6165237632   0.984486396

4         sp   0.945000000  0.9268653430   0.959756088

5    diag.ac   0.943627451  0.9255197857   0.958436211

6    diag.or 120.272727273 26.5044169172 545.778047894

7       nndx   1.219512195  1.0590499968   1.840301892

8     youden   0.820000000  0.5433891062   0.944242484

9     pv.pos   0.241379310  0.1387012000   0.371689716

10    pv.neg   0.997361478  0.9905015414   0.999680303

11    lr.pos  15.909090909 11.3036587611  22.390907130

12    lr.neg   0.132275132  0.0361754961   0.483661941

13    p.rout   0.928921569  0.9090806715   0.945588096

14     p.rin   0.071078431  0.0544119039   0.090919329

15    p.tpdn   0.055000000  0.0402439118   0.073134657

16    p.tndp   0.125000000  0.0155136038   0.383476237

17    p.dntp   0.758620690  0.6283102838   0.861298800

18    p.dptn   0.002638522  0.0003196972   0.009498459

Si nota immediatamente che quando la prevalenza della malattia diminuisce dal 50% al 2% (precisamente diventa 1.96%)

                 est       lower        upper

tp        0.5000000   0.4752070    0.5247930 [primo esempio]

tp       0.01960784  0.01124809   0.03164703 [secondo esempio]

il valore predittivo del test positivo diminuisce passando dal 94.1% al 24.1%

                 est       lower        upper

pv.pos    0.9408602   0.9214177    0.9567038  [primo esempio]

pv.pos    0.24137931  0.13870120   0.37168972 [secondo esempio]

e il valore predittivo del test negativo aumenta passando da 88.3% a 99.7%

                 est       lower        upper

pv.neg    0.8831776   0.8597456    0.9039310  [primo esempio]

pv.neg    0.99736148  0.99050154   0.99968030 [secondo esempio]

Mentre sensibilità e specificità non cambiano al cambiare della popolazione alla quale un test diagnostico viene applicato, in quanto sono caratteristiche intrinseche al test, la prevalenza della malattia non dipende dal test, ma è una caratteristica della popolazione alla quale il test viene applicato, e può cambiare notevolmente. Per esempio la prevalenza della malaria in Italia è diversa da quella che si ha nei paesi dell'Africa equatoriale; la prevalenza della malattia reumatica nei pazienti che accedono all'ambulatorio di un medico generico è differente dalla prevalenza della malattia reumatica nei pazienti ricoverati in un Reparto ospedaliero specializzato nella cura di questa malattia; la prevalenza del cancro al polmone è differente nei soggetti fumatori rispetto ai soggetti non fumatori; e così via. E questo cambia il valore predittivo del test in quanto al diminuire della prevalenza della malattia, cioè della probabilità a priori di malattia, diminuisce il valore predittivo del test positivo (la probabilità a posteriori di malattia in un soggetto cui il test è risultato positivo), e aumenta il valore predittivo del test negativo (la probabilità a posteriori di non-malattia in un soggetto cui il test è risultato negativo) [7].

Questo è un punto cruciale, e rende conto del perché in molte situazioni di bassa prevalenza un test sia utile, e venga impiegato, soprattutto per escludere una malattia. Mentre il teorema di Bayes mostra (e non solo in campo medico) come e quanto l'incertezza delle nostre conclusioni può essere ridotta grazie all'informazione/evidenza fornita dai fatti, e R fornisce opportuni strumenti di calcolo.

----------

[1] Abbiamo anche visto in un altro post come sovrapporre due Curve ROC sovrapposte.

[2] Gerhardt W, Keller H. Evaluation of test data from clinical studies. Scand J Clin Lab Invest, 46, 1986 (supplement 181).

https://www.tandfonline.com/doi/pdf/10.1080/00365518709168461

[3] Package ‘epiR’.

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

[4] Digitate help(c) nella Console di R per la documentazione della funzione c().

[5] Digitate help(matrix) nella Console di R per la documentazione della funzione matrix.

[6] Deeks JJ, Altman DG. Diagnostic tests 4: likelihood ratios. BMJ 2004;329;168-169.

https://www.bmj.com/content/329/7458/168.full

[7] Per i temi trattati nei post sul teorema di Bayes vedere Altman, Federspil, Galen, Gerhardt, Gigerenzer, Gill, Motterlini e Nordenström nella pagina Bibliografia.

APPENDICE

Si vuole dimostrare che i risultati ottenuti sopra con i conteggi dei casi sono identici a quelli ottenuti impiegando il teorema di Bayes con i valori di probabilità e le espressioni (iv) per il valore predittivo del test positivo e (v) per il valore predittivo del test negativo (vedi post Teorema di Bayes e diagnosi medica).

Primo esempio riportato sopra: abbiamo sensibilità del test 0.875, specificità del test 0.945 e prevalenza della malattia 0.50.

Valore predittivo del test positivo

                              P(T+|M+) · P(M+)

P(M+|T+) = ———————————————————— (iv)

                   P(T+|M+) · P(M+) + P(T+|M-) · P(M-)

che si legge come

                                        sensibilità · prevalenza

P(M+|T+) = ——————————————————————————————

                   sensibilità · prevalenza + (1 – specificità) · (1 - prevalenza)

quindi il valore predittivo del test positivo è uguale a

                              0.875 · 0.50

P(M+|T+) = ———————————————— = 0.9408602 

                   (0.875 · 0.50) + (0.055 · 0.50)

Valore predittivo del test negativo

                               P(T-|M-) · P(M-)

P(M-|T-) = ——————————————————— (v)

                 P(T-|M-) · P(M-) + P(T-|M+) · P(M+)

che si legge come

                                   specificità · (1 – prevalenza)

P(M-|T-) = ——————————————————————————————

                 specificità · (1 – prevalenza) + (1 – sensibilità) · prevalenza

quindi il valore predittivo del test negativo è uguale a

                            0.945 · 0.50

P(M-|T-) = ———————————————— = 0.8831776 

                 (0.945 · 0.50) + (0.125 · 0.50)

Secondo esempio riportato sopra: abbiamo sensibilità del test 0.875, specificità del test 0.945 e prevalenza della malattia 0.01969784.

Valore predittivo del test positivo

                                     0.875 · 0.01960784

P(M+|T+) = ———————————————————————— = 0.24137928 

                   (0.875 · 0.01960784) + (0.055 · 0.98039216)

Valore predittivo del test negativo

                                    0.945 · 0.98039216

P(M-|T-) = ———————————————————————— = 0.99736148 

                 (0.945 · 0.98039216) + (0.125 · 0.01960784)

I risultati del valore predittivo ottenuti impiegando i conteggi dei casi sono identici a quelli ottenuti impiegando il teorema di Bayes con i valori di probabilità, come volevasi dimostrare.

----------