とりあえず作図

（データへのアクセス、ヒストグラム、箱ひげ図、散布図）

ここでは、人に説明する用や論文用のためのキレイな作図ではなく、データの傾向をつかむための自分用の作図について簡単に解説する。

きれいな図を作成するための方法については、こちらのサイト（R でプログラミング：データの一括処理とグラフ描き）で大変詳しく解説されている。

さらに、最近はR使いの間ではパッケージ"ggplot2"が標準的になりつつあるようだ。

データへのアクセス

手持ちのデータの傾向をつかむためには、データを様々な角度から図示して眺めるのが有効である。

そのためには、まずデータを様々に切り取る方法を知る必要がある。

データフレームから一部を切り取る方法には、subset関数、カギかっこ[]、ドルマーク$、head関数が挙げられる（データフレームはRの標準的なデータ形式）。

Rには最初から用意されているデータセットがいくつか存在するので、それを使ってひとつずつ説明する。

データ"iris"

"iris"は、アヤメのがく片と花びらの長さ・幅を種ごとに測定したデータセットである。

この項ではこのデータを使ってやり方を説明していく。

subset関数

subset関数は二番目の引数に"条件式"を指定することで、その条件に合うデータを取ってくる関数である。

subset(iris, Species=="setosa") #列Speciesがsetosaであるデータを返す

subset(iris, Sepal.Length==5) #列Sepal.Lengthが5であるデータを返す

subset(iris, Sepal.Length>5) #列Sepal.Lengthが5より大きいデータを返す

subset(iris, Sepal.Length<=5) #列Sepal.Lengthが5以下であるデータを返す

subset(iris, Species!="setosa" & Sepal.Length>5) #列Speciesがsetosaでない、かつ、列Sepal.Lengthが5より大きいデータを返す

subset(iris, Species=="setosa" | Species=="virginica") #列Speciesがsetosaである、または、列Speciesがvirginicaであるデータを返す

カギかっこ[]

データフレームのデータには行番号と列番号がついている。

データフレーム名の後にカギかっこ[]をつなげ、かっこ内に行番号と列番号を指定することで、特定の行または列のデータを取ってこれる。

iris[1:3, ] #1～3行目を返す。1:3は1から3までの公差1の等差数列（1, 2, 3）を示す

iris[, 1:3] #1～3列目を返す

iris[, c(1, 3, 5)] #1, 3, 5列目を返す

iris[, c("Sepal.Length", "Petal.Length")] #列Sepal.LengthとPetal.Lengthを返す。このように列名を指定して呼び出すことも可能

iris[2:6, 1:3] #2～3行目かつ1～3列目を返す

iris[, -1] #1列目以外を返す

ドルマーク$

データフレームには列名がついている。

データフレーム名の後にドルマークをつけ、その直後に列名をつけることで、列をベクトルとして取り出すことができる（ベクトルはRの標準的な変数の形式）。

subset関数と組み合わせることで、取り出したい条件の変数をベクトルとして呼び出せる。

iris$Sepal.Width #列Sepal.Widthをベクトルとして返す

subset(iris, Sepal.Length>5)$Sepal.Width #列Sepal.Lengthが5より大きいデータのSepal.Widthをベクトルとして返す

subset(iris, Species=="setosa")$Sepal.Width #列SpeciesがsetosaであるデータのSepal.Widthをベクトルとして返す

head関数

データフレームの1～6行目を呼び出す関数である。

とりあえずデータフレームの見た目をチェックするのに使える。

head(iris) #データの1~6行目を返す

以上のスクリプトの出力結果はこちら。

ヒストグラム

ヒストグラムは横軸にデータの値の取る範囲、縦軸に値の範囲ごとのデータの個数を取る棒グラフで、データの分布の形を観察するためのグラフである。

データ解析でまず重要なのが、データの分布が一山型になっているかどうかである。

ANOVAなど多くの解析では、データが正規分布、つまり左右対称の一山型であることを仮定している。

したがって、ヒストグラムでデータの分布を観察しておくことは、解析の前提が満たされるかどうかを確認する手段のひとつとなる。

また、データの分布が一山型でないことは、データに影響している他の要因の存在を暗示する場合もある。

以下にirisを使った例を示す。

ヒストグラムを描くには、関数histを用いる。

histのカッコ内に、数値ベクトルを一つ入れることでヒストグラムが描かれる。

まずは、先ほどのirisのデータを用いて、花びらの長さのヒストグラムを描いてみる。

head(iris) #データの確認

hist(iris$Petal.Length) #全種のPetal.Lengthのデータのヒストグラム

iris$Petal.Length #データの確認、値の範囲ごとに数を数えれば、ヒストグラムの高さと同じになっているはず

ここでは、データを呼び出すためにドルマーク$を使っている。

> head(iris) #データの確認

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1 5.1 3.5 1.4 0.2 setosa

2 4.9 3.0 1.4 0.2 setosa

3 4.7 3.2 1.3 0.2 setosa

4 4.6 3.1 1.5 0.2 setosa

5 5.0 3.6 1.4 0.2 setosa

6 5.4 3.9 1.7 0.4 setosa

> hist(iris$Petal.Length) #全種のPetal.Lengthのデータのヒストグラム

> iris$Petal.Length #データの確認、値の範囲ごとに数を数えれば、ヒストグラムの高さと同じになっているはず

[1] 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 1.5 1.6 1.4 1.1 1.2 1.5 1.3 1.4 1.7 1.5 1.7 1.5 1.0 1.7 1.9

[26] 1.6 1.6 1.5 1.4 1.6 1.6 1.5 1.5 1.4 1.5 1.2 1.3 1.4 1.3 1.5 1.3 1.3 1.3 1.6 1.9 1.4 1.6 1.4 1.5 1.4

[51] 4.7 4.5 4.9 4.0 4.6 4.5 4.7 3.3 4.6 3.9 3.5 4.2 4.0 4.7 3.6 4.4 4.5 4.1 4.5 3.9 4.8 4.0 4.9 4.7 4.3

[76] 4.4 4.8 5.0 4.5 3.5 3.8 3.7 3.9 5.1 4.5 4.5 4.7 4.4 4.1 4.0 4.4 4.6 4.0 3.3 4.2 4.2 4.2 4.3 3.0 4.1

[101] 6.0 5.1 5.9 5.6 5.8 6.6 4.5 6.3 5.8 6.1 5.1 5.3 5.5 5.0 5.1 5.3 5.5 6.7 6.9 5.0 5.7 4.9 6.7 4.9 5.7

[126] 6.0 4.8 4.9 5.6 5.8 6.1 6.4 5.6 5.1 5.6 6.1 5.6 5.5 4.8 5.4 5.6 5.1 5.1 5.9 5.7 5.2 5.0 5.2 5.4 5.1

図1. 花びらの長さのヒストグラム。横軸は花びらの長さの値の取る範囲である。たとえば1から1.5の間にはデータが38個くらい、1.5から2の間にはデータが13個くらいある、というように読む。

このヒストグラムを見ると、花びらの長さのデータは1～2と2.5～7の二つの山に明らかに分かれていることがわかる。

これは、データが何らかの分類に基づくいくつかの群に分けられる（花びらの長さに対して何らかの要因が影響している）ことを暗示している。

このヒストグラムは、すべての種の花びらの長さをまとめて一つのデータとして扱って描いたものだった。

種ごとに花びらの長さが違うことが、このような一山型でないヒストグラムをもたらした可能性がある。

次は、種ごとに花びらの長さのヒストグラムを描いてみよう。

par(mfrow=c(3, 1)) #一つの図に複数のパネルを並べる。ここでは3行1列の表形式でパネルを並べる

hist(subset(iris, Species=="setosa")$Petal.Length) #subset関数でsetosaのデータを抽出

hist(subset(iris, Species=="versicolor")$Petal.Length)

hist(subset(iris, Species=="virginica")$Petal.Length)

ここで活躍するのがsubset関数とドルマーク$の合わせ技である。

また、ここで登場したpar関数は、作成する図の見た目に関して様々な調整を行うための関数である。

parのカッコ内でmfrow=c(x, y)（xとyは行数と列数）と指定することで、一つの図に複数のパネルをx行y列の表形式で並べることができる。

図2. 種ごとの花びらの長さのヒストグラム。

種ごとにヒストグラムを描いてみると、どれもおおむね一山型になっていることが確認できる。

setosa（一番上のパネル）のヒストグラムの値の取る範囲は1～2、versicolorとvirginicaはそれぞれ3～5.5、4.5～7であり、図1の左の山はすべてsetosaのデータだったことがわかる。いっぽう、図1の右の山はversicologrとvirginicaのデータの双方が含まれていたこともわかる。

どうやら、種ごとに花びらの長さが異なる（要因"種"が花びらの長さに影響を与えている）ようである。

ただし、図2では横軸のスケールが異なっており、種内でのデータの分布は見やすいが、種間での比較が難しい。

次は、横軸のスケールをパネルごとに合わせてみよう。

par(mfrow=c(3, 1))

hist(subset(iris, Species=="setosa")$Petal.Length, xlim=c(1, 7)) #引数xlimでx軸の値の範囲を指定

hist(subset(iris, Species=="versicolor")$Petal.Length, xlim=c(1, 7))

hist(subset(iris, Species=="virginica")$Petal.Length, xlim=c(1, 7))

ここでは、histのカッコ内にxlim=c(1, 7)という指定を追加している。

これはx limitの略で、x軸の取る範囲を指定する追加のオプションのようなものである。

図3. 種ごとの花びらの長さのヒストグラム。横軸のスケールをパネルごとに揃えたもの。

横軸のスケールをパネルごとに合わせた結果、種ごとに花びらの長さが大きく異なることが浮き彫りになってきた。

setosaは他の2種に比べてずいぶん値が小さく、versicolorとvirginicaでは前者の方が値が小さいがその差は大きくないことがわかる。

もう一つこの図で重要なのが、種間で値のばらつきの大きさも異なることである。

setosaは1cmの範囲内でしかばらついていないが、versicolorもvirginicaもともに2.5cmもばらついている。

このことは、花びらの長さの平均が種間で異なるだけではなく、花びらの長さの分散も種間で異なること（異分散）を示唆する。ANOVAなどの一般線形モデルを用いた解析法は群間で分散が同じであること（等分散）を仮定するため、このようなデータにはそれなりの対処を施す必要がある。

最後に、図3を図1と比べてみよう。図3は図1を種という要因で分割したものであることがわかるだろう。

箱ひげ図

箱ひげ図はヒストグラムを上から見たようなグラフである。

箱ひげ図は、関数boxplotにベクトルを入れることで描くことができる。

a <- c(2, 4, 4, 3, 3, 6, 4, 5, 6, 4, 4, 1, 12, 9, 6, 7) #でたらめな数値ベクトルを作成し、aに保存

boxplot(a)

図4. 箱ひげ図の例。

このように、箱とひげが描かれる。

箱の中にある太い線が中央値（順位が上から、または下から数えて50%の数）、箱の両端が四分位点を示す。

また、ひげの両端はデータの端を示す。ひげの外側にある点は、boxplot関数の内部で外れ値と判断されたデータを示している。

図5. 数列aのヒストグラムと箱ひげ図の比較。

箱ひげ図の意味は、このようにヒストグラムと並べてみるとわかりやすい。

図5を見ると、ヒストグラムの山と箱の大きさおよび中央値がほぼ一致している。データの範囲も同様である。

これを見ると、箱ひげ図はヒストグラムを上から見たようなものだという意味が分かると思う。

このことは、ヒストグラムが一山型である限り、両者が表す情報はほとんど同じだということを示している。

ヒストグラムがデータの分布の"形"を表現するのに長けているいっぽう、箱ひげ図は一山型のデータの真ん中と幅をヒストグラムより少ないスペースで表すことに特長がある。

図3のように、多群のデータを比較したい場合、ヒストグラムではスペースを食うし、群ごとの分布の形に目がいってしまう。実はこのような目的には箱ひげ図の方が適している。

次はirisを使って、図3と同じデータを箱ひげ図で表現してみよう。

ここで導入するのが、モデル式と呼ばれる記法である。

科学における図は多くの場合、x軸(横軸)がy軸(縦軸)に影響を与える、ということを暗黙の前提として描かれる。

これを式にすると、y = x と表現できる。

モデル式の記法では、これを y ~ x のように記述する。

boxplot関数はこの記法を利用して、直感的な作図を行うことができる。

#まずモデル式を記述し、そのあとに用いるデータフレーム名を指定する

#モデル式では、" ~ "の左にy軸の変数を、右にx軸の変数を記述する

boxplot(Petal.Length~Species, iris, ylab="Petal length") #y軸には変数Petal.Lengthを、x軸には変数Speciesを指定、さらに引数ylabでy軸のラベルを指定

図6. 箱ひげ図による花びらの長さの種間比較。

図3と比べると、setosaが他の2種よりSepal.Lengthが短く、分散も小さいというデータが、パネル1枚ですっきりと表現されていることがわかる。

特に群間の比較を行いたい場合は、各群のデータが一山型である限り、ヒストグラムより箱ひげ図の方が適している。

通常、文字列の変数は、因子ベクトルも含めて、アルファベット順に並べられる。

図6でも、箱は左からアルファベット順に並んでいる。

しかし、図示の都合上、別の順序に並べたいときもあるだろう（順序変数のときとか）。

その場合は、以下のように入力すれば解決する。

iris$Species <- factor(iris$Species, levels=c("versicolor", "setosa", "virginica")) #factor関数の中に順序を変更したい変数を入れ、引数levelsに新たな順序を指定する

boxplot(Petal.Length~Species, iris, ylab="Petal length")

図7. 図6の順序を変えたもの。

ヒストグラムのセクションでも説明した異分散性に対処する方法として、変数変換がある。

箱ひげ図を使うことで、変数変換が等分散に近づける効果があったかどうかある程度確認することができる。

par(mfrow=c(1, 3))

boxplot(Petal.Length~Species, iris, ylab="Petal length", main="raw") #引数mainでタイトルを書く

boxplot(log(Petal.Length)~Species, iris, ylab="Petal length", main="log") #y軸の変数を関数logでかこむ。関数logは自然対数

boxplot(sqrt(Petal.Length)~Species, iris, ylab="Petal length", main="sqrt") #y軸の変数を関数sqrtでかこむ。関数sqrtは平方根√のこと

このように、単に変換したい関数でモデル式内の変数yを囲んでしまえばよい。

図7. 変数変換前と後との箱ひげ図の比較。

比較すると、log変換では種間で分散がある程度同じ大きさになっているように見える。

平方根変換はlog変換に比べると、まだsetosaの分散が小さい。

もし、ANOVAによって花びらの長さを種間で比較するなら、log変換が有効そうである。

散布図

散布図は、2つの連続変数の関係を調べるためのグラフである。

plot関数によって描くことができる。

x <- c(1, 2, 3, 4, 5, 6) #x軸のデータを適当に作成

y <- c(1, 3, 4, 6, 4, 6) #y軸のデータを適当に作成

plot(x, y) #最初にx軸のデータを、次にy軸のデータを入れる

図9. データxとyの関係を表した散布図。

言うまでもないかもしれないが、x軸のデータとy軸のデータの対応関係がここには表現されている。

xの1つめの要素とyの1つめの要素はともに1なので、x, y座標1,1の場所に点が打たれている。他の要素も同様である。

次は、irisのデータを使って花びらの長さとがく片の長さの関係をみてみよう。

plot関数もboxplot関数と同様にモデル式が使える。

モデル式を用いて花びらの長さとがく片の長さの関係を散布してみる。

plot(Sepal.Length ~ Petal.Length, iris) #がく片の長さをyに、花びらの長さをxに指定

図10. 花びらの長さとがく片の長さの関係。

全体では、花びらが長いとがく片も長いといった、正の相関関係がありそうである。

しかし、よく見るとはっきりした正の関係があるのはPetal.Lengthが3以上の部分だけで、2以下の飛び離れた点の塊ではそのような正の関係はほとんどないようにも見える。

注目したいのは、横軸のデータは図1のヒストグラムで表したデータとまったく同じであることだ。

図1のヒストグラムは図3のヒストグラムが合わさったものだった。

同様に、この散布図から示唆されるのは、花びらの長さとがく片の長さの関係が種ごとに異なっている可能性である。

そこで、ヒストグラムと同様に、種ごとに分けて散布図を描くということをやってみる。

ここでも活躍するのがsubset関数である。

subset関数によって種ごとのデータを抽出し、それをさらにplot関数に入れればよい。

par(mfrow=c(2, 2))

plot(Sepal.Length ~ Petal.Length, subset(iris, Species=="setosa"), main="setosa")

plot(Sepal.Length ~ Petal.Length, subset(iris, Species=="versicolor"), main="versicolor")

plot(Sepal.Length ~ Petal.Length, subset(iris, Species=="virginica"), main="virginica")

par(mfrow=c(2, 2))

for (i in levels(iris$Species)){ #プログラミングに長けている人ならば、for構文を使えば同じ図をもっと短いスクリプトで描ける

plot(Sepal.Length ~ Petal.Length, subset(iris, Species==i), main=i)

}

図11. 花びらの長さとがく片の長さの関係の種間比較。

図10と図1, 3を合わせて考えた通り、やはり種ごとに花びらの長さとがく片の長さの関係は異なるようである。

versicolorとvirginicaはPetal.Lengthが長くなるほどSepal.Lengthも長くなっているが、setosaにはそのような正の相関関係はほとんどみられない。

このように、ヒストグラムでやったような、いろいろデータを切ってみるアプローチが散布図でも有効である。

小技として、unclass関数を用いた色分け、形わけがある。

uncalss関数は様々なオブジェクトを数字にしてしまう関数なのだが、それはそれとして、以下のように入力すれば種ごとに色を分けて散布図を作成することができる。

plot(Sepal.Length ~ Petal.Length, iris, pch=16, col=c("black", "red", "blue")[unclass(Species)])

legend("topleft", pch=16, legend=c("setosa", "versicolor", "virginica"), col=c("black", "red", "blue"))

なお、legend関数は凡例を描き加える関数で、まず場所を書き、その後に必要な情報を入力する。

図12が出力される図である。

ここまでくると、ほぼ人に見せても差し支えない図になってきている。

このような小技は、ggplot2の普及とともに忘れ去られていく技術かもしれないが…。

図12. 花びらの長さとがく片の長さの関係の種間比較。色分けでやったもの。

’18 10/30

Google Sites

Report abuse