Rを使った統計解析

アナウンス: 「千葉大学留学フェア」開催中(6月9日(月)〜13日(金))

第9回授業の獲得目標: [worried]

Rに慣れ親しむ [smile]

Rは、これまで皆さんが使ってきたMS WordやMS Excelとは全く異なる使い方をするソフトウェアで、ちょっと、心理的な敷居が高いかも知れない。
でも、ぜひ、使って、慣れてみてほしい。今回もゆっくり復習しながらRを使ってみるので、まずは、「こんなことができるんだ」ということを分かって欲しい。
それと、解析に使うコマンドや関数は、必ずしも覚えなくても構わない(私も覚えていない)。大切なのは、「目的の操作」がどこに載っているかを知ること。別に、Googleで解析方法を検索しても良いし、この授業ページに戻って方法を確かめても構わない。まずは、「自分で操作して、結果が表示される」というところを目指そう。

そして、Rに慣れてきたら、一人で使い方を調べられるようになろう。  参考になる情報

Rの使い方をさらに勉強したいのなら、本屋や図書館にある様々な参考書から、どれか自分の気に入ったものを一冊を選んで、随時参照しながら操作をするといいだろう。
下に示す、フリーのオンライン情報だけを使って勉強しても、十分使いこなせる。

Rを使った統計解析 [smile]

生物学を学ぶ限り、自分が実験や観察で得たデータを使うときに、統計解析は必ず使うことになる。生物学のための統計解析の本も、数多く出版されている。一年生の皆さんも、機会をみつけて、できるだけ統計学に慣れ親しんでおいて欲しい。そんなとき、Rは非常に役に立つ。Rは非常に強力な統計解析ソフトであるため、与えられたデータの持つ統計的特徴を瞬時に計算したり、データの持つ意味を推測するための図版を一瞬に作図してくれる。この授業ではRを使って実際にデータを解析することで、統計解析に慣れ親しむ方法を解説する。実際の理論的裏付けは、統計の本を読んだり、授業に出たりして、各自学習してほしい。

統計解析の2つの主要目的

 この授業では統計解析を、次の2つの目的のために使う。

1. ある対象(母集団)から得られた部分的な数値データから、その対象の持っている性質や特徴を知ること

2. ある対象や実験で得られた数値データを用いて、その対象についてどういう判断を下せば良いかを論じること

母集団の性質の推測

"ある対象(集団)から得られた部分的な数値データから、その対象の持っている性質や特徴を知ること"

っていうけれど、「母集団」って何?
「母集団」というのは、興味の対象になっている集団全体のことだ。
例えば、
生物学の大きな目的の一つに、自然界に存在する様々な生物の特徴を知ることがある。しかし、自然界の全ての生物個体を計測して、特徴を知ることは不可能。そこで、対象とする生物から、一部だけ(標本)を偏り無く取り出して(これを任意抽出とかランダムサンプリングと呼ぶ)、特徴を計測し、対象生物全体(母集団)の特徴を推定する。こういう作業を「統計的推測」と言う。

Untitled-1.gif

この、「一部を使って全体を知る」という点が、統計の基本

 生態学などのマクロ系の生物学では、上のような方法で標本を採集し、母集団の性質の解析に統計学的手法を用いることが、非常に多くある。この授業では、統計学全般について語っている時間は無いので、ほんのさわりしか扱わないが、下に示す参考書などを読んで、勉強しておいて欲しい。受験で培った数学力が落ちておらず、学問に対する情熱が非常に高い、1年生の今が統計学習得のチャンスですぞ!

それでは、実際に統計的推測やってみよう。。

【ここをクリック】←ここに、産業技術総合研究所 デジタルヒューマン研究センターから提供して頂いた、1997年の男子大学生110人の身長、体重のデータがある(表示されない人は、moodleページを見てみよう)。このデータを使って、1997年当時の大学生の身長と体重の一般的傾向について議論してみよう。

このとき、

母集団1997年当時の大学生全て

標本1997年の男子大学生110人

標本母集団からランダムに(無作為に)抽出されたとする

今回の統計的推測では、

110人分のデータから、1997年当時の大学生全体の身長と体重の傾向を推測する

ことを目的としている。

  1. まず、Rにデータを取り込む。データは数が多いし空白で区切られているので、c()を使ってデータを取り込むよりは、scan()という関数を使って、h というオブジェクトにデータを入れる方がいい(データの取り込みには他にもいろんな方法があるが、ここではコピー・ペーストで行えるscan()を使って解説する)。
データを一まとめにしてオブジェクトに代入する2番目の方法
コピー・ペーストや手入力でデータをまとめたいときには、scan() という関数を使う
> h=scan()
1:

画面に表示された 1: の後に、上のページから身長データをコピー・ペーストして、最後にenter キーを押す。次のように表示される。

1: 1775  1710 ..
  ..<省略>..
111: 

同様にして、体重データも w というオブジェクトにを入れる。次のように表示される。

1: 79.8 58.0 ..
  ..<省略>..
111: 

確認のために

> h
> w

と入力。データが一覧表示される。

  1. 標本データの視覚化

    計測データが得られたら、まずは、グラフにしてみるよう。グラフとして視覚化することで、データの持つ性質が直観的に理解できる。こういうとき、Rは、いろんなグラフを一瞬にして描画してくれるので、非常に便利。
    まずはヒストグラム(柱状グラフ)を使っみる。Rを使ってヒストグラムを書く場合は

    > hist(h)
    キーボードからたった7文字打ち込むだけで、身長の度数分布がグラフで表示される。
    hist.gif
    このグラフを見ただけで、標本データでは身長170cmぐらいの学生が最も多く、
    データは釣り鐘型の分布をしているということが分かる。
    同様のことを体重についてもやってみよう。
    > hist(w)
  2. 要約統計量の表示

    次に、この標本データの平均とか最大値とか、最小値とか、データの集まりがもっている基本的な性質を表示させて見よう。次のように入力。

    > 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.というのは、全体の4分の1の人数(つまり25%)がどの高さまでの間に含まれているかということを表している。つまり、167.1cm以下の人が25%、173.4cm以上の人が25%、その間の人が全体の50%ということ。
    体重も同様に解析してみよう。
    > summary(w)
  3. 身長と体重の関係を視覚化

    身長も体重も釣り鐘型の分布になっていた。では、この2つの数値の間に何か、関係があるだろうか?これもまた、グラフで表示させてみる。こういう場合、縦軸に体重、横軸に身長をとって、一人一人の持つ身長と体重の散布図を作る。Rの場合、次のように入力するだけ。

    > plot(h,w)
    身長と体重の関係は右上がりの直線関係にありそうだ。
    2変数間の関係は、相関係数で表され、Rで簡単に計算できる(→R-Tipsの説明もみておこう
    >cor(h,w)
    [1] 0.5872979   #値が正なら正の相関
    plot.gif
  4. 統計的推定で分かったこと

    今回、110人分のデータを使って、1997年当時の大学生男子の身長と体重について、

    グラフによる視覚化
    要約統計量の表示
    散布図の作成
    という3つのことを行った。この標本が母集団の性質を表しているとするならば、1997年当時の大学生について、次のような推定ができる。
    1997年当時の男子大学生は、身長の平均が170.6cmで、データは釣り鐘型に分布している
    体重は平均が59.53kgで、身長と同様にデータは釣り鐘型に分布している
    身長と体重の間には、正の相関がある
    これで、1997年当時の男子大学生という母集団の身長・体重について、統計的推定を行うことができた!

    なお、今回は分散や標準偏差の話しは省略した。Rを使った統計解析の基本(データの図示、代表値の計算、データ間の相関)などについてはここでは説明しきれないので、下のリンクにある中澤さんによる解説などを読んで、勉強してみて欲しい。。

演習(時間があったら):エクセルファイルのデータをグラフ化

エクセルデータをそのまま読み込むには、http://cse.naro.affrc.go.jp/takezawa/r-tips/r/41.html を参照。

統計的検定(仮説検定)を直観的に理解する [smile]

この授業で対象とする統計解析の2つめは、
"ある対象や実験で得られた数値データを用いて、その対象についてどういう判断を下せば良いかを論じること"

それでは、次の設問について、皆さんの判断を考えて下さい。

判断基準が重要

皆さんの判断は様々だったが、こういう問題はどのように考えるのが良いのだろうか?問題は、

判断基準

だ。科学的に説明するためには、論理的な「判断基準」が無ければならない。
そして、ある「判断基準」に基づいて判断を下すのが、これから学ぶ「仮説検定」だ。
まず、この授業では、正確さは二の次にして、仮説検定の流れを直観的に理解してもらえるような説明を試みる。
(この説明のどこが、後で出てくる「帰無仮説」、「有意水準」、「帰無仮説の棄却」にあたるのかは、後で考えてみて欲しい。なんとなく理解できたら、どれでもいいから統計の教科書を読んでみて、正確に理解することをオススメする。)

直観的説明

snmpenwl.gif
010605icon075-trans.png

A君とB君が20回ジャンケンをしたら、B君は5回しか勝てなかった。ジャンケンの勝ち負けって、半々の確率で決まるはずだけど、こんなに勝ち負けが偏ることって、たまたま偶然で起こったことなのかなー?それとも、A君が何かズルをしていたと考えた方がいいのかなー?。。


ジャンケンの勝敗が、ズルとかじゃなくて、本当に半々の確率で決まっているのなら、...

「20回中5回しか勝てない」っていうのは、偶々(たまたま)起こったことなんだろうね。でも、その「偶々(たまたま)」というのは、どの程度偶々(たまたま)起こることなんだろうか?


snmpensl.gif

よし。「ジャンケンの勝敗が完全に半々の確率で決まるというのを前提として、20回ジャンケンをする実験」を何度も何度も繰り返したとき、片方が何回勝つかをグラフで表してみよう... もしも、「20回中5回しか勝てない」っていうできごとが、統計的には滅多に起こらないことだとしたら、そもそも「A君とB君のジャンケンの勝敗が半々の確率で決まっている」という前提が間違っていると言って良いんじゃないかな?
そうすると、「滅多に起こらない」をどうやって判断するか、基準が必要だけど、ここでは、5%を基準にして、実験を繰り返したときに全体の5%でしか起きないことは、滅多に起こらないと判断することにしよう。




#09_2.gif
snmpengi.gif

グラフを作ってみると、なんと!! 「20回中5回しか勝てない」状況は、全体の5%未満しか起きていない...。つまり、5%を基準にすると、 「20回中5回しか勝てない」ことは「滅多に起こらない」んだ!


「ジャンケンの勝敗は完全に半々の確率で決まる」という前提の下では、「B君が20回中5回しか勝てない」っていうのは滅多に起こらないというわけだから、

「ジャンケンは完全に半々の確率で決まる」っていう前提が間違っていたと考える方がいいんだろうね。

うん、きっと、どんな方法かは分からないけれど、A君はB君にジャンケンで勝てる方法を見つけてるんだ


snmpenas.gif

でも、待てよ。「滅多に起こらないこと」の判断基準が5%って大きすぎない?
5%っていうと、、、100回中5回。。。 20回中に1回だ!
っていうことは、もしも、「ジャンケンは完全に半々の確率で決まる」にしても、20回に1回は、間違った判断を下してしまう可能性があるってことだね。

もう少しだけ詳しい統計的検定の説明

上の例では、20回のジャンケンの結果を使って統計的検定の直観的な説明を試みた。こんな単純な例であっても、実際の生物学の研究で使う、統計的検定と、推論の進め方は同じだ。(注意: 下の例は、今回のジャンケンの例にあわせるために、かなり説明を端折っている。ここで直観的に理解したら、上記統計の教科書を見てみよう)

まず、統計的検定で、知りたいのは、

「A君とB君が20回ジャンケンをしたとき、5回しか勝てない」という状況は、
「ジャンケンの勝ち負けが、本当に半々の確率で決まる」という前提のもとで、
 単なる偶然で起こるようなことか、
 あるいは、偶然では滅多に起こらないことか」

ということ。
このように、得られたデータが単なる偶然で生じることかどうかを知るときに、統計的検定という作業が必要になる。

 統計的検定では、次のような論証う。

  1. 帰無仮説と対立仮説を作る
    問題にしているのは、 「A君とB君が20回ジャンケンをしたとき、B君が5回しか勝てない」という状況は、「A君とB君のジャンケンの勝率には違いがあるかどうか」ということ
    このとき、「A君とB君のジャンケンの勝敗」の決まり方について、次のような仮説を作る
    A君とB君のジャンケンの勝敗は、半々の確率でランダムに決まる
    この仮説のことを 帰無仮説と言いう。(統計の教科書を見ると、帰無仮説には「違いがない」とか「差がない」という仮説を使うと書かれている。今回の場合「『ジャンケンの勝敗は、半々の確率で決まる』というだけで、『差が無い』とか『違いが無い』とかは言って無いじゃないか?」と思う人がいるかもしれない。しかし、「ジャンケンの勝敗は、半々の確率でランダムに決まる」ということは、「ジャンケンの勝ち・負けそれぞれが生じる確率には違いが無い」ということを意味している。字面だけで判断するとこんがらがる場合があるので、要注意。)
     帰無仮説に対応して、あなたが持っている仮説は、得られたデータは偶然のせいでたまたま偏ったのでは無く、対象の持っている性質そのものが偏っているというもの。この仮説の事を対立仮説と呼ぶ。今回の場合、対立仮説は、
    A君とB君のジャンケンの勝敗は、半々の確率ではきまらない(A君はB君よりもジャンケンに強い)
    というもの。対立仮説は帰無仮説を否定したものになっている。
  2. p値の計算
     次に、帰無仮説が正しいと仮定した場合に、得られたデータが実現する確率を計算する。
    「A君とB君のジャンケンの勝敗は、半々の確率で決まる」(帰無仮説)と仮定しているけれど、サンプリングの際の偶然のバラツキで、20回中5回しか勝てないということだって、よくあることなのかもしれない。
    そこで、「20回ジャンケンをしたとき、B君が5回しか勝てない」という確率を計算して、その確率を用いて、「滅多に起こらないこと」かどうかを判断する。この確率の事をp値(Rではp-valueと表示される)と呼ぶ。

     このp値の計算で、一つ注意しなければならないことがある。それは、
    p値の計算をするときに、「20回中勝ちが5回」の場合だけを計算するだけではなく
    20回中勝ちが5回以下(5, 4, 3, 2, 1, 0)の場合の確率も計算しなくてはならないということ
    ということ。この後の議論に出てくるように、起きた事象がどれだけ希にしか起きないことかを判断して、それが滅多におこらないようなことなら、前提となった帰無仮説が間違っていると判断する。つまり、「勝ちが5回以下」のそれぞれのできごとは「勝ちが5回」という出来事よりも起こりにくいことなので、そう言う場合にも、帰無仮説は棄却される。そのため、p値の計算は「5回以下」の場合の確率が必要になるの。この確率のことを累積確率と呼ぶ。
    hist(xlim=c(0,20),rbinom(10000,20,0.5)) でグラフを表示させ、勝ちが5回以下の範囲を見てみよう

    今回の解説では、片側検定、両側検定の話しは意図的に省略してある

  3. 有意水準との比較
     次に、得られたデータが実現する確率(p値)ある基準とを比較することで、こういう状況が起きることが、「滅多に起こらないこと」かどうかを判断する。多くの場合、基準には0.05(つまり5%)あるいは0.001(つまり0.1%)とか0.0001(つまり0.01%))という値が使われる。この基準のことを有意水準とか危険率と呼ぶ。上で得られたp値と有意水準を比較すると、次の2つの場合が考えられる:
    1. p値が有意水準よりも低い: 「A君とB君のジャンケンの勝敗は半々の確率で決まる」という仮説(帰無仮説)のもとで、「たまたまデータに偏りが生じて5回しか勝っていない」という状況が生じる確率(p値)が、「滅多に起こらないこと」の判断基準(有意水準)よりも小さいので、この状況は滅多に起こらないと考える。このとき、今の考えの前提になった帰無仮説は否定されるということになる。
      つまり、
      「A君とB君のジャンケンの勝敗は半々の確率で決まる」(帰無仮説)
                (この場合、5回しか勝てないというデータの偏りは単なる偶然)
      が棄却されたので、
      「A君とB君のジャンケンの勝敗は、半々の確率ではきまらない」
                  (対立仮説)(この場合、データの偏りは、単なる偶然では無い)

      という仮説が採択される(こういう論理の進め方を背理法と呼ぶ)。
      今回の場合、「得られたデータ(この場合、ジャンケン20回におけるA君B君の勝率)には、...%の有意水準で、有意差がある」と言う。また、「「20回中5回しか勝てないという状況は、偶然で生じることはほとんど無いと考えられるので、A君とB君のジャンケンの勝敗は半々では決まらない」と考える。
    2. p値が有意水準よりも高い:上と同様に考えると、
      「A君とB君のジャンケンの勝敗は、半々の確率で決まる」(帰無仮説)
      が正しいという仮定の下で、「20回中5回しか勝てない」という状況が生じる確率(p値)が、「滅多に起こらないこと」の判断基準(有意水準)よりも大きい。言い換えると、「A君とB君のジャンケンの勝敗は、半々の確率で決まるときに、20回中5回しか勝てないという状況は、別に不思議ではない」ということ。つまり
      「A君とB君のジャンケンの勝敗は、半々の確率で決まる」(帰無仮説)
      という仮説は棄却 されない
      得られたデータからは帰無仮説を否定できなかったので、この場合、得られたデータ(この場合、ジャンケン20回におけるA君B君の勝率い)には...%の有意水準で、有意差が無いと言う。ここで注意しなければならないのは、「有意差が無い」からといって、「A君とB君のジャンケンの勝敗は、半々の確率で決まる」と結論づけてはいけないこと。「有意産が無い」のは、「A君とB君のジャンケンの勝敗は、半々の確率で決まっていてもおかしくはない」ということを意味しているだけで、単なる偶然の結果かもしれないし、データ数が少ないせいかもしれないし、いろんな場合があり得る。

Rを用いた実際の解析:

 では、実際に先ほどの例を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値の値は0.04139なので、有意水準が5%よりも小さい。そこで、帰無仮説「A君とB君のジャンケンの勝敗は、半々の確率で決まる」は棄却される。(もちろん、有意水準を0.001にすると棄却されない)

つまり、「A君とB君のジャンケンの勝敗は、有意水準5%で、半々の確率で決まるとは言えない」ということができる(有意水準0.1%の場合は、「2人の勝率に差があるかどうかはわからない」)。「有意水準5%では、A君とB君のジャンケンの勝敗は、半々の確率で決まるという帰無仮説は棄却された」という言い方もある。

どうやら、やっぱりA君の方がB君よりも強いようだ。もしかすると、B君の癖をよんでいるのかもしれない。

今回の解説では、片側検定、両側検定の話しは意図的に省略してある

Rを使った実際の解析2:2集団の差の検定

 先ほど1997年の男子大学生110の身長データをhというオブジェクトに入れた。【ここをクリック】のもう少し下の方には、1998年当時60歳以上の男性51人の身長のデータものっている。大学生(平均20.52歳)と60以上(平均68.61歳)では、48歳ほどの差がある。同じ成人男性といったって、40年も違う時代に成長期を迎えたわけだから、身長には違いがあってもおかしく無いような気がする。この2つのデータをみて、2つのサンプルの母集団(ランダムサンプリングが行われたとして、1997年当時の男子大学生全てと1998年当時の60歳以上の男性全て)の間で身長には差があると言っても良いだろうか?

60歳以上の男性の身長データをオブジェクトに格納。scan()関数を使って

> h2=scan()
1:  1547  1689 ..
  ..<省略>..
52: 

hとh2のそれぞれの平均値やヒストグラムを比べる。

> summary(h)
> summary(h2)
> hist(h)
> hist(h2)

平均値やヒストグラムを見た限りでは、けっこう違いがありそう。

Untitled-6.gif
もう一つグラフを作ってみよう。2集団のデータを簡単に比較するときに便利な関数が、boxplot()。

> boxplot(h,h2)

さて、いよいいよ検定。考え方は上で示したのと全く同じ。ただし、データの持っている性質が違うので、p値の計算方法は異なる。このような2集団の違いの検定には、Wilcoxonの順位和検定という方法を使う。検定の説明はブラックボックスになりがちだが、上で紹介した粕谷さんの本や、中澤さんのRの使い方pdfには、詳しい考え方の説明があるので、読んでおいて欲しい。今回はRを使った統計的検定がいかに楽かということを示すのが目的なので、詳しくは触れない。

 まず、ジャンケンの例で見たように、帰無仮説と対立仮説を決める。

帰無仮説: 2つのサンプルの母集団の身長に差は無い
対立仮説: 2つのサンプルの母集団の身長に差はある

 2つの母集団に違いがなければ、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%よりもずーっと小さい値なので、2つの身長データの違いが同一母集団から偶然のサンプリングのバラツキでで生まれたということはほとんど無い、ということを示している。統計学っぽく表現すると、「2集団の差をWilcoxonテストで検定した結果、0.01%の有意水準で有意な差があった」と書く。

どうやら、40年も年齢の異なる2つの母集団の身長には差があると言って良さそう。やっぱり、成長期に摂取できた栄養の影響などが、きっと大きいのだろう。

リンク:

第9回授業の課題

moodleページから提出すること。
http://bean.bio.chiba-u.jp/moodle24/