![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
Rを使った統計解析 †アナウンス： 「千葉大学留学フェア」開催中（6月9日(月)〜13日(金)） †
第9回授業の獲得目標：
|
| データをコピー・ペーストや手入力で、オブジェクトに代入したいときには |
| scan() という関数を使う |
> h=scan() 1:
画面に表示された 1: の後に、上のページから身長データをコピー・ペーストして、最後にenter キーを押す。次のように表示される。
1: 1775 1710 .. ..＜省略＞.. 111:
同様にして、体重データも w というオブジェクトにを入れる。次のように表示される。
1: 79.8 58.0 .. ..＜省略＞.. 111:
確認のために
> h > w
と入力。データが一覧表示される。
計測データが得られたら、まずは、グラフにしてみるよう。グラフとして視覚化することで、データの持つ性質が直観的に理解できる。こういうとき、Rは、いろんなグラフを一瞬にして描画してくれるので、非常に便利。
まずはヒストグラム（柱状グラフ）を使っみる。Rを使ってヒストグラムを書く場合は
> hist(h)キーボードからたった７文字打ち込むだけで、身長の度数分布がグラフで表示される。
このグラフを見ただけで、標本データでは身長170cmぐらいの学生が最も多く、 データは釣り鐘型の分布をしているということが分かる。同様のことを体重についてもやってみよう。
> hist(w)
次に、この標本データの平均とか最大値とか、最小値とか、データの集まりがもっている基本的な性質を表示させて見よう。次のように入力。
> summary(h) Min. 1st Qu. Median Mean 3rd Qu. Max. 1578 1671 1700 1706 1734 1839これで、1997年当時の調査した大学生110人の身長は、平均170.6cm（最大:183.9cm, 最小: 157.8cm）ということが分かった。1st Qu.と3rd Qu.というのは、全体の４分の１の人数(つまり25%）がどの高さまでの間に含まれているかということを表している。つまり、167.1cm以下の人が25%、173.4cm以上の人が25%、その間の人が全体の50%ということ。
> summary(w)
身長も体重も釣り鐘型の分布になっていた。では、この２つの数値の間に何か、関係があるだろうか？これもまた、グラフで表示させてみる。こういう場合、縦軸に体重、横軸に身長をとって、一人一人の持つ身長と体重の散布図を作る。Rの場合、次のように入力するだけ（注: 得られた身長・体重データは、得られた順に並んでいる。つまり、同じ位置（最初からの順番）にある身長・体重は１人の人から得られたもの）。
> plot(h,w)身長と体重の関係は右上がりの直線関係にありそうだ。
>cor(h,w) [1] 0.5872979 #値が正なら正の相関
今回、110人分のデータを使って、1997年当時の大学生男子の身長と体重について、
グラフによる視覚化 要約統計量の表示 散布図の作成という３つのことを行った。この標本が母集団の性質を表しているとするならば、1997年当時の大学生について、次のような推定ができる。
1997年当時の男子大学生は、身長の平均が170.6cmで、データは釣り鐘型に分布している 体重は平均が59.53kgで、身長と同様にデータは釣り鐘型に分布している 身長と体重の間には、正の相関があるこれで、1997年当時の男子大学生という母集団の身長・体重について、統計的推定を行うことができた！
なお、今回は分散や標準偏差の話しは省略した。Rを使った統計解析の基本（データの図示、代表値の計算、データ間の相関）などについてはここでは説明しきれないので、下のリンクにある中澤さんによる解説などを読んで、勉強してみて欲しい。。
エクセルデータをそのまま読み込むには、http://cse.naro.affrc.go.jp/takezawa/r-tips/r/41.html を参照。
![[smile]](image/face/smile.png)
†この授業で対象とする統計解析の２つめは、
"ある対象や実験で得られた数値データを用いて、その対象についてどういう判断を下せば良いかを論じること"
それでは、次の設問について、皆さんの判断を考えて下さい。
・A君の方がB君よりもじゃんけんに強いと言っていい ・A君の方がB君よりもじゃんけんに強いとは言えない ・どちらとも判断できない次に、moodleページの演習の指示に従って、グループワークを行おう。
皆さんからは、最初は様々な意見が出ていると思うが、こういう問題はどのように考えるのが良いのだろうか？問題は、
判断基準
だ。科学的に説明するためには、論理的な「判断基準」が無ければならない。
そして、ある「判断基準」に基づいて判断を下すのが、これから学ぶ「仮説検定」だ。
まず、この授業では、正確さは二の次にして、仮説検定の流れを直観的に理解してもらえるような説明を試みる。
（この説明のどこが、後で出てくる「帰無仮説」、「有意水準」、「帰無仮説の棄却」にあたるのかは、後で考えてみて欲しい。なんとなく理解できたら、どれでもいいから統計の教科書を読んでみて、正確に理解することをオススメする。）
A君とB君が20回ジャンケンをしたら、B君は5回しか勝てなかった。ジャンケンの勝ち負けって、半々の確率で決まるはずだけど、こんなに勝ち負けが偏ることって、たまたま偶然で起こったことなのかなー？それとも、A君が何かズルをしていたと考えた方がいいのかなー？。。
ジャンケンの勝敗が、ズルとかじゃなくて、本当に半々の確率で決まっているのなら、...
「20回中5回しか勝てない」っていうのは、偶々（たまたま）起こったことなんだろうね。でも、その「偶々（たまたま）」というのは、どの程度の頻度で起こることなんだろうか？
よし。「ジャンケンの勝敗が完全に半々の確率で決まるというのを前提として、20回ジャンケンをする実験」を何度も何度も繰り返したとき、片方が何回勝つかをグラフで表してみよう... もしも、「20回中5回しか勝てない」っていうできごとが、統計的には滅多に起こらないことだとしたら、そもそも「A君とB君のジャンケンの勝敗が半々の確率で決まっている」という前提が間違っていると言って良いんじゃないかな？
そうすると、「滅多に起こらない」をどうやって判断するか、基準が必要だけど、ここでは、5%を基準にして、実験を繰り返したときに全体の5%でしか起きないことは、滅多に起こらないと判断することにしよう。
グラフを作ってみると、なんと！! 「20回中5回しか勝てない」状況は、全体の5%未満しか起きていない...。つまり、5%を基準にすると、 「20回中5回しか勝てない」ことは「滅多に起こらない」んだ！
「ジャンケンの勝敗は完全に半々の確率で決まる」という前提の下では、「B君が20回中5回しか勝てない」っていうのは滅多に起こらないというわけだから、
「ジャンケンは完全に半々の確率で決まる」っていう前提が間違っていたと考える方がいいんだろうね。
うん、きっと、どんな方法かは分からないけれど、A君はB君にジャンケンで勝てる方法を見つけてるんだ
でも、待てよ。「滅多に起こらないこと」の判断基準が5%って大きすぎない？
5%っていうと、、、100回中5回。。。 20回中に１回だ！
っていうことは、もしも、「ジャンケンは完全に半々の確率で決まる」にしても、20回に1回は、間違った判断を下してしまう可能性があるってことだね。
上の例では、20回のジャンケンの結果を使って統計的検定の直観的な説明を試みた。こんな単純な例であっても、実際の生物学の研究で使う、統計的検定と、推論の進め方は同じだ。（注意： 下の例は、今回のジャンケンの例にあわせるために、かなり説明を端折っている。ここで直観的に理解したら、上記統計の教科書を見てみよう）
まず、統計的検定で、知りたいのは、
「A君とB君が20回ジャンケンをしたとき、5回しか勝てない」という状況は、 「ジャンケンの勝ち負けが、本当に半々の確率で決まる」という前提のもとで、 単なる偶然で起こるようなことか、 あるいは、偶然では滅多に起こらないことか」
ということ。
このように、得られたデータが単なる偶然で生じることかどうかを知るときに、統計的検定という作業が必要になる。
統計的検定では、次のような論証う。
A君とB君のジャンケンの勝敗は、半々の確率でランダムに決まる （つまり「勝率に差は無い」ということ）この仮説のことを 帰無仮説と言いう。（統計の教科書を見ると、帰無仮説には「違いがない」とか「差がない」という仮説を使うと書かれている。今回の場合「『ジャンケンの勝敗は、半々の確率で決まる』というだけで、『差が無い』とか『違いが無い』とかは言って無いじゃないか？」と思う人がいるかもしれない。しかし、「ジャンケンの勝敗は、半々の確率でランダムに決まる」ということは、「ジャンケンの勝率には差が無い」ということを意味している。字面だけで判断するとこんがらがる場合があるので、要注意。）
A君とB君のジャンケンの勝敗は、半々の確率ではきまらない（A君はB君よりもジャンケンに強い）というもので、対立仮説は帰無仮説を否定したものになっている。
p値の計算をするときに、「20回中勝ちが5回」の場合だけを計算するだけではなく 20回中勝ちが5回以下（5, 4, 3, 2, 1, 0）の場合の確率も計算しなくてはならないということということ。この後の議論に出てくるように、起きた事象がどれだけ希にしか起きないことかを判断して、それが滅多におこらないようなことなら、前提となった帰無仮説が間違っていると判断する。つまり、「勝ちが5回以下」のそれぞれのできごとは「勝ちが5回」という出来事よりも起こりにくいことなので、そう言う場合にも、帰無仮説は棄却される。そのため、p値の計算は「5回以下」の場合の確率が必要になるの。この確率のことを累積確率と呼ぶ。
今回の解説では、片側検定、両側検定の話しは意図的に省略してある
![[worried]](image/face/worried.png)
journal.pone.0091836.g001&representation=PNG_L
「A君とB君のジャンケンの勝敗は半々の確率で決まる」（帰無仮説）
（この場合、5回しか勝てないというデータの偏りは単なる偶然）
が棄却されたので、
「A君とB君のジャンケンの勝敗は、半々の確率ではきまらない」
（対立仮説）（この場合、データの偏りは、単なる偶然では無い）
「A君とB君のジャンケンの勝敗は、半々の確率で決まる」（帰無仮説）が正しいという仮定の下で、「20回中5回しか勝てない」という状況が生じる確率(p値）が、「滅多に起こらないこと」の判断基準（有意水準）よりも大きいので、「A君とB君のジャンケンの勝敗が半々の確率で決まるときに、20回中5回しか勝てないという状況は、生じても別に不思議ではない」ということ。このとき、
「A君とB君のジャンケンの勝敗は、半々の確率で決まる」という帰無仮説は、5%の有意水準で棄却 されないという。得られたデータからは帰無仮説を棄却できなかったので、この場合、「A君とB君のジャンケンの勝率には、有意な違いは無かった」と言う。ここで注意しなければならないのは、「有意な違いが無い」からといって、「A君とB君のジャンケンの勝敗は、半々の確率で決まる」と結論づけてはいけないこと。「有意な違いが無い」のは、「A君とB君のジャンケンの勝敗は、半々の確率で決まっていてもおかしくはない」ということを意味しているだけで、単なる偶然の結果かもしれないし、データ数が少ないせいかもしれないし、いろんな場合があり得る。
では、実際に先ほどの例をRを使って解析してみよう。この検定は二項検定と呼ばれている。また、ここでは有意水準を5%としよう。
まず、帰無仮説は
A君とB君のジャンケンの勝敗は、半々の確率で決まる（つまり「勝率に差は無い」ということ）
ということ。この帰無仮説が正しいと仮定すると、 得られたデータの偏りは、サンプリングの際の偶然の偏りであると考えられる。ジャンケンの勝ちという、確率50%で生じる事象が、20回中5回得られたわけだが、検定の際には、これより希な条件も全て含めて確率を計算する。つまり、20回中、5回以下と、さらには勝ちの数が15以上の事象が生じる確率を計算。 手で計算してもいいが、今回はすぐにRを使ってみよう。上の検定を行うにはbinom.test()関数を使う。
> binom.test(x,n)
と入力するだけで、実際のx個のサンプルの中で、データ値がn回実現したという状況で、そのデータ値は（もう一つの対になるデータ値に対して）母集団中で1:1で出現するという帰無仮説をテストできる。
このテストの詳しい解説は、上記、粕谷さんの統計の本に載っている。
上の例だと、20回中5回が勝ちで、帰無仮説は「A君とB君のジャンケンの勝敗は、半々の確率で決まる」だったので、
> binom.test(5,20)
Exact binomial test
data: 5 and 20
number of successes = 5, number of trials = 20, p-value = 0.04139
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.5089541 0.9134285
sample estimates:
probability of success
0.25
p値(p-value)は0.04139なので、有意水準が5%よりも小さい。そこで、帰無仮説「A君とB君のジャンケンの勝敗は、半々の確率で決まる」は棄却される。（もちろん、有意水準を0.001にすると棄却されない）
つまり、「A君とB君のジャンケンの勝敗は、有意水準5%で、半々の確率で決まるとは言えない」ということができる（有意水準0.1%の場合は、「２人の勝率に差があるかどうかはわからない」）。「有意水準5%では、A君とB君のジャンケンの勝敗は、半々の確率で決まるという帰無仮説は棄却された」という言い方もある。
どうやら、やっぱりA君の方がB君よりも強いようだ。もしかすると、B君の癖をよんでいるのかもしれない。
今回の解説では、片側検定、両側検定の話しは意図的に省略してある
先ほど1997年の男子大学生110の身長データをhというオブジェクトに入れた。【ここをクリック】のもう少し下の方には、1998年当時60歳以上の男性51人の身長のデータものっている。大学生(平均20.52歳）と60以上（平均68.61歳）では、48歳ほどの差がある。同じ成人男性といったって、40年も違う時代に成長期を迎えたわけだから、身長には違いがあってもおかしく無いような気がする。この２つのデータをみて、２つのサンプルの母集団（ランダムサンプリングが行われたとして、1997年当時の男子大学生全てと1998年当時の60歳以上の男性全て）の間で身長には差があると言っても良いだろうか？
60歳以上の男性の身長データをオブジェクトに格納。scan()関数を使って
> h2=scan() 1: 1547 1689 .. ..＜省略＞.. 52:
hとh2のそれぞれの平均値やヒストグラムを比べる。
> summary(h) > summary(h2) > hist(h) > hist(h2)
平均値やヒストグラムを見た限りでは、けっこう違いがありそう。
もう一つグラフを作ってみよう。２集団のデータを簡単に比較するときに便利な関数が、boxplot()。
> boxplot(h,h2)
さて、いよいいよ検定。考え方は上で示したのと全く同じ。ただし、データの持っている性質が違うので、p値の計算方法は異なる。このような２集団の違いの検定には、Wilcoxonの順位和検定という方法を使う。検定の説明はブラックボックスになりがちだが、上で紹介した粕谷さんの本や、中澤さんのRの使い方pdfには、詳しい考え方の説明があるので、読んでおいて欲しい。今回はRを使った統計的検定がいかに楽かということを示すのが目的なので、詳しくは触れない。
まず、ジャンケンの例で見たように、帰無仮説と対立仮説を決める。
帰無仮説： ２つのサンプルの母集団の身長に差は無い 対立仮説： ２つのサンプルの母集団の身長に差はある
２つの母集団に違いがなければ、hとh2という集団のデータは、母集団から抽出したときに、たまたま偶然で偏ってしまったものということ。そういう偶然の偏りでこれだけの違いが生じるのは、どのくらいの確率(p値）で生じるのかを計算する。Rではwilcox.test()という関数を使います。
> wilcox.test(h,h2) Wilcoxon rank sum test with continuity correction data: h and h2 W = 4957, p-value = 5.329e-15 alternative hypothesis: true mu is not equal to 0
p値の値は 5.329e-15 で、これは 5.329x10^(-15) を示している。つまり、0.0000000000000005329。これは有意水準の0.01%よりもずーっと小さい値なので、２つの身長データの違いが同一母集団から偶然のサンプリングのバラツキでで生まれたということはほとんど無い、ということを示している。統計学っぽく表現すると、「２集団の差をWilcoxonテストで検定した結果、0.01%の有意水準で有意な差があった」と書く。
どうやら、40年も年齢の異なる２つの母集団の身長には差があると言って良さそう。やっぱり、成長期に摂取できた栄養の影響などが、きっと大きいのだろう。
moodleページから提出すること。
http://bean.bio.chiba-u.jp/moodle24/