生物統計学−Rの使い方

授業で使う範囲でRの使い方を説明します。

データの種類、構造と変数

データの種類（型）：数値、文字列
数値は計算できる数。他のプログラミング言語と違い、整数と実数の使用上の区別はない。文字列は文字を表示するためのデータで、計算はできない。数値や文字列のようなデータの種類をデータの型という。
ベクトル
ベクトルは同じ型のデータを複数同時に取り扱える構造で、Rでは最も基本の構造となる。例えば、一つの数値を変数に代入するとRでは長さ１のベクトルとして扱う。他のプログラミング言語ではこの構造を「配列」と呼ぶが、Rの配列は別のものなので注意。c(1,3,5)、c("月", "火", "水", "木", "金")、c(1:100)のように表記する。
因子ベクトル
統計ではいくつかの水準をもった因子をよく使うが、それを表現するベクトルで、factor(ベクトル)で通常のベクトルを因子に変換できる。因子は分散分析には欠かせない。
リスト
リストはベクトルに似ているが、要素のデータ型がそれぞれ違うものである。list(1, "a", 2, "b")のように表記する。リストの中にリストやベクトルも入れることができる。
行列
行列は同じ長さの行を複数集めた、二次元の表形式のデータで、マトリックスとデータフレームの２種類がある。
マトリックス
行列のデータが全て数値型で、全体を計算に使う時に利用する。
データフレーム
データに数値、文字列、因子を含む場合に使用する。分散分析に適している。
配列
配列は行列を三次元以上に拡張したもの。
変数
プログラミング言語での変数は数学とは違い、データを入れておく容器のようなものと考えると分かりやすい。変数の名前は文字で始まり、文字、数字、ピリオドの組み合わせとする。半角数字から始めることはできない。大文字・小文字の区別があるので、resultとResultは別のもの。日本語（全角文字）も使用可能だが、一般的にプログラム中では文字コードの問題があるため日本語は使わないほうがトラブルが少ない。最近はUTF-8が標準になりつつあるため、以前に比べると気にすることは少なくなった。変数への代入（付値）はx<-c(2,5,9)のように「<-」を使う。

データフレームの扱い方

data.frame(ベクトル１, ベクトル２, ...)
複数のベクトルを合わせてデータフレームを作る。列の名前は元の変数名になる。data.frame(名前1=ベクトル1, 名前2=ベクトル2)のように書くと、列の名前を変更できる。
> x<-c(2,5,9)
> y<-1:3
> z<-data.frame(x,y)
> z
subset(データフレーム, 条件式)
条件式に合う行のみを抽出する。

特殊な数値

NULL
存在しないことを示す。
TRUE
比較演算の結果に使われ、真であることを示す。Tと省略できる。
FALSE
比較演算の結果に使われ、偽であることを示す。Fと省略できる。
NA
欠測値を示す。
Inf
無限大（例えば1÷0の結果など）を示す。–Infで –∞ を表す。
NaN
数学では不定となる値（例えば0÷0の結果など）を示す。Not a Numberの頭文字。

ベクトル、行列の計算

Rではベクトルに対する四則演算はベクトルの要素それぞれに対して四則演算した結果のベクトルとなる。例えば

> c(1,2,3) * 3
[1] 3 6 9

のようになる。同様に行列に対しても各要素の計算を行い、その結果が行列として返る。

基本的な演算子

四則演算
+, -, *, / を使う。
べき乗
^または**を使う。2^3は8になる。2**3と書いても同じ。
等差数列
:は公差1または–1の等差数列を与える。1:5は1, 2, 3, 4, 5になる。
比較演算子
同じであるかどうかは==を使う。その他、>, <, >=, <=があり、いずれも成立すればTRUE, しなければFALSEとなる。
行列のかけ算
%*%を使う。
> rv <- runif(6)
> rv
[1] 0.5458232 0.7941123 0.0848588 0.9254300 0.4667509 0.0112287

> a <- matrix(rv,3,2)
> a
          [,1]       [,2]
[1,] 0.5458232 0.9254300
[2,] 0.7941123 0.4667509
[3,] 0.0848588 0.0112287

> b <- c(3, 5)
> a * b
          [,1]       [,2]
[1,] 1.6374695 4.6271499
[2,] 3.9705613 1.4002528
[3,] 0.2545765 0.0561436

> a %*% b
          [,1]
[1,] 6.2646194
[2,] 4.7160914
[3,] 0.3107201

[2016.6.22 加筆]

関数

Rの関数はあるデータと条件を指定し、何らかの操作をした結果が新しいデータとして返ってくる仕組みである。例えばsum(変数)は変数が数値ベクトルのとき、その合計を計算する。関数は組み込みのもの、ライブラリとして提供されているもの、自分で定義するものがある。

基本的なRの関数について

summary(変数)
変数の型に応じてデータを要約したものを表示する。データが数値ベクトルの場合、最大、最小、平均、四分位数が、因子ベクトルの場合、水準とその頻度が表示される。
attributes(変数)
変数の型に関するデータを表示する。変数の行、列ラベル、型の情報など。
nrow(変数)
変数が行列のとき、その行数を返す。
length(変数)
変数がベクトルのときはデータ数、マトリックスのときは全データ数（行×列）、データフレームのときは列数を返す。
sum(変数)
変数の合計を計算する。
sqrt(変数)
変数の平方根を計算する。変数^0.5としても同様に計算できる。

因子を操作するRの関数について

factor(ベクトル)
数値または文字列のベクトルを因子ベクトルに変換する。
interaction(因子ベクトル, 因子ベクトル, ...)
複数の因子を組み合わせた新しい因子を作る。

作業ディレクトリの変更

Windowsの場合

Fileメニューから「Change dir...」を選び、マイドキュメントなどに作ったR用のフォルダを選択する。

Macintoshの場合

その他メニューから「作業ディレクトリの変更...」を選び、R用のフォルダを選択する。

データの読み込み

dir()と入力すると作業ディレクトリの中のファイル一覧が表示される。

csv形式のファイルは次のようにするとデータフレームとして取り込める。以下の例で#記号の右側はRにおけるコメントで、何を書いても無視される。入力する必要もない。

data0<-read.table("ファイル名", sep=",", header=T) #先頭行がカラム名の場合

data0<-read.table("ファイル名", sep=",", header=F) #先頭行から数値の場合

データの要約

summary関数でデータの分布範囲を把握する。

hist関数で度数分布図を作成し、正規分布なのか特別な分布（二つの正規分布が重なっているなど）なのかを見る。

データが一つの場合

hist(x)

Rだとこれだけでいいので簡単。

二つのデータを一つのグラフで重ねて表示する

# 透過色を使う
hist(a, breaks=seq(0,100,5), col="#FF00007F", xlim=c(0,100), main="", xlab="" )
hist(b, breaks=seq(0,100,5), col="#0000FF7F", add=T)

t検定に関するRの関数について

subset(データフレーム, 条件式)
条件式に合う行のみを抽出する。
mean(変数)
変数の平均を計算する。変数は通常、ベクトルとする。
sd(変数)
変数の不偏標準偏差を計算する。
var(変数)
変数の不偏分散を計算する。
標本分散・標本標準偏差の計算 [2014-06-27修正]
Rには組み込みで標本分散や標本標準偏差を計算する関数はないので次のようにする。
dataが数値ベクトルの場合、
m1<-mean(data)
ss<-sum((data-m1)^2) #偏差平方和
sx2<-ss/length(data) #標本分散
sx<-sqrt(sx2) #標本標準偏差
t値の計算 [2014-05-30修正]
２つの母集団が正規分布に従い、分散が等しいとき、次のように計算できる。ただし、data1, data2はベクトル。
n1<-length(data1)
n2<-length(data2)
m1<-mean(data1)#グループ１の平均
ss1<-sum((data1-m1)^2) #グループ１の偏差平方和
m2<-mean(data2)#グループ２の平均
ss2<-sum((data2-m2)^2) #グループ２の偏差平方和
v<-(ss1+ss2)/(n1+n2-2) #共通の分散
t0<-(m1-m2)/sqrt(v*(1/n1+1/n2))
このt0は自由度n1+n2-2のt分布に従う。
ペチュニアのデータであれば、まずdata1（inflata）とdata2（integrifolia）をデータフレームから次のように取り出す。
inflata<-subset(data0,sp=="inflata")
data1<-inflata$X1 # X1は必要に応じてX4, X30などとする。
integrifolia<-subset(data0,sp=="integrifolia")
data2<-integrifolia$X1 # X1は必要に応じてX4, X30などとする。
qt(面積, 自由度)
指定した自由度のt分布において–∞から積分した面積が指定と一致するようなtの値を計算する。有意水準1%の両側検定なら棄却域は t < qt(0.005, 自由度), t > qt(0.995, 自由度) となる。上側の棄却限界値はqt(0.005, 自由度, lower.tail=F)としても良いし、t分布が左右対称であることから-qt(0.005, 自由度)としても良い。
pt(t値, 自由度)
指定した自由度のt分布において–∞からt値の面積を計算する。pt(t値, 自由度, lower.tail=F)とすると、t値から∞の面積となる。
t.test(変数1, 変数2)
Welchの検定を行う。母分散σ12, σ22とも未知の場合や等分散性が棄却された場合のt検定。p値の解釈（両側なのか、片側なのか）に注意する。何も指定しない時は両側検定として有意水準が示されている。t.test(変数1, 変数2, alternative ="less")またはt.test(変数1, 変数2, alternative ="greater")とすると、それぞれ変数1が2に対して小さい、または大きいという対立仮説のもとでの片側検定における有意水準が表示される。

分散分析に関するRの関数について

by(変数, 因子, 関数)
by関数は指定した変数を因子の水準ごとに分けて関数を適用する。例えば by(iris$Sepel.Length, iris$Species, mean) を実行すると、Sepal.Lengthの平均をSpeciesごと（ここではsetosa, versicolor, virginica）に分けて計算する。結果はリストの１種である、by型の変数で、x <- by(iris$Sepel.Length, iris$Species, mean) とすると x["setosa"]（またはx[1]）にsetosaの平均が、x["virginica"]（またはx[3]）にはvirginicaの平均が格納される。
qf(面積, 自由度1, 自由度2)
指定した自由度のF分布において–∞から積分した面積が指定と一致するようなFの値を計算する。qf(面積, 自由度1, 自由度2, lower.tail=F)とすると、F値から∞の面積となる。分散分析では分母を誤差（級内）変動の平均平方、分子を級間変動の平均平方としたときの分散比を検定統計量F0として使う。帰無仮説のもとでは、水準による変動はないため、F0は0に近くなるはずである。F分布は0から始まるため、分散分析は常に片側検定になる。従って分散分析における有意水準1%の棄却域は F0 > qf(0.99, 自由度1, 自由度2)または、F0 > qf(0.01, 自由度1, 自由度2, lower.tail=F)となる。(自由度1が級間の自由度=水準数–1）
aov(formula, データフレーム)
formulaは観測値~因子A+因子B+因子A*因子Bのように書く。因子A*因子Bは交互作用を示し、因子A*因子B*因子Cのように三因子以上でも良い。誤差はformulaには入れない。これらの名前はデータフレーム中の列の名前を指定し、因子にあたるものは因子ベクトルになっていなければならない。結果は変数に入れておき、summary関数で表示させる。帰無仮説が棄却されたら、次のTukeyHSDで多重比較を行う。
TukeyHSD(aovの結果, 因子名)
aovで計算された結果をもとに、指定した因子の水準ごとの多重比較をTukeyのHSD法によって行う。
glht(aovの結果, mcp=(因子名="Tukey")); cld(glhtの結果)
aovで計算された結果をもとに、指定した因子の水準ごとの多重比較をTukeyのHSD法によって行い、水準間の有意差の有無によってアルファベットをつける。

主成分分析に関するRの関数について

dist(変数)
データは数値行列。次のprincompで使う距離行列を計算する。デフォルトはユークリッド距離。
princomp(変数, 因子, 関数)
データにはあらかじめ計算した距離行列を与える必要がある。また、変数の数がケースの数より多い時はエラーになる。一般的に、次のprcompを使う方がいい。
prcomp(formula,data,scale.=TRUE)
formulaは従属変数を書かずに、分析に使う数値型の独立変数のみを書く。（例　~X1+X2+X3+X4）dataはデータフレーム。計測値の単位が異なるときなど、数値の桁や変動の範囲が変数ごとに大きく異なる場合はscale.=Tで標準化する。scaleのあとのピリオドに注意。デフォルトでは標準化しない。結果はデータフレームとして返され、$xに主成分得点が入っている。

クラスター分析に関するRの関数について

dist(変数)
データは数値行列。次のhclustで使う距離行列を計算する。デフォルトはユークリッド距離。
hclust(距離行列, method="average")
methodはクラスタリングの方法で、指定しなければ"complete" (最遠隣法) になる。その他、"single" (最近隣法)、"average" (群平均法)、"centroid" (重心法)、"median" (メディアン法)、"ward.D" (ウォード法) などが指定できる。

クラスター分析の実行例

> pet<-read.table("petunia.csv",sep=",",header=T)
> pet$X31[37]<-NA
> pet$X31[47]<-NA
> pet2<-subset(pet,complete.cases(pet)) # 欠測値を含む２ケースを除外
> pet_dist<-dist(pet2[,3:23]) # pl, sp, X33, X34を除いたデータで距離行列を計算、[,3:25]とすればpl, sp以外の全データを使う
> cl_res<-hclust(pet_dist,method="average")
> plot(cl_res,labels=pet2$pl,hang=-1)

グラフィックスに関するRの関数について

rglパッケージ
インタラクティブな3Dグラフィックスを提供するパッケージ。このパッケージの関数で描かれたグラフィックスはマウスを使って回転やズームができる。主成分得点の散布図を３つの主成分について同時に見ることができるため、主成分分析で有効に使える。
RSVGTipsDeviceパッケージ
RのグラフィックスをSVG形式で保存するためのパッケージ。パッケージインストーラーでダウンロード、インストール後にパッケージマネージャを利用するかlibrary(RSVGTipsDevice)として読み込む。現在のところ、Windowsバージョンは提供されていないので、linuxもしくはOS Xでのみ利用可能。ただし、RToolsというのを使うとWindowsでも使えるようになるとのこと。（やや古い記事ですが、参照ページ）
devSVGTips(file="ファイル名")
plot関数やhist関数などで通常なら画面に出力されるグラフィックスを指定したファイルにSVG形式で書き込む。
dev.off()
ファイルを閉じて、出力を画面に戻す。