統計3
Tomoki Motegi
2024-09-27
有意差検定
今回も引き続きirisデータセットを使って考えていく。まずは解析にデータとそれを格納したを呼び出す。
## Rに収載されているデータ"iris"を呼び出して"iris_data"に格納
## なお、以降このデータを使うため、一度Rstudioを閉じた場合は一番最初から順々に動かして、各データを取ると良い。
iris_data = iris
## 最初の6行のみを表示
head(iris_data)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa## setosa, versicolor, virginicaの各グループごとに分けたデータを作成
iris_sub_set = subset(x = iris_data, subset = iris_data$Species == "setosa")
iris_sub_ver = subset(x = iris_data, subset = iris_data$Species == "versicolor")
iris_sub_vir = subset(x = iris_data, subset = iris_data$Species == "virginica")では最初にやることは何か？それは仮説の設定である。この時点で仮説を決定していない場合、HARKing (Hypothesizing After the Results are Known)である。ではまず「setosaとversicolorのがくの長さは異なる」と仮説を立ててみる。
t検定（2群の比較）
では上記の仮説の検証を行う。帰無仮説は「両群に差はない」であり、この帰無仮説を棄却するための有意水準を5%未満（P<0.05）とする。現時点での基礎情報が無いため、setosaとversicolorのがくの長さはどちらが大きいあるいは小さいということはわからない。このため両側検定を行っていく。
## setosaとversicolorのがくの長さを入れてt検定を行う
## t.test(x, y)でデータを入力すると実行される
## 内部で気にはウェルチのt検定を行っている(var.equal = FALSE)
## デフォルトでpaired = FALSE, alternative = "two.side"となっているが、わかりやすいように明示している
t.test(iris_sub_set$Sepal.Length, iris_sub_ver$Sepal.Length,
paired = FALSE, alternative = "two.side", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: iris_sub_set$Sepal.Length and iris_sub_ver$Sepal.Length
## t = -10.521, df = 86.538, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.1057074 -0.7542926
## sample estimates:
## mean of x mean of y
## 5.006 5.936実行すると2群のウェルチt検定が実行され、p-value < 2.2e-16となるため帰無仮説は棄却され、2群は有意に異なることが示された。しかしこの結果からはどちらが有意に大きいかはわからない。そこで再現の平均を見るとxよりyのほうが大きい。今回xにsetosa、yにversicolorを入れているため、「versicolorのがくの長さはsetosaに比べて有意に大きい」ことが分かった。
t検定（多群比較）
では仮説が「setosaとversicolorとvirginicaの三品種では全てがくの長さが異なる」と仮説を立ててみる。この場合の帰無仮説は「全群において差はない」だが、2群の検定で行う場合、品種の組み合わせが3通りのため有意水準をP<0.05/3とBonferroniの補正する必要がある。ではこれで行ってみる。
## setosaとversicolorのがくの長さを入れてt検定（上と同じ）
t.test(iris_sub_set$Sepal.Length, iris_sub_ver$Sepal.Length,
paired = FALSE, alternative = "two.side", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: iris_sub_set$Sepal.Length and iris_sub_ver$Sepal.Length
## t = -10.521, df = 86.538, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.1057074 -0.7542926
## sample estimates:
## mean of x mean of y
## 5.006 5.936## setosaとvirginicaのがくの長さを入れてt検定
t.test(iris_sub_set$Sepal.Length, iris_sub_vir$Sepal.Length,
paired = FALSE, alternative = "two.side", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: iris_sub_set$Sepal.Length and iris_sub_vir$Sepal.Length
## t = -15.386, df = 76.516, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.78676 -1.37724
## sample estimates:
## mean of x mean of y
## 5.006 6.588## versicolorとvirginicaのがくの長さを入れてt検定
t.test(iris_sub_ver$Sepal.Length, iris_sub_vir$Sepal.Length,
paired = FALSE, alternative = "two.side", var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: iris_sub_ver$Sepal.Length and iris_sub_vir$Sepal.Length
## t = -5.6292, df = 94.025, p-value = 1.866e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.8819731 -0.4220269
## sample estimates:
## mean of x mean of y
## 5.936 6.588この結果は全てp-value < 0.016666……であるため多重検定を行ってもすべての群で帰無仮説を棄却できる。このため「アヤメの三品種（setosaとversicolorとvirginica）では全てがくの長さが有意に異なる」と主張できる。ではこれを「アヤメでは全ての品種でがくの長さが有意に異なる」と記載できるか？これは当然間違いである。なぜなら他のアヤメを調べていないため、この結論を導くことはできない。しかし「アヤメでは全ての品種でがくの長さが有意に異なる可能性がある」ということは言えるだろうか？これはこの説を支持できる何らかの情報があれば主張できる。（例えばアヤメは歴史的にがくの長さで品種を分けていたなど。）このように有意差検定から出てきたものは結果であり、ここに1段階までの解釈を加えたものをディスカッションと呼ぶ。基本的には2段階まで解釈を加えた場合、overspeculationとして判断される。
Mann–WhitneyのU検定
続いてMann–WhitneyのU検定を行ってみる。irisデータセットはパラメトリックなため、このデータを取る日に台風が来てしまい、10個体ずつしか測定できなかったというストーリーがあり、更に世界で誰もアヤメのがくや花びらのサイズを調べておらず各データがどのような分布に従うかわからないと仮定してみる。ではまずデータを作成してみる。
## 各データから上から10行ずつ取ってくる
iris10_sub_set = head(x = iris_sub_set, n = 10L)
iris10_sub_ver = head(x = iris_sub_ver, n = 10L)
iris10_sub_vir = head(x = iris_sub_vir, n = 10L)
## iris10_sub_setだけ中身を見てみる
iris10_sub_set## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosaでは今回は花びらの幅（Petal.Width)について考える。まずデータを視覚化してみる。
## 視覚化で使いやすいようにデータをまとめてみる
iris10 = rbind(iris10_sub_set, iris10_sub_ver, iris10_sub_vir)
## データをバイオリンプロットで表現
library(ggplot2)## Warning: パッケージ 'ggplot2' はバージョン 4.2.3 の R の下で造られましたbp_gg = ggplot(data = iris10,
aes(x = Species, y = Petal.Width, fill = Species)) +
geom_violin(trim = FALSE) +
labs(x = "", y = "Petal.Width") +
geom_jitter(size = 1.0, width = 0.1, alpha = 0.3) +
scale_fill_manual(values = c("red", "green", "blue"))+
theme_classic()
## グラフを表示
bp_gg 
データが歪んでいて、正規分布するようなデータに見えない。我々の仮説は「3品種で花びらの幅は異なる」である。このため各群を比較するためにU検定を実施して各群が異なっているかを調べてみる。多群比較のt検定と同じく、品種の組み合わせが3通りのため、有意水準をP<0.05/3とBonferroniの補正を行う。
## Mann–WhitneyのU検定はWilcoxonの順位和検定と同じものを見ている
## このためwilcox.test(x, y)で実施できる
## setosaとversicolorの花びらの幅を入れてU検定
wilcox.test(iris10_sub_set$Sepal.Length, iris10_sub_ver$Sepal.Length,
paired = FALSE)## Warning in wilcox.test.default(iris10_sub_set$Sepal.Length,
## iris10_sub_ver$Sepal.Length, : タイがあるため、正確な p
## 値を計算することができません##
## Wilcoxon rank sum test with continuity correction
##
## data: iris10_sub_set$Sepal.Length and iris10_sub_ver$Sepal.Length
## W = 6, p-value = 0.0009817
## alternative hypothesis: true location shift is not equal to 0## setosaとvirginicaの花びらの幅を入れてU検定
wilcox.test(iris10_sub_set$Sepal.Length, iris10_sub_vir$Sepal.Length,
paired = FALSE)## Warning in wilcox.test.default(iris10_sub_set$Sepal.Length,
## iris10_sub_vir$Sepal.Length, : タイがあるため、正確な p
## 値を計算することができません##
## Wilcoxon rank sum test with continuity correction
##
## data: iris10_sub_set$Sepal.Length and iris10_sub_vir$Sepal.Length
## W = 5, p-value = 0.0007442
## alternative hypothesis: true location shift is not equal to 0## versicolorとvirginicaの花びらの幅を入れてU検定
wilcox.test(iris10_sub_ver$Sepal.Length, iris10_sub_vir$Sepal.Length,
paired = FALSE)## Warning in wilcox.test.default(iris10_sub_ver$Sepal.Length,
## iris10_sub_vir$Sepal.Length, : タイがあるため、正確な p
## 値を計算することができません##
## Wilcoxon rank sum test with continuity correction
##
## data: iris10_sub_ver$Sepal.Length and iris10_sub_vir$Sepal.Length
## W = 32, p-value = 0.1849
## alternative hypothesis: true location shift is not equal to 0警告: タイがあるため、正確な p 値を計算することができませんと表示されるが、このエラーは同じ順位になってしまうものがある場合に表示される。U検定は漸化式により正確な
P値を求めているが同じ順位の値がある場合は漸化式は使えないため、正規分布による近似値を用いている。しかしそもそもn=50の場合は自動的に正規分布に切り替わってしまう。このため実質的にはこの警告は無視されることが多い。気になる場合はexactRankTestsパッケージのwilcox.exact関数を用いる。
## 入っていない人はパッケージを以下のコマンドでインストールして使う
## BiocManager::install(c("exactRankTests"), force = TRUE)
library(exactRankTests)## Warning: パッケージ 'exactRankTests' はバージョン 4.2.3 の R の下で造られました## Package 'exactRankTests' is no longer under development.
## Please consider using package 'coin' instead.## setosaとversicolorの花びらの幅を入れてU検定
wilcox.exact(iris10_sub_set$Sepal.Length, iris10_sub_ver$Sepal.Length,
paired = FALSE)##
## Exact Wilcoxon rank sum test
##
## data: iris10_sub_set$Sepal.Length and iris10_sub_ver$Sepal.Length
## W = 6, p-value = 0.0002815
## alternative hypothesis: true mu is not equal to 0## setosaとvirginicaの花びらの幅を入れてU検定
wilcox.exact(iris10_sub_set$Sepal.Length, iris10_sub_vir$Sepal.Length,
paired = FALSE)##
## Exact Wilcoxon rank sum test
##
## data: iris10_sub_set$Sepal.Length and iris10_sub_vir$Sepal.Length
## W = 5, p-value = 0.0002057
## alternative hypothesis: true mu is not equal to 0## versicolorとvirginicaの花びらの幅を入れてU検定
wilcox.exact(iris10_sub_ver$Sepal.Length, iris10_sub_vir$Sepal.Length,
paired = FALSE)##
## Exact Wilcoxon rank sum test
##
## data: iris10_sub_ver$Sepal.Length and iris10_sub_vir$Sepal.Length
## W = 32, p-value = 0.1828
## alternative hypothesis: true mu is not equal to 0やってみるとわかるが、どちらも大した差はない。今回の結果を見るとsetosaはversicolor・virginica両社と差があるが、versicolorとvirginicaはp-value
=
0.1828であり、帰無仮説を棄却できない。このため、「アヤメの三品種のうちsetosaはversicolorとvirginicaの両者と比べて花びらの幅が有意に異なるが、versicolorとvirginicaの花びらの幅は有意に異なるかはわからなかった」となる。
カイ二乗（χ2）検定とフィッシャーの正確確率検定
今度は独立性の検定を行っていく。irisデータセットは連続値のため、本書で述べたデータを用いる。今回の帰無仮説は「膵炎の有無とv-Lipに関係性はない」である。
## 膵炎の有無(Panc_Posi or Nega)とv-Lipの増加(v-Lip_High or Norm)の2×2の分割表を作る。
data_matrix = matrix(c(52, 12, 14, 38), nrow = 2, byrow = TRUE)
df = as.data.frame(data_matrix)
rownames(df) = c("v-Lip_High", "v-Lip_Norm")
colnames(df) = c("Panc_Posi", "Panc_Nega")
df## Panc_Posi Panc_Nega
## v-Lip_High 52 12
## v-Lip_Norm 14 38ではまずχ2二乗検定を行っていく。
chisq.test(df)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: df
## X-squared = 32.347, df = 1, p-value = 1.289e-08p-value = 1.289e-08となり、十分に小さい値を取っているため帰無仮説は棄却され、「膵炎の有無とv-Lipの上昇には関係性がある」という結論になる。では続いてフィッシャーの正確確率検定を行ってみる。
fisher.test(df)##
## Fisher's Exact Test for Count Data
##
## data: df
## p-value = 5.971e-09
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 4.521894 31.187470
## sample estimates:
## odds ratio
## 11.44447こちらもp-value = 5.971e-09であり、帰無仮説は棄却された。今回は期待度数が5未満のものは無いため、当然フィッシャーの正確確率検定を無理に行う理由はない。なお、感度・特異度・陽性尤度比を計算すると下記のようになる。
## 感度
Sensitivity = df$Panc_Posi[1]/(df$Panc_Posi[1] + df$Panc_Posi[2])
Sensitivity## [1] 0.7878788## 特異度
Specificity = df$Panc_Nega[2]/(df$Panc_Nega[1] + df$Panc_Nega[2])
Specificity## [1] 0.76## 陽性尤度比
LR_posi = Sensitivity/(1-Specificity)
LR_posi## [1] 3.282828生存時間分析
カプランマイヤー曲線の引き方
この項ではカプランマイヤー曲線の引き方とその有意差検定のやり方を学んでいく。理論は本書を確認して貰いたい。まず今回はMASSパッケージ内にある42人の白血病患者に対して6-メルカプトプリンを投与したgehanデータ（Remission Times of Leukaemia Patients）を用いる。
## MASSパッケージの呼び出し
library(MASS)## Warning: パッケージ 'MASS' はバージョン 4.2.3 の R の下で造られました## gehanデータをdata_gehanに格納
data_gehan = gehan
## 中身の確認
head(data_gehan)## pair time cens treat
## 1 1 1 1 control
## 2 1 10 1 6-MP
## 3 2 22 1 control
## 4 2 7 1 6-MP
## 5 3 3 1 control
## 6 3 32 0 6-MPこの表は以下のように表現されている。
pair：投薬 (6-メルカプトプリン) 有無の比較のためのペア番号
time：生存時間（週）
cens：打ち切りか否か（0が打ち切り、1が再発）
treat：投薬したかどうかを表す (control = 無，6-MP = 有)
これを用いてまずはカプランマイヤー曲線を引いてみる。
## 必要なパッケージの呼び出し
library(survival)
library(survminer)## Warning: パッケージ 'survminer' はバージョン 4.2.3 の R の下で造られました## 要求されたパッケージ ggpubr をロード中です## Warning: パッケージ 'ggpubr' はバージョン 4.2.3 の R の下で造られました##
## 次のパッケージを付け加えます: 'survminer'## 以下のオブジェクトは 'package:survival' からマスクされています:
##
## myeloma## 生存時間情報を格納
## y = x の関係を生存時間（time）~ 治療（treat）で表現されている
KP_gehan = survfit(Surv(time, cens) ~ treat, data = data_gehan)
ggsurvplot(KP_gehan, ylab = "Probability of Remission", risk.table = FALSE) 
このように打ち切りも記載されたカプランマイヤー曲線を描画することができる。より詳細にカプランマイヤー曲線を引くために作ったデータの中身を見ていく。
KP_gehan## Call: survfit(formula = Surv(time, cens) ~ treat, data = data_gehan)
##
## n events median 0.95LCL 0.95UCL
## treat=6-MP 21 9 23 16 NA
## treat=control 21 21 8 4 12これを見ると寛解期間中央値（median）はは6-MP群で23週、control群で8週である。0.95LCLは寛解期間の95%信頼区間下限値であり、0.95UCLは寛解期間の95%信頼区間の上限値である。信頼区間も併せて表現すると以下のようなカプランマイヤー曲線が引ける。
ggsurvplot(KP_gehan, ylab = "Probability of Remission", risk.table = FALSE,
conf.int = TRUE, pval = FALSE) 
例数や打ち切り、瞬間死亡率の違いによって一部で信頼区間が重なっている場所もあるが、おおむね二群は離れていることがわかる。また比例ハザード性は症例の半分がいなくなると、評価できなくなるため、この情報も加えておく。
ggsurvplot(KP_gehan, ylab = "Probability of Remission", risk.table = TRUE,
conf.int = TRUE, pval = FALSE) 
ログランク検定とCOX比例ハザードモデル
ではログランク検定を行っていく。検定自体はsurvivalパッケージ内にあるsurvdiff関数を用いると2群間のイベント時間の違いを計算できる。今回の帰無仮説は「2群の寛解期間に差はない」である。
survdiff(Surv(time, cens) ~ treat, data = data_gehan)## Call:
## survdiff(formula = Surv(time, cens) ~ treat, data = data_gehan)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## treat=6-MP 21 9 19.3 5.46 16.8
## treat=control 21 21 10.7 9.77 16.8
##
## Chisq= 16.8 on 1 degrees of freedom, p= 4e-05この結果P = 4e-05となり十分に小さいため、帰無仮説は棄却され「白血病患者における6-MP治療は対照群と比べ寛解期間を有意に延長する」ことを示せた。このP値を載せたグラフは通常のRにおける画像描画用のplot関数で行うと煩雑になるが、survminerパッケージのggsurvplot関数では非常に簡単にできる。
ggsurvplot(KP_gehan, ylab = "Remission Time", risk.table = TRUE,
conf.int = TRUE, pval = TRUE) 
またコックス比例ハザードモデルを用いることもできる。単変量で実施した場合はコックス・マンテル検定を行っていること同じになる。実施する場合はsurvivalパッケージに入っているcoxph関数を用いる。
cox_gehan = coxph(Surv(time, cens) ~ treat, data = data_gehan)
cox_gehan## Call:
## coxph(formula = Surv(time, cens) ~ treat, data = data_gehan)
##
## coef exp(coef) se(coef) z p
## treatcontrol 1.5721 4.8169 0.4124 3.812 0.000138
##
## Likelihood ratio test=16.35 on 1 df, p=5.261e-05
## n= 42, number of events= 30pを見ると0.000138で有意であることは変わりないが、コックス比例ハザードモデルはハザード比を算出できる。ハザード比を見るために中身の詳細を見ると以下のようになる。
summary(cox_gehan)## Call:
## coxph(formula = Surv(time, cens) ~ treat, data = data_gehan)
##
## n= 42, number of events= 30
##
## coef exp(coef) se(coef) z Pr(>|z|)
## treatcontrol 1.5721 4.8169 0.4124 3.812 0.000138 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## treatcontrol 4.817 0.2076 2.147 10.81
##
## Concordance= 0.69 (se = 0.041 )
## Likelihood ratio test= 16.35 on 1 df, p=5e-05
## Wald test = 14.53 on 1 df, p=1e-04
## Score (logrank) test = 17.25 on 1 df, p=3e-05中段にハザード比とその信頼区間（lower .95とupper
.95）が表示されるようになり、治療の効果を評価することができるようになる。今回であれば治療群と対照群と比較した場合のハザード（イベントが発生するリスク）が約4.82倍になることを意味する。つまり、治療群は寛解期間に有意な効果を持っており、対照群に比べて治療群のリスクが大幅に低下していることを示唆している。さらにハザード比の95%信頼区間は
[2.147, 10.81]
のため、治療によるハザード比が2.147倍から10.81倍の範囲であると推定されることを示している。今回範囲には1.00を下回っていないため、治療群と対照群の間に統計的に有意な差があると結論付けられる。