並べ替え検定と正規近似
2017-07-11投稿, 2019.7.4更新
正確なｐ値 
特定の確率分布をもとに推定を行うのではなく、母集団の未知のパラメータやサンプリング誤差が入らないため計算上も正しいp値が得られる．
並べ替え検定
例）ｘ群とｙ群を比較します
x <- c(5, 9) #平均 = 7 y <- c(6, 12, 14, 16) #平均 = 12
xとｙは同じ母集団からのサンプリングと考えます．それぞれのグループへの割付の際にたまたま差が生じました．なおｘとｙは正規分布には従いません．
帰無仮説：xとｙに差はない
対立仮説：xよりｙが大きい
並べ替え検定の考え方
もしが正しいと考えるとき、このサンプリングの平均差がどの程度大きいのかを考えます．（ 5 , 6 , 9 , 12 , 14 , 16 ) がマークされている、6個の同質の玉が袋に入っていると考えます．
帰無仮説は、どのように取り出してもx（2個）とｙ（4個）が示す増加量が等しいということになります．
取り出し方は、
choose (6, 2)
の15通りあります．帰無仮説が正しいとすると、15通り全てが「xとｙは等しい」ということになります．
ｘの平均は7、ｙの平均は12なので、ｙ－ｘ＝5となります．つまり15通り中、差が平均の差5より大きくなる確率を正確なｐ値として考えます．
差が5以上になるのは2通りなので、片側検定のｐ値は2／15 =0.1333333 、両側検定のp値は 4/15 = 0.2666667 となります．
Rでは以下のような計算式になります
Package ‘coin’, March 8, 2019, Version 1.3-0, Date 2019-03-04より
https://cran.r-project.org/web/packages/coin/coin.pdf
#サンプルデータセット diffusion <- data.frame( pd = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46, 1.15, 0.88, 0.90, 0.74, 1.21), age = factor(rep(c("At term", "12-26 Weeks"), c(10, 5))) ) # ソート (sortlist <- order(diffusion$pd)) ( diffusion <- diffusion[sortlist,] ) diffusion$no <- c(1:15);diffusion #セット > diffusion pd age no 9 0.73 At term 1 14 0.74 12-26 Weeks 2 1 0.80 At term 3 2 0.83 At term 4 12 0.88 12-26 Weeks 5 13 0.90 12-26 Weeks 6 4 1.04 At term 7 11 1.15 12-26 Weeks 8 15 1.21 12-26 Weeks 9 6 1.38 At term 10 5 1.45 At term 11 10 1.46 At term 12 8 1.64 At term 13 3 1.89 At term 14 7 1.91 At term 15
Rのパッケージ"coin"を使って検定します
#install.packages("coin") library(coin) #正確なP値を求めます：並び替え検定、Exact Wilcoxon-Mann-Whitney test (wt <- wilcox_test(pd ~ age, data = diffusion, distribution = "exact", conf.int = TRUE)) Exact Wilcoxon-Mann-Whitney Test data: pd by age (12-26 Weeks, At term) Z = -1.2247, p-value = 0.2544 alternative hypothesis: true mu is not equal to 0 95 percent confidence interval: -0.76 0.15 sample estimates: difference in location -0.305
期待値や分散など
statistic(wt, type = "linear")# sum of the ranks for age = "12-26 Weeks" 12-26 Weeks 30 expectation(wt)# 平均順位の差の期待値 12-26 Weeks 40 covariance(wt)# 平均順位の差の分散 12-26 Weeks 12-26 Weeks 66.66667 pvalue(wt) [1] 0.2544123 confint(wt) 95 percent confidence interval: -0.76 0.15 sample estimates: difference in location -0.305
Mann-Whitney の U 検定, wilcoxonの順位和検定（Wilcoxon rank-sum test）
正規近似を行って検定する方法．なお，この正規近似は m，n が 7 より大きければかなり正確であることも示されている．
ノンパラメトリック検定
(wt2 <- wilcox_test(pd ~ age, data = diffusion, conf.int = TRUE)) Asymptotic Wilcoxon-Mann-Whitney Test data: pd by age (12-26 Weeks, At term) Z = -1.2247, p-value = 0.2207 alternative hypothesis: true mu is not equal to 0 95 percent confidence interval: -0.7599742 0.1499702 sample estimates: difference in location -0.3038399
コンピュータの発達により、並べ替え検定や正確確率検定を行うのは無理ではなくなったので、わざわざｐ値の近似値を求める従来のパラメトリック検定やノンパラメトリック検定よりも、直観的かつ、わかりやすい結果が得られるといえるだろう． 水本 篤(2010)より引用
参考
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p162-164
水本 篤：統計数理研究所共同研究リポート 238『言語コーパス分析における数理データの統計的処理手法の検討』(2010) pp. 1–14
Rで簡単（二元配置分散分析：データには対応なし、繰り返し数は5）
投稿2016.11.17 最終更新日2019.6.5
検定は帰無仮説が大切です！
商品のの主効果：帰無仮説＝カップ麺でもインスタント麺でも評価の母平均は等しい
スープの主効果：帰無仮説＝スープが違っても評価の母平均は等しい
商品とスープの主効果：帰無仮説＝商品とスープの組合せと評価は関係ない
商品 スープ 評価 カップ とんこつ 10 カップ とんこつ 11 カップ とんこつ 11 カップ とんこつ 9 カップ とんこつ 9 カップ 醤油 11 カップ 醤油 8 カップ 醤油 10 カップ 醤油 8 カップ 醤油 8 カップ みそ 7 カップ みそ 3 カップ みそ 5 カップ みそ 2 カップ みそ 3 インスタント とんこつ 11 インスタント とんこつ 11 インスタント とんこつ 10 インスタント とんこつ 9 インスタント とんこつ 10 インスタント 醤油 9 インスタント 醤油 7 インスタント 醤油 9 インスタント 醤油 8 インスタント 醤油 8 インスタント みそ 4 インスタント みそ 4 インスタント みそ 2 インスタント みそ 3 インスタント みそ 3
コピーしてRに読み込みます
x<-read.table("clipboard",header=T) summary(aov(x$評価~x$商品*x$スープ)) Df Sum Sq Mean Sq F value Pr(>F) x$商品 1 1.63 1.63 1.077 0.310 x$スープ 2 231.67 115.83 76.374 3.93e-11 *** x$商品:x$スープ 2 1.67 0.83 0.549 0.584 Residuals 24 36.40 1.52 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
スープによる主効果あり！
次の式でも同じ結果になります
summary(aov(x$評価~x$商品+x$スープ+x$商品:x$スープ))
ベクトルを使用した方法
評価<-c(x$評価) 商品<-factor(c(rep("カップ",15),rep("インスタント",15))) スープ<-factor(rep(c(rep("とんこつ",5),rep("醤油",5),rep("みそ",5)),2)) summary(aov(x$評価~x$商品*x$スープ))
検証するために図を描きましょう！
par(mfrow=c(1,2)) interaction.plot(商品,スープ,評価) interaction.plot(スープ,商品,評価) par(mfrow=c(1,1))
相関行列と共分散行列
相関行列
x1 <- c(93, 89, 115, 90, 75) x2 <- c(121, 136, 121, 161, 125) x3 <- c(101, 115, 118, 122, 106) #データフレームにすること dat <- data.frame(x1, x2, x3) cor(dat) x1 x2 x3 x1 1.0000000 -0.2220449 0.3981155 x2 -0.2220449 1.0000000 0.6767614 x3 0.3981155 0.6767614 1.0000000
散布図
pairs(dat)
共分散行列
これもデータフレームで var(dat) x1 x2 x3 x1 207.80 -54.15 49.80 x2 -54.15 286.20 99.35 x3 49.80 99.35 75.30
参考
install.packages("psych") library(psych) x <- rnorm(100,5,3) y <- rnorm(100,15,8) z <- rnorm(100,2,0.5) dat2 <- data.frame(x,y,z) pairs.panels(dat2)
対応のある一元配置分散分析（繰り返しのない二元配置分散分析）、多重比較
以下のサイトに移転いたしました
生存時間解析の基本
パッケージのインストール
install.packages("survival") library(survival)
survivalのデータセットheartを使用します
head(heart) start, stop, event: Entry and exit time and status for this interval of time age: age-48 years year: year of acceptance (in years after 1 Nov 1967) surgery: prior bypass surgery 1=yes#バイパス手術あり transplant: received transplant 1=yes#移植を受けている id: patient id
バイパス手術あり＆なしの中央値
pfit <- survfit(Surv(stop-start,event) ~ surgery, heart) Call: survfit(formula = Surv(stop - start, event) ~ surgery, data = heart) n events median 0.95LCL 0.95UCL surgery=0 143 66 102 64 280 surgery=1 29 9 897 322 NA
バイパス手術あり＆なしで生存曲線を描く、赤(あり)vs黒(なし)
plot(pfit, col = 1:2)
バイパス手術あり＆なしの有意差検定 log-rank test
rho = 0 （デフォルト）this is the log-rank or Mantel-Haenszel test
rho = 1 it is equivalent to the Peto & Peto modification of the Gehan-Wilcoxon test.
diff_surg <- survdiff(Surv(stop-start,event)~surgery , heart, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ surgery, data = heart, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V surgery=0 143 66 57.8 1.15 5.17 surgery=1 29 9 17.2 3.87 5.17 Chisq= 5.2 on 1 degrees of freedom, p= 0.02
p値のみを算出
pchisq(diff_surg <- survdiff(Surv(stop-start,event)~surgery , heart, rho = 0)$chisq, 1, lower.tail=FALSE ) [1] 0.02294185
surgeryを1にすれば、全体の曲線と95％信頼区間をプロットする
pfit2 <- survfit(Surv(stop-start,event) ~ 1, heart) plot(pfit2, col = c(1,4,4))
年齢で比較 48歳以上と48歳未満で検定
diff_age1 <- survdiff(Surv(stop-start,event)~ age<0 , heart, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ age < 0, data = heart, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V age < 0=FALSE 81 39 32.5 1.283 2.29 age < 0=TRUE 91 36 42.5 0.983 2.29 Chisq= 2.3 on 1 degrees of freedom, p= 0.1
年齢で比較 50歳以上と40以下
sub <- subset(heart, age>=2 | age <= -8) diff_age2 <- survdiff(Surv(stop-start,event)~ age>=2 , sub, rho = 0) Call: survdiff(formula = Surv(stop - start, event) ~ age >= 2, data = sub, rho = 0) N Observed Expected (O-E)^2/E (O-E)^2/V age >= 2=FALSE 33 10 18.2 3.67 6.82 age >= 2=TRUE 57 31 22.8 2.92 6.82 Chisq= 6.8 on 1 degrees of freedom, p= 0.009
カテゴリー変数がある場合の層別 カテゴリー変数名をＡとします。Ａには1,2,3,4が含まれており、1と3を比較します
x<-1; y<-3 cat <- subset(heart, Ａ==x | Ａ==y, data=heart ) diff_c <- survdiff(Surv(stop-start,event)~"カテゴリ－変数名" , cat)
COX比例ハザードモデル
cox_age_sur <- coxph(Surv(stop-start,event) ~ age + surgery, data=heart) summary(cox_age_sur) Call: coxph(formula = Surv(stop - start, event) ~ age + surgery, data = heart) n= 172, number of events= 75 coef exp(coef) se(coef) z Pr(>|z|) age 0.02941 1.02985 0.01384 2.125 0.0336 * surgery -0.82104 0.43997 0.35872 -2.289 0.0221 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.03 0.971 1.0023 1.0582 surgery 0.44 2.273 0.2178 0.8887 Concordance= 0.603 (se = 0.039 ) Rsquare= 0.062 (max possible= 0.975 ) Likelihood ratio test= 10.96 on 2 df, p=0.004 Wald test = 9.77 on 2 df, p=0.008 Score (logrank) test = 10.17 on 2 df, p=0.006
生存率の予測
fit_surv <- survfit(cox_age_sur) plot(fit_surv,col = c(2,3,3)) haz50 <- summary(fit_surv) Call: survfit(formula = cox_age_sur) time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 172 1 0.995 0.00537 0.9842 1.000 1.0 171 2 0.984 0.00926 0.9659 1.000 2.0 166 3 0.967 0.01314 0.9421 0.994 3.0 160 4 0.945 0.01703 0.9124 0.979 5.0 150 1 0.939 0.01792 0.9050 0.975 6.0 147 2 0.928 0.01962 0.8900 0.967 8.0 144 1 0.922 0.02042 0.8826 0.963 9.0 141 1 0.916 0.02121 0.8750 0.958 10.0 140 1 0.910 0.02197 0.8675 0.954
