読者です 読者をやめる 読者になる 読者になる

山羊の午後

研究関係の備忘録。

一元配置の分散分析

分散分析は3群以上に分けられるデータに差があるかどうかを比較するときに行います。
データを分ける要因が一つである場合、一要因の分散分析を行います。

たとえば、「三年生の国語のテスト結果」をクラスごとに比較する場合、生徒のデータは「クラス」だけで区切られています。よって一要因の分散分析をもちいます。


要因:クラス
A組B組C組
国語の点数 生徒1 : 67点
生徒2 : 63点
生徒3 : 77点
生徒4 : 70点
生徒5 : 80点
生徒6 : 97点
生徒7 : 85点
生徒8 : 45点
生徒9 : 65点
生徒10 : 87点
生徒11 : 46点
生徒12 : 88点
生徒13 : 48点
生徒14 : 55点
生徒15 : 64点
生徒16 : 67点
生徒17 : 62点
生徒18 : 58点
生徒19 : 60点
生徒20 : 87点
生徒21 : 75点
平均点 77.062.067.6

  • 国語の点数は各生徒から一度だけ測定されているので「データにくり返しはない」といいます。
  • データ数 n = 21

等分散性を仮定した一要因の分散分析

class.A <- c(67, 63, 77, 70, 80, 97, 85)
class.B <- c(45, 65, 87, 46, 88, 48, 55)
class.C <- c(64, 67, 62, 58, 60, 87, 75)
#データの入力

dat <- data.frame(class.A,class.B,class.C) 
# 上記表とおなじく行に生徒の点数、列にクラスをとったwide形式の表

dat <- stack(dat)
# Rが扱いやすいstack形式の表に変換する
# ind 列にクラス名(要因の違い)、values 列に点数(従属変数)が並ぶ
# このとき各行は生徒個人を表しています。

oneway.test(data=dat, values~ind, var.equal=T)
# 繰り返しのない一要因の分散分析(等分散性を仮定)を実行する
# データ dat の values を ind の違いによって比較
  • Rでは ~ (チルダ)を使って Y~Xとかくと「YをXによって」という意味を示す
oneway.test()
一要因の分散分析を実行する。var.equal=T で各群の分散は等しいと仮定。Fで等分散性を仮定しないWelchの方法を実行する。
  • 大量のデータを入力する場合、エクセル、テキストファイルで表をつくっておいてread.table()で読み込む。
  • スタック形式の表について

yaginogogo.hatenablog.jp


Oneway.test(var.equal=T)の結果の見方

   One-way analysis of means

 data: values and ind
 F = 2.0441, num df = 2, denom df = 18, p-value = 0.1585

  • F: F value、F値 分散分析の統計量
  • num df: numenator df、分子の自由度、要因の自由度
  • denom df: denominator df、分母の自由度、多くの場合「残差の自由度」
  • p-value:p値、num df, denom dfのF分布表でのF値の有意確率

分散分析では、要因の自由度、残差の自由度から決定されるF分布表をもちいて、算出されたF値の有意確立pを求める。したがって、p値を示すときには F(要因の自由度, 残差の自由度)= F値 もあわせて示す。

p値が有意水準α(通常0.05)より小さいとき、values(点数)はind(クラス)によって異なる、といえる。
ここでは p > 0.05 であるため有意ではない。
したがって「国語の点数はクラスごとに差があるとは言えない」*1


繰り返しのない一要因の分散分析(等分散性を仮定)の結果は、関数 aov()、lm()でも同じものが得られる。

comp.test <- aov(data=dat, values~ind)
summary(comp.test)

comp.test.2 <- lm(data=dat, values~ind)
anova(comp.test.2)

結果は下記のように分散分析表で表示される。

    Df Sum Sq Mean Sq F value Pr(>F)
ind 2 48.1 24.07 0.18 0.837
Residuals 12 1601.6 133.47
  • Df:自由度、Sum Sq 平方和、Mean Sq 平均平方、F値、p値
  • Residuals 残差
  • Mean Sq = Sum Sq/Df
  • F値 = 要因の平均平方/残差の平均平方

等分散性のチェック

等分散でないデータに等分散性を仮定した検定を行うと、誤った結果を得る可能性がある。したがって、Oneway.test(var.equal=T)を使う際には、データが等分散であることを確認しておく必要がある。

データのばらつきを見るにはボックスプロット(箱ひげ図)を作ってみる。

boxplot(data=dat, values~ind)

f:id:yk_uminami:20170225005002p:plain

B組のばらつきがやや大きいように見える。
このデータが統計的に等分散といえるかどうかはバートレット検定を用いて確認できる。

bartlett.test(data=dat, values~ind)

結果はこのように表示される。

   Bartlett test of homogeneity of variances

 data: values by ind
 Bartlett's K-squared = 2.3881, df = 2, p-value = 0.303

バートレット検定では帰無仮説H0は「各群の分散は等しい」であり、p > 有意水準α(0.05)のとき、帰無仮説は棄却できない。したがって、「等しくない、とはいえない」と結論付けられる。*2
このデータでは p = 0.30 > α = 0.05 であるので、等しくないとは言えない。つまり等分散性を仮定した分散分析を用いても悪くない。

いっぽう、p < 0.05のとき帰無仮説は棄却され、「等しくない」という対立仮説H1を採択する。
このとき、等分散性を仮定しないWelchの方法を使うべきである。

一要因の分散分析 Welchの方法

Welchの方法は等分散であってもなくても利用できる。*3

たとえば、下記のようなデータではクラスごとの点数のばらつきは等分散とはいえない。

class.D <- c(67, 70, 77, 70, 80, 70, 85) #平均 77.0
class.E <- c(45, 40, 55, 80, 92, 48, 55) #平均 62.0
class.F <- c(64, 67, 62, 58, 60, 87, 75) #平均 67.6
dat2 <- stack(data.frame(class.D, class.E, class.F))

bartlett.test(data=dat2, values~ind)

   Bartlett test of homogeneity of variances

 data: values by ind
 Bartlett's K-squared = 6.2796, df = 2, p-value = 0.04329

ボックスプロット(箱ひげ図)をみるとE組の点数のばらつきがD、F組より大きいことが分かる。

boxplot(data=dat2, values~ind)

f:id:yk_uminami:20170225003700p:plain

したがって、等分散を仮定した一要因の分散分析ではなく、等分散を仮定しないWelchの方法で一要因の分散分析をおこなうべきである。

oneway.test(data=dat2, values~ind, var.equal=F)

   One-way analysis of means (not assuming equal variances)

 data: values and ind
 F = 2.3196, num df = 2.000, denom df = 10.642, p-value = 0.1458

等分散性を仮定した方法とは、denom dfおよびF値の算出方法が異なる。*4
のでlm()およびaov()では実行できないし、分散分析表も作らない。しかしF値はnum dfとdenom dfに従うのでF(num df, denom df)=F値 と書く。

ノンパラメトリック法による一元配置分散分析:クラスカル・ウォリス法

データの分布が正規分布に当てはまらないときは、ノンパラメトリックな手法を用いることが適切である。
データの分布を確認するには、ヒストグラムと確率密度曲線を描いて判断することが簡便。

hist(dat$values)
par(new=T) #グラフへの上書き許可
plot(density(dat$values),xlim=c(40,100),axes=F,ann=F)
# xlim には表示されたヒストグラムのx軸下限、上限を指定する。
# axes=F 軸の非表示、ann=F ラベルの非表示

hist(dat2$values)
par(new=T)
plot(density(dat2$values),xlim=c(40,100),axes=F,ann=F)

f:id:yk_uminami:20170225003817p:plain

どちらも大まかにいって正規分布にあてはまっているだろう。

日常的なデータの多くは正規分布すると仮定されている。
しかし、たとえば毎日の欠席者数のように0からカウントアップされ、頻度が多くないデータでは、その分布は正規分布に従わないことが多い。

要因:クラス
A組B組C組
十日間の遅刻者数 1日目 : 1人
2日目 : 0人
3日目 : 1人
4日目 : 0人
5日目 : 0人
6日目 : 2人
7日目 : 1人
8日目 : 1人
9日目 : 0人
10日目 : 0人
1日目 : 3人
2日目 : 1人
3日目 : 2人
4日目 : 3人
5日目 : 1人
6日目 : 2人
7日目 : 1人
8日目 : 2人
9日目 : 1人
10日目 : 1人
1日目 : 5人
2日目 : 2人
3日目 : 2人
4日目 : 1人
5日目 : 1人
6日目 : 3人
7日目 : 2人
8日目 : 1人
9日目 : 1人
10日目 : 0人
平均人数 0.61.71.8

  • n = 30 ここでは1~10日目を独立のデータとして扱う。*5
class.A <- c(1, 0, 1, 0, 0, 2, 1, 1, 0, 0)
class.B <- c(3, 1, 2, 3, 1, 2, 1, 2, 1, 1)
class.C <- c(5, 2, 2, 1, 1, 3, 2, 1, 1, 0)

dat3 <- stack(data.frame(class.A, class.B, class.C))

hist(dat3$values)

f:id:yk_uminami:20170225003849p:plain

ヒストグラムから、明らかに正規分布していないことが分かる。
このような場合、正規性を仮定しないノンパラメトリックな一元配置分散分析として、クラスカル・ウォリス検定を用いる。

kruskal.test(data=dat3, values~ind)

   Kruskal-Wallis rank sum test

 data: values by ind
 Kruskal-Wallis chi-squared = 8.4607, df = 2, p-value = 0.01455

p < 0.05 であり、3組の遅刻数には有意な差があると言える。

どれをつかうのかチャート

データは正規分布している 等分散である Oneway.test(var.equal=T)
等分散でない(又はわからない) Oneway.test(var.equal=F)
正規分布していない kruskal.test()

(注)あくまでも目安です。

*1:「差がない」とは言えない。

*2:「等しい」と断言してはいけない。

*3:そのため、最初からWelchの方法を使うべきとする向きもある。

*4:参考:一元配置分散分析 の下部

*5:しかし実際にはこのようなデータは、同じ生徒、同じ曜日、などの繰り返しを考慮すべきだろう。

グラフィックパラメータいろいろ


 Rは図の調整は死ぬほど面倒ですが、その分、自由度は高いです。どんなグラフを描きたいか決めて、一度スクリプトを書いて保存してしまえばあとは楽、なはず。
 調整したい、すべき項目が、描画ウィンドウ(枠、dev()で指定)なのか、図の要素(par()かplot()などのパラメータで指定)なのか、図の付属品(legend()、title()、など)なのかを把握するのが大事です。


par(): グラフィックパラメータの指定

  • 現在のグラフィックウィンドウにパラメータを指定する。また、多くのパラメータはplot()などのグラフ関数内でも指定できる。
  • たいていの項目はpar()のhelpで解説されているので、知りたいことがあったらまずは目を通してみましょう。


描画ウィンドウ (Graphic Device)

dev.new()
新規の描画ウィンドウをつくる。サイズはwidth=5, height=5で指定(単位はinches 1 = 2.54cm, 1.9685 inches = 5 cm)。
dev.size()
現在のウィンドウのサイズを調べる。
dev.list()
現在ある描画ウィンドウ一覧をだす。
dev.set(which= ウィンドウの番号)
描画するウィンドウ(active device)を指定する。

 指定しない限り、グラフは現在のアクティブウィンドウに上書きされます。
 上書きしたくない場合は、dev.new()で新しいウィンドウを作ってからグラフ関数を呼び出します。複数ウィンドウがあるとき、次に描画するグラフは、最も新しいウィンドウ(デバイスナンバーが大きいもの)に上書きされます。

x=c(1,2,3,4,5)
y=c(2,1,3,2,3)
z=c(3,1,3,2,1)

plot(x,y)
dev.new()  #同じ位置に重なって出る

plot(z,y)
a<-dev.list()
a
dev.set(a[1])
abline(lm(y~x))

dev.set(3)
abline(lm(y~z))

描画ウィンドウ内を分割する

par(mfcol=c(列数,行数))
列数×行数に分割する。左上から下方向に順に領域が使用される。
par(mfrow=c(列数,行数))
列数×行数に分割する。左上から横方向に順に領域が使用される。
  • グラフを並べたい順序でどちらか選ぶ。領域は均等に分割される。
  • 枠の数よりも多いグラフが入力されると、最初の枠に戻って描画される。他の枠がリセットされるので、枠の数とグラフの数を合わせておく。
  • 分割した特定の領域だけにグラフを描きたい場合や、任意の大きさに領域を区切りたい場合:layout()やsplit.screen()を使う。
dev.new()
par(mfcol=c(2,2))
for(n in 1:4){plot(1,1,pch=as.character(n))}
#1から4番目の領域に順番にグラフが描かれる。
#5番目のグラフを指定すると1に戻って描かれる。
#(他の枠はリセットされる)

dev.new()
par(mfrow=c(2,3))
for(n in 1:6){plot(1,1,pch=as.character(n))}

dev.set(2)
abline(h=1)
#描き込まれるのは、一つ目のウィンドウの最後の枠

dev.set(3)
plot(1,1,pch="7")
#描き込まれるのは、ウィンドウ3の1番目の枠。
#(他の枠はリセットされる)

図を重ねて書く

par(new=T)
Active deviceにあるグラフに重ね描きを許可する。二つの目のグラフの軸は調整されず、ただ上に重ねられるので、x,y軸の範囲をあらかじめ、xlim=c(),ylim=c()で揃えて指定しておく。
hist(add=T)
描画関数の中には、Add=Tで既に描かれたx,y軸を利用して図を追加して書くことを指定できるものがある。例: hist(), matplot() *1



x軸、y軸

axes()
軸を追加する。side= 1下、2左、3上、4右。at=c(目盛りの場所)。
  • パラメータlabels=T/Fでラベルの有無、labels=c(ラベルの値)で atと異なるラベルを表示。
  • axes()はグラフィックウィンドウ内に設定されている軸を使って線を引くので、軸が設定されていないグラフには設定されていないと
    • *軸の範囲指定
  • 軸を描かない
  • 軸ラベル
  • 軸の0を交差させる(原点の編集): par(xaxs="i", yaxs="i")


マーカー

pch = 0~25
プリセットされた記号を指定。良く使うものは下記。
番号 1 0 2 5 3
×
番号 16 15 17 18 4

文字をマーカーに指定する

pch = "-"
指定した1文字がマーカーになる。複数文字は指定できない。

任意の場所にマーカーを描く

points(x, y, pch=1)
座標(x, y)にマーカーを描く。複数描く場合は、x,yそれぞれにベクトルc()を指定する。




線の種類

任意の線を引く

線の太さ



塗り

色の指定

斜線で塗る


文字のパラメータ

cex = 1
文字の大きさを指定。描画ウィンドウと相対的に決まる。2で二倍サイズ。
ps = 10
文字の大きさをpoint単位で指定。
font= 番号
スタイルを1-4で指定。1標準、2太字、3斜体、4太字斜体。
adj = 数値
文字の揃え位置を指定。0左、0.5中央、1右。c(x, y)で任意のx軸、y軸方向を指定。

タイトル

title(main="title")
メインタイトル(グラフ上)を書く。sub=""でサブタイトル(グラフ下)を書く。xlab,ylabでx,y軸名を書く。

任意の文字を書き込む

text(座標x, 座標y, "text")
グラフの軸内に文字を書く。
mtext("text",side=1~4)
軸の外側に文字を書く。side=1下、2左、3上、4右



凡例 legend(x, y, legend=c(), 描画パラメータ)

  • legendに群名を指定する。各群の描画パラメータと順番が一致することに気をつける。

位置

予約語で指定:相対的に場所を決めてくれるので、便利。
"topleft" "top" "topright"
中央 "left" "center" "right"
"bottomleft" "bottom" "bottomright"
座標で指定: x座標, y座標を指定する。c()でなく、x=X, y=Yと入力する。
  • 凡例の左上角が指定座標になる。
  • 凡例表示のサイズは描画ウィンドウのサイズから相対的に決まるので、座標を固定する場合、注意が必要。
  • グラフの枠外に凡例を描く場合、先に枠外の余白を定義してから、そこに描く指定をする。*2
    • par(oma=c(下, 左, 上, 右)):余白のサイズを指定
    • par(xpd=T):枠外への描画を許可。座標指定でx, y軸の延長線上で枠外を指定できるようになる。(例えばマイナス値)
凡例を2列以上で表示する
  • ncol= 2 :凡例数を指定列数で分割して順に表示する。
  • horiz = T:凡例を列方向に並べる。ncolに優先する。

凡例記号の塗り

  • fill =c(色指定):塗りつぶしが表示される。斜線での塗りつぶし(angle+density)は上手く反映されない。

凡例記号の線

*1:ただし plot() は使えない

*2:一般化しづらいのでお薦めしない

Rの基礎、ベクトルとデータ型

ベクトル

  • Rで使われるベクトルは「複数の値をひとまとめに並べたもの」を意味する。
  • 値には「並べた順番」がつけられ [ ] で示される。番地みたいなもの。これを使って値を呼び出せる。
  • 数値、文字、論理値など色々なデータ型の値をいれることができる。ただし、種類の異なる値を混ぜて入れることは出来ない。
    • たとえば、文字と数字を混ぜて入力すると、数字も文字として扱われ、文字型のベクトルになる。文字型として代入された数字は、計算することが出来ない。
関数 c()
入力された値をひとまとめに並べて、ベクトルを作る。値はベクトル内で順番(例 [1])を与えられる。

f:id:yk_uminami:20160607215302p:plain

変数へのベクトルの代入と、その結果の違いを見る

#変数xに値5を代入する。
#値がひとつの時はc()を使わなくても代入できる。
x <- 5
x <- c(5)
x #変数名をコンソールに入力すると、中身が表示される。

#変数xに値5と10のベクトルを作って代入する
x <- c(5,10)
x

x[1] #xの一番目の値を呼び出す。
x[2] #xの二番目の値を呼び出す。


#変数xに文字「Con」を代入する
x <- "Con"
x <- c("Con")

#変数xに文字ConとTestのベクトルを作って代入する
x <- c("Con","Test")



規則的なベクトルを作る方法いろいろ

x : y :コロン。xからyまでの整数ベクトルを作る。(公差1の等差数列をつくる)
LETTERS[x:y] 変数LETTERSは大文字アルファベットのベクトル。[x:y]でx番目からy番目の文字を抜き出す。
rep(x, y) xをy回くりかえしたベクトルを作る。xにはベクトルも指定できる。
seq(x, y, z) xからyまでをz間隔で区切ったベクトルを作る。等間隔の値が欲しい時に。
  • seq(x, y, 1) = x:y
1 : 10
# 1  2  3  4  5  6  7  8  9 10

LETTERS[1:10]
# "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"

rep(0,10)
# 0 0 0 0 0 0 0 0 0 0

seq(0,1000,250)
#  0  250  500  750 1000



データ型

  • Rはベクトルの中身にいれられた値に「データ型」という分類をつけ、それにしたがって処理を変化させる。そのためデータ型が間違っていると、関数plot()や統計解析から得られる結果が変わってしまう。入力したデータがどんな型として扱われるかは、つねに注意が必要。
    • 調べる関数例 str(), class(), mode() 下記参照

数値 numeric

  • 一般的な「数」の分類。計算できる値として扱われる。他に数を扱う方として、整数型 integralと複素数型 complex がある。

文字 character

  • ""で囲まれて示される。間に半角スペースを入れられる。
  • "1"は文字であり、計算できない。

日付型 Date

  • 年月日を示す値。"2016-04-20"のように""で囲んで表示されるが、文字ではなく、日数の足し引きなどの計算ができる。

論理値 logical

  • TRUE or FALSE の二つの値だけをもつ。T、Fと略して入力することもできる。判定式の答えなど、Yes/Noの答えに使われる。
x<-5 # xに5を代入
x>2  # xは2以上か?と聞く。
# TRUE の答えが返ってくる。

y <- x>2 # 変数yに「xは2以上か?」を代入する。
y    # 式ではなく値だけが入っていることに注意。 

factor型

  • Rならではの型。含まれる要因の番号とその内容(水準、Levels)で構成される。
関数factor()
入力されたベクトルをfactor型にする。値の種類に番号を付け、その中身と対応付ける(例:1番はオレンジ、2番はりんごを意味する、とする)。そのため、データの中身は数になり、数に応じた要因(Levels)が別に保持される。
  • 大きいデータを効率よく保存するために使われる(?)。
  • 水準に大小・優劣などの相対的な順番をつける 関数ordered() もある。
  • 多くの解析関数、作図関数はcharactor型をfactor型にして読み込んで処理するが、データの入力時に混同すると思わぬ結果がでることもある。
例 character型 と factor型の違い

f:id:yk_uminami:20160607221154p:plain

x <- c("orange", "apple", "orange")
x # 値だけが表示される。

y <- factor(x) #xをfactor型にしてyに代入
y # 値とLevels(水準名)

str(y) # 変数yの構造をみる
# Factor w/ 2 levels "apple","orange": 2 1 2
# Levels(水準)は  "apple", "orange" の二つ。 
# データの中身は 2 1 2 つまり「水準2、水準1、水準2」の順に並んでいる。

class(x) # データ型を調べる
class(y) # xは文字型と表示され、yはfactor型と表示される。

mode(x) # データの値の型を調べる
mode(y) # xは文字型と表示されるが、yは数字型と表示される。
注意
水準は指定がないとき、アルファベット順に並べられる。水準に順番があるもの(たとえばSmallとLarge)はfactor(x, levels=c(水準の順番))で指定するか、ordered()をつかうといい。
> x<-factor(c("S","M","S"))
> x
[1] S M S
Levels: M S

> x<-factor(c("S","M","S"),levels=c("S","M"))
> x
[1] S M S
Levels: S M

> x<-ordered(c("S","M","S"),levels=c("S","M"))
> x
[1] S M S
Levels: S < M

データ型の変換

  • as.numeric(x) 数字型にする。文字には適応できない。
  • as.character(x) 文字型にする。
  • as.Date(x) 日付型にする。入力値は"2016-04-20"もしくは"2016/04/20"。他の形式で入力するにはパラメータformatで指定する。
  • as.factor(x) factor型にする。
    • factor型をas.character()に入れると、番号をLevelsに変換した値が、as.numeric()に入れると、水準の番号だけが返される。

グラフの選び方:棒グラフか? 折れ線グラフか?

 データを適切に解釈するには、見やすい、分かりやすいグラフを描くことが必要です。自分の好みでグラフを作っていると、誤った結論に陥ることもあります。
 二つのグラフがどんなものに向いているのかを知っておけば、適切にグラフを使い分けることができるでしょう。

棒グラフ

  • 横軸には連続データ(時間、投与量など)または、非連続データ(カテゴリ、雌雄、条件など)をとる。
  • 棒の面積で各値を量的に判断することが出来る。
    • 積み上げ棒グラフ:各項目の割合と全体量の変化を見ることが出来る。
  • 時系列などの連続データにおいて、長期的な傾向を掴むのに適している。短期的な変化を比較するのには向いていない。
  • 縦軸は省略せず、全体量を表示する。
よく使われる例
年・月ごとの収入・支出の変化、実験条件による効果の違い、毎年の収穫量の変化、雨量や積雪量、候補者ごとの得票数の違い、など。どれだけ多いか、少ないかを見たいときに使う。

折れ線グラフ

  • 横軸には連続データ(時間、日時、年、投与量など)を取る。
    • 非連続データ(カテゴリ、雌雄、条件など)を横軸に使わない。*1
  • 線の傾きにより、各値の連続的な変化を追うことが出来る。とくに短期的な変化を追うのに向いている。
  • 複数項目の折れ線グラフは、各項目の順位や相対的な変化を見ることが出来る。
  • 縦軸は省略しないべきだが、必要であれば変化が見やすい範囲に拡大してもよい。
よく使われる例
毎日の体重の変化、月ごとの気温の変化、時間ごとの反応数の移り変わり、毎分ごとの視聴率、など。直前に比べて上がったか、下がったか、が見たいときに使う。

同じデータを用いた二つのグラフ例

f:id:yk_uminami:20160531010536p:plain

  • 棒グラフ
    • Con群とTest群の反応の量的な違いが見やすい。
      • Test群はCon群の2倍以上反応することが分かる。
    • 投与量における各群の変化は上昇していることは分かるが、投与量ごと変化の違いは分かりにくい。
  • 折れ線グラフ
    • Con群とTest群の投与量に対する反応の変化の違いが見やすい。
      • Test群では10から20の変化のほうが20から30の変化より大きい。*2
    • 群間の量的な違いは分かりにくい。
  • この場合、投与量ごとの効果の変化よりも、群間の効果の違いが見たいので、棒グラフを選ぶ。
group<-c(10,20,30) #投与量
m.Con<-c(0.079,0.102,0.146) #Con群 平均値
sd.Con<-c(0.022,0.051,0.026) #Con群 SD
m.Test<-c(0.12,0.25,0.34) #Test群 平均値
sd.Test<-c(0.055,0.031,0.026) #Test群 SD

yRoof=round(max(m.Con+sd.Con, m.Test+sd.Test)*1.2, 1) 
yFloor=round(min(m.Con-sd.Con, m.Test-sd.Test)*1.2, 2) 

dev.new(width=800, height=400)
par(mfcol=c(1,2))

#棒グラフ
graph1<-barplot(rbind(m.Con, m.Test),beside=T, 
	xlab="Dose of drug (mg/kg)",ylab="Response",
      names=group, ylim=c(0,yRoof), col=c(0,8))

arrows(graph1,rbind(m.Con,m.Test),
       graph1, rbind(m.Con+sd.Con,m.Test+sd.Test),
       angle=90, length=0.1,lwd=2)
arrows(graph1,rbind(m.Con,m.Test),
       graph1,rbind(m.Con-sd.Con, m.Test-sd.Test),
       angle=90,length=0.1,lwd=2)

axis(side=1, c(mean(graph1[,1]),mean(graph1[,2]),mean(graph1[,3])), labels=F) 
legend("topleft", c("Con","Test"), fill=c(0,8)) 
title("Bar graph with two factors")

#折れ線グラフ
matplot(group,cbind(m.Con,m.Test),type="b",ylim=c(yFloor,yRoof),
	lty=1,col=c(1,8), cex=2, pch=17, xlab="Dose of drug (mg/kg)",ylab="Response")

arrows(group,cbind(m.Con,m.Test),group,cbind(m.Con+sd.Con,m.Test+sd.Test),
       angle=90,length=0.1,lwd=2,col=c(1,1,1,8,8,8))
arrows(group,cbind(m.Con,m.Test),group,cbind(m.Con-sd.Con,m.Test-sd.Test),
       angle=90,length=0.1,lwd=2,col=c(1,1,1,8,8,8))

legend("topleft",legend=c("Con","Test"),pch=17,col=c(1,8),cex=1.3)
title("Line graph with two factors")



どっちのグラフにするか判別チャート

  • 横軸にとるデータは
    1. 非連続データである:棒グラフ
    2. 連続データである(大小があり、数直線状に並べられるもの)
    • データから読み取りたいのは
      1. 量的、長期的な変化:棒グラフ
      2. 短期的な変化:折れ線グラフ


 あくまで目安です。データからどんな情報を読み取りたいのか考えて選ぶことが大事です。


yaginogogo.hatenablog.jp
yaginogogo.hatenablog.jp

*1:例外はもちろんある。たとえば群間の交互作用を見たいときなどは、折れ線グラフのほうが変化の違いが見やすい。

*2:この違いに意味があるかは別として。

Rで描く、折れ線グラフ

基本の折れ線グラフ: plot(x,y, type="b")

f:id:yk_uminami:20160529202042p:plain

注意
plot()でマーカー付き折れ線グラフtype="b"を指定するとき、x軸データには数値ベクトルしか指定できない。
#データ入力
Elapsed.days<-c(1,2,3,4,5,8) #観察開始からの経過日数
weight<-c(10.2, 10.5, 10.2, 10.8, 11.9, 12.8) #体重

#マーカー付き折れ線グラフ
plot(Elapsed.days, weight, type="b",
     col="royalblue3", pch=15, cex=2, lwd=2, lty=1)
title("Simple line graph with markers")

plot()のプロパティ

type "b" で線+マーカー、"l"で線のみ、"p"でマーカーのみ
col 線とマーカーの色
pch マーカーの形 0=□、15=■
cex マーカーの大きさ
lwd 線の太さ
lty 線のスタイル 1=実線、2=点線

x軸に日付や文字列を使う

  • 数字に置き換えて入力し、あとからラベル表記を変更する。

f:id:yk_uminami:20160529210749p:plain

#データ入力
Elapsed.days<-c(1,2,3,4,5,8) #観察開始からの経過日数
weight<-c(10.2, 10.5, 10.2, 10.8, 11.9, 12.8) #体重

#折れ線グラフを描く
plot(Elapsed.days, weight, type="b",
     col="royalblue3", pch=15, cex=2, lwd=2, lty=1,
     xaxt= "n", xlab="") #x軸とx軸名を描かない。

#x軸の目盛りのみを描く
axis(side=1, labels=F)

#必要なラベルの数を確認
label.n<-length(axis(side=1,labels=F))
print(label.n)

#日付ラベルをつくる
start.day<-as.Date("2016/4/5")
day.names<-seq(start.day, start.day+label.n-1, 1)
#もしくは
# day.names<-as.Date(c("2016/4/5", "2016/4/6","2016/4/7","2016/4/8",
#          "2016/4/9","2016/4/10","2016/4/11","2016/4/12"))

#x軸を描く
axis(side=1, at=1:8, labels=format(day.names, "%m/%d")) 
title("Line graph with date labels", xlab="Date")
  • 注意 記録がない二日間(day 6,7 = 4/10, 4/11)を忘れないように。

使った関数

as.Date("X")
入力文字をデータ型Dateに変換する。入力形式はformat="%y/%m/%d"で指定できる。
seq(from,to,by)
fromからtoまで間隔byで連続ベクトルを作る。
format(a, "X")
データaを形式Xに変換する。形式例はhelp(strptime)に詳しい。



2群~n群の折れ線グラフ matplot(x, y, type="b")

  • matplot()はyに指定されたデータ表の列を群として描画する
    • plot()とおなじくxには数値データが必要

f:id:yk_uminami:20160529214110p:plain

#データ入力
Elapsed.days<-c(1,2,3,4,5,8) #観察開始からの経過日数
w.male<-c(10.2, 10.5, 10.2, 10.8, 11.9, 12.8) #体重
w.female<-c(9.2, 8.5, 8.2, 9.8, 10.2, 11.3) #体重

#2つのデータを列方向に連結
weight<-cbind(w.male,w.female)

#二群の折れ線グラフ
matplot(Elapsed.days, weight, type="b",
     col=c("royalblue3","brown3"), pch=c(15,0), cex=2, lwd=2, lty=1)

#凡例
legend("topleft", legend=c("Male", "Female"),
     col=c("royalblue3","brown3"), pch=c(15,0), lwd=2, lty=1)

#タイトル
title("Line graph with two groups")

データ表

cbind(w.male,w.female)
ベクトルを列方向に合体する
x y (weight)
Elapsed.days w.male w.female
1 10.2 9.2
2 10.5 8.5
3 10.2 8.2
4 10.8 9.8
5 11.9 10.2
8 12.8 11.3

パラメータ

  • 各パラメータの値1が群1に、値2が群2に適用される。値が群より少ない場合、先頭から値が繰り返される。

凡例

legend(場所x,y, legend=群名, 描画パラメータ)
描画パラメータにはmatplot()で指定した各群のパラメータをコピーする。群名の順番と、パラメータの順番が一致すること。
  • 場所の指定
    • 予約語で指定:相対的に場所を決めてくれるので、便利。
"topleft" "top" "topright"
中央 "left" "right" "center"
"bottomleft" "bottom" "bottomright"
    • 座標で指定: x座標, y座標 凡例の左上角が指定座標になる。凡例表示のサイズは描画ウィンドウのサイズから相対的に決まるので、座標を固定する場合、注意が必要。
  • グラフの枠外に凡例を描く:先に枠外の余白を定義してから、そこに描く指定をする。*1
    • par(oma=c(下, 左, 上, 右)):余白のサイズを指定
    • par(xpd=T):枠外への描画を許可。座標指定でx, y軸の延長線上で枠外を指定できるようになる。(例えばマイナス値)
  • 凡例を2列以上で表示する
    • ncol= 2 :凡例数を指定列数で分割して順に表示する。
    • horiz = T:凡例を列方向に並べる。ncolに優先する。

*1:一般化しづらいのでお薦めしない

Rのはじめ方 3:スクリプト入門

スクリプト

 ここでは「なんらかのコマンドを書いたテキストファイル」をスクリプトとよびます。
 Rの便利なところは、一度かいたコマンドを.txtファイルもしくは.Rファイルとして保存することで、すぐに同じ解析、グラフを再現できるところです。
 例えば、前に作ったグラフを新しいデータで作り直すとき、すでに作ってあるスクリプトの「データの読み込み部分」だけを変更すれば、簡単に同じグラフを作ることが出来ます。

スクリプトファイルをつくる

  • 「ファイル」から「新しいスクリプト」をえらぶ
    • もしくは Ctrl + n
  • Rエディタ ウィンドウが開く
  • 保存は.Rファイル形式で行います(テキストファイルとしてメモ帳などの他のソフトでも開くことができる)
    • .txtファイルで保存してもとりあえず問題はない。好みで。

スクリプトをひらく

  • メニューバーの「ファイル」から「スクリプトを開く」を選択
  • デフォルトではRファイル(.R)のみが表示がされているので、ファイル名の右横の欄で「All files」を選択します。
  • 新しいウィンドウ「Rエディタ」が開き、ファイルの中身が表示されます。


 前回 history()から保存したファイルを開きましょう。同じ様に、アヤメのデータ iris を読み込み、グラフを作ってみます。

data(iris)
head(iris)
plot(iris$Species, iris$Sepal.Length)

 

  • エディタウィンドウの中で、この三行をマウスで選択します。(もしくはシフトキー+矢印で選択)
  • 右クリックから「カーソル行もしくは選択中のRコードを実行」を選択します。
  • コンソールに命令文が入力され、グラフが出ます。
    • 右クリック以外にもメニューバーの「編集」からもコードの実行ができます。
    • また、ショートカットキーとして Ctrl+R が利用できます。
    • すべてを選択(Ctrl+A)して、Ctrl+R でエディタ内全体を繰り返すことができます。

スクリプトの書き方

関数名(パラメータ=数値, パラメータ="文字")

  • 値の間はコンマ , で区切る。()内で改行してもよい。
    • 関数ごとに設定できるパラメータの種類と基本の並び順が決められています。詳細はhelp(関数名)で見ることが出来ます。
    • パラメータ名を省略して値を入力すると、基本の並び順にそって関数が実行されます。
    • 値の順番を変更して入力する際は、パラメータ名をそれぞれ書く必要があります。
    • パラメータによって入力できる値の形式も決まっています。数字のみか、文字のみか、などもマニュアルを見て確認します。
    • パラメータに複数の値を指定するときはベクトル c(値1,値2,値3) と書きます。

コメント #

  • 関数の使い方のメモなどを書き込んでおくと、あとでみてわかりやすくなります。
  • #から改行までの記述は実行されない。()の中で改行と組み合わせて使える
plot(x, y, #コメント  ここは実行されない
  パラメータ=1, パラメータ =2 ) # コメント
# 関数は 終わりを示す ) がくるまでスタートしない。

文字の入力 "Hello world"

  • 文字を値として入力するときは""で囲む。

ベクトル c(1,2,3)

  • 複数の値の集合。関数c()にコンマ,で区切って入力する。

値の判定式 A > B

  • AはBより大きいか? と聞く。答えは真偽値TRUE/FALSEで返ってくる。
    • A
    • A==B AとBは同じか? A=Bは代入なので、間違えないように。

代入 A <- B

  • 変数Aに値Bを代入する。以下の書き方でも同じ。
    • B -> A
    • A = B
変数
""で囲まれていない文字列。なにか値が代入されたもの(箱、ラベルのような概念)。
  • 繰り返し使う値に名前を付けて保存し、名前を使って値を呼び出します。
  • このとき値に付けた名前を変数とよび、変数に値を代入する、と表現します。
  • 変数aに値6を代入するには、代入記号「 <- 」を使って a <- 6 と書きます。
  • 変数名には半角スペースが使えないため、ピリオド「.」で区切ったりします。
  • 変数名は大文字と小文字で区別されます。
  • もともとRが使っている関数名を変数として上書きすることが出来てしまうので、注意。
    • 例 mean <- c(x,x,x,)とすると、関数mean()が入力した値で上書きされる。
    • 直すには、remove(mean) で変数を削除する。変数名が使えるかは、変数名をコンソールに入力してみるといい。

代入の例

 身長と体重を入力してみます。

height <- 172 # cm
weight <- 72 # ㎏

 コンソールに「height」と入力してエンターをおすと下記のように代入された値が表示されます。

> height
[1] 172

  • 左端の[1]は一番目の値であることを示します。

 二つの値をつかって、BMIを計算してみます。

# BMI=体重÷身長(m)÷身長(m)
BMI<-weight/(height*0.01)/(height*0.01)
BMI

 代入した値を変数名で呼び出して計算しています。
 では次に weight に三日分の体重を代入し、それぞれのBMIを計算してみます。

height <- 172 # cm
weight <- c(72.0, 71.6, 71.0) # ㎏
BMI<-weight/(height*0.01)/(height*0.01)
BMI

複数の値を代入するときには、c() をつかい、値はコンマ , で区切ります。
結果はこのように返されます。半角スペースが値の区切りとして表示されます。

[1] 24.33748 24.20227 23.99946

体重weightに代入されたすべての値に対して、BMIの計算式が繰り返し適用されました。

ではBMIの変化をプロットするグラフを書いておきましょう。

weight <- c(72.0, 71.6, 71.0) # ㎏

height <- 172 # cm
BMI<-weight/(height*0.01)/(height*0.01)

plot(BMI, type="b")

このスクリプトを保存しておけば、weightに代入する値を追加するだけで、毎日、BMIの変化をグラフ化できます。簡単ですね。

自分に必要なスクリプトを書けるようになることで、Rは一気に便利になります。
次は、スクリプトを書きこなすために、Rの基本ルールをおさえましょう。

Rのはじめ方 2:使い方 入力と出力

 Rの動作をおぼえるために、とにかく使ってみましょう。
 基本的な動作はコマンドを入力して、結果を出し、保存する、この三つだけです。

入力する

 
 Rを起動して出てくるウィンドウには「R Console」と名前がついています。このコンソールが、Rに命令文(コマンド)を入力する場所です。


 Consoleのウィンドウをクリックして、下の赤い>横にカーソルを合わせます。
 ここで「コマンドを書き、エンターを押す」ことでRに命令文をつたえます。


 まず、下の二行をコピぺして、コンソールに入力しましょう。

data(iris)
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
  • コンソールでは 命令文は赤い文字、 結果は青い文字で表示されます。


 ここでは関数data() でRに内蔵されているデータセット「iris」を読み込み、その先頭部分を 関数head() を使って表示しました。
 irisのデータセットはアヤメの三品種を比較した実際のデータです*1


 では次にグラフを作ってみましょう。
 グラフの描画には関数plot()をつかいます。下記をコピペしてコンソールに入力します。

plot(iris$Species, iris$Sepal.Length)

 円たーを押すと下図のようにグラフができます。

 これはアヤメの三品種(Species)における、がく片の長さ(Sepal.Length)の平均とばらつきを示す箱ひげ図です。

 関数 plot(データ1, データ2) は入力されたデータ1をx軸、データ2をy軸としてグラフを描きます。

 ここでは、x軸に iris$Species 、y軸に iris$Sepal.Length を入力しています。

 これはそれぞれ、表 iris の 列名 Species (品種名)、と 列名 Sepal.Length (がく片の長さ)を指定する書き方です。 表iris の 列Species を iris $ Species と書くわけです。$をつかったデータの指定は頻出しますので、慣れてください。

$マーク
データ表の中で名前の付けられた一列を指定する。表名$列名。


f:id:yk_uminami:20160423210255p:plain



 関数 plot(データ1,データ2) は入力されたデータ1をx軸、データ2をy軸としてグラフを描きます。
 そのため、二つのデータの順番を入れ替えて入力すると、異なるグラフが描かれます。

plot(iris$Sepal.Length, iris$Species)

f:id:yk_uminami:20160424011625p:plain

 x軸にがく片の長さ(iris$Sepal.Length)、y軸にアヤメの品種(iris$Species)が入力されたグラフができました。

 でも前回は箱ひげ図になったのに、今回は点でプロットされています。なぜでしょう?


 これは iris$Species が品種名の入力された「文字型」のデータ*2
であるせいです。


 関数plot()には「x軸に文字型が入力された場合は、それぞれの群ごとに箱ひげ図を描く」と定義されています。しかし、y軸方向では文字の扱いは定義されていません。

 そのため、上のグラフを見ると、y軸には品種の違いに応じて一番目の品種、二番目の品種、三番目の品種、と適当な「数字」に変換されて軸が描かれています。

 これではいったい何のグラフなのかわかりません。。。


 正しいグラフを描くためには、関数plot(x軸データ、y軸データ)がもつ入力ルールにおぼえておかねばいけません。Rに習熟する労力の大半は、この コマンドの書き方を覚える ことです。

 それぞれの関数がもつ入力ルールは help(関数名) を使えばオンライン上のマニュアルが参照できます。(ただし英語)
 また、関数名+R で検索すれば、大抵のものは日本語での解説をみつけることができます。



グラフを出力する

 グラフを保存してみましょう。方法は三つあります。

グラフの上で右クリック

メタファイルに保存
.emf形式で画像を保存します。ベクター画像であるため、拡大縮小してもきれいです。パワポに貼ったりするならこちらで保存します。ウィンドウズ以外での互換性に注意が必要です。
ポストスクリプトファイルに保存
.eps形式で画像を保存します。ラスター画像であるため、ややぼわっとしますが、編集が容易だったり、互換性が確かだったりします。

メニューバーの「ファイル」から「別名で保存」を選択

  • Metafile (.emf)、Postscript(.eps)、PDF、PngBmpTIFFJpeg の各形式で保存できます。
  • 「別名で保存」が見つからないときは、グラフのウィンドウ「R Graphics: Device 2」をクリックしてからメニューバーを開きます。(選択されているウィンドウごとにメニューバーの内容が変わります。)

コンソールから命令する

  • 下記を実行すると、PDF形式で現在のフォルダに保存されます。(初期設定ではマイドキュメント。)*3
dev.copy2pdf(file = "test.pdf")



統計解析結果を出力する

 続いて、このデータで統計解析もしてみましょう。
 がく片の長さが三品種の間で有意に違うかを、分散分析で検討します。

oneway.test(iris$Sepal.Length~iris$Species)


 実行すると下記のようにF検定の結果が出力されます。

        One-way analysis of means (not assuming equal variances)

data:  iris$Sepal.Length and iris$Species
F = 138.9083, num df = 2.000, denom df = 92.211, p-value < 2.2e-16

 P値が十分に小さい(0.001以下*4)ので、これらの三品種の間で、がく片の長さに有意に違いがある、といえます。

 コンソールウィンドウの結果をマウスで選択してメモ帳やワードにコピペして保存しておきます。


 では最後に、ここまで入力した命令文を保存しておきましょう。

history()

 と入力すると、「R History」のウィンドウが開きます。
 ここまで入力したコマンドが一覧されていると思います。

 これをメニューバーの「ファイル」から「ファイルを保存」でテキストファイル(.txt)として保存します*5


 次に使うときは「ファイル」から「スクリプトを開く」でテキストファイルを開きます。これで、一度使ったコマンドを再利用することができます。


 ここまでが、Rでやることの一連の流れです。
 それぞれの関数の使い方をおぼえるのは厄介ですが、しかし、すべてを丸暗記する必要はありません。適宜、テキストファイルにメモを取っておけば便利です。

 つぎはコマンドをメモしたテキストファイル、スクリプトの使い方を説明します。


初歩のTips

  • コンソール入力場所で ↑キーを押す 一つ前に入力したコマンドを表示
  • コンソール入力場所で 変数、関数名を途中まで入力してtabキーをおす 残りの入力が保管される。(複数ある場合は一覧が出る)

*1:irisの正体 (R Advent Calendar 2012 6日目) - どんな鳥も

*2:詳しくはyaginogogo.hatenablog.jp

*3:現在のフォルダを変更するには「ファイル」から現在のディレクトリを変更を選択。

*4:2.2e-16は指数表記で、2.2 × 10のマイナス16乗をしめす。とても小さい、くらいの認識

*5:打ち間違いの行はあとで消しましょう