2015年度生物測定基礎実験

2015.9.24

統計解析２

東京大学大学院農学生命科学研究科　大森宏

この実験の目的

　統計解析ソフトＲを用いて，統計解析の理論と実践を学ぶ．

正規母集団からの標本に基づく推論

独立な正規分布の合成分布

• 和の分布
  X ＋ Y は平均 μ₁ ＋ μ₂，分散 σ₁² ＋ σ₂²，の正規分布に従う．
  X ＋ Y ～ N（ μ₁ ＋ μ₂， σ₁² ＋ σ₂² ）

• 差の分布
  X － Y は平均 μ₁ － μ₂，分散 σ₁² ＋ σ₂²，の正規分布に従う．
  X － Y ～ N（ μ₁ － μ₂， σ₁² ＋ σ₂² ）

• 一般の線形結合の分布
  a と b を任意の実数（スカラー）とすると，X と Y の線形結合 ax ＋ by は，
  aX ＋ bY ～ N（ aμ₁ ＋ bμ₂， a² σ₁² ＋ b² σ₂² ）

• 標本平均の分布
  特に，X₁，X₂，…，X_n を 平均 μ，分散 σ² の正規分布からの無作為標本であるとすると，標本平均 X^- の 平均と分散は，前節で示したようにそれぞれ，
  E[X^-] = E[ΣX_i/n] = ΣE[X_i]/n = (Σμ)/n = μ， Var[X^-] = Var[ΣX_i/n] = ΣVar[X_i]/n² = (Σσ²)/n² = σ²/n
  となる．よって，
  
  X^- ～ N（μ，σ²/n)
  となり，標本平均 X^- は，平均 μ，分散 σ²/n の正規分布に従う． この分散の平方根 σ/√n を標本平均 X^- の標準誤差（Standard Error, SE）と呼ぶ こともある．

正規分布に基づく母数の区間推定

分散既知の場合の母平均 μ の区間推定

• 95％ の意味
  　同じ正規母集団から標本抽出を繰り返すと，毎回標本平均として異なる値がえられ，それに 対応して信頼区間も異なる．この信頼区間の 95％ が真の平均 μ を含む，という意味である．
  つまり，100回の標本抽出により，100 個の信頼区間を作ったら平均的にみて，95 個の信頼区間が 真の平均 μ を含むことが期待できる．
  下の図は，平均 0 分散 2 の正規分布 N( 0, 2 ) から大きさ 10 の標本を取りだし，分散が既知であるとして， 母平均に対する信頼区間を 100 個生成したものである．"×" が標本平均を示す．左の "*" は，信頼区間 が母平均の真値 0 を含まなかった場合である．
  

  # 平均 μ の 95 ％信頼区間 100 回生成の R スクリプト

  v <- 2; n <- 10	# 母分散と標本の大きさ
  d <- qnorm(0.975)*sqrt(v/n)	# 95％信頼区間の幅
  x <- c(-2.5, -2.5, 2.5, 2.5)	# グラフの x の範囲
  y <- c(0, 100, 100, 0)	# グラフの y の範囲
  plot(x, y, type="n", xlab="", ylab="")	# グラフ領域確保
  segments(0, 0, 0, 100, col="red")	# 母平均のライン
  for(i in 0:100){	# 100回の繰り返し
  r <- rnorm(n, mean=0, sd=sqrt(v))	# N(0, v)からの大きさ n の標本平均
  m <- mean(r)	# 標本平均
  segments(m-d, i, m+d, i)	# 信頼区間の表示
  points(m, i, pch=4, col="red", cex=0.8)	# 標本平均の赤×表示
  if(m-d>0 || m+d<0) text(-2.5, i, "*")	# 信頼区間が母平均を含まない（失敗）した場合
  }	#
  title(main="N（ 0, 2 ）からの大きさ 10 の標本から得られた
  平均 μ の 95 ％信頼区間を 100 回作成", cex.main=1.0)	#

分散未知の場合の母平均 μ の区間推定

母平均の 95％ 信頼区間幅 ＝ t(n - 1)_0.975×標準誤差

• 母分散未知の時の母平均 95％ 信頼区間
  上のスクリプトで，

  m <- mean(r) 
  

  の部分の後に

  s <- var(r)
  d <- qt(0.975, df=9)*sqrt(s/n)
  

  を入れればよい．

仮説検定（Test of hypothesis）

帰無仮説（H₀）と対立仮説（H₁）

統計学で扱う仮説（hypothesis）とは，母集団に対する断定や推測．たとえば，

• 母集団は正規分布に従っている．
• 母集団平均は 0 である．
• 母集団 A と母集団 B の平均は等しい．

統計的仮説検定で用いられる仮説は，まず，帰無仮説（null hypothesis）という形式で与えられる．
帰無仮説は棄却（reject）されることに意味がある仮説である．
帰無仮説と反対の仮説を対立仮説（altanative hypothesis）という．

上の3番目の例でみると，

母集団 A と母集団 B は異なる処理（薬の投与など）をしているので，実験の目的 は，母集団 A と母集団 B の平均は異なる（処理効果がある）ことを言いたい （対立仮説が正しいことを望む）のだが，まずは 「等しい（処理効果無し）」と仮定してみようという考え方で， 数学の背理法と似た論理である．

背理法：√2 が無理数であることを証明するため，まず√2 が有理数であると仮定し，矛盾があることを 示す．つまり，有理数であることは絶対ありえない（確率 0 である！）ことを示す． この矛盾は，そもそも√2 を有理数とした仮定が誤っていたからであると考え，有理数という仮定を 棄却して，無理数であることを証明する．

検定の概要

検定統計量（Test statistic）

p 値（p value）と有意水準（Significance level）

　統計的仮説検定では，たとえば２つの母集団平均が等しいという帰無仮説を考えると， この帰無仮説のもとで，検定統計量（標本平均の差に基づく t 値など）以上 （もしくは未満）の値が得られる確率を求める．R ではこの値が p 値で表示される．
　p 値は くだけた言い方をすれば，帰無仮説が正しいとしたときに，標本のようなデータが得られる確率 である．
　これが十分小さい（ほとんどありえない）ときは，平均が等しいと仮定したことが誤りであったと判断して 帰無仮説を棄却し，２つの母集団平均には差があると結論づける．
　この確率がそれほど小さく ない場合は，このような統計量が得られることもありえると考え，帰無仮説を採択し，平均が等しいと考え てもよいとする．
　棄却か採択かの判断の基準となる確率を有意水準といい， 5 % や 1 % がよく用いられる．

　統計のソフトが発達していなかった頃は，検定統計量である t 値や F 値を電卓等で算出し，その値を t 分布や F 分布の 5 ％や 1 ％の有意水準に対応する数表と照らし合わせて検定を行い，5 ％有意とかを 記述していた．
　現在では，ソフトが検定統計量に対する p 値を直接計算してくれるので数表はいらなくなった．この結果， 検定に重要な数字であった t 値や F 値より，より直接的な p 値が重要な指標になってきた．p 値を みれば，何％有意かが一目でわかるので，結果にわざわざ 5 ％有意とかを記述する必要がなくなってきており， 論文の書き方も変わってきている．

片側検定と両側検定

実験状況によっては，薬投与などの処理を行った集団（処理群）平均 μ_A が，薬を投与しない 集団（対照群）の平均 μ_B より小さくなることはないことが事前に わかっているような場合が ある．このようなとき，

両側検定と信頼区間

母集団平均に対する両側検定は，母集団平均に対する信頼区間と大きな関係がある． いま，帰無仮説（H₀）と対立仮説（H₁）が，

 

 

帰無仮説を受諾 ⇔

検定における２種類の過誤

検定は，仮説を棄却するか採択するかのいずれかであるが， 統計量は分布をもつので，この判定には間違いが起こることがある．
以下のように，この過誤には 2 種類がある．

統計的検定における２種類の過誤

	仮説の棄却（reject）	仮説の採択（accept）
仮説が真（true）のとき	第１種の過誤（Type 1 error）	正解
仮説が偽（false）のとき	正解	第２種の過誤（Type 2 error）

第１種の過誤が有意水準である．また，第２種の過誤の確率を β としたとき， 仮説が偽のとき正しく仮説を棄却する確率，1 - β，を検出力，もしくは検定力（power）という． よい検定は，第１種の過誤を固定したもとで検出力の高い検定方式である．

　仮説検定とまったく同一の論理構造を持ったものとして，刑事裁判があげられる．真実は，被告人 が無実か犯人かのどちらかであるが，これは本人以外は知り得ないので，裁判により無罪か有罪のいずれか を判断する．このときの帰無仮説と対立仮説は，

帰無仮説（H₀）：被告人は無罪である（無罪推定）
対立仮説（H₁）：被告人は有罪である．

刑事裁判における２種類の過誤

	有　罪	無　罪
被告が無実のとき	第１種の過誤（Type 1 error）	正解
被告が犯人のとき	正解	第２種の過誤（Type 2 error）

１つの母集団に対する検定（One sample problem）

正規母集団の母平均に対する t 検定

なお，この検定は，対のある標本に適用できる．対のある標本とは，n 組のペアー標本（paired smple），

• t 検定の有意水準
  　t 検定の有意水準をシミュレーションにより確認してみる． 平均 μ = 10，分散 σ² = 4 の正規分布から大きさ n = 10 の標本を抽出し， その標本平均 x^- と標本（不偏）分散 s² から t 値を計算する． これを N = 10000回繰り返すと，t 値のヒストグラムができる．このヒストグラムは自由度 9 の t 分布に従っているはずである．
  
  帰無仮説，H₀：μ = 10， 対立仮説，H₁：μ ≠ 10，
  の検定では， t 値の絶対値 |t| が 自由度 9 の t 分布の97.5％分位点 t(9)_0.975 より大きくなった場合に棄却される． このため，この個数を数えれば第１種の過誤の確率，つまり有意水準がシミュレーションにより推定できる． N = 10000回のシミュレーションでは，これがほぼ 5 ％になっていた．
  

  #t 検定の有意水準 の R スクリプト

  N <- 10000	# シミュレーション回数
  n <- 10	# サンプルサイズ
  n1 <- rnorm(N*n, mean=10, sd=2)	# N(10, 4) から大きさ n の標本を N 回シミュレーション
  n1.mat <- matrix(data=n1, ncol=n)	# データ行列
  n1.mean <- apply(n1.mat, 1, mean)	# 各サンプルの平均
  n1.var <- apply(n1.mat, 1, var)	# 各サンプルの分散
  n1.var3 <- n1.var/n	# 各標本平均の分散
  n1.td <- (n1.mean - 10)/sqrt(n1.var3)	# 各サンプルの t 値
  mean(n1.td)	# t 値の平均
  sd(n1.td)	# t 値の標準偏差
  mt <- ceiling(max(abs(n1.td)))	# t 値の絶対値の最大値
  xq1 <- qt(0.025, df = (n-1))	# t(9) の 2.5% 点
  xq2 <- qt(0.975, df = (n-1))	# t(9) の 97.5% 点
  length( n1.td[abs(n1.td) > xq2])	# 5% 検定で有意となった個数
  hist(n1.td, breaks=seq(-mt, mt, by=0.2), probability=TRUE, xlab="t 値", ylab="密度", main="")
  curve(dt(x, df=(n-1) ), -5,5, col="red", add=T)	# t 分布
  abline(v=xq1, col="red")	# 採択域の下限
  abline(v=xq2, col="red")	# 採択域の上限
  title(main="Histogram of t values")	# タイトル

  問題２のヒント
  

  n1 <- rnorm(N*n, mean=10, sd=2)		# N(10, 4) から大きさ n の標本を N 回シミュレーション
  

  のところを

  n1 <- rnorm(N*n, mean=11, sd=2)		# N(11, 4) から大きさ n の標本を N 回シミュレーション
  

  に変えれば，本当は母集団は N(11, 4) だが，μ= 10 と仮定したときの棄却回数が求められる．

• 非心 t 分布（non-central t distribution）
  　帰無仮説が正しいときに，t 値の従う分布は t 分布に従うことを見た．帰無仮説が正しくない ときに t 値の従う分布を非心 t 分布という．いま，母集団分布を N(μ，σ²) とし たとき，
  
  帰無仮説：H₀：μ = μ₀，対立仮説：H₁：μ ≠ μ₀
  の t 検定を，大きさ n の標本で行う場合を考える． このとき標本平均が x^-，標本標準偏差が s であったとすると， 検定統計量
  t = √n(x^- - μ₀)/s
  の分布は，帰無仮説が正しくないとき，自由度 n - 1， 非心パラメーター（non-central parameter）
  δ = √n(μ - μ₀)/σ
  の非心 t 分布に従う．この分布を用いると，t 検定の検出力を求めることができる．
  

  # 非心 t 分布の R スクリプト

  n <- 10 # サンプルサイズ m0 <- 10 # 帰無仮説の平均 s <- 2 # 標準偏差　 m <- 11 # 真の平均 delta <- sqrt(n)*(m-m0)/s # 非心度 curve(dt(x, df=(n-1)), -4,7) curve(dt(x, df=(n-1), ncp=delta), -4,7, col="red", add=T) title(main="non-central t distribution")

• t 検定の検出力（power）
  　前節では帰無仮説が真のとき，t 検定により仮説が棄却する確率を求めた．では，帰無仮説が偽であるとき， 帰無仮説がどれくらいの確率で正しく棄却できるであろうか．すなわち，t 検定の検出力はどの程度であろうか． これは非心 t 分布で求められるが，シミュレーションなら比較的わかりやすく求めることができる． 後ほど問題２で検出力を求める．
  

ノンパラメトリック検定

　ある連続分布をもつ母集団から大きさ n の無作為標本 X₁，…，X_n が抽出されたとする．ここで，母集団分布の位置母数 （location parameter）に関する両側検定問題

Wilcoxon の符号順位和検定（signed rank sum test）

• Wilcoxon の符号順位和検定の有意水準

  問題１：データ数 n = 4 のとき，タイがないものとしてデータの取りえるすべての 可能性を考えることにより，T の帰無仮説のもとでの分布を出せ． また，T ＝ 0 のとき，帰無仮説を棄却すると決めたとすると，この検定の有意水準を求めよ．

  ヒント

  n = 4 のとき，値の絶対値を大きさの順に並べて，X₍₁₎，X₍₂₎， X₍₃₎，X₍₄₎，としたとき， 各 X_(i) は＋由来か－由来かの２通りがあるので， 2⁴ = 16 通りの場合がある．このおのおのの場合 について T の値を出すことにより，T の確率分布が求められる．

• t 検定と符号順位和検定の検出力の比較

  問題２：t 検定の有意水準の R スクリプトを参考にして，t 検定の検出力を求めよ．

  1. 標本の大きさ n = 10 として，帰無仮説が偽であり， 真の母集団平均が μ = 11，12，13，の３通りについて，帰無仮説；H₀：μ = 10，の検定 を行った場合の 検出力（帰無仮説を棄却する割合）を求めよ．
  2. 真の母集団平均が μ = 11 であるとき，検出力を 90％以上にするために 必要なサンプルサイズを求めよ．
  3. 真の母集団平均が μ = 11 であるとき，符号順位和検定を行ったときの 検出力を 90％以上にするために 必要なサンプルサイズを求めよ．

     ヒント：ウイルコックス検定を1000回やって，p値が0.05以下の個数を出す

     N <- 10000
     n <- 10
     n1 <- rnorm(N*n, mean=11, sd=2) 	# N(11, 4) から大きさ n の標本を N 回シミュレーション 　
     n1.mat <- matrix(data=n1, ncol=n) 	# データ行列 　
     pvalue <- NULL
     for(i in 1:N){
       wt <- wilcox.test(n1.mat[i,]-10)
       pvalue <- c(pvalue, wt$p.value)
     }
     length(pvalue[pvalue<0.05])
     

• ラット卵巣のアスコルビン酸濃度
  　生物統計学入門より抜粋．ある組織抽出物が卵巣のアスコルビン酸の量を変化させるかを調べた． ラットには左右 1 対の卵巣があるので，まず片方の卵巣を取ってアスコルビン酸濃度を測り，次に抽出物 を注射して一定時間後に残り片方の卵巣を取ってアスコルビン酸濃度（μg/1000mg 組織）を測った．8 頭の ラットの実験結果は以下のようであった．

  データダウンロード

  　ここで，組織抽出物はアスコルビン酸を変化させないという帰無仮説
  H₀ : アスコルビン酸減少量 = 0
  の検定を行ってみる．

  #１つの母集団に対する検定の R スクリプト

  rat <- read.csv("rat.csv"); rat #データ入力 d <- rat$before-rat$after #減少量 boxplot(d, xlab="decreased concentration") np <- length(d[d > 0]) #＋の個数 nm <- length(d[d < 0]) #－の個数 binom.test(c(np,nm)) #符号検定 mean(d) #平均 sd(d) #標準偏差 t.test(d) #１標本 t 検定 t.test(rat$before, rat$after, paired=TRUE) #同じ検定 wilcox.test(d) #Wilcoxon の符号順位和検定 wilcox.test(rat$before, rat$after, paired=TRUE) #同じ検定

• 園芸療法の効果
  　園芸療法（horticultural therapy）は第2次大戦後，兵士のリハビリテーション等を通して欧米で発達してきた． 日本でも90年代から身体障害者，精神障害者，高齢者等の生活の質（QOL, Quality Of Life）や セルフエスティーム（self-esteem，自尊感情）を向上させる手段として注目されてきていている．
  　園芸療法の数値的･客観的な効果の計測を行うため，認知症診断に用いられる MSE（Mental State Examination，心理機能検査），CDT（Clock Drawing Test，時計描画テスト）と 老年者や認知症患者の日常生活の遂行能力を測る NM スケール（N 式老年者精神状態評価尺度）， N-ADL（N 式老年者日常生活動作能力評価尺度），および，気分の状態（抑うつ，活気のなさ，怒り，疲労，緊張， 混乱）を計測する POMS（Profile Of Mood State），自己記入式 QOL 質問表 QUIK を用いた． これらの調査を園芸療法の前後でクライアントに対して行い，効果の数値的計測を行った．
  　ここでは，MSE の結果についての解析を行う．MSE は 30 点満点で，点数が高いほど日常生活に適応して いると考えられている．MSE に関しては前後とも回答したクライアントの総数は 19 例であった． データは表２にまとめられている．元データに差と符号と順位も付け加えた．このデータを用いて，園芸療法 を施す前と後では，心理機能に差がないという帰無仮説，すなわち，
  H₀：MSE 得点の差 ＝ 0
  の検定を行う．

  データダウンロード

  文献
  

  • 鈴木修（編），大森宏，児玉良治，渡辺俊之，矢野広，山根健治（2004）． 専修学校における園芸療法士教育育成システムの研究開発 （文部科学省委託平成15年度「専修学校先進的教育研究開発事業」）
  • 鈴木修，渡辺俊之，矢野広，山根健治，大森宏，伊東正信，最上正秀，山下容子，小泉力， 児玉良治，頭士智美，細井薫，水口聡子，遠藤久子，樋田奈穂子，小島ユリ，郡司敏幸（2005）． 福祉サービス提供者に対する園芸療法教育システムの研究開発 （平成16年度文部科学省「専修学校社会人キャリアアップ教育推進事業」）

  問題３：MSE 得点データから園芸療法の効果を，符号順位和検定，t 検定を行って 検定せよ．また，クライアントの年齢を85歳以下に限定した場合の効果の検定を行い，結果の考察を行なえ．

  ヒント：85歳以下を取り出すには以下のようにする．

  mse <- read.csv("mse.csv"); mse		#データ読み込み
  d <- mse$after - mse$before 		#処方後得点－処方前得点
  young <- which(mse$age < 86)
  d2 <- d[young]
  

２つの母集団に対する検定（Two sample problem）

母平均の同等性に対する t 検定

　母集団 A から大きさ n_A，母集団 B から大きさ n_B の標本を抽出した． 母集団 A からの標本の標本平均が x^-_A， 標本分散が s_A² であり，母集団 B の 標本平均が x^-_B， 標本分散が s_B² であった．母集団 A，B が共通の 分散 σ² をもつとすると，その推定値 s² は 以下のように推定できる．

母集団 A からの標本の偏差平方和：	SA＝（nA－1）sA2
母集団 B からの標本の偏差平方和：	SB＝（nB－1）sB2
母集団 A，B 全体での偏差平方和：	S ＝ SA ＋ SB ＝（nA－1）sA2＋ （nB－1）sB2
母集団 A，B 共通の標本分散：

帰無仮説（H₀： μ_A ＝ μ_B）のもとでは，μ_A－μ_B＝0，なので， 標本平均の差は，

母集団A，Bで分散の同等性が疑われるときは，ウェルチ（Welch）の検定を用いる．

ノンパラメトリック検定

Mann-Whitney の U 検定

Kolmogorov-Smirnov の 2 標本検定

表．Kolmogorov-Smirnov 2 標本両側検定の臨界値

確率：γ	0.99	0.95	0.90	0.85	0.80
臨界値：kγ	1.63	1.36	1.22	1.14	1.07

２つの母集団分散の同等性の検定

　母集団 A から大きさ n_A，母集団 B から大きさ n_B の標本を抽出した． 母集団 A からの標本の標本平均が x^-_A， 標本分散が s_A² であり，母集団 B の 標本平均が x^-_B， 標本分散が s_B² であるとする．すると， 標本分散に関係した量はそれぞれ

　ところで，帰無仮説が正しいとする と，σ_A² ＝ σ_B² とおけるので， 母集団の分散比は，γ₀ = σ_A²/σ_B² = 1， となる．このとき，標本分散の分散比の統計量 γ が，

• 淡水性ウナギデータの解析
  淡水性ウナギの汽水域での生理活性の違いのデータ を２標本データとして違いを検定してみよう．ただし，２母集団で分散の違いがあるようなので，まず， 分散の同等性の検定を行い，違いが認められるときは，データの対数変換（lx <- log(x)）を行うなどして，分散安定化 変換を行って t 検定を行ってみる．また，データの対数変換などを行わず．ウェルチ検定やノンパラメトリック 検定を行った場合との比較検討を行う．

  データダウンロード

  分散同等性の検定をやってみる．

  unagi <- read.csv("unagi2.csv")
  attach(unagi)
  boxplot(unagi)
  var.test(noNa, Na)
  

  

  分散同等性検定で，p 値 = 0.047 で有意に分散が異なっていることがわかったので， 通常の t 検定を行うことはできない．

  t.test(noNa, Na)
  wilcox.test(noNa, Na)
  ks.test(noNa, Na)
  plot(ecdf(noNa), do.points=F, verticals=T,xlab="生理活性", ylab="累積確率", 
  xlim=range(noNa, Na), main="") 
  plot(ecdf(Na), do.points=F, verticals=T, add=T, col.h='red', col.v='red') 
  

  検定方法	p 値
  分散同等性検定	0.047	*
  Welch 検定	0.01355	*
  順位和検定	0.01452	*
  Kolmogorov-Smirnov 検定	0.07546

  

  Kolmogorov-Smirnov 検定以外では5％有意となり，ウナギの生理活性は淡水と汽水で有意に異なっていた． なお，Kolmogorov-Smirnov 検定の検出力は低いことがわかる．

  分散安定化変換を試す．

  x <- log(noNa)
  y <- log(Na)
  boxplot(x, y)
  var.test(x, y)
  t.test(x, y, var.equal=T)
  

  検定方法	p 値
  分散同等性検定	0.5071
  t 検定	0.01126	*

  

  分散安定化変換を行ったところ，分散の同等性は棄却されず，普通の t 検定を行った．

  昔の指針：２つの母集団の分散の均質性が疑われるときは分散の同等性の検定を行い，これが棄却された ときは仕方なく Welch 検定（自由度の計算が面倒で，小数点自由度のときは t 分布表の線形補間で97.5％を計算 する）を行う．
  今の指針：最初から Welch 検定でよいのでは．

  問題４：園芸療法の MSE 得点データを，園芸療法を行わなかった対照区 19 名と園芸療法を処方 した処理区 19 名のデータとみなし，t 検定，Mann-Whitney の U 検定， Kolmogorov-Smirnov 検定を行え． また，クライアントの年齢を85歳以下に限定した場合の効果の検定を行い，結果の考察を行なえ．

分割表の解析

成功確率（比率）の同等性の近似検定

• 勝率の同等性の検定
  　A と B で将棋を 30 試合行い，A が 21 勝 9 敗の成績であった． 続いて，A と B で囲碁を 30 試合行い，A が 12 勝 18 敗の成績であったとする．このとき， A の B に対する将棋の強さと囲碁の強さとで違いがあるかの検定を行う．帰無仮説は，
  H₀：A の B に対する将棋の勝率 ＝ A の B に対する囲碁の勝率
  である．

  # 勝率同等性検定の R スクリプト

  r <- c(21, 12)	# A の将棋の勝数と囲碁の勝数
  n <- c(30, 30)	# 将棋の試合回数と囲碁の試合回数
  prop.test(r, n, correct=F)	# 比率の同等性検定（連続性の補正なし）
  prop.test(r, n)	# 比率の同等性検定（連続性補正）

２×２分割表

2 × 2 分割表データ

	成　功	失　敗	計
試行 X	a	b	a + b
試行 Y	c	d	c + d
計	a + c	b + d	N

　周辺度数を固定したときに，帰無仮説のもとでの各セル期待度数は，

帰無仮説のもとでの期待度数

	成　功	失　敗
試行 X	(a + b)(a + c)/N	(a + b)(b + d)/N
試行 Y	(c + d)(a + c)/N	(c + d)(b + d)/N

となる．ここで，ピアソン（Pearson）のχ² 値，

• 勝率同等性の χ² 検定
  

  # 勝率同等性ピアソン χ² 検定の R スクリプト

  shogi <- c(21, 9)	# A の将棋の勝ち負け数
  igo <- c(12, 18)	# A の囲碁の勝ち負け数
  x <- rbind(shogi, igo)	# 行連結による行列化
  chisq.test(x, correct=F)	# 勝率同等性 χ2 検定（連続性の補正なし）
  chisq.test(x)	# 勝率同等性 χ2 検定（連続性の補正）

Fisher の直接確率法

将棋と囲碁 2 × 2 分割表

	勝　ち	負　け	計	場合の数
将　棋	7	3	10	10C7
囲　碁	4	6	10	10C4
計	11	9	20	20C11

パターン１ 勝ち 負け 計   場合の数 将　棋 8  2  10    10C8 囲　碁 3  7  10    10C3 計   11  9  20    20C11   確　率  10C8×10C3 /20C11 = 0.03215		パターン２ 勝ち 負け 計   場合の数 将　棋 9  1  10    10C9 囲　碁 2  8  10    10C2 計   11  9  20    20C11   確　率  10C9×10C2 /20C11 = 0.00268		パターン３ 勝ち 負け 計   場合の数 将　棋 10  0  10    10C10 囲　碁 1  9  10    10C1 計   11  9  20    20C11   確　率  10C10×10C1 /20C11 = 0.00006
パターン４ 勝ち 負け 計   場合の数 将　棋 6  4  10    10C6 囲　碁 5  5  10    10C5 計   11  9  20    20C11   確　率  10C6×10C5 /20C11 = 0.315075		パターン５ 勝ち 負け 計   場合の数 将　棋 5  5  10    10C5 囲　碁 6  4  10    10C6 計   11  9  20    20C11   確　率  10C5×10C6 /20C11 = 0.315075		パターン６ 勝ち 負け 計   場合の数 将　棋 4  6  10    10C4 囲　碁 7  3  10    10C7 計   11  9  20    20C11   確　率  10C4×10C7 /20C11 = 0.15004
パターン７ 勝ち 負け 計   場合の数 将　棋 3  7  10    10C3 囲　碁 8  2  10    10C8 計   11  9  20    20C11   確　率  10C3×10C8 /20C11 = 0.03215		パターン８ 勝ち 負け 計   場合の数 将　棋 2  8  10    10C2 囲　碁 9  1  10    10C9 計   11  9  20    20C11   確　率  10C2×10C9 /20C11 = 0.00268		パターン９ 勝ち 負け 計   場合の数 将　棋 1  9  10    10C1 囲　碁 10  0  10    10C10 計   11  9  20    20C11   確　率  10C1×10C10 /20C11 = 0.00006

• 勝率同等性の Fisher 直接確率法検定
  

  shogi <- c(7, 3)  		# A の将棋の勝ち負け数 　 
  igo <- c(4, 6)  		# A の囲碁の勝ち負け数 　 
  x <- rbind(shogi, igo)  # 行連結による行列化 　 
  m <- sum(shogi)  		# 将棋の勝負回数 　 
  chisq.test(x)  			# 勝率同等性 χ2 検定．警告が出る． 
  fisher.test(x)  		# 勝率同等性 Fisher 直接確率 　 
  
  
  

• 超幾何分布による Fisher 直接確率法検定の p 値
  ところで， 将棋を m = 10 回，囲碁を n = 10 回勝負し，全部で k = 11 回の「勝ち」があったとき，将棋では x 回の「勝ち」が得られる確率分布は，超幾何分布 f(x; m, n, k) に従う．
  　超幾何分布で x = 7 回の「勝ち」が得られる確率以下である確率をすべて加え合わせたのが p 値となる．
  

  問題５：勝敗の比率はそのままで，試行回数を増やした（n，m = 30）ときの 勝率同等性 Fisher 直接法の p 値を，超幾何分布の確率密度関数 dhyper(x, m, n, k) か累積分布関数 phyper(x, m, n, k) 用いて算出せよ．

クロス集計と独立性の検定

2 × 2 分割表データ

	反応１	反応２	計
集団１	n11	n12	n1・
集団２	n21	n22	n2・
計	n・1	n・2	n

　ここで，n_ij は，集団 i の標本の中で，反応 j を取った人数（度数） である．また，n_i・ は集団 i の標本の大きさで， n_・j は標本全体（大きさ n）の中で反応 j を取った人数 を表し，それぞれ周辺度数という．
　このような表において，集団で反応パターンに違いがなければ，集団 i で反応 j を取る確率は， 集団 i である確率 n_i・/n に 集団 j である確率 n_・j/n をかけた n_i・n_・j/ n² となることが期待される．これより， 集団 i で反応 j を取る人数（度数）は， n_i・n_・j/n となることが期待される．これを独立性の仮定という．
　いま，集団 i の標本の中で，反応 j を取る確率を p_ij とし，集団 i の周辺確率を p_i.， 反応 j の周辺確率を p_.j とすると，独立性の仮説は，

オッズ比（odds ratio）

• 心筋梗塞死亡率の χ² 独立性検定とオッズ比
  　ある病院で，心筋梗塞発作後の６ヶ月以内の死亡率と 牛乳に対する抗体を持った患者と持っていない患者でまとめたのが下の表である． 心筋梗塞死亡リスクと牛乳に対する抗体の有無との関連を調べる．

  心筋梗塞死亡率と牛乳に対する抗体の有無

  	６ヶ月以内死亡	６ヶ月生存
  牛乳に対する抗体あり	29	71
  牛乳に対する抗体なし	10	94

  # 心筋梗塞死亡率の R スクリプト

  x <- matrix(c(29, 71, 10, 94), nrow=2, byrow=T); x chisq.test(x) # χ2 値，p 値 fisher.test(x) # オッズ比，p 値

• 食中毒原因食材の χ² 独立性検定とオッズ比
  1940 年のNew York 州 Oswego の協会の夕食会における胃腸炎異常発生の喫食 調査データによると食材と食中毒症状とで以下の関係があった．
  詳しいデータ

  食品名	発症あり		発症なし
  	食べた	食べない	食べた	食べない
  ケーキ	27	19	13	16
  ハムステーキ	29	17	17	12
  バニラアイス	43	3	11	18
  チョコアイス	25	20	22	7
  フルーツサラダ	4	42	2	27

  問題６： 食材ごとの発症データから，χ² 独立性検定とオッズ比を使って 食中毒の原因食材を推定せよ．すなわち，オッズ比が 1 より有意に大きくなった食材が食中毒の 原因食材と特定できる．オッズ比が 1 より小さな食材は食中毒発生を抑制した食材である． このような食材があった場合，どうして抑制効果があったのか考察せよ．

参考文献（古い順）

1. Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
2. 「実験」生産環境生物学，東京大学大学院農学生命科学研究科生産・環境生物学専攻編，1999，朝倉書店
3. 工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004， 数理工学社
4. 実践生物統計学－分子から生態まで－（第 1 章，第 2 章）， 東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
5. The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
6. R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
7. R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社