2012年12月5日水曜日

Rいろは・第三部：グラフ作成編


（cf. Rいろはのトップページ http://nhkuma.blogspot.jp/p/r_5.html）

# （2012.12.27）操作ミスによりアドレスが変わってしまいました…

# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）

# データの準備の仕方からグラフ作成、体裁作りまでを実行します


###### 基本的グラフ関数の紹介 #####################
# まずはシンプルなサンプルデータを使用
x0 <- rep(1:5, each=4) # 1~5を4個ずつ
y0 <- sort(rnorm(20, 10)) # 平均10の乱数を20個、ソートしておく

plot(y0) # 大概のデータは"plot()"に放り込めばとりあえず図示してくれる
boxplot(y0 ~ x0) # 箱ひげ図（~の意味は、y0をx0で分ける感じ → 統計でよく使う）
barplot(y0) # 棒グラフ
hist(y0) # ヒストグラム
plot(y0, type="o") # 折れ線グラフ
plot(y0 ~ x0) # 散布図（plot(x0, y0)と表すこともできる。ただし、先がX軸, 後がY軸になる）
plot(y0 ~ x0, ylim=c(8,12)) # Y軸の上下限を8と12に設定（X軸の場合はxlim）
plot(y0 ~ x0, yaxp=c(8,12,4)) # さらにY軸の目盛数を変更（上下限を合わせて目盛数は+1になる）
plot(y0 ~ x0, log="y") # y軸を対数スケール（log="x": x軸を対数、log="xy": 両対数）
plot(y0 ~ x0, cex=c(1:20)) # 散布図でプロットサイズを変更（cex：サイズ指定）

# 注：統計モデルでよくやるように、plot(y ~ x, data=D)のようにデータフレームDのxとyを用いるように指定できるが、plot(x, y, data=D)とするとデータが見つからずエラーとなるので注意。
# （実は ~ を用いるのはplot.formula()という拡張関数であり、~を用いると自動で呼び出される）

###### Rグラフの取り扱い #################
# "ファイル" → "保存"で取り出せる（PDF形式）
# Win版ではメニューから様々なファイル形式での保存が選べるようになっている

# 保存したグラフをOpen officeで編集するやり方がこちらにあり（Win版の説明）


###### 以下、現実的なデータセットと体裁、グラフ形態にしてみよう ########
## 下記のグラフ作成で用いるサンプルデータを用意する。カテゴリー変数：２要因、連続変数：１要因
f1 <- rep(c("A", "B", "C","D"), each=10)
f2 <- rep(c("a", "b"), 20)
X <- sort(sample(0:100,40,replace=T)) # 0~100の数字をランダムに抜き出しソート（重複あり）
Y <- runif(40,0,100) # random uniformの意味。0~100の一様乱数を40個発生させる
D <- data.frame(X=X, Y=Y, F1=f1, F2=f2) # データフレーム化する

## 棒グラフ、折れ線グラフ用に平均と標準偏差を用意
MEAN1 <- tapply(D$Y, D$F1, mean) # 各グループの平均（"対象データ", "グループ分け", "関数名"）
SD1 <- tapply(D$Y, D$F1, sd) # 各グループの標準偏差
MEAN2 <- tapply(D$Y, list(D$F2, D$F1), mean) # 各グループの平均
SD2 <- tapply(D$Y, list(D$F2, D$F1), sd) # 各グループの標準偏差


###### グラフの体裁を整える準備（グラフ枠とラベルのサイズを固定するなら） ####
# par() に入れたパラメータの効果はウインドウを閉じるまで持続する
par(lwd=2, # line widthの意味。線の太さを標準の2倍にした
cex=1.5, # 文字サイズの指定。標準の1.5倍にした
mgp=c(1.75, 0.75, 0), # 軸タイトル、軸ラベル、軸線の間隔を指定
pin=c(3, 3)) # 作図領域のサイズをインチ指定

# family="Helvetica" # フォントを指定する項を追加可能
# この例は Helvetica に指定した例（推奨）。ほか Courier, Times など
# Windowsでは業界標準の Helvetica がインストールされていないので、擬似フォントの Arial を指定する
# （cf. ご自身で記述する時はこんなにマメに改行する必要はありません）

## モニタの小さなPCでは画面に収まりきらずにエラーになる場合あり。そうしたらpin=c(,)の値を小さくしてください。

###### 箱ひげ図（１要因）の作成例 ###################### 
plot(Y ~ F1, data=D, # YをF1でグループ分け、データはDを用いる、の意味
lwd=2, ylim=c(0,100), # y-limitの意味。上限と下限を指定（計算式も利用可）
col="gray80", # カラーを指定：80%白の灰色。"colors()" で色一覧を見られる
axes=F) # いったん軸を除去（plotは、目盛りの線太だけは変えられないようなので）
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1)
# 軸の付け直し（二度手間だが、これをすることで目盛りも同じ線太にできる）

### さらに、図のラベルを左上部に追加する（論文に掲載する図を想定）
mtext("A) Coral growth", side=3, # side：3で上部を指定
at=0.8, line=0.5, # at：左からの横位置、line：グラフから上端からの縦位置
cex=1.5, font=2) # font=2：太字（普通でいいならfont項は不要）
# atの値はデータに応じ細かく調整要…正直、ドロー系ソフト上で調整する方が手っ取り早い


###### 棒グラフ（１要因の場合）の例 ########################
bar1g <- (barplot(MEAN1, # 図示したいデータ列
ylim=c(0,100), lwd=2, # 線太を2倍にした
las=1, # Y軸のラベルを水平方向に変更
col="wheat", # カラーを指定
xlab="group", ylab="volume")) # ラベル名を自分で書く
### エラーバー追加（barplot自体が各棒の中心のx座標データを持っているので利用する）
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2) # x1,y1,x2,y2：両端のx,y座標を指定
mtext("B) Coral happiness", side=3, at=1.2, line=0.8, cex=1.5, font=2)

# cf. points(x,y) # x,y座標を指定し点をプロットできる（cex=数字などの項を加え装飾可能）


###### 応用：間隔を変えた棒グラフ（control１群に対し、treatment複数群の場合）#
barVSctrl <- (barplot(MEAN1,
space=c(1,1,0.5,0.5), # space：棒の幅1に対する棒の間隔を数値で指定
xlim=c(0,7), # xlim：x軸の範囲。棒の幅、棒の間隔の合計より少し大きくなるように設定
ylim=c(0,100), lwd=2, las=1,
col=c("white", "lightblue", "lightgreen", "mistyrose"),
xlab="group", ylab="amount"))
segments(barVSctrl, MEAN1+SD1, barVSctrl, MEAN1-SD1, lwd=2)
mtext("C) Fish sadness", side=3, at=1.0, line=0.8, cex=1.5, font=2)


####### 棒グラフ（２要因（２x２グループ)）の例 ######################
bar2g <- (barplot(MEAN2, beside=T, # TをFにすると積み重ねグラフになる
ylim=c(0,100), lwd=2, las=1, col=c("white", "salmon"),
xlab="group", ylab="Intensity"))
segments(bar2g, MEAN2+SD2, bar2g, MEAN2-SD2, lwd=2)
mtext("D) Shrimp madness", side=3, at=3.0, line=0.8, cex=1.5, font=2)


###### ヒストグラム ######################
hist(D$Y, breaks=c(0,20,40,60,80,100), # 区切りの指定（しなくてもOK）
main="", # でかでかタイトルを除去
lwd=2, las=1, col="seashell", xlab="velocity")
mtext("E) Mollusca marathon", side=3, at=25, line=0.8, cex=1.5, font=2)


###### 折れ線グラフ ######################
plot(MEAN1, ylim=c(0,100), xlab="factor", ylab="richness",
type="o", # 折れ線グラフ（cf. "p"：点のみ、"l"：線のみ）
lwd=2, axes=F, cex=3, # プロット（文字サイズcex）を倍にする
pch=23, # プロットのシンボルを指定（21~25では輪郭と塗りつぶしの色が指定可能）
col="black", # 輪郭の色
bg="tomato") # 塗りつぶしの色
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
segments(c(1:4), MEAN1+SD1, c(1:4), MEAN1-SD1, lwd=2) # 今度はx軸は連続変数なので数値で指定
mtext("F) Polychaeta beauty", side=3, at=1.7, line=0.8, cex=1.5, font=2)


###### 散布図 + 回帰直線 ######################
plot(Y ~ X, xlim=c(0,100), ylim=c(0,120), xlab="Factor", ylab="Loneliness",
axes=F, cex=2, pch=21, col="black", # シンボルの形、輪郭の色は全系列で統一した
bg=c("seagreen","turquoise","slateblue","plum")[F1], data=D)
# 塗りつぶし色を系列によって変える（"F1"のA~D、4カテゴリーで色分けした）
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
legend("topright", legend=levels(D$F1), pt.cex=2, pch=21, # 凡例を追加
pt.bg=c("seagreen","turquoise","slateblue","plum"),horiz=T)
curve(0.2*x + 30, lwd=4, col="orangered", add=T) # 回帰直線：xの一次関数の例
mtext("G) Damselfish", side=3, at=5, line=0.8, cex=1.5, font=2)


############## 複数の図の重ね合わせ ######################
###### 散布図の複数系列を重ねる（上記とは別のやり方で描いてみる） #########
# まず、枠を用意
matplot(c(0,100), c(0,100), xlab="Factor", ylab="Jealousy", axes=F,
type="n") # type="n"：グラフエリア内の余計な文字を消す
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
# １系列ずつ描く
matpoints(D$X, D$Y, cex=2, pch=21, col="black", bg="turquoise") # １系列目
matpoints(D$X*1.1, D$Y*1.1, cex=2, pch=24, col="black", bg="slateblue") # ２系列目
mtext("H) Sea urchin", side=3, at=5, line=0.8, cex=1.5, font=2)



############## 複数の図の並列 ######################
library(lattice) # latticeというパッケージを呼び出す
xyplot(Y ~ X | F1, cex=1.5, data=D) # ぴったり並列して書き出すので見やすい
# ただし、par()が効いていないようだし、plot()とは使い勝手が異なる模様…便利そうなのですが。


###### 異なるタイプのグラフを重ねる ################
# まず、一つ目のグラフ（棒グラフとする）
bar1g <- (barplot(MEAN1, xlim=c(0,5), # x範囲を二つ目のグラフと合致させる
ylim=c(0,100), lwd=2, las=1, col="khaki",
xlab="group", ylab="conflict"))
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2)
par(new=T) # 上書きするコマンド

# 二つ目のグラフ（折れ線グラフとする）
plot(bar1g, # エラーバーと同じ要領でbar1gをx座標に利用
MEAN1*1.5, xlim=c(0,5), # x範囲を一つ目のグラフと合致させる
type="l", # 線だけのグラフタイプを指定
axes=F, xlab="", ylab="", # 枠やラベルを描かないためのコード
lwd=4, col="royalblue")
axis(side=4, lwd=2, las=1) # 右サイドに軸だけを追加
mtext(expression(paste(factor2(m^2))), side=4, line=2, cex=1.5)
# 二つ目のグラフのラベルを右サイドに追加（expression, pasteを組合わせると上付き等も書ける）
mtext("I) Sea cucumber", side=3, at=0.1, line=0.8, cex=1.5, font=2)


####### 画面を分割して複数グラフをまとめる ########################
# parにmfrowを加える
par(mfrow=c(2,2), # multi frame by row の意味。この例では行方向に2x2=4分割
lwd=2, cex=1, mgp=c(1.75, 0.75, 0), pin=c(1.5, 1.5))
# 適当な4枚のグラフを実行してみよう
#（ウィンドウをドラグしてサイズを調整するとグラフの間隔を変えることができる）

## モニタの小さなPCでは画面に収まりきらずにエラーになる場合あり。そうしたらpin=c(,)の値を小さくしてください。


####### 便利技：全変数を一括で図示してみる ########################
# 散布図形式で、全変数の関係を対戦表に配置する
# （自身vs.自身の位置に変数名が入ります。横方向は自身がY、縦方向は自身がXの場合）
pairs(D, panel=panel.smooth) # panel.smoothで変数間の関係を表す曲線が追加される
pairs(Y ~ X+F1+F2, panel=panel.smooth, data=D) # これでも同じ図になる



####### 便利技：色にこだわる ########################
# 色の種類指定が面倒くさいのは難点。せめて色一覧くらいほしいですね。
# htmlの色コードが使用可能：
# plotなどの関数の項で、col="#00FFFF" のように""でコードを囲うと利用できる（水色の例）。

# ちなみに、RGBのhtml色コード末に AA を付加すると半透明になる（例： col="#00FFFFAA"）
# cf. 透明度10%の場合、#00FFFF1A。命名ルールは今一つよくわからない…

# なお、グレースケールでよければ、gray()と0~1の小数を用いてgray(0) ~ gray(1)で表現可能
# 0が黒、1が白。なお、この場合はgray() の前後に""は付けない（関数になっている）。



# その他、実用的な色セットを用途別に一覧にしてくれるパッケージがある（インストール必要）。
install.packages("RColorBrewer", dependencies=T)  # RColorBrewer のインストール（初回のみ必要なコード）
library(RColorBrewer)  # RColorBrewerの呼び出し（毎回必要。名前が長い…忘れてしまうorz）
display.brewer.all()  # 色一覧をグラフにしてくれる（増加、カテゴリー、正負の増減、用途別）
# 各行の文字は色セット名を表しており、この単位で色を利用できる
brewer.pal(8, "Greens") # Greensを選択する場合。数字は色数（段階数）
# ただし、各セットに用意されている色数を超えて指定することはできない（cf. 上図）
brewer.pal(8, "Greens")[1] # 実際の使用では[1]のように添字を付けて色を取り出す
paste(brewer.pal(8, "Greens"), "AA", sep="") # こうすると透過色にできる

# ColorBrewerのサイトも有用です、いろいろなケースにおいて勧められる色セットなどのガイドもあり

# cf. いわゆる"jet"color、二次元マップのように違いを明確にしたい用途に適している
jet <-  colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))  

# 散布図などで重なりをより見やすくするには、xをjitter() 関数に入れて微小な"ブレ"を加えてやるとよい

# cf.  RGBaを数値で定義するにはrgb()関数。例えば青の透過色の場合、rgb(0, 0, 1, 0.1)。
# rgb(R, G, B, 透明度) の順になっている。いずれも最大が1になるように設定する。
# 逆に色名や色コードからRGBの数値を得るにはcol2rgb()関数。


###### おまけ：EPSで出力する場合 ######################
# 下記、postscript…を実行してから、通常のグラフ作成コマンドを実行、
postscript("graphEPS.eps", height=4, width=4, # サイズを指定
horizontal=F, onefile=F, paper="special") # ここは暗号…
# ファイル名とheight、width（インチ単位…）をお好みで書き換えてください
# 作成終了したらdev.off()を実行（括弧内は空白）
dev.off() # 作業ディレクトリ内へ保存される
投稿者 時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)