サイトトップ | Software TipsR tips

統計処理ソフトウェアRについてのTips/基本操作

最終更新:2017年 2月 10日 (金曜日)

このページはRの基本操作を説明する。

Topics

なぜRを使うべきなのか?

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で同じことをする場合,サンプルサイズが小さければ,変数xy(変数名は何でも構わない)に直接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を使う上で支障はなくなったといえよう。

最も基本的な操作(2010年1月5日追補)

データの入出力関連

計算テクニック,他

追加パッケージを使うには?
●プルダウンメニューから選んでもいいが,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.