Rの最大の利点は,オープンソースなフリーソフトで,かつ拡張性が高い点である。世界中の研究者がGISを含む空間統計解析やゲノム解析などに至るまでさまざまな追加ライブラリを公開しているしCRAN(http://cran.r-project.org)というサイト(筑波大ミラー)に集積される仕組みもあるので,CRAN内で検索すれば大抵の処理は見つけることができる),自分で新しい拡張関数を付け加えることもできる(群馬大学社会情報学部の青木繁伸教授のように,検索するより作ってしまう方が早いと言われる方も多いが)。
オープンソースということは,誰でもその気になれば計算の中身をインプリメンテーションレベルでチェックできることを意味する。これは,商用ソフトにはありえない利点である。商用ソフトでは,利用している計算式はわかっても,コードそのものは公開されないために,インプリメンテーションにバグがないかどうかは,サンプルデータ解析を実行してクロスチェックすることでしか確認できない場合が多い。Rの場合は,世界中の人が使いながらコードチェックもしているので,計算の信頼性もかなり高い。
また,R自体の開発は,メーリングリストで連絡をとりながら行われているが,コアグループが受け入れないとバージョンアップに取り入れられない。コアグループは,かなりはっきりした思想をもってコード開発を行っているので,コードが入り乱れてしまうことはない(これは,裏を返せば,コードをチェックしておかしいと思われる点を報告しても,開発コアグループを説得できないと「仕様だ」ということになって,バージョンアップに取り入れられない場合もあることを意味する)。
プログラムコードが洗練されているため,ソフトウェアのサイズが小さく,動作が軽快なのも利点である。大学等での統計教育を考えると,完全に無料で利用できるため,卒業後も覚えた技術が無駄にならない(SASやSPSSのような高価なソフトウェアで実習を受けた場合,卒業しても使えるような環境に就職できる人はほんの一握りなので,多くの学生はその技術を活用することができなくなってしまう)。教育面では,多くの有用なサンプルデータが含まれているため,いちいちデータ入力をしなくても,分析手法に合わせて適切なデータを利用できる点も利点といえる。Fisherのirisデータのような有名なものは当然として,白血病患者の生存時間データとしてよく使われるGehanのデータも,MASSライブラリ(RecommendedなライブラリなのでWindows版バイナリディストリビューションには標準で含まれている)にgehanとして入っている。
また,国際協力などの場面でもライセンスを気にすることなく共用することができる。英文のみならず,仏文,西文などのマニュアルも公開されている。Windowsだけでなく,MacintoshでもLinuxでもFreeBSDでも動作するので,さまざまな環境で同じ統計解析を行なうことができる。R以上に各国語対応している統計解析のフリーソフトウェアとして,米国CDCが提供しているEPIINFOがあるが(ただしEPIINFOはWindowsのみ),利用できる統計解析手法の種類はRの方がずっと多い。RにはSPSSでさえ実装されていないような新しい分析手法が多く含まれている。結果の信頼性も高く,最近では多くの学術論文が統計解析にRを使っている(例えば,2004年2月にはNatureにもRを使って分析された論文が掲載されている。Morris RJ, Lewis OT, Godfray HCJ: Experimental evidence for apparent competition in a tropical forest food web. Nature 428: 310-313, 2004.)。
もちろん,使いやすさもかなり高い水準にある。例えば,Excelで独立2標本の平均値の差の検定をするには,2群の標本データをワークシート入力し,メニューのツールから分析ツールを選んで(パッケージのインストール時にアドインとして分析ツールをインストールしておかないとメニューに出てこない),等分散でないときの2群の平均値の差の検定を選んで2つの標本の範囲を選んで実行しなくてはならない。結果は別々のシートに出力されるが,余分な統計量が雑然と並んでいて,表の体をなしていない。
Rで同じことをする場合,サンプルサイズが小さければ,変数xとy(変数名は何でも構わない)に直接2つの標本データを付値(代入)してから,t.test(x,y)一発でWelchの検定が完了する。
独立2標本の平均値の差の検定は,古典的には,まずF検定をして2群の分散に差がないかを調べ,差がない場合は通常のt検定,有意差があったらWelchの検定と2段階でなされてきたが,奥村先生や青木先生のシミュレーション結果を見ると常にWelchの方法が最良の結果が得られることが明白なので,この辺りの記述を改めた。なお,古典的な2段階式の場合,Excelではいくつもの手順が必要だが,Rでは関数として
a.t.test <- function(x,y) {
if(var.test(x,y)$p.value < 0.05) {t.test(x,y)} else {t.test(x,y,var.equal=T)}}
と定義しておけば,次からはa.t.test(x,y)とするだけで,自動的に等分散性の結果で条件判定して適切な検定をさせることができる。もっといえば,
t.test(x, y, var.equal=(var.test(x,y)$p.value>=0.05))
とすれば,関数定義さえ必要とせず2段階式の検定ができる(頭の固い雑誌のレビューアがどうしてもやれと言ってこない限り,やる必要はないが)。もちろん,表形式のデータを読み込んでから,分析する関数だけ指定することもできる。
ヘビーユーザーにとっての利点は,RがS言語にほぼ互換な言語のインタープリタであり,それが関数型言語だということから生まれる。つまり,結果を変数に代入して保存したり加工したりできる。xtableというライブラリをインストールしておけば,結果をxtable( )の括弧内に入れるだけで,HTML形式やLaTeX形式に変換でき,きれいな表が得られる。
さらに素晴らしいことに,そうやって実行したすべてのプロセスを,テキストファイルとして記録し,保存しておけるので,後になって,どういう分析をしたかをチェックすることができる。しかも,保存しておいたファイルは(例えばtest.Rというファイル名だとすると),source("test.R")とすると再実行できる。どんなに複雑な作業をしても,それを何度でも簡単に再現できるということである。逆に考えれば,適当なテキストエディタでプログラムとしてRのコマンドを書き連ねておいたものを読み込ませれば,複雑な分析手続きでも1回の操作で終わらせることができる。R Commander(Rcmdr)というライブラリを使うと,メニュー形式で操作することもある程度可能になるが,その場合でも,ログはきちんと関数リストとして保存される。
美しい図を作るのも実に簡単で,しかもその図をpdfとかpostscriptとかpngとかjpegとかWindows拡張メタファイルの形式(emf)で保存でき(ただしemfで保存できるのはWindows版のRのみ),他のソフトに容易に取り込める。例えばwin.metafile()関数を使ってemf形式で保存すれば,Microsoft PowerPointやOpenOffice.orgのDrawなどの中で,ベクトルグラフィックスとして再編集することが可能である。
以前は,多くの日本人にとって最大の難点は,日本語が使えないことだったが(データとしては入っていても大丈夫だったが変数名に使えなかったしグラフ内で使えなかった),中間栄治さんが日本語も扱えるようにするパッチを開発して公開されたのでこの問題は解決した。バージョン2からは本体が国際化対応したので,日本のRユーザ有志の手によってメッセージまで日本語化されたものも使えるようになった。2010年10月16日現在の最新版は2.12.0である。
かつては日本語による解説書があまりなかったが,2003年10月に出版された拙著『Rによる統計解析の基礎』(ピアソン・エデュケーション)を皮切りに,現在では多数出版されているので問題ない。ウェブ上の情報も,群馬大学社会情報学部の青木繁伸教授のサイト内の「Rによる統計処理」や岡田昌史氏による「RjpWiki」を始めとして充実したので,環境としてはRを使う上で支障はなくなったといえよう。
- 起動は,Windowsではデスクトップにできるアイコン(またはスタートメニューのプログラムのRにできるアイコン)をクリックするだけでいい。Windows2000なら,コマンドラインでもRterm --no-saveとして起動できる。LinuxやFreeBSDのシェルをtelnetやsshで使う場合は,Rと打てば良い。いずれの場合でも>というプロンプトが表示されて入力待ちになる。
- 終了は,プロンプトに対してq()と打てば良い。コマンドラインパラメータとして--no-saveなどとつけて起動した場合以外はワークスペースを保存するかどうかの問い合わせがあるので,その回のセッションを記録しておきたいならばyを,そうでなければnと打つ。ワークスペースを保存しすぎると.RDataというファイルが大きくなって起動が遅くなるが,作業中は便利な機能である。
- 基本的に,関数にデータを与えて得られる結果を表示したり,変数に付値したりして使う。アルファベットとドットからなる文字列は変数になりうる。付値とは,ほぼ代入を意味する。例えば,xという変数に3, 5, 7という3つの値からなるベクトルを付値するには,x <- c(3,5,7)とする。これら3つの値の平均値を得るには,mean()という関数を使って,mean(x)とすれば良い。付値せずに関数だけを打てば結果を表示するが,もちろん関数の値を別の変数に付値することもできる。例えばy <- mean(x)として,yをまた別の計算に使うこともできる。変数の情報を見るにはstr()という関数が便利である。
- GUI環境では,HelpメニューからR Manual(html)を選べば,階層構造で説明を参照できる。関数へのインデックスもある。
- 関数の使い方を忘れたときは,help(関数)とかhelp.search("キーワード")で説明が得られる。
- example(関数)で関数の利用例が得られる。
- Rは起動時に作業ディレクトリの.RDataというファイルを読む。これは,Rを終了するときに,"Save workspace image?"(日本語メッセージでは「作業スペースを保存しますか?」)と聞かれたときに「はい」を選んでおけば,そのときに保存される作業内容である(実は,.RDataという拡張子がRに関連付けられるので,別名で保存しておくこともできる。.RDataという拡張子をもつファイルからRを起動すれば,そのファイルに保存されている作業スペースが復元された状態でRが起動する)。つまり,必要な関数を読み込ませた作業後,終了するときに,ワークスペースイメージを保存しておけば,次回同じショートカット(起動アイコン)からRを起動するときには,その関数は使えることになる。欠点は,こうやってワークスペースイメージを保存し続けると,ファイルサイズが大きくなるので,起動にかかる時間が長くなることだろう。
- Rprofileの使い方:Rは起動時にRprofileというファイルを読んで実行するので(Windows版ではRがインストールされているディレクトリ内のetcというサブディレクトリ内にある),その中で,
- 関数を直接定義する
- source()で関数を定義してある外部プログラムファイルを読む
- libraryになっていれば、library()する
とすれば,.RDataを肥大させずに起動時から自作の関数を使えるようにできる(この方法は,当時ハーヴァード大学にいらした林啓一さんにお教えいただいた)。
- foreignというパッケージが標準で入っているので,SPSSやSASのデータも読める。例えば,hogehoge.savというSPSSのデータファイルを読んで,hogehoge.txtというタブ区切りテキストファイルに変換して書き出すには,x<-read.spss("hogehoge.sav")としてxというデータフレームに読み込んでから,write.table(x,"hogehoge.txt",sep="\t")とすれば良い。
- 普通に変数名付きのタブ区切りのテキストデータを読み込むのは,x<-read.delim("hogehoge.txt")で良い(より丁寧な説明)。ただし,read.delim()関数やread.table()関数はオーバーヘッドが大きいので,データサイズが大きいときはscan()関数を使わなければならないかもしれない。
- グラフィックウィンドウで右クリックメニューから保存できるファイル形式はWindowsメタファイルとpostscriptしかないが,FileメニューのSave asからだと,pdfやpngが直接選べる。
- グラフィックを書き出すファイル名をプログラム中で指定することもできる(ディスプレイには画像は出ない)。例えば,作業ディレクトリの"test.pdf"という名前のファイルにpdf形式で書き出すには,まずpdf("./test.pdf",horizontal=FALSE)としてデバイスを開き,次にpar(ps=18)などとして全体のグラフィックオプションを指定し(この例ではフォントのポイントサイズを18ポイントに指定している。グラフィック上での日本語を正しく扱うための注意事項も参照されたい),それからbarplot(x,main="test bar plot")などの描画関数を書き連ね,描画が終わったらdev.off()すればデバイスが閉じられる(pdfデバイスに出力する場合だったら,閉じるまでの間は,複数の図を書き出すこともできる。勝手に改ページされる)。Win32環境で,他のソフトで編集するためにはメタファイル形式で書き出すと便利なので,win.metafile("./test.wmf")などとしてウィンドウズメタファイルデバイスを開くことが多い。画面にも表示するには,windows()デバイスやx11()デバイスに描いてからdev.copy()という方法が使える。
- 追加パッケージを使うには?
- ●プルダウンメニューから選んでもいいが,library(パッケージ名),あるいはrequire(パッケージ名)とする方が簡単である。library()とrequire()の最大の違いは返り値。
- ●未インストールのパッケージをCRANからインストールするには,install.packages("パッケージ名". dep=TRUE)とするのが基本。既にインストールしてあるパッケージを最新版に更新するには,update.packages(ask=FALSE, checkBuilt=TRUE)とするのが簡便である。
- ●CRAN以外に存在するパッケージをインストールするには,レポジトリを適切に指定すること。
- ●既にインストール済のパッケージのバージョンなどを確認するには,packageDescription("パッケージ名")またはlibrary(help=パッケージ名)とする。
- 標準化偏回帰係数を得る
- ●ベクトルとして計算させると楽である。例えば,従属変数をy,独立変数をx1, x2, x3とすれば,
- res <- lm(y~x1+x2+x3)で線型回帰を行った後,
- sdd <- c(0,sd(x1),sd(x2),sd(x3))として各独立変数の不偏標準偏差ベクトルを作り(0は切片用),
- stb <- coef(res)*sdd/sd(y)として,偏回帰係数ベクトルに不偏標準偏差ベクトルを掛けて,従属変数の不偏標準偏差で割ってやれば,stbに標準化偏回帰係数のベクトルが得られる。
- 変数選択
- ●Kendallも書いているように,線型モデルにおける機械的なステップワイズの変数選択は,あまり薦められる方法ではなく,できれば,それ以前に先見的知識から独立変数群に入れられる変数を選んで作ったモデルをそのまま解釈するか,またはSASでいうMAXR法のように,K個の独立変数群候補の中からJ個(J<K)の,R2を最大にする変数の組を総当りで計算し,尤度比検定によって尤度が有意に変わらない最小の独立変数群を選ぶとかすべきと思う。
- ●しかし,いろいろな事情によって,変数選択をしなければならないこともあるだろう。そういう場合,Rでは,標準で組み込まれているstep()関数を使うことで,AICによる変数選択をしてくれる。step()に渡すモデルはlm()でもglm()でもよい。directionというオプションで,増加法(direction="forward"),減少法(direction="backward"),増減法(direction="both")を指定することができるが,scopeとして1つの式(upper扱いされる),あるいはupper(独立変数候補のすべて)やlower(常にモデルに入れる独立変数)からなるlistを指定すれば変数増減法(あるいは減増法)になるし,scopeを指定しなければ(その場合,モデルとして与えた式そのものをupperのscopeと解釈してくれる)変数減少法になる。MASSライブラリに入っているstepAIC関数は,正確なAICを計算してくれるとのこと。
- 分類変数(因子,factor)について
- ●身長heightを150cmから180cmまで5cm刻みのカテゴリ変数hcに代入したいとき,
- hc <- cut(height,seq(150,180,by=5))とすればよい。例えば最小区間が(150,155]となる。ただし,このままでは分類変数(factor)に使えないので,
- hcf <- factor(hc)とすれば,hcfは分類変数になる。また,hco <- ordered(hc)とすれば,hcoは順序のある水準になる。
- ●この場合,lm(y~hc)とlm(y~hcf)とlm(y~hco)はすべて結果が違うし,plot(y~hc)ではドットがプロットされるが,plot(y~hcf)とplot(y~hco)では層別箱ヒゲ図になる。
- ●2つの要因で層別した同時散布図をcoplot(y~x|a*b)によって実行する場合,aやbが分類変数でないと層別数はaについてもbについてもデフォルトでは6である(numbers=で変更可)。aやbがfactorまたはorderedなら,カテゴリごとにplot(y~x)がなされる。
- バッチ実行(Windows2000環境の場合)
- ●Rterm.exeコマンドにリダイレクト入力することが必要。つまり,例えば,"C:\Program Files\R\rw1081\bin\Rterm.exe --no-save < %1"という内容のバッチファイルrun_R.batを作って*.Rに関連付けしておいて*.Rをエクスプローラでクリックするとか,このバッチファイルをパスの通ったディレクトリにおいてrun_R foo.Rとすれば,Windows2000の場合,コマンドラインでfoo.Rが実行される。しかし,バッチ処理が終わると同時にRが終了してしまうので,関連付けからの実行でも画面への出力を残しておくには,バッチファイルの2行目に"pause"を入れておく必要がある。あるいは,画面に残さなくてもいいなら,"C:\Program Files\R\rw1081\bin\Rterm.exe --no-save < %1 > _Rresult.txt"としておけば,結果の出力先は画面でなく,カレントディレクトリの_Rresult.txtというファイルになる。グラフィック出力はそれでも残らないが,そもそもwindows()をせずにグラフィック出力をバッチ実行させるとカレントディレクトリのRplots.psというポストスクリプトファイルに書き込まれるので,それを加工してもいいし,画面でみつつpdfで残したければwindows()してからplot()などの描画コマンドを実行し,描画後にdev.copy(pdf,"_Rresult.pdf"); dev.off()するようにバッチ実行するRスクリプト中に書いておけばいい。画面で見なくてもよければ,例えばpdf("_Result.pdf",horizontal=F)とdev.off()で描画コマンドを挟めばよい。
- データの再コーディング
- ●例えば,9つの地域区分がxというデータフレームのAREAという数値変数(値は1~9)に入っているのだが,それを地域名(A~I)がついた分類変数に変換し,かつ3種類の街区(市街地=A,C,G,農村部=B,F,H,工業地区=D,E,I)に区分し直した新しい分類変数REGも使いたいとする。次のコードで可能である([2004年2月18日訂正]……と書いてずっと公開していたが,実はうまく動かなかったことをお詫びする。AREAは数値変数だから分類変数を1つ付値した段階でそれ以降が欠損になってしまうのだった。別の変数を作る限りは大丈夫なので,新しい変数REGを作るだけにしてAREAはそのままにすれば問題なく動く)。
NREG <- c('工','農','市','農','市','市','市','農','農')
x$REG <- rep(0,length(x$AREA))
for (i in 1:length(x$AREA)) {
if (!is.na(x$AREA[i])&(x$AREA[i]>0)&(x$AREA[i]<10)){x$REG[i] <- NREG[x$AREA[i]]}
else {x$REG[i] <- NA}
}
x$REG <- as.factor(x$REG)
- ●(2004年2月18日追記)実はここに載せていたコードではうまく動かないことは群馬大学の青木先生から御指摘いただいたのだが,ついでにもっとシンプルでうまく動くコードをお教えいただいたので許可をいただいて下に掲載する(1行目はデータフレームxに1から9までの整数値をもつ変数AREAを付値している)。ifelse()関数をこんな風に使えるとは驚いた。
x <- data.frame(AREA=1:9)
NAREA <- c('A','B','C','D','E','F','G','H','I')
NREG <- c('工','農','市','農','市','市','市','農','農')
ifelse(!is.na(x$AREA) & (x$AREA>0) & (x$AREA<10), {x$REG <- NREG[x$AREA]; x$AREA <- NAREA[x$AREA]}, x$REG <- x$AREA <- NA)
- ●(2004年3月14日追記)上のAREAの書き換えのように,水準数が変わらずに名前だけ変えるときは,下記のようにlevels()関数を使うと楽である。
x$AREA <- as.factor(x$AREA)
levels(x$AREA) <- c('A','B','C','D','E','F','G','H','I')
- ●(2017年2月10日追記)これを書いた頃はまだRに慣れてなかったのだなあ,と恥ずかしくなる。今なら下のようにする。1行目はサンプルデータ生成のため,1~9の数字をランダムに100個サンプルしてxというデータフレームのAREAという変数に付値。2行目は同じだが,3行目の右辺は下記で十分だ。
x <- data.frame(AREA=sample(1:9, 100, replace=TRUE)
x$AREA2 <- as.factor(x$AREA)
levels(x$AREA2) <- LETTERS[1:9]
str(x)
- ●(2017年2月10日追記2)でも実はもっと簡単で,下記で十分。いきなりlevelsではNAになってしまうので,labelsを使うのが鍵。
x <- data.frame(AREA=sample(1:9, 100, replace=TRUE)
x$AREA2 <- factor(x$AREA, labels=LETTERS[1:9])
str(x)
- Rの文字列操作について(2014年6月1日)
- ●R本体が最も苦手とする処理の1つが文字列操作である。
- ●read.delim()関数などで,文字列をファイルから読み込むとき,通常は自動的に因子型になる。この自動変換をさせないグローバルオプションがoptions(stringsAsFactor=FALSE)である。ちなみに,データフレームでは列ごとに変数の型が異なっていてもいいが,matrix()ではすべての要素が同じ型でなければならない。
- ●既にデータフレームにしてしまった後で因子型の変数を一括で文字列型に変えたい場合は,Stack Overflowのスレッドにあった(bobがデータフレームだとして),
i <- sapply(bob, is.factor)
bob[i] <- lapply(bob[i], as.character)
というコードが便利だと思う。
- ●数値を書式付きで文字列にするのは,C言語と同様の仕様で,sprintf()という関数が使える。表示桁長を見やすく揃えるときも便利。例えば,sprintf("%09d", 4)の結果は,
[1] "000000004"と表示される。
- ●文字列処理関数としてよく使われるのは,paste(),substr(),strsplit()などであるが,あまり機能は充実していない。stringr()パッケージを使うと,例えば,ある文字列に含まれる別の文字列の個数を返すstr_count()関数などが使える。例えば,str_count("abc1234def5432", "4")は,第二引数の文字列が第一引数に2回出現するので2を返す。
リンクと引用について
Correspondence to: minato-nakazawa[at]umin.net.