ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2014-08-20 判定条件の検討

ryamada222014-08-20

[][][][][][][]判定条件の検討

  • スクリーニングという方法がある
  • 「完璧」はありえない
  • 感度を立てれば特異度を失う、そしてその逆も
  • 量的検査値でこれを考えるときにはrocカーブが使える
  • スクリーニングのときには、スクリーニング方法・条件をパラメタ化して、それに応じて感度特異度が出ることがあり、その場合、パラメタの大小について、個々のサンプルのIN/OUTが順序対応するとは限らないので、いわゆるrocカーブ関数をそのまま使うのには適さない
  • そんな話

f:id:ryamada22:20140820091600j:image

library(pROC)
N.true.positive <- 500
N.true.negative <- 4500
N.all <- N.true.positive + N.true.negative
truth <- c(rep(1,N.true.positive),rep(0,N.true.negative))

# 各サンプルについて、値が得られて、その大小で判定する場合
V <- rep(0,N.all)
for(i in 1:N.all){
	V[i] <- rpois(1,20+truth[i]*5)
}
par(mfrow=c(2,3))
boxplot(V~truth)

roc.out <- roc(truth,V)
plot(roc.out)

# 各サンプルについて、条件パラメタごとに、スクリーニングIN/OUTが決まる場合
param <- 0:6
Vs <- matrix(0,N.all,length(param))
for(i in 1:N.all){
	for(j in 1:length(param)){
		if(truth[i]==1){
			tmp.prob <- 0.8^(j)*0.8
			Vs[i,j] <- sample(0:1,1,prob=c(1-tmp.prob,tmp.prob))
		}else{
			tmp.prob <- 0.6^(j)*0.4
			Vs[i,j] <- sample(0:1,1,prob=c(1-tmp.prob,tmp.prob))
		}
	}
}

Sens <- Spec <- PPV <- NPV <- rep(0,length(param))
for(i in 1:length(param)){
	Sens[i] <- length(which(Vs[1:N.true.positive,i]==1))/N.true.positive
	Spec[i] <- length(which(Vs[(N.true.positive+1):N.all,i]==0))/N.true.negative
	PPV[i] <- length(which(Vs[1:N.true.positive,i]==1))/length(which(Vs[,i]==1))
	NPV[i] <- length(which(Vs[(N.true.positive+1):N.all,i]==0))/length(which(Vs[,i]==0))
}
plot(1-Spec,Sens,xlim=c(0,1),ylim=c(0,1),type="b")
text(1-Spec+0.05,Sens,paste("",param))

plot(1-NPV,PPV,xlim=c(0,1),ylim=c(0,1),type="b")
text(1-NPV+0.05,PPV,paste("",param))
matplot(cbind(Sens,Spec,PPV,NPV),type="l")
legend(5,0.5,c("Sens","Spec","PPV","NPV"),col=1:4,lty=1:4)

2012-10-13 あれか、これか、と感度・特異度

ryamada222012-10-13

[][][][]あれか、これか、と感度・特異度

  • あれか、これか、型の検査や情報というものを考える
    • 「吐物の色が赤か茶色か」
    • 「アレルがAかaか」
    • 「光学異性体のD体かL体か」など
  • その前に、検査学の復習をしておこう
    • ある疾患の検査前事前確率がpであり、ある検査があって、その感度がu、特異度がv であるとき
     | test(+)          | test(-)        |
---------------------------------------------
D(+) | pu               | p(1-u)         | p
D(-) | (1-p)(1-v)       | (1-p)v         | 1-p
---------------------------------------------
     | 1-(p+v) + p(u+v) | (p+v) - p(u+v) | 1
    • 検査陽性(test(+))であったなら、疾患である事後確率は pu/(1-(p+v) + p(u+v)) になり、疾患でない事後確率は (1-p)(1-v)/ (1-(p+v) + p(u+v)) になる
    • 検査陽性(test(-))であったなら、疾患である事後確率は p(1-u)/( (p+v) -p(u+v) ) になり、疾患でない事後確率は (1-p)v/( (p+v) - p(u+v) )になる
  • あれか、これか、の検査
    • 光学異性体のD体かL体かの検査のような場合を考える
    • D体らしさの検査をして、「らしかったら」D体と判断して「らしくなかったら」L体と判断するというやり方もあるだろう。この場合は、上の病気の検査と同じ
    • そうではなくて「D体らしかったら」Dに相当するシグナルが返り、「L体らしかったら」Lにそうとうする シグナルが返る、というような検査も可能だろう。ここでDとLはごく対称的なものであるから、その検査特性もよく似ているとする
    • このとき、u=v
     | test(+)          | test(-)            |
-------------------------------------------------
D(+) | pu               | p(1-u)             | p
D(-) | (1-p)(1-u)       | (1-p)u             | 1-p
-------------------------------------------------
     | pu + (1-p)(1-u)  |1-(pu + (1-p)(1-u)) | 1
    • このようなとき、事前確率pが検査後、pu/(pu+(1-p)(1-u)) に変わるか、p(1-u)/(1-(pu+(1-p)(1-u)))に変わる
ps <- seq(from=0,to=1,length=51)
ps <- ps[-1]
ps <- ps[-length(ps)]

us <- ps

pus <- as.matrix(expand.grid(ps,us))

p <- pus[,1]
u <- pus[,2]
post.1 <- p*u/(p*u+(1-p)*(1-u))
post.2 <- p*(1-u)/(1-(p*u+(1-p)*(1-u)))

fold.1 <- post.1/p
fold.2 <- post.2/p

image(matrix(fold.1,length(ps),length(us)))

library(rgl)
plot3d(cbind(pus,log(fold.1)))
open3d()
plot3d(cbind(pus,log(fold.2)))
  • uが高くないと採用しないとすると
  • 確率を高める方は、たかだか、1/p 倍にしかならないので、1 〜 1/p倍を変化するだけだが
  • 確率を低める方は、1 倍から、確率0へと指数関数的に確率を下げるので、
  • uを採用する基準を厳しくすると、それだけ確率を小さくする機会を取り除いてしまうことになる

f:id:ryamada22:20121013115625j:image

fold.change <- function(p,u){
	return(list(positive.fold = u/(p*u+(1-p)*(1-u)),negative.fold = (1-u)/(1-(p*u+(1-p)*(1-u))),inv.negative.fold = 1/((1-u)/(1-(p*u+(1-p)*(1-u)))) ))
}

fold.change(0.01,0.95)

# genotype freq = 0.01
p <- 0.01
# cut.off
u <- 1-(0.1^(1:10))
fold.outs <- fold.change(p,u)


par(mfcol=c(1,2))
plot(log10(1-u),log10(fold.outs$positive.fold))
plot(log10(1-u),log10(fold.outs$inv.negative.fold))
par(mfcol=c(1,1))

2012-03-23 疾患が一様でないときの感度特異度PPVNPV

ryamada222012-03-23

[][][][][]疾患が一様でないとき

  • 関節リウマチを例に取ろう
  • 以下の議論は「RF(+)リウマチとRF(-)リウマチ」という2種類のリウマチがあるので、厳密には、「診断確定のための検査」とは議論が異なるのだが、そのことはまた別の機会に書くとして。(気になる場合には、何かほかの検査なり情報なりに置き換えて考えるのがよいと思います)
  • 関節リウマチに比較的特異的な検査にリウマトイド因子(RF)がある
  • このRFをあるコミュニティの個人の全員に実施したとする
  • このコミュニティには、関節リウマチの患者が含まれ、その割合は有病率(a+b)
  • このコミュニティの関節リウマチ患者のうち、RF陽性の者(a)は、RF(+)リウマチ患者と呼ばれ、RF(-)の者(b)はRF(-)リウマチ患者と呼ばれる。2種類のリウマチは、成因・治療反応性・予後などについて、違いがあるものと考えられている
  • 他方、関節リウマチでなくとも、ある一定の割合のものはRF(+)となる(c)
疾患の有無RF(+)RF(-)
+aba+b
-cdc+d
a+cb+da+b+c+d=1
  • このテーブルを元にすると、有病率はa+b、感度は¥frac{a}{a+b}、特異度は¥frac{d}{c+d}
  • このような状況でRFを用いて、診断の助けにしようとするならば、リウマチであるという事前確率は有病率そのものであって、a+b、RF(+)であれば、PPVは¥frac{a}{a+c}、NPVは¥frac{d}{b+d}
  • 今、あるリウマチ専門の医療機関で同様の集計を取ったとする
疾患の有無RF(+)RF(-)
+a'b'a'+b'
-c'd'c'+d'
a'+c'b'+d'a'+b'+c'+d'=1
  • このような表から、この医療機関における、RF検査の感度・特異度・PPV・NPVは算出できるだろう
  • よくある感度・特異度・PPV・NPVの説明では、感度・特異度は適用対象集団の構成によらずに一定で、有病率(事前確率)だけが変化するように考える
  • しかしながら、このクリニックでの表は、このクリニックでの集計をすれば得られるわけだから、表の自由度は3だろう(a,b,c,dの4つの変数が、a+b+c+d=1という1つの制約を持っているから、4-1=3で、3)
  • 他方、コミュニティ全体の方も自由度は3(有病率と感度と特異度を決めたら、それで表全体もPPVもNPVも決まる)
  • このように、自由度3同士の2つの表で、ちょうどうまい具合に、感度・特異度が揃うとは限らない
  • では、感度・特異度が変化する(もちろんPPV・NPVも変化する)のはどうして?
  • それはこのクリニックの受診集団が、母集団とずれているから
  • そのずれの理由は
    • リウマチの疑いが高い人が訪れる
    • 1次医療機関で診断がつかない人が訪れる
    • 検診で引っかかって、念のため、としてやってくる
  • これらの要因は、事前確率(a'+b')を高めたり低めたりする
  • それ以外に、RF(+)リウマチとRF(-)リウマチの構成比も変化する
  • では、医療機関で、「コミュニティ」での感度・特異度(とそこから決まるPPV・NPV)を使うことの、精度とか是非とかは、どうやって決まるのだろうか?
  • また、年齢・性別・問診・身体所見などによって、リウマチである事前確率が上下し、その上で検査結果を適用しPPV・NPVを考えることは常道だが、そもそも、「事前情報」でリウマチらしい人と、らしくない人とを別々に集計してテーブルを作れば、そこから決まる、感度・特異度・PPV・NPVも変化してくるが、それはどうする?
  • 以下は、こんなことを考えるときのRの落書き
u<-0.9
v<-0.8
x<-seq(from=0,to=1,length=20)
p<-u*x/(u*x+(1-v)*(1-x))
plot(x,p,type="l")
q<-v*(1-x)/((1-u)*x+v*(1-x))
plot(x,q)
plot(x,q,type="l")

matplot(t(cbind(1-q,x,p)),type="l")

N<-2

P<-c(0.6,0.4)

Table<-list()
Table[[1]]<-matrix(c(0.2,0.1,0.3,0.4),byrow=TRUE,2,2)
Table[[2]]<-matrix(c(0.6,0.2,0.1,0.1),byrow=TRUE,2,2)

Whole<-matrix(0,2,2)
for(i in 1:N){
	Whole<-Whole+P[i]*Table[[i]]
}
Whole
Table

mySSPN<-function(t){
	list(sens=t[1,1]/(t[1,1]+t[1,2]),
	spec=t[2,2]/(t[2,1]+t[2,2]),
	PPV=t[1,1]/(t[1,1]+t[2,1]),
	NPV=t[2,2]/(t[1,2]+t[2,2]))
}

mySSPN(Table[[1]])
mySSPN(Table[[2]])
mySSPN(Whole)

2012-01-28 尤度比検定とROCカーブ?

[][][][][][][][][]カットオフ

  • こちらの企画
  • 親子鑑定
    • 親子かそうじゃないか(親子ではなくて、『ただの関係』)
    • 尤度の比をとる
  • ルール・インする尤度比の基準がある
  • ルール・アウトする尤度比の基準がある
  • 尤度比検定
    • パラメタがあって、仮説がパラメタで表現されている
    • 仮説ごとに尤度が計算できる
    • 尤度の比の対数を取ったものと、カイ自乗統計量の間には、パラメタ数の差〜自由度 な関係がある
  • カットオフを決める
    • 感度、特異度、PPV、NPVが決まる
  • この3者の関係を上手に説明してみたい
    • 仮説はいくつある?
    • 仮説を分類する(家系図上の距離で分類)?
    • 標本に観察された情報によって、仮説空間が変わる?
    • 感度・特異度がある?
    • PPV・NPVがある?
    • 感度・特異度・PPV・NPVのすべてがある?すべてはない?
    • 尤度比だけが、「基準」?
  • 尤度分布がパラメタの作る「ユークリッドな空間」にあるわけではなくて、「グラフ的な空間」にあったり、非ユークリッド空間にあったりするってことはない?
  • DNAを用いて、「多型」以外に情報が取れない?
    • メチル化などのエピゲノム、テロメアなどの構造情報・IBDの広がり

2011-10-26 診断基準と言う診断テスト

ryamada222011-10-26

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

  • こちらから
  • 診断基準というものを使って、診断することがある
  • もっとも単純には、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")