第一種の誤りとは

第一種の誤りとは

「帰無仮説が真のとき、これを棄却してしまう」

誤りのことです。

そして、第一種の誤りを犯す確率をαとして、有意水準または危険率と一般に呼びます。

第一種の誤りについてもう少し直感的に理解するには、条件付き確率を扱っていることを認識する必要があります。

すなわち、危険率とは


  1. 帰無仮説が真のとき(こういう条件の下だと)
  2. これを棄却してしまう確率
ということなのです。
ですので、あくまで「帰無仮説が正しいとき」を想定しなければ、危険率αのニュアンスはわからないのです。

正規分布からのサンプルの平均値の検定を題材にして、危険率αとはどのようなものなのかということをシミュレーションを通じて考えてみたいと思います。


今、平均値μ、分散σ^2/μの正規分布(標準正規分布)から、サンプリング(n)を行い、標本平均の分布について考えたいと思います。

一般に平均μ、分散σ^2の正規分布から抽出されたサンプル(n)の平均値の分布は平均値μ、分散σ^2/nの正規分布に従うとされます。
サンプルの平均値を標準化して、
を考えます。Zは標準正規分布に従います。

それでは、平均値12、分散9の正規分布の集団からn=5のサンプリングを行い、その平均値の分布について計算してみます。

一回の試行は以下のようにして行うことができます。
> rnorm(mean=12, sd=3, n=5)
[1] 15.039414  9.264087 11.036560 16.033761 11.038038

Zの値は以下のようにして計算します。
> (mean(rnorm(mean=12, sd=3, n=5) - 12)/(3/sqrt(5)))
[1] 13.86206

ところで、標準正規分布を描写すると以下のようになります。
> curve(dnorm(x), -3, 3)
ここで、上側確率、下側確率を計算してグラフに記入
> qnorm(0.025)
[1] -1.959964
> qnorm(0.975)
[1] 1.959964

abline(v=qnorm(0.025), col="red")
> abline(v=qnorm(0.975), col="red")
それでは今回想定している平均値12、分散9の正規分布の集団からn=5のサンプリングを繰り返し(300回)行い、その都度Z値を計算しています。これらのZ値のうち、上側下側合わせて5%の区間にどの程度入ってくるか調べます。

curve(dnorm(x), -3, 3)
abline(v=qnorm(0.025), col="green")
abline(v=qnorm(0.975), col="green")
for(i in 1 : 100){
abline(v=(mean(rnorm(mean=12, sd=3, n=5) - 12)/(3/sqrt(5))), col = "magenta")
}
確かに、このように繰り返し平均値をもとめていると、
「帰無仮説は正しい」
のにも関わらず、一定の割合でサンプルの平均値は5%以下の確率で生じるとされる区間に入ってきてしまいます。

従って、ある検定を行いたいと考えたときに、サンプルの統計量が偶然には起こりそうにも無い値であったとしても、帰無仮説が誤っている可能性の他に、危険率αの確率で帰無仮説が正しい場合にも同様の統計量が計算される可能性があるといえるのです。

【参考文献】
山田剛史ほか『Rによるやさしい統計学』オーム社 2008 109-119pp

発生学のおすすめサイト

The University of New South Walesのサイトがとてもまとまっています。
#UNSW Embryology
http://php.med.unsw.edu.au/embryology/index.php?title=Main_Page


発生学を動画で学べます。
#UNSW Embryology - Movie
http://php.med.unsw.edu.au/embryology/index.php?title=Movies

第一種の誤りとは ー危険率α=5%の意味ー


第一種の誤りとは

「帰無仮説が真のとき、これを棄却してしまう」

誤りのことです。

そして、第一種の誤りを犯す確率をαとして、有意水準または危険率と一般に呼びます。

第一種の誤りについてもう少し直感的に理解するには、条件付き確率を扱っていることを認識する必要があります。

すなわち、危険率とは


  1. 帰無仮説が真のとき(こういう条件の下だと)
  2. これを棄却してしまう確率
ということなのです。
ですので、あくまで「帰無仮説が正しいとき」を想定しなければ、危険率αのニュアンスはわからないのです。

正規分布からのサンプルの平均値の検定を題材にして、危険率αとはどのようなものなのかということをシミュレーションを通じて考えてみたいと思います。


今、平均値μ、分散σ^2/μの正規分布(標準正規分布)から、サンプリング(n)を行い、標本平均の分布について考えたいと思います。

一般に平均μ、分散σ^2の正規分布から抽出されたサンプル(n)の平均値の分布は平均値μ、分散σ^2/nの正規分布に従うとされます。
サンプルの平均値を標準化して、


を考えます。Zは標準正規分布に従います。

それでは、平均値12、分散9の正規分布の集団からn=5のサンプリングを行い、その平均値の分布について計算してみます。

一回の試行は以下のようにして行うことができます。
rnorm(mean=12, sd=3, n=5)
[1] 15.039414  9.264087 11.036560 16.033761 11.038038

Zの値は以下のようにして計算します。
> (mean(rnorm(mean=12, sd=3, n=5) - 12)/(3/sqrt(5)))
[1] 13.86206

ところで、標準正規分布を描写すると以下のようになります。
> curve(dnorm(x), -3, 3)
ここで、上側確率、下側確率を計算してグラフに記入
qnorm(0.025)
[1] -1.959964
qnorm(0.975)
[1] 1.959964

abline(v=qnorm(0.025), col="red")
> abline(v=qnorm(0.975), col="red")
それでは今回想定している平均値12、分散9の正規分布の集団からn=5のサンプリングを繰り返し(300回)行い、その都度Z値を計算しています。これらのZ値のうち、上側下側合わせて5%の区間にどの程度入ってくるか調べます。

curve(dnorm(x), -3, 3)
abline(v=qnorm(0.025), col="green")
abline(v=qnorm(0.975), col="green")
for(i in 1 : 100){
abline(v=(mean(rnorm(mean=12, sd=3, n=5) - 12)/(3/sqrt(5))), col = "magenta")
+ }
確かに、このように繰り返し平均値をもとめていると、
「帰無仮説は正しい」
のにも関わらず、一定の割合でサンプルの平均値は5%以下の確率で生じるとされる区間に入ってきてしまいます。

従って、ある検定を行いたいと考えたときに、サンプルの統計量が偶然には起こりそうにも無い値であったとしても、帰無仮説が誤っている可能性の他に、危険率αの確率で帰無仮説が正しい場合にも同様の統計量が計算される可能性があるといえるのです。

【参考文献】
山田剛史ほか『Rによるやさしい統計学』オーム社 2008 109-119pp

統計学的検定の根本的な概念

統計学的検定の根本的な概念は

母集団から抽出したデータから計算した統計量が、次のいずれに属するかを検討することだと思います。


  1. 偶然誤差の範囲にある(それくらいのばらつきは生じても不思議ではない)
  2. 偶然誤差の範囲を逸脱している(きっと想定した母集団モデルから逸脱させる何か、すなわち要因が働いているはずである)
以上のように考えると、統計学的検定の意味が少しずつ見えてくる気がします。
成書を読むと、1,2と同様の内容が必ず"さらっと"書かれていることに気がつきます。

Rを用いてdot plotを描写

Rを用いてdot plotを描写する方法を勉強しました。
dot plotには平均値と標準偏差を同時に記述することが多いので、今回もその慣例に従いました。

参考にさせていたいたのは下のリンクです。

http://minato.sip21c.org/swtips/MedicalStat-Rev.pdf

まずは、以下のデータをEXELで用意して、tab区切り形式で保存します。

あとはRを開いて以下のコマンドを実行すればよいです。

$ R


#データの読み込み
> dat <- read.table("body_height.txt", header=T)
dat <- as.data.frame(dat)

> dat
   Male Sex
1   170   1
2   175   1
3   181   1
4   176   1
5   167   1
6   180   1
7   156   2
8   150   2
9   153   2
10  148   2
11  157   2
12  149   2

#Sex列を要因(factor)にする
dat$Sex <- as.factor(dat$Sex)
levels(dat$Sex) <- c('M','F')

dat
   Male Sex
1   170   M
2   175   M
3   181   M
4   176   M
5   167   M
6   180   M
7   156   F
8   150   F
9   153   F
10  148   F
11  157   F
12  149   F

#各列名単独で列の要素にアクセスできるようにする
attach(dat)


#平均値、標準偏差を計算
mean <- tapply(Height,Sex,mean)
sd <- tapply(Height,Sex,sd)
is <- c(1,2)+0.15


#ドットプロットを描写
stripchart(dat$Height~Sex,method="jitter",vert=T,ylab="body height(cm)") 
points(is,mean,pch=18)
arrows(is,mean-sd,is,mean+sd,code=3,angle=90,length=.1)


勝率と二項検定


サッカーのあるチームAがあるチームに対して、年間を通じて、Bチームと30試合を行い25勝5敗だったとします。
このデータから、チームAはチームBより強いと言えるでしょうか。
言い換えると、チームAのチームBに対する真の勝率は5割以上なのでしょうか。それとも、25勝5敗は誤差の範囲で勝ち越しているだけなのでしょうか。
この問題を解くには、二項検定という手法が用いられます。コインの裏表の出る確率の検定のときにも用いられる手法です。

#勝率が5割だったとした時の、勝利数の確率分布を描写する
> win <- 0:30
> plot(dbinom(win, 30, 0.5), type="h", xlab-"number of win")


#累積分布関数を描写する
> plot(pbinom(win, 30, 0.5), type="h", xlab="number of win")
#95%の水準に赤線を引く
abline(h=0.95, col="red")
#5%の危険率で二項検定を行う
#方法1
> 1 - pbinom(25, 30, 0.5)
[1] 2.973806e-05

#方法2
> pbinom(25, 30, 0.5, lower.tail=FALSE)
[1] 2.973806e-05

> pbinom(25, 30, 0.5, lower.tail=FALSE) < 0.05
[1] TRUE
> pbinom(25, 30, 0.5, lower.tail=FALSE) < 0.01
[1] TRUE

> binom.test(25,30, 0.5)

 Exact binomial test

data: 25 and 30 
number of successes = 25, number of trials = 30, p-value = 0.0003249
alternative hypothesis: true probability of success is not equal to 0.5 
95 percent confidence interval:
 0.6527883 0.9435783 
sample estimates:
probability of success 
             0.8333333 

以上の計算結果から、チームAのチームBに対する勝率は有意に50%よりも大きいことが言えます。また、その信頼区間は0.65以上0.95以下です。

今回行った計算方法は、帰無仮説であるp=0.5を他の値に変更することによって、二項分布に従う様々な現象(打率、くじ引き etc)を説明することができます。

【参考文献】
山田剛史ほか『Rによるやさしい統計学』オーム社 2008 第12章

重回帰分析におけるモデルの検討


重回帰分析の手法について分子生物学的な事例をもとにして考えてみます。
ある時刻の遺伝子Mの発現量exp_mが、一時間前の遺伝子n及び遺伝子lの発現量(それぞれexp_n, exp_l)によりどのよう説明されるかというモデルを立てます。
具体的には、




というモデルを想定します。
aを切片、b_1、b_2を偏回帰係数、eを残差と呼びます。
それでは、Rを用いて重回帰分析を行ってみます。

#データの入力
> exp_m <- c(45, 46, 47, 57, 78, 89, 99, 89, 101, 108)
exp_n <- c(102, 106, 109, 120, 124, 130, 146, 178, 190, 196)
exp_l <- c(56, 57, 58, 50, 52, 61, 43, 63, 64, 60)

#exp_mを目的変数、exp_n, exp_lを説明変数として重回帰分析を行う
> summary(lm(exp_m~exp_n+exp_l))

Call:
lm(formula = exp_m ~ exp_n + exp_l)

Residuals:
    Min 1Q Median 3Q Max 
-11.530 -5.954 -3.900 3.842 24.637 

Coefficients:
                       Estimate   Std. Error     t value    Pr(>|t|)    
(Intercept)        36.3470      34.7881      1.045    0.330842    
exp_n                0.6858         0.1221     5.618     0.000801 ***
exp_l                -1.0022         0.6710    -1.494     0.178918    
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 11.91 on 7 degrees of freedom
Multiple R-squared: 0.822,Adjusted R-squared: 0.7712 
F-statistic: 16.16 on 2 and 7 DF, p-value: 0.002379 

exp_nのp値が有意となっているということは、冒頭の式のb_1が有意に0よりも大きい値であるということができます。また、exp_nを用いれば、各値の全体平均からのばらつきを有意に説明するとも言い換えることができます(ANOVAの考え方と同じです)。一方で、exp_lのp値は有意ではありませんでした。
R-squaredとは、決定係数と呼ばれ、「説明変数を用いてどの程度、目的変数の分散を説明できるか」という説明力を意味しています。のR-squaredは0.8と1に近いため、まずまずのモデルと言えるでしょう。
次に、exp_lの寄与は有意ではなかったため、exp_lをはずして、exp_mを説明することを試みます。



#exp_mを目的変数、exp_n, exp_lを説明変数として重回帰分析を行う

> summary(lm(exp_m~exp_n))

Call:
lm(formula = exp_m ~ exp_n)

Residuals:
    Min 1Q Median 3Q Max 
-10.069 -8.693 -6.009 8.438 19.493 

Coefficients:
                      Estimate  Std. Error    t value Pr(>|t|)    
(Intercept)      -9.7452   17.2506     -0.565   0.587616    
exp_n              0.6113   0.1197        5.107   0.000921 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 12.8 on 8 degrees of freedom
Multiple R-squared: 0.7653,Adjusted R-squared: 0.736 
F-statistic: 26.08 on 1 and 8 DF, p-value: 0.0009215 

#散布図に回帰直線を重ねて描写
> plot(exp_n, exp_m)
> abline(lm(exp_m~exp_n)

R-squredは0.06程度低下しましたが、依然として良いモデルと言えます。
実際の研究ではもっと多くの変数の組み合わせを検討してモデル探索を行い、決定係数をもとに独立変数の数を絞ります。

【参考文献】
山田剛史ほか『Rによるやさしい統計学』オーム社 2008 191 - 198pp

LaTeXItを使って簡単に数式を作成する!!!

ブログやプレゼンテーションスライドを作っているときに、一行だけ数式が欲しいということがあります。texを使うのはちょっと面倒だなあというときにおすすめなのが、LaTexItというソフトです。MacTexをインストールすると自動的にインストールされます。

MacTex
http://www.tug.org/mactex/
MacTeX.pkg 


下の画面の中で、白い部分に数式を入力して、LaTex it!をクリックすれば、数式が出力されます。「File」->「Export Image」で保存できます。
これが、LaTeXItで作成した数式です。png形式で出力しています。



対応のあるt検定とANOVA


今回は、「対応のあるt検定」について考えます。対応のあるt 検定を意識したデータ構造は以下のようになります。下の図は、ある薬剤を投与前と投与後の体重の変化を示しています。検証したい命題は「ある薬剤を投与した前後で体重は変化するか」です。

対応のないt検定と比較する目的で、下の表を提示します。各群のデータは順番は入れ替えてあるものの、上の表と同じです。薬剤投与を受けた人と薬剤投与を受けなかった人は別の人を想定します。つまり、上の表の2倍の人数が下の表では動員されていることとなります。



これから、Rを用いて対応のある場合と対応のない場合で、薬剤投与による効果
> pre <- c(67, 70, 78, 65, 75, 77, 70, 76, 79, 83)
> post <- c(66, 68, 78, 63, 74, 76, 68, 75, 78, 82)

> negative <- c(78,67,75,65,70,77,70,83,79,76)
> positive <- c(68,66,78,75,82,78,76,63,68,74)

#ドットプロット
> source("http://aoki2.si.gunma-u.ac.jp/R/src/dot_plot.R", encoding="euc-jp")
#対応あり
> annot <- factor(c(rep("pre", 10), rep("post", 10), rep("negative", 10), rep("positive", 10)))
> dat <- c(pre, post, positive, negative)
> dot.plot(annot, dat)
#対応のあるt検定
> t.test(pre, post, paired=T)

 Paired t-test

data: pre and post 
t = 6, df = 9, p-value = 0.0002025
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.7475686 1.6524314 
sample estimates:
mean of the differences 
                    1.2 

#対応のないt検定
> t.test(negative, positive, paired=F)

 Welch Two Sample t-test

data: negative and positive 
t = 0.4494, df = 17.91, p-value = 0.6585
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -4.411489 6.811489 
sample estimates:
mean of x mean of y 
     74.0 72.8 

驚くことに対応のあるt検定では有意差が検出されたのに、対応のないt検定では有意差が検出されませんでした。
この理由は、個人によるばらつきを対応のあるt検定によって制御できたことにあります。対応のないt検定の場合、個人差によるばらつきが残差に入ってしまうため、残差の分散が大きくなり、その結果薬剤の投与という要因による効果を検出するには至らなかったのです。これはまさにANOVAの考え方そのものです。というか、一般線形モデルの本質的な部分です。

対応のある一元配置分散分析を行うことで、薬剤と人という二つの要因による効果の大きさについて検証することができます。

> dat
 [1] 67 70 78 65 75 77 70 76 79 83 66 68 78 63 74 76 68 75 78 82


> drug <- factor(c(rep("negative", 10), rep("positive",10)))
> drug
 [1] negative negative negative negative negative negative negative negative
 [9] negative negative positive positive positive positive positive positive
[17] positive positive positive positive
Levels: negative positive

> person <- factor(rep(c("Hiroshi", "Takeshi", "Yamato", "Keita", "Tetsuya", "Shinta", "Yohei", "Daichi", "Yusuke", "Daisuke"),2))
> person 
[1] Hiroshi Takeshi Yamato Keita Tetsuya Shinta Yohei Daichi Yusuke 
[10] Daisuke Hiroshi Takeshi Yamato Keita Tetsuya Shinta Yohei Daichi 
[19] Yusuke Daisuke
10 Levels: Daichi Daisuke Hiroshi Keita Shinta Takeshi Tetsuya Yamato ... Yusuke

#drugとpersonの二つの要因を考慮してANOVAを行う。ただし、drugとpersonの間に交互作用はないものとする。
> summary(aov(dat~drug+person)) 
               Df  Sum Sq     Mean   Sq      F value    Pr(>F)    
drug         1           7.2              7.20           36.0    0.000202 ***
person      9        639.8            71.09         355.4    2.16e-10 ***
Residuals 9             1.8             0.20                     
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


#drugだけを要因として考えてANOVA
> summary(aov(dat~drug))
               Df  Sum Sq     Mean   Sq      F value    Pr(>F)   
drug         1           7.2              7.20        0.202     0.658
Residuals  18      641.6            35.64     


drugとpersonを考慮してANOVAを行うと残差(Residuals)の平方和(Mean Sq)が小さくなり、その結果F valueが大きくなることがわかります。F valueを求めるには、

F value = (drugのMean Sq) / (ResidualsのMean Sq)

です。
このように、対応のあるt検定はanovaにおいて、要因を追加することにより系統誤差を制御していることになるのです。その結果、説明できない誤差(すなわち残差)が小さくなるため、F値は大きくなり、その結果、小さなp値が出やすくなります。

以上のことから、実験結果に影響を及ぼす要因はできるだけ見いだしておいて、制御できるものは実験デザインによって制御することが大事であると言えます。

【参考文献】
山田剛史ほか『Rによるやさしい統計学』オーム社 2008 150-156, 175-179pp

フィッシャーの三原則

統計学の父であるフィシャー:Ronald A. Fisherは、農業試験の適切な実施のための方法論として実験計画法という学問分野を確立しました。
Fisherは局所管理、反復、無作為化の3つが大切であると述べています。これはフィッシャーの3原則と呼ばれ、農学以外の幅広い分野において利用されています。
以下に、肥料の種類による作物の生育具合の評価を目的とした農場試験を例に取って説明します。

1. 局所管理(小分け)の原則
実験計画の最も原始的な形態は以下のようなものです。複数水準の間の効果を比較したいのですが、同じ場所で異なる年度で調査を行おうとしています。これでは、年度ごとの気候の差が系統誤差として入ってきます。
そこで以下のように局所管理(小分け)すると、年度ごとの気候の差を調整することができます。


2. 繰り返し(反復)の原則
小分けをしただけではまだ問題が残ります。ある水準で生育がとても良かったときに、それが果たして偶然誤差によるものなのか、水準の違い(肥料の違い)によるものなのかがはっきりしません。この問題は水準ごとに実験を反復することで、解決することができます。同じ水準内で比較することにより、偶然誤差による差を調整することができるためです。

3. 無作為化の原則
2のように水準ごとに繰り返しを作ることで、偶然誤差と系統誤差(水準によるものとそうでないものを含む)に分解することが可能となりました。最後は、実験の割り付けを無作為化することが必要となります。2.のデザインでは、水はけは日当りなど、目的要因以外の要因が、一定の方向で偏りを持って実験結果に影響を与えている可能性が高くなるからです。つまり、水準以外の系統誤差が調整できていないのです。そこで、下の図のように無作為に実験を割り付けることにより、それらの系統誤差を調整することが可能となります。

以上、フィッシャーの3原則でした。
基本的な発想は、「目的となる要因以外からの影響が極力排除すること」です。肥料以外の要因は一定になるように努力するとともに、もしそれを制御しきれないのならば、せめて肥料以外の要因が実験結果に偏って現れないようにデータ収集することが大事なのです。

この発想は臨床試験や生物学実験においても大切な発想であると思います。


【参考文献】
栗原 伸一『入門統計学』Ohmsha, 2011, 88-89, 130-138pp