感度・特異度の計算や比較をRで実施する

はじめに

今回は、以前書いた感度・特異度に関する記事のRでの実施例を紹介する。
mstour.hatenablog.com
例題として、Alonzo et al.(2002)[1]で引用されている前立腺がんの早期発見に用いる診断法を比較した研究結果(Smith et al.(1997)[2])RACIAL DIFFERENCES IN OPERATING CHARACTERISTICS OF PROSTATE CANCER SCREENING TESTS - ScienceDirectを用いることにする。なおこの研究の詳細な内容には立ち入らない。

感度と特異度の計算

ある集団にてPSA検査とDigital rectal examination(DRE)の2つの診断法を実施した結果は以下の通りであった（現論文の表を見るとなぜか合計人数の足し算が合っていないのだけれど、理由がよくわからないので、ここでは合計人数を合わせたものを説明用の例とする）。

前立腺がんあり

	DREで陽性	DREで陰性
PSAで陽性	179	264	443
PSAで陰性	137	0	137
	316	264	580

前立腺がんなし

	DREで陽性	DREで陰性
PSAで陽性	138	717	855
PSAで陰性	976	0	976
	1114	717	1831

感度とは、疾患がある人のうち、陽性と診断できた人の割合である。したがってPSAの感度を $Se_1$ 、DREの感度を $Se_2$ とすると

$\displaystyle Se_1 = \frac{443}{580} = 0.763, \\ Se_2 = \frac{316}{580} = 0.544$

である。
一方、特異度は、疾患がない人のうち陰性になった人の割合となる。PSAの特異度を $Sp_1$ 、DREの特異度を $Sp_2$ とすると

$\displaystyle Sp_1 = \frac{976}{1831} = 0.533, \\ Sp_2 = \frac{717}{1831} = 0.391$

である。
手計算でも簡単にできるが、今回の趣旨はRでの計算なので、以下のようにして実行しよう。

#----- データの作成 -----
disease = c(rep("Y", 580), rep("N", 1831))
PSA = c(rep("+", 443), rep("-", 137), rep("+", 855), rep("-", 976))
DRE = c(rep("+", 179), rep("-", 264), rep("+", 137),
        rep("+", 138), rep("-", 717), rep("+", 976))

all = data.frame(disease = disease, PSA = PSA, DRE = DRE)
pc_y = all[disease == "Y", ]
pc_n = all[disease == "N", ]

一人一人の検査結果を模したデータを作り、疾患ありとなしのデータセットに分けた。
集計表は例えば次のようにして作成できる。

#----- 集計表の作成 -----
addmargins(table(pc_y$PSA, pc_y$DRE))
addmargins(table(pc_n$PSA, pc_n$DRE))

以下の結果が出力される。

> addmargins(table(pc_y$PSA, pc_y$DRE))
     
        -   + Sum
  -     0 137 137
  +   264 179 443
  Sum 264 316 580
> addmargins(table(pc_n$PSA, pc_n$DRE))
     
         -    +  Sum
  -      0  976  976
  +    717  138  855
  Sum  717 1114 1831

陽性は+、陰性は-としている。順番が最初の表と逆になっているが、変更する方法がわからなかった。申し訳ない。
感度と特異度の計算は、tableで集計した結果を使っても、元のデータを使ってもすぐに算出できる。
元のデータを使うなら、以下のようにする。

#PSAの感度
> nrow(pc_y[pc_y$PSA == "+", ])/nrow(pc_y)
[1] 0.7637931
#DREの感度
> nrow(pc_y[pc_y$DRE == "+", ])/nrow(pc_y)
[1] 0.5448276
#PSAの特異度
> nrow(pc_n[pc_n$PSA == "-", ])/nrow(pc_n)
[1] 0.5330421
#DREの特異度
> nrow(pc_n[pc_n$DRE == "-", ])/nrow(pc_n)
[1] 0.3915893

2つの診断方法の比較

f:id:mstour:20201003205421j:plain
例えば、2つの診断方法の感度を比較する、上の例で言えばがんがあることを正しく診断できる性能に違いがあるかどうかを検討するためには、よく知られている分割表の解析方法を用いることになる。特異度の場合も同様である（以下、感度の場合を考える）。
ここでの例では、同じ人に対して2つの診断方法を用いた結果をまとめている。したがって、対応のあるデータがとられている。このような場合には、感度（がんがある人のうち、陽性になる割合）が2つの方法で異なるかどうか、McNemar検定と呼ばれる方法で評価することができる。
ここで、「前立腺がんあり」の $2 \times 2$ 分割表の各セルに該当する確率を次のように表そう：

	DREで陽性	DREで陰性
PSAで陽性	$\pi_{11}$	$\pi_{12}$	$\pi_{1 \cdot}$
PSAで陰性	$\pi_{21}$	$\pi_{22}$	$\pi_{2 \cdot}$
	$\pi_{\cdot 1}$	$\pi_{\cdot 2}$	1

この時、PSAとDREとで感度が同じであるということは、それぞれで陽性になる確率が等しい、つまり $\pi_{1 \cdot} = \pi_{\cdot 1}$ と表現できる。ここで

$\displaystyle \pi_{1 \cdot} = \pi_{11} + \pi_{12}, \\ \pi_{\cdot 1} = \pi_{11} + \pi_{21}$

なので、 $\pi_{12} = \pi_{21}$ と書き直せる。
したがって、検定によって帰無仮説

$\displaystyle H_0 : \pi_{12} = \pi_{21}$

を否定できれば、2つの診断法の感度は異なると主張することができる。
このタイプの代表的な検定がMcNemar検定である。長くなるので、詳細はまた改めて書きたいと思う。
Rで行うには、まず集計表を作っておいて、それに対してmcnemar.test関数を使えばよい。

pc_y_table = table(pc_y$PSA, pc_y$DRE)

mcnemar.test(pc_y_table)

次のような結果が出力される。

> mcnemar.test(pc_y_table)

	McNemar's Chi-squared test with continuity correction

data:  pc_y_table
McNemar's chi-squared = 39.591, df = 1, p-value = 3.131e-10

検定統計量（漸近的にカイ二乗分布に従う）、自由度（ $2 \times 2$ 表のときは1）、p値の3つしか出力されないようだ。原論文と同様の結果になった（現論文では $p < 0.0001$ と記載）。PSAの感度は0.763、DREの感度は0.544だったので、PSAのほうが感度が高いという主張ができるだろう。
ところで、検定で有意になることだけでなく、推定された差が意味のあるレベルの差なのかを考えることがとても大事である。この例では $0.763 - 0.544 = 0.219$ と実際的にも意味のあるほどの差がついていると思われるが（本当にそうなのかは、対象としている分野の専門知識が必要）、サンプルサイズを増やせばいくらでもp値は小さくできるので、統計的に有意かどうかだけで物事を判断することに大して意味がないことには気をつけておかないといけない（実験的研究では、意味のある差がありそうな場合のみ検定で有意になるようにサンプルサイズをあらかじめ設計できるが、そうはいかないことも多い）。

まとめ

今回は感度・特異度についてRで計算する例を紹介した。感度については、「疾患のある人たち」を全サンプルとみなして、よく使われる分割表の解析を行えばよい（特異度については「疾患のない人たち」）。

参考文献

[1] Alonzo, T.A. et al. (2002). "Sample size calculations for comparative studies of medical tests for detecting presence of disease". Statist.Med.21:835–852.
[2] Smith, D. et al. (1997). "Racial differences in operating characteristics of prostate cancer screening tests". Journal of Urology.158:1861–1866.