「データ解析のための統計モデリング入門」を読んで引っかかった点

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

データダウンロード

https://kuboweb.github.io/-kubo/ce/IwanamiBook.html

```
#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の範囲からは、先に「基礎からのベイズ統計学」を読むことをおすすめします。

https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6-%E3%83%8F%E3%83%9F%E3%83%AB%E3%83%88%E3%83%8B%E3%82%A2%E3%83%B3%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E5%AE%9F%E8%B7%B5%E7%9A%84%E5%85%A5%E9%96%80-%E8%B1%8A%E7%94%B0-%E7%A7%80%E6%A8%B9/dp/4254122128/ref=sr_1_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=2AC54MWTVYGZD&keywords=%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6&qid=1696005723&sprefix=%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6%2Caps%2C208&sr=8-1

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

# 2章
グラフの設定

https://symfoware.blog.fc2.com/blog-entry-1501.html

```R
plot(lambda,sapply(lambda,logL),type = "l")
```
`type = "l"`は`l`です。イチではありません、エルです。

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

`R`の練習やね

# 3章

https://zenn.dev/rinc5/articles/c3d6d64af4c46f

```R
plot(d$x,d$y,pch = c(21,19)[d$f])
```
が通らないときは

```R
d$f <- factor(d$f)
```
`d$f`を`factor`型に変更する。文字を数値に変換してくれる。

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

https://qiita.com/stkdev/items/6aba2c1db2fa056170ae


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

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

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

```math
 \log L^* = \max  \log L
```
と考えるみたいな感じか？

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

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

https://www.amazon.co.jp/%E6%83%85%E5%A0%B1%E9%87%8F%E7%B5%B1%E8%A8%88%E5%AD%A6-%E6%83%85%E5%A0%B1%E7%A7%91%E5%AD%A6%E8%AC%9B%E5%BA%A7-%E3%83%BB5%E3%83%BB4/dp/4320021711/ref=sr_1_1?crid=33VIGE7TCRT48&keywords=%E6%83%85%E5%A0%B1%E9%87%8F%E7%B5%B1%E8%A8%88%E5%AD%A6&qid=1695984969&sprefix=jyouhouryoutoukei%2Caps%2C243&sr=8-1

https://www.amazon.co.jp/%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96-%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA%E3%83%BB%E4%BA%88%E6%B8%AC%E3%81%A8%E7%99%BA%E8%A6%8B%E3%81%AE%E7%A7%91%E5%AD%A6-%E5%B0%8F%E8%A5%BF-%E8%B2%9E%E5%89%87/dp/4254127820/ref=sr_1_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=1S1C1WY7MN4TU&keywords=%E6%83%85%E5%A0%B1%E9%87%8F%E7%B5%B1%E8%A6%8F%E6%BA%96&qid=1695984986&sprefix=%E6%83%85%E5%A0%B1%E9%87%8F%E7%B5%B1%E8%A6%8F%E6%BA%96%2Caps%2C185&sr=8-1


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

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

```math
AIC = -2 \times (\log L^* -k)
```
これの意味を理解しようとする。先程の関係式から

```math
AIC = -2 \times E[\log L]
```
つまり、

```math
AIC = -2 \times (得られる対数尤度の期待値)
```
よって、AICは尤度の期待値のことであります。$-2$がかかっているので、
`AIC`高い->尤度**期待値**低い->使えないモデル
`AIC`低い->尤度**期待値**高い->使えるモデル！

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


他の解釈では

```math
AIC = -2 \times (最大尤度 - パラメータ増加による尤度の減少)
```
`AIC`が変わらない -> パラメータを増やした意味がない
`AIC`が増加       ->　パラメータで本質的尤度が減った -> パラメータを増やしたのは逆効果
`AIC`が減少       ->　パラメータで本質的尤度が増えた -> パラメータを増やしたのは意味がある

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


**霊夢知ってるか、`AIC`が<font color="Red">減少する</font>ことはとにかく<font color="Red">良い</font>ことなんだぜ!**

## AICの名前について
赤池情報量規準：An Information Criterion, Akaike's Information Criterion

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

# 5章

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

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

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

# 6章

```R
kubobook_2012 // binomial // data4a.csv

d<-read_csv("data4a.csv")
summary(d)
```
をロード



```R
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)
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/7f7544f9-6d1a-969e-6a05-d84fa848c64e.png)


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

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

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

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

```math
\frac{タバコを吸う人が死ぬ確率}{タバコを吸わない人が死ぬ確率} = 7
```
と考えるか


```math
\frac{死んだ人がタバコを吸っていた確率}{死んだ人がタバコを吸っていない確率} = 7
```
どちらが正しいでしょうか。


https://bellcurve.jp/statistics/course/26726.html#:~:text=%E3%82%AA%E3%83%83%E3%82%BA%E6%AF%94%E3%81%AE%E4%BD%BF%E3%81%84%E6%96%B9,%2D4%E7%AB%A0%E3%82%92%E5%8F%82%E7%85%A7%EF%BC%89%E3%80%82

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

更に、

```math
\frac{死んだ人がタバコを吸っていた確率}{死んだ人がタバコを吸っていない確率} 
/
\frac{生きている人がタバコを吸っていた確率}{生きている人がタバコを吸っていない確率} 
```
をオッズ比と呼びます。解釈をすると

```math
(煙草によって死んだ確率)/(煙草によって生きている確率) = オッズ比
```
オッズ比は原因の比較という感じがします。

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


例

```math
タバコ \rightarrow 死　
```
という関係について、

-------------

```math
前向き研究：タバコを吸う・吸わない集団を集めた
```
その後、ある年代で生きているか・死んでいるか確認した。

```math
リスク比 = \frac{タバコを吸う人が死ぬ確率}{タバコを吸わない人が死ぬ確率} 
```
が求められる

--------------

```math
後ろ向き研究：ある年代で死んだ人・生きている集団を集めた
```
タバコを吸っているか確認した。

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

## AIC

```R
library(MASS)
stepAIC(fit.xf)

fit.null <- glm(cbind(y, N -y )~1,data = d,family = binomial)
-2*(logLik(fit.null)-1)   # 一定のAIC
```
で求められる。

## 交互作用

```R
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を実際に試す
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/6f60d08b-da94-3e43-cf3b-f8b79f79635f.png)

```R
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)

```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/83aa2e9b-28dc-1b97-2281-ca70d8418a10.png)

```R
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のアイデア

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

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

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


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

```math
r_iを最尤推定で最適化する \rightarrow　フルモデルで最適化する
```
よって、$r_i$を推定するのは得策ではない。

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


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

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

```math
\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$を調整するだけ。

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


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

## glmmMLのインストール

```
utils:::menuInstallPkgs()
```
を実行して、インストール

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/bae53876-2642-8c51-6dac-07dcc31bbfcd.png)

```
install.packages("glmmML")
```
でいいって書いてあった。

# 8章：メトロポリスヘイスティングMH法
「基礎からのベイズ統計学」を読んでいたので軽く流した。

https://www.amazon.co.jp/%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6-%E3%83%8F%E3%83%9F%E3%83%AB%E3%83%88%E3%83%8B%E3%82%A2%E3%83%B3%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E5%AE%9F%E8%B7%B5%E7%9A%84%E5%85%A5%E9%96%80-%E8%B1%8A%E7%94%B0-%E7%A7%80%E6%A8%B9/dp/4254122128/ref=sr_1_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=2AC54MWTVYGZD&keywords=%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6&qid=1696005723&sprefix=%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E3%81%AE%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E5%AD%A6%2Caps%2C208&sr=8-1

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

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

```math
GLMの統計モデル＋観測値y_i + 事前分布 \beta　\rightarrow \betaの事後分布
```
を行うことがこの章の目的である。
これまで最尤推定で求めてきた$\beta$であったが、この章から求め方を変えるのである。

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



## WinBUGSの導入

https://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/

から、WinBUGSをダウンロード

```
C://Program Files//WinBUGS14
```
になるよう配置する。



"R2WBWrapper.R"は

https://kuboweb.github.io/-kubo/ce/IwanamiBook.html

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

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

## 実行コード

```WinBUGS: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)
}
```

```R
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")])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/db619db5-bd57-ea1b-7dbc-2cb2d8cf1b40.png)

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/6901c62b-0f96-a392-5415-de3c6e748ef3.png)

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

## ギブスサンプリング

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

# 10章：階層ベイズモデル
7章でGLMに$r_i$を追加したGLMMを考えた。

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

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

```math
r_i \sim 適当な分布関数
```
**$r_i$は未知で何らかの分布に従うと仮定した。**

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

```math

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

# 11章：空間構造のある階層ベイズモデル
個体差だけではなく、隣の農場・隣町の農場のような**グループ分け**が可能という状況は往々にしてある。
これをGLMMモデルに導入する


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

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

```math
r_i \sim 適当な分布関数
```

```math
m_j \sim 適当な分布関数
```

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

![beizu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/41e3e67c-80e0-80fb-39d8-b480d85f1ef9.png)

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

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