世界で一番わかりやすいROC解析

MikuHatsune2012-12-27

いや、世界で一番わかりやすいかはどうだろう…

library(pROC)
n <- 30000 #判定されるものの数
a1 <- rnorm(n, mean=0) #陰性のもの
a2 <- rnorm(n, mean=2, sd=2) #陽性のもの
diagnosis <- rep(0:1, each=n) #各々の判定フラグ

roc1 <- roc(diagnosis, c(a1, a2))
sen_spe <- rbind(roc1$sensitivities, roc1$specificities) #各々の閾値での感度と特異度
ss0 <- which.max(colSums(sen_spe)) #ROC解析上もっとも感度と特異度が良い感じになるカットオフの番地
thres <- roc1$thresholds[ss0 + 1] #そのときのカットオフ

plot(roc1)
points(sen_spe[2, ss0], sen_spe[1, ss0], pch=16, col=2) #もっとも感度と特異度が良い感じになるときのカットオフでやったときの感度と特異度

ROCplot <- function(ss0, thres){
	d1 <- density(a1)
	d2 <- density(a2)
	xlim <- range(d1$x, d2$x)
	ylim <- range(d1$y, d2$y)
	ylim <- range(ylim, -ylim)
	plot(1, type="n", xlim=xlim, ylim=ylim, xlab="value", ylab="density")
	lines(d1$x, d1$y)
	lines(d2$x, -d2$y)
	abline(h=0, lty=2)
	cf <- 1
	legend("topleft", legend=c("TN", round(sen_spe[2, ss0], 3)), bty="n", cex=cf)
	legend("topright", legend=c("FP", round(1-sen_spe[2, ss0], 3)), bty="n", cex=cf)
	legend("bottomleft", legend=c("FN", round(1-sen_spe[1, ss0], 3)), bty="n", cex=cf)
	legend("bottomright", legend=c("TP", round(sen_spe[1, ss0], 3)), bty="n", cex=cf)
	abline(v=thres, lty=2)
}

par(mfrow=c(1, 2))
#for(i in head(seq(roc1$thresholds), -1L)){
for(i in floor(seq(1, 60000, length=100))){
	par(mfrow=c(1, 2))
	ROCplot(ss0=i, thres=roc1$thresholds[i + 1])
	abline(v=thres, lty=2, col=4)
	plot(roc1)
	points(sen_spe[2, i], sen_spe[1, i], pch=16, col=2, cex=2)
	points(sen_spe[2, ss0], sen_spe[1, ss0], pch=16, col=4, cex=2)
}