診断基準の項目数と必要項目数

  • こちらから
  • 診断基準というものを使って、診断することがある
  • もっとも単純には、N個の項目を立て、そのうちn項目が陽性ならば、その疾患であると診断しよう、というようなものである
  • ここで、さらなる単純化をして、N項目について、「真に疾患であるときの陽性率p」が全項目にわたって等しく、「真に疾患でないときの陽性率q」が全項目にわたって等しいという場合を考えてみる
  • N項目のうち、i項目が陽性でN-i項目が陰性である確率は、
    • 疾患の場合\begin{pmatrix} N \\i \end{pmatrix}p^i (1-p)^{N-i}
    • 非疾患の場合\begin{pmatrix} N \\i\end{pmatrix} q^i (1-q)^{N-i}
  • プレテスト・プロバリティは診断基準の感度・特異度には影響しないが、PPV,NPVには影響する
  • いくつかのことが言える
    • 単独のテストとしては感度・特異度・PPV・NPVからして、大したものではないが、それらを合わせることで、
    • 感度・特異度の両方がよい場合が出現すること
    • 特に、Nがかなり大きい場合には、nの値がNの半分の前後で、感度・特異度ともにほぼ1であるようなnが多数出現しうることもわかる
    • 感度・特異度のよくなる場合は、n/Nの値が0.5くらいの場合であること
    • PPV・NPVは、プレテスト・プロバビリティが低ければ、n/Nを大きくする方がよいし、高ければn/Nを小さくする方がよい
# 複数項目式診断の感度・特異度・PPV・NPV

# 総項目数
N<-11
# 診断基準必要項目数
n<-4

# 病気ありのときの個別項目の陽性率
p<-0.7
# 病気なしのときの個別項目の陽性率
q<-0.2

# 病気ありなしの事前確率
P<-0.5
Q<-1-P

M<-matrix(0,2,2)
for(i in n:N){
	M[1,1]<-M[1,1]+choose(N,i)*p^i*(1-p)^(N-i)
	M[2,1]<-M[2,1]+choose(N,i)*q^i*(1-q)^(N-i)
}
M[1,2]<-1-M[1,1]
M[2,2]<-1-M[2,1]

M

make2x2ClinicalTest<-function(N,n,p,q){
	M<-matrix(0,2,2)
	for(i in n:N){
		M[1,1]<-M[1,1]+choose(N,i)*p^i*(1-p)^(N-i)
		M[2,1]<-M[2,1]+choose(N,i)*q^i*(1-q)^(N-i)
	}
	M[1,2]<-1-M[1,1]
	M[2,2]<-1-M[2,1]

	M
}

sensSpecPPVNPV<-function(t,P=0.5,Q=0.5){
	sens<-t[1,1]/(t[1,1]+t[1,2])
	spec<-t[2,2]/(t[2,1]+t[2,2])
	PPV<-t[1,1]*P/(t[1,1]*P+t[2,1]*Q)
	NPV<-t[2,2]*Q/(t[1,2]*P+t[2,2]*Q)
	list(sens=sens,spec=spec,PPV=PPV,NPV=NPV)
}

t<-make2x2ClinicalTest(N,n,p,q)

sensSpecPPVNPV(t,P,Q)

p<-0.6
q<-0.3
N<-11
ns<-1:N
P<-0.2
Q<-1-P

out<-matrix(0,length(ns),4)

for(i in 1:length(ns)){
	t<-make2x2ClinicalTest(N,ns[i],p,q)
	print(t)
	out[i,]<-(unlist(sensSpecPPVNPV(t,P,Q)))
}

par(mfcol=c(2,2))
matplot(out,type="l")
plot(out[,1],out[,2],xlab="Sensitivity",ylab="Specificity",type="b")
plot(out[,3],out[,4],xlab="PPV",ylab="NPV",type="b")

p<-0.6
q<-0.3
N<-11
n<-5
Ps<-seq(from=0.1,to=0.9,by=0.1)
Qs<-1-Ps

out<-matrix(0,length(Ps),4)

for(i in 1:length(Ps)){
	t<-make2x2ClinicalTest(N,n,p,q)
	print(t)
	out[i,]<-(unlist(sensSpecPPVNPV(t,Ps[i],Qs[i])))
}

par(mfcol=c(2,2))
matplot(out,type="l")
plot(out[,1],out[,2],xlab="Sensitivity",ylab="Specificity",type="b")
plot(out[,3],out[,4],xlab="PPV",ylab="NPV",type="b")