2012-11-26

検査を使って判断する

PPV NPV ベイズ事前確率事後確率 R

こちらで臨床検査の感度特異度の話、そのベイズ的使用のためにきおくするべきことの話が出ている
興味深い
現実医療現場では、「状況(事前確率の高低カテゴリ)」「使用検査」「判断したいこと＝確定診断したい・ルールアウトしたい」の組み合わせに対して、「役に立つ」「役に立たない」の２値のどちらであるか、を記憶しているのだと思う
検診で腫瘍マーカーを測って、確定診断する、というのは「役に立たない」
検診で腫瘍マーカーを測って、ルールアウトする、というのも「役に立たない」
不明熱で抗核抗体を測って、SLEの確定診断する、というのは「大して役に立たない」
不明熱で降格固体を測って、SLEのルールアウトをする、というのは「かなり役に立つ」
など
以下の図は、
- 事前確率が0.1のときには確定診断できるような感度特異度の組み合わせはないようなものであり、ルールアウトすることは感度がそこそこあればよく
- 事前確率が0.5のときには、感度のよさ、特異度のよさがそれぞれ、効くが、その値が結構よいことが必要
- 事前確率が0.9のときには、特異度がそこそこあれば確定しやすい
を表している
検診・腫瘍マーカー・確定診断は図の左上の白いところ
検診・腫瘍マーカー・ルールアウトは右上の白いところ
不明熱・ANA・確定診断は、左中段の白いところ
不明熱・ANA・ルールアウトは右中段の黒いところ
と、覚えている…という感じ

# my clinical tests
sens <- seq(from = 0.01,to = 0.99, by = 0.01)
spec <- seq(from = 0.01,to = 0.99, by = 0.01)

pre <- c(0.1,0.5,0.9)
post <- c(0.1,0.5,0.9)

ss <- expand.grid(sens,spec)

PPV <- NPV <- matrix(0,length(ss[,1]),length(pre))
for(i in 1:length(ss[,1])){
	for(j in 1:length(pre)){
		x <- pre[j] * ss[i,1]
		y <- pre[j] * (1-ss[i,1])
		z <- (1-pre[j]) * (1-ss[i,2])
		w <- (1-pre[j]) * ss[i,2]
		PPV[i,j] <- x/(x+z)
		NPV[i,j] <- w/(y+w)
	}
}
par(mfcol=c(1,2))
matplot(PPV,type="l")
matplot(NPV,type="l")

library(rgl)


par(mfcol=c(3,2))


# Pre 0.1 でPPV>=0.95のSens,Specは
ok <- which(PPV[,1] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.1,PPV")
# Pre 0.5 でPPV>=0.95のSens,Specは
ok <- which(PPV[,2] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.5,PPV")
# Pre 0.9 でPPV>=0.95のSens,Specは
ok <- which(PPV[,3] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.9,PPV")



# Pre 0.1 でNPV>=0.95のSens,Specは
ok <- which(NPV[,1] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.1,NPV")

# Pre 0.5 でNPV>=0.95のSens,Specは
ok <- which(NPV[,2] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.5,NPV")


# Pre 0.9 でNPV>=0.95のSens,Specは
ok <- which(NPV[,3] >= 0.95)
plot(ss[ok,],xlim = c(0,1),ylim=c(0,1),xlab="sens",ylab="spec",main="Pre=0.9,NPV")


######## 3d plot

# Pre が0.1のとき、PPVを高めるには
col <- rep(1,length(PPV[,1]))
col[which(PPV[,1]>0.5)] <- 2
plot3d(cbind(ss,NPV[,1]),col = col,main = "pre=0.1",xlab="sens",ylab="spec",zlab="post")

# Pre が0.1のとき、NPVを高めるには
col <- rep(1,length(NPV[,1]))
col[which(NPV[,1]>0.95)] <- 2
plot3d(cbind(ss,NPV[,1]),col = col,main = "pre=0.1",xlab="sens",ylab="spec",zlab="post")

# Pre が0.5のとき、PPVを高めるには、Specの高さが重要
col <- rep(1,length(PPV[,2]))
col[which(PPV[,2]>0.95)] <- 2
plot3d(cbind(ss,PPV[,2]),col = col,main = "pre=0.5",xlab="sens",ylab="spec",zlab="post")

# Pre が0.5のとき、NPVを高めるには、Sensの高さが重要
col <- rep(1,length(NPV[,2]))
col[which(NPV[,2]>0.95)] <- 2
plot3d(cbind(ss,NPV[,2]),col = col,main = "pre=0.5",xlab="sens",ylab="spec",zlab="post")


# Pre が0.9のとき、PPVを高めるには
col <- rep(1,length(PPV[,3]))
col[which(PPV[,3]>0.5)] <- 2
plot3d(cbind(ss,NPV[,3]),col = col,main = "pre=0.9",xlab="sens",ylab="spec",zlab="post")

# Pre が0.9のとき、NPVを高めるには
col <- rep(1,length(NPV[,3]))
col[which(NPV[,3]>0.95)] <- 2
plot3d(cbind(ss,NPV[,3]),col = col,main = "pre=0.9",xlab="sens",ylab="spec",zlab="post")