BLOG

Rのboxplot()コマンドを使うと、箱ひげ図が描画されることが以前にも取り上げました。こんな感じです。

このように軸ラベルが自動的に数字で入るだけです。ここでは、軸ラベルに任意の文字列を入れる方法を説明します。

①ラベル表示をさせないでboxplotを出力

  boxplot(要素1, 要素2, 要素3, …, axisnames=F) と入力

  これで,ラベルが消えて出力される。

②軸ラベル名の指定

  axis(1, at=c(1, 2, 3, 4…), labels=c(“A”, “B”, “C”, …))

   1: x軸への操作、at=c(1, 2, 3): デフォルトで付けられる軸名(数値)

例えば,上図のboxplotの項目名を左から"KNM_Hit", "KNM_FA", "KJKNM_Hit", "KJKNM_FA"ときは,以下のようになります。

axis(1, at=c(1,2,3,4),labels=c("KNM_Hit", "KNM_FA", "KJKNM_Hit", "KJKNM_FA"))

と入力し[Enter]すると下図のように横軸に任意のラベル名が入ります。

6-2まで、二項分布について説明し、いくつかの場面での使用法を説明しました。

 その際には、二項分布曲線を作成し、0.01や0.05や0.04などのいくつかの起きる確率に線を引いて、議論を行いました。

 しかし、一々曲線を描いて、線を引いてという作業はわずらわしいため、ほとんどの場合、二項検定を行います。

 Rにおいては、試行数、起こった回数、期待値を入力することで、二項検定を行うことができます。

 例えば、コイントス50試行中、30回表が出た場合、このコインは表が出やすいかについて二項検定を行う場合は、この場合の否定される仮説(帰無仮説)は、「コインの表と裏がでる確率は等しい」となります。

 Rでは、このように書きます。  binom.test(30, 50, 0.5)

続いて、例3について説明します。

あなたは、人口10万人のA市の保健担当者です。Bという病気(B症)の発症率は人口の5%とされています(話を簡単にするために、年齢性別など関係なく同じ率とします)。A市でもB症の発症率を調べて、一般的な発症率よりも明らかに多ければ、原因を突き止め、予防策などを講じる必要があります。ただ、全住民を調べるには費用と時間がかかります。

 そこで、A市の人口構成と同じような住民100名に協力を得て調査した結果、100名中B症を発症している住民は8名(8%)でした。一般的な発症率に比べて3%多いので、A市のB症発症率は高いので、何らかの対策が必要と結論づけていいでしょうか。


 この例は、前2つと違う感じがしますが、発症した-発症していないの二項に分けられますので、二項分布が使えます。異なるのは、期待値(発症率)が0.05(5%)であることです。また、試行数は、100名を調査したので100となります。この例の場合、100回調べて、発症者が8名ということです。

 これについて、二項分布を描いてみます(下図)。前2例に比べ、かなり分布が原点側に寄っていることが分かります。実際は、横軸は100までありますが、15辺りからずっとほぼ0ですので、20以上は省略しています。

 緑線は、今回の結果である発症者8名のラインです。これを見ると、確率は0.06程度ですので、100名中で発症率8%はレアなケースのように思えます。ただ、一般的な統計的に有意とすることが多い0.05以下ではないので、その点から見ると多いとは言えないので、対策を講講じる必要はないと結論づけることもできます。

 では、このような場合どう考えたらいいでしょうか。

 現実には、他の事柄も判断材料として考慮すべきと考えます。

 例えば、B症はあることをすることで予防が容易であることが明らかであれば、対策を講じる大きな理由となります。あるいは、予防可能であるが、発症した場合の一人当たりの医療費や介護費用が大きい場合も予防策を講じる大きな動機付けとなります。

 一方、予防に費用がかかったり予防が困難であるが、発症した場合でも治療が容易であったり、生命に大きな影響を与えずに治癒するようなものであれば、広報誌などで注意や発症した場合の早期受診を呼びかける程度にとどめるという決断もあり得ます。

 学術研究においては統計的有意差として0.05以下や0.01以下をもってこの例であれば一般的な発症率に比べて明らかに多いと結論づけるかと思います。つまり、この調査結果からは、有意に多いとは言えないと結論づけられます。

 しかし、この結果をもって行政がどのような対策を講じるのかあるいは積極的には講じないのかは、上記のように統計解析とは別の軸で議論される必要があります。統計解析はあくまでも意思決定のためのツールの一つ(判断材料の一つ)であるべきだからです。

ここでは、6. 二項分布で取り上げた例2について考えてみたいと思います。

例2は、こんな疑問です。

 あなたは学校の先生で、生徒に単元終了時に○×クイズを出します。生徒は、○×で回答しますので、あてずっぽうで回答しても1/2は正解する可能性があります。では、20問問題を出したとして、いったい、何問以上正解できたら、生徒はあてずっぽうでなく、分かって答えているといえるでしょうか?

 あてずっぽうに答えると正答率の期待値は1/2=0.5です。つまり、20問であれば、あてずっぽうで答えても10問は正解できる可能性が高いわけです。では、ここで20問(試行)した場合の二項分布をみてみましょう(下図)。

 このように期待値が0.5となる確率が一番高いのはもちろん、20試行中10回正解する場合です。緑線は確率0.05(5%)のラインで14問正解から15問正解の間です。何十回か20問のテストを繰り返した時、あてずっぽうで回答して15問以上正解できる可能性は5%以下ですので、15問以上正解したら、分かって答えているとしてもいいかもしれません。もっと厳しくするとしたら、赤線の0.01(1%)ラインを基準として、17問以上正解としてもいいかもしれません。

 この二項分布を使って考えられるのは、「各試行がお互いに"独立"」していることが前提です。例えば、試験問題の問1で算出された数値を問2でも使って回答するような場合は、問1で正しい答えが出ると問2で正答できる可能性が高くなるといった関係性がでてくるので、"試行が独立していない"と言えるため、二項分布を使うことはできません。

 二項分布とは、結果が二値(二項)である独立した試行をn回行った時の、結果の確率分布のことです。

 結果が二項というのは、yes-no, (コインの)表と裏, ある病気である-ない, 成功-失敗といったものをイメージしてください。

以下に、いくつか例を挙げます。

例1

 あなたは、コイントスをして裏が出たら10円もらえ、表が出たら10円引かれるというゲームをしています。

 コイントスをして、表が出るか裏が出るかは五分五分(確率1/2)です。しかし、このゲームの主催者は、お金をできるだけ節約するためイカサマをしているかもしれません。すなわち、重心が偏ったコインを使っていて、例えば10回投げた(試行した)時に、表が7回も出るかもしれません。しかし、試行数が10回だったのでたまたま表がちょっとだけ多く出て、200回投げたら確率1/2になったかも知れません。二項分布は、10回試行して7回表が出るのがどの程度レアな(起こりにくい)ことなのか(イカサマをしていないか)について統計学的に教えてくれます。

例2

 あなたは学校の先生で、生徒に単元終了時に○×クイズを出します。生徒は、○×で回答しますので、あてずっぽうで回答しても1/2は正解する可能性があります。では、20問問題を出したとして、いったい、何問以上正解できたら、生徒はあてずっぽうでなく、分かって答えているといえるでしょうか?11問でしょうか?12問でしょうか?それとも、18問でしょうか?二項分布は、こういった疑問についても統計的に教えてくれます。

例3

 あなたは、人口10万人のA市の保健担当者です。Bという病気(B症)の発症率は人口の5%とされています(話を簡単にするために、年齢性別など関係なく同じ率とします)。A市でもB症の発症率を調べて、一般的な発症率よりも明らかに多ければ、原因を突き止め、予防策などを講じる必要があります。ただ、全住民を調べるには費用と時間がかかります。

 そこで、A市の人口構成と同じような住民100名に協力を得て調査した結果、100名中B症を発症している住民は8名(8%)でした。一般的な発症率に比べて3%多いので、A市のB症発症率は高いので、何らかの対策が必要と結論づけていいでしょうか。二項分布は、こういった疑問についても統計的に教えてくれます。


では、これらの例に答えるために二項分布について詳しく説明していきます。

 二項分布では、期待値において一番起きる確率が高く、そこから離れるに従って起きる可能性が低くなります。期待値とは、コインの表裏や○×であれば0.5ですし、例3のような発症率であればその発症率(例3の場合は、0.05)です。そこを中心としてどのように可能性が低くなるかというと、正規分布曲線に従って左右対称に減少します。

 言葉での説明では分かりにくいので、二項分布の図をRで描きました。描き方は別のところで触れます。

 下の図は、期待値0.5(コイントスや○×での正誤)の時の、50試行(コイントス50回とか50問とか)行った時の確率分布です。横軸(x軸)が起きる回数で縦軸(y軸)がその時の確率です。例えば、25回表がでる確率は赤線と○が交差する点ですので、約0.11となります。表が50回でる確率はほぼゼロで、表が30回出る確率は0.04くらいです。ここで注意してもらいたいのは、この「確率」はある一回の実験ではなく(50回コイントスをする実験を1回の実験と呼びます)、同じ実験を何十回、何百回行った時に、「表が30回出た」という結果になる実験の確率が0.04(4%)ということです。つまり、50回コイントスをする実験を100回行ったら、表が30回出る結果であるのは、4%=4回ということです。こうしてみると、50回投げて、30/50*100=60%表が出ることは、偏りのないコインであれば、4%とかなりレアケースであることが分かり、イカサマコインである疑いが濃厚であると言えます。

 本来であれば、50回のコイントス×100回の実験=5000回のコイントスをしないとかなりレアケースであることが分からないはずが、50回のコイントス×1回の実験=50回のコイントスで統計学的にかなりレアケースであることが分かる(推測できる)のです。

 ここで、お気づきの方もいるかもしれませんが、上の説明でかなりしつこく"50回"と試行回数を連呼していましたが、これには訳があります。

 下の図を見てください。青点は、上図と同じように50試行の時の二項分布です。グラフ左に寄っている黒点の分布は10試行の時の二項分布を表しています。

 x軸に平行な黒線は、確率0.04のラインを示しています。50試行では、50試行中30回つまり60%ですが、10試行では、10試行中8回つまり80%の辺りに線が引かれています。つまり、試行回数が少ない実験では、より回数多く表が出ないと、「表が出やすい」と結論づけられないということになります。

 このように一回の実験における試行回数が多ければ多いほど、期待値からのずれが少なくなるということが言えます。これを大数の法則といいます。

 一方、試行回数を多くすればするほど、実験に時間がかかります。コイントスのように自分だけが大変な思いをするのであれば、許容できるかもしれませんが、被験者に判断させるような課題ですと、実験時間が数時間に及ぶような実験は現実的ではありません。これについては、また別の場所で詳しく解説したいと思います。

 ここでは、例1に基づいて説明をしたので、6-1では例2、3について説明します。

Rにおいては、データを入れる箱に好きな名前(今回の例では"dummy")を付けることができました。データが箱に入っているおかげで、データの数値を一括して処理することができます。

例えば、加減乗除を行えます。

"KanaM_Hit"のデータ全てに"1"を加算するには、KanaM_Hit+1と入力し[Enter]

下図のようにたちどころに、答えが表示されます。

なお答えの各列の頭にある[1] [17]は、その右隣りが何番目のデータであるかを表しています。

掛け算は"*"(アスタリスク)、割算は"/"(スラッシュ)を使います。

これは、KanjiKanaK_Hitのデータに10を掛けて、5で割るという計算です。

なお、先ほどの計算結果の表示と違って[1] [12] [23]とデータの順番表示が多いのは、Console画面の幅に合わせて自動的に改行されるためです。


変数同士を加減乗除することも可能です。この場合、同じID同士の変数が操作されます。

 例えば、今Rの中にデータとして保持されているのは、Hit率とFalse Alarm率だけですが、正答率(Hit率+Correct Rejection率)×0.5がデータとして欲しいなどの場合です。

 Correct Rejection率=1-False Alarm率ですから、上記の計算をRで行わせることができます。こんな感じです。

 計算の優先順位は、通常の数学での優先順位と同じです。ただし二重かっこの場合は数学の時のように{KanaM_Hit+(1-KanaM_FA)}*0.5のように{ }ではなく普通の( )です。

 また、このままでは、正答率をグラフに描きたいとか、平均値を出したいといったときに、毎回上記の計算式を書く必要があります。それを避ける方法として代入 <-  を使います。これは、データ全体を箱に入れた時に行った処理と同じで、計算した結果を入れる新たな箱を定義してあげるわけです。

 例えばこんな感じです。

"KanaM_Corr"という新たな箱を作り、そこに正答率のデータを入れたという処理です。こうすることで、いつでもKanaM_Corrという箱の名前でデータを呼び出して下図のように処理することができました。

①ヒストグラムの作成

 ヒストグラムの作成は、hist()コマンドを使います。

 ()内にヒストグラムを作成したい変数を入力します。

 hits(KanaKanjiL_Hit)のように・・・

(ここからは、attach()により参照データを明示しなくてもいい方法で説明しています)

 [Enter]すると、すぐにグラフが描画されます。

目盛りの間隔を変えたり、バーの色を変えたりできますが、別のところで触れます。


②2変数をそれぞれXY軸に配置して、散布図を描きたいときはplot()を使います。

plot(x軸, y軸)のように・・・

plot(KanjiKanaL_Hit, KanjiKanaM_Hit)     [Enter]とすると

ただちに、散布図が描画されます。

③箱ひげ図を描きたい場合は、boxplot()を使います。

変数は","カンマで区切ることで何変数でも並べられます。

boxplot(KanaM_Hit, KanaM_FA, KanjiKanaM_Hit, KanjiKanaM_FA, KanjiKanaL_Hit, KanjiKanaL_FA)  [Enter]によりただちに、下図のような箱ひげ図が描画されます。

変数名が数字ですが、()内の変数の順と同じ順番で並びます。

番号を変数名などに変換する方法もありますが、別のところで触れます。

Rのコマンドでは、たとえ、読み込まれたデータが一種類だけであっても、どのデータを参照するか毎回宣言する必要があります。

例えばこんな感じ・・・

mean(dummy$KanaM_Hit)

 しかし、dummy$の部分を省略する方法が用意されています。それがattach()です。

こんな感じに使います。

 一回、attach()を使うと、次にdetach()を使って解除するまで、参照するデータを毎回入力する必要がなくなります。

 dettach()を使って解除するとその後からは、参照するデータを宣言しないとエラーになります。

本格的な統計解析に入る前に、データの操作に慣れておきましょう。

今、読み込まれている、各変数の数値データの平均値や標準偏差を算出したい場合

①平均値

 dummyデータの中の変数KanaM_Hitの平均値は、コマンド mean()を使います。

 mean ( dummy $ KanaM_Hit)   [Enter]

    "dummy"  データが入っている箱の名前

    "$ KanaM_Hit"  dummyの箱から変数"KanaM_Hit"のデータを参照する。

     ($マークは必須)

実行すると下図のように結果がでます。

練習①

dummyデータの"KanjiKanaL_Hit"の平均値を求めましょう。

   答)0.68224

②標準偏差

コマンド sd()を使います。

()内は、①の時と同じように書きます。

例えば、"KanjiKanaM_Hit"の標準偏差を求めたければ

sd(dummy$KanjiKanaM_Hit)のように・・・・

番外編

ここまで、きて毎回dummy$と入力するのは面倒くさいと感じませんか?

Rには、これを省略できる便利なコマンドが用意されています。

それが、attach()コマンドです。その説明はこちらで


③その他

その他によく使うコマンドとして

var() 分散

median() 中央値

quantile() 四分位数

などがあり、()内の書き方は同じです。


また、summary()コマンドでデータの要約をすることもできます。

こんな感じで各変数についての要約が表示されます。

Min(最小値)、1st Qu. (第一四分位)、Median(中央値)、Mean(平均値)、3rd Qu. (第三四分位)、Max. (最大値)です。

 データができcsv形式で保存ができたら、いよいよRを使っての作業になります。

 これからは、1-1で説明した練習用ダミーデータ[D:/Rdata/test_data.csv]を使って説明をしていきます。

 ※ダミーデータは、Dドライブに作られた[Rdata]フォルダの中に[test_data.csv]というファイル名で保存されていることとします。

①Rのアイコンをクリックし、起動します。

 こんな画面が立ち上がったと思います。

内側の四角い領域が、Consoleでこの中で作業を行います。

②データの読み込み

 まずは、これから解析をしたいデータを読み込みます。(今回は、test_data.csv)

 読み込む前に、このデータを入れる箱の呼び名を決めましょう。

 呼び名は、何でもいいです(x, test, data, dummy などなど)が、自分自身が分かりやすく、英数半角文字と_で構成されているものを推奨します。

 ここでは、"dummy"と名付けます。

 そうしたら、consoleのカーソルに、

     dummy <- read.table ("D:/Rdata/test_data.csv", sep=",", header=T)と入力し[Enter]

 それぞれの意味や注意点は、下図を見てください。

 これで、一見、何も起きませんが、データはRの中に読み込まれました。

  なお、見やすさのために、ところどころ半角スペースを入れていますが、Rにとっては、これらのスペースはあってもなくても同じです。(詰めて入力しても問題ない)

本当に、読み込まれたかを確認してみます。

consoleに箱の名前である"dummy"と入力し[Enter]

このように、読み込まれたデータ全て表示されます。

次に、このデータの構造を確認しましょう。

コマンドは、str()です。()内に箱の名前を入れます。

今回の例では  str(dummy)となります。

25の被験者(obj)と7種の変数(variables)で構成されているデータで、ID変数はfactorでそれ以外の変数は数値(num)であることが確認できます。

 手元に具体的なデータがない人のために、練習用ダミーデータを載せておきます。

 被験者25名(ID=fa1**で表記)に、3種類の画像(KanaM, KanjiKanaM, KanjiKanaL)を見てもらい、その中に目的の画像が「ある」か「ない」かを判断させる課題を行った結果のデータです。

 目的の画像の「あり」「なし」という実際に被験者が見た刺激が2種類あり、被験者の反応が2種類あるので、刺激-反応の組み合わせは、2×2=4通りあります。それを信号検出理論の用語で記載し表に整理すると下表のようになります。

 下のデータは、上記の3種類の画像に対する、被験者の「ある」「なし」判断を、Hit, Miss, False Alarm, Correct Rejectionに分け、その内、Hit率(_Hit)とFalse Alarm率(_FA)を被験者別に入力したデータです。

 信号検出理論に関連した諸指標を計算するためには、Miss率, Correct Rejection率は必要ないため、記載されていませんが、Miss率=1-Hit率、Correct Rejection率=1-False Alarm率で計算ができます。

 また、3種類の各画像における全体の正答率(つまり、「あり」画像の時に「あり」と判断し「なし」画像の時に「なし」と判断した率の合計)を算出する場合には、以下の計算をすることにより算出可能です。

正答率={Hit率+(1-False Alarm率)}*0.5

 これら計算で算出可能な数値は、被験者一人ずつ計算をすることももちろんできますが、Rを使うことで、全被験者の計算を一括してすることができます。このデータでは25人の被験者を扱っているので、一人ずつ計算も手間ではありますが、実行不可能なことではありません、しかし、100名規模のデータを扱う場合は、一括して処理できることは大変なメリットになります。


----以下をコピーし、メモ帳などに貼りつけ、ファイル名を適当につけてcsv形式で保存----

ID,KanaM_Hit,KanaM_FA,KanjiKanaM_Hit,KanjiKanaM_FA,KanjiKanaL_Hit,KanjiKanaL_FA

fa101,0.999 ,0.001 ,0.999 ,0.001 ,0.556 ,0.444

fa102,0.999 ,0.001 ,0.999 ,0.001 ,0.833 ,0.167

fa103,0.999 ,0.001 ,0.944 ,0.056 ,0.778 ,0.222

fa104,0.999 ,0.001 ,0.999 ,0.001 ,0.722 ,0.278

fa105,0.944 ,0.056 ,0.889 ,0.111 ,0.611 ,0.389

fa106,0.944 ,0.056 ,0.833 ,0.167 ,0.667 ,0.333

fa107,0.999 ,0.001 ,0.944 ,0.056 ,0.611 ,0.389

fa108,0.778 ,0.222 ,0.833 ,0.167 ,0.722 ,0.278

fa109,0.944 ,0.056 ,0.944 ,0.056 ,0.389 ,0.611

fa110,0.833 ,0.167 ,0.778 ,0.222 ,0.611 ,0.389

fa111,0.999 ,0.001 ,0.778 ,0.222 ,0.667 ,0.333

fa112,0.944 ,0.056 ,0.833 ,0.167 ,0.333 ,0.667

fa113,0.944 ,0.056 ,0.944 ,0.056 ,0.833 ,0.167

fa114,0.944 ,0.056 ,0.778 ,0.222 ,0.722 ,0.278

fa115,0.999 ,0.001 ,0.889 ,0.111 ,0.722 ,0.278

fa116,0.999 ,0.001 ,0.833 ,0.167 ,0.889 ,0.111

fa117,0.999 ,0.001 ,0.778 ,0.222 ,0.889 ,0.111

fa118,0.999 ,0.001 ,0.833 ,0.167 ,0.611 ,0.389

fa119,0.889 ,0.111 ,0.889 ,0.111 ,0.556 ,0.444

fa120,0.778 ,0.222 ,0.889 ,0.111 ,0.556 ,0.444

fa121,0.999 ,0.001 ,0.944 ,0.056 ,0.556 ,0.444

fa122,0.889 ,0.111 ,0.833 ,0.167 ,0.778 ,0.222

fa123,0.778 ,0.222 ,0.833 ,0.167 ,0.889 ,0.111

fa124,0.999 ,0.001 ,0.778 ,0.167 ,0.722 ,0.278

fa125,0.889 ,0.111 ,0.833 ,0.167 ,0.833 ,0.167 


 多くの場合、excelなどの表計算ソフトにデータを入力していくと思いますので、excelにデータを入力する場合で説明します。

 右図の例では、被験者が25名いて、データの種類(変数)が被験者No(ID)も含め7種あるとします。

 データ名は、列(横方向)、被験者は、行(縦方向)に入力します。データ名などは、英数半角でハイフンは使わずに、アンダーバー"_"を使うと、問題が起きにくいです。(使えない訳ではないようですが、私は使ったことがありません)

 入力が終わったら、csv形式で保存します。ファイル名も英数半角を使う方がいいです。

例えば[test_data.csv]

 ファイルの置き場は、後のコマンド入力の手間を省くために、CドライブやDドライブ直下に[Rdata]などのフォルダを作ってそこに入れるようにすると便利です。