Rを使ったプログラミング †

↑

【トピックス】帰無仮説・有意水準・有意差・仮説検定。。。が分からない †

前回授業でやった仮説検定のところがよくわからなかったという意見が多かったです。そこで、正確さは二の次にして、仮説検定で行っている作業をイメージで理解してもらえるようまとめてみました。それぞれどこの部分が、帰無仮説、有意水準、帰無仮説の棄却にあたるのか考えてみてください。また、なんとなく理解できたら、統計の教科書を読んでみることを強くお勧めします。

snmpenwl.gif
WS000002.JPG

1枚のコインを20回投げたら、オモテが5回しかでなかった。このコインはいかさまコインかなー？


もし、このコインがいかさまコインじゃない（正常なコイン）とすると...

20回中5回しかオモテが出ないのっていうのも、わりとフツーに起こることなんだろーね。


snmpensl.gif
pengraph.gif

わりとフツーに起こるこってどうやって数字で示す？確率で 95パーセントぐらいかな?...
正常なコインを20回投げた時のオモテの数ごとに、確率をグラフで表してみると...




snmpengi.gif

なんと！!95パーセント以上っていうのは、オモテの数が6枚から14枚までの場合だっ！

オモテの数が5枚の場合は、5%の範囲というフツーじゃ起こらないとした範囲に含まれちゃってる...


「20回中5回しかオモテが出ない」っていうのが、「正常なコイン」ではフツーに起こらないっていうわけだから、

このコインが「正常なコイン」っていう前提が間違っていたと考える方がいいんだろうね。

うん、きっとこのコインは、いかさまコインなんだ


snmpenas.gif

でも、待てよ。フツーじゃ起こらないことの判断基準が5%って大きすぎない？
5%っていうと、、、100回中5回。。。 20回中に１回だ！
っていうことは、もしも「正常なコイン」という前提が正しくても、20回に1回は、間違った判断を下してしまう可能性があるってことだね。

↑

仮説検定：用語の整理： †

↑

前回やった２項検定の補足：どうしてp値が0.04139になるか？ †

> binom.test(15,20)
	Exact binomial test
data:  15 and 20 
number of successes = 15, 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.75

Untitled-8.gifUntitled-7.gif
file期待値.xls

↑

Rのおさらい †

　前回やったRを使った計算のおさらいを、サンプルデータをつかってやってみましょう。

↑

サンプルデータ：　仲岡先生の実習から †

下のデータは仲岡先生の実習でやった、調査コドラートの面積と、出現種数です

コドラート1    10*10    2種
コドラート2    15*15    4種
コドラート3    20*20    7種
コドラート4    25*25    7種
コドラート5    30*30    7種
コドラート6    35*35    9種
コドラート7    40*40    8種
コドラート8    45*45   12種
コドラート9   50*50   10種
コドラート10   55*55   10種
コドラート11   60*60   11種
コドラート12   65*65   11種
コドラート13   70*70   9種
コドラート14   80*80   10種
コドラート15  100*100  12種
（データ提供、してくれた方に感謝！）

　Rを使ってこのデータの、面積と種数の相関関係を知りたいとき、必要なのは、コドラートの１辺の長さと種数のベクトルです。

　今から何をするか、もうお分かりですよね？

予想通り、K2Editorを立ち上げて、要らないところを消し、

10   2
15   4
.    .

というデータに整形して、エクセルにコピペして、１列ごとのデータを scan() 関数を使ってオブジェクトに読み込みます。では、やってみましょう。（そろそろミミタコだろうから、説明は簡単に書いておきます）

↑

1. K2Editorでの操作 †

  1. 上の囲みのデータをコピペ
  2. 正規表現置換 
    検索文字列：　^.*\*　（説明：行頭から最初のアスタリスクまで）　　
置換文字列：＜無し＞（つまり削除するということ）

検索文字列：　種　　　　　　　　　　　
置換文字列：＜無し＞

検索文字列：  +　（説明：左に書かれているのは半角のスペースと半角のプラス）　
置換文字列： \t　（説明：タブ）
  3. 終わったら全体を選択してコピー。　エクセルを起動
↑

2. エクセルでの操作 †

  1. K2Editorで整形したデータをペースト
  2. ２列に分かれてデータが表示される 
    10  2
15  4
.   .
  3. 1列のデータをコピーして、Rに移動
↑

3. Rでの操作 †

  1. scan()関数を使って、コドラートの一辺の長さを qに、種数をyに代入 
    > q=scan()
[1]　　　　←　ここでペースト。yについても同様
> y=scan()
  2. （確認のため、それぞれ表示しておきます） 
    > q
[1]  10  15  20  25  30  35  40  45  50  55  60  65  70  80 100
> y
[1]  2  4  7  7  7  9  8 12 10 10 11 11  9 10 12
  3. 欲しいのは面積なので、xにqの2乗を入れる 
    x=q^2
  4. それぞれの常用対数を使ってグラフにプロット（対数グラフ用紙必要無し！ [smile]) 
    > plot(log10(x), log10(y))
  5. 回帰分析 
    > result=lm(log10(y) ~ log10(x))
> result
Coefficients:
(Intercept)     log10(x)  
   -0.1375       0.3250  
ここに表示された結果は、log(y)=0.325*log10(x)-0.1375 という回帰直線の式
  6. 回帰直線の描画 
    > abline(result)
  7. 解析のサマリー表示 
    > summary(result)
Call:
lm(formula = log10(y) ~ log10(x))
Residuals:
     Min        1Q    Median        3Q       Max 
-0.211429 -0.054045  0.006334  0.053616  0.142160 
Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.13750    0.15254  -0.901    0.384    
log10(x)     0.32498    0.04713   6.896 1.09e-05 ***　　←有意水準0.1%で有意
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 　←これ、有意水準
Residual standard error: 0.09927 on 13 degrees of freedom
Multiple R-Squared: 0.7853,	Adjusted R-squared: 0.7688 
F-statistic: 47.55 on 1 and 13 DF,  p-value: 1.092e-05
  8. 相関の検定 
    > cor.test(log10(y),log10(x))
Pearson's product-moment correlation
data:  log10(x) and log10(y) 
t = 6.8957, df = 13, p-value = 1.092e-05
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.6847590 0.9618161 
sample estimates:
     cor 
0.886173 　　←　相関係数

　この授業では、それぞれの分析の意味については説明しません。仲岡先生に聞いて下さい。（グラフ用紙で提出しなくてもいいと仰っていたような、いなかったような。。。）

＊大事なこと＊

「こんな解析、自分でやるのは絶対無理！」と思うかもしれませんが、
 そんなことは無いです！

　上のような解析の方法、私だってRでできるとは知りませんでした。 では今回、私はどうやったのでしょう？「回帰分析行いなさい」なんて言われたら、どうすればいいのでしょうか？もしかしてRでできるかなと思ったら、

Adobe Readerを使って、マニュアルを検索します
↑

仮説検定：補足 †

前回アンケートで、仮説検定について「帰無仮説や有意水準などが分からない。」、「Rの操作は分かったけど、その意味が分からない」という意見があります。でも、私は、今の段階では、

全員がRを使って仮説検定の操作を行うことができただけで十分！

と考えています。生物学では収集したデータをもとに何らかの判断を下す場合、ほとんどの場合仮説検定が必要になります。これは別に系統学や生態学に限ったことではありません。最新号のNatureに載っている染色体不安定性を示すマウスの腫瘍に関する論文でもt検定が使われています。

↑

Rを使ったプログラミング演習 †

　さて、今回で３週目になるRという計算ソフト。コマンドを入力するのはちょと難しいけど、エクセルみたいにマウスを使って面倒な操作をしなくても、一瞬にしてグラフが描けてしまう優れものですね。これなら、良いレポートを仕上げるために、何度も解析をやりなおして、グラフを描き直すっていう作業も面倒ではなさそう。。。 　初めて使ったソフトなのに、皆さん、すでにRの操作には慣れてきているし（全員、課題ができていた）、興味も感じてくれたよう。そこで、今回は要望の多かったプログラミングをRを使って勉強します。

↑

プログラミングって何？ †

　ところで、プログラミングって何をすることでしょうか？IT用語辞典で調べてみると、、、

　　載っていません

　唯一近いところでは、「プログラミング言語」というのがあります。説明文からプログラミングに関わるところを抜き出して、「プログラミング」の説明文を作ってみると、

　今日の授業でコンピュータにあれこれ命令を与えるために使うRは、計算ソフトですが、インタプリタ型のプログラミング言語ということもできます。

　つまり、プログラミングとは、簡単に言うと、「言葉を使ってコンピュータが理解できる命令を作ること」　です。

↑

コンピュータにどんな命令ができるのか？　or　どんな命令がしたいか？ †

　「Wordを立ち上げて文章を作成し、印刷する」とか、「Excelで家計簿を管理する」とかは、市販のソフトウェアを使った作業ですが、いずれも、非常に高度で複雑な命令が、使用者からコンピュータに向かって発せられています。マウスを動かしたらポインタが動くなんていうのは、けっこう複雑な命令ですよ。でも、ほとんどの場合、一般の人がコンピュータでやりたいと思うほとんどの作業は、わざわざ「命令」なんていうものを意識する必要は無いはずです。利用者が理解しやすく・使いやすいように、専用のソフトウェアが開発され、販売されているでしょう。

　生物学の研究者でも、プログラミングなどを経験したことの無い人は大勢います。だって、今や、DNAシーケンスの決定や、整列や、系統樹作成にだって、専用のソフトウェアがあるのですから。。。それぞれのソフトウェアの使い方を覚えれば、それですむわけです。私たちも新しい解析方法が開発されるたびに、そのソフトウェアの使い方を覚えるのに苦労しています。

　でも、それだけで、本当にいいのでしょうか？良くないという点が、まず、２つ思いつきます。

1. 自分の目的にあったソフトウェアがいつも存在するというわけではない
2. ソフトが存在する場合、目的の数だけ、ソフトの使い方を覚えなければならない

　ということです。特に１番目は致命的ですね。ソフトウェアが開発されてなかったら、目的が達成できないのですから。例えば、集団遺伝学で学ぶ、遺伝的浮動の解析ソフトなんていうのは、パソコンショップに行っても売ってないです（でも、じつは、このくらいメジャーなものなら、ネットで探せばフリーのものがあります）。

　そこで、この授業でプログラミングを勉強する理由は次の３点です。

  1. 市販のソフトではカバーできないような目的を、コンピュータに命令を与えることで達成する
    • あなたのコンピュータは実は、すごい性能をもっています。その性能はできる限りつかわないと損
  2. プログラミングを行うことで、目的とする生物学的現象の理解が深まる
    • コンピュータは石頭なので、物事を本当に論理立てて説明してやらないと、命令を聞いてくれません
      　　プログラミングの過程で、生物学的現象の意味を論理立てて考える必要に迫られます
  3. プログラミングを行って、命令を実行させることで、初めて、コンピュータに作業をさせていることが実感できる
    • ワープロや表計算とかは、どちらかというと、道具を使って作業をさせられている感が強いですよね
↑

初めて(?)のプログラミング："Hello World!" †

「Hello World!」だなんて、なんか変なタイトルだなーと思いますよね。誰が使い始めたのかは知りませんが、様々なプログラミング言語の教科書で、最初に出てくるのがこのプログラミングです（プログラムを勉強したことのある人なら、たいてい誰でも知っています）。では、早速やってみましょう。Rを立ち上げて、次のように入力してください。

print("Hello World!")

うまくご挨拶できましたか？

ここでやったのは、print()という関数を使って"Hello World!"と画面に表示させただけですが、これも

画面にHello Worldと表示させる

ということを目的としたプログラミングです。

えっ？それなら、Rの最初でやったオブジェクトの内容を画面に出すのもプログラミング？

と聞かれるかもしれませんが、その通りです。次の囲みの中を全てコピーして、Rのコンソールにペーストしてみましょう。

x=c(1,2,3,4,5,6,7,8,9,10)
sum(x)

上でやったプログラミングは、

画面に1から10までの整数の合計を表示させる

というプログラミングです。皆さんがこれまでにやってきたRを使った操作も、何らかの処理の結果を画面に表示させるプログラミングだったわけです。なんとなく、プログラミングが身近になったでしょ？

↑

複数行の命令文 †

上でやった２つのプログラミングですが、違いがありますね。"Hello World!"の方は１行だけの命令だし、合計値の方は2行になっています。プログラムは通常、いろんな処理を組み合わせて作るので、複数行になることが多いです。次に３行からなるプログラミングをやってみましょう。

q=c(10,15,20,25,30,35,40,45,50,55,60,65,70,80,100)
y=c(2,3,4,5,6,7,8,9,10,11,10,11,9,10,12)
plot(log10(q^2), log10(y))

このプログラムは、

２つの数値ベクトルについて計算を行って、対数グラフを描く

というものでした。さっき、上でやったのと同じことを、複数行にまとめて書いて、まとめてRに命令したわけです。

↑

Rで使う命令 †

　授業では一つ一つ解説しますが、復習課題をやるときとか、自分で新たな解析に挑戦するときは、どういう命令があるかが書かれた説明書が必要になるでしょう。

↑

プログラミングの基礎: 繰り返しと代入と条件分岐 †

上で複数の行で様々な命令を一度にＲに与えられるということを解説しました。でもやったことは、それぞれ個々の計算を順番に並べて、一度にペーストしただけですよね。そこで、次にもっとプログラムらしい命令をRに与えてみましょう。それは、繰り返しと代入と条件分岐というもので、プログラミングの基本中の基本です。

↑

繰り返し †

人が不得意でコンピュータが得意なことの一つは、単純な繰り返し作業を際限なく（文句も言わずに）行うことです。要するに、プログラミングの一つの目的は、そういう面倒な繰り返し計算をコンピュータにやらせることです。では、何回繰り返すかというのをどうやってコンピュータに指示するかというと、forという命令を使います。

↑

代入 †

上の for(i in 1:10) { print(i) } というプログラムでは、iの値は繰り返しの度に違う値になります。

1回目  iの値は1
2回目　iの値は2
3回目　iの値は3
....   ........
10回目 iの値は10

では、一つ前の繰り返しの時のiの値を使いたい場合はどうすれば良いでしょうか？例えば、1から10までの合計をfor命令を使ってプログラミングすることを考えると、

今の繰り返し回数のiの値に１つ前回数までの合計を加えて表示させたい
1回目  iの値は1   １つ前までの合計は0  1 を表示
2回目　iの値は2   １つ前までの合計は1  3 を表示
3回目　iの値は3   １つ前までの合計は3  6 を表示
....   ........
10回目 iの値は10  １つ前までの合計は45 55 を表示

こういう場合に代入という命令を使います。その命令は

=

です。通常の計算では'='というのは「3+5=」というように、計算を行わせる意味で用いられますが、プログラミングにおいては、

左辺の変数に右辺の計算結果を代入する

という意味で用いられます。上の計算のプログラムは次のようになります。

goukei=0           #goukeiという変数に初期値0を代入
for(i in 1:10) {   #以下を10回繰り返す
  goukei=goukei+i  #goukeiに前回までのgoukeiの値にiの値を足したもの代入
  print(goukei)    #goukeiの内容を画面に表示
}

goukei=goukei+iという命令文では、右辺の「goukei+i」を先に計算して左辺の「goukei」に代入しているところが、分かりにくいかもしれません。

練習問題：　1から10までの数字を１つずつ順々にかけ合わせたた結果を全て表示させるプログラムを作りなさい
            1x1=1, 1x2=2, 2x3=6, 6x4=24, 24x5=......　と計算してゆくということ
↑

条件分岐 †

プログラミングの基本技の最後は１つは条件分岐です。if()を使って表現します

i=3                     #iに3という数値を代入
if(i==3){print("三")}  #iの値が3ならば画面に三と表示

この例の２つめのif命令で、iの値を評価しています。評価に使われるのは比較演算子というもので、

==　    等しければ
!=      等しく無ければ
>, >=   左辺が右辺よりも大きければ、左辺が右辺以上ならば
<, <=   左辺が右辺よりも小さければ、左辺が右辺以下ならば

ということを意味しています

↑

繰り返しと条件分岐 †

　では、上の繰り返しと条件分岐を組み合わせて見ましょう

練習問題：for命令を使って1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
↑

第10回授業の課題 †

↑

課題１．アンケート調査 †

  1. http://bean.bio.chiba-u.jp/joho19/ に、「自分のID」/10　という新しいページを作成し、下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。
↑

課題１．意見調査 †

　下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。

*第10回授業アンケート
**氏名：
**課題への回答
-今日（6月22日）の授業の進み方は？（はやい、丁度いい、おそい）
--回答：
-今日の授業の難しさはどう感じましたか（簡単すぎ　簡単　丁度いい　難しい　難しすぎ）：
--回答：
-難しいと答えた人は、特にどの点が難しかったですか？：
--回答：
-今日の授業は（よく分かった　分かった　分からなかった）：
--回答：
-分からないと答えた人は、特にどの点が分からなかったですか？：
--回答：
-今日の講義で理解できなかった用語があったら挙げてください：
--回答：
-Ｒで何ができるのか、だんだん分かってきたと思います。今後自分の勉強にＲを使いたいですか？
--回答：
-Ｒでも、Ｒ以外の言語でも、何かやってみたいプログラミングはありますか？
--回答：
-今回やったプログラミングはおもしろかったですか？
--回答：
↑

復習課題：Rによるプログラミング：繰り返しと条件分岐 †

次の問題文の指定にあるプログラムを作成しなさい。作成したプログラムは課題提出ページに貼り付けなさい。 行頭には必ず半角の空白を入れること。

  1. 1から20までの総和を計算するプログラム（?のこと）（注：繰り返し文(for())を使ってプログラムを書いてください。sum()関数は使わないこと）
  2. 　乱数を20個発生させて、0.5未満のものだけを画面に表示させるプログラム
↑

予習課題：Rを使った関数の定義 †

「関数の定義」なんていうと難しそうですが、ようするに、作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにすることです。例えば、

1から入力した数値までの全てを横一列に表示させるプログラムを作り
displayという名前の関数として定義する

どうやるかというと、

display=function(a){       #関数定義の始まり
 kekka=c()                 #kekkaに空ベクトルを代入して初期化
 for(i in 1:a){              #a回（iの値を1からaまで変化させる間）繰り返し
   kekka=append(kekka,i)     #kekkaというベクトルにiを要素として代入
 }
 print(kekka)              #kekkaの内容を表示
}                          #関数定義の終わり

上の方で作ったプログラムの上と下を関数をfunction(){}で囲んで、オブジェクトに代入するだけですから、簡単です。この関数を実行するには、最初aの値を()の中に入力して、

display(10)

とします。 ↑を使って何回も実行してみると、結果がいろいろ変わるのがわかりますよね。

そこで問題です。上の関数の定義方法に従って、入力した数までの合計値を計算するsumupという名前の関数を作成します。下の囲みの中の_の部分（１文字に対応するとは限らない）を埋めて、プログラムを完成して提出しなさい。

 sumup = ______       #関数sumpuを定義
   kotae = _          #kotaeを初期化
   _ (i in ___){      #1からaまで繰り返し
     kotae = _____    #kotaeにiの値を足したものとkotaeに代入
   }                  #繰り返し終了
  _____               #kotaeを表示
  }                   #関数定義の終了