ヒストグラムの階級を細かくして正規分布の密度関数を近似的に描写する

#標準正規分布から乱数を取得
dat <- rnorm(1000000)
#レンジを調べる
range(dat)
png("120812_hist.png")
par(mfrow=c(2,2))
hist(dat, breaks=seq(-10,10,10))
hist(dat, breaks=seq(-10,10,1))
hist(dat, breaks=seq(-10,10,0.1))
hist(dat, breaks=seq(-10,10,0.01))
dev.off()

dat <- rnorm(1000000)

png("120812_hist_2.png")

par(mfrow=c(1,2))

hist(dat,breaks=seq(-6,6,0.1))

curve(length(dat)*dnorm(x))

dev.off()

カイ二乗検定で独立性の検定を行う

いま、20人の大学生がいます。
それぞれに数学と物理学が好きか嫌いかを二択で答えてもらうとします。
結果が以下のようになったときに、数学と物理の好き嫌いの傾向は独立かどうかを検定します。

$ R

> 観測度数 <- c(10, 2 , 4, 4) 

> 期待度数 <- c(8.4, 3.6, 5.6, 2.4)

> カイ二乗要素 <- (観測度数-期待度数)^2/期待度数

> カイ二乗要素
[1] 0.3047619 0.7111111 0.4571429 1.0666667
> カイ二乗値 <- sum(カイ二乗要素)
> カイ二乗値
[1] 2.539683


> pchisq(カイ二乗値, 1, lower.tail=FALSE)
[1] 0.1110171

p > 0.05であるため、このアンケートの結果は偶然でも起こりえる程度のものであるという結論となります。

ちなみに、df= 1のカイ二乗分布は以下のようになります。

png("120807_chi_squere.png")

curve(dchisq(x,1), 0, 30)

abline(v=qchisq(0.05,1,lower.tail=FALSE), col = "red")

abline(v=カイ二乗値)

dev.off()

Mac OS XにEMBOSSをインストール

UNIX環境における配列解析の定番ソフト「EMBOSS」をビルドしてみます。

#ソースコードのダウンロード
ftp ftp://emboss.open-bio.org
> cd pub
> cd EMBOSS
> get EMBOSS-6.5.7.tar.gz
> quit

#解凍
$ tar xzvf EMBOSS-6.5.7.tar.gz
$ cd EMBOSS-6.5.7

#configureスクリプトを実行
$ ./configure
#コンパイル
$ make
#インストール
$ sudo make install

#試しに、embleからxlrhodopの塩基配列(Xenopus laevis rhodopsin mRNA, complete cds)を取得する。

$ seqret
Read and write (return) sequences
Input (gapped) sequence(s): embl:xlrhodop
output sequence(s) [l07770.fasta]:
kappa-no-MacBook:120806 kappa$ ls
l07770.fasta

>L07770 L07770.1 Xenopus laevis rhodopsin mRNA, complete cds.
ggtagaacagcttcagttgggatcacaggcttctagggatcctttgggcaaaaaagaaac
acagaaggcattctttctatacaagaaaggactttatagagctgctaccatgaacggaac
agaaggtccaaatttttatgtccccatgtccaacaaaactggggtggtacgaagcccatt
cgattaccctcagtattacttagcagagccatggcaatattcagcactggctgcttacat
gttcctgctcatcctgcttgggttaccaatcaacttcatgaccttgtttgttaccatcca
gcacaagaaactcagaacacccctaaactacatcctgctgaacctggtatttgccaatca
cttcatggtcctgtgtgggttcacggtgacaatgtacacctcaatgcacggctacttcat
ctttggccaaactggttgctacattgaaggcttctttgctacacttggtggtgaagtggc
cctctggtcactggtagtattggccgttgaaagatatatggtggtctgcaagcccatggc
caacttccgattcggggagaaccatgctattatgggtgtagccttcacatggatcatggc
tttgtcttgtgctgctcctcctctcttcggatggtccagatacatcccagagggaatgca
atgctcatgcggagtagactactacacactgaagcctgaggtcaacaatgaatcctttgt
tatctacatgttcattgtccacttcaccattcccctgattgtcatcttcttctgctatgg
tcgcctgctctgcactgtcaaagaggctgcagcccagcaacaggaatctgctaccaccca
gaaggctgagaaagaggtcaccagaatggttgttatcatggtcgttttcttcctgatctg
ttgggtgccctatgcctatgtggcattctacatcttcacccaccagggctctaactttgg
cccagtcttcatgaccgtcccagctttctttgccaagagctctgctatctacaatcctgt
catctacattgtcttgaacaaacagttccgtaactgcttgatcaccaccctgtgctgtgg
aaagaatccattcggtgatgaagatggctcctctgcagccacctccaagacagaagcttc
ttctgtctcttccagccaggtgtctcctgcataagagcttcaccagggctgtctcagggt
ccgctgcctcacacaattcccatcacttaagccctgtctacttgttgcgaaggcaaagaa
ttccacagttttaatatttacccccattctgcccaaccttggacactgtaagagctgacc
ccattactgctgggaaggcccaagctttgttgcattctgatgtgatcctttcagcagaaa
atgggtggattcaatgaatttcaccaaggctgtacataacaataacattagtctgaaggc
acctcccacccagagaatgcaacacttatttatctctgtcttttcttgacatattgatgc
tgcttctattcatggtcactaacaaaaagtcccattttacaatgcaactgaaagtaatgt
atttttgtaatataataacatatttcatgcaatctcctctgcttattggcaaggtctgat
atagtgaggatagacagccagaccccttgcattaaaatcctgtattaaaaatttctttgc
aagt

サイコロの出る目の度数分布

各数字が1/6の確率で得られる理想的なサイコロがあるとします。
ここで、このサイコロを繰り返し振って、出た目の度数分布を作成します。


必要な関数(ceiling, runif)の仕様を勉強。

$ R
> ? ceiling
Description:
     ‘ceiling’ takes a single numeric argument ‘x’ and returns a
     numeric vector containing the smallest integers not less than the
     corresponding elements of ‘x’.

与えたnumeric vector(要素は小数でも可)の各要素の小数点を切り捨てた整数値を返します。

>?runif

Description:

     These functions provide information about the uniform distribution

     on the interval from ‘min’ to ‘max’.

与えられた最小値から最大値の間で、一様分布に従う値を返すとのこと。

$ vim command.in

png("120805_dice.png")

par(mfrow=c(2,2))

dice1 <- ceiling(runif(n=6, min = 0, max = 6)) 

dice2 <- ceiling(runif(n=60, min = 0, max = 6)) 

dice3 <- ceiling(runif(n=600, min = 0, max = 6)) 

dice4 <- ceiling(runif(n=6000, min = 0, max = 6)) 

hist(dice1,main="n=6",breaks = seq(0,6,1))    

hist(dice2,main="n=60",breaks = seq(0,6,1))    

hist(dice3,main="n=600",breaks = seq(0,6,1))    

hist(dice4,main="n=6000",breaks = seq(0,6,1))   

dev.off()

$ R

> source("commnad.in")

n数を多くすれば母集団の分布を反映したものが得られます。

検出力と例数設計

Rには、t検定の検出力を計算するための関数が用意されています。

それが、power.t.test()です。

今回はこの関数を用いて、検出力と例数設計について検証してみたいと思います。

先日、身長平均身長175cm(±5cm）と他の集団との身長の平均値に関しての比較を扱いました。その時、有意差を検出するためのn数について直感的な検討を行いました。

今回は、パワーテストなので、より定量的に議論したいと思います。

以下のスクリプトは、1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差　2)想定される平均値の差、標準偏差、α危険率、検出力をある値に想定したときに要請されるn数　3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値　について着目しています。

今回はベースに「平均身長175cm(±5cm)の集団と平均身長180cm(±5cm)の集団の平均値の比較を意識していて、サンプル数は１０、p値は0.05を想定している」状況があります。


$ vim command.in


#以下の条件を満たすのに必要なn数を求める。
n_null = power.t.test(delta=5,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件で有意差が検出可能な母平均の差を求める。
delta_null = power.t.test(n=10,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件での検出力を求める。
power_null = power.t.test(delta=5,n=10,sd=5,sig.level=0.05,type=c("two.sample"),alternative="two.sided")


$ R
> source("command.in")

#特定の項目を取り出す方法を調べる。

> names(delta_null)
[1] "n"           "delta"       "sd"          "sig.level"   "power"    
[6] "alternative" "note"        "method"    

#以下、1)-3)のケースの計算結果を表示


> delta_null$delta
[1] 6.624728
> n_null$n
[1] 16.71477
> power_null$power
[1] 0.5619846

1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差


> delta_null$delta
[1] 6.624728

サンプル数は１０、p値は0.05、検出力0.8、標準偏差5を想定しているで検出可能な平均値の差は6.6cmであるということです。

2)標準偏差、α危険率、検出力をある値に想定したときに要請されるn数


> n_null$n
[1] 16.71477

平均値の差が5、p値は0.05、検出力0.8、標準偏差5の条件を満足させるために必要なn数は16.7である。

#実験計画を立てる上で特に重要である。

3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値


> power_null$power
[1] 0.5619846

n数が10、平均値の差が5、p値は0.05、標準偏差5の条件の時の検出力は0.56である。

2種類の過誤を明確に意識し、さらには母集団の平均値の差を実験によりある程度想定しながら、サンプル数を決定したいと思いました。

・第一種過誤（α過誤、偽陽性）は、帰無仮説が実際には真であるのに棄却してしまう過誤

・第二種過誤（β過誤、偽陰性）は、対立仮説が実際には真であるのに帰無仮説を採用してしまう過誤

最後に実際に実験を行うことを意識して

「想定する母集団（標準偏差５）の差を1から0.5cmずつ10cmまで増加させていったときに、p値は0.05、検出力0.8の条件を満足させるために必要なn数を２次元プロットに描写してみます。

$ vim command_2.in

n <- c(1:100)

for(i in 1:100)

{

n[i] <- (power.t.test(delta=0.1*i,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided"))$n

}

delta <- c(1:100)*0.1

png("120801_delta_n_plot.png")

plot(delta,n)

dev.off()

$ fg

> source("command_2.in")

このシミュレーション結果は非常に深い意味があります。

平均値の差が小さい時(2ぐらいまで）は、非常に多くのサンプル数が必要となるのに対して、それよりも大きい平均値の差であるときには、サンプル数がほぼ一定（5程度）で済んでしまうというわけです。

deltaの値を0.11から10までの間に関してプロットを書かせるとよりわかりやすいです。

> plot(delta[11:100],n[11:100])

実験を計画するときはこのような事実を念頭に置いてサンプル数を設定したいものです。