Efficienza e accuratezza diagnostica [2]

Questo codice riprende esattamente i contenuti del post Efficienza e accuratezza diagnostica [1] al quale si rimanda per i commenti impiegando esclusivamente le funzioni del pacchetto base di R e riportando nella grafica i kernel density plot in sostituzione degli istogrammi. 

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

# EFFICIENZA E ACCURATEZZA DIAGNOSTICA con le funzioni del pacchetto base di R

#

mydata <- read.table("c:/Rdati/bayes.txt", header=TRUE, col.names=c("classe","valore"), sep=";", dec=",", allowEscapes=TRUE) # importa i dati

#

# effettua i conteggi dei casi e la loro classificazione in base ai valori soglia

#

t <- data.frame(table(mydata)) # crosstabulazione dei dati

vs <- strtoi(subset(t$valore, t$classe==0)) # valori soglia

s.vs <- subset(t$Freq, t$classe==0) # numero di soggetti sani per ciascun valore soglia

m.vs <- subset(t$Freq, t$classe==1) # numero di soggetti malati per ciascun valore soglia

n <- length(vs) # numero di valori soglia

# inizializza le variabili per le successive sommatorie

vp=rep(0, n)

fp=rep(0, n)

fn=rep(0, n)

vn=rep(0, n)

# calcola veri positivi, falsi positivi, falsi negativi e veri negativi per ciascun valore soglia

for (i in 1:n) {

vp[i] <- sum(m.vs[i:n])

fp[i] <- sum(s.vs[i:n])

fn[i] <- sum(m.vs[1:i-1])

vn[i] <- sum(s.vs[1:i-1])

}

#

t.cont <- data.frame(vs, s.vs, m.vs, vp, fp, fn, vn) # genera la tabella con i valori soglia e i relativi conteggi

names(t.cont) <- c("soglia", "sani", "malati", "VP", "FP", "FN", "VN") # assegna i nomi alle variabili/colonne

t.cont # mostra la tabella

#

# calcola le statistiche per ciascun valore soglia

#

sens <- vp/(vp+fn) # calcola sensibilità

spec <- vn/(vn+fp) # calcola specificità

effi <- sens*spec # calcola efficienza diagnostica

vptp <- vp/(vp+fp) # calcola valore predittivo del test positivo

vptn <- vn/(vn+fn) # calcola valore predittivo del test negativo

accu <- (vp+vn)/(vp+fp+fn+vn) # calcola accuratezza diagnostica

#

t.stat <- data.frame(vs, sens, spec, effi, vptp, vptn, accu) # genera la tabella con i valori soglia e le corrispondenti statistiche

names(t.stat) <- c("soglia", "sensibilità", "specificità", "efficienza", "vptp", "vptn", "accuratezza") # assegna i nomi alle variabili/colonne

t.stat # mostra la tabella

#

paste("Massima efficienza diagnostica ", round(max(effi), digits=3), " per un valore soglia di ", vs[which.max(effi)], sep="") # riporta la massima efficienza diagnostica e il corrispondente valore soglia

paste("Massima accuratezza diagnostica ", round(max(accu), digits=3), " per un valore soglia di ", vs[which.max(accu)], sep="") # riporta la massima accuratezza diagnostica e il corrispondente valore soglia

#

# traccia kernel density-plot e sovrappone curve di sensibilità, specificità ed efficienza diagnostica

#

sani <- subset(mydata$valore, mydata$classe==0) # sottoinsieme con i valori dei soggetti sani per tracciare i kernel density-plot

malati <- subset(mydata$valore, mydata$classe==1) # sottoinsieme con i valori dei soggetti malati per tracciare i kernel density-plot

#

windows() # apre una nuova finestra

par(mar=c(5,4,4,5)+.1) # margini per fare spazio alla scala dell'asse delle y sulla destra

#

plot(density(sani), xlim=c(0,250), xaxp = c(0, 250, 5), ylim=c(0,0.04), yaxp=c(0, 0.04, 4), lty=1, lwd=1, col="lightgreen", main="Sensibilità, specificità, efficienza diagnostica", xlab="Risultato del test", ylab="Frequenza") # density-plot sani

par(new=TRUE) # il prossimo plot viene sovrapposto

plot(density(malati), xlim=c(0,250), xaxp = c(0, 250, 5), ylim=c(0,0.04), yaxp=c(0, 0.04, 4), lty=1, lwd=1, col="red", main="Sensibilità, specificità, efficienza diagnostica", xlab="Risultato del test", ylab="Frequenza") # density-plot malati

#

par(new=TRUE) # sovrappone curva sensibilità

plot(vs, sens, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch= 20, cex=0.5, col="blue", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva

axis(4) # riporta un secondo asse delle y sulla destra per i valori di probabilità

mtext("Probabilità p", side=4, line=3) # aggiunge la legenda all'asse delle y di destra

#

par(new=TRUE) # sovrappone curva specificità

plot(vs, spec, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch=20, cex=0.5, col="mediumorchid", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva

#

par(new=TRUE) # sovrappone curva efficienza diagnostica

plot(vs, effi, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch=20, cex=0.5, col="orange", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva

#

legend(140, 0.97, legend=c("Sani", "Malati", "Sensibilità", "Specificità", "Efficienza diagnostica"), col=c("green", "red", "blue", "mediumorchid", "orange"), lty=c(1,1,2,2,2), lwd=c(1,1,1,1,1), cex=0.7) # aggiunge al grafico la legenda

#

# traccia kernel density plot e sovrappone curve di valore predittivo ed accuratezza diagnostica

#

windows() # apre una nuova finestra

par(mar=c(5,4,4,5)+.1) # margini per fare spazio alla scala dell'asse delle y sulla destra

#

plot(density(sani), xlim=c(0,250), xaxp = c(0, 250, 5), ylim=c(0,0.04), yaxp=c(0, 0.04, 4), lty=1, lwd=1, col="lightgreen", main="Valore predittivo e accuratezza diagnostica", xlab="Risultato del test", ylab="Frequenza") # density-plot sani

par(new=TRUE) # il prossimo plot viene sovrapposto

plot(density(malati), xlim=c(0,250), xaxp = c(0, 250, 5), ylim=c(0,0.04), yaxp=c(0, 0.04, 4), lty=1, lwd=1, col="red", main="Valore predittivo e accuratezza diagnostica", xlab="Risultato del test", ylab="Frequenza") # density-plot malati

#

par(new=TRUE) # sovrappone curva valore predittivo del test positivo

plot(vs, vptp, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch= 20, cex=0.5, col="blue", xaxt="n", yaxt="n", xlab="",ylab="") # traccia la curva

axis(4) # riporta l'asse delle y sulla destra

mtext("Probabilità p", side=4, line=3) # aggiunge legenda asse y

#

par(new=TRUE) # sovrappone curva valore predittivo del test negativo

plot(vs, vptn, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch=20, cex=0.5, col="mediumorchid", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva

#

par(new=TRUE) # sovrappone curva accuratezza diagnostica

plot(vs, accu, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=2, lwd=1, pch=20, cex=0.5, col="orange", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva

#

legend(135, 0.98, legend=c("Sani", "Malati", "Valore predittivo del test +", "Valore predittivo del test -", "Accuratezza diagnostica"), col=c("green", "red", "blue", "mediumorchid", "orange"), lty=c(1,1,2,2,2), lwd=c(1,1,1,1,1), cex=0.7) # aggiunge al grafico la legenda

#