群馬大学 | 医学部 | サイトトップ | Software Tips
最終更新: 2010年 3月 11日 (木曜日) 09時14分
このページでは,国際共同研究のオープンソースなプロジェクトで開発され,GNU GPLに従って公開,配布されている高機能な統計ソフトであるRについてのTipsを扱う。
TOPICS:
参照情報(リンク) | なぜR? | インストール | 基本操作 | 入出力 | グラフィック | 人口ピラミッド | 計算テクニック他 | 数学
保管庫1(2004年1月まで) | 保管庫2(2004年2月から)
保管庫内の主なトピック:
平方和(SS)|「Rによる統計解析の基礎」(保管庫外だがサポート掲示板|正誤表)|オッズ比|レーダーチャート
|1.6.0|1.6.1|1.6.2|1.7.0|1.7.1|1.8.0 | 1.8.1|1.9.0|1.9.1
|2.0.0|2.0.1|2.1.0|2.1.1|2.2.0|2.2.1|2.3.0|2.3.1|2.4.0|2.4.1
|2.5.0|2.5.1|2.6.0|2.6.1|2.6.2|2.7.0|2.7.1|2.7.2|2.8.0|2.8.1|2.9.0|2.9.1|2.9.2|2.10.0|2.10.1
alpha <- function(X) {
dim(X)[2]/(dim(X)[2]-1)*(1-sum(apply(X,2,var))/var(rowSums(X)))
}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ユーザ有志の手によってメッセージまで日本語化されたものも使えるようになった。2008年8月25日に最新版2.7.2がリリースされている。会津大学や筑波大学など国内のCRANミラーサーバから入手することができる。
かつては日本語による解説書があまりなかったが,2003年10月に出版された拙著『Rによる統計解析の基礎』(ピアソン・エデュケーション)を皮切りに,現在では多数出版されているので問題ない。ウェブ上の情報も,群馬大学社会情報学部の青木繁伸教授のサイト内の「Rによる統計処理」や岡田昌史氏による「RjpWiki」を始めとして充実したので,環境としてはRを使う上で支障はなくなったといえよう。
※下記の説明はバージョン1.8.1時点のものなので,別のバージョンをインストールする場合は,適宜読み替えられたい。
all <- list.files(".","\.zip")
install.packages(all, .libPaths()[1], CRAN = NULL)
update.packages(ask=F,destdir=".")
という3行のファイルを作っておいて,そのディレクトリで,"C:\Program Files\R\rw1081\bin\Rterm" --no-save < instlibs.Rを実行すれば,そのディレクトリにあるすべてのパッケージがインストールされ,それより新しいバージョンがCRANにないか確認され,あれば勝手にそのディレクトリにダウンロードされ,アップデートもされるので,古いバージョンを消しておけばいい。
\documentclass{jarticle}
\usepackage{pictex}
\usepackage[dvipdfm]{graphicx}
\begin{document}
\begin{figure}[h]
\centerline{\input{jtext3.tex}}
\caption{3都市の緯度と経度の2次元表示}
\end{figure}
\end{document}
のようにして読み込むファイルtest.texを作って,platex test,dvipdfm testとすれば,test.pdf (16,518 bytes)ができあがる。こちらの方が使えるフォントが多いし,美しくて実用的である(時折そのままでは通らない文字があるのだが,出力されるpictexファイルを直接エディタで編集すれば問題ない。textだけではなく,xlabとかmainでも日本語が使えた)。
sbarplot <- function(freq,xnames,xlabel,title) {
maxc <- max(as.single(dimnames(freq)[[1]]))
z <- rep(0,NROW(xnames))
for (i in c(1:NROW(freq))) {
z[as.single(dimnames(freq)[[1]][i])+1] <- freq[[i]]/sum(freq) }
barplot(z,names.arg=xnames,main=title,xlab=xlabel)}



Package: pyramid Version: 1.1 Date: 2010/1/6 Title: Functions to draw population pyramid Author: Minato Nakazawa <minato-nakazawa@umin.net> Maintainer: Minato Nakazawa <minato-nakazawa@umin.net> Depends: R (>= 2.2.0) Description: Drawing population pyramid using (1) data.frame or (2) vectors. The former is named as pyramid() and the latter pyramids(), as lapper function of pyramid(). License: GPL (>= 2) URL: http://phi.med.gunma-u.ac.jp/swtips/R.html#PYRAMID
\name{pyramid}
\alias{pyramid}
\title{Drawing population pyramid using data.frame}
\description{
Drawing population pyramid using data.frame.
Detailed explanation is given in Japanese
at http://phi.med.gunma-u.ac.jp/swtips/R.html#PYRAMID.
}
\usage{ pyramid(data, Laxis=NULL, Raxis=NULL, Cgap=0.3, Cstep=1,
Llab="Males", Rlab="Females", Clab="Ages", GL=TRUE, Cadj=-0.03
Lcol="Cyan", Rcol="Pink", Ldens=-1, Rdens=-1, main="", ...)
}
\arguments{
\item{data}{A data.frame including left pyramid numbers in the 1st column and
and right pyramid numbers in the 2nd column, where the numbers of males in
each age-class are usually given to left numbers and those of females are to
right numbers. If the data.frame includes 3rd column, it is used as age-class
labels, otherwise the row.names(data) is used as age-class labels.}
\item{Laxis}{A vector of axis for left pyramid. If missing, automatically given.}
\item{Raxis}{A vector of axis for right pyramid. If missing, Laxis is used.}
\item{Cgap}{The width of center gap (as ratio to each panel) to draw age-class. Default is 0.3}
\item{Cstep}{The interval to write the labels of age classes. Default is 1}
\item{Cadj}{The vertical adjustment factor for the labels of age classes. Default is -0.03}
\item{Llab}{The label of the left pyramid. Default is "Males".}
\item{Rlab}{The label of the right pyramid. Default is "Females".}
\item{Clab}{The label of the center age-class. Default is "Ages".}
\item{GL}{Logical value to draw the vertical dotted lines. Default is TRUE.}
\item{Lcol}{The color of the left pyramid. Default is "Cyan".}
\item{Ldens}{The density of hatching lines (/inch) for left pyramid.
Default is -1, when the pyramid will be filled.}
\item{Rcol}{The color of the right pyramid. Default is "Pink".}
\item{Rdens}{The density of hatching lines (/inch) for right pyramid.
Default is -1, when the pyramid will be filled.}
\item{main}{The main title of the pyramid.}
\item{...}{Other options.}
}
\author{Minato Nakazawa \email{minato-nakazawa@umin.net} \url{http://phi.med.gunma-u.ac.jp/}}
\examples{
ages <- c('0-9','10-19','20-29','30-39','40-49','50-59','60-')
males <- c(34,19,11,11,8,7,5)
females <- c(26,25,16,11,7,5,1)
data <- data.frame(males,females,ages)
pyramid(data)
# another example
py.Males <- c(80,40,30,20,10)
names(py.Males) <- c('0-9','10-19','20-29','30-39','40-')
py.Females <- c(60,50,40,30,5)
names(py.Females) <- names(py.Males)
py.df <- data.frame(py.Females,py.Males)
pyramid(py.df,Llab="Females",Rlab="Males",Lcol="navy", Ldens=5, Rcol="red",
Rdens=10, GL=FALSE,main="An example of population pyramid\n with auto-axis")
# The census 2005 for Gunma prefecture data obtained from
# http://www.e-stat.go.jp/SG1/estat/List.do?bid=000001005042&cycode=0
# as "a00411.xls"
GunmaPop2005 <- data.frame(
Males=c(8872, 9144, 9528, 9812, 9817, 10049, 10234, 10047, 10222,
10187, 10319, 10420, 10135, 10473, 10219, 10492, 10943, 11243, 10579,
9653, 9941, 10252, 10396, 10649, 11343, 12031, 12420, 13015, 13354,
14132, 14624, 15873, 16013, 15809, 15338, 14876, 14568, 14536, 14580,
10713, 13844, 12911, 12303, 12158, 12086, 12117, 12620, 12331, 12376,
13138, 14144, 13553, 14353, 15194, 16064, 17286, 18561, 18347, 18230,
12001, 11938, 14363, 13684, 14122, 13475, 12241, 10279, 10866, 10939,
10853, 10093, 9996, 9575, 9334, 9147, 8643, 8249, 7760, 7353, 6812,
6170, 5041, 4099, 3206, 2711, 2778, 2113, 1895, 1630, 1431, 1146, 966,
675, 559, 384, 263, 212, 133, 84, 38, 24, 15, 11, 6, 1, 2, 1, 0, 0),
Females=c(8323, 8750, 8964, 9359, 9559, 9605, 9511, 9800, 9790, 9848,
9939, 9755, 9696, 9884, 9734, 9767, 10293, 10616, 10040, 9383, 9731,
9748, 10211, 10204, 10500, 11086, 11758, 12248, 12548, 13524, 13907,
14930, 15129, 15057, 14344, 13982, 13942, 13587, 13972, 10359, 13212,
12435, 11934, 11673, 11668, 11583, 12100, 11867, 11917, 12746, 13364,
13170, 13968, 15318, 16251, 17125, 18253, 18042, 17927, 11981, 11773,
14450, 14124, 14438, 13502, 12960, 10729, 11710, 11697, 11884, 11413,
11442, 11087, 11035, 11209, 10646, 10482, 9784, 9777, 9491, 8891, 8188,
7636, 7034, 6123, 6103, 4577, 4415, 3861, 3426, 3035, 2571, 2064, 1683,
1209, 878, 739, 536, 333, 193, 134, 86, 56, 28, 13, 8, 5, 3, 1),
Ages=0:108)
pyramid(GunmaPop2005,Llab="Males",Rlab="Females",Clab="",Laxis=seq(0,20000,len=5),
Cstep=10,main="Population pyramid of Gunma Prefecture\n (Data: Census 2005, total by gender)")
}
\keyword{hplot}Dear Sirs, I have uploaded the package named `pyramid' including two functions to draw population pyramid. It's very simple and small. But, as far as I know, there is no such package ever, so that I submitted. I hope that this package can pass the review process and go to main archive. Thank you very much. With Best Wishes, Minato Nakazawa, Ph.D. Associate Professor, Gunma University Graduate School of Medicine

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)
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)
x$AREA <- as.factor(x$AREA)
levels(x$AREA) <- c('A','B','C','D','E','F','G','H','I')とすればいい。出力結果は,[1] 1.644834となる。一方,円周率はpiとして組み込まれているので,[式1]の右辺は,sum(1/(1:10000)^2)
でいい。出力結果は,[1] 1.644934となる。つまり,ちょうど一万分の一の誤差があるわけである。誤差が見えなくなるまで項数を1桁ずつ増やしていったら,千万項が必要なことがわかった。つまり,pi^2/6
の結果は,[1] 1.644934となった。なお,後で青木先生からご指摘を受けた通り,コンピュータでの数値計算の常識として,小さい方から足さないと計算誤差が出るので,足し算は小さい方からやるべきと考えたら,sum(1/(1:10000000)^2)
と書く方がいい。けれども,実はRではどちらでも小さい方から計算されているようだ。つまり,明示的に繰り返しループを書くと,sum(1/(10000000:1)^2)
とした方がいいということだ。この場合はoptions(digits=20)くらいにしないと違いがわからないのだけれども,確かにa <- 0; for (i in 10000000:1) a <- a+1/i^2; a
とすると,結果に微妙な誤差が出る。a <- 0; for (i in 1:10000000) a <- a+1/i^2; a
の出力結果は,options(digits=20)のとき[1] 1.644934066848226となって,完璧にpi^2/6と一致した。sum(1/((1:50)^2*2^(0:49)))+log(2)^2
approxzetatwo <- function(x) { sum(1/((1:x)^2*2^(1:x-1)))+log(2)^2 }と関数定義してから,sapply(1:50,approxzetatwo)とすると,実は41項目で同じ結果が得られていることがわかる。とすれば下図が描ける。ここでは,見やすくするため項数を対数軸表示にした。なお,pi^2/6は,library(gsl)の後でなら,zeta(2)と置き換え可能である。参考までに書いておくと,ここまでの内容は,いつの間にかMaximaと融合し,Gnuplotも内蔵した高機能なものになったEulerというソフトでも,同じような感じで実行できる。plot(cumsum(1/(1:1000)^2),type="l",log="x") lines(c(1,1000),c(pi^2/6,pi^2/6),col="red")

pow<-function(x,y) x^y
plot(2:100,rowSums(sapply(1:10000,pow,-(2:100))),type="l",log="x",xlab="ζ(n)の引数n",ylab="ζ(n)")
require(gsl)
lines(2:100,zeta(2:100),col="red",lty=2)
legend("topright",lty=c(1,2),col=c("black","red"),c("近似","ζ"))
source("http://phi.med.gunma-u.ac.jp/swtips/Zeta.R")
AbsZeta <- function(x) Mod(Zeta2(1/2+x*1i))
curve(AbsZeta(x),0,50)
lines(c(0,50),c(0,0))とすると,条件判定のところで警告メッセージが出る部分があるが,下図ができる。残念ながら計算精度が低いのか,非自明なゼロ点がきちんとゼロになってくれないのだが,だいたい近い形のグラフにはなる。
Correspondence to: minato@phi.med.gunma-u.ac.jp.