東京農工大学

2020.10.23

生物統計学２

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

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

独立な正規分布の合成分布（正規分布の再生性）

• 和の分布
  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 を 平均 μ，分散 σ² の正規分布からの大きさ n 無作為標本であるとすると，標本平均 の平均と分散は，前節で示したようにそれぞれ，
  
  
  となる．よって，
  
  
  となり，標本平均 は，平均 μ，分散 σ²/n の正規分布に従う． この分散の平方根 を標本平均 の標準誤差（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)	#

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

仮説検定（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）

有意水準（Significance level）

p 値（p value）

tv <- 1.58			# t 値
p <- 1 - pt(tv, df=19)		# t 値より大きい確率
2*p				# p値
curve(dt(x, df=19), -4,4)		# 自由度19のt分布
abline(h=0)
x1 <- seq(tv, 4, by=0.1)
y1 <- dt(x1, df=19)
x2 <- seq(-4, -tv, by=0.1)
y2 <- dt(x2, df=19)
segments(tv,0, tv, dt(tv, df=19), col="red")	# t 値
# p 値の部分を色づけ
polygon(c(x1,3,tv), c(y1,0,0), col="red", density=20)
polygon(c(x2,-tv,-3),c(y2,0,0), col="red", density=20)

片側検定と両側検定

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

両側検定と信頼区間

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

 

 

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

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

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

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

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

統計的検定と刑事裁判

　統計的検定と同じロジックであるのが、刑事裁判である。統計的検定における「処理効果が無い」 という帰無仮説は、刑事裁判での「被告人は無罪である（無罪推定）」に対応する。 統計的検定では「帰無仮説」を疑うような、「処理効果が有意に存在することを示すデータ」があれば、 帰無仮説は棄却され、処理効果が認められる。刑事裁判においては、「無罪推定」を覆すような 「客観的証拠」があれば、無罪推定は棄却され有罪となる。

刑事裁判

	有罪（regal）	無罪（irregal）
被告が無実のとき	えん罪（第１種の過誤）	正当
被告が犯人のとき	正当	無罪放免（第２種の過誤）

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

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

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

• t 検定の例
  　いま、平均 μ= 11、分散 σ² = 4 の正規母集団 N(11, 4) から 大きさn = 10の標本を抽出したとき、
  帰無仮説，H₀：μ = 10， 対立仮説，H₁：μ ≠ 10，
  の検定を平均、分散が未知として行ってみる。母集団の平均が11なので、この帰無仮説は偽（false）であり、棄却してくれるとありがたい。

  n <- 10				# サンプルサイズ
  m0 <- 10 ; m1 <- 11		# 帰無仮説における平均と真の平均
  x <- rnorm(n, mean=m1, sd=2) ; x	# N(m1, 4)からn個のサンプリング
  t.test(x, mu=m0)			# μ= m0としたt検定
  

  なお、Rでは t 検定は t.test()で行うことができる。出力は以下のようになった。

  One Sample t-test
  
  data:  x
  t = 1.1159, df = 9, p-value = 0.2934
  alternative hypothesis: true mean is not equal to 10
  95 percent confidence interval:
    	9.360708 11.884105
  sample estimates:
  mean of x 
   	10.62241 
  

  　この標本のもとで、t 値= 1.1159、自由度= 9 = n - 1、p値= 0.2934、 であったことがわかる。p 値が0.05より大きいので 5％有意ではなく、 帰無仮説は棄却できなかった。人によっては、有意になった人もいる。また、平均 μ の95％信頼区間は、 9.360708＜μ＜11.884105、であり、帰無仮説の μ＝10を含むので棄却することはできなかった。なお、標本平均 = 10.62241、であった。
  　これらの値を具体的にデータから計算してみると、

  m <- mean(x); m		# 標本平均
  df <- n-1; df		# t分布の自由度
  se <- sd(x)/sqrt(n)	# 標本平均の標準誤差
  tv <- (m - m0)/se; tv	# t値
  pv <- 2*(1 - pt(tv, df)); pv	# p値
  tq <- qt(0.975, df)	# t(n -1)の97.5％分位点（両側5％点）
  m - tq*se		# 母平均に対する95％信頼区間の下限
  m + tq*se		# 母平均に対する95％信頼区間の上限
  

  となる。
• t 検定の検出力
  　先ほどの t 検定では、帰無仮説が正しくなかったにもかかわらず棄却できなかった。 これは、この検定の検出力が低かったためであると考えられる。 そこで、先ほどの検定の検出力を計算してみよう。 これは、N(11, 4) から大きさ n = 10 の標本を N=10000回抽出して 帰無仮説（H₀：μ= 10）の t 検定を行い、 5％有意になった回数を数えればわかる。
  

  n <- 10; N <- 10000		# 標本の大きさとシミュレーション回数
  pvalue <- tvalue <- numeric(0)
  for(i in 1:N){
  	x <- rnorm(n, mean=m1, sd=2)	# N(m1, 4)からn個のサンプリング
  	res <- t.test(x, mu=m0)
  	pvalue <- c(pvalue, res$p.value)	# t検定のp値
  	tvalue <- c(tvalue, res$statistic)	# t検定のt値
  }
  length(pvalue[pvalue<0.05])	# p値で0.05未満（5％有意）の個数
  tq <- qt(0.975, n-1)
  hist(tvalue, nclass=50)
  abline(v=tq, col="red")		# 棄却域の閾値
  

  これによると、この検定の検出力は30%程度であることがわかった。検出力を上げるには標本の大きさを大きくする必要がある。

ノンパラメトリック検定

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

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

表１．順位のあり得る組み合わせ（n＝4）

１	２	３	４	R+	R-	T
＋	＋	＋	＋	10	0	0
－	＋	＋	＋	9	1	1
＋	－	＋	＋	8	2	2
＋	＋	－	＋	7	3	3
＋	＋	＋	－	6	4	4
－	－	＋	＋	7	3	3
－	＋	－	＋	6	4	4
－	＋	＋	－	5	5	5
＋	－	－	＋	5	5	5
＋	－	＋	－	4	6	4
＋	＋	－	－	3	7	3
－	－	－	＋	4	6	4
－	－	＋	－	3	7	3
－	＋	－	－	2	8	2
＋	－	－	－	1	9	1
－	－	－	－	0	10	0

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

  データダウンロード

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

  rat <- read.csv("rat.csv"); rat  	#データ入力
  d <- rat$before-rat$after			#減少量
  boxplot(d)
  stripchart(d, method="jitter")		# jitterは標本を揺らす
  t.test(d)					#１標本 t 検定
  t.test(rat$before, rat$after, paired=TRUE)	#同じ検定
  wilcox.test(d)				#Wilcoxon の符号順位和検定 
  wilcox.test(rat$before, rat$after, paired=TRUE)    #同じ検定
  

２つの母集団に対する検定（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，なので， 標本平均の差は，

ウェルチ（Welch）検定

ノンパラメトリック検定

Mann-Whitney の U 検定

• 淡水性ウナギの汽水域での生理活性の違い
  　淡水性ウナギの汽水域での生理活性の違いのデータ、ウナギデータ，を２標本データとして違いを検定してみよう。

データダウンロード

unagi <- read.csv("unagi2.csv"); unagi  	#データ入力
boxplot(unagi)
stripchart(unagi)
x <- unagi$noNa; y <- unagi$Na
t.test(x, y) 		# ２標本 Welch 検定 
wilcox.test(x, y)			# Mann-Whitney の U 検定

適合度検定（goodness-of-fit tests）

Peason のχ² 値

Mendelの分離の法則（the law of segregation）

　種子の形について、優性：劣性 = 3:1となっているかみてみよう。観測度数と期待度数は以下のようにまとめられる。

ni <- c(5474, 1850)
n <- sum(ni)	
p <- c(3/4, 1/4)	# 分離の法則から想定される分離確率
npi <- n*p; npi		# 分離の法則から想定される期待度数
v <- (ni-npi)^2/npi; v		# (O - E)^2 / E の計算
ch2 <- sum(v); ch2		# Peason カイ二乗値
1 - pchisq(ch2, df = length(ni) - 1)	# カイ二乗分布上側確率
chisq.test(ni, p = p)   # Rの備え付け関数でもできる。

Mendelの独立の法則（the law of independent assortment）

ni <- c(315, 108, 101, 32)
p <- c(9/16, 3/16, 3/16, 1/16)	# 独立の法則から想定される分離確率
chisq.test(ni, p = p)

想定分布のパラメータを推定する場合

工場労働者一人当たりの遭遇事故数（前回の課題）

　想定分布の第一候補として、ポアソン分布を考えると、分布パラメータλは、観測度数分布の平均 m で推定（最尤推定）される。観測度数分布の分散 v を計算する。m < v ならば、ポアソン分布では無いと考えられるので、想定分布の第二候補として、負の二項分布を考える。負の二項分布のパラメータは n、p の2つであり、負の二項分布の平均は n(1－p)/p、分散は n(1－p)/p² なので、モーメント法により、n、p が推定できる。
　すなわち、
　　　　m = n(1－p)/p,　v = n(1－p)/p² = m /p
とおける。これより、
　　　　p = m /v, 　n = mp / (1－p)
となり、観測度数分布の平均と分散から負の二項分布のパラメータが推定できる。

# ポアソン分布へのあてはめ
x <- 0:5
y <- c(447, 132, 42, 21, 3, 2)	# 5以上は5と仮定
n <- sum(y)
p <- y/n			# 観測度数分布
m <- sum(x*p); m	# 0.465224、観測度数分布平均
q <- dpois(x, m)		# ポアソン分布期待確率
v <- sum((x - m)^2*p); v	# 0.6908308、観測度数分布分散
# m < v なので、ポアソン分布では無さそう。
# 負の二項分布へのあてはめ
pb <- m/v; pb		# 0.673427
nb <- m*pb/(1 - pb); nb	# 0.9593397
r <- dnbinom(x, nb, pb)	# 負の二項分布期待確率
# 5以上は余事象の確率
q[6] <- 1 - sum(q[1:5])
r[6] <- 1 - sum(r[1:5])
nq <- n*q; nq		# ポアソン分布での期待度数
nr <- n*r; nr		# 負の二項分布での期待度数
# グラフ表示
data <- data.frame(y, nq, nr, row.names=c("0", "1", "2", "3", "4", "5+"))
barplot(t(data), beside=T, legend=c("Observed", "Poisson", "neg. bin."))
abline(h=0)

y2 <- c(y[1:4], y[5]+y[6]); q2 <- c(q[1:4], q[5]+q[6])
nq2 <- n*q2
chi2_po <- sum((y2 - nq2)^2/nq2); chi2_po	# ポアソン分布でのχ2 値、
chi2_nb <- sum((y - nr)^2/nr); chi2_nb	# 負の二項分布でのχ2 値、
1 - pchisq(chi2_po, df=3)		# 
1 - pchisq(chi2_nb, df=3)		# 

分割表の解析

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

• 勝率の同等性の検定
  　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 値となる．
  

  m <- 10; n <- 10; k <- 11
  dhyper(7:11, m, n, k)
  pu <- sum(dhyper(7:11, m, n, k)); pu		# 上側確率
  dhyper(0:4, m, n, k)
  pl <- sum(dhyper(0:4, m, n, k)); pl		# 下側確率
  pu + pl					# p 値
  

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

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 x[1,1]*x[2,2]/(x[1,2]*x[2,1]) # 通常のオッズ比 chisq.test(x) # χ2 値，p 値 fisher.test(x) # オッズ比，p 値

参考文献（古い順）

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，九天社