井出草平の研究ノート

トップ > R > prevalenceパッケージまわりの補足[R]

prevalenceパッケージまわりの補足[R]

R 計量

こちらのエントリの補足。

ides.hatenablog.com

感度・特異度がわかっていれば、スクリーニング調査から真の有病率の推定ができる、という手法である。
もちろん、統計的推計なので、診断の有無をしっかりと調べていくようなものを再現できる訳ではなく、いろいろと欠点は存在するが、その点はまた別の機会に。

神戸の新型コロナウィルス抗体検査からの推計

検査人数(神戸): 1000
抗体陽性者数(神戸): 33
感度(クラボウ): 0.7638
特異度(クラボウ): 1

library("prevalence")
Prev <- truePrev(33, 1000, SE = 0.7638, SP = 1, nchains = 4, 
                            burnin = 10000, update = 10000, verbose = FALSE)
summary(Prev)

以前に小数点3位以下が表示できない方法しか示していなかったがsummary()を使えば、普通に使えることがわかった。

$TP
                 mean     median       mode          sd          var       2.5%      97.5% samples
chain 1    0.04433895 0.04404722 0.04423811 0.007389166 5.459978e-05 0.03095667 0.05996641   10000
chain 2    0.04433545 0.04396456 0.04359424 0.007482319 5.598510e-05 0.03073782 0.06040244   10000
chain 3    0.04469223 0.04432558 0.04262175 0.007613095 5.795921e-05 0.03085897 0.06043127   10000
chain 4    0.04431771 0.04396270 0.04456388 0.007538102 5.682298e-05 0.03075770 0.06056437   10000
all chains 0.04442108 0.04407672 0.04416214 0.007507470 5.636211e-05 0.03081621 0.06034275   40000

$SE
             mean median   mode sd var   2.5%  97.5% samples
chain 1    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 2    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 3    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
chain 4    0.7638 0.7638 0.7638  0   0 0.7638 0.7638   10000
all chains 0.7638 0.7638 0.7638  0   0 0.7638 0.7638   40000

$SP
           mean median mode sd var 2.5% 97.5% samples
chain 1       1      1    1  0   0    1     1   10000
chain 2       1      1    1  0   0    1     1   10000
chain 3       1      1    1  0   0    1     1   10000
chain 4       1      1    1  0   0    1     1   10000
all chains    1      1    1  0   0    1     1   40000

　感度・特異度の信頼区間の計算

クラボウのパンフレット(https://www.kurabo.co.jp/bio/English/support/download.php?M=PL&CID=9#catalog111)には感度・特異度が信頼区間の計算をしていないのでやってみよう。

f:id:iDES:20210212150417p:plain

まず分割表を作成する。

dat.kurabo <- as.table(matrix(c(97,0,30,394), nrow = 2, byrow = TRUE))
colnames(dat.kurabo) <- c("PCR+","PCR-")
rownames(dat.kurabo) <- c("Test+","Test-")
dat.kurabo

クラボウの抗体検査キットのテストの分割表である。

      PCR+ PCR-
Test+   97    0
Test-   30  394

感度読緯度はepiRパッケージのepi.testsを使用する。

www.rdocumentation.org

library(epiR)
dat.kurabo <- epi.tests(dat.kurabo, conf.level = 0.95)
print(dat.kurabo); summary(dat.kurabo)

結果は以下。

Point estimates and 95 % CIs:
---------------------------------------------------------
Apparent prevalence                    0.19 (0.15, 0.22)
True prevalence                        0.24 (0.21, 0.28)
Sensitivity                            0.76 (0.68, 0.83)
Specificity                            1.00 (0.99, 1.00)
Positive predictive value              1.00 (0.96, 1.00)
Negative predictive value              0.93 (0.90, 0.95)
Positive likelihood ratio              Inf (NaN, Inf)
Negative likelihood ratio              0.24 (0.17, 0.32)
---------------------------------------------------------

95%信頼区間を考慮したtruePrev()の推定

Prev2 <- truePrev(33, 1000, SE = ~dunif(0.68, 0.83), SP = ~dunif(0.99, 1.00), nchains = 4, 
                            burnin = 10000, update = 10000, verbose = FALSE)
summary(Prev2)

少しまともな推計になっている気がする。

$TP
                 mean     median       mode          sd          var       2.5%      97.5% samples
chain 1    0.03855326 0.03814766 0.03755704 0.008713162 7.591920e-05 0.02252821 0.05702676   10000
chain 2    0.03878964 0.03831414 0.03807106 0.008852792 7.837193e-05 0.02273063 0.05734680   10000
chain 3    0.03872776 0.03827282 0.03691512 0.008820912 7.780849e-05 0.02310364 0.05757549   10000
chain 4    0.03874573 0.03819929 0.03651368 0.008811758 7.764709e-05 0.02292265 0.05766032   10000
all chains 0.03870410 0.03822552 0.03722087 0.008799941 7.743896e-05 0.02279400 0.05738077   40000

$SE
                mean    median      mode         sd         var      2.5%     97.5% samples
chain 1    0.7535388 0.7531481 0.6970268 0.04316914 0.001863575 0.6839307 0.8256182   10000
chain 2    0.7526508 0.7514917 0.7112665 0.04338263 0.001882053 0.6834678 0.8259837   10000
chain 3    0.7521326 0.7505534 0.6963293 0.04311100 0.001858558 0.6833736 0.8257220   10000
chain 4    0.7529574 0.7512809 0.6955498 0.04319945 0.001866192 0.6838645 0.8262247   10000
all chains 0.7528199 0.7516303 0.6949127 0.04321705 0.001867714 0.6836405 0.8258802   40000

$SP
                mean    median      mode          sd          var      2.5%     97.5% samples
chain 1    0.9949075 0.9948642 0.9913735 0.002900040 8.410230e-06 0.9902096 0.9997145   10000
chain 2    0.9950091 0.9949967 0.9986009 0.002892036 8.363873e-06 0.9902785 0.9997547   10000
chain 3    0.9949031 0.9948460 0.9913536 0.002881859 8.305113e-06 0.9902280 0.9997297   10000
chain 4    0.9949990 0.9949875 0.9911999 0.002906795 8.449454e-06 0.9902602 0.9997894   10000
all chains 0.9949547 0.9949262 0.9914553 0.002895512 8.383993e-06 0.9902446 0.9997498   40000

別の書式で

SEとSPを別々に指定し、分布をdistで指定することもできるようだ。

SE <- list(dist = "uniform", min = 0.68, max = 0.83)
SP <- list(dist = "uniform", min = 0.99, max = 1.00)
TP <- truePrev(x = 33, n = 1000, SE = SE, SP = SP, nchains = 4, 
               burnin = 10000, update = 10000, verbose = FALSE)
summary(TP)

プロットも下記のように記述して出力できる。

par(mfrow = c(2, 2))
plot(TP)

f:id:iDES:20210212150441p:plain

密度プロットとトレースプロットだけを表示させるには以下のコードでOKである。

densplot(Prev, exclude_fixed = TRUE)
traceplot(Prev, exclude_fixed = TRUE)