2012年12月7日金曜日

複数のRやWinBUGSを同時に使用する(多重起動)

RやBUGSで計算時間の掛かる作業をしていると、Rが塞がってしまい、その間にRを用いた他の作業ができなくなります。とくにBUGSは数時間〜数日掛かりの計算が当たり前なので、時間的制約が大きな問題でした…。

ところが、じつはひとつのマシンで複数のRやWinBUGSを同時並行で使用(多重起動)が可能なんですね…今日はじめて知りました。私にとっては衝撃的な事実で、早く知っていれば、いったいどれだけの時間を節約できたのだろうかと愕然としています。パンドラの箱を開けた気分です…

ここで紹介するのはMacOSXでのやり方です。
(2013.07.16 追記:Windows版では単にRのアイコンを何度もクリックすれば多重起動できることを知って愕然…Winユーザにとってはこんなに簡単なことだったのですね)
cf. Mac上でWinBUGSを使用する方法はこちら


以下に手順を示します(原則、どんなアプリでも多重起動できるはずです)。

*常套手段ではないかもしれないので、使用の際には自己責任でお願いします。
*多重起動により"落ちやすくなる"というご指摘もいただいています。

1)アプリケーション → ユーティリティ→ ターミナル.app を選択し立ち上げる。

2)ターミナルの画面上で、$(ドルマーク)のすぐ後ろに下記のコードを打ち込み、リターンキー。

open -n /Applications/R64.app    # 64bit版の場合
open -n /Applications/R.app        # 32bit版の場合

*注:-n の後ろには半角スペースを入れるのをお忘れなく
*ちなみに、単にopen /Applications/R64.app とやると、Rが普通に起動する。

3)Rを立ち上げたい数だけ、同じコマンドを繰り返し打ちリターンキー。これで多重起動します。

* -n を付加すると、二回目以降が多重起動になるということです。

4)さらにWinBUGSを多重起動するには、多重起動させた各RからそれぞれR2WinBUGS経由でWinBUGSを呼び出すだけでOKです。MacOSX版なのでX11上のwine(WindowsアプリをUNIX上で動作させるのに使用しているアプリ)を通じて起動しますが、X11上のウィンドウもRの数だけ開かれます。

# なんだか変な感じです…(戦隊物のようだ!)。下に白のライトが当たっているのが起動中のアプリを示します。ちなみに、XマークはUNIXのGUI環境(?)、X11、こちらもRの数だけ立ち上がっていますが、アイコンは一つだけでした。複数立ち上がっているのは、X11でなくwineだからなのでしょう。

(2013.02.14 追記)ターミナルで単に r と打つとRを起動させることができますが、ターミナルを複数ウィンドウ(または複数タブ)開いて、それぞれでRを起動させれば同じように多重起動できました。


次に気になるのが、複数起動による計算速度の低下ですが、簡単な例で調べてみました。
こちらに挙げているサポート関数を用いて、全く同一のWinBUGSの計算を同時並行で実行した際の計算時間を system.time() 関数で調べたものです。WinBUGSとR2WinBUGS使用のサポート関数wbugsにsystem.time()を内蔵させているので、これを利用しました(cf. 作って放ったらかしていたら昨日プログラミングのミスを発見し慌てて修正、現在どんどん実用開始していますが今度こそノートラブルになったと思います^^;)。

# 使用マシン
iMac 2.93 GHz Core i7 (Quad Core) 、メモリ 8GB 1333 MHz DDR3

# 使用コード
source("WBUGS.R")

X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
data <- list(X=X, Y=Y)
parms <- list(a=0, b=0, tau=1, F, sigma=NA) 

model <- function() {
   a ~ dnorm(0, 1.0E-6)
   b ~ dnorm(0, 1.0E-6)
   tau ~ dgamma(1.0E-2, 1.0E-2)
      for (i in 1:100) {
         Y[i] ~ dnorm(mean[i], tau)
         mean[i] <- a + b*X[i] }
   sigma <- 1/sqrt(tau)
                                  }

wbugs(data, parms, model, n.iter=1100000)

# 同時並行で計算したRとWinBUGSの数をX軸に、system.timeの返り値(秒、userのところの値を使用)をY軸に取り、glm(Y ~ X, family=gausiaan("log")) による推定曲線を描いたものです。計算回数がえらく不揃いですがご容赦、生物データでないのでテキトーです。
# Quad Core なので、4と8で階段上に変化するかと思いきや、意外と直線的。


とはいえ、この程度のデータならメモリ使用量も少ないからなのでしょう。データや推定するパラメータ数によってこの速度は相応に変化すると思われます。


ちなみに、このテスト中のCPUとメモリの使用量はこんな感じでした。
(「ユーティリティ」→「アクティビティモニタ」)
まだまだ余裕がある感じです。ちなみに、うちの環境ではwineのCPU負荷は、計算内容に関わらずいつも約90%台です。Rの使用しているスレッドが8というのは、Core i7だからQuad x 2 = 8という意味でしょうかね。wineがなぜか7ですが(残りの1は??)。

なお、いくつかのマシンで試して気づいたのですが、X11はXQuartzにアップグレードした方がwineのCPU負荷が少し抑えられているように見えます。いや、気のせいか…

ちなみに、アクティビティモニタを見る限り、wineってまだ64bit化されていないようですね。(2013.01.31追記)wine wikiによると、大半のWindowsアプリが依然32bitなので、wineも32bitを標準としているとのこと。それに32bit版のWindows OSってメモリをたったの4GBまでしか認識できないのですね…

MacOSX10.8 (Mountain Lion) でWinBUGSを動かす (Wine, R2WinBUGS使用)

すっぴんのMac(OSX10.8プリインストール)に入れることを想定して手順を挙げておきます。まず、X11とMacPortsを入れ、MacPortsを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。

(OSX10.7、10.6に入れる場合はこちら http://nhkuma.blogspot.jp/2012/12/macosx106-107winbugswiner2winbugs.html

なお、インストール作業はターミナルからUNIXコマンドを打ちながらのもの。sudoなどのコマンドは注意深く扱う必要があるようなので、チャレンジする際には慎重に。またXcodeもソフトウェア開発に使うような類のツールなので取り扱い注意です。参考にする際は、この辺りを理解の上、自己責任でお願いします…。


(2013.5.10 追記 再挑戦してみました。OSX10.8リリースから時間が経過して一部の工程が変化ので改定しておきました)


(下記、Rコードは緑、ターミナルのコードは紫にしてみます)


**********************
0)実行環境
・OSX10.8 (Mountain Lion) 搭載のMac

1)X11 (XQuartz) のインストール:
アプリケーション → ユーティリティ → X11 を開く。じつはこのX11アイコンはダミーである。あほらしいですがこのステップは必要です。

開こうとすると、Mountain LionではX11がプリインストールされていない、というメッセージが出て、インストールする場合のリンク先が示されるので、そこにある XQuartz2.7.2 をダウンロードし、インストールする(XQuartzはX11と同等の機能を持つもののようだ)。

★インストールしたら、いったん再起動しましょう。これをしないと、なぜか後の winecfg がうまく通らずに悩むことになりました。


2)App StoreのApple IDを設定し、App StoreからXcodeをダウンロードする(Xcode は無料)。1.4GBくらいあったので、けっこう時間はかかる。

3)Command Line Tools(gcc)のインストール:Xcodeメニューから環境設定(Preference)を開く。環境設定のパネル上部に並んでいるダウンロードボタンをクリック。そこにCommand Line Toolsがあるのでインストールする。

4) MacPortsをインストールする。OSXのバージョンによって細かくインストールファイルが分かれている。Mountain Lion用の"pkg"インストーラをダウンロード。
http://www.macports.org/install.php # ダウンロードサイト
MacPortsとは、紐付けの多いUNIXのアプリのインストールを補助してくれる便利なツール。


5)インストールしたら「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。
すると、冒頭にこのように出ている。

コンピュータ名:~ ユーザ名$

このドルマーク $ の後にコマンドを打っていく。
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。
(以降、パスワードを入れる作業は説明を省略)

念のため、MacPortsを最新版にしておきましょう(さっそくターミナルで次のコマンドを打つ)。

sudo port -v selfupdate
(もしうまくいかないようならば、/opt/local/bin/port -v 、さらにターミナルを再起動させたらコマンドが通ったという情報をいただいています)

6)下記の1行を打つ(コピペでOK)。Web上では改行されて見えているだろうけれど、この長いコードは1行で打つべきもののようだ。

echo export PATH=/opt/local/bin:/opt/local/sbin:\$PATH$'\n'export MANPATH=/opt/local/man:\$MANPATH | sudo tee -a /etc/profile

これはMacPortsのパスを通す作業であるようだ。

このパスが無事に通ると、下記のメッセージが出る。
export PATH=/opt/local/bin:/opt/local/sbin:$PATH
export MANPATH=/opt/local/man:$MANPATH


7)このコードを通す。こちらも改行無しで。

if [ `sysctl -n hw.cpu64bit_capable` -eq 1 ] ; then echo "+universal" | sudo tee -a /opt/local/etc/macports/variants.conf; else echo "not 64bit capable"; fi

マシンが64ビットなのか、そうでないのかをMacPortsに伝える作業らしい。
たぶんダウンロードするインストーラの選択に必要なのだろう。
+universal というメッセージが帰ってくればOKのはず。

# 5、6のコマンドはこちらのサイト(http://www.davidbaumgold.com/tutorials/wine-mac/)が非常に参考になりました。


8)wineのインストールを試みる
sudo port install wine 


cf. 途中で「"javac"を開くには、Java SE 6 ランタイムが必要です。今すぐインストールしますか?」というウィンドウが現れるので「インストール」をクリック(66.7 MB)


9)依存パッケージのインストール作業が開始するが、数十分程度進んだところで途中で止まってしまった。そういえば、MacPortsのインストールページの下の方に、これもやっておけと書いてあったので実行(Xcodeのライセンス確認)

xcodebuild -license

# "space"keyを押していくと最後までスクロールできる、最後のところで同意を求めてくるので、読んで確認できたら、

agree


10)もういちど sudo port install wineを試みるも、またエラーメッセージが出て進まないので、インストールできなかったというヤツ(libunwind-headers)を名指しでインストールしてみる。

sudo port install libunwind-headers

(2013.5.10 追記 今度は dbus ほにゃららがすでにあるからダメだったといって途中で停止。エラーメッセージに port -f activate dbus を実行してみろとあったので、
sudo port -f activate dbus を実行したらうまくいった。その都度のバージョンや環境に応じて、このように異なるエラーが起こるかもしれない)

# 無事インストールできた。では再び、
sudo port install wine

今度こそ途中で止まらず、依存パッケージのインストールを再開できたようだ。Macの前を離れ、3時間くらい静置しておこう。

# wineがインストール完了した時のメッセージ:
Installin wine @1.4.1_0
Activating wine @1.4.1_0
(中略)
No broken files found.
最後にWineがインストールされ、無事に成功したという意味なのだろう。


11)よく見るとWineのインストール中に、"これはデフォルトでは入らないが、入れる時のコマンドはこれだと目立つように記されていたので入れてみる:
sudo launchctl load -w /Library/LaunchDaemons/org.freedesktop.dbus-system.plist


12)以下のコマンドを打つ。

winecfg

しばらく処理音が聞こえた後、X11からWineの環境設定のようなウインドウが表示される。
単に一番下のOKをクリックすればよい。
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよい。


13)ホームディレクトリ(OSXでいわゆる各ユーザのホーム階層)にWinBUGS14.exeを置いた上で…

wine WinBUGS14.exe

と打つと、ブルースクリーン上に通常のWinBUGSのインストール画面が表れるので、順を追ってインストール作業を進める。
普通に、はい、はい、と次々に進めていけばよい(Windowsでのインストール手順と同様ということ)。
インストールが終了すると画面が閉じる。
(cf. 私は最初、まちがえて.zipバージョンを展開しようとして、中途半端に展開され、patchやkeyを当てることができずに困惑した。皆さまはお間違いなく)


14)WinBUGSをアップデートするためにWineから起動する。

wine ~/.wine/drive_c/Program\ Files/WinBUGS14/WinBUGS14.exe


15)WinBUGSのアップデート手順は、Windowsの場合と同じ。
Fileからパッチやkeyのファイルを開き(ホームディレクトリに置くのが一番わかりやすいだろう。ファイルの種類をTextにすると出現する)、Tool → Decode → Decode all。Decode allしても黙って通されたりする(なんかのディレクトリがないので作ってもいいか?と聞かれる場合もある)。
Key適用後、WinBUGSを終了させる。

cf.
open ~/.wine # WinBUGSのインストールされたディレクトリは不可視フォルダになっているが、このコマンドで開くことができる


16)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。

# R2WinBUGSのインストールをお忘れなく
library(R2WinBUGS)
# 真の値は、a=3, b=2, sd=1
X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
data <- list(X=X, Y=Y)
inits <- function() list(a=0, b=0, tau=1) 
parameters <- c("a", "b", "sigma")

model <- function() {
a ~ dnorm(0, 1.0E-6)
b ~ dnorm(0, 1.0E-6)
tau ~ dgamma(1.0E-2, 1.0E-2)
for (i in 1:100) {
Y[i] ~ dnorm(mean[i], tau)
mean[i] <- a + b*X[i] }
sigma <- 1/sqrt(tau)
}
modelpath <- file.path(tempdir(), "model.bug")
write.model(model, modelpath)

mcmc <- bugs(
data=data, inits=inits, parameters=parameters, model.file=modelpath, 
n.chains=3, n.iter=5000, debug=T,
working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,
WINE="/opt/local/bin/wine", WINEPATH="/opt/local/bin/winepath")

print(mcmc) # ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう

# bugs()の中身はいろいろと設定が必要。
# 下記の伊藤さまの記事が非常に参考になりました、感謝!
http://www001.upp.so-net.ne.jp/ito-hi/stat/winbugs.html

# WINE、WINEPATHを設定しないと下記のようなエラーが出た。
以下にエラー findUnixBinary(x = "wine") : couldn't find wine binary file
追加情報: 警告メッセージ:
1: 命令 'which wine' の実行は状態 1 を持ちました
2: 命令 'locate wine | grep bin/wine$' の実行は状態 1 を持ちました


17)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。今度こそ成功!
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3



P.S.
なお、2のステップで、職場のマシンなので新規に専用のApple IDを作ろうとしたら、なぜか認証メールがGmailに届かず。検索してみると、同様の例が最近よく起こっている模様。iOS vs. Androidバトルの余波か!?それはいいとして、アドレスを別のドメインのものにしたら(yahooにしてみました)、瞬間的に認証メールが来たので解決。

ちなみに、Apple IDを作る際に「支払方法」というのが出てくるが、「なし」を選択するとフリーのままで行ける。私的使用でないPCなどではこの選択肢を選べばよいでしょう。

おまけ:wineは全ユーザ用にインストールされるようだが、wine内にインストールするWinBUGSなどのアプリケーションは各ユーザ毎にインストールされる模様。他ユーザでも使いたい場合はWinBUGSだけはもう一度インストールする必要がある。

MacOSX(10.6, 10.7)上でWinBUGS(Wine、R2WinBUGS使用)


WinBUGSによるMCMC推定用に使っていたWindowsマシンがお亡くなりになり、急遽必要に迫られて(台数も確保したかったので)Mac上で動かすことに挑戦した(GeoBUGSを使っているのでJAGSではダメ。それに8/2現在R2.15.1ではJAGSが動作しない?エラーメッセージを読むとR2.14で使用してくれといっているようだ)。

ネット上でよく見る手順は省略があるのか、実際はもっと沢山のステップを要したし、かなり厄介な作業だった。OSX10.7、10.6の両方で計3台でのインストールに成功して手順も大体分かってきました(OSX10.8でのインストール手順はこちら)。

手順としては、まず、MacPortsを入れ、MacPortsを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。

インストール作業はターミナルからUNIXコマンドを打ちながらのもの。sudoなどのコマンドは注意深く扱う必要があるようなので、チャレンジする際には慎重に。またXcodeもソフトウェア開発に使うような類のツールなので取り扱い注意です。参考にする際は、この辺りを理解の上、自己責任でお願いします…。

(下記、Rコードは緑、ターミナルのコードは紫にしてみます)


****************************
0)準備。今回の使用環境(いずれか):
・2010年頃購入のiMac、MacOSX10.7.4、X11 または XQuartz、Xcode4.3(Apple)、Command Line Tool(Xcodeの追加機能? Xcodeの環境設定 → Downloadsパネルから入手可能)、RパッケージのR2WinBUGS
・2010年頃購入のiMacとMacBookAir、OSX10.6.8、X11 または XQuartz、Xcode3.2.6(Xcode3.2はOSX10.6のinstall DVDのオプションインストールの中にあるはず(アプリのDVDではなく)。さらにソフトウェアアップデートで最新にしています)、なお10.6用のCommand Line Toolは無い(10.7ではXcodeの機能をCommand Line Toolに分割している模様)。

1)まずはMacPortsをインストールする。OSXのバージョンによって細かくインストールファイルが分かれている。
自身のMacのOSに合わせてインストーラを選ぶ。
http://www.macports.org/install.php # ダウンロードサイト
MacPortsとは、紐付けの多いUNIXのアプリのインストールを補助してくれる便利なツール。

2)インストールしたら「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。
すると、冒頭にこのように出ている。

コンピュータ名:~ ユーザ名$

このドルマーク $ の後にコマンドを打っていく。
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。
(以降、パスワードを入れる作業は説明を省略)

念のため、MacPortsを最新版にしておきましょう(さっそくターミナルで次のコマンドを打つ)。

sudo port -v selfupdate

(もしうまくいかないようならば、/opt/local/bin/port -v 、さらにターミナルを再起動させたらコマンドが通ったという情報をいただいています)

3)では、下記の1行を打つ(コピペでOK)。Web上では改行されて見えているだろうけれど、この長いコードは1行で打つべきもののようだ。

echo export PATH=/opt/local/bin:/opt/local/sbin:\$PATH$'\n'export MANPATH=/opt/local/man:\$MANPATH | sudo tee -a /etc/profile

これはMacPortsのパスを通す作業であるようだ。

このパスが無事に通ると、下記のメッセージが出る。
export PATH=/opt/local/bin:/opt/local/sbin:$PATH
export MANPATH=/opt/local/man:$MANPATH


4)その次は、このコードを通す。こちらも改行無しで。

if [ `sysctl -n hw.cpu64bit_capable` -eq 1 ] ; then echo "+universal" | sudo tee -a /opt/local/etc/macports/variants.conf; else echo "not 64bit capable"; fi

マシンが64ビットなのか、そうでないのかをMacPortsに伝える作業らしい。
たぶんダウンロードするインストーラの選択に必要なのだろう。
+universal または not 64bit capable というメッセージが帰ってくればOKのはず。

# ここまでの作業はこちらのサイト(http://www.davidbaumgold.com/tutorials/wine-mac/)が非常に参考になりました。


5)次、MacPorts経由でWineと依存パッケージのインストールを行う。

sudo port install wine 

実行すると、2時間くらい本体+依存パッケージをインストールしているので席を離れて待つべし。

たぶん、一回目はインストールが完全には成功しない(OSX10.7では一発で完了するかもです)。
こんなエラーメッセージが出た。
Error:The following dependencies were not installed:…..
(中略)
Error:Processing of port wine failed
wine本体のインストールまで到達せず。

6)インストール順による食い合わせの問題だと邪推して、再度…

sudo port install wine

…と打ってみたところ、さっきインストールできなかったファイルが一通りインストールできた模様。
今度は10分くらいでインストール終了した。
Installin wine @1.4.1_0
Activating wine @1.4.1_0
(中略)
No broken files found.
最後にWineがインストールされ、無事に成功したという意味なのだろう。


7)(取りこぼし?)下記のコマンドを通す。
sudo launchctl load -w /System/Library/LaunchDaemons/com.apple.locate.plist


8)もうひとつ、取りこぼしを除去(よく見るとWineのインストール中に、これはデフォルトでは入らないから入れたければこのコマンドで入れてねみたいなメッセージが出ている)
sudo launchctl load -w /Library/LaunchDaemons/org.freedesktop.dbus-system.plist
launchctl load -w /Library/LaunchAgents/org.freedesktop.dbus-session.plist


9)次、以下のコマンドを打つ。

winecfg

しばらく処理音が聞こえた後、X11からWineの環境設定のようなウインドウが表示される。
単に一番下のOKをクリックすればよい。
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよい。


10)次、ホームディレクトリ(OSXでいわゆる各ユーザのホーム階層)にWinBUGS14.exeを置いた上で…

wine WinBUGS14.exe

と打つと、ブルースクリーン上に通常のWinBUGSのインストール画面が表れるので、順を追ってインストール作業を進める。
普通に、はい、はい、と次々に進めていけばよい(Windowsでのインストール手順と同様ということ)。
インストールが終了すると画面が閉じる。
(cf. 私は最初、まちがえて.zipバージョンを展開しようとして、中途半端に展開され、patchやkeyを当てることができずに困惑した。皆さまはお間違いなく)


11)次、WinBUGSをアップデートするためにWineから起動する。

wine ~/.wine/drive_c/Program\ Files/WinBUGS14/WinBUGS14.exe


12)WinBUGSのアップデート手順は、Windowsの場合と同じ。
Fileからパッチやkeyのファイルを開き(ホームディレクトリに置くのが一番わかりやすいだろう。ファイルの種類をTextにすると出現する)、Tool → Decode → Decode all。Decode allしても黙って通されたりする(なんかのディレクトリがないので作ってもいいか?と聞かれる場合もある)。
Key適用後、WinBUGSを終了させる。

cf.
open ~/.wine # WinBUGSのインストールされたディレクトリは不可視フォルダになっているが、このコマンドで開くことができる


13)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。

# R2WinBUGSのインストールをお忘れなく
library(R2WinBUGS)
# 真の値は、a=3, b=2, sd=1
X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
data <- list(X=X, Y=Y)
inits <- function() list(a=0, b=0, tau=1) 
parameters <- c("a", "b", "sigma")

model <- function() {
a ~ dnorm(0, 1.0E-6)
b ~ dnorm(0, 1.0E-6)
tau ~ dgamma(1.0E-2, 1.0E-2)
for (i in 1:100) {
Y[i] ~ dnorm(mean[i], tau)
mean[i] <- a + b*X[i] }
sigma <- 1/sqrt(tau)
}
modelpath <- file.path(tempdir(), "model.bug")
write.model(model, modelpath)

mcmc <- bugs(
data=data, inits=inits, parameters=parameters, model.file=modelpath, 
n.chains=3, n.iter=5000, debug=T,
working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,
WINE="/opt/local/bin/wine", WINEPATH="/opt/local/bin/winepath")

print(mcmc) # ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう

# bugsの中身、いろいろと設定が必要。
# 下記の伊藤さまの記事が非常に参考になりました、感謝!
http://www001.upp.so-net.ne.jp/ito-hi/stat/winbugs.html

# WINE、WINEPATHを設定しないと下記のようなエラーが出た。
以下にエラー findUnixBinary(x = "wine") : couldn't find wine binary file
追加情報: 警告メッセージ:
1: 命令 'which wine' の実行は状態 1 を持ちました
2: 命令 'locate wine | grep bin/wine$' の実行は状態 1 を持ちました


14)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。今度こそ成功!
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3


cf. 数時間〜数日かかるような重い計算をさせるときに、X11上で動作中のWinBUGSの上にカーソールをかざすとカーソールが消えたりする。一見落ちているのかと疑ったが、放って置いたらちゃんと動作していて無事に計算終了した。私のところの環境での一例に過ぎないですが、悲観的にならないよう参考情報として残しておきます。


cf. 初回の立ち上げ時(初めて、RからR2WinBUGSを用いてWinBUGSを立ち上げようとする時)などに、X11からのメッセージで"Wine Monoをインストールさせろ"とか、メッセージが出るかもしれません。基本的にそういうのには素直に従っておけばOKと思います。

cf. なお、.dmgで配布しているバイナリ版のWine最新版(Wine1.4.1)を用いたインストールは私は断念しました。R2WinBUGSから立ち上げようとするものの、以下のようなエラーが出て、解決方法も見つけられず仕舞いです。
wine ports cannot install for the arch(s) 'i386' because dependency apple-gcc42 does not build for the required arch(s) by default and does not have a universal variant


cf. 最初、試行錯誤していた時に、OSX10.6.8のマシンにソース版Wine1.4.1をインストールした場合で生じたケース(たぶん、上の手順通りやれば、この問題は生じない):
err:x11drv:process_attach failed to load libX11.6.dylib: dlopen(libX11….というようなエラーが出た。

Sys.setenv(DYLD_FALLBACK_LIBRARY_PATH="/usr/x11/lib")

という、コードを推定実行の前に通すとよい(こちらはRのコードです)。
これは、伊藤さまのこちらのページを参考にさせてもらいました
(該当の説明自体はバイナリ版の箇所です)
http://www001.upp.so-net.ne.jp/ito-hi/stat/winbugs.html
Sys.setenv() という関数など思いつくのは私には不可能でした…


cf. Wineのアンインストール(バージョンを上げるときなどに必要かと)
sudo port uninstall Wine

その際には、いったんWineでインストールしたWindowアプリもディレクトリごと全部消して、再度入れ直した方がよいでしょう。
rm -rf .wine

2012年12月5日水曜日

変更履歴

新ブログへ引越し初日、さしあたりRいろは等、R基本操作&統計講座関連のページを移行しました。

cf. 「こちらのブログに引っ越してきました」http://nhkuma.blogspot.jp/2012/12/blog-post_5.html

Rいろは・第五部:現代統計への導入(GLM、GLMMの基礎)

(cf. Rいろはのトップページ http://nhkuma.blogspot.jp/p/r_5.html

# 注:作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

(注:使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません)

### GLM準備編 ###########
# まずは古典統計を統計モデリング風に解いてみよう
#(統計モデルの構築とパラメータ推定の練習) 

# 回帰の例(yへのxの影響が有意か、ではなく、yへのxの影響の仕方を推定する)
x00 <- c(1:100)
y00 <- rnorm(100, mean=x00) # x00の値を元に正規乱数を発生(平均 x00 → 回帰係数1、切片0)
m00 <- lm(y00 ~ x00) # 古典統計編でも使ったlmという推定関数
summary(m00) # y00を実行する度に値は変わります(母集団からのサンプリングを繰り返しているので)

# Coefficients: # 切片、回帰係数がそれぞれゼロから有意に異なっているかを示す(いわゆる検定ではない)
#                  Estimate   Std. Error t value Pr(>|t|) # 推定値、SE、t値、有意確率
# (Intercept) 0.325584 0.186879     1.742 0.0846 . # 切片
# x00             0.996002 0.003213 310.016 <2e-16 *** # x00の回帰係数(約1)

plot(y00 ~ x00)
curve(coef(m00)[1] + coef(m00)[2]*x, add=T)
# y00 = x00 + α くらいの一次直線になるはず(切片の値は正負にばらつく)
# こうやってデータ(Y)と説明変数(X)の関係を説明するために、係数を求めるのが統計モデリングの基本形(パラメータ推定という作業)


# ANOVAの例
F01 <- factor(rep(c("A","B"), each=20))
y01 <- rnorm(40, mean=rep(c(10,20), each=20)) # A, Bそれぞれ平均10, 20になるように設定する
m01 <- lm(y01 ~ F01)
summary(m01) # 実行する度に値は変わります(母集団からのサンプリングを繰り返しているので)

# Coefficients:
#                  Estimate Std. Error t value Pr(>|t|)
# (Intercept) 10.2747 0.2316 44.36 <2e-16 *** # 切片=F01のAの推定値(約10)
# F01B           10.0316 0.3276 30.62 <2e-16 *** # 切片にこれを足すとF01のBの推定値(約20)
# F01という要因がAかBかによって、このぐらいyの値が異なるという関係を推定した
#(Aのとき10、Bのとき20という、真の値に近い推定が得られた)
# こういう、推定というプロセスは普通のANOVAではやらないので馴染みがないでしょうね。

plot(y01 ~ F01) # 図示する
points(1, coef(m01)[1], cex=3) # 推定値をプロットしてみる
points(2, coef(m01)[1] + coef(m01)[2], cex=3) # Bの方はこういう足し合わせた値になる



######## 乱数で遊ぼう(GLM準備) ######
# 一般化線形モデルの内部で計算してくれていることを理解するための演習、
# けっこう理解は難しいかも…ここは雰囲気を掴んでもらえるだけでも十分かなと思います
sample(1:10, replace=T) # 1~10から適当に数を選択(replace=T→同じ値を繰り返し取り得る)
rpois(10, 10) # 平均10の整数 x 10個(前の数字は個数、後ろの数字は平均lambda(ラムダ)を表す)
rnorm(10, 0, 1) # 平均0、標準偏差1の実数 x 10個
# 適当に値を変えて、コピペ&リターンを繰り返してみよう。
# パラメータは同じ値のままでも、毎回の結果は変わります(それが乱数。母集団からのサンプリングを繰り返していることになる)

### 今度はサイコロ
# まず、各目の出る確率はいずれも 1/6 = 0.1666...
# 1回振って、ある目("6"とします)が出たら"1"その他なら"0"となる
rbinom(1, 1, 1/6)
# これを10回繰り返し、平均を取ってみる = 出る確率
mean(rbinom(10, 1, 1/6)) # 1/6にはならない
# 繰り返しを増やしてみる(サイコロ10万 笑)
mean(rbinom(100000, 1, 1/6)) # ほぼ 0.166 (1/6)
#(繰り返しを増やすと真の値に近くなる:大数の法則、という性質)


### 図示してみよう
hist(rnorm(10000, mean=5, sd=2)) # 正規分布(実数)、平均と分散で分布が決まる分布
hist(rnorm(10000, mean=100, sd=2)) # 平均の大きさにかかわらず、いつも左右対称な性質を持っている
hist(rpois(10000, lambda=5)) # ポアソン分布(正の整数)、lambda: 平均かつ分散(平均=分散、な分布)
hist(rpois(10000, lambda=100)) # 平均値を大きくすると、左右対称に近くなる(正規分布に近くなる)
hist(rbinom(10000, size=10, prob=0.3)) # 二項分布(0~上限の整数)、size:上限数(上限と確率で分布が決まる)
hist(rbinom(10000, size=100, prob=0.3)) # 上限数(試行回数)を大きくすると、やはり左右対称(
正規分布に近くなる)


##### 尤度(そのデータが生じる確率)を計算するツール ########
po <- rpois(10000, 1) # 平均1のポアソン分布に従う乱数を10000個発生させ、poに格納する
length(po[po==1])/10000 # 上記のうち、実際に値1を取った乱数の個数を全個数10000で割る
dpois(1, 1) # 1という値が、平均1のポアソン分布で生じる確率
# 両者の値は、ほぼ同じになるはず。rpois、dpoisの計算していること、分かりますか??


# もう少し進めていきます
dpois(1:10, 5) # 1~10というデータが、それぞれ平均5のポアソン分布で生じる確率を計算
# dpois(1,5); dpois(2,5); dpois(3,5) ;;;; dpois(9,5); dpois(10,5)# を一括で計算するのと同じ意味
hist(rpois(10000, 5), breaks=c(0:16)) # 平均5のポアソン分布を図示
lines(1:10, dpois(1:10, 5)*10000, type="b") # 確率を重ねて図示(平均である約5で最大となる)
sum(dpois(1:10, 5, log=T)) # 扱いやすいよう対数を取り合計するのが一般的(掛け算のlogは足し算になる)
# これが対数尤度、これを最大化させるようなパラメータ(平均値など)を推定する
#(=観察されたデータが生じる確率が最大となる確率分布を推定するということ)


###### パラメータ推定をシンプルな例で試してみます ########
# まず、ツールの紹介
curve(-(x^2 - 4*x + 4), xlim=c(-10,10)) # 高校1年くらいで習うあの曲線
abline(v=2)
abline(h=0) # 見やすくするよう、縦横のラインを追加
Ys <- function(x) -(x^2 - 4*x + 4) # -(x - 2)^2 = 0 を解くと、x = 2 で最大でしたね?
optimize(Ys, interval=c(-5, 5), maximum=T) # これを計算してくれるツール(区間は指定要)
# X=2のとき、最大値0を取ることが確認できましたか?(二行上の計算と同じ)

# これを用いて推定してみます(尤度計算に即した、もう少し複雑な例へ拡張)
# ちょっと単純化して Y = b*X を考えてみる(切片なしの回帰式)
b <- 5 # bの真の値は5とする(5の時に尤度が最大になるはず)
X <- c(1:100) # Xは1~100
Y <- rpois(100, exp(b*X)) # 平均 exp(b*X)、なぜ指数を取るのかは今はスルーします
poi <- function(b) sum(dpois(Y, exp(b*X), log=T))
# bの値を変化させたときの反応を見るため、関数にしてみる(上記のx^2 - 4*x + 4の状態)
# ちなみにlog=TはXの各データの確率の対数を取るという意味
# なぜ?って、同時に各データが生じる確率は掛け算になるわけだけど(Π、パイの大文字)、
# logを取れば足し算に変わるから扱いやすくなるから(対数の性質、大丈夫ですか?)
optimize(poi, interval=c(0,10), maximum=T) # 対数尤度が最大となるbの値を推定
# 便宜上、推定するためのxの範囲は指定しなければならない(0~10にしてみた)
# $maximumの値:poiが最大値となるxの値。ちゃんと真の値5にほぼ一致しているはず
# $objectiveの値:こちらもほぼ0に近い値になっているはず(確率値の対数は0で最大)




######(余力があったら)residual deviance(残差逸脱度)の導出過程 ####
Yp <- rpois(100, 5) # Ypは平均5のポアソン分布、100個分のデータ
sum(dpois(Yp, Yp, log=T)) # "各値を同じ値で表現"モデルの対数尤度(カテゴリーが100個あるANOVAと同じこと笑)
# とはいえ、当然ながら、これが対数尤度が最大になるモデル
# ただし、一切の単純化がされていないので、モデルとしての意味が無い
sum(dpois(Yp, 5, log=T)) # 平均5のポアソン="真の母集団"モデル、真のモデルなのに"各値…モデル"より尤度は低い

# Deviance(逸脱度) = -2*対数尤度
# residual.deviance = そのモデルのDevianceから、"各値…"モデル(対数尤度が最大のモデル)のDevianceを引いたもの
# この例では、
-2*sum(dpois(Yp, 5, log=T)) - (-2*sum(dpois(Yp, Yp, log=T)))



############ GLM 実行に用いるサンプルデータを用意する
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 10
X <- rep(c(0:9), each=10)
f1 <- rep(rep(c(0:3), each=5), 5) # 予備操作
F1 <- rep(rep(c("C","T1","T2","T3"), each=5), 5) # f1の0~3をC~T1,2,3に読み替え
Y1 <- rpois(100, lambda=exp(-2 + f1 + 0.5*X)) # 切片-2(=C), Tは1ずつ増加、傾き0.5、の乱数発生
Y2 <- rbinom(100, N, prob=logistic(-10 + f1 + 2*X)) # logisticの中身を基にして、二項乱数を発生
vol <- 10*(X+50) + abs(rnorm(length(X))) # 下記のoffset編で用います。数字に特に意味はなし
Dg <- data.frame(X, F1, Y1, Y2, vol)
# 実行する度に値は変わります(母集団からのサンプリングを繰り返しているので)



############ GLM(ポアソン、ANCOVAタイプ) #######
curve(exp(3*x+1), lwd=4) # ポアソン分布のGLMで推定される曲線を図示
# exp()内の値を変えていじって、曲線の性質に慣れてみよう

# 上で定義した値を用い、切片-2, treat+1, 傾き0.5、に戻れるか否かを試してみる
m1 <- glm(Y1 ~ X + F1, family=poisson(link="log"), data=Dg)
summary(m1) # 解析結果の全貌

# Call:
# glm(formula = Y1 ~ X + F1, family = poisson(link = "log"), data = Dg)
#
# Deviance Residuals:
#  Min        1Q       Median   3Q     Max
# -2.0815 -0.8012 -0.3862 0.6185 3.0358
#
# Coefficients:
 # 切片と回帰係数の推定値、SE、z値(Wald統計量)、P値
#                  Estimate Std. Error z value Pr(>|z|)
# (Intercept) -2.32878 0.16594 -14.034 < 2e-16 *** # 切片の真の値は-2
# X                0.51697 0.01138   45.436 < 2e-16 *** # Xの係数の真の値は0.5
# F1T1          1.18029 0.16677     7.077 1.47e-12 *** # 真の値は1
# F1T2          2.15645 0.15124   14.258 < 2e-16 *** # 真の値は2
# F1T3          3.14441 0.14816   21.222 < 2e-16 *** # 真の値は3…に近い値となりましたか?
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # 各要因が独立の時、ゼロと異なる有意確率
#
# (Dispersion parameter for poisson family taken to be 1) # ポアソン分布なので平均=分散になるべきことを意味
#
# Null deviance: 6470.285 on 99 degrees of freedom
# Residual deviance: 96.548 on 95 degrees of freedom # ここが1対1に近ければ推定はたぶん大丈夫
# AIC: 442.39 # このモデルのAIC値
#
# Number of Fisher Scoring iterations: 5
#
# どうだろう?おおよそ上記の値を推定できているのではないだろうか?
# 上記の値は個別に取り出すことが可能
coef(m1) # 回帰係数。coef(model1)[1]は切片、[]内の数字で抜き出すものを指定可能
logLik(m1) # 対数尤度
deviance(m1) # residual deviance、小さいほどデータにfitしてることを示す
AIC(m1) # 460.1022。deviance + 2*パラメータ数 ←複雑なモデルほど罰則が追加される
# 乱数を使っているので、実行するたびに値は多少変わるはず


# 比較対照:切片だけのモデル
m0 <- glm(Y1 ~ 1, family=poisson(link="log"), data=Dg)
# ヌルモデル(null model)という。何も要因を入れていない切片のみモデル
AIC(m0) - AIC(m1) # ΔAIC = 6544.604、m1による説明力の向上
# Deviance explained:
1 - deviance(m1)/deviance(m0) # 0.9842949
# R2決定係数に相当。ただしR2ほど明快な理屈はないのでR2ほど市民権はない
# "今回得られたデータ"のバラツキのうち98.4%がモデルで説明できた
#(今回のデータだけでなく同じ条件で得られる未来のデータにも等しく適合すべき)

# グラフの図示
plot(Y1 ~ X, cex=2, pch=21, bg=c(1:4)[F1], data=Dg)
curve(exp(coef(m1)[1] + coef(m1)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示(F1のCは切片の値と等しくなっている。T1~3はCからの増分の値)
curve(exp(coef(m1)[1] + coef(m1)[2]*x + coef(m1)[3]), lwd=5, col=2, add=T)
# 二本目、図示(T2)
for (g in 2:4) # 二本目以降を一括で図示
curve(exp(coef(m1)[1] + coef(m1)[2]*x + coef(m1)[g+1]), lwd=5, col=g, add=T)



############# 注:GLMでの、ポアソン分布や二項分布を使用上の注意 #######
# これらの分布では平均と分散が等しいが、実データはそうはなっていないことが多い
# 平均と分散の大きさの目安は下記から参照することができる:
# summary(m1)の一番下の方、Residual deviance: 104.55 on 95 degrees of freedom
# 両者の数値の乖離が大きいと推定が怪しいという目安になる(Crawley 2005)*。
# 左 > 右の場合、overdispersion(過分散)という。
# *しかし、この判定方法は標本数もかなり多いときにのみ有効
# Crawley MJ (2005) Statistics: an introduction using R. John Wiley & Sons, West Sussex

## 過分散は後述するGLMMやGLM負の二項分布型(glm.nb)などの利用で対処できる
# GLMMによる解消ではサンプル毎に異なるランダム切片を取るモデルにするのが汎用性のある方法(Gelman A & Hill J 2006 Data Analysis Using Regression and Multilevel/HierarchicalModels (Analytical Methods for Social Research) Cambridge Univ. Press Pp 648)。

# cf. 過分散する例
data(CO2)
m11 <- glm(conc ~ uptake + Plant, family=poisson, data=CO2)
# CO2濃度が植物の吸収量によって変化するという温暖化対策になりそうな解析例(否、"Plant"違いだったようです 苦笑)
summary(m11)
# かなりの過分散になる Residual deviance: 5094.7 on 71 degrees of freedom



########## GLM(offset:単位あたりの量の推定)#######
# 単位量あたりのy1の解析、密度のようなもの。整数値データを割り算するべからず。
m2 <- glm(Y1 ~ X + F1 + offset(log(vol)), family=poisson(link="log"), data=Dg)
summary(m2)
plot(Y1/vol ~ X, cex=2, pch=21, bg=c(1:4)[F1], ylab="Density", data=Dg)
curve(exp(coef(m2)[1] + coef(m2)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示
for (g in 2:4) # 二本目以降を一括で図示
curve(exp(coef(m2)[1]+ coef(m2)[2]*x + coef(m2)[g+1]), lwd=5, col=g, add=T)
# ちゃんとプロットを貫くジャストの曲線が描けましたね?

# cf. 成長量(after/before)のような連続量の割り算データは、さしあたり正規分布の対数連結(gaussian(link="log")で解析すればよいだろう。連続量の場合はoffset項を使用せず割り算しても推定の違いは小さく見える
# m21 <- glm(after/before ~ X + F1, family=gaussian(link="log"), data)
# m22 <- glm(after ~ X + F1 + offset(log(before)), family=gaussian(link="log"), data)
# あとのやり方はポアソンと同様、いずれもexp(coef(m)[1] + coef(m)[2]*x + ...で図示できる

# 連続量もきちんとやろうとするといろいろあります。ガンマ分布(family=Gamma)なども選択肢として考えたほうがいいかもしれません。


# 注:正規分布/正規分布のような割り算データが大きな問題になるようなケースもあるようです。しかし、今のところ連続量の比率の一般的な対処法は無い模様。




########### GLM(二項分布)ロジスティック回帰型 ########
logistic <- function(xx) 1 / (1 + exp(-xx))
curve(logistic(1*x-5),xlim=c(0,10)) # logisticの数値を変えて図示してみよう

m3 <- glm(cbind(Y2, N-Y2) ~ X + F1, family=binomial(link="logit"), data=Dg)
# cbind()の部分はこう書くしかないようです、ちょっと煩雑ですが
summary(m3)
# Residual deviance: 54.02 on 95 degrees of freedom
# under dispersionのようです(実例ではあまり出ない)。少し推定がよくないはず
# 二項分布ではなるべくGLMMがお勧め(とくにデータ数が二桁程度の少数の場合は)
# cf. Warton & Hui (2011) Ecology 92: 3-11
# N=20以下の過分散サンプルでGLMの推定がひどいことを示す例あり

plot(Y2/10 ~ X, cex=2, pch=21, bg=c(1:4)[F1], ylim=c(0,1), data=Dg)
curve(logistic(coef(m3)[1] + coef(m3)[2]*x), lwd=5, col=1, add=T)
# 一本目、図示
for (g in 2:4) # 二本目以降を一括で図示
curve(logistic(coef(m3)[1] + coef(m3)[2]*x + coef(m3)[g+1]), lwd=5, col=g, add=T)
# 図示すると、当てはまりは悪くないように見える



############ GLM(変数が多い重回帰分析タイプ) ######
# まずは使用するデータの下準備
require(car)
data(Chile) # carパッケージに入ってるChileというデータを使用
head(Chile) # Chileの中身を眺める

### 各連続変数の関係を図示してみる(対戦表のような配置、各変数をx, yとしたときの散布図)
pairs(statusquo ~ income + population + age, panel=panel.smooth, data=Chile)
# 相関のある説明変数があるか?(右上がり:incomeとpopulationは怪しい)
# 交互作用のある説明変数があるか?


#### 複数の連続変数が含まれるモデルの場合、相関のある説明変数同士は片方を除く必要がある(多重共線性の排除)
### 基準を用いて対処してみる
# vif (variance inflation factors) の使用(1 / (1 - cor(でーた)) で計算可能)
 1 / (1 - cor(na.omit(Chile[,c("population","income","statusquo")]))^2)

# population income statusquo
# population Inf 1.282471 0.8309652
# income 1.2824706 Inf 1.0403129
# statusquo 0.8309652 1.040313 Inf
# この値が"大きい"変数は多重共線性を生じている。一方を除くと残された方の値も改善する(関連し合う変数は図示して予測する)
# この値が、甘めの基準では、10未満の変数は使用OK。厳し目の基準では1.5未満なら使用OK…というように文献によって基準には大きな幅がある。しかし、この例の場合は厳しめ基準をもクリアしているのでこのまま使用出来る(注:斜めのInfは同じ変数同士なので当然、相関max)

# 説明要因をペアにしたときの相関係数rの絶対値が0.7より小さければ、どの統計手法でもおおよそ大丈夫(Dorman et al. 2013 Ecography 36: 27–46
→(2016.03.11追記)Dorman論文の解釈を再度修正。|r| < 0.7ならば手法による差は小さい一方、それより相関が大きい強い多重共線性においても頑健な手法はない。樹木法なRandom ForestやBRTでさえもしっかり影響を受ける。RidgeやLasso、Elastic netはパラメータチューニングをちゃんとやればいい方法と評価されている。


#### statusquo(満足度?)をincome, population, ageで回帰してみよう
# その前にちょっと説明変数に細工します
# 基準化(scale: 平均ゼロ、分散1)←解析の是非には影響しないが、
# このデータで解析すれば、回帰係数が標準化される、
# すると、絶対値の大きさでその要因による影響の大きさが比較可能になる
Chile2 <- na.omit(Chile)
Chile2[,c(2,4,6)] <- scale(Chile2[,c(2,4,6)]) # income, population, ageだけscale化


### やっとモデリングします
reg2 <- glm(statusquo ~ income + population + age, family=gaussian, data=Chile2)
summary(reg2)

library(MASS) # 雑多な基本パッケージ
stepAIC(reg2) # 簡易的なモデル選択の実行(減少法)
# stepAICは高速で便利だが、下記のglmerには使用不能だし、選択ミスも生じる弱点が指摘されている


library(MuMIn) # もうひとつのモデル選択の…というより…モデル評価のためのパッケージ
options(na.action=na.fail) # これを通しておかないと動作しないことが多い
dredge(reg2, rank="AIC") # 鬼のようなモデル選択(総当たり…)
# rank: AIC 値をよりどころにする; , fixed=c("X1", "X2",...): 必ず含める要因
# 変数の数が多いと組み合わせが多くなって、かなりの計算時間となります。

# AICの値の順に上から下に向かってモデルを列挙します
# (Int): 切片; X(各要因):回帰係数; k:パラメータ数
# delta:ΔAIC、 ベストのモデルからの AIC の増分
# weigtht:Akaike weight
# Akaike weight:第三者に説明するときに分かりやすく便利だが、理論的背景が明確でない指標。0~1の値を取り、値が大きいほどそのモデルの良さのウエイトが大きいことを意味する。ただし"値が1に近いモデルは選択される確率が高い"と言うのは誤りである
# cf. AICc:N数が少ないときのAICとして好む人がいるが、これは正規分布用なので注意(誤用は沢山)
#(cf. Akaike weightやAICcを推進してる文献 Burnham & Anderson 2002


### 主成分分析の成分に変換すると相関関係(多重共線性)は無くなる
Chile3 <- na.omit(Chile) # NAを含む行を削除
pca <- prcomp(Chile3[,c(2,4,6)], scale=T) # 規格化する項が付いている
summary(pca)
# Importance of components:
#                                       PC1     PC2     PC3
# Standard deviation        1.1070 1.0009 0.8790
# Proportion of Variance  0.4085 0.3339 0.2576
# Cumulative Proportion  0.4085 0.7424 1.0000

# PC2まででも、すでに74.24%の情報量がある→PC3は除いてもよさそうだ
PC1 <- pca$x[,1] # PC1
PC2 <- pca$x[,2] # PC2

plot(pca$rotation[,1], pca$rotation[,2], type="n", xlim=c(-2,2), ylim=c(-2,2)) # limは適宜調整要
text(pca$rotation[,1], pca$rotation[,2], colnames(Chile3[,c(2,4,6)]))
# やはりpopulとincomeは類似→より第一軸、age→より第二軸

reg22 <- glm(statusquo ~ PC1 + PC2, family=gaussian, data=Chile3)
summary(reg22)
AIC(reg2) - AIC(reg22) # -109.2854 このケースではPCAの使用でAICが悪化してしまった
# しかし、reg2が多重共線性を含んでいたことを考えると、一概にreg22がダメとも言えない
# けっきょく説明変数が複数のモデリングは一義的に切れないところがある、と思います



### cf. glm(model, family=gaussian) 、正規分布の場合の分散パラメータの取り出し方
summary(reg22)$dispersion # 分散のスケールで出てくる(ルートを取ると標準偏差)



############### GLMM #######
# 一般化線形混合効果モデル(Generalized linear mixed effects model: GLMM)

### あらためてサンプルデータを用意します
# noise(平均0、分散4の正規分布の乱数)を追加
# noiseのせいでバラツキがひどくなるが、平均(切片や回帰係数の値)は変わらない
# 過分散データなどもノイズが多い状態と言えるだろう
# こんなノイズの多いデータでもGLMMなら"比較的"うまく推定してくれる

logistic <- function(xx) 1 / (1 + exp(-xx)) # ロジスティック関数を定義

N <- 10
noise <- rnorm(100, 0, 4^0.5) # SD = √4
x <- rep(c(0:9), each=10)
y1 <- rpois(100, lambda=exp(-2 + 0.5*x + noise)) # 普通のポアソン分布にノイズを追加
y2 <- rbinom(100, N, prob=logistic(-10 + 2*x + noise)) # logisticの中身を基にして、二項乱数を発生
ID <- 1:length(y1)
Dmm <- data.frame(N,x,y1,y2,ID)

library(lme4) # GLMMを実行する汎用パッケージの読み込み


######### GLMM データが正の整数(ポアソン分布)の場合 ###
Introduction to WinBUGS for EcologistsでPoisson-lognormal modelとして紹介されていました。
(引用文献 Millar RB (2009) Comparison of hierarchical Bayesian models for over dispersed count data using DIC and Bayes factors. Biometrics 65: 962–969
# まずは普通のGLMで推定
plot(y1 ~ x, cex=2, ylab="Frequency", data=Dmm)
pois1 <- glm(y1 ~ x, family=poisson, data=Dmm)
summary(pois1) # 推定値は真の値(-2, 0.5)から大きくずれてしまう
# 案の定ひどい過分散 Residual deviance: 5998.1 on 98 degrees of freedom
# 注:乱数を使用しているので値は実行するたびに多少異なるはず

# 図示してみよう:
curve(exp(coef(pois1)[1] + coef(pois1)[2]*x), lwd=8, lty=3, add=T)
# (太い点線)当てはまってるんだか、どうなんだか、図からはよく分からない

### 次に、GLMMによる推定を試みる
# glmerという関数。IDによるノイズを分散値として推定
# ID毎に推定すると、無数の切片を推定することになる→複雑すぎるモデルなので回避する
# 複雑なデータに対処しつつモデルは単純化できる。これがGLMMの使いどころ。だと思う。
pois2 <- glmer(y1 ~ x + (1|ID), family=poisson, data=Dmm) # 警告は出るが推定はうまくいってる
summary(pois2)

# Generalized linear mixed model fit by the Laplace approximation
# Formula: y1 ~ x + (1 | ID)
# Data: Dmm
# AIC   BIC   logLik deviance # AIC, deviance以外にもいろいろ出してくれる
# 315.3 323.1 -154.6 309.3
# Random effects: # ランダム要因の推定結果
# Groups Name Variance Std.Dev. # ここでは分散とSDだけ
# ID (Intercept) 4.4614 2.1122 # IDごとに切片を設け、そのバラツキの推定値(真の値4, 2)
# Number of obs: 100, groups: ID, 100 # IDの数
#
# Fixed effects: # 固定要因の切片、回帰係数の推定
#                  Estimate Std. Error z value Pr(>|z|)
# (Intercept) -2.29482 0.53798 -4.266 1.99e-05 ** # 切片の真の値は-2
# x                  0.51653 0.09214  5.606 2.07e-08 *** # 回帰係数の真の値は0.5
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Correlation of Fixed Effects: # 固定要因間の相関
# (Intr)
# x -0.888

# 今度はけっこう真の値に近い推定となったはず
#(データDmmに乱数を用いているので、データを定義するたびに値は変わります)

# 再び図示してみよう
curve(exp(fixef(pois2)[1] + fixef(pois2)[2]*x), lwd=8, lty=1, add=T)
# (太い実線)xについての回帰曲線を図示
for (i in 1:nrow(Dmm))
curve(exp(fixef(pois2)[1] + fixef(pois2)[2]*x + ranef(pois2)$ID[[1]][i]), lwd=2, lty=3, add=T)
# さらに、あえてランダム効果の数値(ID毎の無数の切片の値)を取り出し図示してみる
# 普通はこういうことはやらない。だからこそ分散の値だけで割愛している。

# 注:データの指定にattach()関数を用いていると、glmerでは読み込みが失敗する模様。面倒でもattach()は使わずにその都度指定しましょう。


##### 負の二項分布型(もう一つのノイズ対処法):ポアソンよりもより歪んだ整数値を取る分布
library(MASS) # に含まれる glm.nb、歪みの程度を変化させて対処するパラメータを持つ
glmb <- glm.nb(y1 ~ x, data=Dmm) # family指定不要
summary(glmb)
# GLMMとの比較記事を書いています。この数値実験結果によれば、どちらを使うのがよいかをAICで選択可能なようです


####### GLMM、二項分布での例 ###
# サンプル毎に異なるランダム切片を取るモデルにするのが汎用性のある方法だと思いますGelman & Hill (2006)でBinomial-normal modelと記述されていました(Gelman A & Hill J 2006 Data Analysis Using Regression and Multilevel/HierarchicalModels (Analytical Methods for Social Research) Cambridge Univ. Press Pp 648)。help(glmer) の実行例でも、このやり方が紹介されています。

# 図示してみる
plot(y2/N ~ x, cex=1.5, ylab="Probability", data=Dmm)
# まずは普通のGLMで推定を試みる
binom1 <- glm(cbind(y2, N-y2) ~ x, family=binomial, data=Dmm)
summary(binom1) # 真の値(-10, 2)から大きくずれた推定となる
# やはり過分散 Residual deviance: 171.34 on 98 degrees of freedom

# 図示してみる:
curve(logistic(coef(binom1)[1] + coef(binom1)[2]*x), lwd=8, lty=3, add=T)
# (太い点線)今度も、ずれていることが見た目では分かりにくいが…

### GLMM、二項分布の実行 ###
binom2 <- glmer(cbind(y2, N-y2) ~ x + (1|ID), family=binomial, data=Dmm) # 警告は出るが推定はうまくいってる

summary(binom2) # 真の値(-10, 2)により近くなったはず
curve(logistic(fixef(binom2)[1] + fixef(binom2)[2]*x), lwd=8, lty=1, add=T)
for (i in 1:nrow(Dmm))
curve(logistic(fixef(binom2)[1] + fixef(binom2)[2]*x + ranef(binom2)$ID[[1]][i]), lwd=2, lty=3, add=T)




# cf. もうひとつのGLMM関数 glmmML
# メリット:高速なモデル選択関数 stepAICが適用可能
# デメリット:binomialとpoissonにしか使えない
# ランダム項は一つしか指定できない(glmerでもなるべく一つにした方がよい)
# 推定計算自体が実行できずエラーになることもしばしば。
(2013.10.18 加筆:設定をいじることでけっこう計算実行は改善される模様、Bootstrapのためのオプションなどもデフォルトで付いているようで良いかもしれません)

Rいろは・第四部:古典統計編


(cf. Rいろはのトップページ http://nhkuma.blogspot.jp/p/r_5.html

# 注:作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

(注:使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません)


###### "古典"統計 or 統計モデリング ######
# 古典統計(ネイマン・ピアソン流?)のメリットは強力な客観性ですが、とにかく「有意差:P < 0.05」さえ出ればよいといった"悪用"が常態化しています(この論文などを参照)。
# この有意確率Pは真の確率を示すものではなく、また帰無仮説と対立仮説が1対1でない限りは、帰無仮説が棄却されたからと言って、対立仮説が採択されることを意味しません。
# Rは古典統計に"冷たく"できているように感じます。例えば、安易な使用すると痛い目をみるanova()関数など(Anova()関数を推奨、後述)。
# 予測を目的とする場合には、GLMなど統計モデリングと情報量規準(AICなど)、あるいはマシンラーニングを用いた方がよいでしょう。
# 応答変数の予測に興味があれば統計モデリング、outputに影響を与えている説明要因の方に興味があれば古典統計を選択するのがよいかと考えています。
# 統計モデリングに進むには→(cf. RいろはのGLM編)

# しかし、古典統計の使用が根強い分野では、たとえ間違っていても"〜ANOVA"という名前さえ付いていれば受け入れられるが"〜モデル"とかいう耳慣れない方法は門前払いというのはやむを得ないことだ思います…"長いものには巻かれろ"というのも無難な選択肢だと思います。
# 本項はそういった"現状"のために最小限必要な手法をまとめておきました。


#### 今回はRに内蔵のデータを使用します ###
# データCO2を呼び出す
data(CO2)


############ ANOVAや回帰を使用する上での前提条件を満たししているかチェックする ######
# 前提条件のチェック→メインの解析、の手順がすでに多重比較にあたるという議論があります。前提条件の検定さえやっておけば、という思考停止にならないためにも、最近は図示をして確かめるようにしたほうが良い気がしています。回帰分析の項に例示したようにplot(lm(モデル))のQ-Qプロットで判断したと記述するのがよい方法だと考えています。

############ 正規性の検定(例:Shapiro-Wilk test)######
# 応答変数(= 従属変数、目的変数)が正規性を満たすかをチェック

# もっとも市民権のありそうな kolmogorov-smirnov 検定
# CO2のuptake、をドルマークで繋ぐことで表現できる
ks.test(CO2$uptake, "pnorm", mean=mean(CO2$uptake), sd=sd(CO2$uptake))
# あまりにもめんどくさい…

# この関数をコピペして利用すると、データ名の記述は一回で済ませられますね。
KS.test <- function(data) ks.test(data, "pnorm", mean=mean(data), sd=sd(data))
KS.test(CO2$uptake) # 使用例

# 代替法:Rで推奨している?、論文でも近年使用例が増えつつあるShapiro-Wilk検定
shapiro.test(CO2$uptake)

# Shapiro-Wilk normality test
#
# data: CO2$uptake
# W = 0.941, p-value = 0.0007908
# 正規分布に対し有意差有り→正規性を満たしていないので、変数変換をする

CO2$uptake.L <- log(CO2$uptake + 1) # 対数変換した列をCO2に追加
CO2$uptake.R <- sqrt(CO2$uptake + 0.5) # 平方根変換した列をCO2に追加

# 再度トライする
shapiro.test(CO2$uptake.L)
shapiro.test(CO2$uptake.R)

# 結局、むしろ正規性が悪化したという珍しいケース…
# 「改善を試みたがダメだった」と、次善の策として元のまま行くしかないだろう
# ちなみに、なるべく図示して、どういうデータなのか視覚的にチェックしてみよう:
par(mfrow=c(1,3)) # 複数の図を行数、列数を指定してまとめるコード
plot(CO2$uptake)
plot(CO2$uptake.L)
plot(CO2$uptake.R)
# 片山状に下がっていっている感じ、たしかに正規分布ではない…


############ 等分散性の検定(例:Levene's test)######
# 残差に対しANOVAや回帰をする→残差に差や傾向がないことをチェックという検定。ただし等分散性のチェックは特にやらないことも多いです。そんなに影響はしません。

library(car) # car というパッケージに入ってるので呼び出す
# 注:初めから入っているパッケージではないので、インストールする必要があります
(cf. 機能パッケージの追加方法)

leveneTest(uptake ~ Plant, data=CO2)
# CO2というデータのuptakeをPlantの群で分割して、という意味

# Levene's Test for Homogeneity of Variance (center = median)
#            Df F value Pr(>F)
# group 11 0.5877 0.8328
# 72
# 今度は有意差無し→分散が異なるとは言えなかった→等分散性は否定されなかった(舌噛みそう…)



############ 回帰分析 ######
# 量的変化に対する量的反応は、ANOVA+多重比較ではなく回帰分析
(xの変数は不等間隔でも等間隔でもOK)

reg1 <- lm(uptake ~ conc, data=CO2)
summary(reg1)
# 切片や回帰係数(Estimateの部分)、
# それらがゼロと異なっているか否かの検定結果、
# R^2値、自由度調整R^2値などを出力

# Call:
# lm(formula = uptake ~ conc, data = CO2)
#
# Residuals:
# Min 1Q Median 3Q Max
# -22.831 -7.729 1.483 7.748 16.394
#
# Coefficients:
#                      Estimate  Std. Error   t value  Pr(>|t|)
# (Intercept) 19.500290 1.853080 10.523 < 2e-16 ***
# conc            0.017731 0.003529   5.024 2.91e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 9.514 on 82 degrees of freedom
# Multiple R-squared: 0.2354, Adjusted R-squared: 0.2261
# F-statistic: 25.25 on 1 and 82 DF, p-value: 2.906e-06


library(car) # carは事前にインストールが必要なパッケージ
Anova(reg1) # 検定にかけてみる
# 回帰係数の分散分析結果が得られる

# Anova Table (Type II tests)
#
# Response: uptake
#              Sum Sq Df F value Pr(>F)
# conc         2285 1 25.245 2.906e-06 ***
# Residuals 7422 82
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


回帰診断をしてみる。
par(mfrow=c(2,2)) # 4枚のグラフが出てくるが、これを並べて表示するためのコード
plot(reg1) # lmモデルをplotすると回帰診断の結果が得られる。

・左上:残差(Residuals)と当てはめ値(Fitted values)、当てはめ値に対して特に傾向がなければOK。この例では少し山なりになっているのが気になる。
・左下:標準化残差との関係、標準化することで残差の傾向の有無がハッキリする。標準化すると左上での傾向は見られなくなった。たぶんOKだろう。
・右上:Q-Qプロット(quantile versus quantile plot)、残差が正規分布に従っているならば直線的に増加する。中央付近はたいてい直線に乗ってくるが、ズレは上下の端付近に生じる。この例は判断に悩むところだが、多少正規性を満たしていなくても大きな問題ではない。
・右下:影響の強い点の存在を示す図。この例では見えていないが、Cock's distanceの点線の外側にはみ出る点があると外れ値の影響が大きいことを示す。

############ 1 way-ANOVA ######
library(car) # ←上で一度呼び出していれば不要なコードです(以下でも同様)
Anova(lm(uptake ~ Plant, data=CO2)) # 上記の回帰分析と中身が同じ構造!

# Anova Table (Type II tests)
#
# Response: uptake
#                 Sum Sq Df F value Pr(>F)
# Plant        4862.2 11 6.569 1.842e-07 ***
# Residuals 4844.8 72
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


############ 多重比較…の前に、ANOVA+多重比較の乱用に対する批判 ######
# 私はポジティブに多重比較を推奨することはできません。そもそも、仮説があれば事前対比をすべきです(でも、逆に事前対比の方が市民権がないような…)。

# ANOVAは質的に異なる群間の違いを解析する検定です(対照vs処理A群、B群…。あるいは特定の生物種間の性質や反応の違いに着目する必要がある場合のような)。

# よく、X軸に段階的に何かの濃度や環境要因を取り、それに伴うYの変化(量的変化に対する量的反応)をANOVA+多重比較で解析している例を見掛けますが、データの情報量の無駄遣いです。間違いとまでは言わないですが。

# ANOVA+多重比較ではよほどの差がない限り隣合うペアの間に有意差は出ません。データ分布の山の裾野が触れるか触れないかぐらいかまたはそれ以上の差が必要です。単に全体としての増加減少傾向を見たい場合は、そんな極端な差があるかどうかを検定したいわけではないはずです。

# 増加減少傾向を示したいのであれば回帰分析でやりましょう。こちらは回帰直線の傾きがゼロと有意に異なっているか否かを検定した上で、傾きが増加なのか減少なのかを回帰係数の正負で示すことができます。多重比較で差が全然出ない場合でも、回帰分析では問題なく傾きが有意となるケースが多々あります。


############ post-hoc test(多重比較)######
# 多重比較:同じデータについて、何回も検定を繰り返すとそのうち有意になる確率は、繰り返す度に増大する。例えば、コイントスを1回やって表が出なくとも、もう1回やったら表が出そう、というようなこと。同様に、検定をもう1回行う場合には、p値の判定の厳しさを倍にしなければならない。つまり、0.05/2 で0.025。3グループあるときは、多重比較を3回する必要があるので、0.05/2/2/2 = 0.00625を判定基準としなければならない。というのが基本的な多重比較の考え方(Bonferroni法)。

# 多重比較については、こちらのページが分かりやすいです
http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html

# Bonferroni法はあまりに厳しすぎるというので(検出力が落ちすぎる)、より緩めたのが下記に述べるような方法。

# Holm法(sequential-Bonferroni法)
# ☆手計算でも簡単にできるし、ベースとなる統計手法を選ばないので適用範囲が広い方法
p.adjust(probabilities, "holm") # 検定手法を選ばない、出力されている確率値を用いて計算する関数
# 使い方:probabilitiesの部分に比較したい確率を入れる、c(0.01, 0.023, 0.031)のような確率値のベクトルを適用
# 返り値は、修正後のより厳しくなった確率値。

pairwise.t.test(CO2$uptake, CO2$Plant, method="holm") # 分散分析の多重比較の場合
# 群のペア間での差異の有意確率を対戦行列で出力する(P値は修正後の値で出力される)

# TukeyHSD法
TukeyHSD(aov(uptake ~ Plant, data=CO2))
# こちらはペアで有意確率(p adj)を算出、差異の推定と95%信頼区間も出してくれる
# ANOVAとともに使用(パラメトリックな方法)

# cf. multcompでTukey法
library(multcomp) # multcompは要インストール
summary(confint(glht(lm(uptake ~ Plant, data=CO2), linfct=mcp(Plant="Tukey"))))
# linfct=mcp(変数名="Tukey")、ここでは変数名がPlantである

# Dunnett法
# 最初の群(普通は対照群)に対する差異を検定する多重比較
library(multcomp) # multcompは要インストール
summary(confint(glht(lm(uptake ~ Plant, data=CO2), linfct=mcp(Plant="Dunnett"))))
# glhtはかなり高機能な関数の模様、help(glht)で説明を見てみるとよいです


############ 多元ANOVA ######
library(car) # carは要インストール
Anova(lm(uptake ~ Type*Treatment, data=CO2))
# "*":単効果の項 Type, Treatment、さらに交互作用の項 Type:Treatment、の全ての組み合わせを表せる

# Anova Table (Type II tests)
#
# Response: uptake
#                      Sum Sq Df  F value  Pr(>F)
# Type               3365.5 1 52.5086 2.378e-10 ***
# Treatment         988.1 1 15.4164 0.0001817 ***
# Type:Treatment 225.7 1   3.5218 0.0642128 .
# Residuals        5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# このAnova関数は多元の時、標準ではType II 平方和で検定する
# →単効果を優先する平方和であり、交互作用は特殊ケースであるという考え方
# 交互作用も等評価する場合は下記のようにtype="III"と明記するとType III 平方和で検定できる
options(contrasts = c("contr.sum", "contr.sum")) # Type IIIの場合はこのコードも通しておく
Anova(lm(uptake ~ Type*Treatment, data=CO2), type="III")

# Anova Table (Type III tests)
#
# Response: uptake
#                      Sum Sq Df    F value   Pr(>F)
# (Intercept)     26217.3 1 409.0389 < 2.2e-16 ***
# Type                  924.0 1    14.4165 0.0002839 ***
# Treatment         134.6 1       2.1007 0.1511413
# Type:Treatment 225.7 1       3.5218 0.0642128 .
# Residuals        5127.6 80
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Type II の時とは異なり、Treatmentが有意でなくなった

# cf. 多元ANOVAは多くとも3元で留めましょう。多すぎると摩訶不思議な交互作用が出ます。CCAやRDAなどの多変量解析が向いていると思います。

# 注:他に、全部小文字の"anova"関数があるが、こちらはType I 平方和に基づく計算をする。多元ANOVAでN 数が不揃いなときには、一般的なType II, Type III 平方和を利用した解析と異なる結果となるので、使用に注意が必要です。
# Type I はRコードで先に入れた要因から順に優先的に有意にする(Y ~ A + B + A:B なら、A, B, A:Bの順に有意になりやすくなる)。なので恣意的な解析結果となってしまう。
# 断りなく"anova"の方を解説している初心者向けテキストもよく見掛けるが注意されたし。そもそも N数が揃っていない解析がANOVAのポリシーに反しているのかもしれませんが…


############ 多元ANOVAの後で多重比較をする場合 ######
# まず、そうすることが本当に必要なのか、研究の目的と仮説を振り返って見直してください。

# 交互作用が有意でなければ、各要因それぞれにつき1wayANOVAと同様に多重比較する
# (今回は交互作用が有意でなかったが、有意だった場合のやり方も示しておく↓)
# まず、両変数をドッキングさせた新変数を作る
CO2$TyTr <- factor(paste(CO2$Type, CO2$Treatment, sep=""))
#(ただの文字列になってしまったので、要因として認識させる:"factor"という関数)
TukeyHSD(aov(uptake ~ TyTr, data=CO2)) # 新変数でTukeyHSD法を実行する


################ repeated-measures ANOVA #############
# 同一サンプルに対し時間的に条件を変化させつつ複数回測定し、グループ間の差を検定したい場合に使用されてきた方法。言い換えると、全体が実験操作など同一条件変化を経験した時のグループによる応答の違いを見たい場合。単に実験前と後との比較のような2段階比較ならば対応のあるt-検定で十分、この方法は2段階より多い場合に使用。
Rで用意された関数を知りませんが検索すると自作関数が見つかります。単にサンプルIDを要因のひとつに加えた2way-ANOVAをrepeated-measures ANOVAとして紹介してる論文も見ましたが、本当にそれでいいのかどうか確信が持てません。
# 必要な仮定などややこしいので、私ならGLMMで対処します。
# repeated measures ANOVA周辺については、こちらのテキストが詳しいです(簡単ではないが)

*追記:repeated-measures ANOVAと下記の対応のあるANOVAはどうやら同物異名らしい。一方で交互作用のないtwo-way ANOVAと同じコードを適用している例もあり、混乱があるようにも思える。いずれにしても反復が多ければ検出力が落ちていくはずで、利用可能な場合にはやはりGLMMを用いるのがよいと思う。

################ 対応のある ANOVA #############
# 同一サンプルに対し複数回測定し、全測定のデータを用いてグループ間の差を検定したい場合。単に同一サンプルに対する複数回観測。
# 一回あたりの反応の誤差や計測の誤差が大きいと思われる場合には有効かもしれません。
# 通常は、ある一時点のデータのみを用いるようにし、「統計解析には○○時間後の値を使用」すると言い切り、普通のANOVAなどでよいと思います。

data(iris) # データirisを呼び出すのに使用(前に実行していれば不要)
iris$id <- rep(1:15,each=10) # サンプル番号を追加(各10回測定)
# 同じSpeciesを繰り返し使用すると見なす。
summary(aov(Sepal.Length ~ Species + Error(id/Species), data=iris))

# 一番下の、Error: Within のところが"見たい"検定結果。

# Error: id
#             Df Sum Sq Mean Sq
# Species 1  52.83   52.83
#
# Error: id:Species
#            Df  Sum Sq Mean Sq
# Species 2 6.679   3.34
#
# Error: Within
#                 Df    Sum Sq   Mean Sq  F value Pr(>F)
# Species        2   4.60      2.3007     8.705  0.00027 ***
# Residuals 144 38.06      0.2643
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# さらにこのあとに多重比較というのは、TukeyHSD()やglht()はコードが通らなかった。seuquential-Bonferroni法を関数を使わずに手動で実行するしかないかもしれません。
# aovはサンプル数が均等であることを前提としている(Type I 平方和)。不揃いの場合は不可と考えた方がよい。
# やはり、第5部で紹介するGLMMを用いた方がスッキリする。サンプルにid付けをし、glmer(.... + (1|id), ...) のように使用すると、idによるバラツキをノイズとして切り分けることができる。


################ nested ANOVA #############
# メインの要因の効果が、バックグラウンドの条件によっても影響されるかどうかを検定したい場合に使用。
# 多元ANOVAとの違いは、バックグラウンドの要因の単効果による影響は見ない点です。
# つまり、メインの要因の単効果と、メインの要因とバックグラウンドの要因との交互作用、の影響を見るのがnested ANOVAです。

#(以下の例では repeated measures…のところのirisとidを使用します)
model.nA <- lm(Sepal.Length ~ Species/id, data=iris)
# スラッシュを使うのがnestedの指定方法
# これはSpecies + Species:id と記述するのと同じことです…メインの要因であるSpecies単効果と、Speciesとidとの交互作用を見ていることになります。

library(car) # car は事前にインストール要
Anova(model.nA) # 今回は単効果はSpeciesだけなので、anova()によるtype I平方和でも構わない。

# Anova Table (Type II tests)
#
# Response: Sepal.Length
#                 Sum Sq  Df    F value  Pr(>F)  
# Species     63.212     2 119.5812 <2e-16 ***
# Species:id   0.896     3     1.1301  0.339  
# Residuals  38.060 144                  

# これでもしSpecies:id に有意差があった場合「Speciesによる違いはidによっても異なる」ことになる。

# repeated-measures ANOVAは、階層的な誤差を用いてANOVAをする。nested ANOVAは交互作用として階層の要因の影響を評価する。という違いと言えばよいのだろうか…。

# 階層が増えるとか、要因が増えるとか、ちょっと構造が複雑になると、モデル式の書き方がややこしくなるし、計算自体がコケることもあります。また、それ以前に、そんな複雑な計算がまともに出来ているとは思えない。
# そんな時は、やはりGLMM。とくにバックグラウンドの要因が複数あるような場合(時間と場所、または場所が沢山あるなど)にはGLMMでないとまともな計算はできないと思います。また、複数の階層にしたい場合は、WinBUGSなどによるMCMC推定をオススメします。

# nested ANOVAは、一段階まではglhtによる多重比較が"計算可能"のようです。それ以上は手動でsequential-Bonferroni法を使用するしかないかもしれません。

# ちなみに私は、複雑な構造のnested ANOVAを頑張って実行しようとして試行錯誤した結果、どうやってもこんな複雑なデータは検定ベースでは解析できない!と断念して、統計モデリングの世界に入った経緯があります…。


################ cf. 共分散分析 ANCOVA #############

################ cf. 重回帰分析 #############
# 重回帰は、前提条件のチェックなどけっこう厄介なので、何でもかんでも変数を放り込んで"レンジでチン"的な安易な適用は危険だと思います。


################ cf. 多変量分散分析 MANOVA #############
# 1つのサンプルに対し、複数項目の応答(呼吸と光合成など)があり、説明変数に対する反応をトータルとして見る場合(変数間の相関(共分散構造)をも組み込んだ解析をしてくれる)。個別の応答タイプに対し個々の検定をするのではなく、検定を一括で1回で済ますことで、第一種の過誤を抑える狙いがある。また逆に、複数変数を同時に解析することにより、単一変数では検出できない効果を検出できる場合もある。

# manova()関数を使用。
# 前提条件は基本的にはANOVAと同様:
# ・応答変数同士の等分散性、個々の変数の正規性が必要(ただし正規性に対しては頑健)
# ・サンプル数が少ないと検定力が落ちる←サンプル数が少ないからMANOVAで複数の応答変数をトータルで見て"ゆ〜い"にしようというスケベ心には答えてくれない

# 関連(相関)のありそうな複数の指標を用いてグループ間の差を見たいケースなどに有効。全体の傾向を見るのに適している。
# MANOVA以外のやり方としては、主成分分析で応答変数を減らした後、成分ごとに解析を行うのも手だろう。共分散構造を用いて変数をまとめるので、主成分と元の変数との相関関係もちゃんと分かります。
# PerMANOVAもオススメ、library(vegan)のadonis関数で提供されています。生物群集の構造(各種の個体数)を応答変数とした解析が主目的のようですが、method="euclid"を用いれば連続変数にも適用できると思います。

Rいろは・第三部:グラフ作成編


(cf. Rいろはのトップページ http://nhkuma.blogspot.jp/p/r_5.html

# (2012.12.27)操作ミスによりアドレスが変わってしまいました…

# 注:作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。

(注:使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません)

# データの準備の仕方からグラフ作成、体裁作りまでを実行します


###### 基本的グラフ関数の紹介 #####################
# まずはシンプルなサンプルデータを使用
x0 <- rep(1:5, each=4) # 1~5を4個ずつ
y0 <- sort(rnorm(20, 10)) # 平均10の乱数を20個、ソートしておく

plot(y0) # 大概のデータは"plot()"に放り込めばとりあえず図示してくれる
boxplot(y0 ~ x0) # 箱ひげ図(~の意味は、y0をx0で分ける感じ → 統計でよく使う)
barplot(y0) # 棒グラフ
hist(y0) # ヒストグラム
plot(y0, type="o") # 折れ線グラフ
plot(y0 ~ x0) # 散布図(plot(x0, y0)と表すこともできる。ただし、先がX軸, 後がY軸になる)
plot(y0 ~ x0, ylim=c(8,12)) # Y軸の上下限を8と12に設定(X軸の場合はxlim)
plot(y0 ~ x0, yaxp=c(8,12,4)) # さらにY軸の目盛数を変更(上下限を合わせて目盛数は+1になる)
plot(y0 ~ x0, log="y") # y軸を対数スケール(log="x": x軸を対数、log="xy": 両対数)
plot(y0 ~ x0, cex=c(1:20)) # 散布図でプロットサイズを変更(cex:サイズ指定)

# 注:統計モデルでよくやるように、plot(y ~ x, data=D)のようにデータフレームDのxとyを用いるように指定できるが、plot(x, y, data=D)とするとデータが見つからずエラーとなるので注意。
# (実は ~ を用いるのはplot.formula()という拡張関数であり、~を用いると自動で呼び出される)

###### Rグラフの取り扱い #################
# "ファイル" → "保存"で取り出せる(PDF形式)
# Win版ではメニューから様々なファイル形式での保存が選べるようになっている

# 保存したグラフをOpen officeで編集するやり方がこちらにあり(Win版の説明)


###### 以下、現実的なデータセットと体裁、グラフ形態にしてみよう ########
## 下記のグラフ作成で用いるサンプルデータを用意する。カテゴリー変数:2要因、連続変数:1要因
f1 <- rep(c("A", "B", "C","D"), each=10)
f2 <- rep(c("a", "b"), 20)
X <- sort(sample(0:100,40,replace=T)) # 0~100の数字をランダムに抜き出しソート(重複あり)
Y <- runif(40,0,100) # random uniformの意味。0~100の一様乱数を40個発生させる
D <- data.frame(X=X, Y=Y, F1=f1, F2=f2) # データフレーム化する

## 棒グラフ、折れ線グラフ用に平均と標準偏差を用意
MEAN1 <- tapply(D$Y, D$F1, mean) # 各グループの平均("対象データ", "グループ分け", "関数名")
SD1 <- tapply(D$Y, D$F1, sd) # 各グループの標準偏差
MEAN2 <- tapply(D$Y, list(D$F2, D$F1), mean) # 各グループの平均
SD2 <- tapply(D$Y, list(D$F2, D$F1), sd) # 各グループの標準偏差


###### グラフの体裁を整える準備(グラフ枠とラベルのサイズを固定するなら) ####
# par() に入れたパラメータの効果はウインドウを閉じるまで持続する
par(lwd=2, # line widthの意味。線の太さを標準の2倍にした
cex=1.5, # 文字サイズの指定。標準の1.5倍にした
mgp=c(1.75, 0.75, 0), # 軸タイトル、軸ラベル、軸線の間隔を指定
pin=c(3, 3)) # 作図領域のサイズをインチ指定

# family="Helvetica" # フォントを指定する項を追加可能
# この例は Helvetica に指定した例(推奨)。ほか Courier, Times など
# Windowsでは業界標準の Helvetica がインストールされていないので、擬似フォントの Arial を指定する
# (cf. ご自身で記述する時はこんなにマメに改行する必要はありません)

## モニタの小さなPCでは画面に収まりきらずにエラーになる場合あり。そうしたらpin=c(,)の値を小さくしてください。

###### 箱ひげ図(1要因)の作成例 ###################### 
plot(Y ~ F1, data=D, # YをF1でグループ分け、データはDを用いる、の意味
lwd=2, ylim=c(0,100), # y-limitの意味。上限と下限を指定(計算式も利用可)
col="gray80", # カラーを指定:80%白の灰色。"colors()" で色一覧を見られる
axes=F) # いったん軸を除去(plotは、目盛りの線太だけは変えられないようなので)
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1)
# 軸の付け直し(二度手間だが、これをすることで目盛りも同じ線太にできる)

### さらに、図のラベルを左上部に追加する(論文に掲載する図を想定)
mtext("A) Coral growth", side=3, # side:3で上部を指定
at=0.8, line=0.5, # at:左からの横位置、line:グラフから上端からの縦位置
cex=1.5, font=2) # font=2:太字(普通でいいならfont項は不要)
# atの値はデータに応じ細かく調整要…正直、ドロー系ソフト上で調整する方が手っ取り早い


###### 棒グラフ(1要因の場合)の例 ########################
bar1g <- (barplot(MEAN1, # 図示したいデータ列
ylim=c(0,100), lwd=2, # 線太を2倍にした
las=1, # Y軸のラベルを水平方向に変更
col="wheat", # カラーを指定
xlab="group", ylab="volume")) # ラベル名を自分で書く
### エラーバー追加(barplot自体が各棒の中心のx座標データを持っているので利用する)
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2) # x1,y1,x2,y2:両端のx,y座標を指定
mtext("B) Coral happiness", side=3, at=1.2, line=0.8, cex=1.5, font=2)

# cf. points(x,y) # x,y座標を指定し点をプロットできる(cex=数字などの項を加え装飾可能)


###### 応用:間隔を変えた棒グラフ(control1群に対し、treatment複数群の場合)#
barVSctrl <- (barplot(MEAN1,
space=c(1,1,0.5,0.5), # space:棒の幅1に対する棒の間隔を数値で指定
xlim=c(0,7), # xlim:x軸の範囲。棒の幅、棒の間隔の合計より少し大きくなるように設定
ylim=c(0,100), lwd=2, las=1,
col=c("white", "lightblue", "lightgreen", "mistyrose"),
xlab="group", ylab="amount"))
segments(barVSctrl, MEAN1+SD1, barVSctrl, MEAN1-SD1, lwd=2)
mtext("C) Fish sadness", side=3, at=1.0, line=0.8, cex=1.5, font=2)


####### 棒グラフ(2要因(2x2グループ))の例 ######################
bar2g <- (barplot(MEAN2, beside=T, # TをFにすると積み重ねグラフになる
ylim=c(0,100), lwd=2, las=1, col=c("white", "salmon"),
xlab="group", ylab="Intensity"))
segments(bar2g, MEAN2+SD2, bar2g, MEAN2-SD2, lwd=2)
mtext("D) Shrimp madness", side=3, at=3.0, line=0.8, cex=1.5, font=2)


###### ヒストグラム ######################
hist(D$Y, breaks=c(0,20,40,60,80,100), # 区切りの指定(しなくてもOK)
main="", # でかでかタイトルを除去
lwd=2, las=1, col="seashell", xlab="velocity")
mtext("E) Mollusca marathon", side=3, at=25, line=0.8, cex=1.5, font=2)


###### 折れ線グラフ ######################
plot(MEAN1, ylim=c(0,100), xlab="factor", ylab="richness",
type="o", # 折れ線グラフ(cf. "p":点のみ、"l":線のみ)
lwd=2, axes=F, cex=3, # プロット(文字サイズcex)を倍にする
pch=23, # プロットのシンボルを指定(21~25では輪郭と塗りつぶしの色が指定可能)
col="black", # 輪郭の色
bg="tomato") # 塗りつぶしの色
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
segments(c(1:4), MEAN1+SD1, c(1:4), MEAN1-SD1, lwd=2) # 今度はx軸は連続変数なので数値で指定
mtext("F) Polychaeta beauty", side=3, at=1.7, line=0.8, cex=1.5, font=2)


###### 散布図 + 回帰直線 ######################
plot(Y ~ X, xlim=c(0,100), ylim=c(0,120), xlab="Factor", ylab="Loneliness",
axes=F, cex=2, pch=21, col="black", # シンボルの形、輪郭の色は全系列で統一した
bg=c("seagreen","turquoise","slateblue","plum")[F1], data=D)
# 塗りつぶし色を系列によって変える("F1"のA~D、4カテゴリーで色分けした)
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
legend("topright", legend=levels(D$F1), pt.cex=2, pch=21, # 凡例を追加
pt.bg=c("seagreen","turquoise","slateblue","plum"),horiz=T)
curve(0.2*x + 30, lwd=4, col="orangered", add=T) # 回帰直線:xの一次関数の例
mtext("G) Damselfish", side=3, at=5, line=0.8, cex=1.5, font=2)


############## 複数の図の重ね合わせ ######################
###### 散布図の複数系列を重ねる(上記とは別のやり方で描いてみる) #########
# まず、枠を用意
matplot(c(0,100), c(0,100), xlab="Factor", ylab="Jealousy", axes=F,
type="n") # type="n":グラフエリア内の余計な文字を消す
box(lwd=2); axis(side=1, lwd=2); axis(side=2, lwd=2, las=1) # 軸の付け直し
# 1系列ずつ描く
matpoints(D$X, D$Y, cex=2, pch=21, col="black", bg="turquoise") # 1系列目
matpoints(D$X*1.1, D$Y*1.1, cex=2, pch=24, col="black", bg="slateblue") # 2系列目
mtext("H) Sea urchin", side=3, at=5, line=0.8, cex=1.5, font=2)



############## 複数の図の並列 ######################
library(lattice) # latticeというパッケージを呼び出す
xyplot(Y ~ X | F1, cex=1.5, data=D) # ぴったり並列して書き出すので見やすい
# ただし、par()が効いていないようだし、plot()とは使い勝手が異なる模様…便利そうなのですが。


###### 異なるタイプのグラフを重ねる ################
# まず、一つ目のグラフ(棒グラフとする)
bar1g <- (barplot(MEAN1, xlim=c(0,5), # x範囲を二つ目のグラフと合致させる
ylim=c(0,100), lwd=2, las=1, col="khaki",
xlab="group", ylab="conflict"))
segments(bar1g, MEAN1+SD1, bar1g, MEAN1-SD1, lwd=2)
par(new=T) # 上書きするコマンド

# 二つ目のグラフ(折れ線グラフとする)
plot(bar1g, # エラーバーと同じ要領でbar1gをx座標に利用
MEAN1*1.5, xlim=c(0,5), # x範囲を一つ目のグラフと合致させる
type="l", # 線だけのグラフタイプを指定
axes=F, xlab="", ylab="", # 枠やラベルを描かないためのコード
lwd=4, col="royalblue")
axis(side=4, lwd=2, las=1) # 右サイドに軸だけを追加
mtext(expression(paste(factor2(m^2))), side=4, line=2, cex=1.5)
# 二つ目のグラフのラベルを右サイドに追加(expression, pasteを組合わせると上付き等も書ける)
mtext("I) Sea cucumber", side=3, at=0.1, line=0.8, cex=1.5, font=2)


####### 画面を分割して複数グラフをまとめる ########################
# parにmfrowを加える
par(mfrow=c(2,2), # multi frame by row の意味。この例では行方向に2x2=4分割
lwd=2, cex=1, mgp=c(1.75, 0.75, 0), pin=c(1.5, 1.5))
# 適当な4枚のグラフを実行してみよう
#(ウィンドウをドラグしてサイズを調整するとグラフの間隔を変えることができる)

## モニタの小さなPCでは画面に収まりきらずにエラーになる場合あり。そうしたらpin=c(,)の値を小さくしてください。


####### 便利技:全変数を一括で図示してみる ########################
# 散布図形式で、全変数の関係を対戦表に配置する
# (自身vs.自身の位置に変数名が入ります。横方向は自身がY、縦方向は自身がXの場合)
pairs(D, panel=panel.smooth) # panel.smoothで変数間の関係を表す曲線が追加される
pairs(Y ~ X+F1+F2, panel=panel.smooth, data=D) # これでも同じ図になる



####### 便利技:色にこだわる ########################
# 色の種類指定が面倒くさいのは難点。せめて色一覧くらいほしいですね。
# htmlの色コードが使用可能:
# plotなどの関数の項で、col="#00FFFF" のように""でコードを囲うと利用できる(水色の例)。

# ちなみに、RGBのhtml色コード末に AA を付加すると半透明になる(例: col="#00FFFFAA")
# cf. 透明度10%の場合、#00FFFF1A。命名ルールは今一つよくわからない…

# なお、グレースケールでよければ、gray()と0~1の小数を用いてgray(0) ~ gray(1)で表現可能
# 0が黒、1が白。なお、この場合はgray() の前後に""は付けない(関数になっている)。



# その他、実用的な色セットを用途別に一覧にしてくれるパッケージがある(インストール必要)。
install.packages("RColorBrewer", dependencies=T)  # RColorBrewer のインストール(初回のみ必要なコード)
library(RColorBrewer)  # RColorBrewerの呼び出し(毎回必要。名前が長い…忘れてしまうorz)
display.brewer.all()  # 色一覧をグラフにしてくれる(増加、カテゴリー、正負の増減、用途別)
# 各行の文字は色セット名を表しており、この単位で色を利用できる
brewer.pal(8, "Greens") # Greensを選択する場合。数字は色数(段階数)
# ただし、各セットに用意されている色数を超えて指定することはできない(cf. 上図)
brewer.pal(8, "Greens")[1] # 実際の使用では[1]のように添字を付けて色を取り出す
paste(brewer.pal(8, "Greens"), "AA", sep="") # こうすると透過色にできる

# ColorBrewerのサイトも有用です、いろいろなケースにおいて勧められる色セットなどのガイドもあり

# cf. いわゆる"jet"color、二次元マップのように違いを明確にしたい用途に適している
jet <-  colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))  

# 散布図などで重なりをより見やすくするには、xをjitter() 関数に入れて微小な"ブレ"を加えてやるとよい

# cf.  RGBaを数値で定義するにはrgb()関数。例えば青の透過色の場合、rgb(0, 0, 1, 0.1)。
# rgb(R, G, B, 透明度) の順になっている。いずれも最大が1になるように設定する。
# 逆に色名や色コードからRGBの数値を得るにはcol2rgb()関数。


###### おまけ:EPSで出力する場合 ######################
# 下記、postscript…を実行してから、通常のグラフ作成コマンドを実行、
postscript("graphEPS.eps", height=4, width=4, # サイズを指定
horizontal=F, onefile=F, paper="special") # ここは暗号…
# ファイル名とheight、width(インチ単位…)をお好みで書き換えてください
# 作成終了したらdev.off()を実行(括弧内は空白)
dev.off() # 作業ディレクトリ内へ保存される