通称、緑本を読んで突っかかった点と解決をまとめます。

データダウンロード

#bookchapter: folder
ch01: introduction
ch02: distribution
ch03: poisson
ch04: aic
ch05: lrtest
ch06: binomial, gamma
ch07: glmm
ch08: mcmc
ch09: gibbs
ch10: hbm
ch11: spatial

ファイル名にチャプター番号がないので、こちら（Readmeの内容）を参照。
data.RDataがdistributionとaic両方にあるので引っかからないように。
データロードで突っかからないようにね。

（先に）学習後の感想

この本はできれば、Rで実際に試しながらやるべきです。
なぜRであるかといえば、計算に必要な最適化処理がデフォルトで入っているからです。

この本は統計モデルの考え方を教えてくれる本です。その統計モデルになるように同定する手法は載っていません。というより、解法からは統計ではなく数学の範囲になるためです。
我々の学習範囲は既にMCMCなど手計算では不可能なコンピュータ依存の問題に入ってきています。
これを機にRを触ってみてはいかがでしょうか。

内容としては非常にわかりやすく、サラッと読めてしまいます。数式的で具体的な計算はRにまかしていることもわかりやすさに繋がっています。結果として、RかPythonで試さなければあまり頭に入らない気がします。後半はWinBUGSにお任せなので、まぁブラックボックス化してますが、目的が頭に入ればいいかもです。

後半の8章、MCMCの範囲からは、先に「基礎からのベイズ統計学」を読むことをおすすめします。

学習時間：2週間、40hourぐらい？

2章

グラフの設定

plot(lambda,sapply(lambda,logL),type = "l")

type = "l"はlです。イチではありません、エルです。

seq(start,end,step)はrangeみたいに使える。

Rの練習やね

3章

plot(d$x,d$y,pch = c(21,19)[d$f])

が通らないときは

d$f <- factor(d$f)

d$fをfactor型に変更する。文字を数値に変換してくれる。

Rのコロン.は意味を持たず、文字列として扱われる

4章：AIC

AICの解説の章。ただし、理論的な解説はしていない。
個人的にまとめると

データをサンプルして最大尤度となるようにパラメータを推定する。パラメータ数を$k$とする。
結果得られる、対数尤度$\log L$
この対数尤度は、サンプルを変えれば別の値になる
たくさんのサンプル群を用意して、対数尤度の分布を生成
これを対数尤度$\log L$の確率分布と考える。
何らかの仮定をすれば

E[ \log L] = \log L^* - k

が成立するらしい。$\log L^*$が期待値でないのはなぜだろうか・・・全サンプル群での最大尤度としているからだと予想する。つまり、

 \log L^* = \max  \log L

と考えるみたいな感じか？

パラメータ数$k$が増えるほど、最大対数尤度$\log L$の期待値は減少していく。
パラメータ増えると、当てはまりが良くなり、尤度が増えそうと考えるが、過学習となって、状況適応性が下がり、サンプルセットを変えると急に尤度が下がるみたいなことかな？

正確な数理統計学的導出は「情報量統計学」または「情報量規準 (シリーズ・予測と発見の科学) 」にあるらしい。

知りたいが、かなり難解で仮定も厳しかったり、それを知ることでなにかできるわけではなく、裏が取れるというものなので先送りにする。誰かやさしく教えてほしい

とりあえず、AICを以下のように定義する。

AIC = -2 \times (\log L^* -k)

これの意味を理解しようとする。先程の関係式から

AIC = -2 \times E[\log L]

つまり、

AIC = -2 \times (得られる対数尤度の期待値)

よって、AICは尤度の期待値のことであります。$-2$がかかっているので、
AIC高い->尤度期待値低い->使えないモデル
AIC低い->尤度期待値高い->使えるモデル！

AICは期待値であることを言うために、パラメータ数を引いているんですね。

他の解釈では

AIC = -2 \times (最大尤度 - パラメータ増加による尤度の減少)

AICが変わらない -> パラメータを増やした意味がない
AICが増加 ->　パラメータで本質的尤度が減った -> パラメータを増やしたのは逆効果
AICが減少 ->　パラメータで本質的尤度が増えた -> パラメータを増やしたのは意味がある

ということが言えるぐらいの理解ができます。

霊夢知ってるか、AICが減少することはとにかく良いことなんだぜ!

AICの名前について

赤池情報量規準：An Information Criterion, Akaike's Information Criterion

赤池って名字あまりいないですよね。
個人名があるとなんか一般的なものでないと感じてしまいますが、多分赤池さんはAnの方を意図したと思います。赤池さんの論文が広がってAkaikeという解釈になったのかもですね。

5章

第二種の過誤はネイマン検定では起きない。帰無仮説$H0$を棄却できなければ、$H0$が正しいとも・誤っているとも言えないためである。

第一種の過誤は$H0$を棄却したうえで、それが誤りである状態である。P値の薄い確率を引いた状態である。

パラメトリックブートストラップ法は得られているサンプルの代表量と帰無仮説統計モデルを元にサンプルを新たに複数生成して、尤度分布を用意する。対立仮説統計モデルの尤度から確率P値が得られる。

6章

kubobook_2012 // binomial // data4a.csv

d<-read_csv("data4a.csv")
summary(d)

をロード

d<-read_csv("data4a.csv")
summary(d)

fit.xf <- glm(cbind(y, N -y )~x+f,data = d,family = binomial)

d$f <- factor(d$f)
#plot(d$x,d$y/8,xlim = c(8,13),ylim = c(0,1),pch = c(21,19)[d$f])
plot(d$x,d$y/8,xlim = c(8,13),pch = c(21,19)[d$f],xlab="x", ylab="y")
legend("topleft",legend = c("C","T"),pch = c(21,19))

logistic <- function(z){
  1/(1+exp(-z))
}
z <- seq(8,13,0.1)
par(new=T)
plot(z,logistic(-19.536+1.952*z),type="l",xlim = c(8,13),ylim = c(0,1),col = 2,ann = F)
par(new=T)
plot(z,logistic(-19.536+1.952*z+2.022),type="l",xlim = c(8,13),ylim = c(0,1),col = 3,ann = F)

library(MASS)
stepAIC(fit.xf)

当てはまりが良さそうな曲線がかけてます。

追記（Nで割っていますが、確率変数同士の割り算になるので良くないようです）

（脱線）リスクとオッズ比

タバコによる、肺がん死亡リスクは7倍と言われたら

\frac{タバコを吸う人が死ぬ確率}{タバコを吸わない人が死ぬ確率} = 7

と考えるか

\frac{死んだ人がタバコを吸っていた確率}{死んだ人がタバコを吸っていない確率} = 7

どちらが正しいでしょうか。

前者はリスク比、後者はオッズと呼ばれます。

更に、

\frac{死んだ人がタバコを吸っていた確率}{死んだ人がタバコを吸っていない確率} 
/
\frac{生きている人がタバコを吸っていた確率}{生きている人がタバコを吸っていない確率}

をオッズ比と呼びます。解釈をすると

(煙草によって死んだ確率)/(煙草によって生きている確率) = オッズ比

オッズ比は原因の比較という感じがします。

リンクのサイトがわかりやすいですが、後ろ向き調査（結果からサンプルを集める）では、リスク比を求めることはできません。リスク比を求めるためには前向き研究である必要があります。

例

タバコ \rightarrow 死

という関係について、

前向き研究：タバコを吸う・吸わない集団を集めた

その後、ある年代で生きているか・死んでいるか確認した。

リスク比 = \frac{タバコを吸う人が死ぬ確率}{タバコを吸わない人が死ぬ確率}

が求められる

後ろ向き研究：ある年代で死んだ人・生きている集団を集めた

タバコを吸っているか確認した。

オッズ比 = 
\frac{死んだ人がタバコを吸っていた確率}{死んだ人がタバコを吸っていない確率} 
/
\frac{生きている人がタバコを吸っていた確率}{生きている人がタバコを吸っていない確率}

を求めた。

AIC

library(MASS)
stepAIC(fit.xf)

fit.null <- glm(cbind(y, N -y )~1,data = d,family = binomial)
-2*(logLik(fit.null)-1)   # 一定のAIC

で求められる。

交互作用

d<-read_csv("data4a.csv")
summary(d)

fit.xf <- glm(cbind(y, N -y )~x*f,data = d,family = binomial)
fit.xf

7章：GLMM

7.1を実際に試す

data <- read_csv("data.csv")

plot(jitter(data$x,0.2),data$y)

fit.x <- glm(cbind(y, N -y )~x,data = data,family = binomial)

plot(jitter(data$x,0.4),data$y,xlab="x", ylab="y",ann=F)

data <- read_csv("data.csv")

plot(jitter(data$x,0.2),data$y)

fit.x <- glm(cbind(y, N -y )~x,data = data,family = binomial)

xi <- 4
yi <- table(data$y[data$x == xi])

plot(yi,type="p", xlab="x", xlim = c(0,8),ylim = c(0,6),ylab="観測された個体数")
par(new=T)

yj <- seq(0,8,1)
zj <- -2.1487 + 0.5104 * xi

calc.prob <- function(x) { 
  p <- exp(x) / (1 + exp(-1*x))
  return (p)
}
pj <-  calc.prob(zj)
Nj <- dbinom(yj,8,pj)*sum(yi)
plot(yj,Nj,xlim = c(0,8),ylim = c(0,6),type="b",ann=F)

GLMMのアイデア

リンク関数(\theta_i) = \beta_1 + \beta_2 x_i  + r_i

f(y_i) = 分布関数(\theta_i)

r_i \sim 適当な分布関数

GLMに$r_i$を追加した。$r_i$は未知で何らかの分布に従うと仮定する。
$r_i$を追加することで、観測できていない要素の影響を吸収する。

$r$の分布関数と、$y$の分布関数が混合するような形になることから、GLMMにはMixというMの文字が入ってくる。

r_iを最尤推定で最適化する \rightarrow　フルモデルで最適化する

よって、$r_i$を推定するのは得策ではない。

独立した反復（$y_i$をグループ化できない、個体差・場所差の効果が除去できる）であれば、GLMで推定する。
個体・場所などでサンプルをグループ化できれば、階層ベイズモデルを使う

p.156の式

前提
$\beta_1,\beta_2,s$は確率変数ではない。$y_i,r_i$のみが確率変数である。

\beta_1,\beta_2,s にある数値を与えた状態で以下を計算する。

\int p(y_i|r_i) p(r_i) dr_i = \int \frac{p(y_i \bigcap  r_i)}{p(r_i)} dr_i  p(r_i) 
= \int p(y_i \bigcap r_i) dr_i  = p(y_i)

観測された$y_i$に対応する$p(y_i)$が尤度である。

あとは尤度関数が最大化するようにパラメータ$\beta_1,\beta_2,s$を調整するだけ。

\max_{\beta_1,\beta_2,s}  \prod p(y_1) \cdots p(y_n)

具体的な計算はRがやってくれる。

glmmMLのインストール

utils:::menuInstallPkgs()

を実行して、インストール

install.packages("glmmML")

でいいって書いてあった。

8章：メトロポリスヘイスティングMH法

「基礎からのベイズ統計学」を読んでいたので軽く流した。

言っていることはMH法
MH法は初期値の影響が強いので試す人は注意しよう

9章：GLMのパラメータをMCMC法で推定する

これまでパラメタ$\beta_1,\beta_2$は定数を仮定していた。
ベイズ学派では定数を確率変数として考える。
シュレディンガーの猫と同じで、わかってない定数$\beta$を$1$または$0$どちらも同じ確率でとりえるみたいに考え、定数に確率分布をもたせるのである。

GLMの統計モデル＋観測値y_i + 事前分布 \beta　\rightarrow \betaの事後分布

を行うことがこの章の目的である。
これまで最尤推定で求めてきた$\beta$であったが、この章から求め方を変えるのである。

MCMC法は事後分布を計算するアルゴリズムの総称である。

WinBUGSの導入

から、WinBUGSをダウンロード

C://Program Files//WinBUGS14

になるよう配置する。

"R2WBWrapper.R"は

からダウンロードして、RのWorkingディレクトリにコピー

WinBUGSはMCMCの厄介な部分を代わりにやってくれる。

実行コード

model.bug.txt

model{
for(i in 1:N){
	Y[i] ~ dpois(lambda[i])
	log(lambda[i]) <- beta1 + beta2 * (X[i] - Mean.X)
}
beta1 ~ dnorm(0,1.0E-4)
beta2 ~ dnorm(0,1.0E-4)
}

source("R2WBwrapper.R")
load("d.RData")
head(d)

clear.data.param()

set.data("N",nrow(d))
set.data("Y",d$y)
set.data("X",d$x)
set.data("Mean.X",mean(d$x))

set.param("beta1",0)
set.param("beta2",0)

post.bugs <- call.bugs(
  file = "model.bug.txt",
  n.iter = 1600,n.burnin = 100, n.thin = 3
)

post.list <- to.list(post.bugs)
post.mcmc <- to.mcmc(post.bugs)
s <- colnames(post.mcmc) %in% c("beta1","beta2")
plot(post.list[,s,])

plot(as.matrix(post.mcmc)[,c("beta1","beta2")])

パラメータ$\beta_1,\beta_2$の事後分布を求めることができています。

パラメータ$\beta_1,\beta_2$に相関がないこともわかります。

ギブスサンプリング

複数パラメータ版のMH法？

10章：階層ベイズモデル

7章でGLMに$r_i$を追加したGLMMを考えた。

リンク関数(\theta_i) = \beta_1 + \beta_2 x_i  + r_i

f(y_i) = 分布関数(\theta_i)

r_i \sim 適当な分布関数

$r_i$は未知で何らかの分布に従うと仮定した。

この個体差$r_i$を何らかに分布があり


r_i \sim N(0,s^2)

みたいな、新たな代表パラメタ$s$で$r_i$を記述しようとしてみよう！
$r_i \rightarrow s$のような確率変数の階層構造を階層ベイズモデルという。

11章：空間構造のある階層ベイズモデル

個体差だけではなく、隣の農場・隣町の農場のようなグループ分けが可能という状況は往々にしてある。
これをGLMMモデルに導入する

リンク関数(\theta_i) = \beta_1 + \beta_2 x_i  + r_i + m_j

f(y_i) = 分布関数(\theta_i)

r_i \sim 適当な分布関数

m_j \sim 適当な分布関数

$m_j$には空間的構造がある。

$m_1$と$m_2,m_4$は似通った値を取るはずである。近所であることで環境が近い事になり、似通った要素があるはずである。これが空間的構造である。

m_j \sim 適当な分布関数

この分布関数の取り方はわからないが、合理性のあるモデルをとれればよい。