【R言語】RStudioでt検定　95%信頼区間を可視化する

こんにちは。プログラミング超初心者のえいこです。

前回はRStudioを使って、二標本のt検定をしてみました。

今回はそのおまけで、t検定で「平均値に差がある」と言った根拠である95%信頼区間がどれくらい違うのか？RStudioを使って可視化してみようと思います。

Excelを使っていたらここまで突っ込んでデータを見たりしないとは思いますが、勉強がてらやってみようと思います。

使ったデータはX投与群とプラセボ投与群のタンパク質Aの発現比較のデータです。

95%信頼区間をグラフにする

今回作ったグラフがこちら。

横軸にタンパク質Aの発現レベルを、縦軸にそれぞれプラセボ群とX投与群を取って、95%信頼区間がどれくらいの幅を持っているものなのかを示しています。

これがRではサクッと作れてしまうんだからすごい！

どのように作ったのか、早速コードを見ながら振り返ってみます。

> #各群の95%信頼区間を比較するグラフをつくる
> cixtreat=t.test(xtreat)$conf.int
> ciplacebo=t.test(placebo)$conf.int
> dotchart(c(mean(xtreat),mean(placebo)),pch = 16,
　xlim=range(cixtreat,ciplacebo),xlab="タンパク質Aの発現(a.u.)")
> arrows(cixtreat[1],1,cixtreat[2],1,length=0.05,angle=90,code=3)
> arrows(ciplacebo[1],2,ciplacebo[2],2,length=0.05,angle=90,code=3)
> mtext(c("Xtreat","placebo"),side=2,at=1:2,line=0.5,las=1,)

コードの行数が多くなっています。

一つずつ見ていきます。

95%信頼区間の定義づけ

コードの一行目と二行目はX投与群とプラセボ群の95%信頼区間をそれぞれ"cixtreat"と"ciplacebo"という変数で定義しています。

> cixtreat=t.test(xtreat)$conf.int
> ciplacebo=t.test(placebo)$conf.int

一行目と二行目はほぼ同じことをしているので、二行目の説明は割愛します。

あれ？変数に"代入"するときに使うのって"<-"では？？

演算子"="と"<-"の違いは、正直なところよくわかりません。おそらく"<-"でもうまくいくとは思いますが、今回参照したサイトでは"="を採用していたのでこちらを採用します。

追記："="の代わりに"<-"を使用しても問題なくグラフをかけました。

次に、"t.test(xtreat)"では一標本のt検定をそれぞれについて行っています。ここでの検定は、「真の平均値が０と等しい」という帰無仮説について検定しています。

どういう検定なのか一応、結果を出しておきます↓

> t.test(xtreat)

	One Sample t-test

data:  xtreat
t = 21.994, df = 2, p-value = 0.002061
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
4.595655 6.831011
sample estimates:
mean of x 
5.713333 

その後ろに"$conf.int"と書いています。

conf.intはconfidence intervalつまり、信頼区間のことです。

オブジェクト（今回の場合は"t.test()"）が持っているデータを、"$"を使って指定することができます。

今回の場合は、"t.test(xtreat)"が持っている信頼区間"conf.int"を"cixtreat"に代入するということをしています。

ドットチャートをかく

次に、定義した数値を元にグラフをかいています。

三行目のコードはこんな感じでした。

> dotchart(c(mean(xtreat),mean(placebo)),pch = 16,
　xlim=range(cixtreat,ciplacebo),xlab="タンパク質Aの発現(a.u.)")

"dotchart()"の中身は、データのベクトル、プロットのマーカー、X軸のプロットの範囲、X軸のラベルという順番に並んでいます。

どうやらドットチャートでは、グループごとにY軸の高さを変えて結果を表示してくれるようです。（ドットチャート素人...）

なので、今回は"xtreat"と"placebo"という二つのグループがあることをベクトルを使って言っていて、"mean()"を使ってそれぞれのグループの平均値をプロットするようにしています。

"pch"を使ってプロットのマーカーを指定しています。マーカーは全部で26種類で、0から25の数値で指定できます。どんなマーカーがあるのか試してみても楽しいかもしれません。（私はgoogle先生に聞いてしまいましたが...）

"xlim"を使ってx軸のプロットの範囲を指定します。今回は"range()"を使って範囲を指定していますが、具体的に数値を使って範囲を決めることができます。

最後に"xlab"を使ってX軸のラベルを書いています。

最後の二つは"ylim""ylab"という形で箱ひげグラフをかいたときにも出てきました。

ここまでを実行するとこんなグラフがかけます。

ちなみに"pch=0"にしてかいてみると...

プロットのマーカーが変わりました！

95%信頼区間の幅を表示させる

ドットが書けたので、95%信頼区間の幅を"arrow()"関数を使って書きます。

arrow関数は名前の通り図の中に矢印を書き加えるための関数です。矢印の「始点と終点の座標」と「矢印の形状」を指定することでエラーバーなども書くことができます。

今回書いたコードはこちら。

> arrows(cixtreat[1],1,cixtreat[2],1,length=0.05,angle=90,code=3)

"citreat[1],1"で矢印の始点のX座標とY座標を指定しています。

"citreat"という変数の中に95%信頼区間の値が[4.595655,6.831011]という一行2列のベクトルの形で格納されています。

ベクトルの要素を取り出すのには[]を使うんでした！（【R言語】初めての変数・ベクトルを参照）

ということで、"citreat[1],1"は「"citreat"変数のなかの一番目の要素をx、1をyとする座標を矢印の始点とする」という意味です。

同じように"citreat[2],1"で終点の座標を指定しました。

次の"length=0.05"では、矢じりのサイズを指定しています。大きくしたい場合は個々の数値をいじればOKです。

"angle=90"と"code=3"はエラーバーを書くときの決まりの数値のようです。angleの値を変えると矢印っぽくなるようです。（angle=90というのは矢じりの角度が90°ってこと？）

codeの値を変えると、矢印の方向を「始点方向のみ、終点方向のみ、両方」に変えることができます。基本的には両方しか使わないので"code=3"で覚えてしまいます。

ここまでのコードを実行すると、

ここまでできます。

最後にY軸のラベルを書いて完成

このままでは、どっちが薬剤投与群でプラセボ投与群なのかわかりません。

文字を書き込む場合は"mtext()"関数を使うようです。実際に書いたコードはこちらです。

> mtext(c("Xtreat","placebo"),side=2,at=1:2,line=0.5,las=1,)

まず、書き込む文字列を指定します。今回はベクトルを使って指定しています。（複数の要素を扱うときはベクトルが便利ですね）

次に文字列を書き込む位置を"side"で指定します。sideは、1が下、2が左、3が上、4が右で指定できます。今回は左側なので"side=2"で指定しています。

atで文字列を書き込む座標を指定できます。at=1:2で二行に分けるように指定しています。

次のlineで、図形から何行離すかを指定して、"las"で文字列の方向を指定しています。

las=0で各軸に並行して、las=1ですべて水平に、las=2は軸に対して垂直に、las=3ですべて垂直に書くとなります。今回はすべて水平になので"las=1"で指定しました。

これで、書いたコードの説明は終わり！

おまけのつもりで書いたけど、一番長いnoteになってしまった...

それでは、また！