# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。
（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）
### 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のためのオプションなどもデフォルトで付いているようで良いかもしれません）
0 件のコメント:
コメントを投稿