All possible casesを対象としたバリデーション研究における感度とPPV -2
概要
前回の記事の続き。
「日本における傷病名を中心とするレセプト情報から得られる指標のバリデーションに関するタス クフォース」報告書のp20にあるall possible casesの場合の表を正しい表として、本当に感度にバイアスが生じ、PPVにバイアスが生じないのかをシミュレーションで検証する。
データ生成の準備
当該表の値を真の値としてデータを生成するための準備を行う。PPVはWang et al. 2021に記載の式を用いる。
# Set-up
N <- 2000 # total size
true_sens <- 180/240 # true sensitivity
true_fpr <- 1 - ((2000-240)-20)/(2000-240) # true false positive rate
p_true <- 240/N
## calculate ppv from the parameters above
calc_ppv <-
function(.prevalence, .sens, .spe) {
ppv <- .prevalence*.sens/(.prevalence*.sens + (1-.prevalence)*(1-.spe))
return(ppv)
} # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7949469/
true_ppv <- calc_ppv(p_true, true_sens, 1-true_fpr)
> c(true_sensitivity = true_sens, true_PPV = true_ppv, true_specificity = 1 - true_fpr)
true_sensitivity true_PPV true_specificity
0.7500000 0.9000000 0.9886364
データ生成
もっとかっこいい方法があるかもしれない(ChatGPTに聞くべきだったかも)が、上で設定した値を元に2by2テーブルを生成する。
# Generate sample data
set.seed(1)
D_true <- rbinom(N, 1, prob = p_true)
D_def <- rep(0,N)
D_def[which(D_true==1)] <- rbinom(sum(D_true==1), 1, prob = true_sens)
D_def[which(D_true==0)] <- rbinom(sum(D_true==0), 1, prob = true_fpr)
D_trueがD_fefが
> (twobytwo <- table(def = D_def,true = D_true)[c(2,1),c(2,1)])
true
def 1 0
1 185 19
0 58 1738
上で設定した値から生成したデータにおける感度・PPV・特異度は下記の通りであった。こちらはランダムサンプリングして抽出した集団に対する数値なので、この元となった推定量は不偏で、サンプリング誤差が乗っているのみである。
# Estimates
sens_ppv_spec <-
function(mat) {
sens <- mat[1,1]/sum(mat[,1])
ppv <- mat[1,1]/sum(mat[1,])
spec <- mat[2,2]/sum(mat[,2])
return(c(sens = sens, ppv = ppv, spec = spec))
}
> (correct_estimates <- sens_ppv_spec(twobytwo))
sens ppv spec
0.7613169 0.9068627 0.9891861
シミュレーションの実行
上で生成したD_defとD_trueをベースに、q (Zは
create_2by2_apc <-
function(D_def,D_true,q) {
Z <- rep(0,N) # APC flag vector
Z[which(D_def==1)] <- rbinom(sum(D_def==1), size = 1, prob = q)
# Arbitrary sampling probability of APC from non-definition positive sample
Z[which(D_def==0)] <- rbinom(sum(D_def==0), size = 1, prob = min(q*runif(1),0.01))
apc_size <- sum(Z)
# Create APC sample
D_def_apc <- rep(0,apc_size)
D_true_apc <- rep(0,apc_size)
D_def_apc <- D_def[Z==1]
D_true_apc <- D_true[Z==1]
# APC
twobytwo_apc <- table(def = D_def_apc,true = D_true_apc)[c(2,1),c(2,1)]
return(list(tabdata=data.frame(D_def,D_true,Z),
twobytwo_apc)
)
}
上で作成した関数を引数Uとして指定して、感度・PPV・特異度の理論値からの差分を計算する関数を作成する（引数Uを使うことで関数自体を指定できるという事実を最近知って使ってみたが、特にUを変える必要はなかったので特に意味はなかった）。
get_results <-
function(D_def,D_true,q, U) {
twobytwo_apc <- U(D_def,D_true,q)[[2]]
# Estimates
sens_apc <- sens_ppv_spec(twobytwo_apc)[1]
ppv_apc <- sens_ppv_spec(twobytwo_apc)[2]
spec_apc <- sens_ppv_spec(twobytwo_apc)[3]
return(
c(sen_bias = sens_apc - true_sens,
ppv_bias = ppv_apc - true_ppv, spec = spec_apc - (1-true_fpr))
)
}
シミュレーションの結果
q=1を指定した結果、感度は0.25大きい方向にバイアスがかかっている一方、PPVはほぼバイアスがなかったと言えるだろう。特異度はそもそもall possible casesで推定する対象とはされていないが、非常に小さく出てしまいそうである。（特異度とNPVもall possible caseを用いた方法で推定できるが、apcのみを使っては正しく推定できない。詳しくは、バリデーション研究の報告書を参照されたい）
> q <- 1
> (sim_res <-
+ rowMeans(replicate(100,
+ get_results(D_def,D_true,q,create_2by2_apc)
+ ))
+ )
sen_bias.sens ppv_bias.ppv spec.spec
0.247265530 0.006862745 -0.529494731
試しに、一つ2by2テーブルを作成してall possible casesとして選択された集団とそうでない集団に層別して見てみよう。
sample_data <- create_2by2_apc(D_def,D_true,q)
mat <-
table(
def = sample_data[[1]]$D_def,
true = sample_data[[1]]$D_true,
Z = sample_data[[1]]$Z
)[c(2,1),c(2,1),]
下記から分かるように、Z=0の表では真陽性(
> mat
, , Z = 0
true
def 1 0
1 0 0
0 58 1720
, , Z = 1
true
def 1 0
1 185 19
0 0 18
このサンプルテーブルにおける、all possible casesとnon-all possible casesでの感度などを見てみる。予想通り、non-all possible casesでは感度は0で、そのためにPPVが不定になっている。そして、all possible casesでは感度が1と過大評価されている。
> (results_non_apc <- sens_ppv_spec(mat[,,1]))
sens ppv spec
0 NaN 1
> (results_apc <- sens_ppv_spec(mat[,,2]))
sens ppv spec
1.0000000 0.9068627 0.4864865
例が悪くて申し訳ないが、
> (p_empirical <- sum(sample_data[[1]]$Z*D_true)/sum((D_true)))
[1] 0.7613169
> results_apc[1]*p_empirical
sens
0.7613169
> correct_estimates[1]
sens
0.7613169
Discussion