kubo2017sep16a p.1 ( 1 ) : : :55 kubo ( ( 1 ) / 10

Size: px

Start display at page:

Download "kubo2017sep16a p.1 ( 1 ) : : :55 kubo ( ( 1 ) / 10"

つねたけかみいしづ
4 years ago
Views:

1 kubo2017sep16a p.1 ( 1 ) kubo@ees.hokudai.ac.jp : : :55 kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 :! 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル : GLM GLIM General Linear Model kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 () 1 2 :? 3 GLM: GLM R GLM 4 GLM GLM 1. kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

2 kubo2017sep16a p.2 : :?! ( ) method section p < 0.05 OK kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106? (): kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : 0, 1, 2 (y {0, 1, 2, 3, } ) y x y? kubo ( ( 1 ) / 106 x y ? y 0? x NO! kubo ( ( 1 ) / 106

3 kubo2017sep16a p.3? y x YES! kubo ( ( 1 ) / 106 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル kubo ( ( 1 ) / 106 : : (GLM) 3. MCMC 2. : : 1, 2, 3 4. GLM kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : : () : : kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

4 kubo2017sep16a p.4 : () : R i 50 i {1, 2, 3,, 50} y i {y i }! {y i } = {y 1, y 2,, y 50 } {y i } R > data [1] [26] ! kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : R : table() > table(data) ( ) > hist(data, breaks = seq(-0.5, 9.5, 1))! kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : > mean(data) [1] 3.56 > abline(v = mean(data), col = "red") kubo ( ( 1 ) / 106 : : > var(data) [1] (SD = variance) > sd(data) [1] > sqrt(var(data)) [1] : : kubo ( ( 1 ) / 106

5 kubo2017sep16a p.5 :? :?? kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 :? :? > data.table <- table(factor(data, levels = 0:10)) > cbind(y = data.table, prob = data.table / 50) y prob kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 :? :?? () y p(y λ) = λy exp( λ) y! y! y 4! exp( λ) = e λ (e = ) kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

6 kubo2017sep16a p.6 :? :?? R > y <- 0:9 # () > prob <- dpois(y, lambda = 3.56) # > plot(y, prob, type = "b", lty = 2) > # cbind > cbind(y, prob) (λ) 3.56 Poisson distribution y prob kubo ( ( 1 ) / 106 > hist(data, seq(-0.5, 8.5, 0.5)) # > lines(y, prob, type = "b", lty = 2) # kubo ( ( 1 ) / 106 :? :? λ? λ λ λ 0 : λ = = > # cbind > cbind(y, prob) y prob y {0, 1, 2,, } y 1 p(y λ) = 1 y=0 kubo ( ( 1 ) / 106 : y i {0, 1, 2, } y i : kubo ( ( 1 ) / 106 :? : λ p(y λ) = λy exp( λ) y! λ kubo ( ( 1 ) / 106!? kubo ( ( 1 ) / 106

7 kubo2017sep16a p.7 : (likelihood)? : L(λ) λ λ λ 3 {y 1, y 2, y 3 } = {2, 2, 4} = (the likelihood definition for the example): L(λ) = (y 1 2 ) (y 2 2 ) (y 50 3 ) = p(y 1 λ) p(y 2 λ) p(y 3 λ) p(y 50 λ) = p(y i λ) = λ yi exp( λ), y i i i! kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : : λ () (!) (log likelihood function) log L(λ) = i ( ) y i y i log λ λ log k k log L(λ) L(λ) λ kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : : ˆλ log L(λ) = ( i yi log λ λ y i k log k) log likelihood * ˆλ = (ML estimator): i y i/50 d log L dλ (ML estimate): ˆλ = 3.56 = 0! λ λ ˆλ ˆλ λ 50 kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

8 kubo2017sep16a p.8 : : : λ = y ˆλ = y y kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : : λ = y ˆλ = y : ( : y {0, 1, 2, 3, } y : y {0, 1, 2,, N} N y : < y < kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 : : GLMM 線形モデルの発展階層ベイズモデルもっと自由な統計モデリングを! (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル? (GLM) (Poisson regression) (logistic regression) (linear regression) kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

9 kubo2017sep16a p.9 GLM: GLM: i 3. GLM:? : {y i } : {x i } {f i } f i C: T: y i x i (f i = C): 50 sample (i {1, 2, 50}) (f i = T): 50 sample (i {51, 52, 100}) kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: GLM: data frame d : data: data3a.csv CSV (comma separated value) format file R : > d <- read.csv("data3a.csv") d data frame () data frame d > d y x f C C C T T > d$x [1] [9] [97] > d$y [1] [17] [97] kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: data frame d : GLM: R f > d$f [1] C C C C C C C C C C C C C C C C C C C C C C C C C [26] C C C C C C C C C C C C C C C C C C C C C C C C C [51] T T T T T T T T T T T T T T T T T T T T T T T T T [76] T T T T T T T T T T T T T T T T T T T T T T T T T Levels: C T : C T 2 > class(d) # d data.frame [1] "data.frame" > class(d$y) # y integer [1] "integer" > class(d$x) # x numeric [1] "numeric" > class(d$f) # f factor [1] "factor" kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

10 kubo2017sep16a p.10 GLM: data frame summary() GLM:! Generate Data Plots! Always! > summary(d) y x f Min. : 2.00 Min. : C:50 1st Qu.: st Qu.: T:50 Median : 8.00 Median : Mean : 7.83 Mean : rd Qu.: rd Qu.: Max. :15.00 Max. : kubo ( ( 1 ) / 106 > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) kubo ( ( 1 ) / 106 GLM: GLM: GLM f (box-whisker plot) > plot(d$f, d$y) # note that d$f is factor type! GLM kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: GLM GLM: GLM GLM (GLM)??? : : e.g., β 1 + β 2 x i : () + () x i : kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

11 kubo2017sep16a p.11 GLM: GLM GLM: GLM (?) GLM : (response variable) : (explanatory variable) (linear predictor): () = (, intercept) + ( 1) ( 1) + ( 2) ( 2) : : e.g., β 1 + β 2 x i : ( 3) ( 3) + kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: GLM GLM: GLM GLM logistic R (GLM) GLM () rbinom() glm(family = binomial) : : e.g., β 1 + β 2 x i : logit yi x i rbinom() glm(family = binomial) rpois() glm(family = poisson) rnbinom() glm.nb() in library(mass) () rgamma() glm(family = gamma) rnorm() glm(family = gaussian) glm() GLM kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: GLM Go buck to the seed example y i λ i p(y i λ i ) = λyi i exp( λ i ) y i! i λ i? λ i = exp(β 1 + β 2 x i ) β 1 β 2 () x i i f i kubo ( ( 1 ) / 106 GLM:? i λi GLM λ i = exp(β 1 + β 2 x i ) {β 1, β 2} = { 2, 0.8} {β 1, β 2} = { 1, 0.4} i x i kubo ( ( 1 ) / 106

12 kubo2017sep16a p.12 GLM: GLM GLM: GLM GLM () i λ i λ i = exp(β 1 + β 2 x i ) log(λ i ) = β 1 + β 2 x i log() = : : β 1 + β 2 x i : log kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: R GLM GLM: R GLM glm() R GLM > d y x f C C C T T! > fit <- glm(y ~ x, data = d, family = poisson) kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: R GLM GLM: R GLM glm() glm() > fit <- glm(y ~ x, data = d, family = poisson) all: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x ( z):? link : z (y)? family:? Degrees of Freedom: 99 Total (i.e. Null); Null Deviance: 89.5 Residual Deviance: 85 AIC: Residual kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

13 kubo2017sep16a p.13 GLM: R GLM GLM: R GLM glm() > summary(fit) Call: glm(formula = y ~ x, family = poisson, data = d) Deviance Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) x () kubo ( ( 1 ) / 106 () β 2 (Estimate , SE ) (Estimate 1.29, SE 0.364) β p p ˆβ p 0.5 ˆβ (: ) kubo ( ( 1 ) / 106 GLM: R GLM GLM: R GLM (?) β 2 (Estimate , SE ) (Estimate 1.29, SE 0.364) β %? kubo ( ( 1 ) / 106 > fit <- glm(y ~ x, data = d, family = poisson)... Coefficients: (Intercept) x > plot(d$x, d$y, pch = c(21, 19)[d$f]) # data > xp <- seq(min(d$x), max(d$x), length = 100) > lines(xp, exp( * xp)) kubo ( ( 1 ) / 106 GLM: GLM: f i GLM + y i λ i i λ i β 3 f i p(y i λ i ) = λyi i exp( λ i ) y i! λ i = exp(β 1 + β 2 x i + β 3 d i ) d i = { 0 (f i = C ) 1 (f i = T ) kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

14 kubo2017sep16a p.14 GLM: GLM: glm(y x + f,...) > summary(glm(y ~ x + f, data = d, family = poisson))...()... Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) x ft () x + f > plot(d$x, d$y, pch = c(21, 19)[d$f]) # data > xp <- seq(min(d$x), max(d$x), length = 100) > lines(xp, exp( * xp), col = "blue", lwd = 3) # C > lines(xp, exp( * xp ), col = "red", lwd = 3) # T kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM: GLM: λi f i = C: λ i = exp( x i ) f i = T: λ i = exp( x i 0.032) = exp( x i ) exp( 0.032) x i exp( 0.032)! kubo ( ( 1 ) / 106 λi (A) (B) λ = exp(β 1 + β 2 x + ) λ = β 1 + β 2 x + x i x i kubo ( ( 1 ) / 106 GLM: GLM: GLM: y x log y x 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

15 kubo2017sep16a p.15 GLM GLM : 4. GLM N i = 8 : 8 : i y i = 3 : kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM example data: dead or alive of seeds y i (no individual difference!) kubo ( ( 1 ) / 106 GLM q q i N i y i ) p(y i q) = In this example... ( Ni y i q y i (1 q) N i y i, no individual difference q kubo ( ( 1 ) / 106 GLM : 20 {y i } q 20 q L(q {y i }) = 20 i=1 p(y i q) kubo ( ( 1 ) / 106 GLM L(q ) ˆq 20 log L(q ) = log i=1 ( Ni 20 + {y i log(q) + (N i y i ) log(1 q)} i=1 q kubo ( ( 1 ) / 106 y i )

16 kubo2017sep16a p.16 GLM (MLE) do you remember? GLM 8 y i L(q ) q log L(q ) q 0 ˆq log L(q )/ q = 0 q log likelihood ˆq = = = * ˆq = 0.46 ( ) 8 y 0.46 y y y i kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM GLM GLM GLM GLM logistic GLM : : e.g., β 1 + β 2 x i yi : logit x i kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM GLM GLM GLM N y 8 y! f i C: T: i N i = 8 y i = 3 (alive) (dead) x i data4a.csv CSV (comma separated value) format file R : > d <- read.csv("data4a.csv") or > d <- read.csv( + " d data frame () kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

17 kubo2017sep16a p.17 GLM GLM GLM GLM data frame d > summary(d) N y x f Min. :8 Min. :0.00 Min. : C:50 1st Qu.:8 1st Qu.:3.00 1st Qu.: T:50 Median :8 Median :6.00 Median : Mean :8 Mean :5.08 Mean : rd Qu.:8 3rd Qu.:8.00 3rd Qu.: Max. :8 Max. :8.00 Max. : > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) yi x i? kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM GLM : N y p(y N, q) = ( ) N q y (1 q) N y y ( N ) y N y logit p(y i 8, q) q = 0.1 q = 0.8 q = y i kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM GLM (z i: e.g. z i = β 1 + β 2x i) 1 q i = logistic(z i ) = 1 + exp( z i ) > logistic <- function(z) 1 / (1 + exp(-z)) # > z <- seq(-6, 6, 0.1) > plot(z, logistic(z), type = "l") q q = 1 1+exp( z) z kubo ( ( 1 ) / 106 {β 1, β 2 } = {0, 2}(A) β 2 = 2 β 1 (B) β 1 = 0 β 2 q (A) β 2 = 2 β 1 = 2 β 1 = x β 1 = (B) β 1 = 0 β 2 = 4 β 2 = x β 2 = 1 {β 1, β 2 } x q 0 q 1 kubo ( ( 1 ) / 106

18 kubo2017sep16a p.18 GLM logit link function GLM R β 1 β 2 logistic 1 q = 1 + exp( (β 1 + β 2 x)) = logistic(β 1 + β 2 x) logit q logit(q) = log 1 q = β 1 + β 2 x logit logistic logistic logit logit is the inverse function of logistic function, vice versa y (A) f i =C x (B) > glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)... Coefficients: (Intercept) x ft x kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106 GLM GLM : : yi (A) f i =C x i (B) f i =T x i (GLM) 3. MCMC 4. GLM kubo ( ( 1 ) / 106 kubo ( ( 1 ) / 106

19 kubo2017sep16b p.1 ( 2 ) GLM kubo@ees.hokudai.ac.jp : : :58 kubo ( ( 2 ) / 90 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル kubo ( ( 2 ) / 90 : (GLM) 3. MCMC () 1 MCMC MCMC: 2 MCMC GLM kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 MCMC MCMC : () 1. MCMC : () y i kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

20 kubo2017sep16b p.2 MCMC q () MCMC (MLE) () i N i y i ) p(y i q) = ( Ni y i q y i (1 q) N i y i, q L(q ) q log L(q ) q 0 ˆq log L(q )/ q = 0 q log likelihood ˆq = = = * kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 MCMC 8 y i () MCMC ˆq = 0.46 ( ) 8 y 0.46 y y y i kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 MCMC MCMC : p(q Y) = () p(y q) p(q) p(y) () p(y q): q Y () p(q): q () p(q Y): Y q () p(y): Y (q!) 1 p(y q) 2 p(q) 3 p(q Y ) 1 p(y q) 2 p(q) 3 p(q Y ) p(y q) 4 p(q)? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

kubo2017sep16b p.3 MCMC MCMC log likelihood -50-45 -40 log L(q) * 0.3 0.4 0.5 0.6 q L(q) p(y q) L(q) q? p(q Y) Y p(y q) L(q) p(y q) q p(q) OK? 0.0 0.2 0.4 0.6 0.8 1.0? : q q?

gl/ufq2) ( 2 ) 2017 09 16 16 / 90 MCMC MCMC: MCMC MCMC: : MCMC () Markov chain Monte Carlo (MCMC) Metropolis Method (Metropolis method) :?? MCMC Metropolis Method kubo (http://goo.

21 kubo2017sep16b p.3 MCMC MCMC log likelihood log L(q) * q L(q) p(y q) L(q) q? p(q Y) Y p(y q) L(q) p(y q) q p(q) OK? ? : q q? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 MCMC MCMC MCMC: : GLM MCMC:? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 MCMC MCMC: MCMC MCMC: : MCMC () Markov chain Monte Carlo (MCMC) Metropolis Method (Metropolis method) :?? MCMC Metropolis Method kubo ( ( 2 ) / 90 MCMC log likelihood log L(q) * : q log likelihood q () kubo ( ( 2 ) / 90

22 kubo2017sep16b p.4 MCMC MCMC: MCMC MCMC: q log likelihood q: seed sruvivorship 1 q 2 q 3 (1), (2) kubo ( ( 2 ) / 90 Metropolis Method : 1 q ( q 0.3) 2 q ( q q new ) 3 q new L(q new ) L(q) L(q new ) L(q) (): q q new L(q new ) < L(q) (): r = L(q new)/l(q) q qnew 1 r q 4 2. (q = 0.01 q = 0.99 ) kubo ( ( 2 ) / 90 MCMC MCMC: MCMC MCMC: Metropolis Method q log likelihood log likelihood (MCMC) Metropolis Method kubo ( ( 2 ) / 90 q Metropolis Method ( MCMC) log likelihood ? q kubo ( ( 2 ) / 90 MCMC MCMC: MCMC MCMC: q q? q q? q MCMC MCMC?? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

23 kubo2017sep16b p.5 MCMC MCMC: MCMC MCMC: q MCMC q q? kubo ( ( 2 ) / 90 MCMC? log L(q) log likelihood * L(q) MCMC () kubo ( ( 2 ) / 90 MCMC MCMC: MCMC MCMC: MCMC q () MCMC p(q) q 95%! kubo ( ( 2 ) / 90 q Metropolis Method p(q Y) L(q) p(q) q kubo ( ( 2 ) / 90 MCMC MCMC MCMC? 2. MCMC Gibbs sampling 1 : : MCMC 2 R package : package : 3 BUGS Gibbs sampler : Gibbs sampler? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

24 kubo2017sep16b p.6 MCMC のためのソフトウェア MCMC のためのソフトウェアさまざまな MCMC アルゴリズム Gibbs sampling とは何か? MCMC アルゴリズムのひとついろいろな MCMC Metropolis Method: 試行錯誤で値を変化させていく MCMC Metropolis-Hastings: その改良版複数のパラメーターの MCMC サンプリングに使う例: パラメーター β1 と β2 の Gibbs sampling 1 2 Gibbs sampling: 条件つき確率分布を使った MCMC β2 に何か適当な値を与える β2 の値はそのままにしてその条件のもとでの β1 の MCMC sampling をする (条件つき事後分布) 3 複数の変数 (パラメーター状態) を効率よくサンプリング an β1 の値はそのままにしてその条件のもとでの β2 の MCMC sampling をする (条件つき事後分布) eﬃcient method for sampling of parameter values をくりかえす教科書の第 9 章の例題で説明 kubo ( 統計モデリング入門 (第 2 部) / 90 kubo ( MCMC のためのソフトウェア図解: Gibbs sampling β2 のサンプリング BUGS 言語 (+ っぽいもの) でベイズモデルを 8 10 記述できるソフトウェア β β β step β WinBUGS 歴史を変えてさようなら? OpenBUGS 予算が足りなくて停滞? JAGS お手軽で良いどんな OS でも動く Stan いま一番の注目今日は紹介しませんが 8 10 リンク集: β1 3 kubo ( step 3 32 / 90 便利な BUGS 汎用 Gibbs sampler たち (統計モデリング入門の第 9 章) step MCMC のためのソフトウェア β1 のサンプリング MCMC 統計モデリング入門 (第 2 部) β 統計モデリング入門 (第 2 部) えーと BUGS 言語って何? / 90 MCMC のためのソフトウェア kubo ( 統計モデリング入門 (第 2 部) / 90 MCMC のためのソフトウェアいろいろな OS で使える JAGS このベイズモデルを BUGS 言語で記述したいデータ Y[i] 種子数８個のうちの生存数二項分布 dbin(q,8) BUGS 言語コード R core team のひとり Martyn Plummer さんが開発 Just Another Gibbs Sampler C++ で実装されている R がインストールされていることが必要生存確率 q Linux, Windows, Mac OS X バイナリ版もある無情報事前分布開発進行中 R からの使う: library(rjags) 矢印は手順ではなく依存関係をあらわしている BUGS 言語: ベイズモデルを記述する言語 Spiegelhalter et al BUGS: Bayesian Using Gibbs Sampling version kubo ( 統計モデリング入門 (第 2 部) / 90 kubo ( 統計モデリング入門 (第 2 部) / 90

25 kubo2017sep16b p.7 MCMC JAGS R モデルの構造データとパラメーターの初期値サンプリングの詳細 Input BUGS 言語 JAGS MCMC sampling from posterior distributions 事後分布からのランダムサンプル Trace of beta[1] Iterations Trace of beta[2] Iterations Output Density of beta[1] N = 3000 Bandwidth = Density of beta[2] N = 3000 Bandwidth = kubo ( ( 2 ) / 90 MCMC R JAGS (1 / 3) library(rjags) library(r2winbugs) # to use write.model() model.bugs <- function() { for (i in 1:N.data) { Y[i] ~ dbin(q, 8) # } q ~ dunif(0.0, 1.0) # q } file.model <- "model.bug.txt" write.model(model.bugs, file.model) # # kubo ( ( 2 ) / 90 MCMC R JAGS (2 / 3) load("mcmc.rdata") # (data.rdata mcmc.rdata!!) list.data <- list(y = data, N.data = length(data)) inits <- list(q = 0.5) n.burnin < n.chain <- 3 n.thin <- 1 n.iter <- n.thin * 1000 model <- jags.model( file = file.model, data = list.data, inits = inits, n.chain = n.chain ) # kubo ( ( 2 ) / 90 MCMC R JAGS (3 / 3) # burn-in update(model, n.burnin) # burn in # post.mcmc.list post.mcmc.list <- coda.samples( model = model, variable.names = names(inits), n.iter = n.iter, thin = n.thin ) # kubo ( ( 2 ) / 90 MCMC burn in? MCMC step kubo ( ( 2 ) / 90 MCMC MCMC ˆR = MCMC ˆR = MCMC MCMC step! kubo ( ( 2 ) / 90

26 kubo2017sep16b p.8 MCMC ˆR MCMC Gibbs sampling plot(post.mcmc.list) gelman.diag(post.mcmc.list) R-hat Gelman-Rubin ˆR var ˆ + (ψ y) = W var ˆ + (ψ y) = n 1 n W + 1 n B W : variance B : variance Gelman et al Bayesian Data Analysis. Chapman & Hall/CRC Trace of q Density of q Iterations N = 1000 Bandwidth = kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90! GLM ! y i? ( 10 ) kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 (overdispersion) y i 0.5 : overdispersion : Modeling of individual difference i N i y i ) p(y i q i ) = ( Ni y i q y i i (1 q i) N i y i, q i kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

27 kubo2017sep16b p.9 Survivorship: logistic regression (GLM) q i = q(z i ) q(z) = 1/{1 + exp( z)} q(z) z i = a + r i a: r i : i () z r i > (a {r 1, r 2,, r 100 }) /! () kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 {r i } s = 1.0 s = 1.5 s = 3.0 r i ) 1 p(r i s) = ( exp r2 i 2πs 2 2s 2 p(r i s) r i r i r i kubo ( ( 2 ) / 90 : r i (A) I III I II I I y i (B) p(r i s) s = {r i} s = 3.0 r i 1 q i = 1+exp( r i) II I II I II I I I III I II I I I I I p(y i q i) y i kubo ( ( 2 ) / 90 r i r i {r i } 100 r i s = 1.0 s = 1.5 r i p(r i s) = s = πs 2 exp ( r2 i 2s 2 ) (A) (!) (B) (C) s kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

28 kubo2017sep16b p.10 r i! s = 1.0 s = 1.5 r i s = 3.0 p(r i s) = 1 2πs 2 exp ( r2 i 2s 2 kubo ( ( 2 ) / 90 ) 全データ個体個体 3 3 のデータのデータ個体 1 のデータ個体個体 3 3 のデータのデータ個体 2 のデータ {r 1, r 2, r 3,..., r 100 } s local parameter random effects a global parameter fixed effects? kubo ( ( 2 ) / 90 {r i } s (B) (C) s = 1.0 a, s {r i } s s = 1.5 s = 3.0 s s (non-informative prior) 0 < s < 10 4 kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 a : Hierarchical and non-informative priors ( 0; 1) ( 0; 100) (logit) a a 種子 8 個のうち Y[i] が生存生存 q[i] 全個体共通の平均 a r[i] s hyper parameter kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

29 kubo2017sep16b p.11 JAGS R JAGS BUGS model { for (i in 1:N.data) { Y[i] ~ dbin(q[i], 8) logit(q[i]) <- a + r[i] } a ~ dnorm(0, 1.0E-4) for (i in 1:N.data) { r[i] ~ dnorm(0, tau) } tau <- 1 / (s * s) s ~ dunif(0, 1.0E+4) } 種子 8 個のうち Y[i] が生存生存 q[i] 全個体共通の平均 a r[i] s hyper parameter kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 JAGS > source("mcmc.list2bugs.r") # > post.bugs <- mcmc.list2bugs(post.mcmc.list) # bugs 80% 10 interval 5 0 for each 5 10 chain 1 R hat a * r[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] s [76] * array truncated for lack of space 3 chains, each with 4000 iterations (first 2000 discarded) a * r s medians and 80% intervals kubo ( ( 2 ) / 90 bugs post.bugs print(post.bugs, digits.summary = 3) 95% 3 chains, each with 4000 iterations (first 2000 discarded), n.thin = 2 n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a s r[1] r[2] r[3] r[4] r[5] () r[99] r[100] kubo ( ( 2 ) / 90 R Trace of a Density of a post.mcmc <- to.mcmc(post.bugs) Iterations Trace of s Iterations N = 1000 Bandwidth = Density of s N = 1000 Bandwidth = matrix kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

30 kubo2017sep16b p.12 GLMM? 線形モデルの発展推定計算方法階層ベイズモデル (HBM) MCMC もっと自由な統計モデリングを! 一般化線形混合モデル (GLMM) 最尤推定法個体差場所差といった変量効果をあつかいたい一般化線形モデル (GLM) 正規分布以外の確率分布をあつかいたい最小二乗法線形モデル (Generalized Linear Mixed Model) 4. + GLMM local parameter GLMM kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 :? pot G pot E pot B pot H pot C pot I pot A pot D pot J pot F y i ( 100 ) (f j = C) 5 ( 50 ) (f j = T) 5 ( 50 ) kubo ( ( 2 ) / 90 pot A pot G pot J pot H pot E pot D pot C pot I pot F pot B y i ( 100 ) (f j = C) 5 ( 50 ) (f j = T) 5 ( 50 ) kubo ( ( 2 ) / 90!! I > d <- read.csv("d1.csv") > head(d) id pot f y id : {1, 2, 3,, 100} 1 1 A C 6 pot : {A, B, C, 2 2 A C 3, J} 3 3 A C 19 f : : C, 4 4 A C 5 T 5 5 A C 0 y : () 6 6 A C 19 d$y D D J B AA D G I D E I A G II A AA A A B DE E G C D E H B C A B C B C CC D EE D E EF F F F I G H H I II J J JJ J d$id plot(d$id, d$y, pch = as.character(d$pot),...)? kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90

31 kubo2017sep16b p.13? I d$y D D J B AA D G I D E I A G II A AA A A B DE E G C D E H B C A B C B C CC D EE D E EF F F F I G H H I II J J JJ J d$id A B C D E F G H I J? () plot(d$pot, d$y, col = rep(c("blue", "red"), each = 5)) random effects kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 () 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル GLM: > summary(glm(y ~ f, data = d, family = poisson))...()... Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) < 2e-16 ft e-06...()... (f)? AIC kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 GLMM:! > library(glmmml) > summary(glmmml(y ~ f, data = d, family = poisson, + cluster = id))...()... coef se(coef) z Pr(> z ) (Intercept) e-12 ft e-03...()...?? kubo ( ( 2 ) / 90,! kubo ( ( 2 ) / 90

32 kubo2017sep16b p.14 R! Use JAGS! + R GLMM library(glmmml) glmmml() library(lme4) lmer() library(nlme) nlme() () GLMM + () log log(λ i ) = a + bf i + () + () a f i b () ( σ 1, σ 2 ) σ ([0, 10 4 ] ) kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 データ各個体の種子数 Y[i] 個ポアソン分布平均 lambda[i] 全個体共通 beta[1] 無情報事前分布データ ( 説明変数 ) 施肥処理 F[i] 植木鉢差 rp[pot[i]] 個体差 r[i] 階層事前分布階層事前分布 s[1] 個体差の s[2] 植木鉢差の beta[2] ばらつきばらつき無情報事前分布 ( 超事前分布 ) kubo ( ( 2 ) / 90 + BUGS code (1) model { for (i in 1:N.sample) { Y[i] ~ dpois(lambda[i]) log(lambda[i]) <- a + b * F[i] + r[i] + rp[pot[i]] } # BUGS coding f i {C, T} F[i] 0, 1 Pot[i] 1, 2,..., 10 rp[...] kubo ( ( 2 ) / 90 + BUGS code (2) # a ~ dnorm(0, 1.0E-4) # b ~ dnorm(0, 1.0E-4) # for (i in 1:N.sample) { r[i] ~ dnorm(0, tau[1]) # } for (j in 1:N.pot) { rp[j] ~ dnorm(0, tau[2]) # () } for (k in 1:N.tau) { tau[k] <- 1.0 / (sigma[k] * sigma[k]) # sigma[k] ~ dunif(0, 1.0E+4) } } kubo ( ( 2 ) / 90 ( b)? mean sd 2.5% 25% 50% 75% 97.5% Rhat n. a b sigma[1] ()... kubo ( ( 2 ) / 90

33 kubo2017sep16b p.15 () rp[j]! random effects random effects fixed effects! kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 : (GLM) 3. MCMC 4. GLM kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 kubo ( ( 2 ) / 90 階層ベイズモデルもっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル個体差場所差といった変量効果をあつかいたい一般化線形モデル正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル () kubo ( ( 2 ) / 90

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht kubo2015ngt6 p.1 2015 (6 MCMC kubo@ees.hokudai.ac.jp, @KuboBook http://goo.gl/m8hsbm 1 ( 2 3 4 5 JAGS : 2015 05 18 16:48 kubo (http://goo.gl/m8hsbm 2015 (6 1 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 2 /