Rを使ったデータ解析 †このページを見たら、まずやっておいて欲しいこと †男子学生のみ!やってください †【ここをクリック】←左のリンクのページに入り、画面の説明に従って、データを入力。今回の復習課題に使う(間違った形式で入力すると、自分でやるデータ整形処理が複雑になるので注意)。入力が同じ時間に集中すると「更新の衝突」がおきるので。「コメントの挿入」ボタンをクリックする直前に、一度、ページをリロード(再読込)すること。 第9回授業の獲得目標: †
Rの日本語説明 †
Rの基本操作を理解しよう: †Rを一人で使う場合: 参考になる情報 †Rには様々な参考書が出版されているので、どれか一冊を選んで買ってきて、随時参照しながら操作をするといいだろう。フリーのオンライン情報だけを使って使いたい人は、次のようにするといい。 1. Rの命令一覧が載っている説明書のPDF(上記)を画面に表示させる 2. Adobe Readerの検索機能を使って、行いたい処理を検索 3. 説明を読んで入力すべき命令文をコピー 4. Rにペーストして、値を編集して実行 こうやって、常に説明書を参照しながら操作すれば、命令を覚える必要は無い。 Rを起動する †まず、インストールしたRのフォルダで、Rのアイコンをダブルクリックする。 > が表示される。この">"をプロンプトと言う。プロンプトの右側には通常、カーソルがあり、「このマークの後に何かコマンド(命令)を打ち込んでくださいよ」とあなたに促している。式を入力して、最後にenter(or returnキー)を押せば計算される。 Rによる簡単な数値計算: Rは関数電卓代わりに使える! †以下、Rを使った簡単な計算。プロンプト(>)から右の部分をコピーペーストすれば、計算できる。 > 4/3*pi*5^3 #半径5の球の体積 > pi #円周率が表示される(デフォルトでは小数点以下6桁まで) > (1+ 2 + 3 + 4 + 5 )/5 #式の間に半角スペースがいくつか入っても大丈夫 > 【ここでキーボードの上矢印(↑)を押す: 前に入力した命令を表示させることができる】 > (10 + 2 + 3 + 4 + 15 )/5 #前に入力した式が表示されるので、左矢印(←)でカーソルを動かし修正
オブジェクトへの数値の代入 †Rでは、好きな名前をつけたオブジェクト("もの")に、数値や文字列などのデータを代入できます。オブジェクトには大文字小文字の区別があります。 > x = 3 #xに3という数字を代入した > X = 4 #Xに4という数字を代入した > x + X #xとXの足し算 (なお、#の後はコメント文と呼ばれ、処理には関係の無い説明を書いておける) なお、上では"="を使って数値を代入しましたが、 "<-"を使っても同じことができます。 > x <- 3 > X <- 4 > x+X (Rの代入はもともと "<- " だけだったが、比較的最近、 " = " も使えるようになった) Rで「沢山のデータを一まとめにして扱う方法」を理解しよう †みなさんは、沢山のデータを1まとめにして扱う方法をすでに知っている。 例えばエクセルを使った計算では、
という表を使ってデータを一まとめにし、sum()関数をつかったり、1つのセルに5を足すという計算結果を全部のセルにコピペして、全ての値に5ずつ0足した結果を得たりした。 (こんなやり方だと面倒くさい) > 10.4+5.6+3.1+6.4+21.7 > 10.4+5 > 5.6+5 Rでは、データを横1行に並べて扱う(これをベクトルという) †そこで、Rでは、たくさんのデータを一まとめにして扱うときには、カンマで区切って横1行に並べて使う。このようなデータの構造をベクトルと呼ぶ。ベクトルとうと、高校の数学でやった「大きさと方向をもった量」を思い出すと思うが、 Rで使うベクトルは、同じ種類のデータに順番をつけて並べたもの と定義しておこう。これは、Rが沢山のデータをまとめて処理するために用いる、データの構造だ。例えば、 (10.4, 5.6, 3.1, 6.4, 21.7) というデータの集まりは、同じタイプのデータに順番をつけてひとまとめにしたものであり、Rでいうところのベクトルだ。
x=c(10.4, 5.6, 3.1, 6.4, 21.7) という命令を使う。上の命令を入力した後、「 x 」とだけ入力すると、xの内容が表示される。 > x [1] 10.4 5.6 3.1 6.4 21.7 上の結果に表示された[1]は何だろうか?これは、表示された行の最初の数字が、ベクトルの1番目の要素だということを意味している。では、次のようなベクトルをyに代入して、内容を表示させてみよう。 > y=c(1984,1985,1986,1987,1988,1989,1990,1991,1992,1993, 1994,1995,1996,1997,1998,1999,2000,2001) ウィンドウの大きさによって表示のされ方は異なるが、もし1行に収まらなければ、 > y [1] 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 [15] 1998 1999 2000 2001 というように、2行目の始まりに[15]というような数字が表示される。これは、2行目が15番目の要素から始まっていることを示している。 > y[5] [1] 1988 と、オブジェクトの名前に順番を[]に入れたものをくっつける事で、目的の要素を取り出すことができる。 オブジェクトへのベクトルの代入 †では、Rで、沢山のデータを一まとめにして、いろいろ計算してみよう。 まずはおさらい。上で出てきた、 c() という関数を使って、オブジェクトにデータを代入でする。 > x=c(1,2,3,4,5,6,7) > entterキーを押しても何もおこらない?いえいえ。エラーメッセージが何も表示されず、新しいプロンプトが表示されたら大成功。 では、プロンプトの後に"x"と入力して、結果を見てみよう。 >x [1] 1 2 3 4 5 6 7 上の様に表示されれば「xというオブジェクトに(1,2,3,4,5,6,7)が代入されたことになる。では、このデータを対象に、統計解析の基本である平均値の計算を計算してみよう。 mean()という関数(meanは英語で平均値の意味)で一まとめにしたデータ(ベクトル)の平均値が計算できる。 命令の入力に自信の無い人は、下の囲みの中をコピー・ペースト。 >mean(x) [1] 4 平均値"4"が表示されれば大成功(これで、あなたがRを使って行う、初めての統計解析は、無事に成功した) > x=c(1,2,3,4,5,6,7) と表示された。そしたら"←"と"Back space"などを使って、数値を好きに変更して、最後に"Enter"キーを押そう。 > x=c(1232,223.33,3 ,4 , 5666) 半角文字であるかぎり、別にスペースがいくつ入っていても大丈夫。入力したら、上と同じことをしてみると、 > x [1] 1232.00 223.33 3.00 4.00 5666.00 > mean(x) [1] 1425.666 これで、自分の好きな値を入力して、平均を計算することができた。 Rを使った統計解析 †生物学を学ぶ限り、自分が実験や観察で得たデータを使うときに、統計解析は必ず使うことになる。生物学のための統計解析の本も、数多く出版されている。一年生の皆さんも、機会をみつけて、できるだけ統計学に慣れ親しんでおいて欲しい。そんなとき、Rは非常に役に立つ。Rは非常に強力な統計解析ソフトであるため、与えられたデータの持つ統計的特徴を瞬時に計算したり、データの持つ意味を推測するための図版を一瞬に作図してくれる。この授業ではRを使って実際にデータを解析することで、統計解析に慣れ親しむ方法を解説する。実際の理論的裏付けは、統計の本を読んだり、授業に出たりして、各自学習してほしい。 統計解析の2つの主要目的 †この授業では統計解析を、次の2つの目的のために使う。 1. ある対象(母集団)から得られた部分的な数値データから、その対象の持っている性質や特徴を知ること 2. ある対象や実験で得られた数値データを用いて、その対象についてどういう判断を下せば良いかを論じること 母集団の性質の推測 †"ある対象(集団)から得られた部分的な数値データから、その対象の持っている性質や特徴を知ること" っていうけれど、「母集団」って何? この、「一部を使って全体を知る」という点が、統計の基本! 生態学などのマクロ系の生物学では、上のような方法で標本を採集し、母集団の性質の解析に統計学的手法を用いることが、非常に多くある。この授業では、統計学全般について語っている時間は無いので、ほんのさわりしか扱わないが、下に示す参考書などを読んで、勉強しておいて欲しい。受験で培った数学力が落ちておらず、学問に対する情熱が非常に高い、1年生の今が統計学習得のチャンスですぞ!
それでは、実際に統計的推測やってみよう。。 【ここをクリック】←ここに、産業技術総合研究所 デジタルヒューマン研究センターから提供して頂いた、1997年の男子大学生110人の身長、体重のデータがある。このデータを使って、1997年当時の大学生の身長と体重の一般的傾向について議論してみよう。 このとき、 母集団は1997年当時の大学生全て 標本は1997年の男子大学生110人 今回の統計的推測では、 110人分のデータから、1997年当時の大学生全体の身長と体重の傾向を推測する ことを目的としている。
> h=scan() 1: 画面に表示された 1: の後に、上のページから身長データをコピー・ペーストして、最後にenter キーを押す。次のように表示される。 1: 1775 1710 .. ..<省略>.. 111: 同様にして、体重データも w というオブジェクトにを入れる。次のように表示される。 1: 79.8 58.0 .. ..<省略>.. 111: 確認のために > h > w と入力。データが一覧表示される。
演習(時間があったら):課題8の結果の分析 †TAの木内君が、課題8のOzawakenによるタイプ練習の結果をRでグラフにして分析してくれた。 X=c(4,17,1,14,3,1,1,1,1) #Excelで集計した結果 Xlab=c("G〜D","C","C++++","B","A","A+","A++","S+","S++") #横軸ラベルの名前 barplot(X,names=Xlab,ylim=c(0,20)) #棒グラフを書く命令 title(main="OZAWA-KENの結果") #グラフタイトルを書く命令 統計的検定(仮説検定)を直観的に理解する †この授業で対象とする統計解析の2つめは、 それでは、次の設問について、皆さんの判断を教えて下さい。
判断基準が重要 †皆さんの判断は様々だったが、こういう問題はどのように考えるのが良いのだろうか?問題は、 判断基準 だ。科学的に説明するためには、論理的な「判断基準」が無ければならない。 直観的説明 †
グラフを作ってみると、なんと!! 「20回中5回しか勝てない」というのは、全体の5%未満しか起こっていない...。5%の未満の範囲っていうのは、「滅多に起こらない」っていうことにした範囲だよなー。
「ジャンケンは完全に半々の確率で決まる」っていう前提が間違っていたと考える方がいいんだろうね。 うん、きっと、どんな方法かは分からないけれど、A君はB君にジャンケンで勝てる方法を見つけてるんだ
でも、待てよ。「滅多に起こらないこと」の判断基準が5%って大きすぎない?
もう少しだけ詳しい統計的検定の説明 †上の例では、20回のジャンケンの結果を使って統計的検定の直観的な説明を試みた。こんな単純な例であっても、実際の生物学の研究で使う、統計的検定と、推論の進め方は同じだ。(注意: 下の例は、今回のジャンケンの例にあわせるために、かなり説明を端折っている。ここで直観的に理解したら、上記統計の教科書を見てみよう) まず、統計的検定で、知りたいのは、 「A君とB君が20回ジャンケンをしたとき、5回しか勝てない」という状況は、 「ジャンケンの勝ち負けが、本当に半々の確率で決まる」という前提のもとで、 単なる偶然で起こるようなことか、 あるいは、偶然では滅多に起こらないことか」 ということ。 統計的検定では、次のような論証う。
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) 平均値やヒストグラムを見た限りでは、けっこう違いがありそう。
> 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回授業の課題 †
課題1.アンケート調査 †
課題2、3の提出方法 †
復習課題2. Rを使った統計解析:ヒストグラムの表示 †
復習課題3.Rを使った統計解析 †
発展課題. Rを使った統計的検定 †
|