理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

Rで簡単 多重ロジスティック回帰分析

サンプルはRのmtcarsを以下のように加工して、応答変数をY1, 説明変数をx2,x5,x6とします.
車に詳しくないので回帰分析に相応しいサンプルか分かりませんが、学習のため数値のみ使用させていただきます.
f:id:yoshida931:20180517164535p:plain:w300

クリップボードの取り込みから

df <- read.table("clipboard",header=T,row.names=1)
dat <- df[,c(8,2,5,6)]    #ここがポイント8列目を応答変数として、2列目・5列目・6列目を説明変数として設定
ans <- glm(dat$Y1~.,data=dat,family=binomial) 

ans <- glm(dat$Y1~dat$x2+dat$x5+dat$x6,family=binomial)   #もちろんこの式も同じ答えになります

これで終了です.

summary(ans)

Call:
glm(formula = dat$Y1 ~ dat$x2 + dat$x5 + dat$x6, family = binomial)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.43462  -0.15612  -0.06808   0.22595   1.35668  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   22.302     11.577   1.926   0.0541 .
dat$x2        -3.375      1.544  -2.186   0.0288 *
dat$x5        -2.189      1.874  -1.168   0.2428  
dat$x6         1.777      1.660   1.070   0.2845  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.86  on 31  degrees of freedom
Residual deviance: 13.92  on 28  degrees of freedom
AIC: 21.92

Number of Fisher Scoring iterations: 7

to-kei.net 上記サイトで勉強させていただきました.

ロジスティック回帰分析(説明変数が単一かつ連続の場合)

投稿日:2018.2.13 最終更新日:2018.5.17

ロジット関数とロジスティック関数

準備として関数の特徴を押さえておきます

ロジット関数 \ \  y=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

(標準)ロジスティック関数 \ \ ロジット関数の逆関数=P = \frac{1}{1+exp(-y)}

mathwords.net

サンプルirisより

品種"virginica=1"、"別の品種=0"という2値を目的変数、”Sepal.Length(がく片の長さ)(cm)”を説明変数とします. ロジスティック回帰で求めることは、”Sepal.Length(がく片の長さ)(cm)”から予測される、品種が"virginica"になる確率です.

class(iris$Species)  #属性を確認
is <- as.numeric(iris$Species)   #  iris$Speciesが要因になっているので数字に変換(as.integerでも同じ)
y <- gsub(1, 0, is)  # is の1を0に変換
y <- gsub(2, 0, y)  # is の2を0に変換
y <- gsub(3, 1, y)  # is の3を1に変換

x <- iris$Sepal.Length

#実際のデータは以下のようになっています
> y
  [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [24] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [47] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [70] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [93] "0" "0" "0" "0" "0" "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[116] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[139] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
> x
  [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4 5.1 5.7 5.1 5.4 5.1 4.6
 [24] 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8
 [47] 5.1 4.6 5.3 5.0 7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8 6.2
 [70] 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4 6.0 6.7 6.3 5.6 5.5 5.5 6.1
 [93] 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8
[116] 6.4 6.5 7.7 7.7 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7 6.3 6.4
[139] 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9

ちなみに、そのまま単回帰分析を行ってみます.

summary(lm(formula=y~x) )
plot(x,y)
abline(lm(y~x)) #回帰直線
f:id:yoshida931:20180208184829p:plain:w300

目的変数は0、1なのですが、上図ではかなりはみ出してしまいます.0より少なくても、1より大きくても、良い予測とは言えません.したがって、この回帰直線は適当なモデルとは言えません.そこで指数を用いたロジスティックモデルを使用します.

ロジスティック変換 logistic transformation

目的変数が2値変数(品種=virginica=1、品種=別の品種=0)の回帰分析をおこなうためには、2値変数に与えられる上限と下限を取り除けば解決します.ロジスティック変換は、0~1の区間制限を取り除くことで、統計的推測に柔軟性を持たせた方法といえます.

品種が"virginica(y=1)"になる確率をPとした場合のロジスティック(ロジット)変換 を考えてみます.

オッズ=\frac{品種=virginica=1 である確率}{品種=他の品種=0 である確率}=\frac{p}{1-p}

ロジット(logit)とは、0から1の値をとる確率p に対し、そのオッズ\frac{p}{1-p}の対数から計算されます.

Pのロジスティック変換\ \  logit(p)=log\frac{P}{1-P}

f(P)=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

グラフで確認してみます

pl <- function(p){ return(log(p/(1-p))) }
plot(pl,xlab="P",ylab = "f(P)")
abline(0,0)
ロジットf(p)がどのような値になっても、0<p< 1となっています.
f:id:yoshida931:20180213133241p:plain:w300

\ \  f(P) = b0+b1*x

ロジスティック回帰モデル

2値の変数(virginica=1、別の品種=0)が1になる確率のロジットを目的変数にしたものが、ロジスティック回帰分析になります。

ロジスティック回帰モデル

\ \  f(P) = b0+b1*x

P = \frac{1}{1+exp(-(b0+b1*x))}

パラメータの推定(最尤推定法) P(yi=1|b0,b1)

尤度関数
L(θ)= \prod_{i=0}^{150} (\frac{1}{1+exp(-(b0+b1x))})^{yi})*(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}

このままでは大変なので対数をとります
対数尤度関数 L(θ)= \sum_{i=0}^{150} log(\frac{1}{1+exp(-(b0+b1x))})^{yi}log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi} 

yiは0か1です.例えば例題のy1と y150を見てみます.それぞれy1=1, y150=0となっています.
y1のときは
f\theta (y_1)=log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}
y150のときは
f\theta (y_{150})=log(\frac{1}{1+exp(-(b0+b1x))})^{yi}
となります.これを、1~150全て掛け合わせることになります.

ここからNewton法により最尤法を解いていくことになります.
私には解説が困難ので以下のサイトをご覧ください

https://qiita.com/NaoyaOura/items/6ad5142b0306476d9293qiita.com

Rを使ったロジスティック回帰分析

ここではRを使って楽させていただきます.一般化線形モデル(glm:Generalized Linear Models)は最尤法によりパラメータを推定する関数です! Rの関数はglmそのままです.

# glm関数で、応答変数の分布は二項分布なので引数をbinomialにすればOKのはず・・・
(ans <- glm(y~x,family=binomial))

Error in weights * y : non-numeric argument to binary operator
というエラーが出たので、例題で使用したxyの属性を調べます

typeof(x);typeof(y)
[1] "double"
[1] "character"

xは数で、yが文字になっていました.ということでyを数に置きなおして再トライ

y <- as.numeric(y)
(ans <- glm(y~x,family=binomial))

Call:  glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
    -16.320        2.592  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3   AIC: 121.3

#やっぱりRは便利です

b0=-16.32, b1=2.596という解になりました

#対数尤度を求めてみます
x1 <- x[101:150]
x2 <- x[1:100]
e1 <- exp(-(-16.3198+2.592*x1))
e2 <- exp(-(-16.3198+2.592*x2))
(ly <- sum(log(1/(1+e1)))+sum(log(1-1/(1+e2))))

[1] -58.67273

視覚的にイメージしておきます.

par(mfrow = c(1,2))
plot(x,y,yaxt="n",
     ylab = "目的変数(0,1)",
     xlab = "Sepal.Length(がく片の長さ)(cm)")     #まずは散布図
name<-c("0","1")
axis(side=2,at=c(0,1),labels=name) 
abline(lm(y~x))    #回帰直線
p <- function(x){ return(1/(1+exp(-(-16.320+2.592*x)))) } #pの関数を作ります
plot(p,-16.320+2.592*x,ylim=c(0,1),xlim = c(3.3,9.3),
     ylab = "virginicaである確率",
     xlab="b0+b1*x") 
f:id:yoshida931:20180517150938p:plain:w700

Rを使用して最尤法で推定したパラメータの有意性も確認することが可能です. 関数summaryを使用します.

summary(ans)

Call:
glm(formula = y ~ x, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9870  -0.5520  -0.2614   0.5832   2.7001  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -16.3198     2.6581  -6.140 8.27e-10 ***
x             2.5921     0.4316   6.006 1.90e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 190.95  on 149  degrees of freedom
Residual deviance: 117.35  on 148  degrees of freedom
AIC: 121.35     # -2*対数尤度+2*パラメータ数

Number of Fisher Scoring iterations: 5

標準誤差の算出については、またいつか・・・

オッズについて

ロジスティック変換より

がく片の長さxcmのオッズ=\frac{P}{1-P} = exp(-16.3198+2.5921x)

がく片の長さxcmから1cm増えた場合のオッズ =\frac{P}{1-P} = exp(-16.3198+2.5921(x+1))

オッズ比= \frac{がく片の長さがx+1cmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}{ がく片の長さがxcmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}= exp(-16.3198+2.5921(x+1))-exp(-16.3198+2.5921x)=exp(2.5921)=13.35646

xが1増えた場合の対数オッズ比は2.59、オッズ比は13.36になります.

 summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900 

xには大きな変動がないので、xが0.5増えた場合のオッズ比を求めてみます

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}=exp(2.5921×0.5)=3.654832

つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます.オッズ比にはxが含まれていないので、どの幅(xcm)のに対しても同じことが言えます.


参考・引用

回帰分析入門 (Rで学ぶ最新データ解析)

回帰分析入門 (Rで学ぶ最新データ解析)

  • 作者: 豊田秀樹,尾崎幸謙,川端一光,岩間徳兼,鈴川由美,久保沙織,池原一哉,阿部昌利,大橋洸太郎,秋山 隆,大久保治信
  • 出版社/メーカー: 東京図書
  • 発売日: 2012/01/13
  • メディア: 単行本(ソフトカバー)
  • 購入: 2人 クリック: 36回
  • この商品を含むブログ (6件) を見る

相関係数のイメージ

パッケージmvtnormを使用して相関係数0.0, 0.2, 0.5, 0.7, 0.8, 0.9のグラフを作成してみます

install.packages("mvtnorm")
library(mvtnorm)

共分散行列.分散を全て1に設定しているので共分散=相関係数となります.

sigma00 <- matrix(c(1,0,0,1), ncol=2)  
sigma02 <- matrix(c(1,0.2,0.2,1), ncol=2)
sigma05 <- matrix(c(1,0.5,0.5,1), ncol=2)
sigma07 <- matrix(c(1,0.7,0.7,1), ncol=2)
sigma08 <- matrix(c(1,0.8,0.8,1), ncol=2)
sigma09 <- matrix(c(1,0.9,0.9,1), ncol=2)

乱数から各相関図に300個プロットします

n00 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma00)  
n02 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma02)  
n05 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma05)  
n07 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma07)  
n08 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma08)  
n09 <- rmvnorm(n=300, mean=c(0,0), sigma=sigma09)  

図を作成

par(mfrow = c(2,3))
plot(n00,xlab = "ρ=0.0",ylab="")  
plot(n02,xlab = "ρ=0.2",ylab="") 
plot(n05,xlab = "ρ=0.5",ylab="")    
plot(n07,xlab = "ρ=0.7",ylab="")
plot(n08,xlab = "ρ=0.8",ylab="")
plot(n09,xlab = "ρ=0.9",ylab="")

f:id:yoshida931:20180514142213p:plain:w800

変数の呼称について(目的変数と説明変数)

それぞれの研究界のご意見はあると思うのですが・・・
ややこしや

目的変数 Yは以下のように呼ばれています

目的変数 objective variable
応答変数 response variable 反応変数 reaction variable(response variable ) 結果変数 outcome variable
従属変数 dependent variable
基準変数 criterion variable
外的基準 external criterion
被説明変数 explained variable

説明変数 Xは以下のように呼ばれています

説明変数 explanatory variable
予測変数 predictor variable
独立変数 independent variable

パターンとしては

「目的変数&説明変数」
「従属変数&独立変数」
「被説明変数&説明変数」

個人的には…
応答変数&説明変数が理解しやすいかな?

2変量の正規分布をグラフでイメージ(persp)

また、ここで勉強させていただきました.
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html
忘れないように要点のみ転記させていただます.まさに備忘録.

今回はRの関数perspを使用して、密度関数の数式から3Dのグラフを描いてみます
確率変数x1、x2が正規分布に従い、無相関であることを仮定して進めていきます.

x1 <- seq(-3, 3, length=50)      # -3~3の範囲を50分割

head(x1)  #先頭部分だけ確認してみます
[1] -3.000000 -2.877551 -2.755102 -2.632653 -2.510204 -2.387755

x2 <- x1  

これで変数の設定は完了です.次に、分散1共分散0のマトリックスを作成します.

sigma.zero <- matrix(c(1,0,0,1), ncol=2)

     [,1] [,2]
[1,]    1    0
[2,]    0    1

Rの関数functionにx1、x2の同時確率密度を定義します. 上記sigma.zeroを分散共分散行列とする密度関数です  

f <- function(x1,x2) { dmvnorm(matrix(c(x1,x2), ncol=2), mean=c(0,0), sigma=sigma.zero) }

#6ポイントのみ確認してみます
f(-3,-3);f(0,0);f(3,3)
[1] 1.964128e-05
[1] 0.1591549
[1] 1.964128e-05

f(-1,-1);f(0,0);f(1,1)
[1] 0.05854983
[1] 0.1591549
[1] 0.05854983

関数outerを使用して、(x2,x2)であらわされる座標に対しFUN(f <- function(x1,x2))を適用します.
つまり、(x2,x2)に該当する値(f)が定まることになり3次元の描画を可能にします.

z <- outer(x1, x2, f)  
#50×50=2500個のデータが生成されました(length(z)で確認)
#10個だけ見てみましょう
z[1:10]
 [1] 1.964128e-05 2.814820e-05 3.973927e-05 5.526846e-05
 [5] 7.572219e-05 1.022015e-04 1.358875e-04 1.779878e-04
 [9] 2.296621e-04 2.919285e-04

http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html に書かれあります,
z の欠損値の置換は省略しております.
色塗りはRのヘルプ persp {graphics} をそのまま使用しました.

nrz <- nrow(z)
ncz <- ncol(z)
# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(x1, x2, z, col = color[facetcol], phi = 30, theta = -30)
f:id:yoshida931:20180512100244p:plain:w500

2変量の正規分布をグラフでイメージ(scatterplot3d)

ここで勉強させていただきました.
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html
忘れないように要点のみ転記させていただます.

必要なパッケージをインストールします

install.packages("mvtnorm")
library(mvtnorm)
install.packages("scatterplot3d")
library(scatterplot3d)

パッケージmvtnormは以下を参照
yoshida931.hatenablog.com

次に散布図を描きます

#共分散が0なので2変数には相関がありません
sigma.zero <- matrix(c(1,0,0,1), ncol=2)

     [,1] [,2]
[1,]    1    0
[2,]    0    1

#ともに平均=0,共分散=0となる2変数(乱数)を3組生成(n=100, n=1000, n=10000)
x100 <- rmvnorm(n=100, mean=c(0,0), sigma=sigma.zero)  
x1000 <- rmvnorm(n=1000, mean=c(0,0), sigma=sigma.zero)  
x10000 <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero)  

#3組並べて散布図
par(mfrow = c(1,3))
plot(x100)  
plot(x1000) 
plot(x10000)  

f:id:yoshida931:20180511164429p:plain:w500
確かに共分散はゼロです.

3Dで2変量の正規分布の同時確率密度関数を描いてみましょう

par(mfrow = c(1,3))
scatterplot3d(x100[,1], x100[,2], dmvnorm(x100, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)
scatterplot3d(x1000[,1], x1000[,2], dmvnorm(x1000, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)
scatterplot3d(x10000[,1], x10000[,2], dmvnorm(x10000, mean=c(0,0), sigma=sigma.zero), highlight=TRUE)

f:id:yoshida931:20180511165033p:plain:w600

Fisherの直接法

投稿日:2017-10-25、最終更新日:2018-05-02

Fisherの正確確率検定やFisherの直接確率検定、他にFisherの正確検定などと呼ばれています(統一してくれれば良いのにといつも思います).もともとカイ二乗検定は近似法でP値を求めています.一つのセルに度数が4以下が存在する場合には近似法の精度が落ちるため、直接P値を計算します.

超幾何分布

非復元抽出
袋の中にA+B個の玉が入っています.赤玉A個、白玉B個.その袋からn個取り出したときの赤玉の数をxとします.
f:id:yoshida931:20180502165147p:plain:w300

個数 残り 総数
赤玉 x A-x A
白玉 n-x B-(n-x) B
n A+B-n A+B

赤玉数A、白玉数B、個数nを与えることで、周辺和がすべて固定される
4つのセルのなかで、一つを決めると、他のセルも決まる(自由度は1)

赤玉の数がx個となる確率 P(X=x)=\frac{_AC_x\times_BC_{n-x}}{_{A+B}C_n}

上記クロス表のxが従う分布を超幾何分布といいます.
つまり、周辺和を与えたときの分割表が得られる確率分布のことです.

Fisherの直接法ではこの確率を使用します. それでは参考・引用文献の例題を参考にして勉強していきます. 上記のような2×2表が生起する確率を直接求めることになります.

Fisherの直接法

例) 小規模な二群並行ランダム化試験

有効 無効
対照群 7 8 15
試験群 12 2 14
合計 19 10 29

帰無仮説H_0: 対照群が有効となる確率 = 試験群が有効となる確率
対立仮説H_1: 対照群が有効となる確率 < 試験群が有効となる確率

上記のような実験結果から試験群の有効率が対照群の有効率より大きいことを証明します.まずは何も考えずRを使用してカイ二乗検定を行ってみます.

da <- c(7,12,8,2)
md <- matrix(da,2,2)

 md
     [,1] [,2]
[1,]    7    8
[2,]   12    2

ch <- chisq.test ( md, correct=F )
  
    Pearson's Chi-squared test

data:  md
X-squared = 4.8871, df = 1, p-value = 0.02706

Warning message:
In chisq.test(md, correct = F) :  カイ自乗近似は不正確かもしれません 

yoshida931.hatenablog.com

次に直接法で計算してみます.
対立仮説H_1: 対照群が有効となる確率 < 試験群が有効となる確率の方向を示すデータの確率を求めます.

     [有効] [無効]
[対照群]    a   b
[試験群]    c   d

c>12 となる確率が対立仮説の方向を強く示すデータとなります

a,b,c,d = 7,8,12,2
a,b,c,d = 6,9,13,1
a,b,c,d = 5,10,14,0

それぞれが生起する確率を求めます

c1 <- choose(15,7)
c2 <- choose(14,12)
c3 <- choose(29,19)
c1*c2/c3
 = 0.02923538

c21 <- choose(15,6)
c22 <- choose(14,13)
c23 <- choose(29,19)
c21*c22/c23
= 0.003498251

c31 <- choose(15,5)
c32 <- choose(14,14)
c33 <- choose(29,19)
 = 0.000149925

p-value = c31*c32/c33+c21*c22/c23+c1*c2/c3
        = 0.03288356

#Rの関数を使用
mx=matrix(c(7, 8, 12, 2), nrow=2, byrow=T)
fisher.test(mx,alternative = "less")  #ここでは割合の小さい対照群が上段にあるので"less"

    Fisher's Exact Test for Count Data

data:  mx
p-value = 0.03288
alternative hypothesis: true odds ratio is less than 1
95 percent confidence interval:
 0.00000 0.86222
sample estimates:
odds ratio 
 0.1565941 

p-value=0.03288356
対照群の有効数が12以上になる確率は <0.05 ということになりました.したがって、「a,b,c,d = 7,8,12,2」というサンプルは、帰無仮説のもとで起こりえない(有意水準0.05)組合せで、対照群が有効となる確率 < 試験群が有効となる確率という結果になります.

参考・引用

バイオ統計の基礎―医薬統計入門 (バイオ統計シリーズ)

バイオ統計の基礎―医薬統計入門 (バイオ統計シリーズ)