random dispersal

"時間が掛かるからRでforループを使うな"、"applyファミリーを使え（ほかにsapply, lapplyなど）"とよく言われる。かといって無理に使うと難解なコードになる場合もあるし、せっかく実現しても計測してみたら、むしろforループの方が早かったなんてこともあった（今回の動機）。並列化（複数のCPUコアを用いる）するという手も考えたが、大きなー繋がりのタスクを1回やるだけならよいが、入れ子になっていると呼び出す時間の方が律速になったりもする。

applyファミリーやforeachなどで置換えにくい計算例として、値を毎回更新するようなプログラムがあるだろう。例えば以下のような計算例（"0.9* "さえ無ければapply(x, 2, cumsum)で一発なのだが、sapplyとapplyを組み合わせて書いたらforより遅くなった…）。

NN <- 10^6
x <- matrix(0, NN, 10)
# NN回分の結果を格納するための入れ物（1回あたりは長さ10のベクトル）
system.time( for(r in 2:(nrow(x))) x[r,] <- 0.9*x[r-1,] + rnorm(10) )
# ひとつ前の結果に0.9を掛け、正規乱数を足していくマルコフ連鎖
# ユーザ システム 経過
# 5.457 0.363 5.778

手元の環境では6秒弱掛かった。

いろいろ調べて行き着いたのが Rcpp というパッケージ（要インストール）。RからC++（"しーぷらぷら"と呼ぶみたい、なのでcppなのだろう）で作成した関数を呼び出せるというもの。Cで作られているというパッケージは最近よく見かけるようになったが、自分でカスタマイズできるのはすごい。もっとも、C++を書けるならばですが。それでも部分的にだけでも記述できればだいぶ速くすることができるはず。目的のforループ作成を達成するのに丸一日掛かりましたが、最終目的の計算に掛かる予想時間が数ヶ月だったことを考えれば大幅な時間短縮になりました。

手順として、まず開発環境のセットアップが必要です。Macの場合はXcode（AppStoreからフリーでダウンロード）とX11（ユーティリティ内にあるX11のアイコンをダブルクリックして立ち上げてライセンスにOKすれば大丈夫なはず）。Windows版はすみませんよく分かりません…コマンドプロンプトでC++をコンパイルする必要があると思います。
次に肝心の、C++で記述したコードが必要です。ただし随所に純粋のC++と異なる書き方を要する箇所があります（私はここでつまづきました）。C++で記述したファイルをRの作業ディレクトリに置いて読み込むと、Rの関数として使用できるようになります。





まずは最小限のコード例：（以下、C++のコードは紫にしておきます）
require(Rcpp) # Rcppパッケージを呼び出す（この行はRに直打ち）

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double testF(double x) {
return(x*x);
}
紫で書いた6行分をコードエディタなどで書いて test.cpp というファイル名で作業ディレクトリに保存します（ .cpp ってC++ファイルの拡張子なのだと思います）。
C++用コードの初めの3行は決まり文句のようなものです。3行目は用途によって変える場合があるようですが、とりあえず忘れてよさそう。
まずRとの大きな違いはすべての変数型を逐一定義すること。この性質は嫌いじゃないです。むしろ勝手に変数型が変わるのはRが嫌われる点でもありますね。

double は実数の意味、整数なら int 、のように、使う変数の前に半角スペースを置いて書いて定義する必要があります。

testF は、いま定義しようとしている自作関数の名前です。括弧の中で x という変数が与えられた時に { } の処理を行う変数を作成することを宣言しています。

return(x*x); これは x*x の結果を返すことを示しています。もうひとつRと大きな違いとして各行の最後に ;（セミコロン）が入ることです（中括弧 { の後ろ以外）。


次にR上で以下のようにして、この自作関数を呼び出します（初回の呼び出しには数秒掛かるかもしれない）。
sourceCpp("test.cpp")

すると、testF 関数が使用できるようになります。例えば…
testF(1.41421356) # 一夜一夜に人見頃
結果はちゃんと 2 が返ってくるはず。


ちなみにRの上でC++コードを走らせることも可能です。""で挟んで、code=code.testで指定してやるだけです。
sourceCpp(code="
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double testF(double x) {
return(x*x);
} ")

ただし、Rcpp.h以外のファイルを読み込む必要がある場合のやり方が今ひとつ分からず。inlineパッケージのcxxfunctionを使ってできそうなのですが。






次の例では段階を上げて行きます。
rnorm 関数をC++で書いてみましょう。なお、乱数の発生アルゴリズムとしてRの場合はMersenne Twister法が使用されていますが、C++でやろうとすると一工夫必要でした。これについては、詳しい解説とコードも配布されているサイトがありましたのでリンクを貼っておきます。MT.hファイルと内包されるinit_genrand関数、genrand_real3関数はこちらから取得しました。（追記：他にも、匿名さんから頂いたコメントで乱数発生アルゴリズムのまとめサイトがあるとのことです。）

r_normal という名前でC++版rnormを作成します。
下準備としてMt.hと名付けたファイルをRの作業ディレクトリに置きます。下記がRcppから読み込む用のC++コードになります（もはや純粋なC++コードではないので、こう呼びます）。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericVector r_normal( int NN, double mu, double sigma ) {
NumericVector x(NN);
init_genrand((unsigned)time(NULL));

for(int i=0; i<NN; i++) {
x[i] = mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
return(x);
}
新たに加わった include "MT.h" で上記のファイルを詠みこみます

#include <stdio.h>
#include "MT.h"

次に、この行ですが、これがクセモノです。
NumericVector r_normal( int NN, double mu, double sigma ) {

NumericVector は返り値がベクトルになるように定義しますという意味です（Rcpp特有の表現：cf. Rcpp Quick reference）。そして r_normal という名前の自作関数を作ろうとします。

int NN, double mu, double sigma は、ベクトルの要素数 NN（整数）、平均が mu（実数）、標準偏差が sigma（実数）の値を用いることを宣言しています。

NumericVector x(NN);
次にこれですが、xという、要素数が NN のベクトル（初期値は全部0になっている）を用意することを意味しています（Rcpp特有の定義のようだ）。

init_genrand((unsigned)time(NULL));
これは乱数のタネの指定、現在時刻を秒精度で使用します。

for(int i=0; i<NN; i++) {
xのi番目の要素に正規乱数を格納するコード。右辺はごちゃごちゃしていますが、genrand_real3が一様乱数を返す関数です。詳しい説明が参照元にあります。
return(x); で x を返り値で返すようにすれば、r_normal関数が使えるようになります。

rnormと速度比較をしてみます。

NN <- 10^7
system.time(rnorm(n, 0, 1))
# ユーザ システム 経過
# 0.677 0.006 0.679
system.time(r_normal(n, 0, 1))
# ユーザ システム 経過
# 0.363 0.005 0.367

速度は2倍弱といったところ。もう少し期待したかったですが、それでもベクトル化されているRコードよりもさらに高速です。


一方、Rの関数をC++側に渡すことも可能です。Rの関数をC++で計算させることができます。やり方としては、cppFunction関数を用いてrnormのC++版、rnormCppを定義します。
cppFunction( "
NumericVector rnormCpp(int N, double me, double sd) {
return(rnorm(N, me, sd));
}" )

定義の仕方などはだいたい同じですが、rnormが使われていることに注意です。これはC++の関数ではなく、たしかにRのrnorm関数です。
これを通すと、rnormCppという関数が使えるようになります。気になる速度ですが…

system.time(rnormCpp(10^7, 0, 1))
# ユーザ システム 経過
# 0.577 0.019 0.594

Rのrnormより少し速くはなりましたが、C++で1から作るより遅いです。それでも手軽に計算速度を向上させることができるので、何かと試したくなる技です。


ではいよいよ、本題のforループを r_markov 関数として記述してみます（同様にファイル保存します）。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericMatrix r_markov( int repl, int length, double mu, double sigma ) {
NumericMatrix x(repl, length);
int i,j;

init_genrand((unsigned)time(NULL));

for(i=1; i<repl; i++) {
for(j=0; j<length; j++) {
x(i,j) = 0.9*x(i-1,j) + mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
}
return(x);
}

繰り返しを各行に割り当てるため、int repl、各回のベクトル要素数は int length です。
ひとつ前のr_normal関数との違いは、まず毎回の結果を行列に格納している点です。
NumericMatrixが行列型を示しています。
一行目でもxをrepl行、length列の行列として定義します。NumericMatrix x(repl, length);
毎回の繰り返しでは、ひとつ前の結果に0.9を掛けて正規乱数を足していくので、iのforループは1から始まるようにしています（1行目は0のままで放置、Rコードで除去します）。
もう一つ、分かりにくい違いが添字の書き方です。さっきはx[i] = ...と肩のある括弧を使いましたが、行列の場合はx(i,j)と普通の括弧を使用するようです。C++の解説を見るとx[i][j]とありますが、こうするとエラーになりました。x[i,j]もダメです。



それでは冒頭のR版forループと速度比較してみます。
NN <- 10^6
system.time(r_markov(NN+1, 10, 0, 1)[-1,])
# ユーザ システム 経過
# 0.439 0.024 0.460

R版では5.778秒だったので10倍以上早いです。C++の速さを思い知らされました。

Rcpp、少しずつ部分的にC++化できるので、練習にも持って来いだと思いました。たぶん純粋C++に移行することはなくハイブリッドでやっていくことになりそうな気がします。

OSX10.8以降、長らくアップデートしていなかったので、更新することにしました。もういい加減WinBUGSは卒業したくはあるのですが、依然StanやJAGSでは通らないコードも残っておりまだ手放せずにいます。早くStanが離散パラメータを直接扱えるようになってほしい…

クリーンインストールのEl Capitanに入れることを想定して手順を挙げておきます。まず、X11とHomebrew（パッケージ化されていないアプリを容易にインストールするための補助ツール？）を入れ、Homebrewを用いてWine（非Windows OS上でWin専用アプリケーションを実行する環境）をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。以前はMacPortsを使用していましたが、Homebrewの方がはるかに簡単のようです。以前利用できていたユーザの上の階層がEl Capitanでは使いにくくなったというのも大きな理由です。

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

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

cf. usr/local階層に変更を加えるための認証は、OSのマイナーアップグレードの度に必要になるかもしれない。その場合、ターミナルで以下のコードを打ち込む
sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local



**********************
0）実行環境
・OSX10.11 (El Capitan) 搭載のMac

1）App StoreのApple IDを設定しておく


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

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

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


3) 下準備、/usr/local/フォルダを作る、ロックをいったん外して操作をするという動作のようなので、推奨されていない動作だということをお忘れなく。

リカバリモード（⌘+R を押しながら起動）で起動し、ターミナルを立ち上げる


4) ターミナルの $ マークの直後に、下記のコードを打ち込む（コピペでOK、以降も同様）

5) 通常の再起動をする


6) ターミナルに下記を打ち込む（改行されて見えているだろうが、改行無しで打ち込む）

sudo mkdir /usr/local && sudo chflags norestricted /usr/local && sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local


7) 再度、リカバリモードで起動


8) 下記のコードを打ち込む

csrutil enable


9) 今一度、通常の再起動をする


10) homebrewをインストールする

ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

（せいぜい数分でインストールは終了するはず）


11) Wineのインストールを試みる、下記のコードをターミナルに打ち込む

brew install wine

数分ほどでエラーメッセージが出る、メッセージの最後に下記のように書かれている
To continue, you must install Xcode from the App Store, or the CLT by running:
xcode-select --install

先にXcodeをインストールしておけば回避できるのだろうが、このやり方でインストールするほうがむしろ手間が省けるだろう。


12) 上記のメッセージ通り、下記のコードを打ち込む
xcode-select --install

すると、AppStoreからXcodeをインストールしてもよいかと聞かれるのでOKをする。


13) インストールが終わったら、再びWineのインストールを試みる

brew install wine


14) 回線速度にもよるが、良好な環境では30分程度でインストールは終わるだろう。これが終わったら、次のコマンドを打つ。

winecfg

しばらく処理音が聞こえた後、X11からWineの環境設定のようなウインドウが表示される。単に一番下のOKをクリックすればよい。
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよさそうだ。
15) WinBUGSのインストール
パッチなどを当ててある展開済みのWinBUGSフォルダを用意する。Windows7、Windows8へのインストールも現在はこのやり方で行くしかないことを考えれば、Macでも同様にすればいいだろう。適用済みのWinBUGSも公開されている。

下記のコードで不可視フォルダにあるProgram Filesフォルダを開く

open ~/.wine/drive_c/Program\ Files/

ここへWinBUGSフォルダを入れればインストール完了


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

# R2WinBUGSのインストールをお忘れなく
require(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="/usr/local/bin/wine", WINEPATH="/usr/local/bin/winepath")

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

# 今回、opt/localではなくusr/localにパスを通すよう変更する必要が出た。以前のWINE、WINEPATHは/opt/local/bin/になっていたが、ここは/usr/local/bin/に変更していることに注意。


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

RでのGIS操作、いずれまとめようと思いつつ放ったらかしてました。だれにでも有用そうなものから少しずつアップしていく予定です。

基本の関数の備忘録、とくにshapefileの読み込みと書き出しの関数が長くて忘れてしまいがちです。。

（うっかり、同様の検証記事を消してしまったので、ついでにアップデートします）
前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット（cumulative logit）、
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。

しかし、悪い＝0、ふつう＝1、よい＝2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか？たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。

テストデータを用いて、両者の推定と推定値の求め方を比較する。

# まず解析用データの設定
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 3 # 最大3、つまり 0, 1, 2, 3 の値を取りうる
X <- rep(c(1:10), each=10)
Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) # logisticの中身を基にして、二項乱数を発生
Y2 <- factor(Y1, ordered=T) # 累積ロジット用に、ランク化した応答変数を用意
D <- data.frame(X, Y1, Y2)
require(ordinal)
require(VGAM)
M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) # 二項分布のGLM
M2 <- clm(Y2 ~ X, data=D) # 累積ロジット
M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)
# 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定（切片のみランク毎の場合、比例オッズモデルともいうようだ）

累積ロジット（cumulative logit model）を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。
（cf. 累積ロジット：よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ）

例えば、vglmではランダム効果が入れられない。

mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない（手計算すれば良い話ではあるが…）。

最近、ordinalというパッケージを見つけた。1つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。

ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。

require(ordinal)
example(ordinal) # パッケージで用意されている例を出力

# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)

ordinl> summary(fm1)
formula: rating ~ contact + temp
data: wine

link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01

Coefficients:
Estimate Std. Error z value Pr(>|z|) # 回帰係数部分（この例ではカテゴリカルだが）
contactyes 1.5278 0.4766 3.205 0.00135 **
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients: # 各ランクの切片
Estimate Std. Error z value
1|2 -1.3444 0.5171 -2.600
2|3 1.2508 0.4379 2.857
3|4 3.4669 0.5978 5.800
4|5 5.0064 0.7309 6.850


# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)

Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052

広島の生態学会に来ています。発表のネタの準備中に大きなミスを見つけて取り返したりで時間に追い詰められグッタリしてますorz
発表は明後日17日のポスター、国内温帯域の藻場の衰退とサンゴの分布拡大の話です。温暖化影響や生物多様性データベースのGIS化などに興味のある方には興味の持てる内容かと思います。来聴歓迎します。
http://www.esj.ne.jp/meeting/abst/61/PB3-083.html

# 2014年02月20日追記：下記のlmer()は現在はglmer()に相当します

かつて、数値実験をしてみて「できないようだ」と書いたのですが、GLMM側（lmerやglmmML）のAIC計算方法をGLM側（glmやglm.nb）へと合わせるように修正すれば可能であることを確かめました。じつは、両者でAICの計算方法が異なっていたことが分かりました。道理でいつもGLMMのAICが小さ過ぎるわけです…。

★以前アップしていた数値実験を大幅に修正・加筆しておきました（cf.「過分散データ：GLM負の二項分布、GLMMによる解析の比較」）

なお、glmやglm.nb と、lmerやglmmMLのAICの計算方法の違いを検証するところまでは、検索すると複数のブログで見つかります。

ところがさらに検証してみたところ、逆にglmやglm.nbのAICの方をlmerやglmmMLのAIC計算方法へと修正すると、不当にglm、とくにglm.nbのAICが小さくなり、誤ったモデル選択となる確率が大幅に増大することも分かりました（ほぼ逆転してしまう！）。もっとも、ちゃんと数値実験をやっているとはいえ、なぜそうなるのか理論的裏付けがよく分からないので、なにか思いついたら検証したいと思います。

ちなみに、まったくのポアソン分布データのパラメータ推定でも、glm、glm.nb、lmerはいずれも真の値をほぼ推定できていました。過分散問題がけっこう厄介なことを思うと、もはや glm(, family=poisson) を紹介する意義はほとんど無いのかもしれません…。



数値実験では、乱数発生させた、ポアソン分布データと、負の二項分布データ、ポアソン分布に正規乱数ノイズを加えたデータの三者で、GLM（ポアソン分布：glm）とGLM（負の二項分布：glm.nb）、GLMM（ポアソン分布：lmer）が真の値を推定できているかをクロスチェックし、同時にAICが正しい推定をちゃんと反映できているかをもチェックしました。lmerの方は本来のGLMMの使い方とは異なり、各データを異なるid（1データ1id）としているので過分散の対処に用途は限定しています（こういう使い方は、help(lmer)や、Warton & Hui (2011) Ecology 92:3–10、和文では粕谷先生のGLM本「一般化線形モデル」で紹介されています）。

混乱しやすいのでメモ的に整理します。
まず、AIC()を用いた時の、glmやglm.nbと、lmerやglmmMLの計算方法の違い：
glmやglm.nbでは、
AIC = deviance + 2*推定するパラメータ数
lmerやglmmMLでは、
AIC = residual deviance + 2*推定するパラメータ数

また、deviance()は、いずれの場合もresidual devianceを計算する。
一方、logLik()では、glmやglm.nbではdeviance/(-2)、lmerやglmmMLではresidual deviance/(-2)


lmerやglmmMLのAICをglmの方へ合わせるには、下記のようにする（上述の理由により、逆はやめておきましょう）。

# 例えば、こんな2つのモデルがあるとき（id=1~N数）
pois <- glm(Y ~ X, family=poisson)
pois.m <- lmer(Y ~ X + (1|id), family=poisson)

glmのAICはふつうに、AIC(pois)、で求められる。
GLMMのAICは、上述の違いに基づくと、こうすればglmのと同じように求められる：
AIC(pois.m) - deviance(pois) -2*logLik(pois)
# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
#（ランダム効果を抜く以外は同等であるglmモデルを用意してやる必要がある）

# ちなみに、AICがそのように変換できることの検証として、GLMMでidを全部 1 にしてやると（1グループのみ）、グループ間の分散は 0 と推定されます。切片と回帰係数の推定も、単なるGLMのそれと同一の推定結果となります。ただし、分散値は推定するパラメータ数としてカウントされるので、変換後のGLMMのAICはGLMのAICよりもキッカリ2だけ増加した値となります。これは、AICの計算式にある、2*パラメータ数（ここでは、2*1）、の部分に相当します。

しかし、なぜGLMMのパッケージがresidual devianceに基づいたAIC計算にしているのか、その意図が何なのかわかりません。こういう計算方法の修正をしないでくれと言っているようにも見えます。さしあたり、過分散データへの対処という限定的な使用法について数値実験をしてみた限りでは問題は無いようです。

こういう煩わしさを思うと、もはや個人的にはBUGSに全面移行してしまおうかなどと考えたりします…Stanの発展も待ち遠しいです。


# （当時いただいたコメント）
2012年11月15日(木) 23:39:11 | URL | 高橋
興味深く読ませていただきました。
ちょっと気になったのですが、
＞# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
この式の2つのresidual devianceは別物だと思います。
residual deviance = (当該モデルのdeviance) - (飽和モデルのdeviance)
なので、
＞- deviance(pois) -2*logLik(pois)
は、
- ((ポアソンモデルのdeviance) - (ポアソンモデルの飽和モデルのdeviance)) + (ポアソンモデルのdeviance)
となります。
したがって、
＞AIC(pois.m) - deviance(pois) -2*logLik(pois)
は、
(混合モデルのdeviance) - (混合モデルの飽和モデルのdeviance) + 2*パラメータ数 + (ポアソンモデルの飽和モデルのdeviance)
となり、「混合モデルの飽和モデル」が「ポアソンモデルの飽和モデル」と等しい場合（多分そうなのだと思いますが）、混合モデルのdevianceに基づいてAICを求めているということになると思います。
glm.nb()の関数定義を見ると、推定されたThetaのもとで飽和モデルのdevianceを計算しているようです。この値は、ポアソンモデルにおける飽和モデルのdevianceよりずっと大きくなるので、residual devianceを用いてAICを計算すると、負の二項分布の常勝となるのでしょう。
ちなみに飽和モデルの対数尤度は、応答変数をyとすれば、
sum(dpois(y, y, log = TRUE))
で求められます。負の二項分布の場合は、glm.nb()が返すThetaパラメータを使って、
sum(dnbinom(y, mu = y, size = Theta, log = TRUE))