## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = TRUE, fig.width = 5, fig.height = 5)
options(width = 116, scipen = 10)
setwd("~/statistics/Rmedstats/")
## ROC曲線を描いてみます。
## Introductory Statistics with R 2nd ed. (Peter Dalgaard)のデータを使っています。
## 1版は和訳 Rによる医療統計学 (岡田 昌史) もあります。
## 使い方
## この記載をすべてコピーして手元のRのエディタに貼付ける。
## Windows版では左上 [ファイル] → [新しいスクリプト]でエディタ。
## Mac版では真っ白い紙の形のアイコンでエディタ。
## 実行する行、あるいは範囲を選択して
## WindowsであればCtrl＋R、MacであればCommand＋Returnで実行。
## インストール用の関数定義
## rownames(installed.packages()でインストールされているパッケージの一覧
## PKG %in% ... という表現は...の中にPKGが含まれればTRUEを返します。!をつけて評価を逆転。
inst <- function (PKG) {
if (!(PKG %in% rownames(installed.packages()))) {
install.packages(pkgs = PKG, dependencies = TRUE)
}
}
## Introductory Statistics with Rの本のデータ集がなければ勝手にインストールします。
inst("ISwR")
## library()はパッケージを読み込むコマンド。ISwRが上記の本のデータ集パッケージです。
library(ISwR)
## データフレーム作り。要は月経発来(menarche)を年齢(age)で予測するためのデータ。
## juulというデータの読み込み。どんなデータかは ?juul で読めます。
data(juul)
## menarcheという変数が0, 1になっていますがわかりにくいのでNo, Yesにかえます
juul$menarche <- factor(juul$menarche, labels = c("No","Yes"))
## subset()はデータフレームから条件でsubgroupを抜き取ります。
## ここではage<8かつage<20かつ complete.cases()をつかってmenarche欠損値でないもの。
juul.girl <- subset(juul, age > 8 & age < 20 & complete.cases(menarche))
## ROC解析のパッケージの一つpROCをインストール。
## 制作元Swiss Institute of Bioinformatics: http://web.expasy.org/pROC/
## 説明: http://www.biomedcentral.com/1471-2105/12/77
inst("pROC")
## pROCを読み込みます。
library(pROC)
## ROC解析を行います。
## まずROCオブジェクトをつくります。
## outcome ~ predictorの書式を使います。
roc.menarche.by.age <- roc(formula = menarche ~ age, data = juul.girl)
## そのまま実行すると結果が見られます。
roc.menarche.by.age
Call:
roc.formula(formula = menarche ~ age, data = juul.girl)
Data: age in 256 controls (menarche No) < 263 cases (menarche Yes).
Area under the curve: 0.977
## ci()で信頼区間をつけられます。
ci(roc.menarche.by.age, of = "auc", method = "delong") # DeLong
95% CI: 0.967-0.986 (DeLong)
ci(roc.menarche.by.age, of = "auc", method = "boot") # bootstrap
95% CI: 0.967-0.985 (2000 stratified bootstrap replicates)
## of = "thresholds" をつけると任意のカットオフでの感度特異度の95%信頼区間がだせます。
ci(roc.menarche.by.age, of = "thresholds")
95% CI (2000 stratified bootstrap replicates):
thresholds sp.low sp.median sp.high se.low se.median se.high
11.21 0.566 0.625 0.684 1.000 1.000 1.000
11.54 0.609 0.668 0.723 0.989 0.996 1.000
11.66 0.629 0.688 0.742 0.981 0.992 1.000
11.88 0.688 0.742 0.793 0.973 0.989 1.000
11.98 0.719 0.770 0.820 0.970 0.985 0.996
12.14 0.730 0.785 0.832 0.958 0.977 0.992
12.18 0.734 0.789 0.836 0.951 0.973 0.992
12.23 0.746 0.797 0.844 0.947 0.970 0.989
12.32 0.754 0.805 0.852 0.943 0.966 0.989
12.41 0.766 0.816 0.863 0.939 0.962 0.985
12.54 0.789 0.836 0.879 0.932 0.958 0.981
12.57 0.797 0.844 0.887 0.928 0.954 0.977
12.72 0.812 0.859 0.902 0.920 0.947 0.973
12.80 0.816 0.863 0.902 0.909 0.939 0.966
12.82 0.828 0.871 0.910 0.905 0.935 0.966
12.84 0.832 0.875 0.914 0.894 0.924 0.954
12.88 0.844 0.883 0.918 0.886 0.920 0.951
12.93 0.844 0.887 0.922 0.878 0.913 0.947
13.09 0.863 0.902 0.938 0.871 0.905 0.939
13.11 0.867 0.906 0.938 0.867 0.901 0.935
13.29 0.875 0.910 0.942 0.852 0.890 0.924
13.38 0.895 0.930 0.957 0.840 0.882 0.920
13.48 0.902 0.934 0.961 0.833 0.875 0.913
13.57 0.906 0.938 0.965 0.825 0.867 0.905
13.68 0.922 0.949 0.973 0.817 0.859 0.901
13.75 0.926 0.953 0.977 0.802 0.848 0.890
13.82 0.941 0.965 0.984 0.798 0.844 0.886
13.89 0.961 0.980 0.996 0.795 0.840 0.882
14.04 0.973 0.988 1.000 0.779 0.825 0.871
14.35 0.980 0.992 1.000 0.749 0.795 0.844
14.71 0.988 0.996 1.000 0.688 0.741 0.795
14.95 1.000 1.000 1.000 0.646 0.703 0.760
## plot()するとグラフになります。
## 何も考えずに書くとpROCはX軸を一般的ではない(けど確かにわかりやすい)方法でかきます。
## legacy.axes = TRUE を与えるとX軸が 1 - specificityになります。
plot(roc.menarche.by.age)
Call:
roc.formula(formula = menarche ~ age, data = juul.girl)
Data: age in 256 controls (menarche No) < 263 cases (menarche Yes).
Area under the curve: 0.977
plot(roc.menarche.by.age, legacy.axes = TRUE)
Call:
roc.formula(formula = menarche ~ age, data = juul.girl)
Data: age in 256 controls (menarche No) < 263 cases (menarche Yes).
Area under the curve: 0.977
## 複数のROC曲線の比較を行ってみましょう。
## pROCの中にaSAHという113例の動脈瘤破裂のSAHのデータがあるので使用します。
## Intensive Care Med 2010;36:107-15 http://www.ncbi.nlm.nih.gov/pubmed/19760205
data(aSAH)
## データ構造をみます。outcomeがGood vs Poorのdichotomous outcomeです。
## gos6: Glasgow Outcome Score (GOS) at 6 months は高いほどよいアウトカム
## outcome: gos6が4,5をGoodと定義している。
## wfns: World Federation of Neurological Surgeons (WFNS) at admission なんだろう?
## s100b: S100beta at admission 血清マーカーのようです。brain injury-related biomarkersとpubmedに。
## ndka: Nucleoside DiPhosphate Kinase A at admission。これも同じような血清マーカー。
summary(aSAH)
gos6 outcome gender age wfns s100b ndka
1:28 Good:72 Male :42 Min. :18.0 1:39 Min. :0.030 Min. : 3.0
2: 0 Poor:41 Female:71 1st Qu.:42.0 2:32 1st Qu.:0.090 1st Qu.: 9.0
3:13 Median :51.0 3: 4 Median :0.140 Median : 12.2
4: 6 Mean :51.1 4:16 Mean :0.247 Mean : 19.7
5:66 3rd Qu.:61.0 5:22 3rd Qu.:0.330 3rd Qu.: 17.3
Max. :81.0 Max. :2.070 Max. :419.2
## responseがアウトカム、predictorが予測する変数です。3つの変数すべてで行います。
## formulaはoutcome ~ predictorで書きます。
roc.wfns <- roc(formula = outcome ~ wfns, data = aSAH)
roc.s100b <- roc(formula = outcome ~ s100b, data = aSAH)
roc.ndka <- roc(formula = outcome ~ ndka, data = aSAH)
## そのまま表示するとAUCがみられます。
roc.wfns
Call:
roc.formula(formula = outcome ~ wfns, data = aSAH)
Data: wfns in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.824
roc.s100b
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH)
Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.731
roc.ndka
Call:
roc.formula(formula = outcome ~ ndka, data = aSAH)
Data: ndka in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.612
## グラフ化しましょう
plot(roc.wfns, lty=1, legacy.axes = TRUE)
Call:
roc.formula(formula = outcome ~ wfns, data = aSAH)
Data: wfns in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.824
## s100bとndkaを書き込みます。add = TRUEでさっきのグラフに書き入れます。Line TYpeで点線にします。
plot(roc.s100b, lty = 2, add = TRUE)
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH)
Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.731
plot(roc.ndka, lty = 3, add = TRUE)
Call:
roc.formula(formula = outcome ~ ndka, data = aSAH)
Data: ndka in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.612
## 凡例はやや煩雑ですが、Box TYpe="None"で凡例の枠なし。Line WiDthで太さ、Line TYpeはさっきの1,2,3にあわせて
## legendでさっきの順番でwfns, s100b, ndkaと入れています。
legend(x = 0.7, y = 0.3, bty = "n", lwd = 2, lty = 1:3,
legend = c("wfns","s100b","ndka"))
## AUCを書き加えてみます。legendで線の太さを0にして線を消してAUCの数字だけ少しずらして入れています。
legend(x = 0.5, y = 0.3, bty = "n", lwd = 0,
legend =
c(
round(roc.wfns$auc, digits = 2),
round(roc.s100b$auc, digits = 2),
round(roc.ndka$auc, digits = 2)
)
)
## 2つのROCの比較を行います。書式は下記の二つ。
## roc.test(roc1=, roc2=, method="delong") 上記のようにrocオブジェクトを作ってあればそれを使用。
roc.test(roc1 = roc.wfns, roc2 = roc.s100b, method = "delong")
DeLong's test for two correlated ROC curves
data: roc.wfns and roc.s100b
Z = 2.209, p-value = 0.02718
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
0.8237 0.7314
##
## DeLong's test for two correlated ROC curves # correlatedtつまり同じ患者に二種類の
## # 検査をしたpairedの検定です。
## data: roc.wfns and roc.s100b
## Z = 2.209, p-value = 0.02718
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8236789 0.7313686
## roc.test(response =, predictor1 =, predictor2 =) オブジェクトがなければ変数を直に入れることもできる。
roc.test(response = aSAH$outcome,
predictor1 = aSAH$wfns,
predictor2 = aSAH$s100b,
method = "delong")
DeLong's test for two correlated ROC curves
data: aSAH$wfns and aSAH$s100b by aSAH$outcome (Good, Poor)
Z = 2.209, p-value = 0.02718
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
0.8237 0.7314
## Bootstrap法もサポート。だいたいにたようなp-valueになるかと思います。
roc.test(roc1 = roc.wfns, roc2 = roc.s100b, method = "bootstrap", boot.n = 2000)
Bootstrap test for two correlated ROC curves
data: roc.wfns and roc.s100b
D = 2.219, boot.n = 2000, boot.stratified = 1, p-value = 0.02647
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
0.8237 0.7314
## unpairedの検定、つまり同じ検査を二群の患者に行ってROCを作って比較することもできます。
## 男女で患者を分けて
aSAH.male <- aSAH[aSAH$gender == "Male",] # 男性群
aSAH.female <- aSAH[aSAH$gender == "Female",] # 女性群
roc.wfns.male <- roc(formula = outcome ~ wfns, data = aSAH.male)
roc.wfns.female <- roc(formula = outcome ~ wfns, data = aSAH.female)
roc.test(roc.wfns.male, roc.wfns.female)
DeLong's test for two ROC curves
data: roc.wfns.male and roc.wfns.female
D = 1.277, df = 106, p-value = 0.2043
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2
0.8761 0.7786
##
## DeLong's test for two ROC curves # correlatedがなくなっていてunpairedであることが分かります。
##
## data: roc.wfns.male and roc.wfns.female
## D = 1.2772, df = 106.014, p-value = 0.2043
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.8761364 0.7785714
if (FALSE) {
### 細かい数値化はDiagnosisMed::ROC()が得意です。グラフはいまいちですが。# CRANから削除されました。
## DiagnosisMedがなければインストール
inst("DiagnosisMed") ; library(DiagnosisMed)
Diag.Med.ROC <-
DiagnosisMed::ROC(gold = as.numeric(aSAH$outcome), # 0,1の数値を与える必要あり。
test = as.numeric(aSAH$wfns),
Plot = FALSE, Print = FALSE )
plot(Diag.Med.ROC) # グラフはみにくいです Diag.Med.ROC
## 結果を表示するといろいろ出ます。
## 感度、特異度、PPV, NPV、positive LR, negative LRなどがcut offごとに取り出せます。
Diag.Med.ROC$test.diag.table[,c("test.values", "TP", "FN", "FP", "TN", "Sensitivity","Specificity", "PPV", "NPV", "PLR", "NLR")]
}