統計解析ログ/R作図 - 引越作業中

R作図

このページでは、Rで作図する基本的な手順を紹介したいと思います。いい図はデータ解析を助けてくれますので、ぜひ自分の思い通りの図を描けるようになってください。

Rで作図するメリット・デメリット

メリット

• データを解析しながら作図できる
• 複数の図を描くときに自動化できる
• 様々な形式で図が出力可能（png, jpg, メタファイル, pdf, epsなどなど）
• 低水準作図を駆使すれば、フォーマットが決まっていない図も描ける（地図など）

デメリット

• マウスで何か操作する訳ではないので、関数を知らないとどうにもならない
• 複雑な作図だとコードが長くなる
• 低水準作図で次々と要素を足していったときに、途中で間違えても戻れない

個人的には、複数の図が自動化ですぐにかける点と、低水準作図を駆使することでフォーマットにとらわれない図がかける点が特に便利だと思います。これらの点は、データを見る上で（ひいては統計解析する上で）データをよく見るために欠かせません。

ということで、私も昔はExcelで作図していましたが、Rで作図する方が研究にとっては100倍ぐらい楽だと思っています。

なお、筆者の力量不足のため、このページではlatticeによる作図は取り扱いません。

Rにおける作図手順

1. どんな図を描くかよく考えます。既成のフォーマットにとらわれる必要はありません。ノートにラフスケッチするのもいいでしょう。
2. 図を描くにあたり、データの形（並び方）が作図に適しているか考えます。何かしらの要因が横方向に展開しているデータの形式は好ましくありません。また、自動でいくつもの図を描くのであれば、split()などでデータを分割しておく必要があるでしょう。
3. 出来上がる図が、高水準作図で用意されている図で作成できそうか考えます。どうにも当てはまらない場合、高水準作図で空の空間だけ発生させます。
4. 低水準作図で細かいところまで仕上げます。

とにかく大事なのは、どういう図を描きたいのかよく考えることです。そしてその図は、既に見たことがある図でもいいですし、全く新しい図でもかまわないのです。自分が言いたいことをもっともよく表現してくれそうな図がどんな図か、よく考えましょう。

高水準作図

高水準作図は、「すでに形ができあがっている図を出力する作図」です。図の種類としては、

• plot(): 散布図
• pairs(): 対散布図（全ての2変数の組み合わせについて散布図を作成）
• boxplot(): 箱ひげ図（カテゴリーごとのデータのばらつきを見るのに有効）
• barplot(): 棒グラフ。Rコアチームの意図か、デフォルトではひどく見栄えの悪い関数
• pie(): 円グラフ
• scatterplot3d(): 三次元散布図。library(scatterplot3d)で呼び出せる

などがあります。

低水準作図関数

低水準作図は、図に細かく装飾をしていく作図です。これを駆使すると、一から新しい図を作図することもできます。

凡例を表示

legend()を使う。

例えば

> d <- rnorm(100, 10, 5); e <- d + rnorm(100, 5, 2)
> plot(y ~ x) #まず図の生成
> legend("topleft", legend="x") #凡例を図に表示させる

legend("出力位置", legend="つけたい名前")という指定をしている。

2組の凡例（例えばxとy）を同時発生させる場合は、

> plot(y ~ x)
> legend("topleft", legend=c("x", "y"))

とする。単に名前の部分がベクトルで複数になっただけ。

さらに、2組の凡例で、かつ（シンボル テキスト）というような見せ方をしたい場合。

> plot(y ~ x)
> legend("topleft", legend=c("x", "y"), pch=c(15,19), col=c("yellow","purple"))

legend("出力位置", legend=c("xのラベル","yのラベル"), pch=c(xのシンボル,yのシンボル), col=c("xのシンボルの色","yのシンボルの色"))という指定をしている。

デフォルトでは凡例に枠がついてしまうが、これを消したい場合は、legend()の引数で

> legend(......, bty="n")

とすればよい。

特殊なフォントを使う

expression()を使う。

例えばplot()中でイタリック文字を使いたいとき、

> plot(y ~ x, xlab=expression(italic("イタリック文字ににしたい文字")))

とすればよい。

数式も記述できる。

> plot(y ~ x, xlab=expression(paste("RGR (mg ",mg^-1,year^-1,")"), ylab="")

expression内のpaste内で記述する。普通の文字にしたい部分は""でかこみ、数式処理をしたいところはそのまま記述する。今回は^で上付き処理をしている（下付きは[]）。通常部分と数式部分は,で区切る。

棒グラフにエラーバーを発生させる

• barplotとarrows関数を使う
• barplot2関数を使う

barplot()の場合

> y <- c(1,3,11,8,15) #値のベクトル
> error <- c(0.1,0.3,1.1,0.8,1.5) #エラーバーの長さ（標準偏差とか）
> location <- (0.2 + 1) * 1:5 - 1/2 #エラーバーを発生させる位置（x座標）。

最後のエラーバーの位置指定は一見よくわからない式だが、Rの棒グラフは棒グラフ間は0.2、棒グラフの太さは1がデフォルトである。棒グラフの中心にエラーバーを発生させるのであれば、x座標は最初必ず0.2あいて、それに棒グラフ1本が来て1が足され、そこから半分の0.5が引かれている。そして1:5で2本目、3本目...とどんどん位置を指定している。ちなみに名前は何でもよい。

というように、値そのもの、エラーバーの長さ、エラーバーのx軸の位置を指定するオブジェクトを生成しておいて、

> barplot(y) #ここでまず棒グラフが生成される
> arrows(location, y, location, y + error, angle=90)

（エラーバーの基部のx座標,基部のy座標,エラーバーの先端のx座標,先端のy座標,エラーバーの角度）という指定の仕方

これで棒グラフに上向きのエラーバーが生成される。もし下向きのも付け加えたいときは

> arrows(location, y, location, y - error, angle=90)

とすればよい。

barplot2()の場合

barplot2はgplotsというパッケージ内にある。

> y <- c(1,3,11,8,15)
> error <- c(0.1,0.3,1.1,0.8,1.5)
> barplot2(y, plot.ci=T, ci.l=y-error, ci.u=y+error)

というように、barplot2 (値そのもの、エラーバーを作るか作らないか（Tで作るにしている)、エラーバーの下限値、エラーバーの上限値）
とすればよい。

散布図に直線を当てはめる

lsfit()で回帰直線を生成し、abline()で散布図に結合する。
例えば

> plot(x,y)  #まず散布図を生成
> abline(lsfit(data$x, data$y)) #lsfit内で回帰直線を生成、#ablineで散布図に結合

とすればよい。

または、

> model <- lm(y ~ x) #回帰分析の結果をmodelに結合
> plot(x,y) #散布図の生成
> abline(model) #modelに結合された回帰直線をablineによって散布図に結合

対数軸で自由に軸幅設定

axis()を使う。

通常図の軸幅は

> par(lab=c(*, *, 7))

によって指定する。前の二つで、x軸とy軸につける目盛の数を指定できる。最後の数字はラベルのサイズを指定するはずなのだが、どうやらRには移植されていないコマンドらしい。

ところが対数軸にすると、このグラフィックスパラメータは機能しない。そこで以下のようにする。

> plot(x, y, xaxt="n", ...) #xaxt="n"で軸を書かないようにする。
> axis(1, c(0.1,1.0,10.0,100.0)) #自分で軸を設定。

axis(軸の場所（1は図の下に軸をおけという意味（3が上、右が4、左が2））, c(に区切りを挿入する場所))である。

軸目盛ラベルの角度変更

軸目盛ラベルの角度は通常グラフィックスパラメータlas=*を高水準作図関数内の引数として記述する。

• las=0 #各軸（x,y）に対して平行にラベルが出現
• las=1 #軸に関わらず全て水平にラベルが出現
• las=2 #各軸に対して垂直にラベルが出現
• las=3 #軸に関わらず全て垂直にラベルが出現

さらに、角度を詳細に指定して描きたい場合がある。この場合、一度ラベルなしの軸を描き、その上に低水準作図関数text()を利用してラベルを挿入する。ただしtext()で今回のように欄外に文字を発生させる場合、あらかじめ欄外に文字を書いてもよいという指定をする必要がある。以下ではboxplotでの作成例。

> par(xpd=T) #図の外に文字を書くことを許可する指定
> boxplot(column1,column2,column3, xaxt="n") #xaxt="n"で軸を書かないように
> axis(1,at=1:*,labels=F) # *は欲しい目盛の数（何個目盛を打つか）、labels=Fは目盛ラベルを表示しない指定
> x <- 1:* # 目盛ラベルのx座標の指定、*は欲しい目盛の数
> y <- rep(-5,*) #目盛ラベルのy座標の指定、-5はグラフの値によって適時修正、*は欲しい目盛の数
> text(x,y,c("A","B",...), srt=**, adj=1)

text()内の説明：x,yはそれぞれ先ほどの位置を指定する変数、c()は実際に入力したい目盛ラベルの名称が入ったベクトル、srt=**が本題の、目盛ラベルの角度を指定するパラメータ、adj=1はラベルを右揃えにするパラメータ（ラベルの角度によってadj=0左揃え、adj=0.5中央揃えのほうが見栄えがよいときがあります）

カテゴリーごとにシンボルの異なる図を描く

カテゴリーを数字として扱うこと、数字で色やシンボルの形が指定できるのがポイントです。具体例を以下で示します。

> e <- rnorm(100, 10, 5) #架空データの生成
> f <- e + rnorm(100, 5, 2) #同上
> g <- rep(LETTERS[1:4], each=25) #同上（カテゴリーデータ）。
#LETTERSはRに最初から用意されている文字列ベクトル
> d <- data.frame(e, f, g) #データフレーム化
> plot(f ~ e, d) #普通に描画すると、全てのカテゴリーが同じシンボルになってしまう

> plot(f ~ e, pch=as.numeric(g), d) #カテゴリーごとにシンボルが異なる

as.numeric()を使いこなす

as.numeric()は文字列だろうが無理矢理数字にしてしまう関数です（数字はアルファベットの速い順から1、2、......とつきます）。これをカテゴリー変数にすることで、各カテゴリーごとに異なる数字がシンボルの形を示す引数pchや色を示す引数colに投入されます。

なので、色も同様に変えられます

> plot(f ~ e, pch=as.numeric(g), col=as.numeric(g), d)

描画するシンボルや色を指定したい！

as.numeric()を使うと確かにデータにあわせてシンボルや色は変えられますが、基本的には1から順に番号が振られるので、必ずしも思い通りのシンボルや色になりません。こういう時には、事前にシンボルや色を指定しておきます。

> iro <- c("black", "red", "blue", "green")
> katachi <- c(0, 1, 15, 16) #番号によってどのようなシンボルが出力されるかは
#library(Hmisc) ; show.pch()を見てください。
> plot(f ~ e, pch=katachi[as.numeric(g)], col=iro[as.numeric(g)], d)

さて、上記のコマンドは何を意味しているのでしょうか？ベクトルの性質を思い出してみましょう。

ベクトルはある一つながりのデータです。先ほど作ったiroやkatachiもベクトルです。ベクトルの特定の要素（何番目）を取り出すには、[]を使うことができました。例えば、

> iro[3]
[1] "blue"

というように、上記の通りiroの3番目は"blue"でしたが、そのiroの3番目の要素を取り出しなさいということを、iro[3]で実行したわけです。

ここで上記の作図コマンドに戻ります。iroやkatachiの[]内にはいずれもas.numeric(g)が入ってますね。as.numeric(g)は前述の通り、カテゴリーを数字化したものです。すなわち、アルファベットで最初にくるカテゴリーについては1が、アルファベット順で2番目のカテゴリーについては2が、......とそれぞれ値が返ります。

これによって、iroあるいはkatachiの何番目の要素を取り出すかを、as.numeric()でカテゴリーごとに指定することで、カテゴリーごとに別な色やシンボルを表示させています。

これらの性質を考慮すると、例えば1番目のカテゴリーの色は何色でシンボルはこの形、という希望があれば、iroやkatachiの1番目の要素をその希望する内容に変えてしまえばいいのです。これによって、カテゴリーごとに色やシンボルを変更でき、かつ希望のパタンに変更する、ということが可能になります。

複数の図を自動で描く

先ほどは異なるカテゴリーのデータを同一の図に描画しました。今度は、カテゴリーごとに図を描きます。

といっても、一つ一つの図のコマンドを書いていたのではコードが膨大になり、また間違いが発覚した際に直す箇所も複数生じてしまいます。こういう時は、for()を使って自動化することを考えます。

このように、複数の図を自動で描く場合の基本的な手順は以下のようになります。

• 図を何行×何列で表示させるのか決める
• データフレームをsplit()でカテゴリーごとに分割する
• for()で、カテゴリーごとの図の描画を全カテゴリー分繰り返す

では、以下で実際例を見ていきましょう。

図のレイアウト

図を何行×何列で表示させるのか指定する方法としては、

• par(mfrow=c(行数, 列数))
• layout()

の2つの方法があります。図を描く順番等を特に考える必要がないなら前者でいいでしょう。

layout()で自由な図の配置

layout()は、何番目の図をどこに描くかを数字で指定します。そのため、例えばある図だけ2行にわたって表示させたり、4カテゴリー目の図を1番最初に表示させたいなど、不規則な図の表示をさせるときに使います。例えば、

1番目の図   4番目の図
2番目の図   3番目の図

というような配置にしたければ、

> layout(matrix(c(1, 4, 2, 3), ncoL=2, byrow=TRUE))

と指定します。matrix()はその名の通り行列を作る関数で、ncolは行列の列数、byrow=TRUEにすると、c()で投入した文字が行方向から埋めるように指定することになります。

別の例を見てみましょう。

1番目の図   4番目の図   3番目の図
2        番        目        の        図

というような配置にしたければ、

> layout(matrix(c(1, 4, 3, 2, 2, 2), ncol=3, byrow=TRUE))

と指定します。2番目の図が他の図3つ分のスペースに広がって表示させるために、2を3回指定しています。

相性抜群のsplit()とfor()

図のレイアウトが決まったら、カテゴリーごとに図を描画します。その際の基本的な考え方は、

データをカテゴリーごとに分割し、カテゴリーごとのデータフレームを使って
図を描く作業を全カテゴリー分繰り返す

となります。

split()の動作

split()は、データフレームをあるカテゴリーごとに分割することができます。実行方法は、

> d2 <- split(d, d$g)

となります。では、実際にはどのように分割され、またどのようにして取り出すことができるのでしょうか？

今回のカテゴリーですが、

> levels(d$g)
[1] "A" "B" "C" "D"

という4カテゴリーです。ちなみに、

> nlevels(d$g)
[1] 4

とすると、カテゴリー数を取り出すことができます。

つまり、データは4つに分割されているはずです。この、分割されたデータを取り出す方法は2つあります。

• カテゴリー名を指定する
• カテゴリーの順番を番号で指定する

実際に見てみましょう。

#カテゴリー名指定
> d2$A

#カテゴリーの順番を番号で指定
> d2[[1]]

この、番号でも分割したデータをカテゴリーごとに取り出せる点が、繰り返し動作をさせる関数for()との相性の上で重要です。

for()の動作

for()は、ある動作を繰り返し実行させる関数です。最初から言葉で説明するとわかりにくいと思うので、実行例を見てから説明したいと思います。

> par(mfrow=c(2,2)) #図を2×2で配置すると指定
> for (i in 1:4) {
+      plot(f ~ e, d2[[i]])
+ }

すると、図が4つ生成されたと思います。これを順に説明します。

まず1行目です。

for (i in 1:4) {

とありますね。

• iは、繰り返しを指定する「文字」です。別にiである必要はなく、どんな文字で構いません。
• inは、先ほどのiという文字に、inの右にある内容を代入する、という意味です。
• 1:4は1,2,3,4という数字の集まりです。
• 一番左の{はこれ以降が繰り返す内容だよ、としきるための括弧です。

つまり、iという文字に1、2、3、4が次々と代入されるという意味になります。

2行目です。

plot(f ~ e, d2[[i]])

とありますね。この行は先ほどの{の内側なので、繰り返される内容です。ポイントはplot()で使うデータフレームで、

d2[[i]]

とあります。ここで1行目で指定したiが出てきましたね。

ここがポイントで、iには1行目で指定した数字が代入されます。すなわち、最初は

d2[[1]]

となり、次は

d2[[2]]

となり、......と、全カテゴリーについて繰り返されることになるのです。

split()で分割したデータフレームは数字でも取り出すことができ、for()によって数字を次々に代入することで、分割した各カテゴリーごとに図を描く作業を、全カテゴリー分繰り返すことが可能になるのです。

3行目です。単に}があるだけですね。これは、繰り返す内容がここで終わりだよ、と示しています。

例外処理

さて、このように複数の図を描くとき、全ての図が同じ形でいいということはあまりありません。例えば左側の図だけy軸を表示したいなど、例外処理が必ず発生します。例外処理は、if (条件式) {動作内容}やif (条件式) {動作内容} else {条件式がFALSEの場合の動作内容}などで制御します。

図全体に対する装飾

完成コード

以上の手順をまとめると、以下のようになります。

#データの準備
> e <- rnorm(100, 10, 5) #架空データの生成
> f <- d + rnorm(100, 5, 2) #同上
> g <- rep(LETTERS[1:4], each=25) #同上（カテゴリーデータ）
> d <- data.frame(e, f, g) #データフレーム化
> d2 <- split(d, d$g) #カテゴリーごとにデータを分割

#図のレイアウト
> par(mfrow=c(2,2))

#forによる自動作図
> for (i in 1:nlevels(d$g)) {
+      plot(f ~ e, d2[[i]])
+      legend("topleft", legend=unique(d2[[i]]$g), bty="n")
+ }

#図全体の装飾
> mtext(1, line=2, outer=TRUE, text="X")
> mtext(2, line=2, outer=TRUE, text="Y")